text:
# differentiated version of functions in flux.jl

#------------------------------------------------------------------------------
# calcFaceIntegral

function calcFaceIntegral_nopre_diff(
                                mesh::AbstractDGMesh{Tmsh},
                                sbp::AbstractOperator,
                                eqn::EulerData{Tsol, Tres, Tdim, :conservative},
                                opts,
                                functor::FluxType_diff,
                                interfaces::AbstractArray{Interface, 1},
                                assembler::AssembleElementData) where {Tmsh, Tsol, Tres, Tdim}

  nfaces = length(interfaces)
  params = eqn.params

  data = params.calc_face_integrals_data
  @unpack data q_faceL q_faceR flux_dotL flux_dotR res_jacLL res_jacLR res_jacRL res_jacRR

  for i=1:nfaces
    iface_i = interfaces[i]
    fill!(flux_dotL, 0)
    fill!(flux_dotR, 0)

    qL = ro_sview(eqn.q, :, :, iface_i.elementL)
    qR = ro_sview(eqn.q, :, :, iface_i.elementR)
    interiorFaceInterpolate!(mesh.sbpface, iface_i, qL, qR, q_faceL,
                             q_faceR)

    # compute dF/dqface
    for j=1:mesh.numNodesPerFace
      qL_j = ro_sview(q_faceL, :, j)
      qR_j = ro_sview(q_faceR, :, j)

      eqn.aux_vars_face[1, j, i] = calcPressure(params, qL_j)
      aux_vars = ro_sview(eqn.aux_vars_face, :, j, i)

      nrm_xy = ro_sview(mesh.nrm_face, :, j, i)
#      flux_j = sview(flux_face, :, j)
      flux_dotL_j = sview(flux_dotL, :, :, j)
      flux_dotR_j = sview(flux_dotR, :, :, j)

      functor(params, qL_j, qR_j, aux_vars, nrm_xy, flux_dotL_j, flux_dotR_j)
    end  # end loop j

    # compute dR/dq
    #TODO: for sparse faces, only zero out the needed entries, to avoid loading
    #      the entire array into cache
    # Tests show this doesn't make a difference
    fill!(res_jacLL, 0.0)
    fill!(res_jacLR, 0.0)
    fill!(res_jacRL, 0.0)
    fill!(res_jacRR, 0.0)

    interiorFaceCombined_jac!(mesh.sbpface, iface_i, flux_dotL, flux_dotR,
                             res_jacLL, res_jacLR, res_jacRL, res_jacRR,
                             SummationByParts.Subtract())

    # multiply by Minv if needed
    if params.use_Minv == 1
      for q=1:mesh.numNodesPerElement
        for p=1:mesh.numNodesPerElement
          valL = mesh.jac[p, iface_i.elementL]/sbp.w[p]  # entry in Minv
          valR = mesh.jac[p, iface_i.elementR]/sbp.w[p]
          @simd for m=1:mesh.numDofPerNode
            @simd for n=1:mesh.numDofPerNode
              res_jacLL[n, m, p, q] *= valL
              res_jacLR[n, m, p, q] *= valL
              res_jacRL[n, m, p, q] *= valR
              res_jacRR[n, m, p, q] *= valR
            end
          end
        end
      end
    end

    # assemble into the Jacobian
    assembleInterface(assembler, mesh.sbpface, mesh, iface_i, res_jacLL,
                      res_jacLR, res_jacRL, res_jacRR)
  end  # end loop i

  return nothing
end

"""
  Reverse mode wrt metrics of [`calcFaceIntegral_nopre`](@ref)

  **Inputs**

   * mesh: bar fields are updated
   * sbp
   * eqn: res_bar fields in input
   * opts
   * functor_revm: reverse mode wrt metrics functor (`FluxType_revm`)
   * interfaces
"""
function calcFaceIntegral_nopre_revm(
        mesh::AbstractDGMesh{Tmsh},
        sbp::AbstractOperator,
        eqn::EulerData{Tsol, Tres, Tdim, :conservative},
        opts,
        functor_revm::FluxType_revm,
        interfaces::AbstractArray{Interface, 1}) where {Tmsh, Tsol, Tres, Tdim}

  nfaces = length(interfaces)
  params = eqn.params

  q_faceL = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceR = zeros(q_faceL)

  flux_face_bar = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)

  # reverse sweep
  for i=1:nfaces
    iface_i = interfaces[i]

    qL = ro_sview(eqn.q, :, :, iface_i.elementL)
    qR = ro_sview(eqn.q, :, :, iface_i.elementR)
    interiorFaceInterpolate!(mesh.sbpface, iface_i, qL, qR, q_faceL,
                             q_faceR)

    fill!(flux_face_bar, 0)
    resL_bar = sview(eqn.res_bar, :, :, iface_i.elementL)
    resR_bar = sview(eqn.res_bar, :, :, iface_i.elementR)
    interiorFaceIntegrate_rev!(mesh.sbpface, iface_i, flux_face_bar, resL_bar,
                               resR_bar, SummationByParts.Subtract())

    for j=1:mesh.numNodesPerFace
      qL_j = ro_sview(q_faceL, :, j)
      qR_j = ro_sview(q_faceR, :, j)

      eqn.aux_vars_face[1, j, i] = calcPressure(params, qL_j)
      aux_vars = ro_sview(eqn.aux_vars_face, :, j, i)

      nrm_xy = sview(mesh.nrm_face, :, j, i)
      nrm_bar = sview(mesh.nrm_face_bar, :, j, i)
      flux_bar_j = sview(flux_face_bar, :, j)

      functor_revm(params, qL_j, qR_j, aux_vars, nrm_xy, nrm_bar, flux_bar_j)
    end  # end loop j

  end  # end loop i

  return nothing
end

"""
  Reverse mode wrt q of [`calcFaceIntegral_nopre`](@ref)

  **Inputs**

   * mesh
   * sbp
   * eqn: `res_bar` field as input, `q_bar` as output
   * opts
   * functor_revq: reverse mode wrt metrics functor (`FluxType_revq`)
   * interfaces
"""
function calcFaceIntegral_nopre_revq(
        mesh::AbstractDGMesh{Tmsh},
        sbp::AbstractOperator,
        eqn::EulerData{Tsol, Tres, Tdim, :conservative},
        opts,
        functor_revq::FluxType_revq,
        interfaces::AbstractArray{Interface, 1}) where {Tmsh, Tsol, Tres, Tdim}

  nfaces = length(interfaces)
  params = eqn.params

  q_faceL = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceR = zeros(q_faceL)
  q_faceL_bar = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceR_bar = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)

  flux_face_bar = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)

  # reverse sweep
  for i=1:nfaces
    iface_i = interfaces[i]

    qL = ro_sview(eqn.q, :, :, iface_i.elementL)
    qR = ro_sview(eqn.q, :, :, iface_i.elementR)
    qL_bar = sview(eqn.q_bar, :, :, iface_i.elementL)
    qR_bar = sview(eqn.q_bar, :, :, iface_i.elementR)
    interiorFaceInterpolate!(mesh.sbpface, iface_i, qL, qR, q_faceL,
                             q_faceR)

    fill!(flux_face_bar, 0)
    resL_bar = sview(eqn.res_bar, :, :, iface_i.elementL)
    resR_bar = sview(eqn.res_bar, :, :, iface_i.elementR)
    interiorFaceIntegrate_rev!(mesh.sbpface, iface_i, flux_face_bar, resL_bar,
                               resR_bar, SummationByParts.Subtract())

    fill!(q_faceL_bar, 0); fill!(q_faceR_bar, 0)
    for j=1:mesh.numNodesPerFace
      qL_j = ro_sview(q_faceL, :, j); qL_bar_j = sview(q_faceL_bar, :, j)
      qR_j = ro_sview(q_faceR, :, j); qR_bar_j = sview(q_faceR_bar, :, j)

      eqn.aux_vars_face[1, j, i] = calcPressure(params, qL_j)
      aux_vars = ro_sview(eqn.aux_vars_face, :, j, i)

      nrm_xy = sview(mesh.nrm_face, :, j, i)
      flux_bar_j = sview(flux_face_bar, :, j)

      functor_revq(params, qL_j, qL_bar_j, qR_j, qR_bar_j, aux_vars, nrm_xy, flux_bar_j)
    end  # end loop j

    interiorFaceInterpolate_rev!(mesh.sbpface, iface_i, qL_bar, qR_bar,
                                 q_faceL_bar, q_faceR_bar)

  end  # end loop i

  return nothing
end

#------------------------------------------------------------------------------
# getFaceElemenIntegral

function getFaceElementIntegral_diff(
                           mesh::AbstractDGMesh{Tmsh},
                           sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
                           face_integral_functor::FaceElementIntegralType,
                           flux_functor::FluxType_diff,
                           sbpface::AbstractFace,
                           interfaces::AbstractArray{Interface, 1},
                           assembler::AssembleElementData) where {Tmsh, Tsol, Tres, Tdim}

  params = eqn.params
#  sbpface = mesh.sbpface
  nfaces = length(interfaces)
#  resL2 = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerElement)
#  resR2 = zeros(resL2)

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

  for i=1:nfaces
    iface = interfaces[i]
    elL = iface.elementL
    elR = iface.elementR
    qL = ro_sview(eqn.q, :, :, elL)
    qR = ro_sview(eqn.q, :, :, elR)
    aux_vars = ro_sview(eqn.aux_vars, :, :, elL)
    nrm_face = ro_sview(mesh.nrm_face, :, :, i)
    resL = sview(eqn.res, :, :, elL)
    resR = sview(eqn.res, :, :, elR)

    fill!(res_jacLL, 0); fill!(res_jacLR, 0)
    fill!(res_jacRL, 0); fill!(res_jacRR, 0)
    calcFaceElementIntegral_diff(face_integral_functor, params, sbpface, iface,
                       qL, qR, aux_vars,
                       nrm_face, flux_functor, res_jacLL, res_jacLR, res_jacRL, res_jacRR)

    assembleInterface(assembler, sbpface, mesh, iface, res_jacLL,
                      res_jacLR, res_jacRL, res_jacRR)

    # multiply by Minv if needed
    if params.use_Minv == 1
      for q=1:mesh.numNodesPerElement
        for p=1:mesh.numNodesPerElement
          valL = mesh.jac[p, iface.elementL]/sbp.w[p]  # entry in Minv
          valR = mesh.jac[p, iface.elementR]/sbp.w[p]
          @simd for m=1:mesh.numDofPerNode
            @simd for n=1:mesh.numDofPerNode
              res_jacLL[n, m, p, q] *= valL
              res_jacLR[n, m, p, q] *= valL
              res_jacRL[n, m, p, q] *= valR
              res_jacRR[n, m, p, q] *= valR
            end
          end
        end
      end
    end

  end

  return nothing
end

"""
  Reverse mode wrt metrics of [`getFaceElementIntegral`](@ref)
"""
function getFaceElementIntegral_revm(
                           mesh::AbstractDGMesh{Tmsh},
                           sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
                           face_integral_functor::FaceElementIntegralType,
                           flux_functor::FluxType,
                           sbpface::AbstractFace,
                           interfaces::AbstractArray{Interface, 1}) where {Tmsh, Tsol, Tres, Tdim}

  params = eqn.params
  nfaces = length(interfaces)
  for i=1:nfaces
    iface = interfaces[i]
    elL = iface.elementL
    elR = iface.elementR
    qL = ro_sview(eqn.q, :, :, elL)
    qR = ro_sview(eqn.q, :, :, elR)
    aux_vars = ro_sview(eqn.aux_vars, :, :, elL)
    nrm_face = ro_sview(mesh.nrm_face, :, :, i)
    nrm_face_bar = sview(mesh.nrm_face_bar, :, :, i)
    resL_bar = ro_sview(eqn.res_bar, :, :, elL)
    resR_bar = ro_sview(eqn.res_bar, :, :, elR)

    calcFaceElementIntegral_revm(face_integral_functor, params, sbpface, iface,
                        qL, qR, aux_vars, nrm_face, nrm_face_bar, flux_functor,
                        resL_bar, resR_bar)

  end

  return nothing
end

