text:
# differentiated version of shock sensors

#------------------------------------------------------------------------------
# ShockSensorPP

"""
  Differentiated versionn of `getShockSensor` for [`ShockSensorPP`](@ref)
"""
function getShockSensor_diff(params::ParamType{Tdim}, sbp::AbstractOperator,
                      sensor::ShockSensorPP,
                      q::AbstractMatrix{Tsol}, elnum::Integer,
                      coords::AbstractMatrix, dxidx::Abstract3DArray,
                      jac::AbstractVector{Tmsh},
                      Se_jac::Abstract4DArray{Tres},
                      ee_jac::Abstract4DArray{Tres}
                     ) where {Tsol, Tmsh, Tres, Tdim}
# computes the Jacobian of Se and ee wrt q. Does not take in q_dot because
# that would be way more expensive
# Se_jac and ee_jac are overwritten
# The third output argument tells if ee_jac is all zeros (hopefully a common
# case).

  numDofPerNode, numNodesPerElement = size(q)
  dim = size(coords, 1)

  @unpack sensor up up_tilde up1_tilde s0 kappa e0 num_dot den_dot ee_dot lambda_max_dot
  fill!(num_dot, 0); fill!(den_dot, 0); fill!(Se_jac, 0); fill!(ee_jac, 0)

  @simd for i=1:numNodesPerElement
    up[i] = q[1, i]
  end

  # for getFiltered solution, the matrix itself is the Jacobian
  getFilteredSolutions(params, sensor.Vp, up, up_tilde, up1_tilde)
#  getFilteredSolution(params, sensor.Vp, up, up_tilde)
#  getFilteredSolution(params, sensor.Vp1, up, up1_tilde)
  up_tilde_dotT = sensor.Vp.filtT  # use transposed because of memory order
  up1_tilde_dotT = sensor.Vp.filt1T
  #up1_tilde_dotT = sensor.Vp1.filtT
  #up1_tilde_dotT = sensor.Vp1.filt1T

  # compute the inner product
  num = zero(Tres)
  den = zero(Tres)

  if sensor.use_filtered
    @simd for i=1:numNodesPerElement
      fac = sbp.w[i]/jac[i]
      delta_u = up_tilde[i] - up1_tilde[i]

      num += delta_u*fac*delta_u
      @simd for j=1:numNodesPerElement
        delta_u_dot = up_tilde_dotT[j, i] - up1_tilde_dotT[j, i]
        num_dot[j] += 2*fac*delta_u*delta_u_dot
      end

      den += up_tilde[i]*fac*up_tilde[i]
      @simd for j=1:numNodesPerElement
        den_dot[j] += 2*fac*up_tilde[i]*up_tilde_dotT[j, i]
      end
    end
  else
    @simd for i=1:numNodesPerElement
      fac = sbp.w[i]/jac[i]
      delta_u = up[i] - up1_tilde[i]

      num += delta_u*fac*delta_u

      @simd for j=1:numNodesPerElement
        delta_u_dot = -up1_tilde_dotT[j, i]
        if j == i
          delta_u_dot += 1
        end

        num_dot[j] += 2*fac*delta_u*delta_u_dot
      end

      den += up[i]*fac*up[i]
      den_dot[i] += 2*fac*up[i]
    end
  end

  Se = num/den
  fac2 = 1/(den*den)
#  fac3 = 1/(2*Se)
  @simd for i=1:numNodesPerElement
    Se_jac[1, 1, 1, i] = (num_dot[i]*den - den_dot[i]*num)*fac2#*fac3
  end

  se = log10(Se)
  eejac_zero = true
  ee_dot = zeros(Tres, numNodesPerElement)

  if se < s0 - kappa
    #ee = zero(Tres)
    ee = Tres(0.0)
    fill!(ee_dot, 0)
  elseif se > s0 - kappa && se < s0 + kappa
    ee = 0.5*e0*(1 + sinpi( (se - s0)/(2*kappa)))

    # derivative of ee wrt Se (not se)
    fac3 = 0.5*e0*cospi( (se - s0)/(2*kappa) ) * (Float64(pi)/(2*kappa*log(10)*Se))
    fill!(ee_dot, 0)
    @simd for i=1:numNodesPerElement
      ee_dot[i] = fac3*Se_jac[1, 1, 1, i]
    end
    eejac_zero = false
  else
    ee = Tres(e0)
    fill!(ee_dot, 0)
    eejac_zero = false
  end

  # multiply by lambda_max * h/p to get subcell resolution
  lambda_max = zero(Tsol)
  lambda_max_dot = zeros(Tres, numDofPerNode, numNodesPerElement)  #TODO: use preallocated version
  h_avg = zero(Tmsh)
  @simd for i=1:numNodesPerElement
    q_i = sview(q, :, i)
    lambda_max_dot_i = sview(lambda_max_dot, :, i)
    lambda_max += getLambdaMax_diff(params, q_i, lambda_max_dot_i)
    h_avg += sbp.w[i]/jac[i]
  end

  lambda_max /= numNodesPerElement
  scale!(lambda_max_dot, 1/numNodesPerElement)
  h_avg = h_avg^(1/Tdim)

  fill!(ee_jac, 0)
  @simd for i=1:numNodesPerElement
    ee_jac[1, 1, 1, i] = lambda_max*h_avg*ee_dot[i]/sbp.degree
    @simd for j=1:numDofPerNode
      ee_jac[1, j, 1, i] += lambda_max_dot[j, i]*ee*h_avg/sbp.degree
    end
  end
  ee *= lambda_max*h_avg/sbp.degree

  # fill in the rest of the jac arrays
  @simd for q=1:numNodesPerElement
    @simd for p=1:numNodesPerElement
      @simd for j=1:numDofPerNode
        @simd for i=1:dim
          Se_jac[i, j, p, q] = Se_jac[1, j, 1, q]
          ee_jac[i, j, p, q] = ee_jac[1, j, 1, q]
        end
      end
    end
  end

  return eejac_zero
end

function getShockSensor_revq(params::ParamType{Tdim}, sbp::AbstractOperator,
                        sensor::ShockSensorPP{Tsol, Tres},
                        q::AbstractMatrix, q_bar::AbstractMatrix,
                        elnum::Integer,
                        coords::AbstractMatrix,
                        dxidx::Abstract3DArray, jac::AbstractVector{Tmsh},
                        ee_mat::AbstractMatrix, ee_mat_bar::AbstractMatrix
                        ) where {Tsol, Tres, Tmsh, Tdim}

  numNodesPerElement = size(q, 2)

  @unpack sensor up up_tilde up1_tilde s0 kappa e0 up_tilde_bar up1_tilde_bar
  @unpack sensor up_bar
  fill!(up_tilde_bar, 0); fill!(up1_tilde_bar, 0); fill!(up_bar, 0)

  # use density as the key variable
  for i=1:numNodesPerElement
    up[i] = q[1, i]
  end

  getFilteredSolutions(params, sensor.Vp, up, up_tilde, up1_tilde)
  #getFilteredSolution(params, sensor.Vp, up, up_tilde)
  #getFilteredSolution(params, sensor.Vp1, up, up1_tilde)

  # compute the inner product
  num = zero(Tres)
  den = zero(Tres)
  if sensor.use_filtered
    for i=1:numNodesPerElement
      fac = sbp.w[i]/jac[i]
      delta_u = up_tilde[i] - up1_tilde[i]

      num += delta_u*fac*delta_u
      den += up_tilde[i]*fac*up_tilde[i]
    end
  else
    for i=1:numNodesPerElement
      fac = sbp.w[i]/jac[i]
      delta_u = up[i] - up1_tilde[i]

      num += delta_u*fac*delta_u
      den += up[i]*fac*up[i]
    end
  end

  Se = num/den
  se = log10(Se)

  # should this be a separate function from computing Se?
  if se < s0 - kappa
    ee = Tres(0.0)
    is_nonzero = false
  elseif se > s0 - kappa && se < s0 + kappa
    ee = 0.5*e0*(1 + sinpi( (se - s0)/(2*kappa)))
    is_nonzero = true
  else
    ee = Tres(e0)
    is_nonzero = true
  end

  # scale by lambda_max * h/p to get subcell resolution
  lambda_max = zero(Tsol)
  h_avg = zero(Tmsh)
  for i=1:numNodesPerElement
    q_i = sview(q, :, i)
    lambda_max += getLambdaMax(params, q_i)
    h_avg += sbp.w[i]/jac[i]
  end

  lambda_max /= numNodesPerElement
  h_avg = h_avg^(1/Tdim)

  ee_final = ee*lambda_max*h_avg/sbp.degree

  # use elementwise constant
  fill!(ee_mat, ee)

  #---------------------
  # reverse sweep
  ee_final_bar = zero(Tres)
  for i=1:numNodesPerElement
    for d=1:Tdim
      ee_final_bar += ee_mat_bar[d, i]
    end
  end

  ee_bar = zero(Tres); lambda_max_bar = zero(Tres)
  ee_bar         += ee_final_bar*lambda_max*h_avg/sbp.degree
  lambda_max_bar += ee_final_bar*ee*h_avg/sbp.degree

  lambda_max_bar /= numNodesPerElement
  for i=1:numNodesPerElement
    q_i = ro_sview(q, :, i)
    q_bar_i = sview(q_bar, :, i)
    getLambdaMax_revq(params, q_i, q_bar_i, lambda_max_bar)
  end

  # if s < s0 - kappa or > s0 + kappa, then ee is either 1 or 0, so there is
  # no q dependence and no need for further computation
  if se > s0 - kappa && se < s0 + kappa

    # derivative of ee wrt Se (not se)
    Se_bar = ee_bar*0.5*e0*cospi( (se - s0)/(2*kappa) ) * (Float64(pi)/(2*kappa*log(10)*Se))

    num_bar = Se_bar/den
    den_bar = -num*Se_bar/(den*den)

#    up_tilde_bar = zeros(Tsol, numNodesPerElement)
#    up1_tilde_bar = zeros(Tsol, numNodesPerElement)a

    if sensor.use_filtered
      for i=1:numNodesPerElement
        fac = sbp.w[i]/jac[i]
        delta_u = up_tilde[i] - up1_tilde[i]

        delta_u_bar       = 2*delta_u*fac*num_bar
        up_tilde_bar[i]  += 2*up_tilde[i]*fac*den_bar
        up_tilde_bar[i]  += delta_u_bar
        up1_tilde_bar[i] -= delta_u_bar
      end
    else
      for i=1:numNodesPerElement
        fac = sbp.w[i]/jac[i]
        delta_u = up[i] - up1_tilde[i]

        delta_u_bar       = 2*delta_u*fac*num_bar
        up_bar[i]        += 2*up[i]*fac*den_bar
        up_bar[i]        += delta_u_bar
        up1_tilde_bar[i] -= delta_u_bar
      end
    end

    getFilteredSolutions_rev(params, sensor.Vp, up_bar, up_tilde_bar, up1_tilde_bar)

    for i=1:numNodesPerElement
      q_bar[1, i] += up_bar[i]
    end
  end

  return is_nonzero
end

