text:
function solve_lin_system!(A_mat, b_vec)
    # passing Val(false) to the second argument turns off pivoting. it seems to
    # improve speed for the small matrices that we typically consider. we can
    # revisit this (or make it a parameter) if needed.
    A_fact = lu!(A_mat, Val(false))
    ldiv!(A_fact, b_vec)        # b_vec is overwritten to store the result
    nothing
end

function solve_lin_system_pivot!(A_mat, b_vec)
    # for the larger matrices, in the coupled equations, it's better to pivot
    A_fact = lu!(A_mat, Val(true))
    ldiv!(A_fact, b_vec)        # b_vec is overwritten to store the result
    nothing
end

struct Aux{T<:Real}
    A_mat   :: Matrix{T}
    b_vec   :: Vector{T}
    ABCS    :: Vector{T}
    A_mat2  :: Matrix{T}
    b_vec2  :: Vector{T}
    AA      :: Matrix{T}
    BB      :: Matrix{T}
    CC      :: Matrix{T}
    SS      :: Vector{T}
    function Aux{T}(N::Int) where {T<:Real}
        A_mat   = zeros(T, N, N)
        b_vec   = zeros(T, N)
        ABCS    = zeros(T, 4)
        A_mat2  = zeros(T, 2*N, 2*N)
        b_vec2  = zeros(T, 2*N)
        AA      = zeros(T, 2,2)
        BB      = zeros(T, 2,2)
        CC      = zeros(T, 2,2)
        SS      = zeros(T, 2)
        new(A_mat, b_vec, ABCS, A_mat2, b_vec2, AA, BB, CC, SS)
    end
end
function Aux{T}(sys::System) where {T}
    Nu, Nx, Ny = size(sys)
    nt = Threads.nthreads()
    # pre-allocate thread-local aux quantities
    [Aux{T}(Nu) for _ in 1:nt]
end

struct BC{T}
    S    :: Array{T,2}
    Fx   :: Array{T,2}
    Fy   :: Array{T,2}
    B1d  :: Array{T,2}
    B2d  :: Array{T,2}
    Gd   :: Array{T,2}
    phid :: Array{T,2}
    Sd   :: Array{T,2}
    A    :: Array{T,2}
    S_u  :: Array{T,2}
    Fx_u :: Array{T,2}
    Fy_u :: Array{T,2}
    A_u  :: Array{T,2}
end
function BC{T}(Nx::Int, Ny::Int) where {T<:Real}
    S    = Array{T}(undef, Nx, Ny)
    Fx   = Array{T}(undef, Nx, Ny)
    Fy   = Array{T}(undef, Nx, Ny)
    B1d  = Array{T}(undef, Nx, Ny)
    B2d  = Array{T}(undef, Nx, Ny)
    Gd   = Array{T}(undef, Nx, Ny)
    phid = Array{T}(undef, Nx, Ny)
    Sd   = Array{T}(undef, Nx, Ny)
    A    = Array{T}(undef, Nx, Ny)
    S_u  = Array{T}(undef, Nx, Ny)
    Fx_u = Array{T}(undef, Nx, Ny)
    Fy_u = Array{T}(undef, Nx, Ny)
    A_u  = Array{T}(undef, Nx, Ny)
    BC{T}(S, Fx, Fy, B1d, B2d, Gd, phid, Sd, A, S_u, Fx_u, Fy_u, A_u)
end

function BCs(systems::SystemPartition)
    sys1  = systems[1]
    T     = Jecco.coord_eltype(sys1.ucoord)
    _, Nx, Ny = size(sys1)

    [BC{T}(Nx, Ny) for sys in systems]
end

#= Notes

for each metric function there are radial ODEs at each z point. since they are
x,y-independent, they can all be solved independently and simultaneously. this is
achieved via the trivial Threads.@thread parallelisation below.

the matrix A_mat is obtained, for each equation, through

  A_mat = A D_uu + B D_u + C Id

which builds the differential operator, and then the top and bottom lines are
replaced to enforce the boundary conditions, as shown below. the vector b_vec is
simply built through

  b_vec = -S