"""
  Reverse mode wrt q of [`getFaceElementIntegral`](@ref)
"""
function getFaceElementIntegral_revq(
                           mesh::AbstractDGMesh{Tmsh},
                           sbp::AbstractOperator, eqn::EulerData{Tsol, Tres, Tdim},
                           face_integral_functor::FaceElementIntegralType,
                           flux_functor_revq::FluxType_revq,
                           sbpface::AbstractFace,
                           interfaces::AbstractArray{Interface, 1}) where {Tmsh, Tsol, Tres, Tdim}

  params = eqn.params
  nfaces = length(interfaces)
  for i=1:nfaces
    iface = interfaces[i]
    elL = iface.elementL
    elR = iface.elementR
    qL = ro_sview(eqn.q, :, :, elL); qL_bar = sview(eqn.q_bar, :, :, elL)
    qR = ro_sview(eqn.q, :, :, elR); qR_bar = sview(eqn.q_bar, :, :, elR)
    aux_vars = ro_sview(eqn.aux_vars, :, :, elL)
    nrm_face = ro_sview(mesh.nrm_face, :, :, i)
    resL_bar = ro_sview(eqn.res_bar, :, :, elL)
    resR_bar = ro_sview(eqn.res_bar, :, :, elR)

    calcFaceElementIntegral_revq(face_integral_functor, params, sbpface, iface,
                        qL, qL_bar, qR, qR_bar, aux_vars, nrm_face,
                        flux_functor_revq, resL_bar, resR_bar)

  end

  return nothing
end

#------------------------------------------------------------------------------
# calcSharedFaceIntegrals_element (and _inner)

"""
  Differentiated version of
  [`calcSharedFaceIntegrals_element_inner`](@ref).
  It presents the interface required by [`finishExchangeData`](@ref)
"""
function calcSharedFaceIntegrals_element_diff(
                            mesh::AbstractDGMesh{Tmsh},
                            sbp::AbstractOperator, eqn::EulerData{Tsol},
                            opts, data::SharedFaceData) where {Tmsh, Tsol}

    calcSharedFaceIntegrals_nopre_element_inner_diff(mesh, sbp, eqn, opts, data, eqn.flux_func_diff, eqn.assembler)

  return nothing
end

"""
  Differentiated version of [`calcSharedFaceIntegrals_nopre_element_inner`](@ref),
"""
function calcSharedFaceIntegrals_nopre_element_inner_diff(
                            mesh::AbstractDGMesh{Tmsh},
                            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres},
                            opts, data::SharedFaceData, functor::FluxType_diff,
                            assembler::AssembleElementData) where {Tmsh, Tsol, Tres}

  params = eqn.params
  fdata = params.calc_face_integrals_data
  @unpack fdata q_faceL q_faceR flux_dotL flux_dotR res_jacLL res_jacLR res_jacRL res_jacRR

  # get data
  idx = data.peeridx
  interfaces = data.interfaces
  bndries_local = data.bndries_local
  bndries_remote = data.bndries_remote
#    qL_arr = eqn.q_face_send[i]
  qR_arr = data.q_recv
  nrm_arr = mesh.nrm_sharedface[idx]
  aux_vars_arr = eqn.aux_vars_sharedface[idx]
#  flux_face = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)

  start_elnum = mesh.shared_element_offsets[idx]

  for j=1:length(interfaces)
    iface_j = interfaces[j]
    bndryL_j = bndries_local[j]
    bndryR_j = bndries_remote[j]
    fL = bndryL_j.face

    fill!(flux_dotL, 0); fill!(flux_dotR, 0)

    # interpolate to face
    qL = ro_sview(eqn.q, :, :, iface_j.elementL)
    el_r = iface_j.elementR - start_elnum + 1
    qR = ro_sview(qR_arr, :, :, el_r)
    interiorFaceInterpolate!(mesh.sbpface, iface_j, qL, qR, q_faceL, q_faceR)

    # calculate flux jacobian
    for k=1:mesh.numNodesPerFace
      qL_k = ro_sview(q_faceL, :, k)
      qR_k = ro_sview(q_faceR, :, k)

      aux_vars = ro_sview(aux_vars_arr, :, k, j)
      nrm_xy = ro_sview(nrm_arr, :, k, j)
#      flux_k = sview(flux_face, :, k)
      flux_dotL_k = sview(flux_dotL, :, :, k)
      flux_dotR_k = sview(flux_dotR, :, :, k)

      parent(aux_vars)[1] = calcPressure(params, qL_k)

      functor(params, qL_k, qR_k, aux_vars, nrm_xy, flux_dotL_k, flux_dotR_k)
     end

     # this is excessive because we don't need jacRL, jacRR, but
     # boundaryFaceCombined_jac can't handle jacLR
     fill!(res_jacLL, 0.0)
     fill!(res_jacLR, 0.0)
     interiorFaceCombined_jac!(mesh.sbpface, iface_j, flux_dotL, flux_dotR,
                                res_jacLL, res_jacLR, res_jacRL, res_jacRR,
                                SummationByParts.Subtract())

    # multiply by Minv if needed
    if params.use_Minv == 1
      for q=1:mesh.numNodesPerElement
        for p=1:mesh.numNodesPerElement
          valL = mesh.jac[p, iface_j.elementL]/sbp.w[p]  # entry in Minv
          @simd for m=1:mesh.numDofPerNode
            @simd for n=1:mesh.numDofPerNode
              res_jacLL[n, m, p, q] *= valL
              res_jacLR[n, m, p, q] *= valL
            end
          end
        end
      end
    end

    assembleSharedFace(assembler, mesh.sbpface, mesh, iface_j, res_jacLL, res_jacLR)
   end  # end loop over interfaces

  return nothing
end

"""
  Reverse mode wrt metrics of [`calcSharedFaceIntegrals_element`](@ref)
"""
 function calcSharedFaceIntegrals_element_revm(mesh::AbstractDGMesh{Tmsh},
                            sbp::AbstractOperator, eqn::EulerData{Tsol},
                            opts, data::SharedFaceData) where {Tmsh, Tsol}

  calcSharedFaceIntegrals_nopre_element_inner_revm(mesh, sbp, eqn, opts, data, eqn.flux_func_revm)

  return nothing
end

"""
  Reverse mode of [`calcSharedFaceIntegrals_nopre_element_inner`](@ref) wrt
  the metrics
"""
function calcSharedFaceIntegrals_nopre_element_inner_revm(
                            mesh::AbstractDGMesh{Tmsh},
                            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres},
                            opts, data::SharedFaceData,
                            functor_revm::FluxType_revm
                           ) where {Tmsh, Tsol, Tres}

  if getParallelData(data) != PARALLEL_DATA_ELEMENT
    error("element parallel data required")
  end

  q = eqn.q
  params = eqn.params

  # TODO: make these fields of params
  q_faceL = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceR = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)
  flux_face_bar = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)

  idx = data.peeridx
  interfaces = data.interfaces
  bndries_local = data.bndries_local
  bndries_remote = data.bndries_remote

  # eqn.q_face_send should actually contain el
  # data, but we want to use eqn.q because that's the canonical source
  # qL_arr = eqn.q_face_send[i]
  qR_arr = data.q_recv
#  resR_bar_arr = data_bar.q_recv
  nrm_arr = mesh.nrm_sharedface[idx]
  nrm_arr_bar = mesh.nrm_sharedface_bar[idx]
  aux_vars_arr = eqn.aux_vars_sharedface[idx]

  start_elnum = mesh.shared_element_offsets[idx]

  for j=1:length(interfaces)
    iface_j = interfaces[j]
    bndryL_j = bndries_local[j]
    bndryR_j = bndries_remote[j]
    fL = bndryL_j.face

    # interpolate to face
    qL = ro_sview(q, :, :, iface_j.elementL)
    el_r = iface_j.elementR - start_elnum + 1
    qR = ro_sview(qR_arr, :, :, el_r)
    interiorFaceInterpolate!(mesh.sbpface, iface_j, qL, qR, q_faceL, q_faceR)

    fill!(flux_face_bar, 0)
    resL_bar = sview(eqn.res_bar, :, :, iface_j.elementL)
    #resR_bar = sview(resR_bar_arr, :, :, el_r)
    boundaryFaceIntegrate_rev!(mesh.sbpface, fL, flux_face_bar, resL_bar, SummationByParts.Subtract())

    # calculate flux
    for k=1:mesh.numNodesPerFace
      qL_k = ro_sview(q_faceL, :, k)
      qR_k = ro_sview(q_faceR, :, k)
      aux_vars = ro_sview(aux_vars_arr, :, k, j)
      nrm_xy = ro_sview(nrm_arr, :, k, j)
      nrm_bar = sview(nrm_arr_bar, :, k, j)
      flux_bar_k = sview(flux_face_bar, :, k)

      parent(aux_vars)[1] = calcPressure(params, qL_k)

      functor_revm(params, qL_k, qR_k, aux_vars, nrm_xy, nrm_bar, flux_bar_k)
     end

   end  # end loop over interfaces

  return nothing
end

"""
  Reverse mode wrt q of [`calcSharedFaceIntegrals_element`](@ref)
"""
 function calcSharedFaceIntegrals_element_revq(mesh::AbstractDGMesh{Tmsh},
                            sbp::AbstractOperator, eqn::EulerData{Tsol},
                            opts, data::SharedFaceData,
                            data_bar::SharedFaceData) where {Tmsh, Tsol}

  calcSharedFaceIntegrals_nopre_element_inner_revq(mesh, sbp, eqn, opts, data,
                                                   data_bar, eqn.flux_func_revq)

  return nothing
end

#=
"""
  Reverse mode of [`calcSharedFaceIntegrals_nopre_element_inner`](@ref) wrt q
"""
function calcSharedFaceIntegrals_nopre_element_inner_revq(
                            mesh::AbstractDGMesh{Tmsh},
                            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres},
                            opts, data::SharedFaceData,
                            data_bar::SharedFaceData,
                            functor_revq::FluxType_revq
                           ) where {Tmsh, Tsol, Tres}

  if getParallelData(data) != PARALLEL_DATA_ELEMENT
    error("element parallel data required")
  end

  q = eqn.q
  params = eqn.params

  # TODO: make these fields of params
  q_faceL = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceR = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceL_bar = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceR_bar = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)
  flux_face_bar = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)

  idx = data.peeridx
  interfaces = data.interfaces
  bndries_local = data.bndries_local
  bndries_remote = data.bndries_remote

  # eqn.q_face_send should actually contain el
  # data, but we want to use eqn.q because that's the canonical source
  # qL_arr = eqn.q_face_send[i]
  qR_arr = data.q_recv
  qR_arr_bar = data_bar.q_recv
#  resR_bar_arr = data_bar.q_recv
  nrm_arr = mesh.nrm_sharedface[idx]
  aux_vars_arr = eqn.aux_vars_sharedface[idx]

  start_elnum = mesh.shared_element_offsets[idx]

  for j=1:length(interfaces)
    iface_j = interfaces[j]
    bndryL_j = bndries_local[j]
    bndryR_j = bndries_remote[j]
    elL = iface_j.elementL
    elR = iface_j.elementR - start_elnum + 1
    fL = bndryL_j.face

    # interpolate to face
    qL = ro_sview(q, :, :, elL); qL_bar = sview(eqn.q_bar, :, :, elL)
    qR = ro_sview(qR_arr, :, :, elR); qR_bar = sview(qR_arr_bar, :, :, elR)
    fill!(q_faceL_bar, 0); fill!(q_faceR_bar, 0)
    interiorFaceInterpolate!(mesh.sbpface, iface_j, qL, qR, q_faceL, q_faceR)

    fill!(flux_face_bar, 0)
    resL_bar = sview(eqn.res_bar, :, :, iface_j.elementL)
    #resR_bar = sview(resR_bar_arr, :, :, el_r)
    boundaryFaceIntegrate_rev!(mesh.sbpface, fL, flux_face_bar, resL_bar, SummationByParts.Subtract())

    # calculate flux
    for k=1:mesh.numNodesPerFace
      qL_k = ro_sview(q_faceL, :, k); qL_bar_k = sview(q_faceL_bar, :, k)
      qR_k = ro_sview(q_faceR, :, k); qR_bar_k = sview(q_faceR_bar, :, k)
      aux_vars = ro_sview(aux_vars_arr, :, k, j)
      nrm_xy = ro_sview(nrm_arr, :, k, j)
      flux_bar_k = sview(flux_face_bar, :, k)

      parent(aux_vars)[1] = calcPressure(params, qL_k)

      functor_revq(params, qL_k, qL_bar_k, qR_k, qR_bar_k, aux_vars, nrm_xy,
                   flux_bar_k)
     end

     interiorFaceInterpolate_rev!(mesh.sbpface, iface_j, qL_bar, qR_bar,
                                  q_faceL_bar, q_faceR_bar)

   end  # end loop over interfaces

  return nothing