function getShockSensor_revm(params::ParamType{Tdim}, sbp::AbstractOperator,
                        sensor::ShockSensorPP{Tsol, Tres},
                        q::AbstractMatrix,
                        elnum::Integer,
                        coords::AbstractMatrix, coords_bar::AbstractMatrix,
                        dxidx::Abstract3DArray, dxidx_bar::Abstract3DArray,
                        jac::AbstractVector{Tmsh}, jac_bar::AbstractVector,
                        ee_mat_bar::AbstractMatrix
                       ) where {Tsol, Tres, Tmsh, Tdim}

  numNodesPerElement = size(q, 2)

  @unpack sensor up up_tilde up1_tilde s0 kappa e0
  #fill!(up_tilde, 0); fill!(up1_tilde, 0)

  # use density as the key variable
  for i=1:numNodesPerElement
    up[i] = q[1, i]
  end

  getFilteredSolutions(params, sensor.Vp, up, up_tilde, up1_tilde)
  #getFilteredSolution(params, sensor.Vp, up, up_tilde)
  #getFilteredSolution(params, sensor.Vp1, up, up1_tilde)

  # compute the inner product
  num = zero(Tres)
  den = zero(Tres)
  if sensor.use_filtered
    for i=1:numNodesPerElement
      fac = sbp.w[i]/jac[i]
      delta_u = up_tilde[i] - up1_tilde[i]

      num += delta_u*fac*delta_u
      den += up_tilde[i]*fac*up_tilde[i]

      # see if this makes the sensor less sensitive
      #den += up[i]*fac*up[i]
    end
  else
    for i=1:numNodesPerElement
      fac = sbp.w[i]/jac[i]
      delta_u = up[i] - up1_tilde[i]

      num += delta_u*fac*delta_u
      den += up[i]*fac*up[i]
    end

  end

  Se = num/den
  se = log10(Se)

  # should this be a separate function from computing Se?
  if se < s0 - kappa
    ee = Tres(0.0)
    is_nonzero = false
  elseif se > s0 - kappa && se < s0 + kappa
    ee = 0.5*e0*(1 + sinpi( (se - s0)/(2*kappa)))
    is_nonzero = true
  else
    ee = Tres(e0)
    is_nonzero = true
  end

  # scale by lambda_max * h/p to get subcell resolution
  lambda_max = zero(Tsol)
  h_avg = zero(Tmsh)
  for i=1:numNodesPerElement
    q_i = sview(q, :, i)
    lambda_max += getLambdaMax(params, q_i)
    h_avg += sbp.w[i]/jac[i]
  end

  lambda_max /= numNodesPerElement
  h_avg2 = h_avg^(1/Tdim)

  ee_final = ee*lambda_max*h_avg2/sbp.degree

  # use elementwise constant
  #fill!(ee_mat, ee)

  #---------------------
  # reverse sweep
  ee_final_bar = zero(Tres)
  for i=1:numNodesPerElement
    for d=1:Tdim
      ee_final_bar += ee_mat_bar[d, i]
    end
  end

  ee_bar = zero(Tres); lambda_max_bar = zero(Tres); h_avg2_bar = zero(Tres)
  ee_bar         += ee_final_bar*lambda_max*h_avg2/sbp.degree
  h_avg2_bar += ee*lambda_max*ee_final_bar/sbp.degree

  h_avg_bar = (1/Tdim)*h_avg2_bar*h_avg^((1/Tdim) - 1)

  for i=1:numNodesPerElement
    jac_bar[i] += -sbp.w[i]*h_avg_bar/(jac[i]*jac[i])
  end
#=
  lambda_max_bar /= numNodesPerElement
  for i=1:numNodesPerElement
    q_i = ro_sview(q, :, i)
    q_bar_i = sview(q_bar, :, i)
    getLambdaMax_revq(params, q_i, q_bar_i, lambda_max_bar)
  end
=#
  # if s < s0 - kappa or > s0 + kappa, then ee is either 1 or 0, so there is
  # no metric dependence and no need for further computation
  if se > s0 - kappa && se < s0 + kappa

    # derivative of ee wrt Se (not se)
    Se_bar = ee_bar*0.5*e0*cospi( (se - s0)/(2*kappa) ) * (Float64(pi)/(2*kappa*log(10)*Se))

    num_bar = Se_bar/den
    den_bar = -num*Se_bar/(den*den)

    #TODO: why didn't this error?
    if sensor.use_filtered
      for i=1:numNodesPerElement
        fac = sbp.w[i]/jac[i]
        delta_u = up_tilde[i] - up1_tilde[i]

        fac_bar     = delta_u*delta_u*num_bar
        fac_bar    += up_tilde[i]*up_tilde[i]*den_bar
        jac_bar[i] += -sbp.w[i]*fac_bar/(jac[i]*jac[i])
      end
    else
      for i=1:numNodesPerElement
        fac = sbp.w[i]/jac[i]
        delta_u = up[i] - up1_tilde[i]

        fac_bar     = delta_u*delta_u*num_bar
        fac_bar    += up[i]*up[i]*den_bar
        jac_bar[i] += -sbp.w[i]*fac_bar/(jac[i]*jac[i])
      end

    end
  end

  return is_nonzero
end

"""
  Computed the two filtered solutions required by the shock sensor (well, only
  one of them is required, the other one may or may not be useful)

  **Inputs**

   * params
   * vand: `VandermondeData`
   * u_nodal_bar:

  **Inputs/Outputs**

   * u_filt1_bar:
   * u_filt2_bar
"""
function getFilteredSolutions_rev(params::ParamType, vand::VandermondeData,
                              u_nodal, u_filt1, u_filt2)

  smallmatTvec_kernel!(vand.filt, u_filt1, u_nodal, 1, 1)
  smallmatTvec_kernel!(vand.filt1, u_filt2,  u_nodal, 1, 1)

  return nothing
end

#------------------------------------------------------------------------------
# ShockSensorNone

function getShockSensor_diff(params::ParamType, sbp::AbstractOperator,
                      sensor::ShockSensorNone,
                      q::AbstractMatrix{Tsol}, elnum::Integer,
                      coords::AbstractMatrix, dxidx::Abstract3DArray,
                      jac::AbstractVector{Tmsh},
                      Se_jac::Abstract4DArray{Tres},
                      ee_jac::Abstract4DArray{Tres}) where {Tsol, Tmsh, Tres}

  error("getShockSensor_diff called for ShockSensorNone: did you forget to specify the shock capturing scheme?")
end

#------------------------------------------------------------------------------
# ShockSensorEverywhere

function getShockSensor_diff(params::ParamType, sbp::AbstractOperator,
                      sensor::ShockSensorEverywhere,
                      q::AbstractMatrix{Tsol}, elnum::Integer,
                      coords::AbstractMatrix, dxidx::Abstract3DArray,
                      jac::AbstractVector{Tmsh},
                      Se_jac::Abstract4DArray{Tres},
                      ee_jac::Abstract4DArray{Tres}) where {Tsol, Tmsh, Tres}

  fill!(Se_jac, 0)
  fill!(ee_jac, 0)

  return true
end

function getShockSensor_revq(params::ParamType, sbp::AbstractOperator,
                        sensor::ShockSensorEverywhere{Tsol, Tres},
                        q::AbstractMatrix, q_bar::AbstractMatrix,
                        elnum::Integer,
                        coords::AbstractMatrix,
                        dxidx::Abstract3DArray, jac::AbstractVector{Tmsh},
                        ee_mat::AbstractMatrix, ee_bar::AbstractMatrix
                        ) where {Tsol, Tres, Tmsh}

  # ee is constant, so q_bar += 0*ee_bar

  return true
end

#------------------------------------------------------------------------------
# ShockSensorVelocity

function getShockSensor_diff(params::ParamType, sbp::AbstractOperator,
                      sensor::ShockSensorVelocity,
                      q::AbstractMatrix{Tsol}, elnum::Integer,
                      coords::AbstractMatrix, dxidx::Abstract3DArray,
                      jac::AbstractVector{Tmsh},
                      Se_jac::Abstract4DArray{Tres},
                      ee_jac::Abstract4DArray{Tres}) where {Tsol, Tmsh, Tres}

  numNodesPerElement = size(q, 2)
  dim = size(coords, 1)

  fill!(Se_jac, 0)
  fill!(ee_jac, 0)

  # derivative with momentum in direction d = d, all others zero
  alpha = sensor.alpha
  @simd for p=1:numNodesPerElement
    @simd for d=1:dim
      Se_jac[d, d+1, p, p] = d*alpha
      ee_jac[d, d+1, p, p] = d*alpha
    end
  end

  return false
end

function getShockSensor_revq(params::ParamType, sbp::AbstractOperator,
                        sensor::ShockSensorVelocity{Tsol, Tres},
                        q::AbstractMatrix, q_bar::AbstractMatrix,
                        elnum::Integer,
                        coords::AbstractMatrix,
                        dxidx::Abstract3DArray, jac::AbstractVector{Tmsh},
                        ee_mat::AbstractMatrix, ee_bar::AbstractMatrix
                        ) where {Tsol, Tres, Tmsh}

  # ee is constant, so q_bar += 0*ee_bar
  numNodesPerElement = size(q, 2)
  dim = size(coords, 1)

  numNodesPerElement = size(q, 2)
  dim = size(coords, 1)

  @simd for i=1:numNodesPerElement
    @simd for d=1:dim
      ee_mat[d, i] = d*q[d+1, i]
      q_bar[d+1, i] += ee_bar[d, i]*d
    end
  end

  return true
end

function getShockSensor_revm(params::ParamType{Tdim}, sbp::AbstractOperator,
                        sensor::ShockSensorVelocity{Tsol, Tres},
                        q::AbstractMatrix,
                        elnum::Integer,
                        coords::AbstractMatrix, coords_bar::AbstractMatrix,
                        dxidx::Abstract3DArray, dxidx_bar::Abstract3DArray,
                        jac::AbstractVector{Tmsh}, jac_bar::AbstractVector,
                        ee_bar::AbstractMatrix
                       ) where {Tsol, Tres, Tmsh, Tdim}

  # nothing to do here

  return nothing
end

#------------------------------------------------------------------------------
# ShockSensorHIso

function getShockSensor_diff(params::ParamType{Tdim}, sbp::AbstractOperator,
                      sensor::ShockSensorHIso,
                      q::AbstractMatrix{Tsol}, elnum::Integer,
                      coords::AbstractMatrix, dxidx::Abstract3DArray,
                      jac::AbstractVector{Tmsh},
                      Se_jac::Abstract4DArray{Tres},
                      ee_jac::Abstract4DArray{Tres}) where {Tsol, Tmsh, Tres, Tdim}

  # compute | div(F) |
  # Do this in Cartesian coordinates because it makes the differentiation easier
  numDofPerNode, numNodesPerElement = size(q)
  dim = size(coords, 1)
  @unpack sensor res res_jac

  _Se_jac = zeros(Tres, numDofPerNode, numNodesPerElement)
  computeStrongResidual_diff(params, sbp, sensor.strongdata,
                             q, coords, dxidx, jac, res, res_jac)
  val = computeL2Norm_diff(params, sbp, jac, res, res_jac, _Se_jac)

  h_avg = computeElementVolume(params, sbp, jac)
  h_fac = h_avg^((2 - sensor.beta)/Tdim)

  ee = Tres(sensor.C_eps*h_fac*val)
  @simd for q=1:numNodesPerElement
    @simd for j=1:numDofPerNode
      ee_jac[1, j, 1, q] = sensor.C_eps*h_fac*_Se_jac[j, q]
    end
  end

  # fill in the rest of the jac arrays
  @simd for q=1:numNodesPerElement
    @simd for p=1:numNodesPerElement
      @simd for j=1:numDofPerNode
        @simd for i=1:dim
          Se_jac[i, j, p, q] = _Se_jac[j, q]
          ee_jac[i, j, p, q] = ee_jac[1, j, 1, q]
        end
      end
    end
  end

  # calling | div(f) | as the shock sensor, which is somewhat arbitrary
  return false
end