the first and last entries are replaced to enforce the boundary conditions,
leading to the schematic form

  (  1   0   0  ...  0  ) ( x0 )   (  u0  )
  ( a10 a11 a12 ... a1N ) ( x1 )   (  b1  )
  ( ................... ) ( x2 ) = (  b2  )                               (*)
  ( ................... ) ( .. )   (  ..  )
  ( d00 d01 d02 ... d0N ) ( xN )   (  u'0 )

where N = Nu and d00, d01, d02, ..., d0N are the coefficients of the first
line in the first derivative operator D_u:

         ( d00 d01 d02 d03 ... d0N )
         ( d10 d11 d12 d13 ... d1N )
  D_u =  ( d20 d21 d22 d23 ... d2N )
         ( ....................... )
         ( dN0 dN1 dN2 dN3 ... dNN )

the first equation from the system (*) then enforces x0 = u0 while the last
equation enforces dx/du_{u=u0} = u'0, where u0 and u'0 need to be given for each
equation. the remaining equations from (*) enforce the differential equations
themselves.

note that this example is valid only for the second order ODEs. for the first
order ones, obviously, just one BC is needed. thus, we accordingly skip the step
of replacing the last line of the A_mat matrix and last entry of b_vec vector.

=#

function solve_S!(bulk::Bulk, bc::BC, gauge::Gauge, deriv::BulkDeriv, aux_acc,
                  sys::System, evoleq::AffineNull)
    Nu, Nx, Ny = size(sys)

    Du_B1   = deriv.Du_B1
    Du_B2   = deriv.Du_B2
    Du_G    = deriv.Du_G
    Du_phi  = deriv.Du_phi
    # Duu_B1  = deriv.Duu_B1
    # Duu_B2  = deriv.Duu_B2
    # Duu_G   = deriv.Duu_G

    Du  = sys.Du
    Duu = sys.Duu
    # Dx  = sys.Dx
    # Dxx = sys.Dxx
    # Dy  = sys.Dy
    # Dyy = sys.Dyy

    phi0  = evoleq.phi0

    @fastmath @inbounds @threads for j in 1:Ny
        @inbounds for i in 1:Nx
            id  = Threads.threadid()
            aux = aux_acc[id]

            xi  = gauge.xi[1,i,j]

            @inbounds @simd for a in 1:Nu
                u     = sys.ucoord[a]

                B1    = bulk.B1[a,i,j]
                B2    = bulk.B2[a,i,j]
                G     = bulk.G[a,i,j]
                phi   = bulk.phi[a,i,j]

                B1p   = -u*u * Du_B1[a,i,j]
                B2p   = -u*u * Du_B2[a,i,j]
                Gp    = -u*u * Du_G[a,i,j]
                phip  = -u*u * Du_phi[a,i,j]

                vars = (phi0, u, xi, B1, B1p, B2, B2p, G, Gp, phi, phip)

                S_eq_coeff!(aux.ABCS, vars, sys.gridtype)

                aux.b_vec[a]   = -aux.ABCS[4]
                @inbounds @simd for aa in 1:Nu
                    aux.A_mat[a,aa] = aux.ABCS[1] * Duu[a,aa] + aux.ABCS[2] * Du[a,aa]
                end
                aux.A_mat[a,a] += aux.ABCS[3]
            end

            # BC

            aux.b_vec[1]    = bc.S[i,j]
            aux.A_mat[1,:] .= 0.0
            aux.A_mat[1,1]  = 1.0

            aux.b_vec[end]    = bc.S_u[i,j]
            @inbounds @simd for aa in 1:Nu
                aux.A_mat[end,aa]  = Du[1,aa]
            end

            solve_lin_system!(aux.A_mat, aux.b_vec)

            @inbounds @simd for aa in 1:Nu
                bulk.S[aa,i,j] = aux.b_vec[aa]
            end

        end
    end

    nothing
end

function solve_Fxy!(bulk::Bulk, bc::BC, gauge::Gauge, deriv::BulkDeriv, aux_acc,
                    sys::System, evoleq::AffineNull)
    Nu, Nx, Ny = size(sys)

    Du_B1   = deriv.Du_B1
    Du_B2   = deriv.Du_B2
    Du_G    = deriv.Du_G
    Du_phi  = deriv.Du_phi
    Du_S    = deriv.Du_S
    Duu_B1  = deriv.Duu_B1
    Duu_B2  = deriv.Duu_B2
    Duu_G   = deriv.Duu_G
    Duu_phi = deriv.Duu_phi
    Duu_S   = deriv.Duu_S

    Du  = sys.Du
    Duu = sys.Duu
    Dx  = sys.Dx
    Dxx = sys.Dxx
    Dy  = sys.Dy
    Dyy = sys.Dyy

    phi0  = evoleq.phi0

    @fastmath @inbounds @threads for j in 1:Ny
        @inbounds for i in 1:Nx
            id  = Threads.threadid()
            aux = aux_acc[id]

            xi    = gauge.xi[1,i,j]
            xi_x  = Dx(gauge.xi, 1,i,j)
            xi_y  = Dy(gauge.xi, 1,i,j)

            @inbounds @simd for a in 1:Nu
                u     = sys.ucoord[a]
                u2    = u * u
                u3    = u * u2
                u4    = u2 * u2

                B1    = bulk.B1[a,i,j]
                B1p   = -u2 * Du_B1[a,i,j]
                B1_x  = Dx(bulk.B1, a,i,j)
                B1_y  = Dy(bulk.B1, a,i,j)
                B1pp  = 2*u3 * Du_B1[a,i,j] + u4 * Duu_B1[a,i,j]
                B1p_x = -u2 * Dx(Du_B1, a,i,j)
                B1p_y = -u2 * Dy(Du_B1, a,i,j)

                B2    = bulk.B2[a,i,j]
                B2p   = -u2 * Du_B2[a,i,j]
                B2_x  = Dx(bulk.B2, a,i,j)
                B2_y  = Dy(bulk.B2, a,i,j)
                B2pp  = 2*u3 * Du_B2[a,i,j] + u4 * Duu_B2[a,i,j]
                B2p_x = -u2 * Dx(Du_B2, a,i,j)
                B2p_y = -u2 * Dy(Du_B2, a,i,j)

                G     = bulk.G[a,i,j]
                Gp    = -u2 * Du_G[a,i,j]
                G_x   = Dx(bulk.G, a,i,j)
                G_y   = Dy(bulk.G, a,i,j)
                Gpp   = 2*u3 * Du_G[a,i,j] + u4 * Duu_G[a,i,j]
                Gp_x  = -u2 * Dx(Du_G, a,i,j)
                Gp_y  = -u2 * Dy(Du_G, a,i,j)

                phi   = bulk.phi[a,i,j]
                phip  = -u2 * Du_phi[a,i,j]
                phi_x = Dx(bulk.phi, a,i,j)
                phi_y = Dy(bulk.phi, a,i,j)
                # phipp   = 2*u3 * Du_phi[a,i,j] + u4 * Duu_phi[a,i,j]
                # phip_x  = -u2 * Dx(Du_phi, a,i,j)
                # phip_y  = -u2 * Dy(Du_phi, a,i,j)

                S     = bulk.S[a,i,j]
                Sp    = -u2 * Du_S[a,i,j]
                S_x   = Dx(bulk.S, a,i,j)
                S_y   = Dy(bulk.S, a,i,j)
                Spp   = 2*u3 * Du_S[a,i,j] + u4 * Duu_S[a,i,j]
                Sp_x  = -u2 * Dx(Du_S, a,i,j)
                Sp_y  = -u2 * Dy(Du_S, a,i,j)

                vars = (
                    phi0, u, xi, xi_x, xi_y,
                    B1     ,    B2     ,    G      ,    phi    ,    S      ,
                    B1p    ,    B2p    ,    Gp     ,    phip   ,    Sp     ,
                    B1pp   ,    B2pp   ,    Gpp    ,                Spp    ,
                    B1_x   ,    B2_x   ,    G_x    ,    phi_x  ,    S_x    ,
                    B1_y   ,    B2_y   ,    G_y    ,    phi_y  ,    S_y    ,
                    B1p_x  ,    B2p_x  ,    Gp_x   ,                Sp_x   ,
                    B1p_y  ,    B2p_y  ,    Gp_y   ,                Sp_y
                )

                Fxy_eq_coeff!(aux.AA, aux.BB, aux.CC, aux.SS, vars, sys.gridtype)

                aux.b_vec2[a]    = -aux.SS[1]
                aux.b_vec2[a+Nu] = -aux.SS[2]
                @inbounds @simd for aa in 1:Nu
                    aux.A_mat2[a,aa]         = aux.AA[1,1] * Duu[a,aa] + aux.BB[1,1] * Du[a,aa]
                    aux.A_mat2[a,aa+Nu]      = aux.AA[1,2] * Duu[a,aa] + aux.BB[1,2] * Du[a,aa]
                    aux.A_mat2[a+Nu,aa]      = aux.AA[2,1] * Duu[a,aa] + aux.BB[2,1] * Du[a,aa]
                    aux.A_mat2[a+Nu,aa+Nu]   = aux.AA[2,2] * Duu[a,aa] + aux.BB[2,2] * Du[a,aa]
                end
                aux.A_mat2[a,a]       += aux.CC[1,1]
                aux.A_mat2[a,a+Nu]    += aux.CC[1,2]
                aux.A_mat2[a+Nu,a]    += aux.CC[2,1]
                aux.A_mat2[a+Nu,a+Nu] += aux.CC[2,2]
            end

            # BC

            aux.b_vec2[1]     = bc.Fx[i,j]
            aux.A_mat2[1,:]  .= 0.0
            aux.A_mat2[1,1]   = 1.0

            aux.b_vec2[1+Nu]      = bc.Fy[i,j]
            aux.A_mat2[1+Nu,:]   .= 0.0
            aux.A_mat2[1+Nu,1+Nu] = 1.0

            aux.b_vec2[Nu]   = bc.Fx_u[i,j]
            aux.b_vec2[2*Nu] = bc.Fy_u[i,j]

            aux.A_mat2[Nu,:]   .= 0.0
            aux.A_mat2[2*Nu,:] .= 0.0
            @inbounds @simd for aa in 1:Nu
                aux.A_mat2[Nu,aa]      = Du[1,aa]
                aux.A_mat2[2*Nu,aa+Nu] = Du[1,aa]
            end

            solve_lin_system_pivot!(aux.A_mat2, aux.b_vec2)

            @inbounds @simd for aa in 1:Nu
                bulk.Fx[aa,i,j] = aux.b_vec2[aa]
                bulk.Fy[aa,i,j] = aux.b_vec2[aa+Nu]
            end

        end
    end

    nothing
end

function solve_Sd!(bulk::Bulk, bc::BC, gauge::Gauge, deriv::BulkDeriv, aux_acc,
                   sys::System, evoleq::AffineNull)
    Nu, Nx, Ny = size(sys)

    Du_B1   = deriv.Du_B1
    Du_B2   = deriv.Du_B2
    Du_G    = deriv.Du_G
    Du_phi  = deriv.Du_phi
    Du_S    = deriv.Du_S
    Du_Fx   = deriv.Du_Fx
    Du_Fy   = deriv.Du_Fy
    Duu_B1  = deriv.Duu_B1
    Duu_B2  = deriv.Duu_B2
    Duu_G   = deriv.Duu_G
    Duu_phi = deriv.Duu_phi
    Duu_S   = deriv.Duu_S
    Duu_Fx  = deriv.Duu_Fx
    Duu_Fy  = deriv.Duu_Fy

    Du  = sys.Du
    Duu = sys.Duu
    Dx  = sys.Dx
    Dxx = sys.Dxx
    Dy  = sys.Dy
    Dyy = sys.Dyy

    potential = evoleq.potential
    phi0      = evoleq.phi0

    @fastmath @inbounds @threads for j in 1:Ny
        @inbounds for i in 1:Nx
            id   = Threads.threadid()
            aux  = aux_acc[id]

            xi    = gauge.xi[1,i,j]
            xi_x  = Dx(gauge.xi, 1,i,j)
            xi_y  = Dy(gauge.xi, 1,i,j)
            xi_xx = Dxx(gauge.xi, 1,i,j)
            xi_yy = Dyy(gauge.xi, 1,i,j)
            xi_xy = Dx(Dy, gauge.xi, 1,i,j)

            @inbounds @simd for a in 1:Nu
                u     = sys.ucoord[a]
                u2    = u * u
                u3    = u * u2
                u4    = u2 * u2

                B1    = bulk.B1[a,i,j]
                B2    = bulk.B2[a,i,j]
                G     = bulk.G[a,i,j]
                phi   = bulk.phi[a,i,j]
                S     = bulk.S[a,i,j]
                Fx    = bulk.Fx[a,i,j]
                Fy    = bulk.Fy[a,i,j]

                # r derivatives

                B1p        = -u2 * Du_B1[a,i,j]
                B2p        = -u2 * Du_B2[a,i,j]
                Gp         = -u2 * Du_G[a,i,j]
                phip       = -u2 * Du_phi[a,i,j]
                Sp         = -u2 * Du_S[a,i,j]
                Fxp        = -u2 * Du_Fx[a,i,j]
                Fyp        = -u2 * Du_Fy[a,i,j]

                B1pp       = 2*u3 * Du_B1[a,i,j]  + u4 * Duu_B1[a,i,j]
                B2pp       = 2*u3 * Du_B2[a,i,j]  + u4 * Duu_B2[a,i,j]
                Gpp        = 2*u3 * Du_G[a,i,j]   + u4 * Duu_G[a,i,j]
                phipp      = 2*u3 * Du_phi[a,i,j] + u4 * Duu_phi[a,i,j]
                Spp        = 2*u3 * Du_S[a,i,j]   + u4 * Duu_S[a,i,j]
                Fxpp       = 2*u3 * Du_Fx[a,i,j]  + u4 * Duu_Fx[a,i,j]
                Fypp       = 2*u3 * Du_Fy[a,i,j]  + u4 * Duu_Fy[a,i,j]

                # x and y derivatives

                B1_x       = Dx(bulk.B1, a,i,j)
                B2_x       = Dx(bulk.B2, a,i,j)
                G_x        = Dx(bulk.G,  a,i,j)
                phi_x      = Dx(bulk.phi,a,i,j)
                S_x        = Dx(bulk.S,  a,i,j)
                Fx_x       = Dx(bulk.Fx, a,i,j)
                Fy_x       = Dx(bulk.Fy, a,i,j)

                B1_y       = Dy(bulk.B1, a,i,j)
                B2_y       = Dy(bulk.B2, a,i,j)
                G_y        = Dy(bulk.G,  a,i,j)
                phi_y      = Dy(bulk.phi,a,i,j)
                S_y        = Dy(bulk.S,  a,i,j)
                Fx_y       = Dy(bulk.Fx, a,i,j)
                Fy_y       = Dy(bulk.Fy, a,i,j)

                B1p_x      = -u2 * Dx(Du_B1, a,i,j)
                B2p_x      = -u2 * Dx(Du_B2, a,i,j)
                Gp_x       = -u2 * Dx(Du_G,  a,i,j)
                phip_x     = -u2 * Dx(Du_phi,a,i,j)
                Sp_x       = -u2 * Dx(Du_S,  a,i,j)
                Fxp_x      = -u2 * Dx(Du_Fx,  a,i,j)
                Fyp_x      = -u2 * Dx(Du_Fy,  a,i,j)

                B1p_y      = -u2 * Dy(Du_B1, a,i,j)
                B2p_y      = -u2 * Dy(Du_B2, a,i,j)
                Gp_y       = -u2 * Dy(Du_G,  a,i,j)
                phip_y     = -u2 * Dy(Du_phi,a,i,j)
                Sp_y       = -u2 * Dy(Du_S,  a,i,j)
                Fxp_y      = -u2 * Dy(Du_Fx, a,i,j)
                Fyp_y      = -u2 * Dy(Du_Fy, a,i,j)

                B1_xx      = Dxx(bulk.B1, a,i,j)
                B2_xx      = Dxx(bulk.B2, a,i,j)
                G_xx       = Dxx(bulk.G,  a,i,j)
                phi_xx     = Dxx(bulk.phi,a,i,j)
                S_xx       = Dxx(bulk.S,  a,i,j)

                B1_yy      = Dyy(bulk.B1, a,i,j)
                B2_yy      = Dyy(bulk.B2, a,i,j)
                G_yy       = Dyy(bulk.G,  a,i,j)
                phi_yy     = Dyy(bulk.phi,a,i,j)
                S_yy       = Dyy(bulk.S,  a,i,j)

                B2_xy      = Dx(Dy, bulk.B2, a,i,j)
                G_xy       = Dx(Dy, bulk.G,  a,i,j)
                S_xy       = Dx(Dy, bulk.S,  a,i,j)

                vars = (
                    potential, phi0, u, xi, xi_x, xi_y, xi_xx, xi_yy, xi_xy,
                    B1     ,    B2     ,    G      ,    phi    ,    S      ,    Fx     ,    Fy     ,
                    B1p    ,    B2p    ,    Gp     ,    phip   ,    Sp     ,    Fxp    ,    Fyp    ,
                    B1pp   ,    B2pp   ,    Gpp    ,    phipp  ,    Spp    ,    Fxpp   ,    Fypp   ,
                    B1_x   ,    B2_x   ,    G_x    ,    phi_x  ,    S_x    ,    Fx_x   ,    Fy_x   ,
                    B1_y   ,    B2_y   ,    G_y    ,    phi_y  ,    S_y    ,    Fx_y   ,    Fy_y   ,
                    B1p_x  ,    B2p_x  ,    Gp_x   ,    phip_x ,    Sp_x   ,    Fxp_x  ,    Fyp_x  ,
                    B1p_y  ,    B2p_y  ,    Gp_y   ,    phip_y ,    Sp_y   ,    Fxp_y  ,    Fyp_y  ,
                    B1_xx  ,    B2_xx  ,    G_xx   ,    phi_xx ,    S_xx   ,
                    B1_yy  ,    B2_yy  ,    G_yy   ,    phi_yy ,    S_yy   ,
                                B2_xy  ,    G_xy   ,                S_xy
                )

                Sd_eq_coeff!(aux.ABCS, vars, sys.gridtype)

                aux.b_vec[a]   = -aux.ABCS[4]
                @inbounds @simd for aa in 1:Nu
                    aux.A_mat[a,aa] = aux.ABCS[1] * Duu[a,aa] + aux.ABCS[2] * Du[a,aa]
                end
                aux.A_mat[a,a] += aux.ABCS[3]
            end

            # BC (first order equation)

            aux.b_vec[1]    = bc.Sd[i,j]
            aux.A_mat[1,:] .= 0.0
            aux.A_mat[1,1]  = 1.0

            solve_lin_system!(aux.A_mat, aux.b_vec)

            @inbounds @simd for aa in 1:Nu
                bulk.Sd[aa,i,j] = aux.b_vec[aa]
            end

        end
    end

    nothing
end

function solve_B2d!(bulk::Bulk, bc::BC, gauge::Gauge, deriv::BulkDeriv, aux_acc,
                    sys::System, evoleq::AffineNull)
    Nu, Nx, Ny = size(sys)

    Du_B1   = deriv.Du_B1
    Du_B2   = deriv.Du_B2
    Du_G    = deriv.Du_G
    Du_phi  = deriv.Du_phi
    Du_S    = deriv.Du_S
    Du_Fx   = deriv.Du_Fx
    Du_Fy   = deriv.Du_Fy
    Duu_B1  = deriv.Duu_B1
    Duu_B2  = deriv.Duu_B2
    Duu_G   = deriv.Duu_G
    Duu_phi = deriv.Duu_phi
    Duu_S   = deriv.Duu_S
    Duu_Fx  = deriv.Duu_Fx
    Duu_Fy  = deriv.Duu_Fy

    Du  = sys.Du
    Duu = sys.Duu
    Dx  = sys.Dx
    Dxx = sys.Dxx
    Dy  = sys.Dy
    Dyy = sys.Dyy

    phi0  = evoleq.phi0

    @fastmath @inbounds @threads for j in 1:Ny
        @inbounds for i in 1:Nx
            id  = Threads.threadid()
            aux = aux_acc[id]

            xi    = gauge.xi[1,i,j]
            xi_x  = Dx(gauge.xi, 1,i,j)
            xi_y  = Dy(gauge.xi, 1,i,j)
            xi_xx = Dxx(gauge.xi, 1,i,j)
            xi_yy = Dyy(gauge.xi, 1,i,j)
            xi_xy = Dx(Dy, gauge.xi, 1,i,j)

            @inbounds @simd for a in 1:Nu
                u     = sys.ucoord[a]
                u2    = u * u
                u3    = u * u2
                u4    = u2 * u2

                B1    = bulk.B1[a,i,j]
                B2    = bulk.B2[a,i,j]
                G     = bulk.G[a,i,j]
                phi   = bulk.phi[a,i,j]
                S     = bulk.S[a,i,j]
                Fx    = bulk.Fx[a,i,j]
                Fy    = bulk.Fy[a,i,j]
                Sd    = bulk.Sd[a,i,j]

                # r derivatives

                B1p        = -u2 * Du_B1[a,i,j]
                B2p        = -u2 * Du_B2[a,i,j]
                Gp         = -u2 * Du_G[a,i,j]
                phip       = -u2 * Du_phi[a,i,j]
                Sp         = -u2 * Du_S[a,i,j]
                Fxp        = -u2 * Du_Fx[a,i,j]
                Fyp        = -u2 * Du_Fy[a,i,j]

                B1pp       = 2*u3 * Du_B1[a,i,j]  + u4 * Duu_B1[a,i,j]
                B2pp       = 2*u3 * Du_B2[a,i,j]  + u4 * Duu_B2[a,i,j]
                Gpp        = 2*u3 * Du_G[a,i,j]   + u4 * Duu_G[a,i,j]
                phipp      = 2*u3 * Du_phi[a,i,j] + u4 * Duu_phi[a,i,j]
                Spp        = 2*u3 * Du_S[a,i,j]   + u4 * Duu_S[a,i,j]
                Fxpp       = 2*u3 * Du_Fx[a,i,j]  + u4 * Duu_Fx[a,i,j]
                Fypp       = 2*u3 * Du_Fy[a,i,j]  + u4 * Duu_Fy[a,i,j]

                # x and y derivatives

                B1_x       = Dx(bulk.B1, a,i,j)
                B2_x       = Dx(bulk.B2, a,i,j)
                G_x        = Dx(bulk.G,  a,i,j)
                phi_x      = Dx(bulk.phi,a,i,j)
                S_x        = Dx(bulk.S,  a,i,j)
                Fx_x       = Dx(bulk.Fx, a,i,j)
                Fy_x       = Dx(bulk.Fy, a,i,j)

                B1_y       = Dy(bulk.B1, a,i,j)
                B2_y       = Dy(bulk.B2, a,i,j)
                G_y        = Dy(bulk.G,  a,i,j)
                phi_y      = Dy(bulk.phi,a,i,j)
                S_y        = Dy(bulk.S,  a,i,j)
                Fx_y       = Dy(bulk.Fx, a,i,j)
                Fy_y       = Dy(bulk.Fy, a,i,j)

                B1p_x      = -u2 * Dx(Du_B1, a,i,j)
                B2p_x      = -u2 * Dx(Du_B2, a,i,j)
                Gp_x       = -u2 * Dx(Du_G,  a,i,j)
                phip_x     = -u2 * Dx(Du_phi,a,i,j)
                Sp_x       = -u2 * Dx(Du_S,  a,i,j)
                Fxp_x      = -u2 * Dx(Du_Fx,  a,i,j)
                Fyp_x      = -u2 * Dx(Du_Fy,  a,i,j)

                B1p_y      = -u2 * Dy(Du_B1, a,i,j)
                B2p_y      = -u2 * Dy(Du_B2, a,i,j)
                Gp_y       = -u2 * Dy(Du_G,  a,i,j)
                phip_y     = -u2 * Dy(Du_phi,a,i,j)
                Sp_y       = -u2 * Dy(Du_S,  a,i,j)
                Fxp_y      = -u2 * Dy(Du_Fx, a,i,j)
                Fyp_y      = -u2 * Dy(Du_Fy, a,i,j)

                B1_xx      = Dxx(bulk.B1, a,i,j)
                B2_xx      = Dxx(bulk.B2, a,i,j)
                G_xx       = Dxx(bulk.G,  a,i,j)
                phi_xx     = Dxx(bulk.phi,a,i,j)
                S_xx       = Dxx(bulk.S,  a,i,j)

                B1_yy      = Dyy(bulk.B1, a,i,j)
                B2_yy      = Dyy(bulk.B2, a,i,j)
                G_yy       = Dyy(bulk.G,  a,i,j)
                phi_yy     = Dyy(bulk.phi,a,i,j)
                S_yy       = Dyy(bulk.S,  a,i,j)

                B2_xy      = Dx(Dy, bulk.B2, a,i,j)
                G_xy       = Dx(Dy, bulk.G,  a,i,j)
                S_xy       = Dx(Dy, bulk.S,  a,i,j)

                vars = (
                    phi0, u, xi, xi_x, xi_y, xi_xx, xi_yy, xi_xy,
                    B1     ,    B2     ,    G      ,    phi    ,    S      ,    Fx     ,    Fy     ,  Sd,
                    B1p    ,    B2p    ,    Gp     ,    phip   ,    Sp     ,    Fxp    ,    Fyp    ,
                    B1pp   ,    B2pp   ,    Gpp    ,    phipp  ,    Spp    ,    Fxpp   ,    Fypp   ,
                    B1_x   ,    B2_x   ,    G_x    ,    phi_x  ,    S_x    ,    Fx_x   ,    Fy_x   ,
                    B1_y   ,    B2_y   ,    G_y    ,    phi_y  ,    S_y    ,    Fx_y   ,    Fy_y   ,
                    B1p_x  ,    B2p_x  ,    Gp_x   ,    phip_x ,    Sp_x   ,    Fxp_x  ,    Fyp_x  ,
                    B1p_y  ,    B2p_y  ,    Gp_y   ,    phip_y ,    Sp_y   ,    Fxp_y  ,    Fyp_y  ,
                    B1_xx  ,    B2_xx  ,    G_xx   ,    phi_xx ,    S_xx   ,
                    B1_yy  ,    B2_yy  ,    G_yy   ,    phi_yy ,    S_yy   ,
                                B2_xy  ,    G_xy   ,                S_xy
                )

                B2d_eq_coeff!(aux.ABCS, vars, sys.gridtype)

                aux.b_vec[a]   = -aux.ABCS[4]
                @inbounds @simd for aa in 1:Nu
                    aux.A_mat[a,aa] = aux.ABCS[1] * Duu[a,aa] + aux.ABCS[2] * Du[a,aa]
                end
                aux.A_mat[a,a] += aux.ABCS[3]
            end

            # BC (first order equation)

            aux.b_vec[1]    = bc.B2d[i,j]
            aux.A_mat[1,:] .= 0.0
            aux.A_mat[1,1]  = 1.0

            solve_lin_system!(aux.A_mat, aux.b_vec)

            @inbounds @simd for aa in 1:Nu
                bulk.B2d[aa,i,j] = aux.b_vec[aa]
            end

        end
    end

    nothing
end

function solve_B1dGd!(bulk::Bulk, bc::BC, gauge::Gauge, deriv::BulkDeriv, aux_acc,
                      sys::System, evoleq::AffineNull)
    Nu, Nx, Ny = size(sys)

    Du_B1   = deriv.Du_B1
    Du_B2   = deriv.Du_B2
    Du_G    = deriv.Du_G
    Du_phi  = deriv.Du_phi
    Du_S    = deriv.Du_S
    Du_Fx   = deriv.Du_Fx
    Du_Fy   = deriv.Du_Fy
    Duu_B1  = deriv.Duu_B1
    Duu_B2  = deriv.Duu_B2
    Duu_G   = deriv.Duu_G
    Duu_phi = deriv.Duu_phi
    Duu_S   = deriv.Duu_S
    Duu_Fx  = deriv.Duu_Fx
    Duu_Fy  = deriv.Duu_Fy

    Du  = sys.Du
    Duu = sys.Duu
    Dx  = sys.Dx
    Dxx = sys.Dxx
    Dy  = sys.Dy
    Dyy = sys.Dyy

    phi0  = evoleq.phi0

    @fastmath @inbounds @threads for j in 1:Ny
        @inbounds for i in 1:Nx
            id  = Threads.threadid()
            aux = aux_acc[id]

            xi    = gauge.xi[1,i,j]
            xi_x  = Dx(gauge.xi, 1,i,j)
            xi_y  = Dy(gauge.xi, 1,i,j)
            xi_xx = Dxx(gauge.xi, 1,i,j)
            xi_yy = Dyy(gauge.xi, 1,i,j)
            xi_xy = Dx(Dy, gauge.xi, 1,i,j)

            @inbounds @simd for a in 1:Nu
                u     = sys.ucoord[a]
                u2    = u * u
                u3    = u * u2
                u4    = u2 * u2

                B1    = bulk.B1[a,i,j]
                B2    = bulk.B2[a,i,j]
                G     = bulk.G[a,i,j]
                phi   = bulk.phi[a,i,j]
                S     = bulk.S[a,i,j]
                Fx    = bulk.Fx[a,i,j]
                Fy    = bulk.Fy[a,i,j]
                Sd    = bulk.Sd[a,i,j]

                # r derivatives

                B1p        = -u2 * Du_B1[a,i,j]
                B2p        = -u2 * Du_B2[a,i,j]
                Gp         = -u2 * Du_G[a,i,j]
                phip       = -u2 * Du_phi[a,i,j]
                Sp         = -u2 * Du_S[a,i,j]
                Fxp        = -u2 * Du_Fx[a,i,j]
                Fyp        = -u2 * Du_Fy[a,i,j]

                B1pp       = 2*u3 * Du_B1[a,i,j]  + u4 * Duu_B1[a,i,j]
                B2pp       = 2*u3 * Du_B2[a,i,j]  + u4 * Duu_B2[a,i,j]
                Gpp        = 2*u3 * Du_G[a,i,j]   + u4 * Duu_G[a,i,j]
                phipp      = 2*u3 * Du_phi[a,i,j] + u4 * Duu_phi[a,i,j]
                Spp        = 2*u3 * Du_S[a,i,j]   + u4 * Duu_S[a,i,j]
                Fxpp       = 2*u3 * Du_Fx[a,i,j]  + u4 * Duu_Fx[a,i,j]
                Fypp       = 2*u3 * Du_Fy[a,i,j]  + u4 * Duu_Fy[a,i,j]

                # x and y derivatives

                B1_x       = Dx(bulk.B1, a,i,j)
                B2_x       = Dx(bulk.B2, a,i,j)
                G_x        = Dx(bulk.G,  a,i,j)
                phi_x      = Dx(bulk.phi,a,i,j)
                S_x        = Dx(bulk.S,  a,i,j)
                Fx_x       = Dx(bulk.Fx, a,i,j)
                Fy_x       = Dx(bulk.Fy, a,i,j)

                B1_y       = Dy(bulk.B1, a,i,j)
                B2_y       = Dy(bulk.B2, a,i,j)
                G_y        = Dy(bulk.G,  a,i,j)
                phi_y      = Dy(bulk.phi,a,i,j)
                S_y        = Dy(bulk.S,  a,i,j)
                Fx_y       = Dy(bulk.Fx, a,i,j)
                Fy_y       = Dy(bulk.Fy, a,i,j)

                B1p_x      = -u2 * Dx(Du_B1, a,i,j)
                B2p_x      = -u2 * Dx(Du_B2, a,i,j)
                Gp_x       = -u2 * Dx(Du_G,  a,i,j)
                phip_x     = -u2 * Dx(Du_phi,a,i,j)
                Sp_x       = -u2 * Dx(Du_S,  a,i,j)
                Fxp_x      = -u2 * Dx(Du_Fx,  a,i,j)
                Fyp_x      = -u2 * Dx(Du_Fy,  a,i,j)

                B1p_y      = -u2 * Dy(Du_B1, a,i,j)
                B2p_y      = -u2 * Dy(Du_B2, a,i,j)
                Gp_y       = -u2 * Dy(Du_G,  a,i,j)
                phip_y     = -u2 * Dy(Du_phi,a,i,j)
                Sp_y       = -u2 * Dy(Du_S,  a,i,j)
                Fxp_y      = -u2 * Dy(Du_Fx, a,i,j)
                Fyp_y      = -u2 * Dy(Du_Fy, a,i,j)

                B1_xx      = Dxx(bulk.B1, a,i,j)
                B2_xx      = Dxx(bulk.B2, a,i,j)
                G_xx       = Dxx(bulk.G,  a,i,j)
                phi_xx     = Dxx(bulk.phi,a,i,j)
                S_xx       = Dxx(bulk.S,  a,i,j)

                B1_yy      = Dyy(bulk.B1, a,i,j)
                B2_yy      = Dyy(bulk.B2, a,i,j)
                G_yy       = Dyy(bulk.G,  a,i,j)
                phi_yy     = Dyy(bulk.phi,a,i,j)
                S_yy       = Dyy(bulk.S,  a,i,j)

                B2_xy      = Dx(Dy, bulk.B2, a,i,j)
                G_xy       = Dx(Dy, bulk.G,  a,i,j)
                S_xy       = Dx(Dy, bulk.S,  a,i,j)

                vars = (
                    phi0, u, xi, xi_x, xi_y, xi_xx, xi_yy, xi_xy,
                    B1     ,    B2     ,    G      ,    phi    ,    S      ,    Fx     ,    Fy     ,  Sd,
                    B1p    ,    B2p    ,    Gp     ,    phip   ,    Sp     ,    Fxp    ,    Fyp    ,
                    B1pp   ,    B2pp   ,    Gpp    ,    phipp  ,    Spp    ,    Fxpp   ,    Fypp   ,
                    B1_x   ,    B2_x   ,    G_x    ,    phi_x  ,    S_x    ,    Fx_x   ,    Fy_x   ,
                    B1_y   ,    B2_y   ,    G_y    ,    phi_y  ,    S_y    ,    Fx_y   ,    Fy_y   ,
                    B1p_x  ,    B2p_x  ,    Gp_x   ,    phip_x ,    Sp_x   ,    Fxp_x  ,    Fyp_x  ,
                    B1p_y  ,    B2p_y  ,    Gp_y   ,    phip_y ,    Sp_y   ,    Fxp_y  ,    Fyp_y  ,
                    B1_xx  ,    B2_xx  ,    G_xx   ,    phi_xx ,    S_xx   ,
                    B1_yy  ,    B2_yy  ,    G_yy   ,    phi_yy ,    S_yy   ,
                                B2_xy  ,    G_xy   ,                S_xy
                )

                B1dGd_eq_coeff!(aux.AA, aux.BB, aux.CC, aux.SS, vars, sys.gridtype)

                aux.b_vec2[a]    = -aux.SS[1]
                aux.b_vec2[a+Nu] = -aux.SS[2]
                @inbounds @simd for aa in 1:Nu
                    aux.A_mat2[a,aa]         = aux.AA[1,1] * Duu[a,aa] + aux.BB[1,1] * Du[a,aa]
                    aux.A_mat2[a,aa+Nu]      = aux.AA[1,2] * Duu[a,aa] + aux.BB[1,2] * Du[a,aa]
                    aux.A_mat2[a+Nu,aa]      = aux.AA[2,1] * Duu[a,aa] + aux.BB[2,1] * Du[a,aa]
                    aux.A_mat2[a+Nu,aa+Nu]   = aux.AA[2,2] * Duu[a,aa] + aux.BB[2,2] * Du[a,aa]
                end
                aux.A_mat2[a,a]       += aux.CC[1,1]
                aux.A_mat2[a,a+Nu]    += aux.CC[1,2]
                aux.A_mat2[a+Nu,a]    += aux.CC[2,1]
                aux.A_mat2[a+Nu,a+Nu] += aux.CC[2,2]
            end

            # BC (first order system)

            aux.b_vec2[1]     = bc.B1d[i,j]
            aux.A_mat2[1,:]  .= 0.0
            aux.A_mat2[1,1]   = 1.0

            aux.b_vec2[1+Nu]      = bc.Gd[i,j]
            aux.A_mat2[1+Nu,:]   .= 0.0
            aux.A_mat2[1+Nu,1+Nu] = 1.0

            solve_lin_system_pivot!(aux.A_mat2, aux.b_vec2)

            @inbounds @simd for aa in 1:Nu
                bulk.B1d[aa,i,j] = aux.b_vec2[aa]
                bulk.Gd[aa,i,j]  = aux.b_vec2[aa+Nu]
            end

        end
    end

    nothing
end

function solve_phid!(bulk::Bulk, bc::BC, gauge::Gauge, deriv::BulkDeriv, aux_acc,
                     sys::System, evoleq::AffineNull)
    Nu, Nx, Ny = size(sys)

    Du_B1   = deriv.Du_B1
    Du_B2   = deriv.Du_B2
    Du_G    = deriv.Du_G
    Du_phi  = deriv.Du_phi
    Du_S    = deriv.Du_S
    Du_Fx   = deriv.Du_Fx
    Du_Fy   = deriv.Du_Fy
    Duu_B1  = deriv.Duu_B1
    Duu_B2  = deriv.Duu_B2
    Duu_G   = deriv.Duu_G
    Duu_phi = deriv.Duu_phi
    Duu_S   = deriv.Duu_S
    Duu_Fx  = deriv.Duu_Fx
    Duu_Fy  = deriv.Duu_Fy

    Du  = sys.Du
    Duu = sys.Duu
    Dx  = sys.Dx
    Dxx = sys.Dxx
    Dy  = sys.Dy
    Dyy = sys.Dyy

    potential = evoleq.potential
    phi0      = evoleq.phi0

    # if phi0 = 0 set phid to zero and return
    if abs(phi0) < 1e-9
        fill!(bulk.phid, 0)
        return
    end

    @fastmath @inbounds @threads for j in 1:Ny
        @inbounds for i in 1:Nx
            id  = Threads.threadid()
            aux = aux_acc[id]

            xi    = gauge.xi[1,i,j]
            xi_x  = Dx(gauge.xi, 1,i,j)
            xi_y  = Dy(gauge.xi, 1,i,j)
            xi_xx = Dxx(gauge.xi, 1,i,j)
            xi_yy = Dyy(gauge.xi, 1,i,j)
            xi_xy = Dx(Dy, gauge.xi, 1,i,j)

            @inbounds @simd for a in 1:Nu
                u     = sys.ucoord[a]
                u2    = u * u
                u3    = u * u2
                u4    = u2 * u2

                B1    = bulk.B1[a,i,j]
                B2    = bulk.B2[a,i,j]
                G     = bulk.G[a,i,j]
                phi   = bulk.phi[a,i,j]
                S     = bulk.S[a,i,j]
                Fx    = bulk.Fx[a,i,j]
                Fy    = bulk.Fy[a,i,j]
                Sd    = bulk.Sd[a,i,j]

                # r derivatives

                B1p        = -u2 * Du_B1[a,i,j]
                B2p        = -u2 * Du_B2[a,i,j]
                Gp         = -u2 * Du_G[a,i,j]
                phip       = -u2 * Du_phi[a,i,j]
                Sp         = -u2 * Du_S[a,i,j]
                Fxp        = -u2 * Du_Fx[a,i,j]
                Fyp        = -u2 * Du_Fy[a,i,j]

                B1pp       = 2*u3 * Du_B1[a,i,j]  + u4 * Duu_B1[a,i,j]
                B2pp       = 2*u3 * Du_B2[a,i,j]  + u4 * Duu_B2[a,i,j]
                Gpp        = 2*u3 * Du_G[a,i,j]   + u4 * Duu_G[a,i,j]
                phipp      = 2*u3 * Du_phi[a,i,j] + u4 * Duu_phi[a,i,j]
                Spp        = 2*u3 * Du_S[a,i,j]   + u4 * Duu_S[a,i,j]
                Fxpp       = 2*u3 * Du_Fx[a,i,j]  + u4 * Duu_Fx[a,i,j]
                Fypp       = 2*u3 * Du_Fy[a,i,j]  + u4 * Duu_Fy[a,i,j]

                # x and y derivatives

                B1_x       = Dx(bulk.B1, a,i,j)
                B2_x       = Dx(bulk.B2, a,i,j)
                G_x        = Dx(bulk.G,  a,i,j)
                phi_x      = Dx(bulk.phi,a,i,j)
                S_x        = Dx(bulk.S,  a,i,j)
                Fx_x       = Dx(bulk.Fx, a,i,j)
                Fy_x       = Dx(bulk.Fy, a,i,j)

                B1_y       = Dy(bulk.B1, a,i,j)
                B2_y       = Dy(bulk.B2, a,i,j)
                G_y        = Dy(bulk.G,  a,i,j)
                phi_y      = Dy(bulk.phi,a,i,j)
                S_y        = Dy(bulk.S,  a,i,j)
                Fx_y       = Dy(bulk.Fx, a,i,j)
                Fy_y       = Dy(bulk.Fy, a,i,j)

                B1p_x      = -u2 * Dx(Du_B1, a,i,j)
                B2p_x      = -u2 * Dx(Du_B2, a,i,j)
                Gp_x       = -u2 * Dx(Du_G,  a,i,j)
                phip_x     = -u2 * Dx(Du_phi,a,i,j)
                Sp_x       = -u2 * Dx(Du_S,  a,i,j)
                Fxp_x      = -u2 * Dx(Du_Fx,  a,i,j)
                Fyp_x      = -u2 * Dx(Du_Fy,  a,i,j)

                B1p_y      = -u2 * Dy(Du_B1, a,i,j)
                B2p_y      = -u2 * Dy(Du_B2, a,i,j)
                Gp_y       = -u2 * Dy(Du_G,  a,i,j)
                phip_y     = -u2 * Dy(Du_phi,a,i,j)
                Sp_y       = -u2 * Dy(Du_S,  a,i,j)
                Fxp_y      = -u2 * Dy(Du_Fx, a,i,j)
                Fyp_y      = -u2 * Dy(Du_Fy, a,i,j)

                B1_xx      = Dxx(bulk.B1, a,i,j)
                B2_xx      = Dxx(bulk.B2, a,i,j)
                G_xx       = Dxx(bulk.G,  a,i,j)
                phi_xx     = Dxx(bulk.phi,a,i,j)
                S_xx       = Dxx(bulk.S,  a,i,j)

                B1_yy      = Dyy(bulk.B1, a,i,j)
                B2_yy      = Dyy(bulk.B2, a,i,j)
                G_yy       = Dyy(bulk.G,  a,i,j)
                phi_yy     = Dyy(bulk.phi,a,i,j)
                S_yy       = Dyy(bulk.S,  a,i,j)

                B2_xy      = Dx(Dy, bulk.B2, a,i,j)
                G_xy       = Dx(Dy, bulk.G,  a,i,j)
                phi_xy     = Dx(Dy, bulk.phi,a,i,j)
                S_xy       = Dx(Dy, bulk.S,  a,i,j)

                vars = (
                    potential, phi0, u, xi, xi_x, xi_y, xi_xx, xi_yy, xi_xy,
                    B1     ,    B2     ,    G      ,    phi    ,    S      ,    Fx     ,    Fy     ,  Sd,
                    B1p    ,    B2p    ,    Gp     ,    phip   ,    Sp     ,    Fxp    ,    Fyp    ,
                    B1pp   ,    B2pp   ,    Gpp    ,    phipp  ,    Spp    ,    Fxpp   ,    Fypp   ,
                    B1_x   ,    B2_x   ,    G_x    ,    phi_x  ,    S_x    ,    Fx_x   ,    Fy_x   ,
                    B1_y   ,    B2_y   ,    G_y    ,    phi_y  ,    S_y    ,    Fx_y   ,    Fy_y   ,
                    B1p_x  ,    B2p_x  ,    Gp_x   ,    phip_x ,    Sp_x   ,    Fxp_x  ,    Fyp_x  ,
                    B1p_y  ,    B2p_y  ,    Gp_y   ,    phip_y ,    Sp_y   ,    Fxp_y  ,    Fyp_y  ,
                    B1_xx  ,    B2_xx  ,    G_xx   ,    phi_xx ,    S_xx   ,
                    B1_yy  ,    B2_yy  ,    G_yy   ,    phi_yy ,    S_yy   ,
                                B2_xy  ,    G_xy   ,    phi_xy,     S_xy
                )

                phid_eq_coeff!(aux.ABCS, vars, sys.gridtype)

                aux.b_vec[a]   = -aux.ABCS[4]
                @inbounds @simd for aa in 1:Nu
                    aux.A_mat[a,aa] = aux.ABCS[1] * Duu[a,aa] + aux.ABCS[2] * Du[a,aa]
                end
                aux.A_mat[a,a] += aux.ABCS[3]
            end

            # BC (first order equation)

            aux.b_vec[1]    = bc.phid[i,j]
            aux.A_mat[1,:] .= 0.0
            aux.A_mat[1,1]  = 1.0

            solve_lin_system!(aux.A_mat, aux.b_vec)

            @inbounds @simd for aa in 1:Nu
                bulk.phid[aa,i,j] = aux.b_vec[aa]
            end

        end
    end

    nothing
end

function solve_A!(bulk::Bulk, bc::BC, gauge::Gauge, deriv::BulkDeriv, aux_acc,
                  sys::System, evoleq::AffineNull)
    Nu, Nx, Ny = size(sys)

    Du_B1   = deriv.Du_B1
    Du_B2   = deriv.Du_B2
    Du_G    = deriv.Du_G
    Du_phi  = deriv.Du_phi
    Du_S    = deriv.Du_S
    Du_Fx   = deriv.Du_Fx
    Du_Fy   = deriv.Du_Fy
    Duu_B1  = deriv.Duu_B1
    Duu_B2  = deriv.Duu_B2
    Duu_G   = deriv.Duu_G
    Duu_phi = deriv.Duu_phi
    Duu_S   = deriv.Duu_S
    Duu_Fx  = deriv.Duu_Fx
    Duu_Fy  = deriv.Duu_Fy

    Du  = sys.Du
    Duu = sys.Duu
    Dx  = sys.Dx
    Dxx = sys.Dxx
    Dy  = sys.Dy
    Dyy = sys.Dyy

    potential = evoleq.potential
    phi0      = evoleq.phi0

    @fastmath @inbounds @threads for j in 1:Ny
        @inbounds for i in 1:Nx
            id  = Threads.threadid()
            aux = aux_acc[id]

            xi    = gauge.xi[1,i,j]
            xi_x  = Dx(gauge.xi, 1,i,j)
            xi_y  = Dy(gauge.xi, 1,i,j)
            xi_xx = Dxx(gauge.xi, 1,i,j)
            xi_yy = Dyy(gauge.xi, 1,i,j)
            xi_xy = Dx(Dy, gauge.xi, 1,i,j)

            @inbounds @simd for a in 1:Nu
                u     = sys.ucoord[a]
                u2    = u * u
                u3    = u * u2
                u4    = u2 * u2

                B1    = bulk.B1[a,i,j]
                B2    = bulk.B2[a,i,j]
                G     = bulk.G[a,i,j]
                phi   = bulk.phi[a,i,j]
                S     = bulk.S[a,i,j]
                Fx    = bulk.Fx[a,i,j]
                Fy    = bulk.Fy[a,i,j]
                Sd    = bulk.Sd[a,i,j]
                B1d   = bulk.B1d[a,i,j]
                B2d   = bulk.B2d[a,i,j]
                Gd    = bulk.Gd[a,i,j]
                phid  = bulk.phid[a,i,j]

                # r derivatives

                B1p        = -u2 * Du_B1[a,i,j]
                B2p        = -u2 * Du_B2[a,i,j]
                Gp         = -u2 * Du_G[a,i,j]
                phip       = -u2 * Du_phi[a,i,j]
                Sp         = -u2 * Du_S[a,i,j]
                Fxp        = -u2 * Du_Fx[a,i,j]
                Fyp        = -u2 * Du_Fy[a,i,j]

                B1pp       = 2*u3 * Du_B1[a,i,j]  + u4 * Duu_B1[a,i,j]
                B2pp       = 2*u3 * Du_B2[a,i,j]  + u4 * Duu_B2[a,i,j]
                Gpp        = 2*u3 * Du_G[a,i,j]   + u4 * Duu_G[a,i,j]
                phipp      = 2*u3 * Du_phi[a,i,j] + u4 * Duu_phi[a,i,j]
                Spp        = 2*u3 * Du_S[a,i,j]   + u4 * Duu_S[a,i,j]
                Fxpp       = 2*u3 * Du_Fx[a,i,j]  + u4 * Duu_Fx[a,i,j]
                Fypp       = 2*u3 * Du_Fy[a,i,j]  + u4 * Duu_Fy[a,i,j]

                # x and y derivatives

                B1_x       = Dx(bulk.B1, a,i,j)
                B2_x       = Dx(bulk.B2, a,i,j)
                G_x        = Dx(bulk.G,  a,i,j)
                phi_x      = Dx(bulk.phi,a,i,j)
                S_x        = Dx(bulk.S,  a,i,j)
                Fx_x       = Dx(bulk.Fx, a,i,j)
                Fy_x       = Dx(bulk.Fy, a,i,j)

                B1_y       = Dy(bulk.B1, a,i,j)
                B2_y       = Dy(bulk.B2, a,i,j)
                G_y        = Dy(bulk.G,  a,i,j)
                phi_y      = Dy(bulk.phi,a,i,j)
                S_y        = Dy(bulk.S,  a,i,j)
                Fx_y       = Dy(bulk.Fx, a,i,j)
                Fy_y       = Dy(bulk.Fy, a,i,j)

                B1p_x      = -u2 * Dx(Du_B1, a,i,j)
                B2p_x      = -u2 * Dx(Du_B2, a,i,j)
                Gp_x       = -u2 * Dx(Du_G,  a,i,j)
                phip_x     = -u2 * Dx(Du_phi,a,i,j)
                Sp_x       = -u2 * Dx(Du_S,  a,i,j)
                Fxp_x      = -u2 * Dx(Du_Fx,  a,i,j)
                Fyp_x      = -u2 * Dx(Du_Fy,  a,i,j)

                B1p_y      = -u2 * Dy(Du_B1, a,i,j)
                B2p_y      = -u2 * Dy(Du_B2, a,i,j)
                Gp_y       = -u2 * Dy(Du_G,  a,i,j)
                phip_y     = -u2 * Dy(Du_phi,a,i,j)
                Sp_y       = -u2 * Dy(Du_S,  a,i,j)
                Fxp_y      = -u2 * Dy(Du_Fx, a,i,j)
                Fyp_y      = -u2 * Dy(Du_Fy, a,i,j)

                B1_xx      = Dxx(bulk.B1, a,i,j)
                B2_xx      = Dxx(bulk.B2, a,i,j)
                G_xx       = Dxx(bulk.G,  a,i,j)
                phi_xx     = Dxx(bulk.phi,a,i,j)
                S_xx       = Dxx(bulk.S,  a,i,j)

                B1_yy      = Dyy(bulk.B1, a,i,j)
                B2_yy      = Dyy(bulk.B2, a,i,j)
                G_yy       = Dyy(bulk.G,  a,i,j)
                phi_yy     = Dyy(bulk.phi,a,i,j)
                S_yy       = Dyy(bulk.S,  a,i,j)

                B2_xy      = Dx(Dy, bulk.B2, a,i,j)
                G_xy       = Dx(Dy, bulk.G,  a,i,j)
                phi_xy     = Dx(Dy, bulk.phi,a,i,j)
                S_xy       = Dx(Dy, bulk.S,  a,i,j)

                vars = (
                    potential, phi0, u, xi, xi_x, xi_y, xi_xx, xi_yy, xi_xy,
                    B1   , B2   , G   , phi   , S    , Fx    , Fy    , Sd, B1d, B2d, Gd, phid,
                    B1p  , B2p  , Gp  , phip  , Sp   , Fxp   , Fyp   ,
                    B1pp , B2pp , Gpp , phipp , Spp  , Fxpp  , Fypp  ,
                    B1_x , B2_x , G_x , phi_x , S_x  , Fx_x  , Fy_x  ,
                    B1_y , B2_y , G_y , phi_y , S_y  , Fx_y  , Fy_y  ,
                    B1p_x, B2p_x, Gp_x, phip_x, Sp_x , Fxp_x , Fyp_x ,
                    B1p_y, B2p_y, Gp_y, phip_y, Sp_y , Fxp_y , Fyp_y ,
                    B1_xx, B2_xx, G_xx, phi_xx, S_xx ,
                    B1_yy, B2_yy, G_yy, phi_yy, S_yy ,
                           B2_xy, G_xy, phi_xy, S_xy
                )

                A_eq_coeff!(aux.ABCS, vars, sys.gridtype)

                aux.b_vec[a]   = -aux.ABCS[4]
                @inbounds @simd for aa in 1:Nu
                    aux.A_mat[a,aa] = aux.ABCS[1] * Duu[a,aa] + aux.ABCS[2] * Du[a,aa]
                end
                aux.A_mat[a,a] += aux.ABCS[3]
            end

            # BC

            aux.b_vec[1]    = bc.A[i,j]
            aux.A_mat[1,:] .= 0.0
            aux.A_mat[1,1]  = 1.0

            aux.b_vec[end]    = bc.A_u[i,j]
            @inbounds @simd for aa in 1:Nu
                aux.A_mat[end,aa]  = Du[1,aa]
            end

            solve_lin_system!(aux.A_mat, aux.b_vec)

            @inbounds @simd for aa in 1:Nu
                bulk.A[aa,i,j] = aux.b_vec[aa]
            end

        end
    end

    nothing
end

function solve_nested!(bulkconstrain::BulkConstrained, bulkevol::BulkEvolved, bc::BC,
                       gauge::Gauge, deriv::BulkDeriv, aux_acc,
                       sys::System, evoleq::AffineNull)
    Du_B1   = deriv.Du_B1
    Du_B2   = deriv.Du_B2
    Du_G    = deriv.Du_G
    Du_phi  = deriv.Du_phi
    Du_S    = deriv.Du_S
    Du_Fx   = deriv.Du_Fx
    Du_Fy   = deriv.Du_Fy
    Du_A    = deriv.Du_A
    Duu_B1  = deriv.Duu_B1
    Duu_B2  = deriv.Duu_B2
    Duu_G   = deriv.Duu_G
    Duu_phi = deriv.Duu_phi
    Duu_S   = deriv.Duu_S
    Duu_Fx  = deriv.Duu_Fx
    Duu_Fy  = deriv.Duu_Fy
    Duu_A   = deriv.Duu_A

    Du  = sys.Du
    Duu = sys.Duu
    Dx  = sys.Dx
    Dxx = sys.Dxx
    Dy  = sys.Dy
    Dyy = sys.Dyy

    bulk = Bulk(bulkevol, bulkconstrain)

    # solve for S
    vprint("INFO: solve_S")
    solve_S!(bulk, bc, gauge, deriv, aux_acc, sys, evoleq)

    # take u-derivatives of S
    vprint("INFO: S derivatives")
    @sync begin
        @spawn mul!(Du_S,   Du,  bulk.S)
        @spawn mul!(Duu_S,  Duu, bulk.S)
    end

    # solve for Fx and Fy
    vprint("INFO: solve_Fxy")
    solve_Fxy!(bulk, bc, gauge, deriv, aux_acc, sys, evoleq)

    # take u-derivatives of Fx and Fy
    vprint("INFO: Fxy derivatives")
    @sync begin
        @spawn mul!(Du_Fx,   Du,  bulk.Fx)
        @spawn mul!(Du_Fy,   Du,  bulk.Fy)
        @spawn mul!(Duu_Fx,  Duu, bulk.Fx)
        @spawn mul!(Duu_Fy,  Duu, bulk.Fy)
    end

    # solve for Sd
    vprint("INFO: solve_Sd")
    solve_Sd!(bulk, bc, gauge, deriv, aux_acc, sys, evoleq)

    # equations for B2d, (B1d, Gd) and phid are actually independent of each
    # other and could be solved in parallel. however, it seems that use of
    # @spawn here is actually harmful for scaling. since the loops in each
    # function are already threaded, it seems better not to @spawn.

    vprint("INFO: solve_B2d")
    solve_B2d!(bulk, bc, gauge, deriv, aux_acc, sys, evoleq)

    vprint("INFO: solve_B1dGd")
    solve_B1dGd!(bulk, bc, gauge, deriv, aux_acc, sys, evoleq)

    vprint("INFO: solve_phid")
    solve_phid!(bulk, bc, gauge, deriv, aux_acc, sys, evoleq)

    # solve for A
    vprint("INFO: solve_A")
    solve_A!(bulk, bc, gauge, deriv, aux_acc, sys, evoleq)

    # take u-derivatives of A. they will be needed for syncing the domains and
    # also for the xi_t function
    vprint("INFO: A derivatives")
    mul!(Du_A,   Du,  bulk.A)
    mul!(Duu_A,  Duu, bulk.A)

    nothing
end

function syncBCs!(bc::BC, bulk::BulkConstrained, deriv::BulkDeriv)
    Nu, Nx, Ny = size(bulk.S)

    @fastmath @inbounds @threads for j in 1:Ny
        @inbounds @simd for i in 1:Nx
            bc.S[i,j]    = bulk.S[end,i,j]
            bc.Fx[i,j]   = bulk.Fx[end,i,j]
            bc.Fy[i,j]   = bulk.Fy[end,i,j]
            bc.Sd[i,j]   = bulk.Sd[end,i,j]
            bc.B1d[i,j]  = bulk.B1d[end,i,j]
            bc.B2d[i,j]  = bulk.B2d[end,i,j]
            bc.Gd[i,j]   = bulk.Gd[end,i,j]
            bc.phid[i,j] = bulk.phid[end,i,j]
            bc.A[i,j]    = bulk.A[end,i,j]

            bc.S_u[i,j]   = deriv.Du_S[end,i,j]
            bc.Fx_u[i,j]  = deriv.Du_Fx[end,i,j]
            bc.Fy_u[i,j]  = deriv.Du_Fy[end,i,j]
            bc.A_u[i,j]   = deriv.Du_A[end,i,j]
        end
    end

    nothing
end

function set_innerBCs!(bc::BC, bulk::BulkEvolved, boundary::Boundary,
                       gauge::Gauge, deriv::BulkDeriv, sys::System{Inner},
                       evoleq::AffineNull)
    _, Nx, Ny = size(sys)

    Dx  = sys.Dx
    Dy  = sys.Dy

    phi0  = evoleq.phi0
    phi02 = phi0 * phi0
    phi03 = phi0 * phi02
    phi04 = phi02 * phi02

    @fastmath @inbounds @threads for j in 1:Ny
        @inbounds @simd for i in 1:Nx
            xi      = gauge.xi[1,i,j]
            xi3     = xi*xi*xi

            phi     = bulk.phi[1,i,j]
            phi_u   = deriv.Du_phi[1,i,j]

            phi_x   = Dx(bulk.phi, 1,i,j)
            phi_y   = Dy(bulk.phi, 1,i,j)

            xi_x    = Dx(gauge.xi, 1,i,j)
            xi_y    = Dy(gauge.xi, 1,i,j)

            b14     = bulk.B1[1,i,j]
            b24     = bulk.B2[1,i,j]
            g4      = bulk.G[1,i,j]

            a4      = boundary.a4[1,i,j]

            fx2     = boundary.fx2[1,i,j]
            fy2     = boundary.fy2[1,i,j]

            fx2_x   = Dx(boundary.fx2, 1,i,j)
            fy2_y   = Dy(boundary.fy2, 1,i,j)

            b14_x   = Dx(bulk.B1,1,i,j)
            b24_x   = Dx(bulk.B2,1,i,j)
            phi_x   = Dx(bulk.phi,1,i,j)
            g4_x    = Dx(bulk.G,1,i,j)

            b14_y   = Dy(bulk.B1,1,i,j)
            b24_y   = Dy(bulk.B2,1,i,j)
            phi_y   = Dy(bulk.phi,1,i,j)
            g4_y    = Dy(bulk.G,1,i,j)

            phi2    = phi03 * phi - phi0 * xi * xi
            phi2_x  = phi03 * phi_x - 2 * phi0 * xi * xi_x
            phi2_y  = phi03 * phi_y - 2 * phi0 * xi * xi_y

            bc.S[i,j]   = phi04 * (1 - 18 * phi) / 54
            bc.S_u[i,j] = phi02 * (-12 * xi3 +
                                   phi02 * xi * (18 * phi - 5) -
                                   24 * phi02 * phi_u) / 90

            bc.Fx[i,j]   = fx2
            bc.Fx_u[i,j] = -2 * fx2 * xi - 12 / 15 * (b14_x + b24_x - g4_y) +
                4/15 * phi0 * phi2_x

            bc.Fy[i,j]   = fy2
            bc.Fy_u[i,j] = -2 * fy2 * xi - 12 / 15 * (-b14_y + b24_y - g4_x) +
                4/15 * phi0 * phi2_y

            bc.Sd[i,j] = a4 / 2 - 5 * phi04/36 + phi0 * phi2/2

            bc.B2d[i,j] = -2 * b24
            bc.B1d[i,j] = -2 * b14
            bc.Gd[i,j]  = -2 * g4

            bc.A[i,j]   = a4
            bc.A_u[i,j] = -2 * xi * a4 - 2 * phi0 * xi * phi2 -
                2/3 * (phi02 * xi3 + fx2_x + fy2_y + phi04 * phi_u)
        end
    end

    # separate the phid case, since the if statement may prevent loop
    # vectorization

    if abs(phi0) < 1e-9
        fill!(bc.phid, 0)
        return
    end

    @fastmath @inbounds for j in 1:Ny
        @inbounds @simd for i in 1:Nx
            xi      = gauge.xi[1,i,j]
            phi     = bulk.phi[1,i,j]
            phi2    = phi03 * phi - phi0 * xi * xi

            bc.phid[i,j] = 1/3 - 3/2 * phi2 / phi03
        end
    end

    nothing
end

function set_outerBCs!(bc::BC, bulk::BulkConstrained, gauge::Gauge,
                       deriv::BulkDeriv, sys::System, evoleq::AffineNull)
    _, Nx, Ny = size(sys)

    phi0 = evoleq.phi0

    # we are here assuming that the inner and outer grids merely touch at the
    # interface, so we pass the values at this point without any interpolation
    u0 = sys.ucoord[end]

    @fastmath @inbounds @threads for j in 1:Ny
        @inbounds @simd for i in 1:Nx
            xi     = gauge.xi[1,i,j]

            S      = bulk.S[end,i,j]
            Fx     = bulk.Fx[end,i,j]
            Fy     = bulk.Fy[end,i,j]
            Sd     = bulk.Sd[end,i,j]
            B1d    = bulk.B1d[end,i,j]
            B2d    = bulk.B2d[end,i,j]
            Gd     = bulk.Gd[end,i,j]
            phid   = bulk.phid[end,i,j]
            A      = bulk.A[end,i,j]

            S_u    = deriv.Du_S[end,i,j]
            Fx_u   = deriv.Du_Fx[end,i,j]
            Fy_u   = deriv.Du_Fy[end,i,j]
            A_u    = deriv.Du_A[end,i,j]

            bc.S[i,j]   = S_inner_to_outer(S, u0, xi, phi0)
            bc.S_u[i,j] = S_u_inner_to_outer(S_u, S, u0, xi, phi0)

            bc.Fx[i,j]   = F_inner_to_outer(Fx, u0)
            bc.Fy[i,j]   = F_inner_to_outer(Fy, u0)
            bc.Fx_u[i,j] = F_u_inner_to_outer(Fx_u, Fx, u0)
            bc.Fy_u[i,j] = F_u_inner_to_outer(Fy_u, Fy, u0)

            bc.Sd[i,j]  = Sd_inner_to_outer(Sd, u0, xi, phi0)

            # B1d, B2d, and Gd transform in the same way
            bc.B2d[i,j] = Bd_inner_to_outer(B2d, u0)
            bc.B1d[i,j] = Bd_inner_to_outer(B1d, u0)
            bc.Gd[i,j]  = Bd_inner_to_outer(Gd, u0)

            bc.phid[i,j] = phid_inner_to_outer(phid, u0, phi0)

            bc.A[i,j]   = A_inner_to_outer(A, u0, xi, phi0)
            bc.A_u[i,j] = A_u_inner_to_outer(A_u, A, u0, xi, phi0)
        end
    end

    nothing
end

# We assume that the first entry on these arrays is the inner grid, and that
# there is only one domain spanning this grid. If we ever change this
# construction we must remember to make the appropriate changes here.
function solve_nesteds!(bulkconstrains, bulkevols, boundary::Boundary, gauge::Gauge,
                        bcs, derivs, aux_accs,
                        systems::SystemPartition, evoleq::AffineNull)
    Nsys = length(systems)

    # take all u-derivatives of the bulkevols functions
    vprint("INFO: bulkevols derivatives")
    @sync begin
        @inbounds for i in 1:Nsys
            @spawn mul!(derivs[i].Du_B1,  systems[i].Du,  bulkevols[i].B1)
            @spawn mul!(derivs[i].Du_B2,  systems[i].Du,  bulkevols[i].B2)
            @spawn mul!(derivs[i].Du_G,   systems[i].Du,  bulkevols[i].G)
            @spawn mul!(derivs[i].Du_phi, systems[i].Du,  bulkevols[i].phi)
            @spawn mul!(derivs[i].Duu_B1, systems[i].Duu, bulkevols[i].B1)
            @spawn mul!(derivs[i].Duu_B2, systems[i].Duu, bulkevols[i].B2)
            @spawn mul!(derivs[i].Duu_G,  systems[i].Duu, bulkevols[i].G)
            @spawn mul!(derivs[i].Duu_phi,systems[i].Duu, bulkevols[i].phi)
        end
    end

    vprint("INFO: innerBCs")
    set_innerBCs!(bcs[1], bulkevols[1], boundary, gauge, derivs[1], systems[1], evoleq)

    vprint("INFO: solve_nested 1")
    solve_nested!(bulkconstrains[1], bulkevols[1], bcs[1], gauge,
                  derivs[1], aux_accs[1], systems[1], evoleq)

    vprint("INFO: outerBCs")
    set_outerBCs!(bcs[2], bulkconstrains[1], gauge, derivs[1], systems[1], evoleq)

    @inbounds for i in 2:Nsys-1
        vprint("INFO: solve_nested $i")
        solve_nested!(bulkconstrains[i], bulkevols[i], bcs[i], gauge,
                      derivs[i], aux_accs[i], systems[i], evoleq)
        syncBCs!(bcs[i+1], bulkconstrains[i], derivs[i])
    end
    vprint("INFO: solve_nested $Nsys")
    solve_nested!(bulkconstrains[Nsys], bulkevols[Nsys], bcs[Nsys], gauge,
                  derivs[Nsys], aux_accs[Nsys], systems[Nsys], evoleq)

    nothing
end

# for testing only: set all constrained variables to zero
function solve_nesteds!(bulkconstrains, bulkevols, boundary::Boundary, gauge::Gauge,
                        bcs, derivs, aux_accs,
                        systems, evoleq::EvolTest0)
    Nsys = length(systems)

    for i in 1:Nsys
        fill!(bulkconstrains[i], 0)
    end
    nothing
end

struct Nested{S,C,D,B,A}
    systems        :: S
    bulkconstrains :: C
    derivs         :: D
    bcs            :: B
    aux_accs       :: A
end

function Nested(systems::SystemPartition, bulkconstrains::BulkPartition,
                derivs::BulkPartition)
    T        = Jecco.coord_eltype(systems[1].ucoord)
    bcs      = BCs(systems)
    aux_accs = [Aux{T}(sys) for sys in systems]

    Nested{typeof(systems),typeof(bulkconstrains),typeof(derivs),
           typeof(bcs),typeof(aux_accs)}(systems, bulkconstrains, derivs, bcs,
                                         aux_accs)
end

function Nested(systems::SystemPartition, bulkconstrains::BulkPartition)
    T        = Jecco.coord_eltype(systems[1].ucoord)
    bcs      = BCs(systems)
    aux_accs = [Aux{T}(sys) for sys in systems]
    derivs   = BulkDerivPartition(systems)

    Nested(systems, bulkconstrains, derivs)
end

function Nested(systems::SystemPartition)
    T        = Jecco.coord_eltype(systems[1].ucoord)
    bcs      = BCs(systems)
    aux_accs = [Aux{T}(sys) for sys in systems]

    bulkconstrains = BulkConstrainedPartition(systems)
    derivs         = BulkDerivPartition(systems)

    Nested(systems, bulkconstrains, derivs)
end

function (nested::Nested)(bulkevols::BulkPartition, boundary::Boundary,
                          gauge::Gauge, evoleq::EvolutionEquations)
    solve_nesteds!(nested.bulkconstrains, bulkevols, boundary, gauge, nested.bcs,
                   nested.derivs, nested.aux_accs, nested.systems, evoleq)
end


text:
truncated(d::LogUniform, lo::T, hi::T) where {T<:Real} = LogUniform(max(d.a, lo), min(d.b, hi))
truncated(d::LogUniform, lo::Real, ::Nothing) = LogUniform(max(d.a, lo), d.b)
truncated(d::LogUniform, ::Nothing, hi::Real) = LogUniform(d.a, min(d.b, hi))

ext: jl


text:
#------------------------------------------------------------------------------
"""
    OlsGMFn(Y,X)

LS of Y on X; for one dependent variable, Gauss-Markov assumptions

# Usage
(b,u,Yhat,V,R2) = OlsGMFn(Y,X)

# Input
- `Y::Vector`:    T-vector, the dependent variable
- `X::Matrix`:    Txk matrix of regressors (including deterministic ones)

# Output
- `b::Vector`:    k-vector, regression coefficients
- `u::Vector`:    T-vector, residuals Y - yhat
- `Yhat::Vector`: T-vector, fitted values X*b
- `V::Matrix`:    kxk matrix, covariance matrix of b
- `R2::Number`:   scalar, R2 value

"""
function OlsGMFn(Y,X)

    T    = size(Y,1)

    b    = X\Y
    Yhat = X*b
    u    = Y - Yhat

    σ2   = var(u)
    V    = inv(X'X)*σ2
    R2   = 1 - σ2/var(Y)

    return b, u, Yhat, V, R2

end
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
"""
    OlsNWFn(Y,X,m=0)

LS of Y on X; for one dependent variable, using Newey-West covariance matrix

# Usage
(b,u,Yhat,V,R2) = OlsNWFn(Y,X,m)

# Input
- `Y::Array`:     Tx1, the dependent variable
- `X::Array`:     Txk matrix of regressors (including deterministic ones)
- `m::Int`:       scalar, bandwidth in Newey-West

# Output
- `b::Array`:     kx1, regression coefficients
- `u::Array`:     Tx1, residuals Y - Yhat
- `Yhat::Vector`: Tx1, fitted values X*b
- `V::Array`:     kxk matrix, covariance matrix of b
- `R2::Number`:   scalar, R2 value

"""
function OlsNWFn(Y,X,m=0)

    T    = size(Y,1)

    b    = X\Y
    Yhat = X*b
    u    = Y - Yhat

    S    = CovNWFn(X.*u,m)         #Newey-West covariance matrix
    Sxx  = X'X
    V    = inv(Sxx)'S*inv(Sxx)     #covariance matrix of b
    R2   = 1 - var(u)/var(Y)

    return b, u, Yhat, V, R2

end
#------------------------------------------------------------------------------

max_stars_repo_name: duzaichuan/FinancialEconometrics

max_issues_repo_name: duzaichuan/FinancialEconometrics

max_forks_repo_name: duzaichuan/FinancialEconometrics


text:
include("../models/baseline.jl")

model = Baseline()

function run_step(img, speed, command)
    img = convert.(Float32, img[:,:,1:3]) / 255
    img = reshape(img, (size(img)..., 1))
    img = KnetArray{Float32, 4}(img)

    speed = KnetArray{Float32}([speed])

    preds = model(img)
    throttle = preds[1,1]
    steer = preds[2,1]
    return return throttle, steer
end


text:
for n = 1:10
	if eu / factorial(n+1) / 2^n <= 1./10.^5
		@show n
		break
	end
end

max_stars_repo_name: Magikis/Uniwersity

max_issues_repo_name: Magikis/Uniwersity

max_forks_repo_name: Magikis/Uniwersity


text:
__precompile__(true)
module JDF

using Blosc: Blosc
using Missings: Missings
using BufferedStreams
#using RLEVectors
using WeakRefStrings, TimeZones

using StatsBase: rle, inverse_rle, countmap, sample

import Base: show, getindex, setindex!, eltype, names

using Base: size#, @v_str, >=, include, VERSION

using Serialization: serialize, deserialize

if VERSION >= v"1.3.0"
    import Base.Threads: @spawn
else
    macro spawn(_)
        println("JDF: parallel save/load do not work in < Julia 1.3")
    end
end

function __init__()
    if VERSION >= v"1.3.0"
        Blosc.set_num_threads(Threads.nthreads())
    else
        Blosc.set_num_threads(isdefined(Sys, :CPU_CORES) ? Sys.CPU_CORES : Sys.CPU_THREADS)
    end
end

export savejdf, loadjdf, ssavejdf, sloadjdf#, save, load
export column_loader, column_loader!
export type_compress!, type_compress
export compress_then_write
export JDFFile, @jdf_str, jdfmetadata, metadata, size, names
export IsBitsType, eachcol, some_elm, getindex, istable

include("JDFFile.jl")
include("type-writer-loader/Bool.jl")
include("type-writer-loader/Char.jl")
include("type-writer-loader/categorical-arrays.jl")
include("type-writer-loader/pooled-arrays.jl")
include("type-writer-loader/Missing.jl")
include("type-writer-loader/Nothing.jl")
include("type-writer-loader/String.jl")
include("type-writer-loader/StringArray.jl")
include("type-writer-loader/Symbol.jl")
include("type-writer-loader/ZonedDateTime.jl")
include("type-writer-loader/substring.jl")

include("column_loader.jl")
include("compress_then_write.jl")

include("load-columns.jl")
include("loadjdf.jl")
include("savejdf.jl")
include("type_compress.jl")

include("metadata.jl")
include("eachcol.jl")
include("Tables.jl")

# Blosc.set_num_threads(6)

end # module


text:
struct Path
    nodes::Vector{Node}
    edges::Vector{Edge}
end

function get_node(p::Path, idx::Int)
    return p.nodes[idx]
end

function get_relationship(p::Path, idx::Int)
    return p.edges[idx]
end

function first_node(p::Path)
    return p.nodes[1]
end

function last_node(p::Path)
    return p.nodes[nodes_count(p)]
end

function edge_count(p::Path)
    return length(p.edges)
end

function nodes_count(p::Path)
    return length(p.nodes)
end


text:
### The functions in this file extend functions of the same names in `Base` and `LinearAlgebra`.
### It may be necessary to instead replace the `Base` and `LinearAlgebra` functions
### if `IsApprox` is widely adopted.

### isone, iszero

isone(x, approx_test::AbstractApprox=Equal()) = isapprox(approx_test, x, one(x))
iszero(x, approx_test::AbstractApprox=Equal()) = isapprox(approx_test, x, zero(x))
iszero(x::AbstractArray, approx_test::AbstractApprox=Equal()) =
    all(y -> iszero(y, approx_test), x)
isone(x::BigInt, ::Equal) = Base.isone(x)
iszero(x::BigInt, ::Equal) = Base.iszero(x)

isone(A::StridedMatrix, approx_test::Approx) = isapprox(approx_test, A, one(A))
iszero(A::StridedMatrix, approx_test::Approx) = isapprox(approx_test, A, zero(A))

# dense.jl
const ISONE_CUTOFF = 2^21 # 2M

function isone(A::StridedMatrix, approx_test::AbstractApprox=Equal())
    m, n = size(A)
    m != n && return false # only square matrices can satisfy x == one(x)
    if sizeof(A) < ISONE_CUTOFF
        _isone_triacheck(A, m, approx_test)
    else
        _isone_cachefriendly(A, m, approx_test)
    end
end

@inline function _isone_triacheck(A::StridedMatrix, m::Int, approx_test::AbstractApprox)
    @inbounds for i in 1:m, j in i:m
        if i == j
            isone(A[i,i], approx_test) || return false
        else
            iszero(A[i,j], approx_test) && iszero(A[j,i], approx_test) || return false
        end
    end
    return true
end

# Inner loop over rows to be friendly to the CPU cache
@inline function _isone_cachefriendly(A::StridedMatrix, m::Int, approx_test::AbstractApprox)
    @inbounds for i in 1:m, j in 1:m
        if i == j
            isone(A[i,i], approx_test) || return false
        else
            iszero(A[j,i], approx_test) || return false
        end
    end
    return true
end

### ishermitian, issymmetric

# The call `ishermitian(A::AbstractMatrix, B::AbstractMatrix) lowers
# to exactly the same code as that in LinearAlgebra
function ishermitian(A::AbstractMatrix, approx_test::AbstractApprox=Equal())
    indsm, indsn = axes(A)
    if indsm != indsn
        return false
    end
    for i = indsn, j = i:last(indsn)
        if ! isapprox(approx_test, A[i,j], adjoint(A[j,i]))
            return false
        end
    end
    return true
end

# This uses the isapprox interface, which compares using a norm
ishermitian(A::AbstractMatrix, approx_test::Approx) = isapprox(approx_test, A, adjoint(A))

ishermitian(x::Number, approx_test::AbstractApprox=Equal()) = isapprox(approx_test, x, conj(x))

issymmetric(A::AbstractMatrix{<:Real}, approx_test::AbstractApprox=Equal()) =
    ishermitian(A, approx_test)

# Copied from LinearAlgebra. Why does the iteration over i differ slightly from ishermitian ?
function issymmetric(A::AbstractMatrix, approx_test::AbstractApprox=Equal())
    indsm, indsn = axes(A)
    if indsm != indsn
        return false
    end
    for i = first(indsn):last(indsn), j = (i):last(indsn)
        if ! isapprox(approx_test, A[i,j], transpose(A[j,i]))
            return false
        end
    end
    return true
end

issymmetric(x::Number, approx_test::AbstractApprox=Equal()) = isapprox(approx_test, x, x)

### isreal

# complex.jl
isreal(x::Real, approx_test::AbstractApprox=Equal()) = true
isreal(z::Complex, approx_test::AbstractApprox=Equal()) = isapprox(approx_test, real(z), z)
# Old way
#isreal(z::Complex, approx_test::AbstractApprox=Equal()) = iszero(imag(z), approx_test)

### isinteger

# This must be changed if this is integrated into Base
isinteger(x) = isinteger(x, Equal())
# We use union or explicit types to avoid method ambiguity. Is there a way around this ?
isinteger(x::BigFloat, ::Equal) = Base.isinteger(x)
isinteger(x::Rational, ::Equal) = Base.isinteger(x)

# number.jl
isinteger(x::Integer, ::AbstractApprox=Equal()) = true

# floatfuncs.jl
## The original is x - trunc(x) == 0. So this implementation might differ
## from the stock version of isinteger for some user's subtype of AbstractFloat.
## We choose to this implementation because the default relative tolerance is
## reasonable. That is, `isapprox(approx_test, x - trunc(x), 0)` requires
## specifying `atol`.
isinteger(x::AbstractFloat, approx_test::AbstractApprox=Equal()) = isapprox(approx_test, x, trunc(x))

# complex.jl
isinteger(z::Complex, approx_test::AbstractApprox=Equal()) =
    isreal(z, approx_test) & isinteger(real(z), approx_test)

isinteger(x::Rational, approx_test::AbstractApprox) = isinteger(float(x), approx_test)

### istriu, istril, isbanded, isdiag

function istriu(A::AbstractMatrix, k::Integer = 0; approx::AbstractApprox=Equal())
    Base.require_one_based_indexing(A)
    m, n = size(A)
    for j in 1:min(n, m + k - 1)
        for i in max(1, j - k + 1):m
            iszero(A[i, j], approx) || return false
        end
    end
    return true
end
istriu(x::Number, ::AbstractApprox) = true

function istril(A::AbstractMatrix, k::Integer = 0; approx::AbstractApprox=Equal())
    Base.require_one_based_indexing(A)
    m, n = size(A)
    for j in max(1, k + 2):n
        for i in 1:min(j - k - 1, m)
            iszero(A[i, j], approx) || return false
        end
    end
    return true
end
istril(x::Number, ::AbstractApprox) = true

isbanded(A::AbstractMatrix, kl::Integer, ku::Integer; approx::AbstractApprox=Equal()) =
    istriu(A, kl; approx=approx) && istril(A, ku; approx=approx)

isdiag(A::AbstractMatrix; approx::AbstractApprox=Equal()) = isbanded(A, 0, 0; approx=approx)
isdiag(x::Number; approx::AbstractApprox=Equal()) = true


text:
module Fourier1D
using ..FourierBasis: inverse_unidimensional_index, ϕ, purge

using LinearAlgebra, SparseArrays, IntervalArithmetic
using FastTransforms, FFTW

export assemble_matrix, noise_matrix, truncation_error, rigorous_norm
export evaluate_trig_polynomial, abscissas

function restrictfft!(new::Vector, orig::Vector, Nx)
    FFTNx = length(orig)
    new[1:Nx+1] = @view orig[1:Nx+1]
    new[Nx+1:2*Nx+1] = @view orig[FFTNx-Nx:FFTNx]
    return new
end

# function arb_restrictfft(x,y, Nx)
#     FFTNx = length(y)
#     for i in 1:Nx+1
#         x[i] =  y[i]
#         #x[Nx+1:2*Nx+1] = y[FFTNx-Nx:FFTNx]
#         x[Nx+i] = y[FFTNx-Nx+i-1]
#     end
#     return x
# end

function extendfft(orig::Vector, FFTNx)
    N = length(orig)
    Nx = (N-1) ÷ 2
    return [orig[1:Nx]; zeros(Complex{Float64}, FFTNx-N); orig[Nx+1:end]]
end

"""
assemble_matrix(T, Nx; FFTNx, x_0, x_1, ϵ)

This function assembles the matrix of the transfer operator of T using the FFT

Arguments:
- T : the dynamic
- Nx : the truncation frequence
- FFTNx : the number of points used to compute the FFT
- x_0 : left endpoint of the domain of T, defaults to 0
- x_1 : right endpoint of the domain of T, defaults to 1
- ϵ : the purge threshold, under this threshold the entries are set to 0
"""
#import IntervalArithmetic: radius
#radius(x::Complex{Interval{T}}) where T = radius(real(x))

using FastRounding

function  norm_2_upper_bound(v::Vector{Float64})
    sum = 0.0
    for x in v
        sum = sum ⊕₊ square_round(x, RoundUp)
    end
    return sqrt_round(sum, RoundUp)
end

function  norm_2_upper_bound(v::Vector{BigFloat})
    old_rnd = rounding(BigFloat)
    setrounding(BigFloat, RoundUp)
    sum = 0.0
    for x in v
        sum = sum + x^2
    end
    sum = sqrt(sum)
    setrounding(BigFloat, old_rnd)
    return sum
end

function abs_upper_square(x::Complex{Float64})
    absval = square_round(real(x), RoundUp)
    absval = absval ⊕₊ square_round(imag(x), RoundUp)
    return absval
end

function  norm_2_upper_bound(v::Vector{Complex{Float64}})
    sum = 0.0
    for x in v
        sum = sum ⊕₊ abs_upper_square(x)
    end
    return sqrt_round(sum, RoundUp)
end

function  norm_2_upper_bound(v::Vector{Complex{BigFloat}})
    old_rnd = rounding(BigFloat)
    setrounding(BigFloat, RoundUp)
    sum = 0.0
    for x in v
        sum = sum + real(x)^2+imag(x)^2
    end
    sum = sqrt(sum)
    setrounding(BigFloat, old_rnd)
    return sum
end

"""
Assembles the truncated Fourier matrix for a dynamic D,
with frequencies -Nx≤ i≤ Nx

Input:
    D: Dynamic
    Nx: size of the basis
    FFTNx: size of the FFT
    x_0 : left interval of the dynamic interval
    x_1 : right interval of the dynamic interval
    T: Floating point type, supports Float64, BigFloat

    Output:
    M: matrix
    n: ||M-R||_2 where R is the interval of matrix that contains M
"""
function assemble_matrix(D, Nx; FFTNx = 2048, x_0 = 0, x_1 = 1, T = Float64)
    I = Interval{T}
    dx = [I(x_0)+I(i)*(I(x_1)-I(x_0))/FFTNx for i in 0:FFTNx-1];

    Lx = abs(I(x_1)-I(x_0))

    one = [T(1) for x in dx]
    P = plan_fft(one)
    # we take the adjoint since we are computing the adjoint operator

    Dx = [D(x) for x in dx]

    N = (2*Nx+1) # we are taking Nx positive and negative frequencies and the 0 frequency

    M = zeros(Complex{T},(2*Nx+1, 2*Nx+1))

    observablevalue = zeros(Complex{T}, FFTNx)
    observablerad = zeros(Float64, FFTNx)

    onedtransform = zeros(Complex{T}, FFTNx)
    new = zeros(Complex{T}, 2*Nx+1)
    l2error = 0.0
    # a priori estimate of the FFT error
    # we use the fact that
    # ||FFT(v)/FFTNx||₂ ≤ 1/√FFTNX
    # and
    # ||Fl(FFT(v))- FFT(v)||₂ ≤ tη/(1-η)||v||₂
    # where η = μ + γ₄(√2+μ), t > log₂(FFTNx)
    # and μ is the absolute error in the computation of the wiggle factors
    # (we assume them to be precise to machine precision)
    # and γ₄ = 4u/(1-4u) and u is the unit roundoff
    # this bound is from Higham N. J. - Accuracy and Stability of Numerical Algorithms
    # Second Edition - SIAM

    u = Float64(eps(T), RoundUp)
    γ₄ = (4.0 ⊗₊ u)⊘₊(1.0 ⊖₋ 4.0 ⊗₋ u)
    #@info γ₄
    μ = u
    η = μ ⊕₊ γ₄ ⊗₊ (sqrt_round(2.0, RoundUp) ⊕₊ μ)
    #@info μ
    t = ceil(log2(FFTNx))
    rel_err_fft = (t ⊗₊ η) ⊘₊(1.0 ⊖₋ η)
    f_FFTNx = Float64(FFTNx)
    norm_FFT_normalized_2 = 1.0 ⊘₊(sqrt(f_FFTNx, RoundUp))

    for i in 1:2*Nx+1
        l = inverse_unidimensional_index(i, Nx) # the index in the form [0, ..., Nx, -Nx, ..., -1]

        for (ind, val) in pairs(Dx)
            obsval = ϕ(l, val; L=Lx)
            real_m, real_r = midpoint_radius(real(obsval))
            imag_m, imag_r = midpoint_radius(imag(obsval))
            observablevalue[ind] =  real_m+im*imag_m
            observablerad[ind] = sqrt(real_r^2+imag_r^2)
        end
        norm_obs = Float64(norm_2_upper_bound(observablevalue), RoundUp)
        norm_rad = Float64(norm_2_upper_bound(observablerad), RoundUp)
        err_fft = norm_FFT_normalized_2⊗₊(rel_err_fft ⊗₊ norm_obs)⊕₊ norm_FFT_normalized_2⊗₊norm_rad
        l2error= max(l2error, err_fft)

        mul!(onedtransform, P, observablevalue)
        restrictfft!(new, onedtransform, Nx)

        for (ind, val) in pairs(new)
                if abs(val)!= 0
                    # this is the adjoint matrix of the Koopman operator
                    M[i, ind] = conj(val)/FFTNx
                end
        end
    end
    return M, l2error ⊗₊ sqrt_round(Float64(N), RoundUp)
end

# using ArbNumerics

# """
# assemble_matrix_high_precision(T, Nx; FFTNx, x_0, x_1, ϵ, prec)

# This function assembles the matrix of the transfer operator of T using the FFT in high precision, using ArbNumerics

# Arguments:
# - T : the dynamic
# - Nx : the truncation frequence
# - FFTNx : the number of points used to compute the FFT
# - x_0 : left endpoint of the domain of T, defaults to 0
# - x_1 : right endpoint of the domain of T, defaults to 1
# - ϵ : the purge threshold, under this threshold the entries are set to 0
# - prec : the internal working precision (which is the same as the displayed precision with setextrabits(0)) to a given number of bits, defaults to 1000
# """
# function assemble_matrix_high_precision(T, Nx; FFTNx = 2048, x_0 = 0, x_1 = 1, ϵ = 2^-30, prec = 300)
#     setextrabits(0)
#     setprecision(ArbReal, bits = prec)
#     setprecision(ArbComplex, bits = prec)

#     dx = [ArbReal(x_0+i*(x_1-x_0)/FFTNx) for i in 0:FFTNx-1];

#     Lx = ArbReal(abs(x_1-x_0))
#     Tx = [ArbComplex((T(x))) for x in dx]

#     N = (2*Nx+1) # we are taking Nx positive and negative frequencies and the 0 frequency

#     M = zeros(ArbComplex,(N, N))
#     for i in 1:N
#         l = inverse_unidimensional_index(i, Nx) # the index in the form [0, ..., Nx, -Nx, ..., -1]
#         observablevalue = ϕ.(l, Tx; L=Lx)
#         @info norm(observablevalue, 2)
#         onedtransform = dft(observablevalue)

#         for j in 1:Nx
#             y = conj(onedtransform[j])/FFTNx
#             @info y
#             #M[i, j] = y #conj()/FFTNx
#         end

#         #new = [ArbComplex(0) for s in 1:2*Nx+1]
#         #new = arb_restrictfft(new, onedtransform, Nx)
#         #@info new

#         #for (ind, val) in pairs(new)
#         #        if abs(val)!= 0
#                     # this is the adjoint matrix of the Koopman operator
#         #            M[i, ind] = conj(val)/FFTNx
#         #        end
#         #end
#     end
#     return M
# end
"""
noise_matrix(σ, Nx)

returns the matrix of the (periodic) convolution operator with a Gaussian of average 0
and variance σ in the Fourier basis, truncated at frequence Nx
(Check what happens when we change x_0 and x_1, probably rescaling σ)

"""
noise_matrix(σ, Nx) = Diagonal([[exp(-(k*σ)^2/2) for k in 0:Nx]; [exp(-(k*σ)^2/2) for k in -Nx:-1]])

"""
truncation_error(σ, Nx)

This function returns the truncation error when truncating frequencies at Nx
when the Gaussian has variance σ
(check the rescaling)
"""
truncation_error(σ, Nx) = exp(-(Nx*σ)^2/2)/(Nx*σ*sqrt(2*pi))

"""
This function computes a candidate for the top singular value
of M|_{V_0} square, i.e., it returns λ^2 and the associated singular vector
"""
function norm_2_estimator_square_nonrig(M; n_it= 10)
    n, m = size(M)
    @assert n==m
    v = ones(n-1)
    N = M[2:end, 2:end]

    start_norm =  norm(v, 2)

    for i in 1:n_it
        v = N'*(N*v)
        v = v/norm(v,2)
    end
    w = v

    v = N*v

    return abs(v'*v), w/norm(w,2)
end

using IntervalArithmetic, IntervalRootFinding

widen(x::T, fatten) where {T} = x+2.0^(-fatten)*Interval{T}(-1, 1)
widen(x::Complex{T}, fatten) where {T} = x+2.0^(-fatten)*(Interval{T}(-1, 1)+im*Interval{T}(-1, 1))

function Yamamoto_certify(M, λ, v; fatten = 10)
    N = M[2:end, 2:end]
    F(λ, v) = [N'*N*v-λ*v; v'*v-1]
    v_fat = widen.(v, fatten)
    lam = widen(λ, fatten)

    DF(λ, v) = [N'*N-λ*I -v;
                v'          0]

    w = [lam; v_fat]
    w0 = [λ; v] - DF(lam, v_fat)\F(λ, v)

    w = [λ; v] - DF(w0[1], w0[2:end])\F(λ, v)

    #t = intersect( abs(w0[1]), abs(w[1]))

    return (abs(w[1])).hi
end

"""
This function uses the classical bound, cited in

RESIDUAL BOUNDS ON APPROXIMATE EIGENSYSTEMS OF NONNORMAL MATRICES*
W. KAHAN, B. N. PARLETT AND E. JIANGt

Theorem 1

"""
function faster_certify(A, λ, v)
    N = A[2:end, 2:end]
    w = Interval.(v)
    ν = Interval(λ)
    #@info λ
    ρ = norm(N'*N*w-ν*w, 2)/norm(w, 2)
    ν = hull(-ρ, ρ) + λ
    return ν
end

using IntervalArithmetic

upper(x::Interval) = x.hi

"""
rigorous_norm(M; k = 10)

Computes the L² and L∞ norm of the powers of M|_V₀ up to the power k.
The algorithm first computes a numeric estimate of the top Singular Value
and then certifies it by an Interval Newton step

Outputs: (v2, v∞)
- v₂ : it is the vector with the rigorous bounds on the L² norm
- v∞ : it is the vector with the rigorous bounds on the L∞ norm;
       these are obtained by observing that
       ||f||₂ = ||w||₂,
       where w is the vector of the Fourier coefficients of the trigonometric polynomial f
       ||f||₂ <= ||f||∞ <= ∑ |w_i| <= √n*||f||₂
       so
       ||M||∞ <= √n *||M||₂
"""
function rigorous_norm(M; k = 1000)
   A = M
   n, m = size(M)
   norms = zeros(Interval{Float64}, k)
   j = 0
   for i in 1:k
       λ, v = norm_2_estimator_square_nonrig(A)
       #@info λ
       #norms[i] = Yamamoto_certify(A, λ, v; fatten = 100)
       norms[i] = sqrt(faster_certify(A, λ, v))
       A = M*A
       j = j+1
       @info(j)
   end
   rescale = sqrt(Interval(n-1))
   @info "rescale", rescale
   return norms, rescale*norms #check to make it rigorous
end

"""
evaluate_trig_polynomial(v; FFTNx = 16384)

Extends a vector v to size FFTNx (padding with zeros), and then takes the real part of the
inverse fourier transform, this gives us the value of IFFT(v) at FFTNx points.
Useful for plots... beware the fact that when plotting you need to give abscissas
"""
evaluate_trig_polynomial(v; FFTNx = 16384) = real.(ifft(extendfft(v, FFTNx)))

"""
Returns the abscissas for plotting, i.e., FFTNx equispaced points between x_0 and x_1
"""
abscissas(x_0, x_1, FFTNx) = [x_0+(i-1)*(x_1-x_0)/FFTNx for i in 1:FFTNx]

end


text:
using Test, LinearAlgebra
using JuliaBLAS: Block, addmul!

function matmul_test(name, datagen, generic_only = false)
    @testset "$name" begin
        siz = (3, 5, 9, 13, 32, 200)
        for m in siz, n in siz, k in siz
            C = datagen(m, n)
            A = datagen(m, k)
            B = datagen(k, n)
            @test true # ensure base functionality doesn't error out before here

            if !generic_only
                ABC = A*B + C
                addmul!(C, A, B, generic=Val(false), mr=Val(12), nr=Val(4))
                @test C ≈ ABC
                ABC = A*B + C
                addmul!(C, A, B, generic=Val(false), mr=Val(8), nr=Val(6))
                @test C ≈ ABC
                ABC = A*B + C
                addmul!(C, A, B, generic=Val(false), mr=Val(4), nr=Val(4))
                @test C ≈ ABC
            end

            ABC = A*B + C
            addmul!(C, A, B, generic=Val(true), mr=Val(12), nr=Val(4))
            @test C ≈ ABC
            ABC = A*B + C
            addmul!(C, A, B, generic=Val(true), mr=Val(8), nr=Val(6))
            @test C ≈ ABC
            ABC = A*B + C
            addmul!(C, A, B, generic=Val(true), mr=Val(4), nr=Val(4))
            @test C ≈ ABC
        end
    end
end

matmul_test.([
              Int64,
              Int32,
              UInt64,
              UInt32,
              Float64,
              Float32,
             ],
              randn)

# Weird type, with weird math
#struct Stringish
#    str::String
#end
#
#Stringish(x) = Stringish(string(x))
#
#Base.zero(::Stringish) = Stringish("")
#Base.one(::Stringish) = Stringish("")
#Base.:+(s1::Stringish, s2::Stringish) = Stringish(s1.str * s2.str)
#Base.:*(s1::Stringish, s2::Stringish) = Stringish("[$(s1.str),$(s2.str)]")
#
#Base.rand(::Type{Stringish}, (m, n)::Tuple{Int64, Int64}) = [Stringish(rand(0:9)) for ii in 1:m, jj in 1:n]
#Base.:≈(x::Stringish, y::Stringish) = x.str == y.str
##################################################################
#
#
#datagen(T) = (m,n) -> rand(T, (m,n))
#
#@testset "Matrix Multiplication Tests" begin
#    for T in [
#              #Int128,
#              #Int16,
#              #Int32,
#              Int64,
#              #Int8,
#              #UInt128,
#              #UInt16,
#              #UInt32,
#              UInt64,
#              #UInt8,
#              #Float16,
#              #Float32,
#              Float64,
#              #Bool
#             ]
#
#        matmul_test(string(T), datagen(T))
#        #T = Complex{T}
#        #matmul_test(string(T), datagen(T))
#    end
#
#
#    #for T in [BigFloat,]
#    #    matmul_test(string(T), datagen(T), true)
#    #    T = Complex{T}
#    #    matmul_test(string(T), datagen(T), true)
#    #end
#
#    matmul_test(string(Stringish), datagen(Stringish), true)
#
#    matmul_test(string(BigInt), (m,n) -> rand(-big"2"^256 : big"2"^256, (m,n)), true)
#end


text:
"""
    IndicatorSOS1Bridge{T,S<:MOI.AbstractScalarSet}

The `IndicatorSOS1Bridge` replaces an indicator constraint of the following
form:
``z \\in \\mathbb{B}, z == 1 \\implies f(x) \\in S`` with a SOS1 constraint:
``z \\in \\mathbb{B}, slack \\text{ free}, f(x) + slack \\in S, SOS1(slack, z)``.
"""
struct IndicatorSOS1Bridge{T,S<:MOI.AbstractScalarSet} <: AbstractBridge
    slack::MOI.VariableIndex
    z::MOI.VariableIndex
    sos_index::MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.SOS1{T}}
    affine_index::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},S}
end

function bridge_constraint(
    ::Type{IndicatorSOS1Bridge{T,S}},
    model::MOI.ModelLike,
    f::MOI.VectorAffineFunction{T},
    s::MOI.Indicator{MOI.ACTIVATE_ON_ONE,S},
) where {T<:Real,S}
    f_scalars = MOI.Utilities.eachscalar(f)
    z = convert(MOI.VariableIndex, f_scalars[1])
    slack = MOI.add_variable(model)
    sos_index = MOI.add_constraint(
        model,
        MOI.VectorOfVariables([slack, z]),
        MOI.SOS1{T}([0.4, 0.6]),  # This weight vector is arbitrary!
    )
    new_f = MOI.Utilities.operate(+, T, f_scalars[2], slack)
    affine_index = MOI.add_constraint(model, new_f, s.set)
    return IndicatorSOS1Bridge{T,S}(slack, z, sos_index, affine_index)
end

function MOI.supports_constraint(
    ::Type{<:IndicatorSOS1Bridge},
    ::Type{<:MOI.AbstractVectorFunction},
    ::Type{<:MOI.Indicator{MOI.ACTIVATE_ON_ONE,<:MOI.AbstractScalarSet}},
)
    return true
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintSet,
    b::IndicatorSOS1Bridge,
)
    set = MOI.get(model, attr, b.affine_index)
    return MOI.Indicator{MOI.ACTIVATE_ON_ONE}(set)
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintFunction,
    b::IndicatorSOS1Bridge{T},
) where {T}
    f = MOI.get(model, attr, b.affine_index)
    terms = MOI.VectorAffineTerm{T}[
        MOI.VectorAffineTerm(2, t) for t in f.terms if t.variable != b.slack
    ]
    push!(terms, MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(one(T), b.z)))
    return MOI.VectorAffineFunction(terms, [zero(T), f.constant])
