text:
@views function StrainRate!( Eps::Tensor2D, Vx::Matrix{Float64}, Vy::Matrix{Float64}, dom::ModelDomain )
    dx, dy = dom.dx, dom.dy
    @tturbo Eps.div    .= diff(Vx[:,2:end-1],dims=1)./dx .+ diff(Vy[2:end-1,:],dims=2)./dy
    # ncx, ncy = size(div,1), size(div,2)
    # @tturbo for i=1:ncx
    #     for j=1:ncy
    #         div[i,j] = 1.0/dx * (Vx[i+1,j+1] - Vx[i,j+1]) + 1.0/dy * (Vy[i+1,j+1] - Vy[i+1,j])
    #     end
    # end
    @tturbo Eps.xx .= diff(Vx[:,2:end-1],dims=1)./dx .- 1.0/3.0 .* Eps.div
    @tturbo Eps.yy .= diff(Vy[2:end-1,:],dims=2)./dy .- 1.0/3.0 .* Eps.div
    @tturbo Eps.zz .= .-(Eps.xx .+ Eps.yy)
    @tturbo Eps.xy .= 0.5.*( diff(Vx,dims=2)./dy .+ diff(Vy,dims=1)/dx )
    VerticesToCentroids!( Eps.xy_c, Eps.xy )
    @tturbo Eps.II .= sqrt.(0.5*(Eps.xx.^2 .+ Eps.xx.^2 .+ Eps.zz.^2) .+ Eps.xy_c.^2)
end
export StrainRate!

########

@views function Stress!( Tau::Tensor2D, Eps::Tensor2D, f::Fields2D )
    @tturbo Tau.xx .= 2.0.*f.etac.*Eps.xx
    @tturbo Tau.yy .= 2.0.*f.etac.*Eps.yy
    @tturbo Tau.zz .= 2.0.*f.etac.*Eps.zz
    @tturbo Tau.xy .= 2.0.*f.etav.*Eps.xy
    VerticesToCentroids!( Tau.xy_c, Tau.xy )
    @tturbo Tau.II .= sqrt.(0.5*(Tau.xx.^2 .+ Tau.xx.^2 .+ Tau.zz.^2) .+ Tau.xy_c.^2)
end
export Stress!

########

@views function StressVE!( Tau::Tensor2D, Tau0::Tensor2D, Eps::Tensor2D, f::Fields2D, dt::Float64 )
    @tturbo Tau.xx .= 2f.etac.* ( Eps.xx .+ Tau0.xx ./ ( 2f.Gc.*dt ) )
    @tturbo Tau.yy .= 2f.etac.* ( Eps.yy .+ Tau0.yy ./ ( 2f.Gc.*dt ) )
    @tturbo Tau.zz .= 2f.etac.* ( Eps.zz .+ Tau0.zz ./ ( 2f.Gc.*dt ) )
    @tturbo Tau.xy .= 2f.etav.* ( Eps.xy .+ Tau0.xy ./ ( 2f.Gv.*dt ) )
    VerticesToCentroids!( Tau.xy_c, Tau.xy )
    @tturbo Tau.II .= sqrt.(0.5*((f.Damc.*Tau.xx).^2 .+ (f.Damc.*Tau.yy).^2 .+ (f.Damc.*Tau.zz).^2) .+ (f.Damc.*Tau.xy_c).^2)
end
export StressVE!

########
function StrainEnergy!( We::Matrix{Float64}, Tau::Tensor2D, Eps::Tensor2D, Str0::Tensor2D, Str::Tensor2D, P::Matrix{Float64}, dt::Float64, f::Fields2D, dx::Float64, dy::Float64 )
    Str.xx   .= Str0.xx   .+ dt.*(Eps.xx .+ 1.0/3.0.*Eps.div)
    Str.yy   .= Str0.yy   .+ dt.*(Eps.yy .+ 1.0/3.0.*Eps.div)
    Str.zz   .= Str0.zz   .+ dt.*(Eps.zz .+ 1.0/3.0.*Eps.div)
    Str.xy_c .= Str0.xy_c .+ dt.*(Eps.xy_c)
    @tturbo Str.II .= sqrt.(0.5*(Str.xx.^2 .+ Str.xx.^2 .+ Str.zz.^2) .+ Str.xy_c.^2)
    We  .= (Tau.xx .- P) .* Str.xx .+ (Tau.yy .- P) .* Str.yy .+ (Tau.zz .- P) .* Str.zz .+ 2.0.*Tau.xy_c .* Str.xy_c
    We .*= 0.5
    # We ./= f.Damc
    We[f.We0.>We] .= f.We0[f.We0.>We]

    # Exx = diff(f.Vx[:,2:end-1],dims=1)./dx * dt
    # Eyy = diff(f.Vy[2:end-1,:],dims=2)./dy * dt
    # Exyv = 0.5.*( diff(f.Vx,dims=2)./dy .+ diff(f.Vy,dims=1)/dx ) * dt
    # Exyc = zeros(Float64, size(Exx))
    # VerticesToCentroids!( Exyc, Exyv )

    # Sxx = (f.Kc.+4.0./3.0*f.Gc).*Exx + (f.Kc.-2.0./3.0*f.Gc).*Eyy
    # Syy = (f.Kc.-2.0./3.0*f.Gc).*Exx + (f.Kc.+4.0./3.0*f.Gc).*Eyy
    # Szz = (f.Kc.-4.0./3.0*f.Gc).*Exx + (f.Kc.-2.0./3.0*f.Gc).*Eyy
    # Sxy = 2.0*f.Gc.*Exyc

    # We .= 2.0*Sxy.*Exyc .+ Sxx.*Exx .+ Syy.*Eyy
end
export StrainEnergy!
########

@views function Residuals!( f::Fields2D, Tau::Tensor2D, Eps::Tensor2D, BC::BoundaryConditions, dom::ModelDomain, params::ModelParameters )
    dx, dy = dom.dx, dom.dy
    @tturbo f.Fx .= 0.0
    @tturbo f.Fy .= 0.0
    @tturbo f.Fp .= 0.0
    if BC.periodix==0
        f.Fx[2:end-1,:] .= ( diff(Tau.xx .- f.Pc, dims=1)./dx .+ diff(Tau.xy[2:end-1,:], dims=2)/dy )
    else
        Sxx_ex             = zeros(ncx+1,ncy)
        @tturbo Sxx_ex[2:end-0,:] .= -f.Pc .+ Tau.xx
        Sxx_ex[      1,:] .= -f.Pc[end,:] .+ Tau.xx[end,:]
        #Sxx_ex[    end,:] .= -Pc[  1,:] .+ Tau.xx[  1,:] # Do not assemble last column
        f.Fx .= ( diff(Sxx_ex, dims=1)./dx .+ diff(Tau.xy[1:end-1,:], dims=2)/dy )
    end
    f.Fy[:,2:end-1] .= ( diff(Tau.yy .- f.Pc, dims=2)./dy .+ diff(Tau.xy[:,2:end-1], dims=1)/dx )
    @tturbo f.Fp            .= -Eps.div
    # # For periodic
    # if BC.periodix==1
    #     Fx = Fx[1:end-1,:]
    # end
end
export Residuals!

########

@views function ResidualsComp!( f::Fields2D, Tau::Tensor2D, Eps::Tensor2D, BC::BoundaryConditions, dom::ModelDomain, params::ModelParameters )
    dx, dy, ncx, ncy = dom.dx, dom.dy, dom.ncx, dom.ncy
    @tturbo f.Fx .= 0.0
    @tturbo f.Fy .= 0.0
    @tturbo f.Fp .= 0.0
    if BC.periodix==0
        f.Fx[2:end-1,:] .= ( diff(Tau.xx .- f.Pc, dims=1)./dx .+ diff(Tau.xy[2:end-1,:], dims=2)/dy )
    else
        Sxx_ex             = zeros(ncx+1,ncy)
        @tturbo Sxx_ex[2:end-0,:] .= -f.Pc .+ Tau.xx
        Sxx_ex[      1,:] .= -f.Pc[end,:] .+ Tau.xx[end,:]
        #Sxx_ex[    end,:] .= -Pc[  1,:] .+ Tau.xx[  1,:] # Do not assemble last column
        f.Fx .= ( diff(Sxx_ex, dims=1)./dx .+ diff(Tau.xy[1:end-1,:], dims=2)/dy )
    end
    f.Fy[:,2:end-1] .= ( diff(Tau.yy .- f.Pc, dims=2)./dy .+ diff(Tau.xy[:,2:end-1], dims=1)/dx )
    if params.comp == 1
        @tturbo f.Fp .= -Eps.div .- (f.Pc .- f.Pc0) ./ (f.Kc.*params.dt)
    else
        @tturbo f.Fp .= -Eps.div
    end
end
export ResidualsComp!

########

@views function ResidualsCompDam!( f::Fields2D, Tau::Tensor2D, Eps::Tensor2D, BC::BoundaryConditions, dom::ModelDomain, params::ModelParameters )
    dx, dy, ncx, ncy = dom.dx, dom.dy, dom.ncx, dom.ncy
    @tturbo f.Fx .= 0.0
    @tturbo f.Fy .= 0.0
    @tturbo f.Fp .= 0.0
    if BC.periodix==0
        f.Fx[2:end-1,:] .= ( diff(f.Damc.*(Tau.xx .- f.Pc), dims=1)./dx .+ diff(f.Damv[2:end-1,:].*Tau.xy[2:end-1,:], dims=2)/dy )
    else
        Sxx_ex             = zeros(ncx+1,ncy)
        @tturbo Sxx_ex[2:end-0,:] .= -f.Pc .+ Tau.xx
        Sxx_ex[      1,:] .= -f.Pc[end,:] .+ Tau.xx[end,:]
        #Sxx_ex[    end,:] .= -Pc[  1,:] .+ Tau.xx[  1,:] # Do not assemble last column
        f.Fx .= ( diff(Sxx_ex, dims=1)./dx .+ diff(Tau.xy[1:end-1,:], dims=2)/dy )
    end
    f.Fy[:,2:end-1] .= ( diff(f.Damc.*(Tau.yy .- f.Pc), dims=2)./dy .+ diff(f.Damv[:,2:end-1].*Tau.xy[:,2:end-1], dims=1)/dx )
    if params.comp == 1
        @tturbo f.Fp .= -Eps.div .- (f.Pc .- f.Pc0) ./ (f.Kc.*params.dt)
    else
        @tturbo f.Fp .= -Eps.div
    end
end
export ResidualsCompDam!

########

@views function SetBCs( Vx::Matrix{Float64}, Vy::Matrix{Float64}, BC::BoundaryConditions )
    # Boundaries
    if BC.Vx.type_S == 22 Vx[:,  1] .= Vx[:,    2] end
    if BC.Vx.type_S == 11 Vx[:,  1] .= 2.0.*BC.Vx.Dir_S .- Vx[:,    2] end
    if BC.Vx.type_N == 22 Vx[:,end] .= Vx[:,end-1] end
    if BC.Vx.type_N == 11 Vx[:,end] .= 2.0.*BC.Vx.Dir_N .- Vx[:,end-1] end
    if BC.periodix == 1
        Vy[1,  :] .= Vy[end-1,:]
        Vy[end,:] .= Vy[    2,:]
    else
        if BC.Vy.type_W == 22  Vy[1,  :] .= Vy[2,    :] end
        if BC.Vy.type_W == 11  Vy[1,  :] .= 2.0.*BC.Vy.Dir_W .- Vy[2,    :] end
        if BC.Vy.type_E == 22  Vy[end,:] .= Vy[end-1,:] end
        if BC.Vy.type_E == 11  Vy[end,:] .= 2.0.*BC.Vy.Dir_E .- Vy[end-1,:] end
    end
end
export SetBCs

########

function NumberingStokes!( f::Fields2D, BC::BoundaryConditions, dom::ModelDomain )
    ncx, ncy = dom.ncx, dom.ncy
    # Numbering
    if BC.periodix==0
        f.NumVx     = collect(reshape(1:(ncx+1)*ncy,ncx+1,ncy))
    else
        f.NumVx             = zeros(Int64,(ncx+1),ncy)
        f.NumVx[1:end-1,:] .= collect(reshape(1:(ncx+0)*ncy,ncx,ncy))
        f.NumVx[end,:]     .= f.NumVx[1,:]
    end
    f.NumVy     = collect(reshape(1:ncx*(ncy+1),ncx,ncy+1) .+ maximum(f.NumVx))
    f.NumP      = collect(reshape(1:(ncx)*ncy,ncx,ncy))
    return
end
export NumberingStokes!

########