"""
  Differentiated version

  **Inputs**

   * params
   * sbp
   * data: `StrongFormData` to hold temporary arrays
   * q: solution on the entire element, `numDofPerNode` x `numNodesPerElement`
   * dxidx: scaled mapping jacobian for entire element, `dim` x `dim` x
            `numNodesPerElement`
   * jac: mapping jacobian determinant for entire element, `numNodesPerElement`

  **Inputs/Outputs**

   * val_jac: jacobian of `val` wrt `q`, same size as `q`

  **Outputs**

   * val
"""
function computeStrongResidual_diff(params::ParamType{Tdim},
                      sbp::AbstractOperator,
                      data::StrongFormData{Tsol, Tres},
                      q::AbstractMatrix,
                      coords::AbstractMatrix, dxidx::Abstract3DArray,
                      jac::AbstractVector{Tmsh},
                      res::AbstractMatrix, res_jac::Abstract4DArray) where {Tsol, Tres, Tmsh, Tdim}

  numDofPerNode, numNodesPerElement = size(q)
  @unpack data flux nrm aux_vars work flux_jac Dx
  fill!(res, 0); fill!(nrm, 0); fill!(flux_jac, 0); fill!(res_jac, 0)

  @simd for i=1:numNodesPerElement
    q_i = sview(q, :, i)
    @simd for d=1:Tdim
      nrm[d] = 1
      flux_d = sview(flux, :, i, d)
      fluxjac_d = sview(flux_jac, :, :, d, i, i)

      calcEulerFlux(params, q_i, aux_vars, nrm, flux_d)
      calcEulerFlux_diff(params, q_i, aux_vars, nrm, fluxjac_d)

      nrm[d] = 0
    end
  end

  applyDx(sbp, flux, dxidx, jac, work, res)

  calcDx(sbp, dxidx, jac, Dx)
  applyOperatorJac(Dx, flux_jac, res_jac)

  return nothing
end

function computeStrongResidual_revq(params::ParamType{Tdim},
                               sbp::AbstractOperator,
                               data::StrongFormData{Tsol, Tres},
                               q::AbstractMatrix, q_bar::AbstractMatrix,
                               dxidx::Abstract3DArray,
                               jac::AbstractVector,
                               res_bar::AbstractMatrix
                              ) where {Tdim, Tsol, Tres}

  @unpack data flux nrm aux_vars work flux_bar
  fill!(nrm, 0); fill!(flux_bar, 0)

  numDofPerNode, numNodesPerElement = size(q)
#=
  for i=1:numNodesPerElement
    q_i = sview(q, :, i)
    for d=1:Tdim
      nrm[d] = 1
      flux_d = sview(flux, :, i, d)

      calcEulerFlux(params, q_i, aux_vars, nrm, flux_d)

      nrm[d] = 0
    end
  end

  applyDx(sbp, flux, dxidx, jac, work, res)
=#
  # reverse sweep
  applyDx_revq(sbp, flux_bar, dxidx, jac, work, res_bar)

  @simd for i=1:numNodesPerElement
    q_i = sview(q, :, i)
    q_bar_i = sview(q_bar, :, i)
    @simd for d=1:Tdim
      nrm[d] = 1
      flux_bar_d = ro_sview(flux_bar, :, i, d)

      calcEulerFlux_revq(params, q_i, q_bar_i, aux_vars, nrm, flux_bar_d)

      nrm[d] = 0
    end
  end

  return nothing
end

function computeStrongResidual_revm(params::ParamType{Tdim},
                             sbp::AbstractOperator,
                             data::StrongFormData{Tsol, Tres},
                             q::AbstractMatrix,
                             dxidx::Abstract3DArray, dxidx_bar::Abstract3DArray,
                             jac::AbstractVector, jac_bar::AbstractVector,
                             res_bar::AbstractMatrix
                            ) where {Tdim, Tsol, Tres}

  @unpack data flux nrm aux_vars work flux_bar
  fill!(nrm, 0); fill!(flux_bar, 0)

  numDofPerNode, numNodesPerElement = size(q)

  @simd for i=1:numNodesPerElement
    q_i = sview(q, :, i)
    @simd for d=1:Tdim
      nrm[d] = 1
      flux_d = sview(flux, :, i, d)

      calcEulerFlux(params, q_i, aux_vars, nrm, flux_d)

      nrm[d] = 0
    end
  end

  #applyDx(sbp, flux, dxidx, jac, work, res)

  # reverse sweep
  applyDx_revm(sbp, flux, flux_bar, dxidx, dxidx_bar, jac, jac_bar, work, res_bar)

  return nothing
end

"""
  **Inputs**

   * params
   * sbp
   * jac: mapping jacobian determinant for the element, `numNodesPerElement`
   * res: `numDofPerNode` x `numNodesPerElement` array to compute the norm
          of.
   * res_jac: 4D jacobian of `res` wrt `q`

  **Inputs/Outputs**

   * val_jac: array to be overwritten with jacobian of `val` wrt `q`

  **Outputs**

   * val: L2 norm
"""
function computeL2Norm_diff(params::ParamType{Tdim}, sbp::AbstractOperator,
                            jac::AbstractVector,
                            res::AbstractMatrix{Tres}, res_jac::Abstract4DArray,
                            val_jac::AbstractMatrix) where {Tres, Tdim}

  numDofPerNode, numNodesPerElement = size(res)

  # compute norm
  val = zero(Tres)
  fill!(val_jac, 0)
  @simd for i=1:numNodesPerElement
    fac = sbp.w[i]/jac[i]
    @simd for j=1:numDofPerNode
      val += res[j, i]*fac*res[j, i]
    end
  end

  val = sqrt(val)

  # jacobian of val (including the sqrt)
  fac_val = 1/(2*val)
  @simd for q=1:numNodesPerElement
    @simd for p=1:numNodesPerElement
      fac = sbp.w[p]/jac[p]
      @simd for j=1:size(res_jac, 2)
        @simd for i=1:size(res_jac, 1)
          val_jac[j, q] += fac_val*2*res[i, p]*fac*res_jac[i, j, p, q]
        end
      end
    end
  end

  return val
end

#------------------------------------------------------------------------------
# ShockSensorBO

function getShockSensor_diff(params::ParamType{Tdim}, sbp::AbstractOperator,
                      sensor::ShockSensorBO,
                      q::AbstractMatrix{Tsol}, elnum::Integer,
                      coords::AbstractMatrix, dxidx::Abstract3DArray,
                      jac::AbstractVector{Tmsh},
                      Se_jac::Abstract4DArray{Tres},
                      ee_jac::Abstract4DArray{Tres}) where {Tsol, Tmsh, Tres, Tdim}

  numDofPerNode, numNodesPerElement = size(q)
  dim = size(coords, 1)

  lambda_max = zero(Tsol)
  _Se_jac = zeros(Tres, numDofPerNode, numNodesPerElement)
  #lambda_max_dot = zeros(Tres, numDofPerNode, numNodesPerElement)
  h_avg = zero(Tmsh)
  @simd for i=1:numNodesPerElement
    q_i = sview(q, :, i)
    lambda_max_dot_i = sview(_Se_jac, :, i)
    lambda_max += getLambdaMax_diff(params, q_i, lambda_max_dot_i)
    h_avg += sbp.w[i]/jac[i]
  end

  lambda_max /= numNodesPerElement
  scale!(_Se_jac, 1/numNodesPerElement)
  h_avg = h_avg^(1/Tdim)

  @simd for i=1:numNodesPerElement
    @simd for j=1:numDofPerNode
      ee_jac[1, j, 1, i] = sensor.alpha*h_avg*_Se_jac[j, i]/sbp.degree
    end
  end

  # fill in the rest of the jac arrays
  @simd for q=1:numNodesPerElement
    @simd for p=1:numNodesPerElement
      @simd for j=1:numDofPerNode
        @simd for i=1:dim
          Se_jac[i, j, p, q] = _Se_jac[j, q]
          ee_jac[i, j, p, q] = ee_jac[1, j, 1, q]
        end
      end
    end
  end

  return false
end

function getShockSensor_revq(params::ParamType{Tdim}, sbp::AbstractOperator,
                        sensor::ShockSensorBO{Tsol, Tres},
                        q::AbstractMatrix, q_bar::AbstractMatrix,
                        elnum::Integer,
                        coords::AbstractMatrix,
                        dxidx::Abstract3DArray, jac::AbstractVector{Tmsh},
                        ee_mat::AbstractMatrix, ee_mat_bar::AbstractMatrix
                        ) where {Tsol, Tres, Tmsh, Tdim}

  numDofPerNode, numNodesPerElement = size(q)

#  lambda_max = zero(Tsol)
  h_avg = zero(Tmsh)
  for i=1:numNodesPerElement
    q_i = sview(q, :, i)
#    lambda_max += getLambdaMax(params, q_i)
    h_avg += sbp.w[i]/jac[i]
  end

#  lambda_max /= numNodesPerElement
  h_avg = h_avg^(1/Tdim)

#  ee = sensor.alpha*lambda_max*h_avg/sbp.degree
#  fill!(Se_mat, lambda_max)
#  fill!(ee_mat, sensor.alpha*lambda_max*h_avg/sbp.degree)

  #---------------------
  # reverse sweep
  ee_bar = zero(Tres)
  for i=1:numNodesPerElement
    for d=1:Tdim
      ee_bar += ee_mat_bar[d, i]
    end
  end

  lambda_max_bar = zero(Tres)
  lambda_max_bar += ee_bar*sensor.alpha*h_avg/sbp.degree

  lambda_max_bar /= numNodesPerElement
  for i=1:numNodesPerElement
    q_i = ro_sview(q, :, i)
    q_bar_i = sview(q_bar, :, i)
    getLambdaMax_revq(params, q_i, q_bar_i, lambda_max_bar)
  end

  return nothing
end

function getShockSensor_revm(params::ParamType{Tdim}, sbp::AbstractOperator,
                        sensor::ShockSensorBO{Tsol, Tres},
                        q::AbstractMatrix,
                        elnum::Integer,
                        coords::AbstractMatrix, coords_bar::AbstractMatrix,
                        dxidx::Abstract3DArray, dxidx_bar::Abstract3DArray,
                        jac::AbstractVector{Tmsh}, jac_bar::AbstractVector,
                        ee_mat_bar::AbstractMatrix
                       ) where {Tsol, Tres, Tmsh, Tdim}

 numDofPerNode, numNodesPerElement = size(q)

  lambda_max = zero(Tsol)
  h_avg = zero(Tmsh)
  for i=1:numNodesPerElement
    q_i = sview(q, :, i)
    lambda_max += getLambdaMax(params, q_i)
    h_avg += sbp.w[i]/jac[i]
  end

  lambda_max /= numNodesPerElement
  #h_avg2 = h_avg^(1/Tdim)

#  ee = sensor.alpha*lambda_max*h_avg/sbp.degree
#  fill!(Se_mat, lambda_max)
#  fill!(ee_mat, sensor.alpha*lambda_max*h_avg/sbp.degree)

  #---------------------
  # reverse sweep
  ee_bar = zero(Tres)
  for i=1:numNodesPerElement
    for d=1:Tdim
      ee_bar += ee_mat_bar[d, i]
    end
  end

  h_avg2_bar = zero(Tres)
  h_avg2_bar += sensor.alpha*lambda_max*ee_bar/sbp.degree

  h_avg_bar = (1/Tdim)*h_avg2_bar*h_avg^((1/Tdim) - 1)

  for i=1:numNodesPerElement
    jac_bar[i] += -sbp.w[i]*h_avg_bar/(jac[i]*jac[i])
  end

  return nothing
end

#------------------------------------------------------------------------------
# ShockSensorHHO

