text:
function reconstruct!(edge_u, u, balance_law, meshx::UniformPeriodicMesh1D,
                      reconstruction!, stencil_width::Val{1}, parallel)
    # no stencils with overlap; set mean values as reconstructions
    @inbounds for cell in 1:length(u)
        edge_u[1,cell] = edge_u[2,cell] = u[cell]
    end
    nothing
end

function reconstruct!(edge_u, u, balance_law, meshx::UniformPeriodicMesh1D,
                      reconstruction!, stencil_width::Val{3}, parallel)
    # stencils with overlap
    @inbounds begin
        cell = 1
        u_m1 = u[end]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        reconstruction!(edge_u, cell, u_m1, u_0, u_p1, balance_law, meshx)

        cell = length(u)
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[1]
        reconstruction!(edge_u, cell, u_m1, u_0, u_p1, balance_law, meshx)
    end

    reconstruct_inner_loop!(edge_u, u, balance_law, meshx, reconstruction!, stencil_width, parallel)
end

@inline function reconstruct_inner_loop!(edge_u, u, balance_law, meshx::UniformPeriodicMesh1D,
                                         reconstruction!, stencil_width::Val{3}, parallel)
    # serial
    @inbounds for cell in 2:length(u)-1
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        reconstruction!(edge_u, cell, u_m1, u_0, u_p1, balance_law, meshx)
    end
    nothing
end

@inline function reconstruct_inner_loop!(edge_u, u, balance_law, meshx::UniformPeriodicMesh1D,
                                         reconstruction!, stencil_width::Val{3}, ::Val{:threads})
    # threaded
    @inbounds Threads.@threads for cell in 2:length(u)-1
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        reconstruction!(edge_u, cell, u_m1, u_0, u_p1, balance_law, meshx)
    end
    nothing
end

function reconstruct!(edge_u, u, balance_law, meshx::UniformPeriodicMesh1D,
                      reconstruction!, stencil_width::Val{5}, parallel)
    # stencils with overlap
    @inbounds begin
        cell = 1
        u_m2 = u[end-1]
        u_m1 = u[end]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        reconstruction!(edge_u, cell, u_m2, u_m1, u_0, u_p1, u_p2, balance_law, meshx)

        cell = 2
        u_m2 = u[end]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        reconstruction!(edge_u, cell, u_m2, u_m1, u_0, u_p1, u_p2, balance_law, meshx)

        cell = length(u)-1
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[1]
        reconstruction!(edge_u, cell, u_m2, u_m1, u_0, u_p1, u_p2, balance_law, meshx)

        cell = length(u)
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[1]
        u_p2 = u[2]
        reconstruction!(edge_u, cell, u_m2, u_m1, u_0, u_p1, u_p2, balance_law, meshx)
    end

    reconstruct_inner_loop!(edge_u, u, balance_law, meshx, reconstruction!, stencil_width, parallel)
end

@inline function reconstruct_inner_loop!(edge_u, u, balance_law, meshx::UniformPeriodicMesh1D,
                                         reconstruction!, stencil_width::Val{5}, parallel)
    # serial
    @inbounds for cell in 3:length(u)-2
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        reconstruction!(edge_u, cell, u_m2, u_m1, u_0, u_p1, u_p2, balance_law, meshx)
    end
    nothing
end

@inline function reconstruct_inner_loop!(edge_u, u, balance_law, meshx::UniformPeriodicMesh1D,
                                         reconstruction!, stencil_width::Val{5}, ::Val{:threads})
    # threaded
    @inbounds Threads.@threads for cell in 3:length(u)-2
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        reconstruction!(edge_u, cell, u_m2, u_m1, u_0, u_p1, u_p2, balance_law, meshx)
    end
    nothing
end

function reconstruct!(edge_u, u, balance_law, meshx::UniformPeriodicMesh1D,
                      reconstruction!, stencil_width::Val{7}, parallel)
    # stencils with overlap
    @inbounds begin
        cell = 1
        u_m3 = u[end-2]
        u_m2 = u[end-1]
        u_m1 = u[end]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        reconstruction!(edge_u, cell, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, balance_law, meshx)

        cell = 2
        u_m3 = u[end-1]
        u_m2 = u[end]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        reconstruction!(edge_u, cell, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, balance_law, meshx)

        cell = 3
        u_m3 = u[end]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        reconstruction!(edge_u, cell, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, balance_law, meshx)

        cell = length(u)-2
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[1]
        reconstruction!(edge_u, cell, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, balance_law, meshx)

        cell = length(u)-1
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[1]
        u_p3 = u[2]
        reconstruction!(edge_u, cell, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, balance_law, meshx)

        cell = length(u)
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[1]
        u_p2 = u[2]
        u_p3 = u[3]
        reconstruction!(edge_u, cell, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, balance_law, meshx)
    end

    reconstruct_inner_loop!(edge_u, u, balance_law, meshx, reconstruction!, stencil_width, parallel)
end

@inline function reconstruct_inner_loop!(edge_u, u, balance_law, meshx::UniformPeriodicMesh1D,
                                         reconstruction!, stencil_width::Val{7}, parallel)
    # serial
    @inbounds for cell in 4:length(u)-3
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        reconstruction!(edge_u, cell, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, balance_law, meshx)
    end
    nothing
end

@inline function reconstruct_inner_loop!(edge_u, u, balance_law, meshx::UniformPeriodicMesh1D,
                                         reconstruction!, stencil_width::Val{7}, ::Val{:threads})
    # threaded
    @inbounds Threads.@threads for cell in 4:length(u)-3
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        reconstruction!(edge_u, cell, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, balance_law, meshx)
    end
    nothing
end

function reconstruct!(edge_u, u, balance_law, meshx::UniformPeriodicMesh1D,
                      reconstruction!, stencil_width::Val{9}, parallel)
    # stencils with overlap
    @inbounds begin
        cell = 1
        u_m4 = u[end-3]
        u_m3 = u[end-2]
        u_m2 = u[end-1]
        u_m1 = u[end]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        reconstruction!(edge_u, cell, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, balance_law, meshx)

        cell = 2
        u_m4 = u[end-2]
        u_m3 = u[end-1]
        u_m2 = u[end]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        reconstruction!(edge_u, cell, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, balance_law, meshx)

        cell = 3
        u_m4 = u[end-1]
        u_m3 = u[end]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        reconstruction!(edge_u, cell, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, balance_law, meshx)

        cell = 4
        u_m4 = u[end]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        reconstruction!(edge_u, cell, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, balance_law, meshx)

        cell = length(u)-3
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[1]
        reconstruction!(edge_u, cell, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, balance_law, meshx)

        cell = length(u)-2
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[1]
        u_p4 = u[2]
        reconstruction!(edge_u, cell, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, balance_law, meshx)

        cell = length(u)-1
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[1]
        u_p3 = u[2]
        u_p4 = u[3]
        reconstruction!(edge_u, cell, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, balance_law, meshx)

        cell = length(u)
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[1]
        u_p2 = u[2]
        u_p3 = u[3]
        u_p4 = u[4]
        reconstruction!(edge_u, cell, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, balance_law, meshx)
    end

    reconstruct_inner_loop!(edge_u, u, balance_law, meshx, reconstruction!, stencil_width, parallel)
end

@inline function reconstruct_inner_loop!(edge_u, u, balance_law, meshx::UniformPeriodicMesh1D,
                                         reconstruction!, stencil_width::Val{9}, parallel)
    # serial
    @inbounds for cell in 5:length(u)-4
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        reconstruction!(edge_u, cell, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, balance_law, meshx)
    end
    nothing
end

@inline function reconstruct_inner_loop!(edge_u, u, balance_law, meshx::UniformPeriodicMesh1D,
                                         reconstruction!, stencil_width::Val{9}, ::Val{:threads})
    # threaded
    @inbounds Threads.@threads for cell in 5:length(u)-4
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        reconstruction!(edge_u, cell, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, balance_law, meshx)
    end
    nothing
end

function reconstruct!(edge_u, u, balance_law, meshx::UniformPeriodicMesh1D,
                      reconstruction!, stencil_width::Val{11}, parallel)
    # stencils with overlap
    @inbounds begin
        cell = 1
        u_m5 = u[end-4]
        u_m4 = u[end-3]
        u_m3 = u[end-2]
        u_m2 = u[end-1]
        u_m1 = u[end]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        reconstruction!(edge_u, cell, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, balance_law, meshx)

        cell = 2
        u_m5 = u[end-3]
        u_m4 = u[end-2]
        u_m3 = u[end-1]
        u_m2 = u[end]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        reconstruction!(edge_u, cell, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, balance_law, meshx)

        cell = 3
        u_m5 = u[end-2]
        u_m4 = u[end-1]
        u_m3 = u[end]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        reconstruction!(edge_u, cell, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, balance_law, meshx)

        cell = 4
        u_m5 = u[end-1]
        u_m4 = u[end]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        reconstruction!(edge_u, cell, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, balance_law, meshx)

        cell = 5
        u_m5 = u[end]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        reconstruction!(edge_u, cell, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, balance_law, meshx)

        cell = length(u)-4
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[1]
        reconstruction!(edge_u, cell, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, balance_law, meshx)

        cell = length(u)-3
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[1]
        u_p5 = u[2]
        reconstruction!(edge_u, cell, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, balance_law, meshx)

        cell = length(u)-2
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[1]
        u_p4 = u[2]
        u_p5 = u[3]
        reconstruction!(edge_u, cell, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, balance_law, meshx)

        cell = length(u)-1
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[1]
        u_p3 = u[2]
        u_p4 = u[3]
        u_p5 = u[4]
        reconstruction!(edge_u, cell, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, balance_law, meshx)

        cell = length(u)
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[1]
        u_p2 = u[2]
        u_p3 = u[3]
        u_p4 = u[4]
        u_p5 = u[5]
        reconstruction!(edge_u, cell, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, balance_law, meshx)
    end

    reconstruct_inner_loop!(edge_u, u, balance_law, meshx, reconstruction!, stencil_width, parallel)
