text:
# ===============================
# Written by AAE, TU Delft
# simulkade.com
# ===============================

# ================== Linear averaging scheme ==================
function linearMean(phi::CellValue)
# calculates the average values of a cell variable. The output is a
# face variable
    phi_face = createFaceVariable(phi.domain, 0.0)
    linearMean!(phi, phi_face)
    return phi_face

# d=phi.domain.dimension
# if d==1 || d==1.5
#   dx = phi.domain.cellsize.x
#   FaceValue(phi.domain,
#     (dx[2:end].*phi.value[1:end-1]+dx[1:end-1].*phi.value[2:end])./(dx[2:end]+dx[1:end-1]),
#     [1.0],
#     [1.0])
# elseif d==2 || d==2.5 || d==2.8
#   dx = phi.domain.cellsize.x
#   Ny = phi.domain.dims[2]
#   dy = zeros( 1, Ny+2)
#   dy[:] = phi.domain.cellsize.y
#   FaceValue(phi.domain,
#     (dx[2:end].*phi.value[1:end-1,2:end-1]+dx[1:end-1].*phi.value[2:end,2:end-1])./(dx[2:end]+dx[1:end-1]),
#     (dy[:,2:end].*phi.value[2:end-1,1:end-1]+dy[:,1:end-1].*phi.value[2:end-1,2:end])./(dy[:,2:end]+dy[:,1:end-1]),
#     [1.0])
# elseif d==3 || d==3.2
#   Ny = phi.domain.dims[2]
#   Nz = phi.domain.dims[3]
#   dx = phi.domain.cellsize.x
#   dy= zeros( 1, Ny+2)
#   dy[:] = phi.domain.cellsize.y
#   dz= zeros( 1, 1, Nz+2)
#   dz[:] = phi.domain.cellsize.z
#   FaceValue(phi.domain,
#     (dx[2:end].*phi.value[1:end-1,2:end-1,2:end-1]+dx[1:end-1].*phi.value[2:end,2:end-1,2:end-1])./(dx[2:end]+dx[1:end-1]),
#     (dy[:,2:end].*phi.value[2:end-1,1:end-1,2:end-1]+dy[:,1:end-1].*phi.value[2:end-1,2:end,2:end-1])./(dy[:,1:end-1]+dy[:,2:end]),
#     (dz[:,:,2:end].*phi.value[2:end-1,2:end-1,1:end-1]+dz[:,:,1:end-1].*phi.value[2:end-1,2:end-1,2:end])./(dz[:,:,1:end-1]+dz[:,:,2:end]))
# end
end

function linearMean!(phi::CellValue, phi_face::FaceValue)
    # calculates the average values of a cell variable. The output is a
    # face variable
    d=phi.domain.dimension
    if d==1 || d==1.5
        dx = phi.domain.cellsize.x

        phi_face.xvalue .= (dx[2:end].*phi.value[1:end-1]+dx[1:end-1].*phi.value[2:end])./(dx[2:end]+dx[1:end-1])
        phi_face.yvalue .= [1.0]
        phi_face.zvalue .= [1.0]
    elseif d==2 || d==2.5 || d==2.8
        dx = phi.domain.cellsize.x
        Ny = phi.domain.dims[2]
        dy = zeros( 1, Ny+2)
        dy[:] = phi.domain.cellsize.y

        phi_face.xvalue .= (dx[2:end].*phi.value[1:end-1,2:end-1]+dx[1:end-1].*phi.value[2:end,2:end-1])./(dx[2:end]+dx[1:end-1])
        phi_face.yvalue .= (dy[:,2:end].*phi.value[2:end-1,1:end-1]+dy[:,1:end-1].*phi.value[2:end-1,2:end])./(dy[:,2:end]+dy[:,1:end-1])
        phi_face.zvalue .= [1.0]
    elseif d==3 || d==3.2
        Ny = phi.domain.dims[2]
        Nz = phi.domain.dims[3]
        dx = phi.domain.cellsize.x
        dy= zeros( 1, Ny+2)
        dy[:] = phi.domain.cellsize.y
        dz= zeros( 1, 1, Nz+2)
        dz[:] = phi.domain.cellsize.z

        phi_face.xvalue .= (dx[2:end].*phi.value[1:end-1,2:end-1,2:end-1].+dx[1:end-1].*phi.value[2:end,2:end-1,2:end-1])./(dx[2:end]+dx[1:end-1])
        phi_face.yvalue .= (dy[:,2:end].*phi.value[2:end-1,1:end-1,2:end-1]+dy[:,1:end-1].*phi.value[2:end-1,2:end,2:end-1])./(dy[:,1:end-1]+dy[:,2:end])
        phi_face.zvalue .= (dz[:,:,2:end].*phi.value[2:end-1,2:end-1,1:end-1]+dz[:,:,1:end-1].*phi.value[2:end-1,2:end-1,2:end])./(dz[:,:,1:end-1]+dz[:,:,2:end])
    end
end

# ================== Aritmetic averaging scheme ==================
function arithmeticMean(phi::CellValue)
# calculates the average values of a cell variable. The output is a
# face variable
d=phi.domain.dimension
if d==1 || d==1.5
  dx = phi.domain.cellsize.x
  FaceValue(phi.domain,
    (dx[1:end-1].*phi.value[1:end-1]+dx[2:end].*phi.value[2:end])./(dx[2:end]+dx[1:end-1]),
    [1.0],
    [1.0])
elseif d==2 || d==2.5 || d==2.8
  dx = phi.domain.cellsize.x
  Ny = phi.domain.dims[2]
  dy = zeros( 1, Ny+2)
  dy[:] = phi.domain.cellsize.y
  FaceValue(phi.domain,
    (dx[1:end-1].*phi.value[1:end-1,2:end-1]+dx[2:end].*phi.value[2:end,2:end-1])./(dx[2:end]+dx[1:end-1]),
    (dy[:,1:end-1].*phi.value[2:end-1,1:end-1]+dy[:,2:end].*phi.value[2:end-1,2:end])./(dy[:,2:end]+dy[:,1:end-1]),
    [1.0])
elseif d==3 || d==3.2
  Ny = phi.domain.dims[2]
  Nz = phi.domain.dims[3]
  dx = phi.domain.cellsize.x
  dy= zeros( 1, Ny+2)
  dy[:] = phi.domain.cellsize.y
  dz= zeros( 1, 1, Nz+2)
  dz[:] = phi.domain.cellsize.z
  FaceValue(phi.domain,
    (dx[1:end-1].*phi.value[1:end-1,2:end-1,2:end-1]+dx[2:end].*phi.value[2:end,2:end-1,2:end-1])./(dx[2:end]+dx[1:end-1]),
    (dy[:,1:end-1].*phi.value[2:end-1,1:end-1,2:end-1]+dy[:,2:end].*phi.value[2:end-1,2:end,2:end-1])./(dy[:,1:end-1]+dy[:,2:end]),
    (dz[:,:,1:end-1].*phi.value[2:end-1,2:end-1,1:end-1]+dz[:,:,2:end].*phi.value[2:end-1,2:end-1,2:end])./(dz[:,:,1:end-1]+dz[:,:,2:end]))
end
end

# ================== Geometric averaging scheme ==================
function geometricMean(phi::CellValue)
# calculates the average values of a cell variable. The output is a
# face variable
d=phi.domain.dimension
if d==1 || d==1.5
  dx = phi.domain.cellsize.x
  n= phi.domain.dims[1]
  phix=zeros(n+1)
  for i=1:n+1
      if phi.value[i]==0.0 || phi.value[i+1]==0.0
          phix[i]=0.0
      else
          phix[i]=exp.((dx[i]*log.(phi.value[i])+dx[i+1]*log.(phi.value[i+1]))/(dx[i+1]+dx[i]))
      end
  end
  FaceValue(phi.domain,
    phix,
    [1.0],
    [1.0])
elseif d==2 || d==2.5 || d==2.8
  dx = phi.domain.cellsize.x
  Ny = phi.domain.dims[2]
  dy = zeros( 1, Ny+2)
  dy[:] = phi.domain.cellsize.y
  FaceValue(phi.domain,
    exp.((dx[1:end-1].*log.(phi.value[1:end-1,2:end-1])+dx[2:end].*log.(phi.value[2:end,2:end-1]))./(dx[2:end]+dx[1:end-1])),
    exp.((dy[:,1:end-1].*log.(phi.value[2:end-1,1:end-1])+dy[:,2:end].*log.(phi.value[2:end-1,2:end]))./(dy[:,2:end]+dy[:,1:end-1])),
    [1.0])
elseif d==3 || d==3.2
  Ny = phi.domain.dims[2]
  Nz = phi.domain.dims[3]
  dx = phi.domain.cellsize.x
  dy= zeros( 1, Ny+2)
  dy[:] = phi.domain.cellsize.y
  dz= zeros( 1, 1, Nz+2)
  dz[:] = phi.domain.cellsize.z
  FaceValue(phi.domain,
    exp.((dx[1:end-1].*log.(phi.value[1:end-1,2:end-1,2:end-1])+dx[2:end].*log.(phi.value[2:end,2:end-1,2:end-1]))./(dx[2:end]+dx[1:end-1])),
    exp.((dy[:,1:end-1].*log.(phi.value[2:end-1,1:end-1,2:end-1])+dy[:,2:end].*log.(phi.value[2:end-1,2:end,2:end-1]))./(dy[:,1:end-1]+dy[:,2:end])),
    exp.((dz[:,:,1:end-1].*log.(phi.value[2:end-1,2:end-1,1:end-1])+dz[:,:,2:end].*log.(phi.value[2:end-1,2:end-1,2:end]))./(dz[:,:,1:end-1]+dz[:,:,2:end])))
end
end

# ================== Harmonic averaging scheme ==================
function harmonicMean(phi::CellValue)
# calculates the average values of a cell variable. The output is a
# face variable
d=phi.domain.dimension
if d==1 || d==1.5
  dx = phi.domain.cellsize.x
  n=phi.domain.dims[1]
  phix=zeros(n+1)
  for i=1:n+1
      if phi.value[i]==0.0 || phi.value[i+1]==0.0
          phix[i]=0.0
      else
          phix[i]=(dx[i+1]+dx[i])/(dx[i+1]/phi.value[i+1]+dx[i]/phi.value[i])
      end
  end
  FaceValue(phi.domain,
    phix,
    [1.0],
    [1.0])
elseif d==2 || d==2.5 || d==2.8
  dx = phi.domain.cellsize.x
  Ny = phi.domain.dims[2]
  dy = zeros( 1, Ny+2)
  dy[:] = phi.domain.cellsize.y
  FaceValue(phi.domain,
    phi.value[2:end,2:end-1].*phi.value[1:end-1,2:end-1].*(dx[2:end]+dx[1:end-1])./(dx[2:end].*phi.value[1:end-1,2:end-1]+dx[1:end-1].*phi.value[2:end,2:end-1]),
    phi.value[2:end-1,2:end].*phi.value[2:end-1,1:end-1].*(dy[:,2:end]+dy[:,1:end-1])./(dy[:,2:end].*phi.value[2:end-1,1:end-1]+dy[:,1:end-1].*phi.value[2:end-1,2:end]),
    [1.0])
elseif d==3 || d==3.2
  Ny = phi.domain.dims[2]
  Nz = phi.domain.dims[3]
  dx = phi.domain.cellsize.x
  dy= zeros( 1, Ny+2)
  dy[:] = phi.domain.cellsize.y
  dz= zeros( 1, 1, Nz+2)
  dz[:] = phi.domain.cellsize.z
  FaceValue(phi.domain,
    phi.value[2:end,2:end-1,2:end-1].*phi.value[1:end-1,2:end-1,2:end-1].*(dx[2:end]+dx[1:end-1])./(dx[2:end].*phi.value[1:end-1,2:end-1,2:end-1]+dx[1:end-1].*phi.value[2:end,2:end-1,2:end-1]),
    phi.value[2:end-1,2:end,2:end-1].*phi.value[2:end-1,1:end-1,2:end-1].*(dy[:,1:end-1]+dy[:,2:end])./(dy[:,2:end].*phi.value[2:end-1,1:end-1,2:end-1]+dy[:,1:end-1].*phi.value[2:end-1,2:end,2:end-1]),
    phi.value[2:end-1,2:end-1,2:end].*phi.value[2:end-1,2:end-1,1:end-1].*(dz[:,:,1:end-1]+dz[:,:,2:end])./(dz[:,:,2:end].*phi.value[2:end-1,2:end-1,1:end-1]+dz[:,:,1:end-1].*phi.value[2:end-1,2:end-1,2:end]))
end
end