function getShockSensor_diff(params::ParamType{Tdim}, sbp::AbstractOperator,
                      sensor::ShockSensorHHO,
                      q_el::AbstractMatrix{Tsol}, elnum::Integer,
                      coords::AbstractMatrix, dxidx::Abstract3DArray,
                      jac::AbstractVector{Tmsh},
                      Se_jac::Abstract4DArray{Tres},
                      ee_jac::Abstract4DArray{Tres}) where {Tsol, Tmsh, Tres, Tdim}

  numDofPerNode, numNodesPerElement = size(q_el)
  dim = size(coords, 1)

  @unpack sensor p_jac press_el press_dx work res res_jac p_hess Dx px_jac val_dot Rp fp Rp_jac fp_jac

  fill!(Rp, 0)
  fill!(Rp_jac, 0); fill!(fp_jac, 0)
  fill!(press_dx, 0); fill!(res, 0); fill!(px_jac, 0)

  computeStrongResidual_diff(params, sbp, sensor.strongdata,
                             q_el, coords, dxidx, jac, res, res_jac)

  # compute pressure and pressure_jac @simd for entire element
  @simd for p=1:numNodesPerElement
    q_p = ro_sview(q_el, :, p)
    p_dot = sview(p_jac, 1, :, p)
    press_el[1, p] = calcPressure_diff(params, q_p, p_dot)
  end

  # compute Rp
  # jacobian first: p_dot/p * dR/dq + R*d/dq (p_dot/p)
  # first term

  @simd for q=1:numNodesPerElement
    @simd for p=1:numNodesPerElement
      press = press_el[1, p]
      p_dot = sview(p_jac, 1, :, p)
      @simd for j=1:numDofPerNode
        @simd for i=1:numDofPerNode
          #res_jac[i, j, p, q] *= p_dot[i]/press
          Rp_jac[p, j, q] += res_jac[i, j, p, q]*p_dot[i]/press
        end
      end
    end
  end

  # second term
  @simd for p=1:numNodesPerElement
    q_p = sview(q_el, :, p)
    press = press_el[1, p]
    p_dot = sview(p_jac, 1, :, p)
    calcPressure_hess(params, q_p, p_hess)
    @simd for j=1:numDofPerNode
      @simd for i=1:numDofPerNode
        Rp_jac[p, j, p] += res[i, p]*(-p_dot[i]*p_dot[j]/(press*press) + p_hess[i, j]/press)
      end
    end
  end

  # residual itself
  @simd for i=1:numNodesPerElement
    q_i = ro_sview(q_el, :, i)
    press = press_el[1, i]
    p_dot = sview(p_jac, 1, :, i)

    @simd for j=1:numDofPerNode
      Rp[i] += p_dot[j]*res[j, i]/press
    end
  end

  # compute pressure switch fp
  calcDx(sbp, dxidx, jac, Dx)
  applyOperatorJac(Dx, p_jac, px_jac)

  applyDx(sbp, press_el, dxidx, jac, work, press_dx)
  @simd for p=1:numNodesPerElement
    val = zero(Tres)
    fill!(val_dot, 0)
    @simd for d=1:Tdim
      val += press_dx[1, p, d]*press_dx[1, p, d]
      @simd for q=1:numNodesPerElement
        @simd for j=1:numDofPerNode
          val_dot[j, q] += 2*press_dx[1, p, d]*px_jac[1, j, d, p, q]
        end
      end
    end  # end d
    press = press_el[1, p]

    val2 = sqrt(val)
    scale!(val_dot, 1/(2*val2))

    t1 = (press + 1e-12)
    val3 = val2/t1

    # contribution of derivative of val2
    @simd for q=1:numNodesPerElement
      @simd for j=1:numDofPerNode
        val_dot[j, q] = val_dot[j, q]/t1
      end
    end

    # contribution of derivative of p
    @simd for j=1:numDofPerNode
      val_dot[j, p] += val2*(-p_jac[1, j, p]/(t1*t1))
    end

    fp[p] = val3
    @simd for q=1:numNodesPerElement
      @simd for j=1:numDofPerNode
        fp_jac[j, p, q] = val_dot[j, q]
      end
    end
  end

  @simd for q=1:numNodesPerElement
    @simd for p=1:numNodesPerElement
      @simd for j=1:numDofPerNode
        @simd for d=1:Tdim
          h_fac = sensor.h_k_tilde[d, elnum]
          Se_jac[d, j, p, q] = h_fac*(fp[p]*sign_c(Rp[p])*Rp_jac[p, j, q] +
                                      fp_jac[j, p, q]*absvalue(Rp[p]))
          ee_jac[d, j, p, q] = Se_jac[d, j, p, q]*h_fac*h_fac*sensor.C_eps
        end
      end
    end
  end

  return false
end

function getShockSensor_revq(params::ParamType{Tdim}, sbp::AbstractOperator,
                        sensor::ShockSensorHHO{Tsol, Tres},
                        q::AbstractMatrix, q_bar::AbstractMatrix,
                        elnum::Integer,
                        coords::AbstractMatrix,
                        dxidx::Abstract3DArray, jac::AbstractVector{Tmsh},
                        ee_mat::AbstractMatrix, ee_bar::AbstractMatrix
                       ) where {Tsol, Tres, Tmsh, Tdim}

  numDofPerNode, numNodesPerElement = size(q)

  @unpack sensor p_dot press_el press_dx work res Rp fp
  fill!(press_dx, 0); fill!(res, 0); fill!(Rp, 0)

  # compute Rm
  computeStrongResidual(params, sbp, sensor.strongdata, q, dxidx, jac, res)

  # compute Rp
  @simd for i=1:numNodesPerElement
    q_i = ro_sview(q, :, i)
    press = calcPressure_diff(params, q_i, p_dot)
    press_el[1, i] = press

    @simd for j=1:numDofPerNode
      Rp[i] += p_dot[j]*res[j, i]/press
    end
  end

  # compute pressure switch fp (excluding the h_k factor)
  applyDx(sbp, press_el, dxidx, jac, work, press_dx)
  @simd for i=1:numNodesPerElement
    val = zero(Tres)
    @simd for d=1:Tdim
      val += press_dx[1, i, d]*press_dx[1, i, d]
    end
    fp[i] = sqrt(val)/(press_el[1, i] + 1e-12)
  end

  @simd for i=1:numNodesPerElement
    @simd for d=1:Tdim
      h_fac = sensor.h_k_tilde[d, elnum]
      Se = h_fac*fp[i]*absvalue(Rp[i])
      ee_mat[d, i] = Se*h_fac*h_fac*sensor.C_eps
    end
  end

  # reverse sweep
  @unpack sensor fp_bar Rp_bar press_el_bar press_dx_bar p_dot_bar res_bar
  fill!(fp_bar, 0); fill!(Rp_bar, 0); fill!(press_el_bar, 0)
  fill!(press_dx_bar, 0); fill!(res_bar, 0)

  @simd for i=1:numNodesPerElement
    @simd for d=1:Tdim
      h_fac = sensor.h_k_tilde[d, elnum]
      Se_bar = h_fac*h_fac*sensor.C_eps*ee_bar[d, i]
      fp_bar[i] += h_fac*absvalue(Rp[i])*Se_bar
      Rp_bar[i] += h_fac*fp[i]*sign_c(Rp[i])*Se_bar
    end
  end

  # fp_bar
  @simd for i=1:numNodesPerElement
    val = zero(Tres)
    val_bar = zero(Tres)
    @simd for d=1:Tdim
      val += press_dx[1, i, d]*press_dx[1, i, d]
    end
    val_bar += fp_bar[i]/(2*sqrt(val)*(press_el[1, i] + 1e-12))
    press_el_bar[1, i] += -sqrt(val)*fp_bar[i]/(
                          (press_el[1, i] + 1e-12)*(press_el[1, i] + 1e-12) )
    @simd for d=1:Tdim
      press_dx_bar[1, i, d] += 2*press_dx[1, i, d]*val_bar
    end
  end
  applyDx_revq(sbp, press_el_bar, dxidx, jac, work, press_dx_bar)

  # Rp_bar
   @simd for i=1:numNodesPerElement
    q_i = ro_sview(q, :, i)
    q_bar_i = sview(q_bar, :, i)
    calcPressure_diff(params, q_i, p_dot)
    press = press_el[1, i]
    fill!(p_dot_bar, 0)

    @simd for j=1:numDofPerNode
      #Rp[i] += p_dot[j]*res[j, i]/press
      p_dot_bar[j]       += res[j, i]*Rp_bar[i]/press
      res_bar[j, i]      += p_dot[j]*Rp_bar[i]/press
      press_el_bar[1, i] += -p_dot[j]*res[j, i]*Rp_bar[i]/(press*press)
    end

    calcPressure_diff_revq(params, q_i, q_bar_i, p_dot_bar, press_el_bar[1, i])
  end

  computeStrongResidual_revq(params, sbp, sensor.strongdata, q, q_bar, dxidx,
                             jac, res_bar)

  return true
end

function getShockSensor_revm(params::ParamType{Tdim}, sbp::AbstractOperator,
                        sensor::ShockSensorHHO{Tsol, Tres},
                        q::AbstractMatrix,
                        elnum::Integer,
                        coords::AbstractMatrix, coords_bar::AbstractMatrix,
                        dxidx::Abstract3DArray, dxidx_bar::Abstract3DArray,
                        jac::AbstractVector{Tmsh}, jac_bar::AbstractVector,
                        ee_bar::AbstractMatrix
                       ) where {Tsol, Tres, Tmsh, Tdim}

  numDofPerNode, numNodesPerElement = size(q)

  @unpack sensor p_dot press_el press_dx work res Rp fp
  fill!(press_dx, 0); fill!(res, 0); fill!(Rp, 0)

  # compute Rm
  computeStrongResidual(params, sbp, sensor.strongdata, q, dxidx, jac, res)

  # compute Rp
  @simd for i=1:numNodesPerElement
    q_i = ro_sview(q, :, i)
    press = calcPressure_diff(params, q_i, p_dot)
    press_el[1, i] = press

    @simd for j=1:numDofPerNode
      Rp[i] += p_dot[j]*res[j, i]/press
    end
  end

  # compute pressure switch fp (excluding the h_k factor)
  applyDx(sbp, press_el, dxidx, jac, work, press_dx)
  @simd for i=1:numNodesPerElement
    val = zero(Tres)
    @simd for d=1:Tdim
      val += press_dx[1, i, d]*press_dx[1, i, d]
    end
    fp[i] = sqrt(val)/(press_el[1, i] + 1e-12)
  end
#=
  @simd for i=1:numNodesPerElement
    @simd for d=1:Tdim
      h_fac = sensor.h_k_tilde[d, elnum]
      Se = h_fac*fp[i]*absvalue(Rp[i])
      ee_mat[d, i] = Se*h_fac*h_fac*sensor.C_eps
    end
  end