end

function MOI.delete(model::MOI.ModelLike, bridge::IndicatorSOS1Bridge)
    MOI.delete(model, bridge.sos_index)
    MOI.delete(model, bridge.affine_index)
    MOI.delete(model, bridge.slack)
    return
end

function MOI.Bridges.added_constrained_variable_types(
    ::Type{<:IndicatorSOS1Bridge},
)
    return Tuple{Type}[]
end

function MOI.Bridges.added_constraint_types(
    ::Type{<:IndicatorSOS1Bridge{T,S}},
) where {T,S}
    return Tuple{Type,Type}[
        (MOI.VectorOfVariables, MOI.SOS1{T}),
        (MOI.ScalarAffineFunction{T}, S),
    ]
end

function concrete_bridge_type(
    ::Type{<:IndicatorSOS1Bridge{T}},
    ::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.Indicator{MOI.ACTIVATE_ON_ONE,S}},
) where {T,S}
    return IndicatorSOS1Bridge{T,S}
end

MOI.get(::IndicatorSOS1Bridge, ::MOI.NumberOfVariables)::Int64 = 1

function MOI.get(b::IndicatorSOS1Bridge, ::MOI.ListOfVariableIndices)
    return [b.slack]
end

function MOI.get(
    ::IndicatorSOS1Bridge,
    ::MOI.NumberOfConstraints{MOI.VectorOfVariables,<:MOI.SOS1},
)::Int64
    return 1
