text:
"""
  Does basic testing of explicit jacobian calculation in parallel
"""
function test_jac_parallel()

  # SBPOmega
  fname = "input_vals_jac3dp.jl"
  fname2 = "input_vals_jac_tmp.jl"
  mesh, sbp, eqn, opts = run_solver(fname)

  MPI.Barrier(mesh.comm)
  test_jac_parallel_inner(mesh, sbp, eqn, opts)
  opts["preallocate_jacobian_coloring"] = true
  test_jac_parallel_inner(mesh, sbp, eqn, opts, is_prealloc_exact=true, set_prealloc=false)

  return nothing
end

add_func1!(EulerTests, test_jac_parallel, [TAG_SHORTTEST, TAG_JAC])

"""
  Does more thorough testing of jacobian calculation in parallel
"""
function test_jac_parallel_long()

  @testset "----- Testing jacobian assembly long -----" begin
    fname = "input_vals_jac3dp.jl"
    fname2 = "input_vals_jac_tmp.jl"

    myrank = MPI.Comm_rank(MPI.COMM_WORLD)

    # SBPGamma
    if myrank == 0
      opts_tmp = read_input_file(fname)
      opts_tmp["operator_type"] = "SBPGamma"
      make_input(opts_tmp, fname2)
    end
    MPI.Barrier(MPI.COMM_WORLD)
    mesh4, sbp4, eqn4, opts4 = run_solver(fname2)

    # SBPDiagonalE
    if myrank == 0
      opts_tmp = read_input_file(fname)
      opts_tmp["operator_type"] = "SBPDiagonalE"
      make_input(opts_tmp, fname2)
    end
    MPI.Barrier(MPI.COMM_WORLD)
    mesh5, sbp5, eqn5, opts5 = run_solver(fname2)

    if myrank == 0
      opts_tmp = read_input_file(fname)
      opts_tmp["operator_type"] = "SBPDiagonalE"
      opts_tmp["use_Minv"] = true
      make_input(opts_tmp, fname2)
    end
    MPI.Barrier(MPI.COMM_WORLD)
    mesh6, sbp6, eqn6, opts6 = run_solver(fname2)

    # SBPGamma ES
    if myrank == 0
      opts_tmp = read_input_file(fname)
      opts_tmp["operator_type"] = "SBPGamma"
      opts_tmp["volume_integral_type"] = 2
      opts_tmp["Volume_flux_name"] = "IRFlux"
      opts_tmp["Flux_name"] = "IRFlux"
      opts_tmp["face_integral_type"] = 2
      opts_tmp["FaceElementIntegral_name"] = "ESLFFaceIntegral"
      make_input(opts_tmp, fname2)
    end
    MPI.Barrier(MPI.COMM_WORLD)
    mesh7, sbp7, eqn7, opts7 = run_solver(fname2)

    # SBPDiagonalE ES
    if myrank == 0
      opts_tmp = read_input_file(fname)
      opts_tmp["operator_type"] = "SBPDiagonalE"
      opts_tmp["volume_integral_type"] = 2
      opts_tmp["Volume_flux_name"] = "IRFlux"
      opts_tmp["Flux_name"] = "IRSLFFlux"
      make_input(opts_tmp, fname2)
    end
    MPI.Barrier(MPI.COMM_WORLD)
    mesh8, sbp8, eqn8, opts8 = run_solver(fname2)

    # SBPDiagonalE ES shock capturing
    if myrank == 0
      opts_tmp = read_input_file(fname)
      opts_tmp["operator_type"] = "SBPDiagonalE"
      opts_tmp["volume_integral_type"] = 2
      opts_tmp["Volume_flux_name"] = "IRFlux"
      opts_tmp["Flux_name"] = "IRSLFFlux"
      opts_tmp["addShockCapturing"] = true
      opts_tmp["shock_capturing_name"] = "SBPParabolicReduced"
      opts_tmp["shock_sensor_name"] = "SensorVelocity"
      make_input(opts_tmp, fname2)
    end
    MPI.Barrier(MPI.COMM_WORLD)
    mesh9, sbp9, eqn9, opts9 = run_solver(fname2)

    opts4_tmp = copy(opts4)
    test_jac_parallel_inner(mesh4, sbp4, eqn4, opts4)
    test_jac_homotopy(mesh4, sbp4, eqn4, opts4_tmp)
    test_revm_product(mesh4, sbp4, eqn4, opts4)
    test_revq_product(mesh4, sbp4, eqn4, opts4)

    test_jac_parallel_inner(mesh5, sbp5, eqn5, opts5)

    # run test twice to make sure arrays are zeroed out correctly
    test_jac_parallel_inner(mesh5, sbp5, eqn5, opts5)

    test_revm_product(mesh5, sbp5, eqn5, opts5)

    test_jac_parallel_inner(mesh6, sbp6, eqn6, opts6)

    test_jac_parallel_inner(mesh7, sbp7, eqn7, opts7)
    test_jac_parallel_inner(mesh7, sbp7, eqn7, opts7)

    test_revm_product(mesh7, sbp7, eqn7, opts7)

    test_revq_product(mesh7, sbp7, eqn7, opts7)

    test_revm_product2(mesh9, sbp9, eqn9, opts9)
    test_revq_product(mesh9, sbp9, eqn9, opts9)
    test_jac_parallel_inner(mesh9, sbp9, eqn9, opts9)

    # test functional that require parallel communication
    func_test = ["lpsdissipation", "neglpsdissipation"]
    for (name, func_ctor) in EulerEquationMod.FunctionalDict
      func = func_ctor(Complex128, mesh4, sbp4, eqn4, opts4, [1])
      if getParallelData(func) != PARALLEL_DATA_NONE || name in func_test
        if mesh4.myrank == 0
          println("testing functional ", func_ctor)
        end
        test_functional_comm(mesh4, sbp4, eqn4, opts4, func)
        test_functional_deriv_m(mesh4, sbp4, eqn4, opts4, func)

        test_functional_comm(mesh7, sbp7, eqn7, opts7, func)
        test_functional_deriv_q(mesh7, sbp7, eqn7, opts7, func)
        test_functional_deriv_m(mesh7, sbp7, eqn7, opts7, func)
      end
    end

    testEntropyDissFunctional2(mesh8, sbp8, eqn8, opts8)

  end

  return nothing
end

add_func1!(EulerTests, test_jac_parallel_long, [TAG_LONGTEST, TAG_JAC])

#------------------------------------------------------------------------------
# functions that run individual tests

function rand_realpart(dims...)

  a = rand(Complex128, dims...)
  for i=1:length(a)
    a[i] = real(a[i])
  end

  return a
end

function test_jac_parallel_inner(mesh, sbp, eqn, opts; is_prealloc_exact=true, set_prealloc=true)

  # use a spatially varying solution
  icfunc = EulerEquationMod.ICDict["ICExp"]
  icfunc(mesh, sbp, eqn, opts, eqn.q_vec)
  array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)

  # get the correct differentiated flux function (this is needed because the
  # input file set calc_jac_explicit = false
  eqn.flux_func_diff = EulerEquationMod.FluxDict_diff[opts["Flux_name"]]
  eqn.volume_flux_func_diff = EulerEquationMod.FluxDict_diff[opts["Volume_flux_name"]]

  startSolutionExchange(mesh, sbp, eqn, opts)

  opts["calc_jac_explicit"] = false
  pc1, lo1 = NonlinearSolvers.getNewtonPCandLO(mesh, sbp, eqn, opts)
  opts["calc_jac_explicit"] = true
  val_orig = opts["preallocate_jacobian_coloring"]
  if set_prealloc
    opts["preallocate_jacobian_coloring"] = false
  end
  pc2, lo2 = NonlinearSolvers.getNewtonPCandLO(mesh, sbp, eqn, opts)
  opts["preallocate_jacobian_coloring"] = val_orig

  jac1 = getBaseLO(lo1).A
  jac2 = getBaseLO(lo2).A

  assembler = Jacobian._AssembleElementData(getBaseLO(lo2).A, mesh, sbp, eqn, opts)

  # compute jacobian via coloring
  println("\n\nComputing Jacobian via coloring")
  opts["calc_jac_explicit"] = false
  ctx_residual = (evalResidual,)
  NonlinearSolvers.physicsJac(mesh, sbp, eqn, opts, jac1, ctx_residual)

  # compute jacobian explicitly
  println("\n\nComputing Jacobian explicitly")
  opts["calc_jac_explicit"] = true
  evalJacobian(mesh, sbp, eqn, opts, assembler)
  println("finished computing jacobian explicitly")
  assembly_begin(jac1, MAT_FINAL_ASSEMBLY)
  assembly_begin(jac2, MAT_FINAL_ASSEMBLY)

  assembly_end(jac1, MAT_FINAL_ASSEMBLY)
  assembly_end(jac2, MAT_FINAL_ASSEMBLY)

  # multiply against a random vector to make sure the jacobian is
  # the same
  for i=1:1
#    x = rand(PetscScalar, mesh.numDof)
    x = zeros(PetscScalar, mesh.numDof)
    if mesh.myrank == 0
      x[111] = 1
    end
    b1 = zeros(PetscScalar, mesh.numDof)
    b2 = zeros(PetscScalar, mesh.numDof)
    b3 = zeros(PetscScalar, mesh.numDof)

    t = 0.0
    applyLinearOperator(lo1, mesh, sbp, eqn, opts, ctx_residual, t, x, b1)
    applyLinearOperator(lo2, mesh, sbp, eqn, opts, ctx_residual, t, x, b2)
    evaldRdqProduct(mesh, sbp, eqn, opts, x, b3)

    @test isapprox( norm(b1 - b2), 0.0) atol=1e-12
    @test isapprox( norm(b1 - b3), 0.0) atol=1e-12
    println("norm(b1 - b2) = ", norm(b1 - b2))
  end

  A = getBaseLO(lo2).A
  if typeof(A) <: PetscMat
    matinfo = MatGetInfo(A, PETSc2.MAT_LOCAL)
    if is_prealloc_exact
      @test ( matinfo.nz_unneeded )== 0
    else
      @test  matinfo.nz_unneeded  > 0
    end
  end

  free(lo1)
  free(lo2)
  free(pc1)
  free(pc2)

  return nothing
end

import EulerEquationMod: EulerData

function test_jac_homotopy(mesh, sbp, eqn::EulerData{Tsol, Tres}, opts) where {Tsol, Tres}

  # use a spatially varying solution
  icfunc = EulerEquationMod.ICDict["ICExp"]
  icfunc(mesh, sbp, eqn, opts, eqn.q_vec)
  array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)

  # get the correct differentiated flux function (this is needed because the
  # input file set calc_jac_explicit = false
  eqn.flux_func_diff = EulerEquationMod.FluxDict_diff["RoeFlux"]

  res1 = zeros(eqn.res)
  res2 = zeros(eqn.res)
#  EulerEquationMod.calcHomotopyDiss(mesh, sbp, eqn, opts, res1)
  h = 1e-20
  pert = Complex128(0, h)
  eqn.q[1] += pert
  EulerEquationMod.calcHomotopyDiss(mesh, sbp, eqn, opts, res2)
  eqn.q[1] -= pert
#=
  println("diffnorm = ", vecnorm(res1 - res2))
  println("res1 = \n", res1)
  println("res2 = \n", res2)
  println("diff = \n", res1 - res2)
  @assert vecnorm(res1 - res2) < 1e-13
=#
  startSolutionExchange(mesh, sbp, eqn, opts, wait=true)

  hdata = NonlinearSolvers.HomotopyData{Tsol, Float64}(mesh, sbp, eqn, opts,
                                                       evalResidual, evalHomotopy)
  hdata.lambda = 1
  opts["calc_jac_explicit"] = false
  pc1, lo1 = NonlinearSolvers.getHomotopyPCandLO(mesh, sbp, eqn, opts)
  lo1.hdata = hdata

  opts["calc_jac_explicit"] = true
  pc2, lo2 = NonlinearSolvers.getHomotopyPCandLO(mesh, sbp, eqn, opts)
  lo2.hdata = hdata

  jac1 = getBaseLO(lo1).A
  jac2 = getBaseLO(lo2).A

  assembler = Jacobian._AssembleElementData(getBaseLO(lo2).A, mesh, sbp, eqn, opts)

  function _evalHomotopy(mesh, sbp, eqn, opts, t)
    evalHomotopy(mesh, sbp, eqn, opts, eqn.res, t)
  end

  ctx_residual = (_evalHomotopy,)

  opts["calc_jac_explicit"] = false
  NonlinearSolvers.physicsJac(mesh, sbp, eqn, opts, jac1, ctx_residual)

  # compute jacobian explicitly
  opts["calc_jac_explicit"] = true

  evalHomotopyJacobian(mesh, sbp, eqn, opts, assembler, lo2.hdata.lambda)

  assembly_begin(jac1, MAT_FINAL_ASSEMBLY)
  assembly_begin(jac2, MAT_FINAL_ASSEMBLY)

  # multiply against a random vector to make sure the jacobian is
  # the same
  for i=1:10
    x = rand(PetscScalar, mesh.numDof)
    b1 = zeros(PetscScalar, mesh.numDof)
    b2 = zeros(PetscScalar, mesh.numDof)

    t = 0.0
    applyLinearOperator(lo1, mesh, sbp, eqn, opts, ctx_residual, t, x, b1)
    applyLinearOperator(lo2, mesh, sbp, eqn, opts, ctx_residual, t, x, b2)

    @test isapprox( norm(b1 - b2), 0.0) atol=1e-12
  end

  free(lo1)
  free(lo2)
  free(pc1)
  free(pc2)

  return nothing
end

#DEBUGGING
function zeroSharedElements(mesh, sbp, eqn, opts, qvec)

  for i=1:mesh.npeers
    for el in mesh.local_element_lists[i]
      for j=1:mesh.numNodesPerElement
        for k=1:mesh.numDofPerNode
          qvec[mesh.dofs[k, j, el]] = 0
        end
      end
    end
  end

  return nothing
end

