text:
# functions for calculating node coordinates and metrics for curvilinear
# meshes

#------------------------------------------------------------------------------
# allocators for curvilinear meshes

"""
  Allocates mesh.vert_coords
"""
function allocateMeshCoordinateArray(mesh::PumiMeshDG{Tmsh}, sbp::AbstractSBP) where Tmsh

  num_coord_nodes = mesh.coord_numNodesPerElement

  if !isFieldDefined(mesh, :vert_coords)
    mesh.vert_coords = Array{Tmsh}(mesh.dim, num_coord_nodes, mesh.numEl)
    mesh.vert_coords_bar = Array{Tmsh}(mesh.dim, num_coord_nodes, mesh.numEl)
  else
    fill!(mesh.vert_coords, 0.0)
  end

  return nothing
end

"""
  Allocates the following fields of mesh:
    coords
    dxidx
    jac

  Also allocates the following fields with size zero (because they are not used):
   dxidx_face
   jac_face
   dxidx_bndry
   jac_bndry
   dxidx_sharedface
   jac_sharedface
"""
function allocateCurvilinearCoordinateAndMetricArrays(mesh::PumiMeshDG{Tmsh},                                                         sbp::AbstractSBP) where Tmsh

  dim = mesh.dim
  sbpface = mesh.sbpface

  if !isFieldDefined(mesh, :coords, :dxidx, :jac)
    mesh.coords = zeros(Tmsh, mesh.dim, mesh.numNodesPerElement, mesh.numEl)
    mesh.dxidx = zeros(Tmsh, mesh.dim, mesh.dim, mesh.numNodesPerElement, mesh.numEl)
    mesh.jac = zeros(Tmsh, mesh.numNodesPerElement, mesh.numEl)

    mesh.dxidx_bar = zeros(mesh.dxidx)
    mesh.jac_bar = zeros(mesh.jac)

    # these arrays are not used for curvilinear meshes

    # interior faces
    mesh.dxidx_face = zeros(Tmsh, 0, 0, 0, 0)
    mesh.jac_face = zeros(Tmsh, 0, 0)

    mesh.dxidx_face_bar = zeros(mesh.dxidx_face)
    mesh.jac_face_bar = zeros(mesh.jac_face)

    # boundary faces
    mesh.dxidx_bndry = zeros(Tmsh, 0, 0, 0, 0)
    mesh.jac_bndry = zeros(Tmsh, 0, 0)

    mesh.dxidx_bndry_bar = zeros(mesh.dxidx_bndry)
    mesh.jac_bndry_bar = zeros(mesh.jac_bndry)

    # parallel shared faces
    mesh.dxidx_sharedface = Array{Array{Tmsh, 4}}(0)
    mesh.jac_sharedface = Array{Array{Tmsh, 2}}(0)

    mesh.dxidx_sharedface_bar = Array{Array{Tmsh, 4}}(0)
    mesh.jac_sharedface_bar = Array{Array{Tmsh, 2}}(0)
  else
    fill!(mesh.coords, 0.0)
    fill!(mesh.dxidx, 0.0)
    fill!(mesh.jac, 0.0)
  end

  return nothing
end

#------------------------------------------------------------------------------
# curvilinear functions

"""
  This function gets the coordinates that define the mesh
  (not the solution node coordinates) and puts them the mesh.vert_coords field
  of the mesh object.  For each element, the ordering of the coordinates
  is verts, then edges, then faces, then regions.  This function uses
  smart allocators to allocate the array if needed.

  Input:
    mesh: a DG mesh
    sbp: an SBP operators
"""
function getMeshCoordinates(mesh::PumiMeshDG{Tmsh}, sbp::AbstractSBP) where Tmsh

  allocateMeshCoordinateArray(mesh, sbp)

  coords_tmp = zeros(Float64, mesh.dim, mesh.coord_numNodesPerElement)
  for i=1:mesh.numEl
    el_i = mesh.elements[i]
    coords_i = sview(mesh.vert_coords, :, :, i)
    apf.getAllEntityCoords(mesh.m_ptr, el_i, coords_tmp)
    copy!(coords_i, coords_tmp)
  end

  return nothing
end

"""
  This function calculates the fields of the mesh that hold coordinates of the
  face nodes for boundaries, interfaces, and sharedfaces.  This function uses
  smart allocators to allocate the arrays if needed
"""
function getFaceCoordinatesAndNormals(mesh::PumiMeshDG{Tmsh}, sbp::AbstractSBP) where Tmsh

  allocateFaceCoordinates(mesh)

  allocateNormals(mesh, sbp)

  if length(mesh.bndryfaces) > 0  # debugging: don't call if unneeded
    calcFaceCoordinatesAndNormals(mesh, sbp, mesh.bndryfaces,
                                  mesh.coords_bndry, mesh.nrm_bndry)
  end

  calcFaceCoordinatesAndNormals(mesh, sbp, mesh.interfaces,
                               mesh.coords_interface, mesh.nrm_face)
  for i=1:mesh.npeers
    calcFaceCoordinatesAndNormals(mesh, sbp, mesh.bndries_local[i],
                               mesh.coords_sharedface[i], mesh.nrm_sharedface[i])
  end

  return nothing
end

"""
  This function back propigates mesh.nrm_*_bar and mesh.coords_*_bar to
  mesh.vert_coords_bar.  mesh.vert_coords_bar is updated (not overwritten)
  with the results.
"""
function getFaceCoordinatesAndNormals_rev(mesh::PumiMeshDG{Tmsh},
                                          sbp::AbstractSBP) where Tmsh
  if length(mesh.bndryfaces) > 0  # debugging: don't call if unneeded
    calcFaceCoordinatesAndNormals_rev(mesh, sbp, mesh.bndryfaces,
                                      mesh.coords_bndry,
                                      mesh.coords_bndry_bar,
                                      mesh.nrm_bndry,
                                      mesh.nrm_bndry_bar)
  end

  coords_face_bar = zeros(Tmsh, mesh.dim, mesh.numNodesPerFace, mesh.numInterfaces)
  calcFaceCoordinatesAndNormals_rev(mesh, sbp, mesh.interfaces,
                                mesh.coords_interface, coords_face_bar,
                                mesh.nrm_face, mesh.nrm_face_bar)

  for i=1:mesh.npeers
    coords_sharedface_bar = zeros(Tmsh, mesh.dim, mesh.numNodesPerFace, length(mesh.bndries_local[i]))
    calcFaceCoordinatesAndNormals_rev(mesh, sbp, mesh.bndries_local[i],
                                  mesh.coords_sharedface[i],
                                  coords_sharedface_bar,
                                  mesh.nrm_sharedface[i],
                                  mesh.nrm_sharedface_bar[i])
  end

  return nothing
end

"""
  This is the inner function used by getFaceCoordinatesandNormals.  It
  gets the coordinates of the face nodes and their normal vectors for
  a given array of Interfaces/Boundaries

  Inputs:
    mesh: a DG mesh object
    faces: an array of Boundaries or Interfaces

  Inputs/Outputs
    coords_face: array Tdim x numNodesPerFace x length(faces)
    nrm_face: array, same size as coords_face

  Aliasing restrictions: none, although it would be weird if coords_face and
                         nrm_face aliased
"""
function calcFaceCoordinatesAndNormals(
                    mesh::PumiMeshDG{Tmsh}, sbp::AbstractSBP,
                    faces::AbstractArray{I, 1},
                    coords_face::AbstractArray{Tmsh, 3},
                    nrm_face::AbstractArray{Tmsh, 3}) where {Tmsh, I <: Union{Boundary, Interface}}

  blocksize = 1000  # magic parameter: number of faces to do in a group
  nfaces = length(faces)

  # calculate number of blocks
  nblocks_full = div(nfaces, blocksize)
  nrem = nfaces % blocksize

  numNodesPerElement = mesh.coord_numNodesPerElement
  numNodesPerFace = mesh.coord_numNodesPerType[mesh.dim]

  # some temporary arrays
  down_faces = Array{Ptr{Void}}(12)
  coords_lag_face = Array{Tmsh}(mesh.dim, mesh.coord_numNodesPerFace, blocksize)

  # get the parametic coordinates of the face nodes
  face_xi = mesh.coord_facexi
  ref_verts = baryToXY(face_xi, mesh.sbpface.vtx)

  face_idx = 1  # index in faces array
  for block=1:nblocks_full
    start_idx = (block - 1)*blocksize + 1
    end_idx = block*blocksize
    coords_face_block = sview(coords_face, :, :, start_idx:end_idx)
    nrm_face_block = sview(nrm_face, :, :, start_idx:end_idx)

    # load up data for current block
    for i=1:blocksize
      el_i = getElementL(faces[face_idx])
      face_i = getFaceL(faces[face_idx])
      el_ptr = mesh.elements[el_i]

      coords_i = sview(coords_lag_face, :, :, i)
      getMeshFaceCoordinates(mesh, el_i, face_i, coords_i)
      face_idx += 1
    end

    # populate output array for current block
    calcFaceNormals!(mesh.sbpface, mesh.coord_order, ref_verts, coords_lag_face,
                     coords_face_block, nrm_face_block)
    fill!(coords_lag_face, 0.0)

    fixOutwardNormal(mesh, sview(faces, start_idx:end_idx), nrm_face_block)
  end  # end loop over full blocks

  # do remainder loop
  start_idx = nblocks_full*blocksize + 1
  end_idx = nfaces
  @assert end_idx - start_idx + 1 <= blocksize
  @assert end_idx - start_idx + 1 == nrem

  coords_face_block = sview(coords_face, :, :, start_idx:end_idx)
  nrm_face_block = sview(nrm_face, :, :, start_idx:end_idx)
  coords_lag_face_block = sview(coords_lag_face, :, :, 1:nrem)

  for i=1:nrem  # loop over remaining faces
    el_i = getElementL(faces[face_idx])
    face_i = getFaceL(faces[face_idx])
    el_ptr = mesh.elements[el_i]

    coords_i = sview(coords_lag_face, :, :, i)
    getMeshFaceCoordinates(mesh, el_i, face_i, coords_i)
    face_idx += 1

  end

  calcFaceNormals!(mesh.sbpface, mesh.coord_order, ref_verts,
                   coords_lag_face_block, coords_face_block, nrm_face_block)

  # make sure the normal vectors point outwards
  fixOutwardNormal(mesh, sview(faces, start_idx:end_idx), nrm_face_block)

  return nothing
end

"""
  Reverse mode of calcFaceCoordinatesAndNormals_rev, back propigates
  coords_face_bar and nrm_face_bar to mesh.vertcoords_bar.

  This function also recalculates coords_face and nrm_face in the course
  of doing the reverse mode.
"""
function calcFaceCoordinatesAndNormals_rev(
                    mesh::PumiMeshDG{Tmsh}, sbp::AbstractSBP,
                    faces::AbstractArray{I, 1},
                    coords_face::AbstractArray{Tmsh, 3},
                    coords_face_bar::AbstractArray{Tmsh, 3},
                    nrm_face::AbstractArray{Tmsh, 3},
                    nrm_face_bar::AbstractArray{Tmsh, 3}; print=false) where {Tmsh, I <: Union{Boundary, Interface}}

  blocksize = 1000  # magic parameter: number of faces to do in a group
  nfaces = length(faces)

  # calculate number of blocks
  nblocks_full = div(nfaces, blocksize)
  nrem = nfaces % blocksize

  numNodesPerElement = mesh.coord_numNodesPerElement
  numNodesPerFace = mesh.coord_numNodesPerType[mesh.dim]

  # some temporary arrays
  down_faces = Array{Ptr{Void}}(12)
  coords_lag_face = Array{Tmsh}(mesh.dim, mesh.coord_numNodesPerFace, blocksize)
  coords_lag_face_bar = zeros(coords_lag_face)

  # get the parametic coordinates of the face nodes
  face_xi = mesh.coord_facexi
  ref_verts = baryToXY(face_xi, mesh.sbpface.vtx)

  for block=1:nblocks_full
    start_idx = (block - 1)*blocksize + 1
    end_idx = block*blocksize
    faces_block = sview(faces, start_idx:end_idx)
    coords_face_block = sview(coords_face, :, :, start_idx:end_idx)
    coords_face_bar_block = sview(coords_face_bar, :, :, start_idx:end_idx)
    nrm_face_block = sview(nrm_face, :, :, start_idx:end_idx)
    nrm_face_bar_block = sview(nrm_face_bar, :, :, start_idx:end_idx)

    # load up data for current block
    for i=1:blocksize
      el_i = getElementL(faces_block[i])
      face_i = getFaceL(faces_block[i])
      el_ptr = mesh.elements[el_i]

      coords_i = sview(coords_lag_face, :, :, i)
      getMeshFaceCoordinates(mesh, el_i, face_i, coords_i)
    end

    # forward sweep
    # need the *original* face normals for fix_outward_normal, so recalculate
    # them here
    calcFaceNormals!(mesh.sbpface, mesh.coord_order, ref_verts, coords_lag_face,
                     coords_face_block, nrm_face_block)

    # reverse sweep
    fixOutwardNormal_rev(mesh, sview(faces, start_idx:end_idx), nrm_face_block,
                         nrm_face_bar_block)

    fill!(coords_lag_face_bar, 0.0)
    calcFaceNormals_rev!(mesh.sbpface, mesh.coord_order, ref_verts,
                         coords_lag_face, coords_lag_face_bar,
                         coords_face_bar_block, nrm_face_bar_block)

    #TODO: unecessary?
    fill!(coords_lag_face, 0.0)

    for i=1:blocksize
      el_i = getElementL(faces_block[i])
      face_i = getFaceL(faces_block[i])
      el_ptr = mesh.elements[el_i]

      coords_bar_i = sview(coords_lag_face_bar, :, :, i)
      getMeshFaceCoordinates_rev(mesh, el_i, face_i, coords_bar_i)
    end
  end  # end loop over full blocks

  # do remainder loop
  start_idx = nblocks_full*blocksize + 1
  end_idx = nfaces
  @assert end_idx - start_idx + 1 <= blocksize
  @assert end_idx - start_idx + 1 == nrem

  faces_block = sview(faces, start_idx:end_idx)
  coords_face_block = sview(coords_face, :, :, start_idx:end_idx)
  coords_face_bar_block = sview(coords_face_bar, :, :, start_idx:end_idx)
  nrm_face_block = sview(nrm_face, :, :, start_idx:end_idx)
  nrm_face_bar_block = sview(nrm_face_bar, :, :, start_idx:end_idx)
  coords_lag_face_block = sview(coords_lag_face, :, :, 1:nrem)
  coords_lag_face_bar_block = sview(coords_lag_face_bar, :, :, 1:nrem)

  for i=1:nrem  # loop over remaining faces
    el_i = getElementL(faces_block[i])
    face_i = getFaceL(faces_block[i])

    coords_i = sview(coords_lag_face, :, :, i)
    getMeshFaceCoordinates(mesh, el_i, face_i, coords_i)
  end

  # forward sweep
  # we need to use the *original* face normals, not the already reversed
  # ones for fixOutwardNormal_rev
  calcFaceNormals!(mesh.sbpface, mesh.coord_order, ref_verts,
                   coords_lag_face_block, coords_face_block, nrm_face_block)

  # reverse sweep
  fixOutwardNormal_rev(mesh, sview(faces, start_idx:end_idx), nrm_face_block,
                    nrm_face_bar_block)

  fill!(coords_lag_face_bar, 0.0)

  calcFaceNormals_rev!(mesh.sbpface, mesh.coord_order, ref_verts,
                   coords_lag_face_block, coords_lag_face_bar_block,
                   coords_face_bar_block, nrm_face_bar_block)

  for i=1:nrem
    el_i = getElementL(faces_block[i])
    face_i = getFaceL(faces_block[i])

    coords_bar_i = sview(coords_lag_face_bar_block, :, :, i)
    getMeshFaceCoordinates_rev(mesh, el_i, face_i, coords_bar_i)
  end

  return nothing