# ================== Upwind averaging scheme ==================
function upwindMean(phi::CellValue, u::FaceValue)
# calculates the average values of a cell variable. The output is a
# face variable
d=phi.domain.dimension
phi_tmp = Base.copy(phi.value)
if d==1 || d==1.5
  ux=u.xvalue
  # assign the value of the left boundary to the left ghost cell
  phi_tmp[1] = 0.5*(phi.value[1]+phi.value[2])
  # assign the value of the right boundary to the right ghost cell
  phi_tmp[end] = 0.5*(phi.value[end]+phi.value[end-1])
  FaceValue(phi.domain,
    (ux.>0.0).*phi_tmp[1:end-1]+(ux.<0.0).*phi_tmp[2:end]+
    0.5*(ux.==0).*(phi.value[1:end-1]+phi.value[2:end]),
    [1.0],
    [1.0])
elseif d==2 || d==2.5 || d==2.8
  ux = u.xvalue
  uy = u.yvalue
  # assign the value of the left boundary to the left ghost cells
  phi_tmp[1,:] = 0.5*(phi.value[1,:]+phi.value[2,:])
  # assign the value of the right boundary to the right ghost cells
  phi_tmp[end,:] = 0.5*(phi.value[end,:]+phi.value[end-1,:])
  # assign the value of the bottom boundary to the bottom ghost cells
  phi_tmp[:,1] = 0.5*(phi.value[:,1]+phi.value[:,2])
  # assign the value of the top boundary to the top ghost cells
  phi_tmp[:,end] = 0.5*(phi.value[:,end]+phi.value[:,end-1])
  FaceValue(phi.domain,
    (ux.>0.0).*phi_tmp[1:end-1,2:end-1]+
    (ux.<0.0).*phi_tmp[2:end,2:end-1]+
    0.5*(ux.==0.0).*(phi.value[1:end-1,2:end-1]+phi.value[2:end,2:end-1]),
    (uy.>0.0).*phi_tmp[2:end-1,1:end-1]+
    (uy.<0.0).*phi_tmp[2:end-1,2:end]+
    0.5*(uy.==0.0).*(phi.value[2:end-1,1:end-1]+phi.value[2:end-1,2:end]),
    [1.0])
elseif d==3 || d==3.2
  ux = u.xvalue
  uy = u.yvalue
  uz = u.zvalue
  # assign the value of the left boundary to the left ghost cells
  phi_tmp[1,:,:] = 0.5*(phi.value[1,:,:]+phi.value[2,:,:])
  # assign the value of the right boundary to the right ghost cells
  phi_tmp[end,:,:] = 0.5*(phi.value[end,:,:]+phi.value[end-1,:,:])
  # assign the value of the bottom boundary to the bottom ghost cells
  phi_tmp[:,1,:] = 0.5*(phi.value[:,1,:]+phi.value[:,2,:])
  # assign the value of the top boundary to the top ghost cells
  phi_tmp[:,end,:] = 0.5*(phi.value[:,end,:]+phi.value[:,end-1,:])
  # assign the value of the back boundary to the back ghost cells
  phi_tmp[:,:,1] = 0.5*(phi.value[:,:,1]+phi.value[:,:,2])
  # assign the value of the front boundary to the front ghost cells
  phi_tmp[:,:,end] = 0.5*(phi.value[:,:,end]+phi.value[:,:,end-1])
  FaceValue(phi.domain,
    (ux.>0.0).*phi_tmp[1:end-1,2:end-1,2:end-1]+
    (ux.<0.0).*phi_tmp[2:end,2:end-1,2:end-1]+
    0.5*(ux.==0.0).*(phi.value[1:end-1,2:end-1,2:end-1]+phi.value[2:end,2:end-1,2:end-1]),
    (uy.>0.0).*phi_tmp[2:end-1,1:end-1,2:end-1]+
    (uy.<0.0).*phi_tmp[2:end-1,2:end,2:end-1]+
    0.5*(uy.==0.0).*(phi.value[2:end-1,1:end-1,2:end-1]+phi.value[2:end-1,2:end,2:end-1]),
    (uz.>0.0).*phi_tmp[2:end-1,2:end-1,1:end-1]+
    (uz.<0.0).*phi_tmp[2:end-1,2:end-1,2:end]+
    0.5*(uz.==0.0).*(phi.value[2:end-1,2:end-1,1:end-1]+phi.value[2:end-1,2:end-1,2:end]))
end
end

# ================== TVD averaging scheme ==================
function tvdMean(phi::CellValue, u::FaceValue, FL::Function)
# u is a face variable
# phi is a cell variable

# a function to avoid division by zero
eps1 = 1.0e-20
fsign(phi_in) = (abs.(phi_in).>=eps1).*phi_in+eps1*(phi_in.==0.0)+eps1*(abs.(phi_in).<eps1).*sign.(phi_in)

d=phi.domain.dimension

if d == 1 || d == 1.5
  # extract data from the mesh structure
  Nx = u.domain.dims[1]
  dx = 0.5*(u.domain.cellsize.x[1:end-1]+u.domain.cellsize.x[2:end])
  phi_p = zeros(Float64, Nx+1)
  phi_m = zeros(Float64, Nx+1)

  # extract the velocity data
  ux = u.xvalue

  # calculate the upstream to downstream gradient ratios for u>0 (+ ratio)
  dphi_p = (phi.value[2:Nx+2]-phi.value[1:Nx+1])./dx
  rp = dphi_p[1:end-1]./fsign(dphi_p[2:end])
  phi_p[2:Nx+1] = phi.value[2:Nx+1]+0.5*FL.(rp).*(phi.value[3:Nx+2]-phi.value[2:Nx+1])
  phi_p[1] = (phi.value[1]+phi.value[2])/2.0 # left boundary

  # calculate the upstream to downstream gradient ratios for u<0 (- ratio)
  rm = dphi_p[2:end]./fsign(dphi_p[1:end-1])
  phi_m[1:Nx] = phi.value[2:Nx+1]+0.5*FL.(rm).*(phi.value[1:Nx]-phi.value[2:Nx+1])
  phi_m[Nx+1] = (phi.value[end]+phi.value[end-1])/2.0 # right boundary

  FaceValue(phi.domain,
      (ux.>0.0).*phi_p+(ux.<0.0).*phi_m+
      0.5*(ux.==0).*(phi.value[1:end-1]+phi.value[2:end]),
      [1.0],
      [1.0])
elseif d==2 || d==2.5 || d==2.8
  # extract data from the mesh structure
  Nx = u.domain.dims[1]
  Ny = u.domain.dims[2]
  dx=0.5*(u.domain.cellsize.x[1:end-1]+u.domain.cellsize.x[2:end])
  dy=zeros( 1, Ny+1)
  dy[:]=0.5*(u.domain.cellsize.y[1:end-1]+u.domain.cellsize.y[2:end])
  phi_p = zeros(Float64, Nx+1)
  phi_m = zeros(Float64, Nx+1)
  phiX_p = zeros(Float64, Nx+1, Ny)
  phiX_m = zeros(Float64, Nx+1,Ny)
  phiY_p = zeros(Float64, Nx,Ny+1)
  phiY_m = zeros(Float64, Nx,Ny+1)

  # extract the velocity data
  ux = u.xvalue
  uy = u.yvalue

  # calculate the upstream to downstream gradient ratios for u>0 (+ ratio)
  # x direction
  dphiX_p = (phi.value[2:Nx+2, 2:Ny+1]-phi.value[1:Nx+1, 2:Ny+1])./dx
  rX_p = dphiX_p[1:end-1,:]./fsign(dphiX_p[2:end,:])
  phiX_p[2:Nx+1,:] = phi.value[2:Nx+1, 2:Ny+1]+0.5*FL.(rX_p).*
  (phi.value[3:Nx+2,2:Ny+1]-phi.value[2:Nx+1, 2:Ny+1])
  phiX_p[1, :] = (phi.value[1, 2:Ny+1]+phi.value[2, 2:Ny+1])/2.0  # left boundary
  # y direction
  dphiY_p = (phi.value[2:Nx+1, 2:Ny+2]-phi.value[2:Nx+1, 1:Ny+1])./dy
  rY_p = dphiY_p[:,1:end-1]./fsign(dphiY_p[:,2:end])
  phiY_p[:,2:Ny+1] = phi.value[2:Nx+1, 2:Ny+1]+0.5*FL.(rY_p).*
      (phi.value[2:Nx+1,3:Ny+2]-phi.value[2:Nx+1, 2:Ny+1])
  phiY_p[:,1] = (phi.value[2:Nx+1,1]+phi.value[2:Nx+1,2])/2.0  # Bottom boundary

  # calculate the upstream to downstream gradient ratios for u<0 (- ratio)
  # x direction
  rX_m = dphiX_p[2:end,:]./fsign(dphiX_p[1:end-1,:])
  phiX_m[1:Nx,:] = phi.value[2:Nx+1, 2:Ny+1]+0.5*FL.(rX_m).*
      (phi.value[1:Nx, 2:Ny+1]-phi.value[2:Nx+1, 2:Ny+1])
  phiX_m[Nx+1,:] = (phi.value[end, 2:Ny+1]+phi.value[end-1, 2:Ny+1])/2.0  # right boundary
  # y direction
  rY_m = dphiY_p[:,2:end]./fsign(dphiY_p[:,1:end-1])
  phiY_m[:,1:Ny] = phi.value[2:Nx+1, 2:Ny+1]+0.5*FL.(rY_m).*
      (phi.value[2:Nx+1, 1:Ny]-phi.value[2:Nx+1, 2:Ny+1])
  phiY_m[:, Ny+1] = (phi.value[2:Nx+1, end]+phi.value[2:Nx+1, end-1])/2.0  # top boundary

  FaceValue(phi.domain,
      (ux.>0.0).*phiX_p+(ux.<0.0).*phiX_m+
         0.5*(ux.==0.0).*(phi.value[1:Nx+1,2:Ny+1]+phi.value[2:Nx+2,2:Ny+1]),
      (uy.>0.0).*phiY_p+(uy.<0.0).*phiY_m+
         0.5*(uy.==0.0).*(phi.value[2:Nx+1,1:Ny+1]+phi.value[2:Nx+1,2:Ny+2]),
      [1.0])