function test_revm_product(mesh, sbp, eqn, opts)

  h = 1e-20
  pert = Complex128(0, h)

  icfunc = EulerEquationMod.ICDict["ICExp"]
  icfunc(mesh, sbp, eqn, opts, eqn.q_vec)
  eqn.q_vec .+= 0.01.*rand(size(eqn.q_vec))  # add a little noise, to make jump across
  array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)
                                     # interfaces non-zero

  # fields: dxidx, jac, nrm_bndry, nrm_face, coords_bndry

  res_bar = rand_realpart(mesh.numDof)

  dxidx_dot       = rand_realpart(size(mesh.dxidx))
  jac_dot         = rand_realpart(size(mesh.jac))
  nrm_bndry_dot   = rand_realpart(size(mesh.nrm_bndry))
  nrm_face_dot    = rand_realpart(size(mesh.nrm_face_bar))
  coords_bndry_dot = rand_realpart(size(mesh.coords_bndry))
  nrm_sharedface_dot = Array{Array{Complex128, 3}}(mesh.npeers)
  for i=1:mesh.npeers
    nrm_sharedface_dot[i] = rand_realpart(size(mesh.nrm_sharedface[i]))
    #nrm_sharedface_dot[i] = zeros(Complex128, size(mesh.nrm_sharedface[i]))
  end
  nrm_sharedface_dot[1][1] = 1

  zeroBarArrays(mesh)

  mesh.dxidx        .+= pert*dxidx_dot
  mesh.jac          .+= pert*jac_dot
  mesh.nrm_bndry    .+= pert*nrm_bndry_dot
  mesh.nrm_face     .+= pert*nrm_face_dot
  mesh.coords_bndry .+= pert*coords_bndry_dot
  for i=1:mesh.npeers
    mesh.nrm_sharedface[i] .+= pert*nrm_sharedface_dot[i]
  end

  startSolutionExchange(mesh, sbp, eqn, opts)

  fill!(eqn.res, 0)
  evalResidual(mesh, sbp, eqn, opts)
  array3DTo1D(mesh, sbp, eqn, opts, eqn.res, eqn.res_vec)
  val = sum(imag(eqn.res_vec)/h .* res_bar)

#  val = MPI.Allreduce(val, MPI.SUM, eqn.comm)

  mesh.dxidx        .-= pert*dxidx_dot
  mesh.jac          .-= pert*jac_dot
  mesh.nrm_bndry    .-= pert*nrm_bndry_dot
  mesh.nrm_face     .-= pert*nrm_face_dot
  mesh.coords_bndry .-= pert*coords_bndry_dot
  for i=1:mesh.npeers
    mesh.nrm_sharedface[i] .-= pert*nrm_sharedface_dot[i]
  end

  evalResidual_revm(mesh, sbp, eqn, opts, res_bar)
  val2 = 0
  val2 = sum(mesh.dxidx_bar .* dxidx_dot)              +
         sum(mesh.jac_bar .* jac_dot)                  +
         sum(mesh.nrm_bndry_bar .* nrm_bndry_dot)      +
         sum(mesh.nrm_face_bar .* nrm_face_dot)        +
         sum(mesh.coords_bndry_bar .* coords_bndry_dot)

  for i=1:mesh.npeers
    val2 += sum(mesh.nrm_sharedface_bar[i] .* nrm_sharedface_dot[i])
  end

#  val2 = MPI.Allreduce(val2, MPI.SUM, eqn.comm)

  @test abs(val - val2) < 1e-12

  return nothing
end

"""
  Tests the reverse mode product by perturbing mesh.vert_coords and
  recalculating the metrics from that.  This has the benefit of testing
  the `mesh.remote_metrics_bar`
"""
function test_revm_product2(mesh, sbp, eqn, opts)

  h = 1e-20
  pert = Complex128(0, h)

  icfunc = EulerEquationMod.ICDict["ICExp"]
  icfunc(mesh, sbp, eqn, opts, eqn.q_vec)
  array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)
  eqn.q .+= 0.01.*rand(size(eqn.q))  # add a little noise, to make jump across
                                     # interfaces non-zero

  vert_coords_dot = rand_realpart(size(mesh.vert_coords))

  res_bar = rand_realpart(mesh.numDof)

  # complex step
  mesh.vert_coords .+= pert*vert_coords_dot
  recalcCoordinatesAndMetrics(mesh, sbp, opts)
  println("after recalculating metrics:")
  println("has nan dxidx = ", any(isnan.(mesh.dxidx)))
  println("has nan jac = ", any(isnan.(mesh.jac)))
  for peer=1:mesh.npeers
    println("  peer ", peer)
    obj = mesh.remote_metrics[peer]
    println("  has nan dxidx = ", any(isnan.(obj.dxidx)))
    println("  has nan jac = ", any(isnan.(obj.jac)))
  end

  startSolutionExchange(mesh, sbp, eqn, opts)
  evalResidual(mesh, sbp, eqn, opts)
  println("after evalResidual, isnan eqn.res = ", any(isnan.(eqn.res)))
  array3DTo1D(mesh, sbp, eqn, opts, eqn.res, eqn.res_vec)
  val1 = sum(imag(eqn.res_vec)/h .* res_bar)

  # undo perturbation
  for i=1:length(mesh.vert_coords)
    mesh.vert_coords[i] = real(mesh.vert_coords[i])
  end
  recalcCoordinatesAndMetrics(mesh, sbp, opts)

  # reverse mode

  zeroBarArrays(mesh)
  evalResidual_revm(mesh, sbp, eqn, opts, res_bar)
  getAllCoordinatesAndMetrics_rev(mesh, sbp, opts)

  val2 = sum(mesh.vert_coords_bar .* vert_coords_dot)

  val1 = MPI.Allreduce(val1, MPI.SUM, eqn.comm)
  val2 = MPI.Allreduce(val2, MPI.SUM, eqn.comm)
  println("val1 = ", val1)
  println("val2 = ", val2)

  @test abs(val1 - val2) < max( min(abs(val1), abs(val2))*1e-12, 1e-12)

  return nothing
end

function test_revq_product(mesh, sbp, eqn, opts)

  srand(1234)
  h = 1e-20
  pert = Complex128(0, h)

  icfunc = EulerEquationMod.ICDict["ICExp"]
  icfunc(mesh, sbp, eqn, opts, eqn.q_vec)
  eqn.q_vec .+= 0.01.*rand(size(eqn.q_vec))  # add a little noise, to make jump across
                                     # interfaces non-zero

  array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)

  # fields: dxidx, jac, nrm_bndry, nrm_face, coords_bndry

  res_vec_bar = rand_realpart(mesh.numDof)
  q_vec_dot = rand_realpart(mesh.numDof)
  q_vec_bar = zeros(Complex128, mesh.numDof)

  fill!(eqn.res, 0)
  eqn.q_vec .+= pert*q_vec_dot
  array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)

  if opts["parallel_type"] != 1
    startSolutionExchange(mesh, sbp, eqn, opts)
  end

  evalResidual(mesh, sbp, eqn, opts)
  eqn.q_vec .-= pert*q_vec_dot
  array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)

  array3DTo1D(mesh, sbp, eqn, opts, eqn.res, eqn.res_vec)
  val = sum(res_vec_bar .* imag(eqn.res_vec)/h)
  val = MPI.Allreduce(val, MPI.SUM, mesh.comm)

  evalResidual_revq(mesh, sbp, eqn, opts, res_vec_bar, q_vec_bar)
  val2 = sum(q_vec_bar .* q_vec_dot)
  val2 = MPI.Allreduce(val2, MPI.SUM, mesh.comm)

  @test abs(val - val2) < 1e-12

  # test accumulation
  q_vec_bar_orig = copy(q_vec_bar)
  evalResidual_revq(mesh, sbp, eqn, opts, res_vec_bar, q_vec_bar, zero_output=false)

  @test maximum(abs.(2.*q_vec_bar_orig - q_vec_bar)) < 1e-13

  # test start_comm = false
  fill!(q_vec_bar, 0)
  evalResidual_revq(mesh, sbp, eqn, opts, res_vec_bar, q_vec_bar, start_comm=false)
  val2 = sum(q_vec_bar .* q_vec_dot)
  val2 = MPI.Allreduce(val2, MPI.SUM, mesh.comm)

  @test abs(val - val2) < 1e-12

  return nothing
end

"""
  Test functional with start_comm=true and false
"""
function test_functional_comm(mesh, sbp, eqn, opts, func)

  h = 1e-60
  pert = Complex128(0, h)

  # use a spatially varying solution
  icfunc = EulerEquationMod.ICDict["ICRho1E2U3"]
  icfunc(mesh, sbp, eqn, opts, eqn.q_vec)
  eqn.q_vec .+= 0.1*rand(length(eqn.q_vec))
  array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)

  f = evalFunctional(mesh, sbp, eqn, opts, func)
  f2 = evalFunctional(mesh, sbp, eqn, opts, func, start_comm=false)

  @test abs(f - f2) < 1e-13

  return nothing
end

"""
  Test functional derivative wrt q against complex step
"""
function test_functional_deriv_q(mesh, sbp, eqn, opts, func)

  h = 1e-60
  pert = Complex128(0, h)

  # use a spatially varying solution
  icfunc = EulerEquationMod.ICDict["ICRho1E2U3"]
  icfunc(mesh, sbp, eqn, opts, eqn.q_vec)
  eqn.q_vec .+= 0.1*rand(length(eqn.q_vec))
  array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)

  q_dot = rand(size(eqn.q))
  q_bar = zeros(eqn.q)

  eqn.q_vec .+= pert*vec(q_dot)
  f = evalFunctional(mesh, sbp, eqn, opts, func)
  val = imag(f)/h
  # evalFunctional does the Allreduce, don't duplicate
  eqn.q_vec .-= pert*vec(q_dot)

  evalFunctionalDeriv_q(mesh, sbp, eqn, opts, func, q_bar)
  val2 = sum(q_bar .* q_dot)
  val2 = MPI.Allreduce(val2, MPI.SUM, eqn.comm)

   @test abs(val - val2) < 1e-12

  # test start_comm = false
  q_bar_orig = copy(q_bar)
  evalFunctionalDeriv_q(mesh, sbp, eqn, opts, func, q_bar, start_comm=false)

  @test maximum(abs.(q_bar - q_bar_orig)) < 1e-13

  return nothing
end

function test_functional_deriv_m(mesh, sbp, eqn, opts, func)

  h = 1e-20
  pert = Complex128(0, h)

  icfunc = EulerEquationMod.ICDict["ICExp"]
  icfunc(mesh, sbp, eqn, opts, eqn.q_vec)
  eqn.q_vec .+= 0.01.*rand(size(eqn.q_vec))  # add a little noise, to make jump across
  array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)
                                     # interfaces non-zero

  # fields: dxidx, jac, nrm_bndry, nrm_face, coords_bndry

  res_bar = rand_realpart(mesh.numDof)

  dxidx_dot       = rand_realpart(size(mesh.dxidx))
  jac_dot         = rand_realpart(size(mesh.jac))
  nrm_bndry_dot   = rand_realpart(size(mesh.nrm_bndry))
  nrm_face_dot    = rand_realpart(size(mesh.nrm_face_bar))
  coords_bndry_dot = rand_realpart(size(mesh.coords_bndry))
  nrm_sharedface_dot = Array{Array{Complex128, 3}}(mesh.npeers)
  for i=1:mesh.npeers
    nrm_sharedface_dot[i] = rand_realpart(size(mesh.nrm_sharedface[i]))
  end

  zeroBarArrays(mesh)

  mesh.dxidx        .+= pert*dxidx_dot
  mesh.jac          .+= pert*jac_dot
  mesh.nrm_bndry    .+= pert*nrm_bndry_dot
  mesh.nrm_face     .+= pert*nrm_face_dot
  mesh.coords_bndry .+= pert*coords_bndry_dot
  for i=1:mesh.npeers
    mesh.nrm_sharedface[i] .+= pert*nrm_sharedface_dot[i]
  end

  f = evalFunctional(mesh, sbp, eqn, opts, func)
  val = imag(f)/h

  mesh.dxidx        .-= pert*dxidx_dot
  mesh.jac          .-= pert*jac_dot
  mesh.nrm_bndry    .-= pert*nrm_bndry_dot
  mesh.nrm_face     .-= pert*nrm_face_dot
  mesh.coords_bndry .-= pert*coords_bndry_dot
  for i=1:mesh.npeers
    mesh.nrm_sharedface[i] .-= pert*nrm_sharedface_dot[i]
  end

  evalFunctionalDeriv_m(mesh, sbp, eqn, opts, func)

  val2 = sum(mesh.dxidx_bar .* dxidx_dot)              +
         sum(mesh.jac_bar .* jac_dot)                  +
         sum(mesh.nrm_bndry_bar .* nrm_bndry_dot)      +
         sum(mesh.nrm_face_bar .* nrm_face_dot)        +
         sum(mesh.coords_bndry_bar .* coords_bndry_dot)

  for i=1:mesh.npeers
    val2 += sum(mesh.nrm_sharedface_bar[i] .* nrm_sharedface_dot[i])
  end

  # the functional does an allreduce, so do it here too
  val2 = MPI.Allreduce(val2, MPI.SUM, eqn.comm)

  @test abs(val - val2) < 1e-12

  return nothing
end

function testEntropyDissFunctional2(mesh, sbp, eqn, opts)

  # test that the two methods of computing the entropy dissipation give the
  # same answer
  # This depends on the eqn object using the IRSLF flux with diagonal E
  # operators

  eqn.q .+= 0.01*rand(size(eqn.q))
  startSolutionExchange(mesh, sbp, eqn, opts, wait=true)

  func1 = createFunctional(mesh, sbp, eqn, opts, "entropydissipation", [1])
  func2 = createFunctional(mesh, sbp, eqn, opts, "entropydissipation2", [2])

  val1 = evalFunctional(mesh, sbp, eqn, opts, func1)
  val2 = evalFunctional(mesh, sbp, eqn, opts, func2)

  @test abs(val1 - val2) < 1e-13

  setParallelData(eqn.shared_data, PARALLEL_DATA_FACE)
  val1 = evalFunctional(mesh, sbp, eqn, opts, func1)
  val2 = evalFunctional(mesh, sbp, eqn, opts, func2)

  @test abs(val1 - val2) < 1e-13

  setParallelData(eqn.shared_data, PARALLEL_DATA_ELEMENT)

  return nothing
end


text:
"""
Base type for objects wrapping measurement data.
"""
abstract AbstractData

