text:
# functions for calculating the node coordinates and metrics of straight
# sided elements (ie. meshes with linear coordinate fields)

#------------------------------------------------------------------------------
# allocators for straight sides meshes
"""
  This function allocates the arrays that store the mesh node coordinates,
  the mesh vertex coordinates and assigns them to the corresponding field of
  the mesh object.
"""
function allocateCoordinateArrays(mesh::PumiMeshDG{Tmsh},
                                           sbp::AbstractSBP) where Tmsh

  num_coord_nodes = mesh.coord_numNodesPerElement

  if !isFieldDefined(mesh, :coords, :vert_coords)
    mesh.coords = Array{Float64}(mesh.dim, sbp.numnodes, mesh.numEl)
    mesh.vert_coords = Array{Float64}(mesh.dim, num_coord_nodes, mesh.numEl)
    mesh.vert_coords_bar = zeros(mesh.vert_coords)

  else
    fill!(mesh.coords, 0.0)
    fill!(mesh.vert_coords, 0.0)
  end

  return nothing
end

"""
  This function allocates the arrays that store the dxidx and jac and
  assigns them to the corresponding field of the mesh object.
"""
function allocateMetricsArrays(mesh::PumiMeshDG{Tmsh}, sbp::AbstractSBP) where Tmsh

  if !isFieldDefined(mesh, :jac, :dxidx)
    mesh.dxidx = Array{Tmsh}(mesh.dim, mesh.dim, sbp.numnodes, mesh.numEl)
    mesh.jac = Array{Tmsh}(sbp.numnodes, mesh.numEl)
  else
    fill!(mesh.dxidx, 0.0)
    fill!(mesh.jac, 0.0)
  end
  mesh.dxidx_bar = zeros(mesh.dxidx)
  mesh.jac_bar = zeros(mesh.jac)

  return nothing
end

#------------------------------------------------------------------------------
# straight-sided coordinate and metric functions

"""
  This function gets the coordinates of both the mesh coordinate field and the
  node field and stores them into the member fields of the mesh object.
  It uses a smart allocator to allocate
  these arrays needed, or avoid doing so if they have already been allocated.
  See allocateCoordinateArrays for which arrays are populated by
  this function.  Non curvilinear meshes only.
"""
#TODO: stop using slice notation
# can be generalized with numVertsPerElement
function getCoordinates(mesh::PumiMeshDG, sbp::AbstractSBP)
# populate the coords array of the mesh object

#  mesh.coords = Array{Float64}(mesh.dim, sbp.numnodes, mesh.numEl)
#  mesh.vert_coords = Array{Float64}(mesh.dim, nvert_per_el, mesh.numEl)

  allocateCoordinateArrays(mesh, sbp)

  nvert_per_el = mesh.numTypePerElement[1]
  @assert size(mesh.coords) == (mesh.dim, sbp.numnodes, mesh.numEl)
  @assert size(mesh.vert_coords) == (mesh.dim, nvert_per_el, mesh.numEl)
  @assert mesh.coord_order == 1

  #println("entered getCoordinates")
  numVertsPerElement = mesh.numTypePerElement[1]
  coords_i = zeros(mesh.dim ,numVertsPerElement)
  coords_it = zeros(numVertsPerElement, mesh.dim)
  for i=1:mesh.numEl  # loop over elements

    el_i = mesh.elements[i]
    apf.getAllEntityCoords(mesh.m_ptr, el_i, coords_i)
    mesh.vert_coords[:, :, i] = coords_i[:, :]
    mesh.coords[:, :, i] = vertToVolumeCoords(mesh, sbp, coords_i)
#    coords_it[:,:] = coords_i[1:mesh.dim, :].'
#    mesh.coords[:, :, i] = SummationByParts.SymCubatures.calcnodes(sbp.cub, coords_it)
  end

  return nothing

end

function vertToVolumeCoords(mesh, sbp, vert_coords)
  coords_it = vert_coords[1:mesh.dim, :].'
  coords = SummationByParts.SymCubatures.calcnodes(sbp.cub, coords_it)
  return coords
end

"""
  This function calculates the dxidx and jac values for the entire mesh.  It
  uses a smart allocator to allocate the arrays if needed.  Non-curvilinear
  meshes only.
"""
function getMetrics(mesh::PumiMeshDG, sbp::AbstractSBP)

  allocateMetricsArrays(mesh, sbp)

#  if mesh.dim == 2
#    fill!(mesh.coords, 0.0)  # when using calcMappingJacobian, this isn't needed
#    @assert size(mesh.vert_coords, 2) == 3
#    coord_order = 1
#    calcMappingJacobian!(sbp, coord_order, sbp.vtx.', mesh.vert_coords,
#                         mesh.coords, mesh.dxidx, mesh.jac)
#  else
    mappingjacobian!(sbp, mesh.coords, mesh.dxidx, mesh.jac)
#  end

  return nothing
end

"""
  This function uses reverse mode differentiation to back propigate dxidx_bar
  to vert_coords_bar

  Inputs/Outputs:
    mesh: a DG mesh object, must be 2D.  The mesh.vert_coords_bar field is
          incremented (ie. not overwritten) by this function.

    sbp: an SBP operator
"""
function getVertCoords_rev(mesh::PumiMeshDG{Tmsh}, sbp) where Tmsh
#  @assert mesh.dim == 2
  @assert mesh.coord_order == 1

#  if mesh.dim == 2
#    #TODO: go back to old mappingjacobian!
#    coord_order = 1
#    coords_bar = zeros(Tmsh, size(mesh.coords)...)  # we don't allow perturbing the mesh coords
#    Eone_bar = zeros(Tmsh, mesh.numNodesPerElement, mesh.dim, mesh.numEl)
#    calcMappingJacobianRevDiff!(sbp, coord_order, sbp.vtx.',
#                            mesh.vert_coords_bar, coords_bar, mesh.dxidx,
#                            mesh.dxidx_bar,
#                            mesh.jac, mesh.jac_bar,
#                            Eone_bar)
#
#    else

      coords_bar = zeros(Tmsh, size(mesh.coords)...)
      mappingjacobian_rev!(sbp, mesh.coords, coords_bar, mesh.dxidx_bar,
                           mesh.jac_bar)
      volumeToVertCoords_rev(mesh, sbp, mesh.vert_coords_bar, coords_bar)
#    end

  return coords_bar
end

"""
  Back propigates the adjoint part of the volume coordinates to the adjoint
  part of the vertex coordinates

  Inputs/Outputs:
    mesh: a DG mesh object. mesh.vert_coords_bar is updated with the result
    sbp: an SBP operator
    vertcoords_bar: updated with the results
    coords_bar: the adjoint part of the volume node coordinates (input)
"""
function volumeToVertCoords_rev(mesh::PumiMeshDG, sbp, vertcoords_bar::Abstract3DArray,  coords_bar::Abstract3DArray)

  #TODO; check that face/boundary/shareface coordinates are back propigated
  #      to the volume coordinates
  # because the mapping is linear, compute the jacobian using complex step
  # and apply it (transposed) to each element

  numVertPerElement = mesh.numTypePerElement[1]
  jac = zeros(mesh.dim*mesh.numNodesPerElement, numVertPerElement*mesh.dim)
  coords_i = zeros(Complex128, mesh.dim, numVertPerElement)
#  coords_it = zeros(Complex128, numVertPerElement, mesh.dim)
  # coords_out = mesh.dim x mesh.numNodesPerElement
  h = 1e-20
  pert = Complex128(0, h)

  # make a complex SBP operator because we need it for calcnodes
  #TODO: make this work for 3D
#  sbp_complex = TriSBP{Complex128}(degree=sbp.degree, reorder=sbp.reorder, internal=sbp.internal, vertices=sbp.vertices)

  for i=1:length(coords_i)
    coords_i[i] += pert
    coords_it = coords_i.'
    coords_out = SummationByParts.SymCubatures.calcnodes(sbp.cub, coords_it)
    for j=1:length(coords_out)
      jac[j, i] = imag(coords_out[j])/h
    end
    coords_i[i] -= pert
  end

  # now apply it
  coords_bar_reshaped = zeros(numVertPerElement* mesh.dim)
  vertcoords_bar_reshaped = zeros(numVertPerElement*mesh.dim)
  for i=1:mesh.numEl
    coords_bar_i = sview(coords_bar, :, :, i)
    vertcoords_bar_i = sview(vertcoords_bar, :, :, i)
    # make coords_bar into a vector
    for j=1:length(coords_bar_i)
      coords_bar_reshaped[j] = coords_bar_i[j]
    end

    smallmatTvec!(jac, coords_bar_reshaped, vertcoords_bar_reshaped)

    # move back to vertcoords_bar
    for j=1:length(vertcoords_bar_i)
      vertcoords_bar_i[j] += vertcoords_bar_reshaped[j]
    end
  end  # end loop i

  return nothing
end

"""
  This function dispatches on the type of the mesh to get the coordinates of
  a mesh element (a face in 2d or a region in 3d)
"""
function getElementCoords(mesh::PumiMesh2D, entity::Ptr{Void}, coords::AbstractMatrix)
  # coords must be 3 x numVertsPerElement
  sx, sy = size(coords)
  apf.getFaceCoords(mesh.m_ptr, entity, coords, sx, sy)
end

function getElementCoords(mesh::PumiMesh3D, entity::Ptr{Void}, coords::AbstractMatrix)

  sx, sy = size(coords)
  apf.getElCoords(mesh.m_ptr, entity, coords, sx, sy)
end

#TODO: stop using slice notation
"""
  This function calculates the node coordinates, dxidx, jac, for a
  CG mesh
"""
function getCoordinatesAndMetrics(mesh::PumiMeshCG{Tmsh}, sbp::AbstractSBP) where Tmsh
# populate the coords array of the mesh object
mesh.coords = Array{Float64}(mesh.dim, sbp.numnodes, mesh.numEl)

#println("entered getCoordinates")

numVertsPerElement = mesh.numTypePerElement[1]
coords_i = zeros(3,numVertsPerElement)
coords_it = zeros(numVertsPerElement, mesh.dim)
for i=1:mesh.numEl  # loop over elements
  el_i = mesh.elements[i]
  getElementCoords(mesh, el_i, coords_i)

  coords_it[:,:] = coords_i[1:mesh.dim, :].'
  mesh.coords[:, :, i] = calcnodes(sbp, coords_it)
end

  mesh.dxidx = Array{Tmsh}(mesh.dim, mesh.dim, sbp.numnodes, mesh.numEl)
  mesh.jac = Array{Tmsh}(sbp.numnodes, mesh.numEl)
  mappingjacobian!(sbp, mesh.coords, mesh.dxidx, mesh.jac)

return nothing

end

"""
  Populates the input array with the coordinates of the nodes on list of
  boundary faces.  Non-curvilinear meshes only.  In 3D it uses mesh.topo
  to get the vertices of the correct order.

  Inputs:
    mesh: a mesh object
    bndryfaces: an array of Boundary objects to calculate the coordinates of

  Inputs/Outputs:
    coords_bndry: an array dim x numfacenodes x length(bndryfaces) to populate
                  with the coordinates of the nodes on the faces
"""
function getBndryCoordinates(mesh::PumiMeshDG2{Tmsh},
                       bndryfaces::Array{Boundary},
                       coords_bndry::Array{Tmsh, 3}) where Tmsh
# calculate the coordinates on the boundary for the specified faces
# and store in coords_bndry

#  println("----- Entered getBndryCoordinates -----")
  sbpface = mesh.sbpface

  coords_i = zeros(3, 3)
  coords_it = zeros(3, 2)
  coords_edge = zeros(2, 2)

  for i=1:length(bndryfaces)
    bndry_i = bndryfaces[i]

    el = bndry_i.element
    el_ptr = mesh.elements[el]
    face = bndry_i.face

    sizex, sizey = size(coords_i)
    apf.getFaceCoords(mesh.m_ptr, el_ptr, coords_i, sizex, sizey)

    coords_it[:, :] = coords_i[1:2, :].'

    # extract the needed vertex coords
#    v1 = facemap[1, face]
#    v2 = facemap[2, face]
    v1 = face
    v2 = mod(face,3) + 1
    coords_edge[1, 1] = coords_it[v1, 1]
    coords_edge[1, 2] = coords_it[v1, 2]
    coords_edge[2, 1] = coords_it[v2, 1]
    coords_edge[2, 2] = coords_it[v2 ,2]

    coords_bndry[:, :, i] = SummationByParts.SymCubatures.calcnodes(sbpface.cub, coords_edge)

  end

end

function getBndryCoordinates(mesh::PumiMesh3DG{Tmsh}, bndryfaces::Array{Boundary}, coords_bndry::Array{Tmsh, 3}) where Tmsh

#  println("----- entered getBndryCoordinates -----")
  sbpface = mesh.sbpface
  coords_i = zeros(3, mesh.numEntitiesPerType[1])
#  coords_i = zeros(3, mesh.numTypePerElement[1])
  coords_face = zeros(3, 3)
  vertmap = mesh.topo.face_verts
  for i=1:length(bndryfaces)
    bndry_i = bndryfaces[i]
    el = bndry_i.element
    el_ptr = mesh.elements[el]
    face = bndry_i.face

    getElementCoords(mesh, el_ptr, coords_i)

    # extract the vertices of the face
    for v=1:3
      vidx = vertmap[v, face]
      for dim=1:3
        coords_face[v, dim] = coords_i[dim, vidx]
      end
    end

    coords_bndry[:, :, i] = SummationByParts.SymCubatures.calcnodes(sbpface.cub, coords_face)
  end

  return nothing
end

"""
  This is basically the same as getBndryCoordinates, except it operates on
  interfaces.  It calculates all quantities from the perspective of elementL

  I think this can be combined with getBndryCoordinates
"""
function getInterfaceCoordinates(mesh::PumiMeshDG2{Tmsh},
                             bndryfaces::Array{Interface},
                             coords_bndry::Array{Tmsh, 3}) where Tmsh
# calculate the coordinates on the boundary for the specified faces
# and store in coords_bndry