elseif d==3 || d==3.2
  # extract data from the mesh structure
  Nx = u.domain.dims[1]
  Ny = u.domain.dims[2]
  Nz = u.domain.dims[3]
  dx=0.5*(u.domain.cellsize.x[1:end-1]+u.domain.cellsize.x[2:end])
  dy=zeros( 1, Ny+1)
  dy[:]=0.5*(u.domain.cellsize.y[1:end-1]+u.domain.cellsize.y[2:end])
  dz=zeros( 1, 1, Nz+1)
  dz[:]=0.5*(u.domain.cellsize.z[1:end-1]+u.domain.cellsize.z[2:end])
  # extract the velocity data
  ux = u.xvalue
  uy = u.yvalue
  uz = u.zvalue

  # define the tvd face vectors
  phiX_p = zeros(Float64, Nx+1,Ny,Nz)
  phiX_m = zeros(Float64, Nx+1,Ny,Nz)
  phiY_p = zeros(Float64, Nx,Ny+1,Nz)
  phiY_m = zeros(Float64, Nx,Ny+1,Nz)
  phiZ_p = zeros(Float64, Nx,Ny,Nz+1)
  phiZ_m = zeros(Float64, Nx,Ny,Nz+1)

  # calculate the upstream to downstream gradient ratios for u>0 (+ ratio)
  # x direction
  dphiX_p = (phi.value[2:Nx+2, 2:Ny+1, 2:Nz+1]-phi.value[1:Nx+1, 2:Ny+1, 2:Nz+1])./dx
  rX_p = dphiX_p[1:end-1,:,:]./fsign(dphiX_p[2:end,:,:])
  phiX_p[2:Nx+1,:,:] = phi.value[2:Nx+1, 2:Ny+1, 2:Nz+1]+0.5*FL.(rX_p).*
      (phi.value[3:Nx+2,2:Ny+1,2:Nz+1]-phi.value[2:Nx+1,2:Ny+1,2:Nz+1])
  phiX_p[1,:,:] = (phi.value[1,2:Ny+1,2:Nz+1]+phi.value[2,2:Ny+1,2:Nz+1])/2.0  # left boundary
  # y direction
  dphiY_p = (phi.value[2:Nx+1, 2:Ny+2, 2:Nz+1]-phi.value[2:Nx+1, 1:Ny+1, 2:Nz+1])./dy
  rY_p = dphiY_p[:,1:end-1,:]./fsign(dphiY_p[:,2:end,:])
  phiY_p[:,2:Ny+1,:] = phi.value[2:Nx+1, 2:Ny+1, 2:Nz+1]+0.5*FL.(rY_p).*
      (phi.value[2:Nx+1,3:Ny+2,2:Nz+1]-phi.value[2:Nx+1, 2:Ny+1,2:Nz+1])
  phiY_p[:,1,:] = (phi.value[2:Nx+1,1,2:Nz+1]+phi.value[2:Nx+1,2,2:Nz+1])/2.0  # Bottom boundary
  # z direction
  dphiZ_p = (phi.value[2:Nx+1, 2:Ny+1, 2:Nz+2]-phi.value[2:Nx+1, 2:Ny+1, 1:Nz+1])./dz
  rZ_p = dphiZ_p[:,:,1:end-1]./fsign(dphiZ_p[:,:,2:end])
  phiZ_p[:,:,2:Nz+1] = phi.value[2:Nx+1, 2:Ny+1, 2:Nz+1]+0.5*FL.(rZ_p).*
      (phi.value[2:Nx+1,2:Ny+1,3:Nz+2]-phi.value[2:Nx+1,2:Ny+1,2:Nz+1])
  phiZ_p[:,:,1] = (phi.value[2:Nx+1,2:Ny+1,1]+phi.value[2:Nx+1,2:Ny+1,2])/2.0  # Back boundary

  # calculate the upstream to downstream gradient ratios for u<0 (- ratio)
  # x direction
  rX_m = dphiX_p[2:end,:,:]./fsign(dphiX_p[1:end-1,:,:])
  phiX_m[1:Nx,:,:] = phi.value[2:Nx+1, 2:Ny+1, 2:Nz+1]+0.5*FL.(rX_m).*
      (phi.value[1:Nx, 2:Ny+1, 2:Nz+1]-phi.value[2:Nx+1, 2:Ny+1, 2:Nz+1])
  phiX_m[Nx+1,:,:] = (phi.value[end,2:Ny+1,2:Nz+1]+phi.value[end-1,2:Ny+1,2:Nz+1])/2.0  # right boundary
  # y direction
  rY_m = dphiY_p[:,2:end,:]./fsign(dphiY_p[:,1:end-1,:])
  phiY_m[:,1:Ny,:] = phi.value[2:Nx+1,2:Ny+1,2:Nz+1]+0.5*FL.(rY_m).*
      (phi.value[2:Nx+1,1:Ny,2:Nz+1]-phi.value[2:Nx+1,2:Ny+1,2:Nz+1])
  phiY_m[:,Ny+1,:] = (phi.value[2:Nx+1, end,2:Nz+1]+phi.value[2:Nx+1, end-1,2:Nz+1])/2.0  # top boundary
  # z direction
  rZ_m = dphiZ_p[:,:,2:end]./fsign(dphiZ_p[:,:,1:end-1])
  phiZ_m[:,:,1:Nz] = phi.value[2:Nx+1,2:Ny+1,2:Nz+1]+0.5*FL.(rZ_m).*
      (phi.value[2:Nx+1,2:Ny+1,1:Nz]-phi.value[2:Nx+1,2:Ny+1,2:Nz+1])
  phiZ_m[:,:,Nz+1] = (phi.value[2:Nx+1,2:Ny+1,end]+phi.value[2:Nx+1,2:Ny+1,end-1])/2.0  # front boundary

  FaceValue(phi.domain,
      (ux.>0.0).*phiX_p+(ux.<0.0).*phiX_m+
         0.5*(ux.==0.0).*(phi.value[1:Nx+1,2:Ny+1,2:Nz+1]+phi.value[2:Nx+2,2:Ny+1,2:Nz+1]),
      (uy.>0.0).*phiY_p+(uy.<0).*phiY_m+
         0.5*(uy.==0.0).*(phi.value[2:Nx+1,1:Ny+1,2:Nz+1]+phi.value[2:Nx+1,2:Ny+2,2:Nz+1]),
      (uz.>0.0).*phiZ_p+(uz.<0).*phiZ_m+
         0.5*(uz.==0.0).*(phi.value[2:Nx+1,2:Ny+1,1:Nz+1]+phi.value[2:Nx+1,2:Ny+1,2:Nz+2]))
end

end

max_stars_repo_name: simulkade/JFVM

max_forks_repo_name: simulkade/JFVM


text:
module aoc_13

using LinearAlgebra

export parse_input, problem13, time_until_next_bus, find_next_sync
export solve_p2

struct problem13
    timestamp::Int
    bus_ids::Vector{Int}
    positions::Vector{Int}
end

function parse_input(filename, keepx=false)
    lines = readlines(filename)
    timestamp = parse(Int,lines[1])
    strings = split(lines[2],",")
#    if keepx
    bus_ids = [parse(Int,replace(s, "x" => "0")) for s in strings]
    positions = findall(bus_ids .!= 0 )
    bus_ids = bus_ids[ positions ]

#    else
#        bus_ids = [parse(Int,s) for s in strings if s != "x"]
#    end
    return problem13(timestamp, bus_ids, positions)
end

time_until_next_bus(p) = p.bus_ids .- mod.(p.timestamp, p.bus_ids)

function get_product(p)
    time_until = time_until_next_bus(p)
    min_index = findmin(time_until)
    return p.bus_ids[ min_index[2] ] * min_index[1]
end

function find_next_sync(bus_ids, positions, start=0)
    i = 0
    t = start
    product = prod(bus_ids[1:end-1])
    while true
        if all( mod(t + positions[end] - positions[1],bus_ids[end]) == 0 )
            println(t)
            break
        end
        t += product
    end
    return t
end

function solve_p2(p::problem13)
    perm = sortperm(p.bus_ids, rev=true)
    p_sorted = problem13( 0 , p.bus_ids[perm], p.positions[perm] )
    prev_sync = 0
    for l = 2:length(perm)
        prev_sync = find_next_sync(p_sorted.bus_ids[1:l], p_sorted.positions[1:l], prev_sync)
    end
    return prev_sync - p_sorted.positions[1] + p.positions[1]
end

end

ext: jl


text:
#__precompile__()

module Ravana

# Module Initializer
function __init__()
    init_raft()
end

# Source Files
include("Utils.jl")
include("CoreOps.jl")
include("Raft.jl")
include("ProtocolServer.jl")
include("Server.jl")

# Exported
export raft_client, init_cluster

# Deprecated

end # module Ravana


text:
const stext = ["one", "two", "three", "four", "five", "six", "seven", "eight", "nine"]
const teentext = ["eleven", "twelve", "thirteen", "fourteen", "fifteen", "sixteen",
    "seventeen", "eighteen", "nineteen"]
const tenstext = ["ten", "twenty", "thirty", "forty", "fifty", "sixty", "seventy",
    "eighty", "ninety"]
const ordstext = ["million", "billion", "trillion", "quadrillion", "quintillion", "sextillion",
    "septillion", "octillion", "nonillion", "decillion", "undecillion", "duodecillion",
    "tredecillion", "quattuordecillion", "quindecillion", "sexdecillion", "septendecillion",
    "octodecillion", "novemdecillion", "vigintillion"]

function normalize_digits!(a::Array{T,1}) where T<:Integer
    while 0 < length(a) && a[end] == 0
        pop!(a)
    end
    length(a)
end

function digits2text!(d::Array{T,1}, use_short_scale=true) where T<:Integer
    ndig = normalize_digits!(d)
    0 < ndig || return ""
    if ndig < 7
        s = ""
        if 3 < ndig
            t = digits2text!(d[1:3])
            s = digits2text!(d[4:end]) * " thousand"
            0 < length(t) || return s
            return occursin(t, "and") ? s * " " * t : s * " and " * t
        end
        if ndig == 3
            s *= stext[pop!(d)] * " hundred"
            ndig = normalize_digits!(d)
            0 < ndig || return s
            s *= " and "
        end
        1 < ndig || return s*stext[pop!(d)]
        j, i = d
        j ≠ 0 || return s*tenstext[i]
        i ≠ 1 || return s*teentext[j]
        return s*tenstext[i] * "-" * stext[j]
    end
    s = digits2text!(d[1:6])
    d = d[7:end]
    dgrp = use_short_scale ? 3 : 6
    ord = 0
    while dgrp < length(d)
        ord += 1
        t = digits2text!(d[1:dgrp])
        d = d[(dgrp+1):end]
        0 < length(t) || continue
        t = t * " " * ordstext[ord]
        s = length(s) == 0 ? t : t * " " * s
    end
    ord += 1
    t = digits2text!(d) * " " * ordstext[ord]
    0 < length(s) || return t
    t * " " * s
end

function num2text(n::T, use_short_scale=true) where T<:Integer
    -1 < n || return "minus "*num2text(-n, use_short_scale)
    0 < n || return "zero"
    toobig = use_short_scale ? big(10)^66 : big(10)^126
    n < toobig || return "too big to say"
    digits2text!(digits(n, base=10), use_short_scale)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
import Base: convert
import Base: ==

const BeliefArray{T} = Union{Array{T,2}, Adjoint{T, Array{T,2}} }

"""
$(TYPEDEF)
"""
struct SingleThreaded
end
"""
$(TYPEDEF)
"""
struct MultiThreaded
end

"""
$(TYPEDEF)

Solver parameters for the DistributedFactoGraph.

Dev Notes
- FIXME change to using kwargs from Parameters.jl
- TODO remove NothingUnion
- TODO Upgrade to common @kwargs struct approach
"""
mutable struct SolverParams <: DFG.AbstractParams
  dimID::Int
  # TODO remove NothingUnion
  registeredModuleFunctions::NothingUnion{Dict{Symbol, Function}} # remove from
  reference::NothingUnion{Dict{Symbol, Tuple{Symbol, Vector{Float64}}}}
  stateless::Bool
  qfl::Int # Quasi fixed length
  isfixedlag::Bool # true when adhering to qfl window size for solves
  limitfixeddown::Bool # if true, then fixed lag will also not update marginalized during down (default false)
  # new functions
  incremental::Bool
  useMsgLikelihoods::Bool
  upsolve::Bool
  downsolve::Bool
  drawtree::Bool
  drawCSMIters::Bool
  showtree::Bool
  drawtreerate::Float64
  dbg::Bool
  async::Bool
  limititers::Int
  N::Int
  multiproc::Bool
  logpath::String
  graphinit::Bool
  treeinit::Bool # still experimental with known errors
  limittreeinit_iters::Int
  algorithms::Vector{Symbol} # list of algorithms to run [:default] is mmisam
  spreadNH::Float64 # experimental, entropy spread adjustment used for both null hypo cases.
  maxincidence::Int # maximum incidence to a variable in an effort to enhance sparsity
  alwaysFreshMeasurements::Bool
  devParams::Dict{Symbol,String}
  #
end

SolverParams(;dimID::Int=0,
              registeredModuleFunctions=nothing,
              reference=nothing,
              stateless::Bool=false,
              qfl::Int=99999999999,
              isfixedlag::Bool=false,
              limitfixeddown::Bool=false,
              incremental::Bool=true,
              useMsgLikelihoods::Bool=true,
              upsolve::Bool=true,
              downsolve::Bool=true,
              drawtree::Bool=false,
              drawCSMIters::Bool=true,
              showtree::Bool=false,
              drawtreerate::Float64=0.5,
              dbg::Bool=false,
              async::Bool=false,
              limititers::Int=500,
              N::Int=100,
              multiproc::Bool=1 < nprocs(),
              logpath::String="/tmp/caesar/$(now())",
              graphinit::Bool=true,
              treeinit::Bool=false,
              limittreeinit_iters::Int=10,
              algorithms::Vector{Symbol}=[:default],
              spreadNH::Float64=3.0,
              maxincidence::Int=500,
              alwaysFreshMeasurements::Bool=true,
              devParams::Dict{Symbol,String}=Dict{Symbol,String}()
            ) = SolverParams( dimID,
                              registeredModuleFunctions,
                              reference,
                              stateless,
                              qfl,
                              isfixedlag,
                              limitfixeddown,
                              incremental,
                              useMsgLikelihoods,
                              upsolve,
                              downsolve,
                              drawtree,
                              drawCSMIters,
                              showtree,
                              drawtreerate,
                              dbg,
                              async,
                              limititers,
                              N,
                              multiproc,
                              logpath,
                              graphinit,
                              treeinit,
                              limittreeinit_iters,
                              algorithms,
                              spreadNH,
                              maxincidence,
                              alwaysFreshMeasurements,
                              devParams )
#

"""
$(TYPEDEF)

Notes
- type-stable `usercache`, see https://github.com/JuliaRobotics/IncrementalInference.jl/issues/783#issuecomment-665080114

DevNotes
- TODO why not just a NamedTuple? Perhaps part of #467
- TODO better consolidate with CCW, CPT, CalcFactor
- TODO standardize -- #927, #1025, #784, #692, #640
- TODO make immutable #825
"""
mutable struct FactorMetadata{FV<:AbstractVector{<:DFGVariable},
                              VL<:AbstractVector{Symbol},
                              AR<:AbstractVector{<:AbstractArray},
                              CD}
  # full list of Vector{DFGVariable} connected to the factor
  fullvariables::FV # Vector{<:DFGVariable}
  #TODO full variable can perhaps replace this
  variablelist::VL # Vector{Symbol}
  # TODO consolidate, same as ARR used in CCW,
  arrRef::AR # Vector{Matrix{Float64}}
  # label of which variable is being solved for
  solvefor::Symbol
  # for type specific user data, see (? #784)
  cachedata::CD
