text:
# crank_nicolson.jl
# Crank-Nicolson implicit solver for PDEs

export crank_nicolson, cnResidual

#TODO: stop doing this
push!(LOAD_PATH, joinpath(Pkg.dir("PumiInterface"), "src"))
push!(LOAD_PATH, joinpath(Pkg.dir("PDESolver"), "src/solver/euler"))
push!(LOAD_PATH, joinpath(Pkg.dir("PDESolver"), "src/NonlinearSolvers"))
push!(LOAD_PATH, joinpath(Pkg.dir("PDESolver"), "src/Debugging"))
push!(LOAD_PATH, joinpath(Pkg.dir("PDESolver"), "src/Utils"))

@doc """
crank_nicolson

  This function performs Crank-Nicolson implicit time solution, using a function
  of the form du/dt = f(u, t)

  Arguments:
    * f  : function evaluation, must have signature (ctx..., opts, t)
    * h  : time step size
    * t_max: time value to stop time stepping (time starts at 0)
    * q_vec: vector of the u values
    * res_vec: vector of du/dt values (the output of the function f)
    * pre_func: function to to be called after the new u values are put into
                q_vec but before the function f is evaluated.  Must have
                signature: pre_func(ctx..., opts)
    * post_func: function called immediately after f is called.  The function
                 must have the signature res_norm = post_func(ctx..., opts,
                 calc_norm=true),
                 where res_norm is a norm of res_vec, and calc_norm determines
                 whether or not to calculate the norm.
    * ctx: a tuple (or any iterable container) of the objects needed by
           f, pre_func, and post func.  The tuple is splatted before being
           passed to the functions.
    * opts : options dictionary

    Keyword Arguments:
    * majorIterationCallback: a callback function called after the first
                              stage, useful to do output or logging
    * res_tol : keyword arg, residual topping tolerance
    * real_time : do actual time marching, not pseudo-time marching

   The eqn.q_vec should hold the whichever variables (conservative or
   entropy) that the simulation should use.

   The idea behind pre_func, post func, and ctx is that they abstract what kind
   of system rk4 is timestepping.  rk4 only needs to know about q_vec and
   res_vec.

   For physics modules, ctx should be (mesh, sbp, eqn) and q_vec and res_vec
   should be eqn.q_vec and eqn.res_vec.

   This function supported jacobian/preconditioner freezing with the prefix
   "CN", with the default setting to never recalculate either.  newtonInner
   will use its recalculation policy to recalculate the PC and jacobian.
"""->
function crank_nicolson(f::Function, h::AbstractFloat, t_max::AbstractFloat,
                        mesh::AbstractMesh, sbp::AbstractOperator, eqn::AbstractSolutionData,
                        opts, res_tol=-1.0, real_time=true)

  myrank = eqn.myrank
  if myrank == 0
    println(BSTDOUT, "\nEntered Crank-Nicolson")
    println(BSTDOUT, "res_tol = ", res_tol)
  end
  flush(BSTDOUT)

  output_freq = opts["output_freq"]::Int
  write_vis = opts["write_vis"]::Bool
  use_itermax = opts["use_itermax"]::Bool
  jac_type = opts["jac_type"]
  if use_itermax
    itermax = opts["itermax"]
  end

  use_checkpointing = opts["use_checkpointing"]::Bool
  chkpoint_freq = opts["checkpoint_freq"]::Int
  ncheckpoints = opts["ncheckpoints"]::Int

#  if jac_type == 4
#    throw(ErrorException("CN not implemented for matrix-free ops. (jac_type cannot be 4)"))
#  end

  @mpi_master if myrank == 0
    _f1 = open("convergence.dat", "a")
    f1 = BufferedIO(_f1)
  end

  t = 0.0
  t_steps = round(Int, t_max/h)

  # eqn_nextstep = deepcopy(eqn)
  eqn_nextstep = eqn_deepcopy(mesh, sbp, eqn, opts)

  # TODO: copyForMultistage! does not give correct values.
  #     deepcopy works for now, but uses more memory than copyForMultistage!, if it worked
  # eqn_nextstep = copyForMultistage!(eqn)
  eqn_nextstep.q = reshape(eqn_nextstep.q_vec, mesh.numDofPerNode,
                           mesh.numNodesPerElement, mesh.numEl)
  eqn_nextstep.res = reshape(eqn_nextstep.res_vec, mesh.numDofPerNode,
                             mesh.numNodesPerElement, mesh.numEl)

  @debug1 println("============ In CN ============")

   # setup all the checkpointing related data
  chkpointer, chkpointdata, skip_checkpoint = explicit_checkpoint_setup(opts, myrank)
  istart = chkpointdata.i

  #-------------------------------------------------------------------------------
  # allocate Jac outside of time-stepping loop

  # NOTE: eqn_nextstep changed to eqn 20161013
  pc, lo = getCNPCandLO(mesh, sbp, eqn, opts)
  ls = StandardLinearSolver(pc, lo, eqn.comm, opts)
  newton_data, rhs_vec = setupNewton(mesh, mesh, sbp, eqn, opts, ls)
  newton_data.itermax = 30
  recalc_policy = getRecalculationPolicy(opts, "CN")

  # this loop is 2:(t_steps+1) when not restarting
  for i = istart:(t_steps + 1)

    t = (i-2)*h
    @mpi_master println(BSTDOUT, "\ni = ", i, ", t = ", t)
    @debug1 println(eqn.params.f, "====== CN: at the top of time-stepping loop, t = $t, i = $i")
    @debug1 flush(eqn.params.f)

    if use_checkpointing && i % chkpoint_freq == 0
      if skip_checkpoint
        skip_checkpoint = false
      else
        @mpi_master println(BSTDOUT, "Saving checkpoint at timestep ", i)
        skip_checkpoint = false
        # save all needed variables to the chkpointdata
        chkpointdata.i = i

        if countFreeCheckpoints(chkpointer) == 0
          freeOldestCheckpoint(chkpointer)  # make room for a new checkpoint
        end

        # save the checkpoint
        saveNextFreeCheckpoint(chkpointer, mesh, sbp, eqn, opts, chkpointdata)
      end   # end of if skip_checkpoint check
    end   # end of if use_checkpointing check

#=
    #----------------------------
    # zero out Jac
    #   this works for both PETSc and Julia matrices.
    #   when jac is a Julia matrix, this is effectively wrapping: fill!(jac, 0.0)
    if jac_type != 4
      MatZeroEntries(jac)
    end
=#

    # TODO: Allow for some kind of stage loop: ES-Dirk

    # TODO: output freq

    # NOTE: Must include a comma in the ctx tuple to indicate tuple
    # f is the physics function, like evalEuler

    t_nextstep = t + h

    # allow for user to select CN's internal Newton's method. Only supports dense FD Jacs, so only for debugging
    if opts["cleansheet_CN_newton"]
      cnNewton(mesh, sbp, opts, h, f, eqn, eqn_nextstep, t_nextstep)
    else

      ctx_residual = (f, eqn, h, newton_data)

      # recalculate PC and LO if needed
      doRecalculation(recalc_policy, i,
                    ls, mesh, sbp, eqn_nextstep, opts, ctx_residual, t_nextstep)

      newtonInner(newton_data, mesh, sbp, eqn_nextstep, opts, cnRhs, ls,
                  rhs_vec, ctx_residual, t_nextstep)
    end

    # do the callback using the current eqn object at time t
    eqn.majorIterationCallback(i, mesh, sbp, eqn, opts, BSTDOUT)

    # need to assemble solution into res_vec?
    res_norm = calcNorm(eqn, eqn.res_vec)
    # logging
    @mpi_master if i % 1 == 0
      println(f1, i, " ", res_norm)
    end

    @mpi_master if i % output_freq == 0
      println(BSTDOUT, "flushing convergence.dat to disk")
      flush(f1)
    end

    # check stopping conditions
    if (res_norm < res_tol)  # ???
      if myrank == 0
        println(BSTDOUT, "breaking due to res_tol, res norm = $res_norm")
        close(f1)
        flush(BSTDOUT)
      end
      break
    end

    if use_itermax && i > itermax
      if myrank == 0
        println(BSTDOUT, "breaking due to itermax")
        close(f1)
        flush(BSTDOUT)
      end
      break
    end

    # This allows the solution to be updated from _nextstep without a deepcopy.
    # There are two memory locations used by eqn & eqn_nextstep,
    #   and this flips the location of eqn & eqn_nextstep every time step
    eqn_temp = eqn
    eqn = eqn_nextstep
    eqn_nextstep = eqn_temp

    # Note: we now need to copy the updated q over for the initial newton guess
    for i = 1:mesh.numDof
      eqn_nextstep.q_vec[i] = eqn.q_vec[i]
    end
    array1DTo3D(mesh, sbp, eqn_nextstep, opts, eqn_nextstep.q_vec, eqn_nextstep.q)

#    t = t_nextstep
    flush(BSTDOUT)

  end   # end of t step loop

  # final time update
  t += h

  # depending on how many timesteps we do, this may or may not be necessary
  #   usage: copyForMultistage!(dest, src)
  copyForMultistage!(eqn, eqn_nextstep)

  #TODO: return the NewtonData?
  free(newton_data)

  @debug1 println("============= end of CN: t = $t ===============")
  return t

  flush(BSTDOUT)

end   # end of crank_nicolson function

"""
  Construct CN preconditioner and linear operator based on options dictionary

  **Inputs**

   * mesh
   * sbp
   * eqn
   * opts
"""
function getCNPCandLO(mesh, sbp, eqn, opts)

  jac_type = opts["jac_type"]

  if jac_type <= 2
    pc = PCNone(mesh, sbp, eqn, opts)
  elseif opts["use_volume_preconditioner"]
    pc = CNBDiagPC(mesh, sbp, eqn, opts)
  else
    pc = CNMatPC(mesh, sbp, eqn, opts)
  end

  if jac_type == 1
    lo = CNDenseLO(pc, mesh, sbp, eqn, opts)
  elseif jac_type == 2
    lo = CNSparseDirectLO(pc, mesh, sbp, eqn, opts)
  elseif jac_type == 3
    lo = CNPetscMatLO(pc, mesh, sbp, eqn, opts)
  elseif jac_type == 4
    lo = CNPetscMatFreeLO(pc, mesh, sbp, eqn, opts)
  end

  return pc, lo
end

#------------------------------------------------------------------------------
# Regular (matrix-based) PC

mutable struct CNMatPC <: AbstractPetscMatPC
  pc_inner::NewtonMatPC
end

function CNMatPC(mesh::AbstractMesh, sbp::AbstractOperator,
                    eqn::AbstractSolutionData, opts::Dict)

  pc_inner = NewtonMatPC(mesh, sbp, eqn, opts)

  return CNMatPC(pc_inner)
end

function calcPC(pc::CNMatPC, mesh::AbstractMesh, sbp::AbstractOperator,
                eqn::AbstractSolutionData, opts::Dict, ctx_residual, t)

  calcPC(pc.pc_inner, mesh, sbp, eqn, opts, ctx_residual, t)
  modifyJacCN(pc, mesh, sbp, eqn, opts, ctx_residual, t)

  return nothing