end

function MOI.get(
    ::IndicatorSOS1Bridge{T,S},
    ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T},S},
)::Int64 where {T,S}
    return 1
end

function MOI.get(
    b::IndicatorSOS1Bridge{T},
    ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables,<:MOI.SOS1},
) where {T}
    return [b.sos_index]
end

function MOI.get(
    b::IndicatorSOS1Bridge{T,S},
    ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T},S},
) where {T,S}
    return [b.affine_index]
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimal,
    bridge::IndicatorSOS1Bridge,
)
    z = MOI.get(model, MOI.VariablePrimal(attr.result_index), bridge.z)
    w = MOI.get(model, MOI.VariablePrimal(attr.result_index), bridge.slack)
    f = MOI.get(model, attr, bridge.affine_index)
    return [z, f - w]
end

function MOI.supports(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimalStart,
    ::Type{IndicatorSOS1Bridge{T,S}},
) where {T,S}
    ci = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},S}
    return MOI.supports(model, MOI.VariablePrimalStart(), MOI.VariableIndex) &&
           MOI.supports(model, attr, ci)
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimalStart,
    bridge::IndicatorSOS1Bridge,
)
    z = MOI.get(model, MOI.VariablePrimalStart(), bridge.z)
    w = MOI.get(model, MOI.VariablePrimalStart(), bridge.slack)
    f = MOI.get(model, attr, bridge.affine_index)
    return [z, f - w]