@views function StokesAssembly( BC::BoundaryConditions, f::Fields2D, DirScale::Float64, dom::ModelDomain )

    ncx, ncy = dom.ncx, dom.ncy
    dx, dy = dom.dx, dom.dy
    NumVx, NumVy, NumP = f.NumVx, f.NumVy, f.NumP
    etac, etav = f.Damc.*f.etac, f.Damv.*f.etav

    if BC.periodix==0
        nxvx = ncx+1
        sx   = 0
    else
        nxvx = ncx+1#ncx # remove last column
        sx   = 0#1   # shift
    end

    # Connectivity
    iVxC      =  ones(Int64, nxvx, ncy);
    iVxC     .=  NumVx[1:end-sx,:]
    iVxW      =  ones(Int64, nxvx, ncy); iVxW[2:end-0,: ] = NumVx[1:end-1,:]
    iVxE      =  ones(Int64, nxvx, ncy); iVxE[1:end-1,: ] = NumVx[2:end-0,:]
    iVxS      =  ones(Int64, nxvx, ncy); iVxS[: ,2:end-0] = NumVx[:,1:end-1]
    iVxN      =  ones(Int64, nxvx, ncy); iVxN[: ,1:end-1] = NumVx[:,2:end-0]
    iVySW     =  ones(Int64, nxvx, ncy); iVySW[2:end-0,:] = NumVy[:,1:end-1]
    iVySE     =  ones(Int64, nxvx, ncy); iVySE[1:end-1,:] = NumVy[:,1:end-1]
    iVyNW     =  ones(Int64, nxvx, ncy); iVyNW[2:end-0,:] = NumVy[:,2:end-0]
    iVyNE     =  ones(Int64, nxvx, ncy); iVyNE[1:end-1,:] = NumVy[:,2:end-0]
    iPW       =  ones(Int64, nxvx, ncy); iPW[2:end-0,:]   = NumP[:,:]
    iPE       =  ones(Int64, nxvx, ncy); iPE[1:end-1,:]   = NumP[:,:]

    # Viscosity coefficients
    etaW      = zeros(Float64, nxvx, ncy); etaW[2:end-0,:] = etac[1:end-0,:]
    etaE      = zeros(Float64, nxvx, ncy); etaE[1:end-1,:] = etac[1:end-0,:]
    etaS      = zeros(Float64, nxvx, ncy); etaS[:,1:end-0] = etav[:,1:end-1]
    etaN      = zeros(Float64, nxvx, ncy); etaN[:,1:end-0] = etav[:,2:end-0]

    DamW =  ones(Float64, size(NumVx)); DamW[2:end-0,:]   .= f.Damc
    DamE =  ones(Float64, size(NumVx)); DamE[1:end-1,:]   .= f.Damc

    if BC.periodix==1
        DamW[  1,:] = Damc[end,:]
        DamE[end,:] = Damc[1,:]
        etaW[  1,:] = etac[end,:]
        etaE[end,:] = etac[1,:]
        iVxW[  1,:] = NumVx[end-1,:]
        iVxE[end,:] = NumVx[    2,:]
        iVySW[  1,:]  = NumVy[end,1:end-1]
        iVyNW[  1,:]  = NumVy[end,2:end-0]
        iVySE[end,:]  = NumVy[  1,1:end-1]
        iVyNE[end,:]  = NumVy[  1,2:end-0]
        iPW[  1,:] = NumP[end,:]
        iPE[end,:] = NumP[  1,:]
    end

    # Finite difference coefficients
    @tturbo cVxC  = -(-1.0.*etaN./dy - 1.0.*etaS./dy)./dy - (-4/3*etaE./dx - 4/3*etaW./dx)./dx
    @tturbo cVxW  = -4/3*etaW./dx.^2
    @tturbo cVxE  = -4/3*etaE./dx.^2
    @tturbo cVxS  = -1.0*etaS./dy.^2
    @tturbo cVxN  = -1.0*etaN./dy.^2
    @tturbo cVySW = -1.0*etaS./(dx.*dy) + 2/3*etaW./(dx.*dy)
    @tturbo cVySE = -2/3*etaE./(dx.*dy) + 1.0*etaS./(dx.*dy)
    @tturbo cVyNW = 1.0*etaN./(dx.*dy) - 2/3*etaW./(dx.*dy)
    @tturbo cVyNE = 2/3*etaE./(dx.*dy) - 1.0*etaN./(dx.*dy)
    @tturbo cPW   = -DamW./dx .*  ones(Float64, nxvx, ncy)
    @tturbo cPE   =  DamE./dx .*  ones(Float64, nxvx, ncy)

    if BC.Vx.type_S==11
        cVxC[:,  1] .-= cVxS[:,  1]
    end
    if BC.Vx.type_N==11
        cVxC[:,end] .-= cVxN[:,end]
    end
    if BC.Vx.type_S==22
        cVxC[:,  1] .+= cVxS[:,  1]
    end
    if BC.Vx.type_N==22
        cVxC[:,end] .+= cVxN[:,end]
    end
    cVxS[:,  1] .= 0.0
    cVxN[:,end] .= 0.0

    # Symmetry - kill Dirichlet connections
    cVySW[  :,  1] .= 0.0
    cVySE[  :,  1] .= 0.0
    cVyNW[  :,end] .= 0.0
    cVyNE[  :,end] .= 0.0

    if BC.periodix==0
        cVxS[:,  1] .= 0.0;
        cVxN[:,end] .= 0.0;
        cVxC[1,:]  .= DirScale; cVxC[end,:]  .= DirScale
        cVxW[1,:]  .= 0.0;      cVxW[end,:]  .= 0.0
        cVxE[1,:]  .= 0.0;      cVxE[end,:]  .= 0.0
        cVxS[1,:]  .= 0.0;      cVxS[end,:]  .= 0.0
        cVxN[1,:]  .= 0.0;      cVxN[end,:]  .= 0.0
        cVySW[1,:] .= 0.0;      cVySW[end,:] .= 0.0
        cVyNW[1,:] .= 0.0;      cVyNW[end,:] .= 0.0
        cVySE[1,:] .= 0.0;      cVySE[end,:] .= 0.0
        cVyNE[1,:] .= 0.0;      cVyNE[end,:] .= 0.0
        cPW[1,:]   .= 0.0;      cPW[end,:]   .= 0.0
        cPE[1,:]   .= 0.0;      cPE[end,:]   .= 0.0
        # Symmetry - kill Dirichlet connections
        cVxW[    2,:] .= 0.0
        cVxE[end-1,:] .= 0.0
    end

    ###################

    # Connectivity
    iVyC      = NumVy
    iVyW      =  ones(Int64, size(NumVy)); iVyW[2:end-0,: ] = NumVy[1:end-1,:]
    iVyE      =  ones(Int64, size(NumVy)); iVyE[1:end-1,: ] = NumVy[2:end-0,:]
    iVyS      =  ones(Int64, size(NumVy)); iVyS[: ,2:end-0] = NumVy[:,1:end-1]
    iVyN      =  ones(Int64, size(NumVy)); iVyN[: ,1:end-1] = NumVy[:,2:end-0]
    iVxSW     =  ones(Int64, size(NumVy)); iVxSW[:,2:end-0] = NumVx[1:end-1,:]
    iVxSE     =  ones(Int64, size(NumVy)); iVxSE[:,2:end-0] = NumVx[2:end-0,:]
    iVxNW     =  ones(Int64, size(NumVy)); iVxNW[:,1:end-1] = NumVx[1:end-1,:]
    iVxNE     =  ones(Int64, size(NumVy)); iVxNE[:,1:end-1] = NumVx[2:end-0,:]
    iPS       =  ones(Int64, size(NumVy)); iPS[:,2:end-0]   = NumP
    iPN       =  ones(Int64, size(NumVy)); iPN[:,1:end-1]   = NumP

    if BC.periodix==1
        iVyW[  1,:] = NumVy[end,:]
        iVyE[end,:] = NumVy[  1,:]
    end

    # Viscosity coefficients
    etaS      = zeros(size(NumVy)); etaS[:,2:end-0] = etac[:,1:end-0]
    etaN      = zeros(size(NumVy)); etaN[:,1:end-1] = etac[:,1:end-0]
    etaW      = zeros(size(NumVy)); etaW[1:end-0,:] = etav[1:end-1,:]
    etaE      = zeros(size(NumVy)); etaE[1:end-0,:] = etav[2:end-0,:]
    # Finite difference coefficients
    @tturbo cVyC  = -(-4/3*etaN./dy - 4/3*etaS./dy)./dy - (-1.0*etaE./dx - 1.0*etaW./dx)./dx
    @tturbo cVyW  = -1.0*etaW./dx.^2
    @tturbo cVyE  = -1.0*etaE./dx.^2
    @tturbo cVyS  = -4/3*etaS./dy.^2
    @tturbo cVyN  = -4/3*etaN./dy.^2
    @tturbo cVxSW = 2/3*etaS./(dx.*dy) - 1.0*etaW./(dx.*dy)
    @tturbo cVxSE = 1.0*etaE./(dx.*dy) - 2/3*etaS./(dx.*dy)
    @tturbo cVxNW = -2/3*etaN./(dx.*dy) + 1.0*etaW./(dx.*dy)
    @tturbo cVxNE = -1.0*etaE./(dx.*dy) + 2/3*etaN./(dx.*dy)
    DamS =  ones(Float64, size(NumVy)); DamS[:,2:end-0]   .= f.Damc
    DamN =  ones(Float64, size(NumVy)); DamN[:,1:end-1]   .= f.Damc
    @tturbo cPS   = -DamS./dy .* ones(size(NumVy)); cPS[:,  1] .= 0.0;  cPS[:,end] .= 0.0
    @tturbo cPN   =  DamN./dy .* ones(size(NumVy)); cPN[:,  1] .= 0.0;  cPN[:,end] .= 0.0

    if BC.periodix==0
        if BC.Vy.type_W==11
            cVyC[  1,:] .-= cVyW[  1,:]
        end
        if BC.Vy.type_W==22
            cVyC[  1,:] .+= cVyW[  1,:]
        end
        if BC.Vy.type_E==11
            cVyC[end,:] .-= cVyE[end,:]
        end
        if BC.Vy.type_E==22
            cVyC[end,:] .+= cVyE[end,:]
        end
        cVyW[  1,:] .= 0.0
        cVyE[end,:] .= 0.0
    end

    # N-S Dirichlet nodes
    cVyC[:,1]  .= DirScale; cVyC[:,end]  .= DirScale
    cVyW[:,1]  .= 0.0;      cVyW[:,end]  .= 0.0
    cVyE[:,1]  .= 0.0;      cVyE[:,end]  .= 0.0
    cVyS[:,1]  .= 0.0;      cVyS[:,end]  .= 0.0
    cVyN[:,1]  .= 0.0;      cVyN[:,end]  .= 0.0
    cVxSW[:,1] .= 0.0;      cVxSW[:,end] .= 0.0
    cVxNW[:,1] .= 0.0;      cVxNW[:,end] .= 0.0
    cVxSE[:,1] .= 0.0;      cVxSE[:,end] .= 0.0
    cVxNE[:,1] .= 0.0;      cVxNE[:,end] .= 0.0
    # Symmetry - kill Dirichlet connections
    cVyS[:,     2] .= 0.0
    cVyN[:, end-1] .= 0.0

    if BC.periodix==0
        # cVyW[  1,:] .= 0.0;
        # cVyE[end,:] .= 0.0;

        # cVyC[:,1] .= 1e3; cVyC[:,end] .= 1e3
        # cVyW[:,1] .= 0.0; cVyW[:,end] .= 0.0
        # cVyE[:,1] .= 0.0; cVyE[:,end] .= 0.0
        # cVyS[:,1] .= 0.0; cVyS[:,end] .= 0.0
        # cVyN[:,1] .= 0.0; cVyN[:,end] .= 0.0
        # cVxSW[:,1] .= 0.0; cVxSW[:,end] .= 0.0
        # cVxNW[:,1] .= 0.0; cVxNW[:,end] .= 0.0
        # cVxSE[:,1] .= 0.0; cVxSE[:,end] .= 0.0
        # cVxNE[:,1] .= 0.0; cVxNE[:,end] .= 0.0
        # # Symmetry - kill Dirichlet connections
        # cVyS[:,   2]   .= 0.0
        # cVyN[:,end-1]  .= 0.0
        # cVxSW[  1,  :] .= 0.0
        # cVxSE[end,  :] .= 0.0
        # cVxNW[  1,  :] .= 0.0
        # cVxNE[end,  :] .= 0.0
        cVyW[  1,:] .= 0.0;
        cVyE[end,:] .= 0.0;
        # Symmetry - kill Dirichlet connections
        cVxSW[  1,  :] .= 0.0
        cVxSE[end,  :] .= 0.0
        cVxNW[  1,  :] .= 0.0
        cVxNE[end,  :] .= 0.0
    end

    if BC.periodix==1
        # Remove redundant Vx equation on the right side to make Kuu matrix symmetric positive definite
        iVxC  = collect(iVxC[1:end-1,:]);  cVxC  = collect(cVxC[1:end-1,:])
        iVxW  = collect(iVxW[1:end-1,:]);  cVxW  = collect(cVxW[1:end-1,:])
        iVxE  = collect(iVxE[1:end-1,:]);  cVxE  = collect(cVxE[1:end-1,:])
        iVxS  = collect(iVxS[1:end-1,:]);  cVxS  = collect(cVxS[1:end-1,:])
        iVxN  = collect(iVxN[1:end-1,:]);  cVxN  = collect(cVxN[1:end-1,:])
        iVySW = collect(iVySW[1:end-1,:]); cVySW = collect(cVySW[1:end-1,:])
        iVySE = collect(iVySE[1:end-1,:]); cVySE = collect(cVySE[1:end-1,:])
        iVyNW = collect(iVyNW[1:end-1,:]); cVyNW = collect(cVyNW[1:end-1,:])
        iVyNE = collect(iVyNE[1:end-1,:]); cVyNE = collect(cVyNE[1:end-1,:])
        iPW   = collect(iPW[1:end-1,:]);   cPW   = collect(cPW[1:end-1,:])
        iPE   = collect(iPE[1:end-1,:]);   cPE   = collect(cPE[1:end-1,:])
        # ncx = size(cVxC,1)-1
        # ncy = size(cVxC,2)
        # iVxC1  = zeros(Int64,(ncx,ncy)); cVxC1  = zeros(Float64,(ncx,ncy))
        # iVxW1  = zeros(Int64,(ncx,ncy)); cVxW1  = zeros(Float64,(ncx,ncy))
        # iVxE1  = zeros(Int64,(ncx,ncy)); cVxE1  = zeros(Float64,(ncx,ncy))
        # iVxS1  = zeros(Int64,(ncx,ncy)); cVxS1  = zeros(Float64,(ncx,ncy))
        # iVxN1  = zeros(Int64,(ncx,ncy)); cVxN1  = zeros(Float64,(ncx,ncy));
        # iVySW1  = zeros(Int64,(ncx,ncy)); cVySW1  = zeros(Float64,(ncx,ncy));
        # iVySE1  = zeros(Int64,(ncx,ncy)); cVySE1  = zeros(Float64,(ncx,ncy));
        # iVyNW1  = zeros(Int64,(ncx,ncy)); cVyNW1  = zeros(Float64,(ncx,ncy));
        # iVyNE1  = zeros(Int64,(ncx,ncy)); cVyNE1  = zeros(Float64,(ncx,ncy));
        # iPW1  = zeros(Int64,(ncx,ncy)); cPW1  = zeros(Float64,(ncx,ncy));
        # iPE1  = zeros(Int64,(ncx,ncy)); cPE1  = zeros(Float64,(ncx,ncy));
        # iVxC1  .= iVxC[1:end-1,:];  cVxC1  .= cVxC[1:end-1,:]
        # iVxW1  .= iVxW[1:end-1,:];  cVxW1  .= cVxW[1:end-1,:]
        # iVxE1  .= iVxE[1:end-1,:];  cVxE1  .= cVxE[1:end-1,:]
        # iVxS1  .= iVxS[1:end-1,:];  cVxS1  .= cVxS[1:end-1,:]
        # iVxN1  .= iVxN[1:end-1,:];  cVxN1  .= cVxN[1:end-1,:]
        # iVySW1 .= iVySW[1:end-1,:]; cVySW1 .= cVySW[1:end-1,:]
        # iVySE1 .= iVySE[1:end-1,:]; cVySE1 .= cVySE[1:end-1,:]
        # iVyNW1 .= iVyNW[1:end-1,:]; cVyNW1 .= cVyNW[1:end-1,:]
        # iVyNE1 .= iVyNE[1:end-1,:]; cVyNE1 .= cVyNE[1:end-1,:]
        # iPW1   .= iPW[1:end-1,:];   cPW1   .= cPW[1:end-1,:]
        # iPE1   .= iPE[1:end-1,:];   cPE1   .= cPE[1:end-1,:]
    end

    # Sparse matrix Kuu
    nVx = size(cVxC,1)*size(cVxC,2)
    nVy = size(cVyC,1)*size(cVyC,2)
    I   = zeros(  Int64, 9*(nVx + nVy) )
    J   = zeros(  Int64, 9*(nVx + nVy) )
    V   = zeros(Float64, 9*(nVx + nVy) )

    #------------------- Vx
    FillCoefficients!(I, J, V, 0*nVx, iVxC[:], iVxC[:], cVxC[:])
    FillCoefficients!(I, J, V, 1*nVx, iVxC[:], iVxW[:], cVxW[:])
    FillCoefficients!(I, J, V, 2*nVx, iVxC[:], iVxE[:], cVxE[:])
    FillCoefficients!(I, J, V, 3*nVx, iVxC[:], iVxS[:], cVxS[:])
    FillCoefficients!(I, J, V, 4*nVx, iVxC[:], iVxN[:], cVxN[:])
    FillCoefficients!(I, J, V, 5*nVx, iVxC[:], iVySW[:], cVySW[:])
    FillCoefficients!(I, J, V, 6*nVx, iVxC[:], iVySE[:], cVySE[:])
    FillCoefficients!(I, J, V, 7*nVx, iVxC[:], iVyNW[:], cVyNW[:])
    FillCoefficients!(I, J, V, 8*nVx, iVxC[:], iVyNE[:], cVyNE[:])
    #------------------- Vy
    FillCoefficients!(I, J, V, 9*nVx+0*nVy, iVyC[:], iVyC[:], cVyC[:])
    FillCoefficients!(I, J, V, 9*nVx+1*nVy, iVyC[:], iVyW[:], cVyW[:])
    FillCoefficients!(I, J, V, 9*nVx+2*nVy, iVyC[:], iVyE[:], cVyE[:])
    FillCoefficients!(I, J, V, 9*nVx+3*nVy, iVyC[:], iVyS[:], cVyS[:])
    FillCoefficients!(I, J, V, 9*nVx+4*nVy, iVyC[:], iVyN[:], cVyN[:])
    FillCoefficients!(I, J, V, 9*nVx+5*nVy, iVyC[:], iVxSW[:], cVxSW[:])
    FillCoefficients!(I, J, V, 9*nVx+6*nVy, iVyC[:], iVxSE[:], cVxSE[:])
    FillCoefficients!(I, J, V, 9*nVx+7*nVy, iVyC[:], iVxNW[:], cVxNW[:])
    FillCoefficients!(I, J, V, 9*nVx+8*nVy, iVyC[:], iVxNE[:], cVxNE[:])
    #------------------- Assemble
    Kuu = sparse( I, J, V)
    droptol!(Kuu, 1e-9)
    # display(UnicodePlots.spy(K))

    # Sparse matrix Kup
    I   = zeros(  Int64, 2*(nVx + nVy) )
    J   = zeros(  Int64, 2*(nVx + nVy) )
    V   = zeros(Float64, 2*(nVx + nVy) )
    FillCoefficients!(I, J, V, 0*nVx+0*nVy, iVxC[:], iPW[:], cPW[:])
    FillCoefficients!(I, J, V, 1*nVx+0*nVy, iVxC[:], iPE[:], cPE[:])
    FillCoefficients!(I, J, V, 2*nVx+0*nVy, iVyC[:], iPS[:], cPS[:])
    FillCoefficients!(I, J, V, 2*nVx+1*nVy, iVyC[:], iPN[:], cPN[:])
    Kup = sparse( I, J, V)
    droptol!(Kup, 1e-9)

    # Sparse matrix Kpu: Coefficients
    iVxW  =  ones(  Int64, size(NumP)); iVxW .= NumVx[1:end-1,:]
    iVxE  =  ones(  Int64, size(NumP)); iVxE .= NumVx[2:end-0,:]
    iVyS  =  ones(  Int64, size(NumP)); iVyS .= NumVy[:,1:end-1]
    iVyN  =  ones(  Int64, size(NumP)); iVyN .= NumVy[:,2:end-0]
    cVxW  = zeros(Float64, size(NumP)); cVxW .= -1.0 ./ dx
    cVxE  = zeros(Float64, size(NumP)); cVxE .=  1.0 ./ dx
    cVyS  = zeros(Float64, size(NumP)); cVyS .= -1.0 ./ dy
    cVyN  = zeros(Float64, size(NumP)); cVyN .=  1.0 ./ dy

    # Kill Dirichlet connections
    cVyS[:,1] .= 0.0; cVyN[:,end] .= 0.0
    if BC.periodix==0
        cVxW[1,:] .= 0.0; cVxE[end,:] .= 0.0
    end
    # Sparse matrix Kpu: Triplets
    nP = length(NumP)
    I  = zeros(  Int64, 4*nP)
    J  = zeros(  Int64, 4*nP)
    V  = zeros(Float64, 4*nP)
    FillCoefficients!(I, J, V, 0*nP, NumP[:], iVxW[:], cVxW[:])
    FillCoefficients!(I, J, V, 1*nP, NumP[:], iVxE[:], cVxE[:])
    FillCoefficients!(I, J, V, 2*nP, NumP[:], iVyS[:], cVyS[:])
    FillCoefficients!(I, J, V, 3*nP, NumP[:], iVyN[:], cVyN[:])
    Kpu = sparse( I, J, V)
    droptol!(Kpu, 1e-9)
    # display(UnicodePlots.spy(Div))
    return Kuu, Kup, Kpu