end
=#

"""
  Reverse mode of [`calcSharedFaceIntegrals_nopre_element_inner`](@ref) wrt q
"""
function calcSharedFaceIntegrals_nopre_element_inner_revq(
                            mesh::AbstractDGMesh{Tmsh},
                            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres},
                            opts, data::SharedFaceData,
                            data_bar::SharedFaceData,
                            functor_revq::FluxType_revq
                           ) where {Tmsh, Tsol, Tres}

  if getParallelData(data) != PARALLEL_DATA_ELEMENT
    error("element parallel data required")
  end

  q = eqn.q
  params = eqn.params

  # TODO: make these fields of params
  q_faceL = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceR = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceL_bar = Array{Tres}(mesh.numDofPerNode, mesh.numNodesPerFace)
#  q_faceR_bar = Array{Tres}(mesh.numDofPerNode, mesh.numNodesPerFace)
  qR_bar_k = Array{Tres}(mesh.numDofPerNode)  # throwaway array
#  q_faceR_bar = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)
  flux_face_bar = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerFace)

  idx = data.peeridx
  interfaces = data.interfaces
  bndries_local = data.bndries_local
  bndries_remote = data.bndries_remote

  # eqn.q_face_send should actually contain el
  # data, but we want to use eqn.q because that's the canonical source
  # qL_arr = eqn.q_face_send[i]
  qR_arr = data.q_recv
  resR_bar_arr = data_bar.q_recv
  nrm_arr = mesh.nrm_sharedface[idx]
  aux_vars_arr = eqn.aux_vars_sharedface[idx]

  start_elnum = mesh.shared_element_offsets[idx]

  for j=1:length(interfaces)
    iface_j = interfaces[j]
    bndryL_j = bndries_local[j]
    bndryR_j = bndries_remote[j]
    elL = iface_j.elementL
    elR = iface_j.elementR - start_elnum + 1
    fL = bndryL_j.face

    # interpolate to face
    qL = ro_sview(q, :, :, elL); qL_bar = sview(eqn.q_bar, :, :, elL)
    qR = ro_sview(qR_arr, :, :, elR);
    fill!(q_faceL_bar, 0);
    interiorFaceInterpolate!(mesh.sbpface, iface_j, qL, qR, q_faceL, q_faceR)

    fill!(flux_face_bar, 0)
    resL_bar = ro_sview(eqn.res_bar, :, :,  elL)
    resR_bar = ro_sview(resR_bar_arr, :, :, elR)
    interiorFaceIntegrate_rev!(mesh.sbpface, iface_j, flux_face_bar, resL_bar,
                               resR_bar, SummationByParts.Subtract())

    # calculate flux
    for k=1:mesh.numNodesPerFace
      qL_k = ro_sview(q_faceL, :, k); qL_bar_k = sview(q_faceL_bar, :, k)
      qR_k = ro_sview(q_faceR, :, k);
      aux_vars = ro_sview(aux_vars_arr, :, k, j)
      nrm_xy = ro_sview(nrm_arr, :, k, j)
      flux_bar_k = sview(flux_face_bar, :, k)

      parent(aux_vars)[1] = calcPressure(params, qL_k)

      functor_revq(params, qL_k, qL_bar_k, qR_k, qR_bar_k, aux_vars, nrm_xy,
                   flux_bar_k)
     end

     boundaryFaceInterpolate_rev!(mesh.sbpface, fL, qL_bar, q_faceL_bar)

   end  # end loop over interfaces

  return nothing
end

#------------------------------------------------------------------------------
# calcSharedFaceElementIntegrals_element

"""
  Differentiated version of [`calcSharedFaceElementIntegrals_element`](@ref)
"""
function calcSharedFaceElementIntegrals_element_diff(
                            mesh::AbstractDGMesh{Tmsh},
                            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres},
                            opts, data::SharedFaceData) where {Tmsh, Tsol, Tres}

  if opts["use_staggered_grid"]
    error("explicit computaton of jacobian not supported for staggered grids")

  else
    calcSharedFaceElementIntegrals_element_inner_diff(mesh, sbp, eqn, opts,
                            data, eqn.face_element_integral_func,  eqn.flux_func_diff, eqn.assembler)
  end

  return nothing
end

"""
  This function loops over given set of shared faces and computes a face
  integral that
  uses data from all volume nodes.  See [`FaceElementIntegralType`](@ref)
  for details on the integral performed.

  **Inputs**:

   * mesh
   * sbp
   * eqn
   * opts
   * data: a SharedFaceData specifying which shared faces to compute
   * face_integral_functor
   * flux_functor

"""
function calcSharedFaceElementIntegrals_element_inner_diff(
                            mesh::AbstractDGMesh{Tmsh},
                            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres},
                            opts, data::SharedFaceData,
                            face_integral_functor::FaceElementIntegralType,
                            flux_functor::FluxType_diff,
                            assembler::AssembleElementData) where {Tmsh, Tsol, Tres}

  if getParallelData(data) != PARALLEL_DATA_ELEMENT
    throw(ErrorException("cannot use calcSharedFaceIntegrals_element_diff without parallel element data"))
  end

  #TODO: consider writing specialized 1-sided flux functions that only do
  #      the local half.
  fdata = eqn.params.calc_face_integrals_data
  @unpack fdata res_jacLL res_jacLR res_jacRL res_jacRR

  q = eqn.q
  params = eqn.params
  # we don't care about elementR here, so use this throwaway array
  resR = Array{Tres}(mesh.numDofPerNode, mesh.numNodesPerElement)

  # get the data for the parallel interface
  idx = data.peeridx
  interfaces = data.interfaces
  bndries_local = data.bndries_local
  bndries_remote = data.bndries_remote
#    qL_arr = eqn.q_face_send[i]
  qR_arr = data.q_recv
  nrm_face_arr = mesh.nrm_sharedface[idx]

  start_elnum = mesh.shared_element_offsets[idx]

  # compute the integrals
  for j=1:length(interfaces)
    iface_j = interfaces[j]
    elL = iface_j.elementL
    elR = iface_j.elementR - start_elnum + 1  # is this always equal to j?
    qL = ro_sview(q, :, :, elL)
    qR = ro_sview(qR_arr, :, :, elR)
    aux_vars = ro_sview(eqn.aux_vars, :, :, elL)
    nrm_face = ro_sview(nrm_face_arr, :, :, j)

    fill!(res_jacLL, 0); fill!(res_jacLR, 0)
    fill!(res_jacRL, 0); fill!(res_jacRR, 0)
    calcFaceElementIntegral_diff(face_integral_functor, eqn.params, mesh.sbpface,
                            iface_j, qL, qR, aux_vars,
                            nrm_face, flux_functor, res_jacLL, res_jacLR, res_jacRL, res_jacRR)

    assembleSharedFace(assembler, mesh.sbpface, mesh, iface_j, res_jacLL, res_jacLR)
    # multiply by Minv if needed
    if params.use_Minv == 1
      for q=1:mesh.numNodesPerElement
        for p=1:mesh.numNodesPerElement
          valL = mesh.jac[p, iface_j.elementL]/sbp.w[p]  # entry in Minv
          @simd for m=1:mesh.numDofPerNode
            @simd for n=1:mesh.numDofPerNode
              res_jacLL[n, m, p, q] *= valL
              res_jacLR[n, m, p, q] *= valL
            end
          end
        end
      end
    end

  end  # end loop j

  return nothing
end

"""
  Reverse mode of [`calcSharedFaceElementIntegrals_element`](@ref) wrt metrics
"""
function calcSharedFaceElementIntegrals_element_revm(
                            mesh::AbstractDGMesh{Tmsh},
                            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres},
                            opts, data::SharedFaceData) where {Tmsh, Tsol, Tres}

  if opts["use_staggered_grid"]
    error("staggered grid not supported for revm product")
  else
    calcSharedFaceElementIntegrals_element_inner_revm(mesh, sbp, eqn, opts, data,
                            eqn.face_element_integral_func,  eqn.flux_func)
  end

  return nothing
end

"""
  Reverse mode of [calcSharedFaceElementIntegrals_element_inner`](@ref)

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * data
   * face_integral_functor
   * flux_functor: a `FluxType` object (not `FluxType_revm`)
"""
function calcSharedFaceElementIntegrals_element_inner_revm(
                            mesh::AbstractDGMesh{Tmsh},
                            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres},
                            opts, data::SharedFaceData,
                            face_integral_functor::FaceElementIntegralType,
                            flux_functor::FluxType) where {Tmsh, Tsol, Tres}

  if getParallelData(data) != PARALLEL_DATA_ELEMENT
    throw(ErrorException("cannot use calcSharedFaceIntegrals_element without parallel element data"))
  end

  q = eqn.q
  params = eqn.params
  # we don't care about elementR here, so use this throwaway array
  resR_bar = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerElement)

  # get the data for the parallel interface
  idx = data.peeridx
  interfaces = data.interfaces
#    qL_arr = eqn.q_face_send[i]
  qR_arr = data.q_recv
  nrm_face_arr = mesh.nrm_sharedface[idx]
  nrm_face_arr_bar = mesh.nrm_sharedface_bar[idx]

  start_elnum = mesh.shared_element_offsets[idx]

  # compute the integrals
  for j=1:length(interfaces)
    iface_j = interfaces[j]
    elL = iface_j.elementL
    elR = iface_j.elementR - start_elnum + 1
    qL = ro_sview(q, :, :, elL)
    qR = ro_sview(qR_arr, :, :, elR)
    aux_vars = ro_sview(eqn.aux_vars, :, :, elL)
    nrm_face = ro_sview(nrm_face_arr, :, :, j)
    nrm_bar = sview(nrm_face_arr_bar, :, :, j)
    resL_bar = ro_sview(eqn.res_bar, :, :, elL)

    calcFaceElementIntegral_revm(face_integral_functor, eqn.params, mesh.sbpface,
                            iface_j, qL, qR, aux_vars,
                            nrm_face, nrm_bar, flux_functor, resL_bar, resR_bar)
  end  # end loop j

  return nothing
end

"""
  Reverse mode of [`calcSharedFaceElementIntegrals_element`](@ref) wrt q
"""
function calcSharedFaceElementIntegrals_element_revq(
                            mesh::AbstractDGMesh{Tmsh},
                            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres},
                            opts, data::SharedFaceData,
                            data_bar::SharedFaceData) where {Tmsh, Tsol, Tres}

  if opts["use_staggered_grid"]
    error("staggered grid not supported for revq product")
  else
    calcSharedFaceElementIntegrals_element_inner_revq(mesh, sbp, eqn, opts,
          data, data_bar,  eqn.face_element_integral_func,  eqn.flux_func_revq)
  end

  return nothing
end

#=
"""
  Reverse mode of [calcSharedFaceElementIntegrals_element_inner`](@ref)

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
   * data
   * data_q_bar: a `SharedFaceData`.  `q_recv`, not `q_send` will be updated
                 with the back-propigation of `eqn.res_bar`
   * face_integral_functor
   * flux_functor: a `FluxType_revq` object
"""
function calcSharedFaceElementIntegrals_element_inner_revq(
                            mesh::AbstractDGMesh{Tmsh},
                            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres},
                            opts,
                            data::SharedFaceData, data_q_bar::SharedFaceData,
                            face_integral_functor::FaceElementIntegralType,
                            flux_functor_revq::FluxType_revq
                            ) where {Tmsh, Tsol, Tres}

  if getParallelData(data) != PARALLEL_DATA_ELEMENT || getParallelData(data_q_bar) != PARALLEL_DATA_ELEMENT
    throw(ErrorException("cannot use calcSharedFaceIntegrals_element without parallel element data"))
  end

  q = eqn.q
  params = eqn.params
  # we don't care about elementR here, so use this throwaway array
  resR_bar = zeros(Tres, mesh.numDofPerNode, mesh.numNodesPerElement)

  # get the data for the parallel interface
  idx = data.peeridx
  interfaces = data.interfaces