#  println("----- Entered getBndryCoordinates -----")
  sbpface = mesh.sbpface

  coords_i = zeros(3, 3)
  coords_it = zeros(3, 2)
  coords_edge = zeros(2, 2)  # one point per row

  for i=1:length(bndryfaces)
    bndry_i = bndryfaces[i]

    el = bndry_i.elementL
    el_ptr = mesh.elements[el]
    face = bndry_i.faceL

    sizex, sizey = size(coords_i)
    apf.getFaceCoords(mesh.m_ptr, el_ptr, coords_i, sizex, sizey)

    coords_it[:, :] = coords_i[1:2, :].'

    # extract the needed vertex coords
#    v1 = facemap[1, face]
#    v2 = facemap[2, face]
    v1 = face
    v2 = mod(face,3) + 1
    coords_edge[1, 1] = coords_it[v1, 1]
    coords_edge[1, 2] = coords_it[v1, 2]
    coords_edge[2, 1] = coords_it[v2, 1]
    coords_edge[2, 2] = coords_it[v2 ,2]

    coords_bndry[:, :, i] = SummationByParts.SymCubatures.calcnodes(sbpface.cub, coords_edge)

  end

end

function getInterfaceCoordinates(mesh::PumiMesh3DG{Tmsh}, bndryfaces::Array{Interface}, coords_bndry::Array{Tmsh, 3}) where Tmsh

#  println("----- entered getInterfaceCoordinates -----")
  sbpface = mesh.sbpface
  coords_i = zeros(3, mesh.numTypePerElement[1])
  coords_face = zeros(3, 3)
  vertmap = mesh.topo.face_verts
  for i=1:length(bndryfaces)
    bndry_i = bndryfaces[i]
    el = bndry_i.elementL
    el_ptr = mesh.elements[el]
    face = bndry_i.faceL

    getElementCoords(mesh, el_ptr, coords_i)

    # extract the vertices of the face
    for v=1:3
      vidx = vertmap[v, face]
      for dim=1:3
        coords_face[v, dim] = coords_i[dim, vidx]
      end
    end

    coords_bndry[:, :, i] = SummationByParts.SymCubatures.calcnodes(sbpface.cub, coords_face)
  end

  return nothing
end

"""
  This function calculates the face normal vectors for all interfaces,
  boundary faces, and shared faces.  A smart allocator is used to allocate
  the arrays if needed.  Non curvilinear meshes only.

  The interpolated dxidx arrays must be populated before this function is
  called.

  Inputs:
    mesh: a PumiMeshdG
    sbp: an SBP operator
"""
function getFaceNormals(mesh::PumiMeshDG, sbp)

  allocateNormals(mesh, sbp)

  calcFaceNormal(mesh, sbp, mesh.bndryfaces, mesh.dxidx_bndry, mesh.nrm_bndry)
  calcFaceNormal(mesh, sbp, mesh.interfaces, mesh.dxidx_face, mesh.nrm_face)
  for i=1:mesh.npeers
    calcFaceNormal(mesh, sbp, mesh.bndries_local[i], mesh.dxidx_sharedface[i],
                   mesh.nrm_sharedface[i])
  end

  return nothing
end

"""
  This function calculates the face normal for a set of faces given the
  dxidx values for each face node and the facenormal in parametric space

  Inputs:
    mesh: a DG mesh
    sbp: an SBP operator containing the face normals in parametric space
    faces: list of Interfaces or Boundaries, should be dim x numFacesPerElement
    dxidx: d(xi)/d(x) at the face nodes, should be dim x dim x numfaceNodes x
           length(faces)

  Inputs/Outputs:
    nrm: array, dim x numfacenodes x length(faces) containing the face normal
         vector in x-y space at each face node of each face

  Aliasing restrictions: none
"""
function calcFaceNormal(
mesh::PumiMeshDG, sbp, faces::AbstractArray{Tbndry, 1},
dxidx::AbstractArray{Tmsh, 4},
nrm::AbstractArray{Tmsh, 3}) where {Tmsh, Tbndry <: Union{Boundary, Interface}}

  nfaces = length(faces)
  Tdim = mesh.dim
  for i=1:nfaces
    faceL = getFaceL(faces[i])

    for j=1:mesh.numNodesPerFace
      for dim=1:Tdim
        nrm_d = zero(Tmsh)
        for d=1:Tdim
          nrm_d += mesh.sbpface.normal[d, faceL]*dxidx[d, dim, j, i]
        end  # end loop d
        nrm[dim, j, i] = nrm_d
      end  # end loop dim
    end  # end loop j

  end  # end loop i

  return nothing
end


text:
# ---
# title: 69. Sqrt(x)
# id: problem69
# author: Indigo
# date: 2022-03-03
# difficulty: Easy
# categories: Math, Binary Search
# link: <https://leetcode.com/problems/sqrtx/description/>
# hidden: true
# ---
#
# Given a non-negative integer `x`, compute and return _the square root of_ `x`.
#
# Since the return type is an integer, the decimal digits are **truncated** ,
# and only **the integer part** of the result is returned.
#
#
#
# **Example 1:**
#
#
#
#     Input: x = 4
#     Output: 2
#
#
# **Example 2:**
#
#
#
#     Input: x = 8
#     Output: 2
#     Explanation: The square root of 8 is 2.82842..., and since the decimal part is truncated, 2 is returned.
#
#
#
# **Constraints:**
#
#   * `0 <= x <= 231 - 1`
#
#
## @lc code=start
using LeetCode

my_sqrt = isqrt
function mysqrt(x::Int)::Int
    x < 0 && throw(DomainError(x))
    left, right = 0, x
    while left <= right
        mid = left + ((right - left) >> 1)
        if mid ^ 2 > x
            right = mid - 1
        else
            left = mid + 1
        end
    end
    right
end

## @lc code=end

ext: jl


text:
# Get show compatible function with var set
printpolyfun(var) = (io, p, mimetype = MIME"text/plain"()) -> printpolydesc(io, p, var, mimetype)

# Specialized function, copied and adapted from Polynomials.jl
# Ignores var as set in Poly
""" Prints polynomial in descending order, with variable `var`
"""
function printpolydesc(io::IO, p::Poly{T}, var, mimetype = MIME"text/plain"()) where {T}
    printpoly2(io, p, var, mimetype)
end

# Specialized function, copied and adapted from Polynomials.jl
# Ignores var as set in Poly
function printpoly2(io::IO, p::Poly{T}, var, mimetype) where {T}
    first = true
    printed_anything = false
    for i in reverse(eachindex(p))
        printed = showterm2(io, p, i, first, var, mimetype)
        first &= !printed
        printed_anything |= printed
    end
    printed_anything || print(io, zero(T))
    return nothing
end

#Specialized function, copied and adapted from Polynomials.jl
# Ignores var as set in Poly
function showterm2(io::IO, p::Poly{T}, j, first, var, mimetype) where {T}
    pj = p[j]

    pj == zero(T) && return false

    pj = Polynomials.printsign(io, pj, first, mimetype)
    Polynomials.printcoefficient(io, pj, j, mimetype)
    Polynomials.printproductsign(io, pj, j, mimetype)
    Polynomials.printexponent(io, var, j, mimetype)
    true
end


text:
function LiesWithin(project,activeTile,e,n)
    iActiveTile = findall(project.tileRegister.tileRefs.==activeTile)
    gridRefLimits = Dict(
        "easting_min" => minimum(project.tileRegister.eastings_min[iActiveTile]),
        "easting_max" => maximum(project.tileRegister.eastings_max[iActiveTile]),
        "northing_min" => minimum(project.tileRegister.northings_min[iActiveTile]),
        "northing_max" => maximum(project.tileRegister.northings_max[iActiveTile]))
    iKeep = intersect(
        findall(e .>= gridRefLimits["easting_min"]),
        findall(e .<  gridRefLimits["easting_max"]),
        findall(n .>= gridRefLimits["northing_min"]),
        findall(n .<  gridRefLimits["northing_max"])
    )
end

#-------------------------------------------------------------------------------
function ExtractData(project,DataSource,activeTiles,selected_cols;method="gridref")

	# Find file names (from "<ProjectDir>/.DataFileLinks/<DataSource>.txt")
	fnames = DataFileName(DataSource)
	fnames = unique([fnames[tile] for tile in activeTiles])

	# Extract csv data for all returned files, linked to active tiles
	tmp = [CSV.read(
		joinpath(env.paths[:data_sensitive],string(fname)),
		DataFrame, normalizenames=true)
	for fname in fnames]

	# Merge all returned dataframes
	tmp = vcat(tmp...)

	# Find points within active tiles
	if method == "gridref"
		iKeep = sort(vcat(LiesWithin.((project,),activeTiles,(tmp.Easting,),(tmp.Northing,))...))
	elseif method == "uprn"
		iKeep = 1:nrow(tmp)
	end

	# Push to project struct (selected cols only)
	push!(project.dat,DataSource=>tmp[iKeep,selected_cols])

	return project
end

#-------------------------------------------------------------------------------
"""
`addFields` should be a dictionary, with keys named as target dataframe columns
from GML, values as a Tuple of lenght 2: first value DataType, second value
default placeholders (for `fill!(Array{}...)`).

Example:

```
addFields = Dict(
    "eastings"=>(Float64,[]),
    "northings"=>(Float64,[]),
)
```
"""
function AddFields(project,addFields)
    fields = string.(keys(addFields))

    for i in 1:length(fields)
        # Template array
        tmp = fill!(Array{
            addFields[fields[i]][1]
        }(undef,nrow(project.dat["master"])),
            addFields[fields[i]][2]
        )

        for k in 1:length(project.dat["gml"])
            # Find data rows in master dataframe, relating to current tile
            iTile = intersect(
				findall(project.dat["master"].osgb_tile .== project.dat["gml"][k].tile),
				findall(project.dat["master"].iGml.!=0)
			)
            # Find corresponding entries in GML
            iGml = project.dat["master"].iGml[iTile]

			coldata = project.dat["gml"][k].summary[!,fields[i]][iGml]
            if fields[i] == "eastings" || fields[i] == "northings"
				arraypos1(coldata) = coldata[1]
				coldata = arraypos1.(coldata)
			end
			tmp[iTile] .= coldata
        end

        project.dat["master"][!,"$(fields[i])_GML"] = tmp
    end
end

#-------------------------------------------------------------------------------
function LoadStockData(env,project)
	# Send compiled stock data to JSON (by tile)
	function WriteStockJSON(project,tile)
		data_byTile = project.dat["master"][project.dat["master"].osgb_tile.==tile,:]
		json_string = JSON.json(data_byTile)

		writeJSON(json_string,"BuildingStock",tile)
	end

	# List active tiles
	activetiles(tile) = tile.tile
	activeTiles = activetiles.(project.dat["gml"])

	# List saved stock JSONs
	BuildingStockDir = readdir(joinpath(env.paths[:projects],"BuildingStock"))

	# Find relevant saved JSONs
	ExistingTileJSONs(activeTile,BuildingStockDir) = sum(occursin.(activeTile,BuildingStockDir))>0
	iKeep = findall(vcat(ExistingTileJSONs.(activeTiles,(BuildingStockDir,))...))

	# Restore existing JSONs (if any JSON are found)
	stockData = []
	if !isempty(iKeep)
		# Try/catch rejects the saved JSONs when multiple tiles are used AND the
		# ... DataFrame columns don't match. Rebuild carried out at next step
		try
			push!(stockData,vcat(readJSON.("BuildingStock",activeTiles[iKeep])...))
		catch
			iKeep = []
		end
	end

	# Evaluate new stock data, where existing JSONs are not found/don't match
	if length(iKeep)!=length(activeTiles)
		project = LinkHaData(env,project,activeTiles[Not(iKeep)])
		project = link_epcdata(env,project,activeTiles[Not(iKeep)])
		push!(stockData,project.dat["master"])
		WriteStockJSON.((project,),activeTiles[Not(iKeep)])
	end

	# Check that column names match, rerun if necessary, then merge all stock data
	if length(stockData).==2
		MatchingDfCols(col₁,cols₂) = sum(cols₂.==col₁)
		cols = sort.(names.(stockData))
		if sum(MatchingDfCols.(cols[1],(cols[2],)))!=length(cols[1])
			project = LinkHaData(env,project,activeTiles[iKeep])
			project = link_epcdata(env,project,activeTiles[iKeep])
			stockData[1] = deepcopy(project.dat["master"])
			project.dat["master"] = vcat(stockData...)
			WriteStockJSON.((project,),activeTiles[iKeep])
		end
	end

	project.dat["master"] = vcat(stockData...)

	return project
end

#-------------------------------------------------------------------------------
function LinkHaData(env,project,activeTiles)

	# Link HA data
	DataSource = "HA"
	selected_cols = [:UPRN,:Easting,:Northing,:Full_postal_address,
		:Building_block_ID,:Data_zone_code,:Data_zone,:Property_type,]
	project = ExtractData(project,DataSource,activeTiles,selected_cols;method="gridref")

	# Create new master dataframe from HA data
	push!(project.dat,"master"=>deepcopy(project.dat["HA"]))
	rename!(project.dat["master"],:Easting=>:Easting_HA)
	rename!(project.dat["master"],:Northing=>:Northing_HA)

	# Geo-reference master data against GIS data
	project.dat["master"][!,"osgb_tile"] = fill!(Array{String}(undef,size(project.dat["master"],1)),"")
	iTiles = LiesWithin.((project,),activeTiles,(project.dat["HA"].Easting,),(project.dat["HA"].Northing,))
	IMap(iTile,activeTile) = project.dat["master"].osgb_tile[iTile].=activeTile
	IMap.(iTiles,activeTiles)
	project.dat["master"][!,"iGml"] = fill!(Array{Int64}(undef,size(project.dat["master"],1)),0)
	project.dat["master"][!,"osgb"] = fill!(Array{Int64}(undef,size(project.dat["master"],1)),0)
	[GeoRef(project,k) for k in 1:length(project.dat["gml"])]

	# Introduce other HA data fields (Tuple of (DataType,<placeholder value>))
	addFields = Dict(
	    "eastings"=>(Array{Float64},[]),
	    "northings"=>(Array{Float64},[]),
	)
	AddFields(project,addFields)

	# Remove redunant vertices and scan for adjacencies
	project = ProcessGeom(project)

	# Eval polygon areas
	i_linked = findall(isempty.(project.dat["master"].eastings_GML).==false)
	shoelacearea_closedpoly(x,y) = shoelacearea(x[1:end-1],y[1:end-1])
	project.dat["master"][!,:OS_area_eval] = fill!(Array{Union{Missing,Float32}}(undef,nrow(project.dat["master"])),missing)
	project.dat["master"].OS_area_eval[i_linked] = shoelacearea_closedpoly.(
		project.dat["master"].eastings_GML[i_linked],
		project.dat["master"].northings_GML[i_linked]
    )

	# Integrate Ordnance Survey 'Building Height Attribute' data
	project = linkbuildingheights(env,project,activeTiles)

	# Find plans that have more than one UPRN under same roof, generate unique string of all UPRNs
	n_uprns_undersameroof, uprns_undersameroof = multidwellplans(project.dat["master"][:,[:UPRN,:osgb,:iGml]])
	project.dat["master"][!,:nUPRNs_undersameroof] = n_uprns_undersameroof
	project.dat["master"][!,:UPRNs_undersameroof] = uprns_undersameroof

	return project