end

"""
$(TYPEDEF)

DevNotes
- FIXME rm ::Union on `.activehypo`
- TODO consider renaming `.p` to `.decisionDims`
  - TODO `.decisionDims::DD where DD <: Union{<:AbstractVector{Int},Colon}` -- ensure type stability
- TODO consolidate with CCW, FMd, CalcFactor
- TODO figure out if we want static parameter THRID
- TODO make static params {XDIM, ZDIM, P}
- TODO make immutable
"""
mutable struct ConvPerThread{T}
  thrid_::Int
  # the actual particle being solved at this moment
  particleidx::Int
  # additional/optional data passed to user function
  factormetadata::FactorMetadata
  # subsection indices to select which params should be used for this hypothesis evaluation
  activehypo::Union{UnitRange{Int},Vector{Int}}
  # Select which decision variables to include in a particular optimization run
  p::Vector{Int}
  # slight numerical perturbation for degenerate solver cases such as division by zero
  perturb::Vector{Float64}
  # working memory location for optimization routines on target decision variables
  X::Array{Float64,2}
  # working memory to store residual for optimization routines
  res::T # was Vector{Float64}
end

"""
$(TYPEDEF)

DevNotes
- FIXME remove inner constructor
"""
mutable struct CommonConvWrapper{T<:FunctorInferenceType} <: FactorOperationalMemory
  ### Values consistent across all threads during approx convolution
  usrfnc!::T # user factor / function
  # general setup
  xDim::Int
  zDim::Int
  # special case settings
  specialzDim::Bool # is there a special zDim requirement -- defined by user
  partial::Bool # is this a partial constraint -- defined by user
  # multi hypothesis settings
  hypotheses::Union{Nothing, Distributions.Categorical} # categorical to select which hypothesis is being considered during convolution operation
  certainhypo::Union{Nothing, Vector{Int}}
  nullhypo::Float64
  # values specific to one complete convolution operation
  params::Vector{Array{Float64,2}} # parameters passed to each hypothesis evaluation event on user function
  varidx::Int # which index is being solved for in params?
  measurement::Tuple # user defined measurement values for each approxConv operation
  threadmodel::Union{Type{SingleThreaded}, Type{MultiThreaded}}
  ### particular convolution computation values per particle idx (varies by thread)
  cpt::Vector{ConvPerThread}

  # remove
  CommonConvWrapper{T}() where {T<:FunctorInferenceType} = new{T}()
end

function ConvPerThread( X::Array{Float64,2},
                        zDim::Int,
                        factormetadata::FactorMetadata;
                        particleidx::Int=1,
                        activehypo= 1:length(params),
                        p=collect(1:size(X,1)),
                        perturb=zeros(zDim),
                        # Y=zeros(zDim), #zeros(size(X,1)),
                        res=zeros(zDim),
                        thrid_ = 0  )
  #
  # cpt = ConvPerThread()
  return ConvPerThread( thrid_,
                        particleidx,
                        factormetadata,
                        activehypo,
                        p,
                        perturb,
                        X,
                        res)
end

function CommonConvWrapper( fnc::T,
                            X::Array{Float64,2},
                            zDim::Int,
                            params::Vector{Array{Float64,2}},
                            factormetadata::FactorMetadata;
                            specialzDim::Bool=false,
                            partial::Bool=false,
                            hypotheses=nothing,
                            certainhypo=nothing,
                            activehypo= 1:length(params),
                            nullhypo::Real=0,
                            varidx::Int=1,
                            measurement::Tuple=(zeros(0,1),),
                            particleidx::Int=1,
                            p=collect(1:size(X,1)),
                            perturb=zeros(zDim),
                            # Y=zeros(zDim), # zeros(size(X,1)),
                            xDim=size(X,1),
                            res=zeros(zDim),
                            threadmodel=MultiThreaded  ) where {T<:FunctorInferenceType}
  #
  ccw = CommonConvWrapper{T}()

  ccw.usrfnc! = fnc
  ccw.xDim = xDim
  ccw.zDim = zDim
  ccw.specialzDim = specialzDim
  ccw.partial = partial
  ccw.hypotheses = hypotheses
  ccw.certainhypo=certainhypo
  ccw.nullhypo=nullhypo
  ccw.params = params
  ccw.varidx = varidx
  ccw.measurement = measurement
  ccw.threadmodel = threadmodel

  # thread specific elements
  ccw.cpt = Vector{ConvPerThread}(undef, Threads.nthreads())
  for i in 1:Threads.nthreads()
    ccw.cpt[i] = ConvPerThread(X, zDim,
                    factormetadata,
                    particleidx=particleidx,
                    activehypo=activehypo,
                    p=p,
                    perturb=perturb,
                    # Y=Y,
                    res=res )
  end

  return ccw
end

#


text:
for (d, lattype) in enumerate((
        MonteCarlo.Chain,
        MonteCarlo.SquareLattice,
        MonteCarlo.CubicLattice
    ))
    constructor = d < 3 ? L -> lattype(L) : L -> lattype(d, L)

    @testset "$(d)D $lattype" begin
        for L in (3, 4)
            l = constructor(L)
            @test length(l) == L^d

            bonds = collect(MonteCarlo.neighbors(l, Val(true)))
            @test length(bonds) == 2*d*L^d
            @test allunique(bonds)
            for i in 0:length(l)-1
                # same source
                @test all(bonds[2d*i + 1][1] == x[1] for x in bonds[2d*i .+ (2:2d)])
                # different target
                @test allunique(x[2] for x in bonds[2d*i .+ (1:2d)])
            end

            reduced_bonds = collect(MonteCarlo.neighbors(l, Val(false)))
            @test length(reduced_bonds) == d*L^d
            # If directed is false, only one of (i, j) and (j, i)
            # should be kept.
            mirrored_bonds = [[trg, src] for (src, trg) in reduced_bonds]
            all_bonds = vcat(reduced_bonds, mirrored_bonds)
            @test sort(all_bonds) == sort(bonds)
        end
    end
end

@testset "2D Honeycomb" begin
    for L in (2, 3)
        l = MonteCarlo.HoneycombLattice(L)
        @test length(l) == (2L)^2

        bonds = collect(MonteCarlo.neighbors(l, Val(true)))
        @test length(bonds) == 3 * (2L)^2
        @test allunique(bonds)
        for i in 0:length(l)-1
            # same source
            @test all(bonds[3i+1][1] == x[1] for x in bonds[3i .+ (2:3)])
            # different target
            @test allunique(x[2] for x in bonds[3i .+ (1:3)])
        end

        reduced_bonds = collect(MonteCarlo.neighbors(l, Val(false)))
        @test length(reduced_bonds) == round(Int64, 1.5 * (2L)^2)
        # If directed is false, only one of (i, j) and (j, i)
        # should be kept.
        mirrored_bonds = [[trg, src] for (src, trg) in reduced_bonds]
        all_bonds = vcat(reduced_bonds, mirrored_bonds)
        @test sort(all_bonds) == sort(bonds)

        # For L < 3 NNN wrap around the lattice
        L < 3 && continue
        l = MonteCarlo.HoneycombLattice(L, true)
        @test size(l.NNNs) == (6, length(l))
        for i in 1:length(l)
            @test allunique(l.NNNs[:, i])
        end
    end
end


text:
# Copyright 2015, 2016, 2017, 2018 Martin Holters
# See accompanying license file.

export Circuit, add!, connect!, disconnect!, @circuit

import Base: delete!

const Net = Vector{Tuple{Symbol,Symbol}} # pairs of element designator and pin name

struct Circuit
    elements::OrderedDict{Symbol, Element}
    nets :: Vector{Net}
    net_names :: Dict{Symbol, Net}
    Circuit() = new(OrderedDict{Symbol, Element}(), [], Dict{Symbol, Net}())
end

elements(c::Circuit) = values(c.elements)

for n in [:nb; :nx; :nq; :nu; :nl; :ny; :nn]
    @eval ($n)(c::Circuit) = Compat.reduce(+, init=0, $n(elem) for elem in elements(c))
end

for mat in [:mv; :mi; :mx; :mxd; :mq; :mu; :pv; :pi; :px; :pxd; :pq]
    # blockdiag() does not work, so include an empty matrix of desired type in
    # case elements(c) is empty
    # as blockdiag for unknown number of arguments cannot be inferred properly,
    # add type-assertion
    @eval ($mat)(c::Circuit) =
         blockdiag(spzeros(Rational{BigInt}, 0, 0),
                 (convert(SparseMatrixCSC{Rational{BigInt}}, elem.$mat)
                  for elem in elements(c))...
                )::SparseMatrixCSC{Rational{BigInt},Int}
end

u0(c::Circuit) = vcat((elem.u0 for elem in elements(c))...)

function incidence(c::Circuit)
    i = sizehint!(Int[], 2nb(c))
    j = sizehint!(Int[], 2nb(c))
    v = sizehint!(Int[], 2nb(c))
    for (row, pins) in enumerate(c.nets), (elemname, pinname) in pins
        elem = c.elements[elemname]
        offset = branch_offset(c, elem)
        bps = elem.pins[pinname]
        for (branch, polarity) in bps
            push!(i, row)
            push!(j, offset + branch)
            push!(v, polarity)
        end
    end
    # ensure zeros due to short-circuited branches are removed with dropzeros!
    dropzeros!(sparse(i, j, v, length(c.nets), nb(c)))
end

function nonlinear_eq(c::Circuit, elem_idxs=1:length(elements(c)))
    # construct a block expression containing all element's expressions after
    # offsetting their indexes into q, J and res

    row_offset = 0
    col_offset = 0
    nl_expr = Expr(:block)
    for elem in collect(elements(c))[elem_idxs]
        index_offsets = Dict( :q => (col_offset,),
                              :J => (row_offset, col_offset),
                              :res => (row_offset,) )

        function offset_indexes(expr::Expr)
            ret = Expr(expr.head)
            if expr.head == :ref && haskey(index_offsets, expr.args[1])
                push!(ret.args, expr.args[1])
                offsets = index_offsets[expr.args[1]]
                length(expr.args) == length(offsets) + 1 ||
                    throw(ArgumentError(string(expr.args[1], " must be indexed with exactly ", length(offsets), " index(es)")))
                for i in 1:length(offsets)
                    push!(ret.args,
                          :($(offsets[i]) + $(offset_indexes(expr.args[i+1]))))
                end
            else
                append!(ret.args, offset_indexes.(expr.args))
            end
            ret
        end

        function offset_indexes(s::Symbol)
            haskey(index_offsets, s) && throw(ArgumentError(string(s, " used without indexing expression")))
            s
        end

        offset_indexes(x::Any) = x

        # wrap into a let to keep variables local
        push!(nl_expr.args, :( let; $(offset_indexes(elem.nonlinear_eq)) end))

        row_offset += nn(elem)
        col_offset += nq(elem)
    end
    nl_expr
end

"""
    add!(c::Circuit, elem::Element)

Adds the element `elem` to the circuit `c`, creating and returning a new, unique
reference designator, leaving its pins unconnected.
"""
function add!(c::Circuit, elem::Element)
    for (k, v) in c.elements
        if v == elem
            return k
        end
    end
    designator = gensym()
    add!(c, designator, elem)
    return designator
end

add!(c::Circuit, elems::Element...) = ((add!(c, elem) for elem in elems)...,)

"""
    add!(c::Circuit, designator::Symbol, elem::Element)

Adds the element `elem` to the circuit `c` with the reference designator
`designator`, leaving its pins unconnected. If the circuit already contained
an element named `designator`, it is removed first.
"""
function add!(c::Circuit, designator::Symbol, elem::Element)
    if haskey(c.elements, designator)
        delete!(c, designator)
    end
    for pin in keys(elem.pins)
        push!(c.nets, [(designator, pin)])
    end
    c.elements[designator] = elem
end

"""
    delete!(c::Circuit, designator::Symbol)

Deletes the element named `designator` from the circuit `c` (disconnecting all
its pins).
"""
function delete!(c::Circuit, designator::Symbol)
    for net in c.nets
        filter!(elempin -> elempin[1] != designator, net)
    end
    delete!(c.elements, designator)
end

function branch_offset(c::Circuit, elem::Element)
    offset = 0
    for el in elements(c)
        el == elem && return offset
        offset += nb(el)
    end
    throw(ArgumentError("Element not found in circuit"))
end

function netfor!(c::Circuit, p::Tuple{Symbol,Symbol})
    for net in c.nets
        p ∈ net && return net
    end
    throw(ArgumentError("Unknown pin $p"))
end
netfor!(c::Circuit, p::Tuple{Symbol,Any}) = netfor!(c, (p[1], Symbol(p[2])))