end
export StokesAssembly

########

@views function SetInitialVelocity!( Vx::Matrix{Float64}, Vy::Matrix{Float64}, BC::BoundaryConditions, xv::LinRange{Float64}, yv::LinRange{Float64}, xce::LinRange{Float64}, yce::LinRange{Float64}, xmin::Float64, xmax::Float64, ymin::Float64, ymax::Float64, ncx::Int64, ncy::Int64 )

    # Initial guess - Vx
    Vx[1    ,:] .= BC.Vx.Dir_W
    Vx[end  ,:] .= BC.Vx.Dir_E
    Vx[:,    2] .= BC.Vx.Dir_S
    Vx[:,end-1] .= BC.Vx.Dir_N
    Vy[2    ,:] .= BC.Vy.Dir_W
    Vy[end-1,:] .= BC.Vy.Dir_E
    Vy[:,    1] .= BC.Vy.Dir_S
    Vy[:,  end] .= BC.Vy.Dir_N

    Vx_ini = zeros(ncx+1,ncy+2)
    Vx_ini .= Vx
    Vy_ini = zeros(ncx+2,ncy+1)
    Vy_ini .= Vy

    for i=1:size(Vx,1)
        for j=1:size(Vx,2)
            wW = 1.0 - (xv[i]-xmin)/(xmax-xmin)
            Vx[i,j] = wW * BC.Vx.Dir_W[j] + (1.0-wW) * BC.Vx.Dir_E[j]
            if i>1 && i<size(Vx,1)
                wS = 1.0 - (yce[j]-ymin)/(ymax-ymin)
                Vx[i,j] += wS * BC.Vx.Dir_S[i] + (1.0-wS) * BC.Vx.Dir_N[i]
            end
        end
    end

    # Initial guess - Vy
    for i=1:size(Vy,1)
        for j=1:size(Vy,2)
            wS = 1.0 - (yv[j]-ymin)/(ymax-ymin)
            Vy[i,j] = wS * BC.Vy.Dir_S[i] + (1.0-wS) * BC.Vy.Dir_N[i]
            if j>1 && j<size(Vy,2)
                wW = 1.0 - (xce[i]-xmin)/(xmax-xmin)
                Vy[i,j] += wW * BC.Vy.Dir_W[j] + (1.0-wW) * BC.Vy.Dir_E[j]
            end
        end
    end
    return
end
export SetInitialVelocity!

########

@views function SetInitialVelocity2!( Vx::Matrix{Float64}, Vy::Matrix{Float64}, BC::BoundaryConditions, dom::ModelDomain )

    ncx, ncy = dom.ncx, dom.ncy

    # Initial guess - Vx
    Vx[1    ,:] .= BC.Vx.Dir_W
    Vx[end  ,:] .= BC.Vx.Dir_E
    Vx[:,    2] .= BC.Vx.Dir_S
    Vx[:,end-1] .= BC.Vx.Dir_N
    Vy[2    ,:] .= BC.Vy.Dir_W
    Vy[end-1,:] .= BC.Vy.Dir_E
    Vy[:,    1] .= BC.Vy.Dir_S
    Vy[:,  end] .= BC.Vy.Dir_N

    Vx_ini = zeros(ncx+1,ncy+2)
    Vx_ini .= Vx
    Vy_ini = zeros(ncx+2,ncy+1)
    Vy_ini .= Vy

    for i=1:size(Vx,1)
        for j=1:size(Vx,2)
            wW = 1.0 - (dom.xv[i]-dom.xmin)/(dom.xmax-dom.xmin)
            Vx[i,j] = wW * BC.Vx.Dir_W[j] + (1.0-wW) * BC.Vx.Dir_E[j]
            if i>1 && i<size(Vx,1)
                wS = 1.0 - (dom.yce[j]-dom.ymin)/(dom.ymax-dom.ymin)
                Vx[i,j] += wS * BC.Vx.Dir_S[i] + (1.0-wS) * BC.Vx.Dir_N[i]
            end
        end
    end

    # Initial guess - Vy
    for i=1:size(Vy,1)
        for j=1:size(Vy,2)
            wS = 1.0 - (dom.yv[j]-dom.ymin)/(dom.ymax-dom.ymin)
            Vy[i,j] = wS * BC.Vy.Dir_S[i] + (1.0-wS) * BC.Vy.Dir_N[i]
            if j>1 && j<size(Vy,2)
                wW = 1.0 - (dom.xce[i]-dom.xmin)/(dom.xmax-dom.xmin)
                Vy[i,j] += wW * BC.Vy.Dir_W[j] + (1.0-wW) * BC.Vy.Dir_E[j]
            end
        end
    end
    return
end
export SetInitialVelocity2!


text:
fortran_double(s) = parse(Float64, replace(s, "D", "E"))

function parse_zf_block(block)
    lower = split(strip(block[1][5:end]))
    upper = split(strip(block[2][5:end]))
    len = split(block[4])[[4,7,10]]
    vel = split(block[5])
    Dict(:lower => lower[2],
         :upper => upper[2],
         :E1 => parse(Float64, lower[1]),
         :E2 => parse(Float64, upper[1]),
         :SL => fortran_double(len[1]),
         :SV => fortran_double(vel[1]),
         :GFL => fortran_double(len[2]),
         :GFV => fortran_double(vel[2]),
         :AKIL => fortran_double(len[3]),
         :AKIV => fortran_double(vel[3]))
end

function read_zf_res(directory::AbstractString,
                     name::AbstractString = "zf_res",
                     dd::AbstractString = "")
    blocks = []
    open("$(directory)/$(name)") do file
        readline(file)
        block = []
        for line in eachline(file)
            if strip(line) == ""
                push!(blocks, parse_zf_block(block))
                block = []
            else
                push!(block, line)
            end
        end
        push!(blocks, parse_zf_block(block))
    end

    if dd != ""
        mat_els = readdlm(dd)
        m,n = maximum(mat_els[:,1]),maximum(mat_els[:,2])
        m*n == length(blocks) || error("Incongruent transition matrix size")
        L = [fortran_double(s) for s in mat_els[:,3]]
        V = [fortran_double(s) for s in mat_els[:,4]]
        for i in eachindex(blocks)
            blocks[i][:DL] = L[i]
            blocks[i][:DV] = V[i]
        end
    end

    blocks
end

read_zf_res(directory::AbstractString,
            i::Integer, j::Integer) =
                read_zf_res(directory,
                            @sprintf("trans.%d_%d", i, j),
                            @sprintf("dd.%03d_%03d", i, j))

export read_zf_res

ext: jl


text:
using Test
using AutomotiveDrivingModels
using Distributions

include("test_roadways.jl")
include("test_agent_definitions.jl")
include("test_states.jl")
include("test_collision_checkers.jl")
include("test_actions.jl")
include("test_features.jl")
include("test_behaviors.jl")
include("test_simulation.jl")


text:
import CUDAnative

struct Single{T}
    u::T
end

import Base: promote_rule, convert
import Base.+, Base.-, Base.*, Base./, Base.^
import Base.>, Base.<, Base.>=, Base.<=, Base.==, Base.!=
import Base.isapprox
import Base.zero, Base.one

zero(::Single{T}) where T = Single(zero(T))
one(::Single{T}) where T = Single(one(T))
convert(::Type{Single{T}}, c::T) where T<:Number = Single(c)
promote_rule(::Type{Single{T}}, ::Type{T}) where T<:Number = Single{T}

(+)(x::Single) = Single(x.u)
(+)(x::Single{T}, y::Single{T}) where T = Single(x.u+y.u)
(+)(x::Single{T}, c::Number) where T = (+)(promote(x,T(c))...)
(+)(c::Number, y::Single{T}) where T = (+)(promote(T(c),y)...)

(-)(x::Single) = Single(-x.u)
(-)(x::Single{T}, y::Single{T}) where T = Single(x.u-y.u)
(-)(x::Single{T}, c::Number) where T = (-)(promote(x,T(c))...)
(-)(c::Number, y::Single{T}) where T = (-)(promote(T(c),y)...)