end

"""
  This function check to make sure each face normal vector is oriented
  outwards, and flips it if needed.

  Inputs:
    mesh
    faces: array of Boundary of Interfaces to check for normal orientation

  Inputs/Outputs:
    nrm_face: array of size mesh.dim x mesh.numNodesPerFace x length(faces)
              containing the normal vector at each node of each face.
              Updated in place to make normal vector point outward
"""
function fixOutwardNormal(mesh,
faces::AbstractArray{I, 1},
nrm_face::AbstractArray{Tmsh, 3}) where {I <: Union{Boundary, Interface}, Tmsh}

  nfaces = length(faces)
  should_flip_node = Array{Bool}(mesh.numNodesPerFace)
  for i=1:nfaces
    is_inward_normal(mesh, faces[i], sview(nrm_face, :, :, i), should_flip_node)

    for j=1:mesh.numNodesPerFace
      if should_flip_node[j]
        for p=1:mesh.dim
          nrm_face[p, j, i] = -nrm_face[p, j, i]
        end
      end
    end  # end j

  end  # end loop i

   return nothing
end

"""
  Reverse mode of fixOutwardNormal, reverses the primal normal vectors.
  On entry nrm_face should have the normal vectors as calculated by
  SBP.  On exit, they will point outwards.

  Inputs:
    mesh
    faces

  Inputs/Outputs:
    nrm_face
    nrm_face_bar: adjoint part of nrm_face
"""
function fixOutwardNormal_rev(mesh,
                          faces::AbstractArray{I, 1},
                          nrm_face::AbstractArray{Tmsh, 3},
                          nrm_face_bar::AbstractArray{Tmsh, 3}) where {Tmsh, I <: Union{Boundary, Interface}}

  nfaces = length(faces)
  should_flip_node = Array{Bool}(mesh.numNodesPerFace)
  for i=1:nfaces
    is_inward_normal(mesh, faces[i], sview(nrm_face, :, :, i), should_flip_node)

    for j=1:mesh.numNodesPerFace
      if should_flip_node[j]
        for p=1:mesh.dim
          nrm_face_bar[p, j, i] = -nrm_face_bar[p, j, i]
          nrm_face[p, j, i] = -nrm_face[p, j, i]
        end
      end
    end  # end j

  end  # end loop i

  return nothing
end

"""
  Check if the normal vector on each node of a face is pointing inward

  The algorithm determines orientation by comparing the normal vector against
  the vector along an edge of the simplex (from the vertex not on the face to
  the vertex on the face).  For highly skewed elements, using different
  vertices on the face to compute the vector can give different results, so
  only the vector with the largest magnitude dot product with the normal vector
  is considered.

  Inputs:
    mesh: the mesh
    iface: either a Boundary or an Interface
    nrm_face: mesh.dim x mesh.numNodesPerFace array containing the normal
              vectors for each face node

  Inputs/Outputs:
    should_flip_node: an Bool array of length numNodesPerFace specifying whether
                      the normal vector at each node is pointing inward
"""
function is_inward_normal(mesh, iface::Union{Boundary, Interface},
                    nrm_face::AbstractMatrix{Tmsh},
                    should_flip_node::AbstractVector{Bool}) where Tmsh

  tmp = zeros(3)  # temporary vector to hold coordinates
  topo = mesh.topo
  numVertPerElement = mesh.numTypePerElement[1]
  numVertPerFace = numVertPerElement - 1

  # temporary arrays
#  el_verts = Array{Ptr{Void}}(numVertPerElement)
  other_vert_coords = zeros(Tmsh, mesh.dim)
#  face_verts = Array{Ptr{Void}}(numVertPerElement - 1)
  face_vert_coords = zeros(Tmsh, mesh.dim, numVertPerFace)

  elnum = getElementL(iface)
  facenum_local = getFaceL(iface)

#  el_i = mesh.elements[elnum]
#  apf.getDownward(mesh.m_ptr, el_i, 0, el_verts)

  for j=1:numVertPerFace
    v_j = topo.face_verts[j, facenum_local]

    for p=1:mesh.dim
      face_vert_coords[p, j] = mesh.vert_coords[p, v_j, elnum]
    end
  end

  # get the vert not on the face
  other_vert = 0
  for j=1:numVertPerElement
     if !(j in sview(topo.face_verts, :, facenum_local))
       other_vert = j
     end
  end

  for p=1:mesh.dim
    other_vert_coords[p] = mesh.vert_coords[p, other_vert, elnum]
  end
  # check that the face normal is in the opposite direction as the
  # vectors from a vertex on the face to the vertex not on the face

  # in some cases the normal vector can be nearly orthoogonal to vertex
  # vectors, so use the one with the greatest magnitude dot product
  should_flip = false
  max_mag = 0.0  # maximum dot product
  for j=1:mesh.numNodesPerFace
    outward_count = 0  # count number of calculations that showed outward
    for k=1:numVertPerFace
      val = zero(Tmsh)
      for p=1:mesh.dim
        r1_p = other_vert_coords[p] - face_vert_coords[p, k]
        val += nrm_face[p, j]*r1_p  # accumulate dot product
      end

      if abs(val) > max_mag
        max_mag = abs(val)
        # flip if the value is greater than 0
        should_flip = real(val) > 0
        #=
        if val < 0
          should_flip = false
        else
          should_flip = true
        end  # end if
        =#

      end
    end  # end k

    should_flip_node[j] = should_flip
  end  # end loop j

  return nothing
end

"""
  This function calculates the node coordinates, scaled mapping jacobian, and
  mapping jacobian determinant for a curvilinear mesh and stores them to
  the fields of the mesh object.  This function uses smart allocators to
  allocate the arrays if needed

  Inputs:
    mesh: a DG mesh object
    sbp: an SBP operator
"""
function getCurvilinearCoordinatesAndMetrics(mesh::PumiMeshDG{Tmsh},
                                            sbp::AbstractSBP) where Tmsh

  allocateCurvilinearCoordinateAndMetricArrays(mesh, sbp)

  ref_vtx = baryToXY(mesh.coord_xi, sbp.vtx)
#  if mesh.dim == 2
#    calcMappingJacobian!(sbp, mesh.coord_order, ref_vtx, mesh.vert_coords,
#                         mesh.coords, mesh.dxidx, mesh.jac)
#  else  # need to calculate Eone

    # block format
    blocksize = 1000  # number of elements per block
    nblocks_full = div(mesh.numEl, blocksize)
    nrem = mesh.numEl % blocksize

    Eone = zeros(Tmsh, mesh.numNodesPerElement, mesh.dim, blocksize)

    for block=1:nblocks_full
      start_idx = (block - 1)*blocksize + 1
      end_idx = block*blocksize
      element_range = start_idx:end_idx

      getCurvilinearMetricsAndCoordinates_inner(mesh, sbp, element_range, Eone)

    end  # end loop over blocks

    if nrem != 0
      # remainder block
      start_idx = nblocks_full*blocksize + 1
      end_idx = mesh.numEl
      @assert end_idx - start_idx + 1 <= blocksize
      @assert end_idx - start_idx + 1 == nrem

      element_range = start_idx:end_idx

      # make sure dimensions of Eone is correct (number of elements might be
      # less than before)
      Eone_rem = sview(Eone, :, :, 1:nrem)

      getCurvilinearMetricsAndCoordinates_inner(mesh, sbp, element_range, Eone_rem)
    end
#  end  # end if dim == 2

  return nothing
end

"""
  Back propigate dxidx to mesh.vert_coords and mesh.nrm.  See the primal method
  for details.

  Inputs:
    mesh: vert_coords_bar, nrm_face_bar, nrm_bndry_bar, nrm_sharedface are
          updated with the results
    sbp:
"""
function getCurvilinearCoordinatesAndMetrics_rev(mesh::PumiMeshDG{Tmsh},
                                                 sbp::AbstractSBP) where Tmsh

  # we have to compute E1_bar for both 2d and 3D
  blocksize = 1000  # number of elements per block
  nblocks_full = div(mesh.numEl, blocksize)
  nrem = mesh.numEl % blocksize

  Eone_bar = zeros(Tmsh, mesh.numNodesPerElement, mesh.dim, blocksize)

  for block=1:nblocks_full
    start_idx = (block - 1)*blocksize + 1
    end_idx = block*blocksize
    element_range = start_idx:end_idx

    fill!(Eone_bar, 0.0)
    getCurvilinearMetricsAndCoordinates_inner_rev(mesh, sbp, element_range,
                                                  Eone_bar)
  end

 if nrem != 0
    # remainder block
    start_idx = nblocks_full*blocksize + 1
    end_idx = mesh.numEl
    @assert end_idx - start_idx + 1 <= blocksize
    @assert end_idx - start_idx + 1 == nrem

    element_range = start_idx:end_idx

    # make sure dimensions of Eone is correct (number of elements might be
    # less than before)
    Eone_bar_rem = sview(Eone_bar, :, :, 1:nrem)

    fill!(Eone_bar_rem, 0.0)
    getCurvilinearMetricsAndCoordinates_inner_rev(mesh, sbp, element_range,
                                                  Eone_bar_rem)
  end

  return Eone_bar
end

"""
  This function calculates the metrics and coordinates for one block of
  elements.  Used by getCurvilinearMetricsAndCoordinates
"""
function getCurvilinearMetricsAndCoordinates_inner(mesh, sbp,
                             element_range::UnitRange, Eone::AbstractArray{T, 3}) where T

  # calculate Eone for current range

  vert_coords_block = sview(mesh.vert_coords, :, :, element_range)
  coords_block = sview(mesh.coords, :, :, element_range)
  dxidx_block = sview(mesh.dxidx, :, :, :, element_range)
  jac_block = sview(mesh.jac, :, element_range)
  ref_vtx = baryToXY(mesh.coord_xi, sbp.vtx)

  calcEone(mesh, sbp, element_range, Eone)

  #=
  println("checking Eone")
  for el=1:size(Eone, 3)
    println("el = ", el)
    println("Eone = \n", Eone[:, :, el])
    for d=1:mesh.dim
      println("  d = ", d)
      val = abs(sum(Eone[:, d, el]))
      println("  val = ", val)
      @assert val < 1e-6
    end
  end
  =#

#=
  println("sbp.numnodes = ", sbp.numnodes)
  println("size(coords) = ", size(coords_block))
  println("size(dxidx) = ", size(dxidx_block))
  println("size(jac) = ", size(jac_block))
  println("size(Eone) = ", size(Eone))
=#
  calcMappingJacobian!(sbp, mesh.coord_order, ref_vtx, vert_coords_block,
                       coords_block, dxidx_block, jac_block, Eone)

  fill!(Eone, 0.0)
  return nothing
end

"""
  This function calculates the metrics and coordinates for one block of
  elements.  Used by getCurvilinearMetricsAndCoordinates
"""
function getCurvilinearMetricsAndCoordinates_inner_rev(mesh, sbp,
                             element_range::UnitRange, Eone_bar::AbstractArray{T, 3}) where T

  vert_coords_block = sview(mesh.vert_coords, :, :, element_range)
  # we don't allow perturbing mesh coordinate directly (yet)
  coords_bar_block = zeros(T, mesh.dim, mesh.numNodesPerElement, length(element_range))
  dxidx_block = sview(mesh.dxidx, :, :, :, element_range)
  dxidx_bar_block = sview(mesh.dxidx_bar, :, :, :, element_range)
#  jac_block = sview(mesh.jac, :, element_range)
  jac_bar_block = sview(mesh.jac_bar, :, element_range)
  #@assert vecnorm(jac_bar_block) < 1e-13  # this is currently broken in SBP

  # outputs
  vert_coords_bar_block = sview(mesh.vert_coords_bar, :, :, element_range)
  fill!(Eone_bar, 0.0)

  ref_vtx = baryToXY(mesh.coord_xi, sbp.vtx)

  # back propigate dxidx to vert_coords, E1
  calcMappingJacobian_rev!(sbp, mesh.coord_order, ref_vtx, vert_coords_block,
                           vert_coords_bar_block, coords_bar_block,
                           dxidx_bar_block, jac_bar_block, Eone_bar)

  # back propigate E1 to the face normals
  calcEone_rev(mesh, sbp, element_range, Eone_bar)

  return nothing
end

function calcEone(mesh::PumiMeshDG{Tmsh}, sbp, element_range,
                  Eone::AbstractArray{Tmsh, 3}) where Tmsh

  # search interfaces and bndryfaces

  first_el = element_range[1]
  last_el = element_range[end]

  offset = first_el - 1

  # R times vector of ones
  # the permutation doesn't matter because it is being multiplied by a constant
  # vector.
  # also R1 = 1 by definition
  Rone = ones(Float64, mesh.numNodesPerFace)