end

@inline function reconstruct_inner_loop!(edge_u, u, balance_law, meshx::UniformPeriodicMesh1D,
                                         reconstruction!, stencil_width::Val{11}, parallel)
    # serial
    @inbounds for cell in 6:length(u)-5
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        reconstruction!(edge_u, cell, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, balance_law, meshx)
    end
    nothing
end

@inline function reconstruct_inner_loop!(edge_u, u, balance_law, meshx::UniformPeriodicMesh1D,
                                         reconstruction!, stencil_width::Val{11}, ::Val{:threads})
    # threaded
    @inbounds Threads.@threads for cell in 6:length(u)-5
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        reconstruction!(edge_u, cell, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, balance_law, meshx)
    end
    nothing
end

function reconstruct!(edge_u, u, balance_law, meshx::UniformPeriodicMesh1D,
                      reconstruction!, stencil_width::Val{13}, parallel)
    # stencils with overlap
    @inbounds begin
        cell = 1
        u_m6 = u[end-5]
        u_m5 = u[end-4]
        u_m4 = u[end-3]
        u_m3 = u[end-2]
        u_m2 = u[end-1]
        u_m1 = u[end]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        reconstruction!(edge_u, cell, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, balance_law, meshx)

        cell = 2
        u_m6 = u[end-4]
        u_m5 = u[end-3]
        u_m4 = u[end-2]
        u_m3 = u[end-1]
        u_m2 = u[end]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        reconstruction!(edge_u, cell, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, balance_law, meshx)

        cell = 3
        u_m6 = u[end-3]
        u_m5 = u[end-2]
        u_m4 = u[end-1]
        u_m3 = u[end]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        reconstruction!(edge_u, cell, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, balance_law, meshx)

        cell = 4
        u_m6 = u[end-2]
        u_m5 = u[end-1]
        u_m4 = u[end]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        reconstruction!(edge_u, cell, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, balance_law, meshx)

        cell = 5
        u_m6 = u[end-1]
        u_m5 = u[end]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        reconstruction!(edge_u, cell, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, balance_law, meshx)

        cell = 6
        u_m6 = u[end]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        reconstruction!(edge_u, cell, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, balance_law, meshx)

        cell = length(u)-5
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[1]
        reconstruction!(edge_u, cell, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, balance_law, meshx)

        cell = length(u)-4
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[1]
        u_p6 = u[2]
        reconstruction!(edge_u, cell, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, balance_law, meshx)

        cell = length(u)-3
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[1]
        u_p5 = u[2]
        u_p6 = u[3]
        reconstruction!(edge_u, cell, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, balance_law, meshx)

        cell = length(u)-2
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[1]
        u_p4 = u[2]
        u_p5 = u[3]
        u_p6 = u[4]
        reconstruction!(edge_u, cell, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, balance_law, meshx)

        cell = length(u)-1
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[1]
        u_p3 = u[2]
        u_p4 = u[3]
        u_p5 = u[4]
        u_p6 = u[5]
        reconstruction!(edge_u, cell, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, balance_law, meshx)

        cell = length(u)
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[1]
        u_p2 = u[2]
        u_p3 = u[3]
        u_p4 = u[4]
        u_p5 = u[5]
        u_p6 = u[6]
        reconstruction!(edge_u, cell, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, balance_law, meshx)
    end

    reconstruct_inner_loop!(edge_u, u, balance_law, meshx, reconstruction!, stencil_width, parallel)
end

@inline function reconstruct_inner_loop!(edge_u, u, balance_law, meshx::UniformPeriodicMesh1D,
                                         reconstruction!, stencil_width::Val{13}, parallel)
    # serial
    @inbounds for cell in 7:length(u)-6
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        reconstruction!(edge_u, cell, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, balance_law, meshx)
    end
    nothing
end

@inline function reconstruct_inner_loop!(edge_u, u, balance_law, meshx::UniformPeriodicMesh1D,
                                         reconstruction!, stencil_width::Val{13}, ::Val{:threads})
    # threaded
    @inbounds Threads.@threads for cell in 7:length(u)-6
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        reconstruction!(edge_u, cell, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, balance_law, meshx)
    end
    nothing
end

function reconstruct!(edge_u, u, balance_law, meshx::UniformPeriodicMesh1D,
                      reconstruction!, stencil_width::Val{15}, parallel)
    # stencils with overlap
    @inbounds begin
        cell = 1
        u_m7 = u[end-6]
        u_m6 = u[end-5]
        u_m5 = u[end-4]
        u_m4 = u[end-3]
        u_m3 = u[end-2]
        u_m2 = u[end-1]
        u_m1 = u[end]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        u_p7 = u[cell+7]
        reconstruction!(edge_u, cell, u_m7, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, u_p7, balance_law, meshx)

        cell = 2
        u_m7 = u[end-5]
        u_m6 = u[end-4]
        u_m5 = u[end-3]
        u_m4 = u[end-2]
        u_m3 = u[end-1]
        u_m2 = u[end]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        u_p7 = u[cell+7]
        reconstruction!(edge_u, cell, u_m7, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, u_p7, balance_law, meshx)

        cell = 3
        u_m7 = u[end-4]
        u_m6 = u[end-3]
        u_m5 = u[end-2]
        u_m4 = u[end-1]
        u_m3 = u[end]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        u_p7 = u[cell+7]
        reconstruction!(edge_u, cell, u_m7, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, u_p7, balance_law, meshx)

        cell = 4
        u_m7 = u[end-3]
        u_m6 = u[end-2]
        u_m5 = u[end-1]
        u_m4 = u[end]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        u_p7 = u[cell+7]
        reconstruction!(edge_u, cell, u_m7, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, u_p7, balance_law, meshx)

        cell = 5
        u_m7 = u[end-2]
        u_m6 = u[end-1]
        u_m5 = u[end]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        u_p7 = u[cell+7]
        reconstruction!(edge_u, cell, u_m7, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, u_p7, balance_law, meshx)

        cell = 6
        u_m7 = u[end-1]
        u_m6 = u[end]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        u_p7 = u[cell+7]
        reconstruction!(edge_u, cell, u_m7, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, u_p7, balance_law, meshx)

        cell = 7
        u_m7 = u[end]
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        u_p7 = u[cell+7]
        reconstruction!(edge_u, cell, u_m7, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, u_p7, balance_law, meshx)

        cell = length(u)-6
        u_m7 = u[cell-7]
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        u_p7 = u[1]
        reconstruction!(edge_u, cell, u_m7, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, u_p7, balance_law, meshx)

        cell = length(u)-5
        u_m7 = u[cell-7]
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[1]
        u_p7 = u[2]
        reconstruction!(edge_u, cell, u_m7, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, u_p7, balance_law, meshx)

        cell = length(u)-4
        u_m7 = u[cell-7]
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[1]
        u_p6 = u[2]
        u_p7 = u[3]
        reconstruction!(edge_u, cell, u_m7, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, u_p7, balance_law, meshx)

        cell = length(u)-3
        u_m7 = u[cell-7]
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[1]
        u_p5 = u[2]
        u_p6 = u[3]
        u_p7 = u[4]
        reconstruction!(edge_u, cell, u_m7, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, u_p7, balance_law, meshx)

        cell = length(u)-2
        u_m7 = u[cell-7]
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[1]
        u_p4 = u[2]
        u_p5 = u[3]
        u_p6 = u[4]
        u_p7 = u[5]
        reconstruction!(edge_u, cell, u_m7, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, u_p7, balance_law, meshx)

        cell = length(u)-1
        u_m7 = u[cell-7]
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[1]
        u_p3 = u[2]
        u_p4 = u[3]
        u_p5 = u[4]
        u_p6 = u[5]
        u_p7 = u[6]
        reconstruction!(edge_u, cell, u_m7, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, u_p7, balance_law, meshx)

        cell = length(u)
        u_m7 = u[cell-7]
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[1]
        u_p2 = u[2]
        u_p3 = u[3]
        u_p4 = u[4]
        u_p5 = u[5]
        u_p6 = u[6]
        u_p7 = u[7]
        reconstruction!(edge_u, cell, u_m7, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, u_p7, balance_law, meshx)
    end

    reconstruct_inner_loop!(edge_u, u, balance_law, meshx, reconstruction!, stencil_width, parallel)
end

@inline function reconstruct_inner_loop!(edge_u, u, balance_law, meshx::UniformPeriodicMesh1D,
                                         reconstruction!, stencil_width::Val{15}, parallel)
    # serial
    @inbounds for cell in 8:length(u)-7
        u_m7 = u[cell-7]
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        u_p7 = u[cell+7]
        reconstruction!(edge_u, cell, u_m7, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, u_p7, balance_law, meshx)
    end
    nothing
end

@inline function reconstruct_inner_loop!(edge_u, u, balance_law, meshx::UniformPeriodicMesh1D,
                                         reconstruction!, stencil_width::Val{15}, ::Val{:threads})
    # threaded
    @inbounds Threads.@threads for cell in 8:length(u)-7
        u_m7 = u[cell-7]
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        u_p7 = u[cell+7]
        reconstruction!(edge_u, cell, u_m7, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, u_p7, balance_law, meshx)
    end
    nothing