end


text:
module ETOPO

using OceanGrids
using Dates
using Unitful
using DataDeps              # For storage location of data
using NCDatasets
using Distances

function fallback_download(remotepath, localdir)
    @assert(isdir(localdir))
    filename = basename(remotepath)  # only works for URLs with filename as last part of name
    localpath = joinpath(localdir, filename)
    Base.download(remotepath, localpath)
    return localpath
end

# Create registry entry for ETOPO
function register_ETOPO()
    register(
        DataDep(
            "ETOPO",
            """
            Reference:
            $(citation()) [accessed $(Dates.today())]
            """,
            "https://www.ngdc.noaa.gov/mgg/global/relief/ETOPO1/data/bedrock/grid_registered/netcdf/ETOPO1_Bed_g_gdal.grd.gz",
            sha2_256,
            fetch_method = fallback_download,
            post_fetch_method = unpack
        )
    )
    return nothing
end

"""
    Z, lats, lons = load()

Returns the fine resolution topography from ETOPO.
"""
function load()
    register_ETOPO()
    nc_file = @datadep_str string("ETOPO/ETOPO1_Bed_g_gdal.grd")
    ds = Dataset(nc_file)
    Z, lats, lons, dims =
    Dataset(nc_file,"r") do ds
        dims = Tuple(ds["dimension"][:])
        ds["z"][:], range(ds["y_range"]..., length=dims[2]),
        range(ds["x_range"]..., length=dims[1]), dims
    end # ds is closed
    @info """You are about to use the ETOPO data set.
          If you use it for research, please cite:

          - $(citation())

          You can find the corresponding BibTeX entries in the CITATION.bib file
          at the root of the AIBECS.jl package repository.
          (Look for the "Amante_Eakins_2009" key.)
          """
    return reverse(reshape(Z, dims), dims=2), lats, lons # TODO
end

citation() = "Amante, C. and B.W. Eakins, 2009. ETOPO1 1 Arc-Minute Global Relief Model: Procedures, Data Sources and Analysis. NOAA Technical Memorandum NESDIS NGDC-24. National Geophysical Data Center, NOAA. doi:10.7289/V5C8276M [access date]."

"""
    bintopo(Z, lat, lon, grd)

Bins the depths `Z` into at locations `(lat,lon)` onto `grd`.
To be used with the ETOPO dataset to dertermine the amount of subgrid topography in each box area
"""
function bintopo(Z::AbstractArray{T,2}, lats, lons, grd) where T
    bin3D = zeros(size(grd))
    # Use the top edge to invert the functions lon[i], lat[j], and depth[k]
    # I.e., to find the index in grd for a given lon, lat, depth
    edepth = ustrip.(grd.depth .+ 0.5grd.δdepth)
    elon = ustrip.(grd.lon .+ 0.5grd.δlon)
    elat = ustrip.(grd.lat .+ 0.5grd.δlat)
    kmax = dropdims(count(iswet(grd), dims=3), dims=3)
    for (j,lat) in enumerate(lats)
        jgrd = searchsortedfirst(elat, lat)
        for (i,lon) in enumerate(lons)
            igrd = searchsortedfirst(elon, mod(lon, 360))
            kmax[jgrd,igrd] == 0 && continue
            kgrd = searchsortedfirst(view(edepth, 1:kmax[jgrd, igrd]-1), -Z[i,j])
            bin3D[jgrd, igrd, kgrd] += 1
        end
    end
    return bin3D
end

"""
    fractiontopo(grd)

Returns the fraction of subgrid sediments in each wet box of grd.

Note that (i) points located above sea level are counted in the top layer,
and (ii) points located in a dry box or below are counted in the deepest
wet box.
"""
function fractiontopo(grd)
    @warn "Binning ETOPO to grd. This may take a few seconds"
    Z, lats, lons = load()
    return fractiontopo(bintopo(Z, lats, lons, grd), grd)
end
fractiontopo(bin3D, grd) = vectorize(bin3D ./ repeat(sum(bin3D, dims=3), outer=(1, 1, size(grd)[3])), grd)
export fractiontopo

function roughnesstopo(grd)
    @warn "Estimating roughness from ETOPO. This may take a few seconds"
    Z, lats, lons = load()
    R_earth = ustrip(u"m", 6371.0u"km")
    δlon = step(lons)
    δlat = step(lons)
    midlats = 0.5(lats[1:end-1] + lats[2:end])
    midlons = 0.5(lons[1:end-1] + lons[2:end])
    δxs = [Haversine(R_earth)((0.0, lat), (δlon, lat)) for lat in midlats]
    δy = Haversine(R_earth)((lons[1], 0.0), (lons[1], δlat))
    A_lon = abs.(diff(Z, dims=1)) .* δy
    Z_lon = 0.5(Z[1:end-1,:] + Z[2:end,:])
    A_lat = abs.(diff(Z, dims=2)) .* repeat(reshape(δxs, (1, length(δxs))), outer=(size(Z,1), 1))
    Z_lat = 0.5(Z[:,1:end-1] + Z[:,2:end])
    return sumroughness(A_lon, Z_lon, lats, midlons, grd) + sumroughness(A_lat, Z_lat, midlats, lons, grd)
end

function sumroughness(A::AbstractArray{T,2}, Z, lats, lons, grd) where T
    roughness = zeros(size(grd))
    edepth = ustrip.(grd.depth .+ 0.5grd.δdepth)
    elon = ustrip.(grd.lon .+ 0.5grd.δlon)
    elat = ustrip.(grd.lat .+ 0.5grd.δlat)
    kmax = dropdims(count(iswet(grd), dims=3), dims=3)
    for (j,lat) in enumerate(lats)
        jgrd = searchsortedfirst(elat, lat)
        for (i,lon) in enumerate(lons)
            igrd = searchsortedfirst(elon, mod(lon, 360))
            kmax[jgrd,igrd] == 0 && continue
            kgrd = searchsortedfirst(view(edepth, 1:kmax[jgrd, igrd]-1), -Z[i,j])
            roughness[jgrd, igrd, kgrd] += A[i,j]
        end
    end
    return roughness[iswet(grd)]
end

end # module

export ETOPO


text:
using Test
import MadNLP: jac_structure!, hess_structure!, obj, grad!, cons!, jac_coord!, hess_coord!, jac_dense!, hess_dense!
using NLPModels
using LinearAlgebra
using MadNLPTests
using SparseArrays
using Random

function _compare_dense_with_sparse(kkt_system, n, m, ind_fixed, ind_eq)
    sparse_options = Dict{Symbol, Any}(
        :kkt_system=>MadNLP.SPARSE_KKT_SYSTEM,
        :linear_solver=>MadNLPLapackCPU,
        :print_level=>MadNLP.ERROR,
    )
    dense_options = Dict{Symbol, Any}(
        :kkt_system=>kkt_system,
        :linear_solver=>MadNLPLapackCPU,
        :print_level=>MadNLP.ERROR,
    )

    nlp = MadNLPTests.DenseDummyQP(; n=n, m=m, fixed_variables=ind_fixed, equality_cons=ind_eq)

    ips = MadNLP.InteriorPointSolver(nlp, option_dict=sparse_options)
    ipd = MadNLP.InteriorPointSolver(nlp, option_dict=dense_options)

    MadNLP.optimize!(ips)
    MadNLP.optimize!(ipd)

    # Check that dense formulation matches exactly sparse formulation
    @test ips.cnt.k == ipd.cnt.k
    @test ips.obj_val ≈ ipd.obj_val atol=1e-10
    @test ips.x ≈ ipd.x atol=1e-10
    @test ips.l ≈ ipd.l atol=1e-10
    @test ips.kkt.jac_com[:, 1:n] == ipd.kkt.jac
    if isa(ipd.kkt, MadNLP.AbstractReducedKKTSystem)
        @test Symmetric(ips.kkt.aug_com, :L) ≈ ipd.kkt.aug_com atol=1e-10
    end
end

@testset "MadNLP: API $(kkt_type)" for (kkt_type, kkt_options) in [
        (MadNLP.DenseKKTSystem, MadNLP.DENSE_KKT_SYSTEM),
        (MadNLP.DenseCondensedKKTSystem, MadNLP.DENSE_CONDENSED_KKT_SYSTEM),
    ]

    n = 10 # number of variables
    @testset "Unconstrained" begin
        dense_options = Dict{Symbol, Any}(
            :kkt_system=>kkt_options,
            :linear_solver=>MadNLPLapackCPU,
        )
        m = 0
        nlp = MadNLPTests.DenseDummyQP(; n=n, m=m)
        ipd = MadNLP.InteriorPointSolver(nlp, option_dict=dense_options)

        kkt = ipd.kkt
        @test isa(kkt, kkt_type)
        @test isempty(kkt.jac)
        # Special test for DenseKKTSystem
        if kkt_type <: MadNLP.DenseKKTSystem
            @test kkt.hess === kkt.aug_com
        end
        @test ipd.linear_solver.dense === kkt.aug_com
        @test size(kkt.hess) == (n, n)
        @test length(kkt.pr_diag) == n
        @test length(kkt.du_diag) == m

        # Test that using a sparse solver is forbidden in dense mode
        dense_options_error = Dict{Symbol, Any}(
            :kkt_system=>kkt_options,
            :linear_solver=>MadNLPUmfpack,
        )
        @test_throws Exception MadNLP.InteriorPointSolver(nlp, dense_options_error)
    end
    @testset "Constrained" begin
        dense_options = Dict{Symbol, Any}(
            :kkt_system=>MadNLP.DENSE_KKT_SYSTEM,
            :linear_solver=>MadNLPLapackCPU,
        )
        m = 5
        nlp = MadNLPTests.DenseDummyQP(; n=n, m=m)
        ipd = MadNLP.InteriorPointSolver(nlp, option_dict=dense_options)
        ns = length(ipd.ind_ineq)

        kkt = ipd.kkt
        @test isa(kkt, MadNLP.DenseKKTSystem)
        @test size(kkt.jac) == (m, n)
        @test ipd.linear_solver.dense === kkt.aug_com
        @test size(kkt.hess) == (n, n)
        @test length(kkt.pr_diag) == n + ns
        @test length(kkt.du_diag) == m
    end
end

@testset "MadNLP: option kkt_system=$(kkt_system)" for kkt_system in [MadNLP.DENSE_KKT_SYSTEM, MadNLP.DENSE_CONDENSED_KKT_SYSTEM]
    @testset "Size: ($n, $m)" for (n, m) in [(10, 0), (10, 5), (50, 10)]
        _compare_dense_with_sparse(kkt_system, n, m, Int[], Int[])
    end
    @testset "Equality constraints" begin
        n, m = 20, 15
        _compare_dense_with_sparse(kkt_system, n, m, Int[], Int[1, 8]) # test with non-trivial equality constraints
    end
    @testset "Fixed variables" begin
        n, m = 10, 5
        _compare_dense_with_sparse(kkt_system, n, m, Int[1, 2], Int[])
    end
end

@testset "MadNLP: restart (PR #113)" begin
    n, m = 10, 5
    nlp = MadNLPTests.DenseDummyQP(; n=n, m=m)
    sparse_options = Dict{Symbol, Any}(
        :kkt_system=>MadNLP.SPARSE_KKT_SYSTEM,
        :linear_solver=>MadNLPLapackCPU,
        :print_level=>MadNLP.ERROR,
    )
    ips = MadNLP.InteriorPointSolver(nlp, option_dict=sparse_options)
    MadNLP.optimize!(ips)
    # Restart (should hit MadNLP.reinitialize function)
    res = MadNLP.optimize!(ips)
    @test ips.status == MadNLP.SOLVE_SUCCEEDED
end


text:
####
# Rosenbrock's
@counted function rosenbrock(x::Vector)
    return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end