end

#------------------------------------------------------------------------------
# Volume preconditioner

"""
  Volume integral jacobian preconditioner for CN.

"""
mutable struct CNBDiagPC <: AbstractPetscMatFreePC
  pc_inner::PetscMatFreePC
  other_pc::NewtonBDiagPC
end  # end type definition

function CNBDiagPC(mesh::AbstractMesh, sbp,
                              eqn::AbstractSolutionData, opts)

  pc_inner = PetscMatFreePC(mesh, sbp, eqn, opts)
  other_pc = NewtonBDiagPC(mesh, sbp, eqn, opts)

  return CNBDiagPC(pc_inner, other_pc)
end

function calcPC(pc::CNBDiagPC, mesh::AbstractMesh, sbp::AbstractOperator,
                eqn::AbstractSolutionData, opts::Dict, ctx_residual, t)

# eqn should be eqn_nextstep from the CN function

  cnCalcBDiagPC(pc, mesh, sbp, eqn, opts, ctx_residual, t)
  factorBDiagPC(pc, mesh, sbp, eqn, opts)

  # this needs to come last (because calcBDiagPC calls it too)
  setPCCtx(pc, mesh, sbp, eqn, opts, ctx_residual, t)

  return nothing
end

function applyPC(pc::CNBDiagPC, mesh::AbstractMesh, sbp::AbstractOperator,
                 eqn::AbstractSolutionData, opts::Dict, t, b::AbstractVector,
                 x::AbstractVector)

  applyBDiagPC(pc.other_pc, mesh, sbp, eqn, opts, b, x)

  return nothing
end

#------------------------------------------------------------------------------
# Linear operators

#TODO: make a macro to generate these definitions

mutable struct CNDenseLO <: AbstractDenseLO
  lo_inner::NewtonDenseLO
end

function CNDenseLO(pc::PCNone, mesh::AbstractMesh,
                    sbp::AbstractOperator, eqn::AbstractSolutionData, opts::Dict)

  lo_inner = NewtonDenseLO(pc, mesh, sbp, eqn, opts)

  return CNDenseLO(lo_inner)
end

mutable struct CNSparseDirectLO <: AbstractSparseDirectLO
  lo_inner::NewtonSparseDirectLO
end

function CNSparseDirectLO(pc::PCNone, mesh::AbstractMesh,
                    sbp::AbstractOperator, eqn::AbstractSolutionData, opts::Dict)

  lo_inner = NewtonSparseDirectLO(pc, mesh, sbp, eqn, opts)

  return CNSparseDirectLO(lo_inner)
end

mutable struct CNPetscMatLO <: AbstractPetscMatLO
  lo_inner::NewtonPetscMatLO
end

function CNPetscMatLO(pc::AbstractPetscPC, mesh::AbstractMesh,
                    sbp::AbstractOperator, eqn::AbstractSolutionData, opts::Dict)

  lo_inner = NewtonPetscMatLO(pc, mesh, sbp, eqn, opts)

  return CNPetscMatLO(lo_inner)
end

"""
  All CN LOs that have matrices
"""
const CNMatLO = Union{CNDenseLO, CNSparseDirectLO, CNPetscMatLO}

function calcLinearOperator(lo::CNMatLO, mesh::AbstractMesh,
                            sbp::AbstractOperator, eqn::AbstractSolutionData,
                            opts::Dict, ctx_residual, t)

  calcLinearOperator(lo.lo_inner, mesh, sbp, eqn, opts, ctx_residual, t)

  modifyJacCN(lo, mesh, sbp, eqn, opts, ctx_residual, t)

  return nothing
end

# matrix-free
mutable struct CNPetscMatFreeLO <: AbstractPetscMatFreeLO
  lo_inner::NewtonPetscMatFreeLO
end

function CNPetscMatFreeLO(pc::AbstractPetscPC, mesh::AbstractMesh,
                      sbp::AbstractOperator, eqn::AbstractSolutionData, opts::Dict)

  lo_inner = NewtonPetscMatFreeLO(pc, mesh, sbp, eqn, opts)

  return CNPetscMatFreeLO(lo_inner, rhs_func)
end

"""
  Any PC or LO that has a matrix in the field `A`
"""
const CNHasMat = Union{CNMatPC, CNDenseLO, CNSparseDirectLO, CNPetscMatLO}

function calcLinearOperator(lo::CNPetscMatFreeLO, mesh::AbstractMesh,
                            sbp::AbstractOperator, eqn::AbstractSolutionData,
                            opts::Dict, ctx_residual, t)

  calcLinearOperator(lo.lo_inner, mesh, sbp, eqn, opts, ctx_residual, t)

  setLOCtx(lo, mesh, sbp, eqn, opts, ctx_residual, t)

  return nothing
end

function applyLinearOperator(lo::CNPetscMatFreeLO, mesh::AbstractMesh,
                       sbp::AbstractOperator, eqn::AbstractSolutionData{Tsol},
                       opts::Dict, ctx_residual, t, x::AbstractVector,
                       b::AbstractVector) where Tsol

  h = ctx[3]
  # the CN residual has the form: I - 0.5*delta_t*jac, where jac is the physics
  # Jacobian.
  # thus a matrix vector product is: v - 0.5*delta_t*jac*v
  applyLinearOperator(lo, mesh, sbp, eqn, opts, ctx_residual, t, x, b)

  scale!(b, -0.5*h)
  @simd for i=1:length(b)
    b[i] += x[i]
  end

  return nothing
end

function applyLinearOperatorTranspose(lo::CNPetscMatFreeLO,
                             mesh::AbstractMesh,
                             sbp::AbstractOperator, eqn::AbstractSolutionData{Tsol},
                             opts::Dict, ctx_residual, t, x::AbstractVector,
                             b::AbstractVector) where Tsol

  # see the non-transpose method for an explanation of the math

  h = ctx[3]
  applyLinearOperatorTranspose(lo, mesh, sbp, eqn, opts, ctx_residual, t, x, b)

  scale!(b, -0.5*h)
  @simd for i=1:length(b)
    b[i] += x[i]
  end

  return nothing
end

"""
  Takes the Jacobian of the physics and modifies it to be the Crank-Nicolson
  Jacobian.

  CN_Jac = I - dt/2 * physicsJac

  **Inputs**

   * lo: a CNHasMat
   * mesh
   * sbp
   * eqn
   * opts
   * ctx_residual: delta_t must be the 3rd element
   * t
"""
function modifyJacCN(lo::CNHasMat, mesh, sbp, eqn, opts, ctx_residual, t)

  lo2 = getBaseLO(lo)
  h = ctx_residual[3]

  assembly_begin(lo2.A, MAT_FINAL_ASSEMBLY)
  assembly_end(lo2.A, MAT_FINAL_ASSEMBLY)

  # scale jac by -delta_t/2
#  scale_factor = h*-0.5
  petsc_scale_factor = PetscScalar(-h*0.5)
  scale!(lo2.A, petsc_scale_factor)

  # add the identity
  diagonal_shift!(lo2.A, 1)

  return nothing
end

"""
  This function uses [`calcBDiagPC`](@ref) to compute a
  a preconditioner for the matrix used by Crank-Nicolson.  This function
  has the same signature as [`cnJac`](@ref) to facilitate matrix-free
  operations.
"""
function cnCalcBDiagPC(pc::CNBDiagPC, mesh::AbstractDGMesh,
               sbp::AbstractOperator, eqn_nextstep::AbstractSolutionData, opts,
               ctx_residual, t)
  #TODO: is this ctx_residual?
  physics_func = ctx[1]
#  eqn = ctx[2]
  h = ctx[3]

  epsilon = opts["epsilon"]
  pert = Complex128(0, epsilon)

  other_pc = pc.other_pc

#  NonlinearSolvers.physicsJac(newton_data, mesh, sbp, eqn_nextstep, opts, jac, ctx, t_nextstep)
  calcBDiagPC(other_pc, mesh, sbp, eqn_nextstep, opts, pert,
                           physics_func, t)

  volume_prec = other_pc.vol_prec
  scale!(other_pc.volume_jac, -0.5*h)

  # add the identity matrix
  jac_size = mesh.numDofPerNode*mesh.numNodesPerElement

  for i=1:mesh.numEl
    for j=1:jac_size
      other_pc.volume_jac[j, j, i] += 1
    end
  end

  return nothing
end

@doc """
###NonlinearSolvers.cnRhs

  RHS of the CN calculation

  ctx:
    physics func must be the first element, i.e. evalEuler
    eqn_nextstep must be the second element
    h must be the third element

"""->
function cnRhs(mesh::AbstractMesh, sbp::AbstractOperator, eqn_nextstep::AbstractSolutionData, opts, rhs_vec, ctx, t)

  # eqn comes in through ctx_residual, which is set up in CN before the newtonInner call

  physics_func = ctx[1]
  eqn = ctx[2]
  h = ctx[3]

  # evalute residual at t_nextstep
  # q_vec -> q
  array1DTo3D(mesh, sbp, eqn_nextstep, opts, eqn_nextstep.q_vec)

  # start parallel communication if needed
  time = eqn.params.time
  time.t_send += @elapsed if opts["parallel_type"] == 2 && mesh.npeers > 0
    startSolutionExchange(mesh, sbp, eqn_nextstep, opts)
  end

  physics_func(mesh, sbp, eqn_nextstep, opts, t)
  array3DTo1D(mesh, sbp, eqn_nextstep, opts, eqn_nextstep.res,
                   eqn_nextstep.res_vec)

  # evalute residual at t - h
  if opts["parallel_type"] == 2 && mesh.npeers > 0
    startSolutionExchange(mesh, sbp, eqn, opts)
  end
  physics_func(mesh, sbp, eqn, opts, t - h)
  array3DTo1D(mesh, sbp, eqn, opts, eqn.res, eqn.res_vec)

  # compute rhs

  #   what this is doing:
  #   u_(n+1) - 0.5*dt* (del dot G_(n+1)) - u_n - 0.5*dt* (del dot G_n)
  for i = 1:mesh.numDof

    temp1 = eqn_nextstep.q_vec[i] - 0.5*h*eqn_nextstep.res_vec[i]
    temp2 = eqn.q_vec[i] + 0.5*h*eqn.res_vec[i]

    rhs_vec[i] = temp1 - temp2

    # NOTE: question: is there a sign problem here? should rhs_vec = -rhs_vec ?
    #     NO. this negative gets applied in newton.jl, where res_0[i] = -res_0[i]

  end

  # calculate the norm
  rhs_0_norm = calcNorm(eqn, rhs_vec, strongres=true)

  return rhs_0_norm

end     # end of function cnRhs