#    qL_arr = eqn.q_face_send[i]
  qR_arr = data.q_recv
  qR_arr_bar = data_q_bar.q_recv
  nrm_face_arr = mesh.nrm_sharedface[idx]

  start_elnum = mesh.shared_element_offsets[idx]

  # compute the integrals
  for j=1:length(interfaces)
    iface_j = interfaces[j]
    elL = iface_j.elementL
    elR = iface_j.elementR - start_elnum + 1
    qL = ro_sview(q, :, :, elL);      qL_bar = sview(eqn.q_bar, :, :, elL)
    qR = ro_sview(qR_arr, :, :, elR); qR_bar = sview(qR_arr_bar, :, :, elR)
    aux_vars = ro_sview(eqn.aux_vars, :, :, elL)
    nrm_face = ro_sview(nrm_face_arr, :, :, j)
    resL_bar = ro_sview(eqn.res_bar, :, :, elL)

    calcFaceElementIntegral_revq(face_integral_functor, params, mesh.sbpface,
                            iface_j, qL, qL_bar, qR,  qR_bar, aux_vars,
                            nrm_face, flux_functor_revq, resL_bar, resR_bar)
  end  # end loop j

  return nothing
end
=#

"""
  Reverse mode of [`calcSharedFaceElementIntegrals_element`](@ref) wrt q
"""
function calcSharedFaceElementIntegrals_element_inner_revq(
                            mesh::AbstractDGMesh{Tmsh},
                            sbp::AbstractOperator, eqn::EulerData{Tsol, Tres},
                            opts,
                            data::SharedFaceData, data_q_bar::SharedFaceData,
                            face_integral_functor::FaceElementIntegralType,
                            flux_functor_revq::FluxType_revq
                            ) where {Tmsh, Tsol, Tres}

  if getParallelData(data) != PARALLEL_DATA_ELEMENT || getParallelData(data_q_bar) != PARALLEL_DATA_ELEMENT
    throw(ErrorException("cannot use calcSharedFaceIntegrals_element without parallel element data"))
  end

  q = eqn.q
  params = eqn.params
  qR_bar = Array{Tres}(mesh.numDofPerNode, mesh.numNodesPerElement) # throwaway

  # get the data for the parallel interface
  idx = data.peeridx
  interfaces = data.interfaces
#    qL_arr = eqn.q_face_send[i]
  qR_arr = data.q_recv
  resR_arr_bar = data_q_bar.q_recv
  nrm_face_arr = mesh.nrm_sharedface[idx]

  start_elnum = mesh.shared_element_offsets[idx]

  # compute the integrals
  for j=1:length(interfaces)
    iface_j = interfaces[j]
    elL = iface_j.elementL
    elR = iface_j.elementR - start_elnum + 1
    qL = ro_sview(q, :, :, elL);      qL_bar = sview(eqn.q_bar, :, :, elL)
    qR = ro_sview(qR_arr, :, :, elR);
    aux_vars = ro_sview(eqn.aux_vars, :, :, elL)
    nrm_face = ro_sview(nrm_face_arr, :, :, j)
    resL_bar = ro_sview(eqn.res_bar, :, :, elL)
    resR_bar = ro_sview(resR_arr_bar, :, :, elR)

    calcFaceElementIntegral_revq(face_integral_functor, params, mesh.sbpface,
                            iface_j, qL, qL_bar, qR,  qR_bar, aux_vars,
                            nrm_face, flux_functor_revq, resL_bar, resR_bar)
  end  # end loop j

  return nothing
end

#------------------------------------------------------------------------------
# Functors

"""
  Calls the [`RoeSolver_diff`](@ref)
"""
mutable struct RoeFlux_diff <: FluxType_diff
end

function (obj::RoeFlux_diff)(params::ParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractVector{Tmsh},
              F_dotL::AbstractArray{Tres},
              F_dotR::AbstractArray{Tres}) where {Tsol, Tres, Tmsh}

  RoeSolver_diff(params, uL, uR, aux_vars, nrm, F_dotL, F_dotR)
  return nothing
end

mutable struct LFFlux_diff <: FluxType_diff
end

function (obj::LFFlux_diff)(params::ParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractVector{Tmsh},
              F_dotL::AbstractArray{Tres},
              F_dotR::AbstractArray{Tres}) where {Tsol, Tres, Tmsh}

  calcLFFlux_diff(params, uL, uR, aux_vars, nrm, F_dotL, F_dotR)

  return nothing
end

"""
  Calls the [`calcEulerFlux_IR_diff`](@ref)
"""
mutable struct IRFlux_diff <: FluxType_diff
end

function (obj::IRFlux_diff)(params::ParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractArray{Tmsh},
              F_dotL::AbstractArray{Tres},
              F_dotR::AbstractArray{Tres}) where {Tsol, Tres, Tmsh}

  calcEulerFlux_IR_diff(params, uL, uR, aux_vars, nrm, F_dotL, F_dotR)

  return nothing
end

mutable struct StandardFlux_diff <: FluxType_diff
end

function (obj::StandardFlux_diff)(params::ParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractVector{Tmsh},
              F_dotL::AbstractArray{Tres},
              F_dotR::AbstractArray{Tres}) where {Tsol, Tres, Tmsh}

  # this should be implemented eventually, but for now this method needs to
  # exist as a stub
  error("StandardFlux() called, something unexpected has occured")

end

mutable struct IRSLFFlux_diff <: FluxType_diff
end

function (obj::IRSLFFlux_diff)(params::ParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractVector{Tmsh},
              F_dotL::AbstractArray{Tres},
              F_dotR::AbstractArray{Tres}) where {Tsol, Tres, Tmsh}

  calcEulerFlux_IRSLF_diff(params, uL, uR, aux_vars, nrm, F_dotL, F_dotR)
end

mutable struct ErrorFlux_diff <: FluxType_diff
end

function (obj::ErrorFlux_diff)(params::ParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractVector{Tmsh},
              F_dotL::AbstractArray{Tres},
              F_dotR::AbstractArray{Tres}) where {Tsol, Tres, Tmsh}

  error("ErrorFlux() called, something unexpected has occured")

end

mutable struct HLLFlux_diff <: FluxType_diff
end

function (obj::HLLFlux_diff)(params::ParamType,
              uL::AbstractArray{Tsol,1},
              uR::AbstractArray{Tsol,1},
              aux_vars::AbstractVector{Tres},
              nrm::AbstractVector{Tmsh},
              F_dotL::AbstractArray{Tres},
              F_dotR::AbstractArray{Tres}) where {Tsol, Tres, Tmsh}

  calcHLLFlux_diff(params, uL, uR, aux_vars, nrm, F_dotL, F_dotR)
end

"""
  Container for all differentiated flux functors.  Maps name to object.
  The names are exactly the same as the non-differentiated functor.
"""
global const FluxDict_diff = Dict{String, FluxType_diff}(
"RoeFlux" => RoeFlux_diff(),
"LFFlux" => LFFlux_diff(),
"IRFlux" => IRFlux_diff(),
"IRSLFFlux" => IRSLFFlux_diff(),
"StandardFlux" => StandardFlux_diff(),
"HLLFlux" => HLLFlux_diff(),
"ErrorFlux" => ErrorFlux_diff(),
)

"""
  Gets the functor for the differentiated version of the flux function and
  stores is to `eqn.flux_func_diff`

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts

  **Options Keys**

  Uses the `Flux_name` key to get the functor (same as for the non-differentiated
  version).

  This function only gets the functor if the jacobian will be computed explicitly
  as specified by `calc_jac_explicit`, otherwise it uses a functor that
  will throw an error if used.
"""
function getFluxFunctors_diff(mesh::AbstractDGMesh, sbp, eqn, opts)

  if opts["calc_jac_explicit"]
    name = opts["Flux_name"]
    name2 = opts["Volume_flux_name"]
  else
    name = "ErrorFlux"
    name2 = "ErrorFlux"
  end

  eqn.flux_func_diff = FluxDict_diff[name]
  eqn.volume_flux_func_diff = FluxDict_diff[name2]

  return nothing
end


text:
"""
General Interface for Computation of variance interval of price

		variance_=variance(mcProcess,rfCurve,mcBaseData,payoff_,)

Where:\n
		mcProcess          = Process to be simulated.
		rfCurve  = Zero Rate Data.
		mcBaseData = Basic properties of MonteCarlo simulation
		payoff_ = Payoff(s) to be priced

		variance_     = variance of the payoff of the derivative

"""
function variance(mcProcess::BaseProcess, rfCurve::AbstractZeroRateCurve, mcConfig::MonteCarloConfiguration, abstractPayoff::AbstractPayoff)
    set_seed(mcConfig)
    T = maturity(abstractPayoff)
    S = simulate(mcProcess, rfCurve, mcConfig, T)
    Payoff = payoff(S, abstractPayoff, rfCurve)
    variance_ = var(Payoff)
    return variance_
end

function variance(mcProcess::BaseProcess, rfCurve::AbstractZeroRateCurve, mcConfig::MonteCarloConfiguration, abstractPayoffs::Array{abstractPayoff_}) where {abstractPayoff_ <: AbstractPayoff}
    set_seed(mcConfig)
    maxT = maximum([maturity(abstractPayoff) for abstractPayoff in abstractPayoffs])
    S = simulate(mcProcess, rfCurve, mcConfig, maxT)
    variance_ = [var(payoff(S, abstractPayoff, rfCurve, maxT)) for abstractPayoff in abstractPayoffs]

    return variance_
end

ext: jl


text:
using Baytes
using Documenter

DocMeta.setdocmeta!(Baytes, :DocTestSetup, :(using Baytes); recursive=true)

makedocs(;
    modules=[Baytes],
    authors="Patrick Aschermayr <p.aschermayr@gmail.com>",
    repo="https://github.com/paschermayr/Baytes.jl/blob/{commit}{path}#{line}",
    sitename="Baytes.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://paschermayr.github.io/Baytes.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Introduction" => "intro.md",
    ],
)

deploydocs(;
    repo="github.com/paschermayr/Baytes.jl",
    devbranch="main",
)


text:
struct TightenedNMDT <: AbstractRelaxation
    num_layers::Int
end

function reformulate_unary_quadratic_term!(
    model::MOI.ModelLike,
    reformulation::Reformulation{TightenedNMDT},
    _x::VI,
    _y::VI,
    _x_l::Float64,
    _x_u::Float64,
)
    x, y = _map_to_interval_if_needed(model, _x, _y, _x_l, _x_u)
    x_l, x_u = 0.0, 1.0

    L = reformulation.relaxation.num_layers
    α_v, α_c =
        MOI.add_constrained_variables(model, [MOI.ZeroOne() for i in 1:L])
    Δx_v, Δx_c =
        MOI.add_constrained_variable(model, MOI.Interval(0.0, 2.0^(-L)))
    u_v = MOI.add_variables(model, L)
    Δz_v = MOI.add_variable(model)

    α = [SV(α_v[i]) for i in 1:L]
    u = [SV(u_v[i]) for i in 1:L]
    x = 1.0 * SV(x)
    y = 1.0 * SV(y)
    Δx = SV(Δx_v)
    Δz = SV(Δz_v)

    MOI.add_constraint(
        model,
        x - Δx - sum(2.0^(-i) * α[i] for i in 1:L),
        MOI.EqualTo(0.0),
    )
    MOI.add_constraint(
        model,
        y - Δz - sum(2.0^(-i) * u[i] for i in 1:L),
        MOI.EqualTo(0.0),
    )
    for i in 1:L
        a, b, c = x + Δx, α[i], u[i]
        a_l, a_u = x_l, x_u + 2.0^(-L)
        b_l, b_u = 0.0, 1.0

        MOIU.normalize_and_add_constraint(
            model,
            c - (a_u * b + a * b_u - a_u * b_u),
            MOI.GreaterThan(0.0),
        )
        MOIU.normalize_and_add_constraint(
            model,
            c - (a_l * b + a * b_l - a_l * b_l),
            MOI.GreaterThan(0.0),
        )
        MOIU.normalize_and_add_constraint(
            model,
            c - (a_u * b + a * b_l - a_u * b_l),
            MOI.LessThan(0.0),
        )
        MOIU.normalize_and_add_constraint(
            model,
            c - (a_l * b + a * b_u - a_l * b_u),
            MOI.LessThan(0.0),
        )
    end
    MOI.add_constraint(model, Δz - 2.0^(-L) * Δx, MOI.LessThan(0.0))
    MOI.add_constraint(model, Δz, MOI.GreaterThan(0.0))
    MOIU.normalize_and_add_constraint(
        model,
        Δz - 2.0^(-L + 1) * (Δx - 2.0^(-L)) - 2.0^(-2L),
        MOI.GreaterThan(0.0),
    )
    # Only adding one cut corresponding to x̃_1 = 2^(-L-1)
    xt_1 = 2.0^(-L - 1)
    MOIU.normalize_and_add_constraint(
        model,
        Δz - 2.0 * xt_1 * (Δx - xt_1) - xt_1^2,
        MOI.GreaterThan(0.0),
    )
    return α_v
end