@counted function rosenbrock_gradient(x::Vector)
    storage = zeros(2)
    storage[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
    storage[2] = 200.0 * (x[2] - x[1]^2)
    return storage
end

function rosenbrock_init()
    return clamp.(randn(2), -3.0, 3.0)
end

####
# Powell's
@counted function powell(x::Vector)
    return (x[1] + 10.0 * x[2])^2 + 5.0 * (x[3] - x[4])^2 +
        (x[2] - 2.0 * x[3])^4 + 10.0 * (x[1] - x[4])^4
end

@counted function powell_gradient(x::Vector)
    storage = zeros(4)
    storage[1] = 2.0 * (x[1] + 10.0 * x[2]) + 40.0 * (x[1] - x[4])^3
    storage[2] = 20.0 * (x[1] + 10.0 * x[2]) + 4.0 * (x[2] - 2.0 * x[3])^3
    storage[3] = 10.0 * (x[3] - x[4]) - 8.0 * (x[2] - 2.0 * x[3])^3
    storage[4] = -10.0 * (x[3] - x[4]) - 40.0 * (x[1] - x[4])^3
    return storage
end

function powell_init()
    return clamp.(randn(4), -3.0, 3.0)
end

####
# Himmelblau's
@counted function himmelblau(x::Vector)
    return (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
end

@counted function himmelblau_gradient(x::Vector)
    storage = zeros(2)
    storage[1] = 4.0 * x[1]^3 + 4.0 * x[1] * x[2] -
        44.0 * x[1] + 2.0 * x[1] + 2.0 * x[2]^2 - 14.0
    storage[2] = 2.0 * x[1]^2 + 2.0 * x[2] - 22.0 +
        4.0 * x[1] * x[2] + 4.0 * x[2]^3 - 28.0 * x[2]
    return storage
end

function himmelblau_init()
    return clamp.(randn(2), -3.0, 3.0)
end

const simple_problems = Dict(1=>(f=rosenbrock, g=rosenbrock_gradient, x0=rosenbrock_init, n=20),
                             2=>(f=himmelblau, g=himmelblau_gradient, x0=himmelblau_init, n=40),
                             3=>(f=powell, g=powell_gradient, x0=powell_init, n=100))


text:
module YaoFlux

export Circuit

include("circuit.jl")

end # YaoFlux

max_stars_repo_name: lmondada/YaoFlux

max_issues_repo_name: lmondada/YaoFlux

max_forks_repo_name: lmondada/YaoFlux


text:
PositOrComplex{N,ES} = Union{Complex{Posit{N,ES}}, Posit{N,ES}}

const __real_quire = Quire(Posit{32,2})
const __imag_quire = Quire(Posit{32,2})

#this is a temporary shim to achieve some BLAS calculations

Base.sqrt{N,ES}(x::Posit{N,ES}) = Posit{N,ES}(Float64(x))

#this function is going to be called a lot.

@doc """

  calculates an exact dot product of real posits without quires, but does NOT do
  bounds checking.  This function should be called, but only when that check has
  been performed.

"""
function naked_dot{N,ES}(x::AbstractArray{Posit{N,ES}}, y::AbstractArray{Posit{N,ES}})
    accumulator = zero(Posit{N,ES})
    zero!(__real_quire)
    for idx = 1:length(x)
        @inbounds accumulator = fdp!(__real_quire, x[idx], y[idx])
    end
    accumulator
end

function naked_dot{N,ES}(x::AbstractArray{Complex{Posit{N,ES}}}, y::AbstractArray{Complex{Posit{N,ES}}})
  zero!(__real_quire)
  zero!(__imag_quire)
  real_accumulator = zero(Posit{N,ES})
  imag_accumulator = zero(Posit{N,ES})
  for idx in 1:length(x)
      #note that for complex values, the first vector is conjugated.
      @inbounds fdp!(__real_quire, real(x[idx]), real(y[idx]))
      @inbounds real_accumulator = fdp!(__real_quire, imag(x[idx]), imag(y[idx]))
      @inbounds fdp!(__imag_quire, real(x[idx]), imag(y[idx]))
      @inbounds imag_accumulator = fdp!(__imag_quire, -imag(x[idx]), real(y[idx]))
  end

  real_accumulator + imag_accumulator * im
end

@doc """

  calculates an scaled naked dot product of real posits without quires, but does NOT do
  bounds checking.  This function should be called, but only when that check has
  been performed.  The scale applies to the second vector.

"""
function naked_dot{T<:Posit}(x::AbstractArray{T}, y::AbstractArray{T}, sc::T, seed_1::T = zero(T), seed_2::T = zero(T))
    zero!(__real_quire)
    accumulator = fdp!(__real_quire, seed_1, seed_2)
    for idx = 1:length(x)
        @inbounds accumulator = fdp!(__real_quire, x[idx], y[idx] * sc)
    end
    accumulator
end

@generated function naked_dot{T<:Posit, C}(
    x::AbstractArray{Complex{T}},
    y::AbstractArray{Complex{T}},
    sc::T,
    seed_1::Complex{T} = zero(Complex{T}),
    seed_2::Complex{T} = zero(Complex{T}),
    ::Type{Val{C}} = Val{:noconj})

  if C == :conj
    main_code = quote
      for idx in 1:length(x)
        #note that for complex values, the first vector is conjugated.
        @inbounds                    fdp!(__real_quire,  real(x[idx]), real(y[idx]) * sc)
        @inbounds real_accumulator = fdp!(__real_quire,  imag(x[idx]), imag(y[idx]) * sc)
        @inbounds                    fdp!(__imag_quire,  real(x[idx]), imag(y[idx]) * sc)
        @inbounds imag_accumulator = fdp!(__imag_quire, -imag(x[idx]), real(y[idx]) * sc)
      end
    end
  elseif C == :noconj
    main_code = quote
      for idx in 1:length(x)
        @inbounds                    fdp!(__real_quire,  real(x[idx]), real(y[idx]) * sc)
        @inbounds real_accumulator = fdp!(__real_quire, -imag(x[idx]), imag(y[idx]) * sc)
        @inbounds                    fdp!(__imag_quire,  real(x[idx]), imag(y[idx]) * sc)
        @inbounds imag_accumulator = fdp!(__imag_quire,  imag(x[idx]), real(y[idx]) * sc)
      end
    end
  end

  quote
    real_accumulator = zero(T)
    imag_accumulator = zero(T)
    zero!(__real_quire)
    zero!(__imag_quire)

    fdp!(__real_quire, real(seed_1), sc)
    fdp!(__imag_quire, imag(seed_1), sc)

    $main_code

    real_accumulator + imag_accumulator * im
  end
end

function naked_dot{T<:Posit}(x::AbstractArray{Complex{T}}, y::AbstractArray{Complex{T}}, sc::Complex{T})
  zero!(__real_quire)
  zero!(__imag_quire)
  real_accumulator = zero(T)
  imag_accumulator = zero(T)
  for idx in 1:length(x)
      @inbounds                    fdp!(__real_quire,  real(x[idx]), real(y[idx]) * real(sc))
      @inbounds                    fdp!(__real_quire,  imag(x[idx]), imag(y[idx]) * real(sc))
      @inbounds                    fdp!(__real_quire,  imag(x[idx]), real(y[idx]) * imag(sc))
      @inbounds real_accumulator = fdp!(__real_quire,  real(x[idx]), imag(y[idx]) * imag(sc))

      @inbounds                    fdp!(__imag_quire,  real(x[idx]), real(y[idx]) * imag(sc))
      @inbounds                    fdp!(__imag_quire,  real(x[idx]), imag(y[idx]) * real(sc))
      @inbounds                    fdp!(__imag_quire,  imag(x[idx]), real(y[idx]) * real(sc))
      @inbounds imag_accumulator = fdp!(__imag_quire,  imag(x[idx]), imag(y[idx]) * imag(sc))
  end

  real_accumulator + imag_accumulator * im
end

max_stars_repo_name: ityonemo/SigmoidNumbers

max_issues_repo_name: ityonemo/SigmoidNumbers

max_forks_repo_name: ityonemo/SigmoidNumbers


text:
## Arithmetics

Base.:+(x::DiffCheck{T, S, N}, y::U) where {T<:ReComp, S, N, U<:ReComp} = DiffCheck{T, S}(value(x) + T(y), derivatives(x))
Base.:+(y::U, x::DiffCheck{T, S, N}) where {T<:ReComp, S, N, U<:ReComp} = Base.:+(x, y)

@generated function Base.:+(x::DiffCheck{T₁, S, N}, y::DiffCheck{T₂, Q, M}) where {T₁, T₂<:ReComp, S, Q, N, M}
   Z = Tuple(union(Set(S), Set(Q)))
   N = length(Z)
   str = "("
   for z in Z
        if (z in S && z in Q) # Since Z is the union both cannot be false
           str*="x.der[:$z]+y.der[:$z], "
        elseif z in S
           str*="x.der[:$z],"
        else
           str*="y.der[:$z],"
        end
    end
    T = promote_type(T₁,T₂)
    DATA = StaticTuple{T, Z, N}
    DC = DiffCheck{T, Z}
    str*= ")"
    expr = Meta.parse(str)
    return :($DC(value(x) + value(y),$DATA($expr)))
end

Base.:-(x::DiffCheck{T, S, N}, y::U) where {T<:ReComp, S, N, U<:ReComp} = DiffCheck{T, S}(value(x) - T(y), derivatives(x))
Base.:-(y::U, x::DiffCheck{T, S, N}) where {T<:ReComp, S, N, U<:ReComp} =
    DiffCheck{T, S}(T(y) - value(x), StaticTuple{T, S, N}(-1 .* values(derivatives(x))))

@generated function Base.:-(x::DiffCheck{T₁, S, N}, y::DiffCheck{T₂, Q, M}) where {T₁, T₂<:ReComp, S, Q, N, M}
   Z = Tuple(union(Set(S), Set(Q)))
   N = length(Z)
   str = "("
   for z in Z
        if (z in S && z in Q) # Since Z is the union both cannot be false
           str*="x.der[:$z]-y.der[:$z], "
        elseif z in S
           str*="x.der[:$z],"
        else
           str*="-y.der[:$z],"
        end
    end
    str*= ")"
    expr = Meta.parse(str)
    return :(DiffCheck{promote_type($T₁,$T₂), $Z}(value(x) - value(y), StaticTuple{promote_type($T₁,$T₂),$Z, $N}($expr)))
end

Base.:*(x::DiffCheck{T, S, N}, y::U) where {T<:ReComp, S, N, U<:ReComp} =
    DiffCheck{T, S}(y*value(x), StaticTuple{T, S, N}(y .* values(derivatives(x))))

Base.:*(y::U, x::DiffCheck{T, S, N}) where {T<:ReComp, S, N, U<:ReComp} = x*y

@generated function Base.:*(x::DiffCheck{T₁, S, N}, y::DiffCheck{T₂, Q, M}) where {T₁, T₂<:ReComp, S, Q, N, M}
   Z = Tuple(union(Set(S), Set(Q)))
   N = length(Z)
   str = "("
   for z in Z
        if (z in S && z in Q) # Since Z is the union both cannot be false
           str*="x.value*y.der[:$z]+x.der[:$z]*y.value, "
        elseif z in S
           str*="y.value*x.der[:$z],"
        else
           str*="x.value*y.der[:$z],"
        end
    end
    str*= ")"
    expr = Meta.parse(str)
    return :(DiffCheck{promote_type($T₁,$T₂), $Z}(value(x) * value(y), StaticTuple{promote_type($T₁,$T₂),$Z, $N}($expr)))
end

Base.:/(x::DiffCheck{T, S, N}, y::U) where {T<:ReComp, S, N, U<:ReComp} =
    DiffCheck{T, S}(value(x)/y, StaticTuple{T, S, N}(values(derivatives(x))./ y))

Base.:/(y::U, x::DiffCheck{T, S, N}) where {T<:ReComp, S, N, U<:ReComp} = DiffCheck{T, S}(y/value(x), StaticTuple{T, S, N}( -y/value(x) .* (values(derivatives(x)) ./ value(x))))

@generated function Base.:/(x::DiffCheck{T₁, S, N}, y::DiffCheck{T₂, Q, M}) where {T₁, T₂<:ReComp, S, Q, N, M}
   Z = Tuple(union(Set(S), Set(Q)))
   N = length(Z)
   str = "("
   for z in Z
        if (z in S && z in Q) # Since Z is the union both cannot be false
           str*="-1*(x.value/y.value)*y.der[:$z]/y.value+x.der[:$z]*(x.value/y.value), "
        elseif z in S
           str*="x.der[:$z]/y.value,"
        else
           str*="-1*(x.value/y.value)*y.der[:$z]/y.value,"
        end
    end
    str*= ")"
    expr = Meta.parse(str)
    return :(DiffCheck{promote_type($T₁,$T₂), $Z}(value(x) / value(y), StaticTuple{promote_type($T₁,$T₂),$Z, $N}($expr)))
end


text:
using Statistics, LineSearches, Optim

function buildRecord(stack::stFlux, nEpochs)
    p = params(stack).order
    return [zeros(Float32, size(x)..., nEpochs + 1) for x in p]
end

function buildRecord(stack, nEpochs)
    zeros(eltype(stack), size(stack)..., nEpochs)
end

function expandRecord!(record::Array{Array}, nEpochs)
    for ii in 1:length(record)
        adding = zeros(Float32, size(record[ii])[1:end - 1]..., nEpochs)
        record[ii] = cat(record[ii], adding, dims=ndims(record[ii]))
        size(record[ii])
    end
end

function expandRecord!(record::Array{<:Number}, nEpochs)
    adding = zeros(eltype(record), size(record)[1:end - 1]..., nEpochs)
    record[:] = cat(record, adding, dims=ndims(record))
end

function addCurrent!(record, stack::stFlux, ii)
    p = params(stack).order

    for (kk, x) in enumerate(p)
        ax = axes(x)
        record[kk][ax..., ii + 1] = x
    end
end

function addCurrent!(record, result, ii)
    p = params(stack).order
    result = 3
    record[]
    for (kk, x) in enumerate(p)
        ax = axes(x)
        record[kk][ax..., ii + 1] = x
    end
end

"""
    obj = makeObjFun(target, path::pathLocs, st, normalize=st.normalize,
                     λ=1e-10, jointObj=false)

Make an objective function `obj(x)` that evaluates the difference between a
target stOutput and the input `x`.

    makeObjFun(target::AbstractArray, st, jointObj=false, kwargs...)
jointObj makes a function `obj(x,y)` where `x` is as otherwise and `y` has the
shape of target.

Note that the joint case is defined for y as an array. In that case, the path is
ignored and the difference across the whole output is evaluted.
"""
function makeObjFun(target::ScatteredOut, path::pathLocs, st,
                    normalize=st.normalize, λ=1e-10, jointObj=false, innerProd=false)
    if jointObj # return a function of both x and y?
        if normalize # normalize the input
            if innerProd # just maximize the dot product?
                jointNormIPObjxy(x, y::ScatteredOut) = sum(mse(a[1], a[2]) for a in zip(st(x)[path], y[path])) + λ * norm(x)
                function jointNormIPObjxy(x, y::AbstractArray)
                    -sum(y' * ScatteringTransform.flatten(st(x))) + λ * norm(x)
                end
                return jointNormIPObjxy
            else
                jointNormObjxy(x, y::ScatteredOut) = sum(mse(a[1], a[2]) for a in zip(st(x)[path], y[path])) + λ * norm(x)
                function jointNormObjxy(x, y::AbstractArray)
                    sum((ScatteringTransform.flatten(st(x)) .- y).^2) + λ * norm(x)
                end
                return jointNormObjxy
            end
        else
            jointObjxy(x, y::ScatteredOut) = sum(mse(a[1], a[2]) for a in
                                                 zip(st(x)[path], y[path]))
            function jointObjxy(x, y::AbstractArray)
                sum((ScatteringTransform.flatten(st(x)) .- y).^2)
            end
            return jointObjxy
        end
    else
        if normalize
            return normObjx(x) = sum(mse(a[1], a[2]) for a in zip(st(x)[path], target[path])) + λ * norm(x)
        else
            return unnormObjx(x) = sum(mse(a[1], a[2]) for a in zip(st(x)[path], target[path]))
        end
    end
end

"""
    obj(x) = makeNormMatchObjFun(target, st, normTarget, λ)
fit the target in the scattering domain and try to keep the input domain norm near `normTarget`.
"""
function makeNormMatchObjFun(target, st, normTarget, λ)
    jointNormObj(x) = sum(mse(a[1], a[2]) for a in zip(st(x), target)) + λ * (normTarget - norm(x))^2
    function jointNormObjxy(x)
        sum((ScatteringTransform.flatten(st(x)) .- target).^2) + λ * (normTarget - norm(x))^2
    end
end
"""
    makeObjFun(target, path::pathLocs, st, normalize=st.normalize)
if we want to fit the whole path for the target
"""
function makeObjFun(target::ScatteredOut, st, kwargs...)
    path = pathLocs(0, :, 1, :, 2, :)
    makeObjFun(target, path, st, kwargs...)
end

function makeObjFun(target::AbstractArray, st, kwargs...)
    makeObjFun(roll(target, st), st, kwargs...)
end

function fitReverseSt(N, initGuess; opt=Momentum(), ifGpu=identity, obj=nothing, keepWithin=-1, stochastic=false, errTol=1e5, timeLimit=Year(10), adjust=false, kwargs...)
    ongoing = copy(initGuess |> ifGpu);
    pathOfDescent = zeros(Float64, size(ongoing, 1), 1, size(initGuess, 3), N + 1) |> ifGpu;
    println("initial objective function is $(obj(ongoing))")
    pathOfDescent[:,:,:,1] = ongoing |> ifGpu;
    err = zeros(Float64, N + 1); err[1] = obj(pathOfDescent[:,:,:,1])
    return justTrain(N, pathOfDescent, err, 1, ongoing; opt=opt, ifGpu=ifGpu,
                      obj=obj, keepWithin=keepWithin, stochastic=stochastic,
                      errTol=errTol, timeLimit=timeLimit, adjust=adjust, kwargs...)
end

function adjustStepSizeBasedOnVariance!(err, opt, ii)
    aveErrChange = (mean(err[ii - 22:ii - 18]) - mean(err[ii - 3:ii])) / err[ii]
    # if the err is flat (less than a .1% change over the past 20 steps), change the step size
    if (ii > 30) && (ii % 30 == 0) &&  (aveErrChange < 1e-6)
        # if it varies too much, or has been overall increasing, decrease the step, otherwise, increase
        if (count(abs.(err[ii - 20:ii] .- mean(err[ii - 20:ii])) ./ err[ii] .> .0003) > 5 ||
            aveErrChange < 0 )
            # println("$(mean(err[ii - 22:ii - 18]))  $(mean(err[ii - 3:ii]))")
            opt.eta = opt.eta * .9
        end
    end
end

function justTrain(N, pathOfDescent, err, prevLen, ongoing; opt=Momentum(),
                   ifGpu=identity, obj=obj, keepWithin=-1, stochastic=false,
                   errTol=1e5, timeLimit=Year(10), adjust=false, record=true,
                   relErrEnd=1e-7,relGradNorm=0)
    endTime = now() + timeLimit
    for ii in (prevLen + 1):(prevLen + N)
        err[ii] = obj(ongoing)
        if err[ii] > errTol || isnan(err[ii])
            println("way too big at $(ii), aborting")
            break
        elseif err[ii] / err[1] < relErrEnd
            println("oh, we actually finished?")
            err = err[1:ii]
            pathOfDescent = pathOfDescent[:,:,:,1:(ii - 1)]
            return (pathOfDescent, err, opt)
        end
        if ii % 10 == 0
            println("$ii rel err: $(err[ii] / err[1]), abs err: $(err[ii]), η=$(opt.eta)")
        end
        ∇ = gradient(obj, ongoing)
        if norm(∇[1]) <= relGradNorm
            println("oh, we actually finished?")
            err = err[1:ii]
            pathOfDescent = pathOfDescent[:,:,:,1:(ii - 1)]
            return (pathOfDescent, err, opt)
        end

        if adjust && ii > 30
            adjustStepSizeBasedOnVariance!(err, opt, ii)
        end

        update!(opt, ongoing, ∇[1])
        sz = norm(ongoing)
        if keepWithin > 0 && sz > keepWithin
            ongoing = ongoing / norm(ongoing) * keepWithin
        end

        if stochastic
            ongoing += opt.eta * .01 * norm(∇, 1) * randn(size(ongoing)...)
        end

        if record
            pathOfDescent[:,:,:,ii] = ongoing
        else
            pathOfDescent[:,:,:,end] = ongoing
        end

        if endTime ≤ now()
            println("out of time")
            err = err[1:end - (N - ii + 1)]
            pathOfDescent = pathOfDescent[:,:,:,1:end - (N - ii + 1)]
            return (pathOfDescent, err, opt)
        end
    end
    return (pathOfDescent, err, opt)
end

function continueTrain(N, pathOfDescent, err; opt=Momentum, ifGpu=identity,
                       obj=obj, keepWithin=-1,stochastic=false, errTol=1e5,
                       timeLimit=Year(10), adjust=false, record=true,kwargs...)
    prevLen = size(pathOfDescent, 4)
    currentGuess = pathOfDescent[:,:,:,end] |> ifGpu;
    if record
        pathOfDescent = cat(pathOfDescent,
                            zeros(Float64, size(pathOfDescent, 1), 1, size(pathOfDescent, 3), N),
                            dims=4) |> ifGpu;
    else
        pathOfDescent = cat(pathOfDescent, zeros(Float64, size(pathOfDescent,
                                                               1), 1, size(pathOfDescent, 3), 1), dims=4)
    end
    err = cat(err, zeros(Float64, N), dims=1)
    return justTrain(N, pathOfDescent, err, prevLen, currentGuess; opt=opt, ifGpu=ifGpu,
                      obj=obj, keepWithin=keepWithin, stochastic=stochastic,
                     errTol=errTol, timeLimit=timeLimit, adjust=adjust,kwargs...)
end

# Everything  after this point is built on using Optim.jl frameworks

function maximizeSingleCoordinate(N, initGuess, p,target, st; ifGpu=identity,
                                  obj=nothing, keepWithin=-1,
                                  stochastic=false, spacing=:linear,
                                  adjustStepSize=false, allowedTime=-1,
                                  λ=1e-10, tryCatch=true,kwargs...)
    ongoing = copy(initGuess) |> ifGpu;
    if obj != nothing
        println("initial objective function is $(obj(ongoing))")
    end
    if allowedTime > 0
        totalUsed = 0
        N = 500
    else
        totalUsed = NaN
    end
    paths = Array{Optim.MultivariateOptimizationResults,1}(undef, 20)
    pathErr = Array{Array{Float64,1},1}(undef, 20)
    m = findfirst(p.indices .!= nothing) - 1 # which layer are we working in?

    if spacing == :exponential
        targetVals = exp.(range(0, log(prod(size(target[m])[2:end - 1])), length=20))
    elseif spacing == :linear
        targetVals = range(1, prod(size(target[m])[2:end - 1]), length=20)
    else
        error("no spacing type $(spacing)")
    end

    lastI = -1
    for (ii, t) in enumerate(targetVals)
        println("-----------------------------------------------------------------------")
        println("target value $t")
        println("-----------------------------------------------------------------------")
        target[p] = t
        divvyAllowedTime = (allowedTime - totalUsed) / (20 - ii + 1)
        println("allowed time = $(divvyAllowedTime)")
        tmp, λtmp, err, netTime = fitByPerturbing(N, ongoing, p, target, st;
                                    λ=λ, timeAl=divvyAllowedTime,
                                    tryCatch=tryCatch, kwargs...)
        pathErr[ii] = err
        if tmp != nothing
            paths[ii] = tmp; λ = λtmp
        else
            println("this broke, doubtful it will converge after")
            break
        end
        totalUsed += netTime
        ongoing = paths[ii].minimizer
        lastI = ii
        if totalUsed > allowedTime
            break
        end
    end
    return paths[1:lastI], pathErr[1:lastI]
end

"""
Given a collection of vectors β used in multinomial regression, fit the positive
and negative parts separately
"""
function fitPositiveNegativeSeperately(β, inputSize; Niter=100, nExEach=1,
                                       initGen=randn, trainThese=1:size(β, 2),
                                       varargs...)
    βplus = max.(β, 0)
    βminus = max.(-β, 0)
    initPlus = initGen(inputSize, 1, nExEach, length(trainThese))
    initMinus = initGen(inputSize, 1, nExEach, length(trainThese))
    for i = trainThese
        plusObj = makeObjFun(βplus)
    end

end

"""
    perturbedRes, λ, errGrouped, totalUsed = fitByPerturbing(N, initGuess,
                         p,target, st; obj=nothing, allowedTime = -1,
                         rate = .1, rateFraction=.9, missesInARow=5,
                         NSamples=1000, ifGpu= identity, noiseType=:pink,
                         NAttempts=1e6, varargs...)

implementing the perturbed BFGS idea. Differences: perturbing using pink noise
whose norm is dependent on the size of the signal, whereas theirs is a uniform.
+ `N` means the number of iterations
+ `timeAl` is the amount of time allowed (in seconds)
+ `rate` is the step distance we take during a perturbation. Decreases as the
steps go too far
+ `missesInARow`: the number of misses before adjusting `rate`
"""
function fitByPerturbing(N, initGuess, p,target, st;  timeAl=-1,
                         rate=.1f0, rateFraction=.8f0, λ=1, missesInARow=1,
                         obj=makeObjFun(target, p, st, true, λ),
                         NSamples=500, ifGpu=identity, noiseType=:pink,
                         NAttempts=1e3, varargs...)

    if timeAl > 0
        totalUsed = 0
        N = 10000000000000
    else
        totalUsed = NaN
    end
    pert = copy(initGuess) |> ifGpu
    tooMany = 0
    errGrouped = zeros(0)
    # the first algorithm is fixed, so don't pass it through varargs
    fixedAlgo = filter(x -> x[1] != :algo, varargs)
    (perturbedRes, λ) = fitUsingOptim(N, pert, p, target, st;
                                    timeAl=timeAl - totalUsed, λ=λ,
                                    ifGpu=ifGpu, algo=BFGS, obj=obj, fixedAlgo...)
    fprev = perturbedRes.minimum # get the initial fit value
    errThisRun = [x.value for x in perturbedRes.trace]
    append!(errGrouped, errThisRun)
    totalUsed += perturbedRes.time_run
    for j = 1:NAttempts
        println("perturbing for the $(j)th time. rate = $(rate)")
        # generate the next example
        (pert, minObj, tooMany), timeUsed, _ = @timed perturb(perturbedRes.minimizer, obj, rate, NSamples, noiseType, tooMany)
        println(minObj)
        totalUsed += timeUsed
        oldRes = perturbedRes
        perturbedRes, λ = fitUsingOptim(N, pert, p, target, st;
                                        timeAl=timeAl - totalUsed, obj=obj,
                                        ifGpu=ifGpu, λ=λ, varargs...)
        fcur = perturbedRes.minimum
        # have we stopped progressing during our perturbations?
        if abs(fprev - fcur) / abs(fcur) < eps(typeof(fcur))
            println("there was effectively no change in the objective function, aborting")
            break
        end
        # add condition for if the norm of the function is the source of
        # most of the error, as this is only a fitting concern
        println("----------------------------------------------")
        trueErr = abs(λ * norm(perturbedRes.minimizer) - fcur)
        println("The true error, without regularization: ", trueErr)
        println("----------------------------------------------")
        if trueErr < 1e-1
            println("most of the error is from the norm, breaking")
            break
        end
        fprev = fcur

        # increase the time used
        totalUsed += perturbedRes.time_run

        # record the error
        errThisRun = [x.value for x in perturbedRes.trace]
        append!(errGrouped, errThisRun)

        #  have we had tooMany missesInARow?, then decrease the distance
        #  we're stepping away
        if tooMany >= missesInARow
            rate *= rateFraction
            tooMany = 0
        end

        # are we perturbing by too small of a step?
        if rate < 1e-4
            println("had to decrease the rate too many times")
            break
        end
        # are we over time?
        if totalUsed > timeAl
            println("Out of time")
            break
        end
        # are we sufficiently converged?
        if perturbedRes.minimum < 1e-7
            break
        end
        if 1.01 * oldRes.minimum < perturbedRes.minimum
            println("the new one is more than 10% worse, starting again from the old one")
            perturbedRes = oldRes
        end
    end
    return perturbedRes, λ, errGrouped, totalUsed
end

function perturb(x, obj, stepSize=.01f0, K=100, noiseType=:pink, tooMany=0)
    tmp = genNoise(size(x), K, noiseType) # get your favorite kind of noise
    tmp = stepSize * norm(x) * tmp ./
        reshape([norm(w) for w in eachslice(tmp, dims=4)], (1, 1, 1, K)) # normalize
    # it to match stepSize times the norm of x
    perturbed = x .+ tmp
    errs = [obj(w) for w in eachslice(perturbed, dims=4)]
    μ = argmin(errs)            # find the best
    if errs[μ] > obj(x)
        println("couldn't find a better location. Adjust stepSize or number of samples K")
        tooMany += 1
    else
        tooMany = max(0, tooMany - 1)
    end
    return perturbed[:,:,:,μ], errs[μ], tooMany
end

function genNoise(sz, K, noiseType)
    if noiseType == :white
        tmp = randn(Float32, sz..., K)
        return tmp / norm(tmp)
    elseif noiseType == :pink
        N = sz[1]
        pink = 1 ./ Float32.(1 .+ (0:N >> 1)) # note: not 1:N>>1+1
        phase = exp.(Float32(2π) * im) .* rand(Float32, N >> 1 + 1, sz[2:end]..., K)
        res = irfft(pink .* phase, N, 1)
        return res ./ norm(res)
end
end

"""
given a list of optim solutions, choose the one with the largest entry at
location p
"""
function chooseLargest(p, res, st)
    ii = argmax([st(x.minimizer)[p][1] for x in res])
    return res[ii]
end

function makeOptimGradient(obj)
    function obj∇!(F, ∇, x)
        y, back = pullback(obj, x)
        if ∇ != nothing
            ∇[:] = back(y)[1]
        end
        if F != nothing
            return y
        end
    end
    return obj∇!
end
"""
    fitUsingOptim(N, ongoing, p, target, st; ifGpu=identity, λ=1e-10,
                        timeAl=NaN, tryCatch=true, allowIncrease=false,
                        lineSearch=HagerZhang(), algo=BFGS,
                        obj=makeObjFun(target, p, st, true, λ), kwargs...)

fit the arbitrary one variable objective function obj in a way that
adjusts the normalization penalty if there's an `ArgumentError` and the type of
linesearch if the isn't finite or if B > A.
"""
function fitUsingOptim(N, ongoing, p, target, st; ifGpu=identity, λ=1e-10,
                        timeAl=NaN, tryCatch=true, allowIncrease=false,
                        lineSearch=HagerZhang(), algo=BFGS,
                        obj=makeObjFun(target, p, st, true, λ), kwargs...)
    obj∇! = makeOptimGradient(obj)
    if tryCatch
        try
            res = defaultOptim(obj∇!, ongoing, N, timeAl, allowIncrease, algo,
                               lineSearch; kwargs...)
        catch e
            if isa(e, AssertionError) &&  (occursin("isfinite", e.msg) || occursin("B > A", e.msg))
                res, lineSearch, λ = catchAndSwitchLineSearch(e, obj∇!, ongoing,
                       N, target, p, st, λ, timeAl, allowIncrease,algo,
                       lineSearch; kwargs... )
            elseif isa(e, ArgumentError)
                # too little weight placed on the normalization
                res, lineSearch, λ = catchAndUpλ(e, ongoing, N, target, p, st, λ,
                       timeAl, allowIncrease,algo, lineSearch; kwargs...)
            else
                rethrow(e)
            end
        end
    else
        res = defaultOptim(obj∇!, ongoing, N, timeAl, allowIncrease,
                           algo, lineSearch; kwargs...)
    end
    return res, λ
end
function defaultOptim(obj∇!,ongoing,N,timeAl, allowIncrease=false, algo=BFGS,
                      lineSearch=HagerZhang(); show_every=20, kwargs...)
    constAlgo = algo(;linesearch=lineSearch, kwargs...)
    optimOption = Optim.Options(store_trace=true, show_trace=true, iterations=N,
                           show_every=show_every, time_limit=timeAl, g_tol=1e-8,
                           f_tol=0.0, allow_f_increases=allowIncrease)
    optimize(Optim.only_fg!(obj∇!), ongoing, constAlgo, optimOption)
end

function catchAndSwitchLineSearch(e, obj∇!, ongoing, N, target, p, st, λ, timeAl, allowIncrease, algo, lineSearch; kwargs...)
    if typeof(lineSearch) <: HagerZhang
        lineSearch = BackTracking()
    elseif typeof(lineSearch) <: BackTracking
        lineSearch = MoreThuente()
    else
        println("well, we tried HagerZhang, BackTracking and MoreThente, and they all broke")
        return nothing, nothing, nothing
        # rethrow(e)
    end
    println("the current line search is $(lineSearch)")
    try
        res = defaultOptim(obj∇!, ongoing, N, timeAl, allowIncrease, algo,
                           lineSearch; kwargs...)
        return res, lineSearch, λ
    catch e2
        if isa(e, AssertionError) && (occursin("isfinite", e.msg) || occursin("B > A", e.msg))
            res, lineSearch, λ = catchAndSwitchLineSearch(e2, obj∇!, ongoing, N,
                                                       target, p, st, λ,
                                                       timeAl, allowIncrease,
                                                       algo, lineSearch;
                                                       kwargs...)
            return res, lineSearch, λ
        elseif isa(e, ArgumentError) && occursin("Value and slope", e.msg)
            res, λ = catchAndUpλ(e, ongoing, N, target, p, st, λ, timeAl, allowIncrease, algo, lineSearch; kwargs...)
        else
            rethrow(e2)
        end
    end
end

function catchAndUpλ(e, ongoing, N, target, p, st, λ, timeAl, allowIncrease, algo, lineSearch; kwargs...)
    println("too little weight placed on the normalization, upping to $(λ * 1e2)")
    while λ < 1e8
        λ = λ * 1e2
        obj = makeObjFun(target, p, st, true, λ)
        function obj∇!(F, ∇, x)
            y, back = pullback(obj, x)
            if ∇ != nothing
                ∇[:] = back(y)[1]
            end
            if F != nothing
                return y
            end
        end
        try
            res = defaultOptim(obj∇!,ongoing,N,timeAl, allowIncrease, algo,
                               lineSearch; kwargs...)
            return res, lineSearch, λ
        catch e2
            if isa(e2, AssertionError)
                res, lineSearch, λ = catchAndSwitchLineSearch(e2, obj∇!,
                                                              ongoing, N,
                                                              target, p, st, λ,
                                                              timeAl,
                                                              allowIncrease,
                                                              algo,
                                                              lineSearch;
                                                              kwargs...)
                return res, lineSearch, λ
            elseif isa(e2, ArgumentError)
                if λ < 1e3
                    println("still broken, upping to $(λ * 1e2)")
                else
                    println("well, it seems even $λ isn't large enough.")
                    return nothing, nothing, nothing
                end
            else
                rethrow(e2)
            end
        end
    end
end

distFromI(ii::CartesianIndex, jj, normType) = norm(ii.I .- jj, normType)
distFromI(ii, jj, normType) = norm(ii - jj, normType)
distFromI(ii, jj::Tuple{<:Integer}, normType) = norm(ii - jj[1], normType)
function computeSingleWeightMatrix(target, p::pathLocs, Pnorm, Nd=1)
    outDims = map(x -> x[1 + Nd:end], size(target))
    weights = map(x -> fill(Inf, x...), outDims)
    n = ndims(target)
    for (layerI, jj) in enumerate(p.indices)
        if jj !== nothing
            layInds = eachindex(view(weights[layerI], axes(weights[layerI])...))
            weights[layerI][:] = map(ii -> distFromI(ii, jj[1 + n:end - 1], Pnorm), layInds)
        end
    end
    return weights
end
function computeWeightMatrices(target, pathCollection, distanceNorm)
    allWeights = map(p -> computeSingleWeightMatrix(target, p, distanceNorm), pathCollection)
    # find the minimal distance to each location
    netWeight = map(x -> zeros(size(x)), allWeights[1]) # set weights to zero initally
    for layI in 1:length(netWeight) # iterate over layers
        nDimsWeight = ndims(allWeights[1][layI]) # the ndims of this layer
        joined = cat([allWeights[pathInd][layI] for pathInd = 1:length(allWeights)]..., dims=nDimsWeight + 1)
        shortestDistance = minimum(joined, dims=nDimsWeight + 1) # closest pathLoc in this layer
        netWeight[layI][:] = shortestDistance
    end
    return netWeight
end
function weightByDistance(target, pathCollection, St, distanceNorm=1.0, elementNorm=2.0)
    netWeight = computeWeightMatrices(target, pathCollection, distanceNorm)
    to = target.output
    usedLayers = [i - 1 for i = 1:length(to) if netWeight[i][1] != Inf]
    function objFun(x)
        Stx = St(x)
        total = 0.0
        # minimize the irrelevant locations
        for ii in usedLayers
            w = netWeight[ii + 1] # just a matrix
            total += sum(reshape(w, (1, size(w)...)) .* (Stx[ii] .- target[ii]).^elementNorm)
        end
        # maximize the relevant ones
        onPath = mapreduce(p -> sum((Stx[p] .- target[p]).^elementNorm), +,  pathCollection)
        return total + onPath
    end
end

"""
    adjustPopulation!(setProbW, popMat)
given a BlackBoxOptim `OptController` with some kind of population based method, set the population to be the matrix given in `popMat`; each column of `popMat` is a different individual.
"""
function adjustPopulation!(setProbW, popMat)
    for i = 1:size(popMat, 2)
        setProbW.optimizer.population[i] = copy(popMat[:,i])
    end
end

"""
    bandpassed,objValue, freqInd = BandpassStayInFourier(x,FourierObj)
given a dct representation of a signal `x` and an objective function `FourierObj` which operates on the dct representation of a signal, find the frequency index a super naive way of finding a bandpass threshold so that x (which is represented as dct of the target)  FourierObj
"""
function BandpassStayInFourier(x, FourierObj)
    flattenFrom(x, i) = cat(x[1:i,axes(x)[2:end]...], zeros(size(x, 1) - i), dims=1)
    evalFilters = [FourierObj(flattenFrom(x, ncut)) for ncut = 2:size(x, 1)]
    cutFreq = argmin(evalFilters)
    return flattenFrom(x, cutFreq + 1), evalFilters[cutFreq], cutFreq + 1
end

"""
if sortIndArray isn't defined, it comes from calling sortperm on the list of
scores, e.g.
    scores = [FourierObj(x) for x in eachslice(spacePop, dims=2)]
    sortIndArray = sortperm(scores)
"""
function dampenAboveBestHardbandpass(freqPop, sortIndArray, FourierObj, dampenAmount=.1)
    flattened, val, freqI = BandpassStayInFourier(freqPop[:,sortIndArray[1]], FourierObj)
    newFreqPop = copy(freqPop)
    newFreqPop[freqI + 1:end,:] .*= dampenAmount
    return newFreqPop
end

function saveSerial(name, object...)
    fh = open(name, "w")
    serialize(fh, object...)
    close(fh)
end


text:
module Modules
export AST, Cond
include("AST.jl")
include("Cond.jl")
end


text:
import ACE
using LinearAlgebra, StaticArrays, BenchmarkTools, Test, Printf
using ACE.SphericalHarmonics
using ACE.SphericalHarmonics: dspher_to_dcart, SphericalCoords,
               cart2spher, spher2cart, ALPolynomials
using ACE: evaluate, evaluate_d, evaluate_ed, alloc_B, release!,
           evaluate!, acquire!, acquire_B!, acquire_dB!
using ACE.Testing
using BenchmarkTools

##

L = 10

alp = ALPolynomials(L, Float64)

S = SphericalCoords(0.0, 0.0, rand() * π)

const pl = ACE.ObjectPools.ArrayPool()
const hpl = ACE.ObjectPools.HomogeneousVectorPool{Float64}()
const stpl = ACE.ObjectPools.StaticVectorPool{Float64}()

function _ev(b, x)
   B = acquire!(pl, Float64, (length(b),))
   evaluate!(B, b, x)
   release!(pl, B)
end

function _hev(b, x)
   B = acquire!(hpl, length(b))
   evaluate!(B, b, x)
   release!(hpl, B)
end

function _stev(b, x)
   B = acquire!(stpl, length(b))
   evaluate!(B, b, x)
   release!(stpl, B)
end

B = ACE.acquire_B!(alp)

##

@info("ALP")
@info("non-allocating (for reference)")
@btime evaluate!($B, $alp, $S)
@info("ArrayPool")
@btime _ev($alp, $S)
@info("HomogeneousVectorPool")
@btime _hev($alp, $S)
@info("StaticVectorPool")
@btime _stev($alp, $S)

##

@info("Spherical Harmonics")
sh = SHBasis(L)
x = randn(SVector{3, Float64})
B_sh = acquire_B!(sh)
P_sh = acquire_B!(sh.alp)

@info("non-allocating")
@btime ACE.SphericalHarmonics.__evaluate!($B_sh, $sh, $P_sh, $x)
@info("ALP allocates via pool, Bsh is passed in")
@btime evaluate!($B_sh, $sh, $x)

##

@info("Spherical Harmonics gradients")
dB_sh = ACE.SphericalHarmonics.acquire_dB!(sh, x)
dP_sh = acquire_dB!(sh.alp)

@info("non-allocating")
@btime ACE.SphericalHarmonics.__evaluate_ed!($B_sh, $dB_sh, $sh, $P_sh, $dP_sh, $x)
@info("ALP allocates via pool, Bsh is passed in")
@btime ACE.evaluate_ed!($B_sh, $dB_sh, $sh, $x)


text:
using Statistics
using ProgressMeter
using MATLAB
include("atn_module.jl")

N = 10_000          # number of nodes
m = 5               # number of active links per active node
mu = 0.1            # probability I -> S
tmax = 2000         # maximum duration of the simulation
Ni = Int(0.01 * N)  # initial number of infects
gamma = -2.1      # exponent of the power law
NTrials = 1       # number of trials for each value of lambda
TWindow = 500.0      # time window for evaluating steady state
TMinConvergence = 500

LambdaMin = 0.0
LambdaMax = 1.0
LambdaStep = 0.05
Lambda = LambdaMin:LambdaStep:LambdaMax
LambdaSize = length(Lambda)

EtaMin = 0.0
EtaMax = 5
EtaStep = 0.1
Eta = EtaMin:EtaStep:EtaMax
Eta = 5
EtaSize = length(Eta)
AvgInfectedLambda = zeros(LambdaSize)
Ni = Int(0.01 * N)  # initial number of infects
gamma = -2.1       # exponent of the power law
TWindow = 500       # time window for evaluating steady state
TMinConvergence = 500
ITrend = zeros(Int, 0)
ActMin = 0.001      # minimum activation probability
ActMax = 1
AgentAct =
    (
        (ActMax.^(gamma + 1) - ActMin.^(gamma + 1)) .* rand(N) .+
        ActMin.^(gamma + 1)
    ).^(1.0 / (gamma + 1))
# AvgDegree = 2 * m * mean(AgentAct)
# AgentState = zeros(Bool, N)
# AgentState[1:Ni] .= 1
etaS = 5
eta = 5

# grid=Reshape([(eta, lambda) for eta in Eta, lambda in Lambda],:)
# @showprogress for (ietaI, eta) in enumerate(Eta)
InfectedOverATrial = zeros(0)
InfectedOverTrials = zeros(0)

@showprogress for (ilambda, lambda) in enumerate(Lambda)
# lambda = 1
    Activitylist = zeros(0)

    AgentState = zeros(Bool, N)
    AgentState[1:Int(0.01 * N)] .= 1
    InfectedOverATrial = 0.0
    for i_time = 1:tmax
        AgentStateNext = AgentState
        for i = 1:N
            if AgentState[i] == 1
                r = rand()
                if (r < mu)
                    AgentStateNext[i] = 0
                end
            end
            Activity = AgentAct[i]
            r = rand()
            if (Activity > r)
                for j = 1:m
                    rInt = rand(1:N)
                    if (AgentState[i] == 1)
                        r = rand()
                        if ((AgentState[rInt] == 0) && (r < lambda))
                        # println(rInt)
                            AgentStateNext[rInt] = 1
                        end
                    else
                        r = rand()
                        if ((AgentState[rInt] == 1) && (r < lambda))
                        # println(rInt)
                            AgentStateNext[i] = 1
                        end
                    end
                end
            end
        end
        NumInf = 0
        for k = 1:N
            AgentState[k] = AgentStateNext[k]
            if (AgentState[k] == 1)
                NumInf = NumInf + 1
            end
        end
        push!(InfectedOverTrials, NumInf)
    end
    InfectedOverATrial = InfectedOverATrial / TWindow
    push!(InfectedOverTrials, InfectedOverATrial)
end

mat"figure();set(gcf, 'Position',  [0, 0, 2500, 1500]); hold on;"
# mat"axis([ -0,1,-0,10 ])"
# mat"imagesc([0,1],[0,10],$AvgInfectedEtaLambda)"
# mat"colormap jet"
#
# mat"plot($InfectedOverTrials,'.-','MarkerSize',20)"
mat"[X,Y]=meshgrid($Eta,$Lambda);contourf(X,Y,$AvgInfectedEtaLambda,100)"

mat"colorbar;colormap('inferno')"


text:
# Julia wrapper for header: sc.h
# Automatically generated using Clang.jl

function sc_extern_c_hack_3()
    ccall((:sc_extern_c_hack_3, p4est_lib), Cvoid, ())
end

function sc_malloc(package, size)
    ccall((:sc_malloc, p4est_lib), Ptr{Cvoid}, (Cint, Cint), package, size)
end

function sc_calloc(package, nmemb, size)
    ccall((:sc_calloc, p4est_lib), Ptr{Cvoid}, (Cint, Cint, Cint), package, nmemb, size)
end

function sc_realloc(package, ptr, size)
    ccall((:sc_realloc, p4est_lib), Ptr{Cvoid}, (Cint, Ptr{Cvoid}, Cint), package, ptr, size)
end

function sc_strdup(package, s)
    ccall((:sc_strdup, p4est_lib), Cstring, (Cint, Cstring), package, s)
end

function sc_free(package, ptr)
    ccall((:sc_free, p4est_lib), Cvoid, (Cint, Ptr{Cvoid}), package, ptr)
end

function sc_memory_status(package)
    ccall((:sc_memory_status, p4est_lib), Cint, (Cint,), package)
end

function sc_memory_check(package)
    ccall((:sc_memory_check, p4est_lib), Cvoid, (Cint,), package)
end

function sc_int_compare(v1, v2)
    ccall((:sc_int_compare, p4est_lib), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), v1, v2)
end

function sc_int8_compare(v1, v2)
    ccall((:sc_int8_compare, p4est_lib), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), v1, v2)
end

function sc_int16_compare(v1, v2)
    ccall((:sc_int16_compare, p4est_lib), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), v1, v2)