function netfor!(c::Circuit, p::Pin)
    Base.depwarn("pin specification $p is deprecated, use (refdes, pinname) instead", :netfor!)
    element = p[1]
    designator = add!(c, element)
    local pinname
    for (pname, pbps) in element.pins
        if pbps == p[2]
            pinname = pname
            break
        end
    end
    return netfor!(c, (designator, pinname))
end

function netfor!(c::Circuit, name::Symbol)
    haskey(c.net_names, name) || push!(c.nets, get!(c.net_names, name, []))
    c.net_names[name]
end

"""
    connect!(c::Circuit, pins::Union{Symbol,Tuple{Symbol,Any}}...)

Connects the given pins (or named nets) to each other in the circuit `c`. Named
nets are given as `Symbol`s, pins are given as `Tuple{Symbols,Any}`s, where the
first entry is the reference designator of an element in `c`, and the second
entry is the pin name. For convenience, the latter is automatically converted to
a `Symbol` as needed.

# Example

```jldoctest; output = false, setup = :(using ACME)
circ = Circuit()
add!(circ, :r, resistor(1e3))
add!(circ, :src, voltagesource(5))
connect!(circ, (:src, -), (:r, 2), :gnd) # connect to gnd net

# output

```
"""
function connect!(c::Circuit, pins::Union{Pin,Symbol,Tuple{Symbol,Any}}...)
    nets = unique([netfor!(c, pin) for pin in pins])
    for net in nets[2:end]
        append!(nets[1], net)
        deleteat!(c.nets, findfirst(n -> n == net, c.nets))
        for (name, named_net) in c.net_names
            if named_net == net
                c.net_names[name] = nets[1]
            end
        end
    end
end

function disconnect!(c::Circuit, pin::Tuple{Symbol,Symbol})
    net = netfor!(c, pin)
    filter!(p -> p != pin, net)
    push!(c.nets, [pin])
end
"""
    disconnect!(c::Circuit, p::Tuple{Symbol,Any})

Disconnects the given pin `p` from anything else in the circuit `c`. The pin is
given as a`Tuple{Symbols,Any}`, where the first entry is the reference
designator of an element in `c`, and the second entry is the pin name. For
convenience, the latter is automatically converted to a `Symbol` as needed. Note
that if e.g. three pin `p1`, `p2`, and `p3` are connected then
`disconnect!(c, p1)` will disconnect `p1` from `p2` and `p3`, but leave `p2` and
`p3` connected to each other.
"""
disconnect!(c::Circuit, p::Tuple{Symbol,Any}) = disconnect!(c, (p[1], Symbol(p[2])))

function disconnect!(c::Circuit, pin::Pin)
    Base.depwarn("disconnect!(::Circuit, $p) is deprecated, use (refdes, pinname) to specify the the pin", :disconnect!)
    element = pin[1]
    designator = add!(c, element)
    local pinname
    for (pname, pbps) in element.pins
        if pbps == pin[2]
            pinname = pname
            break
        end
    end
    disconnect!(c, (designator, pinname))
end

function topomat!(incidence::SparseMatrixCSC{T}) where {T<:Integer}
    @assert all(x -> abs(x) == 1, nonzeros(incidence))
    @static if VERSION ≥ v"0.7.0-DEV.4064"
        @assert all(sum(incidence, dims=1) .== 0)
    else
        @assert all(sum(incidence, 1) .== 0)
    end

    t = falses(size(incidence)[2]);

    row = 1;
    for col = 1:size(incidence)[2]
        rows = filter(r -> r ≥ row, findall(!iszero, incidence[:, col]))
        @assert length(rows) ≤ 2

        isempty(rows) && continue
        t[col] = true;

        if rows[1] ≠ row
            incidence[[rows[1], row], :] = incidence[[row, rows[1]], :]
        end
        if length(rows) == 2
            @assert incidence[row, col] + incidence[rows[2], col] == 0
            incidence[rows[2],:] = incidence[rows[2],:] + incidence[row,:]
        end
        if incidence[row, col] < 0
            cols = findall(!iszero, incidence[row, :])
            incidence[row,cols] = -incidence[row,cols]
        end
        rows = findall(==(1), incidence[1:row-1, col])
        incidence[rows, :] .-= incidence[row, :]'
        rows = findall(==(-1), incidence[1:row-1, col])
        incidence[rows, :] .+= incidence[row, :]'
        row += 1
    end

    ti = incidence[1:row-1, :]

    dl = ti[:, (!).(t)]
    tv = spzeros(T, size(dl, 2), size(incidence, 2))
    # findall here works around JuliaLang/julia#27013 (introdcued in 0.7.0-DEV.4983)
    tv[:, findall(t)] = -dl'
    tv[:, findall((!).(t))] = SparseMatrixCSC{T}(I, size(dl, 2), size(dl, 2))

    tv, ti
end

topomat(incidence::SparseMatrixCSC{<:Integer}) = topomat!(copy(incidence))
topomat(c::Circuit) = topomat!(incidence(c))

@doc doc"""
    @circuit begin #= ... =# end

Provides a simple domain-specific language to decribe circuits. The
`begin`/`end` block can hold element definitions of the form
`refdes = elementfunc(params)` and connection specifications of the form
`refdes1[pin1] ⟷ refdes2[pin2]`.

# Example

To create a circuit with a voltage source connected to a resistor:
```jldoctest; output = false, setup = :(using ACME), filter = r"(ACME\.)?Circuit\(.*"s
@circuit begin
    src = voltagesource(5)
    r = resistor(1000)
    src[+] ⟷ r[1]
    src[-] ⟷ r[2]
end

# output

Circuit(...)
```

Alternatively, connection specifications can be given after an element
specification, separated by commas. In that case, the `refdes` may be omitted,
defaulting to the current element.

# Example

```jldoctest; output = false, setup = :(using ACME), filter = r"(ACME\.)?Circuit\(.*"s
@circuit begin
    src = voltagesource(5)
    r = resistor(1000), src[+] ⟷ [1], src[-] ⟷ [2]
end

# output

Circuit(...)
```

Finally, a connection endpoint may simply be of the form `netname`, to connect
to a named net. (Such named nets are created as needed.)

# Example

```jldoctest; output = false, setup = :(using ACME), filter = r"(ACME\.)?Circuit\(.*"s
@circuit begin
    src = voltagesource(5), [-] ⟷ gnd
    r = resistor(1000), [1] ⟷ src[+], [2] ⟷ gnd
end

# output

Circuit(...)
```

If a net or pin specification is not just a single symbol or number, and has to
be put in quotes (e.g. `"in+"`, `"9V"`)

!!! note
    Instead of `⟷` (`\\longleftrightarrow`), one can also use `==`.
""" macro circuit(cdef)
    is_conn_spec(expr::Expr) =
        (expr.head === :call && (expr.args[1] === :(⟷) || expr.args[1] === :(↔) || expr.args[1] === :(==))) ||
        (expr.head === :comparison && all(c -> c === :(==), expr.args[2:2:end]))
    is_conn_spec(::Any) = false

    function elem_spec(expr)
        if !isa(expr, Expr) || expr.head !== :(=)
            error("invalid element specification$locinfo: $(expr)")
        end
        if !isa(expr.args[1], Symbol)
            error("invalid element identifier$locinfo: $(expr.args[1])")
        end
        if isa(expr.args[2], Expr) && expr.args[2].head === :tuple
            if isempty(expr.args[2].args)
                error("invalid element specification$locinfo: $(expr.args[2])")
            end
            elemspec = expr.args[2].args[1]
            conn_exprs = expr.args[2].args[2:end]
        else
            elemspec = expr.args[2]
            conn_exprs = []
        end
        push!(ccode.args, :(add!(circ, $(QuoteNode(expr.args[1])), $(esc(elemspec)))))
        for conn_expr in conn_exprs
            if !is_conn_spec(conn_expr)
                error("invalid connection specification$locinfo: $conn_expr")
            end
            push!(ccode.args, Expr(:call, :connect!, :circ, extractpins(conn_expr, expr.args[1])...))
        end
    end

    function extractpins(expr::Expr, default_element=nothing)
        if expr.head === :call && (expr.args[1] === :(⟷) || expr.args[1] === :(↔) || expr.args[1] === :(==))
            return vcat((extractpins(a, default_element) for a in expr.args[2:end])...)
        elseif expr.head === :comparison && all(c -> c === :(==), expr.args[2:2:end])
            return vcat((extractpins(a, default_element) for a in expr.args[1:2:end])...)
        elseif expr.head === :ref
            return [:(($(QuoteNode(expr.args[1])), $(QuoteNode(expr.args[2]))))]
        elseif expr.head === :vect && length(expr.args) == 1
            if default_element === nothing
                error("missing element$(locinfo): $expr")
            end
            return [:(($(QuoteNode(default_element)), $(QuoteNode(expr.args[1]))))]
        else
            error("invalid pin specification$(locinfo): $expr")
        end
    end

    function extractpins(netname::Symbol, default_element=nothing)
        return [QuoteNode(netname)]
    end

    extractpins(netname::String, default_element=nothing) =
        extractpins(Symbol(netname), default_element)

    if !isa(cdef, Expr) || cdef.head !== :block
        error("@circuit must be followed by a begin/end block")
    end
    ccode = Expr(:block)
    push!(ccode.args, :(circ = Circuit()))
    locinfo = ""
    for expr in cdef.args
        if isa(expr, LineNumberNode)
            locinfo = " at $(expr.file):$(expr.line)"
            continue
        end
        if !isa(expr, Expr)
            error("invalid statement in circuit definition$locinfo: $expr")
        end
        if expr.head === :line
            locinfo = " at $(expr.args[2]):$(expr.args[1])"
        elseif expr.head === :(=)
            elem_spec(expr)
        elseif is_conn_spec(expr)
            push!(ccode.args, Expr(:call, :connect!, :circ, extractpins(expr)...))
        else
            error("invalid statement in circuit definition$locinfo: $expr")
        end
    end

    push!(ccode.args, :(circ))
    return ccode
end


text:
using StochasticDiffEq, DiffEqNoiseProcess, Test, DiffEqDevTools, Random, LinearAlgebra
const μ = 1.01
const σ_const = 0.87

f(u,p,t) = μ * u + μ * u
f1_μ(u,p,t) = μ
f1_μ_analytic(u0,p,t,W) = u0.*exp.((2μ-(σ_const^2)/2)t+σ_const*W)
f2(u,p,t) = μ * u
σ(u,p,t) = σ_const*u
no_noise(u,p,t) = 0.0
f1_no_noise(u,p,t) = μ
f1_no_noise_analytic(u0,p,t,W) = u0.*exp.(2μ*t)

ff1_μ = SplitSDEFunction(f1_μ,f2,σ,analytic=f1_μ_analytic)
prob = SDEProblem(ff1_μ,σ,1/2,(0.0,1.0))
ff1_no_noise = SplitSDEFunction(f1_no_noise,f2,no_noise,analytic=f1_no_noise_analytic)
no_noise_prob = SDEProblem(ff1_no_noise,no_noise,1/2,(0.0,1.0))

sol = solve(prob,IIF1M(),dt=1/10)

prob2 = SDEProblem{false}(f,σ,1/2,(0.0,1.0),noise = NoiseWrapper(sol.W))

sol2 = solve(prob2,EM(),dt=1/10)

sol = solve(no_noise_prob,IIF1M(),dt=1/10)

Random.seed!(100)
dts = (1/2) .^ (7:-1:4) #14->7 good plot
println("IIF scalar")
sim  = test_convergence(dts,prob,IIF1M(),trajectories=Int(1e2))
@test abs(sim.𝒪est[:l2]-0.5) < 0.2 # closer to 1 at this part
sim  = test_convergence(dts,no_noise_prob,IIF1M(),trajectories=Int(2e1))
@test abs(sim.𝒪est[:l2]-1.0) < 0.2 # closer to 1 at this part

dts = (1/2) .^ (7:-1:4) #14->7 good plot
println("IIF no noise scalar")
Random.seed!(100)
sim  = test_convergence(dts,prob,IIF2M(),trajectories=Int(1e2))
@test abs(sim.𝒪est[:l2]-0.5) < 0.2 # closer to 1 at this part
sim  = test_convergence(dts,no_noise_prob,IIF2M(),trajectories=Int(1e1))
@test abs(sim.𝒪est[:l2]-2) < 0.2 # closer to 1 at this part

#=
Random.seed!(200)
sim  = test_convergence(dts,prob,IIF1Mil(),trajectories=Int(2e1))
@test abs(sim.𝒪est[:l2]-1) < 0.3
=#

u0 = rand(2)
A = [-2.0 1.0; 1.0 -2.0]
B = [σ_const 0
    0 σ_const]

function f(du,u,p,t)
  mul!(du,A,u)
  du .+= 1.01u