#  Rone = vec(sum(mesh.sbpface.interp.', 2))
  sbpface = mesh.sbpface
  tmp = zeros(Tmsh, length(Rone))
  nrmL = zeros(Tmsh, mesh.dim, sbpface.numnodes)

  # accumulate E1 for a given element
  Eone_el = zeros(Tmsh, size(sbpface.perm, 1), mesh.dim)

  for i=1:mesh.numInterfaces
    iface_i = mesh.interfaces[i]

    if iface_i.elementL >= first_el && iface_i.elementL <= last_el
      elnum = iface_i.elementL
      facenum_local = iface_i.faceL
      nrm = sview(mesh.nrm_face, :, :, i)

      # call inner functions
      calcEoneElement(sbpface, nrm, Rone, tmp, Eone_el)
      assembleEone(sbpface, elnum - offset, facenum_local, Eone_el, Eone)

    end

    if iface_i.elementR >= first_el && iface_i.elementR <= last_el
      # do same as above, negating and permuting nrm
      elnum = iface_i.elementR
      facenum_local = iface_i.faceR
      orient = iface_i.orient
      for j=1:sbpface.numnodes
        for d=1:mesh.dim
          nrmL[d, j] = -mesh.nrm_face[d, sbpface.nbrperm[j, orient], i]
        end
      end

      calcEoneElement(sbpface, nrmL, Rone, tmp, Eone_el)
      assembleEone(sbpface, elnum - offset, facenum_local, Eone_el, Eone)
    end  # end if/else

  end  # end loop over interfaces

  for i=1:mesh.numBoundaryFaces
    bface_i = mesh.bndryfaces[i]

    if bface_i.element >= first_el && bface_i.element <= last_el
      elnum = bface_i.element
      facenum_local = bface_i.face
      nrm = sview(mesh.nrm_bndry, :, :, i)

      # call inner functions
      calcEoneElement(sbpface, nrm, Rone, tmp, Eone_el)
      assembleEone(sbpface, elnum - offset, facenum_local, Eone_el, Eone)
    end
  end  # end loop over boundary faces

  # check shared faces
  for peer=1:mesh.npeers
    bndryfaces_peer = mesh.bndries_local[peer]
    nrm_peer = mesh.nrm_sharedface[peer]
    for i=1:mesh.peer_face_counts[peer]
      bface_i = bndryfaces_peer[i]

      if bface_i.element >= first_el && bface_i.element <= last_el
        elnum = bface_i.element
        facenum_local = bface_i.face
        nrm = sview(nrm_peer, :, :, i)

        calcEoneElement(sbpface, nrm, Rone, tmp, Eone_el)
        assembleEone(sbpface, elnum - offset, facenum_local, Eone_el, Eone)
      end
    end
  end

  return nothing
end

"""
  Back propigates Eone_bar to the various normal vectors stored in the
  mesh objects (bndry, face, sharedface)

  Inputs:
    mesh: mesh.nrm_face_bar, mesh.nrm_bndry_bar, mesh.nrm_sharedface_bar are
          updated
    sbp
    element range: range of elements to compute
    Eone_bar: the adjoint part of Eone, see the primal method for details
"""
function calcEone_rev(mesh::PumiMeshDG{Tmsh}, sbp, element_range,
                      Eone_bar::AbstractArray{Tmsh, 3}) where Tmsh

  # search interfaces and bndryfaces

  first_el = element_range[1]
  last_el = element_range[end]

  offset = first_el - 1

  # R times vector of ones
  # the permutation doesn't matter because it is being multiplied by a constant
  # vector.
  # also R1 = 1 by def.inition
  Rone = ones(Float64, mesh.numNodesPerFace)
#  Rone = vec(sum(mesh.sbpface.interp.', 2))
  sbpface = mesh.sbpface
  tmp = zeros(Rone)
  nrmL_bar = zeros(Tmsh, mesh.dim, sbpface.numnodes)

  # accumulate E1 for a given element
  Eone_el_bar = zeros(Tmsh, size(sbpface.perm, 1), mesh.dim)

  for i=1:mesh.numInterfaces
    iface_i = mesh.interfaces[i]

    if iface_i.elementL >= first_el && iface_i.elementL <= last_el
      elnum = iface_i.elementL
      facenum_local = iface_i.faceL
      nrm_bar = sview(mesh.nrm_face_bar, :, :, i)

      # call inner functions
      fill!(Eone_el_bar, 0.0)
      assembleEone_rev(sbpface, elnum - offset, facenum_local, Eone_el_bar,
                       Eone_bar)
      calcEoneElement_rev(sbpface, nrm_bar, Rone, tmp, Eone_el_bar)

    end

    if iface_i.elementR >= first_el && iface_i.elementR <= last_el
      # do same as above, negating and permuting nrm
      elnum = iface_i.elementR
      facenum_local = iface_i.faceR
      orient = iface_i.orient

      fill!(Eone_el_bar, 0.0)
      fill!(nrmL_bar, 0.0)
      assembleEone_rev(sbpface, elnum - offset, facenum_local, Eone_el_bar,
                       Eone_bar)
      calcEoneElement_rev(sbpface, nrmL_bar, Rone, tmp, Eone_el_bar)

      for j=1:sbpface.numnodes
        for d=1:mesh.dim
          mesh.nrm_face_bar[d, sbpface.nbrperm[j, orient], i] -= nrmL_bar[d, j]
        end
      end

    end  # end if/else

  end  # end loop over interfaces

  for i=1:mesh.numBoundaryFaces
    bface_i = mesh.bndryfaces[i]

    if bface_i.element >= first_el && bface_i.element <= last_el
      elnum = bface_i.element
      facenum_local = bface_i.face
      nrm_bar = sview(mesh.nrm_bndry_bar, :, :, i)

      # call inner functions
      fill!(Eone_el_bar, 0.0)
      assembleEone_rev(sbpface, elnum - offset, facenum_local, Eone_el_bar,
                       Eone_bar)
      calcEoneElement_rev(sbpface, nrm_bar, Rone, tmp, Eone_el_bar)
    end
  end  # end loop over boundary faces

  # check shared faces
  for peer=1:mesh.npeers
    bndryfaces_peer = mesh.bndries_local[peer]
    nrm_peer_bar = mesh.nrm_sharedface_bar[peer]
    for i=1:mesh.peer_face_counts[peer]
      bface_i = bndryfaces_peer[i]

      if bface_i.element >= first_el && bface_i.element <= last_el
        elnum = bface_i.element
        facenum_local = bface_i.face
        nrm_bar = sview(nrm_peer_bar, :, :, i)

        fill!(Eone_el_bar, 0.0)
        assembleEone_rev(sbpface, elnum - offset, facenum_local, Eone_el_bar,
                         Eone_bar)
        calcEoneElement_rev(sbpface, nrm_bar, Rone, tmp, Eone_el_bar)
      end
    end
  end

  return nothing
end

"""
  Calculates E1 for a given interface.  Used by calcEone.

  Actually what it computes is:
      (R^T)*N*B*R*P*1

  Note that there should be a factor of P^T on the left.  That factor is
  applied in assembleEone().

  Inputs:
    sbpface: an AbstractFace
    nrm: the normal vectors for each node of the face, dim x numFaceNodes
    Rone: the interpolation operator R times the vector of ones
    tmp: a temporary vector of length numFaceNodes, overwritten

  Inputs/Outputs:
    Eone_el: sbp.stencilsize x dim matrix to be populated with the E
             contribution for this face.  Overwritten.

"""
function calcEoneElement(sbpface::AbstractFace, nrm::AbstractMatrix,
                         Rone::AbstractVector, tmp::AbstractVector,
                         Eone_el::AbstractMatrix)

  dim = size(Eone_el, 2)
  numFaceNodes = length(Rone)
  for d=1:dim
    for i=1:numFaceNodes
      tmp[i] = Rone[i]*nrm[d, i]*sbpface.wface[i]
    end

    Eone_dim = sview(Eone_el, :, d)
    smallmatvec!(sbpface.interp, tmp, Eone_dim)
  end

  return nothing
end

"""
  This methods works for SparseFaces (ex. diagonlE operators).
  It takes in all the same arguments as the other method for compatability
  but does not need them
"""
function calcEoneElement(sbpface::SparseFace, nrm::AbstractMatrix,
                         Rone::AbstractVector, tmp::AbstractVector,
                         Eone_el::AbstractMatrix)

  dim = size(Eone_el, 2)
  numFaceNodes = length(Rone)
  for d=1:dim
    for i=1:numFaceNodes
      Eone_el[i, d] = Rone[i]*nrm[d, i]*sbpface.wface[i]
    end
  end

  return nothing
end

# back propigate Eone_el_bar to nrm_bar
"""
  This function uses reverse mode to back-propigate Eone_el_bar to nrm_bar

  Input:
    sbpface: an AbstractFace
    Rone: the interpolation operator R times a vector of ones
    tmp: a temporary vector, overwritten
    Eone_el_bar: the adjoint part of E1

  Inputs/Outputs:
    nrm_bar: the adjoint part of the scaled face normal vector in x-y space
             updated with the contribution from this function

"""
function calcEoneElement_rev(sbpface::AbstractFace,
                         nrm_bar::AbstractMatrix,
                         Rone::AbstractVector, tmp_bar::AbstractVector,
                         Eone_el_bar::AbstractMatrix)

  dim = size(Eone_el_bar, 2)
  numFaceNodes = length(Rone)

  for d=1:dim
    # this function is linear, so no need for a forward sweep

    # reverse sweep
    Eone_dim_bar = sview(Eone_el_bar, :, d)
    # only back propigate to tmp_bar (sbpface.interp is unimportant)
    fill!(tmp_bar, 0.0)
    smallmatvec_revv!(sbpface.interp, tmp_bar, Eone_dim_bar)

    for i=1:numFaceNodes
      nrm_bar[d, i] += tmp_bar[i]*Rone[i]*sbpface.wface[i]
    end
  end

  return nothing
end

"""
  This method works for SparseFace sbp face objects.
"""
function calcEoneElement_rev(sbpface::SparseFace,
                         nrm_bar::AbstractMatrix,
                         Rone::AbstractVector, tmp_bar::AbstractVector,
                         Eone_el_bar::AbstractMatrix)

  dim = size(Eone_el_bar, 2)
  numFaceNodes = length(Rone)

  for d=1:dim
    for i=1:numFaceNodes
      nrm_bar[d, i] += Eone_el_bar[i, d]*Rone[i]*sbpface.wface[i]
    end
  end

  return nothing
end

"""
  Takes an Eone_el matrix from calcEoneElement and assemble it into the big
  Eone.

  Inputs:
    sbpface: an SBP face
    elnum: the element index of Eone to put the values into.  Note that if
           Eone is for the entire mesh, then this is the element number.  If
           Eone is for a range of elements, then elnum is the index of the
           element within the range.
    Eone_el: the E1 contribution of an element

  Inputs/Output
    Eone: a numNodesPerElement x dim x blocksize array to store E1 for a
          range of elements, updated with the new contribution
"""
function assembleEone(sbpface::AbstractFace, elnum::Integer,
                facenum_local::Integer, Eone_el::AbstractMatrix{Tmsh},
                Eone::AbstractArray{Tmsh, 3}) where Tmsh

  dim = size(Eone, 2)
  for d=1:dim
    for i=1:size(sbpface.perm, 1)
      p_i = sbpface.perm[i, facenum_local]
      Eone[p_i, d, elnum] += Eone_el[i, d]
    end
  end

  return nothing
end

"""
  Back propigates Eone_bar to Eone_el_bar.

  Inputs:
    sbpface: an AbstractFace
    elnum: the element number
    facenum_local: the local number of the face that Eone_el is calculated for
    Eone_bar: the adjoint part

  Inputs/Outputs:
    Eone_el_bar: the adjoint part
"""
function assembleEone_rev(sbpface::AbstractFace, elnum::Integer,
                      facenum_local::Integer,
                      Eone_el_bar::AbstractMatrix{Tmsh},
                      Eone_bar::AbstractArray{Tmsh, 3}) where Tmsh

  dim = size(Eone_bar, 2)
  for d=1:dim  # TODO: switch loops: turn an indexed store into a strided store
    for i=1:size(sbpface.perm, 1)
      p_i = sbpface.perm[i, facenum_local]
      # reverse mode step
      Eone_el_bar[i, d] += Eone_bar[p_i, d, elnum]
    end
  end

  return nothing
end

# get the coordinates of the nodes on a given face (in the coordinate field, not
# the solution field)
# this only works up to second order, and uses the ElementTopology to figure
# out the orientation of the edge
# elnum is the global element number
# facenum is the local face number
"""
  This function gets the coordinates of the coordinate field nodes of the
  face of an element, in the orientation specified by mesh.topo.  This only
  works up to second order.  It uses the SBP element topology to determine
  the ordering of the nodes.  I guess this makes sense because calcFaceNormals
  is an SBP function and requires its definition of the vertices that define the
  face.  Because ref_vtx is passed into calcFaceNormals, it should be ok to
  use the Pumi definition of the coordinates of the reference nodes (precisely
  because we use the SBP definition of which vertices define the face, the
  face node coordinates are relative to those vertices).

  Inputs:
    mesh: a 2D mesh
    elnum: the global element number
    facenum: the local face number

  Inputs/Outputs:
    coords: array to be populated with the coordinates, Tdim x
            number of coordinate nodes on a face (2 for linear 2D,
            3 for quadratic 2D)
"""
function getMeshFaceCoordinates(mesh::PumiMesh2DG, elnum::Integer,
                                facenum::Integer,
                                coords::AbstractMatrix)

# coords = a 2 x number of nodes on the face array to be populated with the
# coordinates

  topo = mesh.topo

  # equivalent to topo.face_verts[:, facenum]
  v1_idx = topo.face_verts[1, facenum]
  v2_idx = topo.face_verts[2, facenum]

  coords[1, 1] = mesh.vert_coords[1, v1_idx, elnum]
  coords[2, 1] = mesh.vert_coords[2, v1_idx, elnum]

  coords[1, 2] = mesh.vert_coords[1, v2_idx, elnum]
  coords[2, 2] = mesh.vert_coords[2, v2_idx, elnum]

  if mesh.coord_numNodesPerType[2] > 0
    edge_idx = 3 + topo.face_edges[1, facenum]
    coords[1, 3] = mesh.vert_coords[1, edge_idx, elnum]
    coords[2, 3] = mesh.vert_coords[2, edge_idx, elnum]
  end

  return nothing
end

"""
  This function is the reverse mode of getMeshFaceCoordinates
  This function stores the adjoint part of the vert_coords to mesh.vert_coords

  2nd order coordinate fields only.

  This function *accumulates* into mesh.vertcoords_bar.

  Inputs
    mesh::PumiMesh2DG
    elnum: the element number
    facenum: the local face number
    coords_bar:  the adjoint part of the coordinate field for the given face,
                 2 x coord_numNodesPerFace.
                 The data should be ordered vertices, then mid edge nodes.
"""
function getMeshFaceCoordinates_rev(mesh::PumiMesh2DG, elnum::Integer,
                                facenum::Integer,
                                coords_bar::AbstractMatrix)