=#

  # reverse sweep
  @unpack sensor fp_bar Rp_bar press_dx_bar p_dot_bar press_el_bar res_bar
  fill!(fp_bar, 0); fill!(Rp_bar, 0); fill!(press_el_bar, 0)
  fill!(press_dx_bar, 0); fill!(res_bar, 0)

  @simd for i=1:numNodesPerElement
    @simd for d=1:Tdim
      h_fac = sensor.h_k_tilde[d, elnum]
      Se_bar = h_fac*h_fac*sensor.C_eps*ee_bar[d, i]
      fp_bar[i] += h_fac*absvalue(Rp[i])*Se_bar
      Rp_bar[i] += h_fac*fp[i]*sign_c(Rp[i])*Se_bar
      sensor.h_k_tilde_bar[d, elnum] += 3*fp[i]*absvalue(Rp[i])*Se_bar
    end
  end

  # fp_bar
  @simd for i=1:numNodesPerElement
    val = zero(Tres)
    val_bar = zero(Tres)
    @simd for d=1:Tdim
      val += press_dx[1, i, d]*press_dx[1, i, d]
    end
    val_bar += fp_bar[i]/(2*sqrt(val)*(press_el[1, i] + 1e-12))
    #press_el_bar[1, i] += -sqrt(val)*fp_bar[i]/(
    #                      (press_el[1, i] + 1e-12)*(press_el[1, i] + 1e-12) )
    @simd for d=1:Tdim
      press_dx_bar[1, i, d] += 2*press_dx[1, i, d]*val_bar
    end
  end
  # press_el_bar unneeded, but required by applyDx_revm
  applyDx_revm(sbp, press_el, press_el_bar, dxidx, dxidx_bar, jac, jac_bar,
               work, press_dx_bar)

  # Rp_bar
   @simd for i=1:numNodesPerElement
    q_i = ro_sview(q, :, i)
    calcPressure_diff(params, q_i, p_dot)
    press = press_el[1, i]
    fill!(p_dot_bar, 0)

    @simd for j=1:numDofPerNode
      #Rp[i] += p_dot[j]*res[j, i]/press
      #p_dot_bar[j]       += res[j, i]*Rp_bar[i]/press
      res_bar[j, i]      += p_dot[j]*Rp_bar[i]/press
      #press_el_bar[1, i] += -p_dot[j]*res[j, i]*Rp_bar[i]/(press*press)
    end

    #calcPressure_diff_revq(params, q_i, q_bar_i, p_dot_bar, press_el_bar[1, i])
  end

  computeStrongResidual_revm(params, sbp, sensor.strongdata, q, dxidx,
                             dxidx_bar, jac, jac_bar, res_bar)

  return true
end

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

function getShockSensor_diff(params::ParamType{Tdim}, sbp::AbstractOperator,
                      sensor::ShockSensorHHOConst,
                      q_el::AbstractMatrix{Tsol}, elnum::Integer,
                      coords::AbstractMatrix, dxidx::Abstract3DArray,
                      jac::AbstractVector{Tmsh},
                      Se_jac::Abstract4DArray{Tres},
                      ee_jac::Abstract4DArray{Tres}) where {Tsol, Tmsh, Tres, Tdim}

  numDofPerNode, numNodesPerElement = size(q_el)
  dim = size(coords, 1)

  @unpack sensor p_jac press_el press_dx work res res_jac p_hess Dx px_jac val_dot Rp fp  epsilon Rp_jac fp_jac epsilon_dot val2_dot

  fill!(Rp, 0)
  fill!(Rp_jac, 0); fill!(fp_jac, 0)
  fill!(press_dx, 0); fill!(res, 0); fill!(px_jac, 0)

  computeStrongResidual_diff(params, sbp, sensor.strongdata,
                             q_el, coords, dxidx, jac, res, res_jac)

  # compute pressure and pressure_jac @simd for entire element
  @simd for p=1:numNodesPerElement
    q_p = ro_sview(q_el, :, p)
    p_dot = sview(p_jac, 1, :, p)
    press_el[1, p] = calcPressure_diff(params, q_p, p_dot)
  end

  # compute Rp
  # jacobian first: p_dot/p * dR/dq + R*d/dq (p_dot/p)
  # first term

  @simd for q=1:numNodesPerElement
    @simd for p=1:numNodesPerElement
      press = press_el[1, p]
      p_dot = sview(p_jac, 1, :, p)
      @simd for j=1:numDofPerNode
        @simd for i=1:numDofPerNode
          #res_jac[i, j, p, q] *= p_dot[i]/press
          Rp_jac[p, j, q] += res_jac[i, j, p, q]*p_dot[i]/press
        end
      end
    end
  end

  # second term
  @simd for p=1:numNodesPerElement
    q_p = sview(q_el, :, p)
    press = press_el[1, p]
    p_dot = sview(p_jac, 1, :, p)
    calcPressure_hess(params, q_p, p_hess)
    @simd for j=1:numDofPerNode
      @simd for i=1:numDofPerNode
        Rp_jac[p, j, p] += res[i, p]*(-p_dot[i]*p_dot[j]/(press*press) + p_hess[i, j]/press)
      end
    end
  end

  # residual itself
  @simd for i=1:numNodesPerElement
    q_i = ro_sview(q_el, :, i)
    press = press_el[1, i]
    p_dot = sview(p_jac, 1, :, i)

    @simd for j=1:numDofPerNode
      Rp[i] += p_dot[j]*res[j, i]/press
    end
  end

  # compute pressure switch fp
  calcDx(sbp, dxidx, jac, Dx)
  applyOperatorJac(Dx, p_jac, px_jac)

  applyDx(sbp, press_el, dxidx, jac, work, press_dx)
  @simd for p=1:numNodesPerElement
    val = zero(Tres)
    fill!(val_dot, 0)
    @simd for d=1:Tdim
      val += press_dx[1, p, d]*press_dx[1, p, d]
      @simd for q=1:numNodesPerElement
        @simd for j=1:numDofPerNode
          val_dot[j, q] += 2*press_dx[1, p, d]*px_jac[1, j, d, p, q]
        end
      end
    end  # end d
    press = press_el[1, p]

    val2 = sqrt(val)
    scale!(val_dot, 1/(2*val2))

    t1 = (press + 1e-12)
    val3 = val2/t1

    # contribution of derivative of val2
    @simd for q=1:numNodesPerElement
      @simd for j=1:numDofPerNode
        val_dot[j, q] = val_dot[j, q]/t1
      end
    end

    # contribution of derivative of p
    @simd for j=1:numDofPerNode
      val_dot[j, p] += val2*(-p_jac[1, j, p]/(t1*t1))
    end

    fp[p] = val3
    @simd for q=1:numNodesPerElement
      @simd for j=1:numDofPerNode
        fp_jac[j, p, q] = val_dot[j, q]
      end
    end
  end

  @simd for i=1:numNodesPerElement
    epsilon[1, i] = fp[i]*Rp[i]
  end

  @simd for q=1:numNodesPerElement
    @simd for p=1:numNodesPerElement
      @simd for j=1:numDofPerNode
        epsilon_dot[1, j, p, q] = fp[p]*Rp_jac[p, j, q] + Rp[p]*fp_jac[j, p, q]
      end
    end
  end

  #val = computeL2Norm_diff(params, sbp, jac, epsilon, epsilon_dot, val2_dot)
  fill!(val2_dot, 0)
  @simd for q=1:numNodesPerElement
    @simd for p=1:numNodesPerElement
      @simd for j=1:numDofPerNode
        val2_dot[j, q] += absvalue_deriv(epsilon[1, p])*sbp.w[p]*epsilon_dot[1, j, p, q]
      end
    end
  end

  @simd for q=1:numNodesPerElement
    @simd for p=1:numNodesPerElement
      @simd for j=1:numDofPerNode
        @simd for d=1:Tdim
          h_fac = sensor.h_k_tilde[d, elnum]
          Se_jac[d, j, p, q] = h_fac*val2_dot[j, q]
          ee_jac[d, j, p, q] = Se_jac[d, j, p, q]*h_fac*h_fac*sensor.C_eps
        end
      end
    end
  end

  return false
end

function calcAnisoFactors_revm(mesh::AbstractMesh, sbp, opts,
                          hk_all::AbstractMatrix{T}, hk_all_bar) where {T}

  # compute h_k tilde = h_k/(p + 1), where h_k takes into account the
  # anisotropy of the element in directory k
  # compute the p_k first, then compute h_k, where p_k is the integral of the
  # face normal
  @assert size(hk_all) == (mesh.dim, mesh.numEl)

  fill!(hk_all, 0)

  # do all faces (interior, boundary, shared)
  @simd for i=1:mesh.numInterfaces
    iface_i = mesh.interfaces[i]
    nrm_i = ro_sview(mesh.nrm_face, :, :, i)

    @simd for j=1:mesh.numNodesPerFace
      @simd for k=1:mesh.dim
        fac = mesh.sbpface.wface[j]*absvalue2(nrm_i[k, j])
        hk_all[k, iface_i.elementL] += fac
        hk_all[k, iface_i.elementR] += fac
      end
    end
  end

  @simd for i=1:mesh.numBoundaryFaces
    bndry_i = mesh.bndryfaces[i]
    nrm_i = ro_sview(mesh.nrm_bndry, :, :, i)

    @simd for j=1:mesh.numNodesPerFace
      @simd for k=1:mesh.dim
        fac = mesh.sbpface.wface[j]*absvalue2(nrm_i[k, j])
        hk_all[k, bndry_i.element] += fac
      end
    end
  end

  @simd for peer=1:mesh.npeers
    @simd for i=1:length(mesh.shared_interfaces[peer])
      iface_i = mesh.shared_interfaces[peer]
      nrm_i = ro_sview(mesh.nrm_sharedface[peer], :, :, i)

      @simd for j=1:mesh.numNodesPerFace
        @simd for k=1:mesh.dim
          fac = mesh.sbpface.wface[j]*absvalue2(nrm_i[k, j])
          hk_all[k, iface_i.elementL] += fac
        end
      end
    end
  end

# don't do this now because the reverse sweep needs hk_all to have the
# same values as during the same time during the forward calculation
#=
  # hk_all now contains the p_i in each direction
  # Now compute h_k tilde from p_i
  for i=1:mesh.numEl

    jac_i = ro_sview(mesh.jac, :, i)
    vol = zero(T)
    for j=1:mesh.numNodesPerElement
      vol += sbp.w[j]/jac_i[j]
    end

    h_k = vol^(1/mesh.dim)
    pk_prod = one(T)
    for d=1:mesh.dim
      pk_prod *= hk_all[d, i]
    end

    for d=1:mesh.dim
      hk_all[d, i] = h_k*pk_prod/(hk_all[d, i]*(sbp.degree + 1))
    end
  end