# the goal is to replace newton.jl.
# this will go into CN in the time-stepping loop
function cnNewton(mesh, sbp, opts, h, physics_func, eqn, eqn_nextstep, t)
  println("++++++++++++++++ clean sheet Newton being run ++++++++++")

  println("---- physics_func: ",physics_func)

  # Jac on eqn or eqn_nextstep?

  epsilon = 1e-8
  t_prev = t - h

  jac = zeros(mesh.numDof, mesh.numDof)

  # emulates physicsJac
  # so we need to step through the jacobian column wise.
  #   d res[1]/d q[1]   d res[1]/d q[2]   d res[1]/d q[3] ...
  #   d res[2]/d q[1]   d res[2]/d q[2]   d res[2]/d q[3] ...
  #   d res[3]/d q[1]   d res[3]/d q[2]   d res[3]/d q[3] ...
  #   ...               ...               ...

  newton_itermax = 2
  delta_q_vec = zeros(eqn_nextstep.q_vec)

  # newton_loop starting here?
  for newton_i = 1:newton_itermax

    #--------------------------
    # emulates physicsJac
    unperturbed_q_vec = copy(eqn_nextstep.q_vec)

    physics_func(mesh, sbp, eqn_nextstep, opts, t)
    # needed b/c physics_func only updates eqn.res
    array3DTo1D(mesh, sbp, eqn_nextstep, opts, eqn_nextstep.res, eqn_nextstep.res_vec)
    # Comment here about mass matrix inv multiplication TODO
    applyMassMatrixInv(mesh, eqn_nextstep, eqn_nextstep.res_vec)
    unperturbed_res_vec = copy(eqn_nextstep.res_vec)

    for i = 1:mesh.numDof
      eqn_nextstep.q_vec[i] = eqn_nextstep.q_vec[i] + epsilon

      physics_func(mesh, sbp, eqn_nextstep, opts, t)
      array3DTo1D(mesh, sbp, eqn_nextstep, opts, eqn_nextstep.res, eqn_nextstep.res_vec)
      applyMassMatrixInv(mesh, eqn, eqn_nextstep.res_vec)

      jac[:,i] = (eqn_nextstep.res_vec - unperturbed_res_vec)/epsilon

      eqn_nextstep.q_vec[i] = unperturbed_q_vec[i]

    end

    #--------------------------
    # emulates cnJac
    scale!(jac, -0.5*h)
    for i = 1:mesh.numDof
      jac[i,i] += 1
    end

    #--------------------------
    # emulates cnRhs
    #   what this is doing:
    #   u_(n+1) - 0.5*dt* (del dot G_(n+1)) - u_n - 0.5*dt* (del dot G_n)
    #=
    physics_func(mesh, sbp, eqn, opts, t)
    array3DTo1D(mesh, sbp, eqn, opts, eqn.res, eqn.res_vec)
    physics_func(mesh, sbp, eqn_nextstep, opts, t_nextstep)
    array3DTo1D(mesh, sbp, eqn_nextstep, opts, eqn_nextstep.res, eqn_nextstep.res_vec)

    rhs_vec = zeros(eqn.q_vec)

    for i = 1:mesh.numDof
      rhs_vec[i] = eqn_nextstep.q_vec[i] - h*0.5*eqn_nextstep.res_vec[i] - eqn.q_vec[i] - h*0.5*eqn.res_vec[i]
    end
    =#
    physics_func(mesh, sbp, eqn, opts, t_prev)
    array3DTo1D(mesh, sbp, eqn, opts, eqn.res, eqn.res_vec)
    applyMassMatrixInv(mesh, eqn, eqn.res_vec)
    current_t_step_contribution = zeros(eqn.q_vec)
    for i = 1:mesh.numDof
      current_t_step_contribution[i] = - eqn.q_vec[i] - h*0.5*eqn.res_vec[i]
    end

    # Test for 3D Minv results
    # this works!
#     res_vec_control = deepcopy(eqn.res_vec)
#     res_vec_test = deepcopy(eqn.res_vec)
#     res_control = deepcopy(eqn.res)
#     res_test = deepcopy(eqn.res)
#     applyMassMatrixInv3D(mesh, sbp, eqn, res_test)
#     array3DTo1D(mesh, sbp, eqn, opts, res_test, res_vec_test)
#     println("=+=+=+ norm of diff btwn res_vec_test & res_vec_control: ", norm(res_vec_test - res_vec_control))

    physics_func(mesh, sbp, eqn_nextstep, opts, t)
    array3DTo1D(mesh, sbp, eqn_nextstep, opts, eqn_nextstep.res, eqn_nextstep.res_vec)
    applyMassMatrixInv(mesh, eqn_nextstep, eqn_nextstep.res_vec)
    next_t_step_contribution = zeros(eqn.q_vec)
    for i = 1:mesh.numDof
      next_t_step_contribution[i] = eqn_nextstep.q_vec[i] - h*0.5*eqn_nextstep.res_vec[i]
    end

    rhs_vec = zeros(eqn.q_vec)

    for i = 1:mesh.numDof
      rhs_vec[i] = current_t_step_contribution[i] + next_t_step_contribution[i]
    end

    # TODO: check these args
    rhs_norm = calcNorm(eqn, rhs_vec, strongres=true)

    #--------------------------
    # start of actual Newton
    neg_rhs = scale(rhs_vec, -1.0)

    fill!(delta_q_vec, 0.0)
    delta_q_vec = jac\neg_rhs
    fill!(jac, 0.0)

    for i = 1:mesh.numDof
      eqn_nextstep.q_vec[i] += delta_q_vec[i]
    end

    rhs_norm_tol = 1e-6
    if rhs_norm < rhs_norm_tol
      println("=== cnNewton converged with rhs_norm under $rhs_norm_tol -- newton iters: $newton_i ===")
      return nothing
    end

  end   # end of newton iterations

  println("=== cnNewton did not converge ===")
  return nothing

end

#TODO: why not use the versions in Utils?
# TODO: comment here
function applyMassMatrixInv(mesh, eqn, vec)

  for k = 1:mesh.numDof
    vec[k] = eqn.Minv[k] * vec[k]
  end

  return vec
end

# TODO: comment here
function applyMassMatrixInv3D(mesh, sbp, eqn, arr)

  for i = 1:mesh.numEl
    for j = 1:sbp.numnodes
      for k = 1:mesh.numDofPerNode
        arr[k, j, i] = eqn.Minv3D[k, j, i] * arr[k, j, i]
      end
    end
  end

  return arr
end


text:
#=
Define a transition model for the pedestrian MDP problem
=#

function POMDPs.transition(mdp::PedMDP, s::PedMDPState, a::PedMDPAction)
    # ego transition
    acc = LonAccelDirection(a.acc, 4) #TODO paremeterize
    ego_p = propagate(s.ego, acc, mdp.env.roadway, mdp.ΔT)
    # interpolate ego_p in discrete space
    ego_v_space = get_car_vspace(mdp.env, mdp.vel_res)
    ego_ps, ego_probs = interpolate_state(mdp, ego_p, ego_v_space)
    ped_probs = Float64[]
    # pedestrian transition
    if s.ped.posG == mdp.off_grid
        ped_ps = pedestrian_starting_states(mdp)
        push!(ped_ps, get_off_the_grid(mdp))
        ped_probs = ones(length(ped_ps))
        ped_probs[end] = 1-mdp.ped_birth
        ped_probs[1:end-1] = mdp.ped_birth/(length(ped_ps)-1)
        @assert sum(ped_probs) ≈ 1.0
    else
        ped_actions, ped_actions_probs = get_distribution(mdp.ped_model)
        ped_ps = VehicleState[]
        for (i, ped_a) in enumerate(ped_actions)
            p_a = ped_actions_probs[i]
            ped_p = propagate(s.ped, ped_a, mdp.env.roadway, mdp.ΔT)
            if (ped_p.posF.s >= get_end(get_lane(mdp.env.roadway, ped_p)) && isapprox(ped_p.posF.ϕ, 0.)) ||
               (isapprox(ped_p.posF.s, 0., atol=0.01) && isapprox(ped_p.posF.ϕ, float(pi)))
                ped_pp = get_off_the_grid(mdp)
                index_ped_pp = findall(x -> x==ped_pp, ped_ps)
                if isempty(index_ped_pp)
                    push!(ped_ps, get_off_the_grid(mdp))
                    push!(ped_probs, p_a)
                else
                    ped_probs[index_ped_pp] += p_a
                end
                continue
            end
            # interpolate ped_p in discrete space
            ped_v_space = get_ped_vspace(mdp.env, mdp.vel_ped_res)
            itp_ped_ps, itp_ped_weights = interpolate_pedestrian(mdp, ped_p, ped_v_space)
            for (j, ped_pss) in enumerate(itp_ped_ps)
                index_itp_state = findall(x -> x==ped_pss, ped_ps)
                w = itp_ped_weights[j]*p_a
                # if !(w ≈  0.)
                    if isempty(index_itp_state)
                        push!(ped_ps, ped_pss)
                        push!(ped_probs, w)
                    else
                        ped_probs[index_itp_state] += w
                    end
                # end
            end
        end
        @assert length(ped_probs) == length(ped_ps)
        normalize!(ped_probs, 1)
    end

    # combine ped and ego states
    states_p = Vector{PedMDPState}(length(ego_ps)*length(ped_ps)) # future states
    probs_p = Vector{Float64}(length(ego_ps)*length(ped_ps)) # proba of future states
    k = 1
    for (i, e) in enumerate(ego_ps)
        for (j, p) in enumerate(ped_ps)
            crash =  is_colliding(Vehicle(e, mdp.ego_type, EGO_ID), Vehicle(p, mdp.ped_type, PED_ID))
            states_p[k] = PedMDPState(crash, e, p)
            probs_p[k] = ego_probs[i]*ped_probs[j]
            k += 1
        end
    end
    normalize!(probs_p, 1)
    # remove duplicates
    # sort by si
    p = sortperm(states_p, by=s->stateindex(mdp, s))
    return SparseCat(states_p[p], probs_p[p])
end

ext: jl

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule OpenMPI_jll
using Base
using Base: UUID
using LazyArtifacts
import JLLWrappers

JLLWrappers.@generate_main_file_header("OpenMPI")
JLLWrappers.@generate_main_file("OpenMPI", UUID("fe0851c0-eecd-5654-98d4-656369965a5c"))
end  # module OpenMPI_jll


text:
"""
    randomwalk(g, s, niter)

Perform a random walk on graph `g` starting at vertex `s` and continuing for
a maximum of `niter` steps. Return a vector of vertices visited in order.
"""
function randomwalk(g::AGraphOrDiGraph, s::Integer, niter::Integer)
    s in vertices(g) || throw(BoundsError())
    T = vertextype(g)
    visited = Vector{T}()
    sizehint!(visited, niter)
    currs = s
    i = 1
    while i <= niter
        push!(visited, currs)
        i += 1
        nbrs = out_neighbors(g, currs)
        length(nbrs) == 0 && break
        currs = myrand(nbrs)
    end
    return visited[1:(i - 1)]
end