end

function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimalStart,
    bridge::IndicatorSOS1Bridge{T},
    value::AbstractVector,
) where {T}
    @assert length(value) == 2
    MOI.set(model, MOI.VariablePrimalStart(), bridge.z, value[1])
    w = something(
        MOI.get(model, MOI.VariablePrimalStart(), bridge.slack),
        zero(T),
    )
    MOI.set(model, attr, bridge.affine_index, value[2] + w)
    return
end


text:
module PanelCuttingReadmeExamples

using Markdown
using PanelCutting
using Unitful
using UnitfulUS
using Test

md"""
Suppose we want to make an open-topped box whose outside dimensions
are 12 inches long, 6 inches wide and 4 inches deep out of 1/4 inch
plywood.
"""

box_length = 12u"inch"
box_width = 6u"inch"
box_depth = 4u"inch"

stock_thickness = (1/4)u"inch"

md"""

To use `PanelCutting` we must first describe the target
dimensions of sheet stock after it is cut.  For each panel we must
make a `WantedPanel`.

Note that these dimensions don't account for any joinery.  For this
simple example all panels are butt jointed, with the sides and ends
sitting on top of the bottom and the ends inset bwetween the sides.

Also note that we do not yet include the thickness of stock when
selecting material to cut.  so far my projects have been simple enough
that I've only needed one thickness of stock.  I the furure we should
also model stock thickness and material type.
"""