=#

  # reverse sweep
  @simd for i=1:mesh.numEl

    jac_i = ro_sview(mesh.jac, :, i)
    vol = zero(T)
    @simd for j=1:mesh.numNodesPerElement
      vol += sbp.w[j]/jac_i[j]
    end

    h_k = vol^(1/mesh.dim)
    pk_prod = one(T)
    @simd for d=1:mesh.dim
      pk_prod *= hk_all[d, i]
    end
    pk_prod_d = pk_prod^(1/mesh.dim)

    pk_prod_d_bar = zero(T)
    h_k_bar = zero(T)
    @simd for d=1:mesh.dim
      #hk_all[d, i] = h_k*pk_prod_d/(hk_all[d, i]*(sbp.degree + 1))
      tmp = hk_all[d, i]*(sbp.degree + 1)
      h_k_bar += pk_prod_d*hk_all_bar[d, i]/tmp
      pk_prod_d_bar += h_k*hk_all_bar[d, i]/tmp
      hk_all_bar[d, i] = -h_k*pk_prod_d*hk_all_bar[d, i]/(hk_all[d, i]*tmp)
    end

    pk_prod_bar = (pk_prod_d_bar/mesh.dim)*(pk_prod)^((1/mesh.dim) - 1)
    @simd for d=mesh.dim:-1:1
      pk_prod /= hk_all[d, i]  # reset primal value
      hk_all_bar[d, i] += pk_prod*pk_prod_bar
      pk_prod_bar = hk_all[d, i]*pk_prod_bar
    end

    jac_bar_i = sview(mesh.jac_bar, :, i)
    vol_bar = (h_k_bar/mesh.dim)*vol^((1/mesh.dim) - 1)
    @simd for j=1:mesh.numNodesPerElement
      jac_bar_i[j] += -sbp.w[j]*vol_bar/(jac_i[j]*jac_i[j])
    end

  end  # end i

  @simd for i=1:mesh.numInterfaces
    iface_i = mesh.interfaces[i]
    nrm_i = ro_sview(mesh.nrm_face, :, :, i)
    nrm_bar_i = sview(mesh.nrm_face_bar, :, :, i)

    @simd for j=1:mesh.numNodesPerFace
      @simd for k=1:mesh.dim
        fac = absvalue2_deriv(nrm_i[k, j])*mesh.sbpface.wface[j]
        nrm_bar_i[k, j] += fac*hk_all_bar[k, iface_i.elementL]
        nrm_bar_i[k, j] += fac*hk_all_bar[k, iface_i.elementR]
      end
    end
  end

  @simd for i=1:mesh.numBoundaryFaces
    bndry_i = mesh.bndryfaces[i]
    nrm_i = ro_sview(mesh.nrm_bndry, :, :, i)
    nrm_bar_i = sview(mesh.nrm_bndry_bar, :, :, i)

    @simd for j=1:mesh.numNodesPerFace
      @simd for k=1:mesh.dim
        fac = mesh.sbpface.wface[j]*absvalue2_deriv(nrm_i[k, j])
        nrm_bar_i[k, j] += fac*hk_all_bar[k, bndry_i.element]
      end
    end
  end

  @simd for peer=1:mesh.npeers
    @simd for i=1:length(mesh.shared_interfaces[peer])
      iface_i = mesh.shared_interfaces[peer]
      nrm_i = ro_sview(mesh.nrm_sharedface[peer], :, :, i)
      nrm_bar_i = sview(mesh.nrm_sharedface_bar[peer], :, :, i)

      @simd for j=1:mesh.numNodesPerFace
        @simd for k=1:mesh.dim
          fac = mesh.sbpface.wface[j]*absvalue2_deriv(nrm_i[k, j])
          nrm_bar_i[k, j] += fac*hk_all_bar[k, iface_i.elementL]
        end
      end
    end
  end

  # finish the forward computation, so hk_all has the right values in it on
  # exit.
  # hk_all now contains the p_i in each direction
  # Now compute h_k tilde from p_i
  @simd for i=1:mesh.numEl

    jac_i = ro_sview(mesh.jac, :, i)
    vol = zero(T)
    @simd for j=1:mesh.numNodesPerElement
      vol += sbp.w[j]/jac_i[j]
    end

    h_k = vol^(1/mesh.dim)
    pk_prod = one(T)
    @simd for d=1:mesh.dim
      pk_prod *= hk_all[d, i]
    end
    pk_prod = pk_prod^(1/mesh.dim)

    @simd for d=1:mesh.dim
      hk_all[d, i] = h_k*pk_prod/(hk_all[d, i]*(sbp.degree + 1))
    end
  end

  return nothing
end

#------------------------------------------------------------------------------
# SensorHApprox

# this shock sensor is not a function of q, so the _diff and _revq methods
# are not required

# even though this shared the name with all the other getShockSensor methods,
# it signature and outputs are different.  This is sort-of ok because this
# is a special shock sensor and should never be used as the main shock sensor.
function getShockSensor_revm(params::ParamType{Tdim}, sbp::AbstractOperator,
                        sensor::ShockSensorHApprox{Tsol, Tres},
                        elnum::Integer, jac::AbstractVector{Tmsh},
                        jac_bar::AbstractVector, epsL_bar::Number
                       ) where {Tsol, Tres, Tmsh, Tdim}

  numNodesPerElement = length(jac)

  shared_idx = elnum - first(sensor.shared_els) + 1  # index of shared element
  if (elnum in sensor.local_els) ||
     (elnum in sensor.shared_els && sensor.shared_isShocked[shared_idx])

    # compute element size h for both elements
    h = zero(Tmsh)
    @simd for i=1:numNodesPerElement
      h += sbp.w[i]/jac[i]
    end

    # compute estimates viscosity
    lambda_max = sensor.lambda_max
    fac = sensor.alpha*lambda_max/sbp.degree
    #eps_L = fac*(h^(1/Tdim))

    #----------------------------------------
    # reverse sweep

    h_bar = fac*epsL_bar*(h^(1/Tdim - 1))/Tdim

    @simd for i=1:numNodesPerElement
      jac_bar[i] += -sbp.w[i]*h_bar/(jac[i]*jac[i])
    end
  end

end

#------------------------------------------------------------------------------
# ShockSensorOddBO

function getShockSensor_diff(params::ParamType{Tdim}, sbp::AbstractOperator,
                      sensor::ShockSensorOddBO,
                      q_el::AbstractMatrix{Tsol}, elnum::Integer,
                      coords::AbstractMatrix, dxidx::Abstract3DArray,
                      jac::AbstractVector{Tmsh},
                      Se_jac::Abstract4DArray{Tres},
                      ee_jac::Abstract4DArray{Tres}) where {Tsol, Tmsh, Tres, Tdim}

  if elnum % 2 == 1
    val = getShockSensor_diff(params, sbp, sensor.sensor, q_el, elnum, coords,
                              dxidx, jac, Se_jac, ee_jac)
  else
    fill!(Se_jac, 0)
    fill!(ee_jac, 0)
    val = true
  end

  return val
end

function getShockSensor_revq(params::ParamType{Tdim}, sbp::AbstractOperator,
                        sensor::ShockSensorOddBO{Tsol, Tres},
                        q::AbstractMatrix, q_bar::AbstractMatrix,
                        elnum::Integer,
                        coords::AbstractMatrix,
                        dxidx::Abstract3DArray, jac::AbstractVector{Tmsh},
                        ee_mat::AbstractMatrix, ee_mat_bar::AbstractMatrix
                        ) where {Tsol, Tres, Tmsh, Tdim}

  if elnum % 2 == 1
    getShockSensor_revq(params, sbp, sensor.sensor, q, q_bar, elnum, coords, dxidx,
                        jac, ee_mat, ee_mat_bar)
  end

  return nothing
end

function getShockSensor_revm(params::ParamType{Tdim}, sbp::AbstractOperator,
                        sensor::ShockSensorOddBO{Tsol, Tres},
                        q::AbstractMatrix,
                        elnum::Integer,
                        coords::AbstractMatrix, coords_bar::AbstractMatrix,
                        dxidx::Abstract3DArray, dxidx_bar::Abstract3DArray,
                        jac::AbstractVector{Tmsh}, jac_bar::AbstractVector,
                        ee_bar::AbstractMatrix
                       ) where {Tsol, Tres, Tmsh, Tdim}

  if elnum % 2 == 1
    getShockSensor_revm(params, sbp, sensor.sensor, q, elnum, coords, coords_bar,
                        dxidx, dxidx_bar, jac, jac_bar, ee_bar)
  end

  return nothing
end


text:
abstract type Day5 <: AocDay end
day5 = Day5

function readinput(day::Type{Day5})
    input = readlines("input/day5.txt")

end

function getrowbybp(bp)
    rowinfo = bp[1:7]
    range_start = 0
    range_end = 127
    mid = convert(Int, round(length(range_start:range_end) / 2))
    for r in rowinfo
        if r == 'F'
            range_end -= mid
        end
        if r == 'B'
            range_start += mid
        end
        mid = convert(Int, round(length(range_start:range_end) / 2))
    end
    range_start
end

function getcolbybp(bp)
    colinfo = bp[8:10]
    range_start = 0
    range_end = 7
    mid = convert(Int, round(length(range_start:range_end) / 2))
    for c in colinfo
        if c == 'L'
            range_end -= mid
        end
        if c == 'R'
            range_start += mid
        end
        mid = convert(Int, round(length(range_start:range_end) / 2))
    end
    range_start
end

function getseat(bp)
    row = getrowbybp(bp)
    col = getcolbybp(bp)
    seat = row * 8 + col
    Dict(:row => row, :col => col, :seat => seat)
end

function solve(day::Type{Day5})
    input = readinput(day)
    rcs = getseat.(input)
    maximum(map(x->x[:seat], rcs))
end

function solve_part2(day::Type{day5})
    input = readinput(day)
    rcs = getseat.(input)
    plane_seats = sort([r[:seat] for r in rcs])
    min_plane_seat, max_plane_seat = minimum(plane_seats), maximum(plane_seats)
    full_range = collect(min_plane_seat:max_plane_seat)

    for seat in full_range
        if !(seat in plane_seats)
            return seat
        end
    end
end

ext: jl


text:
using LinearAlgebra
using Plots

function main()

    # vector_ops()
    matrix_ops()

end

function vector_ops()
    x = [1, 2, 3]
    y = [1; 2; 3]
    z = [1 2 3]
    q = [1 2 3; 4 5 6]

    # display(x)
    # display(y)
    # display(z)
    # display(q)

    # norm
    # println(norm(x))
    # println(x / norm(x))
    # println(normalize(x))

    # dot product
    # println(dot(x, y))
    # println(x ⋅ y)

    # cross product
    println(cross(x,y))
    println(x × y)

end

function matrix_ops()

    A = [0 2 4; 2 4 2; 3 3 1]
    B = ones(3,3)
    C = [1 1; 1 1; 1 1]

    # matrix multiplication
    # display(A * B)
    # display(A * C)

    # elementwise multiplication
    # display(A .* B)

    # special operations
    # println(tr(A))
    # println(det(A))
    # display(transpose(A))
    # display(eigvals(A))
    # display(eigvecs(A))
    # display(inv(A))

    # special matrices
    # D = [1 0 0; 0 1 0; 0 0 1]
    # display(D)
    # D = Diagonal(D)

    # linear system Ax = b
    A = [1 2; 1 -1]
    b = [5, -1]

    # A = factorize(A)
    # display(A)
    # display(A.L * A.U)

    x0 = A \ b

    xlin = range(0, 5, length=50)
    f(x) = 5 - 2x
    g(x) = -1 + x

    println(x0)

    plot(xlin, f.(xlin), label="f(x)", xaxis="x", yaxis="y")
    plot!(xlin, g.(xlin), label="g(x)", title="Showing Solved Values")
    scatter!((x0[2], x0[1]), label="solved values", dpi=150)
end

main()

max_stars_repo_name: nunesmelo/djs-office-hours

max_issues_repo_name: nunesmelo/djs-office-hours

max_forks_repo_name: nunesmelo/djs-office-hours


text:
export GphysChannel

abstract type GphysChannel end

function isequal(C::T, D::T) where {T<:GphysChannel}
  q::Bool = true
  F = fieldnames(T)
  for f in F
    if f != :notes
      q = min(q, getfield(C,f) == getfield(D,f))
    end
  end
  return q
end
==(C::T, D::T) where {T<:GphysChannel} = isequal(C,D)

in(s::String, C::GphysChannel) = getfield(C, :id)==s

@doc (@doc namestrip)
namestrip!(C::T) where {T<:GphysChannel} = setfield!(C, :name, namestrip(getfield(C, :name)))


text:
using Distributions

"""
internal type of a random walk (for one cell)
"""
type Random_Walk

    steps_length::Array{Float64,1}
    steps_angles::Array{Float64,1}
    steps_points::Array{Tuple{Float64,Float64},1}

    #save bias and persistance parameters (relevant only for random_walk_gradient model)
    bias::Array{Float64,1}
    persistance::Array{Float64,1}

    """constructor for random walks without bias and persistance arrays"""
    function Random_Walk(steps_length::Array{Float64,1},steps_angles::Array{Float64,1},steps_points::Array{Tuple{Float64,Float64},1})
        new(steps_length,steps_angles,steps_points,zeros(1),zeros(1))
    end

    function Random_Walk(steps_length::Array{Float64,1},steps_angles::Array{Float64,1},steps_points::Array{Tuple{Float64,Float64},1},bias::Array{Float64,1},persistance::Array{Float64,1})
        new(steps_length,steps_angles,steps_points,bias,persistance)
    end

    function Random_Walk()
        new()
    end

end