end
function σ(du,u,p,t)
  mul!(@view(du[:,1]),B,u)
  mul!(@view(du[:,2]),B,u)
end

function f_analytic(u0,p,t,W)
 tmp = (A+1.01I-(B^2))*t + B*sum(W)
 exp(tmp)*u0
end

f1_A(du,u,p,t) = A
function f1_A_analytic(u0,p,t,W)
 tmp = (A+1.01I-(B^2))*t + B*sum(W)
 exp(tmp)*u0
end
f2(du,u,p,t) = du .= μ .* u

ff1_A = SplitSDEFunction(f1_A,f2,σ,analytic=f1_A_analytic)
prob = SDEProblem(ff1_A,σ,u0,(0.0,1.0),noise_rate_prototype=rand(2,2))

f1_no_noise(du,u,p,t) = A
f2(du,u,p,t) = (du .= μ .* u)
function σ22(du,u,p,t)
  du .= 0
end
function f1_no_noise_analytic(u0,p,t,W)
 tmp = (A+1.01I)*t
 exp(tmp)*u0
end
ff1_A = SplitSDEFunction(f1_no_noise,f2,σ22,analytic=f1_no_noise_analytic)
prob_no_noise = SDEProblem(ff1_A,σ22,u0,(0.0,1.0),noise_rate_prototype=rand(2,2))

sol = solve(prob,IIF1M(),dt=1/10)

dts = (1/2) .^ (8:-1:4) #14->7 good plot

Random.seed!(250)
println("IIF")
sim  = test_convergence(dts,prob,IIF1M(),trajectories=Int(5e1))
@test abs(sim.𝒪est[:l2]-0.5) < 0.2

sim  = test_convergence(dts,prob,IIF2M(),trajectories=Int(5e1))
@test abs(sim.𝒪est[:l2]-0.5) < 0.2

println("IIF no noise")
sim  = test_convergence(dts,prob_no_noise,IIF1M(),trajectories=Int(1e1))
@test abs(sim.𝒪est[:l2]-1) < 0.2

sim  = test_convergence(dts,prob_no_noise,IIF2M(),trajectories=Int(1e1))
@test abs(sim.𝒪est[:l2]-2) < 0.1


text:
module ScalarRadau

using StaticArrays: SMatrix, SVector

export radau, radau!

#-------------------------------------------------------------------------------
#RK tableau for RadauIIA 5th order

const c₁ = (4 - √6)/10
const c₂ = (4 + √6)/10
const c₃ = 1.0

const a₁₁ = (88 - 7*√6)/360
const a₁₂ = (296 - 169*√6)/1800
const a₁₃ = (-2 + 3*√6)/225

const a₂₁ = (296 + 169*√6)/1800
const a₂₂ = (88 + 7*√6)/360
const a₂₃ = (-2 - 3*√6)/225

const a₃₁ = (16 - √6)/36
const a₃₂ = (16 + √6)/36
const a₃₃ = 1/9

const e₁ = (-13 - 7*√6)/3
const e₂ = (-13 + 7*√6)/3
const e₃ = -1/3

#-------------------------------------------------------------------------------
#support functions

function ∂f∂y(𝒇::T, x, y, param, f₀, h, ϵ) where {T}
    #don't use a step size that risks roundoff error
    ∂y = max(ϵ*h, sqrt(eps(y)))
    #compute a regular old forward diff
    (𝒇(x, y + ∂y, param) - f₀)/∂y
end

function Jacobian(h, ∂)::SMatrix{3,3}
    #temporary
    q = h*∂
    #column-major storage
    SMatrix{3,3}(
        1.0 - a₁₁*q,
        -a₂₁*q,
        -a₃₁*q,
        -a₁₂*q,
        1.0 - a₂₂*q,
        -a₃₂*q,
        -a₁₃*q,
        -a₂₃*q,
        1.0 - a₃₃*q
    )
end

function xinit(x, h)
    x₁ = x + h*c₁
    x₂ = x + h*c₂
    x₃ = x + h*c₃
    return x₁, x₂, x₃
end

function hinit(x₀, xₙ, f, atol, rtol)
    x = max(abs(x₀), abs(xₙ))
    d = (1/x)^6 + abs(f)^6
    h = ((atol + rtol)/d)^(1/6)
    return min(h, xₙ - x₀)
end

#-------------------------------------------------------------------------------
# wrappers

function radau(𝒇::F,
               y₀::Real,
               x₀::Real,
               xₙ::Real,
               param=nothing;
               kwargs...
               ) where {F}
    radau!((), (), 𝒇, y₀, x₀, xₙ, param; kwargs...)
end

function radau(𝒇::F,
               y₀::Real,
               x₀::Real,
               xₙ::Real,
               nout::Int,
               param=nothing;
               kwargs...
               ) where {F}
    @assert nout > 1 "number of output points should be greater than 1"
    #make y float
    y₀ = float(y₀)
    #evenly spaced output points
    x = LinRange(x₀, xₙ, nout)
    #space for results
    y = zeros(typeof(y₀), nout)
    #integrate!
    radau!(y, x, 𝒇, y₀, x₀, xₙ, param; kwargs...)
    return x, y
end

#-------------------------------------------------------------------------------
#main function

function radau!(yout::Union{AbstractVector{<:Real},Tuple{}},
                xout::Union{AbstractVector{<:Real},Tuple{}},
                𝒇::F,
                y₀::Real,
                x₀::Real,
                xₙ::Real,
                param=nothing;
                rtol::Real=1e-6,
                atol::Real=1e-6,
                facmax::Real=100.0,
                facmin::Real=0.01,
                κ::Real=1e-3,
                ϵ::Real=0.25,
                maxnewt::Int=7,
                maxstep::Int=1000000,
                maxfail::Int=10) where {F}
    #basic checks
    @assert xₙ >= x₀
    @assert rtol < 1
    @assert facmax > 1
    @assert 0 < facmin < 1
    @assert 0 < κ < 1
    @assert 0 < ϵ < 1
    #set initial coordinates
    x, y = float(x₀), float(y₀)
    #initial function eval at x0
    f₀ = 𝒇(x, y, param)
    #output points
    nout = length(xout)
    jout = 1 #tracking index
    #initial step size selection
    h₁ = hinit(x, xₙ, f₀, atol, rtol)
    h₂ = hinit(x, xₙ, 𝒇(x + h₁, y + h₁*f₀, param), atol, rtol)
    h = min(h₁, h₂)
    #allocation, essentially, to keep f₃ in scope
    f₃ = zero(f₀)
    #counter
    nstep = 0
    while x < xₙ
        #don't overshoot the end of the integration interval
        h = min(h, xₙ - x)
        #finite diff ∂f/∂y, precision not necessary in practice, can also hurt
        ∂ = ∂f∂y(𝒇, x, y, param, f₀, h, ϵ)
        #jacobian matrix
        J = Jacobian(h, ∂)
        #x coordinates for function evaluations inside interval
        x₁, x₂, x₃ = xinit(x, h)
        #initial newton guesses, extrapolation appears to make things slower
        z₁, z₂, z₃ = zero(y), zero(y), zero(y)
        #newton iterations
        ΔZ = Inf # ∞ norm of changes to solution
        η = κ*(rtol*abs(y) + atol) #termination threshold
        nnewt = 0
        nfail = 0
        while ΔZ > η
            if nnewt == maxnewt
                #count the convergence failure
                nfail += 1
                #cut off unending disasters
                (nfail == maxfail) && error("repeated Newton convergence failures ($nfail) in Radau")
                #steeply reduce step size
                h /= 10.0
                #wipe the iteration counter
                nnewt = 0
                #reinitialize with the new step size
                J = Jacobian(h, ∂f∂y(𝒇, x, y, param, f₀, h, ϵ))
                x₁, x₂, x₃ = xinit(x, h)
                z₁, z₂, z₃ = zero(y), zero(y), zero(y)
            end
            #function evaluations
            f₁ = 𝒇(x₁, y + z₁, param)
            f₂ = 𝒇(x₂, y + z₂, param)
            f₃ = 𝒇(x₃, y + z₃, param)
            #newton system evaluation β = (h * Af) - z
            β = SVector{3}(
                h*(a₁₁*f₁ + a₁₂*f₂ + a₁₃*f₃) - z₁,
                h*(a₂₁*f₁ + a₂₂*f₂ + a₂₃*f₃) - z₂,
                h*(a₃₁*f₁ + a₃₂*f₂ + a₃₃*f₃) - z₃
            )
            #solve the linear system J*δ = β
            δ₁, δ₂, δ₃ = J\β
            #update
            z₁ += δ₁
            z₂ += δ₂
            z₃ += δ₃
            #norm of updates
            ΔZ = abs(δ₁) + abs(δ₂) + abs(δ₃)
            #count
            nnewt += 1
        end
        #scaled error estimate
        ze = (f₀*h + z₁*e₁ + z₂*e₂ + z₃*e₃)/5
        sc = rtol*abs(max(y, y + z₃)) + atol
        err = abs(ze)/sc
        #accept the step?
        if err < 1.0
            #reserve previous point for interpolation
            xₚ = x
            yₚ = y
            #advance the solution
            x += h
            y += z₃
            #dense output
            @inbounds if (jout <= nout) && ((x > xout[jout]) | (x ≈ xout[jout]))
                #set up cubic Hermite
                u = h*f₀
                v = h*f₃
                H₃ =  2yₚ  + u - 2y + v
                H₂ = -3yₚ - 2u + 3y - v
                #interpolate at all points that have been passed or met
                while (jout <= nout) && (x >= xout[jout])
                    ξ = (xout[jout] - xₚ)/h
                    yout[jout] += yₚ + ξ*(u + ξ*(H₂ + ξ*H₃))
                    jout += 1
                end
            end
            #f₃ is now at the beginning of the next interval
            f₀ = f₃
        end
        #count
        nstep += 1
        if nstep == maxstep
            error("maximum number of steps/attempts ($maxstep) reached in radau @ nx=$x y=$y h=$h")
        end
        #safety factor, loosely dependent on number of newton iterations
        facsaf = 0.9*(maxnewt + 1)/(maxnewt + nnewt)
        #step size selection, double sqrt faster than ^(1/4)
        h *= min(facmax, max(facmin, facsaf*sqrt(sqrt(1/err))))
    end
    return y
end

end


text:
#!/bin/env julia3

using HDF5
using ArgParse

function main()
    s = ArgParseSettings()
    @add_arg_table s begin
        "file1"
        required = true
        "file2"
        required = true
        "--all", "-a"
        help = "Print out fields that are in one file but not the other."
        action = :store_true
    end
    args = parse_args(s)
    f1 = h5open(args["file1"], "r")
    f2 = h5open(args["file2"], "r")
    gdiff(f1["/"], f2["/"], args)
    close(f1)
    close(f2)
end

function gdiff(o1, o2, args)
    for a in names(attrs(o1))
        if !in(a, names(attrs(o2)))
            args["all"] && println("Only in $(filename(o1))[$(name(o1))/$a]")
        end
    end
    for a in names(attrs(o2))
        if !in(a, names(attrs(o1)))
            args["all"] && println("Only in $(filename(o2))[$(name(o2))/$a]")
        else
            v1 = read(attrs(o1)[a])
            v2 = read(attrs(o2)[a])
            if v1 != v2
                println("Attr mismatch: $(name(o1))/$a]=$v1 vs $v2")
            end
        end
    end
    for a in names(o1)
        if !in(a, names(o2))
            args["all"] && println("Only in $(filename(o1))[$(name(o1))/$a]")
        end
    end
    for a in names(o2)
        if !in(a, names(o1))
            args["all"] && println("Only in $(filename(o2))[$(name(o2))/$a]")
        else
            c1 = o1[a]
            c2 = o2[a]
            if typeof(c1) != typeof(c2)
                println("Type mismatch: $(filename(c1))[$(name(c1))]=$(typeof(c1)) $(filename(c2))[$(name(c2))]=$(typeof(c2))")
            elseif isa(c1, HDF5Dataset)
                d1 = read(c1)
                d2 = read(c2)
                if size(d1) != size(d2)
                    println("Size mismatch: $(filename(c1))[$(name(c1))]=$(size(d1)) $(filename(c2))[$(name(c2))]=$(size(d2))")
                else
                    diff = maximum(abs(d1-d2))
                    println("maxdiff($(name(c1))) = $diff")
                end
            elseif isa(c1, HDF5Group)
                gdiff(c1, c2, args)
            else
                println("Don't know how to handle type $typeof(c1)")
            end
        end
    end
end

main()


text:
using Weave, Pkg, InteractiveUtils, IJulia

repo_directory = joinpath(@__DIR__)
cssfile = joinpath(@__DIR__, "templates", "skeleton_css.css")
latexfile = joinpath(@__DIR__, "templates", "julia_tex.tpl")