text:
module problem103

export SpecialSet, can_add, add

# Special set is a set of values, maintaining the two properties:
# 1) Every disjoint subset has a unique sum
# 2) every subset size N has greater sum than all subsets size < N
# The values member contains the items in the special set
# The sums member tracks the unique sums for all possible subsets
type SpecialSet
  values :: Vector{Int}
  sums :: Vector{Int}
end

# add the value x to SpecialSet s, assumes can_add(x, s) is true
function add(x :: Int, s :: SpecialSet)
  return SpecialSet([x;s.values], [x;s.sums;x+s.sums])
end

# Check if value x can be added to SpecialSet s, maintaining the SpecialSet invariants
function can_add(x :: Int, s :: SpecialSet)

  # must be a new minimum value
  if x >= s.values[1]
    #println("can't add $x to $(s.values'), not a minimum")
    return false
  end

  # x + any existing sum should be unique
  # implies that the sum of any subset containing x would be unique
  for y in s.sums
    if in(x+y, s.sums)
      #println("can't add $x to $(s.values'), sum $(x+y) already exists")
      return false
    end
  end

  # larger sized subsets must sum to greater than smaller subsets
  for k = 1:ceil(length(s.values)/2)
    if x+sum(s.values[1:k]) <= sum(s.values[end-k+1:end])
      return false
    end
  end

  return true
end

# Generate all SpecialSets having minimum value=xmin, max value=xmax, length=n
function specials(xmin :: Int, xmax :: Int, n :: Int)

  # base case: all length 2 sets are SpecialSets
  if n == 2
    return [SpecialSet([xmin, xmax], [xmin, xmax, xmin+xmax])]
  end

  # recursive case: Find all special sets from mid to max, length n-1
  # Then check if xmin can be added to the specialset
  results = SpecialSet[]
  mid = max(xmin + 1, xmax - xmin + 1)
  for x = mid : xmax-1
    for s in specials(x, xmax, n-1)
      if can_add(xmin, s)
        push!(results, add(xmin, s))
      end
    end
  end

  return results
end

# returns the values vector of the SpecialSet with minimum sum
# having min, max and length specified
function optimum_special_set(xmin :: Int, xmax :: Int, n :: Int)
  sets = specials(xmin, xmax, n)
  totals = [sum(s.values) for s in sets]
  t, i = findmin(totals)
  return sets[i].values
end

using Base.Test

@test optimum_special_set(3, 7, 4) == [3, 5, 6, 7]
@test optimum_special_set(6, 13, 5) == [6, 9, 11, 12, 13]
@test optimum_special_set(11, 25, 6) == [11, 18, 19, 20, 22, 25]

# based on the hints in the question, the bounds of the solution should be 20 and 45
#println(optimum_special_set(20, 45, 7)')

end

max_stars_repo_name: mbuhot/mbuhot-euler-solutions

max_issues_repo_name: mbuhot/mbuhot-euler-solutions

max_forks_repo_name: mbuhot/mbuhot-euler-solutions


text:
# fourier series
# Pkg.add("Javis")

using Colors, Luxor
function ground(args...)
    background("black")
    sethue("white")
    Javis.translate(-args[1].width/2+50, 0)
    Javis.scale(700, 150)
end

function wave(xs, ys, opacity, color)
    setline(1.5)
    setopacity(opacity)
    sethue(color)
    points = Vector{Luxor.Point}(undef, length(xs))
    for (x, y, i) in zip(xs, ys, 1:length(xs))
        points[i] = Luxor.Point(x, y)
    end
    Javis.prettypoly(points, :stroke, ()->())
end

function term(x, i)
    return 4/π * sin.(2π*(2i-1)*x)/(2i-1)
end

function sum_term(x, k)
    k == 0 && return zeros(length(x))
    return sum(term(x, i) for i in 1:k)
end

nframes = 300
frames_per_wave = 40
video = Video(800, 600)
Background(1:nframes, ground)

x = 0.0:0.001:1.0
k = 6
colors = [RGB(0.0, 1.0, 0.4), RGB(0, 1.0, 1.0), RGB(1.0, 0.25, 0.25), RGB(1.0, 1.0, 0.0), RGB(1.0, 0.5, 1.0), RGB(0.75, 0.75, 1.0)]
waves = Object[]
for i = 1:k
    frames = frames_per_wave*(i-1)+1:nframes
    push!(waves, Object(frames, (args...; y=term(x,i)) -> wave(x, y, 0.5, colors[i])))
    act!(waves[end], Action(5:frames_per_wave,
                    change(:y, term(x,i) => sum_term(x, i))))
end

sum_wave = Object(1:nframes, (args...; y=zeros(length(x)))->wave(x, y, 1.0, "white"))
for i = 1:k
    act!(sum_wave, Action(frames_per_wave*(i-1)+1:frames_per_wave*i,
                         change(:y, sum_term(x, i-1) => sum_term(x, i))))
end

render(video; pathname="fourier_1D.gif")

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
using LinearAlgebra

export fromroots,
       truncate!,
       chop!,
       coeffs,
       degree,
       domain,
       mapdomain,
       order,
       hasnan,
       roots,
       companion,
       vander,
       fit,
       integrate,
       derivative,
       variable,
       isintegral,
       ismonic

"""
    fromroots(::AbstractVector{<:Number}; var=:x)
    fromroots(::Type{<:AbstractPolynomial}, ::AbstractVector{<:Number}; var=:x)

Construct a polynomial of the given type given the roots. If no type is given, defaults to `Polynomial`.

# Examples
```jldoctest common
julia> using Polynomials

julia> r = [3, 2]; # (x - 3)(x - 2)

julia> fromroots(r)
Polynomial(6 - 5*x + x^2)
```
"""
function fromroots(P::Type{<:AbstractPolynomial}, roots::AbstractVector; var::SymbolLike = :x)
    x = variable(P, var)
    p =  prod(x - r for r in roots)
    return truncate!(p)
end
fromroots(r::AbstractVector{<:Number}; var::SymbolLike = :x) =
    fromroots(Polynomial, r, var = var)

"""
    fromroots(::AbstractMatrix{<:Number}; var=:x)
    fromroots(::Type{<:AbstractPolynomial}, ::AbstractMatrix{<:Number}; var=:x)

Construct a polynomial of the given type using the eigenvalues of the given matrix as the roots. If no type is given, defaults to `Polynomial`.

# Examples
```jldoctest common
julia> using Polynomials

julia> A = [1 2; 3 4]; # (x - 5.37228)(x + 0.37228)

julia> fromroots(A)
Polynomial(-1.9999999999999998 - 5.0*x + 1.0*x^2)
```
"""
fromroots(P::Type{<:AbstractPolynomial},
    A::AbstractMatrix{T};
    var::SymbolLike = :x,) where {T <: Number} = fromroots(P, eigvals(A), var = var)
fromroots(A::AbstractMatrix{T}; var::SymbolLike = :x) where {T <: Number} =
    fromroots(Polynomial, eigvals(A), var = var)

"""
    fit(x, y, deg=length(x) - 1; [weights], var=:x)
    fit(::Type{<:AbstractPolynomial}, x, y, deg=length(x)-1; [weights], var=:x)

Fit the given data as a polynomial type with the given degree. Uses linear least squares. When weights are given, as either a `Number`, `Vector` or `Matrix`, will use weighted linear least squares. The default polynomial type is [`Polynomial`](@ref). This will automatically scale your data to the [`domain`](@ref) of the polynomial type using [`mapdomain`](@ref)
"""
function fit(P::Type{<:AbstractPolynomial},
             x::AbstractVector{T},
             y::AbstractVector{T},
             deg::Integer = length(x) - 1;
    weights = nothing,
    var = :x,) where {T}
    x = mapdomain(P, x)
    vand = vander(P, x, deg)
    if weights !== nothing
        coeffs = _wlstsq(vand, y, weights)
    else
        coeffs = pinv(vand) * y
    end
    return P(T.(coeffs), var)
end

fit(P::Type{<:AbstractPolynomial},
    x,
    y,
    deg::Integer = length(x) - 1;
    weights = nothing,
    var = :x,) = fit′(P, promote(collect(x), collect(y))..., deg; weights = weights, var = var)

#  avoid issue  214
fit′(P::Type{<:AbstractPolynomial}, x, y, args...;kwargs...) = throw(MethodError("x and y do not produce abstract   vectors"))
fit′(P::Type{<:AbstractPolynomial},
     x::AbstractVector{T},
     y::AbstractVector{T},
     args...; kwargs...) where {T} = fit(P,x,y,args...;  kwargs...)

fit(x::AbstractVector,
    y::AbstractVector,
    deg::Integer = length(x) - 1;
    weights = nothing,
    var = :x,) = fit(Polynomial, x, y, deg; weights = weights, var = var)

# Weighted linear least squares
_wlstsq(vand, y, W::Number) = _wlstsq(vand, y, fill!(similar(y), W))
_wlstsq(vand, y, W::AbstractVector) = _wlstsq(vand, y, diagm(0 => W))
_wlstsq(vand, y, W::AbstractMatrix) = (vand' * W * vand) \ (vand' * W * y)

"""
    roots(::AbstractPolynomial; kwargs...)

Returns the roots of the given polynomial. This is calculated via the eigenvalues of the companion matrix. The `kwargs` are passed to the `LinearAlgeebra.eigvals` call.

!!! note

        The [PolynomialRoots.jl](https://github.com/giordano/PolynomialRoots.jl) package provides an alternative that is a bit faster and a bit more accurate; the [FastPolynomialRoots](https://github.com/andreasnoack/FastPolynomialRoots.jl) provides an interface to FORTRAN code implementing an algorithm that can handle very large polynomials (it is  `O(n^2)` not `O(n^3)`. the [AMRVW.jl](https://github.com/jverzani/AMRVW.jl) package implements the algorithm in Julia, allowing the use of other  number types.

"""
function roots(q::AbstractPolynomial{T}; kwargs...) where {T <: Number}

    p = convert(Polynomial{T},  q)
    roots(p; kwargs...)

end

"""
    companion(::AbstractPolynomial)

Return the companion matrix for the given polynomial.

# References
[Companion Matrix](https://en.wikipedia.org/wiki/Companion_matrix)
"""
companion(::AbstractPolynomial)

"""
    vander(::Type{AbstractPolynomial}, x::AbstractVector, deg::Integer)

Calculate the psuedo-Vandermonde matrix of the given polynomial type with the given degree.

# References
[Vandermonde Matrix](https://en.wikipedia.org/wiki/Vandermonde_matrix)
"""
vander(::Type{<:AbstractPolynomial}, x::AbstractVector, deg::Integer)

"""
    integrate(::AbstractPolynomial, C=0)

Returns the indefinite integral of the polynomial with constant `C`.
"""
integrate(p::AbstractPolynomial, C::Number = 0) = integrate(p, C)

"""
    integrate(::AbstractPolynomial, a, b)

Compute the definite integral of the given polynomial from `a` to `b`. Will throw an error if either `a` or `b` are out of the polynomial's domain.
"""
function integrate(p::AbstractPolynomial, a::Number, b::Number)
    P = integrate(p)
    return P(b) - P(a)
end

"""
    derivative(::AbstractPolynomial, order::Int = 1)

Returns a polynomial that is the `order`th derivative of the given polynomial. `order` must be non-negative.
"""
derivative(::AbstractPolynomial, ::Int)

"""
    truncate!(::AbstractPolynomial{T};
        rtol::Real = Base.rtoldefault(real(T)), atol::Real = 0)

In-place version of [`truncate`](@ref)
"""
function truncate!(p::AbstractPolynomial{T};
    rtol::Real = Base.rtoldefault(real(T)),
                   atol::Real = 0,) where {T}
    max_coeff = maximum(abs, coeffs(p))
    thresh = max_coeff * rtol + atol
    map!(c->abs(c) <= thresh ? zero(T) : c, coeffs(p), coeffs(p))
    return chop!(p, rtol = rtol, atol = atol)
end

"""
    truncate(::AbstractPolynomial{T};
        rtol::Real = Base.rtoldefault(real(T)), atol::Real = 0)

Rounds off coefficients close to zero, as determined by `rtol` and `atol`, and then chops any leading zeros. Returns a new polynomial.
"""
function Base.truncate(p::AbstractPolynomial{T};
    rtol::Real = Base.rtoldefault(real(T)),
    atol::Real = 0,) where {T}
    truncate!(deepcopy(p), rtol = rtol, atol = atol)
end