abstract ABC_Model

""" internal representation of a model

#### Fields

#### Constructor
** Create a new model type **

"""
type Deterministic_Model <: ABC_Model

    #model informaiton
    name::ASCIIString
    number::Int

    #species and parameters
    num_species::Int
    num_params::Int

    #prior distributions for parameters
    priors::Array{Union{Float64,Distributions.Uniform,Distributions.DiscreteUniform,Distributions.Normal,Distributions.LogNormal},1}

    #integration properties
    integration_mode::ASCIIString
    integration_function::Function
    initial_conditions::Array{Float64,1}

    """this is the model constructor"""
    function Deterministic_Model()

        new()

    end

end

type Stochastic_Model <: ABC_Model

    #model informaiton
    name::ASCIIString
    number::Int

    #species and parameters
    num_species::Int
    num_params::Int

    #prior distributions for parameters
    priors::Array{Union{Float64,Distributions.Uniform,Distributions.DiscreteUniform,Distributions.Normal,Distributions.LogNormal},1}
    #priors::Array{Any,1}

    #integration properties
    integration_mode::ASCIIString
    integration_repeats::Int64
    integration_matrix::Array{Array{Int64,1},1}
    integration_hazards::Array{Function,1}
    initial_conditions::Array{Float64,1}

    """this is the model constructor"""
    function Stochastic_Model()

        new()

    end

end

type Random_Walk_model <: ABC_Model

    #model informaiton
    name::ASCIIString
    number::Int

    #parameters
    num_params::Int

    #prior distributions for parameters
    priors::Array{Union{Float64,Distributions.Uniform,Distributions.DiscreteUniform,Distributions.Normal,Distributions.LogNormal},1}
    #priors::Array{Any,1}

    #integration properties
    integration_mode::ASCIIString

    #walk properties - can be set so that there is only one type of walks, or so that there are different clusters of walks
    n_steps::Union{Int64,Array{Int64,1},Array{Array{Int64,1},1}} #number of steps for each walk - can all be the same, or each walk can have a different length
    n_walks::Int64 #number of total walks (cells)
    n_clusters::Int64 #number of spatio-temporal clusters
    start_point::Union{Array{Tuple{Float64,Float64},1},Array{Array{Tuple{Float64,Float64},1},1}}
    start_time::Union{Array{Int64,1},Array{Array{Int64,1},1}}
    cell_radius::Float64
    bias_angle::Float64

    """this is the model constructor"""
    function Random_Walk_model()

        new()

    end

end

""" internal representation of an Algorithm
"""
type ABC_Algorithm

    #data info
    time_points::FloatRange{Float64}
    data::Union{Array{Array{Float64,1},1},Array{Random_Walk,1},Array{Array{Random_Walk,1},1}}

    #epsilon
    epsilon::Array{Float64,1}

    #number of runs for each step
    particles::Int

    #perturbation kernel
    kernels::Array{Distributions.Uniform,1}

    #distance function
    distance_func::ASCIIString

    """this is the algorithm constructor"""
    function ABC_Algorithm()

        new()

    end

end

"""
internal type of a population from an ABC SMC simulation
"""
type ABC_population

    index::Int64

    #simulations information
    num_steps::Int64
    epsilon::Float64

    #results
    data::Array{Array{Union{Array{Array{Float64,1},1},Array{Random_Walk,1}},1},1}
    particles::Array{Array{Float64,1},1}
    weights::Array{Float64,1}
    distances::Array{Float64,1}

end


text:
mutable struct MaximinModel <: AbstractModel
    settings::Settings
    constraints::Vector{Constraint}
    obj::MaximinObjective
    output::Output
    MaximinModel(settings, constraints, obj, output) =
        new(settings, constraints, obj, output)
    MaximinModel() = new(Settings(), Constraint[], MaximinObjective(), Output())
end


text:
function _make_table_gui()

    #Table of Values Popup
    table_grid=Grid()

    table_list=ListStore(Int32,Int32,Int32,Int32,Bool)
    for i=1:0
        push!(table_list,(i,125,0,0,true))
    end

    table_tv=TreeView(TreeModel(table_list))
    table_rtext1=CellRendererText()
    table_rtext2=CellRendererText()
    setproperty!(table_rtext2, :editable, true)
    table_rtext3=CellRendererText()
    setproperty!(table_rtext3, :editable, true)
    table_rtext4=CellRendererText()
    setproperty!(table_rtext4, :editable, true)
    table_rtog=CellRendererToggle()
    setproperty!(table_rtog, :activatable, true)

    table_c1 = TreeViewColumn("Channel",table_rtext1,Dict([("text",0)]))
    table_c2 = TreeViewColumn("Gain", table_rtext2, Dict([("text",1)]))
    table_c3 = TreeViewColumn("Threshold", table_rtext3, Dict([("text",2)]))
    table_c4 = TreeViewColumn("Reference", table_rtext4, Dict([("text",3)]))
    table_c5 = TreeViewColumn("Enabled",table_rtog,Dict([("active",4)]))

    push!(table_tv,table_c1)
    push!(table_tv,table_c2)
    push!(table_tv,table_c3)
    push!(table_tv,table_c4)
    push!(table_tv,table_c5)

    table_scroll=ScrolledWindow()
    Gtk.GAccessor.min_content_height(table_scroll,500)
    Gtk.GAccessor.min_content_width(table_scroll,500)
    push!(table_scroll,table_tv)

    table_grid[1,1]=table_scroll

    table_win=Window(table_grid)
    setproperty!(table_win, :title, "Parameter List")

    Gtk.showall(table_win)
    visible(table_win,false)

    table_widgets=Table_Widgets(table_win,table_tv,table_list,table_rtext2,table_rtext3,table_rtext4,table_rtog)
end

function add_parameter_callbacks(w,handles,r,fpga)

    signal_connect(w.win, :delete_event) do widget, event
        visible(w.win,false)
        true
    end

    signal_connect(table_col_cb, w.gain_text,"edited",Void,(Ptr{UInt8},Ptr{UInt8}),false,(handles,r,2))
    signal_connect(table_col_cb, w.thres_text,"edited",Void,(Ptr{UInt8},Ptr{UInt8}),false,(handles,r,3))
    signal_connect(table_col_cb, w.ref_text,"edited",Void,(Ptr{UInt8},Ptr{UInt8}),false,(handles,r,4))
    signal_connect(table_en_cb, w.enabled_toggle, "toggled",Void,(Ptr{UInt8},),false,(handles,r))

    nothing
end

function table_col_cb(widget::Ptr, path,new_text,user_data::Tuple{Gui_Handles,RHD2000,Int64})

    han,rhd,col = user_data

    selmodel = Gtk.GAccessor.selection(han.table_widgets.tv)

    iter=Gtk.selected(selmodel)

    num=parse(Int64,unsafe_string(new_text))

    setindex!(han.table_widgets.list,num,iter,col)

    nothing
end

function table_en_cb(widget::Ptr,path,user_data::Tuple{Gui_Handles,RHD2000})

    han,rhd = user_data

    num=parse(Int64,unsafe_string(path))+1

    val=getindex(han.table_widgets.list,num,5)

    setindex!(han.table_widgets.list,!val,num,5)

    nothing
end

function table_cb(widget::Ptr, user_data::Tuple{Gui_Handles,RHD2000})

    han, rhd = user_data

    for i=1:size(rhd.v,2)

        setindex!(han.table_widgets.list,round(Int64,han.scale[i,1]*-1000),i,2)
        setindex!(han.table_widgets.list,han.sc.thres,i,3)
        setindex!(han.table_widgets.list,rhd.refs[i],i,4)
        setindex!(han.table_widgets.list,han.enabled[i],i,5)

    end

    visible(han.table_widgets.win,true)
    nothing
end


text:
using RandomSeparations
using Test

@testset "RandomSeparations.jl" begin

    X_pos = ones(Float32, 10, 3)     # known positive samples
    X_neg = rand(Float32, 1000, 3)   # negative samples
    X_xxx = [0.99f0 0.99f0 0.99f0]   # unknown positive sample
    X_unl = vcat(X_neg, X_xxx)       # unlabeled samples

    X = vcat(X_unl, X_pos)
    Y = vcat(falses(size(X_unl, 1)), trues(size(X_pos, 1)))

    scores = RandomSeparations.run(X, Y)

    @test findmax(scores)[2] == 1001 # index of the unknown positive sample

end


text:
#const Evt{T} = NamedTuple{(:pos, :val),Tuple{Int,T}}
#Evt(pos::Int, val::T) where T = (pos = pos, val = val)

"""
`flt = MaxMinFilter{T}(w)`

Moving maximum-minimum filter of window `w` that maintains state between calls.
"""
mutable struct MaxMinFilter{T}
    window::Int
    x0::T
    need_restart::Bool
    U_buf::RingBuffer{Evt{T}}
    L_buf::RingBuffer{Evt{T}}
    MaxMinFilter{T}(window::Int) where T =
        new{T}(window, 0, true, RingBuffer{Evt{T}}(window), RingBuffer{Evt{T}}(window))
end

function reset!(state::MaxMinFilter)
    reset!(state.L_buf, state.window)
    reset!(state.U_buf, state.window)
    state.need_restart = true
    state
end

function _restart!(state::MaxMinFilter, x0)
    state.x0 = x0
    state.need_restart = false
    state
end

function run(state::MaxMinFilter{T}, x::AbstractArray{T}) where T
    xmax = similar(x)
    xmin = similar(x)
    run!(xmax, xmin, state, x)
    xmax, xmin
end
function run!(xmax::AbstractArray{T}, xmin::AbstractArray{T}, # output
              state::MaxMinFilter{T},
              x::AbstractArray{T}) where T # input

    if state.need_restart && !isempty(x)
        _restart!(state, first(x))
    end
    w = state.window
    x0 = state.x0
    U_buf = state.U_buf
    L_buf = state.L_buf

    @inbounds for i = eachindex(x)
        xi = x[i]
        if xi > x0
            push!(L_buf, Evt(i-1, x0))
            if i == w + first(L_buf).pos
                popfirst!(L_buf)
            end
            while length(U_buf) > 0
                if xi <= last(U_buf).val
                    if i == w + first(U_buf).pos
                        popfirst!(U_buf)
                    end
                    break
                end
                pop!(U_buf)
            end
        else
            push!(U_buf, Evt(i-1, x0))
            if i == w + first(U_buf).pos
                popfirst!(U_buf)
            end
            while length(L_buf) > 0
                if xi >= last(L_buf).val
                    if i == w + first(L_buf).pos
                        popfirst!(L_buf)
                    end
                    break
                end
                pop!(L_buf)
            end
        end
        xmax[i] = ifelse(length(U_buf) > 0, first(U_buf).val, xi) ### x[get_front[U]]
        xmin[i] = ifelse(length(L_buf) > 0, first(L_buf).val, xi) ### x[get_front[L]]
        x0 = xi
    end

    # iterate positions
    len = length(x)
    for i in eachindex(L_buf)
        L_buf[i] = Evt(L_buf[i].pos - len, L_buf[i].val)
    end
    for i in eachindex(U_buf)
        U_buf[i] = Evt(U_buf[i].pos - len, U_buf[i].val)
    end

    state.x0 = x0

    xmax, xmin
end
filt(state::MaxMinFilter, x) = run(state, x)
filt!(xmax, xmin, state::MaxMinFilter, x) = run!(xmax, xmin, state, x)

delay(state::MaxMinFilter) = state.window ÷ 2