"""
    data(s::Symbol,args...)

Create data objects.

Currently implemented variants:
---

"""
data(s::Symbol,args...) = data(Val{s},args...)

"""
    datatypename(d)
Return the name of the data type (used by `show()`).
"""
function datatypename(d::AbstractData) end

"""
    numvalues(d)
Return the number of data values.
"""
function numvalues(d::AbstractData) end

"""
    numvars(d)
Return the number of variables.
"""
function numvars(d::AbstractData) end

"""
    eltype(d)
Return the type of the data values.
"""
function eltype(d::AbstractData) end

"""
    generate!(d)
Generate a set of data values. The function should return its argument d.
"""
function generate!(d::AbstractData) end

"""
    dataindex(d)
Retrieve the data index.
"""
function dataindex(d::AbstractData) end

"""
    datavalues(d)
Retrieve the data values.
"""
function datavalues(d::AbstractData) end

function show(io::IO,d::AbstractData)
    nvals = (n = numvalues(d) ; n>1?string(n," values"):string(n," value"))
    nvars = (n = numvars(d) ; n>1?string(n," variables"):string(n," variable"))
    println(io,"$(datatypename(d)) with $(nvars) and $(nvals)")
    println(io," index: ",dataindex(d))
    println(io," values: ",datavalues(d))
    nothing
end

ext: jl


text:
module FluorophoreMarkers
using Unitful

export FluorophoreMarker, Fluorophore

# TODO: Instead of an enum, this should be a symbol, which is more natural in julia. <18-05-21, kunzaatko> #
# @enum Fluorophore {{{
@enum Fluorophore begin
    # [P4-3]
    # BFP is a basic (constitutively fluorescent) blue fluorescent protein published in 1996, derived from Aequorea victoria.
    BFP
    # EBFP is a basic (constitutively fluorescent) blue fluorescent protein published in 1998, derived from Aequorea victoria.
    EBFP
    # [Y66W]
    # CFP is a basic (constitutively fluorescent) cyan fluorescent protein published in 1994, derived from Aequorea victoria.
    CFP
    # [enhanced GFP, GFPmut1]
    # EGFP is a basic (constitutively fluorescent) green fluorescent protein published in 1996, derived from Aequorea victoria. It is reported to be a rapidly-maturing weak dimer with moderate acid sensitivity.
    EGFP
    # [wtGFP, GFP, gfp10, Green Fluorescent Protein]
    # avGFP is a basic (constitutively fluorescent) long stokes shift fluorescent protein published in 1992, derived from Aequorea victoria. It has low acid sensitivity.
    avGFP
    # dTomato is a basic (constitutively fluorescent) orange fluorescent protein published in 2004, derived from Discosoma sp.. It is reported to be a somewhat slowly-maturing dimer with low acid sensitivity.     dTomato
    dTomato
    # tdTomato is a basic (constitutively fluorescent) orange fluorescent protein published in 2004, derived from Discosoma sp.. It is reported to be a somewhat slowly-maturing tandem dimer with low acid sensitivity.
    tdTomato
    # [drFP583, discRFP, RFP, DsRed1]
    # DsRed is a basic (constitutively fluorescent) red fluorescent protein published in 1999, derived from Discosoma sp..
    DsRed
    # TagRFP is a basic (constitutively fluorescent) red fluorescent protein published in 2007, derived from Entacmaea quadricolor. It has very low acid sensitivity.
    TagRFP
    # [mRFP]
    # mRFP1 is a basic (constitutively fluorescent) red fluorescent protein published in 2002, derived from Discosoma sp.. It is reported to be a somewhat slowly-maturing monomer with low acid sensitivity.
    mRFP1
    # mCherry is a basic (constitutively fluorescent) red fluorescent protein published in 2004, derived from Discosoma sp.. It is reported to be a very rapidly-maturing monomer with low acid sensitivity.
    mCherry
    # TurboRFP is a basic (constitutively fluorescent) orange fluorescent protein published in 2007, derived from Entacmaea quadricolor. It has low acid sensitivity.
    TurboRFP
    # mRuby2 is a basic (constitutively fluorescent) red fluorescent protein published in 2012, derived from Entacmaea quadricolor. It has moderate acid sensitivity.
    mRuby2
    # [mNG]
    # mNeonGreen is a basic (constitutively fluorescent) green/yellow fluorescent protein published in 2013, derived from Branchiostoma lanceolatum. It is reported to be a very rapidly-maturing monomer with moderate acid sensitivity.
    mNeonGreen
    # CyOFP1 is a basic (constitutively fluorescent) long stokes shift fluorescent protein published in 2016, derived from Entacmaea quadricolor. It is reported to be a very rapidly-maturing monomer with moderate acid sensitivity.
    CyOFP1
    # mScarlet-I is a basic (constitutively fluorescent) red fluorescent protein published in 2016, derived from synthetic construct. It is reported to be a rapidly-maturing monomer with moderate acid sensitivity.
    mScarlet_I # FCHO2
    # mScarlet is a basic (constitutively fluorescent) red fluorescent protein published in 2016, derived from synthetic construct. It has moderate acid sensitivity.
    mScarlet
end
# }}}

# struct FluorophoreMarker {{{
struct FluorophoreMarker
    fluorophore::Fluorophore
    exλ::Unitful.Length
    emλ::Unitful.Length
end
# }}}

# FluorophoreMarker(fluorophore::Fluorophore) {{{
function FluorophoreMarker(fluorophore::Fluorophore)
    if fluorophore == BFP::Fluorophore
        exλ = 381u"nm"
        emλ = 440u"nm"
        return FluorophoreMarker(fluorophore, exλ, emλ)
    end
    if fluorophore == EBFP::Fluorophore
        exλ = 380u"nm"
        emλ = 440u"nm"
        return FluorophoreMarker(fluorophore, exλ, emλ)
    end
    if fluorophore == CFP::Fluorophore
        exλ = 456u"nm"
        emλ = 475u"nm"
        return FluorophoreMarker(fluorophore, exλ, emλ)
    end
    if fluorophore == EGFP::Fluorophore
        exλ = 488u"nm"
        emλ = 507u"nm"
        return FluorophoreMarker(fluorophore, exλ, emλ)
    end
    if fluorophore == avGFP::Fluorophore
        exλ = 488u"nm"
        emλ = 509u"nm"
        return FluorophoreMarker(fluorophore, exλ, emλ)
    end
    if fluorophore == dTomato::Fluorophore
        exλ = 554u"nm"
        emλ = 581u"nm"
        return FluorophoreMarker(fluorophore, exλ, emλ)
    end
    if fluorophore == tdTomato::Fluorophore
        exλ = 554u"nm"
        emλ = 581u"nm"
        return FluorophoreMarker(fluorophore, exλ, emλ)
    end
    if fluorophore == DsRed::Fluorophore
        exλ = 558u"nm"
        emλ = 583u"nm"
        return FluorophoreMarker(fluorophore, exλ, emλ)
    end
    if fluorophore == TagRFP::Fluorophore
        exλ = 555u"nm"
        emλ = 584u"nm"
        return FluorophoreMarker(fluorophore, exλ, emλ)
    end
    if fluorophore == mRFP1::Fluorophore
        exλ = 584u"nm"
        emλ = 607u"nm"
        return FluorophoreMarker(fluorophore, exλ, emλ)
    end
    if fluorophore == mCherry::Fluorophore
        exλ = 587u"nm"
        emλ = 610u"nm"
        return FluorophoreMarker(fluorophore, exλ, emλ)
    end
    if fluorophore == TurboRFP::Fluorophore
        exλ = 553u"nm"
        emλ = 574u"nm"
        return FluorophoreMarker(fluorophore, exλ, emλ)
    end
    if fluorophore == mRuby2::Fluorophore
        exλ = 559u"nm"
        emλ = 600u"nm"
        return FluorophoreMarker(fluorophore, exλ, emλ)
    end
    if fluorophore == mNeonGreen::Fluorophore
        exλ = 506u"nm"
        emλ = 517u"nm"
        return FluorophoreMarker(fluorophore, exλ, emλ)
    end
    if fluorophore == CyOFP1::Fluorophore
        exλ = 497u"nm"
        emλ = 589u"nm"
        return FluorophoreMarker(fluorophore, exλ, emλ)
    end
    if fluorophore == mScarlet_I::Fluorophore
        exλ = 569u"nm"
        emλ = 593u"nm"
        return FluorophoreMarker(fluorophore, exλ, emλ)
    end
    if fluorophore == mScarlet::Fluorophore
        exλ = 569u"nm"
        emλ = 594u"nm"
        return FluorophoreMarker(fluorophore, exλ, emλ)
    end
end
# }}}

end


text:
import CSV
using DataFrames
include("plot_utils.jl")

function plot_accuracy_vs_coefficient(dataframe, output_path::String; lower_bound = 1e-7)
    general_setup()
    plot(
        dataframe[!, "coeff"],
        dataframe[!, "error_subg"],
        color = 2,
        line = (3, :solid),
        label = "Subgradient",
    )

    plot!(
        dataframe[!, "coeff"],
        dataframe[!, "error_bundle"],
        color = 1,
        line = (3, :solid),
        label = "Proximal bundle",
    )
    yaxis!("Objective gap")
    xaxis!("Regularizer coefficient")
    savefig(output_path)
end

function plot_log_accuracy_vs_coefficient(
    dataframe,
    output_path::String;
    lower_bound = 1e-7,
)
    general_setup()
    plot(
        dataframe[!, "coeff"],
        dataframe[!, "error_subg"],
        yaxis = (:log10, [lower_bound, :auto]),
        color = 2,
        line = (3, :solid),
        label = "Subgradient",
        legend = :bottomright,
    )

    plot!(
        dataframe[!, "coeff"],
        dataframe[!, "error_bundle"],
        color = 1,
        line = (3, :solid),
        label = "Proximal bundle",
    )
    yaxis!("Objective gap")
    xaxis!("Regularizer coefficient")
    savefig(output_path)
end

function main()
    results = DataFrame(CSV.File("results/svm/results_svm.csv"))
    instances = unique(results[!, "instance"])
    coefficients = unique(results[!, "coeff"])

    for name in instances
        plot_accuracy_vs_coefficient(
            results[(results.instance.==name).&(results.coeff.<10.0), :],
            "results/svm/svm_accuracy_vs_coefficient_" * string(name) * ".pdf",
        )
        plot_log_accuracy_vs_coefficient(
            results[(results.instance.==name).&(results.coeff.<10.0), :],
            "results/svm/svm_log_accuracy_vs_coefficient_" * string(name) * ".pdf",
        )
    end
end

main()


text:
using MieScattering:
    mie_scattering, mie_scattering_n, mie_ab, mie_ab_n, mie_pi_tau, mie_pi_tau_n, mie_S1_S2, mie_S1_S2_n, scattering_function, scattering_function_n
using Test

@testset "scattering values" begin
    @test isapprox(mie_scattering(1.0, 1.0)[2], 0.0, atol = 1e-6)
    @test isapprox(abs(mie_scattering(4.0, 1.0)[2]), 6.062172, atol = 1e-6)
end