"""
    chop!(::AbstractPolynomial{T};
        rtol::Real = Base.rtoldefault(real(T)), atol::Real = 0))

In-place version of [`chop`](@ref)
"""
function chop!(p::AbstractPolynomial{T};
    rtol::Real = Base.rtoldefault(real(T)),
               atol::Real = 0,) where {T}
    isempty(coeffs(p)) && return p
    tol = norm(coeffs(p)) * rtol + atol
    for i = lastindex(p):-1:0
        val = p[i]
        if abs(val) > tol #!isapprox(val, zero(T); rtol = rtol, atol = atol)
            resize!(p.coeffs, i + 1);
            return p
        end
    end
    resize!(p.coeffs, 1)
    return p
end

"""
    chop(::AbstractPolynomial{T};
        rtol::Real = Base.rtoldefault(real(T)), atol::Real = 0))

Removes any leading coefficients that are approximately 0 (using `rtol` and `atol`). Returns a polynomial whose degree will guaranteed to be equal to or less than the given polynomial's.
"""
function Base.chop(p::AbstractPolynomial{T};
    rtol::Real = Base.rtoldefault(real(T)),
    atol::Real = 0,) where {T}
    chop!(deepcopy(p), rtol = rtol, atol = atol)
end

"""
    check_same_variable(p::AbstractPolynomial, q::AbstractPolynomial)

Check if either `p` or `q` is constant or if `p` and `q` share the same variable
"""
check_same_variable(p::AbstractPolynomial, q::AbstractPolynomial) =
    (Polynomials.isconstant(p) || Polynomials.isconstant(q)) || p.var ==  q.var

#=
Linear Algebra =#
"""
    norm(::AbstractPolynomial, p=2)

Calculates the p-norm of the polynomial's coefficients
"""
LinearAlgebra.norm(q::AbstractPolynomial, p::Real = 2) = norm(coeffs(q), p)

"""
    conj(::AbstractPolynomial)

Returns the complex conjugate of the polynomial
"""
LinearAlgebra.conj(p::P) where {P <: AbstractPolynomial} = ⟒(P)(conj(coeffs(p)), p.var)
LinearAlgebra.adjoint(p::P) where {P <: AbstractPolynomial} = ⟒(P)(adjoint.(coeffs(p)), p.var)
LinearAlgebra.transpose(p::AbstractPolynomial) = p
LinearAlgebra.transpose!(p::AbstractPolynomial) = p

#=
Conversions =#
Base.convert(::Type{P}, p::P) where {P <: AbstractPolynomial} = p
Base.convert(P::Type{<:AbstractPolynomial}, x) = P(x)
Base.promote_rule(::Type{<:AbstractPolynomial{T}},
    ::Type{<:AbstractPolynomial{S}},
) where {T,S} = Polynomial{promote_type(T, S)}

#=
Inspection =#
"""
    length(::AbstractPolynomial)

The length of the polynomial.
"""
Base.length(p::AbstractPolynomial) = length(coeffs(p))
"""
    size(::AbstractPolynomial, [i])

Returns the size of the polynomials coefficients, along axis `i` if provided.
"""
Base.size(p::AbstractPolynomial) = size(coeffs(p))
Base.size(p::AbstractPolynomial, i::Integer) = size(coeffs(p), i)
Base.eltype(p::AbstractPolynomial{T}) where {T} = T
# in  analogy  with  polynomial as a Vector{T} with different operations defined.
Base.eltype(::Type{<:AbstractPolynomial}) = Float64
Base.eltype(::Type{<:AbstractPolynomial{T}}) where {T} = T
#Base.eltype(::Type{P}) where {P <: AbstractPolynomial} = P # changed  in v1.1.0
function Base.iszero(p::AbstractPolynomial)
    if length(p) == 0
        return true
    end
    return all(iszero.(coeffs(p))) && p[0] == 0
end

# See discussions in https://github.com/JuliaMath/Polynomials.jl/issues/258
"""
    all(pred, poly::AbstractPolynomial)

Test whether all coefficients of an `AbstractPolynomial` satisfy predicate `pred`.

You can implement `isreal`, etc., to a `Polynomial` by using `all`.
"""
Base.all(pred, poly::AbstractPolynomial) = all(pred, poly[:])
"""
    any(pred, poly::AbstractPolynomial)

Test whether any coefficient of an `AbstractPolynomial` satisfies predicate `pred`.
"""
Base.any(pred, poly::AbstractPolynomial) = any(pred, poly[:])
"""
    map(fn, p::AbstractPolynomial)

Transform coefficients of `p` by applying a function (or other callables) `fn` to each of them.

You can implement `real`, etc., to a `Polynomial` by using `map`.
"""
Base.map(fn, p::P) where {P<:AbstractPolynomial} = ⟒(P)(map(fn, coeffs(p)), p.var)

"""
    isreal(p::AbstractPolynomial)

Determine whether a polynomial is a real polynomial, i.e., having only real numbers as coefficients.

See also: [`real`](@ref)
"""
Base.isreal(p::AbstractPolynomial) = all(isreal, p)
"""
    real(p::AbstractPolynomial)

Construct a real polynomial from the real parts of the coefficients of `p`.

See also: [`isreal`](@ref)

!!! note
    This could cause losing terms in `p`. This method is usually called on polynomials like `p = Polynomial([1, 2 + 0im, 3.0, 4.0 + 0.0im])` where you want to chop the imaginary parts of the coefficients of `p`.
"""
Base.real(p::AbstractPolynomial) = map(real, p)

"""
    isintegral(p::AbstractPolynomial)

Determine whether a polynomial is an integer polynomial, i.e., having only integers as coefficients.
"""
isintegral(p::AbstractPolynomial) = all(isinteger, p)

"""
    ismonic(p::AbstractPolynomial)

Determine whether a polynomial is a monic polynomial, i.e., its leading coefficient is one.
"""
ismonic(p::AbstractPolynomial) = isone(p[end])

"""
    coeffs(::AbstractPolynomial)

Return the coefficient vector `[a_0, a_1, ..., a_n]` of a polynomial.
"""
coeffs(p::AbstractPolynomial) = p.coeffs

"""
    degree(::AbstractPolynomial)

Return the degree of the polynomial, i.e. the highest exponent in the polynomial that
has a nonzero coefficient. The degree of the zero polynomial is defined to be -1.
"""
degree(p::AbstractPolynomial) = iszero(p) ? -1 : length(p) - 1

"""
    isconstant(::AbstractPolynomial)

Is the polynomial  `p` a constant.
"""
isconstant(p::AbstractPolynomial) = degree(p) <= 0

hasnan(p::AbstractPolynomial) = any(isnan.(coeffs(p)))

"""
    domain(::Type{<:AbstractPolynomial})

Returns the domain of the polynomial.
"""
domain(::Type{<:AbstractPolynomial})
domain(::P) where {P <: AbstractPolynomial} = domain(P)

"""
    mapdomain(::Type{<:AbstractPolynomial}, x::AbstractArray)
    mapdomain(::AbstractPolynomial, x::AbstractArray)

Given values of x that are assumed to be unbounded (-∞, ∞), return values rescaled to the domain of the given polynomial.

# Examples
```jldoctest  common
julia> using Polynomials

julia> x = -10:10
-10:10

julia> extrema(mapdomain(ChebyshevT, x))
(-1.0, 1.0)

```
"""
function mapdomain(P::Type{<:AbstractPolynomial}, x::AbstractArray)
    d = domain(P)
    x = collect(x)
    x_zerod = x .- minimum(x)
    x_scaled = x_zerod .* (last(d) - first(d)) ./ maximum(x_zerod)
    x_scaled .+= first(d)
    return x_scaled
end
mapdomain(::P, x::AbstractArray) where {P <: AbstractPolynomial} = mapdomain(P, x)
#=
indexing =#
Base.firstindex(p::AbstractPolynomial) = 0
Base.lastindex(p::AbstractPolynomial) = length(p) - 1
Base.eachindex(p::AbstractPolynomial) = 0:length(p) - 1
Base.broadcastable(p::AbstractPolynomial) = Ref(p)

# iteration
# iteration occurs over the basis polynomials
Base.iterate(p::AbstractPolynomial) = (p[0] * one(typeof(p)), 1)
function Base.iterate(p::AbstractPolynomial, state)
    state <= length(p) - 1 ? (p[state] * basis(p, state), state + 1) : nothing
end

Base.collect(p::P) where {P <: AbstractPolynomial} = collect(P, p)

# getindex
function Base.getindex(p::AbstractPolynomial{T}, idx::Int) where {T <: Number}
    idx < 0 && throw(BoundsError(p, idx))
    idx ≥ length(p) && return zero(T)
    return coeffs(p)[idx + 1]
end
Base.getindex(p::AbstractPolynomial, idx::Number) = getindex(p, convert(Int, idx))
Base.getindex(p::AbstractPolynomial, indices) = [getindex(p, i) for i in indices]
Base.getindex(p::AbstractPolynomial, ::Colon) = coeffs(p)

# setindex
function Base.setindex!(p::AbstractPolynomial, value::Number, idx::Int)
    n = length(coeffs(p))
    if n ≤ idx
        resize!(p.coeffs, idx + 1)
        p.coeffs[n + 1:idx] .= 0
    end
    p.coeffs[idx + 1] = value
    return p
end

Base.setindex!(p::AbstractPolynomial, value::Number, idx::Number) =
    setindex!(p, value, convert(Int, idx))
Base.setindex!(p::AbstractPolynomial, value::Number, indices) =
    [setindex!(p, value, i) for i in indices]
Base.setindex!(p::AbstractPolynomial, values, indices) =
    [setindex!(p, v, i) for (v, i) in zip(values, indices)]
Base.setindex!(p::AbstractPolynomial, value::Number, ::Colon) =
    setindex!(p, value, eachindex(p))
Base.setindex!(p::AbstractPolynomial, values, ::Colon) =
    [setindex!(p, v, i) for (v, i) in zip(values, eachindex(p))]

#=
identity =#
Base.copy(p::P) where {P <: AbstractPolynomial} = P(copy(coeffs(p)), p.var)
Base.hash(p::AbstractPolynomial, h::UInt) = hash(p.var, hash(coeffs(p), h))

#=
zero, one, variable, basis =#
"""
    zero(::Type{<:AbstractPolynomial})
    zero(::AbstractPolynomial)

Returns a representation of 0 as the given polynomial.
"""
Base.zero(::Type{P}, var=:x) where {P <: AbstractPolynomial} = ⟒(P)(zeros(eltype(P), 1), var)
Base.zero(p::P) where {P <: AbstractPolynomial} = zero(P, p.var)
"""
    one(::Type{<:AbstractPolynomial})
    one(::AbstractPolynomial)

Returns a representation of 1 as the given polynomial.
"""
Base.one(::Type{P}, var=:x) where {P <: AbstractPolynomial} = ⟒(P)(ones(eltype(P),1), var)  # assumes  p₀ = 1
Base.one(p::P) where {P <: AbstractPolynomial} = one(P, p.var)

Base.oneunit(::Type{P}, args...) where {P <: AbstractPolynomial} = one(P, args...)
Base.oneunit(p::P, args...) where {P <: AbstractPolynomial} = one(p, args...)

"""
    variable(var=:x)
    variable(::Type{<:AbstractPolynomial}, var=:x)
    variable(p::AbstractPolynomial, var=p.var)

Return the monomial `x` in the indicated polynomial basis.  If no type is give, will default to [`Polynomial`](@ref). Equivalent  to  `P(var)`.

# Examples
```jldoctest  common
julia> using Polynomials

julia> x = variable()
Polynomial(x)

julia> p = 100 + 24x - 3x^2
Polynomial(100 + 24*x - 3*x^2)

julia> roots((x - 3) * (x + 2))
2-element Array{Float64,1}:
 -2.0
  3.0

```
"""
variable(::Type{P}, var::SymbolLike = :x) where {P <: AbstractPolynomial} = MethodError()
variable(p::AbstractPolynomial, var::SymbolLike = p.var) = variable(typeof(p), var)
variable(var::SymbolLike = :x) = variable(Polynomial{Int}, var)

# basis
# var is a positional argument, not a keyword; can't deprecate so we do `_var; var=_var`
#@deprecate basis(p::P, k::Int; var=:x)  where {P<:AbstractPolynomial}  basis(p, k, var)
#@deprecate basis(::Type{P}, k::Int; var=:x) where {P <: AbstractPolynomial} basis(P, k,var)
# return the kth basis polynomial for the given polynomial type, e.g. x^k for Polynomial{T}
function basis(::Type{P}, k::Int, _var::SymbolLike=:x; var=_var) where {P <: AbstractPolynomial}
    zs = zeros(Int, k+1)
    zs[end] = 1
    ⟒(P){eltype(P)}(zs, var)