(*)(x::Single{T}, y::Single{T}) where T = Single(x.u*y.u)
(*)(x::Single{T}, c::Number) where T = (*)(promote(x,T(c))...)
(*)(c::Number, y::Single{T}) where T = (*)(promote(T(c),y)...)

(/)(x::Single{T}, y::Single{T}) where T = Single(x.u/y.u)
(/)(x::Single{T}, c::Number) where T = (/)(promote(x,T(c))...)
(/)(c::Number, y::Single{T}) where T = (/)(promote(T(c),y)...)

CUDAnative.sin(x::Single) = Single(CUDAnative.sin(x.u))
CUDAnative.sin_fast(x::Single) = Single(CUDAnative.sin_fast(x.u))
CUDAnative.sinpi(x::Single) = Single(CUDAnative.sinpi(x.u))

CUDAnative.cos(x::Single) = Single(CUDAnative.cos(x.u))
CUDAnative.cos_fast(x::Single) = Single(CUDAnative.cos_fast(x.u))
CUDAnative.cospi(x::Single) = Single(CUDAnative.cospi(x.u))

CUDAnative.tan(x::Single) = Single(CUDAnative.tan(x.u))
CUDAnative.tan_fast(x::Single) = Single(CUDAnative.tan_fast(x.u))

CUDAnative.acos(x::Single) = Single(CUDAnative.acos(x.u))
CUDAnative.asin(x::Single) = Single(CUDAnative.asin(x.u))
CUDAnative.atan(x::Single) = Single(CUDAnative.atan(x.u))

CUDAnative.atan2(x::Single, y::Single) = Single(CUDAnative.atan2(x.u,y.u))
CUDAnative.atan2(x::Single{T}, c::Number) where T = CUDAnative.atan2(promote(x,T(c))...)
CUDAnative.atan2(c::Number, x::Single{T}) where T = CUDAnative.atan2(promote(T(c),x)...)

CUDAnative.cosh(x::Single) = Single(CUDAnative.cosh(x.u))
CUDAnative.sinh(x::Single) = Single(CUDAnative.sinh(x.u))
CUDAnative.tanh(x::Single) = Single(CUDAnative.tanh(x.u))

CUDAnative.acosh(x::Single) = Single(CUDAnative.acosh(x.u))
CUDAnative.asinh(x::Single) = Single(CUDAnative.asinh(x.u))
CUDAnative.atanh(x::Single) = Single(CUDAnative.atanh(x.u))

CUDAnative.exp(x::Single) = Single(CUDAnative.exp(x.u))
CUDAnative.exp_fast(x::Single) = Single(CUDAnative.exp_fast(x.u))
CUDAnative.exp2(x::Single) = Single(CUDAnative.exp2(x.u))
CUDAnative.exp10(x::Single) = Single(CUDAnative.exp10(x.u))
CUDAnative.expm1(x::Single) = Single(CUDAnative.expm1(x.u))

CUDAnative.log(x::Single) = Single(CUDAnative.log(x.u))
CUDAnative.log_fast(x::Single) = Single(CUDAnative.log_fast(x.u))
CUDAnative.log10(x::Single) = Single(CUDAnative.log10(x.u))
CUDAnative.log10_fast(x::Single) = Single(CUDAnative.log10_fast(x.u))
CUDAnative.log2(x::Single) = Single(CUDAnative.log2(x.u))
CUDAnative.log2_fast(x::Single) = Single(CUDAnative.log2_fast(x.u))
CUDAnative.log1p(x::Single) = Single(CUDAnative.log1p(x.u))

CUDAnative.abs(x::Single) = Single(CUDAnative.abs(x.u))
CUDAnative.saturate(x::Single) = Single(CUDAnative.saturate(x.u))
CUDAnative.sqrt(x::Single) = Single(CUDAnative.sqrt(x.u))
CUDAnative.rsqrt(x::Single) = Single(CUDAnative.rsqrt(x.u))
CUDAnative.cbrt(x::Single) = Single(CUDAnative.cbrt(x.u))
CUDAnative.rcbrt(x::Single) = Single(CUDAnative.rcbrt(x.u))

CUDAnative.hypot(x::Single{T}, y::Single{T}) where T = Single(CUDAnative.hypot(x.u, y.u))
CUDAnative.hypot(x::Single{T}, c::Number) where T = CUDAnative.hypot(promote(x,T(c))...)
CUDAnative.hypot(c::Number, x::Single{T}) where T = CUDAnative.hypot(promote(T(c),x)...)

CUDAnative.erf(x::Single) = Single(CUDAnative.erf(x.u))
CUDAnative.erfc(x::Single) = Single(CUDAnative.erfc(x.u))
CUDAnative.erfinv(x::Single) = Single(CUDAnative.erfinv(x.u))
CUDAnative.erfcinv(x::Single) = Single(CUDAnative.erfcinv(x.u))
CUDAnative.pow(x::Single, p) = Single(CUDAnative.pow(x.u, p))
(^)(x::Single, p) = Single(x.u ^ p)

#############################################################################

CUDAnative.ceil(x::Single) = Single(CUDAnative.ceil(x.u))
CUDAnative.floor(x::Single) = Single(CUDAnative.floor(x.u))

CUDAnative.min(x::Single{T}, y::Single{T}) where T = CUDAnative.min(x.u, y.u)
CUDAnative.min(x::Single{T}, c::Number) where T = CUDAnative.min(x.u, c)
CUDAnative.min(c::Number, x::Single{T}) where T = CUDAnative.min(c, x.u)
CUDAnative.max(x::Single{T}, y::Single{T}) where T = CUDAnative.max(x.u, y.u)
CUDAnative.max(x::Single{T}, c::Number) where T = CUDAnative.max(x.u, c)
CUDAnative.max(c::Number, x::Single{T}) where T = CUDAnative.max(c, x.u)
CUDAnative.isfinite(x::Single) = CUDAnative.isfinite(x.u)
CUDAnative.isinf(x::Single) = CUDAnative.isinf(x.u)
CUDAnative.isnan(x::Single) = CUDAnative.isnan(x.u)
CUDAnative.nearbyint(x::Single) = CUDAnative.nearbyint(x.u)
CUDAnative.nextafter(x::Single) = CUDAnative.nextafter(x.u)
CUDAnative.signbit(x::Single) = CUDAnative.signbit(x.u)

(<)(x::Single{T}, y::Single{T}) where T = (x.u < y.u)
(<)(x::Single{T}, c::Number) where T = (x.u < c)
(<)(c::Number, x::Single{T}) where T = (c < x.u)

(>)(x::Single{T}, y::Single{T}) where T = (x.u > y.u)
(>)(x::Single{T}, c::Number) where T = (x.u > c)
(>)(c::Number, x::Single{T}) where T = (c > x.u)

(<=)(x::Single{T}, y::Single{T}) where T = (x.u <= y.u)
(<=)(x::Single{T}, c::Number) where T = (x.u <= c)
(<=)(c::Number, x::Single{T}) where T = (c <= x.u)

(>=)(x::Single{T}, y::Single{T}) where T = (x.u >= y.u)
(>=)(x::Single{T}, c::Number) where T = (x.u >= c)
(>=)(c::Number, x::Single{T}) where T = (c >= x.u)

(==)(x::Single{T}, y::Single{T}) where T = (x.u == y.u)
(==)(x::Single{T}, c::Number) where T = (x.u == c)
(==)(c::Number, x::Single{T}) where T = (c == x.u)

(!=)(x::Single{T}, y::Single{T}) where T = (x.u != y.u)
(!=)(x::Single{T}, c::Number) where T = (x.u != c)
(!=)(c::Number, x::Single{T}) where T = (c != x.u)

isapprox(x::Single{T}, y::Single{T}) where T = isapprox(x.u, y.u)
isapprox(x::Single{T}, c::Number) where T = isapprox(x.u, c)
isapprox(c::Number, x::Single{T}) where T = isapprox(c, x.u)

#############################################################################


text:
module WinEncoding
#     const cp950=include("cp950.jl")
    const cp1252b=include("cp1252b.jl")
#     const cp936=include("cp936.jl")
#     const cp932=include("cp932.jl")
#     const cp932b=include("cp932b.jl")
#     const _x80=[0xc2,0x80]             # '\u80'
#     const _x80_936=[0xe2,0x82,0xac]     # '\u20ac' '€'
#     const _xff=[0xef, 0xa3, 0xb8]      # '\uf8f8'
#     const _xff_936=[0xef, 0xa3, 0xb5]      # '\uf8f5'
    const _invalid=[0xef, 0xbf, 0xbd]  # '�'
    export decode1252, decode950, encode1252

    """
        decode950(a::Vector{UInt8})
    Convert an array of bytes `a` representing text in encoding `cp950/big5/hkscs` to a string.
    - fallback to big5-hkscs or '\ufffd';no invalid sequence error
    - non-blocking

        decode950(f::Filename)
    Convert file `f` content to Vector{UInt8} from cp950 to utf-8 Vector{UInt8}
    ## Examples
    ```jldoctest
    julia> decode950([0xa4,0x48])
    "人"
    ```
    """
    const decode950(a::Vector{UInt8})=Cp950.decode(a)
    function decode950(f::String,type=Vector{UInt8})
        @assert isfile(f) "$f must be a valid file name"
        bom=read(f,3)
        if bom == [0xef,0xbb,0xbf]
            @info "Detect utf-8 bom(byte order mark) ; just read"
            read(f)
        elseif length(bom) >= 2 && bom[1:2]==[0xff,0xfe]
            @info "Detect UTF-16LE bom ; using built-in transcode"
            transcode(UInt8,reinterpret(read(f),UInt16))
        else
            Cp950.decode(read(f),type)
        end
    end
    const decode936(a)=Cp936.decode(a)
    const decode932(a)=Cp932.decode(a)

    module Cp950
        import .._invalid
        const cp950=include("cp950.jl")
        const _x80=[0xc2,0x80]             # '\u80'
        const _xff=[0xef, 0xa3, 0xb8]      # '\uf8f8'

        function decode(ss::Vector{UInt8},type::Type=String)
            blength=length(ss)
            o=Array{UInt8}(undef,blength*2+1)
            b=1
            skip1=false
            @inbounds for (i,c) in enumerate(ss)
                skip1==true && (skip1=false;continue)
                c<0x80 && (o[b]=c; b+=1;continue)
                c==0x80 && (copyto!(o,b,_x80,1,2); b+=2;continue)
                c==0xff && (copyto!(o,b,_xff,1,3); b+=3;continue)
                i % 1_000_000 == 1 && yield()
                skip1=true
                hh=c-0x80
                ll=((i1=i+1) <= blength ) ? ss[i1] : 0x00
                cc=(ll < 0x40 || ll==0xff || 0x7e<ll<0xa1) ? _invalid : (ll-=0x3f;ll > 0x3f && (ll-=0x22);cp950[hh][ll])
                ccl=length(cc)
                unsafe_copyto!(o,b,cc,1,ccl)
                b+=ccl
            end
            type===String ? String(@view(o[1:b-1])) : o[1:b-1]
        end
    end

    module Cp936
        import .._invalid
        const cp936=include("cp936.jl")
        const _x80=[0xe2,0x82,0xac]     # '\u20ac' '€'
        const _xff=[0xef, 0xa3, 0xb5]      # '\uf8f5'

        """
            decode936(a::Vector{UInt8})
        Convert an array of bytes `a` representing text in encoding ***cp936/gbk*** to a string.
        - fallback to '\ufffd';no invalid sequence error
        - non-blocking
        """
        function decode(ss::Vector{UInt8})
            blength=length(ss)
            o=Array{UInt8}(undef,blength*2+1)
            b=1
            skip1=false
            @inbounds for (i,c) in enumerate(ss)
                skip1==true && (skip1=false;continue)
                c<0x80 && (o[b]=c; b+=1;continue)
                c==0x80 && (copyto!(o,b,_x80,1,3); b+=3;continue)
                c==0xff && (copyto!(o,b,_xff,1,3); b+=3;continue)
                i % 1024*1024 == 0 && yield()
                skip1=true
                hh=c-0x80
                ll=((i1=i+1) <= blength) ? ss[i1] : 0x00
                cc=(ll < 0x40 || ll==0xff) ? _invalid : cp936[hh][ll-0x3f]
                ccl=length(cc)
                unsafe_copyto!(o,b,cc,1,ccl)
                b+=ccl
            end
            @view(o[1:b-1]) |> String
        end
    end #cp936

    module Cp932
        const _invalid=[0xe3, 0x83, 0xbb] #'・'
        const cp932=include("cp932.jl")
        const cp932b=include("cp932b.jl")
        """
            decode932(a::Vector{UInt8})
        Convert an array of bytes `a` representing text in encoding ***cp932/sjis*** to a string.
        - fallback to '・';no invalid sequence error
        - non-blocking
        """
        function decode(ss::Vector{UInt8})
            blength=length(ss)
            o=Array{UInt8}(undef,blength*2+1)
            b=1
            skip1=false
            @inbounds for (i, c) in enumerate(ss)
                skip1==true && (skip1=false;continue)
                c<0x80 && (o[b]=c; b+=1;continue)
                if (c == 0x80 || 0xa0 <= c <0xe0 || c > 0xef)
                    let cc=cp932b[c-0x7f]; ccl=length(cc)
                        copyto!(o,b,cc,1,ccl); b+=ccl;continue
                    end
                end
                i % 1024*1024 == 0 && yield()
                skip1=true
                hh=c-0x80
                hh > 0x1f && (hh-=0x40)
                ll=((i1=i+1) <= blength) ? ss[i1] : 0x00
                cc=(ll < 0x40 || ll > 0xfc) ? _invalid : cp932[hh][ll-0x3f]
                ccl=length(cc)
                unsafe_copyto!(o,b,cc,1,ccl)
                b+=ccl
            end
            @view(o[1:b-1]) |> String
        end
    end #cp932

    #internal use
    function _7bit(ss)
        for c in ss
            c < 0x80 || return false
        end
        true
    end
    #internal use
    function _7bit2(ss)
        for (i,c) in enumerate(ss)
            c < 0x80 || return (false,i-one(i))
        end
        true,length(ss)
    end
    """
        decode1252(a::Vector{UInt8})
    Convert an array of bytes `a` representing text in encoding `cp1252/windows-1252` to a string.
    - [0x8d] => '\\u8d' as Windows API does
    - no invalid sequence error
    - non-blocking
    ## Examples
    ```jldoctest
    julia> decode1252([0x80])
    "€"

    julia> decode1252([0xa9])
    "©"

    ```
    """
    function decode1252(ss::Vector{UInt8},type::Type=String)
        o=Array{UInt8}(undef,length(ss)*3)
        is7bit,len = _7bit2(ss)