function weave_file(folder,file,build_list=(:script,:html,:github,:notebook); kwargs...)
  tmp = joinpath(repo_directory,"tutorials",folder,file)
  args = Dict{Symbol,String}(:folder=>folder,:file=>file)
  if :script ∈ build_list
    println("Building Script")
    dir = joinpath(repo_directory,"script",folder)
    isdir(dir) || mkdir(dir)
    args[:doctype] = "script"
    tangle(tmp;out_path=dir)
  end
  if :html ∈ build_list
    println("Building HTML")
    dir = joinpath(repo_directory,"html",folder)
    isdir(dir) || mkdir(dir)
    args[:doctype] = "html"
    weave(tmp,doctype = "md2html",out_path=dir,args=args; fig_ext=".svg", css=cssfile, kwargs...)
  end
  if :pdf ∈ build_list
    println("Building PDF")
    dir = joinpath(repo_directory,"pdf",folder)
    isdir(dir) || mkdir(dir)
    args[:doctype] = "pdf"
    weave(tmp,doctype="md2pdf",out_path=dir,args=args; template=latexfile, kwargs...)
  end
  if :github ∈ build_list
    println("Building Github Markdown")
    dir = joinpath(repo_directory,"markdown",folder)
    isdir(dir) || mkdir(dir)
    args[:doctype] = "github"
    weave(tmp,doctype = "github",out_path=dir,args=args; kwargs...)
  end
  if :notebook ∈ build_list
    println("Building Notebook")
    dir = joinpath(repo_directory,"notebook",folder)
    isdir(dir) || mkdir(dir)
    args[:doctype] = "notebook"
    weave(tmp,doctype = "notebook",out_path=dir,args=args; kwargs...)
    #Weave.convert_doc(tmp,joinpath(dir,file[1:end-4]*".ipynb"))
  end
end

function weave_all()
  for folder in readdir(joinpath(repo_directory,"tutorials"))
    folder == "appendix.jl" && continue
    weave_folder(folder)
  end
end

function weave_folder(folder)
  for file in readdir(joinpath(repo_directory,"tutorials",folder))
    println("Building $(joinpath(folder,file)))")
    try
      weave_file(folder,file)
    catch
    end
  end
end

#weave_all()

max_stars_repo_name: FSUcilab/sir-julia

max_issues_repo_name: FSUcilab/sir-julia

max_forks_repo_name: FSUcilab/sir-julia


text:
"""
    common.jl

Description:
    All common types, aliases, structs, and methods for the ClusterValidityIndices.jl package.
"""

# -------------------------------------------
# Abstract types
# -------------------------------------------

# Type for all CVIs
abstract type CVI end

# -------------------------------------------
# Aliases
# -------------------------------------------
#   **Taken from StatsBase.jl**
#
#  These types signficantly reduces the need of using
#  type parameters in functions (which are often just
#  for the purpose of restricting the arrays to real)
#
# These could be removed when the Base supports
# covariant type notation, i.e. AbstractVector{<:Real}

# Real-numbered aliases
const RealArray{T<:Real, N} = AbstractArray{T, N}
const RealVector{T<:Real} = AbstractArray{T, 1}
const RealMatrix{T<:Real} = AbstractArray{T, 2}

# Integered aliases
const IntegerArray{T<:Integer, N} = AbstractArray{T, N}
const IntegerVector{T<:Integer} = AbstractArray{T, 1}
const IntegerMatrix{T<:Integer} = AbstractArray{T, 2}

# Specifically floating-point aliases
const RealFP = Union{Float32, Float64}

# System's largest native floating point variable
const Float = (Sys.WORD_SIZE == 64 ? Float64 : Float32)

# Internal label mapping for incremental CVIs
const LabelMap = Dict{Int, Int}

# -------------------------------------------
# Methods
# -------------------------------------------

"""
    get_icvi!(cvi::CVI, sample::RealVector, label::Integer)

Compute and return the criterion value incrementally.

# Arguments
- `cvi::CVI`: the stateful information of the ICVI providing the criterion value.
- `sample::RealVector`: the sample provided to the external. clustering algorithm.
- `label::Integer`: the label prescribed to the `sample` by the external clustering algorithm.

# Examples
```julia
# Create a new CVI object
my_cvi = CH()
# Load in data from some external source
data = load_some_data()
# Cluster the data into a set of labels as an integer vector
labels = my_cluster_algorithm(data)
# Iteratively compute and extract the criterion value at every step
n_samples = length(labels)
criterion_values = zeros(n_samples)
for ix = 1:n_samples
    sample = data[:, ix]
    label = labels[ix]
    criterion_values[ix] = get_icvi!(my_cvi, sample, label)
end
```
"""
function get_icvi!(cvi::CVI, sample::RealVector, label::Integer)
    # Update the ICVI parameters
    param_inc!(cvi, sample, label)

    # Compute the criterion value
    evaluate!(cvi)

    # Return that value
    return cvi.criterion_value
end # get_icvi!(cvi::CVI, sample::RealVector, label::Integer)

"""
    get_cvi!(cvi::CVI, data::RealMatrix, labels::IntegerVector)

Compute and return the criterion value in batch mode.

# Arguments
- `cvi::CVI`: the stateful information of the CVI providing the criterion value.
- `data::RealMatrix`: a matrix of data, columns as samples and rows as features, used in the external clustering process.
- `labels::IntegerVector`: a vector of integers representing labels prescribed to the `data` by the external clustering algorithm.

# Examples
```julia
# Create a new CVI object
my_cvi = CH()
# Load in data from some external source
data = load_some_data()
# Cluster the data into a set of labels as an integer vector
labels = my_cluster_algorithm(data)
# Compute the final criterion value in batch mode
criterion_value = get_cvi!(cvi, data, labels)
```
"""
function get_cvi!(cvi::CVI, data::RealMatrix, labels::IntegerVector)
    # Update the CVI parameters in batch
    param_batch!(cvi, data, labels)

    # Compute the criterion value
    evaluate!(cvi)

    # Return that value
    return cvi.criterion_value
end # get_cvi!(cvi::CVI, data::RealMatrix, labels::IntegerVector)

"""
    get_internal_label!(label_map::LabelMap, label::Int)

Get the internal label and update the label map if the label is new.
"""
function get_internal_label!(label_map::LabelMap, label::Int)
    # If the label map contains the key, return that internal label
    if haskey(label_map, label)
        internal_label = label
    # Otherwise, increment the internal label to preserve monotonicity and store
    else
        internal_label = length(label_map) + 1
        label_map[label] = internal_label
    end

    return internal_label
end # get_internal_label!(label_map::LabelMap, label::Int)

# -------------------------------------------
# Common Documentation
# -------------------------------------------

@doc raw"""
    param_inc!(cvi::CVI, sample::RealVector, label::Integer)

Compute the CVI parameters incrementally.

This method updates only internal parameters of the ICVI algorithm incrementally.
When the criterion value itself is needed, use `evaluate!` and extract it from `cvi.criterion_value`.

# Arguments
- `cvi::CVI`: the stateful information of the CVI/ICVI algorithm.
- `sample::RealVector`: a vector of features used in the external clustering algorithm.
- `label::Integer`: the label that the external clustering algorithm prescribed to the `sample`.

# Examples
```julia-repl
julia> my_cvi = CH()
julia> data = load_some_data()
julia> labels = my_cluster_algorithm(data)
julia> param_inc!(my_cvi, data[:, 1], labels[1])
```
"""
param_inc!(cvi::CVI, sample::RealVector, label::Integer)

@doc raw"""
    param_batch!(cvi::CVI, data::RealMatrix, labels::IntegerVector)

Compute the CVI parameters in batch.

This method updates only the internal parameters of the CVI algorithm in batch.
When the criterion value itself is needed, use `evaluate!` and extract it from `cvi.criterion_value`.

# Arguments
- `cvi::CVI`: the stateful information of the CVI/ICVI algorithm.
- `data::RealMatrix`: a matrix of data where rows are features and columns are samples, used in the external clustering algorithm.
- `labels::IntegerVector`: a vector of labels that the external clustering algorithm prescribed to each column in `data`.

# Examples
```julia-repl
julia> my_cvi = CH()
julia> data = load_some_data()
julia> labels = my_cluster_algorithm(data)
julia> param_batch!(my_cvi, data, labels)
```
"""
param_batch!(cvi::CVI, data::RealMatrix, labels::IntegerVector)

@doc raw"""
    evaluate!(cvi::CVI)

Compute the criterion value of the CVI.

After computation, the resulting criterion value can be extracted from `cvi.criterion_value`.
The criterion value is a function of the CVI/ICVI internal parameters, so at least two classes (i.e., unique labels) must be presented to the CVI in `param_inc!` or `param_batch!` before a non-zero value is returned.

# Arguments
- `cvi::CVI`: the stateful information of the CVI/ICVI to use for computing the criterion value.

# Examples
```julia-repl
julia> my_cvi = CH()
julia> data = load_some_data()
julia> labels = my_cluster_algorithm(data)
julia> param_batch!(my_cvi, data, labels)
julia> evaluate!(my_cvi)
julia> my_criterion_value = my_cvi.criterion_value
```
"""
evaluate!(cvi::CVI)


text:
using HDF5,JLD,KUparser,KUnet

@date net = loadnet("foo11nnet3.jld")
#@date d = load("acl11.trn.jld4")
@date d = load("acl11.dev.jld4")
@date co = d["corpus"][1:100]
@show (ndeps,nbeam,nbatch,ncpu) = (length(d["deprel"]),64,5,nworkers())
@show ft = Flist.zn11cpv
@show pt = ArcEager13

# @date Main.restartmachines()
# @date Main.restartcuda()
# @everywhere using KUparser

@show (pt, length(co), ndeps, length(ft), length(net), nbeam, nbatch, ncpu)
@date p1 = bparse(pt, co, ndeps, ft, net, nbeam, nbatch, ncpu; xy=true)

eqparse(a,b)=(isequal(a[1],b[1]) && isequal(sortcols(vcat(a[2],a[3])), sortcols(vcat(b[2],b[3]))))
@date p2 = bparse(pt, co, ndeps, ft, net, nbeam, nbatch; xy=true)
@show eqparse(p1,p2)

# @date (p0,x0,y0) = bparse(pt, corpus[1:10], ndeps, ft, net, nbeam; xy=true)
# @date (p1,x1,y1) = bparse(pt, corpus, ndeps, ft, net, nbeam; xy=true)
# @date (p2,x2,y2) = bparse(pt, corpus, ndeps, ft, net, nbeam, nbatch; xy=true)
# @date (p3,x3,y3) = bparse(pt, corpus, ndeps, ft, net, nbeam, nbatch, ncpu; xy=true)

# @date bparse(pt, co, ndeps, ft, net, nbeam, 128; xy=true)
# @date bparse(pt, co, ndeps, ft, net, nbeam, 64; xy=true)
# @date bparse(pt, co, ndeps, ft, net, nbeam, 32; xy=true)
# @date bparse(pt, co, ndeps, ft, net, nbeam, 16; xy=true)
# @date bparse(pt, co, ndeps, ft, net, nbeam, 8; xy=true)
# @date bparse(pt, co, ndeps, ft, net, nbeam, 4; xy=true)
# @date bparse(pt, co, ndeps, ft, net, nbeam, 2; xy=true)
# @date bparse(pt, co, ndeps, ft, net, nbeam; xy=true)

:ok


text:
# ---
# title: "`closing`"
# id: op_closing
# cover: assets/closing.png
# ---

# `closing` operator is defined as `dilate(erode(img))`. Intuitively, closing operation
# fills the black holes in the image.

using ImageMorphology
using TestImages
using ImageBase
using ImageShow

img = restrict(testimage_dip3e("fig0940")) # rice
img01 = @. Gray(img > 0.5)
mosaic(img, img01; nrow=1)

#-

out1 = closing(img) # default: all spatial dimensions, r=1, a box-shape SE
out2 = closing(img; dims=(2,)) # only apply to the second dimension
out3 = closing(img; r=5) # half-size r=5

## also to the binary version
out1_01 = closing(img01)
out2_01 = closing(img01; dims=(2,))
out3_01 = closing(img01; r=5)

mosaic(out1, out2, out3, out1_01, out2_01, out3_01; nrow=2, rowmajor=true)

# ## See also

# [`closing`](@ref op_closing) is the dual operator of `closing` in the following sense:

complement.(closing(img)) == closing(complement.(img))

# For a comprehensive and more accurate documentation, please check the [`closing`](@ref)
# reference page.

# save cover image #src
using FileIO #src
mkpath("assets") #src
img = Gray.(restrict(testimage("blobs"))) #src
save("assets/closing.png", closing(img; r=2)) #src


text:
using CRRao, Test, StableRNGs, Logging, RDatasets

Logging.disable_logging(Logging.Warn)

CRRao.setprogress!(false)
CRRao.set_rng(StableRNG(123))

