text:
type TwoSeasonalArms <: BanditEnvironment
  season_length::Int64
end

function arms(env::TwoSeasonalArms, trial::Int64)
  if div(trial, env.season_length) % 2 == 0
    [BernoulliArm(0.1), BernoulliArm(0.9)]
  else
    [BernoulliArm(0.9), BernoulliArm(0.1)]
  end
end

function n_arms(env::TwoSeasonalArms)
  2
end

max_stars_repo_name: rmminusrslash/BanditsBook

max_issues_repo_name: rmminusrslash/BanditsBook

max_forks_repo_name: rmminusrslash/BanditsBook


text:
const _pkg_root = dirname(dirname(@__FILE__))
const _pkg_assets = joinpath(_pkg_root,"assets")

!isdir(_pkg_assets) && mkdir(_pkg_assets)

deps = [
    "https://www.gitcdn.xyz/repo/piever/InteractResources/v0.3.0/bulma/main.min.css" => "bulma.min.css",
    "https://www.gitcdn.xyz/repo/piever/InteractResources/v0.3.0/bulma/main_confined.min.css" => "bulma_confined.min.css",
]

for (dep, name) in deps
    download(dep, joinpath(_pkg_assets, name))
end

ext: jl


text:
#ДАНО: Робот - в произвольной клетке ограниченного прямоугольной рамкой поля без внутренних перегородок и маркеров.
#РЕЗУЛЬТАТ: Робот - в исходном положении в центре косого креста (в форме X) из маркеров.
function mark_x(r)
    for i in 0:3
        mark_line(r,i)
    end
    putmarker!(r)
end
function mark_line(r,i)
    while isborder(r,HorizonSide(i % 4))==0 && isborder(r,HorizonSide((i+1) % 4))==0
        move!(r,HorizonSide(i % 4))
        move!(r,HorizonSide((i+1) % 4))
        putmarker!(r)
    end
    while ismarker(r)==1
        move!(r,HorizonSide((i+2) % 4))
        move!(r,HorizonSide((i+3) % 4))
    end
end
mark_x(r)


text:
function _basic_routing_model_unitary_flow_conservation_flow_formulation(m::Model, rd::RoutingData, routing, f)
    # Write the flow-conservation constraints, the rhs being given by f. This generalises both unitary and unscaled models.
    cstrs_src_in = Dict{Edge{Int}, Any}()
    cstrs_src_out = Dict{Edge{Int}, Any}()
    cstrs_dst_in = Dict{Edge{Int}, Any}()
    cstrs_dst_out = Dict{Edge{Int}, Any}()
    cstrs_bal = Dict{Edge{Int}, Dict{Int, Any}}()

    for d in demands(rd)
        @assert length(outedges(rd.g, d.src)) > 0
        cstrs_src_out[d] = @constraint(m, sum(routing[d, e] for e in outedges(rd.g, d.src)) == f(d))
        set_name(cstrs_src_out[d], "c_conservation_source_out[$(d)]")

        if length(inedges(rd.g, d.src)) > 0
            cstrs_src_in[d] = @constraint(m, sum(routing[d, e] for e in inedges(rd.g, d.src)) == 0)
            set_name(cstrs_src_in[d], "c_conservation_source_in[$(d)]")
        end

        @assert length(inedges(rd.g, d.dst)) > 0
        cstrs_dst_in[d] = @constraint(m, sum(routing[d, e] for e in inedges(rd.g, d.dst)) == f(d))
        set_name(cstrs_dst_in[d], "c_conservation_target_in[$(d)]")

        if length(outedges(rd.g, d.dst)) > 0
            cstrs_dst_out[d] = @constraint(m, sum(routing[d, e] for e in outedges(rd.g, d.dst)) == 0)
            set_name(cstrs_dst_out[d], "c_conservation_target_out[$(d)]")
        end

        cstrs_bal[d] = Dict{Int, Any}()
        for v in vertices(rd.g)
            if v == d.src || v == d.dst
                continue
            end

            cstrs_bal[d][v] = @constraint(m, sum(routing[d, e] for e in inedges(rd.g, v)) == sum(routing[d, e] for e in outedges(rd.g, v)))
            set_name(cstrs_bal[d][v], "c_conservation[$(d), $(v)]")
        end
    end

    return cstrs_src_in, cstrs_src_out, cstrs_dst_in, cstrs_dst_out, cstrs_bal
end

function basic_routing_model_unitary(m::Model, rd::RoutingData, ::FlowFormulation)
    @variable(m, 0 <= routing[d in demands(rd), e in edges(rd)] <= 1)

    cstrs_src_in, cstrs_src_out, cstrs_dst_in, cstrs_dst_out, cstrs_bal =
        _basic_routing_model_unitary_flow_conservation_flow_formulation(m, rd, routing, d -> 1)

    return RoutingModel(rd, m, UnitaryFlows, routing,
        constraints_source_in=cstrs_src_in, constraints_source_out=cstrs_src_out,
        constraints_target_in=cstrs_dst_in, constraints_target_out=cstrs_dst_out,
        constraints_balance=cstrs_bal)
end

function basic_routing_model_unscaled(m::Model, rd::RoutingData, dm, ::FlowFormulation)
    # dm is either a traffic matrix or variables! The only guarantee is that it is indexable by demands.

    @variable(m, routing[d in demands(rd), e in edges(rd)] >= 0)

    cstrs_src_in, cstrs_src_out, cstrs_dst_in, cstrs_dst_out, cstrs_bal =
        _basic_routing_model_unitary_flow_conservation_flow_formulation(m, rd, routing, d -> dm[d])

    return RoutingModel(rd, m, UnscaledFlows, routing,
        constraints_source_in=cstrs_src_in, constraints_source_out=cstrs_src_out,
        constraints_target_in=cstrs_dst_in, constraints_target_out=cstrs_dst_out,
        constraints_balance=cstrs_bal)
end

function total_flow_in_edge(rm::RoutingModel, e::Edge, dm::Dict{Edge{Int}, Float64}, ::FlowFormulation, ::Val{UnitaryFlows})
    flow = AffExpr(0.0)
    for d in keys(dm)
        if ! iszero(dm[d])
            add_to_expression!(flow, dm[d] * rm.routing[d, e])
        end
    end
    return flow
end

function total_flow_in_edge(rm::RoutingModel, e::Edge, ::FlowFormulation, ::Val{UnscaledFlows})
    return sum(rm.routing[d, e] for d in demands(rm.data))
end