end

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

function interpolate!(uval, ξ, cell, u, balance_law, mesh::UniformPeriodicMesh1D,
                      reconstruction!, stencil_width::Val{1})
    # no stencils with overlap; set mean values as reconstructions
    @inbounds for idx in eachindex(uval)
        uval[idx] = u[cell]
    end
    nothing
end

function interpolate!(uval, ξ, cell, u, balance_law, mesh::UniformPeriodicMesh1D,
                      reconstruction, stencil_width::Val{3})
    # stencils with overlap
    @inbounds if cell == 1
        u_m1 = u[end]
        u_0  = u[cell]
        u_p1 = u[cell+1]
    elseif cell == length(u)
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[1]
    else
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
    end

    interpolate!(uval, ξ, u_m1, u_0, u_p1, balance_law, mesh, reconstruction)

    nothing
end

function interpolate!(uval, ξ, cell, u, balance_law, mesh::UniformPeriodicMesh1D,
                      reconstruction, stencil_width::Val{5})
    # stencils with overlap
    @inbounds if cell == 1
        u_m2 = u[end-1]
        u_m1 = u[end]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
    elseif cell == 2
        u_m2 = u[end]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
    elseif cell == length(u)-1
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[1]
    elseif cell == length(u)
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[1]
        u_p2 = u[2]
    else
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
    end

    interpolate!(uval, ξ, u_m2, u_m1, u_0, u_p1, u_p2, balance_law, mesh, reconstruction)

    nothing
end

function interpolate!(uval, ξ, cell, u, balance_law, mesh::UniformPeriodicMesh1D,
                      reconstruction, stencil_width::Val{7})
    # stencils with overlap
    @inbounds if cell == 1
        u_m3 = u[end-2]
        u_m2 = u[end-1]
        u_m1 = u[end]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
    elseif cell == 2
        u_m3 = u[end-1]
        u_m2 = u[end]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
    elseif cell == 3
        u_m3 = u[end]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
    elseif cell == length(u)-2
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[1]
    elseif cell == length(u)-1
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[1]
        u_p3 = u[2]
    elseif cell == length(u)
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[1]
        u_p2 = u[2]
        u_p3 = u[3]
    else
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
    end

    interpolate!(uval, ξ, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, balance_law, mesh, reconstruction)

    nothing
end

function interpolate!(uval, ξ, cell, u, balance_law, mesh::UniformPeriodicMesh1D,
                      reconstruction, stencil_width::Val{9})
    # stencils with overlap
    @inbounds if cell == 1
        u_m4 = u[end-3]
        u_m3 = u[end-2]
        u_m2 = u[end-1]
        u_m1 = u[end]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
    elseif cell == 2
        u_m4 = u[end-2]
        u_m3 = u[end-1]
        u_m2 = u[end]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
    elseif cell == 3
        u_m4 = u[end-1]
        u_m3 = u[end]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
    elseif cell == 4
        u_m4 = u[end]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
    elseif cell == length(u)-3
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[1]
    elseif cell == length(u)-2
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[1]
        u_p4 = u[2]
    elseif cell == length(u)-1
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[1]
        u_p3 = u[2]
        u_p4 = u[3]
    elseif cell == length(u)
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[1]
        u_p2 = u[2]
        u_p3 = u[3]
        u_p4 = u[4]
    else
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
    end

    interpolate!(uval, ξ, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, balance_law, mesh, reconstruction)

    nothing
end

function interpolate!(uval, ξ, cell, u, balance_law, mesh::UniformPeriodicMesh1D,
                      reconstruction, stencil_width::Val{11})
    # stencils with overlap
    @inbounds if cell == 1
        u_m5 = u[end-4]
        u_m4 = u[end-3]
        u_m3 = u[end-2]
        u_m2 = u[end-1]
        u_m1 = u[end]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
    elseif cell == 2
        u_m5 = u[end-3]
        u_m4 = u[end-2]
        u_m3 = u[end-1]
        u_m2 = u[end]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
    elseif cell == 3
        u_m5 = u[end-2]
        u_m4 = u[end-1]
        u_m3 = u[end]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
    elseif cell == 4
        u_m5 = u[end-1]
        u_m4 = u[end]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
    elseif cell == 5
        u_m5 = u[end]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
    elseif cell == length(u)-4
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[1]
    elseif cell == length(u)-3
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[1]
        u_p5 = u[2]
    elseif cell == length(u)-2
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[1]
        u_p4 = u[2]
        u_p5 = u[3]
    elseif cell == length(u)-1
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[1]
        u_p3 = u[2]
        u_p4 = u[3]
        u_p5 = u[4]
    elseif cell == length(u)
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[1]
        u_p2 = u[2]
        u_p3 = u[3]
        u_p4 = u[4]
        u_p5 = u[5]
    else
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
    end

    interpolate!(uval, ξ, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, balance_law, mesh, reconstruction)

    nothing
end

function interpolate!(uval, ξ, cell, u, balance_law, mesh::UniformPeriodicMesh1D,
                      reconstruction, stencil_width::Val{13})
    # stencils with overlap
    @inbounds if cell == 1
        u_m6 = u[end-5]
        u_m5 = u[end-4]
        u_m4 = u[end-3]
        u_m3 = u[end-2]
        u_m2 = u[end-1]
        u_m1 = u[end]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
    elseif cell == 2
        u_m6 = u[end-4]
        u_m5 = u[end-3]
        u_m4 = u[end-2]
        u_m3 = u[end-1]
        u_m2 = u[end]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
    elseif cell == 3
        u_m6 = u[end-3]
        u_m5 = u[end-2]
        u_m4 = u[end-1]
        u_m3 = u[end]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
    elseif cell == 4
        u_m6 = u[end-2]
        u_m5 = u[end-1]
        u_m4 = u[end]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
    elseif cell == 5
        u_m6 = u[end-1]
        u_m5 = u[end]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
    elseif cell == 6
        u_m6 = u[end]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
    elseif cell == length(u)-5
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[1]
    elseif cell == length(u)-4
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[1]
        u_p6 = u[2]
    elseif cell == length(u)-3
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[1]
        u_p5 = u[2]
        u_p6 = u[3]
    elseif cell == length(u)-2
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[1]
        u_p4 = u[2]
        u_p5 = u[3]
        u_p6 = u[4]
    elseif cell == length(u)-1
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[1]
        u_p3 = u[2]
        u_p4 = u[3]
        u_p5 = u[4]
        u_p6 = u[5]
    elseif cell == length(u)
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[1]
        u_p2 = u[2]
        u_p3 = u[3]
        u_p4 = u[4]
        u_p5 = u[5]
        u_p6 = u[6]
    else
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
    end

    interpolate!(uval, ξ, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, balance_law, mesh, reconstruction)

    nothing
end

function interpolate!(uval, ξ, cell, u, balance_law, mesh::UniformPeriodicMesh1D,
                      reconstruction, stencil_width::Val{15})
    # stencils with overlap
    @inbounds if cell == 1
        u_m7 = u[end-6]
        u_m6 = u[end-5]
        u_m5 = u[end-4]
        u_m4 = u[end-3]
        u_m3 = u[end-2]
        u_m2 = u[end-1]
        u_m1 = u[end]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        u_p7 = u[cell+7]
    elseif cell == 2
        u_m7 = u[end-5]
        u_m6 = u[end-4]
        u_m5 = u[end-3]
        u_m4 = u[end-2]
        u_m3 = u[end-1]
        u_m2 = u[end]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        u_p7 = u[cell+7]
    elseif cell == 3
        u_m7 = u[end-4]
        u_m6 = u[end-3]
        u_m5 = u[end-2]
        u_m4 = u[end-1]
        u_m3 = u[end]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        u_p7 = u[cell+7]
    elseif cell == 4
        u_m7 = u[end-3]
        u_m6 = u[end-2]
        u_m5 = u[end-1]
        u_m4 = u[end]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        u_p7 = u[cell+7]
    elseif cell == 5
        u_m7 = u[end-2]
        u_m6 = u[end-1]
        u_m5 = u[end]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        u_p7 = u[cell+7]
    elseif cell == 6
        u_m7 = u[end-1]
        u_m6 = u[end]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        u_p7 = u[cell+7]
    elseif cell == 7
        u_m7 = u[end]
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        u_p7 = u[cell+7]
    elseif cell == length(u)-6
        u_m7 = u[cell-7]
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        u_p7 = u[1]
    elseif cell == length(u)-5
        u_m7 = u[cell-7]
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[1]
        u_p7 = u[2]
    elseif cell == length(u)-4
        u_m7 = u[cell-7]
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[1]
        u_p6 = u[2]
        u_p7 = u[3]
    elseif cell == length(u)-3
        u_m7 = u[cell-7]
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[1]
        u_p5 = u[2]
        u_p6 = u[3]
        u_p7 = u[4]
    elseif cell == length(u)-2
        u_m7 = u[cell-7]
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[1]
        u_p4 = u[2]
        u_p5 = u[3]
        u_p6 = u[4]
        u_p7 = u[5]
    elseif cell == length(u)-1
        u_m7 = u[cell-7]
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[1]
        u_p3 = u[2]
        u_p4 = u[3]
        u_p5 = u[4]
        u_p6 = u[5]
        u_p7 = u[6]
    elseif cell == length(u)
        u_m7 = u[cell-7]
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[1]
        u_p2 = u[2]
        u_p3 = u[3]
        u_p4 = u[4]
        u_p5 = u[5]
        u_p6 = u[6]
        u_p7 = u[7]
    else
        u_m7 = u[cell-7]
        u_m6 = u[cell-6]
        u_m5 = u[cell-5]
        u_m4 = u[cell-4]
        u_m3 = u[cell-3]
        u_m2 = u[cell-2]
        u_m1 = u[cell-1]
        u_0  = u[cell]
        u_p1 = u[cell+1]
        u_p2 = u[cell+2]
        u_p3 = u[cell+3]
        u_p4 = u[cell+4]
        u_p5 = u[cell+5]
        u_p6 = u[cell+6]
        u_p7 = u[cell+7]
    end

    interpolate!(uval, ξ, u_m7, u_m6, u_m5, u_m4, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, u_p4, u_p5, u_p6, u_p7, balance_law, mesh, reconstruction)

    nothing