end

function sc_int32_compare(v1, v2)
    ccall((:sc_int32_compare, p4est_lib), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), v1, v2)
end

function sc_int64_compare(v1, v2)
    ccall((:sc_int64_compare, p4est_lib), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), v1, v2)
end

function sc_double_compare(v1, v2)
    ccall((:sc_double_compare, p4est_lib), Cint, (Ptr{Cvoid}, Ptr{Cvoid}), v1, v2)
end

function sc_set_log_defaults(log_stream, log_handler, log_thresold)
    ccall((:sc_set_log_defaults, p4est_lib), Cvoid, (Ptr{FILE}, sc_log_handler_t, Cint), log_stream, log_handler, log_thresold)
end

function sc_set_abort_handler(abort_handler)
    ccall((:sc_set_abort_handler, p4est_lib), Cvoid, (sc_abort_handler_t,), abort_handler)
end

function sc_log(filename, lineno, package, category, priority, msg)
    ccall((:sc_log, p4est_lib), Cvoid, (Cstring, Cint, Cint, Cint, Cint, Cstring), filename, lineno, package, category, priority, msg)
end

function sc_log_indent_push_count(package, count)
    ccall((:sc_log_indent_push_count, p4est_lib), Cvoid, (Cint, Cint), package, count)
end

function sc_log_indent_pop_count(package, count)
    ccall((:sc_log_indent_pop_count, p4est_lib), Cvoid, (Cint, Cint), package, count)