wanted = [
    WantedPanel(
        label = "bottom",
        # the bottom will occupy the full footprint of the box with
        # the sides and ends sitting on top of it
        length = box_length,
        width = box_width),
    # Sometimes we need several instances of the same shaped panel.
    # WantedPanel can take a pre-multiplier:
    (2 * WantedPanel(
        label = "side",
        length = box_length - 2 * stock_thickness,
        width = box_depth - stock_thickness)
     )...,  # Pre-multiplying gives a Vector, which we must spread.
    (2 * WantedPanel(
        label = "end",
        length = box_width - 2 * stock_thickness,
        width = box_depth - stock_thickness)
     )...
         ]

md"""
We must also describe what sheet stock is available to cut from.  We
do this with `AvasilablePanel`.  Each size of `AvailablePanel` only
needs to be instantiated once.  It is assumed that the supply of
`AvailablePanel`s is unlimited for our purposes.  These sizes and
prices are taken from my local plywood supplier.  My supplier will
accurately cut down stock to specified dimensions at a reasonable
price, so we also note the cost per cut and the kerf width.
"""

BOULTER_PLYWOOD = Supplier(
    name = "Boulter Plywood",
    kerf = (1/8)u"inch",
    cost_per_cut = money(1.50),
    available_stock = [
        AvailablePanel(
            label = "1/4 5 × 5 Baltic Birch",
            length = 5u"ft",
            width = 5u"ft",
            cost=69.00   # I've not been able to get UnitfulCurrency to work
        ),
        AvailablePanel(
            label = "1/4 30 × 60 Baltic Birch",
            length = 60u"inch",
            width = 30u"inch",
            cost=42.00
        ),
        AvailablePanel(
            label = "1/4 30 × 30 Baltic Birch",
            length = 30u"inch",
            width = 30u"inch",
            cost=26.00)
      ])