text:
function test_derivatives_projections_1d(k::Int, N::Int, el::Int)
    f(x) = sin(2π * x)
    df(x) = 2π * cos(2π * x)
    interval = 0..1

    print("Periodic B-Spline Basis:\n")
    basis = PeriodicBSpline(interval, k, N)

    M = basis ⊙ basis
    M_d = mass_matrix(basis')
    D = der_mat(basis)

    fh = M \ (f ⊙ basis)
    dfh = M_d \ (df ⊙ basis')

    @show maximum(abs.( D * fh - dfh))

    x, colmat = col_mat(basis, 100)
    x, colmat_d = col_mat(basis', 100)

    @show maximum(abs.(colmat * fh - f.(x)))
    @show maximum(abs.(colmat_d * dfh - df.(x)))
    @show maximum(abs.(colmat_d * D * fh - df.(x)))

    S = stiffness_matrix(basis')
    @show maximum(abs.(D' * M_d * D - S))

    print("Dirichlet B-Spline Basis:\n")
    basis = DirichletBSpline(interval, k, N)

    M = mass_matrix(basis)
    M_d = mass_matrix(basis')
    D = SplineFEEC.der_mat(basis)

    fh = M \ (f ⊙ basis)
    dfh = M_d \ (df ⊙ basis')

    @show maximum(abs.( D * fh - dfh))

    x, colmat = col_mat(basis, 100)
    x, colmat_d = col_mat(basis', 100)

    @show maximum(abs.(colmat * fh - f.(x)))
    @show maximum(abs.(colmat_d * dfh - df.(x)))
    @show maximum(abs.(colmat_d * D * fh - df.(x)))

    S = stiffness_matrix(basis')
    @show maximum(abs.(D' * M_d * D - S))

    print("Multi-patch Dirichlet B-Spline Basis: \n")

    Δx = (interval.right - interval.left)/el
    basis = DirichletMultipatchBSpline(interval, k, N, el)

    M = mass_matrix(basis)
    M_d =  mass_matrix_d(basis')
    D = der_mat(basis)
    fh = M \ L2_prod(f, basis)
    dfh = M_d\ L2_prod_d(df, basis')

    @show maximum(abs.( D * fh - dfh))

    x, colmat = col_mat(basis, 100)
    x_d, colmat_d = col_mat_d(basis', 100)

    @show maximum(abs.(colmat * fh - f.(x)))
    @show maximum(abs.(colmat_d * dfh - df.(x_d)))
    @show maximum(abs.(colmat_d * D * fh - df.(x_d)))
    S = stiffness_matrix(basis')

    @show maximum(abs.( D' * M_d * D - S))

    print("Multi-patch Periodic B-Spline Basis: \n")

    Δx = (interval.right - interval.left)/el
    basis = PeriodicMultipatchBSpline(interval, k, N, el)

    M = mass_matrix(basis)
    M_d =  mass_matrix_d(basis')
    D = der_mat(basis)
    fh = M \ L2_prod(f, basis)
    dfh = M_d\ L2_prod_d(df, basis')

    @show maximum(abs.( D * fh - dfh))

    x, colmat = col_mat(basis, 100)
    x_d, colmat_d = col_mat_d(basis', 100)

    @show maximum(abs.(colmat * fh - f.(x)))
    @show maximum(abs.(colmat_d * dfh - df.(x_d)))
    @show maximum(abs.(colmat_d * D * fh - df.(x_d)))
    S = stiffness_matrix(basis')

    @show maximum(abs.( D' * M_d * D - S))
end

function test_derivatives_projections_2d(k::NTuple{2,Int}, N::NTuple{2,Int})
    f(x,y) = sin(2π * x)*sin(2π * y)
    gradf(x,y) = 2π * SVector(cos(2π * x)*sin(2π * y),
                                sin(2π * x)*cos(2π * y))

    g(x,y) = SVector(sin(2π * x)*sin(2π * y),
                        sin(2π * y)*sin(2π * x))
    curlg(x,y) = 2π * (sin(2π * y)*cos(2π * x) -  sin(2π * x)*cos(2π * y))
    domain = (0..1) × (0..1)

    print("Periodic B-Spline Basis: \n")
    X, Y = PeriodicBSplineTensorProductBasis(domain, k, N).B
    V0 = X ⊗ Y
    V1 = SVector(X' ⊗ Y , X ⊗ Y')
    V2 = X' ⊗ Y'

    G = disc_grad(V0)
    C = disc_curl(V0)

    xy0, colmat0 = col_mat(V0, 100)
    (xy1_1, colmat1_1), (xy1_2, colmat1_2) = col_mat.(V1, 100)
    xy2, colmat2 = col_mat(V2, 100)

    fh = L2_proj(f, V0)
    gradfh = L2_proj(gradf, V1)

    @show maximum(abs.(vcat((G * fh - gradfh)...)))

    val_f = (i -> f.(xy0[i][1], xy0[i][2])).(1:length(xy0))
    @show maximum(abs.(colmat0 * fh - val_f))

    val_gradf1 = (i -> gradf.(xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_gradf2 = (i -> gradf.(xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
    #return fh, gradfh, colmat1_1, val_gradf1

    @show maximum(abs.(colmat1_1 * gradfh[1] - val_gradf1))
    @show maximum(abs.(colmat1_2 * gradfh[2] - val_gradf2))

    gh = L2_proj(g, V1)
    curlgh = L2_proj(curlg, V2)
    @show maximum(abs.((C * gh - curlgh)))

    val_g1 = (i -> g.(xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_g2 = (i -> g.(xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
    @show maximum(abs.(colmat1_1 * gh[1] - val_g1))
    @show maximum(abs.(colmat1_2 * gh[2] - val_g2))

    val_curlg = (i -> curlg.(xy2[i][1], xy2[i][2])).(1:length(xy2))
    @show maximum(abs.(colmat2 * curlgh - val_curlg))

    print("Dirichlet B-Spline Basis: \n")
    X, Y = DirichletBSplineTensorProductBasis(domain, k, N).B
    V0 = X ⊗ Y
    V1 = SVector(X' ⊗ Y , X ⊗ Y')
    V2 = X' ⊗ Y'

    G = disc_grad(V0)
    C = disc_curl(V0)

    xy0, colmat0 = col_mat(V0, 100)
    (xy1_1, colmat1_1), (xy1_2, colmat1_2) = col_mat.(V1, 100)
    xy2, colmat2 = col_mat(V2, 100)

    fh = L2_proj(f, V0)
    gradfh = L2_proj(gradf, V1)

    @show maximum(abs.(vcat((G * fh - gradfh)...)))

    val_f = (i -> f.(xy0[i][1], xy0[i][2])).(1:length(xy0))
    @show maximum(abs.(colmat0 * fh - val_f))

    val_gradf1 = (i -> gradf.(xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_gradf2 = (i -> gradf.(xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))

    @show maximum(abs.(colmat1_1 * gradfh[1] - val_gradf1))
    @show maximum(abs.(colmat1_2 * gradfh[2] - val_gradf2))

    gh = L2_proj(g, V1)
    curlgh = L2_proj(curlg, V2)
    @show maximum(abs.((C * gh - curlgh)))

    val_g1 = (i -> g.(xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_g2 = (i -> g.(xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
    @show maximum(abs.(colmat1_1 * gh[1] - val_g1))
    @show maximum(abs.(colmat1_2 * gh[2] - val_g2))

    val_curlg = (i -> curlg.(xy2[i][1], xy2[i][2])).(1:length(xy2))
    @show maximum(abs.(colmat2 * curlgh - val_curlg))

    #test collocation projections
    xy0, colmat0 = col_mat_proj(V0)
    (xy1_1, colmat1_1), (xy1_2, colmat1_2) = col_mat_proj.(V1)
    xy2, colmat2 = col_mat_proj(V2)

    val_f = (i -> f.(xy0[i][1], xy0[i][2])).(1:length(xy0))
    @show maximum(abs.(colmat0\val_f - fh))

    val_gradf1 = (i -> gradf.(xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_gradf2 = (i -> gradf.(xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
    @show maximum(abs.(colmat1_1\val_gradf1 - gradfh[1]))
    @show maximum(abs.(colmat1_2\val_gradf2 - gradfh[2]))

    val_g1 = (i -> g.(xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_g2 = (i -> g.(xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
    @show maximum(abs.(colmat1_1\val_g1 - gh[1]))
    @show maximum(abs.(colmat1_2\val_g2 - gh[2]))

    val_curlg = (i -> curlg.(xy2[i][1], xy2[i][2])).(1:length(xy2))
    @show maximum(abs.(colmat2\val_curlg - curlgh))

end

function test_derivatives_projections_Multipatch_2d(k::NTuple{2,Int}, N::NTuple{2,Int}, el::NTuple{2,Int})
    f(x,y) = sin(2π * x)*sin(2π * y)
    gradf(x,y) = 2π * SVector(cos(2π * x)*sin(2π * y),
                                sin(2π * x)*cos(2π * y))

    g(x,y) = SVector(sin(2π * x)*sin(2π * y),
                        sin(2π * y)*sin(2π * x))
    curlg(x,y) = 2π * (sin(2π * y)*cos(2π * x) -  sin(2π * x)*cos(2π * y))
    domain = (0..1) × (0..1)

    print("Multipatch Dirichlet B-Spline Basis: \n")
    X, Y = DirichletMultipatchTensorProductBSpline(domain, k, N, el).B
    V0 = X ⊗ Y
    V1 = SVector(X' ⊗ Y , X ⊗ Y')
    V2 = X' ⊗ Y'

    G = disc_diff(grad, V0)
    C = disc_diff(curl, V0)

    xy0, colmat0 = col_mat(V0, 100)

    #need a mixed setup
    x, cmx = col_mat_d(X', 100)
    y, cmy = col_mat(Y, 100)
    xy1_1 = [(a,b) for a in x for b in y]
    colmat1_1 = cmx ⊗ cmy

    x, cmx = col_mat(X, 100)
    y, cmy = col_mat_d(Y', 100)
    xy1_2 = [(a,b) for a in x for b in y]
    colmat1_2 = cmx ⊗ cmy

    xy2, colmat2 = col_mat_d(V2, 100)

    fh = L2_proj(f, V0)

    #also mixed setup TODO introduce better functions
    M1x = mass_matrix_d(X')
    M1y = mass_matrix(Y)
    M1_1 = M1x ⊗ M1y
    gradfh1 = M1_1\L2_prod_mixed_x((x,y) -> gradf(x,y)[1], X' ⊗ Y)

    M2x = mass_matrix(X)
    M2y = mass_matrix_d(Y')
    M1_2 = M2x ⊗ M2y
    gradfh2 = M1_2\L2_prod_mixed_y((x,y) -> gradf(x,y)[2], X ⊗ Y')
    gradfh = SVector(gradfh1, gradfh2)
    @show maximum(abs.(vcat((G * fh - gradfh)...)))

    val_f = (i -> f.(xy0[i][1], xy0[i][2])).(1:length(xy0))
    @show maximum(abs.(colmat0 * fh - val_f))

    val_gradf1 = (i -> gradf.(xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_gradf2 = (i -> gradf.(xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))

    @show maximum(abs.(colmat1_1 * gradfh[1] - val_gradf1))
    @show maximum(abs.(colmat1_2 * gradfh[2] - val_gradf2))

    #gh = L2_proj(g, V1)
    gh = SVector( M1_1\L2_prod_mixed_x((x,y) -> g(x,y)[1], X' ⊗ Y),
                  M1_2\L2_prod_mixed_y((x,y) -> g(x,y)[2], X ⊗ Y'))

    curlgh = L2_proj_d(curlg, V2)
    @show maximum(abs.((C * gh - curlgh)))

    val_g1 = (i -> g.(xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_g2 = (i -> g.(xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
    @show maximum(abs.(colmat1_1 * gh[1] - val_g1))
    @show maximum(abs.(colmat1_2 * gh[2] - val_g2))

    val_curlg = (i -> curlg.(xy2[i][1], xy2[i][2])).(1:length(xy2))
    @show maximum(abs.(colmat2 * curlgh - val_curlg))

    #test collocation projections
    xy0, colmat0 = col_mat_proj(V0)
    xy1_1, colmat1_1 = col_mat_proj_mixed_x(V1[1])
    xy1_2, colmat1_2 = col_mat_proj_mixed_y(V1[2])
    xy2, colmat2 = col_mat_proj_d(V2)

    val_f = (i -> f.(xy0[i][1], xy0[i][2])).(1:length(xy0))
    @show maximum(abs.(colmat0\val_f - fh))

    val_gradf1 = (i -> gradf.(xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_gradf2 = (i -> gradf.(xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
    @show maximum(abs.(colmat1_1\val_gradf1 - gradfh[1]))
    @show maximum(abs.(colmat1_2\val_gradf2 - gradfh[2]))

    val_g1 = (i -> g.(xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_g2 = (i -> g.(xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
    @show maximum(abs.(colmat1_1\val_g1 - gh[1]))
    @show maximum(abs.(colmat1_2\val_g2 - gh[2]))

    val_curlg = (i -> curlg.(xy2[i][1], xy2[i][2])).(1:length(xy2))
    @show maximum(abs.(colmat2\val_curlg - curlgh))

    print("\nConga Dirichlet B-Spline Basis: \n")

    Δx = (domain.domains[1].right - domain.domains[1].left)/el[1]
    Δy = (domain.domains[2].right - domain.domains[2].left)/el[2]

    X, Y = DirichletMultipatchTensorProductBSpline(domain, k, N, el).B
    V0 = X ⊗ Y
    V1 = SVector(X' ⊗ Y , X ⊗ Y')
    V2 = X' ⊗ Y'

    G = disc_diff(grad, V0)
    C = disc_diff(curl, V0)
    P = conforming_averaging_projection_matrix(V0)
    Px = conforming_averaging_projection_matrix_x(V0)
    Py = conforming_averaging_projection_matrix_y(V0)

    xy0, colmat0 = col_mat(V0, 100)
    xy1_1, colmat1_1 = col_mat_mixed_x(V1[1], 100)
    xy1_2, colmat1_2 = col_mat_mixed_y(V1[2], 100)
    xy2, colmat2 = col_mat_d(V2, 100)

    fh = L2_proj_d(f, V0)
    fh *= 1/(Δx*Δy)

    gradfh = L2_proj_d(gradf, V1)
    gradfh = SVector(1/Δy *gradfh[1],1/Δx * gradfh[2])
    @show maximum(abs.(vcat(( G * SVector{size(P)[1]}(P * fh) - SVector(Py * gradfh[1],Px * gradfh[2]) )...)))

    val_f = (i -> f.(xy0[i][1], xy0[i][2])).(1:length(xy0))
    @show maximum(abs.(colmat0 * P * fh - val_f))

    val_gradf1 = (i -> gradf.(xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_gradf2 = (i -> gradf.(xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
    @show maximum(abs.(colmat1_1 * Py * gradfh[1] - val_gradf1))
    @show maximum(abs.(colmat1_2 * Px * gradfh[2] - val_gradf2))

    gh = L2_proj_d(g, V1)
    gh = SVector(1/Δy * gh[1], 1/Δx * gh[2])
    curlgh = L2_proj_d(curlg, V2)
    @show maximum(abs.((C * SVector(Py * gh[1], Px * gh[2] ) - curlgh)))

    val_g1 = (i -> g.(xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_g2 = (i -> g.(xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
    @show maximum(abs.(colmat1_1 * Py * gh[1] - val_g1))
    @show maximum(abs.(colmat1_2 * Px * gh[2] - val_g2))

    val_curlg = (i -> curlg.(xy2[i][1], xy2[i][2])).(1:length(xy2))
    @show maximum(abs.(colmat2 * curlgh - val_curlg))

    #testing the colocation projection
    xy0, colmat0 = col_mat_proj_d(V0)
    xy1_1, colmat1_1 = col_mat_proj_d(V1[1])
    xy1_2, colmat1_2 = col_mat_proj_d(V1[2])
    xy2, colmat2 = col_mat_proj_d(V2)

    val_f = (i -> f.(xy0[i][1], xy0[i][2])).(1:length(xy0))
    @show maximum(abs.((Δx*Δy) * colmat0\val_f - fh))

    val_gradf1 = (i -> gradf.(xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_gradf2 = (i -> gradf.(xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
    @show maximum(abs.( Δy * colmat1_1\val_gradf1 -  gradfh[1]))
    @show maximum(abs.( Δx * colmat1_2\val_gradf2 -  gradfh[2]))

    val_g1 = (i -> g.(xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_g2 = (i -> g.(xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
    @show maximum(abs.(Δy *colmat1_1\val_g1 - gh[1]))
    @show maximum(abs.(Δx *colmat1_2\val_g2 - gh[2]))

    val_curlg = (i -> curlg.(xy2[i][1], xy2[i][2])).(1:length(xy2))
    @show maximum(abs.(colmat2\val_curlg - curlgh))

    print("\nMultipatch Periodic B-Spline Basis: \n")
    X, Y = SplineFEEC.PeriodicMultipatchTensorProductBSpline(domain, k, N, el).B
    V0 = X ⊗ Y
    V1 = SVector(X' ⊗ Y , X ⊗ Y')
    V2 = X' ⊗ Y'

    G = disc_diff(grad, V0)
    C = disc_diff(curl, V0)

    xy0, colmat0 = col_mat(V0, 100)

    #need a mixed setup
    x, cmx = col_mat_d(X', 100)
    y, cmy = col_mat(Y, 100)
    xy1_1 = [(a,b) for a in x for b in y]
    colmat1_1 = cmx ⊗ cmy

    x, cmx = col_mat(X, 100)
    y, cmy = col_mat_d(Y', 100)
    xy1_2 = [(a,b) for a in x for b in y]
    colmat1_2 = cmx ⊗ cmy

    xy2, colmat2 = col_mat_d(V2, 100)

    fh = L2_proj(f, V0)

    #also mixed setup TODO introduce better functions
    M1x = mass_matrix_d(X')
    M1y = mass_matrix(Y)
    M1_1 = M1x ⊗ M1y
    gradfh1 = M1_1\SplineFEEC.L2_prod_mixed_x((x,y) -> gradf(x,y)[1], X' ⊗ Y)

    M2x = mass_matrix(X)
    M2y = mass_matrix_d(Y')
    M1_2 = M2x ⊗ M2y
    gradfh2 = M1_2\SplineFEEC.L2_prod_mixed_y((x,y) -> gradf(x,y)[2], X ⊗ Y')
    gradfh = SVector(gradfh1, gradfh2)
    @show maximum(abs.(vcat((G * fh - gradfh)...)))

    val_f = (i -> f.(xy0[i][1], xy0[i][2])).(1:length(xy0))
    @show maximum(abs.(colmat0 * fh - val_f))

    val_gradf1 = (i -> gradf.(xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_gradf2 = (i -> gradf.(xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))

    @show maximum(abs.(colmat1_1 * gradfh[1] - val_gradf1))
    @show maximum(abs.(colmat1_2 * gradfh[2] - val_gradf2))

    gh = SVector( M1_1\SplineFEEC.L2_prod_mixed_x((x,y) -> g(x,y)[1], X' ⊗ Y),
                  M1_2\SplineFEEC.L2_prod_mixed_y((x,y) -> g(x,y)[2], X ⊗ Y'))

    curlgh = L2_proj_d(curlg, V2)
    @show maximum(abs.((C * gh - curlgh)))

    val_g1 = (i -> g.(xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_g2 = (i -> g.(xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
    @show maximum(abs.(colmat1_1 * gh[1] - val_g1))
    @show maximum(abs.(colmat1_2 * gh[2] - val_g2))

    val_curlg = (i -> curlg.(xy2[i][1], xy2[i][2])).(1:length(xy2))
    @show maximum(abs.(colmat2 * curlgh - val_curlg))

    #test collocation projections
    xy0, colmat0 = SplineFEEC.col_mat_proj(V0)
    xy1_1, colmat1_1 = SplineFEEC.col_mat_proj_mixed_x(V1[1])
    xy1_2, colmat1_2 = SplineFEEC.col_mat_proj_mixed_y(V1[2])
    xy2, colmat2 = SplineFEEC.col_mat_proj_d(V2)

    val_f = (i -> f.(xy0[i][1], xy0[i][2])).(1:length(xy0))
    @show maximum(abs.(colmat0\val_f - fh))

    val_gradf1 = (i -> gradf.(xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_gradf2 = (i -> gradf.(xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
    @show maximum(abs.(colmat1_1\val_gradf1 - gradfh[1]))
    @show maximum(abs.(colmat1_2\val_gradf2 - gradfh[2]))

    val_g1 = (i -> g.(xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_g2 = (i -> g.(xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
    @show maximum(abs.(colmat1_1\val_g1 - gh[1]))
    @show maximum(abs.(colmat1_2\val_g2 - gh[2]))

    val_curlg = (i -> curlg.(xy2[i][1], xy2[i][2])).(1:length(xy2))
    @show maximum(abs.(colmat2\val_curlg - curlgh))

    print("\nConga Periodic B-Spline Basis: \n")

    Δx = (domain.domains[1].right - domain.domains[1].left)/el[1]
    Δy = (domain.domains[2].right - domain.domains[2].left)/el[2]

    X, Y = SplineFEEC.PeriodicMultipatchTensorProductBSpline(domain, k, N, el).B
    V0 = X ⊗ Y
    V1 = SVector(X' ⊗ Y , X ⊗ Y')
    V2 = X' ⊗ Y'

    G = SplineFEEC.disc_diff(grad, V0)
    C = SplineFEEC.disc_diff(curl, V0)
    P = SplineFEEC.conforming_averaging_projection_matrix(V0)
    Px = SplineFEEC.conforming_averaging_projection_matrix_x(V0)
    Py = SplineFEEC.conforming_averaging_projection_matrix_y(V0)

    xy0, colmat0 = col_mat(V0, 100)
    xy1_1, colmat1_1 = SplineFEEC.col_mat_mixed_x(V1[1], 100)
    xy1_2, colmat1_2 = SplineFEEC.col_mat_mixed_y(V1[2], 100)
    xy2, colmat2 = col_mat_d(V2, 100)

    fh = L2_proj_d(f, V0)
    fh *= 1/(Δx*Δy)

    gradfh = L2_proj_d(gradf, V1)
    gradfh = SVector(1/Δy *gradfh[1],1/Δx * gradfh[2])
    @show maximum(abs.(vcat(( G * SVector{size(P)[1]}(P * fh) - SVector(Py * gradfh[1],Px * gradfh[2]) )...)))

    val_f = (i -> f.(xy0[i][1], xy0[i][2])).(1:length(xy0))
    @show maximum(abs.(colmat0 * P * fh - val_f))

    val_gradf1 = (i -> gradf.(xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_gradf2 = (i -> gradf.(xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
    @show maximum(abs.(colmat1_1 * Py * gradfh[1] - val_gradf1))
    @show maximum(abs.(colmat1_2 * Px * gradfh[2] - val_gradf2))

    gh = L2_proj_d(g, V1)
    gh = SVector(1/Δy * gh[1], 1/Δx * gh[2])
    curlgh = L2_proj_d(curlg, V2)
    @show maximum(abs.((C * SVector(Py * gh[1], Px * gh[2] ) - curlgh)))

    val_g1 = (i -> g.(xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_g2 = (i -> g.(xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
    @show maximum(abs.(colmat1_1 * Py * gh[1] - val_g1))
    @show maximum(abs.(colmat1_2 * Px * gh[2] - val_g2))

    val_curlg = (i -> curlg.(xy2[i][1], xy2[i][2])).(1:length(xy2))
    @show maximum(abs.(colmat2 * curlgh - val_curlg))

    #testing the colocation projection
    xy0, colmat0 = SplineFEEC.col_mat_proj_d(V0)
    xy1_1, colmat1_1 = SplineFEEC.col_mat_proj_d(V1[1])
    xy1_2, colmat1_2 = SplineFEEC.col_mat_proj_d(V1[2])
    xy2, colmat2 = SplineFEEC.col_mat_proj_d(V2)

    val_f = (i -> f.(xy0[i][1], xy0[i][2])).(1:length(xy0))
    @show maximum(abs.((Δx*Δy) * colmat0\val_f - fh))

    val_gradf1 = (i -> gradf.(xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_gradf2 = (i -> gradf.(xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
    @show maximum(abs.( Δy * colmat1_1\val_gradf1 -  gradfh[1]))
    @show maximum(abs.( Δx * colmat1_2\val_gradf2 -  gradfh[2]))

    val_g1 = (i -> g.(xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_g2 = (i -> g.(xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
    @show maximum(abs.(Δy *colmat1_1\val_g1 - gh[1]))
    @show maximum(abs.(Δx *colmat1_2\val_g2 - gh[2]))

    val_curlg = (i -> curlg.(xy2[i][1], xy2[i][2])).(1:length(xy2))
    @show maximum(abs.(colmat2\val_curlg - curlgh))
end

function test_Poisson_1d(k::Int, N::Int, el::Int)
    u(x) = sin(2π * x)
    f(x) = -4π^2 * sin(2π * x)
    interval = 0..1

    print("Periodic B-Spline Basis: \n")
    basis = PeriodicBSpline(interval, k, N)

    Δ = stiffness_matrix(basis')
    R = 1/N .* circ_ones(N)
    P = circ_id(N) - R
    fh = L2_prod(f, basis)
    uh = -(Δ + R) \ (P * fh)

    x, colmat = col_mat(basis, 200)
    @show maximum( abs.(colmat * uh .- u.(x) ))

    print("Dirichlet B-Spline Basis: \n")
    basis = DirichletBSpline(interval, k, N)

    fh = L2_prod(f, basis)
    # enforce homogenous boundary conditions
    Δ = stiffness_matrix(basis')
    ϵ = 1e-12
    Δ[1,1] *= (1 + 1/ϵ)
    Δ[end,end] *= (1 + 1/ϵ)

    uh = -Δ \ fh

    x, colmat = col_mat(basis, 200)
    @show maximum( abs.(colmat * uh .- u.(x) ))

    print("Multi-patch Dirichlet B-Spline Basis: \n")
    basis = SplineFEEC.DirichletMultipatchBSpline(interval, k, N, el)

    fh = L2_prod(f, basis)
    Δ = stiffness_matrix(basis')
    #dirichlet bc
    ϵ = 1e-12
    Δ[1,1] *= (1 +1/ϵ)
    Δ[end,end] *= (1 +1/ϵ)

    uh = -Δ \ fh

    x, colmat = col_mat(basis, 200)
    @show maximum( abs.(colmat * uh .- u.(x) ))
end

function test_Maxwells_1d(k::Int, N::Int, el::Int)
    # exact solution
    B = (t, x) -> cos(2π * t) * cos(2π * x)
    E = (t, x) -> sin(2π * t) * sin(2π * x)
    J = (t, x) -> 0
    domain = (0..1)
    it = 1234
    Δt = 0.001
    T = Δt * it
    Δx = (domain.right - domain.left)/el

    @show (k, N, el)
    @show Δt, it

    print("\nPeriodic B-Spline Basis: \n")
    # define the bases
    X = PeriodicBSpline(domain, k, N)
    V0 = X
    V1 = X'

    # calculate the mass matrices
    M0 = mass_matrix(V0)
    M0inv = inv(M0) # using the circulant inverse
    M1 = mass_matrix(V1)

    # calculate the discrete derivative
    G = disc_grad(V0)

    # discrete projections
    Bh = L2_proj(x -> B(0, x), V1)
    Eh = L2_proj(x -> E(0, x), V0)

    # simple leap frog scheme
    for i in 1:it
        Eh += Δt/2 * M0inv * (G' * (M1 * Bh))
        Bh -= Δt * (G * Eh)
        Eh += Δt/2 * M0inv * (G' * (M1 * Bh))
    end

    # collocation matrices
    x0, cm_0 = col_mat(V0, 100)
    x1, cm_1 = col_mat(V1, 100)

    @show maximum(abs.(cm_0 * Eh - E.(T, x0)))
    @show maximum(abs.(cm_1 * Bh - B.(T, x1)))

    print("\nDirichlet B-Spline Basis: \n")
    X = DirichletBSpline(domain, k, N)
    V0 = X
    V1 = X'

    # calculate the mass matrices
    M0 = mass_matrix(V0)

    ϵ = 1e-14
    M0[1,1] *= (1 + 1/ϵ)
    M0[end,end] *= (1 + 1/ϵ)

    M0inv = inv(M0) # using the circulant inverse
    M1 = mass_matrix(V1)

    # calculate the discrete derivative
    G = disc_grad(V0)

    # discrete projections
    Bh = L2_proj(x -> B(0, x), V1)
    Eh = L2_proj(x -> E(0, x), V0)

    # simple leap frog scheme
    for i in 1:it
        Eh += Δt/2 * M0inv * (G' * (M1 * Bh))
        Bh -= Δt * (G * Eh)
        Eh += Δt/2 * M0inv * (G' * (M1 * Bh))
    end

    # collocation matrices
    x0, cm_0 = col_mat(V0, 100)
    x1, cm_1 = col_mat(V1, 100)

    @show maximum(abs.(cm_0 * Eh - E.(T, x0)))
    @show maximum(abs.(cm_1 * Bh - B.(T, x1)))

    print("\nMulti-patch Dirichlet B-Spline Basis: \n")
    X = DirichletMultipatchBSpline(domain, k, N, el)
    V0 = X
    V1 = X'

    # calculate the mass matrices
    M0 = mass_matrix(V0)

    ϵ = 1e-14
    M0[1,1] *= (1 + 1/ϵ)
    M0[end,end] *= (1 + 1/ϵ)

    M0inv = inv(M0) # using the circulant inverse
    M1 = mass_matrix_d(V1) #special mass matrix

    # calculate the discrete derivative
    G = disc_grad(V0)

    # discrete projections
    Eh = L2_proj(x -> E(0, x), V0)
    Bh = L2_proj_d(x -> B(0, x), V1)

    # simple leap frog scheme
    for i in 1:it
        Eh += Δt/2 * M0inv * (G' * (M1 * Bh))
        Bh -= Δt * (G * Eh)
        Eh += Δt/2 * M0inv * (G' * (M1 * Bh))
    end

    # collocation matrices
    x0, cm_0 = col_mat(V0, 100)
    x1, cm_1 = col_mat_d(V1, 100)

    @show maximum(abs.(cm_0 * Eh - E.(T, x0)))
    @show maximum(abs.(cm_1 * Bh - B.(T, x1)))

    print("\nConga Dirichlet B-Spline Basis: \n")
    X = DirichletMultipatchBSpline(domain, k, N, el)
    V0 = X
    V1 = X'

    # calculate the mass matrices
    M0 = mass_matrix_d(V0)

    ϵ = 1e-14
    M0[1,1] *= (1 + 1/ϵ)
    M0[end,end] *= (1 + 1/ϵ)

    M0inv = 1/Δx^2 * inv(M0) # using the circulant inverse
    M1 = mass_matrix_d(V1) #special mass matrix

    # calculate the discrete derivative
    G = disc_grad(V0)
    P = conforming_averaging_projection_matrix(V0)

    # discrete projections
    Eh = 1/Δx * L2_proj_d(x -> E(0, x), V0)
    Bh = L2_proj_d(x -> B(0, x), V1)

    # simple leap frog scheme
    for i in 1:it
        Eh += Δt/2 * M0inv * (P' * G' * (M1 * Bh))
        Bh -= Δt * (G * P * Eh)
        Eh += Δt/2 * M0inv * (P' * G' * (M1 * Bh))
    end

    # collocation matrices
    x0, cm_0 = col_mat(V0, 100)
    x1, cm_1 = col_mat_d(V1, 100)

    @show maximum(abs.(cm_0 * P * Eh - E.(T, x0)))
    @show maximum(abs.(cm_1 * Bh - B.(T, x1)))

    print("\nMulti-patch Periodic B-Spline Basis: \n")
    X = PeriodicMultipatchBSpline(domain, k, N, el)
    V0 = X
    V1 = X'

    # calculate the mass matrices
    M0 = mass_matrix(V0)

    M0inv = inv(M0) # using the circulant inverse
    M1 = mass_matrix_d(V1) #special mass matrix

    # calculate the discrete derivative
    G = disc_grad(V0)

    # discrete projections
    Eh = L2_proj(x -> E(0, x), V0)
    Bh = L2_proj_d(x -> B(0, x), V1)

    # simple leap frog scheme
    for i in 1:it
        Eh += Δt/2 * M0inv * (G' * (M1 * Bh))
        Bh -= Δt * (G * Eh)
        Eh += Δt/2 * M0inv * (G' * (M1 * Bh))
    end

    # collocation matrices
    x0, cm_0 = col_mat(V0, 100)
    x1, cm_1 = col_mat_d(V1, 100)

    @show maximum(abs.(cm_0 * Eh - E.(T, x0)))
    @show maximum(abs.(cm_1 * Bh - B.(T, x1)))

    print("\nConga Periodic B-Spline Basis: \n")
    X = SplineFEEC.PeriodicMultipatchBSpline(domain, k, N, el)
    V0 = X
    V1 = X'

    # calculate the mass matrices
    M0 = mass_matrix_d(V0)

    M0inv = 1/Δx^2 * inv(M0) # using the circulant inverse
    M1 = mass_matrix_d(V1) #special mass matrix

    # calculate the discrete derivative
    G = disc_grad(V0)
    P = SplineFEEC.conforming_averaging_projection_matrix(V0)

    # discrete projections
    Eh = 1/Δx * L2_proj_d(x -> E(0, x), V0)
    Bh = L2_proj_d(x -> B(0, x), V1)

    # simple leap frog scheme
    for i in 1:it
        Eh += Δt/2 * M0inv * (P' * G' * (M1 * Bh))
        Bh -= Δt * (G * P * Eh)
        Eh += Δt/2 * M0inv * (P' * G' * (M1 * Bh))
    end

    # collocation matrices
    x0, cm_0 = col_mat(V0, 100)
    x1, cm_1 = col_mat_d(V1, 100)

    @show maximum(abs.(cm_0 * P * Eh - E.(T, x0)))
    @show maximum(abs.(cm_1 * Bh - B.(T, x1)))
end

function test_Issautier_1d(k::Int, N::Int, el::Int)
    # exact solution
    B = (t, x) -> (cos(t) - 1) * π * cos(π * x)
    E = (t, x) -> sin(t) * sin(π * x)
    J = (t, x) -> (cos(t) - 1) * π^2 * sin(π * x) - cos(t) * sin(π * x)
    domain = (0..1)
    it = 1234
    Δt = 0.01
    T = Δt * it
    Δx = (domain.right - domain.left)/el

    @show (k, N, el)
    @show Δt, it

    #exact E L2-norm
    tt = 0
    Enex = zeros(it)
    for i in 1:it
        tt += Δt
        Enex[i] = sqrt(hquadrature(x -> E(tt, x).^2, 0, 1)[1])
    end

    print("\nDirichlet B-Spline Basis: \n")
    X = DirichletBSpline(domain, k, N)
    V0 = X
    V1 = X'

    # calculate the mass matrices
    M0 = mass_matrix(V0)
    M0_old = deepcopy(M0)

    ϵ = 1e-13
    M0[1,1] *= (1 + 1/ϵ)
    M0[end,end] *= (1 + 1/ϵ)

    M0inv = inv(M0) # using the circulant inverse
    M1 = mass_matrix(V1)

    # calculate the discrete derivative
    G = disc_grad(V0)

    # discrete projections
    Bh = L2_proj(x -> B(0, x), V1)
    Eh = L2_proj(x -> E(0, x), V0)
    Jh = L2_prod(x -> J(0,x), V0)

    En1 = zeros(it)
    # simple leap frog scheme
    t = 0
    for i in 1:it
        Eh += Δt/2 * M0inv * (G' * (M1 * Bh) - Jh)

        Bh -= Δt * (G * Eh)
        t += Δt
        Jh =  L2_prod(x -> J(t, x), V0)

        Eh += Δt/2 * M0inv * (G' * (M1 * Bh) - Jh)
        En1[i] = 1/(2*X.N) * sqrt(dot(Eh, M0_old * Eh))
    end

    # collocation matrices
    x0, cm_0 = col_mat(V0, 100)
    x1, cm_1 = col_mat(V1, 100)

    @show maximum(abs.(cm_0 * Eh - E.(T, x0)))
    @show maximum(abs.(cm_1 * Bh - B.(T, x1)))
    @show maximum(abs.(En1 - Enex))

    print("\nMulti-patch Dirichlet B-Spline Basis: \n")
    X = DirichletMultipatchBSpline(domain, k, N, el)
    V0 = X
    V1 = X'

    # calculate the mass matrices
    M0 = mass_matrix(V0)
    M0_old = deepcopy(M0)

    ϵ = 1e-13
    M0[1,1] *= (1 + 1/ϵ)
    M0[end,end] *= (1 + 1/ϵ)

    M0inv = inv(M0) # using the circulant inverse
    M1 = mass_matrix_d(V1) #special mass matrix

    # calculate the discrete derivative
    G = disc_grad(V0)

    # discrete projections
    Eh = L2_proj(x -> E(0, x), V0)
    Bh = L2_proj_d(x -> B(0, x), V1)
    Jh = L2_prod(x -> J(0,x), V0)

    # simple leap frog scheme
    En2 = zeros(it)
    t = 0
    for i in 1:it
        Eh += Δt/2 * M0inv * (G' * (M1 * Bh) - Jh)
        Bh -= Δt * (G * Eh)
        t += Δt
        Jh =  L2_prod(x -> J(t, x), V0)
        Eh += Δt/2 * M0inv * (G' * (M1 * Bh) - Jh)
        En2[i] = 1/(2*X.B.N*el) * sqrt(dot(Eh, M0_old * Eh))
    end

    # collocation matrices
    x0, cm_0 = col_mat(V0, 100)
    x1, cm_1 = col_mat_d(V1, 100)

    @show maximum(abs.(cm_0 * Eh - E.(T, x0)))
    @show maximum(abs.(cm_1 * Bh - B.(T, x1)))
    @show maximum(abs.(En2 - Enex))

    print("\nConga Dirichlet B-Spline Basis: \n")
    X = DirichletMultipatchBSpline(domain, k, N, el)
    V0 = X
    V1 = X'

    # calculate the mass matrices
    M0 = mass_matrix_d(V0)
    M0_old = deepcopy(M0)

    ϵ = 1e-13
    M0[1,1] *= (1 + 1/ϵ)
    M0[end,end] *= (1 + 1/ϵ)

    M0inv = 1/Δx^2 * inv(M0) # using the circulant inverse
    M1 = mass_matrix_d(V1) #special mass matrix

    # calculate the discrete derivative
    G = disc_grad(V0)
    P = conforming_averaging_projection_matrix(V0)

    # discrete projections
    Eh = 1/Δx * L2_proj_d(x -> E(0, x), V0)
    Bh = L2_proj_d(x -> B(0, x), V1)
    #we could also directly project J on V^0_b
    Jh = L2_prod(x -> J(0, x), V0)

    En3 = zeros(it)
    # simple leap frog scheme
    t = 0
    for i in 1:it
        Eh += Δt/2 * M0inv * P' * (G' * (M1 * Bh) -   Jh)
        Bh -= Δt * (G * P * Eh)

        t += Δt
        Jh =  L2_prod(x -> J(t, x), V0)

        Eh += Δt/2 * M0inv * P' * (G' * (M1 * Bh) - Jh)
        En3[i] = 1/(2*X.B.N*el^2) * sqrt(dot(Eh, M0_old * Eh))
    end

    # collocation matrices
    x0, cm_0 = col_mat(V0, 100)
    x1, cm_1 = col_mat_d(V1, 100)

    @show maximum(abs.(cm_0 * P * Eh - E.(T, x0)))
    @show maximum(abs.(cm_1 * Bh - B.(T, x1)))
    @show maximum(abs.(En3 - Enex))

    return it, En1, En2, En3, Enex
end

# -Δu = f
function test_Poisson_2d(k::NTuple{2,Int}, N::NTuple{2,Int})
    u(x, y) = sin(2π * x) * sin(2π * y)
    f(x, y) = -4π^2 * sin(2π * x)* sin(2π * y) - 4π^2 * sin(2π * y) * sin(2π * x)
    domain = (0..1) × (0..1)

    print("Periodic B-Spline Basis: \n")
    basis = PeriodicBSplineTensorProductBasis(domain, k, N)

    Δ = ∇(basis) ⊙ ∇(basis)
    NN = prod(N)
    R = 1/NN .* SplineFEEC.circ_ones(NN)
    P = SplineFEEC.circ_id(NN) - R
    fh = L2_prod(f, basis)

    uh = -(Δ + R) \ (P * fh)

    xy, colmat = col_mat(basis, 100)
    val_u = (i -> u.(xy[i][1], xy[i][2])).(1:length(xy))
    @show maximum(abs.(val_u - colmat * uh))
    #return xy, colmat, fh, Δ, uh

    print("Dirichlet B-Spline Basis: \n")
    basis = DirichletBSplineTensorProductBasis(domain, k, N)

    #homogeneous bc built in
    Δ = ∇(basis) ⊙ ∇(basis)
    fh = L2_prod(f, basis)
    uh = -Δ\ fh

    xy, colmat = col_mat(basis, 100)
    val_u = (i -> u.(xy[i][1], xy[i][2])).(1:length(xy))
    @show maximum(abs.(val_u - colmat * uh))
end

function test_Maxwells_Issautier_2d(k::NTuple{2,Int}, N::NTuple{2,Int}, el::NTuple{2,Int})
    E₂(t, x, y) = sin(t) .* SVector( x * sin(π * y), y *  sin(π * x))
    B₂(t, x, y) = (cos(t) - 1) * π * ( y * cos(x*π) - x * cos(π * y))
    J₂(t, x, y) = (cos(t) - 1) .* SVector(π * cos(π * x) + π^2 * x * sin(π * y), π * cos(π * y) + π^2 * y * sin(π * x)) - cos(t) .* SVector( x * sin(π * y),  y * sin(π * x))
    domain = (0..1) × (0..1)

    @show k
    @show N
    @show el

    it = 1234
    Δt = 0.01
    #it = 1234
    #Δt = 0.00001
    T = Δt * it
    @show Δt, it

    tt = 0
    Enex = zeros(it)
    for i in 1:it
        tt += Δt
        Enex[i] = sqrt(hcubature(x -> dot(E₂(tt, x[1], x[2]), E₂(tt, x[1], x[2])), (0, 0), (1,1))[1])
    end

    print("Dirichlet B-Spline Basis: \n")
    # define bases
    X,Y = DirichletBSplineTensorProductBasis(domain, k, N).B
    V0 = X ⊗ Y
    V1 = SVector(X' ⊗ Y, X ⊗ Y')
    V2 = X' ⊗ Y'

    # calculate mass matrices
    M1 = mass_matrix.(V1) #Vector representing the diagonal of the mass blockmatrix
    M1_old = deepcopy(M1)

    ϵ = 1e-12
    M1[1].B[1,1] *= (1+1/ϵ)
    M1[1].B[end,end] *= (1+1/ϵ)
    M1[2].A[1,1] *= (1+1/ϵ)
    M1[2].A[end,end] *= (1+1/ϵ)

    M1inv = inv.(M1) #using the circulant inverse

    M2 = mass_matrix(V2)

    # calculate discrete curl
    C = disc_curl(V0)

    # or by colocation projection
    (xy1_1, colmat1_1), (xy1_2, colmat1_2) = col_mat_proj.(V1)
    xy2, colmat2 = col_mat_proj(V2)

    cmp_inv1 = inv(colmat1_1)
    cmp_inv2 = inv(colmat1_2)

    val_B = (i -> B₂.(0, xy2[i][1], xy2[i][2])).(1:length(xy2))
    val_E1 = (i -> E₂.(0, xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_E2 = (i -> E₂.(0, xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2));
    val_J1 = (i -> J₂.(0, xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_J2 = (i -> J₂.(0, xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2));

    Eh = SVector(cmp_inv1 * val_E1, cmp_inv2 * val_E2)
    Jh = SVector(cmp_inv1 * val_J1, cmp_inv2 * val_J2)
    Bh = colmat2\val_B;

    #prepare time stepping
    t = 0

    En1 = zeros(it)
    for i in 1:it
        Eh += Δt/2 * (M1inv .* (C' *  SVector{length(Bh)}(M2 * Bh) - Jh))
        Bh -= Δt * (C * Eh)

        t += Δt
        val_J1 = (i -> J₂.(t, xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
        val_J2 = (i -> J₂.(t, xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
        Jh = M1_old .* SVector(cmp_inv1 * val_J1, cmp_inv2 * val_J2)

        Eh += Δt/2 * (M1inv .* (C' *  SVector{length(Bh)}(M2 * Bh) - Jh))
        En1[i] =  1/(2*X.N) * 1/(2*Y.N) * sqrt(dot(Eh, M1_old .* Eh))
    end

    (xy1_1, colmat1_1), (xy1_2, colmat1_2) = col_mat.(V1, 10)
    xy2, colmat2 = col_mat(V2, 10)
    val_B = (i -> B₂.(T, xy2[i][1], xy2[i][2])).(1:length(xy2))
    val_E1 = (i -> E₂.(T, xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_E2 = (i -> E₂.(T, xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))

    @show maximum(abs.(colmat1_1 * Eh[1] - val_E1))
    @show maximum(abs.(colmat1_2 * Eh[2] - val_E2))
    @show maximum(abs.(colmat2 * Bh - val_B))
    @show maximum(abs.(En1 - Enex))

    print("Multipatch BSpline Basis:\n")
    # define bases
    X,Y = DirichletMultipatchTensorProductBSpline(domain, k, N, el).B
    V0 = X ⊗ Y
    V1 = SVector(X' ⊗ Y, X ⊗ Y')
    V2 = X' ⊗ Y'

    ϵ = 1e-12
    M1x = mass_matrix_d(X')
    M1y = mass_matrix(Y)
    M1_1_old = deepcopy(M1x ⊗ M1y)

    M1y[1,1] *= (1+1/ϵ)
    M1y[end,end] *= (1+1/ϵ)
    M1_1 = M1x ⊗ M1y

    M2x = mass_matrix(X)
    M2y = mass_matrix_d(Y')
    M1_2_old = deepcopy(M2x ⊗ M2y)

    M2x[1,1] *= (1+1/ϵ)
    M2x[end,end] *= (1+1/ϵ)
    M1_2 = M2x ⊗ M2y

    M1 = SVector(M1_1, M1_2)
    M1_old = SVector(M1_1_old, M1_2_old)

    M1inv = inv.(M1) #using the circulant inverse

    M2 = mass_matrix_d(V2)

    # calculate discrete curl
    C = disc_diff(curl, V0)

    # colocation projection
    xy1_1, colmat1_1 = SplineFEEC.col_mat_proj_mixed_x(V1[1])
    xy1_2, colmat1_2 = SplineFEEC.col_mat_proj_mixed_y(V1[2])
    xy2, colmat2 = SplineFEEC.col_mat_proj_d(V2)

    cmp_inv1 = inv(colmat1_1)
    cmp_inv2 = inv(colmat1_2)

    val_B = (i -> B₂.(0, xy2[i][1], xy2[i][2])).(1:length(xy2))
    val_E1 = (i -> E₂.(0, xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_E2 = (i -> E₂.(0, xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2));
    val_J1 = (i -> J₂.(0, xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_J2 = (i -> J₂.(0, xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2));

    Eh = SVector(cmp_inv1 * val_E1, cmp_inv2 * val_E2)
    Jh = SVector(cmp_inv1 * val_J1, cmp_inv2 * val_J2)
    Bh = colmat2\val_B;

    #prepare time stepping
    t = 0
    En2 = zeros(it)
    for i in 1:it
        Eh += Δt/2 * (M1inv .* (C' *  SVector{length(Bh)}(M2 * Bh) - Jh))
        Bh -= Δt * (C * Eh)

        t += Δt
        val_J1 = (i -> J₂.(t, xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
        val_J2 = (i -> J₂.(t, xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
        Jh = M1_old .* SVector(cmp_inv1 * val_J1, cmp_inv2 * val_J2)

        Eh += Δt/2 * (M1inv .* (C' *  SVector{length(Bh)}(M2 * Bh) - Jh))
        En2[i] =  1/(2*X.B.N*el[1]) * 1/(2 * Y.B.N*el[2]) * sqrt(dot(Eh, M1_old .* Eh))
    end

    x, cmx = col_mat_d(X', 10)
    y, cmy = col_mat(Y, 10)
    xy1_1 = [(a,b) for a in x for b in y]
    colmat1_1 = cmx ⊗ cmy

    x, cmx = col_mat(X, 10)
    y, cmy = col_mat_d(Y', 10)
    xy1_2 = [(a,b) for a in x for b in y]
    colmat1_2 = cmx ⊗ cmy

    xy2, colmat2 = col_mat_d(V2, 10)

    val_B = (i -> B₂.(T, xy2[i][1], xy2[i][2])).(1:length(xy2))
    val_E1 = (i -> E₂.(T, xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_E2 = (i -> E₂.(T, xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))

    @show maximum(abs.(colmat1_1 * Eh[1] - val_E1))
    @show maximum(abs.(colmat1_2 * Eh[2] - val_E2))
    @show maximum(abs.(colmat2 * Bh - val_B))
    @show maximum(abs.(En2 - Enex))

    print("Conga Multipatch BSpline Basis:\n")
    Δx = (domain.domains[1].right - domain.domains[1].left)/el[1]
    Δy = (domain.domains[2].right - domain.domains[2].left)/el[2]

    # define bases
    X,Y = SplineFEEC.DirichletMultipatchTensorProductBSpline(domain, k, N, el).B
    V0 = X ⊗ Y
    V1 = SVector(X' ⊗ Y, X ⊗ Y')
    V2 = X' ⊗ Y'

    M1 = mass_matrix_d.(V1)
    M1_old = deepcopy(M1)

    ϵ = 1e-12
    M1[1].B[1,1] *= (1+1/ϵ)
    M1[1].B[end,end] *= (1+1/ϵ)
    M1[2].A[1,1] *= (1+1/ϵ)
    M1[2].A[end,end] *= (1+1/ϵ)

    M1inv = inv.(M1)
    M1inv = SVector(1/(Δy)^2, 1/(Δx)^2) .* M1inv

    M2 = mass_matrix_d(V2)

    # calculate discrete curl
    C = disc_diff(curl, V0)

    P = SplineFEEC.conforming_averaging_projection_matrix(V0)
    Px = SplineFEEC.conforming_averaging_projection_matrix_x(V0)
    Py = SplineFEEC.conforming_averaging_projection_matrix_y(V0)

    # or by colocation projection
    (xy1_1, colmat1_1), (xy1_2, colmat1_2) = SplineFEEC.col_mat_proj_d.(V1)
    xy2, colmat2 = SplineFEEC.col_mat_proj_d(V2)

    cmp_inv1 = inv(colmat1_1)
    cmp_inv2 = inv(colmat1_2)

    val_B = (i -> B₂.(0, xy2[i][1], xy2[i][2])).(1:length(xy2))
    val_E1 = (i -> E₂.(0, xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_E2 = (i -> E₂.(0, xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2));
    val_J1 = (i -> J₂.(0, xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_J2 = (i -> J₂.(0, xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2));

    Eh = SVector(Δy *cmp_inv1 * val_E1, Δx *Δy *cmp_inv2 * val_E2)
    Jh = SVector(Δy *cmp_inv1 * val_J1, Δx *cmp_inv2 * val_J2)
    Bh = colmat2\val_B;

    #prepare time stepping
    t = 0
    En3 = zeros(it)
    for i in 1:it
        Eh += Δt/2 *  (M1inv .* (SVector(Py', Px') .* (C' *  SVector{length(Bh)}(M2 * Bh)) - Jh))
        Bh -= Δt *  (C * (SVector(Py , Px) .* Eh))

        t += Δt
        val_J1 = (i -> J₂.(t, xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
        val_J2 = (i -> J₂.(t, xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
        Jh =  M1_old .* SVector(Δy * cmp_inv1 * val_J1, Δx * cmp_inv2 * val_J2)

        Eh += Δt/2 *  (M1inv .* (SVector(Py', Px') .* (C' *  SVector{length(Bh)}(M2 * Bh)) - Jh))
        En3[i] =  1/(2*X.B.N*el[1]^2) * 1/(Y.B.N*el[2]^2) * sqrt(dot(Eh, M1_old .* Eh))
    end

    xy1_1, colmat1_1 = SplineFEEC.col_mat_mixed_x(V1[1], 10)
    xy1_2, colmat1_2 = SplineFEEC.col_mat_mixed_y(V1[2], 10)
    xy2, colmat2 = SplineFEEC.col_mat_d(V2, 10)

    val_B = (i -> B₂.(T, xy2[i][1], xy2[i][2])).(1:length(xy2))
    val_E1 = (i -> E₂.(T, xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_E2 = (i -> E₂.(T, xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))

    @show maximum(abs.(colmat1_1 * Py * Eh[1] - val_E1))
    @show maximum(abs.(colmat1_2 * Px * Eh[2] - val_E2))
    @show maximum(abs.(colmat2 * Bh - val_B))
    @show maximum(abs.(En3 - Enex))
    return it, En1, En2, En3, Enex
end

function test_cold_plasma_xwave_2d(k::NTuple{2,Int}, N::NTuple{2,Int}, el::NTuple{2,Int}, it, Δt)

    #todo fill in for energy conservation
    w_c = 2
    w_p = 1
    w = 0.5

    R = 1 - w_p^2/(w*(w-w_c))
    L = 1 - w_p^2/(w*(w+w_c))
    S = 1/2 * (R + L)

    #Langmuir test
    k2 = 0
    k1 = w * sqrt(R*L/S)
    E(t, x, y) = real.(SVector(1, -im * (R+L)/(R-L)) .* exp(-im*(w*t -k1*x -k2*y)))
    B(t, x, y) = real(-k1/w * im * (R+L)/(R-L) * exp(-im*(w*t -k1*x -k2*y)) )
    J(t, x, y) = real.( w_p^2 .* SVector(im * 1/(w+w_c), 1/(w-w_c) * (R+L)/(R-L)) .* exp(-im*(w*t -k1*x -k2*y)))

    E0(x,y) = E(0, x, y)
    B0(x,y) = B(0, x, y)
    J0(x,y) = J(0, x, y)

    ω_c(x,y) = w_c
    ω_p(x,y) = w_p

    domain = (0..2π/k1) × (0..2π)
    T = Δt * it
    Δx = (domain.domains[1].right - domain.domains[1].left)/el[1]
    Δy = (domain.domains[2].right - domain.domains[2].left)/el[2]

    # Periodic Single Patch
    X,Y = PeriodicBSplineTensorProductBasis(domain, k, N).B
    V0 = X ⊗ Y
    V1 = SVector(X' ⊗ Y, X ⊗ Y')
    V2 = X' ⊗ Y'

    M1 = mass_matrix.(V1) #Vector representing the diagonal of the mass blockmatrix
    M2 = mass_matrix(V2)
    M1inv = inv.(M1)

    C = disc_curl(V0)
    Cc = disc_grad(V2).D[1];
    Cc2 = disc_grad(V1[2]).D[1];

    Eh = SplineFEEC.L2_proj(E0, V1)
    Bh = SplineFEEC.L2_proj(B0, V2)
    Jh = SplineFEEC.L2_proj(J0, V1)

    (xy1_1, colmat1_1), (xy1_2, colmat1_2) = SplineFEEC.col_mat.(V1, 50)
    xy2, colmat2 = SplineFEEC.col_mat(V2, 50);

    t = 0;
    #keep track of energy
    En_E = []
    En_B = []
    En_J = []
    En = []
    Enex_E = []
    Enex_J = []
    Enex_B = []
    Enex = [];

    M_p = w_p .* mass_matrix.(V1)
    M_c = w_c * mass_matrix(X', X ) ⊗ mass_matrix(Y,  Y');

    A1 = [2/Δt .* M1[1]  zero(M1[2]) zero(M1[1]);
    zero(M1[1]) 2/Δt .* M1[2]  -1/2*Cc'*M2;
    zero(M1[1])  1/2*Cc2  2/Δt .* I ]

    A2 = [2/Δt .*M1[1] zero(M1[2]) zero(M1[1]);
    zero(M1[1]) 2/Δt .* M1[2]  1/2*Cc'*M2;
    zero(M1[1]) -1/2*Cc2  2/Δt .* I ];

    B11 = [1/Δt .* M1[1] zero(M1[2]) 1/2 .* M1[1] zero(M1[2]);
            zero(M1[1])  1/Δt .* M1[2] zero(M1[2]) 1/2 .* M1[2];
            -M1[1] * w_p^2/2 zero(M1[2]) 1/Δt*M1[1]  w_c/2*M_c;
            zero(M1[1]) -M1[2]*w_p^2/2  -w_c/2*M_c' 1/Δt*M1[2]]

    B2 = [1/Δt .* M1[1] zero(M1[2]) -1/2 .* M1[1] zero(M1[2]);
            zero(M1[1])  1/Δt .* M1[2] zero(M1[2]) -1/2 .* M1[2];
            M1[1] * w_p^2/2 zero(M1[2]) 1/Δt*M1[1]  -w_c/2*M_c;
            zero(M1[1]) M1[2]*w_p^2/2  w_c/2*M_c' 1/Δt*M1[2]];

    A1inv = inv(A1)
    B1inv = inv(B11);

    Eh1 = Eh[1]
    Eh2 = Eh[2]
    Jh1 = Jh[1]
    Jh2 = Jh[2]
    Bh = Bh

    U1 = [Eh1; Eh2; Bh]
    U2 = [Eh1; Eh2; Jh1; Jh2]
    for i in 1:it
            U1 = [Eh1; Eh2; Bh]
            U1 = A1inv * (A2 * U1)

            Eh1 = U1[1:length(Eh1)]
            Eh2 = U1[length(Eh1)+1:length(Eh1)+length(Eh2)]
            Bh = U1[length(Eh1)+length(Eh2)+1:end]

            U2 = [Eh1; Eh2; Jh1; Jh2]
            U2 = B1inv * (B2 * U2)

            Eh1 = U2[1:length(Eh1)]
            Eh2 = U2[length(Eh1)+1:length(Eh1)+length(Eh2)]
            Jh1 = U2[length(Eh1)+length(Eh2)+1:length(Eh1)+length(Eh2)+length(Jh1)]
            Jh2 = U2[length(Eh1)+length(Eh2)+length(Jh1)+1:end]

            U1 = [Eh1; Eh2; Bh]
            U1 = A1inv * (A2 * U1)

            Eh1 = U1[1:length(Eh1)]
            Eh2 = U1[length(Eh1)+1:length(Eh1)+length(Eh2)]
            Bh = U1[length(Eh1)+length(Eh2)+1:end]

            t += Δt
            push!(En_E, (2π/k1)/(2*X.N) * 2π/(2*Y.N) * sqrt.(Eh1'*M1[1]*Eh1 + Eh2'*M1[2]*Eh2))
            push!(En_J, (2π/k1)/(2*X.N) * 2π/(2*Y.N) * sqrt.(Jh1'*M1[1]*Jh1 + Jh2'*M1[2]*Jh2))
            push!(En_B, (2π/k1)/(2*X.N) * 2π/(2*Y.N) * sqrt.(Bh'*M2*Bh))
            push!(En, 1/2*En_E[end]^2 + 1/2*En_J[end]^2 + 1/2*En_B[end]^2 )

            push!(Enex_E, sqrt(hcubature(x -> dot(E(t, x[1], x[2]), E(t, x[1], x[2])), (0, 0), (2π/k1,2π))[1]))
            push!(Enex_J, sqrt(hcubature(x -> dot(J(t, x[1], x[2]), J(t, x[1], x[2])), (0, 0), (2π/k1,2π))[1]))
            push!(Enex_B, sqrt(hcubature(x -> dot(B(t, x[1], x[2]), B(t, x[1], x[2])), (0, 0), (2π/k1,2π))[1]))
            push!(Enex,  1/2*Enex_E[end]^2 + 1/2*Enex_B[end]^2 + 1/2*Enex_J[end]^2  )
    end

    Eh = SVector(Eh1, Eh2)
    Jh = SVector(Jh1, Jh2)
    Bh = Bh;

    t_n = t
    val_B = (i -> B.(t_n, xy2[i][1], xy2[i][2])).(1:length(xy2))
    val_E1 = (i -> E.(t_n, xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_E2 = (i -> E.(t_n, xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
    val_J1 = (i -> J.(t_n, xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
    val_J2 = (i -> J.(t_n, xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))

    err_per_E = max(maximum(abs.(colmat1_1 * Eh[1] - val_E1)), maximum(abs.(colmat1_2 * Eh[2] - val_E2)))
    err_per_J = max(maximum(abs.(colmat1_1 * Jh[1] - val_J1)), maximum(abs.(colmat1_2 * Jh[2] - val_J2)))
    err_per_B = maximum(abs.(colmat2 * Bh - val_B))
    En_per = En
    Enex_per = Enex

    #Periodic Multipatch
    X,Y = SplineFEEC.PeriodicMultipatchTensorProductBSpline(domain, k, N, el).B
            V0 = X ⊗ Y
            V1 = SVector(X' ⊗ Y, X ⊗ Y')
            V2 = X' ⊗ Y'

            M1x = mass_matrix_d(X')
            M1y = mass_matrix(Y)
            M1_1 = M1x ⊗ M1y

            M2x = mass_matrix(X)
            M2y = mass_matrix_d(Y')
            M1_2 = M2x ⊗ M2y

            M1 = SVector(M1_1, M1_2)
            M1inv = inv.(M1)

            M2 = mass_matrix_d(V2)

            Cc = SplineFEEC.disc_diff(curl, V0).D[1]
            Cc2 = disc_diff(curl, V0).D[1]

            Eh = SVector( M1_1\SplineFEEC.L2_prod_mixed_x((x,y) -> E0(x,y)[1], X' ⊗ Y),
                              M1_2\SplineFEEC.L2_prod_mixed_y((x,y) -> E0(x,y)[2], X ⊗ Y'))
            Jh = SVector( M1_1\SplineFEEC.L2_prod_mixed_x((x,y) -> J0(x,y)[1], X' ⊗ Y),
                              M1_2\SplineFEEC.L2_prod_mixed_y((x,y) -> J0(x,y)[2], X ⊗ Y'))
            Bh = SplineFEEC.L2_proj_d(B0, V2)

            xy0, colmat0 = col_mat(V0, 30)
            xy1_1, colmat1_1 = SplineFEEC.col_mat_mixed_x(V1[1], 30)
            xy1_2, colmat1_2 = SplineFEEC.col_mat_mixed_y(V1[2], 30)
            xy2, colmat2 = col_mat_d(V2, 30)

            t = 0;
            En_E = []
            En_B = []
            En_J = []
            En = []
            #Enex_E = []
            #Enex_J = []
            #Enex_B = []
            #Enex = [];

            M_p = w_p^2 .* M1
            M_c = w_c .* SplineFEEC.mixed_mass_matrix_x(V1[1], V1[2])

            A1 = [2/Δt .* M1[1]  zero(M1[2]) zero(M1[1]);
                    zero(M1[1]) 2/Δt .* M1[2]  -1/2*Cc'*M2;
                    zero(M1[1])  1/2*Cc2  2/Δt .* I ]

            A2 = [2/Δt .*M1[1] zero(M1[2]) zero(M1[1]);
                    zero(M1[1]) 2/Δt .* M1[2]  1/2*Cc'*M2;
                    zero(M1[1]) -1/2*Cc2  2/Δt .* I ];

            Bm1 = [1/Δt .* M1[1] zero(M1[2]) 1/2 .* M1[1] zero(M1[2]);
                    zero(M1[1])  1/Δt .* M1[2] zero(M1[2]) 1/2 .* M1[2];
                    -M1[1] * w_p^2/2 zero(M1[2]) 1/Δt*M1[1]  w_c/2*M_c;
                    zero(M1[1]) -M1[2]*w_p^2/2  -w_c/2*M_c' 1/Δt*M1[2]]

            Bm2 = [1/Δt .* M1[1] zero(M1[2]) -1/2 .* M1[1] zero(M1[2]);
                    zero(M1[1])  1/Δt .* M1[2] zero(M1[2]) -1/2 .* M1[2];
                    M1[1] * w_p^2/2 zero(M1[2]) 1/Δt*M1[1]  -w_c/2*M_c;
                    zero(M1[1]) M1[2]*w_p^2/2  w_c/2*M_c' 1/Δt*M1[2]];

            A1inv = inv(A1)
            Bm1inv = inv(Bm1);

            Eh1 = Eh[1]
            Eh2 = Eh[2]
            Jh1 = Jh[1]
            Jh2 = Jh[2]
            Bh = Bh

            U1 = [Eh1; Eh2; Bh]
            U2 = [Eh1; Eh2; Jh1; Jh2]
            for i in 1:it
                    U1 = [Eh1; Eh2; Bh]
                    U1 = A1inv*(A2 * U1)

                    Eh1 = U1[1:length(Eh1)]
                    Eh2 = U1[length(Eh1)+1:length(Eh1)+length(Eh2)]
                    Bh = U1[length(Eh1)+length(Eh2)+1:end]

                    U2 = [Eh1; Eh2; Jh1; Jh2]
                    U2 = Bm1inv*(Bm2 * U2)

                    Eh1 = U2[1:length(Eh1)]
                    Eh2 = U2[length(Eh1)+1:length(Eh1)+length(Eh2)]
                    Jh1 = U2[length(Eh1)+length(Eh2)+1:length(Eh1)+length(Eh2)+length(Jh1)]
                    Jh2 = U2[length(Eh1)+length(Eh2)+length(Jh1)+1:end]

                    U1 = [Eh1; Eh2; Bh]
                    U1 = A1inv*(A2 * U1)

                    Eh1 = U1[1:length(Eh1)]
                    Eh2 = U1[length(Eh1)+1:length(Eh1)+length(Eh2)]
                    Bh = U1[length(Eh1)+length(Eh2)+1:end]

                    t += Δt
                    push!(En_E, (2π/k1)/(2*X.B.N*el[1]) * 2π/(2*Y.B.N*el[2]) * sqrt.(Eh1'*M1[1]*Eh1 + Eh2'*M1[2]*Eh2))
                    push!(En_J, (2π/k1)/(2*X.B.N*el[1]) * 2π/(2*Y.B.N*el[2]) * sqrt.(Jh1'*M1[1]*Jh1 + Jh2'*M1[2]*Jh2))
                    push!(En_B, (2π/k1)/(2*X.B.N*el[1]) * 2π/(2*Y.B.N*el[2]) * sqrt.(Bh'*M2*Bh))
                    push!(En, 1/2*En_E[end]^2 + 1/2*En_J[end]^2 + 1/2*En_B[end]^2 )

                    #push!(Enex_E, sqrt(hcubature(x -> dot(E(tt, x[1], x[2]), E(tt, x[1], x[2])), (0, 0), (2π/k1,2π))[1]))
                    #push!(Enex_J, sqrt(hcubature(x -> dot(J(tt, x[1], x[2]), J(tt, x[1], x[2])), (0, 0), (2π/k1,2π))[1]))
                    #push!(Enex_B, sqrt(hcubature(x -> dot(B(tt, x[1], x[2]), B(tt, x[1], x[2])), (0, 0), (2π/k1,2π))[1]))
                    #push!(Enex,  1/2*Enex_E[end]^2 + 1/2*Enex_B[end]^2 + 1/2*Enex_J[end]^2  )
            end

            Eh = SVector(Eh1, Eh2)
            Jh = SVector(Jh1, Jh2)
            Bh = Bh;

            t_n = t
            val_B = (i -> B.(t_n, xy2[i][1], xy2[i][2])).(1:length(xy2))
            val_E1 = (i -> E.(t_n, xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
            val_E2 = (i -> E.(t_n, xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
            val_J1 = (i -> J.(t_n, xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
            val_J2 = (i -> J.(t_n, xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))

            En_mp_per = En
            err_mp_E = max(maximum(abs.(colmat1_1 * Eh[1] - val_E1)), maximum(abs.(colmat1_2 * Eh[2] - val_E2)))
            err_mp_J = max(maximum(abs.(colmat1_1 * Jh[1] - val_J1)), maximum(abs.(colmat1_2 * Jh[2] - val_J2)))
            err_mp_B = maximum(abs.(colmat2 * Bh - val_B))

            #Cong Multipatch
            X,Y = SplineFEEC.PeriodicMultipatchTensorProductBSpline(domain, k, N, el).B
            V0 = X ⊗ Y
            V1 = SVector(X' ⊗ Y, X ⊗ Y')
            V2 = X' ⊗ Y'

            M1 = mass_matrix_d.(V1)
            M1inv = SVector(1/(Δy)^2, 1/(Δx)^2) .* inv.(M1)
            M1 = SVector((Δy)^2, (Δx)^2) .* M1
            M2 = mass_matrix_d(V2)

            #C = SplineFEEC.disc_diff(curl, V0)

            P = SplineFEEC.conforming_averaging_projection_matrix(V0)
            Px = SplineFEEC.conforming_averaging_projection_matrix_x(V0)
            Py = SplineFEEC.conforming_averaging_projection_matrix_y(V0)

            Cc = SplineFEEC.disc_diff(curl, V0).D[1] * Px
            Cc2 = disc_diff(curl, V0).D[1] * Px

            Eh = SVector(1/Δy, 1/Δx) .* L2_proj_d(E0, V1)
            Jh = SVector(1/Δy, 1/Δx) .* L2_proj_d(J0, V1)
            #SVector( M1_1\SplineFEEC.L2_prod_mixed_x((x,y) -> E0(x,y)[1], X' ⊗ Y),
            #                  M1_2\SplineFEEC.L2_prod_mixed_y((x,y) -> E0(x,y)[2], X ⊗ Y'))
            #Jh = SVector( M1_1\SplineFEEC.L2_prod_mixed_x((x,y) -> J0(x,y)[1], X' ⊗ Y),
            #                  M1_2\SplineFEEC.L2_prod_mixed_y((x,y) -> J0(x,y)[2], X ⊗ Y'))
            Bh = SplineFEEC.L2_proj_d(B0, V2)

            xy0, colmat0 = col_mat(V0, 10)
            xy1_1, colmat1_1 = SplineFEEC.col_mat_mixed_x(V1[1], 10)
            xy1_2, colmat1_2 = SplineFEEC.col_mat_mixed_y(V1[2], 10)
            xy2, colmat2 = col_mat_d(V2, 10)

            t = 0;
            En_E = []
            En_B = []
            En_J = []
            En = []
            #Enex_E = []
            #Enex_J = []
            #Enex_B = []
            #Enex = [];

            M_p = w_p^2 .* M1
            M_c = w_c .* SplineFEEC.mixed_mass_matrix_x_d(V1[1], V1[2]);

            A1 = [2/Δt .* M1[1]  zero(M1[2]) zero(Cc'*M2);
                    zero(M1[1]) 2/Δt .* M1[2]  -1/2*Cc'*M2;
                    zero(Cc2)  1/2*Cc2  2/Δt .* I ]

            A2 = [2/Δt .*M1[1] zero(M1[2]) zero(Cc'*M2);
                    zero(M1[1]) 2/Δt .* M1[2]  1/2*Cc'*M2;
                    zero(Cc2) -1/2*Cc2  2/Δt .* I ];

            Bm1 = [1/Δt .* M1[1] zero(M1[2]) 1/2 .* M1[1] zero(M1[2]);
                    zero(M1[1])  1/Δt .* M1[2] zero(M1[1]) 1/2 .* M1[2];
                    -M1[1] * w_p^2/2 zero(M1[2]) 1/Δt*M1[1]  w_c/2*M_c;
                    zero(M1[1]) -M1[2]*w_p^2/2  -w_c/2*M_c' 1/Δt*M1[2]]

            Bm2 = [1/Δt .* M1[1] zero(M1[2]) -1/2 .* M1[1] zero(M1[2]);
                    zero(M1[1])  1/Δt .* M1[2] zero(M1[1]) -1/2 .* M1[2];
                    M1[1] * w_p^2/2 zero(M1[2]) 1/Δt*M1[1]  -w_c/2*M_c;
                    zero(M1[1]) M1[2]*w_p^2/2  w_c/2*M_c' 1/Δt*M1[2]];

            A1inv = inv(A1)
            Bm1inv = inv(Bm1);

            Eh1 = Eh[1]
            Eh2 = Eh[2]
            Jh1 = Jh[1]
            Jh2 = Jh[2]
            Bh = Bh

            U1 = [Eh1; Eh2; Bh]
            U2 = [Eh1; Eh2; Jh1; Jh2]
            for i in 1:it
                    U1 = [Eh1; Eh2; Bh]
                    U1 = A1inv*(A2 * U1)

                    Eh1 = U1[1:length(Eh1)]
                    Eh2 = U1[length(Eh1)+1:length(Eh1)+length(Eh2)]
                    Bh = U1[length(Eh1)+length(Eh2)+1:end]

                    U2 = [Eh1; Eh2; Jh1; Jh2]
                    U2 = Bm1inv*(Bm2 * U2)

                    Eh1 = U2[1:length(Eh1)]
                    Eh2 = U2[length(Eh1)+1:length(Eh1)+length(Eh2)]
                    Jh1 = U2[length(Eh1)+length(Eh2)+1:length(Eh1)+length(Eh2)+length(Jh1)]
                    Jh2 = U2[length(Eh1)+length(Eh2)+length(Jh1)+1:end]

                    U1 = [Eh1; Eh2; Bh]
                    U1 = A1inv*(A2 * U1)

                    Eh1 = U1[1:length(Eh1)]
                    Eh2 = U1[length(Eh1)+1:length(Eh1)+length(Eh2)]
                    Bh = U1[length(Eh1)+length(Eh2)+1:end]

                    t += Δt
                    push!(En_E, (2π/k1)/(2*X.B.N*el[1]) * 2π/(2*Y.B.N*el[2]) * sqrt.(Eh1'*M1[1]*Eh1 + Eh2'*M1[2]*Eh2))
                    push!(En_J, (2π/k1)/(2*X.B.N*el[1]) * 2π/(2*Y.B.N*el[2]) * sqrt.(Jh1'*M1[1]*Jh1 + Jh2'*M1[2]*Jh2))
                    push!(En_B, (2π/k1)/(2*X.B.N*el[1]) * 2π/(2*Y.B.N*el[2]) * sqrt.(Bh'*M2*Bh))
                    push!(En, 1/2*En_E[end]^2 + 1/2*En_J[end]^2 + 1/2*En_B[end]^2 )

                   # push!(Enex_E, sqrt(hcubature(x -> dot(E(tt, x[1], x[2]), E(tt, x[1], x[2])), (0, 0), (2π/k1,2π))[1]))
                   # push!(Enex_J, sqrt(hcubature(x -> dot(J(tt, x[1], x[2]), J(tt, x[1], x[2])), (0, 0), (2π/k1,2π))[1]))
                   # push!(Enex_B, sqrt(hcubature(x -> dot(B(tt, x[1], x[2]), B(tt, x[1], x[2])), (0, 0), (2π/k1,2π))[1]))
                   # push!(Enex,  1/2*Enex_E[end]^2 + 1/2*Enex_B[end]^2 + 1/2*Enex_J[end]^2  )
            end

            Eh = SVector(Eh1, Eh2)
            Jh = SVector(Jh1, Jh2)
            Bh = Bh;

            t_n = t
            val_B = (i -> B.(t_n, xy2[i][1], xy2[i][2])).(1:length(xy2))
            val_E1 = (i -> E.(t_n, xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
            val_E2 = (i -> E.(t_n, xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))
            val_J1 = (i -> J.(t_n, xy1_1[i][1], xy1_1[i][2])[1]).(1:length(xy1_1))
            val_J2 = (i -> J.(t_n, xy1_2[i][1], xy1_2[i][2])[2]).(1:length(xy1_2))

            En_conga_per = En
            err_conga_E = max(maximum(abs.(colmat1_1 * Py * Eh[1] - val_E1)), maximum(abs.(colmat1_2 * Px * Eh[2] - val_E2)))
            err_conga_J = max(maximum(abs.(colmat1_1 * Py * Jh[1] - val_J1)), maximum(abs.(colmat1_2 * Px * Jh[2] - val_J2)))
            err_conga_B = maximum(abs.(colmat2 * Bh - val_B))

    return it, En_per, En_mp_per, En_conga_per, Enex_per
end


text:
# __BEGIN_LICENSE__
#
# ThreeDeconv.jl
# Author: Hayato Ikoma (h9koma@stanford.edu)
#
# Copyright (c) 2018, Stanford University
#
# All rights reserved.
#
# Redistribution and use in source and binary forms for academic and other
# non-commercial purposes with or without modification, are permitted provided
# that the following conditions are met:
#
# * Redistributions of source code, including modified source code, must retain
#   the above copyright notice, this list of conditions and the following
#   disclaimer.
#
# * Redistributions in binary form or a modified form of the source code must
#   reproduce the above copyright notice, this list of conditions and the
#   following disclaimer in the documentation and/or other materials provided with
#   the distribution.
#
# * Neither the name of The Leland Stanford Junior University, any of its
#   trademarks, the names of its employees, nor contributors to the source code
#   may be used to endorse or promote products derived from this software without
#   specific prior written permission.
#
# * Where a modified version of the source code is redistributed publicly in
#   source or binary forms, the modified source code must be published in a freely
#   accessible manner, or otherwise redistributed at no charge to anyone
#   requesting a copy of the modified source code, subject to the same terms as
#   this agreement.
#
# THIS SOFTWARE IS PROVIDED BY THE TRUSTEES OF THE LELAND STANFORD JUNIOR
# UNIVERSITY "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE LELAND STANFORD JUNIOR
# UNIVERSITY OR ITS TRUSTEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
# OF SUCH DAMAGE.
#
# __END_LICENSE__

module ThreeDeconv

import Base.GC.gc
using LinearAlgebra, FFTW, CUDA
FFTW.set_num_threads(4)

if CUDA.functional()
    to_gpu_or_not_to_gpu(x::AbstractArray) = CuArray(x)
else
    to_gpu_or_not_to_gpu(x::AbstractArray) = x
end

include("psf/psf.jl")
include("util/linearoperator.jl")
include("util/fft.jl")
include("util/util.jl")
include("noiseestimation/noiseestimation.jl")
include("deconvolution/deconvolve.jl")

end # module


text:
function toPointAxis(aout,ain,loninds,latinds)
  iout = 1
  for (ilon,ilat) in zip(loninds,latinds)
    aout[iout]=ain[ilon,ilat]
    iout+=1
  end
end

"""
    extractLonLats(c,pl::Matrix)

Extracts a list of longitude/latitude coordinates from a data cube. The coordinates
are specified through the matrix `pl` where `size(pl)==(N,2)` and N is the number
of extracted coordinates. Returns a data cube without `LonAxis` and `LatAxis` but with a
`SpatialPointAxis` containing the input locations.
"""
function extractLonLats(c,pl::Matrix;kwargs...)
  size(pl,2)==2 || error("Coordinate list must have exactly 2 columns")
  axlist=caxes(c)
  ilon=findAxis("Lon",axlist)
  ilat=findAxis("Lat",axlist)
  ilon>0 || error("Input cube must contain a LonAxis")
  ilat>0 || error("input cube must contain a LatAxis")
  lonax=axlist[ilon]
  latax=axlist[ilat]
  pointax = CategoricalAxis("SpatialPoint", [(pl[i,1],pl[i,2]) for i in 1:size(pl,1)])
  loninds = map(ll->axVal2Index(lonax,ll[1]),pointax.values)
  latinds = map(ll->axVal2Index(latax,ll[2]),pointax.values)
  incubes=InDims("Lon","Lat")
  outcubes=OutDims(pointax)
  y=mapCube(toPointAxis,c,loninds,latinds;indims=incubes,outdims=outcubes,max_cache=1e8,kwargs...)
end


text:
using ArbRadixFloatingPoints
using Test

let x = 2im, z = convert(ArbRadixFloat{x,2}, x)
    @test z.exponent == 1
    @test z.significand == [1, 0]
    @test float(z) == x
end

let x = 2im, z = convert(ArbRadixFloat{-0.5im,2}, x)
    @test z.exponent == 0
    @test z.significand == [0, 1]
    @test float(z) == x
end


text:
function pmap_fwd(B::LLVM.API.LLVMBuilderRef, OrigCI::LLVM.API.LLVMValueRef, gutils::API.EnzymeGradientUtilsRef, normalR::Ptr{LLVM.API.LLVMValueRef}, shadowR::Ptr{LLVM.API.LLVMValueRef})::Cvoid
    orig = LLVM.Instruction(OrigCI)
    mod = LLVM.parent(LLVM.parent(LLVM.parent(orig)))
    ctx = LLVM.context(orig)
    B = LLVM.Builder(B)
    emit_error("fast pfor not implemented");
    return nothing
end

@generated function callfn(ptr, args...)
    ctx = LLVM.Context()
    Pvoid = convert(LLVMType, Ptr{Cvoid}; ctx)
    llvmtys = LLVMType[convert(LLVMType, Int; ctx)]
    realargs = Any[:ptr]
    realtypes = [ptr]
    for (i, a) in enumerate(args)
        if GPUCompiler.isghosttype(a) || Core.Compiler.isconstType(a)
            continue
        end
        push!(llvmtys, convert(LLVMType, a; ctx, allow_boxed=true))
        push!(realargs, :(args[$i]))
        push!(realtypes, a)
    end
    llvm_f, _ = LLVM.Interop.create_function(Pvoid, llvmtys)

	LLVM.Builder(ctx) do builder
		entry = BasicBlock(llvm_f, "entry"; ctx)
		position!(builder, entry)
		params = collect(LLVM.Value, parameters(llvm_f))
		lfn = @inbounds params[1]
		nparams = LLVM.Value[]

        for (parm, source_typ) in zip(params[2:end], realtypes[2:end])
            push!(nparams, parm)
            # codegen_typ = llvmtype(parm)
            # @show parm, source_typ
            # if codegen_typ isa LLVM.PointerType && !issized(eltype(codegen_typ))
            #     push!(nparams, parm)
            #     #push!(args, (cc=GPUCompiler.MUT_REF, typ=source_typ,
            #     #             codegen=(typ=codegen_typ, i=codegen_i)))
            # elseif codegen_typ isa LLVM.PointerType && issized(eltype(codegen_typ)) &&
            #        !(source_typ <: Ptr) && !(source_typ <: Core.LLVMPtr)
            #     if !GPUCompiler.deserves_argbox(source_typ)
            #       push!(nparams, load!(builder, parm))
            #     else
            #       push!(nparams, parm)
            #     end
            #     # push!(args, (cc=GPUCompiler.BITS_REF, typ=source_typ,
            #     #             codegen=(typ=codegen_typ, i=codegen_i)))
            # else
            #     push!(nparams, parm)
            #     # push!(args, (cc=GPUCompiler.BITS_VALUE, typ=source_typ,
            #     #              codegen=(typ=codegen_typ, i=codegen_i)))
            # end
        end

		lfn = inttoptr!(builder, lfn, LLVM.PointerType(LLVM.FunctionType(Pvoid, [llvmtype(x) for x in nparams])))
		res = call!(builder, lfn, nparams)
        ret!(builder, res)
	end

    mod = LLVM.parent(llvm_f)

	ir = string(mod)
	fn = LLVM.name(llvm_f)

    quote
        Base.@_inline_meta
		Base.llvmcall(($ir, $fn), Ptr{Cvoid}, $(Tuple{realtypes...}), $(realargs...))
	end
end

function runtime_pmap_augfwd(count, forward, args...)::Ptr{Ptr{Cvoid}}
    tapes = Base.unsafe_convert(Ptr{Ptr{Cvoid}}, Libc.malloc(sizeof(Ptr{Cvoid})*count))
    function fwd(idx, tapes, f_func, fargs...)
        st = callfn(f_func, idx, fargs...)
        Base.unsafe_store!(tapes, st, idx)
    end
    Enzyme.pmap(count, fwd, tapes, forward, args...)
    return tapes
end

function runtime_pmap_rev(count, adjoint, tapes, args...)
	function adj(idx, tapes, r_func, rargs...)
        st = unsafe_load(tapes, idx)
        callfn(r_func, idx, rargs..., st)
        nothing
	end
	Enzyme.pmap(count, adj, tapes, adjoint, args...)
    Libc.free(tapes)
    return nothing
end

function julia_activity(source_types, FTs, ops, gutils, tape::Bool)
    source_types = source_types[2:end]
    # count, funcT, funcT
    args = Type[source_types[1]]
    dup_args = Type[Const{source_types[1]}]
    if tape
      push!(args, Ptr{Ptr{Cvoid}})
      push!(dup_args, Ptr{Ptr{Cvoid}})
    end
    for T in FTs
      push!(args, T)
      push!(dup_args, T)
    end
    codegen_i = 2

    for source_typ in source_types[3:end]
        if isghosttype(source_typ) || Core.Compiler.isconstType(source_typ)
            continue
        end

        codegen_typ = llvmtype(ops[codegen_i])
        if codegen_typ isa LLVM.PointerType && !issized(eltype(codegen_typ))
            push!(args, source_typ)
            if API.EnzymeGradientUtilsIsConstantValue(gutils, ops[codegen_i]) == 0
              push!(args, source_typ)
              push!(dup_args, Duplicated{source_typ})
            else
              push!(dup_args, Const{source_typ})
            end
            #push!(args, (cc=GPUCompiler.MUT_REF, typ=source_typ,
            #             codegen=(typ=codegen_typ, i=codegen_i)))
        elseif codegen_typ isa LLVM.PointerType && issized(eltype(codegen_typ)) &&
               !(source_typ <: Ptr) && !(source_typ <: Core.LLVMPtr)
            push!(args, source_typ)
            if API.EnzymeGradientUtilsIsConstantValue(gutils, ops[codegen_i]) == 0
              push!(args, source_typ)
              push!(dup_args, Duplicated{source_typ})
            else
              push!(dup_args, Const{source_typ})
            end
            # push!(args, (cc=GPUCompiler.BITS_REF, typ=source_typ,
            #             codegen=(typ=codegen_typ, i=codegen_i)))
        else
            push!(args, source_typ)
            if API.EnzymeGradientUtilsIsConstantValue(gutils, ops[codegen_i]) == 0
              push!(args, source_typ)
              push!(dup_args, Duplicated{source_typ})
            else
              push!(dup_args, Const{source_typ})
            end
            # push!(args, (cc=GPUCompiler.BITS_VALUE, typ=source_typ,
            #              codegen=(typ=codegen_typ, i=codegen_i)))
        end
        codegen_i += 1
    end
    return args, dup_args
end

function commonInnerCompile(runtime_fn, B, orig, gutils, tape)
    mod = LLVM.parent(LLVM.parent(LLVM.parent(orig)))
    ctx = LLVM.context(orig)

    llvmfn = LLVM.called_value(orig)
    mi = nothing
    adjointnm = nothing
    forwardnm = nothing
    for fattr in collect(function_attributes(llvmfn))
        if isa(fattr, LLVM.StringAttribute)
            if kind(fattr) == "enzymejl_mi"
                ptr = reinterpret(Ptr{Cvoid}, parse(Int, LLVM.value(fattr)))
                mi = Base.unsafe_pointer_to_objref(ptr)
                break
            end
            if kind(fattr) == "enzymejl_forward"
                forwardnm = value(fattr)
            end
            if kind(fattr) == "enzymejl_adjoint"
                adjointnm = value(fattr)
            end
        end
    end

    countT = mi.specTypes.parameters[2]
	funcT = mi.specTypes.parameters[3]

    ops = collect(operands(orig))[1:end-1]

    if forwardnm === nothing
        _, dup = julia_activity(mi.specTypes.parameters, [], ops, gutils, #=tape=#false)
        e_tt = Tuple{dup...}
        @static if VERSION >= v"1.8"
          RT = Core.Compiler.return_type(Tuple{funcT, map(eltype, dup)...})
        else
          RT = Core.Compiler.return_type(Core.Compiler.singleton_type(funcT), Tuple{map(eltype, dup)...})
        end
        eprimal, eadjoint = fspec(Core.Compiler.singleton_type(funcT), e_tt)

        etarget = Compiler.EnzymeTarget()
        eparams = Compiler.EnzymeCompilerParams(eadjoint, API.DEM_ReverseModePrimal, Const{RT}, true, #=shadowfunc=#false, #=abiwrap=#false)
        ejob    = Compiler.CompilerJob(etarget, eprimal, eparams)

        cmod, adjointnm, forwardnm = _thunk(ejob)
        LLVM.link!(mod, cmod)
        attributes = function_attributes(llvmfn)
        push!(attributes, StringAttribute("enzymejl_forward", forwardnm; ctx))
        push!(attributes, StringAttribute("enzymejl_adjoint", adjointnm; ctx))
        attributes = function_attributes(llvmfn)
        push!(function_attributes(functions(mod)[forwardnm]), EnumAttribute("alwaysinline"; ctx))
        push!(function_attributes(functions(mod)[adjointnm]), EnumAttribute("alwaysinline"; ctx))
    end

    splat, _ = julia_activity(mi.specTypes.parameters, tape === nothing ? [Int, funcT, funcT] : [Int, Ptr{Ptr{Cvoid}}, funcT, funcT], ops, gutils, #=tape=#false)
    # splat[1] = eltype(splat[1])
    tt = Tuple{splat...}

    funcspec = FunctionSpec(runtime_fn, tt, #=kernel=# false, #=name=# nothing)

    # 3) Use the MI to create the correct augmented fwd/reverse
    # TODO:
    #  - GPU support
    #  - When OrcV2 only use a MaterializationUnit to avoid mutation of the module here

    target = GPUCompiler.NativeCompilerTarget()
    params = Compiler.PrimalCompilerParams()
    job    = CompilerJob(target, funcspec, params)

    otherMod, meta = GPUCompiler.codegen(:llvm, job; optimize=false, validate=false, ctx)
    entry = name(meta.entry)
    optimize!(otherMod, JIT.get_tm())

    # 4) Link the corresponding module
    LLVM.link!(mod, otherMod)

    # 5) Call the function
    entry = functions(mod)[entry]

    B = LLVM.Builder(B)

    T_int64 = LLVM.Int64Type(ctx)
    T_jlvalue = LLVM.StructType(LLVMType[]; ctx)
    T_prjlvalue = LLVM.PointerType(T_jlvalue, #= AddressSpace::Tracked =# 10)
    T_pprjlvalue = LLVM.PointerType(T_prjlvalue)

    # count
	vals = LLVM.Value[LLVM.Value(API.EnzymeGradientUtilsNewFromOriginal(gutils, ops[1]))]

    # function
    run_fn = functions(mod)[tape === nothing ? forwardnm : adjointnm]
    push!(vals, ptrtoint!(B, run_fn, llvmtype(LLVM.ConstantInt(Int(0); ctx))))

    if tape !== nothing
        push!(vals, tape)
    end

    EB = LLVM.Builder(ctx)
    position!(EB, LLVM.BasicBlock(API.EnzymeGradientUtilsAllocationBlock(gutils)))

    i = 2
    for source_typ in mi.specTypes.parameters[3:end]
        if isghosttype(source_typ) || Core.Compiler.isconstType(source_typ)
            continue
        end

        primal = LLVM.Value(API.EnzymeGradientUtilsNewFromOriginal(gutils, ops[i]))
        shadow = if API.EnzymeGradientUtilsIsConstantValue(gutils, ops[i]) == 0
          LLVM.Value(API.EnzymeGradientUtilsInvertPointer(gutils, ops[i], B))
        else
          nothing
        end

        codegen_typ = llvmtype(parameters(entry)[i+1 + (tape !== nothing)])
        if codegen_typ isa LLVM.PointerType && !issized(eltype(codegen_typ))
            push!(vals, primal)
            if shadow !== nothing
              push!(vals, shadow)
            end
            #push!(args, (cc=GPUCompiler.MUT_REF, typ=source_typ,
            #             codegen=(typ=codegen_typ, i=codegen_i)))
        elseif codegen_typ isa LLVM.PointerType && issized(eltype(codegen_typ)) &&
               !(source_typ <: Ptr) && !(source_typ <: Core.LLVMPtr)
            if !GPUCompiler.deserves_argbox(source_typ)
              primA = alloca!(EB, llvmtype(primal))
              store!(B, primal, primA)
              primal = addrspacecast!(B, primA, codegen_typ)
            end
            push!(vals, primal)
            if shadow !== nothing
              if !GPUCompiler.deserves_argbox(source_typ)
                shadowA = alloca!(EB, llvmtype(shadow))
                store!(B, shadow, shadowA)
                shadow = addrspacecast!(B, shadowA, codegen_typ)
              end
              push!(vals, shadow)
            end
            # push!(args, (cc=GPUCompiler.BITS_REF, typ=source_typ,
            #             codegen=(typ=codegen_typ, i=codegen_i)))
        else
            push!(vals, primal)
            if shadow !== nothing
              push!(vals, shadow)
            end
            # push!(args, (cc=GPUCompiler.BITS_VALUE, typ=source_typ,
            #              codegen=(typ=codegen_typ, i=codegen_i)))
        end
        i += 1
    end

    return LLVM.call!(B, entry, vals)
end

function pmap_augfwd(B::LLVM.API.LLVMBuilderRef, OrigCI::LLVM.API.LLVMValueRef, gutils::API.EnzymeGradientUtilsRef, normalR::Ptr{LLVM.API.LLVMValueRef}, shadowR::Ptr{LLVM.API.LLVMValueRef}, tapeR::Ptr{LLVM.API.LLVMValueRef})::Cvoid
    orig = LLVM.Instruction(OrigCI)
    normal = (unsafe_load(normalR) != C_NULL) ? LLVM.Instruction(unsafe_load(normalR)) : nothing
    shadow = (unsafe_load(shadowR) != C_NULL) ? LLVM.Instruction(unsafe_load(shadowR)) : nothing

    @warn "active variables passed by value to jl_pmap not yet supported"
    tape = commonInnerCompile(runtime_pmap_augfwd, B, orig, gutils, nothing)

    # Delete the primal code
    if normal !== nothing
        unsafe_store!(normalR, C_NULL)
    else
        LLVM.API.LLVMInstructionEraseFromParent(LLVM.Instruction(API.EnzymeGradientUtilsNewFromOriginal(gutils, orig)))
    end

    unsafe_store!(tapeR, tape.ref)

    return nothing
end

function pmap_rev(B::LLVM.API.LLVMBuilderRef, OrigCI::LLVM.API.LLVMValueRef, gutils::API.EnzymeGradientUtilsRef, tape::LLVM.API.LLVMValueRef)::Cvoid
    orig = LLVM.Instruction(OrigCI)
    commonInnerCompile(runtime_pmap_rev, B, orig, gutils, LLVM.Value(tape))
    return nothing
end


text:
module Filters

# builtin filters
export Deflate, Shuffle, Fletcher32, Szip, NBit, ScaleOffset

import ..HDF5: Properties, h5doc, API

"""
    Filter

Abstract type to describe HDF5 Filters.
See the Extended Help for information on implementing a new filter.

# Extended Help

## Filter interface

The Filter interface is implemented upon the Filter subtype.

See API.h5z_register for details.

### Required Methods to Implement
* `filterid` - registered filter ID
* `filter_func` - implement the actual filter

### Optional Methods to Implement
* `filtername` - defaults to "Unnamed Filter"
* `encoder_present` - defaults to true
* `decoder_present` - defaults to true
* `can_apply_func` - defaults to nothing
* `set_local_func` - defaults to nothing

### Advanced Methods to Implement
* `can_apply_cfunc` - Defaults to wrapping @cfunction around the result of `can_apply_func`
* `set_local_cfunc` - Defaults to wrapping @cfunction around the result of `set_local_func`
* `filter_cfunc` - Defaults to wrapping @cfunction around the result of `filter_func`
* `register_filter` - Defaults to using the above functions to register the filter

Implement the Advanced Methods to avoid @cfunction from generating a runtime closure which may not work on all systems.
"""
abstract type Filter end

"""
    FILTERS

Maps filter id to filter type.
"""
const FILTERS = Dict{API.H5Z_filter_t, Type{<: Filter}}()

"""
    filterid(F) where {F <: Filter}

The internal filter id of a filter of type `F`.
"""
filterid

"""
    encoder_present(::Type{F}) where {F<:Filter}

Can the filter have an encode or compress the data?
Defaults to true.
Returns a Bool. See `API.h5z_register`.
"""
encoder_present(::Type{F}) where {F<:Filter} = true

"""
    decoder_present(::Type{F}) where {F<:Filter}

Can the filter decode or decompress the data?
Defaults to true.
Returns a Bool.
See `API.h5z_register`
"""
decoder_present(::Type{F}) where {F<:Filter} = true

"""
    filtername(::Type{F}) where {F<:Filter}

What is the name of a filter?
Defaults to "Unnamed Filter"
Returns a String describing the filter. See `API.h5z_register`
"""
filtername(::Type{F}) where {F<:Filter} = "Unnamed Filter"

"""
    can_apply_func(::Type{F}) where {F<:Filter}

Return a function indicating whether the filter can be applied or `nothing` if no function exists.
The function signature is `func(dcpl_id::API.hid_t, type_id::API.hid_t, space_id::API.hid_t)`.
See `API.h5z_register`
"""
can_apply_func(::Type{F}) where {F<:Filter} = nothing

"""
    can_apply_cfunc(::Type{F}) where {F<:Filter}

Return a C function pointer for the can apply function.
By default, this will return the result of using `@cfunction` on the function
specified by `can_apply_func(F)` or `C_NULL` if `nothing`.

Overriding this will allow `@cfunction` to return a `Ptr{Nothing}` rather
than a `CFunction`` closure which may not work on all systems.
"""
function can_apply_cfunc(::Type{F}) where {F<:Filter}
    func = can_apply_func(F)
    if func === nothing
        return C_NULL
    else
        return @cfunction($func, API.herr_t, (API.hid_t,API.hid_t,API.hid_t))
    end
end

"""
    set_local_func(::Type{F}) where {F<:Filter}

Return a function that sets dataset specific parameters or `nothing` if no function exists.
The function signature is `func(dcpl_id::API.hid_t, type_id::API.hid_t, space_id::API.hid_t)`.
See `API.h5z_register`
"""
set_local_func(::Type{F}) where {F<:Filter} = nothing

"""
    set_local_cfunc(::Type{F}) where {F<:Filter}

Return a C function pointer for the set local function.
By default, this will return the result of using `@cfunction` on the function
specified by `set_local_func(F)` or `C_NULL` if `nothing`.

Overriding this will allow `@cfunction` to return a `Ptr{Nothing}` rather
than a `CFunction`` closure which may not work on all systems.
"""
function set_local_cfunc(::Type{F}) where {F<:Filter}
    func = set_local_func(F)
    if func === nothing
        return C_NULL
    else
        return @cfunction($func, API.herr_t, (API.hid_t,API.hid_t,API.hid_t))
    end
end

"""
    filter_func(::Type{F}) where {F<:Filter}

Returns a function that performs the actual filtering.

See `API.h5z_register`
"""
filter_func(::Type{F}) where {F<:Filter} = nothing

"""
    filter_cfunc(::Type{F}) where {F<:Filter}

Return a C function pointer for the filter function.
By default, this will return the result of using `@cfunction` on the function
specified by `filter_func(F)` or will throw an error if `nothing`.

Overriding this will allow `@cfunction` to return a `Ptr{Nothing}` rather
than a `CFunction`` closure which may not work on all systems.
"""
function filter_cfunc(::Type{F}) where {F<:Filter}
    func = filter_func(F)
    if func === nothing
        error("Filter function for $f must be defined via `filter_func`.")
    end
    c_filter_func = @cfunction($func, Csize_t,
                               (Cuint, Csize_t, Ptr{Cuint}, Csize_t,
                               Ptr{Csize_t}, Ptr{Ptr{Cvoid}}))
    return c_filter_func
end

# Generic implementation of register_filter
"""
    register_filter(::Type{F}) where F <: Filter

Register the filter with the HDF5 library via API.h5z_register.
Also add F to the FILTERS dictionary.
"""
function register_filter(::Type{F}) where F <: Filter
    id = filterid(F)
    encoder = encoder_present(F)
    decoder = decoder_present(F)
    name = filtername(F)
    can_apply = can_apply_cfunc(F)
    set_local = set_local_cfunc(F)
    func = filter_cfunc(F)
    GC.@preserve name begin
        API.h5z_register(API.H5Z_class_t(
            API.H5Z_CLASS_T_VERS,
            id,
            encoder,
            decoder,
            pointer(name),
            can_apply,
            set_local,
            func
        ))
    end
    FILTERS[id] = F
    return nothing
end

struct UnknownFilter <: Filter
    filter_id::API.H5Z_filter_t
    flags::Cuint
    data::Vector{Cuint}
    name::String
    config::Cuint
end
filterid(filter::UnknownFilter) = filter.filter_id
filtername(filter::UnknownFilter) = filter.name
filtername(::Type{UnknownFilter}) = "Unknown Filter"
encoder_present(::Type{UnknownFilter}) = false
decoder_present(::Type{UnknownFilter}) = false

"""
    FilterPipeline(plist::DatasetCreateProperties)

The filter pipeline associated with `plist`. Acts like a `AbstractVector{Filter}`,
supporting the following operations:

- `length(pipeline)`: the number of filters.
- `pipeline[i]` to return the `i`th filter.
- `pipeline[FilterType]` to return a filter of type `FilterType`
- `push!(pipline, filter)` to add an extra filter to the pipeline.
- `append!(pipeline, filters)` to add multiple filters to the pipeline.
- `delete!(pipeline, FilterType)` to remove a filter of type `FilterType` from the pipeline.
- `empty!(pipeline)` to remove all filters from the pipeline.
"""
struct FilterPipeline{P<:Properties} <: AbstractVector{Filter}
    plist::P
end

function Base.length(f::FilterPipeline)
    API.h5p_get_nfilters(f.plist)
end
Base.size(f::FilterPipeline) = (length(f),)

function Base.getindex(f::FilterPipeline, i::Integer)
    id = API.h5p_get_filter(f.plist, i-1, C_NULL, C_NULL, C_NULL, 0, C_NULL, C_NULL)
    F = get(FILTERS, id, UnknownFilter)
    return getindex(f, F, i)
end

function Base.getindex(f::FilterPipeline, ::Type{UnknownFilter}, i::Integer, cd_values::Vector{Cuint} = Cuint[])
    flags = Ref{Cuint}()
    cd_nelmts = Ref{Csize_t}(length(cd_values))
    namebuf = Array{UInt8}(undef, 256)
    config = Ref{Cuint}()
    id = API.h5p_get_filter(f.plist, i-1, flags, cd_nelmts, cd_values, length(namebuf), namebuf, config)
    if cd_nelmts[] < length(cd_values)
        resize!(cd_values, cd_nelmts[])
        return getindex(f, UnknownFilter, i, cd_values)
    end
    resize!(namebuf, findfirst(isequal(0), namebuf)-1)
    resize!(cd_values, cd_nelmts[])
    return UnknownFilter(id, flags[], cd_values, String(namebuf), config[])
end

function Base.getindex(f::FilterPipeline, ::Type{F}, i::Integer) where {F<:Filter}
    @assert isbitstype(F)
    ref = Ref{F}()
    GC.@preserve ref begin
        id = API.h5p_get_filter(f.plist, i-1, C_NULL, div(sizeof(F), sizeof(Cuint)), pointer_from_objref(ref), 0, C_NULL, C_NULL)
    end
    @assert id == filterid(F)
    return ref[]
end
function Base.getindex(f::FilterPipeline, ::Type{F}) where {F<:Filter}
    @assert isbitstype(F)
    ref = Ref{F}()
    GC.@preserve ref begin
        API.h5p_get_filter_by_id(f.plist, filterid(F), C_NULL, div(sizeof(F), sizeof(Cuint)), pointer_from_objref(ref), 0, C_NULL, C_NULL)
    end
    return ref[]
end

function Base.empty!(filters::FilterPipeline)
    API.h5p_remove_filter(filters.plist, API.H5Z_FILTER_ALL)
    return filters
end
function Base.delete!(filters::FilterPipeline, ::Type{F}) where {F<:Filter}
    API.h5p_remove_filter(filters.plist, filterid(F))
    return filters
end
function Base.append!(filters::FilterPipeline, extra)
    for filter in extra
        push!(filters, filter)
    end
    return filters
end
function Base.push!(p::FilterPipeline, f::F) where F <: Filter
    ref = Ref(f)
    GC.@preserve ref begin
        API.h5p_set_filter(p.plist, filterid(F), API.H5Z_FLAG_OPTIONAL, div(sizeof(F), sizeof(Cuint)), pointer_from_objref(ref))
    end
    return p
end
function Base.push!(p::FilterPipeline, f::UnknownFilter)
    GC.@preserve f begin
        API.h5p_set_filter(p.plist, f.filter_id, f.flags, length(f.data), pointer(f.data))
    end
end

include("builtin.jl")

end # module


text:
using Logging

struct ChannelLogger <: AbstractLogger
    c::Channel{Tuple}
    timestamp::Bool
    ChannelLogger(;timestamp=true) = new(Channel{Tuple}(1000), timestamp)
end

const collect_max = 10

function collect_channel(c::AbstractChannel{T}) where T
    v = Vector{T}()
    push!(v, take!(c))
    while !isempty(c) && length(v) < collect_max
        push!(v, take!(c))
    end
    return v
end

Base.collect(l::ChannelLogger) = collect_channel(l.c)
Base.take!(l::ChannelLogger) = take!(l.c)
Base.isready(l::ChannelLogger) = isready(l.c)
Base.isempty(l::ChannelLogger) = isempty(l.c)
drain!(l::ChannelLogger) = while !isempty(l) ; take!(l) ; end

Logging.shouldlog(::ChannelLogger, args...) = true
Logging.min_enabled_level(::ChannelLogger) = Logging.Debug
Logging.catch_exceptions(::ChannelLogger) = false

log_squash(v::Union{Integer,Float64,Bool}) = v
log_squash(v) = string(v)

function Logging.handle_message(l::ChannelLogger,
                                level, message, _module, group, id, file, line;
                                timestamp=nothing, kw...)

    # To avoid feedback loops,
    # don't put TinyRPC Debug logs in the remote logging channel
    if _module == TinyRPC && level < Logging.Info
        return
    end

    # Dump oldest etries if the buffer is full.
    if length(l.c.data) >= l.c.sz_max
        take!(l.c)
    end

    # Convert the kw args into a Vector.
    kwv = Pair{Symbol,Any}[Symbol(k) => log_squash(v) for (k, v) in kw]

    # Add a timestamp.
    if l.timestamp && timestamp == nothing
        timestamp == time()
    end
    if timestamp != nothing
        push!(kwv, :timestamp => timestamp)
    end

    # Convert args to basic types.
    put!(l.c, ((LogLevel(level),
                string(message),
                string(_module),
                string(group),
                Symbol(id),
                string(file),
                Int32(line)),
                kwv))
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

mean(d, v::Symbol, w::WeightingMethod) = mean(d[v], weight(d, w))
mean(d, v::Symbol, s::Number) = mean(d, v, BlockWeighting(ntuple(i->s,embeddim(d))))
mean(d, v::Symbol) = mean(d, v, mode_heuristic(d))
mean(d, w::WeightingMethod) = Dict(v => mean(d, v, w) for v in name.(variables(d)))
mean(d, s::Number) = mean(d, BlockWeighting(ntuple(i->s,embeddim(d))))

"""
    mean(sdata)
    mean(sdata, v)
    mean(sdata, v, s)

Spatial mean of spatial data `sdata`. Optionally,
specify the variable `v` and the block side `s`.
"""
mean(d::Data) = mean(d, mode_heuristic(d))

var(d, v::Symbol, w::WeightingMethod) = var(d[v], weight(d, w), mean=mean(d, v, w), corrected=false)
var(d, v::Symbol, s::Number) = var(d, v, BlockWeighting(ntuple(i->s,embeddim(d))))
var(d, v::Symbol) = var(d, v, mode_heuristic(d))
var(d, w::WeightingMethod) = Dict(v => var(d, v, w) for v in name.(variables(d)))
var(d, s::Number) = var(d, BlockWeighting(ntuple(i->s,embeddim(d))))

"""
    var(sdata)
    var(sdata, v)
    var(sdata, v, s)

Spatial variance of spatial data `sdata`. Optionally,
specify the variable `v` and the block side `s`.
"""
var(d::Data) = var(d, mode_heuristic(d))

quantile(d, v::Symbol, p, w::WeightingMethod) = quantile(d[v], weight(d, w), p)
quantile(d, v::Symbol, p, s::Number) = quantile(d, v, p, BlockWeighting(ntuple(i->s,embeddim(d))))
quantile(d, v::Symbol, p) = quantile(d, v, p, mode_heuristic(d))
quantile(d, p, w::WeightingMethod) = Dict(v => quantile(d, v, p, w) for v in name.(variables(d)))
quantile(d, p::T, s::Number) where {T<:Union{Number,AbstractVector}} = quantile(d, p, BlockWeighting(ntuple(i->s,embeddim(d))))

"""
    quantile(sdata, p)
    quantile(sdata, v, p)
    quantile(sdata, v, p, s)

Spatial quantile of spatial data `sdata` at probability `p`.
Optionally, specify the variable `v` and the block side `s`.
"""
quantile(d::Data, p) = quantile(d, p, mode_heuristic(d))

function dist_matrix_random_sample(d, npoints=1000)
  # select at most 1000 points at random
  nel = nelements(d)
  inds = sample(1:nel, min(nel, npoints), replace=false)
  X = (coordinates(centroid(d, ind)) for ind in inds)
  pairwise(Euclidean(), X)
end

bound_box_constr(d) = minimum(sides(boundingbox(domain(d))))

function median_heuristic(d)
  D = dist_matrix_random_sample(d)
  # median heuristic
  n = size(D, 1)
  m = median(D[i,j] for i in 1:n for j in 1:n if i > j)

  l = bound_box_constr(d)

  min(m, l)
end

"""
    mode_heuristic(d)

Return the estimated mode of the pairwise distances for a set of locations.
"""
function mode_heuristic(d)
  D = dist_matrix_random_sample(d)
  n = size(D, 1)
  m = hsm_mode([D[i,j] for i in 1:n for j in 1:n if i > j])
  l = bound_box_constr(d)
  min(m, l)
end

"""
    hsm_mode(x)

Return the mode of the vector `x`.

## References

* Bickel & Frühwirth, 2005. [On a fast, robust estimator
  of the mode: Comparisons to other robust estimators
  with applications](https://doi.org/10.1016/j.csda.2005.07.011)
"""
hsm_mode(x) = hsm_recursion(sort(x))

function hsm_recursion(x)
  n = length(x)

  # base cases
  n == 1 && return x[1]
  n == 2 && return (x[1] + x[2]) / 2
  if n == 3
    d1 = x[2] - x[1]
    d2 = x[3] - x[2]
    d1 < d2 && return (x[1] + x[2]) / 2
    d1 > d2 && return (x[2] + x[3]) / 2
    d1 == d2 && return x[2]
  end

  # find index of half interval
  k = ceil(Int, n / 2)
  i = argmin([x[j+k] - x[j] for j in 1:n-k])

  # perform recursion
  hsm_recursion(view(x, i:i+k))
end


text:
@inline function apply_nonlinearity(x::Real,::NLRelu)
  return max(x,0.0)
end
struct NLRectifiedQuadratic <: AbstractNonlinearity end
@inline function apply_nonlinearity(x::R,::NLRectifiedQuadratic) where R<:Real
  if x <= 0.0
    return 0.0
  else
    return x*x
  end
end


text:
struct LTTransDist
    terminal::Bool
    rob::Coord
    prev_opp::Coord
    probs::SVector{5, Float64} # corresponds to opp moving in the cardinal directions or staying
                     # sum(probs) is always 1!
end
LTTransDist(rob::Coord, prev_opp::Coord, probs::AbstractVector) = LTTransDist(false, rob, prev_opp, probs)

function Random.rand(rng::AbstractRNG, d::LTTransDist)
    if d.terminal
        return LTState(d.rob, d.prev_opp, true)
    end
    i = sample(rng, Weights(d.probs, 1.0))
    if i < 5
        opp = d.prev_opp + CARDINALS[i]
    else
        opp = d.prev_opp
    end
    return LTState(d.rob, opp, false)
end

function Distributions.pdf(d::LTTransDist, s::LTState)::Float64
    if d.terminal
        return s.terminal ? 1.0 : 0.0
    elseif s.terminal || s.robot != d.rob || sum(abs, s.opponent-d.prev_opp) > 1
        return 0.0
    elseif s.opponent == d.prev_opp
        return d.probs[5]
    else
        dir = s.opponent-d.prev_opp
        if dir[1] == 0
            if dir[2] == 1
                return d.probs[1]
            else
                return d.probs[3]
            end
        elseif dir[1] == 1
            return d.probs[2]
        else
            return d.probs[4]
        end
        # return d.probs[findfirst(CARDINALS, dir)]
    end
end

POMDPs.support(d::LTTransDist) = d

function Base.iterate(d::LTTransDist, i::Int=1)
    if d.terminal
        if i > 1
            return nothing
        else
            return (LTState(d.rob, d.prev_opp, true), i+1)
        end
    end
    while i <= 5 && d.probs[i] == 0.0
        i += 1
    end
    if i > 5
        return nothing
    elseif i <= 4
        return (LTState(d.rob, d.prev_opp+CARDINALS[i], false), i+1)
    else
        return (LTState(d.rob, d.prev_opp, false), i+1)
    end
end

function Base.length(d::LTTransDist)
    if d.terminal
        return 1
    else
        return sum(prob > 0.0 for prob in d.probs)
    end
end

function POMDPs.transition(p::LaserTagPOMDP, s::LTState, a::Int)
    if s.terminal || a == TAG_ACTION && s.robot == s.opponent
        return LTTransDist(true, s.robot, s.opponent, SVector(1., 0., 0., 0., 0.))
    end

    probs = fill!(MVector{5, Float64}(undef), 0.0)

    opp = s.opponent
    rob = s.robot
    f = p.floor
    obst = p.obstacles

    # opponent behavior (see base_tag.cpp line 576)
    # 0.4 chance of moving in x direction
    if opp[1] == rob[1]
        if !opaque(f, obst, opp + Coord(1,0))
            probs[2] += 0.2
        end
        if !opaque(f, obst, opp + Coord(-1,0))
            probs[4] += 0.2
        end
    elseif opp[1] > rob[1] && !opaque(f, obst, opp + Coord(1,0))
        probs[2] += 0.4
    elseif opp[1] < rob[1] && !opaque(f, obst, opp + Coord(-1,0))
        probs[4] += 0.4
    end

    # 0.4 chance of moving in y direction
    if opp[2] == rob[2]
        if !opaque(f, obst, opp + Coord(0,1))
            probs[1] += 0.2
        end
        if !opaque(f, obst, opp + Coord(0,-1))
            probs[3] += 0.2
        end
    elseif opp[2] > rob[2] && !opaque(f, obst, opp + Coord(0,1))
        probs[1] += 0.4
    elseif opp[2] < rob[2] && !opaque(f, obst, opp + Coord(0,-1))
        probs[3] += 0.4
    end

    # 0.2 + all out of bounds mass chance staying the same
    probs[5] = 1.0 - sum(probs)

    next_rob = add_if_clear(f, obst, rob, ACTION_DIRS[a])

    return LTTransDist(next_rob, opp, probs)
end


text:
module firstTask

# Global variable
goldenRatio = (1 + sqrt(5)) / 2

"""
Function optimizes single variable function using golden section method (to find local minimum value in provided range)
"""
function goldenSection(f, a, b, tol = 1e-8)
  # Make sure that entered range is correct
  a = min(a, b)
  b = max(a, b)

  while abs(a - b) >= tol
    # Calculate next 2 points
    x1 = b - (b - a) / goldenRatio
    x2 = a + (b - a) / goldenRatio
    # Compare function values from calculated points and determine where minimum is located
    if f(x1) < f(x2)
      b = x2
    else
      a = x1
    end
  end
  # Get the approximate value and return it
  result = (a + b) / 2
  resultY = f(result)
  return (result, resultY)
end

export goldenSection

end # module

max_stars_repo_name: KK-studia/firstTask

max_issues_repo_name: KK-studia/firstTask

max_forks_repo_name: KK-studia/firstTask


text:
# External packages
using Revise, DelimitedFiles, SymPy

# Add Exasim to Julia search path
cdir = pwd(); ii = findlast("Exasim", cdir);
include(cdir[1:ii[end]] * "/Installation/setpath.jl");

# Exasim packages
using Preprocessing, Mesh, Gencode, Postprocessing

# create pde structure and mesh structure
pde, mesh = Preprocessing.initializeexasim();

# Define PDE model: governing equations, initial solutions, and boundary conditions
pde.model = "ModelD";            # ModelC, ModelD, ModelW
include("pdemodel.jl");          # include the PDE model file

# Set discretization parameters, physical parameters, and solver parameters
pde.porder = 3;          # polynomial degree
pde.torder = 2;          # time-stepping order of accuracy
pde.nstage = 2;          # time-stepping number of stages
pde.dt = 0.05*ones(50);   # time step sizes
pde.soltime = collect(10:10:length(pde.dt)); # steps at which solution are collected
pde.visdt = 0.05; # visualization timestep size

gam = 1.4;                      # specific heat ratio
Re = 1600;                      # Reynolds number
Pr = 0.71;                      # Prandtl number
Minf = 0.2;                     # Mach number
pde.physicsparam = [gam Re Pr Minf];
pde.tau = [2.0];                # DG stabilization parameter
pde.GMRESrestart=40;
pde.linearsolvertol=0.0001;
pde.linearsolveriter=41;
pde.precMatrixType=2;
pde.NLiter=2;

# Choose computing platform and set number of processors
#pde.platform = "gpu";           # choose this option if NVIDIA GPUs are available
pde.mpiprocs = 4;                # number of MPI processors

# create a linear mesh for a cube domain
mesh.p, mesh.t = Mesh.cubemesh(16,16,16,1);
mesh.p = 2*pi*mesh.p;
# expressions for disjoint boundaries
mesh.boundaryexpr = [p -> (p[2,:] .< 1e-3), p -> (p[1,:] .> 2*pi-1e-3), p -> (p[2,:] .> 2*pi-1e-3), p -> (p[1,:] .< 1e-3), p -> (p[3,:] .< 1e-3), p -> (p[3,:] .> 2*pi-1e-3)];
mesh.boundarycondition = [1 1 1 1 1 1]; # Set boundary condition for each disjoint boundary
mesh.periodicexpr = [2 p->p[2:3,:] 4 p->p[2:3,:]; 1 p->p[[1,3],:] 3 p->p[[1,3],:]; 5 p->p[[1,2],:] 6 p->p[[1,2],:]];

# call exasim to generate and run C++ code to solve the PDE model
sol, pde, mesh,~,~,~,~  = Postprocessing.exasim(pde,mesh);