#  vertmap = mesh.topo_pumi.face_verts
  topo = mesh.topo_pumi

  # equivalent to topo.face_verts[:, facenum]
  v1_idx = topo.face_verts[1, facenum]
  v2_idx = topo.face_verts[2, facenum]

  mesh.vert_coords_bar[1, v1_idx, elnum] += coords_bar[1, 1]
  mesh.vert_coords_bar[2, v1_idx, elnum] += coords_bar[2, 1]

  mesh.vert_coords_bar[1, v2_idx, elnum] += coords_bar[1, 2]
  mesh.vert_coords_bar[2, v2_idx, elnum] += coords_bar[2, 2]

  if mesh.coord_numNodesPerType[2] > 0
    edge_idx = facenum + 3  # for simplex elements
    mesh.vert_coords_bar[1, edge_idx, elnum] += coords_bar[1, 3]
    mesh.vert_coords_bar[2, edge_idx, elnum] += coords_bar[2, 3]
  end

  return nothing
end

function getMeshFaceCoordinates(mesh::PumiMesh3DG, elnum::Integer,
                                facenum::Integer, coords::AbstractMatrix)

  topo = mesh.topo
  # get the face vertices
  face_vert_idx = sview(topo.face_verts, :, facenum)

  for i=1:3  # 3 vertices per face
    for j=1:3  # 3 coordinates at each vertex
      coords[j, i] = mesh.vert_coords[j, face_vert_idx[i], elnum]
    end
  end

  if mesh.coord_numNodesPerType[2] > 0
    for i=1:3  # 3 edges per face
      edge_idx = topo.face_edges[i, facenum] + 4 # 4 = numVerts
      for j=1:3  # 3 coordinates per face
        coords[j, 3 + i] = mesh.vert_coords[j, edge_idx, elnum]
      end
    end

  end

  return nothing
end

"""
  Similar to the 2d method.

  coords_bar should be ordered vertices, edge nodes, then face nodes
"""
function getMeshFaceCoordinates_rev(mesh::PumiMesh3DG, elnum::Integer,
                                facenum::Integer,
                                coords_bar::AbstractMatrix)

#  vertmap = mesh.topo.face_verts
  topo = mesh.topo

  # get the face vertices
  face_vert_idx = sview(topo.face_verts, :, facenum)

  for i=1:3  # 3 vertices per face
    for j=1:3  # 3 coordinates at each vertex
      mesh.vert_coords_bar[j, face_vert_idx[i], elnum] += coords_bar[j, i]
    end
  end

  if mesh.coord_numNodesPerType[2] > 0
    offset = 3 + (facenum - 1)*3  # 3 vertices + 1 edge node per face

    for i=1:3  # 3 edges per face
      edge_idx = topo.face_edges[i, facenum] + 4  # 4 = numVerts
      for j=1:3  # 3 coordinates per face
        mesh.vert_coords_bar[j, edge_idx, elnum] += coords_bar[j, 3 + i]  # 3 verts
      end
    end

  end

  return nothing
end


text:
@testset "Defining and evaluating Gaussian SPNs" begin

    import SPNetworks: GaussianDistribution

    GSPN = SumProductNetwork([
        SumNode([2,3,4],[4/20,9/20,7/20]),# 1 Sum P(A, B) 0.2,0.5,0.3
        ProductNode([5,7]),               # 2 Prod P(A, B)=P(A)*P(B)
        ProductNode([5,8]),               # 3 Prod P(A,B)=P(A)*P(B)
        ProductNode([6,8]),               # 4 Prod P(A,B)=P(A)*P(B)
        GaussianDistribution(1, 2, 18),   # 5 Normal(A, mean=2, var=18)
        GaussianDistribution(1, 11, 8),   # 6 Normal(A, mean=11, var=8)
        GaussianDistribution(2, 3, 10),   # 7 Normal(B, mean=3, var=10)
        GaussianDistribution(2, -4, 7)    # 8 Normal(B, mean=-4, var=7)
    ])

    @testset "Evaluation" begin
        res = GSPN([11.0,-4.0])
        @test res ≈ 0.008137858167261642
        @test logpdf(GSPN,[11.0,-4.0]) ≈ -4.811228258006023
    end
    @testset "Sampling" begin
        #N = 500000
        N = 1000
        data = rand(GSPN,N)

        max = -Inf
        amax = 0
        for n=1:N
            v = logpdf(GSPN,view(data,n,:))  # GSPN(data[n,:])
            if v > max
                max = v
                amax = n
            end
        end
        a, b = 11.0, -4.0
        ref = logpdf(GSPN,[a,b]) #GSPN([a,b])
        # println("max ln S($(data[amax,:])) = $max ≈ ln S($a,$b) = $ref")
    end
end # END of test set

ext: jl


text:
function emulate(fname::String)
    myChip8 = Chip()
    loadApplication(myChip8, fname)
end


text:
module problem132

description = """
Large repunit factors

A number consisting entirely of ones is called a repunit. We shall define R(k) to be a repunit of length k.
For example, R(10) = 1111111111 = 11×41×271×9091, and the sum of these prime factors is 9414.

Find the sum of the first forty prime factors of R(10^9).
"""

# Optimized search
# R(n) = (10^n - 1)/9
# R(n) = 0 mod p
# (10^n - 1) / 9 = 0 mod p
# 10^n - 1 = 0 mod 9p
# 10^n = 1 mod 9p
# so search for p satisfying 10^n = 1 mod 9p using Julias powermod function
function search(n = 10^9, facs = 40, maxprime = 200000)
  factors = Int[]
  for p in primes(maxprime)
    if powermod(10, n, 9*p) == 1
      push!(factors, p)
      if length(factors) == facs
        return sum(factors)
      end
    end
  end
end

using Base.Test

@test search(10, 4) == sum([11, 41, 271, 9091])

end

max_stars_repo_name: mbuhot/mbuhot-euler-solutions

max_issues_repo_name: mbuhot/mbuhot-euler-solutions

max_forks_repo_name: mbuhot/mbuhot-euler-solutions


text:
using Test
using Dates

using TOMLX: TOMLX, parsestring, tryparsestring, ParserError, Internals

include("readme.jl")
include("toml_test.jl")
include("values.jl")
include("invalids.jl")
include("error_printing.jl")


text:
module GenDirectionalStats

include("utils.jl")
include("vonmisesfisher.jl")
include("rot3.jl")
include("rot2.jl")
include("direction3.jl")
include("hopf.jl")
include("moves.jl")

end # module


text:
using ImportMacros

@import ChargedParticleDynamics.GuidingCenter4d.TokamakFastDeeplyPassing as problem

# run configuration
const Δt    = 2.5
const ntime = 1250000

# plot configuration
const nplot = 1
const splot = true
const XLIM  = ( 0.80,  2.70)
const YLIM  = (-0.90, +0.90)

# output configuration
const coll_title = "Guiding Center 4D Fast Deeply Passing"


text:
# to test vtk-files
using StableRNGs
OVERWRITE_CHECKSUMS = false
checksums_file = joinpath(dirname(@__FILE__), "checksums.sha1")
checksum_list = read(checksums_file, String)
if OVERWRITE_CHECKSUMS
    csio = open(checksums_file, "w")
else
    csio = open(checksums_file, "r")
end

@testset "Grid, DofHandler, vtk" begin
    for (celltype, dim) in ((Line,                   1),
                            (QuadraticLine,          1),
                            (Quadrilateral,          2),
                            (QuadraticQuadrilateral, 2),
                            (Triangle,               2),
                            (QuadraticTriangle,      2),
                            (Hexahedron,             3),
                            (Cell{3,20,6},           3),
                            (Tetrahedron,            3))

        # create test grid, do some operations on it and then test
        # the resulting sha1 of the stored vtk file
        # after manually checking the exported vtk
        nels = ntuple(x->5, dim)
        right = Vec{dim, Float64}(ntuple(x->1.5, dim))
        left = -right
        grid = generate_grid(celltype, nels, left, right)

        transform!(grid, x-> 2x)

        radius = 2*1.5
        addcellset!(grid, "cell-1", [1,])
        addcellset!(grid, "middle-cells", x -> norm(x) < radius)
        addfaceset!(grid, "middle-faceset", x -> norm(x) < radius)
        addfaceset!(grid, "right-faceset", getfaceset(grid, "right"))
        addnodeset!(grid, "middle-nodes", x -> norm(x) < radius)

        gridfilename = "grid-$(Ferrite.celltypes[celltype])"
        vtk_grid(gridfilename, grid) do vtk
            vtk_cellset(vtk, grid, "cell-1")
            vtk_cellset(vtk, grid, "middle-cells")
            vtk_nodeset(vtk, grid, "middle-nodes")
        end

        # test the sha of the file
        sha = bytes2hex(open(SHA.sha1, gridfilename*".vtu"))
        if OVERWRITE_CHECKSUMS
            write(csio, sha, "\n")
        else
            @test chomp(readline(csio)) == sha
            rm(gridfilename*".vtu")
        end

        # Create a DofHandler, add some things, write to file and
        # then check the resulting sha
        dofhandler = DofHandler(grid)
        push!(dofhandler, :temperature, 1)
        push!(dofhandler, :displacement, 3)
        close!(dofhandler)
        ch = ConstraintHandler(dofhandler)
        dbc = Dirichlet(:temperature, union(getfaceset(grid, "left"), getfaceset(grid, "right-faceset")), (x,t)->1)
        add!(ch, dbc)
        dbc = Dirichlet(:temperature, getfaceset(grid, "middle-faceset"), (x,t)->4)
        add!(ch, dbc)
        for d in 1:dim
            dbc = Dirichlet(:displacement, union(getfaceset(grid, "left")), (x,t) -> d, d)
            add!(ch, dbc)
        end
        close!(ch)
        update!(ch, 0.0)
        rng = StableRNG(1234)
        u = rand(rng, ndofs(dofhandler))
        apply!(u, ch)

        dofhandlerfilename = "dofhandler-$(Ferrite.celltypes[celltype])"
        vtk_grid(dofhandlerfilename, dofhandler) do vtk
            vtk_point_data(vtk, ch)
            vtk_point_data(vtk, dofhandler, u)
        end

        # test the sha of the file
        sha = bytes2hex(open(SHA.sha1, dofhandlerfilename*".vtu"))
        if OVERWRITE_CHECKSUMS
            write(csio, sha, "\n")
        else
            @test chomp(readline(csio)) == sha
            rm(dofhandlerfilename*".vtu")
        end

    end

end # of testset

close(csio)

@testset "vtk tensor export" begin
    # open files
    checksums_file_tensors = joinpath(dirname(@__FILE__), "checksums2.sha1")
    if OVERWRITE_CHECKSUMS
        csio = open(checksums_file_tensors, "w")
    else
        csio = open(checksums_file_tensors, "r")
    end

    # 3D grid
    grid = generate_grid(Hexahedron, (1,1,1))

    sym_tensor_data = [SymmetricTensor{2,3}(ntuple(i->i, 6)) for j=1.0:8.0]
    tensor_data = [Tensor{2,3}(ntuple(i->i, 9)) for j=1.0:8.0]
    vector_data = [Vec{3}(ntuple(i->i, 3)) for j=1:8]

    filename_3d = "test_vtk_3d"
    vtk_grid(filename_3d, grid) do vtk_file
        vtk_point_data(vtk_file, sym_tensor_data, "symmetric tensor")
        vtk_point_data(vtk_file, tensor_data, "tensor")
        vtk_point_data(vtk_file, vector_data, "vector")
    end

    # 2D grid
    grid = generate_grid(Quadrilateral, (1,1))

    sym_tensor_data = [SymmetricTensor{2,2}(ntuple(i->i, 3)) for j=1.0:4.0]
    tensor_data = [Tensor{2,2}(ntuple(i->i, 4)) for j=1.0:4.0]
    tensor_data_1D = [SymmetricTensor{2,1}(ntuple(i->i, 1)) for j=1.0:4.0]
    vector_data = [Vec{2}(ntuple(i->i, 2)) for j=1:4]

    filename_2d = "test_vtk_2d"
    vtk_grid(filename_2d, grid) do vtk_file
        vtk_point_data(vtk_file, sym_tensor_data, "symmetric tensor")
        vtk_point_data(vtk_file, tensor_data, "tensor")
        vtk_point_data(vtk_file, tensor_data_1D, "tensor_1d")
        vtk_point_data(vtk_file, vector_data, "vector")
    end

    # test the shas of the files
    files = [filename_3d, filename_2d]
    for filename in files
        sha = bytes2hex(open(SHA.sha1, filename*".vtu"))
        if OVERWRITE_CHECKSUMS
            write(csio, sha, "\n")
        else
            @test chomp(readline(csio)) == sha
            rm(filename*".vtu")
        end
    end

    close(csio)
end

@testset "Grid utils" begin

    grid = Ferrite.generate_grid(QuadraticQuadrilateral, (1, 1), Vec((0.,0.)), Vec((1.,1.)))

    addcellset!(grid, "cell_set", [1]);
    node_set = Set(1:getnnodes(grid))
    addnodeset!(grid, "node_set", node_set)

    @test getnodesets(grid) == Dict("node_set" => node_set)

    @test getnodes(grid, [1]) == [getnodes(grid, 1)] # untested

    @test length(getnodes(grid, "node_set")) == 9

    @test collect(getcoordinates(getnodes(grid, 5)).data) ≈ [0.5, 0.5]

    @test getcells(grid, "cell_set") == [getcells(grid, 1)]

    f(x) = Tensor{1,1,Float64}((1 + x[1]^2 + 2x[2]^2, ))

    values = compute_vertex_values(grid, f)
    @test f([0.0, 0.0]) == values[1]
    @test f([0.5, 0.5]) == values[5]
    @test f([1.0, 1.0]) == values[9]

    @test compute_vertex_values(grid, collect(1:9), f) == values

    # Can we test this in a better way? The set makes the order random.
    @test length(compute_vertex_values(grid, "node_set", f)) == 9

    # CellIterator on a grid without DofHandler
    grid = generate_grid(Triangle, (4,4))
    n = 0
    ci = CellIterator(grid)
    @test length(ci) == getncells(grid)
    for c in ci
        getcoordinates(c)
        getnodes(c)
        n += cellid(c)
    end
    @test n == div(getncells(grid)*(getncells(grid) + 1), 2)