#         o[1:len]=ss[1:len]
#         unsafe_copyto!(o,1,ss,1,len)
        copyto!(o,1,ss,1,len)
        b=len+1
        @inbounds for c in @view(ss[b:end])
            if c<0x80
                o[b]=c
                b+=1
            else
                b % 1024*1024 == 0 && yield()
                cc=cp1252b[c+1]
                ccl=length(cc)
                unsafe_copyto!(o,b,cc,1,ccl)
                b+=ccl
            end
        end
        type===String ? String(@view(o[1:b-1])) : o[1:b-1]
    end
    function decode1252(f::String,type=Vector{UInt8})
        @assert isfile(f) "Input must be a valid file name"
        decode1252(read(f),type)
    end
    function encode1252(ss::String)
        tbl=Dict('Ú' => 0xda, '\uad' => 0xad, '¢' => 0xa2, 'Í' => 0xcd, 'ç' => 0xe7, 'ÿ' => 0xff, '®' => 0xae, 'ã' => 0xe3, '«' => 0xab, '¼' => 0xbc, 'å' => 0xe5, 'æ' => 0xe6, '€' => 0x80, '†' => 0x86, 'ä' => 0xe4, 'À' => 0xc0, 'ž' => 0x9e, 'é' => 0xe9, 'ò' => 0xf2, '°' => 0xb0, 'Š' => 0x8a, '¥' => 0xa5, '£' => 0xa3, '¦' => 0xa6, 'Á' => 0xc1, '“' => 0x93, 'Þ' => 0xde, '¿' => 0xbf, 'Å' => 0xc5, '\u8d' => 0x8d, 'Ã' => 0xc3, '¨' => 0xa8, 'œ' => 0x9c, '§' => 0xa7, 'á' => 0xe1, '\u81' => 0x81, '´' => 0xb4, 'Ô' => 0xd4, 'Ý' => 0xdd, 'à' => 0xe0, 'Ò' => 0xd2, 'ß' => 0xdf, 'ê' => 0xea, '²' => 0xb2, 'ù' => 0xf9, 'û' => 0xfb, 'þ' => 0xfe, '„' => 0x84, '©' => 0xa9, '‰' => 0x89, 'ð' => 0xf0, 'Ö' => 0xd6, '¬' => 0xac, '‡' => 0x87, 'ü' => 0xfc, '‘' => 0x91, 'ª' => 0xaa, 'ý' => 0xfd, '\u9d' => 0x9d, 'Õ' => 0xd5, 'Ž' => 0x8e, 'Ø' => 0xd8, 'Ù' => 0xd9, '¡' => 0xa1, '¸' => 0xb8, '›' => 0x9b, 'Ê' => 0xca, '×' => 0xd7, 'è' => 0xe8, '¹' => 0xb9, '…' => 0x85, '\ua0' => 0xa0, 'ˆ' => 0x88, '–' => 0x96, 'Ì' => 0xcc, 'ñ' => 0xf1, '÷' => 0xf7, '¾' => 0xbe, 'õ' => 0xf5, '½' => 0xbd, 'Û' => 0xdb, 'ô' => 0xf4, 'ó' => 0xf3, '»' => 0xbb, 'ö' => 0xf6, 'Ï' => 0xcf, '±' => 0xb1, 'Œ' => 0x8c, 'Â' => 0xc2, 'ƒ' => 0x83, 'Î' => 0xce, 'Ü' => 0xdc, '¤' => 0xa4, 'Ä' => 0xc4, '˜' => 0x98, 'µ' => 0xb5, 'Ð' => 0xd0, '³' => 0xb3, 'Ñ' => 0xd1, 'Ÿ' => 0x9f, '’' => 0x92, 'š' => 0x9a, 'ì' => 0xec, '—' => 0x97, 'ë' => 0xeb, 'î' => 0xee, '¶' => 0xb6, 'ø' => 0xf8, 'Ó' => 0xd3, 'Æ' => 0xc6, 'ú' => 0xfa, 'Ç' => 0xc7, '·' => 0xb7, '\u8f' => 0x8f, '™' => 0x99, 'ï' => 0xef, 'í' => 0xed, 'â' => 0xe2, 'È' => 0xc8, 'É' => 0xc9, '‚' => 0x82, '‹' => 0x8b, 'Ë' => 0xcb, 'º' => 0xba, '”' => 0x94, '¯' => 0xaf, '\u90' => 0x90, '•' => 0x95)
        chars=collect(ss)
        o=Array{UInt8}(undef,ncodeunits(ss))
        b=1
        for c in chars
            cc= codepoint(c)
            o[b]= cc < 0x00000080 ? UInt8(cc) : get(tbl,c,0x3f)
            b+=1
        end
        o[1:b-1]
    end #encode1252
end # module

max_stars_repo_name: testercwt/WinEncoding

max_issues_repo_name: testercwt/WinEncoding

max_forks_repo_name: testercwt/WinEncoding


text:
import AdversarialPrediction: CM_TP, CM_FP, CM_FN, CM_TN, CM_AP, CM_AN, CM_PP, CM_PN, CM_ALL

tp = CM_TP()
fp = CM_FP()
fn = CM_FN()
tn = CM_TN()
ap = CM_AP()
an = CM_AN()
pp = CM_PP()
pn = CM_PN()
al = CM_ALL()

@testset "expression" begin
    @test string(2 * tp) == "2.0 TP"
    @test string(tp + fp) == "PP"
    @test string(ap + an) == "ALL"
    @test string(ap + pp) == "(AP) + (PP)"
    @test string(2 * tp / (ap + pp)) == "(2.0 TP) / ((AP) + (PP))"

    kappa = ((tp + tn) / al - (ap * pp + an * pn) / al^2) / (1 - (ap * pp + an * pn) / al^2)
    str_kappa = "((((TP) + (TN)) / (ALL)) + (-((((AP) * (PP)) + ((AN) * (PN))) / ((ALL)^2.0)))) / ((1.0) + (-((((AP) * (PP)) + ((AN) * (PN))) / ((ALL)^2.0))))"
    @test string(kappa) == str_kappa
end

@testset "expression info" begin
    acc = (tp + tn) / (ap + an)
    @test acc.info.needs_adv_sum_marg == false

    prec = tp / pp
    @test prec.info.needs_pred_sum_marg == true
    @test prec.info.needs_adv_sum_marg == false

    gm = tp / (sqrt(pp * ap))
    @test gm.info.needs_adv_sum_marg == true
    @test gm.info.needs_pred_sum_marg == true
    @test gm.info.is_linear_tp_tn == true

    kappa = ((tp + tn) / al - (ap * pp + an * pn) / al^2) / (1 - (ap * pp + an * pn) / al^2)
    @test kappa.info.needs_pred_sum_marg == true
    @test kappa.info.needs_adv_sum_marg == true
    @test kappa.info.is_linear_tp_tn == true

    gm_err = tp / (sqrt(pp * tn))
    @test gm_err.info.needs_pred_sum_marg == true
    @test gm_err.info.needs_adv_sum_marg == false
    @test gm_err.info.is_linear_tp_tn == false
end

@testset "invalid metric" begin
    @metric Metric1
    function define(::Type{Metric1}, C::ConfusionMatrix)
        return C.tp * C.tn / sqrt(C.ap * C.pp)
    end
    @metric Metric2
    function define(::Type{Metric2}, C::ConfusionMatrix)
        return C.ap / C.all
    end
    @metric Metric3
    function define(::Type{Metric3}, C::ConfusionMatrix)
        return (C.tp / C.all) * (1 - C.fp) / (C.pp)
    end
    @metric Metric4
    function define(::Type{Metric4}, C::ConfusionMatrix)
        return C.tp / (C.ap + C.pp + C.fp)
    end

    @test_throws ErrorException Metric1()
    @test_throws ErrorException Metric2()
    @test_throws ErrorException Metric3()
    @test_throws ErrorException Metric4()
end


text:
module ShellInterop
export shell, configure_shell,
       eval_sh, read_sh, end_sh,
       incomplete_msg_sh, err_sh, eval_file_sh,
       FXCRD_PATH

const FXCRD_PATH = chomp(readstring(`fxcrd_path`))

struct Shell{EvalIO<:IO,ReadIO<:IO}
  eval::EvalIO
  read::ReadIO
end

mutable struct ShellConfig
  eval::String
  read::String
end
const SHELL_CONFIG = ShellConfig("","")

function configure_shell(ev,rd)
  SHELL_CONFIG.eval = ev
  SHELL_CONFIG.read = rd
end

function shell(f)
  open("$(SHELL_CONFIG.eval)", "w") do eval_io
    open("$(SHELL_CONFIG.read)", "r") do read_io
      f(Shell(eval_io, read_io))
    end
  end
end

function sh_fxcrd_signal(name)
  "\$\${\$fxcrd_$name\$}\$\$"
end

function sh_send_fxcrd_signal(sh, name)
  println(sh.eval, sh_fxcrd_signal(name))
  flush(sh.eval)
end

function err_sh(sh)
  sh_send_fxcrd_signal(sh, "error")
end

function end_sh(sh)
  sh_send_fxcrd_signal(sh, "end")
end

function read_sh(sh, cmd)
  println(sh.eval, cmd)
  sh_send_fxcrd_signal(sh, "read")

  line=""
  result=""
  while line != sh_fxcrd_signal("return")
    result="$result\n$line"
    line=readline(sh.read)
  end
  return lstrip(result, '\n')
end

function clear_sh(sh)
  println(sh.eval, "")
end

function incomplete_msg_sh(sh, msg)
  println(sh.eval, msg)
end

function eval_sh(sh, cmd)
  println(sh.eval, cmd)
  sh_send_fxcrd_signal(sh, "eval")
end

function eval_file_sh(sh, file)
  println(sh.eval, ". '$file'")
  sh_send_fxcrd_signal(sh, "eval")
end

end

max_stars_repo_name: fcard/Fxcrd

max_issues_repo_name: fcard/Fxcrd

max_forks_repo_name: fcard/Fxcrd


text:
using StateSpaceEcon
using TimeSeriesEcon
using ModelBaseEcon

@info "Loading Model examples. Might take some time to pre-compile."
@using_example E1
@using_example E2
@using_example E3
@using_example E3nl
@using_example E6
@using_example E7
@using_example E7A

using Test
using Suppressor

@testset "1dsolvers" begin
    # f(x) = (x-2)*(x-3) = a x^2 + b x + c with vals = [a, x, b, c]
    let f(v) = v[1] * v[2]^2 + v[3] * v[2] + v[4],
        fdf(v) = (f(v), [v[2]^2, v[1] * v[2] * 2 + v[3], v[2], 1.0]),
        vals = [1.0, NaN, -5.0, 6.0]

        vals[2] = 0.0
        @test StateSpaceEcon.SteadyStateSolver.newton1!(fdf, vals, 2; tol=eps(), maxiter=8)
        @test vals ≈ [1.0, 2.0, -5.0, 6.0] atol = 1e3 * eps()

        vals[2] = 6.0
        @test StateSpaceEcon.SteadyStateSolver.newton1!(fdf, vals, 2; tol=eps(), maxiter=8)
        @test vals ≈ [1.0, 3.0, -5.0, 6.0] atol = 1e3 * eps()

        vals[2] = 0.0
        @test StateSpaceEcon.SteadyStateSolver.bisect!(f, vals, 2, fdf(vals)[2][2]; tol=eps())
        @test vals ≈ [1.0, 2.0, -5.0, 6.0] atol = 1e3 * eps()

        vals[2] = 6.0
        @test StateSpaceEcon.SteadyStateSolver.bisect!(f, vals, 2, fdf(vals)[2][2]; tol=eps())
        @test vals ≈ [1.0, 3.0, -5.0, 6.0] atol = 1e3 * eps()
    end
end

@testset "Plans" begin
    m = E1.model
    p = Plan(m, 1:3)
    @test first(p.range) == ii(0)
    @test last(p.range) == ii(4)
    out = @capture_out print(p)
    @test length(split(out, '\n')) == 2
    @test p[1] == [:y_shk]
    @test p[ii(1)] == [:y_shk]
    endogenize!(p, :y_shk, ii(1))
    @test isempty(p[ii(1)])
    endogenize!(p, :y_shk, ii(1):ii(3))
    exogenize!(p, :y, ii(2))
    exogenize!(p, :y, ii(4))
    # make sure indexing with integers works as well
    @test p[ii(0)] == p[1] == [:y_shk]
    @test p[ii(1)] == p[2] == []
    @test p[ii(2)] == p[3] == [:y]
    @test p[ii(3)] == p[4] == []
    @test p[ii(4)] == p[5] == [:y, :y_shk]
    out = @capture_out print(p)
    @test length(split(out, '\n')) == 6
    out = @capture_out print(IOContext(stdout, :displaysize => (7, 80)), p)
    @test length(split(out, '\n')) == 4
    @test length(split(out, '⋮')) == 2
    let p = Plan(m, 2000Q1:2020Q4)
        endogenize!(p, shocks(m), 2000Q1:2002Q4)
        @test isempty(p[2000Q1])
        out = @capture_out print(p)
        length(split(out, "\n")) == 4
    end
end

include("simdatatests.jl")
include("sstests.jl")