md"""
Now that we have described what we want and what we're starting width,
we can search for an optimal cut pattern.
"""

@testset "README Example 1" begin
    # The cut optimization is performed by a Searcher:
    searcher = Searcher(BOULTER_PLYWOOD, wanted)
    search(searcher)
    # The most efficient cut pattern (least expensive in terms of
    # stock used and number of cuts) is in `searcher.cheapest`.
    @test length(searcher.cheapest.finished) == 5

    open("example_panel_cut_report.html", "w") do io
        write(io, report(searcher;
                         includeCutDiagram=true,
                         includeCutGraph=false).content)
    end
end

end


text:
export SatNeighborhood

struct SatNeighborhood <: NeighborhoodAlgorithm
    k::Int
end

function SatNeighborhood()
    return SatNeighborhood(32)
end

function optimize_neighborhood!(algo::SatNeighborhood, index::SearchGraph{T}, dist::Function, perf, recall) where {T}
end

function neighborhood(algo::SatNeighborhood, index::SearchGraph{T}, dist::Function, item::T) where {T}
    N = Int32[]
    knn = search(index, dist, item, KnnResult(algo.k))
    @inbounds for p in knn
        pobj = index.db[p.objID]
        near = NnResult()
        push!(near, zero(Int32), p.dist)
        for nearID in N
            d = dist(index.db[nearID], pobj)
            push!(near, nearID, d)
        end

        if first(near).objID == 0
            push!(N, p.objID)
        end
    end

    return knn, N