end

@testset "Grid sets" begin

    grid = Ferrite.generate_grid(Hexahedron, (1, 1, 1), Vec((0.,0., 0.)), Vec((1.,1.,1.)))

    #Test manual add
    addcellset!(grid, "cell_set", [1]);
    addnodeset!(grid, "node_set", [1])
    addfaceset!(grid, "face_set", [FaceIndex(1,1)])
    addedgeset!(grid, "edge_set", [EdgeIndex(1,1)])
    addvertexset!(grid, "vert_set", [VertexIndex(1,1)])

    #Test function add
    addfaceset!(grid, "left_face", (x)-> x[1] ≈ 0.0)
    addedgeset!(grid, "left_lower_edge", (x)-> x[1] ≈ 0.0 && x[3] ≈ 0.0)
    addvertexset!(grid, "left_corner", (x)-> x[1] ≈ 0.0 && x[2] ≈ 0.0 && x[3] ≈ 0.0)

    @test 1 in Ferrite.getnodeset(grid, "node_set")
    @test FaceIndex(1,5) in getfaceset(grid, "left_face")
    @test EdgeIndex(1,4) in getedgeset(grid, "left_lower_edge")
    @test VertexIndex(1,1) in getvertexset(grid, "left_corner")

end

@testset "grid coloring" begin
    function test_coloring(grid)
        for alg in (Ferrite.GREEDY, Ferrite.WORKSTREAM)
            color_vectors = create_coloring(grid; alg=alg)
            @test sum(length, color_vectors) == getncells(grid)
            @test union(Set.(color_vectors)...) == Set(1:getncells(grid))
            conn = Ferrite.create_incidence_matrix(grid)
            for color in color_vectors, c1 in color, c2 in color
                @test !conn[c1, c2]
            end
        end
    end
    test_coloring(generate_grid(Line, (5,)))
    test_coloring(generate_grid(QuadraticLine, (5,)))
    test_coloring(generate_grid(Triangle, (5, 5)))
    test_coloring(generate_grid(QuadraticTriangle, (5, 5)))
    test_coloring(generate_grid(Quadrilateral, (5, 5)))
    test_coloring(generate_grid(QuadraticQuadrilateral, (5, 5)))
    test_coloring(generate_grid(Tetrahedron, (5, 5, 5)))
    # test_coloring(generate_grid(QuadraticTetrahedron, (5, 5, 5)))
    test_coloring(generate_grid(Hexahedron, (5, 5, 5)))
    # test_coloring(generate_grid(QuadraticHexahedron, (5, 5, 5)))
end


text:
using CSV, TimeSeries, Dates, Statistics, StatsBase, Measures, Plots; pyplot()

df = CSV.read("../data/oneOnEpsilonBlogs.csv",copycols = true)
tsA = TimeArray(Date.(df.Day,Dates.DateFormat("m/d/y")),df.Users)
tsB = moving(mean,tsA,7,padding = true)
tsC = TimeArray(timestamp(tsA), values(tsA) - values(tsB))

dow = dayofweek.(timestamp(tsA))
dayDiv = [filter((x)->!isnan(x), values(tsC)[dow .== d]) for d in 1:7]
dayMeans = mean.(dayDiv)
dailyCorrection = dayMeans[dow]
errs = filter((x)->!isnan(x),values(tsC) - dailyCorrection)
diffs = diff(errs)

lags = 1:50
acc = autocor(diffs,lags)

default(legend = false)
p1 = plot(diffs, c=:blue,
        xlabel="Day",ylabel="Difference Between Corrected Days")
p2 = plot(lags,acc,line=:stem, c=:blue,
        xlabel="Lag",ylabel="Autocorrelation")
plot(p1,p2,size=(900,400),margin = 5mm)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
module CLI # Bukdu

using ..Bukdu.Routing
using ..Bukdu.Naming

"""
    CLI.routes(io::IO = stdout)

Showing the routing table.
"""
function routes(io::IO = stdout)
    A = Routing.store[:routing_tables]
    isempty(A) && return
    ncols = length(first(A)) # (verb, url, C, action, pipe)
    nrows = length(A)
    paddings = map(1:ncols) do col
        maximum([(length ∘ string)(row[col]) for row in A]) .+ 2
    end
    function f(idx, el, lastcolumn)
        if idx == lastcolumn
            el
        else
            rpad(el, paddings[idx])
        end
    end
    for row in A
        lastcolumn = isempty(row[ncols]) ? ncols-1 : ncols
        print.(Ref(io), [f(idx, el, lastcolumn) for (idx, el) in enumerate(row[1:lastcolumn])])
        println(io)
    end
end

end # module Bukdu.CLI


text:
"""
Graphics context attached to a window. Used to register drawing commands on the window surface.
"""
mutable struct GraphicsContext
    conn::Connection
    id::xcb_window_t
    function GraphicsContext(conn::Connection, win_id::xcb_window_t)
        id = XCB.xcb_generate_id(conn)
        gc = new(conn, id)
        @check :error xcb_create_gc(gc.conn, gc.id, win_id, 0, C_NULL)
        Base.finalizer(x -> @check(:error, xcb_free_gc(gc.conn, gc.id)), gc)
    end
end

function GraphicsContext(conn, win_id; attributes=[], values=[])
    gc = GraphicsContext(convert(Connection, conn), convert(xcb_window_t, win_id))
    set_attributes(gc, attributes, values)
    gc
end

function set_attributes(gc::GraphicsContext, attributes, values)
    values = values[sortperm(collect(attributes))]
    list = zeros(UInt32, 32)
    setindex!.(Ref(list), values, 1:length(values))
    @flush @check xcb_change_gc(gc.conn, gc.id, reduce(|, attributes), list)
end


text:
function runOneCmd(cmd; igs::Bool=false)
    p = (igs == true) ? "ignore" : ""
    println("$(p)>> ", string(cmd))
    if igs
        cmd = ignorestatus(cmd)
    end
    run(cmd)
end

function pleaseRun(cmd; igs::Bool=false)
    if isa(cmd, Array)
        for i = 1:length(cmd)
            runOneCmd(cmd[i]; igs = igs)
        end
    else
        runOneCmd(cmd; igs = igs)
    end
end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
module PolyaUrnSimulator

using ProgressMeter
using StatsBase

export SSW!, WSW!, run_simulation

function poppush!(v::Vector{T}, e::T) where {T}
  pop!(v)
  pushfirst!(v, e)
end

function SSW!(_::Int, buffer::Vector{Int}, urn::Vector{Int}, exchanged::Int)
  if !(exchanged in buffer)
    poppush!(buffer, exchanged)
  end
end

function WSW!(buffer_size::Int, buffer::Vector{Int}, urn::Vector{Int}, exchanged::Int)
  set = Set{Int}()
  while length(set) != buffer_size
    push!(set, rand(urn))
  end
  buffer .= collect(set)
end

"""
アクションの呼び出され側ノードを選ぶ
アクションの呼び出し側ノードの壺内からランダムに1個選び出す
"""
function get_called(caller_node_urns::Vector{Int}, caller::Int)::Int
  while true
    called = rand(caller_node_urns)
    if called != caller
      return called
    end
  end
  throw(BoundsError())
end

"""
アクションの呼び出し側ノードを選ぶ
各々のノードが持つ壺の重さに比例する確率でノードを選ぶ
"""
function get_caller(node_urn_sizes::Vector{Int})::Int
  # 壺のサイズに基づいてインタラクションの起点ノードを選択
  return sample(Weights(node_urn_sizes))
end

function interact!(
  rho::Int,
  nu_plus_one::Int,
  strategy!::Function,
  node_urns::Vector{Vector{Int}},
  node_buffers::Vector{Vector{Int}},
  node_urn_sizes::Vector{Int},
  total_urn_size::Int,
  history::Vector{Tuple{Int,Int}};
  async_strategy::Bool = false
)
  ##### Model Rule (2) >>> #####

  "アクションを起こす起点のノード"
  caller = get_caller(node_urn_sizes)

  "アクションを起こされる終点のノード"
  called = get_called(node_urns[caller], caller)

  append!(history, [(caller, called)])

  ##### <<< Model Rule (2) #####

  ##### Model Rule (5) >>> #####

  # もしcalledノードが今まで呼ばれたことの無いノード(=壺が空のノード)である場合
  if node_urn_sizes[called] == 0

    # nu_plus_one個のノードを生成
    for _ = 1:nu_plus_one
      append!(node_urns, Vector{Int}[Int[]])
      append!(node_buffers, Vector{Int}[Int[]])
      append!(node_urn_sizes, Int[0])
    end

    # 生成したノードをcalledノードの壺とメモリバッファに追加
    for i = 1:nu_plus_one
      idx = length(node_urns) - (i - 1)
      append!(node_urns[called], Int[idx])
      append!(node_buffers[called], Int[idx])
      node_urn_sizes[called] += 1
      total_urn_size += 1
    end
  end

  ##### <<< Model Rule (5) #####

  ##### Model Rule (3) >>> #####
  append!(node_urns[caller], fill(called::Int, rho))
  node_urn_sizes[caller] += rho
  total_urn_size += rho

  append!(node_urns[called], fill(caller::Int, rho))
  node_urn_sizes[called] += rho
  total_urn_size += rho

  ##### <<< Model Rule (3) #####

  ##### Model Rule (4) >>> #####
  # メモリバッファを交換する
  append!(node_urns[caller], node_buffers[called])
  node_urn_sizes[caller] += nu_plus_one
  total_urn_size += nu_plus_one

  append!(node_urns[called], node_buffers[caller])
  node_urn_sizes[called] += nu_plus_one
  total_urn_size += nu_plus_one

  # メモリバッファを更新する
  strategy!(nu_plus_one, node_buffers[caller], node_urns[caller], called)
  if !async_strategy # async_strategyモードのときは、呼び出し側のノードのみがメモリバッファを更新する
    strategy!(nu_plus_one, node_buffers[called], node_urns[called], caller)
  end
  ##### <<< Model Rule (4) #####
end

"""壺モデルを最初から走らせる
## returns
`(history, node_urns, node_buffers)`
"""
function run_simulation(rho::Int, nu::Int, strategy!::Function, steps::Int; async_strategy::Bool = false, show_progress::Bool = true)

  nu_plus_one = nu + 1

  ##### Model Rule (1) >>> #####

  # 初期値の設定
  "壺"
  node_urns::Vector{Vector{Int}} = [[2], [1]]

  "メモリバッファ"
  node_buffers::Vector{Vector{Int}} = [Int[], Int[]]

  "壺のサイズ"
  node_urn_sizes = Int[1, 1]

  "系内に存在する全ての壺のサイズの合計" # 系内のノードの壺サイズを全て毎回計算すると計算時間がかかるので、保持しておく
  total_urn_size::Int = 2 # 最初は2つのノードがそれぞれお互いを持っているので 2

  "履歴"
  history = Tuple{Int,Int}[]

  # 初期ノードが初期状態でメモリバッファに持っているノードを作成
  for i = 1:2*nu_plus_one
    push!(node_urns, Int[])
    push!(node_urn_sizes, 0)
    push!(node_buffers, Int[])
  end

  # 初期ノードが初期状態でメモリバッファに持っているノードを設定
  for i = 1:nu_plus_one
    push!(node_urns[1], 2 + i)
    node_urn_sizes[1] += 1
    total_urn_size += 1
    push!(node_urns[2], 2 + i + nu_plus_one)
    node_urn_sizes[2] += 1
    total_urn_size += 1

    push!(node_buffers[1], 2 + i)
    push!(node_buffers[2], 2 + i + nu_plus_one)
  end

  ##### <<< Model Rule (1) #####

  p = Progress(length(1:steps); showspeed = true, enabled = show_progress)

  @inbounds @simd for _ = 1:steps
    interact!(rho, nu_plus_one, strategy!, node_urns, node_buffers, node_urn_sizes, total_urn_size, history; async_strategy = async_strategy)
    ProgressMeter.next!(p)
  end

  return (history, node_urns, node_buffers)
end

"""壺モデルを途中から走らせる
## returns
`(history, node_urns, node_buffers)`
"""
function run_simulation(
  rho::Int,
  nu_plus_one::Int,
  strategy!::Function,
  node_urns::Vector{Vector{Int}},
  node_buffers::Vector{Vector{Int}},
  steps::Int;
  async_strategy::Bool = false,
  show_progress::Bool = true
)
  node_urn_sizes::Vector{Int} = map(length, node_urns)
  total_urn_size::Int = sum(node_urn_sizes)

  history::Vector{Tuple{Int,Int}} = Tuple{Int,Int}[]

  p = Progress(steps; showspeed = true, enabled = show_progress)

  for _ = 1:steps
    interact!(rho, nu_plus_one, strategy!, node_urns, node_buffers, node_urn_sizes, total_urn_size, history; async_strategy = async_strategy)
    next!(p)
  end

  return (history, node_urns, node_buffers)
end

end

max_stars_repo_name: sudame/PolyaUrnSimulator

max_issues_repo_name: sudame/PolyaUrnSimulator

max_forks_repo_name: sudame/PolyaUrnSimulator


text:
#=
  MATLAB_HOME needs to be set to where Matlab lives for using MATLAB,
  for example you could put the following in .bash_profile
  export MATLAB_HOME="/Applications/MATLAB_R2017a.app"

  and need gtimeout, which can get from brew install coreutils
=#

using MATLAB

function timeLimitCmg(limit, la, b; tol::Real=1e-8, maxits=1000, verbose=false)

    mf = MatFile("fromJulia.mat","w")
    put_variable(mf, "la", la)
    put_variable(mf, "b", b)
    put_variable(mf, "tol", tol)
    put_variable(mf, "maxits", maxits)
    close(mf)

    fn = "$(Pkg.dir("Laplacians"))/matlab/timeCmg.m"
    mat = ENV["MATLAB_HOME"]
    matlab = "$(mat)/bin/matlab"
    cmd = `gtimeout $(limit) $(matlab) -nojvm \< $(fn)`

    t0 = now()

    bt = Inf
    st = Inf
    err = Inf
    iter = Inf
    start = -Inf

    try
        run(cmd)
    catch
        println("Matlab Died")
    end

    try
        mf = MatFile("toJulia.mat")

        bt = get_variable(mf,:bt)
        st = get_variable(mf,:st)
        err = get_variable(mf,:relres)
        iter = get_variable(mf,:iter)
        start = get_variable(mf,:startTime)
        close(mf)
    catch
        println("No .mat file generated")
    end

    if verbose
      println("Build Time: ", bt)
      println("Solve Time: ", st)
      println("Iterations: ", iter)
      println("error: ", err)
      d1 = DateTime(start,"d-u-yyyy H:M:S")
      println("Time to load and start Matlab: $(d1-DateTime(t0))/1000)")
    end

    return (st, bt, iter, err)