end

function sc_log_indent_push()
    ccall((:sc_log_indent_push, p4est_lib), Cvoid, ())
end

function sc_log_indent_pop()
    ccall((:sc_log_indent_pop, p4est_lib), Cvoid, ())
end

function sc_abort()
    ccall((:sc_abort, p4est_lib), Cvoid, ())
end

function sc_abort_verbose(filename, lineno, msg)
    ccall((:sc_abort_verbose, p4est_lib), Cvoid, (Cstring, Cint, Cstring), filename, lineno, msg)
end

function sc_abort_collective(msg)
    ccall((:sc_abort_collective, p4est_lib), Cvoid, (Cstring,), msg)
end

function sc_package_register(log_handler, log_threshold, name, full)
    ccall((:sc_package_register, p4est_lib), Cint, (sc_log_handler_t, Cint, Cstring, Cstring), log_handler, log_threshold, name, full)
end

function sc_package_is_registered(package_id)
    ccall((:sc_package_is_registered, p4est_lib), Cint, (Cint,), package_id)
end

function sc_package_lock(package_id)
    ccall((:sc_package_lock, p4est_lib), Cvoid, (Cint,), package_id)
end

function sc_package_unlock(package_id)
    ccall((:sc_package_unlock, p4est_lib), Cvoid, (Cint,), package_id)