@testset "self-consistency" begin

    m = 1.33 + 0.01im
    x = 3

    an, bn = mie_ab(m, x, 3)

    @test isapprox(an[1], mie_ab_n(m, x, 1)[1], atol=1e-6)
    @test isapprox(an[2], mie_ab_n(m, x, 2)[1], atol=1e-6)
    @test isapprox(an[3], mie_ab_n(m, x, 3)[1], atol=1e-6)
    @test isapprox(bn[1], mie_ab_n(m, x, 1)[2], atol=1e-6)
    @test isapprox(bn[2], mie_ab_n(m, x, 2)[2], atol=1e-6)
    @test isapprox(bn[3], mie_ab_n(m, x, 3)[2], atol=1e-6)

    m = 2 + 1im
    x = 5

    an, bn = mie_ab(m, x, 3)
    @test isapprox(an[1], mie_ab_n(m, x, 1)[1])
    @test isapprox(an[2], mie_ab_n(m, x, 2)[1])
    @test isapprox(an[3], mie_ab_n(m, x, 3)[1])
    @test isapprox(bn[1], mie_ab_n(m, x, 1)[2])
    @test isapprox(bn[2], mie_ab_n(m, x, 2)[2])
    @test isapprox(bn[3], mie_ab_n(m, x, 3)[2])

    # error
    # pin, taun = mie_pi_tau(mu, 1)

    mu = cos(0)
    pn, tn = mie_pi_tau(mu, 3)
    @test isapprox(pn[1], mie_pi_tau_n(mu, 1)[1])
    @test isapprox(pn[2], mie_pi_tau_n(mu, 2)[1])
    @test isapprox(pn[3], mie_pi_tau_n(mu, 3)[1])
    @test isapprox(tn[1], mie_pi_tau_n(mu, 1)[2])
    @test isapprox(tn[2], mie_pi_tau_n(mu, 2)[2])
    @test isapprox(tn[3], mie_pi_tau_n(mu, 3)[2])

    mu = cos(π / 3)
    pn, tn = mie_pi_tau(mu, 3)
    @test isapprox(pn[1], mie_pi_tau_n(mu, 1)[1])
    @test isapprox(pn[2], mie_pi_tau_n(mu, 2)[1])
    @test isapprox(pn[3], mie_pi_tau_n(mu, 3)[1])
    @test isapprox(tn[1], mie_pi_tau_n(mu, 1)[2])
    @test isapprox(tn[2], mie_pi_tau_n(mu, 2)[2])
    @test isapprox(tn[3], mie_pi_tau_n(mu, 3)[2])

    mu = cos(π / 2)
    pn, tn = mie_pi_tau(mu, 3)
    @test isapprox(pn[1], mie_pi_tau_n(mu, 1)[1])
    @test isapprox(pn[2], mie_pi_tau_n(mu, 2)[1])
    @test isapprox(pn[3], mie_pi_tau_n(mu, 3)[1])
    @test isapprox(tn[1], mie_pi_tau_n(mu, 1)[2])
    @test isapprox(tn[2], mie_pi_tau_n(mu, 2)[2])
    @test isapprox(tn[3], mie_pi_tau_n(mu, 3)[2])

    m = 1.33 + 0.01im
    x = 3

    qext1, qsca1, _, _ = mie_scattering(m, x)
    qext2, qsca2, _, _ = mie_scattering_n(m, x)

    @test isapprox(qext1, qext2)
    @test isapprox(qsca1, qsca2)

    m = 2 + 1im
    x = 5

    qext1, qsca1, _, _ = mie_scattering(m, x)
    qext2, qsca2, _, _ = mie_scattering_n(m, x)

    @test isapprox(qext1, qext2)
    @test isapprox(qsca1, qsca2)

    m = 1.33 + 0.01im
    x = 3

    mu = cos(0)
    s1s2 = mie_S1_S2(m, x, mu, 2)
    s1s2_n = mie_S1_S2_n(m, x, mu, 1) .+ mie_S1_S2_n(m, x, mu, 2)
    @test all(isapprox.(s1s2, s1s2_n))

    mu = cos(π / 3)
    s1s2 = mie_S1_S2(m, x, mu, 2)
    s1s2_n = mie_S1_S2_n(m, x, mu, 1) .+ mie_S1_S2_n(m, x, mu, 2)
    @test all(isapprox.(s1s2, s1s2_n))

    mu = cos(π / 2)
    s1s2 = mie_S1_S2(m, x, mu, 2)
    s1s2_n = mie_S1_S2_n(m, x, mu, 1) .+ mie_S1_S2_n(m, x, mu, 2)
    @test all(isapprox.(s1s2, s1s2_n))

    m = 2 + 1im
    x = 5

    mu = cos(0)
    s1s2 = mie_S1_S2(m, x, mu, 2)
    s1s2_n = mie_S1_S2_n(m, x, mu, 1) .+ mie_S1_S2_n(m, x, mu, 2)
    @test all(isapprox.(s1s2, s1s2_n))

    mu = cos(π / 3)
    s1s2 = mie_S1_S2(m, x, mu, 2)
    s1s2_n = mie_S1_S2_n(m, x, mu, 1) .+ mie_S1_S2_n(m, x, mu, 2)
    @test all(isapprox.(s1s2, s1s2_n))

    mu = cos(π / 2)
    s1s2 = mie_S1_S2(m, x, mu, 2)
    s1s2_n = mie_S1_S2_n(m, x, mu, 1) .+ mie_S1_S2_n(m, x, mu, 2)
    @test all(isapprox.(s1s2, s1s2_n))

    m = 1.33 + 0.01im
    x = 3

    theta = 0
    r1 = scattering_function(theta, m, x)
    r2 = scattering_function_n(theta, m, x)
    @test all(isapprox.(r1,r2))

    theta = π/3
    r1 = scattering_function(theta, m, x)
    r2 = scattering_function_n(theta, m, x)
    @test all(isapprox.(r1,r2))

    theta = π/2
    r1 = scattering_function(theta, m, x)
    r2 = scattering_function_n(theta, m, x)
    @test all(isapprox.(r1,r2))

    m = 2 + 1im
    x = 5

    theta = 0
    r1 = scattering_function(theta, m, x)
    r2 = scattering_function_n(theta, m, x)
    @test all(isapprox.(r1,r2))

    theta = π/3
    r1 = scattering_function(theta, m, x)
    r2 = scattering_function_n(theta, m, x)
    @test all(isapprox.(r1,r2))

    theta = π/2
    r1 = scattering_function(theta, m, x)
    r2 = scattering_function_n(theta, m, x)
    @test all(isapprox.(r1,r2))

end

max_stars_repo_name: fekad/jlmie

max_issues_repo_name: fekad/jlmie

max_forks_repo_name: fekad/jlmie


text:
"""
    niggli_reduce(lattice::AbstractMatrix, symprec=1e-5)
    niggli_reduce(cell::Cell, symprec=1e-5)

Apply Niggli reduction to input basis vectors `lattice`.
"""
function niggli_reduce(lattice::AbstractMatrix, symprec = 1e-5)
    clattice = convert(Matrix{Cdouble}, transpose(lattice))
    exitcode = ccall(
        (:spg_niggli_reduce, libsymspg),
        Cint,
        (Ptr{Cdouble}, Cdouble),
        clattice,
        symprec,
    )
    if exitcode == 0
        throw(SpglibError("Niggli reduce failed!"))
    end
    return transpose(clattice)
end
function niggli_reduce(cell::Cell, symprec = 1e-5)
    lattice = cell.lattice
    clattice = niggli_reduce(lattice, symprec)
    return Cell(clattice, cell.positions, cell.types, cell.magmoms)
end

"""
    delaunay_reduce(lattice::AbstractMatrix, symprec=1e-5)
    delaunay_reduce(cell::Cell, symprec=1e-5)

Apply Delaunay reduction to input basis vectors `lattice`.
"""
function delaunay_reduce(lattice::AbstractMatrix, symprec = 1e-5)
    clattice = convert(Matrix{Cdouble}, transpose(lattice))
    exitcode = ccall(
        (:spg_delaunay_reduce, libsymspg),
        Cint,
        (Ptr{Cdouble}, Cdouble),
        clattice,
        symprec,
    )
    if exitcode == 0
        throw(SpglibError("Delaunay reduce failed!"))
    end
    return transpose(clattice)
end
function delaunay_reduce(cell::Cell, symprec = 1e-5)
    lattice = cell.lattice
    clattice = delaunay_reduce(lattice, symprec)
    return Cell(clattice, cell.positions, cell.types, cell.magmoms)
end


text:
using Plots
gr()

Random.seed!(2)

colors = [RGB(.22,.596,.149), RGB(.8,.361,.361), RGB(.702, .322, .8)]
border = [RGB(.133, .541, .133), RGB(.8, .2, .2), RGB(.584, .345, .698)]

y = rand(50)
for i in eachindex(y)
    if rand() < .8
        y[i] = 0.0
    end
end

scatter(y, c = colors, ms = 10, markerstrokewidth = 5, markerstrokecolor = border,
    legend=false, grid=false, axis=false)

savefig(joinpath(@__DIR__(), "logo.png"))


text:
task = @async begin
    condition = Base.AsyncCondition()
    output::Ref{UInt32} = C_NULL
    joinhandle = ccall((:multithreaded, "libccall_with_threads"), Ptr{Cvoid}, (Ref{UInt32}, Ptr{Cvoid}), output, condition.handle)
    wait(condition)
    ccall((:drop_handle, "libccall_with_threads"), Cvoid, (Ptr{Cvoid},), joinhandle)

    output[]
end

task2 = @async begin
    while !istaskdone(task)
        println("Still running")
        sleep(0.1)
    end

    @assert fetch(task) == 127 "Wrong result"
end

wait(task)
wait(task2)

max_stars_repo_name: marhkb/jlrs

max_issues_repo_name: marhkb/jlrs

max_forks_repo_name: marhkb/jlrs


text:
"""
    FENodeSetModule

Module for the finite element node set.
"""
module FENodeSetModule

__precompile__(true)

using ..FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import Base.count

"""
Finite element node set type.
"""
mutable struct FENodeSet
	"""
	Array of node locations. Indexed with the node number.
	The location of node `j` is given by `xyz[j,:]`.
	Clearly, the nodes needs to be numbered between `1` and `size(xyz, 1)`.
	The constructor makes a *copy* of the input `xyz` array for safety.
	"""
	xyz::Array{FFlt, 2}

	function FENodeSet(xyz::FFltMat)
		self = new(deepcopy(xyz)); # Need to make a COPY of the input array!
		return self
	end
end

"""
    spacedim(self::FENodeSet)

Number of dimensions of the space in which the node lives, 1, 2, or 3.
"""
spacedim(self::FENodeSet) =size(self.xyz,2)

"""
    xyz3(self::FENodeSet)

Get the  3-D coordinate that define the location  of the node.
Even if the nodes  were specified in  lower dimension (1-D, 2-D)
this function returns  a 3-D coordinate  by padding with zeros.
"""
function xyz3(self::FENodeSet)
    if (size(self.xyz,2)==1)
        val = [self.xyz zeros(size(self.xyz,1),2)];
    elseif (size(self.xyz,2)==2)
        val = [self.xyz zeros(size(self.xyz,1),1)];
    else
        val = deepcopy(self.xyz);
    end
end

"""
    count(self::FENodeSet)

Get the number of finite element nodes in the node set.
"""
function count(self::FENodeSet)
    return size(self.xyz,1)::FInt
end

end


text:
# Hacky script to download the latest version of all packages registered in the
# General registry for testing the parser.
#
# This uses internal Pkg APIs and seems to work on Julia 1.7

using Pkg
using Downloads

registry = only(filter(r->r.name == "General", Pkg.Registry.reachable_registries()))

packages = []

for (uuid,pkg) in registry
    versions = collect(Pkg.Registry.registry_info(pkg).version_info)
    latest_ver, ver_info = last(sort(versions, by=first))
    if ver_info.yanked
        continue
    end

    push!(packages, (; uuid, pkg.name, version=latest_ver, ver_info.git_tree_sha1))

end

server = Pkg.pkg_server()
output_dir = "pkgs"
mkpath(output_dir)

asyncmap(packages, ntasks=5) do pkg
    url = "$server/package/$(pkg.uuid)/$(pkg.git_tree_sha1)"
    outfile_path = joinpath(output_dir, "$(pkg.name)_$(pkg.version).tgz")
    if isfile(outfile_path)
        @info "Skipping package" pkg
        return outfile_path
    else
        @info "Download package" url outfile_path
        for i=1:5
            try
                Downloads.download(url, outfile_path)
                break
            catch
                @error "Error downloading" pkg exception=current_exceptions()
            end
            sleep(i)
        end
    end
end


text:
using ITensors
using ITensorTDVP
using Random
using Test

@testset "DMRG" begin
  N = 10
  cutoff = 1e-12

  s = siteinds("S=1/2", N)

  os = OpSum()
  for j in 1:(N - 1)
    os += 0.5, "S+", j, "S-", j + 1
    os += 0.5, "S-", j, "S+", j + 1
    os += "Sz", j, "Sz", j + 1
  end

  H = MPO(os, s)

  psi = randomMPS(s; linkdims=2)

  nsweeps = 10
  nsite = 2
  maxdim = [10, 20, 40, 100]
  sweeps = Sweeps(nsweeps) # number of sweeps is 5
  maxdim!(sweeps, 10, 20, 40, 100) # gradually increase states kept
  cutoff!(sweeps, cutoff)
  psi = ITensorTDVP.dmrg(
    H, psi; nsweeps, maxdim, cutoff, nsite, solver_krylovdim=3, solver_maxiter=1
  )

  e2, psi2 = dmrg(H, psi, sweeps; normalize=false, outputlevel=0)

  @test inner(psi', H, psi) ≈ inner(psi2', H, psi2)
end

nothing


text:
using ITensors,
      Test