"""
    nonbacktracking_randomwalk(g, s, niter)

Perform a non-backtracking random walk on directed graph `g` starting at
vertex `s` and continuing for a maximum of `niter` steps. Return a
vector of vertices visited in order.
"""
function nonbacktracking_randomwalk(g::AGraph, s::Integer, niter::Integer)
    s in vertices(g) || throw(BoundsError())
    T = vertextype(g)
    visited = Vector{T}()
    sizehint!(visited, niter)
    currs = s
    prev = -1
    i = 1

    push!(visited, currs)
    i += 1
    nbrs = out_neighbors(g, currs)
    length(nbrs) == 0 && return visited[1:(i - 1)]
    prev = currs
    currs = myrand(nbrs)

    while i <= niter
        push!(visited, currs)
        i += 1
        # Use collect to work with generic iterable
        # becouse we need indexing (alternative is IterTools.nth)
        # e.g. Network returns generators
        # TODO: do something more efficient
        nbrs = out_neighbors(g, currs) |> collect
        length(nbrs) == 1 && break
        idnext = rand(1:(length(nbrs) - 1))
        next = nbrs[idnext]
        if next == prev
            next = nbrs[end]
        end
        prev = currs
        currs = next
    end
    return visited[1:(i - 1)]
end

function nonbacktracking_randomwalk(g::ADiGraph, s::Integer, niter::Integer)
    s in vertices(g) || throw(BoundsError())
    T = vertextype(g)
    visited = Vector{T}()
    sizehint!(visited, niter)
    currs = s
    prev = -1
    i = 1

    while i <= niter
        push!(visited, currs)
        i += 1
        nbrs = out_neighbors(g, currs) |> collect
        length(nbrs) == 0 && break
        next = myrand(nbrs)
        if next == prev
            length(nbrs) == 1 && break
            idnext = rand(1:(length(nbrs) - 1))
            next = nbrs[idnext]
            if next == prev
                next = nbrs[end]
            end
        end
        prev = currs
        currs = next
    end
    return visited[1:(i - 1)]
end

"""
    self_avoiding_randomwalk(g, s, niter)

Perform a [self-avoiding walk](https://en.wikipedia.org/wiki/Self-avoiding_walk)
on graph `g` starting at vertex `s` and continuing for a maximum of `niter` steps.
Return a vector of vertices visited in order.
"""
function self_avoiding_randomwalk(g::AGraphOrDiGraph, s::Integer, niter::Integer)
    s in vertices(g) || throw(BoundsError())
    T = vertextype(g)
    visited = Vector{T}()
    svisited = Set{T}()
    sizehint!(visited, niter)
    sizehint!(svisited, niter)
    currs = s
    i = 1
    while i <= niter
        push!(visited, currs)
        push!(svisited, currs)
        i += 1
        nbrs = setdiff(Set(out_neighbors(g, currs)), svisited)
        length(nbrs) == 0 && break
        currs = myrand(nbrs)
    end
    return visited[1:(i - 1)]
end


text:
isdefined(Main, :__precompile__) && __precompile__()

module KnowedMimeTypes

export get_mime, MIMETYPES, MimeTypes

type MimeTypes
	data

	MimeTypes(data) = new(data)
	MimeTypes() = new(Void)
end

include("mimetypes.jl")

const MIMETYPES = MimeTypes(MIMETYPES_DICT)

Base.get(mt::MimeTypes, ext::AbstractString, default=Void) = Base.get(mt.data, ext, default)
get_mime(mt::MimeTypes, ext::AbstractString, default=Void) = Base.get(mt, ext, default)
get_mime(ext::AbstractString, default=Void) = Base.get(MIMETYPES, ext, default)

end # module


text:
#=
    Used from https://github.com/lanl-ansi/Juniper.jl
=#

#=
    Used from https://github.com/lanl-ansi/Alpine.jl
=#
function expr_dereferencing!(expr, m)
    for i in 2:length(expr.args)
        if isa(expr.args[i], Union{Float64,Int64,JuMP.VariableRef})
            k = 0
        elseif expr.args[i].head == :ref
            @assert isa(expr.args[i].args[2], MOI.VariableIndex)
            expr.args[i] = JuMP.VariableRef(m, expr.args[i].args[2])
        elseif expr.args[i].head == :call
            expr_dereferencing!(expr.args[i], m)
        else
            error("expr_dereferencing :: Unexpected term in expression tree.")
        end
    end
end

function expr_dereferencing(expr, m)
    c_expr = copy(expr)
    expr_dereferencing!(c_expr, m)
    return c_expr
end

function get_type_dict(obj)
    T = typeof(obj)
    type_dict = Dict{Symbol,Type}()
    for (name, typ) in zip(fieldnames(T), T.types)
        type_dict[name] = typ
    end
    return type_dict
end

"""
    is_global_status(state::MOI.TerminationStatusCode)
Returns true if either ALMOST_OPTIMAL, OPTIMAL or INFEASIBLE and false otherwise
"""
function is_global_status(state::MOI.TerminationStatusCode)
    return state == MOI.ALMOST_OPTIMAL || state == MOI.OPTIMAL || state == MOI.INFEASIBLE
end

"""
    only_almost_solved(state::MOI.TerminationStatusCode)
Returns true if either ALMOST_OPTIMAL or ALMOST_LOCALLY_SOLVED
"""
function only_almost_solved(state::MOI.TerminationStatusCode)
    return state == MOI.ALMOST_OPTIMAL || state == MOI.ALMOST_LOCALLY_SOLVED
end

"""
    state_is_optimal(state::MOI.TerminationStatusCode; allow_almost=false)
Returns true if either optimal or locally solved. If allow_almost then check for `ALMOST_LOCALLY_SOLVED` and `ALMOST_OPTIMAL`
"""
function state_is_optimal(state::MOI.TerminationStatusCode; allow_almost=false)
    return state == MOI.OPTIMAL || state == MOI.LOCALLY_SOLVED ||
            (allow_almost && state == MOI.ALMOST_LOCALLY_SOLVED) || (allow_almost && state == MOI.ALMOST_OPTIMAL)
end

"""
    state_is_infeasible(state::MOI.TerminationStatusCode)
Returns true if either infeasible or locally infeasible
"""
function state_is_infeasible(state::MOI.TerminationStatusCode)
    return state == MOI.INFEASIBLE || state == MOI.LOCALLY_INFEASIBLE
end

"""
    add_obj_constraint(jp::HybridProblem, rhs::Float64)
Add a constraint for the objective based on whether the objective is linear/quadratic or non linear.
If the objective sense is :MIN than add objective <= rhs else objective >= rhs
"""
function add_obj_constraint(jp::HybridProblem, rhs::Float64)
    if jp.has_nl_objective
        obj_expr = MOI.objective_expr(jp.nlp_evaluator)
        if jp.obj_sense == :Min
            obj_constr = Expr(:call, :<=, obj_expr, rhs)
        else
            obj_constr = Expr(:call, :>=, obj_expr, rhs)
        end
        HybridMINLPSolver.expr_dereferencing!(obj_constr, jp.model)
        JuMP.add_NL_constraint(jp.model, obj_constr)
    else # linear or quadratic
        backend = JuMP.backend(jp.model);
        if isa(jp.objective, MOI.SingleVariable)
            if jp.obj_sense == :Min
                JuMP.set_upper_bound(jp.x[jp.objective.variable.value], rhs)
            else
                JuMP.set_lower_bound(jp.x[jp.objective.variable.value], rhs)
            end
        else
            if jp.obj_sense == :Min
                MOI.add_constraint(backend, jp.objective, MOI.LessThan(rhs))
            else
                MOI.add_constraint(backend, jp.objective, MOI.GreaterThan(rhs))
            end
        end
    end
end

"""
    evaluate_objective(optimizer::MOI.AbstractOptimizer, jp::HybridProblem, xs::Vector{Float64})
If no objective exists => return 0
Evaluate the objective whether it is non linear, linear or quadratic
"""
function evaluate_objective(optimizer::MOI.AbstractOptimizer, jp::HybridProblem, xs::Vector{Float64})
    if optimizer.nlp_data.has_objective
        return MOI.eval_objective(optimizer.nlp_data.evaluator, xs)
    elseif optimizer.objective !== nothing
        return MOIU.eval_variables(vi -> xs[vi.value], optimizer.objective)
    else
        return 0.0
    end
end

function set_time_limit!(optimizer, time_limit::Union{Nothing,Float64})
    old_time_limit = Inf
    if MOI.supports(optimizer, MOI.TimeLimitSec())
        old_time_limit = MOI.get(optimizer, MOI.TimeLimitSec())
        MOI.set(optimizer, MOI.TimeLimitSec(), time_limit)
    end
    return old_time_limit
end

"""
    optimize_get_status_backend(model::JuMP.Model; solver=nothing)
Run optimize! and get the status and the backend
"""
function optimize_get_status_backend(model::JuMP.Model; solver=nothing)
    if solver === nothing
        JuMP.optimize!(model)
    else
        JuMP.set_optimizer(model, solver)
        JuMP.optimize!(model)
    end
    backend = JuMP.backend(model)
    status = MOI.get(backend, MOI.TerminationStatus())
    return status, backend
end

function register_functions!(model, registered_functions)
    if registered_functions !== nothing
        for reg_f in registered_functions
            if reg_f.gradf === nothing
                JuMP.register(model, reg_f.s, reg_f.dimension, reg_f.f; autodiff=reg_f.autodiff)
            elseif reg_f.grad2f === nothing
                JuMP.register(model, reg_f.s, reg_f.dimension, reg_f.f, reg_f.gradf)
            else
                JuMP.register(model, reg_f.s, reg_f.dimension, reg_f.f, reg_f.gradf, reg_f.grad2f)
            end
        end
    end
end

max_stars_repo_name: MdMahbubarRahman/Duran

max_issues_repo_name: MdMahbubarRahman/Duran

max_forks_repo_name: MdMahbubarRahman/Duran


text:
module Test_KinematicRevoluteJoints

using  Modia3D
using  Modia3D.StaticArrays
import Modia3D.ModiaMath

vgreen  = Modia3D.Material(color="Green", transparency=0.5)
vyellow = Modia3D.Material(color="Yellow", transparency=0.5)
vred    = Modia3D.Material(color="Red", transparency=0.5)
cmat    = "Steel"

Dx = 0.2
Lx = 1.0
Ly = 0.3
Lz = 0.3

# a few signal macros
@signal SinWave(;A=5.0, B=2.0) begin
   y1 = ModiaMath.RealScalar("y", causality=ModiaMath.Output, numericType=ModiaMath.WR)
end
function Modia3D.computeSignal(signal::SinWave, sim::ModiaMath.SimulationState)
    signal.y1.value = signal.A*sin(signal.B*sim.time)
end

@signal CosWave(;C=5.0) begin
   y = ModiaMath.RealScalar("y", causality=ModiaMath.Output, numericType=ModiaMath.WR)
end
function Modia3D.computeSignal(signal::CosWave, sim::ModiaMath.SimulationState)
    signal.y.value = signal.C*cos(sim.time)
end

@signal LinearSig(;D=5.0, E=2.0) begin
   y = ModiaMath.RealScalar("y", causality=ModiaMath.Output, numericType=ModiaMath.WR)
end
function Modia3D.computeSignal(signal::LinearSig, sim::ModiaMath.SimulationState)
    signal.y.value = signal.D*sim.time + signal.E
end

@assembly Pendulum() begin
world    = Modia3D.Object3D()
world_f1 = Modia3D.Object3D(world , r=[0.5*Lx, 0.0, Lz/2])

green = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(Lx,Ly,Lz), "Aluminium", vgreen),[ [-Lx/2,0,0],[ Lx/2,0,0] ])
rev1  = Modia3D.Revolute(world, green.frames[1]; axis = 2)
green2 = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(Lx,Ly,Lz), "Aluminium", vgreen),[ [-Lx/2,0,0],[ Lx/2,0,0] ])
rev12  = Modia3D.Revolute(green.frames[2], green2.frames[1]; axis = 1)
gelb = Modia3D.Part(Modia3D.Solid(Modia3D.SolidBox(Lx,Ly,Lz), "Aluminium", vyellow),[ [-Lx/2,0,0],[ Lx/2,0,0] ])
rev2 = Modia3D.Revolute(world_f1, gelb.frames[1])

sinWave     = SinWave()
cosWave     = CosWave()
linearSig   = LinearSig()
sigGreen    = Modia3D.SignalToFlangeAngle(sinWave.y1)
sigGreen2   = Modia3D.SignalToFlangeAngle(cosWave.y)
sigGelb     = Modia3D.SignalToFlangeAngle(cosWave.y)

Modia3D.connect(sigGreen,  rev1)
Modia3D.connect(sigGreen2, rev12)
Modia3D.connect(sigGelb,   rev2)
end
pendulum = Pendulum(sceneOptions=Modia3D.SceneOptions(visualizeFrames=true,defaultFrameLength=0.3, enableContactDetection=false))

#=
Modia3D.connect(pendulum.sigGreen,  pendulum.rev1)
Modia3D.connect(pendulum.sigGreen2, pendulum.rev12)
Modia3D.connect(pendulum.sigGelb,   pendulum.rev2)
=#

model    = Modia3D.SimulationModel(pendulum, analysis=ModiaMath.KinematicAnalysis)
# ModiaMath.print_ModelVariables(model)
result   = ModiaMath.simulate!(model, stopTime=5.0)

ModiaMath.plot(result, ["sinWave.y1", "rev1.phi", "cosWave.y", "rev12.phi", "linearSig.y", "rev2.phi"] )

println("... success of Test_KinematicRevoluteJoints.jl!")
end


text:
"""
    unitvector(v)

Normalize input vector `v`.
"""
function unitvector(v)
    n = norm(v)
    if n > 0
        return v / n
    else
        # zero vector
        return v
    end
end

"""
    random_in_unit_disk()

Generate a uniformly random point within the unit disk.
"""
function random_in_unit_disk()
    while true
        p = 2*rand(2) .- 1
        if dot(p, p) < 1
            return p
        end
    end
end

"""
    random_in_unit_sphere()

Generate a uniformly random point within the unit sphere.
"""
function random_in_unit_sphere()
    while true
        p = 2*rand(3) .- 1
        if dot(p, p) < 1
            return p
        end
    end
end


text:
"""
Convert DynamcHMC samples to a chain
* `posterior`: an array of NamedTuple consisting of mcmc samples
"""
function nptochain(posterior,tune)
    Np = length(vcat(posterior[1]...))+1 #include lf_eps
    Ns = length(posterior)
    a3d = Array{Float64,3}(undef,Ns,Np,1)
    ϵ=tune.ϵ
    for (i,post) in enumerate(posterior)
        temp = Float64[]
        for p in post
            push!(temp,values(p)...)
        end
        push!(temp,ϵ)
        a3d[i,:,1] = temp'
    end
    parameter_names = getnames(posterior)
    push!(parameter_names,"lf_eps")
    chns = MCMCChains.Chains(a3d,parameter_names,
        Dict(:internals => ["lf_eps"]))
    return chns
end

function getnames(post)
    nt = post[1]
    Np =length(vcat(nt...))
    parm_names = fill("",Np)
    cnt = 0
    for (k,v) in pairs(nt)
        N = length(v)
        if isa(v,Array)
            for i in 1:N
                cnt += 1
                parm_names[cnt] = string(k,"[",i,"]")
            end
        else
            cnt+=1
            parm_names[cnt] = string(k)
        end
    end
    return parm_names
end


text:
push!(LOAD_PATH, @__DIR__)

using TSPD
using Statistics

function read_data(filename, n)
    dir_path = @__DIR__
    inst_dir = "../../ExactSolutions/instances/n$(n)/"
    filepath = joinpath(dir_path, inst_dir, filename)

    f = open(filepath, "r")
    readline(f)
    speed_of_truck = parse(Float64, readline(f))
    readline(f)
    speed_of_drone = parse(Float64, readline(f))
    readline(f)
    n_nodes = parse(Int, readline(f))
    readline(f)
    depot_line = split(readline(f))
    depot_x = parse(Float64, depot_line[1])
    depot_y = parse(Float64, depot_line[2])
    readline(f)
    x = Vector{Float64}(undef, n_nodes)
    y = Vector{Float64}(undef, n_nodes)
    x[1] = depot_x
    y[1] = depot_y
    for i in 2:n_nodes
        cx, cy, _ = split(readline(f))
        x[i] = parse(Float64, cx)
        y[i] = parse(Float64, cy)
    end

    return x, y, speed_of_truck, speed_of_drone
end

exact_opt_obj = Dict{Int, Vector{Float64}}()
# n11
exact_opt_obj[11] = Float64[
    221.18876576478925,
    205.76050725572097,
    192.96313461174037,
    241.25592289521398,
    248.1379946498235,
    217.68894293889753,
    237.34013623078425,
    214.76536428997835,
    256.33972821148967,
    227.90300661076967
]

# n15
exact_opt_obj[15] = Float64[
    260.19649903254805,
    279.8748849328301,
    278.2501294150187,
    227.4927030824479,
    233.0058389136502,
    252.9062262214419,
    288.5045658852292,
    277.4202725739683,
    280.23702686030765,
    271.52779790946164
]

function test_agatz(n)
    println("Testing for n =", n)

    obj_val = Float64[]
    for i in 1:10
        filename = "uniform-$(i)-n$(n).txt"
        tsp_x, tsp_y, truck_cost_factor, drone_cost_factor = read_data(filename, n)

        val, _, _ = solve_tspd(tsp_x, tsp_y, truck_cost_factor, drone_cost_factor, method="TSP-ep-all")

        println(val)
        push!(obj_val, val)
    end

    opt_val = exact_opt_obj[n] :: Vector{Float64}
    gap = (obj_val - opt_val) ./ opt_val * 100
    println("average gap = ", mean(gap), " %")
    println("maximum gap = ", maximum(gap), " %")
end

println("-----------------------------")
t0 = time()
test_agatz(11)
t1 = time()
println("average time = ", (t1-t0) / 10, " seconds")
println("-----------------------------")
t0 = time()
test_agatz(15)
t1 = time()
println("average time = ", (t1-t0) / 10, " seconds")
println("-----------------------------")

# let
#     x = rand(1:100, 10) .* 1.0
#     y = rand(1:100, 10) .* 1.0
#     x[1] = rand()
#     y[1] = rand()
#     @show tspd_heuristic(x, y, 1, 0.5, method="TSP-ep")
#     @show tspd_heuristic(x, y, 1, 0.5, method="TSP-ep-1p")
#     @show tspd_heuristic(x, y, 1, 0.5, method="TSP-ep-2p")
#     @show tspd_heuristic(x, y, 1, 0.5, method="TSP-ep-2opt")
#     @show tspd_heuristic(x, y, 1, 0.5, method="TSP-ep-all", ep="Agatz")
#     @show tspd_heuristic(x, y, 1, 0.5, method="TSP-ep-all", ep="Poikonen")
#     @show branch_and_bound(x, y, 1, 0.5)
# end

# aep_vals = Float64[]
# pep_vals = Float64[]
# bab_vals = Float64[]
# for k in 1:10
#     x = rand(1:100, 11) .* 1.0
#     y = rand(1:100, 11) .* 1.0
#     x[1] = rand()
#     y[1] = rand()
#     aep, _, _ = tspd_heuristic(x, y, 1, 0.5, method="TSP-ep-all", ep="Agatz")
#     pep, _, _ = tspd_heuristic(x, y, 1, 0.5, method="TSP-ep", ep="Poikonen")
#     bab, _, _ = branch_and_bound(x, y, 1, 0.5)

#     push!(aep_vals, aep)
#     push!(pep_vals, pep)
#     push!(bab_vals, bab)
# end
# @show mean(aep_vals), mean(pep_vals), mean(bab_vals)


text:
#I Do not posses the scala compiler, so if anything is wrong please correct it.

println("Hello World")


text:
struct ppFPFState{T, TF1, TF2} <: AbstractFilterState
    ensemble::UnweightedParticleEnsemble{T}
    eq_dN::PoissonEquation{T, TF1}
    eq_dt::PoissonEquation{T, TF2}
    function ppFPFState(ens::UnweightedParticleEnsemble{T}, eq_dN::PoissonEquation{T, TF1}, eq_dt::PoissonEquation{T, TF2}) where {T, TF1, TF2}
        update!(eq_dN, ens)
        update!(eq_dt, ens)
        new{T, TF1, TF2}(ens, eq_dN, eq_dt)
    end
end

struct ppFPF{TS<:HiddenStateModel, TO<:ObservationModel, T1<:GainEstimationMethod, T2<:ParticleFlowMethod} <: AbstractFilteringAlgorithm{ContinuousTime, ContinuousTime}
    state_model::TS
    obs_model::TO
    gain_method::T1
    flow_method::T2
    N::Int
    function ppFPF(st_mod::HiddenStateModel{S1, ContinuousTime}, ob_mod::ObservationModel{S1, S2, ContinuousTime}, gain_method::GainEstimationMethod, flow_method::ParticleFlowMethod, N::Int) where {S1, S2}
        if ob_mod isa CountingObservationModel
            return new{typeof(st_mod), typeof(ob_mod), typeof(gain_method), typeof(flow_method)}(st_mod, ob_mod, gain_method, flow_method, N)
        else
            error("Error: the ppFPF algorithm is only defined for counting process observations.")
        end
    end
end

#####################
###### METHODS ######
#####################

initial_condition(filter::ppFPF) = initial_condition(state_model(filter))
no_of_particles(filter::ppFPF) = filter.N
no_of_particles(st::ppFPFState) = no_of_particles(st.ensemble)
state_model(filter::ppFPF) = filter.state_model
obs_model(filter::ppFPF) = filter.obs_model
gain_estimation_method(filter::ppFPF) = filter.gain_method
flow_method(filter::ppFPF) = filter.flow_method

mean(st::ppFPFState) = mean(st.ensemble)
cov(st::ppFPFState)  = cov(st.ensemble)
var(st::ppFPFState)  = var(st.ensemble)