end
basis(p::P, k::Int, _var::SymbolLike=:x; var=_var) where {P<:AbstractPolynomial} = basis(P, k, var)

#=
arithmetic =#
Base.:-(p::P) where {P <: AbstractPolynomial} = P(-coeffs(p), p.var)
Base.:+(c::Number, p::AbstractPolynomial) = +(p, c)
Base.:-(p::AbstractPolynomial, c::Number) = +(p, -c)
Base.:-(c::Number, p::AbstractPolynomial) = +(-p, c)
Base.:*(c::Number, p::AbstractPolynomial) = *(p, c)

function Base.:*(p::P, c::S) where {P <: AbstractPolynomial,S}
    T = promote_type(P, S)
    return T(coeffs(p) .* c, p.var)
end

function Base.:/(p::P, c::S) where {T,P <: AbstractPolynomial{T},S}
    R = promote_type(P, eltype(one(T) / one(S)))
    return R(coeffs(p) ./ c, p.var)
end

Base.:-(p1::AbstractPolynomial, p2::AbstractPolynomial) = +(p1, -p2)

function Base.:+(p::P, n::Number) where {P <: AbstractPolynomial}
    p1, p2 = promote(p, n)
    return p1 + p2
end

function Base.:+(p1::P, p2::O) where {P <: AbstractPolynomial,O <: AbstractPolynomial}
    p1, p2 = promote(p1, p2)
    return p1 + p2
end

function Base.:*(p1::P, p2::O) where {P <: AbstractPolynomial,O <: AbstractPolynomial}
    p1, p2 = promote(p1, p2)
    return p1 * p2
end

Base.:^(p::AbstractPolynomial, n::Integer) = Base.power_by_squaring(p, n)

function Base.divrem(num::P, den::O) where {P <: AbstractPolynomial,O <: AbstractPolynomial}
    n, d = promote(num, den)
    return divrem(n, d)
end

"""
    gcd(a::AbstractPolynomial, b::AbstractPolynomial; atol::Real=0, rtol::Real=Base.rtoldefault)

Find the greatest common denominator of two polynomials recursively using
[Euclid's algorithm](http://en.wikipedia.org/wiki/Polynomial_greatest_common_divisor#Euclid.27s_algorithm).

# Examples

```jldoctest common
julia> using Polynomials

julia> gcd(fromroots([1, 1, 2]), fromroots([1, 2, 3]))
Polynomial(4.0 - 6.0*x + 2.0*x^2)

```
"""
function Base.gcd(p1::AbstractPolynomial{T}, p2::AbstractPolynomial{S}; kwargs...) where {T,S}
    gcd(promote(p1, p2)...; kwargs...)
end

function Base.gcd(p1::P, p2::Q;
                  atol::Real=zero(real(T)),
                  rtol::Real=Base.rtoldefault(real(T))
                  ) where {T, P <: AbstractPolynomial{T}, Q <: AbstractPolynomial{T}}

    r₀, r₁ = p1, p2
    iter = 1
    itermax = length(r₁)

    while !iszero(r₁) && iter ≤ itermax
        _, rtemp = divrem(r₀, r₁)
        r₀ = r₁
        r₁ = truncate(rtemp; atol=atol, rtol=rtol)
        iter += 1
    end
    return r₀
end

"""
    div(::AbstractPolynomial, ::AbstractPolynomial)
"""
Base.div(n::AbstractPolynomial, d::AbstractPolynomial) = divrem(n, d)[1]

"""
    rem(::AbstractPolynomial, ::AbstractPolynomial)
"""
Base.rem(n::AbstractPolynomial, d::AbstractPolynomial) = divrem(n, d)[2]

#=
Comparisons =#
Base.isequal(p1::P, p2::P) where {P <: AbstractPolynomial} = hash(p1) == hash(p2)
Base.:(==)(p1::AbstractPolynomial, p2::AbstractPolynomial) =
    check_same_variable(p1,p2) && (coeffs(p1) == coeffs(p2))
Base.:(==)(p::AbstractPolynomial, n::Number) = degree(p) <= 0 && p[0] == n
Base.:(==)(n::Number, p::AbstractPolynomial) = p == n

function Base.isapprox(p1::AbstractPolynomial{T},
    p2::AbstractPolynomial{S};
    rtol::Real = (Base.rtoldefault(T, S, 0)),
                       atol::Real = 0,) where {T,S}

    p1, p2 = promote(p1, p2)
    check_same_variable(p1, p2)  || error("p1 and p2 must have same var")

    # copy over from abstractarray.jl
    Δ  = norm(p1-p2)
    if isfinite(Δ)
        return Δ <= max(atol, rtol*max(norm(p1), norm(p2)))
    else
        for i in 0:max(degree(p1), degree(p2))
            isapprox(p1[i], p2[i]; rtol=rtol, atol=atol) || return false
        end
        return true
    end
end

function Base.isapprox(p1::P,
                       n::S;
                       rtol::Real = (Base.rtoldefault(T, S, 0)),
                       atol::Real = 0,) where {T,S, P<:AbstractPolynomial{T}}
    return isapprox(p1, ⟒(P){T}(n,p1.var))
end

Base.isapprox(n::S,
    p1::AbstractPolynomial{T};
    rtol::Real = (Base.rtoldefault(T, S, 0)),
    atol::Real = 0,) where {T,S} = isapprox(p1, n, rtol = rtol, atol = atol)


text:
#Written By Johnathan Bizzano
module CommonStaticFunctions

    module Single
        using Base.MathConstants
        using ...Computation
        import ....@Fun

        export Sin, Cos, Tan, ASin, ACos, ATan, Tanh

        Wrap(sig::Signature, f::Function) = MemoryWrapper(
            sig,
            function (sett)
                IT = sett.inputType
                return arg -> IT(f(arg))
            end
        )

        Sin() = Wrap(Sig(:Sin), sin)
        Cos() = Wrap(Sig(:Cos), cos)
        Tan() = Wrap(Sig(:Tan), tan)
        ASin() = Wrap(Sig(:ASin), asin)
        ACos() = Wrap(Sig(:ACos), acos)
        ATan() = Wrap(Sig(:ATan), atan)
        Tanh() = Wrap(Sig(:Tanh), x -> (e ^ x - e ^ -x) / (e ^ x + e ^ -x))

        export RMax, RMin, LMax, LMin,Sigmoid
        RMax(Lower_Bound) = Wrap(Sig(:RMax), x -> max(Lower_Bound, x))
        RMin(Lower_Bound) = Wrap(Sig(:RMin), x -> min(Lower_Bound, x))
        LMax(Upper_Bound) = Wrap(Sig(:LMax), x -> max(x, Upper_Bound))
        LMin(Upper_Bound) = Wrap(Sig(:LMin), x -> min(x, Upper_Bound))
        Sigmoid() = Wrap(Sig(:Sigmoid), x -> 1 / (1 + e ^ -x))

        export Linear, Quadratic, Cubic, Sqrt
        Add(constant) = Wrap(Sig(:Add), x -> x + constant)
        Multiply(constant) = Wrap(Sig(:Multiply), x -> x * constant)
        Linear(slope, constant) = Wrap(Sig(:Linear), x -> slope * x + constant)
        Quadratic(a, b, c) = Wrap(Sig(:Quadratic), x -> a * x ^ 2 + b * x + c)
        Cubic(a, b, c, d) = Wrap(Sig(:Cubic), x -> a * x ^ 3 + b * x ^ 2 + c * x + d)
        Sqrt() = Wrap(Sig(:Sqrt), sqrt)

        export LeftValuePieceWise, RightValuePieceWise
        function LeftValuePieceWise(functions::AbstractArray{Function}, values)
            sort!(values)
            return Wrap(
                Sig(:LeftValuePieceWise),
                function (x)
                    i = 1
                    while i < length(values) && x < values[i]
                        i += 1
                    end
                    functions[i](x)
                end)
        end

        function RightValuePieceWise(functions::AbstractArray{Function}, values)
            sort!(values)
            return Wrap(
                Sig(:RightValuePieceWise),
                    function (x)
                        i = length(slopes)
                        while i > 1 && x > values[i]
                            i -= 1
                        end
                        functions[i](x)
                    end)
        end

    end
end


text:
using OffsetArrays
using Combinatorics
using Base: RefValue

const program = OffsetVector(parse.(Int, split(readchomp("input"), ',')), -1)

struct Machine
    memory::OffsetArray{Int,1,Vector{Int}}
    ip::RefValue{Int}
end

Machine(program::AbstractVector, ip::Integer=0) = Machine(program, Ref(ip))

retrieve(program::AbstractVector, x::Integer, mode::Integer) =
    mode == 0 ? program[x] :
    mode == 1 ? x :
    error("unrecognized mode $mode")

function (m::Machine)(input::Channel{Int}, output::Channel{Int})
    program = m.memory
    pos = m.ip[]
    while pos < length(program) - 1 && program[pos] != 99
        code = digits(program[pos], pad=5)
        opcode, _, mode1, mode2, mode3 = code
        @assert mode3 == 0
        if opcode == 1
            a, b, c = program[pos+1:pos+3]
            program[c] = retrieve(program, a, mode1) + retrieve(program, b, mode2)
            pos += 4
        elseif opcode == 2
            a, b, c = program[pos+1:pos+3]
            program[c] = retrieve(program, a, mode1) * retrieve(program, b, mode2)
            pos += 4
        elseif opcode == 3
            isready(input) || wait(input)
            program[program[pos+1]] = take!(input)
            pos += 2
        elseif opcode == 4
            put!(output, retrieve(program, program[pos+1], mode1))
            pos += 2
        elseif opcode == 5
            if retrieve(program, program[pos+1], mode1) > 0
                pos = retrieve(program, program[pos+2], mode2)
            else
                pos += 3
            end
        elseif opcode == 6
            if retrieve(program, program[pos+1], mode1) == 0
                pos = retrieve(program, program[pos+2], mode2)
            else
                pos += 3
            end
        elseif opcode == 7
            a, b, c = program[pos+1:pos+3]
            program[c] = Int(retrieve(program, a, mode1) < retrieve(program, b, mode2))
            pos += 4
        elseif opcode == 8
            a, b, c = program[pos+1:pos+3]
            program[c] = Int(retrieve(program, a, mode1) == retrieve(program, b, mode2))
            pos += 4
        else
            error("unexpected opcode $opcode as position $pos")
        end
    end
    m.ip[] = pos
    nothing
end

struct Amplifier
    machine::Machine
    phase::Int
    input::Channel{Int}
    output::Channel{Int}
end

function (amp::Amplifier)(input::Integer)
    put!(amp.input, amp.phase)
    put!(amp.input, input)
    amp.machine(amp.input, amp.output)
    take!(amp.output)
end

function part1(program)
    best = maximum(permutations(0:4)) do phases
        channels = [Channel{Int}(Inf) for _ in 1:6]
        mapfoldr((f, x)->f(x), 1:5, init=0) do i
            Amplifier(Machine(copy(program)), phases[i], channels[i], channels[i+1])
        end
    end
    println("Part 1 answer: ", best)
end

function part2(program)
    best = maximum(permutations(5:9)) do phases
        channels = [Channel{Int}(Inf) for _ in 1:5]
        amps = Vector{Amplifier}(undef, 5)
        for i = 1:5
            amps[i] = Amplifier(Machine(copy(program)), phases[i],
                                channels[i], channels[mod1(i+1,5)])
            put!(channels[i], phases[i])
        end
        put!(channels[1], 0)
        tasks = [@async(amps[i].machine(amps[i].input, amps[i].output)) for i = 1:5]
        wait(tasks[1])
        take!(channels[1])
    end
    println("Part 2 answer: ", best)
end

part1(program)
part2(program)

max_stars_repo_name: ararslan/advent-solutions

max_issues_repo_name: ararslan/advent-solutions

max_forks_repo_name: ararslan/advent-solutions


text:
using Test
using LogicCircuits
using LogicCircuits: Element # test some internals