end

function sc_package_set_verbosity(package_id, log_priority)
    ccall((:sc_package_set_verbosity, p4est_lib), Cvoid, (Cint, Cint), package_id, log_priority)
end

function sc_package_set_abort_alloc_mismatch(package_id, set_abort)
    ccall((:sc_package_set_abort_alloc_mismatch, p4est_lib), Cvoid, (Cint, Cint), package_id, set_abort)
end

function sc_package_unregister(package_id)
    ccall((:sc_package_unregister, p4est_lib), Cvoid, (Cint,), package_id)
end

function sc_package_print_summary(log_priority)
    ccall((:sc_package_print_summary, p4est_lib), Cvoid, (Cint,), log_priority)
end

function sc_init(mpicomm, catch_signals, print_backtrace, log_handler, log_threshold)
    ccall((:sc_init, p4est_lib), Cvoid, (MPI_Comm, Cint, Cint, sc_log_handler_t, Cint), mpicomm, catch_signals, print_backtrace, log_handler, log_threshold)
end

function sc_finalize()
    ccall((:sc_finalize, p4est_lib), Cvoid, ())
end

function sc_is_root()
    ccall((:sc_is_root, p4est_lib), Cint, ())
end

function sc_extern_c_hack_4()
    ccall((:sc_extern_c_hack_4, p4est_lib), Cvoid, ())
end
# Julia wrapper for header: sc_containers.h
# Automatically generated using Clang.jl

function sc_array_memory_used()
    ccall((:sc_array_memory_used, p4est_lib), Cint, ())
end

function sc_array_new(elem_size)
    ccall((:sc_array_new, p4est_lib), Ptr{sc_array_t}, (Cint,), elem_size)
end

function sc_array_new_count(elem_size, elem_count)
    ccall((:sc_array_new_count, p4est_lib), Ptr{sc_array_t}, (Cint, Cint), elem_size, elem_count)
end

function sc_array_new_view(array, offset, length)
    ccall((:sc_array_new_view, p4est_lib), Ptr{sc_array_t}, (Ptr{sc_array_t}, Cint, Cint), array, offset, length)
end

function sc_array_new_data(base, elem_size, elem_count)
    ccall((:sc_array_new_data, p4est_lib), Ptr{sc_array_t}, (Ptr{Cvoid}, Cint, Cint), base, elem_size, elem_count)
end

function sc_array_destroy(array)
    ccall((:sc_array_destroy, p4est_lib), Cvoid, (Ptr{sc_array_t},), array)
end

function sc_array_destroy_null(parray)
    ccall((:sc_array_destroy_null, p4est_lib), Cvoid, (Ptr{Ptr{sc_array_t}},), parray)
end

function sc_array_init(array, elem_size)
    ccall((:sc_array_init, p4est_lib), Cvoid, (Ptr{sc_array_t}, Cint), array, elem_size)
end

function sc_array_init_size(array, elem_size, elem_count)
    ccall((:sc_array_init_size, p4est_lib), Cvoid, (Ptr{sc_array_t}, Cint, Cint), array, elem_size, elem_count)
end

function sc_array_init_count(array, elem_size, elem_count)
    ccall((:sc_array_init_count, p4est_lib), Cvoid, (Ptr{sc_array_t}, Cint, Cint), array, elem_size, elem_count)
end

function sc_array_init_view(view, array, offset, length)
    ccall((:sc_array_init_view, p4est_lib), Cvoid, (Ptr{sc_array_t}, Ptr{sc_array_t}, Cint, Cint), view, array, offset, length)
end

function sc_array_init_data(view, base, elem_size, elem_count)
    ccall((:sc_array_init_data, p4est_lib), Cvoid, (Ptr{sc_array_t}, Ptr{Cvoid}, Cint, Cint), view, base, elem_size, elem_count)
end

function sc_array_memset(array, c)
    ccall((:sc_array_memset, p4est_lib), Cvoid, (Ptr{sc_array_t}, Cint), array, c)
end

function sc_array_reset(array)
    ccall((:sc_array_reset, p4est_lib), Cvoid, (Ptr{sc_array_t},), array)
end

function sc_array_truncate(array)
    ccall((:sc_array_truncate, p4est_lib), Cvoid, (Ptr{sc_array_t},), array)
end

function sc_array_rewind(array, new_count)
    ccall((:sc_array_rewind, p4est_lib), Cvoid, (Ptr{sc_array_t}, Cint), array, new_count)
end

function sc_array_resize(array, new_count)
    ccall((:sc_array_resize, p4est_lib), Cvoid, (Ptr{sc_array_t}, Cint), array, new_count)
end

function sc_array_copy(dest, src)
    ccall((:sc_array_copy, p4est_lib), Cvoid, (Ptr{sc_array_t}, Ptr{sc_array_t}), dest, src)
end

function sc_array_copy_into(dest, dest_offset, src)
    ccall((:sc_array_copy_into, p4est_lib), Cvoid, (Ptr{sc_array_t}, Cint, Ptr{sc_array_t}), dest, dest_offset, src)
end

function sc_array_move_part(dest, dest_offset, src, src_offset, count)
    ccall((:sc_array_move_part, p4est_lib), Cvoid, (Ptr{sc_array_t}, Cint, Ptr{sc_array_t}, Cint, Cint), dest, dest_offset, src, src_offset, count)
end

function sc_array_sort(array, compar)
    ccall((:sc_array_sort, p4est_lib), Cvoid, (Ptr{sc_array_t}, Ptr{Cvoid}), array, compar)
end

function sc_array_is_sorted(array, compar)
    ccall((:sc_array_is_sorted, p4est_lib), Cint, (Ptr{sc_array_t}, Ptr{Cvoid}), array, compar)
end

function sc_array_is_equal(array, other)
    ccall((:sc_array_is_equal, p4est_lib), Cint, (Ptr{sc_array_t}, Ptr{sc_array_t}), array, other)
end

function sc_array_uniq(array, compar)
    ccall((:sc_array_uniq, p4est_lib), Cvoid, (Ptr{sc_array_t}, Ptr{Cvoid}), array, compar)
end

function sc_array_bsearch(array, key, compar)
    ccall((:sc_array_bsearch, p4est_lib), ssize_t, (Ptr{sc_array_t}, Ptr{Cvoid}, Ptr{Cvoid}), array, key, compar)
end

function sc_array_split(array, offsets, num_types, type_fn, data)
    ccall((:sc_array_split, p4est_lib), Cvoid, (Ptr{sc_array_t}, Ptr{sc_array_t}, Csize_t, Cint, Ptr{Cvoid}), array, offsets, num_types, type_fn, data)
end

function sc_array_is_permutation(array)
    ccall((:sc_array_is_permutation, p4est_lib), Cint, (Ptr{sc_array_t},), array)
end

function sc_array_permute(array, newindices, keepperm)
    ccall((:sc_array_permute, p4est_lib), Cvoid, (Ptr{sc_array_t}, Ptr{sc_array_t}, Cint), array, newindices, keepperm)
end

function sc_array_checksum(array)
    ccall((:sc_array_checksum, p4est_lib), UInt32, (Ptr{sc_array_t},), array)
end

function sc_array_pqueue_add()
    ccall((:sc_array_pqueue_add, p4est_lib), Cint, ())
end

function sc_array_pqueue_pop()
    ccall((:sc_array_pqueue_pop, p4est_lib), Cint, ())
end

function sc_mstamp_init(mst, stamp_unit, elem_size)
    ccall((:sc_mstamp_init, p4est_lib), Cvoid, (Ptr{sc_mstamp_t}, Csize_t, Csize_t), mst, stamp_unit, elem_size)
end

function sc_mstamp_reset(mst)
    ccall((:sc_mstamp_reset, p4est_lib), Cvoid, (Ptr{sc_mstamp_t},), mst)
end

function sc_mstamp_truncate(mst)
    ccall((:sc_mstamp_truncate, p4est_lib), Cvoid, (Ptr{sc_mstamp_t},), mst)