"""
`flt = RangeFilter{T}(w)`

Moving max-min range filter of window `w` that maintains state between calls.
"""
mutable struct RangeFilter{T}
    window::Int
    x0::T
    need_restart::Bool
    U_buf::RingBuffer{Evt{T}}
    L_buf::RingBuffer{Evt{T}}
    RangeFilter{T}(window::Int) where T =
        new{T}(window, 0, true, RingBuffer{Evt{T}}(window), RingBuffer{Evt{T}}(window))
end

function reset!(state::RangeFilter)
    reset!(state.L_buf, state.window)
    reset!(state.U_buf, state.window)
    state.need_restart = true
    state
end

function _restart!(state::RangeFilter, x0)
    state.x0 = x0
    state.need_restart = false
    state
end

function run(state::RangeFilter{T}, x::AbstractArray{T}) where T
    xrange = similar(x)
    run!(xrange, state, x)
    xrange
end
function run!(xrange::AbstractArray{T}, # output
              state::RangeFilter{T},
              x::AbstractArray{T}) where T # input

    if state.need_restart && !isempty(x)
        _restart!(state, first(x))
    end
    w = state.window
    x0 = state.x0
    U_buf = state.U_buf
    L_buf = state.L_buf

    @inbounds for i = eachindex(x)
        xi = x[i]
        if xi > x0
            push!(L_buf, Evt(i-1, x0))
            if i == w + first(L_buf).pos
                popfirst!(L_buf)
            end
            while length(U_buf) > 0
                if xi <= last(U_buf).val
                    if i == w + first(U_buf).pos
                        popfirst!(U_buf)
                    end
                    break
                end
                pop!(U_buf)
            end
        else
            push!(U_buf, Evt(i-1, x0))
            if i == w + first(U_buf).pos
                popfirst!(U_buf)
            end
            while length(L_buf) > 0
                if xi >= last(L_buf).val
                    if i == w + first(L_buf).pos
                        popfirst!(L_buf)
                    end
                    break
                end
                pop!(L_buf)
            end
        end
        xmax = ifelse(length(U_buf) > 0, first(U_buf).val, xi) ### x[get_front[U]]
        xmin = ifelse(length(L_buf) > 0, first(L_buf).val, xi) ### x[get_front[L]]
        xrange[i] = xmax - xmin
        x0 = xi
    end

    # iterate positions
    len = length(x)
    for i in eachindex(L_buf)
        L_buf[i] = Evt(L_buf[i].pos - len, L_buf[i].val)
    end
    for i in eachindex(U_buf)
        U_buf[i] = Evt(U_buf[i].pos - len, U_buf[i].val)
    end
    state.x0 = x0

    xrange
end
filt(state::RangeFilter, x) = run(state, x)
filt!(xrange, state::RangeFilter, x) = run!(xrange, state, x)

delay(state::RangeFilter) = state.window ÷ 2

"""
`flt = MaxFilter{T}(w)`

Moving maximum filter of window `w` that maintains state between calls.
"""
mutable struct MaxFilter{T}
    window::Int
    x0::T
    need_restart::Bool
    U_buf::RingBuffer{Evt{T}}
    MaxFilter{T}(window::Int) where T =
        new{T}(window, 0, true, RingBuffer{Evt{T}}(window))
end

function reset!(state::MaxFilter)
    reset!(state.U_buf, state.window)
    state.need_restart = true
    state
end

function _restart!(state::MaxFilter, x0)
    state.x0 = x0
    state.need_restart = false
    state
end

function run(state::MaxFilter{T}, x::AbstractArray{T}) where T
    xmax = similar(x)
    run!(xmax, state, x)
    xmax
end
function run!(xmax::AbstractArray{T}, # output
              state::MaxFilter{T},
              x::AbstractArray{T}) where T # input

    if state.need_restart && !isempty(x)
        _restart!(state, first(x))
    end
    w = state.window
    x0 = state.x0
    U_buf = state.U_buf

    @inbounds for i = eachindex(x)
        xi = x[i]
        if xi > x0
            ### ...
            while length(U_buf) > 0
                if xi <= last(U_buf).val
                    if i == w + first(U_buf).pos
                        popfirst!(U_buf)
                    end
                    break
                end
                pop!(U_buf)
            end
        else
            push!(U_buf, Evt(i-1, x0))
            if i == w + first(U_buf).pos
                popfirst!(U_buf)
            end
            ### ...
        end
        xmax[i] = ifelse(length(U_buf) > 0, first(U_buf).val, xi) ### x[get_front[U]]
        x0 = xi
    end

    # iterate positions
    len = length(x)
    for i in eachindex(U_buf)
        U_buf[i] = Evt(U_buf[i].pos - len, U_buf[i].val)
    end
    state.x0 = x0

    xmax
end
filt(state::MaxFilter, x) = run(state, x)
filt!(xmax, state::MaxFilter, x) = run!(xmax, state, x)

delay(state::MaxFilter) = state.window ÷ 2

"""
`flt = MinFilter{T}(w)`

Moving minimum filter of window `w` that maintains state between calls.
"""
mutable struct MinFilter{T}
    window::Int
    x0::T
    need_restart::Bool
    L_buf::RingBuffer{Evt{T}}
    MinFilter{T}(window::Int) where T =
        new{T}(window, 0, true, RingBuffer{Evt{T}}(window))
end

function reset!(state::MinFilter)
    reset!(state.L_buf, state.window)
    state.need_restart = true
    state
end

function _restart!(state::MinFilter, x0)
    state.x0 = x0
    state.need_restart = false
    state
end

function run(state::MinFilter{T}, x::AbstractArray{T}) where T
    xmin = similar(x)
    run!(xmin, state, x)
    xmin
end
function run!(xmin::AbstractArray{T}, # output
              state::MinFilter{T},
              x::AbstractArray{T}) where T # input

    if state.need_restart && !isempty(x)
        _restart!(state, first(x))
    end
    w = state.window
    x0 = state.x0
    L_buf = state.L_buf

    @inbounds for i = eachindex(x)
        xi = x[i]
        if xi > x0
            push!(L_buf, Evt(i-1, x0))
            if i == w + first(L_buf).pos
                popfirst!(L_buf)
            end
            # ...
        else
            # ...
            while length(L_buf) > 0
                if xi >= last(L_buf).val
                    if i == w + first(L_buf).pos
                        popfirst!(L_buf)
                    end
                    break
                end
                pop!(L_buf)
            end
        end
        xmin[i] = ifelse(length(L_buf) > 0, first(L_buf).val, xi) ### x[get_front[U]]
        x0 = xi
    end

    # iterate positions
    len = length(x)
    for i in eachindex(L_buf)
        L_buf[i] = Evt(L_buf[i].pos - len, L_buf[i].val)
    end
    state.x0 = x0

    xmin
end
filt(state::MinFilter, x) = run(state, x)
filt!(xmin, state::MinFilter, x) = run!(xmin, state, x)

delay(state::MinFilter) = state.window ÷ 2

"""
`flt = EnvelopeFilter{T}(w1, w2)`

Moving max-min envelope filter that maintains state between calls.
Envelope is defined as consecutive min-max RangeFilter of window `w1` (1-d dilation)
and MinFilter of window `w2` (1-d erosion). By default `w2 = w1`.
Filter delay = `(w1 + w2) ÷ 2`
"""
mutable struct EnvelopeFilter{T}
    rangefilter::RangeFilter{T}
    minfilter::MinFilter{T}
    EnvelopeFilter{T}(w1::Int, w2::Int = w1-1) where T =
        new{T}(RangeFilter{T}(w1), MinFilter{T}(w2))
end

function reset!(state::EnvelopeFilter)
    reset!(state.rangefilter)
    reset!(state.minfilter)
    state
end

function run(state::EnvelopeFilter{T}, x::AbstractArray{T}) where T
    xenvelope = similar(x)
    run!(xenvelope, state, x)
    xenvelope
end
function run!(xenvelope::AbstractArray{T},
              state::EnvelopeFilter{T},
              x::AbstractArray{T}) where T
    run!(xenvelope, state.rangefilter, x)
    if (state.minfilter.window > 0)
        run!(xenvelope, state.minfilter, xenvelope)
    end
end
filt(state::EnvelopeFilter, x) = run(state, x)
filt!(xenvelope, state::EnvelopeFilter, x) = run!(xenvelope, state, x)

delay(state::EnvelopeFilter) = (state.rangefilter.window + state.minfilter.window) ÷ 2


text:
ismacro(f::Function) = startswith(string(methods(f).mt.name), "@")

wstype(x) = ""
wstype(::Module) = "module"
wstype(f::Function) = "function"
wstype(::Type) = "type"
wstype(::Expr) = "mixin"
wstype(::Symbol) = "tag"
wstype(::AbstractString) = "property"
wstype(::Number) = "constant"
wstype(::Exception) = "tag"

wsicon(x) = ""
wsicon(f::Function) = ismacro(f) ? "icon-mention" : ""
wsicon(::AbstractArray) = "icon-file-binary"
wsicon(::AbstractVector) = "icon-list-ordered"
wsicon(::AbstractString) = "icon-quote"
wsicon(::Expr) = "icon-code"
wsicon(::Symbol) = "icon-code"
wsicon(::Exception) = "icon-bug"
wsicon(::Number) = "n"

wsnamed(name, val) = false
wsnamed(name, f::Function) = name == methods(f).mt.name
wsnamed(name, m::Module) = name == module_name(m)
wsnamed(name, T::DataType) = name == Symbol(T.name)

function wsitem(name, val)
  d(:name  => name,
    :value => render′(Inline(), val),
    :type  => wstype(val),
    :icon  => wsicon(val))
end

wsitem(mod::Module, name::Symbol) = wsitem(name, getfield(mod, name))

handle("workspace") do mod
  mod = getmodule′(mod)
  ns = filter!(x->!Base.isdeprecated(mod, x), Symbol.(CodeTools.filtervalid(names(mod, all=true))))
  filter!(n -> isdefined(mod, n), ns)
  # TODO: only filter out imported modules
  filter!(n -> !isa(getfield(mod, n), Module), ns)
  contexts = [d(:context => string(mod), :items => map(n -> wsitem(mod, n), ns))]
  if isdebugging()
    prepend!(contexts, Debugger.contexts())
  end
  return contexts
end


text:
# Autogenerated wrapper script for COINBLAS_jll for powerpc64le-linux-gnu-libgfortran3
export libcoinblas

## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"

# Relative path to `libcoinblas`
const libcoinblas_splitpath = ["lib", "libcoinblas.so"]

# This will be filled out by __init__() for all products, as it must be done at runtime
libcoinblas_path = ""

# libcoinblas-specific global declaration
# This will be filled out by __init__()
libcoinblas_handle = C_NULL

# This must be `const` so that we can use it with `ccall()`
const libcoinblas = "libcoinblas.so.0"

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"COINBLAS")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # We first need to add to LIBPATH_list the libraries provided by Julia
    append!(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)])
    global libcoinblas_path = normpath(joinpath(artifact_dir, libcoinblas_splitpath...))

    # Manually `dlopen()` this right now so that future invocations
    # of `ccall` with its `SONAME` will find this path immediately.
    global libcoinblas_handle = dlopen(libcoinblas_path)
    push!(LIBPATH_list, dirname(libcoinblas_path))

    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(LIBPATH_list, ':')

    # Add each element of LIBPATH to our DL_LOAD_PATH (necessary on platforms
    # that don't honor our "already opened" trick)
    #for lp in LIBPATH_list
    #    push!(DL_LOAD_PATH, lp)
    #end
end  # __init__()


text:
using Query
using DataFrames

source = Dict("John"=>34., "Sally"=>56.)

result = @from i in source begin
         @where i.second>36.
         @select {Name=lowercase(i.first)}
         @collect DataFrame
end

println(result)

result = @from i in source begin
         @where i.second>36.
         @select {Name=lowercase(i.first)}
         @collect
end

println(result)