end

function timeLimitIcc(limit, la, b; tol::Real=1e-8, maxits=1000, verbose=false)

    mf = MatFile("fromJulia.mat","w")
    put_variable(mf, "la", la)
    put_variable(mf, "b", b)
    put_variable(mf, "tol", tol)
    put_variable(mf, "maxits", maxits)
    close(mf)

    fn = "$(Pkg.dir("Laplacians"))/matlab/timeIccLap.m"
    mat = ENV["MATLAB_HOME"]
    matlab = "$(mat)/bin/matlab"
    cmd = `gtimeout $(limit) $(matlab) -nojvm \< $(fn)`

    t0 = now()

    bt = Inf
    st = Inf
    err = Inf
    iter = Inf
    start = -Inf

    try
        run(cmd)
    catch
        println("Matlab Died")
    end

    try
        mf = MatFile("toJulia.mat")

        bt = get_variable(mf,:bt)
        st = get_variable(mf,:st)
        err = get_variable(mf,:relres)
        iter = get_variable(mf,:iter)
        start = get_variable(mf,:startTime)
        close(mf)
    catch
        println("No .mat file generated")
    end

    if verbose
      println("Build Time: ", bt)
      println("Solve Time: ", st)
      println("Iterations: ", iter)
      println("error: ", err)
      d1 = DateTime(start,"d-u-yyyy H:M:S")
      println("Time to load and start Matlab: $(d1-DateTime(t0))/1000)")
    end

    return (st, bt, iter, err)

end

function timeLimitLamg(limit, la, b; tol::Real=1e-8, maxits=1000, verbose=false)

    mf = MatFile("fromJulia.mat","w")
    put_variable(mf, "la", la)
    put_variable(mf, "b", b)
    put_variable(mf, "tol", tol)
    put_variable(mf, "maxits", maxits)
    close(mf)

    fn = "$(Pkg.dir("Laplacians"))/matlab/timeLamg.m"
    mat = ENV["MATLAB_HOME"]
    matlab = "$(mat)/bin/matlab"
    cmd = `gtimeout $(limit) $(matlab) -nojvm \< $(fn)`

    t0 = now()

    bt = Inf
    st = Inf
    err = Inf
    iter = Inf
    start = -Inf

    try
        run(cmd)
    catch
        println("Matlab Died")
    end

    try
        mf = MatFile("toJulia.mat")

        bt = get_variable(mf,:bt)
        st = get_variable(mf,:st)
        err = get_variable(mf,:relres)
        iter = get_variable(mf,:iter)
        start = get_variable(mf,:startTime)
        close(mf)
    catch
        println("No .mat file generated")
    end

    if verbose
      println("Build Time: ", bt)
      println("Solve Time: ", st)
      println("Iterations: ", iter)
      println("error: ", err)
      d1 = DateTime(start,"d-u-yyyy H:M:S")
      println("Time to load and start Matlab: $(d1-DateTime(t0))/1000)")
    end

    return (st, bt, iter, err)

end

function timeLimitLamgSddm(limit, la, b; tol::Real=1e-8, maxits=1000, verbose=false)

    mf = MatFile("fromJulia.mat","w")
    put_variable(mf, "la", la)
    put_variable(mf, "b", b)
    put_variable(mf, "tol", tol)
    put_variable(mf, "maxits", maxits)
    close(mf)

    fn = "$(Pkg.dir("Laplacians"))/matlab/timeLamgSddm.m"
    mat = ENV["MATLAB_HOME"]
    matlab = "$(mat)/bin/matlab"
    cmd = `gtimeout $(limit) $(matlab) -nojvm \< $(fn)`

    t0 = now()

    bt = Inf
    st = Inf
    err = Inf
    iter = Inf
    start = -Inf

    try
        run(cmd)
    catch
        println("Matlab Died")
    end

    try
        mf = MatFile("toJulia.mat")

        bt = get_variable(mf,:bt)
        st = get_variable(mf,:st)
        err = get_variable(mf,:relres)
        iter = get_variable(mf,:iter)
        start = get_variable(mf,:startTime)
        close(mf)
    catch
        println("No .mat file generated")
    end

    if verbose
      println("Build Time: ", bt)
      println("Solve Time: ", st)
      println("Iterations: ", iter)
      println("error: ", err)
      d1 = DateTime(start,"d-u-yyyy H:M:S")
      println("Time to load and start Matlab: $(d1-DateTime(t0))/1000)")
    end

    return (st, bt, iter, err)

end


text:
cd("git/SparseInverseProblems.jl")
# Function to test running 2d simple ADCG

# import Pkg;

# this uses a forked version
# git clone https://github.com/marius311/SparseInverseProblems.jl.git
Pkg.clone("git@github.com:albertasdvirnas/SparseInverseProblems.jl.git")
# which includes a toml file and some
# modifications for this to run in Julia properly

# load all necessary functions.
using Pkg
Pkg.activate(".")
using SparseInverseProblems
using Compose
using SpecialFunctions
using Images
using Plots
import SparseInverseProblems: lmo, phi, solveFiniteDimProblem, localDescent,
getStartingPoint, parameterBounds, computeGradient
using SparseInverseProblems.Util
using SparseInverseProblems
using NLopt, SparseInverseProblems.Util
using LinearAlgebra
include("examples/smi/gaussblur.jl")

# functions which are changed from default..

function computeGradient(model :: GaussBlur2D, weights :: Vector{Float64},
  thetas :: Matrix{Float64}, v :: Vector{Float64})
  gradient=0
  v = reshape(v, model.n_pixels, model.n_pixels)
  #
  gradient = zeros(length(thetas))
  # #allocate temporary variables...
  f_x = zeros(model.n_pixels)
  f_y = zeros(model.n_pixels)
  fpy = zeros(model.n_pixels)
  fpx = zeros(model.n_pixels)
  v_x = zeros(model.n_pixels)
  v_y = zeros(model.n_pixels)
  v_yp = zeros(model.n_pixels)

  # #compute gradient
  for i = 1:size(thetas,2)
    point = vec(thetas[:,i])
    computeFG(model, point[1], f_x, fpx)
    computeFG(model, point[2], f_y, fpy)
    v_x = mul!(v_x,v,f_x)
    v_y = mul!(v_y,v,f_y)
    v_yp = mul!(v_yp,v,fpy)
    function_value = dot(f_y, v_x)
    g_x = dot(v_y, fpx)
    g_y = dot(v_x, fpy)

    gradient[:,i] = weights[i]*[g_x; g_y]
  end
  return gradient
end

function gensim()

  noise_mean = 0.0021172377176794793
  sigmasq = let
    lambda = 723.0
    NA = 1.4
    FWHM = lambda/(2*NA)
    sigma = FWHM/(2*log(2.0))
    (sigma/(64*100.0))^2
  end

  gb_sim = GaussBlur2D(sigmasq,64,500)

  return gb_sim
end