@testset "Trimmed SDD test" begin

    num_vars = 7
    mgr = SddMgr(num_vars, :balanced)

    @test num_variables(mgr) == num_vars
    @test num_nodes(mgr) == 2*num_vars-1
    @test num_edges(mgr) == 2*num_vars-2
    @test mgr isa SddMgr

    @test varsubset(left_most_descendent(mgr), mgr)
    @test varsubset(mgr.left, mgr)
    @test varsubset(mgr.right, mgr)
    @test varsubset_left(mgr.left, mgr)
    @test varsubset_left(mgr.left.left, mgr)
    @test varsubset_left(mgr.left.right, mgr)
    @test varsubset_right(mgr.right, mgr)
    @test varsubset_right(mgr.right.right, mgr)
    @test varsubset_right(mgr.right.left, mgr)

    @test !varsubset(mgr, left_most_descendent(mgr))
    @test !varsubset_left(mgr.right, mgr)
    @test !varsubset_left(mgr.right.left, mgr)
    @test !varsubset_left(mgr.right.right, mgr)
    @test !varsubset_left(mgr, mgr)
    @test !varsubset_left(mgr, mgr.left)
    @test !varsubset_left(mgr, mgr.right)
    @test !varsubset_right(mgr.left, mgr)
    @test !varsubset_right(mgr.left.right, mgr)
    @test !varsubset_right(mgr.left.left, mgr)
    @test !varsubset_right(mgr, mgr)
    @test !varsubset_right(mgr, mgr.left)
    @test !varsubset_right(mgr, mgr.right)

    x = Var(1)
    y = Var(2)

    x_c = compile(mgr, var2lit(x))
    y_c = compile(mgr, var2lit(y))

    @test x_c != y_c

    @test variable(x_c) == x
    @test literal(x_c) == var2lit(x)
    @test vtree(x_c) ∈ mgr
    @test ispositive(x_c)
    @test x_c == compile(mgr, var2lit(x))

    @test variable(y_c) == y
    @test literal(y_c) == var2lit(y)
    @test vtree(y_c) ∈ mgr
    @test ispositive(y_c)
    @test y_c == compile(mgr, var2lit(y))

    notx = -var2lit(x)

    notx_c = compile(mgr,notx)

    @test sat_prob(x_c) == 1//2
    @test sat_prob(notx_c) == 1//2
    @test model_count(x_c,num_vars) == BigInt(2)^(num_vars-1)
    @test model_count(notx_c,num_vars) == BigInt(2)^(num_vars-1)
    @test !notx_c(true)
    @test notx_c(false)

    @test variable(notx_c) == x
    @test literal(notx_c) == notx
    @test vtree(notx_c) ∈ mgr
    @test isnegative(notx_c)
    @test notx_c == compile(mgr, notx)

    true_c = compile(mgr,true)

    @test istrue(true_c)
    @test constant(true_c) == true

    false_c = compile(mgr,false)

    @test isfalse(false_c)
    @test constant(false_c) == false

    @test !true_c == false_c
    @test !false_c == true_c
    @test !x_c == notx_c
    @test !notx_c == x_c

    @test model_count(true_c,num_vars) == BigInt(2)^(num_vars)
    @test model_count(false_c,num_vars) == BigInt(0)

    v1 = compile(mgr, Lit(1))
    v2 = compile(mgr, Lit(2))
    v3 = compile(mgr, Lit(3))
    v4 = compile(mgr, Lit(4))
    v5 = compile(mgr, Lit(5))
    v6 = compile(mgr, Lit(6))
    v7 = compile(mgr, Lit(7))
    @test_throws Exception compile(mgr, Lit(8))

    p1 = [Element(true_c,v3)]
    @test canonicalize(p1, mgr.left.right) === v3
    p2 = [Element(v1,true_c), Element(!v1,false_c)]
    @test canonicalize(p2, mgr.left) === v1

    p3 = [Element(v1,v4), Element(!v1,v7)]
    n1 = canonicalize(p3, mgr)
    p4 = [Element(!v1,v7), Element(v1,v4)]
    n2 = canonicalize(p4,mgr)
    @test n1.vtree.left === mgr.left
    @test n1.vtree.right === mgr.right
    @test n1 === n2
    @test isdeterministic(n1)
    @test n1(true, false, false, true, false, false, false)
    @test !n1(false, true, false, true, false, false, false)

end


text:
using Revise
using NNFEM
using ForwardDiff

E11 = 1.0
E22 = 2.0
G12 = 3.0
C1 = 1.0
C2 = 1.0

f = x -> PlaneStressIncompressibleRivlinSaundersJacobian(x[1], x[2], x[3], C1, C2)[1]
J0 = ForwardDiff.jacobian(f, [E11;E22;G12])

J = PlaneStressIncompressibleRivlinSaundersJacobian(x[1], x[2], x[3], C1, C2)[2]
@show J
@show J0

# using SymPy
# E11, E22, G12, C1, C2 = @vars E11 E22 G12 C1 C2
# E33 = (1.0/(1 + 2*E11 + G12 + 4*E11*E22 - G12*G12) - 1.0) /2.0
# E =[ E11 G12 0
#      G12 E22 0
#      0   0   E33]
# lambda = eigvals(2*E+SymPy.I).^2
# I1 = simplify(sum(lambda))
# I2 = simplify(lambda[1]*lambda[2] + lambda[1]*lambda[3] + lambda[2]*lambda[3])

# W = C1*(I1 - 3) + C2*(I2 - 3)
# σ = Array{Sym}(undef,3)
# σ[1] = simplify(diff(W, E11))
# σ[2] = simplify(diff(W, E22))
# σ[3] = simplify(diff(W, G12))
# J = Matrix{Sym}(undef, 3, 3)
# for i = 1:3
#     J[i,1] = simplify(diff(σ[i], E11))
#     J[i,2] = simplify(diff(σ[i], E22))
#     J[i,3] = simplify(diff(σ[i], G12))
# end
# s = "J = zeros(3,3);\n"
# for i = 1:3
#     for j = 1:3
#         global s *= "J[$i,$j]="*sympy.julia_code(J[i,j])*"\n"
#     end
# end
# t = "σ = zeros(3)\n"
# for i = 1:3
#     global t *= "σ[$i]="*sympy.julia_code(σ[i])*"\n"
# end


text:
@testset "li4" begin
    data = open(readdlm, joinpath(@__DIR__, "data", "Li4.txt"))

    for r in 1:size(data, 1)
        row      = data[r, :]
        z        = row[1] + row[2]*1im
        expected = row[3] + row[4]*1im

        @test PolyLog.li4(z) ≈ expected atol=1e-14
    end
end


text:
## Author:  Pearson
## Date: April 2020
## Julia version: 1.3.1
## Purpose: Solve EKNR(2016) Simple 1 City Case

# Precompile Packages
using Parameters
using MATLAB
using NLsolve
using Plots

# Introduce Functions
include("Simple_1Cty_Functions.jl")

# Program Parameters
T_data = 20
T_tail = 500
T = T_data+T_tail

# Initial Conditions and Time Invariant Parameters
K_init_rel_SS = 1/2
βᴸ = 2/3
βᴷ = 1-βᴸ
δ = 0.1
ρ = 0.95
B = βᴸ^-βᴸ*βᴷ^-βᴷ

# Read data
fname = MatFile(joinpath(@__DIR__, "..", "input", "EKNR_Simple_1Cty_Shocks.mat"))
χ_data = get_variable(fname,"chi_data")
Aᴰ_data = get_variable(fname,"AD_data")
L_data = get_variable(fname,"L_data")

# Hypothetical Data on Shocks(Assumed Constant During T_tail)
χ_tail = χ_data[T_data]*ones(T_tail,1)
Aᴰ_tail = Aᴰ_data[T_data]*ones(T_tail,1)
L_tail = L_data[T_data]*ones(T_tail,1)
χ = [χ_data'; χ_tail]
Aᴰ = [Aᴰ_data'; Aᴰ_tail]
L = [L_data'; L_tail]

# Hypothetical Data on Shocks in Change Form
χ̂ = [χ[2:T]./χ[1:T-1]; 1]
Âᴰ = [Aᴰ[2:T]./Aᴰ[1:T-1]; 1]
L̂ = [L[2:T]./L[1:T-1]; 1]

# Steady State Values (FROM ANALYTICAL SOLUTION)
Kₛₛ = L[T]*βᴷ/βᴸ*(ρ*χ[T]*Aᴰ[T]/(1-ρ*(1-δ)))^(1/βᴸ)
Yₛₛ = (1-ρ*(1-δ))/(1-ρ*(1-δ)-δ*ρ*βᴷ)

###########################################################################################
#                         Solve for Solution to the Problem in Levels
###########################################################################################
# Initial Conditions and Guess
K₁ = Kₛₛ*K_init_rel_SS
guess_lev = ones(eltype(Yₛₛ),T,1).*Yₛₛ

# Wrap All Exogenous Variables and Parameters
myinit_lev = @with_kw (K₁ = K₁,)
myexos_lev = @with_kw (χ = χ, Aᴰ = Aᴰ, L = L, K₁ = K₁, Kₛₛ = Kₛₛ)
myparams_lev = @with_kw (βᴷ = βᴷ, βᴸ = βᴸ, δ = δ, ρ = ρ, B = B, T = T)
init_lev = myinit_lev()
exos_lev = myexos_lev()
params_lev = myparams_lev()

# Solve for the Problem in levels
println("\n\nStart to solve the problem in levels.")
println("Run time and memory cost:")
@time results_lev =
    try
        results_lev = nlsolve(
        (res_lev, guess_lev) -> Fun_1Cty_Levels!(res_lev,
            guess_lev, init_lev, exos_lev, params_lev),
        guess_lev,
        ftol=1e-6,
        autodiff=:forward,
        method=:newton,
        show_trace=false,
    )
    catch err
        if isa(err, DomainError)
            error("Failed to solve the problem in levels, please try again.")
        end
    end

# Check Convergence
converged(results_lev) || error("Failed to converge in $(results_lev.iterations) iterations")
println("Successfully solved problem in levels.\n")

# Catch Solutions
res_lev = similar(results_lev.zero)
res_lev, K_lev, Y_lev = Fun_1Cty_Levels!(
    res_lev, results_lev.zero, init_lev, exos_lev, params_lev)

###########################################################################################
#                         Solve for Solution to the Problem in Changes
###########################################################################################
# Guess
guess_hat = ones(T,1)

# Wrap All Exogenous Variables and Parameters
myinit_hat = @with_kw (Y₁ = Y_lev[1],)
myexos_hat = @with_kw (χ̂ = χ̂, Âᴰ = Âᴰ, L̂ = L̂)
myparams_hat = @with_kw (βᴷ = βᴷ, βᴸ = βᴸ, δ = δ, ρ = ρ, B = B, T = T)
init_hat = myinit_hat()
exos_hat = myexos_hat()
params_hat = myparams_hat()

# Solve for Solution to Problem in Changes
println("Start to solve the problem in changes.")
println("Run time and memory cost:")
@time results_hat =
    try
        results_hat = nlsolve(
            (res_hat, guess_hat) -> Fun_1Cty_Changes!(
                res_hat, guess_hat, init_hat, exos_hat, params_hat),
        guess_hat,
        ftol=1e-6,
        autodiff=:forward,
        method=:newton,
        show_trace=false,
    )
    catch err
        if isa(err, DomainError)
            error("Failed to solve the problem in changes, please try again.")
        end
    end

# Check Convergence
converged(results_hat) || error("Failed to converge in $(results_hat.iterations) iterations")
println("Successfully solved problem in changes.\n")

# Catch Solutions
res_hat = similar(results_hat.zero)
res_hat, K̂, Ŷ = Fun_1Cty_Changes!(
    res_hat, results_hat.zero, init_hat, exos_hat, params_hat)
K_hat = similar(K̂)
Y_hat = similar(Ŷ)
K_hat[1] = K₁
Y_hat[1] = Y_lev[1]
for tt in 2:T
    K_hat[tt] = K_hat[tt-1]*K̂[tt-1]
    Y_hat[tt] = Y_hat[tt-1]*Ŷ[tt-1]
end

###########################################################################################
#                                      Plot Figures
###########################################################################################
println("Plotting the figure.")
l = @layout [a; b]
p1 = Plots.plot(K_lev[1:75], line=(:solid, 2), label="Level Solution");
    Plots.plot!(K_hat[1:75], line=(:dash, 2), label="Change Solution")
p2 = Plots.plot(Y_lev[1:75], line=(:solid, 2), label="Level Solution");
    Plots.plot!(Y_hat[1:75], line=(:dash, 2), label="Change Solution")
figures_all = plot(p1, p2, layout = l)
display(figures_all)
savefig(figures_all, joinpath(@__DIR__, "..", "output", "EKNR_1Cty.pdf"))
println("The figure is saved.")