using StateSpaceEcon.StackedTimeSolver: dict2array, dict2data, array2dict, array2data, data2dict, data2array
@testset "dict2array" begin
    m = E3.model
    sim = m.maxlag .+ (1:10)
    p = Plan(m, sim)

    # random data
    d = zerodict(m, p)
    for v in values(d)
        v .= rand(Float64, size(v))
    end

    @test size(dict2array(d, [:pinf, :ygap])) == (length(p.range), 2)
    @test size(dict2array(d, ["pinf", "ygap"])) == (length(p.range), 2)
    @test size(dict2array(d, m.variables)) == (length(p.range), 3)
    a = dict2array(d, m.allvars)
    @test size(a) == (length(p.range), 6)
    @test a == hcat((d[string(v)] for v in m.allvars)...)

    # error variable missing from dictionary
    @test_throws ArgumentError dict2array(d, [:pinf, :ygap, :nosuchvar])
    # error out of range
    @test_throws ArgumentError dict2array(d, [:pinf, :ygap], range=10U:20U)

    # warning variables with different ranges
    d["wrong_var"] = TSeries(3U, rand(10))
    b = @test_logs (:warn, r".*Using\s+intersection\s+range: 3U:12U") dict2array(d, [:pinf, :wrong_var])
    @test size(b) == (10, 2)
    @test b[:,1] == d["pinf"][3U:12U].values
    @test b[:,2] == d["wrong_var"][3U:12U].values

    s = dict2data(d, m.allvars)
    @test s == a

    sa = data2array(s)
    @test sa == s
    s.pinf[3U:5U] = 3:5
    @test sa == s

    sd = data2dict(s)
    @test Set(keys(sd)) == Set(string.(colnames(s)))

    as = array2data(a, m.allvars, first(p.range))
    @test as == a
    a[1,2] = 2.5
    @test as == a
    as = array2data(a, m.allvars, first(p.range), copy=true)
    @test as == a
    a[1,2] = 3.0
    @test as != a
    as[1,2] = 3.0
    @test as == a

    ad = array2dict(a, m.allvars, first(p.range))
    @test length(ad) == size(a, 2)
    @test all(ad[string(v)].values == a[:,i] for (i, v) in enumerate(m.allvars))
end

include("simtests.jl")
include("logsimtests.jl")


text:
AbstractTables.tblrowdim(::Table) = AbstractTables.HasRowDim()
AbstractTables.column_indexable(::Table) = true


text:
"""Returns a sparse boolean adjacency matrix for a graph, indexed by `[u, v]`
vertices. `true` values indicate an edge between `u` and `v`. Users may
specify a direction (`:in`, `:out`, or `:both` are currently supported; `:out`
is default for both directed and undirected graphs) and a data type for the
matrix (defaults to `Int`).

Note: This function is optimized for speed.
"""
function adjacency_matrix(g::SimpleGraph, dir::Symbol=:out, T::DataType=Int)
    n_v = nv(g)
    nz = ne(g) * (is_directed(g)? 1 : 2)
    colpt = ones(Int, n_v + 1)

    if dir == :out
        neighborfn = out_neighbors
    elseif dir == :in
        neighborfn = in_neighbors
    elseif dir == :both
        if is_directed(g)
            neighborfn = all_neighbors
            nz *= 2
        else
            neighborfn = out_neighbors
        end
    else
        error("Not implemented")
    end
    rowval = sizehint!(Vector{Int}(), nz)
    selfloops = Vector{Int}()
    for j in 1:n_v
        if has_edge(g,j,j)
            push!(selfloops, j)
        end
        dsts = neighborfn(g, j)
        colpt[j+1] = colpt[j] + length(dsts)
        append!(rowval, sort!(dsts))
    end
    spmx = SparseMatrixCSC(n_v,n_v,colpt,rowval,ones(T,nz))

    # this is inefficient. There should be a better way of doing this.
    # the issue is that adjacency matrix entries for self-loops are 2,
    # not one(T).
    for i in selfloops
        if !(T <: Bool)
            spmx[i,i] += one(T)
        end
    end
    return spmx
end

"""Returns a sparse [Laplacian matrix](https://en.wikipedia.org/wiki/Laplacian_matrix)
for a graph `g`, indexed by `[u, v]` vertices. For undirected graphs, `dir`
defaults to `:out`; for directed graphs, `dir` defaults to `:both`. `T`
defaults to `Int` for both graph types.
"""
function laplacian_matrix(g::Graph, dir::Symbol=:out, T::DataType=Int)
    A = adjacency_matrix(g, dir, T)
    D = spdiagm(sum(A,2)[:])
    return D - A
end

function laplacian_matrix(g::DiGraph, dir::Symbol=:both, T::DataType=Int)
    A = adjacency_matrix(g, dir, T)
    D = spdiagm(sum(A,2)[:])
    return D - A
end

doc"""Returns the eigenvalues of the Laplacian matrix for a graph `g`, indexed
by vertex. Warning: Converts the matrix to dense with $nv^2$ memory usage. Use
`eigs(laplacian_matrix(g);  kwargs...)` to compute some of the
eigenvalues/eigenvectors. Default values for `dir` and `T` are the same as
`laplacian_matrix`.
"""
laplacian_spectrum(g::Graph, dir::Symbol=:out, T::DataType=Int) = eigvals(full(laplacian_matrix(g, dir, T)))
laplacian_spectrum(g::DiGraph, dir::Symbol=:both, T::DataType=Int) = eigvals(full(laplacian_matrix(g, dir, T)))

doc"""Returns the eigenvalues of the adjacency matrix for a graph `g`, indexed
by vertex. Warning: Converts the matrix to dense with $nv^2$ memory usage. Use
`eigs(adjacency_matrix(g);kwargs...)` to compute some of the
eigenvalues/eigenvectors. Default values for `dir` and `T` are the same as
`adjacency_matrix`.
"""
adjacency_spectrum(g::Graph, dir::Symbol=:out, T::DataType=Int) = eigvals(full(adjacency_matrix(g, dir, T)))
adjacency_spectrum(g::DiGraph, dir::Symbol=:both, T::DataType=Int) = eigvals(full(adjacency_matrix(g, dir, T)))

"""Returns a sparse node-arc incidence matrix for a graph, indexed by
`[v, i]`, where `i` is in `1:ne(g)`, indexing an edge `e`. For
directed graphs, a value of `-1` indicates that `src(e) == v`, while a
value of `1` indicates that `dst(e) == v`. Otherwise, the value is
`0`. For undirected graphs, both entries are `1`.
"""
function incidence_matrix(g::SimpleGraph, T::DataType=Int)
    isdir = is_directed(g)
    n_v = nv(g)
    n_e = ne(g)
    nz = 2 * n_e

    # every col has the same 2 entries
    colpt = collect(1:2:(nz + 1))
    nzval = repmat([isdir ? -one(T) : one(T), one(T)], n_e)

    # iterate over edges for row indices
    rowval = zeros(Int, nz)
    i = 1
    for u in vertices(g)
        for v in out_neighbors(g, u)
            if isdir || u < v # add every edge only once
                rowval[2*i - 1] = u
                rowval[2*i] = v
                i += 1
            end
        end
    end

    spmx = SparseMatrixCSC(n_v,n_e,colpt,rowval,nzval)
    return spmx
end

"""
Given two oriented edges i->j and k->l in g, the
non-backtraking matrix B is defined as

B[i->j, k->l] = δ(j,k)* (1 - δ(i,l))

returns a matrix B, and an edgemap storing the oriented edges' positions in B
"""
function non_backtracking_matrix(g::SimpleGraph)
    # idedgemap = Dict{Int, Edge}()
    edgeidmap = Dict{Edge, Int}()
    m = 0
    for e in edges(g)
        m += 1
        edgeidmap[e] = m
    end

    if !is_directed(g)
        for e in edges(g)
            m += 1
            edgeidmap[reverse(e)] = m
        end
    end

    B = zeros(Float64, m, m)

    for (e,u) in edgeidmap
        i, j = src(e), dst(e)
        for k in in_neighbors(g,i)
            k == j && continue
            v = edgeidmap[Edge(k, i)]
            B[v, u] = 1
        end
    end

    return B, edgeidmap
end

"""Nonbacktracking: a compact representation of the nonbacktracking operator

    g: the underlying graph
    edgeidmap: the association between oriented edges and index into the NBT matrix

The Nonbacktracking operator can be used for community detection.
This representation is compact in that it uses only ne(g) additional storage
and provides an implicit representation of the matrix B_g defined below.

Given two oriented edges i->j and k->l in g, the
non-backtraking matrix B is defined as

B[i->j, k->l] = δ(j,k)* (1 - δ(i,l))

This type is in the style of GraphMatrices.jl and supports the necessary operations
for computed eigenvectors and conducting linear solves.

Additionally the contract!(vertexspace, nbt, edgespace) method takes vectors represented in
the domain of B and represents them in the domain of the adjacency matrix of g.
"""
type Nonbacktracking{G}
    g::G
    edgeidmap::Dict{Edge,Int}
    m::Int
end

function Nonbacktracking(g::SimpleGraph)
    edgeidmap = Dict{Edge, Int}()
    m = 0
    for e in edges(g)
        m += 1
        edgeidmap[e] = m
    end
    if !is_directed(g)
        for e in edges(g)
            m += 1
            edgeidmap[reverse(e)] = m
        end
    end
    return Nonbacktracking(g, edgeidmap, m)
end

size(nbt::Nonbacktracking) = (nbt.m,nbt.m)
eltype(nbt::Nonbacktracking) = Float64
issymmetric(nbt::Nonbacktracking) = false

function *{G, T<:Number}(nbt::Nonbacktracking{G}, x::Vector{T})
    length(x) == nbt.m || error("dimension mismatch")
    y = zeros(T, length(x))
    for (e,u) in nbt.edgeidmap
        i, j = src(e), dst(e)
        for k in in_neighbors(nbt.g,i)
            k == j && continue
            v = nbt.edgeidmap[Edge(k, i)]
            y[v] += x[u]
        end
    end
    return y
end

function coo_sparse(nbt::Nonbacktracking)
    m = nbt.m
    #= I,J = zeros(Int, m), zeros(Int, m) =#
    I,J = zeros(Int, 0), zeros(Int, 0)
    for (e,u) in nbt.edgeidmap
        i, j = src(e), dst(e)
        for k in in_neighbors(nbt.g,i)
            k == j && continue
            v = nbt.edgeidmap[Edge(k, i)]
            #= J[u] = v =#
            #= I[u] = u =#
            push!(I, v)
            push!(J, u)
        end
    end
    return I,J,1.0
end

sparse(nbt::Nonbacktracking) = sparse(coo_sparse(nbt)..., nbt.m,nbt.m)

function *{G, T<:Number}(nbt::Nonbacktracking{G}, x::AbstractMatrix{T})
    y = zeros(x)
    for i in 1:nbt.m
        y[:,i] = nbt * x[:,i]
    end
    return y
end

function A_mul_B!(C, nbt::Nonbacktracking, B)
    # computs C = A*B
    for i in 1:size(B,2)
        C[:,i] = nbt*B[:,i]
    end
    return C
end

"""contract!(vertexspace, nbt, edgespace) in place version of
contract(nbt, edgespace). modifies first argument
"""
function contract!(vertexspace::Vector, nbt::Nonbacktracking, edgespace::Vector)
    for i=1:nv(nbt.g)
        for j in neighbors(nbt.g, i)
            u = nbt.edgeidmap[i > j ? Edge(j,i) : Edge(i,j)]
            vertexspace[i] += edgespace[u]
        end
    end
end

"""contract(nbt, edgespace)
Integrates out the edges by summing over the edges incident to each vertex.
"""
function contract(nbt::Nonbacktracking, edgespace::Vector)
    y = zeros(eltype(edgespace), nv(nbt.g))
    contract!(y,nbt,edgespace)
    return y
end

"""spectral_distance(G₁, G₂ [, k])
Compute the spectral distance between undirected n-vertex
graphs G₁ and G₂ using the top k ≤ n greatest eigenvalues.
If k is ommitted, uses full spectrum.

For further details, please refer to:

JOVANOVIC, I.; STANIC, Z., 2014. Spectral Distances of
Graphs Based on their Different Matrix Representations
"""
function spectral_distance(G₁::Graph, G₂::Graph, k::Integer)
  A₁ = adjacency_matrix(G₁)
  A₂ = adjacency_matrix(G₂)

  λ₁ = k < nv(G₁)-1 ? eigs(A₁, nev=k, which=:LR)[1] : eigvals(full(A₁))[end:-1:end-(k-1)]
  λ₂ = k < nv(G₂)-1 ? eigs(A₂, nev=k, which=:LR)[1] : eigvals(full(A₂))[end:-1:end-(k-1)]

  sumabs(λ₁ - λ₂)
end

function spectral_distance(G₁::Graph, G₂::Graph)
  @assert nv(G₁) == nv(G₂) "spectral distance not defined for |G₁| != |G₂|"
  spectral_distance(G₁, G₂, nv(G₁))
end


text:
@testset "fill_dem" begin
    # TODO: need more fill tests
    dx = 0.1
    xs = -1.5:dx:1
    ys = -0.5:dx:3.0
    dem = dem1.(xs, ys', withpit=true)
    area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem, bnd_as_pits=false)
    demf = WWF.fill_dem(dem, pits, dir)
    @test sum(demf.-dem) ≈ 2.1499674517313414
    @test sum(demf.-dem .> 0) == 5
    @test all([c[pits[cc]]==cc  for cc=1:length(pits)]) # pit in catchment of same color
end

@testset "catchment" begin
    for demfn in [peaks, peaks2, peaks2_nan, peaks2_nan_edge]
        xs, dem = demfn()
        ys = xs

        @test size(dem)==(length(xs), length(ys))
        area, slen, dir, nout, nin, pits, c, bnds = WWF.waterflows(dem);

        for cc=1:length(pits)
            ij = pits[cc]
            @test catchment(dir, ij) == (c.==cc)
        end

        for (cc,dd) in zip(1:length(pits), length(pits):-1:1)
            ii, jj = pits[cc], pits[dd]
            @test catchment(dir, [ii,jj]) == ((c.==cc) .| (c.==dd))
        end

        ci = CartesianIndices((2:4,7:9))
        @test catchment(dir, ci) == catchment(dir, collect(ci)[:])

        @test catchment_flux(ones(size(dir)), catchment(dir, ci)) > 0
        @test catchment_flux(ones(size(dir)), c, 5) > 0
        @test catchment_flux(zeros(size(dir)), catchment(dir, ci)) == 0

        sinks = [CartesianIndices((2:4,7:9)), CartesianIndices((5:8,1:3))]
        ss = [collect(sinks[1])[:]; collect(sinks[2])[:]]
        @test (catchments(dir, sinks) .> 0) == catchment(dir, ss)

        @test length(unique(prune_catchments(c, 10; val=0))) < length(unique(c))
    end