function phi(s :: GaussBlur2D, parameters :: Matrix{Float64},weights :: Matrix{Float64})
  n_pixels = s.n_pixels
  if size(parameters,2) == 0
    return zeros(n_pixels*n_pixels)
  end
  v_x = computeFs(vec(parameters[1,:]),n_pixels,s.sigma).*weights
  v_y = computeFs(vec(parameters[2,:]),n_pixels,s.sigma)
  return vec(v_y*v_x')
end

# Generate simulation
gb_sim = gensim()

parameters = rand(2,4);
weights = 100000*rand(1,4);

# this gives us the output
y = phi(gb_sim,parameters, weights)

# for showing heatmap, we need to reshape to 64x64 (or change how these are chosen in gb_sim)
heatmap(reshape(y,64,64))

# callback function. remove mean  noise from target
target = vec(y).-noise_mean;
function callback(old_thetas,thetas, weights,output,old_obj_val)
  #evalute current OV
  new_obj_val,t = loss(LSLoss(), output - target)
  if old_obj_val - new_obj_val < 7E-5
    return true
  end
  return false
end

# now default parameters for ADCG calculation
lossFn =LSLoss()
tau = 100.0;
max_iters = 50
min_optimality_gap = 1E-5
max_cd_iters = 200
fully_corrective = false
@assert(tau > 0.0)

output = zeros(length(y))

residual = output - y
#evalute the objective value and gradient of the loss
objective_value, grad = loss(lossFn, residual)
#compute the next parameter value to add to the support

# end
# next parameter value and score. Using the default gives us
# score = 0.0..
theta,score = lmo(gb_sim, grad)

# go through lmo methods step by step.
# 1 starting point
initial_x = getStartingPoint(gb_sim, grad)
#
# v = reshape(grad, gb_sim.n_pixels, gb_sim.n_pixels)
#
# # grad2 = reshape(grad, gb_sim.n_pixels, gb_sim.n_pixels)
# ng = length(gb_sim.grid)
# grid_objective_values = vec(gb_sim.grid_f'*v*gb_sim.grid_f)
# best_point_idx = argmin(grid_objective_values)

# parameters
model = gb_sim
v = residual

lb,ub = parameterBounds(model)
initial_x = getStartingPoint(model, v)
println(initial_x)
p = length(lb)

function fg!(point :: Vector{Float64},gradient_storage :: Vector{Float64})
  output = phi(model,reshape(point,p,1), ones(1,1))
  ip = dot(output,v)
  s = sign(ip)
  gradient_storage[:] = -s*computeGradient(model, [1.0],reshape(point,length(point),1), v)
  # gradient_storage[:] = -s*computeGradient(model, [1.0],reshape(point,length(point),1), v)
  return -s*ip
end

 phi(model,reshape(initial_x,2,1),ones(1,1))

 function initializeOptimizer!(model :: BoxConstrainedDifferentiableModel, opt :: Opt)
  ftol_abs!(opt, 1e-6)
  xtol_rel!(opt, 0.0)
  maxeval!(opt, 200)
end

#
# i=1
# dpsi(model, vec(initial_x[:,i]))'*v
# computeGradient(model,ones(1,1))

opt = Opt(:LD_MMA, p)
initializeOptimizer!(model, opt)
min_objective!(opt, fg!)
# opt.min_objective = f_and_g!
lower_bounds!(opt, lb)
upper_bounds!(opt, ub)
(optf,optx,ret) = optimize(opt, initial_x)

# Niw optf is non-zero.
println(optf)
println(optx)

#
#
# weights = vec(ones(1,1))
# thetas = reshape(initial_x,2,1)
# v = residual
# grfff = computeGradient(model, weights,thetas, v)


text:
function coefficients{C<:Chebyshev,DD<:IntervalDomain}(f::Vector,sp1::JacobiWeight{C,DD},sp2::JacobiWeight{C,DD})
    α,β=sp1.α,sp1.β
    c,d=sp2.α,sp2.β
    tol=10E-5  #TODO:tol choice
    @assert isapproxinteger(α-c) && isapproxinteger(β-d)
    if isapprox(c,α) && isapprox(d,β)
        f
    elseif c>α && d>β
        g=Fun(f,sp1.space)
        @assert abs(first(g))<tol && abs(last(g))<tol
        coefficients(divide_singularity(g).coefficients,JacobiWeight(α+1,β+1,sp1.space),sp2)
    elseif c>α
        g=Fun(f,sp1.space)
        @assert abs(first(g))<tol
        coefficients(divide_singularity(false,g).coefficients,JacobiWeight(α+1,β,sp1.space),sp2)
    elseif d>β
        g=Fun(f,sp1.space)
        @assert abs(last(g))<tol
        coefficients(divide_singularity(true,g).coefficients,JacobiWeight(α,β+1,sp1.space),sp2)
    elseif c<α && d<β
        x=Fun()
        coefficients(((1-x^2)*Fun(f,Chebyshev())).coefficients,JacobiWeight(α-1,β-1,sp1.space),sp2)
    elseif d<β
        x=Fun()
        coefficients(((1-x)*Fun(f,Chebyshev())).coefficients,JacobiWeight(α,β-1,sp1.space),sp2)
    elseif c<α
        x=Fun()
        coefficients(((1+x)*Fun(f,Chebyshev())).coefficients,JacobiWeight(α-1,β,sp1.space),sp2)
    else
        error("Something has gone wrong")
    end
end


text:
# Helper functions
function z_subproblem!{T}(z_fft::Array{Complex{T}, 3},
                          t_fft::Array{Complex{T}, 3}, mu_t::T, lambda_t_fft::Array{Complex{T}, 3},
                          DtD_fft::Array{Complex{T}, 3}, Dtx_fft::Array{Complex{T}, 3})
  convLen = size(t_fft, 3);
  dictSize = size(DtD_fft, 1);
  numImgs = size(Dtx_fft, 2);
  invmu_t = 1/mu_t;
  # (Eq. 7)
  invA = zeros(Complex128, dictSize, dictSize); b = zeros(Complex128, dictSize, numImgs);
  for i=1:convLen;
    # z_fft[:, :, i] = ( DtD_fft[:, :, i] + mu_t.*I ) \
    #   ( Dtx_fft[:, :, i] + mu_t.*t_fft[:, :, i] - lambda_t_fft[:, :, i] )
    # using the Sherman-Morrison formula
    a = 1 + invmu_t*sum(diag(DtD_fft[:, :, i]));
    scalar = -1/(a*mu_t*mu_t);
    broadcast!(*, invA, scalar, DtD_fft[:, :, i]);
    for j=1:dictSize; invA[j, j] = invmu_t + invA[j, j]; end
    copy!(b, Dtx_fft[:, :, i]);
    b += mu_t.*t_fft[:, :, i];
    b -= lambda_t_fft[:, :, i];
    z_fft[:, :, i] = invA*b;
  end
  return z_fft::Array{Complex{T}, 3};
end

function t_subproblem!{T}(t::Array{T, 2}, z::Array{T, 2}, mu_t::T, lambda_t::Array{T, 2}, beta::T)
  # (Eq. 9)
  if mu_t>0;
    # shrink!(t, z+lambda_t/mu_t, beta/mu_t);
    for i=1:length(t);
      t[i] = lambda_t[i]/mu_t + z[i];
    end
    shrink!(t, t, beta/mu_t);
  else;
    fill!(t, 0);
  end
  return t::Array{T, 2};
end

function fcsc_coding_fastmulti{T}(data::FastCodingData{T}, optParams::FastCodingOpts{T})
# Alg. parameters:
#   mu_max: 1e5
#   tau: [1.01 1.5] larger means faster convergence but sacrifices primal feasibility,
#        smaller means slower convergence

  const mu_max = 1e5;
  const tau = 1.01;

  const Dtx_fft      = data.Dtx_fft;
  const DtD_fft      = data.DtD_fft;
  const convH        = data.convH;
  const convW        = data.convW;
        z            = optParams.z;
        z_fft        = optParams.z_fft;
        t            = optParams.t;
        t_fft        = optParams.t_fft;
        lambda_t     = optParams.lambda_t;
        lambda_t_fft = optParams.lambda_t_fft;
  const beta::T      = optParams.beta;
        mu_t::T      = optParams.mu_t;
  const max_iter     = optParams.max_iter;
  const tol          = optParams.tol;
  const DEBUG        = optParams.DEBUG;
        debugData    = optParams.debugData;

  numImgs  = size(Dtx_fft, 2);
  dictSize = size(DtD_fft, 1);
  convLen  = convH*convW;

  assert(size(DtD_fft, 3)==convLen, "DtD_fft does not have the correct depth.")
  assert(size(Dtx_fft, 1)==dictSize && size(Dtx_fft, 3)==convLen,
    "Dtx_fft does not have the correct dimensions.")
  z = initcheck(z, T, (convLen*dictSize, numImgs), "z does not have the correct dimensions.");
  z_fft = initcheck(z_fft, Complex{T}, (dictSize, numImgs, convLen), "z_fft does not have the correct dimensions.");
  update_t = false;
  if vecnorm(z)==0;
    # Initialize z to pure reconstruction solution
    fill!(z_fft, 0);
    # Increase from 1e-14 if singular errors still occur.
    I = 1e-14.*eye(dictSize);
    for i=1:convLen;
      z_fft[:, :, i] = ( DtD_fft[:, :, i]+I ) \ Dtx_fft[:, :, i];
    end
    for j=1:numImgs;
      z_tmp = reshape(squeeze(z_fft[:, j, :], 2).', convH, convW, dictSize);
      z[:, j] = reshape(real(ifft2(z_tmp)), dictSize*convLen, 1);
    end
    update_t = true;
  end

  t = initcheck(t, T, (convLen*dictSize, numImgs), "t does not have the correct dimensions.");
  lambda_t = initcheck(lambda_t, T, (convLen*dictSize, numImgs), "lambda_t does not have the correct dimensions.");
  lambda_t_fft = initcheck(lambda_t_fft, Complex{T}, (dictSize, numImgs, convLen), "lambda_t_fft does not have the correct dimensions.");
  if isempty(t_fft);
    t_fft = zeros(Complex{T}, dictSize, numImgs, convLen);
    update_t = true;
  else;
    assert(size(t_fft)==(dictSize, numImgs, convLen), "t_fft does not have the correct dimensions.")
  end
  if update_t;
    copy!(t, z);
    copy!(t_fft, z_fft);
  end

  if DEBUG;
    debugData["z"] = z;
    debugData["t"] = t;
    debugData["lambda_t"] = lambda_t;
    (prev_obj, rerr, spar) = compute_lagrangian(debugData);
  end

  iter = 0;
  converged = false;
  while ~converged;
    iter += 1;

    if max_iter>1; @printf("\b\b\b%03d", iter); end

    # Solve for z* in the frequency domain
    z_subproblem!(z_fft, t_fft, mu_t, lambda_t_fft, DtD_fft, Dtx_fft);
    # TODO parfor this loop
    for j=1:numImgs;
      z_tmp = reshape(squeeze(z_fft[:, j, :], 2).', convH, convW, dictSize);
      z[:, j] = reshape(real(ifft2(z_tmp)), dictSize*convLen, 1);
    end
    if DEBUG;
      debugData["z"] = z;
      (obj, rerr, spar) = compute_lagrangian(debugData);
      is_mono = (prev_obj>=obj) ? 'O' : 'X';
      @printf("          Z: %c obj=%.10f (rerr=%.10f, spar=%.10f)\n", is_mono, obj, rerr, spar);
      prev_obj = obj;
    end

    # Solve for t* in the spatial domain
    t_subproblem!(t, z, mu_t, lambda_t, beta);
    if DEBUG;
      debugData["t"] = t;
      (obj, rerr, spar) = compute_lagrangian(debugData);
      is_mono = (prev_obj>=obj) ? 'O' : 'X';
      @printf("          T: %c obj=%.10f (rerr=%.10f, spar=%.10f)\n", is_mono, obj, rerr, spar);
      prev_obj = obj;
    end

    if max_iter==1; return (z, t); end # Don't do unnecessary work
    # TODO parfor this loop
    for j=1:numImgs;
      t_fft[:, j, :] = reshape(fft2(reshape(t[:, j], convH, convW, dictSize)), convLen, dictSize).';

      # Update the Lagrangian multipliers
      for i=1:convLen*dictSize; lambda_t[i, j] += mu_t*( z[i, j]-t[i, j] ); end
      lambda_t_fft[:, j, :] = reshape(fft2(reshape(lambda_t[:, j], convH, convW, dictSize)), convLen, dictSize).';
    end
    if DEBUG;
      debugData["lambda_t"] = lambda_t;
      (obj, rerr, spar) = compute_lagrangian(debugData);
      is_mono = (prev_obj>=obj) ? 'O' : 'X';
      @printf("          L: %c obj=%.10f (rerr=%.10f, spar=%.10f)\n", is_mono, obj, rerr, spar);
      prev_obj = obj;
    end

    # Update penalty term
    mu_t = min(tau*mu_t, mu_max);
    if DEBUG;
      debugData["mu_t"] = mu_t;
      (obj, rerr, spar) = compute_lagrangian(debugData);
      is_mono = (prev_obj>=obj) ? 'O' : 'X';
      @printf("          M: %c obj=%.10f (rerr=%.10f, spar=%.10f)\n", is_mono, obj, rerr, spar);
      prev_obj = obj;
    end

    codesConverged = true;
    for j=1:numImgs;
      if mu_t>0 && mean((z[:, j]-t[:, j]).^2)>tol;
        codesConverged = false;
        break;
      end
    end
    if iter>=max_iter || codesConverged;
      converged = true;
    end
  end

  return (z, t);
end

max_stars_repo_name: bkong/FCSC


text:
"""
Run the PIMDP synthesis
"""
function run_pimdp_synthesis(pimdp::PIMDP, pimdp_filename::String; add_opt=false, horizon=-1)
    write_pimdp_to_file(pimdp, pimdp_filename)

    @info "Synthesizing a controller... (maximize pessimistic)"
    res_mat = run_unbounded_imdp_synthesis(pimdp_filename)

    if add_opt
        @info "Synthesizing a controller... (maximize optimistic)"
        basename = pimdp_filename[1:end-3]
        opt_pimdp_filename = "$basename-optimistic.txt"
        res_mat_opt = run_imdp_synthesis(pimdp_filename, horizon, mode2="optimistic", save_mats=false)
        for j in 1:length(res_mat[:,1])
            if res_mat[j, 3] == res_mat_opt[j,3] && res_mat_opt[j,4] > res_mat[j,4]
                res_mat[j, 2] = res_mat_opt[j,2]
                res_mat[j, 4] = res_mat_opt[j,4]
            end
        end
    end

    return res_mat
end


text:
# An Initializer is used to initialize network parameters
export Initializer
export ConstantInitializer
export XavierInitializer
export GaussianInitializer
export OrthogonalInitializer

@compat abstract type Initializer end # The root type of all initializer

immutable NullInitializer <: Initializer end
function init(::NullInitializer, blob::Blob)
  # do nothing
end

immutable ConstantInitializer <: Initializer
  value
end

function init(initializer::ConstantInitializer, blob::Blob)
  fill!(blob, initializer.value)
end

################################################################################
# An initializer  based on the paper [Bengio and Glorot 2010]: Understanding
# the difficulty of training deep feedforward neuralnetworks, but does not
# use the fan_out value.
#
# It fills the incoming matrix by randomly sampling uniform data from
# [-scale, scale] where scale = sqrt(3 / fan_in) where fan_in is the number
# of input nodes.
#
# For a ND-tensor blob parameter, the product of the 1 ~ (N-1) dimensions
# are considered as fan-in, and the last dimension is considered as fan-out.
################################################################################
immutable XavierInitializer <: Initializer
end
function init(initializer::XavierInitializer, blob::Blob)
  fan_in = get_fea_size(blob)
  scale = convert(eltype(blob), sqrt(3.0 / fan_in))
  init_val = rand(eltype(blob), size(blob)) * 2scale - scale
  copy!(blob, init_val)
end

immutable GaussianInitializer <: Initializer
  mean :: AbstractFloat
  std  :: AbstractFloat
end
GaussianInitializer(;mean=0.0, std=1.0) = GaussianInitializer(mean, std)

function init(initializer::GaussianInitializer, blob::Blob)
  init_val = randn(size(blob)) * initializer.std + initializer.mean
  init_val = convert(Array{eltype(blob)}, init_val)
  copy!(blob, init_val)
end

####### Orthogonal Initializer ##############################################
#
# Based on https://github.com/Lasagne/Lasagne/blob/master/lasagne/init.py
#
#############################################################################

immutable OrthogonalInitializer <: Initializer
  gain::AbstractFloat
end
OrthogonalInitializer() = OrthogonalInitializer(1.0) # but use OrthogonalInitializer(sqrt(2)) for ReLU units

function init(initializer::OrthogonalInitializer, blob::Blob)
  dims = size(blob)
  if length(dims) < 2
    error("Must have length(size(blob)) >= 2")
  end
  flatshape = (dims[1], prod(dims[2:end]))
  x = randn(flatshape)
  u, _, v = svd(x)
  x = (size(u) == flatshape) ? u : v
  copy!(blob, convert(Array{eltype(blob)}, x*initializer.gain))
end


text:
using LinearFitXYerrors
using Test

@testset "LinearFitXYerrors.jl" begin

    # 0a - TEST DATA WITH NO ERRORS
    # INPUT DATA:  http://pba.ucdavis.edu/files/45007.pdf

    X = [16.4, 17.2, 17.6, 18.0, 18.2, 18.5]
    Y = [2.67, 2.75, 2.99, 3.14, 3.88, 4.23]

    stxy = linearfitxy(X,Y)

    @test   stxy.a    ≈ -9.4628187535736
    @test   stxy.b    ≈  0.7217838765009
    @test   stxy.σa   ≈  3.6596191146816
    @test   stxy.σb   ≈  0.2071823670235
    @test   stxy.σa95 ≈ 10.1607315778461
    @test   stxy.σb95 ≈  0.5752304688056
    @test   stxy.S    ≈  0.3537301433210
    @test   stxy.ρ    ≈  0.8672492974278

    # 0b - TEST DATA WITH NO ERRORS
    # INPUT DATA:  Table I from Altman and Gardner [1988]

    X = [91., 104, 107, 107, 106, 100, 92, 92, 105, 108]
    Y = [9.8, 7.4, 7.9, 8.3, 8.3, 9.0, 9.7, 8.8, 7.6, 6.9]

    stxy = linearfitxy(X,Y)

    @test   stxy.a    ≈ 20.1891143911439
    @test   stxy.b    ≈ -0.1167896678967
    @test   stxy.σa   ≈  2.6713186714425
    @test   stxy.σb   ≈  0.0263407280356
    @test   stxy.σa95 ≈  6.1600719027946
    @test   stxy.σb95 ≈  0.0607418277744
    @test   stxy.S    ≈  0.5484946692667
    @test   stxy.ρ    ≈ -0.8430654646904

    # 1 - TEST UNCORRELATED ERRORS
    # INPUT DATA:
    # Tables I and II from York [1966], data taken from Pearson (1901)
    # also in Cantrell (2008)

    X = [0.0, 0.9, 1.8, 2.6, 3.3, 4.4, 5.2, 6.1, 6.5, 7.4]
    Y = [5.9, 5.4, 4.4, 4.6, 3.5, 3.7, 2.8, 2.8, 2.4, 1.5]

    # the weights ωX, ωY are the inverse of the variances
    σX = 1 ./ sqrt.([1000., 1000, 500, 800, 200, 80,  60, 20, 1.8, 1])
    σY = 1 ./ sqrt.([1., 1.8, 4, 8, 20, 20, 70, 70, 100, 500])
    r = 0*X

    stxy = linearfitxy(X,Y; σX=σX, σY=σY, r=r)

    @test  stxy.a    ≈  5.47991022403286
    @test  stxy.b    ≈ -0.48053340744620
    @test  stxy.σa   ≈  0.35924652255111
    @test  stxy.σb   ≈  0.07062026952877
    @test  stxy.σa95 ≈  0.82842396656058
    @test  stxy.σb95 ≈  0.16285063356258
    @test  stxy.S    ≈  1.21790564053940
    @test  stxy.ρ    ≈ -0.94325300584483

    # 2 - TEST CORRELATED ERRORS

    # INPUT DATA: Table 3.1 from Amen, S. [2012]

    # Num,    ue,    n,   var_ue,   var_n,   cov_ue_n
    M = [1   23.371   5.992   0.106   0.323   0.170;
        2    24.462   8.899   1.692   10.195  3.887;
        3    22.414   5.361   0.124   0.460   0.207;
        4    27.628   10.271  1.088   1.460   1.255;
        5    23.802   7.114   0.107   0.496   0.191;
        6    23.092   7.975   0.281   1.651   0.582;
        7    23.263   7.486   0.316   1.695   0.661;
        8    20.843   4.169   0.389   0.335   0.357;
        9    24.421   5.645   0.381   0.716   0.507;
        10   23.716   9.220   0.733   3.656   1.401;
        11   22.032   5.358   0.094   0.268   0.106;
        12   21.448   3.166   0.167   0.170   0.001;
        13   21.631   4.002   0.127   0.106   0.082;
        14   19.947   2.073   0.015   0.015   0.001;
        15   20.652   3.338   0.040   0.104   0.021;
        16   20.588   2.028   0.007   0.011   0.002;
        17   20.745   1.978   0.010   0.010  -0.001;
        18   21.658   2.526   0.009   0.058  -0.003;
        19   19.539   2.043   0.044   0.070  -0.039;
        20   21.638   3.361   1.909   3.961  -2.557;
        21   19.863   2.452   0.014   0.019  -0.010;
        22   18.398   2.199   0.042   0.051  -0.037;
        23   20.729   1.725   0.013   0.026   0.003;
        24   20.511   1.906   0.003   0.010   0.001;
        25   22.140   3.372   0.047   0.108   0.062;
        26   20.510   2.130   0.004   0.020   0.002;
        27   20.141   1.897   0.006   0.011  -0.005]

    μₑ = M[:,2]
    η = M[:,3]
    σμₑ = sqrt.(M[:,4])   # standard deviation of errors in μₑ
    ση = sqrt.(M[:,5])     # standard deviation of errors in η
    cov_μₑ_η = M[:,6]   # covariance between errors in μₑ and η
    r = cov_μₑ_η ./ (σμₑ .* ση )   # correlation coefficient between errors

    # COMPUTE:
    stxy = linearfitxy(μₑ, η; σX=σμₑ, σY=ση, r=r, isplot=false)

    @test   stxy.a  ≈ -21.9553309045508
    @test   stxy.b  ≈   1.1762288806682
    @test   stxy.σa ≈   1.2846101711080
    @test   stxy.σb ≈   0.0578260394408
    @test   stxy.S  ≈   3.4885539308643
    @test   stxy.ρ  ≈   0.8435087435271

end


text:
# SerialMemCounter
# ================

"""
    SerialMemCounter{M<:AbstractMer,C<:CountMode}

KmerAnalysis' simplest kmer counter.

Build a sorted list (vector) of kmer counts (MerCount), serially and entirely in memory.

Build a SerialMemCounter using the [`serial_mem`](@ref) method. This returns
a SerialMemCounter.

A SerialMemCounter can be treated as a function or functor and passed to other
functions as an argument, and it can be called on other arguments.

It is a general purpose kmer counting method:

It basically collects all the kmers from iterating over the input, before sorting
the collected kmers and collapsing them into counts. To reduce allocations and
GC from many repeated calls of a SerialMemCounter, they keep a few buffers that
are used often for kmer collection and count collapsing.

Whilst general purpose for a variety of argument types, it is most suitable for:

1. Counting the kmers in a small number of long-ish `BioSequences` (e.g. a
   reference genome or low(ish)-complexity genome graph.).

2. Counting the kmers in a ReadDatastore if the dataset size is small enough
   such that all the kmers collected from every dataset can fit in a single
   machines memory. If you have a decent MacBook Pro, for example, datasets like
   _E. coli_ paired end sequencing reads will be no problem.

3. Other larger ReadDatastores IF you have the RAM to throw at it, e.g. you're on
   a HPC machine configured for large memory jobs.
"""
struct SerialMemCounter{M<:AbstractMer,C<:CountMode} <: AbstractKmerCounter{M}
    mer_buffer::Vector{M}
    count_buffer::Vector{MerCount{M}}
end

serial_mem(::Type{M}, count_mode::C) where {M<:AbstractMer,C<:CountMode} = SerialMemCounter{M,C}(Vector{M}(), Vector{MerCount{M}}())

@inline function count!(smc::SerialMemCounter{M,C}, input, args...) where {M<:AbstractMer,C<:CountMode}
    _serial_mem!(smc.count_buffer, smc.mer_buffer, C(), input, args...)
end

@inline function (smc::SerialMemCounter{M,C})(input, args...) where {M<:AbstractMer,C}
    count!(smc, input, args...)
    return copy(smc.count_buffer)
end

@inline function _serial_mem!(count_buffer::Vector{MerCount{M}}, mer_buffer::Vector{M}, count_mode::CountMode, input, args...) where {M<:AbstractMer}
    collect_mers!(mer_buffer, count_mode, input, args...)
    collapse_into_counts!(count_buffer, mer_buffer)
end

@inline function serial_mem(::Type{M}, input::PairedReads{A}, count_mode::CountMode) where {M,A}
    counts = Vector{MerCount{M}}()
    mers = Vector{M}()
    _serial_mem!(counts, mers, count_mode, input)
    return counts
end


text:
abstract AbstractTimestamp

immutable Timestamp{T} <: AbstractTimestamp
    timestamp::Union(Date,Day,DateTime)
    value::T
end

typealias TimestampArray{T} Array{Timestamp{T},1}

function Base.show{T}(io::IO, ts::Timestamp{T})
    print(io, ts.timestamp, " | ", ts.value)
end

# show{T}(io::IO, ts::Timestamp{T}) =
#     try
#        isnull(ts.value) ? print(io, ts.timestamp, " | ", "NA") : nothing
#     catch
#        print(io, ts.timestamp, " | ", ts.value)
# end

# single date

#@function getindex{T<:Timestamp}(tsa::TimestampArray{T}, d::TimeType)
function getindex(tsa::TimestampArray, d::TimeType)
    counter = Int[]
    for i in 1:length(tsa)
        if d == tsa[i].timestamp
            push!(counter, i)
        end
    end
    return tsa[counter]
end

# range of dates
function getindex(tsa::TimestampArray, dates::Union(Vector{Date}, Vector{DateTime}))
    counter   = Int[]
    tsa_dates = TimeType[t.timestamp for t in tsa]
  #  counter = int(zeros(length(dates)))
    for i in 1:length(dates)
        if findfirst(tsa_dates, dates[i]) != 0
        #counter[i] = findfirst(ta.timestamp, dates[i])
            push!(counter, findfirst(tsa_dates, dates[i]))
        end
    end
    tsa[counter]
end

#function getindex{T<:Timestamp}(tsa::TimestampArray{T}, r::Union(StepRange{Date}, StepRange{DateTime}))
function getindex(tsa::TimestampArray, r::Union(StepRange{Date}, StepRange{DateTime}))
    tsa[[r;]]
end


text:
io = IOBuffer(read(`top -b -i -H -w 256 -d 1 -n 10`));  # takes 10 seconds

using TopParser

processes = collect(TopParser.processes(io));

using DataFrames

df = DataFrame(processes);  # `processes` is a table

seekstart(io);

samples = map(s -> collect(s.processes), TopParser.samples(io));


text:
#-----------------------------------------------------------------------------# CTable
"""
A Columnar data table.  It is a simple wrapper around a NamedTuple of AbstractVectors (See
`Tables.ColumnTable`).  By default, ranges (e.g. `1:10`) will be `collect`-ed.

## Examples

    x = 1:10
    y = rand(10)

    # With field names:
    CTable((field1 = x, field2 = y))
    CTable(field1= x , field2 = y)

    # Auto-named:
    CTable((x, y))
    CTable(x, y)
    CTable([x y])
"""
struct CTable{T<:ColumnTable}
    cols::T
    function CTable(t::ColumnTable; allow_ranges=false)
        check(t)
        if allow_ranges
            t2 = t
        else
            t2 = map(t) do x
                x isa AbstractRange ? collect(x) : x
            end
        end
        new{typeof(t2)}(t2)
    end
end
function CTable(t::Tuple{Vararg{<:AbstractVector}}; kw...)
    CTable(namedtuple(tuple((Symbol("x$i") for i in 1:length(t))...), t); kw...)
end
CTable(; kw...) = CTable(namedtuple(kw))
CTable(args...) = CTable(args)
CTable(x::AbstractMatrix) = CTable(eachcol(x)...)

cols(t::CTable) = getfield(t, :cols)

schema(t::CTable) = CTable((
    field = collect(names(t)),
    type = collect(map(typeof, cols(t))),
    eltype = collect(map(eltype, cols(t)))
))

check(t::ColumnTable) = all(x -> length(x) == length(first(t)), t) || error("Columns have different lengths: $(map(length, t))")
nrows(t::CTable) = length(cols(t)[1])
ncols(t::CTable) = length(cols(t))
head(t::CTable, n=10) = t[1:min(length(t), n)]
tail(t::CTable, n=10) = t[max(1, end-n):end]

function rename(t::CTable, old_new::Pair...)
    for old_new in old_new
        old_key = field_name(t, old_new[1])
        new_key = field_name(t, old_new[2])
        names = collect(names(t))
        i = findfirst(==(old_key), names)
        names[i] = new_key
        t = CTable(namedtuple(names, cols(t)))
    end
    return t
end

#-----------------------------------------------------------------------------# Base
Base.summary(t::CTable) = "CTable ($(nrows(t)) × $(ncols(t)))"
Base.length(t::CTable) = isempty(t) ? 0 : length(cols(t)[1])

Base.names(t::CTable) = keys(cols(t))
Base.size(t::CTable) = (nrows(t), ncols(t))
Base.size(t::CTable, i::Integer) = size(t)[i]
Base.copy(t::CTable) = CTable(map(copy, cols(t)))
Base.keys(t::CTable) = Base.LinearIndices(Base.OneTo(length(t)))

Base.deleteat!(t::CTable, i) = (map(x -> deleteat!(x, i), cols(t)); t)
Base.isempty(t::CTable) = all(isempty, cols(t))
Base.hcat(a::CTable, b::CTable) = CTable(merge(cols(a), cols(b)))

function Base.vcat(a::CTable, b::CTable; cols=false)
    all(in.(names(a), Ref(names(b)))) || error("Cannot vcat")
    CTable(map((k,v) -> vcat(v, cols(b)[k]), pairs(cols(a))))
end

Base.push!(t::CTable, row) = (foreach(x -> push!(x...), zip(cols(t), row)); t)

function Base.append!(t::CTable{T}, t2::CTable{T}) where {T}
    map(x -> append!(x...), zip(cols(t), cols(t2)))
    t
end

Base.getindex(t::CTable, i::Integer) = map(x -> x[i], cols(t))
Base.getindex(t::CTable, i::AbstractVector{<:Integer}) = CTable(map(x -> x[i], cols(t)))
Base.getindex(t::CTable, x::Symbol) = cols(t)[x]
Base.getindex(t::CTable, x::AbstractVector{Symbol}) = CTable(select(cols(t), x))
function Base.setindex!(t::CTable, val, i::Integer)
    for (k,v) in pairs(cols(t))
        v[i] = val[k]
    end
    val
end
Base.lastindex(t::CTable) = nrows(t)

Base.getproperty(t::CTable, x::Symbol) = cols(t)[x]
function Base.setproperty!(t::CTable, val, x::Symbol)
    for (k, v) in pairs(cols(t))
        v[i] = getproperty(val, k)
    end
    val
end

Base.iterate(t::CTable, i=1) = i > length(t) ? nothing : (t[i], i + 1)
Base.eltype(t::CTable) = NamedTuple{names(t), Tuple{map(eltype, cols(t))...}}

Base.merge(a::CTable, b::CTable) = CTable(merge(cols(a), cols(b)))

function Base.sort!(a::CTable, sel=All())
end

function fake(::Type{CTable}, n = 10_000)
    data = (
        x1 = rand(1:10, n),
        x2 = rand(n),
        x3 = rand(Bool, n),
        y1 = rand(1f0:10f0, n),
        y2 = rand(["A","B","C"], n),
        y3 = rand(1:10, n),
        z1 = [rand(Bool) ? missing : rand() for i in 1:n]
    )
    CTable(data)
end

#-----------------------------------------------------------------------------# filter
Base.filter(f, t::CTable, sel=All()) = filter!(f, copy(t), sel)

function Base.filter!(f, t::CTable, sel=All())
    idx = map(!f, select(t, sel))
    map(x -> deleteat!(x, idx), cols(t))
    t
end

#-----------------------------------------------------------------------------# drop missing
Base.filter!(::Not{Missing}, t::CTable, sel=All()) = error("No filter! method.  Did you mean filter?")

function Base.filter(f::Not{Missing}, t::CTable, sel=All())
    temp = select(t, sel)
    idx = map(x -> !any(ismissing, x), temp)
    temp2 = temp[idx]
    new_cols = map(cols(temp2)) do col
        T = eltype(col)
        if Missing <: T
            new_col = similar(col, nonmissingtype(T))
            new_col[:] = col
            new_col
        else
            col
        end
    end
    CTable(merge(cols(t[idx]), new_cols))
end

#-----------------------------------------------------------------------------# show
function Base.show(io::IO, t::CTable)
    n = length(t)
    printstyled(io, "CTable ", bold=true)
    printstyled(io, "($n × $(ncols(t)))"; color=:light_black)

    isempty(t) && return

    nr, nc = displaysize(io)
    println(io)

    rows_omitted = length(t) - nr + 6

    cols2show = Vector{String}[]

    width = 0
    for (k,v) in pairs(cols(t))
        vals = string.(vcat(k, v[1:min(n, nr - 6 - (rows_omitted > 0))]))
        width += maximum(length, vals) + 2
        width <= nc - 2 ? push!(cols2show, vals) : break
    end

    for c in cols2show
        w = maximum(length, c) + 1
        h = c[1]
        printstyled(io, "|"; color = :light_black)
        printstyled(io, " $h"; color = :light_green)
        printstyled(io, ' ' ^ (w - length(h)); color = :light_green)
    end
    printstyled(io, " |"; color = :light_black)

    println(io)

    for c in cols2show
        w = maximum(length, c) + 1
        h = c[1]
        printstyled(io, "|" * '-' ^ (w + 1); color = :light_black)
    end
    printstyled(io, "-|"; color = :light_black)

    println(io)

    for i in 2:length(cols2show[1])
        for c in cols2show
            w = maximum(length, c)
            printstyled(io, "|"; color = :light_black)
            print(io, " $(c[i]) ")
            print(io, ' ' ^ (w -= length(c[i])))
        end
        printstyled(io, " |"; color = :light_black)
        println(io)
    end

    rows_omitted > 0 && printstyled(io, "  (Rows Omitted: ", rows_omitted, ")", color=:cyan)

    cols_omitted = length(cols(t)) - length(cols2show)

    cols_omitted > 0 && printstyled(io, "  (Cols Omitted: ", cols_omitted, ")", color=:cyan)
end


text:
# Autogenerated wrapper script for GRAMPC_jll for armv7l-linux-musleabihf
export libgrampc_julia

JLLWrappers.@generate_wrapper_header("GRAMPC")
JLLWrappers.@declare_library_product(libgrampc_julia, "libgrampc_julia.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libgrampc_julia,
        "lib/libgrampc_julia.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Testing part 1

include("../scripts-part1/part1.jl") # modify to include the correct script

# Add unit and reference tests

# Reference test using ReferenceTests.jl
"Compare all dict entries"
comp(d1, d2) = keys(d1) == keys(d2) && all([ isapprox(v1, v2; atol = 1e-8) for (v1,v2) in zip(values(d1), values(d2))])
inds = Int.(ceil.(LinRange(1, length(xc), 12)))
d = Dict(:X=> xc[inds], :H=>H[inds, inds, 15])

@testset "Ref-file" begin
    @test_reference "reftest-files/test_1.bson" d by=comp
end