end


text:
using JLD
using PyPlot
include("utilities.jl")
include("self_confidence.jl")
pygui(false)

searchdir(path,key) = filter(x->contains(x,key), readdir(path))

folder = "gcloud_data"

for file in searchdir(folder,"jld")
    println("Processing: $file")
    fname = split("$(folder)/$file",".")[1]
    fsize = 12

    data = JLD.load("$fname.jld")
    max_steps = data["max_steps"]
    u_vals = data["u_vals"]
    its_axis = data["its_axis"]
    d_axis = data["d_axis"]
    DA = data["DA"]
    R = data["R"]
    rewards = data["rewards"]
    IA = data["IA"]
    ST = data["ST"]

    fig,ax_ary = PyPlot.subplots(1,1,sharex=true)
    fig[:set_size_inches](8.0,4.0)

    (x_vals,x_ticks,x_lbls,x_rng) = return_indices(DA[:])

    #  vline_frac = (mdp.road_net.gprops[:net_stats].diam - x_rng[1])/(x_rng[2]-x_rng[1])
    #  vline_tick = minimum(x_ticks)+vline_frac*(maximum(x_ticks)-minimum(x_ticks))

    # points = # of points where the Kernel Density Estimator is estimated
    # bw_method = bandwidth of the KDE
    trusted_solver_num = 9
    #  ax_ary[:violinplot](R[1:end .!= trusted_solver_num,:]',x_ticks[1:end .!= trusted_solver_num],widths=0.5,points=250,showmeans=false,showmedians=true,showextrema=true)
    tbp = ax_ary[:boxplot](R[trusted_solver_num,:]',positions=[x_ticks[trusted_solver_num]],widths=0.25,notch=false,showfliers=false,autorange=true,patch_artist=true)
    cbp = ax_ary[:boxplot](R[1:end .!= trusted_solver_num,:]',positions=x_ticks[1:end .!= trusted_solver_num],widths=0.25,notch=false,showfliers=false,autorange=true)
    #  for x in x_ticks
        #  ax_ary[:scatter](float(x)*ones(length(R[x,:])),R[x,:],alpha=0.1)
    #  end
    #  trusted_violin = ax_ary[:violinplot](R[trusted_solver_num,:]',[x_ticks[trusted_solver_num]],widths=0.5,points=250,showmeans=false,showmedians=true,showextrema=true)
    for element in ["boxes","whiskers","means","caps"]
        PyPlot.setp(cbp[element],color=:red)
    end
    for element in ["boxes","whiskers","means","caps"]
        PyPlot.setp(tbp[element],color=:blue)
    end
    for patch in tbp["boxes"]
        patch[:set_facecolor](:blue)
        patch[:set_alpha](0.35)
    end

    for i in x_ticks
        #  display(R[i,:])
        SQ = X3(R[i,:],R[trusted_solver_num,:],global_rwd_range=[minimum(R),maximum(R)])
        println("SQ: $SQ")
        #  @info "X3 at $i: $SQ"
        ax_ary[:annotate]("SQ: $(@sprintf("%0.2f",SQ))",xy=(i,mean(R[i,:])),size=8)
    end

    ax_ary[:set_ylabel]("Reward",size=fsize)
    ax_ary[:set_xticks](x_ticks)
    ax_ary[:xaxis][:set_ticklabels](x_lbls)
    ax_ary[:set_xlabel]("MCTS Depth",size=fsize)
    ax_ary[:set_title]("Reward vs. MCTS Depth",size=fsize)

    PyPlot.savefig("$fname.pdf",dpi=300,transparent=true)
end

ext: jl

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
## Define the new measure evaluation method
# Make alias for our new method
struct NewUniEvalMethod <: InfiniteOpt.MeasureToolbox.AbstractUnivariateMethod end
struct NewMultiEvalMethod <: InfiniteOpt.MeasureToolbox.AbstractMultivariateMethod end

# Extend generate_supports_and_coeffs for scalar params
function InfiniteOpt.MeasureToolbox.generate_integral_data(
    pref::InfiniteOpt.GeneralVariableRef,
    lower_bound::Real,
    upper_bound::Real,
    method::Type{NewUniEvalMethod};
    num_supports::Int = InfiniteOpt.DefaultNumSupports,
    weight_func::Function
    )::InfiniteOpt.AbstractMeasureData # REPLACE WITH ACTUAL ALIAS
    # REPLACE WITH ACTUAL FUNCTIONALITY
    # ADD CHECKS IF NECESSARY
    increment = (upper_bound - lower_bound) / (num_supports - 1)
    supports = [lower_bound + (i - 1) * increment for i in 1:num_supports]
    coeffs = ones(num_supports) / num_supports * (upper_bound - lower_bound)
    # MAKE SURE IT RETURNS APPROPRIATE MEASURE DATA
    return InfiniteOpt.DiscreteMeasureData(
        pref, coeffs, supports,
        weight_function = weight_func,
        lower_bound = lower_bound,
        upper_bound = upper_bound)
end

# Extend generate_supports_and_coeffs for multi-dimensional params
function InfiniteOpt.MeasureToolbox.generate_integral_data(
    prefs::Vector{InfiniteOpt.GeneralVariableRef},
    lower_bounds::Vector{<:Real},
    upper_bounds::Vector{<:Real},
    method::Type{NewMultiEvalMethod};
    num_supports::Int = InfiniteOpt.DefaultNumSupports,
    weight_func::Function
    )::InfiniteOpt.AbstractMeasureData # REPLACE WITH ACTUAL ARGUMENTS
    # REPLACE WITH ACTUAL FUNCTIONALITY
    # ADD CHECKS IF NECESSARY
    if !all(i -> (dispatch_variable_ref(i) isa InfiniteOpt.IndependentParameterRef), prefs)
        @warn("The method is implemented for independent multivariate parameters.")
    end
    supports_dict = Dict()
    for i in eachindex(lower_bounds)
        increment = (upper_bounds[i] - lower_bounds[i]) / (num_supports - 1)
        supports_dict[i] = [lower_bounds[i] + (j - 1) * increment for j in 1:num_supports]
    end
    supports = Array{Vector{Float64}, 1}(undef, num_supports)
    for j in 1:num_supports
        supports[j] = [supports_dict[k][j] for k in keys(lower_bounds)]
    end
    coeffs = ones(num_supports) / num_supports * prod(upper_bounds .- lower_bounds)
    # MAKE SURE IT RETURNS APPROPRIATE MEASURE DATA
    return InfiniteOpt.DiscreteMeasureData(
        prefs, coeffs, supports,
        weight_function = weight_func,
        lower_bounds = lower_bounds,
        upper_bounds = upper_bounds)
end


text:
struct ParameterRef
    model
    index::Int
    subindex::Any
    subindex_names::Any
end

struct SectorRef
    model
    index::Int
    subindex::Any
    subindex_names::Any
end

struct CommodityRef
    model
    index::Int
    subindex::Any
    subindex_names::Any
end

struct ConsumerRef
    model
    index::Int
    subindex::Any
    subindex_names::Any
end

abstract type Parameter end;

mutable struct ScalarParameter <: Parameter
    name::Symbol
    value::Float64
    description::String

    function ScalarParameter(name::Symbol; value::Float64=1., description::AbstractString="")
        return new(name, value, description)
    end
end

mutable struct IndexedParameter <: Parameter
    name::Symbol
    indices::Any
    value::DenseAxisArray
    description::String

    function IndexedParameter(name::Symbol, indices; value::Float64=1., description::AbstractString="")
        return new(name, indices, DenseAxisArray(fill(value, length.(indices)...), indices...), description)
    end
end

function Parameter(name; indices=nothing, kwargs...)
    return indices===nothing ? ScalarParameter(name; kwargs...) : IndexedParameter(name, indices; kwargs...)
end

abstract type Sector end;

mutable struct ScalarSector <: Sector
    name::Symbol
    benchmark::Float64
    description::String
    fixed::Bool

    function ScalarSector(name::Symbol; description::AbstractString="", benchmark::Float64=1., fixed=false)
        return new(name, benchmark, description, fixed)
    end
end

mutable struct IndexedSector <: Sector
    name::Symbol
    indices::Any
    benchmark::DenseAxisArray
    description::String
    fixed::DenseAxisArray

    function IndexedSector(name::Symbol, indices; description::AbstractString="", benchmark::Float64=1., fixed=false)
        return new(name, indices, DenseAxisArray(fill(benchmark, length.(indices)...), indices...), description, DenseAxisArray(fill(fixed, length.(indices)...), indices...))
    end
end

function Sector(name; indices=nothing, kwargs...)
    return indices===nothing ? ScalarSector(name; kwargs...) : IndexedSector(name, indices; kwargs...)
end

abstract type Commodity end;

mutable struct ScalarCommodity <: Commodity
    name::Symbol
    benchmark::Float64
    description::String
    fixed::Bool

    function ScalarCommodity(name::Symbol; description::AbstractString="", benchmark::Float64=1., fixed=false)
        return new(name, benchmark, description, fixed)
    end
end

mutable struct IndexedCommodity <: Commodity
    name::Symbol
    indices::Any
    benchmark::DenseAxisArray
    description::String
    fixed::DenseAxisArray

    function IndexedCommodity(name::Symbol, indices; description::AbstractString="", benchmark::Float64=1., fixed=false)
        return new(name, indices, DenseAxisArray(fill(benchmark, length.(indices)...), indices...), description, DenseAxisArray(fill(fixed, length.(indices)...), indices...))
    end
end

function Commodity(name; indices=nothing, kwargs...)
    return indices===nothing ? ScalarCommodity(name; kwargs...) : IndexedCommodity(name, indices; kwargs...)
end

abstract type Consumer end;

mutable struct ScalarConsumer <: Consumer
    name::Symbol
    benchmark::Float64
    description::String
    fixed::Bool

    function ScalarConsumer(name::Symbol; description::AbstractString="", benchmark::Float64=1., fixed=false)
        return new(name, benchmark, description, fixed)
    end
end

mutable struct IndexedConsumer <: Consumer
    name::Symbol
    indices::Any
    benchmark::DenseAxisArray
    description::String
    fixed::DenseAxisArray

    function IndexedConsumer(name::Symbol, indices; description::AbstractString="", benchmark::Float64=1., fixed=false)
        return new(name, indices, DenseAxisArray(fill(benchmark, length.(indices)...), indices...), description, DenseAxisArray(fill(fixed, length.(indices)...), indices...))
    end
end

function Consumer(name; indices=nothing, kwargs...)
    return indices===nothing ? ScalarConsumer(name; kwargs...) : IndexedConsumer(name, indices; kwargs...)
end

mutable struct Input
    commodity::CommodityRef
    quantity::Union{Float64,Expr}
    production_function::Any

    function Input(commodity::CommodityRef, quantity::Union{Float64,Expr})
        return new(commodity, quantity, nothing)
    end
end

mutable struct Output
    commodity::CommodityRef
    quantity::Union{Float64,Expr}
    production_function::Any

    function Output(commodity::CommodityRef, quantity::Union{Float64,Expr})
        return new(commodity, quantity, nothing)
    end
end

struct Production
    sector::SectorRef
    elasticity::Union{Float64,Expr}
    outputs::Vector{Output}
    inputs::Vector{Input}

    function Production(sector::SectorRef, elasticity::Union{Float64,Expr}, outputs::Vector{Output}, inputs::Vector{Input})
        x = new(sector, elasticity, outputs, inputs)

        for output in outputs
            output.production_function = x
        end

        for input in inputs
            input.production_function = x
        end

        return x
    end
end

struct Endowment
    commodity::CommodityRef
    quantity::Union{Float64,Expr}
end

mutable struct Demand
    commodity::CommodityRef
    quantity::Union{Float64,Expr}
    demand_function::Any

    function Demand(commodity::CommodityRef, quantity::Union{Float64,Expr})
        return new(commodity, quantity, nothing)
    end
end

struct DemandFunction
    consumer::ConsumerRef
    demands::Vector{Demand}
    endowments::Vector{Endowment}

    function DemandFunction(consumer::ConsumerRef, demands::Vector{Demand}, endowments::Vector{Endowment})
        x = new(consumer, demands, endowments)

        for demand in demands
            demand.demand_function = x
        end

        return x
    end
end

mutable struct Model
    _parameters::Vector{Parameter}
    _sectors::Vector{Sector}
    _commodities::Vector{Commodity}
    _consumers::Vector{Consumer}

    _productions::Vector{Production}
    _demands::Vector{DemandFunction}

    _jump_model::Union{Nothing,JuMP.Model}
    _status

    _nlexpressions::Vector{Any}

    function Model()
        return new(
            Parameter[],
            Sector[],
            Commodity[],
            Consumer[],
            Production[],
            DemandFunction[],
            nothing,
            nothing,
            []
        )
    end
end

function Base.show(io::IO, m::Model)
    println(io, "MPSGE model with $(length(m._sectors)) sectors, $(length(m._commodities)) commodities and $(length(m._consumers)) consumers.")

    if length(m._sectors) > 0
        print(io, "  Sectors: ")
        print(io, join(["$(s.name) (bm=$(s.benchmark))" for s in m._sectors], ", "))
        println(io)
    end

    if length(m._commodities) > 0
        print(io, "  Commodities: ")
        print(io, join(["$(c.name) (bm=$(c.benchmark))" for c in m._commodities], ", "))
        println(io)
    end

    if length(m._consumers) > 0
        print(io, "  Consumers: ")
        print(io, join(["$(c.name) (bm=$(c.benchmark))" for c in m._consumers], ", "))
        println(io)
    end

    if m._jump_model!==nothing
        if m._status==:Solved
            println(io, "Solution:")

            for n in JuMP.all_variables(m._jump_model)
                println(io, "  $n:\t$(Complementarity.result_value(n))")
            end
        else
            println(io, "Did not solve with error: $(m._status).")
        end
    end
end

function get_name(sector::SectorRef, include_subindex=false)
    if sector.subindex===nothing || include_subindex==false
        return sector.model._sectors[sector.index].name
    else
        return Symbol("$(sector.model._sectors[sector.index].name )[$(sector.subindex_names)]")
    end
end

function get_name(commodity::CommodityRef, include_subindex=false)
    if commodity.subindex===nothing || include_subindex==false
        return commodity.model._commodities[commodity.index].name
    else
        return Symbol("$(commodity.model._commodities[commodity.index].name )[$(commodity.subindex_names)]")
    end
end

function get_name(consumer::ConsumerRef, include_subindex=false)
    if consumer.subindex===nothing || include_subindex===false
        return consumer.model._consumers[consumer.index].name
    else
        return Symbol("$(consumer.model._consumers[consumer.index].name )[$(consumer.subindex_names)]")
    end
end

function get_full(s::SectorRef)
    return s.model._sectors[s.index]
end

function get_full(c::CommodityRef)
    return c.model._commodities[c.index]
end

# Outer constructors

function Input(commodity::CommodityRef, quantity::Number)
    return Input(commodity, convert(Float64, quantity))
end

function Output(commodity::CommodityRef, quantity::Number)
    return Output(commodity, convert(Float64, quantity))
end

function Production(sector::SectorRef, elasticity::Union{Number,Expr}, outputs::Vector{Output}, inputs::Vector{Input})

    if isa(elasticity,Number)
        elasticity = convert(Float64, elasticity)
    end

    return Production(sector, elasticity, outputs, inputs)
end

function Endowment(commodity::CommodityRef, quantity::Number)
    return Endowment(commodity, convert(Float64, quantity))
end

function Demand(commodity::CommodityRef, quantity::Number)
    return Demand(commodity, convert(Float64, quantity))
end

function add!(m::Model, s::ScalarSector)
    m._jump_model = nothing
    push!(m._sectors, s)
    return SectorRef(m, length(m._sectors), nothing, nothing)
end

function add!(m::Model, s::IndexedSector)
    m._jump_model = nothing
    push!(m._sectors, s)

    temp_array = Array{SectorRef}(undef, length.(s.indices)...)

    for i in CartesianIndices(temp_array)
        # TODO Fix the [1] thing here to properly work with n-dimensional data
        temp_array[i] = SectorRef(m, length(m._sectors), i, string(s.indices[1][i]))
    end
    return JuMP.Containers.DenseAxisArray(temp_array, s.indices...)
end

function add!(m::Model, c::ScalarCommodity)
    m._jump_model = nothing
    push!(m._commodities, c)
    return CommodityRef(m, length(m._commodities), nothing, nothing)
end

function add!(m::Model, c::IndexedCommodity)
    m._jump_model = nothing
    push!(m._commodities, c)

    temp_array = Array{CommodityRef}(undef, length.(c.indices)...)

    for i in CartesianIndices(temp_array)
        # TODO Fix the [1] thing here to properly work with n-dimensional data
        temp_array[i] = CommodityRef(m, length(m._commodities), i, string(c.indices[1][i]))
    end
    return JuMP.Containers.DenseAxisArray(temp_array, c.indices...)
end

function add!(m::Model, cn::ScalarConsumer)
    m._jump_model = nothing
    push!(m._consumers, cn)
    return ConsumerRef(m, length(m._consumers), nothing, nothing)
end

function add!(m::Model, cn::IndexedConsumer)
    m._jump_model = nothing
    push!(m._consumers, cn)

    temp_array = Array{ConsumerRef}(undef, length.(cn.indices)...)

    for i in CartesianIndices(temp_array)
        # TODO Fix the [1] thing here to properly work with n-dimensional data
        temp_array[i] = ConsumerRef(m, length(m._consumers), i, string(cn.indices[1][i]))
    end
    return JuMP.Containers.DenseAxisArray(temp_array, cn.indices...)
end

function add!(m::Model, p::Production)
    m._jump_model = nothing
    push!(m._productions, p)
    return m
end

function add!(m::Model, c::DemandFunction)
    m._jump_model = nothing
    push!(m._demands, c)
    return m
end

function add!(m::Model, p::ScalarParameter)
    m._jump_model = nothing
    push!(m._parameters, p)
    return ParameterRef(m, length(m._parameters), nothing, nothing)
end

function add!(m::Model, p::IndexedParameter)
    m._jump_model = nothing
    push!(m._parameters, p)

    temp_array = Array{ParameterRef}(undef, length.(p.indices)...)

    for i in CartesianIndices(temp_array)
        # TODO Fix the [1] thing here to properly work with n-dimensional data
        temp_array[i] = ParameterRef(m, length(m._parameters), i, string(p.indices[1][i]))
    end
    return JuMP.Containers.DenseAxisArray(temp_array, p.indices...)
end

function JuMP.value(m::Model, name::Symbol)
    Complementarity.result_value(m._jump_model[name])
end

function solve!(m::Model; solver::Symbol=:PATH, kwargs...)
    if m._jump_model===nothing
        m._jump_model = build(m)
    end

    set_all_start_values(m)
    set_all_parameters(m)
    set_all_bounds(m)

    m._status = Complementarity.solveMCP(m._jump_model; solver=solver, kwargs...)

    return m
end

function JuMP.set_value(p::ParameterRef, new_value::Float64)
    p.model._parameters[p.index].value = new_value
end

function JuMP.set_value(c::ConsumerRef, new_value::Float64)
    c.model._consumers[c.index].benchmark = new_value
end

function JuMP.set_value(c::CommodityRef, new_value::Float64)
    c.model._commodities[c.index].benchmark = new_value
end

function set_fixed!(commodity::CommodityRef, new_value::Bool)
    c = commodity.model._commodities[commodity.index]
    if c isa ScalarCommodity
        c.fixed = new_value
    else
        c.fixed[commodity.subindex] = new_value
    end
    return nothing
end

function set_fixed!(consumer::ConsumerRef, new_value::Bool)
    c = consumer.model._consumers[consumer.index]
    if c isa ScalarConsumer
        c.fixed = new_value
    else
        c.fixed[consumer.subindex] = new_value
    end
    return nothing
end

function set_fixed!(sector::SectorRef, new_value::Bool)
    s = sector.model._sectors[sector.index]
    if s isa ScalarSector
        s.fixed = new_value
    else
        s.fixed[sector.subindex] = new_value
    end
    return nothing
end


text:
workspace()
using MCTSParser


text:
templates_folder() = joinpath(homedir(), ".azmanagers")

function save_template(templates_filename::AbstractString, name::AbstractString, template::Dict)
    templates = isfile(templates_filename) ? JSON.parse(read(templates_filename, String)) : Dict{String,Any}()
    templates[name] = template
    if !ispath(templates_folder())
        mkdir(templates_folder())
    end
    write(templates_filename, json(templates, 1))
    nothing
end

#
# scale-set templates
#
"""
    AzManagers.build_sstemplate(name; kwargs...)

returns a dictionary that is an Azure scaleset template for use in `addprocs` or for saving
to the `~/.azmanagers` folder.

# required key-word arguments
* `subscriptionid` Azure subscription
* `admin_username` ssh user for the scaleset virtual machines
* `location` Azure data-center location
* `resourcegroup` Azure resource-group
* `imagegallery` Azure image gallery that contains the VM image
* `imagename` Azure image
* `vnet` Azure virtual network for the scaleset
* `subnet` Azure virtual subnet for the scaleset
* `skuname` Azure VM type

# optional key-word arguments
* `subscriptionid_image` Azure subscription corresponding to the image gallery, defaults to `subscriptionid`
* `resourcegroup_vnet` Azure resource group corresponding to the virtual network, defaults to `resourcegroup`
* `resourcegroup_image` Azure resource group correcsponding to the image gallery, defaults to `resourcegroup`
* `skutier = "Standard"` Azure SKU tier.
* `datadisks=[]` list of data disks to create and attach [1]
* `tempdisk = "sudo mkdir -m 777 /mnt/scratch\nln -s /mnt/scratch /scratch"` cloud-init commands used to mount or link to temporary disk

# Notes
[1] Each datadisk is a Dictionary. For example,
```julia
Dict("createOption"=>"Empty", "diskSizeGB"=>1023, "managedDisk"=>Dict("storageAccountType"=>"PremiumSSD_LRS"))
```
or, to accept the defaults,
```julia
Dict("diskSizeGB"=>1023)
```
The above example is populated with the default options.  So, if `datadisks=[Dict()]`, then the default options
will be included.
"""
function build_sstemplate(name;
        subscriptionid,
        subscriptionid_image = "",
        admin_username,
        location,
        resourcegroup,
        resourcegroup_vnet = "",
        resourcegroup_image = "",
        imagegallery,
        imagename,
        vnet,
        subnet,
        skutier="Standard",
        datadisks=[],
        tempdisk="sudo mkdir -m 777 /mnt/scratch\nln -s /mnt/scratch /scratch",
        skuname)
    resourcegroup_vnet == "" && (resourcegroup_vnet = resourcegroup)
    resourcegroup_image == "" && (resourcegroup_image = resourcegroup)
    subscriptionid_image == "" && (subscriptionid_image = subscriptionid)
    subnetid = "/subscriptions/$subscriptionid/resourceGroups/$resourcegroup_vnet/providers/Microsoft.Network/virtualNetworks/$vnet/subnets/$subnet"
    image = "/subscriptions/$subscriptionid_image/resourceGroups/$resourcegroup_image/providers/Microsoft.Compute/galleries/$imagegallery/images/$imagename"

    _datadisks = Dict{String,Any}[]
    ultrassdenabled = false
    for (idisk,datadisk) in enumerate(datadisks)
        datadisk_template = Dict{String,Any}("createOption"=>"Empty", "diskSizeGB"=>1023, "lun"=>1, "managedDisk"=>Dict("storageAccountType"=>"Premium_LRS"))
        merge!(datadisk_template, datadisk)
        merge!(datadisk_template, Dict("lun"=>idisk))
        push!(_datadisks, datadisk_template)
        if datadisk_template["managedDisk"]["storageAccountType"] ∈ ("UltraSSD_LRS",)
            ultrassdenabled = true
        end
    end

    Dict(
        "tempdisk" => tempdisk,
        "value" => Dict(
            "sku" => Dict(
                "tier" => skutier,
                "capacity" => 2,
                "name" => skuname
            ),
            "location" => location,
            "properties" => Dict(
                "overprovision" => true,
                "singlePlacementGroup" => false,
                "additionalCapabilities" => Dict(
                    "ultraSSDEnabled" => ultrassdenabled
                ),
                "virtualMachineProfile" => Dict{String,Any}(
                    "storageProfile" => Dict(
                        "imageReference" => Dict(
                            "id" => image
                        ),
                        "osDisk" => Dict(
                            "caching" => "ReadWrite",
                            "managedDisk" => Dict(
                                "storageAccountType" => "Standard_LRS"
                            ),
                            "createOption" => "FromImage"
                        ),
                        "dataDisks" => _datadisks
                    ),
                    "osProfile" => Dict(
                        "computerNamePrefix" => name,
                        "adminUsername" => admin_username,
                        "linuxConfiguration" => Dict(
                            "ssh" => Dict(
                                "publicKeys" => []
                                ),
                            "disablePasswordAuthentication" => true
                        )
                    ),
                    "networkProfile" => Dict(
                        "networkInterfaceConfigurations" => [
                            Dict(
                                "name" => name,
                                "properties" => Dict(
                                    "primary" => true,
                                    "ipConfigurations" => [
                                        Dict(
                                            "name" => name,
                                            "properties" => Dict(
                                                "subnet" => Dict(
                                                    "id" => subnetid
                                                )
                                            )
                                        )
                                    ] # ipConfigurations
                                )
                            )
                        ] # networkInterfaceConfigurations
                    ), # networkProfile
                ), # virtualMachineProfile
                "upgradePolicy" => Dict(
                    "mode" => "Manual"
                )
            ) # properties
        )
    )
end

templates_filename_scaleset() = joinpath(templates_folder(), "templates_scaleset.json")

"""
    AzManagers.save_template_scaleset(scalesetname, template)

Save `template::Dict` generated by AzManagers.build_sstemplate to $(templates_filename_scaleset()).
"""
save_template_scaleset(name::AbstractString, template::Dict) = save_template(templates_filename_scaleset(), name, template)

#
# NIC templates
#
"""
    AzManagers.build_nictemplate(nic_name; kwargs...)

Returns a dictionary for a NIC template, and that can be passed to the `addproc` method, or written
to AzManagers.jl configuration files.

# Required keyword arguments
* `subscriptionid` Azure subscription
* `resourcegroup_vnet` Azure resource group that holds the virtual network that the NIC is attaching to.
* `vnet` Azure virtual network for the NIC to attach to.
* `subnet` Azure sub-network name.
* `location` location of the Azure data center where the NIC correspond to.

# Optional keyword arguments
* `accelerated=true` use accelerated networking (not all VM sizes support accelerated networking).
"""
function build_nictemplate(name;
        subscriptionid,
        resourcegroup_vnet,
        vnet,
        subnet,
        accelerated = true,
        location)
    subnetid = "/subscriptions/$subscriptionid/resourceGroups/$resourcegroup_vnet/providers/Microsoft.Network/virtualNetworks/$vnet/subnets/$subnet"

    body = Dict(
        "properties" => Dict(
            "enableAcceleratedNetworking" => accelerated,
            "ipConfigurations" => [
                Dict(
                    "name" => "ipConfig1",
                    "properties" => Dict(
                        "subnet" => Dict(
                            "id" => subnetid
                        )
                    )
                )
            ]
        ),
        "location" => location
    )
end

templates_filename_nic() = joinpath(templates_folder(), "templates_nic.json")

"""
    AzManagers.save_template_nic(nic_name, template)

Save `template::Dict` generated by AzManagers.build_nictmplate to $(templates_filename_nic()).
"""
save_template_nic(name::AbstractString, template::Dict) = save_template(templates_filename_nic(), name, template)

#
# VM templates
#
"""
    AzManagers.build_vmtemplate(vm_name; kwargs...)

Returns a dictionary for a virtual machine template, and that can be passed to the `addproc` method
or written to AzManagers.jl configuration files.

# Required keyword arguments
* `subscriptionid` Azure subscription
* `admin_username` ssh user for the scaleset virtual machines
* `location` Azure data center location
* `resourcegroup` Azure resource group where the VM will reside
* `imagegallery` Azure shared image gallery name
* `imagename` Azure image name that is in the shared image gallery
* `vmsize` Azure vm type, e.g. "Standard_D8s_v3"

# Optional keyword arguments
* `resourcegroup_vnet` Azure resource group containing the virtual network, defaults to `resourcegroup`
* `subscriptionid_image` Azure subscription containing the image gallery, defaults to `subscriptionid`
* `resourcegroup_image` Azure resource group containing the image gallery, defaults to `subscriptionid`
* `nicname = "cbox-nic"` Name of the NIC for this VM
* `datadisks=[]` additional data disks to attach
* `tempdisk = "sudo mkdir -m 777 /mnt/scratch\nln -s /mnt/scratch /scratch"`  cloud-init commands used to mount or link to temporary disk

# Notes
[1] Each datadisk is a Dictionary. For example,
```julia
Dict("createOption"=>"Empty", "diskSizeGB"=>1023, "managedDisk"=>Dict("storageAccountType"=>"PremiumSSD_LRS"))
```
The above example is populated with the default options.  So, if `datadisks=[Dict()]`, then the default options
will be included.
"""
function build_vmtemplate(name;
        subscriptionid,
        admin_username,
        subscriptionid_image = "",
        location,
        resourcegroup,
        resourcegroup_vnet = "",
        resourcegroup_image = "",
        imagegallery,
        imagename,
        vmsize,
        datadisks = [],
        tempdisk = "sudo mkdir -m 777 /mnt/scratch\nln -s /mnt/scratch /scratch",
        nicname = "cbox-nic")
    resourcegroup_vnet == "" && (resourcegroup_vnet = resourcegroup)
    resourcegroup_image == "" && (resourcegroup_image = resourcegroup)
    subscriptionid_image == "" && (subscriptionid_image = subscriptionid)

    image = "/subscriptions/$subscriptionid_image/resourceGroups/$resourcegroup_image/providers/Microsoft.Compute/galleries/$imagegallery/images/$imagename"

    ultrassdenabled = false
    _datadisks = Dict{String,Any}[]
    for (idisk,datadisk) in enumerate(datadisks)
        datadisk_template = Dict{String,Any}("createOption"=>"Empty", "diskSizeGB"=>1023, "lun"=>1, "managedDisk"=>Dict("storageAccountType"=>"Premium_LRS"))
        merge!(datadisk_template, datadisk)
        merge!(datadisk_template, Dict("lun"=>idisk, "name"=>"scratch$idisk"))
        push!(_datadisks, datadisk_template)
        if datadisk_template["managedDisk"]["storageAccountType"] ∈ ("UltraSSD_LRS",)
            ultrassdenabled = true
        end
    end

    Dict(
        "tempdisk" => tempdisk,
        "value" => Dict(
            "location" => location,
            "properties" => Dict(
                "additionalCapabilities" => Dict(
                    "ultraSSDEnabled"=>ultrassdenabled
                ),
                "hardwareProfile" => Dict(
                    "vmSize" => vmsize
                ),
                "storageProfile" => Dict(
                    "imageReference" => Dict(
                        "id" => image
                    ),
                    "osDisk" => Dict(
                        "caching" => "ReadWrite",
                        "managedDisk" => Dict(
                            "storageAccountType" => "Standard_LRS"
                        ),
                        "createOption" => "FromImage"
                    ),
                    "dataDisks" => _datadisks
                ),
                "osProfile" => Dict(
                    "computerName" => name,
                    "adminUsername" => admin_username,
                    "linuxConfiguration" => Dict(
                        "ssh" => Dict(
                            "publicKeys" => []
                        ),
                        "disablePasswordAuthentication" => true
                    )
                ),
                "networkProfile" => Dict(
                    "networkInterfaces" => [
                        Dict(
                            "id" => "/subscription/$subscriptionid/resourceGroups/$resourcegroup_vnet/providers/Microsoft.Network/networkInterfaces/$nicname",
                            "properties" => Dict(
                                "primary" => true
                            )
                        )
                    ]
                )
            )
        )
    )
end

templates_filename_vm() = joinpath(templates_folder(), "templates_vm.json")

"""
    AzManagers.save_template_vm(vm_name, template)

Save `template::Dict` generated by AzManagers.build_vmtmplate to $(templates_filename_vm()).
"""
save_template_vm(name::AbstractString, template::Dict) = save_template(templates_filename_vm(), name, template)


text:
# This file is a part of BinaryCommitteeMachineFBP.jl. License is MIT: http://github.com/carlobaldassi/BinaryCommitteeMachineFBP.jl/LICENCE.md

module AtanhErf

export atanherf, batanherf

using StatsFuns
using SpecialFunctions
using JLD2
using Interpolations

const builddir = joinpath(dirname(@__FILE__), "..", "deps", "builds")
batanherf(x::Float64) = Float64(atanh(erf(big(x))))

let
    mm = 16.0
    st = 1e-4
    r = 1.0:st:mm
    rb = big.(first(r):step(r):last(r))

    interp_degree = Quadratic
    interp_boundary = Line
    interp_grid = OnGrid
    # interp_type = Interpolations.BSplineInterpolation{Float64,1,Vector{Float64},
    #                                                   BSpline{interp_degree{interp_boundary}},
    #                                                   OnGrid,0}

    interp_type = Interpolations.BSplineInterpolation{Float64,1,Array{Float64,1},
                                                      BSpline{interp_degree{interp_boundary{interp_grid}}},
                                                      Tuple{Base.Slice{UnitRange{Int64}}}}
    function getinp!()
        isdir(builddir) || mkdir(builddir)
        filename = joinpath(builddir, "atanherf_interp.max_$mm.step_$st.jld")
        if isfile(filename)
            @load filename inp
        else
            @info("Computing atanh(erf(x)) table, this may take a while...")
            inp = setprecision(BigFloat, 512) do
                interpolate!(Float64[atanh(erf(x)) for x in rb], BSpline(interp_degree(interp_boundary(interp_grid()))))
            end
            jldopen(filename, "w") do f
                write(f, "inp", inp)
            end
        end
        # @show interp_type
        # @show typeof(inp)
        # x = 0.0
        # @show inp[(x - first(r)) / step(r) + 1]::Float64
        return inp::interp_type
    end

    inp = getinp!()

    global atanherf_interp
    atanherf_interp(x::Float64) = inp((x - first(r)) / step(r) + 1)::Float64
end

function atanherf_largex(x::Float64)
    x² = x^2
    t = 1/x²

    return sign(x) * (2log(abs(x)) + log4π + 2x² +
                      t * @evalpoly(t,
                                    1,
                                    -1.25,
                                    3.0833333333333335,
                                    -11.03125,
                                    51.0125,
                                    -287.5260416666667,
                                    1906.689732142857,
                                    -14527.3759765625,
                                    125008.12543402778,
                                    -1.1990066259765625e6)
                     ) / 4
end

function atanherf(x::Float64)
    ax = abs(x)
    ax ≤ 2 && return atanh(erf(x))
    ax ≤ 15 && return sign(x) * atanherf_interp(ax)
    return atanherf_largex(x)
end

end


text:
using CSV
include("oeid.jl");
include("idanalysis.jl");

AoP, AoQ = importPQ("invivo_data/Ao_baseline1_Ao",4);
checksampling(["AoP", "AoQ"], AoP, AoQ)
pltPQ = plotPQ(AoP, AoQ, "in vivo")

h = 0.005 # sampling rate

# PV loops
title="in vivo"
for i = 1:4
    local pltPQ = plot(plot(AoP[i].timestamp, AoP[i].vb_va, xticks=:none, ylabel = "AoQ [L/min]"), plot(AoQ[i].timestamp, AoQ[i].vb_va, ylabel = "AoQ [L/min]"), xlims = (AoQ[i][1,:timestamp]|>Dates.value, AoQ[i][end,:timestamp]|>Dates.value), layout = (2,1));
    V, pltPV = pvloop(AoP[i].vb_va,AoQ[i].vb_va,h)

    display(plot(pltPV, pltPQ,  layout = (1, 2), size = (1200,400), legend=:none, title = [title*" $(i) at ~70 bpm" "" ""]))
    # savefig("./results/pv_invivo$(i)_" * Dates.format(now(), "yyyymmdd") * ".png")
end

# θ=[Rp,C,Rc,L]
u = vcat([AoQ[i].vb_va for i in 1:4]...);
y = vcat([AoP[i].vb_va for i in 1:4]...);
Rp = mean(y)/mean(u);
θ0 = [Rp, 1, 1, 1]
x00 = solvex0(G(θ0,h), u)
θid = [1,1,1,1]
x0id = false
x0solve = true
idinvivo = IdConditions(θ0,x00,θid,x0id,x0solve)

# id first attempt
@time phat, res = identify(θ0,x00,u,y,h,θid=θid,x0id=x0id,x0solve=x0solve);
θhat, x0hat = deal(phat, idinvivo)
θ0, x00
x0hatsolve = solvex0(G(θhat,h),u)
hess = Optim.trace(res)[end].metadata["h(x)"]
svd(hess).U
yhat, xhat, t, plt = plot_result(θhat,x0hatsolve,y,u,h)
plt

# multiple init
θ0s, θhats, x0hats, res = id_multiple_init(u, y, x00, θ0, 100, θid=θid,x0id=x0id,x0solve=x0solve, paramdist=[false,true,true,true]);

results4, nanratio = multiinitresults(θhats, x0hats, res);
# equal MSE, so parameters set with more near-physiological L and x1 chosen
results4[99,:hess]
results4[99,:hessSVD].U
bestinvivo4 = results4[99,[1:7...,11,12]]

θbest4 = Array(results4[99, 1:length(θ0)]);
[θ0, round.(θbest4, sigdigits=3)]
x0best4 = Array(results4[99,["x01", "x02"]]);
[x00, round.(x0best4, sigdigits=3)]
yhat4, xhat4, t, plt_best = plot_result(θbest4, x0best4, y, u, h)
plot(plt_best, size=(1200,400))

# savefig("results/best_fit_invivo4_$(Dates.format(now(), "yyyymmdd")).png")
# CSV.write("results/results_invivo4_$(Dates.format(now(), "yyyymmdd")).csv",results4)
plot(t, xhat4, title = "exvivo $(round.(θbest4, sigdigits=3))", xlabel = "Time [s]", label=["x1 (pk-pk ≈ $(round(maximum(xhat4[1:100,1])-minimum(xhat4[1:100,1]), sigdigits=2)))" "x2"], size=(1200,400))

# plot parameter histograms before and after fit
θlabels = ["Rp" "C" "Rc" "L"]
fithistograms(θ0s,θhats,θlabels)
# savefig("./results/param_conv_" * Dates.format(now(), "yyyymmdd-HHMMSS") * ".png")

w = exp10.(range(-3, stop=3, length=500))
mag4, phase4, w = bode(G(θbest4,h,discretize=false),w)
magplt = plot(w,mag4[:],xaxis=:log, yaxis=:log, label="WK4",title="in vivo");
phaseplt = plot(w,phase4[:], xaxis=:log, label="WK4");
plot(magplt, phaseplt, layout=(2,1))

# multiple initialization, identify 3 element windkessel model
θ0 = copy(θbest4[1:3])
x00 = solvex0(G3(θ0,h), u)
θid = [1,1,1]
x0id = false
x0solve = true
θ0s, θhats, x0hats, res = id_multiple_init(u, y, x00, θ0, 100, θid=θid,x0id=x0id,x0solve=x0solve, G_=G3, paramdist=[false,true,true]);

results3, nanratio = multiinitresults(θhats, x0hats, res);
θbest3 = Array(results3[1,1:3])
x0best3 = [results3[1,"x01"]]
results3[1,:hess]
results3[1,:hessSVD].U
bestinvivo3 = results3[1,["Rp", "C", "Rc", "x01", "MSE","hesscond"]]
# CSV.write("results/results_invivo3_$(Dates.format(now(), "yyyymmdd")).csv",results3)
# results3 = CSV.read("results/results_invivo3_20210308.csv", DataFrame)

yhat3, = lsim(G3(θbest3,h),u,h.*(1:length(u)),x0best3)
mag3, phase3, w = bode(G3(results3[1,1:3],h,discretize=false),w)
mag3plt = plot(w,mag3[:],xaxis=:log, yaxis=:log, label="WK3",title="in vivo");
phase3plt = plot(w,phase3[:], xaxis=:log, label="WK3");
plot(mag3plt, phase3plt, layout=(2,1))

# multiple initialization, identify 2 element windkessel model
θ0 = copy(θbest4[1:2])
x00 = solvex0(G2(θ0,h), u)
θid = [1,1]
θ0s, θhats, x0hats, res = id_multiple_init(u, y, x00, θ0, 100, θid=θid,x0id=x0id,x0solve=x0solve, G_=G2, paramdist=[false,true]);

results2, nanratio = multiinitresults(θhats, x0hats, res);
θbest2 = Array(results2[1,1:2])
x0best2 = [results2[1,"x01"]]
results2[1,:hess]
results2[1,:hessSVD].U
bestinvivo2 = results2[1,["Rp", "C", "x01", "MSE","hesscond"]]
# CSV.write("results/results_invivo2_$(Dates.format(now(), "yyyymmdd")).csv",results2)
# results2 = CSV.read("results/results_invivo2_20210308", DataFrame)

yhat2, = lsim(G2(θbest2,h),u,h.*(1:length(u)),x0best2)
mag2, phase2, w = bode(G2(θbest2,h,discretize=false),w)
mag2plt = plot(w,mag2[:],xaxis=:log, yaxis=:log, label="WK2",title="in vivo");
phase2plt = plot(w,phase2[:], xaxis=:log, label="WK2");
plot(mag2plt, phase2plt, layout=(2,1))

v, = pvloop(y,u,h)
# load results
# results4 = CSV.read("results/results_invivo4_20210309.csv", DataFrame)
# θbest4 = Array(results4[99,1:4])
# x0best4 = solvex0(G(θbest4,h),u)
# yhat4 = simulate(θbest4,x0best,u,h)[:,1]
plot([y,yhat4],xlims=(500,1000))
plot([v[1:300],v[1:300]],[y[1:300],yhat4[1:300]])

# export bode data for pgfplots
comment= *(
    "# fit to data: in vivo Ao_baseline1 20201201\n",
    "# mag4,phase4: 4-element parallel Windkessel model\n",
    "# θhat = [Rc,L,C,Rp]\n",
    "# θhat = $(θbest4)\n",
    "# mag3,phase3: 3-element parallel Windkessel model\n",
    "# θhat = [Rc,C,Rp]\n",
    "# θhat = $(θbest3)\n",
    "# mag2,phase2: 2-element parallel Windkessel model\n",
    "# θhat = [C,Rp]\n",
    "# θhat = $(θbest2)\n")
filename = "bodeinvivo.csv"
open(filename, "w") do f
    write(f, comment)
end
CSV.write(filename, DataFrame(w = w,
mag4 = mag4[:], phase4 = phase4[:],
mag3 = mag3[:], phase3 = phase3[:],
mag2 = mag2[:], phase2 = phase2[:]),header=true,append=true)

# export pressure-flow-volume data for pgfplots
beatstart = 512
beatend = 679
comment = *("# data: in vivo Ao_baseline1 20201201, fourth beat only $(beatstart):$(beatend)\n",
"# t: time in seconds\n",
"# q: measured Aortic flow in L/min\n",
"# p: measured Aortic pressure in mmHg\n",
"# phat4: modelled Aortic pressure in mmHg\n",
"# v: measured volume through the aorta in ml\n",
"# epsilon4: error signal of the model p-phat4\n",
"# 4-element parallel Windkessel model\n",
"# θhat = [Rc,L,C,Rp]\n",
"# θhat = $(θbest4); x0hat = $(x0best4)\n",
"# as well as the corresonding phat and epsilon values for\n",
"# 3-element parallel Windkessel model\n",
"# θhat = [Rc,C,Rp]\n",
"# θhat = $(θbest3); x0hat = $(x0best3)\n",
"# 2-element parallel Windkessel model\n",
"# θhat = [C,Rp]\n",
"# θhat = $(θbest2); x0hat = $(x0best2)\n")
filename = "pqvinvivo.csv"
open(filename, "w") do f
    write(f, comment)
end
CSV.write(filename, DataFrame(t = round.(collect(0:beatend-beatstart)*h, digits=5), q = u[beatstart:beatend], p = y[beatstart:beatend], phat4 = yhat4[beatstart:beatend], phat3 = yhat3[beatstart:beatend], phat2 = yhat2[beatstart:beatend], v = v[beatstart:beatend].-v[beatstart], epsilon4 = y[beatstart:beatend].-yhat4[beatstart:beatend],epsilon3 = y[beatstart:beatend].-yhat3[beatstart:beatend],epsilon2 = y[beatstart:beatend].-yhat2[beatstart:beatend]),header=true,append=true)

max_stars_repo_name: hpigot/windkessel-id

max_issues_repo_name: hpigot/windkessel-id

max_forks_repo_name: hpigot/windkessel-id


text:
using OrdinaryDiffEq, DiffEqDevTools, Test, LinearAlgebra
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear, van

testTol = 0.2

for prob in [prob_ode_linear, prob_ode_2Dlinear]
  sim21 = test_convergence(1 .//2 .^(6:-1:3),prob,RadauIIA5())
  @test sim21.𝒪est[:final] ≈ 5 atol=testTol
end

# test adaptivity
for iip in (true, false)
  vanstiff = ODEProblem{iip}(van, [0;sqrt(3)], (0.0,1.0), 1e6)
  @test length(solve(vanstiff, RadauIIA5())) < 110
  @test length(solve(remake(vanstiff, p=1e7), RadauIIA5())) < 150
  @test length(solve(remake(vanstiff, p=1e7), reltol=[1e-4, 1e-6], RadauIIA5())) < 160
  @test length(solve(remake(vanstiff, p=1e7), RadauIIA5(), reltol=1e-9, abstol=1e-9)) < 860
  @test length(solve(remake(vanstiff, p=1e9), RadauIIA5())) < 160
  @test length(solve(remake(vanstiff, p=1e10), RadauIIA5())) < 190
end


text:
module Exemplo

greet() = print("Hello World!")

include("Novo_arquivo.jl")

export soma

end # module