end


text:
using Distributed
using SharedArrays
using Formatting
using Printf
using ArgParse
using Dates

include("ProgramTunnel/src/julia/ProgramTunnel_fs_new.jl")
include("julia_lib/NetCDFIO.jl")
include("julia_lib/parseMsg.jl")
include("julia_lib/BinaryIO.jl")
include("../share/ModelTime.jl")

using JSON
using .ProgramTunnel_fs
using .NetCDFIO
using .BinaryIO
using .ModelTime

include("config.jl")
include("driver.jl")

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
module CellArraysTests

using Test
using Gridap

using ..CellValuesMocks

l = 10
a = rand(3,2,4,2)

cn = TestIterCellValue(a,l)
@test isa(cn,CellArray)
@test isa(cn,IterCellArray)
r = [a for i in 1:l]
test_iter_cell_array(cn,r)
@test !isa(cn,CellNumber)
@test !isa(cn,IterCellNumber)
@test !isa(cn,CellMap)
@test !isa(cn,IterCellMap)

cn = TestIndexCellValue(a,l)
@test isa(cn,CellArray)
@test isa(cn,IndexCellArray)
test_index_cell_array(cn,r)
@test !isa(cn,CellNumber)
@test !isa(cn,IndexCellNumber)
@test !isa(cn,CellMap)
@test !isa(cn,IndexCellMap)

a = rand(3)

cn = TestIterCellValue(a,l)
@test isa(cn,CellVector)
@test isa(cn,IterCellVector)
@test !isa(cn,CellNumber)
@test !isa(cn,IterCellNumber)
@test !isa(cn,CellMap)
@test !isa(cn,IterCellMap)

cn = TestIndexCellValue(a,l)
@test isa(cn,CellVector)
@test isa(cn,IndexCellVector)
@test !isa(cn,CellNumber)
@test !isa(cn,IndexCellNumber)
@test !isa(cn,CellMap)
@test !isa(cn,IndexCellMap)

a = rand(3,5)

cn = TestIterCellValue(a,l)
@test isa(cn,CellMatrix)
@test isa(cn,IterCellMatrix)
@test !isa(cn,CellNumber)
@test !isa(cn,IterCellNumber)
@test !isa(cn,CellMap)
@test !isa(cn,IterCellMap)

cn = TestIndexCellValue(a,l)
@test isa(cn,CellMatrix)
@test isa(cn,IndexCellMatrix)
@test !isa(cn,CellNumber)
@test !isa(cn,IndexCellNumber)
@test !isa(cn,CellMap)
@test !isa(cn,IndexCellMap)

end # module


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: mq
using AWS.Compat
using AWS.UUIDs