@testset "CRRao.jl" begin
    @testset "Basic Tests" begin
        @testset "Linear Regression" begin
            include("basic/LinearRegression.jl")
        end

        @testset "Logistic Regression" begin
            include("basic/LogisticRegression.jl")
        end

        @testset "Poisson Regression" begin
            include("basic/PoissonRegression.jl")
        end

        @testset "Negative Binomial Regression" begin
            include("basic/NegBinomialRegression.jl")
        end
    end
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/AbaqusReader.jl/blob/master/LICENSE

using AbaqusReader: abaqus_read_mesh
using Test

datadir = first(splitext(basename(@__FILE__)))
filename = joinpath(datadir, "mesh.inp")
mesh = abaqus_read_mesh(filename)
@test length(mesh["elements"]) == 7


text:
using Base.Test
using DataFrames

let
	@assert NA | true == true
	@assert isna(NA | false)
	@assert isna(NA | NA)
	@assert true | NA == true
	@assert isna(false | NA)

	@assert isna(NA & true)
	@assert NA & false == false
	@assert isna(NA & NA)
	@assert isna(true & NA)
	@assert false & NA == false

	@assert any(DataVector[1, 2, NA] .== 1) == true
	@assert any(DataVector[NA, 1, 2] .== 1) == true
	@assert isna(any(DataVector[1, 2, NA] .== 3))
	@assert any(DataVector[1, 2, 3] .== 4) == false

	@assert isna(all(DataVector[1, 1, NA] .== 1))
	@assert isna(all(DataVector[NA, 1, 1] .== 1))
	@assert all(DataVector[1, 1, 1] .== 1) == true
	@assert all(DataVector[1, 2, 1] .== 1) == false
end


text:
using Test
using Adapt: adapt
@testset "testing adapt to Float32" begin
	@test typeof(adapt(Float32,Flux.param(randn(3,3,3)))) <: Flux.TrackedArray{Float32,3}
	@test typeof(adapt(Float32,Flux.param(randn(3,3)))) <: Flux.TrackedArray{Float32,2}
	@test typeof(adapt(Float32,Flux.param(randn(3)))) <: Flux.TrackedArray{Float32,1}
	@test typeof(adapt(Float32,Flux.Dense(3,3)).W) <: Flux.TrackedArray{Float32,2}
	@test typeof(adapt(Float32,Flux.Dense(3,3)).b) <: Flux.TrackedArray{Float32,1}
	@test typeof(adapt(Float32, ResDense(3, 3)).a.W) <: Flux.TrackedArray{Float32,2}
end

max_stars_repo_name: pevnak/FluxExtensions

max_issues_repo_name: pevnak/FluxExtensions

max_forks_repo_name: pevnak/FluxExtensions


text:
export dropout

doc"""
    dropout(x::Var, droprate::Float64)

Drops elements randomly with probability ``droprate`` and scales the other elements by factor ``1 / (1 - droprate)``.
"""
function dropout(x::Var, droprate::Float64)
    droprate == 0.0 && return x
    istrain() || return x
    ydata, work = dropout(x.data, droprate)
    Var(ydata, ∇dropout!, (x,droprate,work))
end
dropout(x::Node, droprate) = Node(dropout, (x,droprate))

function dropout(x::Array{T}, droprate::Float64) where T
    work = rand(T, length(x))
    scale = T(1 / (1-droprate))
    y = similar(x)
    @inbounds for i = 1:length(x)
        y[i] = work[i] <= droprate ? T(0) : scale*x[i]
    end
    y, work
end

dropout(x::CuArray, droprate) = CUDNN.dropout(x, droprate)

function ∇dropout!(y::Var, x::Var, droprate::Float64, work)
    isnothing(x.grad) && return
    ∇dropout!(y.grad, x.grad, droprate, work)
end

function ∇dropout!(gy::Array{T}, gx::Array{T}, droprate::Float64, work::Vector{T}) where T
    scale = T(1 / (1-droprate))
    @inbounds for i = 1:length(gx)
        gx[i] += work[i] <= droprate ? T(0) : scale*gy[i]
    end
end

∇dropout!(gy::CuArray, gx::CuArray, droprate, dropdesc) = CUDNN.∇dropout!(gy, gx, dropdesc)


text:
module Model

using CSV
using DataFrames
using Distances
using JuMP
using Gurobi

import DispatchOps

include("interface.jl")

end


text:
macro publish(mod,name)
  quote
    using Gridap.$mod: $name; export $name
  end
end

# Reexport from LinearAlgebra (just for convenience)
using LinearAlgebra:  det, inv, tr, cross, dot, norm, ×, ⋅
export det, inv, tr, cross, dot, norm, ×, ⋅

@publish Helpers operate
@publish Helpers GridapType

@publish Algebra solve
@publish Algebra solve!
@publish Algebra symbolic_setup
@publish Algebra numerical_setup
@publish Algebra numerical_setup!
@publish Algebra LUSolver
@publish Algebra BackslashSolver
@publish Algebra zero_initial_guess
@publish Algebra NLSolver
@publish Algebra get_matrix
@publish Algebra get_vector
@publish Algebra SparseMatrixCSR
@publish Algebra SymSparseMatrixCSR

@publish Arrays array_cache
@publish Arrays getindex!
@publish Arrays get_array
@publish Arrays apply
@publish Arrays reindex

@publish TensorValues VectorValue
@publish TensorValues TensorValue
@publish TensorValues inner
@publish TensorValues outer
@publish TensorValues diagonal_tensor
@publish TensorValues num_components
using Gridap.TensorValues: ⊙; export ⊙
using Gridap.TensorValues: ⊗; export ⊗

@publish Fields gradient
@publish Fields ∇
@publish Fields integrate
@publish Fields Point
@publish Fields evaluate
@publish Fields evaluate!
@publish Fields curl
@publish Fields laplacian
@publish Fields divergence
@publish Fields Δ
@publish Fields ε
@publish Fields symmetric_gradient

@publish Integration get_coordinates
@publish Integration get_weights

@publish ReferenceFEs is_simplex
@publish ReferenceFEs is_n_cube
@publish ReferenceFEs simplexify
@publish ReferenceFEs num_dims
@publish ReferenceFEs num_cell_dims
@publish ReferenceFEs num_point_dims
@publish ReferenceFEs VERTEX
@publish ReferenceFEs SEGMENT
@publish ReferenceFEs TRI
@publish ReferenceFEs QUAD
@publish ReferenceFEs TET
@publish ReferenceFEs HEX
@publish ReferenceFEs WEDGE
@publish ReferenceFEs PYRAMID
@publish ReferenceFEs LagrangianRefFE
@publish ReferenceFEs is_first_order
@publish ReferenceFEs is_affine
@publish ReferenceFEs is_Q
@publish ReferenceFEs is_P
@publish ReferenceFEs is_S
@publish ReferenceFEs VERTEX1
@publish ReferenceFEs SEG2
@publish ReferenceFEs TRI3
@publish ReferenceFEs QUAD4
@publish ReferenceFEs TET4
@publish ReferenceFEs HEX8
@publish ReferenceFEs SerendipityRefFE
@publish ReferenceFEs Polytope

@publish Geometry get_triangulation
@publish Geometry num_cells
@publish Geometry Triangulation
@publish Geometry get_normal_vector
@publish Geometry get_cell_coordinates
@publish Geometry restrict
@publish Geometry get_physical_coordinate
@publish Geometry get_cell_map
@publish Geometry CartesianGrid
@publish Geometry CartesianDiscreteModel
@publish Geometry DiscreteModel
@publish Geometry DiscreteModelFromFile
@publish Geometry num_tags
@publish Geometry num_entities
@publish Geometry get_grid
@publish Geometry get_face_labeling
@publish Geometry add_tag!
@publish Geometry add_tag_from_tags!
@publish Geometry BoundaryTriangulation
@publish Geometry jump
@publish Geometry mean
@publish Geometry SkeletonTriangulation
@publish Geometry RestrictedTriangulation
@publish Geometry InterfaceTriangulation

@publish CellData CellQuadrature
@publish CellData QPointCellField
@publish CellData CellField

@publish FESpaces FESpace
@publish FESpaces TrialFESpace
@publish FESpaces TestFESpace
@publish FESpaces FETerm
@publish FESpaces FEEnergy
@publish FESpaces AffineFETerm
@publish FESpaces LinearFETerm
@publish FESpaces FESource
@publish FESpaces AffineFEOperator
@publish FESpaces LinearFESolver
@publish FESpaces get_free_values
@publish FESpaces get_dirichlet_values
@publish FESpaces num_dirichlet_dofs
@publish FESpaces num_free_dofs
@publish FESpaces num_dirichlet_tags
@publish FESpaces FEFunction
@publish FESpaces interpolate
@publish FESpaces interpolate_everywhere
@publish FESpaces interpolate_dirichlet
@publish FESpaces FEOperator
@publish FESpaces FESolver
@publish FESpaces update_state_variables!

using Gridap.FESpaces: @law; export @law

@publish MultiField MultiFieldFESpace

@publish Visualization writevtk
@publish Visualization createvtk
@publish Visualization print_op_tree


text:
"""A package for graph computations related to graph Laplacians

Graphs are represented by sparse adjacency matrices, etc.
"""
module Laplacians

#=

Started by Dan Spielman.
Other contributors:

This represents graphs as sparse adjacency matrices.
The reasons are twofold:
  # 1. It is much faster than the Graphs.jl library, and
  # 2. It is more compatible with linear algebraic operations

This module begins by importing the packages it requires.
It then defines a few functions for dealing with the graphs:
deg(graph, node),  nbri(graph, node, i)  and weighti(graph, node i).
Using these is slower than what we actually do in the code.
But, you can start with this and then convert.

This then includes code from other files,
and exports the functions for which it seems appropriate.

=#

  function __init__()
    if !isdefined(Main, :LAPLACIANS_NOPLOT)
        eval(Expr(:using, :PyPlot))
    end

    if !isdefined(Main, :LAPLACIANS_NOAMG)
        eval(Expr(:using, :PyAMG))
    end
  end

  using DataStructures

  include("fastCSC.jl")
  export symPermuteCSC
  export symTransposeCSC
  export submatrixCSC

  include("graphUtils.jl")
  export deg
  export nbri
  export weighti
  export nbrs
  export wdeg
  export setValue
  export backIndices
  export flipIndex
  export findEntries
  export compConductance
  export getVolume
  export getObound

  include("graphGenerators.jl")
  export readIJ
  export ringGraph
  export generalizedRing
  export generalizedNecklace
  export randMatching
  export randRegular
  export grownGraph
  export grownGraphD
  export prefAttach
  export hyperCube
  export completeBinaryTree
  export completeGraph
  export pathGraph

  export wGrid2
  export wGrid3

  export grid2
  export grid2coords

  export randGenRing
  export randperm
  export ErdosRenyi
  export ErdosRenyiCluster
  export ErdosRenyiClusterFix
  export pureRandomGraph

  export chimera
  export randWeight
  export wtedChimera, semiWtedChimera

  include("IO.jl")
  export readIJ, readIJV, writeIJV

  include("graphOps.jl")

  export unweight, unweight!
  export mapweight
  export uniformWeight, uniformWeight!

  export edgeVertexMat

  export productGraph
  export subsampleEdges

  export twoLift
  export joinGraphs, disjoin

  export plotGraph

  export shortIntGraph, floatGraph

  export lap
  export adj
  export spectralCoords
  export spectralDrawing

  export toUnitVector

  export diagmat

  include("graphAlgs.jl")

  export components
  export biggestComp
  export vecToComps
  export isConnected

  export shortestPaths, shortestPathTree, pathFromParents
  export kruskal, prim

  include("treeAlgs.jl")

  export RootedTree
  export matToTree
  export matToTreeDepth
  export tarjanStretch
  export compDepth
  export compStretches
  export dfsOrder

  include("pcg.jl")

  export cg, cgSolver
  export pcg, pcgSolver, pcgLapSolver

  include("flow.jl")

  export maxflow

  include("akpw.jl")

  export akpw, akpwU

  include("localClustering.jl")

  export prn
  export apr
  export localImprove

  include("cutHeuristics.jl")

  export refineCut
  export dumb

  include("randTrees.jl")
  export randishKruskal, randishPrim

  include("sampler.jl")
  include("fastSampler.jl")

  include("solvers.jl")
  export lapWrapSolver, lapChol, augmentTree, augTreePrecon, augTreeSolver
  export augTreeLapPrecon, augTreeLapSolver
  export AMGSolver, AMGLapSolver

  include("complexSolvers.jl")
  export SDDSolvers
  export LapSolvers

  include("johnlind.jl")
  export johnlind

  include("toposort.jl")
  export toposort, dirEdgeVertexMat

  # include("isotonicIPM.jl")
  # export isotonicIPM, isotonicIPMrelEps

end # module yinsGraph