end


text:
let
    function rosenbrock{T}(x::Vector{T})
        o = one(T)
        c = convert(T,100)
        return (o - x[1])^2 + c * (x[2] - x[1]^2)^2
    end

    function rosenbrock_gradient!{T}(x::Vector{T}, storage::Vector{T})
        o = one(T)
        c = convert(T,100)
        storage[1] = (-2*o) * (o - x[1]) - (4*c) * (x[2] - x[1]^2) * x[1]
        storage[2] = (2*c) * (x[2] - x[1]^2)
    end

    function rosenbrock_hessian!{T}(x::Vector{T}, storage::Matrix{T})
        o = one(T)
        c = convert(T,100)
        f = 4*c
        storage[1, 1] = (2*o) - f * x[2] + 3 * f * x[1]^2
        storage[1, 2] = -f * x[1]
        storage[2, 1] = -f * x[1]
        storage[2, 2] = 2*c
    end

    d2 = DifferentiableFunction(rosenbrock,
                                rosenbrock_gradient!)
    d3 = TwiceDifferentiableFunction(rosenbrock,
                                     rosenbrock_gradient!,
                                     rosenbrock_hessian!)

    for method in (NelderMead(), SimulatedAnnealing())
        optimize(rosenbrock, [0.0,0,.0], method)
        optimize(rosenbrock, Float32[0.0, 0.0], method)
    end

    for method in (BFGS(),
                   ConjugateGradient(),
                   GradientDescent(),
                   MomentumGradientDescent(),
                   AcceleratedGradientDescent(),
                   LBFGS())
        optimize(d2, [0.0,0,.0], method)
        optimize(d2, Float32[0.0, 0.0], method)
    end

    for method in (Newton(),)# NewtonTrustRegion())
        optimize(d3, [0.0,0.0], method)
        optimize(d3, Float32[0.0, 0.0], method)
    end
end


text:
module field

using QuantumOpticsBase, LinearAlgebra
using ..interaction, ..meanfield, ..mpc, ..quantum, ..reducedspin, ..CollectiveSpins

"""
    field.intensity(r, S, state)

    Arguments:
    * `r`: Position (vector in R3).
    * `S`: SpinCollection System
    * `state`: system state (Independent, Meanfield, MPC or Quantum)
"""
function intensity(r::Vector{Float64}, S::SpinCollection, state::Union{Vector, ProductState, MPCState, StateVector, DenseOpType})
    N = length(S.spins)
    @assert length(r) == 3

    if isa(state, ProductState)
     @assert N == state.N
     SM = 0.5*(meanfield.sx(state) - 1im*meanfield.sy(state))
    return norm(sum(greenstensor(r-S.spins[i].position)*S.polarizations[i]*SM[i] for i=1:N))^2
    elseif isa(state, MPCState)
        @assert N == state.N
        SM = 0.5*(mpc.sx(state) - 1im*mpc.sy(state))
        return norm(sum(GreenTensor(r-S.spins[i].position)*S.polarizations[i]*SM[i] for i=1:N))^2
    elseif (isa(state, StateVector) || isa(state, DenseOpType))
        # Check which Basis is used
        # ReducedSpinBasis and tensor(Spin(1/2), N) are handeled automatically
        reduced = false
        spin12 = false

        if isa(state, StateVector)
            if isa(state.basis, reducedspin.ReducedSpinBasis)
                reduced = true
            elseif isa(state.basis, CompositeBasis)
                B = [ isa(state.basis.bases[i], SpinBasis{1//2}) for i=1:length(state.basis.bases)]
                @assert length(B) == N
                if (sum(B) == N)
                    spin12 = true
                end
            end
        end

        if isa(state, DenseOpType)
            if isa(state.basis_l, reducedspin.ReducedSpinBasis)
                reduced = true
            elseif isa(state.basis_l, CompositeBasis)
                B = [ isa(state.basis_l.bases[i], SpinBasis{1//2}) for i=1:length(state.basis.bases)]
                @assert length(B) == N
                if (sum(B) == N)
                    spin12 = true
                end
            end
        end

        e(r, i) = GreenTensor(r - S.spins[i].position) * S.polarizations[i]

        if reduced
            b = isa(state, StateVector) ? state.basis : state.basis_l
            smr(i) = reducedspin.reducedsigmam(b, i)
            intensity = sum(dot(e(r, i), e(r, j))*dagger(smr(i))*smr(j) for i=1:N, j=1:N)
            return expect(intensity, state)
        elseif spin12
            # Spin-1/2 Ensemnle
            sm(i) = embed(quantum.basis(S), i, sigmam(SpinBasis(1//2)))
            intensity = sum(dot(e(r, i), e(r, j))*dagger(sm(i))*sm(j) for i=1:N, j=1:N)
            return expect(intensity, state)
        else
            throw("Field cannot be calculated automatically. Please supply a function sm(j) defining the sigma-minus operators.")
        end
    end
end

"""
    field.intensity(r, S, state, sm)

    Arguments:
    * `r`: Position (vector in R3)
    * `S`: SpinCollection System
    * `state`: system state (Independent, Meanfield, MPC or Quantum)
    * `sm`: function sm(j::Int) that gives the sigma-minus operator for the j-th particle.
"""
function intensity(r::Vector{Float64}, S::SpinCollection, state::Union{StateVector, DenseOpType}, sm::Function)
    N = length(S.spins)
    @assert length(r) == 3

    e(r, i) = GreenTensor(r - S.spins[i].position) * S.polarizations[i]
    intensity = sum(dot(e(r, i), e(r, j))*dagger(sm(i))*sm(j) for i=1:N, j=1:N)
     return expect(intensity, state)
end

end # module


text:
using Test
using Plots
using OnlineStats: CovMatrix
using LyapunovExponents
using LyapunovExponents: VecVariance, with_input_noise, with_output_noise
using LyapunovExponents.TestTools: @test_nothrow

@time @testset "Smoke test CLV" begin
    @test_nothrow begin
        prob = LyapunovExponents.lorenz_63(t_attr=3).prob :: LEProblem
        solver = CLVSolver(prob)
        solve!(solver)
    end
    @testset "Recording CLV" begin
        prob = LyapunovExponents.lorenz_63(t_attr=3).prob :: LEProblem
        solver = CLVSolver(
            prob;
            record = [:G, :C, :D, :x])
        solve!(solver)
        @test isdefined(solver.sol, :G_history)
        @test isdefined(solver.sol, :C_history)
        @test isdefined(solver.sol, :D_history)
        @test isdefined(solver.sol, :x_history)
        @test unique(size.(solver.sol.G_history)) == [(3, 3)]
        @test unique(size.(solver.sol.C_history)) == [(3, 3)]
        @test unique(size.(solver.sol.D_history)) == [(3,)]
        @test unique(size.(solver.sol.x_history)) == [(3,)]
    end
end

@time @testset "Smoke test demo" begin
    @test begin
        demo = solve!(LyapunovExponents.lorenz_63(t_attr=3))
        plot(demo, show = false)
        true
    end
    @test begin
        demo = solve!(LyapunovExponents.lorenz_63(t_attr=3, dim_lyap=1))
        plot(demo, show = false)
        true
    end
    @testset "main_stat = $ST" for ST in [VecVariance, CovMatrix]
        demo = LyapunovExponents.lorenz_63(t_attr=3)
        @test_nothrow solve!(demo; main_stat = ST)
        sol = demo.solver.sol
        @test isa(sol.main_stat, ST)
        @test_nothrow plot(demo, show = false)
    end
end

@time @testset "Smoke test RODE" begin
    ode = LyapunovExponents.lorenz_63().prob.phase_prob
    rode = with_input_noise(ode, 1)
    @test_nothrow begin
        t_attr = 3
        prob = LEProblem(rode; t_attr=t_attr)
        sol = solve(prob; record=true, integrator_options=[:dt => 1e-3])
        plot(sol)
    end

    rode = with_output_noise(ode, 1)
    @test_nothrow begin
        t_attr = 3
        prob = LEProblem(rode; t_attr=t_attr)
        sol = solve(prob; record=true, integrator_options=[:dt => 1e-3])
        plot(sol)
    end
end


text:
module PiCamera

export preview

using AbstractPlotting, GLMakie, Makie, VideoIO, CoordinateTransformations, Interact, Observables

# cam_rotation = slider(range(0, stop = 2π, length = 100))
# cam_position = slider(range(0, stop = 2π, length = 100))
# cam_transform = PerspectiveMap() ∘ inv(AffineMap(cam_rotation, cam_position))
# screen_points = map(cam_transform, points)
function preview()
    VideoIO.viewcam()
end

end # module


text:
using BenchmarkTools
using Particles, Distributions, ConjugatePriors
using ConjugatePriors: NormalInverseChisq

truth = MixtureModel([Normal(-2.), Normal(2.)])
y = rand(truth, 100)

function f(y)
    ps = ChenLiuParticles(100,
                          NormalInverseChisq(0., 1., 0.1, 1.0),
                          ChineseRestaurantProcess(1.))
    filter!(ps, y, false)
end

f(y)
@btime(f($y))

@profiler f(y)

function g(y)
    ps = FearnheadParticles(100,
                            NormalInverseChisq(0., 1., 0.1, 1.0),
                            ChineseRestaurantProcess(1.))
    filter!(ps, y, false)
end

g(y)
@btime(g($y))

@profiler g(y)