function Base.show(io::IO, ::MIME"text/plain", filter::ppFPF)
    print(io, "Point-process feedback particle filter algorithm
    hidden state model:                     ", state_model(filter),"
    observation model:                      ", obs_model(filter),"
    ensemble size:                          ", no_of_particles(filter),"
    gain estimation method:                 ", gain_estimation_method(filter),"
    particle flow method:                   ", flow_method(filter))
end

function Base.show(io::IO, filter::ppFPF)
    print(io, "ppFPF with ", no_of_particles(filter)," particles,", gain_estimation_method(filter), ", and ", flow_method(filter))
end

######################
### MAIN ALGORITHM ###
######################

function initialize(filter::ppFPF)
    return ppFPFState(filter)
end

function update!(filter_state::ppFPFState, filt_algo::ppFPF, dY, dt)
    propagate!(filter_state, filt_algo, dt)
    assimilate!(dY, filter_state, filt_algo, dt)
end

function propagate!(filter_state::ppFPFState, filt_algo::ppFPF, dt)
    propagate!(filter_state.ensemble, filt_algo.state_model, dt)
end

function assimilate!(dN, filter_state::ppFPFState, filt_algo::ppFPF, dt)
    ensemble  = filter_state.ensemble
    eq_dt     = filter_state.eq_dt
    eq_dN     = filter_state.eq_dN
    method    = gain_estimation_method(filt_algo)
    fl_method = flow_method(filt_algo)

    # dt update
    solve!(eq_dt, method)

    ensemble.positions .+= dt .* view(eq_dt.gain, :, :, 1)

    # event update
    logl(x) = LinearAlgebra.dot(eq_dN.h(x), dN)
    flow!(ensemble, logl, fl_method)

    update!(eq_dt, ensemble)
    return ensemble.positions
end

################################
### CONVENIENCE CONSTRUCTORS ###
################################

function ppFPFState(filter::ppFPF)
    N = no_of_particles(filter)
    return ppFPFState(state_model(filter), obs_model(filter), N)
end

function ppFPFState(filt_prob::AbstractFilteringProblem, N)
    return ppFPFState(state_model(filt_prob), obs_model(filt_prob), N)
end

function ppFPFState(state_model::HiddenStateModel, obs_model::ObservationModel, N)
    ens    = UnweightedParticleEnsemble(hcat([initialize(state_model) for i in 1:N]...))
    eq_dt  = PoissonEquation(x -> -sum(observation_function(obs_model)(x)), ens)
    eq_dN  = PoissonEquation(x -> log.(observation_function(obs_model)(x)), ens)
    return ppFPFState(ens, eq_dN, eq_dt)
end

function ppFPF(filt_prob::AbstractFilteringProblem, method::GainEstimationMethod, fl_method::ParticleFlowMethod, N)
    st_mod = state_model(filt_prob)

    if obs_model(filt_prob) isa CountingObservationModel
        ob_mod = obs_model(filt_prob)
    else
        error("Error: the EKSPF algorithm is only defined for counting process observations.")
    end

    return ppFPF(st_mod, ob_mod, method, fl_method, N)
end

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
using .Hidato

hidat = """
__ 33 35 __ __  .  .  .
__ __ 24 22 __  .  .  .
__ __ __ 21 __ __  .  .
__ 26 __ 13 40 11  .  .
27 __ __ __  9 __  1  .
 .  . __ __ 18 __ __  .
 .  .  .  . __  7 __ __
 .  .  .  .  .  .  5 __"""

const kingmoves = [[-1, -1], [-1, 0], [-1, 1], [0, -1], [0, 1], [1, -1], [1, 0], [1, 1]]

board, maxmoves, fixed, starts = hidatoconfigure(hidat)
printboard(board)
hidatosolve(board, maxmoves, kingmoves, fixed, starts[1][1], starts[1][2], 1)
printboard(board)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Testing
# =======
#
# Utilities to assist testing of BioJulia packages.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioCore.jl/blob/master/LICENSE.md

module Testing

import YAML

const FMT_SPECIMEN_PATH = joinpath(dirname(dirname(@__FILE__)), "BioFmtSpecimens")

"""
    get_bio_fmt_specimens(checkout = "master", auto_checkout = true, fresh = false)

Install/update and return the path of BioJulia's biological data format
specimen archive.

When the BioFmtSpecimens archive is fetched from the web, the branch or tag
specified by  `checkout` is checked out for use. Unless, `auto_checkout` is
true, in which case, the latest tagged release of the BioFmtSpecimens archive
will be checked out.

If `fresh` is set to true, this will force a deletion of any currently installed
BioFmtSpecimens archive repository, and fetch it from the web again. This may
be useful if updating the installed BioFmtSpecimens archive is problematic.
"""
function get_bio_fmt_specimens(checkout = "master", auto_checkout = true, fresh = false)
    if fresh
        rm(FMT_SPECIMEN_PATH, force = true, recursive = true)
    end
    if !isdir(FMT_SPECIMEN_PATH)
        run(`git clone https://github.com/BioJulia/BioFmtSpecimens.git $(FMT_SPECIMEN_PATH)`)
    end
    cd(FMT_SPECIMEN_PATH) do
        if auto_checkout
            checkout = open(readline, `git describe --tags`, "r+")
        end
        run(`git fetch origin`)
        run(`git checkout $(checkout)`)
    end
    return FMT_SPECIMEN_PATH
end

"""
    bio_fmt_specimens(format::String, fn::Function, checkout = "master", auto_checkout = true, fresh = false)

Get the paths for file format specimen files from BioJulia's biological data
format specimen archive.

Will return a vector of paths for specimen files of a given `format` (e.g. FASTA)
which satisfy the filter `fn`. the input for `fn` should be a single argument,
which is a `Dict{Any, Any}`. Each `Dict{Any, Any}` represents a format file
specimen from BioJulia's biological data format specimen archive, and has the
following fields:

* **"filename"**: Specimen filename.
* **"valid"**: `true` or `false`, indicates whether the example conforms to the format.
* **"origin"** The contributor or source from which a specimen was taken.
* **"tags"** Zero or more words used to group specimens by shared features.
* **"comments"** (Optional) Any additional information that might be of interest.

When the BioFmtSpecimens archive is fetched from the web or updated, the branch or tag
specified by  `checkout` is checked out for use. Unless, `auto_checkout` is
true, in which case, the latest tagged release of the BioFmtSpecimens archive
will be checked out.

If `fresh` is set to true, this will force a deletion of any currently installed
BioFmtSpecimens archive repository, and fetch it from the web again. This may
be useful if updating the installed BioFmtSpecimens archive is problematic.

```@example
# Get paths for FASTA format specimens which are examples of a valid file.
bio_fmt_specimens("FASTA", (x) -> x["valid"] == true)
```
"""
function bio_fmt_specimens(format::String, fn::Function, checkout = "master", auto_checkout = true, fresh = false)
    get_bio_fmt_specimens(checkout, auto_checkout, fresh)
    specimens = YAML.load_file(joinpath(FMT_SPECIMEN_PATH, format, "index.yml"))
    filtered_specimens = Vector{String}(undef, length(specimens))
    fsi = 0
    for specimen in specimens
        if fn(specimen)
            fsi += 1
            filtered_specimens[fsi] = joinpath(FMT_SPECIMEN_PATH, format, specimen["filename"])
        end
    end
    resize!(filtered_specimens, fsi)
    return filtered_specimens
end

function random_array(n::Integer, elements, probs)
    cumprobs = cumsum(probs)
    x = Vector{eltype(elements)}(undef, n)
    for i in 1:n
        x[i] = elements[searchsorted(cumprobs, rand()).start]
    end
    return x
end

function random_seq(n::Integer, nts, probs)
    cumprobs = cumsum(probs)
    x = Vector{Char}(undef, n)
    for i in 1:n
        x[i] = nts[searchsorted(cumprobs, rand()).start]
    end
    return String(x)
end

function random_dna(n, probs=[0.24, 0.24, 0.24, 0.24, 0.04])
    return random_seq(n, ['A', 'C', 'G', 'T', 'N'], probs)
end

function random_rna(n, probs=[0.24, 0.24, 0.24, 0.24, 0.04])
    return random_seq(n, ['A', 'C', 'G', 'U', 'N'], probs)
end

function random_aa(len)
    return random_seq(len,
        ['A', 'R', 'N', 'D', 'C', 'Q', 'E', 'G', 'H', 'I',
         'L', 'K', 'M', 'F', 'P', 'S', 'T', 'W', 'Y', 'V', 'X' ],
        push!(fill(0.049, 20), 0.02))
end

function intempdir(fn::Function, parent=tempdir())
    dirname = mktempdir(parent)
    try
        cd(fn, dirname)
    finally
        rm(dirname, recursive=true)
    end
end

function random_interval(minstart, maxstop)
    start = rand(minstart:maxstop)
    return start:rand(start:maxstop)
end

end # Module Testing


text:
# Core debugging logic.
# Hopefully someday much of this will be replaced by Gallium.

const VarnameType = Tuple{Vararg{Union{Symbol,Expr}}}  # Expr covers `foo(x, (a,b), y)` destructured-tuple signatures
struct Stored
    method::Method
    varnames::VarnameType
    varvals

    function Stored(m, names, vals)
        new(m, names, safe_deepcopy(vals...))
    end
end

const stashed     = Ref{Any}(nothing)
const stored      = Dict{UUID,Stored}()              # UUID => store data
const storefunc   = Dict{UUID,Function}()            # UUID => function that puts inputs into `stored`
const storemap    = Dict{Tuple{Method,Bool},UUID}()  # (method, overwrite) => UUID

struct StopException   <: Exception end
struct StashingFailed  <: Exception end   # stashing saves the function and arguments from the caller (transient)
struct StorageFailed   <: Exception end   # storage puts callee values into `stored` (lasts until `stored` is cleared)
struct DefMissing      <: Exception
    method::Method
    exception
end
struct SignatureError  <: Exception
    method::Method
end
struct StepException   <: Exception
    msg::String
end
struct EvalException   <: Exception
    exprstring
    exception
end

"""
    Rebugger.clear()

Clear internal data. This deletes storage associated with stored variables, but also
forces regeneration of capture methods, which can be handy while debugging Rebugger itself.
"""
function clear()
    stashed[] = nothing
    empty!(stored)
    empty!(storefunc)
    empty!(storemap)
    nothing
end

### Stacktraces

"""
    r = linerange(expr, offset=0)

Compute the range of lines occupied by `expr`.
Returns `nothing` if no line statements can be found.
"""
function linerange(def::ExLike, offset=0)
    start, haslinestart = findline(def, identity)
    stop, haslinestop  = findline(def, Iterators.reverse)
    (haslinestart & haslinestop) && return (start+offset):(stop+offset)
    return nothing
end

function findline(ex, order)
    ex.head == :line && return ex.args[1], true
    for a in order(ex.args)
        a isa LineNumberNode && return a.line, true
        if a isa ExLike
            ln, hasline = findline(a, order)
            hasline && return ln, true
        end
    end
    return 0, false
end

"""
    usrtrace, defs = pregenerated_stacktrace(trace, topname=:capture_stacktrace)

Generate a list of methods `usrtrace` and their corresponding definition-expressions `defs`
from a stacktrace.
Not all methods can be looked up, but this attempts to resolve, e.g., keyword-handling methods
and so on.
"""
function pregenerated_stacktrace(trace; topname = :capture_stacktrace)
    usrtrace, defs = Method[], RelocatableExpr[]
    methodsused = Set{Method}()
    for (i, sf) in enumerate(trace)
        sf.func == topname && break  # truncate at the chosen spot
        sf.func ∈ notrace && continue
        mi = sf.linfo
        if mi isa Core.MethodInstance
            method = mi.def
            # Set up tracking, if necessary
            file = String(sf.file)
            if !haskey(Revise.fileinfos, file)
                ret = Revise.find_file(file, method.module)
                ret == nothing && continue
                file, recipemod = ret
                if !haskey(Revise.fileinfos, file)
                    try
                        @info "tracking $recipemod"
                        Revise.track(recipemod)
                    catch err
                        err isa Revise.GitRepoException && continue
                        rethrow(err)
                    end
                end
            end
            haskey(Revise.fileinfos, file) || continue
            fi = Revise.fileinfos[file]
            Revise.maybe_parse_from_cache!(fi, file)
            funcname = String(sf.func)
            if startswith(funcname, '#')
                # This is a generated method, perhaps it's a keyword function handler
                # Look for it by line number
                defmap = fi.fm[method.module].defmap
                for (def, info) in defmap
                    info == nothing && continue
                    sigts, offset = info
                    r = linerange(def, offset)
                    r == nothing && continue
                    if sf.line ∈ r
                        mths = Base._methods_by_ftype(last(sigts), -1, typemax(UInt))
                        if length(mths) == 1
                            m = mths[1][3]
                            if m ∉ methodsused
                                push!(defs, def)
                                push!(usrtrace, m)
                                push!(methodsused, m)
                                break
                            end
                        end
                    end
                end
            else
                method ∈ methodsused && continue
                def = Revise.get_def(method)
                def isa ExLike || continue
                push!(defs, def)
                push!(usrtrace, method)
            end
        end
    end
    return usrtrace, defs
end

"""
    uuids = capture_stacktrace(mod, command)

Execute `command` in module `mod`. `command` must throw an error.
Then instrument the methods in the stacktrace so that their input
variables are stored in `Rebugger.stored`.
After storing the inputs, restore the original methods.

Since this requires two `eval`s of `command`, usage should be limited to
deterministic expressions that always result in the same call chain.
"""
function capture_stacktrace(mod::Module, command::Expr)
    errored = true
    trace = try
        Core.eval(mod, command)
        errored = false
    catch
        stacktrace(catch_backtrace())
    end
    errored || error("$command did not throw an error")
    usrtrace, defs = pregenerated_stacktrace(trace)
    println(stderr, "Captured elements of stacktrace:")
    show(stderr, MIME("text/plain"), usrtrace)
    length(unique(usrtrace)) == length(usrtrace) || @error "the same method appeared twice, not supported. Try stepping into the command."
    uuids = UUID[]
    capture_stacktrace!(uuids, usrtrace, defs) do
        Core.eval(mod, command)
    end
    uuids
end
capture_stacktrace(command::Expr) = capture_stacktrace(Main, command)

function capture_stacktrace!(f::Function, uuids::Vector, usrtrace, defs)
    if isempty(usrtrace)
        # We've finished modifying all the methods, time to run the command
        try
            f()
            @warn "traced method did not throw an error"
        catch
        end
        return
    end
    method, def = usrtrace[end], defs[end]
    if def != nothing
        push!(uuids, method_capture_from_callee(method, def; overwrite=true))
    end
    # Recurse up the stack until we get to the top...
    pop!(usrtrace)
    pop!(defs)
    capture_stacktrace!(f, uuids, usrtrace, defs)
    # ...after which it will call the erroring function and come back here.
    # Having come back, restore the original definition
    if def != nothing
        eval_noinfo(method.module, def)  # unfortunately this doesn't restore the original method as a viable key to storemap
    end
    return uuids
end

### Stepping

function stepin(io)
    @assert Rebugger.stashed[] === nothing
    # Step 1: rewrite the command to stash the call function and its arguments.
    prepare_caller_capture!(io)
    capexpr, stop = Meta.parse(content(io), 1)
    try
        Core.eval(Main, capexpr)
        throw(StashingFailed())
    catch err
        err isa StashingFailed && rethrow(err)
        if !(err isa StopException)
            throw(EvalException(content(io), err))
        end
    end
    f, args, kwargs = Rebugger.stashed[]
    Rebugger.stashed[] = nothing
    # Step 2: determine which method is called, and if need be create a function
    # that captures all of the callee's inputs. (This allows us to capture default arguments,
    # keyword arguments, and type parameters.)
    method = which(f, Base.typesof(args...))
    uuid = get(storemap, (method, false), nothing)
    if uuid === nothing
        uuid = method_capture_from_callee(method; overwrite=false)
    end
    # Step 3: execute the command to store the inputs.
    fcapture = storefunc[uuid]
    tv, decls = Base.arg_decl_parts(method)
    if !isempty(decls[1][1])
        # This is a call-overloaded method, prepend the calling object
        args = (f, args...)
    end
    try
        Base.invokelatest(fcapture, args...; kwargs...)
        throw(StorageFailed())   # this should never happen
    catch err
        err isa StopException || rethrow(err)
    end
    return uuid, generate_let_command(method, uuid)
end

"""
    callexpr = prepare_caller_capture!(io)

Given a buffer `io` representing a string and "point" (the seek position) set at a call expression,
replace the call with one that stashes the function and arguments of the call.

For example, if `io` has contents

    <some code>
    if x > 0.5
        ^fcomplex(x, 2; kw1=1.1)
        <more code>

where in the above `^` indicates `position(s)` ("point"), rewrite this as

    <some code>
    if x > 0.5
        Main.Rebugger.stashed[] = (fcomplex, (x, 2), (kw1=1.1,))
        throw(Rebugger.StopException())
        <more code>

(Keyword arguments do not affect dispatch and hence are not stashed.)
Consequently, if this is `eval`ed and execution reaches "^", it causes the arguments
of the call to be placed in `Rebugger.stashed`.

`callexpr` is the original (unmodified) expression specifying the call, i.e.,
`fcomplex(x, 2; kw1=1.1)` in this case.

This does the buffer-preparation for *caller* capture.
For *callee* capture, see [`method_capture_from_callee`](@ref),
and [`stepin`](@ref) which puts these two together.
"""
function prepare_caller_capture!(io)  # for testing, needs to work on a normal IO object
    start = position(io)
    callstring = content(io, start=>bufend(io))
    callexpr, len = Meta.parse(callstring, 1; raise=false)
    callexpr == nothing && throw(StepException("Got empty expression from $callstring"))
    isa(callexpr, Expr) || throw(StepException("Rebugger can only step into expressions, got $callexpr"))
    if callexpr.head == :error
        iend = len
        for i = 1:2
            iend = prevind(callstring, iend)
        end
        callstring = callstring[1:iend]
        callexpr, len = Meta.parse(callstring, 1)
    end
    if callexpr.head == :tuple && !(startswith(callstring, "tuple") || startswith(callstring, "("))
        # An expression like foo(bar(x)..., 1) where point is positioned at bar
        callexpr = callexpr.args[1]
    end
    if callexpr.head == :ref
        callexpr = Expr(:call, :getindex, callexpr.args...)
    elseif callexpr.head == :(=) && isa(callexpr.args[1], Expr) && callexpr.args[1].head == :ref
        ref, val = callexpr.args
        callexpr = Expr(:call, :setindex!, ref.args[1], val, ref.args[2:end]...)
    elseif (callexpr.head == :&& || callexpr.head == :||) && isa(callexpr.args[1], Expr)
        callexpr = callexpr.args[1]
    elseif callexpr.head == :...
        callexpr = callexpr.args[1]
    end
    # Must be a call or broadcast
    ((callexpr.head == :call) | (callexpr.head == :.)) || throw(Meta.ParseError("point must be at a call expression, got $callexpr"))
    if callexpr.head == :call
        fname, args = callexpr.args[1], callexpr.args[2:end]
    else
        fname, args = :broadcast, [callexpr.args[1], callexpr.args[2].args...]
    end
    # In the edited callstring separate any kwargs now. They don't affect dispatch.
    kwargs = []
    if length(args) >= 1 && isa(args[1], Expr) && args[1].head == :parameters
        # foo(x; kw1=1, ...) syntax    (with the semicolon)
        append!(kwargs, popfirst!(args).args)
    end
    while !isempty(args) && isa(args[end], Expr) && args[end].head == :kw
        # foo(x, kw1=1, ...) syntax    (with a comma, no semicolon)
        push!(kwargs, pop!(args))
    end
    captureexpr = quote
        Main.Rebugger.stashed[] = ($fname, (($(args...)),), Main.Rebugger.kwstasher(; $(kwargs...)))
        throw(StopException())
    end
    # Now insert this in place of the marked call
    # Unfortunately we have to convert to a string and there are scoping issues
    capturestr = string(captureexpr, '\n')
    regexunscoped = r"(?<!\.)StopException"
    capturestr = replace(capturestr, regexunscoped=>(s->"Rebugger."*s))
    regexscoped   = r"(?<!\.)Rebugger\.StopException"
    capturestr = replace(capturestr, regexscoped=>(s->"Main."*s))
    edit_splice!(io, start=>start+len-1, capturestr)
    return callexpr
end

"""
    uuid = method_capture_from_callee(method; overwrite::Bool=false)

Create a version of `method` that stores its inputs in `Main.Rebugger.stored`.
For a method

    function fcomplex(x::A, y=1, z=""; kw1=3.2) where A<:AbstractArray{T} where T
        <body>
    end

if `overwrite=false`, this generates a new method

    function hidden_fcomplex(x::A, y=1, z=""; kw1=3.2) where A<:AbstractArray{T} where T
        Main.Rebugger.stored[uuid] = Main.Rebugger.Stored(fcomplex, (:x, :y, :z, :kw1, :A, :T), deepcopy((x, y, z, kw1, A, T)))
        throw(StopException())
    end

(If a `uuid` already exists for `method` from a previous call to `method_capture_from_callee`,
it will simply be returned.)

With `overwrite=true`, there are two differences:

- it replaces `fcomplex` rather than defining `hidden_fcomplex`
- rather than throwing `StopException`, it re-inserts `<body>` after the line performing storage

The returned `uuid` can be used for accessing the stored data.
"""
function method_capture_from_callee(method, def; overwrite::Bool=false)
    uuid = get(storemap, (method, overwrite), nothing)
    uuid != nothing && return uuid
    sigr, body = get_signature(def), unquote(funcdef_body(def))
    sigr == nothing && throw(SignatureError(method))
    sigex = convert(Expr, sigr)
    if sigex.head == :(::)
        sigex = sigex.args[1]  # return type declaration
    end
    methname, argnames, kwnames, paramnames = signature_names!(sigex)
    # Check for call-overloading method, e.g., (obj::ObjType)(x, y...) = <body>
    callerobj = nothing
    if methname isa Expr && methname.head == :(::)
        @assert length(methname.args) == 2
        callerobj = methname
        argnames = (methname.args[1], argnames...)
        methname = methname.args[2]
        if methname isa Expr
            if methname.head == :curly
                methname = methname.args[1]
            else
                dump(methname)
                error("unexpected call-overloading type")
            end
        end
    end
    allnames = (argnames..., kwnames..., paramnames...)
    qallnames = QuoteNode.(allnames)
    uuid = uuid1()
    uuidstr = string(uuid)
    storeexpr = :(Main.Rebugger.setstored!($uuidstr=>Main.Rebugger.Stored($method, ($(qallnames...),), ($(allnames...),)) ) )
    capture_body = overwrite ? quote
        $storeexpr
        $body
    end : quote
        $storeexpr
        throw(Main.Rebugger.StopException())
    end
    capture_name = try
        _gensym(methname)
    catch
        nothing
    end
    capture_name == nothing && (dump(methname); dump(sigex); error("couldn't gensym"))
    # capture_name = _gensym(methname)
    mod = method.module
    capture_function = Expr(:function, overwrite ? sigex : rename_method(sigex, capture_name, callerobj), capture_body)
    result = Core.eval(mod, capture_function)
    if !overwrite
        storefunc[uuid] = result
    end
    storemap[(method, overwrite)] = uuid
    return uuid
end
function method_capture_from_callee(method; kwargs...)
    # Could use a default arg above but this generates a more understandable error message
    local def
    try
        def = get_def(method)
    catch err
        throw(DefMissing(method, err))
    end
    def == nothing && throw(DefMissing(method, nothing))
    method_capture_from_callee(method, def; kwargs...)
end

function generate_let_command(method::Method, uuid::UUID)
    s = stored[uuid]
    @assert method == s.method
    argstring = '(' * join(s.varnames, ", ") * (length(s.varnames)==1 ? ",)" : ')')
    body = convert(Expr, striplines!(deepcopy(funcdef_body(get_def(method)))))
    return """
        @eval $(method.module) let $argstring = Main.Rebugger.getstored(\"$uuid\")
        $body
        end"""
end
function generate_let_command(uuid::UUID)
    s = stored[uuid]
    generate_let_command(s.method, uuid)
end

"""
    args_and_types = Rebugger.getstored(uuid)

Retrieve the values of stored arguments and type-parameters from the store specified
`uuid`. This makes a copy of values, so as to be safe for repeated execution of methods
that modify their inputs.
"""
getstored(uuidstr::AbstractString) = safe_deepcopy(Main.Rebugger.stored[UUID(uuidstr)].varvals...)

function setstored!(p::Pair{S,Stored}) where S<:AbstractString
    uuidstr, val = p.first, p.second
    Main.Rebugger.stored[UUID(uuidstr)] = val
end

kwstasher(; kwargs...) = kwargs

### Utilities

"""
    fname, argnames, kwnames, parameternames = signature_names!(sigex::Expr)

Return the function name `fname` and names given to its arguments, keyword arguments,
and parameters, as specified by the method signature-expression `sigex`.

`sigex` will be modified if some of the arguments are unnamed.

# Examples

```jldoctest; setup=:(using Rebugger)
julia> Rebugger.signature_names!(:(complexargs(w::Ref{A}, @nospecialize(x::Integer), y, z::String=""; kwarg::Bool=false, kw2::String="", kwargs...) where A <: AbstractArray{T,N} where {T,N}))
(:complexargs, (:w, :x, :y, :z), (:kwarg, :kw2, :kwargs), (:A, :T, :N))

julia> ex = :(myzero(::Float64));     # unnamed argument

julia> Rebugger.signature_names!(ex)
(:myzero, (:__Float64_1,), (), ())

julia> ex
:(myzero(__Float64_1::Float64))
```
"""
function signature_names!(sigex::ExLike)
    # TODO: add parameter names
    argname(s::Symbol) = s
    function argname(ex::ExLike)
        if ex.head == :(...) && length(ex.args) == 1
            # varargs function
            ex = ex.args[1]
            ex isa Symbol && return ex
        end
        (ex.head == :macrocall || ex.head == :meta) && return argname(ex.args[end])  # @nospecialize
        ex.head == :kw && return argname(ex.args[1])  # default arguments
        ex.head == :tuple && return ex    # tuple-destructuring argument
        ex.head == :(::) || throw(ArgumentError(string("expected :(::) expression, got ", ex)))
        arg = ex.args[1]
        if length(ex.args) == 1 && (arg isa Symbol)
            # This argument has a type but no name
            return arg, true
        end
        if isa(arg, Expr) && arg.head == :curly && arg.args[1] == :Type
            # Argument of the form ::Type{T}
            return arg.args[2], false
        end
        return arg
    end
    paramname(s::Symbol) = s
    function paramname(ex::ExLike)
        ex.head == :(<:) && return paramname(ex.args[1])
        throw(ArgumentError(string("expected parameter expression, got ", ex)))
    end

    kwnames, parameternames = (), []
    while sigex.head == :where
        parameternames = [paramname.(sigex.args[2:end])..., parameternames...]
        sigex = sigex.args[1]
    end
    name = sigex.args[1]
    offset = 1
    if length(sigex.args) > 1 && isa(sigex.args[2], ExLike) && sigex.args[2].head == :parameters
        # keyword arguments
        kwnames = tuple(argname.(sigex.args[2].args)...)
        offset += 1
    end
    # Argnames. For any unnamed arguments we have to generate a name.
    empty!(usdict)
    argnames = Union{Symbol,Expr}[]
    for i = offset+1:length(sigex.args)
        arg = sigex.args[i]
        name = argname(arg)
        if name isa Tuple
            name, genname = name
            if genname
                # This argument is missing a real name
                argt = name
                name = genunderscored(argt)
                sigex.args[i] = :($name::$argt)
            else
                # This is a ::Type{T} argument. We should remove this from the list of parameters
                parameternames = filter(!isequal(name), parameternames)
            end
        end
        push!(argnames, name)
    end

    return sigex.args[1], tuple(argnames...), kwnames, tuple(parameternames...)
end

function rename_method!(ex::ExLike, name::Symbol, callerobj)
    sig = ex.head ∈ (:call, :where) ? ex : get_signature(ex)
    while isa(sig, ExLike) && sig.head == :where
        sig = sig.args[1]
    end
    sig.head == :call || (dump(ex); throw(ArgumentError(string("expected call expression, got ", ex))))
    sig.args[1] = name
    if callerobj != nothing
        # Call overloading, add an argument
        sig.args = [sig.args[1]; callerobj; sig.args[2:end]]
    end
    return ex
end
rename_method(ex::ExLike, name::Symbol, callerobj) = rename_method!(copy(ex), name, callerobj)

# Use to re-evaluate an expression without leaving "breadcrumbs" about where
# the eval is coming from. This is used below to prevent the re-evaluaton of an
# original method from being attributed to Rebugger itself in future backtraces.
eval_noinfo(mod::Module, ex::Expr) = ccall(:jl_toplevel_eval, Any, (Any, Any), mod, ex)
eval_noinfo(mod::Module, rex::RelocatableExpr) = eval_noinfo(mod, convert(Expr, rex))

function unquote(ex::Expr)
    if ex.head == :quote
        return Expr(:block, ex.args...)
    end
    ex
end
unquote(rex::RelocatableExpr) = unquote(convert(Expr, rex))

_gensym(sym::Symbol) = gensym(sym)
_gensym(q::QuoteNode) = _gensym(q.value)
_gensym(ex::Expr) = (@assert ex.head == :. && length(ex.args) == 2; _gensym(ex.args[2]))

const usdict = Dict{Symbol,Int}()
function genunderscored(sym::Symbol)
    n = get(usdict, sym, 0) + 1
    usdict[sym] = n
    return Symbol("__"*String(sym)*'_'*string(n))
end

const notrace = (:error, :throw)


text:
function test_relu_neuron(backend::Backend, T)
  println("-- Testing ReLU neuron on $(typeof(backend)){$T}...")

  eps = 1e-10
  data = rand(T, 3,4,5,6) - 0.5
  data_blob = make_blob(backend, data)
  neuron0 = Neurons.ReLU()
  neuron1 = Neurons.ReLU(1e-6)

  println("    > Forward")
  for (neuron, threshval) in [(neuron0, 0), (neuron1, 1e-6)]
    copy!(data_blob, data)
    forward(backend, neuron, data_blob)
    expected_data = max(data, threshval)
    got_data = zeros(T, size(data))
    copy!(got_data, data_blob)
    @test all(-eps .< got_data - expected_data .< eps)
  end

  println("    > Backward")
  grad = rand(T, size(data))
  grad_blob = make_blob(backend, grad)

  for (neuron, threshval) in [(neuron0, 0), (neuron1, 1e-6)]
    copy!(data_blob, data)
    forward(backend, neuron, data_blob)
    backward(backend, neuron, data_blob, grad_blob)
    expected_grad = grad .* (data .> threshval)
    got_grad = zeros(T, size(expected_grad))
    copy!(got_grad, grad_blob)

    @test all(-eps .< got_grad - expected_grad .< eps)
  end
end
function test_relu_neuron(backend::Backend)
  test_relu_neuron(backend, Float32)
  test_relu_neuron(backend, Float64)
end

if test_cpu
  test_relu_neuron(backend_cpu)
end
if test_gpu
  test_relu_neuron(backend_gpu)
end


text:
# # Basic demonstration of IteratedModel wrapper for a tree booster

using Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

using MLJ
using MLJIteration
using Dates

using Plots
pyplot(size = (600, 300*(sqrt(5) - 1)))

using Statistics
using Random
Random.seed!(123)

MLJ.color_off()

X, y = make_moons(1000, rng=123)

# Import an model type:

EvoTreeClassifier = @load EvoTreeClassifier verbosity=0

model = EvoTreeClassifier(rng=123, η=0.005)

iterations = Int[0]
losses = Float64[0]

function update_plot(loss)
    push!(losses, loss)
    plot(iterations[2:end],
         losses[2:end],
         xlim=[1,300],
         ylim=[0,0.5]) |> display
end

imodel = IteratedModel(model=model,
                       resampling=Holdout(rng=123),
                       measures=[brier_loss, log_loss],
                       iteration_parameter=:nrounds,
                       controls=[Step(2),
                                 WithIterationsDo(i->push!(iterations, i)),
                                 WithLossDo(update_plot),
                                 GL(200),
                                 TimeLimit(Second(30))])

mach = machine(imodel, X, y) |> fit!

using Literate #src
Literate.markdown(@__FILE__, @__DIR__, execute=false) #src
Literate.notebook(@__FILE__, @__DIR__, execute=true) #src


text:
module BinaryClassificationTools

include("evaluation/diagnostics.jl")
include("evaluation/roc_auc.jl")

end


text:
function fahrenheit_to_celsius(temp_f)
    temp_c = (temp_f - 32.0) * (5.0/9.0)
    return temp_c
end

temp_c = fahrenheit_to_celsius(100.0)
println(temp_c)


text:
using Circuits
using JuliaFormatter
using Test

@testset "Circuits.jl" begin
    # No tests yet... TODO
end

@testset "Formatting" begin
    @test format(Circuits; overwrite=false)
end