end

function sc_mstamp_alloc(mst)
    ccall((:sc_mstamp_alloc, p4est_lib), Ptr{Cvoid}, (Ptr{sc_mstamp_t},), mst)
end

function sc_mstamp_memory_used()
    ccall((:sc_mstamp_memory_used, p4est_lib), Cint, ())
end

function sc_mempool_memory_used()
    ccall((:sc_mempool_memory_used, p4est_lib), Cint, ())
end

function sc_mempool_new(elem_size)
    ccall((:sc_mempool_new, p4est_lib), Ptr{sc_mempool_t}, (Csize_t,), elem_size)
end

function sc_mempool_new_zero_and_persist(elem_size)
    ccall((:sc_mempool_new_zero_and_persist, p4est_lib), Ptr{sc_mempool_t}, (Csize_t,), elem_size)
end

function sc_mempool_init(mempool, elem_size)
    ccall((:sc_mempool_init, p4est_lib), Cvoid, (Ptr{sc_mempool_t}, Csize_t), mempool, elem_size)
end

function sc_mempool_destroy(mempool)
    ccall((:sc_mempool_destroy, p4est_lib), Cvoid, (Ptr{sc_mempool_t},), mempool)
end

function sc_mempool_destroy_null(pmempool)
    ccall((:sc_mempool_destroy_null, p4est_lib), Cvoid, (Ptr{Ptr{sc_mempool_t}},), pmempool)
end

function sc_mempool_reset(mempool)
    ccall((:sc_mempool_reset, p4est_lib), Cvoid, (Ptr{sc_mempool_t},), mempool)
end

function sc_mempool_truncate(mempool)
    ccall((:sc_mempool_truncate, p4est_lib), Cvoid, (Ptr{sc_mempool_t},), mempool)
end

function sc_list_memory_used()
    ccall((:sc_list_memory_used, p4est_lib), Cint, ())
end

function sc_list_new(allocator)
    ccall((:sc_list_new, p4est_lib), Ptr{sc_list_t}, (Ptr{sc_mempool_t},), allocator)
end

function sc_list_destroy(list)
    ccall((:sc_list_destroy, p4est_lib), Cvoid, (Ptr{sc_list_t},), list)
end

function sc_list_init(list, allocator)
    ccall((:sc_list_init, p4est_lib), Cvoid, (Ptr{sc_list_t}, Ptr{sc_mempool_t}), list, allocator)
end

function sc_list_reset(list)
    ccall((:sc_list_reset, p4est_lib), Cvoid, (Ptr{sc_list_t},), list)
end

function sc_list_unlink(list)
    ccall((:sc_list_unlink, p4est_lib), Cvoid, (Ptr{sc_list_t},), list)
end

function sc_list_prepend(list, data)
    ccall((:sc_list_prepend, p4est_lib), Ptr{sc_link_t}, (Ptr{sc_list_t}, Ptr{Cvoid}), list, data)
end

function sc_list_append(list, data)
    ccall((:sc_list_append, p4est_lib), Ptr{sc_link_t}, (Ptr{sc_list_t}, Ptr{Cvoid}), list, data)
end

function sc_list_insert(list, pred, data)
    ccall((:sc_list_insert, p4est_lib), Ptr{sc_link_t}, (Ptr{sc_list_t}, Ptr{sc_link_t}, Ptr{Cvoid}), list, pred, data)
end

function sc_list_remove(list, pred)
    ccall((:sc_list_remove, p4est_lib), Ptr{Cvoid}, (Ptr{sc_list_t}, Ptr{sc_link_t}), list, pred)
end

function sc_list_pop(list)
    ccall((:sc_list_pop, p4est_lib), Ptr{Cvoid}, (Ptr{sc_list_t},), list)
end

function sc_hash_function_string(s, u)
    ccall((:sc_hash_function_string, p4est_lib), UInt32, (Ptr{Cvoid}, Ptr{Cvoid}), s, u)
end

function sc_hash_memory_used()
    ccall((:sc_hash_memory_used, p4est_lib), Cint, ())
end

function sc_hash_new(hash_fn, equal_fn, user_data, allocator)
    ccall((:sc_hash_new, p4est_lib), Ptr{sc_hash_t}, (sc_hash_function_t, sc_equal_function_t, Ptr{Cvoid}, Ptr{sc_mempool_t}), hash_fn, equal_fn, user_data, allocator)
end

function sc_hash_destroy(hash)
    ccall((:sc_hash_destroy, p4est_lib), Cvoid, (Ptr{sc_hash_t},), hash)
end

function sc_hash_destroy_null(phash)
    ccall((:sc_hash_destroy_null, p4est_lib), Cvoid, (Ptr{Ptr{sc_hash_t}},), phash)
end

function sc_hash_truncate(hash)
    ccall((:sc_hash_truncate, p4est_lib), Cvoid, (Ptr{sc_hash_t},), hash)
end

function sc_hash_unlink(hash)
    ccall((:sc_hash_unlink, p4est_lib), Cvoid, (Ptr{sc_hash_t},), hash)
end

function sc_hash_unlink_destroy(hash)
    ccall((:sc_hash_unlink_destroy, p4est_lib), Cvoid, (Ptr{sc_hash_t},), hash)
end

function sc_hash_lookup(hash, v, found)
    ccall((:sc_hash_lookup, p4est_lib), Cint, (Ptr{sc_hash_t}, Ptr{Cvoid}, Ptr{Ptr{Ptr{Cvoid}}}), hash, v, found)
end

function sc_hash_insert_unique(hash, v, found)
    ccall((:sc_hash_insert_unique, p4est_lib), Cint, (Ptr{sc_hash_t}, Ptr{Cvoid}, Ptr{Ptr{Ptr{Cvoid}}}), hash, v, found)
end

function sc_hash_remove(hash, v, found)
    ccall((:sc_hash_remove, p4est_lib), Cint, (Ptr{sc_hash_t}, Ptr{Cvoid}, Ptr{Ptr{Cvoid}}), hash, v, found)
end

function sc_hash_foreach(hash, fn)
    ccall((:sc_hash_foreach, p4est_lib), Cvoid, (Ptr{sc_hash_t}, sc_hash_foreach_t), hash, fn)
end

function sc_hash_print_statistics(package_id, log_priority, hash)
    ccall((:sc_hash_print_statistics, p4est_lib), Cvoid, (Cint, Cint, Ptr{sc_hash_t}), package_id, log_priority, hash)
end

function sc_hash_array_memory_used()
    ccall((:sc_hash_array_memory_used, p4est_lib), Cint, ())
end

function sc_hash_array_new(elem_size, hash_fn, equal_fn, user_data)
    ccall((:sc_hash_array_new, p4est_lib), Ptr{sc_hash_array_t}, (Csize_t, sc_hash_function_t, sc_equal_function_t, Ptr{Cvoid}), elem_size, hash_fn, equal_fn, user_data)
end

function sc_hash_array_destroy(hash_array)
    ccall((:sc_hash_array_destroy, p4est_lib), Cvoid, (Ptr{sc_hash_array_t},), hash_array)
end

function sc_hash_array_is_valid(hash_array)
    ccall((:sc_hash_array_is_valid, p4est_lib), Cint, (Ptr{sc_hash_array_t},), hash_array)
end

function sc_hash_array_truncate(hash_array)
    ccall((:sc_hash_array_truncate, p4est_lib), Cvoid, (Ptr{sc_hash_array_t},), hash_array)
end

function sc_hash_array_lookup(hash_array, v, position)
    ccall((:sc_hash_array_lookup, p4est_lib), Cint, (Ptr{sc_hash_array_t}, Ptr{Cvoid}, Ptr{Csize_t}), hash_array, v, position)
end

function sc_hash_array_insert_unique(hash_array, v, position)
    ccall((:sc_hash_array_insert_unique, p4est_lib), Ptr{Cvoid}, (Ptr{sc_hash_array_t}, Ptr{Cvoid}, Ptr{Csize_t}), hash_array, v, position)
end

function sc_hash_array_rip(hash_array, rip)
    ccall((:sc_hash_array_rip, p4est_lib), Cvoid, (Ptr{sc_hash_array_t}, Ptr{sc_array_t}), hash_array, rip)
end

function sc_recycle_array_init(rec_array, elem_size)
    ccall((:sc_recycle_array_init, p4est_lib), Cvoid, (Ptr{sc_recycle_array_t}, Csize_t), rec_array, elem_size)
end

function sc_recycle_array_reset(rec_array)
    ccall((:sc_recycle_array_reset, p4est_lib), Cvoid, (Ptr{sc_recycle_array_t},), rec_array)
end

function sc_recycle_array_insert(rec_array, position)
    ccall((:sc_recycle_array_insert, p4est_lib), Ptr{Cvoid}, (Ptr{sc_recycle_array_t}, Ptr{Csize_t}), rec_array, position)
end

function sc_recycle_array_remove(rec_array, position)
    ccall((:sc_recycle_array_remove, p4est_lib), Ptr{Cvoid}, (Ptr{sc_recycle_array_t}, Csize_t), rec_array, position)
end
# Julia wrapper for header: sc_io.h
# Automatically generated using Clang.jl

function sc_io_sink_destroy(sink)
    ccall((:sc_io_sink_destroy, p4est_lib), Cint, (Ptr{sc_io_sink_t},), sink)
end

function sc_io_sink_write(sink, data, bytes_avail)
    ccall((:sc_io_sink_write, p4est_lib), Cint, (Ptr{sc_io_sink_t}, Ptr{Cvoid}, Csize_t), sink, data, bytes_avail)
end

function sc_io_sink_complete(sink, bytes_in, bytes_out)
    ccall((:sc_io_sink_complete, p4est_lib), Cint, (Ptr{sc_io_sink_t}, Ptr{Csize_t}, Ptr{Csize_t}), sink, bytes_in, bytes_out)
end

function sc_io_sink_align(sink, bytes_align)
    ccall((:sc_io_sink_align, p4est_lib), Cint, (Ptr{sc_io_sink_t}, Csize_t), sink, bytes_align)
end

function sc_io_source_destroy(source)
    ccall((:sc_io_source_destroy, p4est_lib), Cint, (Ptr{sc_io_source_t},), source)
end

function sc_io_source_read(source, data, bytes_avail, bytes_out)
    ccall((:sc_io_source_read, p4est_lib), Cint, (Ptr{sc_io_source_t}, Ptr{Cvoid}, Csize_t, Ptr{Csize_t}), source, data, bytes_avail, bytes_out)
end

function sc_io_source_complete(source, bytes_in, bytes_out)
    ccall((:sc_io_source_complete, p4est_lib), Cint, (Ptr{sc_io_source_t}, Ptr{Csize_t}, Ptr{Csize_t}), source, bytes_in, bytes_out)
end

function sc_io_source_align(source, bytes_align)
    ccall((:sc_io_source_align, p4est_lib), Cint, (Ptr{sc_io_source_t}, Csize_t), source, bytes_align)
end

function sc_io_source_activate_mirror(source)
    ccall((:sc_io_source_activate_mirror, p4est_lib), Cint, (Ptr{sc_io_source_t},), source)
end

function sc_io_source_read_mirror(source, data, bytes_avail, bytes_out)
    ccall((:sc_io_source_read_mirror, p4est_lib), Cint, (Ptr{sc_io_source_t}, Ptr{Cvoid}, Csize_t, Ptr{Csize_t}), source, data, bytes_avail, bytes_out)
end

function sc_vtk_write_binary(vtkfile, numeric_data, byte_length)
    ccall((:sc_vtk_write_binary, p4est_lib), Cint, (Ptr{FILE}, Cstring, Csize_t), vtkfile, numeric_data, byte_length)
end

function sc_vtk_write_compressed(vtkfile, numeric_data, byte_length)
    ccall((:sc_vtk_write_compressed, p4est_lib), Cint, (Ptr{FILE}, Cstring, Csize_t), vtkfile, numeric_data, byte_length)
end

function sc_fwrite(ptr, size, nmemb, file, errmsg)
    ccall((:sc_fwrite, p4est_lib), Cvoid, (Ptr{Cvoid}, Csize_t, Csize_t, Ptr{FILE}, Cstring), ptr, size, nmemb, file, errmsg)
end

function sc_fread(ptr, size, nmemb, file, errmsg)
    ccall((:sc_fread, p4est_lib), Cvoid, (Ptr{Cvoid}, Csize_t, Csize_t, Ptr{FILE}, Cstring), ptr, size, nmemb, file, errmsg)
end

function sc_mpi_write(mpifile, ptr, zcount, t, errmsg)
    ccall((:sc_mpi_write, p4est_lib), Cvoid, (MPI_File, Ptr{Cvoid}, Csize_t, MPI_Datatype, Cstring), mpifile, ptr, zcount, t, errmsg)
end


text:
## Bulge chasing for the RDS and CSS case

## A bulge is introduced by a unitary transformation. In the single shift case this is just U^t A U; for
## the double shift case this is (VU)^T A (VU).
## In either case, we have three steps:
##
## * create the unitary transformations (U or (V,U)) -- create_bulge
## * absorb the left side -- absorb_Ut
## * chase the right side down until it is absorbed -- absorb_U
##
function bulge_step(QF::QFactorization{T, S, Vt}, RF, storage, ctr, m=0) where {T, S, Vt}

    create_bulge(QF, RF, storage, ctr)
    absorb_Ut(QF, RF, storage, ctr) # see CSS or RDS
    absorb_U(QF, RF, storage, ctr)

    return nothing
end

##################################################
## Absorb U
## After U absorbtion we have either
## A * U or (W,A) * V * U
## THis chases U or (V*U) through until the value can be absorbed
##
## This just pushes work to passthrough_triu and passthrough_Q
##
## * passthrough_triu is basically passthrough(RF) with some efficiencies for first
## time through
##
## * passthrough_Q passes through chain and return false, until fusion when
## it returns true
##
function absorb_U(QF::QFactorization{T, S, Vt}, RF, storage, ctr) where {T, S, Vt}

    flag = false
    while !flag

        passthrough_triu(QF, RF, storage, ctr, Val(:right))
        flag = passthrough_Q(QF, RF, storage,  ctr, Val(:right))
        ## implicit unitary operation to move U from left to right

    end

    return nothing

end