"""
    create_broker(auto_minor_version_upgrade, broker_name, deployment_mode, engine_type, engine_version, host_instance_type, publicly_accessible, users)
    create_broker(auto_minor_version_upgrade, broker_name, deployment_mode, engine_type, engine_version, host_instance_type, publicly_accessible, users, params::Dict{String,<:Any})

Creates a broker. Note: This API is asynchronous. To create a broker, you must either use
the AmazonMQFullAccess IAM policy or include the following EC2 permissions in your IAM
policy. ec2:CreateNetworkInterface This permission is required to allow Amazon MQ to create
an elastic network interface (ENI) on behalf of your account.
ec2:CreateNetworkInterfacePermission This permission is required to attach the ENI to the
broker instance. ec2:DeleteNetworkInterface ec2:DeleteNetworkInterfacePermission
ec2:DetachNetworkInterface ec2:DescribeInternetGateways ec2:DescribeNetworkInterfaces
ec2:DescribeNetworkInterfacePermissions ec2:DescribeRouteTables ec2:DescribeSecurityGroups
ec2:DescribeSubnets ec2:DescribeVpcs For more information, see Create an IAM User and Get
Your AWS Credentials and Never Modify or Delete the Amazon MQ Elastic Network Interface in
the Amazon MQ Developer Guide.

# Arguments
- `auto_minor_version_upgrade`: Enables automatic upgrades to new minor versions for
  brokers, as new versions are released and supported by Amazon MQ. Automatic upgrades occur
  during the scheduled maintenance window of the broker or after a manual broker reboot. Set
  to true by default, if no value is specified.
- `broker_name`: Required. The broker's name. This value must be unique in your AWS
  account, 1-50 characters long, must contain only letters, numbers, dashes, and underscores,
  and must not contain white spaces, brackets, wildcard characters, or special characters.
- `deployment_mode`: Required. The broker's deployment mode.
- `engine_type`: Required. The type of broker engine. Currently, Amazon MQ supports
  ACTIVEMQ and RABBITMQ.
- `engine_version`: Required. The broker engine's version. For a list of supported engine
  versions, see Supported engines.
- `host_instance_type`: Required. The broker's instance type.
- `publicly_accessible`: Enables connections from applications outside of the VPC that
  hosts the broker's subnets. Set to false by default, if no value is provided.
- `users`: Required. The list of broker users (persons or applications) who can access
  queues and topics. This value can contain only alphanumeric characters, dashes, periods,
  underscores, and tildes (- . _ ~). This value must be 2-100 characters long. Amazon MQ for
  RabbitMQ When you create an Amazon MQ for RabbitMQ broker, one and only one administrative
  user is accepted and created when a broker is first provisioned. All subsequent broker
  users are created by making RabbitMQ API calls directly to brokers or via the RabbitMQ web
  console.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"authenticationStrategy"`: Optional. The authentication strategy used to secure the
  broker. The default is SIMPLE.
- `"configuration"`: A list of information about the configuration.
- `"creatorRequestId"`: The unique ID that the requester receives for the created broker.
  Amazon MQ passes your ID with the API action. Note: We recommend using a Universally Unique
  Identifier (UUID) for the creatorRequestId. You may omit the creatorRequestId if your
  application doesn't require idempotency.
- `"encryptionOptions"`: Encryption options for the broker. Does not apply to RabbitMQ
  brokers.
- `"ldapServerMetadata"`: Optional. The metadata of the LDAP server used to authenticate
  and authorize connections to the broker. Does not apply to RabbitMQ brokers.
- `"logs"`: Enables Amazon CloudWatch logging for brokers.
- `"maintenanceWindowStartTime"`: The parameters that determine the WeeklyStartTime.
- `"securityGroups"`: The list of rules (1 minimum, 125 maximum) that authorize connections
  to brokers.
- `"storageType"`: The broker's storage type.
- `"subnetIds"`: The list of groups that define which subnets and IP ranges the broker can
  use from different Availability Zones. If you specify more than one subnet, the subnets
  must be in different Availability Zones. Amazon MQ will not be able to create VPC endpoints
  for your broker with multiple subnets in the same Availability Zone. A SINGLE_INSTANCE
  deployment requires one subnet (for example, the default subnet). An
  ACTIVE_STANDBY_MULTI_AZ Amazon MQ for ActiveMQ deployment requires two subnets. A
  CLUSTER_MULTI_AZ Amazon MQ for RabbitMQ deployment has no subnet requirements when deployed
  with public accessibility. Deployment without public accessibility requires at least one
  subnet. If you specify subnets in a shared VPC for a RabbitMQ broker, the associated VPC to
  which the specified subnets belong must be owned by your AWS account. Amazon MQ will not be
  able to create VPC endpoints in VPCs that are not owned by your AWS account.
- `"tags"`: Create tags when creating the broker.
"""
function create_broker(
    autoMinorVersionUpgrade,
    brokerName,
    deploymentMode,
    engineType,
    engineVersion,
    hostInstanceType,
    publiclyAccessible,
    users;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mq(
        "POST",
        "/v1/brokers",
        Dict{String,Any}(
            "autoMinorVersionUpgrade" => autoMinorVersionUpgrade,
            "brokerName" => brokerName,
            "deploymentMode" => deploymentMode,
            "engineType" => engineType,
            "engineVersion" => engineVersion,
            "hostInstanceType" => hostInstanceType,
            "publiclyAccessible" => publiclyAccessible,
            "users" => users,
            "creatorRequestId" => string(uuid4()),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_broker(
    autoMinorVersionUpgrade,
    brokerName,
    deploymentMode,
    engineType,
    engineVersion,
    hostInstanceType,
    publiclyAccessible,
    users,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mq(
        "POST",
        "/v1/brokers",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "autoMinorVersionUpgrade" => autoMinorVersionUpgrade,
                    "brokerName" => brokerName,
                    "deploymentMode" => deploymentMode,
                    "engineType" => engineType,
                    "engineVersion" => engineVersion,
                    "hostInstanceType" => hostInstanceType,
                    "publiclyAccessible" => publiclyAccessible,
                    "users" => users,
                    "creatorRequestId" => string(uuid4()),
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_configuration(engine_type, engine_version, name)
    create_configuration(engine_type, engine_version, name, params::Dict{String,<:Any})

Creates a new configuration for the specified configuration name. Amazon MQ uses the
default configuration (the engine type and version).

# Arguments
- `engine_type`: Required. The type of broker engine. Currently, Amazon MQ supports
  ACTIVEMQ and RABBITMQ.
- `engine_version`: Required. The broker engine's version. For a list of supported engine
  versions, see Supported engines.
- `name`: Required. The name of the configuration. This value can contain only alphanumeric
  characters, dashes, periods, underscores, and tildes (- . _ ~). This value must be 1-150
  characters long.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"authenticationStrategy"`: Optional. The authentication strategy associated with the
  configuration. The default is SIMPLE.
- `"tags"`: Create tags when creating the configuration.
"""
function create_configuration(
    engineType, engineVersion, name; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mq(
        "POST",
        "/v1/configurations",
        Dict{String,Any}(
            "engineType" => engineType, "engineVersion" => engineVersion, "name" => name
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_configuration(
    engineType,
    engineVersion,
    name,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mq(
        "POST",
        "/v1/configurations",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "engineType" => engineType,
                    "engineVersion" => engineVersion,
                    "name" => name,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_tags(resource-arn)
    create_tags(resource-arn, params::Dict{String,<:Any})

Add a tag to a resource.

# Arguments
- `resource-arn`: The Amazon Resource Name (ARN) of the resource tag.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"tags"`: The key-value pair for the resource tag.
"""
function create_tags(resource_arn; aws_config::AbstractAWSConfig=global_aws_config())
    return mq(
        "POST",
        "/v1/tags/$(resource-arn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_tags(
    resource_arn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mq(
        "POST",
        "/v1/tags/$(resource-arn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_user(broker-id, password, username)
    create_user(broker-id, password, username, params::Dict{String,<:Any})

Creates an ActiveMQ user.

# Arguments
- `broker-id`: The unique ID that Amazon MQ generates for the broker.
- `password`: Required. The password of the user. This value must be at least 12 characters
  long, must contain at least 4 unique characters, and must not contain commas, colons, or
  equal signs (,:=).
- `username`: The username of the ActiveMQ user. This value can contain only alphanumeric
  characters, dashes, periods, underscores, and tildes (- . _ ~). This value must be 2-100
  characters long.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"consoleAccess"`: Enables access to the ActiveMQ Web Console for the ActiveMQ user.
- `"groups"`: The list of groups (20 maximum) to which the ActiveMQ user belongs. This
  value can contain only alphanumeric characters, dashes, periods, underscores, and tildes (-
  . _ ~). This value must be 2-100 characters long.
"""
function create_user(
    broker_id, password, username; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mq(
        "POST",
        "/v1/brokers/$(broker-id)/users/$(username)",
        Dict{String,Any}("password" => password);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_user(
    broker_id,
    password,
    username,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mq(
        "POST",
        "/v1/brokers/$(broker-id)/users/$(username)",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("password" => password), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_broker(broker-id)
    delete_broker(broker-id, params::Dict{String,<:Any})

Deletes a broker. Note: This API is asynchronous.

# Arguments
- `broker-id`: The unique ID that Amazon MQ generates for the broker.

"""
function delete_broker(broker_id; aws_config::AbstractAWSConfig=global_aws_config())
    return mq(
        "DELETE",
        "/v1/brokers/$(broker-id)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_broker(
    broker_id,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mq(
        "DELETE",
        "/v1/brokers/$(broker-id)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_tags(resource-arn, tag_keys)
    delete_tags(resource-arn, tag_keys, params::Dict{String,<:Any})

Removes a tag from a resource.

# Arguments
- `resource-arn`: The Amazon Resource Name (ARN) of the resource tag.
- `tag_keys`: An array of tag keys to delete

"""
function delete_tags(
    resource_arn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mq(
        "DELETE",
        "/v1/tags/$(resource-arn)",
        Dict{String,Any}("tagKeys" => tagKeys);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_tags(
    resource_arn,
    tagKeys,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mq(
        "DELETE",
        "/v1/tags/$(resource-arn)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("tagKeys" => tagKeys), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_user(broker-id, username)
    delete_user(broker-id, username, params::Dict{String,<:Any})

Deletes an ActiveMQ user.

# Arguments
- `broker-id`: The unique ID that Amazon MQ generates for the broker.
- `username`: The username of the ActiveMQ user. This value can contain only alphanumeric
  characters, dashes, periods, underscores, and tildes (- . _ ~). This value must be 2-100
  characters long.

"""
function delete_user(broker_id, username; aws_config::AbstractAWSConfig=global_aws_config())
    return mq(
        "DELETE",
        "/v1/brokers/$(broker-id)/users/$(username)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_user(
    broker_id,
    username,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mq(
        "DELETE",
        "/v1/brokers/$(broker-id)/users/$(username)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_broker(broker-id)
    describe_broker(broker-id, params::Dict{String,<:Any})

Returns information about the specified broker.

# Arguments
- `broker-id`: The unique ID that Amazon MQ generates for the broker.

"""
function describe_broker(broker_id; aws_config::AbstractAWSConfig=global_aws_config())
    return mq(
        "GET",
        "/v1/brokers/$(broker-id)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_broker(
    broker_id,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mq(
        "GET",
        "/v1/brokers/$(broker-id)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_broker_engine_types()
    describe_broker_engine_types(params::Dict{String,<:Any})

Describe available engine types and versions.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"engineType"`: Filter response by engine type.
- `"maxResults"`: The maximum number of brokers that Amazon MQ can return per page (20 by
  default). This value must be an integer from 5 to 100.
- `"nextToken"`: The token that specifies the next page of results Amazon MQ should return.
  To request the first page, leave nextToken empty.
"""
function describe_broker_engine_types(; aws_config::AbstractAWSConfig=global_aws_config())
    return mq(
        "GET",
        "/v1/broker-engine-types";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_broker_engine_types(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mq(
        "GET",
        "/v1/broker-engine-types",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_broker_instance_options()
    describe_broker_instance_options(params::Dict{String,<:Any})

Describe available broker instance options.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"engineType"`: Filter response by engine type.
- `"hostInstanceType"`: Filter response by host instance type.
- `"maxResults"`: The maximum number of brokers that Amazon MQ can return per page (20 by
  default). This value must be an integer from 5 to 100.
- `"nextToken"`: The token that specifies the next page of results Amazon MQ should return.
  To request the first page, leave nextToken empty.
- `"storageType"`: Filter response by storage type.
"""
function describe_broker_instance_options(;
    aws_config::AbstractAWSConfig=global_aws_config()
)
    return mq(
        "GET",
        "/v1/broker-instance-options";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_broker_instance_options(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mq(
        "GET",
        "/v1/broker-instance-options",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_configuration(configuration-id)
    describe_configuration(configuration-id, params::Dict{String,<:Any})

Returns information about the specified configuration.

# Arguments
- `configuration-id`: The unique ID that Amazon MQ generates for the configuration.

"""
function describe_configuration(
    configuration_id; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mq(
        "GET",
        "/v1/configurations/$(configuration-id)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_configuration(
    configuration_id,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mq(
        "GET",
        "/v1/configurations/$(configuration-id)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_configuration_revision(configuration-id, configuration-revision)
    describe_configuration_revision(configuration-id, configuration-revision, params::Dict{String,<:Any})

Returns the specified configuration revision for the specified configuration.

# Arguments
- `configuration-id`: The unique ID that Amazon MQ generates for the configuration.
- `configuration-revision`: The revision of the configuration.

"""
function describe_configuration_revision(
    configuration_id,
    configuration_revision;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mq(
        "GET",
        "/v1/configurations/$(configuration-id)/revisions/$(configuration-revision)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_configuration_revision(
    configuration_id,
    configuration_revision,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mq(
        "GET",
        "/v1/configurations/$(configuration-id)/revisions/$(configuration-revision)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_user(broker-id, username)
    describe_user(broker-id, username, params::Dict{String,<:Any})

Returns information about an ActiveMQ user.

# Arguments
- `broker-id`: The unique ID that Amazon MQ generates for the broker.
- `username`: The username of the ActiveMQ user. This value can contain only alphanumeric
  characters, dashes, periods, underscores, and tildes (- . _ ~). This value must be 2-100
  characters long.

"""
function describe_user(
    broker_id, username; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mq(
        "GET",
        "/v1/brokers/$(broker-id)/users/$(username)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_user(
    broker_id,
    username,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mq(
        "GET",
        "/v1/brokers/$(broker-id)/users/$(username)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_brokers()
    list_brokers(params::Dict{String,<:Any})

Returns a list of all brokers.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of brokers that Amazon MQ can return per page (20 by
  default). This value must be an integer from 5 to 100.
- `"nextToken"`: The token that specifies the next page of results Amazon MQ should return.
  To request the first page, leave nextToken empty.
"""
function list_brokers(; aws_config::AbstractAWSConfig=global_aws_config())
    return mq("GET", "/v1/brokers"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET)
end
function list_brokers(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mq(
        "GET", "/v1/brokers", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end

"""
    list_configuration_revisions(configuration-id)
    list_configuration_revisions(configuration-id, params::Dict{String,<:Any})

Returns a list of all revisions for the specified configuration.

# Arguments
- `configuration-id`: The unique ID that Amazon MQ generates for the configuration.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of brokers that Amazon MQ can return per page (20 by
  default). This value must be an integer from 5 to 100.
- `"nextToken"`: The token that specifies the next page of results Amazon MQ should return.
  To request the first page, leave nextToken empty.
"""
function list_configuration_revisions(
    configuration_id; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mq(
        "GET",
        "/v1/configurations/$(configuration-id)/revisions";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_configuration_revisions(
    configuration_id,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mq(
        "GET",
        "/v1/configurations/$(configuration-id)/revisions",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_configurations()
    list_configurations(params::Dict{String,<:Any})

Returns a list of all configurations.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of brokers that Amazon MQ can return per page (20 by
  default). This value must be an integer from 5 to 100.
- `"nextToken"`: The token that specifies the next page of results Amazon MQ should return.
  To request the first page, leave nextToken empty.
"""
function list_configurations(; aws_config::AbstractAWSConfig=global_aws_config())
    return mq(
        "GET", "/v1/configurations"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_configurations(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mq(
        "GET",
        "/v1/configurations",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_tags(resource-arn)
    list_tags(resource-arn, params::Dict{String,<:Any})

Lists tags for a resource.

# Arguments
- `resource-arn`: The Amazon Resource Name (ARN) of the resource tag.

"""
function list_tags(resource_arn; aws_config::AbstractAWSConfig=global_aws_config())
    return mq(
        "GET",
        "/v1/tags/$(resource-arn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_tags(
    resource_arn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mq(
        "GET",
        "/v1/tags/$(resource-arn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_users(broker-id)
    list_users(broker-id, params::Dict{String,<:Any})

Returns a list of all ActiveMQ users.

# Arguments
- `broker-id`: The unique ID that Amazon MQ generates for the broker.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of brokers that Amazon MQ can return per page (20 by
  default). This value must be an integer from 5 to 100.
- `"nextToken"`: The token that specifies the next page of results Amazon MQ should return.
  To request the first page, leave nextToken empty.
"""
function list_users(broker_id; aws_config::AbstractAWSConfig=global_aws_config())
    return mq(
        "GET",
        "/v1/brokers/$(broker-id)/users";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_users(
    broker_id,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mq(
        "GET",
        "/v1/brokers/$(broker-id)/users",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    reboot_broker(broker-id)
    reboot_broker(broker-id, params::Dict{String,<:Any})

Reboots a broker. Note: This API is asynchronous.

# Arguments
- `broker-id`: The unique ID that Amazon MQ generates for the broker.

"""
function reboot_broker(broker_id; aws_config::AbstractAWSConfig=global_aws_config())
    return mq(
        "POST",
        "/v1/brokers/$(broker-id)/reboot";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function reboot_broker(
    broker_id,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mq(
        "POST",
        "/v1/brokers/$(broker-id)/reboot",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    update_broker(broker-id)
    update_broker(broker-id, params::Dict{String,<:Any})

Adds a pending configuration change to a broker.

# Arguments
- `broker-id`: The unique ID that Amazon MQ generates for the broker.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"authenticationStrategy"`: Optional. The authentication strategy used to secure the
  broker. The default is SIMPLE.
- `"autoMinorVersionUpgrade"`: Enables automatic upgrades to new minor versions for
  brokers, as new versions are released and supported by Amazon MQ. Automatic upgrades occur
  during the scheduled maintenance window of the broker or after a manual broker reboot.
- `"configuration"`: A list of information about the configuration.
- `"engineVersion"`: The broker engine version. For a list of supported engine versions,
  see Supported engines.
- `"hostInstanceType"`: The broker's host instance type to upgrade to. For a list of
  supported instance types, see Broker instance types.
- `"ldapServerMetadata"`: Optional. The metadata of the LDAP server used to authenticate
  and authorize connections to the broker. Does not apply to RabbitMQ brokers.
- `"logs"`: Enables Amazon CloudWatch logging for brokers.
- `"maintenanceWindowStartTime"`: The parameters that determine the WeeklyStartTime.
- `"securityGroups"`: The list of security groups (1 minimum, 5 maximum) that authorizes
  connections to brokers.
"""
function update_broker(broker_id; aws_config::AbstractAWSConfig=global_aws_config())
    return mq(
        "PUT",
        "/v1/brokers/$(broker-id)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_broker(
    broker_id,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mq(
        "PUT",
        "/v1/brokers/$(broker-id)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    update_configuration(configuration-id, data)
    update_configuration(configuration-id, data, params::Dict{String,<:Any})

Updates the specified configuration.

# Arguments
- `configuration-id`: The unique ID that Amazon MQ generates for the configuration.
- `data`: Required. The base64-encoded XML configuration.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"description"`: The description of the configuration.
"""
function update_configuration(
    configuration_id, data; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mq(
        "PUT",
        "/v1/configurations/$(configuration-id)",
        Dict{String,Any}("data" => data);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_configuration(
    configuration_id,
    data,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mq(
        "PUT",
        "/v1/configurations/$(configuration-id)",
        Dict{String,Any}(mergewith(_merge, Dict{String,Any}("data" => data), params));
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    update_user(broker-id, username)
    update_user(broker-id, username, params::Dict{String,<:Any})

Updates the information for an ActiveMQ user.

# Arguments
- `broker-id`: The unique ID that Amazon MQ generates for the broker.
- `username`: The username of the ActiveMQ user. This value can contain only alphanumeric
  characters, dashes, periods, underscores, and tildes (- . _ ~). This value must be 2-100
  characters long.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"consoleAccess"`: Enables access to the the ActiveMQ Web Console for the ActiveMQ user.
- `"groups"`: The list of groups (20 maximum) to which the ActiveMQ user belongs. This
  value can contain only alphanumeric characters, dashes, periods, underscores, and tildes (-
  . _ ~). This value must be 2-100 characters long.
- `"password"`: The password of the user. This value must be at least 12 characters long,
  must contain at least 4 unique characters, and must not contain commas, colons, or equal
  signs (,:=).
"""
function update_user(broker_id, username; aws_config::AbstractAWSConfig=global_aws_config())
    return mq(
        "PUT",
        "/v1/brokers/$(broker-id)/users/$(username)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_user(
    broker_id,
    username,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mq(
        "PUT",
        "/v1/brokers/$(broker-id)/users/$(username)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
import Pkg
import GitCommand
import RegistryTools
import Test

function gather_stdlib_uuids()
    return Set{Base.UUID}(x for x in keys(RegistryTools.stdlibs()))
end

# For when you have a registry that has packages with dependencies obtained from
# another dependency registry. For example, packages registered at the BioJuliaRegistry
# that have General dependencies. BJW.
function load_registry_dep_uuids(registry_deps_urls::Vector{<:AbstractString} = String[])
    return with_temp_depot() do
        # Get the registries!
        for url in registry_deps_urls
            Pkg.Registry.add(Pkg.RegistrySpec(url = url))
        end
        # Now use the RegistrySpec's to find the Project.toml's. I know
        # .julia/registires/XYZ/ABC is the most likely place, but this way the
        # function never has to assume. BJW.
        extrauuids = Set{Base.UUID}()
        for spec in Pkg.Types.collect_registries()
            if spec.url ∈ registry_deps_urls
                reg = Pkg.TOML.parsefile(joinpath(spec.path, "Registry.toml"))
                for x in keys(reg["packages"])
                    push!(extrauuids, Base.UUID(x))
                end
            end
        end
        return extrauuids
    end
end

#########################
# Testing of registries #
#########################

"""
    test(path)

Run various checks on the registry located at `path`.
Checks for example that all files are parsable and
understandable by Pkg and consistency between Registry.toml
and each Package.toml.

If your registry has packages that have dependencies that are registered in other
registries elsewhere, then you may provide the github urls for those registries
using the `registry_deps` parameter.
"""
function test(path = pwd();
              registry_deps::Vector{<:AbstractString} = String[])
    Test.@testset "(Registry|Package|Versions|Deps|Compat).toml" begin; cd(path) do
        reg = Pkg.TOML.parsefile("Registry.toml")
        reguuids = Set{Base.UUID}(Base.UUID(x) for x in keys(reg["packages"]))
        stdlibuuids = gather_stdlib_uuids()
        registry_dep_uuids = load_registry_dep_uuids(registry_deps)
        alluuids = reguuids ∪ stdlibuuids ∪ registry_dep_uuids

        # Test that each entry in Registry.toml has a corresponding Package.toml
        # at the expected path with the correct uuid and name
        for (uuid, data) in reg["packages"]
            # Package.toml testing
            pkg = Pkg.TOML.parsefile(abspath(data["path"], "Package.toml"))
            Test.@test Base.UUID(uuid) == Base.UUID(pkg["uuid"])
            Test.@test data["name"] == pkg["name"]
            Test.@test Base.isidentifier(data["name"])
            Test.@test haskey(pkg, "repo")

            # Versions.toml testing
            vers = Pkg.TOML.parsefile(abspath(data["path"], "Versions.toml"))
            for (v, data) in vers
                Test.@test VersionNumber(v) isa VersionNumber
                Test.@test haskey(data, "git-tree-sha1")
            end

            # Deps.toml testing
            depsfile = abspath(data["path"], "Deps.toml")
            if isfile(depsfile)
                deps = Pkg.TOML.parsefile(depsfile)
                # Require all deps to exist in the General registry or be a stdlib
                depuuids = Set{Base.UUID}(Base.UUID(x) for (_, d) in deps for (_, x) in d)
                Test.@test depuuids ⊆ alluuids
                # Test that the way Pkg loads this data works
                Test.@test Pkg.Operations.load_package_data_raw(Base.UUID, depsfile) isa
                    Dict{Pkg.Types.VersionRange,Dict{String,Base.UUID}}
            end

            # Compat.toml testing
            compatfile = abspath(data["path"], "Compat.toml")
            if isfile(compatfile)
                compat = Pkg.TOML.parsefile(compatfile)
                # Test that all names with compat is a dependency
                compatnames = Set{String}(x for (_, d) in compat for (x, _) in d)
                if !(isempty(compatnames) || (length(compatnames)== 1 && "julia" in compatnames))
                    depnames = Set{String}(x for (_, d) in Pkg.TOML.parsefile(depsfile) for (x, _) in d)
                    push!(depnames, "julia") # All packages has an implicit dependency on julia
                    @assert compatnames ⊆ depnames
                end
                # Test that the way Pkg loads this data works
                Test.@test Pkg.Operations.load_package_data_raw(Pkg.Types.VersionSpec, compatfile) isa
                    Dict{Pkg.Types.VersionRange,Dict{String,Pkg.Types.VersionSpec}}
            end
        end
    end end
    return
end