@testset "Physics SiteSets" begin

  N = 10

  @testset "Spin Half Site" begin
    s = SpinHalfSite(3)
    @test hastags(ind(s),"Site")
    @test hastags(ind(s),"S=1/2")
    @test hastags(ind(s),"n=3")

    @test val(state(s,2)) == 2
  end

  @testset "Spin Half SiteSet" begin
    sites = spinHalfSites(N)

    Sz_2 = op(sites,"Sz",2)
    @test Sz_2[1,1] ≈ 0.5
    @test Sz_2[2,2] ≈ -0.5

    s3 = sites[3]
    Sp_3 = op(sites,"S+",3)
    @test Sp_3[s3'(1),s3(2)] ≈ 1.0
  end

  @testset "Spin One Site" begin
    s = SpinOneSite(4)
    @test hastags(ind(s),"Site")
    @test hastags(ind(s),"S=1")
    @test hastags(ind(s),"n=4")

    @test val(state(s,3)) == 3
  end

  @testset "Spin One SiteSet" begin
    sites = spinOneSites(N)

    Sz_2 = op(sites,"Sz",2)
    @test Sz_2[1,1] ≈ 1.0
    @test Sz_2[2,2] ≈ 0.0
    @test Sz_2[3,3] ≈ -1.0

    s3 = sites[3]
    Sp_3 = op(sites,"S+",3)
    @test Sp_3[s3'(1),s3(2)] ≈ sqrt(2.)
    @test Sp_3[s3'(2),s3(3)] ≈ sqrt(2.)
  end

  @testset "Electron Site" begin
    s = ElectronSite(5)
    @test hastags(ind(s),"Site")
    @test hastags(ind(s),"Electron")
    @test hastags(ind(s),"n=5")

    @test val(state(s,4)) == 4
  end

  @testset "Spin One SiteSet" begin
    sites = electronSites(N)

    Nup_2 = op(sites,"Nup",2)
    @test Nup_2[2,2] ≈ 1.0
    @test Nup_2[4,4] ≈ 1.0
  end

  @testset "tJ Site" begin
    s = tJSite(5)
    @test hastags(ind(s),"Site")
    @test hastags(ind(s),"tJ")
    @test hastags(ind(s),"n=5")

    @test val(state(s,3)) == 3
  end

  @testset "tJ SiteSet" begin
    sites = tJSites(N)

    Nup_2 = op(sites,"Nup",2)
    @test Nup_2[2,2] ≈ 1.0
  end
end


text:
function extracellular{T}(e::Extracellular{T},n::Neuron,sigma::Float64)

    coeffs=zeros(Float64,0)
    inds=zeros(Int64,0)

    for i=1:length(n.secs)
        for j in n.secs[i].pnode
            if n.nodes[j].internal == true

                push!(inds,j)
                mypt3d=n.secs[i].pt3d[n.nodes[j].pt3d]
                myc=calc_coeffs(e,mypt3d,n.nodes[j].parent)
                push!(coeffs,(1/(4*pi*sigma))*myc)

            end
        end
    end

    Extra_coeffs(coeffs,inds)
end

calc_coeffs{T<:Point}(e::Extracellular{T},p::Array{Pt3d,1},par::Int64)=point_coeffs(p,e.xyz)

calc_coeffs{T<:Line}(e::Extracellular{T},p::Array{Pt3d,1},par::Int64)=line_coeffs(p,e.xyz)

function calc_coeffs{T<:Mixed}(e::Extracellular{T},p::Array{Pt3d,1},par::Int64)
    if par==0
        c=point_coeffs(p,e.xyz)
    else
        c=line_coeffs(p,e.xyz)
    end
    c
end

function extracellular{T<:Neuron}(n::Array{T,1},e::Extracellular)

    mycoeffs=Array(Extra_coeffs,length(n))

    for i=1:length(n)
        mycoeffs[i]=extracellular(e,n[i],0.3)
    end

    mycoeffs
end

function extracellular{T<:Line}(extra::Extracellular{T},xyz::Array{Float64,2},minds::Array{Array{Int64,1},1},sigma::Float64,cinds::Array{Int64,1})

    coeffs=zeros(Float64,div(size(xyz,1),2))

    i=1
    count=1
    while i<size(xyz,1)
        @inbounds coeffs[count]=(1/(4*pi*sigma))*line_coeffs(xyz,extra.xyz,i)
        i+=2
        count+=1
    end
    Extra_coeffs(coeffs,cinds)
end

function line_coeffs(pt3d::Array{Pt3d,1},xyz::Array{Float64,1})

    (unit_ds, delta_s)=pt3d_vec(pt3d[1],pt3d[end])

    d1=vect(pt3d[1],xyz)
    d2=vect(pt3d[end],xyz)

    #ln=dot(unit_ds,d1)
    hn=dot(unit_ds,d2)
    ln=delta_s+hn

    if (ln>0) && (hn<0)

        a=sqrt(sum(d2.^2))+abs(hn)
        b=(ln+sqrt(sum(d1.^2))) / abs(sum(d2.^2) - hn*hn)

        d=log(a*b)

    else

        a=sqrt(sum(d2.^2))+abs(hn)
        b=sqrt(sum(d1.^2))+abs(ln)

        d=log(a/b)

    end

    1/delta_s*d
end

function line_coeffs(xyz::Array{Float64,2},xyze::Array{Float64,1},i::Int64)

    @inbounds (ds_x,ds_y,ds_z,delta_s)=pt3d_vec(xyz[i,1],xyz[i,2],xyz[i,3],xyz[i+1,1],xyz[i+1,2],xyz[i+1,3])

    @inbounds (d1_x,d1_y,d1_z)=vect(xyz[i,1],xyz[i,2],xyz[i,3],xyze[1],xyze[2],xyze[3])
    @inbounds (d2_x,d2_y,d2_z)=vect(xyz[i+1,1],xyz[i+1,2],xyz[i+1,3],xyze[1],xyze[2],xyze[3])

    #ln=pdot(ds_x,ds_y,ds_z,d1_x,d1_y,d1_z)
    hn=pdot(ds_x,ds_y,ds_z,d2_x,d2_y,d2_z)
    ln=delta_s+hn

    if (ln>0) && (hn<0)

        a=sqrt(d2_x*d2_x+d2_y*d2_y+d2_z*d2_z)+hn
        b=(ln+sqrt(d1_x*d1_x+d1_y*d1_y+d1_z*d1_z)) / abs(d2_x*d2_x+d2_y*d2_y+d2_z*d2_z - hn*hn)

        d=log(a*b)

    else

        a=sqrt(d2_x*d2_x+d2_y*d2_y+d2_z*d2_z)+abs(hn)
        b=sqrt(d1_x*d1_x+d1_y*d1_y+d1_z*d1_z)+abs(ln)

        d=log(a/b)

    end

    1/delta_s*d
end

function point_coeffs(p::Array{Pt3d,1},xyz::Array{Float64,1})

    middle=round(Int64, length(p)/2)

    dist1=dist(p[middle],xyz)

    1/(dist1)
end

dist(p1::Pt3d, p2::Pt3d)=sqrt((p1.x-p2.x)^2+(p1.y-p2.y)^2+(p1.z-p2.z)^2)

dist(p1::Pt3d,xyz::Array{Float64,1})=sqrt((p1.x-xyz[1])^2+(p1.y-xyz[2])^2+(p1.z-xyz[3])^2)

dist(x1::Float64,y1::Float64,z1::Float64,x2::Float64,y2::Float64,z2::Float64)=sqrt((x1-x2)^2+(y1-y2)^2+(z1-z2)^2)

vect(p::Pt3d,xyz::Array{Float64,1})=[xyz[1]-p.x,xyz[2]-p.y,xyz[3]-p.z]

vect(p1::Pt3d,p2::Pt3d)=[p2.x-p1.x,p2.y-p1.y,p2.z-p1.z]

vect(x1::Float64,y1::Float64,z1::Float64,x2::Float64,y2::Float64,z2::Float64)=(x2-x1,y2-y1,z2-z1)

pt3d_vec(p1::Pt3d,p2::Pt3d)=(n=vect(p1,p2); mag=dist(p1,p2); n /= mag; (n, mag))

function pt3d_vec(x1::Float64,y1::Float64,z1::Float64,x2::Float64,y2::Float64,z2::Float64)
    (vx,vy,vz)=vect(x1,y1,z1,x2,y2,z2)
    mag=dist(x1,y1,z1,x2,y2,z2)
    (vx/mag,vy/mag,vz/mag,mag)
end

pdot(x1,y1,z1,x2,y2,z2)=x1*x2+y1*y2+z1*z2

function getv(im::Array{Float64,2},e::Extra_coeffs,v::Array{Float64,2},k::Int64)

    mystart=20

    for i=1:size(v,1)
        for j=1:length(e.c)
            @inbounds v[i,k]+=e.c[j]*im[e.inds[j],mystart+i-1]
        end
    end
    nothing
end

function getv(im::Array{Float64,2},e::Extra_coeffs,v::Array{Float64,3},k::Int64,l::Int64)

    mystart=20

    for i=1:size(v,1)
        for j=1:length(e.c)
            @inbounds v[i,k,l]+=e.c[j]*im[e.inds[j],mystart+i-1]
        end
    end
    nothing
end

function nete(n::Neuron,im::Array{Float64,2},ex::Extracellular,num::Int64)

    mystart=20
    v=zeros(Float64,length(mystart:size(im,2)),num)
    (xyz,inds)=getxyz(n);
    cinds=findcind(n)
    xyz1=copyn(xyz)

    for i=1:num
        resetn!(xyz,xyz1)
        randomize_shape!(n,xyz1,inds)
        translate3d!(xyz1,rand(-500.0:500.0),rand(-500.0:500.0),rand(-500.0:500.0))
        c=extracellular(ex,xyz1,inds,.3,cinds)
        getv(im,c,v,i)
    end
    v
end

function nete{T}(n::Neuron,im::Array{Float64,2},ex::Array{T,1},num::Int64)

    mystart=20
    v=zeros(Float64,length(mystart:size(im,2)),num,length(ex))
    (xyz,inds)=getxyz(n);
    cinds=findcind(n)
    xyz1=copyn(xyz)

    for i=1:num
        resetn!(xyz,xyz1)
        randomize_shape!(n,xyz1,inds)
        translate3d!(xyz1,rand(-500.0:500.0),rand(-500.0:500.0),rand(-500.0:500.0))
        for j=1:length(ex)
            c=extracellular(ex[j],xyz1,inds,.3,cinds)
            getv(im,c,v,i,j)
        end
    end
    v
end

function extrap(v::Array{Float64,2},t::Float64)

    v2=zeros(Float64,length(0.0:.025:t))

    spike=find_cspikes(v)
    if length(spike)>0
        cur_spike=spike[1]
        spike_ind=1
    else
        cur_spike=0
    end

    st=[Array(Int64,0) for i=1:length(spike)]

    firing=zeros(Int64,length(0.0:.025:t)-size(v,1))

    @inbounds for i=1:size(v,2)

        d=Poisson(100*rand()/40000)

        rand!(d,firing)

        j=1
        if i==cur_spike
            while j<length(firing)
                if firing[j]>0
                    push!(st[spike_ind],j)
                    for k=1:size(v,1)
                        if abs(v[k,i])>abs(v2[j])
                            v2[j]=v[k,i]+randn()*.01
                        end
                        j+=1
                    end
                end
                j+=1
            end
            spike_ind+=1
            if spike_ind<=length(spike)
                cur_spike=spike[spike_ind]
            end
        else
            while j<length(firing)
                if firing[j]>0
                    for k=1:size(v,1)
                        if abs(v[k,i])>abs(v2[j])
                            v2[j]=v[k,i]+randn()*.01
                        end
                        j+=1
                    end
                end
                j+=1
            end
        end
    end
    (v2,st)
end

function extrap(v::Array{Float64,3},t::Float64)

    v2=zeros(Float64,length(0.0:.025:t),size(v,3))

    spike=find_cspikes(mean(v,3))
    if length(spike)>0
        cur_spike=spike[1]
        spike_ind=1
    else
        cur_spike=0
    end

    st=[Array(Int64,0) for i=1:length(spike)]

    firing=zeros(Int64,length(0.0:.025:t)-size(v,1))

    @inbounds for i=1:size(v,2)

        d=Poisson(100*rand()/40000)

        rand!(d,firing)

        j=1
        if i==cur_spike
            while j<length(firing)
                if firing[j]>0
                    push!(st[spike_ind],j)
                    for k=1:size(v,1)
                        for l=1:size(v2,2)
                            if abs(v[k,i,l])>abs(v2[j,l])
                                v2[j,l]=v[k,i,l]+randn()*.01
                            end
                        end
                        j+=1
                    end
                end
                j+=1
            end
            spike_ind+=1
            if spike_ind<=length(spike)
                cur_spike=spike[spike_ind]
            end
        else
            while j<length(firing)
                if firing[j]>0
                    for k=1:size(v,1)
                        for l=1:size(v2,2)
                            if abs(v[k,i,l])>abs(v2[j,l])
                                v2[j,l]=v[k,i,l]+randn()*.01
                            end
                        end
                        j+=1
                    end
                end
                j+=1
            end
        end
    end
    (v2,st)
end

function find_cspikes(v)

    spike_num=Array(Int64,0)
    l=size(v,1)
    for i=1:size(v,2)
        for j=1:l
            @inbounds if v[j,i]<-0.12
                push!(spike_num,i)
                break
            end
        end
    end
    spike_num
end

function runc(n::NetworkS,init=false)

    myi=zeros(Float64,length(n.neur.N_1[1].i_vm),length(n.t))
    if init==true
        init!(n)
    end
    for i=1:5
	n.neur.N_1[1].v[i]=30.0
    end
    for i=0:2
	n.neur.N_1[1].v[end-i]=30
    end
    for i=1:length(n.t)

        @inbounds main(n.neur)

        myi[:,i]=n.neur.N_1[1].i_vm[:]
    end
    myi
end

max_stars_repo_name: paulmthompson/JNeuron

max_issues_repo_name: paulmthompson/JNeuron

max_forks_repo_name: paulmthompson/JNeuron


text:
module TestLandscapes

# 1D potentials
include("potentials1D.jl")
export SymmetricDoubleWell, AsymmetricDoubleWell, FatSkinnyDoubleWell10

# 2D potentials
include("potentials2D.jl")
export EntropicSwitch, SymmetricTwoChannel, Muller, Rosenbrock, Zpotential, EntropicBox

end

max_stars_repo_name: gideonsimpson/TestLandscapes

max_issues_repo_name: gideonsimpson/TestLandscapes

max_forks_repo_name: gideonsimpson/TestLandscapes


text:
export create_new_net_with_prior_counts

function create_new_net_with_prior_counts(net::Network, r::AbstractVector{Int}, d::AbstractMatrix{Int};
    priorcount :: Int = 1
    )

    nvars = length(r)
    stats = statistics(net, r, d)
    net2 = Network()

    for i = 1 : nvars
        node = get_node(net, i-1)
        name = get_name(node)
        nodeid = add_node(net2, DSL_CPT, name)
        set_number_of_outcomes(definition(get_node(net2, nodeid)), r[i])
        @assert(nodeid == i - 1)
    end
    for i = 1 : nvars
        nodeid = i - 1
        for p_id in get_parents(net, nodeid)
            add_arc(net2, p_id, nodeid)
        end
    end
    for i = 1 : nvars

        nodeid = i - 1

        ri = r[i]
        syscoord = SysCoordinates(definition(get_node(net2, nodeid)))
        counts = stats[i]

        # note that the ordering of parents is different than the ordering for parents in Bayes nets
        # Smile:     parent 1 is the most  significant bit so changes the least quickly
        # BayesNets: parent 1 is the least significant bit so changes the most  quickly

        parent_ids = get_parents(net, nodeid)

        nPa = length(parent_ids)
        if nPa > 0
            q = prod(r[parent_ids+1])
            dims = tuple(r[parent_ids+1]...)
            revdims = tuple(ri, reverse(dims)...)

            for ind_dmat in 1 : q
                sub = ind2sub(dims, ind_dmat)
                ind_dset = sub2ind(revdims, 1, sub...)-1

                # get the new normalized probability
                p = (counts[:,ind_dmat] .+ 1)
                p ./= sum(p)

                for j = 1 : ri
                    go_to(syscoord, ind_dset + j - 1)
                    set_unchecked_value(syscoord, p[j])
                end

            end
        else

            ind_dmat = 1
            ind_dset = 0

            # get the new normalized probability
            p = counts[:,ind_dmat] .+ priorcount
            p ./= sum(p)

            for j = 1 : ri
                go_to(syscoord, ind_dset + j - 1)
                set_unchecked_value(syscoord, p[j])
            end
        end
    end

    net2
end
function create_new_net_with_prior_counts(net::Network, r::AbstractVector{Int}, stats::Vector{Matrix{Float64}};
    priorcount :: Int = 1
    )

    nvars = length(r)
    net2 = Network()

    for i in 1 : nvars
        nodeid = add_node(net2, DSL_CPT, get_name(get_node(net, i-1)))
        set_number_of_outcomes(definition(get_node(net2, nodeid)), r[i])
        @assert(nodeid == i - 1)
    end
    for i in 1 : nvars
        nodeid = i - 1
        for p_id in get_parents(net, nodeid)
            add_arc(net2, p_id, nodeid)
        end
    end
    for i in 1 : nvars

        nodeid = i - 1

        ri = r[i]
        syscoord = SysCoordinates(definition(get_node(net2, nodeid)))
        counts = stats[i]

        # note that the ordering of parents is different than the ordering for parents in Bayes nets
        # Smile:     parent 1 is the most  significant bit so changes the least quickly
        # BayesNets: parent 1 is the least significant bit so changes the most  quickly

        parent_ids = get_parents(net, nodeid)

        nPa = length(parent_ids)
        if nPa > 0
            q = prod(r[parent_ids+1])
            dims = tuple(r[parent_ids+1]...)
            revdims = tuple(ri, reverse(dims)...)

            for ind_dmat in 1 : q
                sub = ind2sub(dims, ind_dmat)
                ind_dset = sub2ind(revdims, 1, sub...)-1

                # get the new normalized probability
                p = (counts[:,ind_dmat] .+ 1)
                p ./= sum(p)

                for j = 1 : ri
                    go_to(syscoord, ind_dset + j - 1)
                    set_unchecked_value(syscoord, p[j])
                end

            end
        else

            ind_dmat = 1
            ind_dset = 0

            # get the new normalized probability
            p = counts[:,ind_dmat] .+ priorcount
            p ./= sum(p)

            for j = 1 : ri
                go_to(syscoord, ind_dset + j - 1)
                set_unchecked_value(syscoord, p[j])
            end
        end
    end

    net2
end


text:
function newSol(y, μ)
	# This function is used to generate new point according to lower and upper
	# and a random factor proportional to current point.
	return(((1.0 .+ μ) .^ abs.(y) .- 1) / μ) .* sign.(y)
end


text:
using NiLang, NiLang.AD
using Test

@testset "check grad" begin
    for opm in [PlusEq, MinusEq]
        @test check_grad(opm(identity), (1.0, 2.0); verbose=true, iloss=1)
        @test check_grad(opm(*), (1.0, 2.0, 2.0); verbose=true, iloss=1)
        @test check_grad(opm(+), (1.0, 2.0, 2.0); verbose=true, iloss=1)
        @test check_grad(opm(-), (1.0, 2.0, 2.0); verbose=true, iloss=1)
        @test check_grad(opm(^), (1.0, 2.0, 2); verbose=true, iloss=1)
        @test check_grad(opm(^), (1.0, 2.0, 2.0); verbose=true, iloss=1)
        @test check_grad(opm(inv), (1.0, 2.0); verbose=true, iloss=1)
        @test check_grad(opm(sqrt), (1.0, 2.0); verbose=true, iloss=1)
        @test check_grad(opm(abs), (1.0, -2.0); verbose=true, iloss=1)
        @test check_grad(opm(abs2), (1.0, -2.0); verbose=true, iloss=1)
        @test check_grad(opm(exp), (1.0, 2.0); verbose=true, iloss=1)
        @test check_grad(opm(log), (1.0, 2.0); verbose=true, iloss=1)
        @test check_grad(opm(sin), (1.0, 2.0); verbose=true, iloss=1)
        @test check_grad(opm(sinh), (1.0, 2.0); verbose=true, iloss=1)
        @test check_grad(opm(asin), (1.0, 0.2); verbose=true, iloss=1)
        @test check_grad(opm(cos), (1.0, 2.0); verbose=true, iloss=1)
        @test check_grad(opm(cosh), (1.0, 2.0); verbose=true, iloss=1)
        @test check_grad(opm(acos), (1.0, 0.2); verbose=true, iloss=1)
        @test check_grad(opm(tan), (1.0, 2.0); verbose=true, iloss=1)
        @test check_grad(opm(tanh), (1.0, 2.0); verbose=true, iloss=1)
        @test check_grad(opm(atan), (1.0, -2.0); verbose=true, iloss=1)
        @test check_grad(opm(atan), (1.0, -2.0, 1.5); verbose=true, iloss=1)
        @test check_grad(opm(convert), (Fixed43(0.5), 2.0); verbose=true, iloss=1)
        @test check_grad(opm(/), (1.0, 2.0, 2.0); verbose=true, iloss=1)
        @test check_grad(opm(min), (1.0, 2.0, 3.0); verbose=true, iloss=1)
        @test check_grad(opm(max), (1.0, 2.0, 3.0); verbose=true, iloss=1)
        @test check_grad(opm(min), (1.0, 3.0, 2.0); verbose=true, iloss=1)
        @test check_grad(opm(max), (1.0, 3.0, 2.0); verbose=true, iloss=1)
        @test_broken check_grad(opm(÷), (1.0, 2.0, 2.0); verbose=true, iloss=1)
        @test gradient(opm(sqrt), (1.0, 0.0); iloss=1)[2] == 0
    end
    @test check_grad(NEG, (1.0,); verbose=true, iloss=1)
    @test check_grad(INV, (3.0,); verbose=true, iloss=1)
    @test check_grad(AddConst(2.0), (3.0,); verbose=true, iloss=1)
    @test check_grad(SubConst(2.0), (3.0,); verbose=true, iloss=1)
    @test check_grad(INC, (1.0,); verbose=true, iloss=1)
    @test check_grad(DEC, (1.0,); verbose=true, iloss=1)
    @test check_grad(ROT, (1.0, 2.0, 2.0); verbose=true, iloss=1)
    @test check_grad(ROT, (1.0, 2.0, 2.0); verbose=true, iloss=2)
    @test check_grad(IROT, (1.0, 2.0, 2.0); verbose=true, iloss=1)
    @test check_grad(IROT, (1.0, 2.0, 2.0); verbose=true, iloss=2)
    @test check_grad(HADAMARD, (3.0, 2.0); verbose=true, iloss=1)
    @test check_grad(HADAMARD, (3.0, 2.0); verbose=true, iloss=2)
end

@testset "partial gvar" begin
    @i function testf1(f, a, b)
	f(a, b, 2.0)
    end
    @i function testf2(f, a, b)
	f(a, 2.0, b)
    end
    for testf in [testf1, testf2]
    	for opm in [PlusEq, MinusEq]
            @test check_grad(testf, (opm(*), 1.0, 2.0); verbose=true, iloss=2)
            @test check_grad(testf, (opm(+), 1.0, 2.0); verbose=true, iloss=2)
            @test check_grad(testf, (opm(-), 1.0, 2.0); verbose=true, iloss=2)
            @test check_grad(testf, (opm(^), 1.0, 2.0); verbose=true, iloss=2)
            @test check_grad(testf, (opm(atan), 1.0, -2.0); verbose=true, iloss=2)
            @test check_grad(testf, (opm(/), 1.0, 2.0); verbose=true, iloss=2)
	end
    end
    @test check_grad(testf1, (ROT, 1.0, 2.0); verbose=true, iloss=2)
    @test check_grad(testf1, (ROT, 1.0, 2.0); verbose=true, iloss=3)
    @test check_grad(testf1, (IROT, 1.0, 2.0); verbose=true, iloss=2)
    @test check_grad(testf1, (IROT, 1.0, 2.0); verbose=true, iloss=3)
    # ROT and HADAMARD does not allow different types of rotation elements
end

@testset "sincos" begin
    @i function f(s, c, x)
        (s, c) += sincos(x)
    end
    @test check_grad(f, (1.0, 2.0, 2.0); verbose=true, iloss=1)
    @test check_grad(f, (1.0, 2.0, 2.0); verbose=true, iloss=2)
end

@testset "AD over pop" begin
    @i function mean(out!::T, x) where T
        anc ← zero(out!)
        for i=1:length(x)
            anc += x[i]
        end
        out! += anc / (@const length(x))
        FLOAT64_STACK[end+1] ↔ anc::T
    end

    @test check_grad(mean, (0.0, [1,2,3.0, 4.0]); iloss=1)
end

@testset "AD over pipe" begin
    @i function mean(out!, anc, x)
        for i=1:length(x)
            PlusEq(identity)(anc, x[i])
            SWAP(anc, x[i])
        end
        out! += anc / (@const length(x))
    end

    @test check_grad(mean, (0.0, 0.0, [1,2,3.0, 4.0]); iloss=1, verbose=true)
end

@testset "push, load data" begin
    stack = []
    val = [1,2,3]
    @instr PUSH!(stack, val)
    @test val == Int[]
    val = 3.0
    @instr PUSH!(stack, val)
    @test val == 0.0
    val = 3.0
    @instr PUSH!(stack, val)
    x = GVar(3.0)
    #@test_throws InvertibilityError @instr POP!(stack, x)
    z = 3.0
    @instr PUSH!(stack, z)
    z = GVar(0.0)
    @instr POP!(stack, z)
    @test z == GVar(3.0)
    x = [1.0, 2.0, 3.0]
    @instr PUSH!(stack, x)
    y = empty(x)
    @instr POP!(stack, y)
    @test y == GVar.([1,2,3.0])
    x = [1.0, 2.0, 3.0]
    @instr PUSH!(stack, x)
    y = empty(x)
    @instr POP!(stack, y)
    @test y == [1,2,3.0]
end

@testset "dataviews" begin
    @i function f(z, y, x)
        y += cos(x |> INV)
        z += tan(y |> AddConst(4.0))
        z += y * (x |> NEG |> SubConst(0.5) |> INV)
        z += sin(x |> INV)
    end
    @test check_grad(f, (0.2, 0.5, 0.8); iloss=1)
end

@testset "additive identity" begin
    struct TestAdd2{T}
        x::T
        y::Vector{T}
    end
    x = TestAdd2(GVar(1.0, 2.0), [GVar(2.0, 1.2)])
    y = TestAdd2(GVar(6.0, 3.0), [GVar(4.0, 4.1)])
    @test getfield.(MinusEq(identity)(x, y), :x) == getfield.((TestAdd2(GVar(-5.0, 2.0), [GVar(-2.0, 1.2)]), TestAdd2(GVar(6.0, 5.0), [GVar(4.0, 5.3)])), :x)
    x = TestAdd2(GVar(1.0, 2.0), [GVar(2.0, 1.2)])
    y = TestAdd2(GVar(6.0, 3.0), [GVar(4.0, 4.1)])
    @test getfield.(MinusEq(identity)(x, y), :y) == getfield.((TestAdd2(GVar(-5.0, 2.0), [GVar(-2.0, 1.2)]), TestAdd2(GVar(6.0, 5.0), [GVar(4.0, 5.3)])), :y)
end


text:
module Bound
using LinearAlgebra

# Estimate constant in Proposition 6 using Frobenius norm.
function estimate_constants(M0, ops)
        c0 = ops.h * LinearAlgebra.opnorm(M0)
        c1 = ops.h * ops.Hxi[1,1]
        return c0, c1
end

end

max_stars_repo_name: ooreilly/sbp

max_issues_repo_name: ooreilly/sbp

max_forks_repo_name: ooreilly/sbp


text:
import Combinatorics

"""
For every right kernel vector of M with nz nonzeros, finds
a vector [i1, ..., inz] of the  indices of the nonzero
entries. Returns an array with at most 'count' index vectors.
If 'count' is negative, finds all (default).
"""
function kernelnz(M::AbstractMatrix, nz::Integer; count::Integer = -1)
    @assert nz > 0
    M = full(M)
    ker = Vector{Vector{Int}}()
    for c in Combinatorics.combinations(1 : size(M, 2), nz)
        @assert length(c) == nz
        0 ≤ count ≤ length(ker) && break
        if rank(M[:,c]) < nz
            push!(ker, c)
        end
    end
    return ker
end

"Returns true if there is a vector v with `nz` non-zeros such that M⋅v = 0."
haskernelnz(M::AbstractMatrix, nz::Integer) = !isempty(kernelnz(M, nz; count = 1))

"""
Matrix flip around a diagonal:
   (a b     (d b
    c d) ->  c a)
"""
function transpoflip(A::Matrix)
    @assert size(A, 1) == size(A, 2)
    N = size(A, 1)
    B = [A[N + 1 - j, N + 1 - i] for i = 1 : N, j = 1 : N]
end

"""
Given the inverse A of a matrix M, returns the inverse matrix
of the matrix M', which differs from M only in n'th diagonal
entry, for which M'[n,n] = M[n,n] + δ. Uses the Sherman-Morrison-Woodbury
formula. See D.S. Bernestein, Matrix Mathematics 2009, Fact 2.16.3.
"""
function invupdate end

function invupdate(A::AbstractMatrix, n::Integer, δ::Real)
    B = copy(A)
    for i = 1 : size(A, 1), j = 1 : size(A, 2)
        B[i,j] = A[i,j] - δ * A[i,n] * A[n,j] / (1 + δ * A[n,n])
    end
    return B
end

function invupdate(A::Hermitian, n::Integer, δ::Real)
    B = copy(A)
    for i = 1 : size(A, 1), j = 1 : i
        Bij = A[i,j] - δ * A[i,n] * A[n,j] / (1 + δ * A[n,n])
        B.data[i,j] = B.data[j,i]= Bij
    end
    return B
end


text:
module Train20

import Base: IteratorSize, IteratorEltype, length, size, iterate, eltype, rand, repeat, summary, show
using Base: @propagate_inbounds, tail, haslength, SizeUnknown
using Base.Iterators: Cycle
using Random: randn, rand, randperm
using Printf: @sprintf
using AutoGrad: AutoGrad, Param, @diff, full
using LinearAlgebra: norm, lmul!, axpy!

include("data.jl"); export minibatch, Data
include("distributions.jl"); export gaussian, xavier, xavier_uniform, xavier_normal, bilinear
include("hyperopt.jl"); export goldensection, hyperband
include("param.jl"); export param, param0
include("progress.jl"); export progress, progress!
include("train.jl"); export minimize, minimize!, converge, converge!, train!
include("update.jl"); export update!, clone, optimizers, SGD, Sgd, sgd, sgd!, Momentum, momentum, momentum!, Nesterov, nesterov, nesterov!, Adagrad, adagrad, adagrad!, Rmsprop, rmsprop, rmsprop!, Adadelta, adadelta, adadelta!, Adam, adam, adam!
include("train_ka.jl") # defines param, update!, _optimizers for KnetArray

end # module


text:
"""
    get_spatial_profile_data_to_save(problem, mode_index)

return (data = ( x=DCM, y=ψ² ), fit=(x=x_fit, y=y_fit) )
"""
function get_spatial_profile_data_to_save(problem, mode_index)
    prepare_spectrum(problem)
    DCM, ψ² = get_spatial_profile_single_mode(problem, mode_index)
    x_fit, y_fit = get_decay_fit(problem, mode_index)
    return (data=(x=DCM, y=ψ²), fit=(x=x_fit, y=y_fit))
end

function plot_atoms_and_mode(problem, mode_index)
    # pyplot()
    # compute eigenvalues if not available
    prepare_spectrum(problem)
    DCM, ψ² = get_spatial_profile_single_mode(problem, mode_index)
    x_fit, y_fit = get_decay_fit(problem, mode_index)

    fig_space = plot_space(problem.atoms, ψ²)
    plot!(; guidefont=16, tickfont=15, size=(800, 600))

    fig_mode = scatter(DCM, ψ²; label="Mode Index = $(mode_index)", yscale=:log10)
    plot!(
        x_fit,
        y_fit;
        xlabel=L"|\mathbf{r}-\mathbf{r}_{cm}|",
        ylabel=L"|\psi_n|^2",
        size=(800, 600),
        guidefont=16,
        tickfont=15,
        legendfont=14,
        label="",
        lw=5,
        linestyle=:dash,
        c=:black,
        grid=false,
        framestyle=:box,
    )

    return fig_space, fig_mode
end

function plot_space(atoms::T where {T<:ThreeD}, ψ²)
    return scatter3d(
        get_X_axes(atoms),
        get_Y_axes(atoms),
        get_Z_axes(atoms);
        label="",
        xlabel="x",
        ylabel="x",
        zlabel="z",
        zcolor=ψ²,
        colorbar_title="|ψ|²",
        camera=(45, 30),
    )
end

prepare_spectrum(problem) = get_spectrum(problem)

function select_matrix_axes(static_array, component)
    N = length(static_array)
    dimensions = length(static_array[1])
    if component ≤ dimensions
        return [static_array[n][component] for n in 1:N]
    else
        return [zero(eltype(static_array[n])) for n in 1:N]
    end
end

get_X_axes(static_array::Vector{StaticArrays.SArray}) = select_matrix_axes(static_array, 1)
get_Y_axes(static_array::Vector{StaticArrays.SArray}) = select_matrix_axes(static_array, 2)
get_Z_axes(static_array::Vector{StaticArrays.SArray}) = select_matrix_axes(static_array, 3)

get_X_axes(atoms::T where {T<:ThreeD}) = select_matrix_axes(atoms.r, 1)
get_Y_axes(atoms::T where {T<:ThreeD}) = select_matrix_axes(atoms.r, 2)
get_Z_axes(atoms::T where {T<:ThreeD}) = select_matrix_axes(atoms.r, 3)

get_X_axes(atoms::T where {T<:TwoD}) = select_matrix_axes(atoms.r, 1)
get_Y_axes(atoms::T where {T<:TwoD}) = select_matrix_axes(atoms.r, 2)
get_Z_axes(atoms::T where {T<:TwoD}) = zeros(length(atoms.r))

get_atoms_matrix(atoms) = convert_StaticArray_to_matrix(atoms.r)

function get_decay_fit(problem, mode_index)
    DCM, ψ² = get_spatial_profile_single_mode(problem, mode_index)

    x_use, y_use = select_points(DCM, ψ²)
    A_L₁, ξ_L₁, R1, y_fit_L₁ = get_coeffs_L1(x_use, y_use)

    if is_R1_localized(R1)
        return x_use, y_fit_L₁
    else
        A_L₂, ξ_L₂, R2, y_fit_L₂ = get_coeffs_L2(x_use, y_use)
        return x_use, y_fit_L₂
    end
end


text:
## -------------------------------------------------------------------
let

    METHODS = [

        :FBA_Z_FIX_MAX_VG_MIN_COST,
        :FBA_Z_FIX_MAX_VG_MAX_COST,
        :FBA_Z_FIX_MIN_VG_MIN_COST,
        :FBA_Z_FIX_MIN_VG_MAX_COST,
        :FBA_Z_FIX_ADJ_MAX_VG_MIN_COST,
        :FBA_Z_FIX_ADJ_MAX_VG_MAX_COST,
        :FBA_Z_FIX_ADJ_MIN_VG_MIN_COST,
        :FBA_Z_FIX_ADJ_MIN_VG_MAX_COST,

        ME_MAX_POL,
        ME_Z_EXPECTED_G_BOUNDED,
        # ME_Z_OPEN_G_OPEN
    ]

    fontsize = 13
    common_params = (;
        thickness_scaling = 1.3,
        xguidefontsize = fontsize, yguidefontsize = fontsize,
        # dpi = 1000
    )

    ps = Plots.Plot[]
    for method in METHODS
        p = plot(;
            xlabel = "vg",
            ylabel = "cost",
            common_params...
        )

        for (iJR, Data) in [
            (ChK.iJR904, ChK.KayserData),
            (ChN.iJR904, ChN.NanchenData),
            (ChF.iJR904, ChF.FolsomData),
            (ChH.iJR904, ChH.HeerdenData),
        ]
            src = nameof(Data)

            marker = (6, source_markers[Data])

            DAT = ChE.load_DAT(src)

            EXPS = DAT[:EXPS]

            @info("Doing", src, method)

            vg_vals = DAT[method, :flx, "GLC", EXPS]
            cost_vals = DAT[method, :flx, "cost", EXPS]

            scatter!(p, abs.(vg_vals), abs.(cost_vals);
                label = "",
                color = :black,
                marker, alpha = 0.8,
            )
        end

        plot!(;legend = :bottomright, ylim = [0.0, Inf])
        mysavefig(p, "vg_cost_study"; method)
        push!(ps, plot!(p; title = string(method)))
    end

    mysavefig(ps, "vg_cost_study")

end


text:
# Constants
const limit_kernel_width_default = 10.0
const speed_of_light = 299792458.0 # m/s

# spectrum.jl:  utils for computing simple spectrum models
include("spectrum.jl")
export AbstractSpectrum, SimpleSpectrum, SimulatedSpectrum, ConvolvedSpectrum
export doppler_shifted_spectrum
#export continuum_model, absorption_lines, gaussian_convolution_kernel
#export std_gaussian_line, absorption_line

# convolution_kernels.jl:  utils for convolution kernels
include("convolution_kernels.jl")
export AbstractConvolutionKernel, GaussianConvolutionKernel, GaussianMixtureConvolutionKernel


text:
include("plain.jl")
include("jld.jl")
include("readremote.jl")
include("extend_test.jl")
include("gc.jl")
if Pkg.installed("DataFrames") != nothing
    include("jld_dataframe.jl")
end


text:
"""
This file contains all the functions that could be commonly used to convert
dates to strings.  This is a very important and necessary tool in climate
observation and reanalysis data, because files are often saved in certain
formats that include the date and time of the observation.

There following types of functions are contained in this .jl module:
1)  single-unit conversion functions (e.g. yr2str) that extracts only a single
    time-unit and converts it to a string
2)  multi-unit conversion functions (e.g. ymd2str) that extracts and combines
    various different time units into a string
3)  conversion of dates into folder directories

"""

# Single-Unit Functions
function yr2str(dnt::TimeType)
    return @sprintf("%04d",Dates.year(dnt));
end

function mo2str(dnt::TimeType)
    return @sprintf("%02d",Dates.month(dnt));
end

function dy2str(dnt::TimeType)
    return @sprintf("%02d",Dates.day(dnt));
end

function hr2str(dnt::TimeType)
    return @sprintf("%02d",Dates.hour(dnt));
end

function mi2str(dnt::TimeType)
    return @sprintf("%02d",Dates.minute(dnt));
end

# Multi-Unit Functions
function yrmo2str(dnt::TimeType)
    yr = Dates.year(dnt); mo = Dates.month(dnt);
    return @sprintf("%04d%02d",yr,mo);
end

function ymd2str(dnt::TimeType)
    yr = Dates.year(dnt); mo = Dates.month(dnt); dy = Dates.day(dnt);
    return @sprintf("%04d%02d%02d",yr,mo,dy);
end

function ymdhm2str(dnt::TimeType)
    yr = Dates.year(dnt); mo = Dates.month(dnt); dy = Dates.day(dnt);
    hr = Dates.hour(dnt); mi = Dates.minute(dnt);
    return @sprintf("%04d%02d%02d%02d%02d00",yr,mo,dy,hr,mi);
end

# Directory Functions
function yrmo2dir(dnt::TimeType)
    yr = Dates.year(dnt); mo = Dates.month(dnt);
    return @sprintf("%04d/%02d",yr,mo);
end

function ymd2dir(dnt::TimeType)
    yr = Dates.year(dnt); mo = Dates.month(dnt); dy = Dates.day(dnt);
    return @sprintf("%04d/%02d/%02d",yr,mo,dy);
end

function yrdy2dir(dnt::TimeType)
    yr = Dates.year(dnt); dy = Dates.dayofyear(dnt);
    return @sprintf("%04d/%02d",yr,dy);
end


text:
# ---
# title: 816. Ambiguous Coordinates
# id: problem816
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String
# link: <https://leetcode.com/problems/ambiguous-coordinates/description/>
# hidden: true
# ---
#
# We had some 2-dimensional coordinates, like `"(1, 3)"` or `"(2, 0.5)"`.  Then,
# we removed all commas, decimal points, and spaces, and ended up with the
# string `S`.  Return a list of strings representing all possibilities for what
# our original coordinates could have been.
#
# Our original representation never had extraneous zeroes, so we never started
# with numbers like "00", "0.0", "0.00", "1.0", "001", "00.01", or any other
# number that can be represented with less digits.  Also, a decimal point within
# a number never occurs without at least one digit occuring before it, so we
# never started with numbers like ".1".
#
# The final answer list can be returned in any order.  Also note that all
# coordinates in the final answer have exactly one space between them (occurring
# after the comma.)
#
#
#
#     **Example 1:**
#     Input: "(123)"
#     Output: ["(1, 23)", "(12, 3)", "(1.2, 3)", "(1, 2.3)"]
#
#
#
#     **Example 2:**
#     Input: "(00011)"
#     Output:  ["(0.001, 1)", "(0, 0.011)"]
#     Explanation:
#     0.0, 00, 0001 or 00.01 are not allowed.
#
#
#
#     **Example 3:**
#     Input: "(0123)"
#     Output: ["(0, 123)", "(0, 12.3)", "(0, 1.23)", "(0.1, 23)", "(0.1, 2.3)", "(0.12, 3)"]
#
#
#
#     **Example 4:**
#     Input: "(100)"
#     Output: [(10, 0)]
#     Explanation:
#     1.0 is not allowed.
#
#
#
#
# **Note:**
#
#   * `4 <= S.length <= 12`.
#   * `S[0]` = "(", `S[S.length - 1]` = ")", and the other elements in `S` are digits.
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# functions: nablasq_3d_grid, return_boundary_nodes,
# return_boundary_nodes_3D, punch_holes_nexus_Cartesian,
# Matern_3d_Grid, Laplace_3D_grid,
# parallel_Matern_3DGrid, parallel_Laplace_3Dgrid

# Internal settings for caching
mutable struct Settings
  A_matrix_STORE_MAX::Int64
end

# There may be up to 15 A_matrices for a 3d grid
const SETTINGS = Settings(15)

# Dict for storing A_matrices
const A_matrix = Dict{Tuple{Int64, Int64, Int64, Int64, Float64, Float64, Float64,
                            Float64}, Matrix{Float64}}()
"""
  nablasq_3d_grid(Nx,Ny)

Construct the 3D Laplace matrix

# Arguments
  - `Nx::Int64`: The number of nodes in the first dimension
  - `Ny::Int64`: The number of nodes in the second dimension
  - `Nz::Int64`: The number of nodes in the third dimension
  - `h::Float64`: Grid spacing in the first dimension
  - `k::Float64`: Grid spacing in the second dimension
  - `l::Float64`: Grid spacing in the third dimension

# Outputs

  - `-nablasq` (discrete Laplacian, real-symmetric positive-definite) on Nx×Ny×Nz grid

"""
function nablasq_3d_grid(Nx, Ny, Nz, h, k, l)
    #haskey(A_matrix, (Nx, Ny, Nz, 1, 0.0, h, k, l)) && return A_matrix[(Nx, Ny, Nz, 1, 0.0, h, k, l)]
    o₁ = ones(Nx) / h
    del1 = spdiagm_nonsquare(Nx + 1, Nx, -1 => -o₁, 0 => o₁)
    o₂ = ones(Ny) / k
    del2 = spdiagm_nonsquare(Ny + 1, Ny, -1 => -o₂,0 => o₂)
    O3 = ones(Nz) / l
    del3 = spdiagm_nonsquare(Nz + 1, Nz, -1 => -O3, 0 => O3)
    A3D = (kron(sparse(I, Nz, Nz), sparse(I, Ny, Ny), del1'*del1) +
            kron(sparse(I, Nz, Nz), del2' * del2, sparse(I, Nx, Nx)) +
            kron(del3' * del3, sparse(I, Ny, Ny), sparse(I, Nx, Nx)))
    BoundaryNodes, xneighbors, yneighbors, zneighbors =
            return_boundary_nodes(Nx, Ny, Nz)
    count = 1
    for i in BoundaryNodes
        A3D[i, i] = 0.0
        A3D[i, i] = A3D[i, i] + xneighbors[count] / h ^ 2 +
                     yneighbors[count] / k ^ 2 + zneighbors[count] / l ^ 2
        count = count + 1
    end
    #length(A_matrix) <  SETTINGS.A_matrix_STORE_MAX && (A_matrix[(Nx, Ny, Nz, 1, 0.0, h, k, l)] = A3D)
    #if length(A_matrix) == SETTINGS.A_matrix_STORE_MAX
    #  @warn "A_matrix cache full, no longer caching Laplace interpolation matrices."
    #end
    return A3D
end

""" Helper function to give the matern matrix """
function _Matern_matrix(Nx, Ny, Nz, m, eps, h, k, l)
    #haskey(A_matrix, (Nx, Ny, Nz, m, eps, h, k, l)) && return A_matrix[(Nx, Ny, Nz, m, eps, h, k, l)]
    A3D = nablasq_3d_grid(Nx, Ny, Nz, h, k, l)
    sizeA = size(A3D, 1)
    for i = 1:sizeA
        A3D[i, i] = A3D[i, i] + eps^2
    end
    A3DMatern = A3D^m
    #length(A_matrix) <  SETTINGS.A_matrix_STORE_MAX && (A_matrix[(Nx, Ny, Nz, m, eps, h, k, l)] = A3DMatern)
    #if length(A_matrix) == SETTINGS.A_matrix_STORE_MAX
    #  @warn "A_matrix cache full, no longer caching Laplace interpolation matrices."
    #end
    A3DMatern
end

"""

  matern_3d_grid(imgg, discard, m, eps, h, k, l)

Interpolates a single punch

...
# Arguments
  - `imgg`: the matrix containing the image
  - `discard::Union{Vector{CartesianIndex{3}}}, Vector{Int64}}`: the linear or
       Cartesian indices of the values to be filled
  - `m::Int64 = 1` : Matern parameter
  - `eps::Float64 = 0.0`: Matern parameter eps
  - `h = 1.0`: Aspect ratio in the first dimension
  - `k = 1.0`: Aspect ratio in the second dimension
  - `l = 1.0`: Aspect ratio in the third dimension

# Outputs
  - array containing the restored image
...

"""
function matern_3d_grid(imgg, discard::Union{Vector{CartesianIndex{3}}, Vector{Int64}},
                m::Int64 = 1,  eps::Float64 = 0.0,
                h = 1.0, k = 1.0, l = 1.0)
    Nx, Ny, Nz = size(imgg)
    A3D = (eps == 0.0)&&(m == 1) ?
                nablasq_3d_grid(Nx, Ny, Nz, h, k, l) :
                _Matern_matrix(Nx, Ny, Nz, m, eps, h, k, l)
    totalsize = Nx * Ny * Nz
    C = sparse(I, totalsize, totalsize)
    rhs_a = copy(imgg)[:]
    for i in discard
        j = (typeof(i) <: CartesianIndex) ? LinearIndices(imgg)[i] : i
        C[j, j] = 0.0
        rhs_a[j] = 0.0
    end
    Id = sparse(I, totalsize, totalsize)
    u = ((C - (Id - C) * A3D)) \ rhs_a
    return reshape(u, Nx, Ny, Nz)
end

# Add m pixels around the punch and then intersect with the size of the full
# image (3D only)
function pad_intersect(discard, m, Nx, Ny, Nz)
    B = bounding_box(discard)
    fi, li = (first(B) - CartesianIndex(m, m, m), last(B) + CartesianIndex(m, m, m))
    return intersect_box(fi:li, CartesianIndices((Nx,Ny,Nz)))
end

function get_Qw_maxmin(wpoints)
    Qw_min, Qw_max = (Int64(floor(minimum(wpoints))), Int64(ceil(maximum(wpoints))))
end

"""

  matern_w_punch(imgg, Qh_min, Qh_max, Qk_min, Qk_max, Ql_min, Ql_max, m, eps, h, k, l, symm)

Interpolate, in serial, multiple punches

...
# Arguments
  - `imgg`: the matrix containing the image
  - ` Qh_min, Qh_max, Qk_min, Qk_max, Ql_min, Ql_max::Int64`: the extents in h,k,l resp
  - `m::Int64 = 1` : Matern parameter
  - `eps::Float64 = 0.0`: Matern parameter eps
  - `h = 1.0`: Aspect ratio in the first dimension
  - `k = 1.0`: Aspect ratio in the second dimension
  - `l = 1.0`: Aspect ratio in the third dimension

# Outputs
  - array containing the interpolated image

# Example

```<julia-repl>
h = k = l = 1.0
symm = 'A'
Nx = Ny = Nz = 61
radius = (0.5, 0.5, 0.5 )
Qh_min = Qk_min = Ql_min = -3.0
Qh_max = Qk_max = Ql_max = 3.0
xpoints = ypoints = zpoints = LinRange(Qh_min, Qh_max, Nx)
imgg = rand(Nx, Ny, Nz)
m = 1
eps = 0.0
interp = matern_w_punch(imgg, #Qh_min, Qh_max, Qk_min, Qk_max, Ql_min, Ql_max, radius,
                      radius, xpoints, ypoints, zpoints, m, eps,
                      h, k, l, symm);
```

...
"""
function matern_w_punch(imgg, radius, xpoints, ypoints, zpoints,
                        m = 1, eps = 0.0, h = 1.0, k = 1.0, l = 1.0, symm = 'G';
                        return_punch_locs = false)
    Qh_min, Qh_max = get_Qw_maxmin(xpoints)
    Qk_min, Qk_max = get_Qw_maxmin(ypoints)
    Ql_min, Ql_max = get_Qw_maxmin(zpoints)
    radius_x, radius_y, radius_z = (typeof(radius) <: Tuple) ? radius :
                                                (radius, radius, radius)
    dx, dy, dz = (xpoints[2] - xpoints[1], ypoints[2] - ypoints[1],
                  zpoints[2] - zpoints[1])
    rpx, rpy, rpz = (Int64(round(radius_x/dx)) + m, Int64(round(radius_y/dy)) + m,
                  Int64(round(radius_z/dz))+ m)
    Nx, Ny, Nz = (length(xpoints), length(ypoints), length(zpoints))
    new_imgg = copy(imgg);
    exclusion_rule = center_check(symm)
    # find the first center
    (cpx, cpy, cpz) = (findmin(abs.(xpoints .- Qh_min))[2],
                                  findmin(abs.(ypoints .- Qk_min))[2],
                                  findmin(abs.(zpoints .- Ql_min))[2])
    stridex, stridey, stridez = (-cpx, -cpy, -cpz) .+ (findmin(abs.(xpoints .- (Qh_min + 1)))[2],
                                  findmin(abs.(ypoints .- (Qk_min + 1)))[2],
                                  findmin(abs.(zpoints .- (Ql_min + 1)))[2])
    punch_locs = return_punch_locs ? ones(Int64, size(imgg)) : nothing
    # Threads.@threads for c in centers
    for (ih, h) in enumerate(Qh_min:Qh_max)
        minpx = maximum([1, cpx + (ih - 1) * stridex - rpx])
        maxpx = minimum([cpx + (ih - 1) * stridex + rpx, Nx])
        for (ik, k) in enumerate(Qk_min:Qk_max)
            minpy = maximum([1, cpy + (ik - 1) * stridey - rpy])
            maxpy = minimum([cpy + (ik - 1) * stridey + rpy, Ny])
            for (il, l) in enumerate(Ql_min:Ql_max)
                minpz = maximum([1, cpz + (il - 1) * stridez - rpz])
                maxpz = minimum([cpz + (il - 1) * stridez + rpz, Nz])
                if exclusion_rule(h, k, l)
                    d = punch_3D_cart((h,k,l), radius, xpoints[minpx:maxpx],
                                      ypoints[minpy:maxpy],
                                      zpoints[minpz:maxpz])
                    # Interpolate
                    new_imgg[minpx:maxpx, minpy:maxpy, minpz:maxpz] = matern_3d_grid(
                                imgg[minpx:maxpx, minpy:maxpy, minpz:maxpz], d, m,
                                eps, h, k, l);
                    if return_punch_locs
                      punch_locs[d] .= 0
                    end
                end
            end
        end
    end
    if return_punch_locs
      return new_imgg, punch_locs
    else
      return new_imgg
    end
end
# function matern_w_punch(imgg, Qh_min, Qh_max, Qk_min, Qk_max, Ql_min, Ql_max, radius,
                      # xpoints, ypoints, zpoints,
                        # m = 1, eps = 0.0, h = 1.0, k = 1.0, l = 1.0, symm = 'G')
    # centers = center_list(symm, Qh_min, Qh_max, Qk_min, Qk_max, Ql_min, Ql_max)
    # radius_x, radius_y, radius_z = (typeof(radius) <: Tuple) ? radius :
                                                # (radius, radius, radius)
    # dx, dy, dz = (xpoints[2] - xpoints[1], ypoints[2] - ypoints[1],
                  # zpoints[2] - zpoints[1])
    # rpx, rpy, rpz = (Int64(round(radius_x/dx)), Int64(round(radius_y/dy)),
                  # Int64(round(radius_z/dz)))
    # Nx, Ny, Nz = (length(xpoints), length(ypoints), length(zpoints))
    # new_imgg = copy(imgg);
    # # Threads.@threads for c in centers
    # for c in centers
        # (cpx, cpy, cpz) = (findmin(abs.(xpoints .- c[1]))[2],
                                      # findmin(abs.(ypoints .- c[2]))[2],
                                      # findmin(abs.(zpoints .- c[3]))[2])
        # minpx, minpy, minpz = (maximum([1, cpx - (rpx + m)]),
                               # maximum([1, cpy - (rpy + m)]),
                               # maximum([1, cpz - (rpz + m)]))
        # maxpx, maxpy, maxpz = (minimum([Nx, cpx + (rpx + m)]),
                               # minimum([Ny, cpy + (rpx + m)]),
                               # minimum([Nz, cpz + (rpz + m)]))
        # d = punch_3D_cart(c, radius, xpoints[minpx:maxpx], ypoints[minpy:maxpy],
                          # zpoints[minpz:maxpz])
        # # Interpolate
        # new_imgg[minpx:maxpx, minpy:maxpy, minpz:maxpz] = matern_3d_grid(
                    # imgg[minpx:maxpx, minpy:maxpy, minpz:maxpz], d, m, eps, h, k, l);
    # end
    # return new_imgg
# end

max_stars_repo_name: rayosborn/nxpeaks

max_issues_repo_name: rayosborn/nxpeaks

max_forks_repo_name: axmas-anl/nxrefine


text:
module TestTypeUtilities
using Test
using EcologicalNetworks
using LinearAlgebra

B = BipartiteProbabilisticNetwork(rand(3, 5))
@test typeof(copy(B)) == typeof(B)
@test size(B) == (3,5)

A = BipartiteNetwork([false true; false false])
@test typeof(A) <: DeterministicNetwork

A = UnipartiteNetwork([false true; false false])
@test typeof(A) <: DeterministicNetwork

# Is there an interaction?
A = UnipartiteNetwork([false true; true false])
@test has_interaction(A, 1, 2)
@test !has_interaction(A, 1, 1)

A = UnipartiteQuantitativeNetwork([0.0 1.0; 1.0 0.0])
@test has_interaction(A, 1, 2)
@test !has_interaction(A, 1, 1)

A = UnipartiteNetwork([false true; true false], [:a, :b])
@test has_interaction(A, :a, :b)
@test has_interaction(A, :b, :a)
@test !has_interaction(A, :b, :b)
@test !has_interaction(A, :a, :a)

A = BipartiteNetwork([false true; true false], [:a, :b], [:c, :d])
@test has_interaction(A, :a, :d)
@test !has_interaction(A, :a, :c)
@test has_interaction(A, :b, :c)
@test !has_interaction(A, :b, :d)

A = UnipartiteNetwork(Matrix(I, (5,5)))
X = nodiagonal(A)
@test !has_interaction(X, 1, 1)
@test !has_interaction(X, 2, 2)
@test !has_interaction(X, 3, 3)

# Accessing an interaction
A = UnipartiteNetwork([true false; true true], [:a, :b])
@test A[:a, :a]
@test A[:b, :a]
@test A[:b, :b]

# Accessing a range of interactions
A = UnipartiteNetwork(Matrix(I, (4,4)), [:a, :b, :c, :d])
B = A[[:a, :b]]
@test typeof(B) <: UnipartiteNetwork
@test richness(B) == 2

# Accessing a range of interactions (bipartite case)
A = BipartiteNetwork([true false true; false true true; true true true], [:A, :B, :C], [:a, :b, :c])

B1 = A[:, [:b, :c]]
@test typeof(B1) == typeof(A)
@test richness(B1) == 5

B2 = A[[:A, :B], :]
@test typeof(B2) == typeof(A)
@test richness(B2) == 5

B3 = A[[:A, :B], [:a, :c]]
@test typeof(B3) == typeof(A)
@test richness(B3) == 4

S1 = A[:A, :]
@test typeof(S1) <: Set
@test eltype(S1) <: Symbol

S2 = A[:, :a]
@test typeof(S2) <: Set
@test eltype(S2) <: Symbol

# Transposition
A = web_of_life("A_HP_001")
s1, s2 = size(A)
B = permutedims(A)
@test size(B) == (s2, s1)

# Get index

A = UnipartiteQuantitativeNetwork([1 2 3; 4 5 6; 7 8 9], [:a, :b, :c])
@test A[:a, :a] == 1
@test A[:a, :b] == 2
@test A[:a, :c] == 3
@test A[:b, :a] == 4
@test A[:b, :b] == 5
@test A[:b, :c] == 6
@test A[:c, :a] == 7
@test A[:c, :b] == 8
@test A[:c, :c] == 9

@test A[[:a, :b]][:a,:b] == 2
@test A[[:a, :b]][:a,:a] == 1
@test A[[:a, :b]][:b,:a] == 4
@test A[[:a, :b]][:b,:b] == 5

# Modify elements
N = UnipartiteNetwork(zeros(Bool, (4,4)))
@test links(N) == 0
N[1,1] = true
@test links(N) == 1
N["s2", "s4"] = true
@test links(N) == 2

# Test `similar`
N = UnipartiteNetwork(ones(Bool, (4,4)))
W = similar(N)
@test !W[1,1,]
@test !W["s2", "s4"]
@test links(W) == 0

N = BipartiteNetwork(zeros(Bool, (4,4)))
@test links(N) == 0
N[1,1] = true
@test links(N) == 1
N["t1", "b2"] = true
@test links(N) == 2
nodiagonal!(N)
@test links(N) == 2

# Interactions
interactions(BipartiteQuantitativeNetwork(rand(Float64, (3,5))))
interactions(BipartiteProbabilisticNetwork(rand(Float64, (3,5))))

# Test throws

@test_throws ArgumentError species(BipartiteQuantitativeNetwork(rand(Float64, (3,5))); dims=5)

# Iteration
for (i, int) in enumerate(unipartitemotifs().S1)
    if i == 1
        @test int.from == :a
        @test int.to == :b
    end
    if i == 2
        @test int.from == :b
        @test int.to == :c
    end
end

for (i, int) in enumerate(null2(unipartitemotifs().S1))
    if i == 1
        @test int.from == :a
        @test int.to == :a
        @test int.probability ≈ 1/6
    end
    if i == 2
        @test int.from == :b
        @test int.to == :a
        @test int.probability ≈ 1/6
    end
end

end