text:
topoint(x::AbstractVector{Point{N,Float32}}) where {N} = x

# GRRR STUPID SubArray, with eltype different from getindex(x, 1)
topoint(x::SubArray) = topoint([el for el in x])

function topoint(x::AbstractArray{<:Point{N,T}}) where {T,N}
    return topoint(Point{N,Float32}.(x))
end

function topoint(x::AbstractArray{<:Tuple{P,P}}) where {P<:Point}
    return topoint(reinterpret(P, x))
end

function create_shader(scene::Scene, plot::Union{Lines,LineSegments})
    # Potentially per instance attributes
    positions = lift(plot[1], transform_func_obs(plot)) do points, trans
        points = apply_transform(trans, topoint(points))
        if plot isa LineSegments
            return points
        else
            # Repeat every second point to connect the lines !
            return topoint(TupleView{2, 1}(points))
        end
        trans
    end
    startr = lift(p -> 1:2:(length(p) - 1), positions)
    endr = lift(p -> 2:2:length(p), positions)
    p_start_end = lift(positions) do positions
        return (positions[startr[]], positions[endr[]])
    end

    per_instance = Dict{Symbol,Any}(:segment_start => Buffer(lift(first, p_start_end)),
                                    :segment_end => Buffer(lift(last, p_start_end)))
    uniforms = Dict{Symbol,Any}()
    for k in (:linewidth, :color)
        attribute = lift(plot[k]) do x
            x = convert_attribute(x, Key{k}(), key"lines"())
            if plot isa LineSegments
                return x
            else
                # Repeat every second point to connect the lines!
                return isscalar(x) ? x : reinterpret(eltype(x), TupleView{2, 1}(x))
            end
        end
        if isscalar(attribute)
            uniforms[k] = attribute
            uniforms[Symbol("$(k)_start")] = attribute
            uniforms[Symbol("$(k)_end")] = attribute
        else
            if attribute[] isa AbstractVector{<:Number} && haskey(plot, :colorrange)
                attribute = lift(array2color, attribute, plot.colormap, plot.colorrange)
            end
            per_instance[Symbol("$(k)_start")] = Buffer(lift(x -> x[startr[]], attribute))
            per_instance[Symbol("$(k)_end")] = Buffer(lift(x -> x[endr[]], attribute))
        end
    end

    uniforms[:resolution] = scene.camera.resolution
    uniforms[:model] = plot.model
    positions = meta(Point2f0[(0, -1), (0, 1), (1, -1), (1, 1)],
                     uv=Vec2f0[(0, 0), (0, 0), (0, 0), (0, 0)])
    instance = GeometryBasics.Mesh(positions, GLTriangleFace[(1, 2, 3), (2, 4, 3)])
    return InstancedProgram(WebGL(), lasset("line_segments.vert"),
                            lasset("line_segments.frag"), instance,
                            VertexArray(; per_instance...); uniforms...)
end


text:
# Defining calculated columns.

mutable struct DefineNode <: TabularNode
    over::Union{SQLNode, Nothing}
    list::Vector{SQLNode}
    label_map::OrderedDict{Symbol, Int}

    function DefineNode(; over = nothing, list = [], label_map = nothing)
        if label_map !== nothing
            return new(over, list, label_map)
        end
        n = new(over, list, OrderedDict{Symbol, Int}())
        for (i, l) in enumerate(n.list)
            name = label(l)
            if name in keys(n.label_map)
                err = DuplicateLabelError(name, path = [l, n])
                throw(err)
            end
            n.label_map[name] = i
        end
        n
    end
end

DefineNode(list...; over = nothing) =
    DefineNode(over = over, list = SQLNode[list...])

"""
    Define(; over; list = [])
    Define(list...; over)

`Define` adds a column to the output.

# Examples

```jldoctest
julia> person = SQLTable(:person, columns = [:person_id, :birth_datetime]);

julia> q = From(person) |>
           Define(:age => Fun.now() .- Get.birth_datetime) |>
           Where(Get.age .> "16 years");

julia> print(render(q))
SELECT "person_1"."person_id", "person_1"."birth_datetime", (NOW() - "person_1"."birth_datetime") AS "age"
FROM "person" AS "person_1"
WHERE ((NOW() - "person_1"."birth_datetime") > '16 years')
```
"""
Define(args...; kws...) =
    DefineNode(args...; kws...) |> SQLNode

dissect(scr::Symbol, ::typeof(Define), pats::Vector{Any}) =
    dissect(scr, DefineNode, pats)

function PrettyPrinting.quoteof(n::DefineNode, ctx::QuoteContext)
    ex = Expr(:call, nameof(Define), quoteof(n.list, ctx)...)
    if n.over !== nothing
        ex = Expr(:call, :|>, quoteof(n.over, ctx), ex)
    end
    ex
end

label(n::DefineNode) =
    label(n.over)

rebase(n::DefineNode, n′) =
    DefineNode(over = rebase(n.over, n′), list = n.list, label_map = n.label_map)

ext: jl


text:
using SQLite
using Query
using DataFrames
using NamedTuples

db = SQLite.DB(joinpath(Pkg.dir("SQLite"), "test", "Chinook_Sqlite.sqlite"))

result = @from i in table(db, "Employee") begin
         @where i.ReportsTo==2
         @select {Name=i.LastName, Adr=i.Address}
         @collect DataFrame
end

println(result)

result = @from i in table(db, "Employee") begin
         @where i.ReportsTo==2
         @select {Name=i.LastName, Adr=i.Address}
         @collect
end

println(result)

# This is an exmaple where the first part gets executed in the DB
# And the second part uses the Enumerable iterator part of Query.jl

result = @from i in result begin
         @select {Mangled=i.Name * i.Adr}
         @collect DataFrame
end

println(result)

# This examples uses SQLite just like any other DataStreams source.

results = @from i in SQLite.Source(db,"SELECT * FROM Employee") begin
    @where i.ReportsTo==2
    @select {Name=i.LastName, Adr=i.Address}
    @collect DataFrame
end

println(result)


text:
# ============================================================================ #
# Variable
# ============================================================================ #

mutable struct Variable <: AbstractVariable
    name::String
    value::Float64

    function Variable(name, value = 0.0)
        return new(string(name), value)
    end
end

macro variable(varname, value = default_variable_value())
    lhs = esc(varname)
    rhs = :(     Variable( $(QuoteNode(varname)), $value)   )
    ex  = :($(lhs) = $(rhs))

    return ex
end

function evaluate(x::Variable)
    return x.value
end

function evaluate(x::A) where {A <: AbstractArray{<: Variable}}
    s = similar(x)
    for ii in eachindex(x)
        s[ii] = x[ii].value
    end
    return s
end

function to_string(x::Variable)
    return x.name
end

function clone(x::Variable)
    y = Variable(x.name, x.value)
    return y
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using Statistics
using HypothesisTests
using Printf
using Distributed

addprocs(4)

@everywhere using WeightedEnsemble

@everywhere include("doublewell_setup.jl");

@everywhere begin
    Δt = 1e-4;  # time step
    tmax = 100.0;    # terminal time
    nΔt = ceil(Int, tmax/Δt);
end

# number of coarse steps in WE
@everywhere n_we_steps = 1000;
# number of time steps during mutation step
@everywhere nΔt_coarse =nΔt ÷ n_we_steps;
# number of samples in coarse matrix
n_samples_per_bin = 10^2;
# ensemble size
n_particles = 10^2;

# define bin structure
@everywhere voronoi_pts = [[x] for x in LinRange(-1.2,1.2,21)];
@everywhere B₀, bin_id, rebin! = setup_Voronoi_bins(voronoi_pts);

# define the mutation mapping
@everywhere opts = MDOptions(n_iters=nΔt_coarse, n_save_iters = nΔt_coarse)
@everywhere mutation! = x-> sample_trajectory!(x, sampler, options=opts);
@everywhere function mutation(x)
    Xvals, _ = sample_trajectory(x, sampler, options=opts);
    return Xvals[end]
end
selection! = (E, B, t)-> WeightedEnsemble.uniform_allocation!(E, B);

f1 = x-> Float64(-0.1 < x[1] < 0.1); # define observables
f2 = x-> Float64( x[1] > 1.1);
obs = (f1, f2);

# set up ensemble
E₀ = WeightedEnsemble.Dirac_to_EnsembleWithBins(x₀, n_particles);
rebin!(E₀, B₀, 0);

Random.seed!(200)
obs_trajectory = prun_we_observables(E₀, B₀, mutation, selection!, rebin!, n_we_steps, obs);


text:
@testset "main.jl" begin
    @testset "unknown ci config" begin
        struct MockCI <: CompatHelper.CIService end

        @test_throws ErrorException("Unknown CI Config: MockCI") CompatHelper.main(
            ENV, MockCI()
        )
    end

    @testset "no pr GitHub" begin
        tmpdir = mktempdir()

        cd(tmpdir) do
            patches = [
                git_clone_patch,
                project_toml_patch,
                clone_all_registries_patch,
                rm_patch,
                pr_titles_mock,
                git_push_patch,
                gh_pr_patch,
                make_clone_https_patch(tmpdir),
                decode_pkey_patch,
                gh_get_repo_patch,
                cd_patch,
                git_gmb_patch,
                git_checkout_patch,
            ]
            apply(patches) do
                withenv(
                    "GITHUB_REPOSITORY" => "CompatHelper.jl", "GITHUB_TOKEN" => "token"
                ) do
                    prs = CompatHelper.main()
                    @test length(prs) == 0
                end
            end
        end
    end

    @testset "successful run GitHub" begin
        tmpdir = mktempdir()

        cd(tmpdir) do
            patches = [
                git_clone_patch,
                project_toml_patch,
                clone_all_registries_patch,
                rm_patch,
                pr_titles_mock,
                git_push_patch,
                gh_pr_patch,
                make_clone_https_patch(tmpdir),
                decode_pkey_patch,
                gh_get_repo_patch,
                cd_patch,
                git_gmb_patch,
                git_checkout_patch,
                gh_make_pr_patch,
            ]
            apply(patches) do
                withenv(
                    "GITHUB_REPOSITORY" => "CompatHelper.jl", "GITHUB_TOKEN" => "token"
                ) do
                    prs = CompatHelper.main()
                    @test length(prs) == 1
                    @test prs[1] isa GitHub.PullRequest
                end
            end
        end
    end

    @testset "successful run GitLab" begin
        mktempdir() do tmpdir
            cd(tmpdir) do
                patches = [
                    git_clone_patch,
                    project_toml_patch,
                    clone_all_registries_patch,
                    rm_patch,
                    pr_titles_mock,
                    git_push_patch,
                    gl_pr_patch,
                    make_clone_https_patch(tmpdir),
                    decode_pkey_patch,
                    gl_get_repo_patch,
                    cd_patch,
                    gl_make_pr_patch,
                ]

                apply(patches) do
                    withenv(
                        "GITLAB_CI" => "true",
                        "CI_PROJECT_PATH" => "CompatHelper.jl",
                        "GITLAB_TOKEN" => "token",
                        "GITHUB_REPOSITORY" => "false",
                    ) do
                        delete!(ENV, "GITHUB_REPOSITORY")
                        prs = CompatHelper.main()
                        @test length(prs) == 1
                        @test prs[1] isa GitLab.MergeRequest
                    end
                end
            end
        end
    end
end


text:
using GridInterpolations
using LocalFunctionApproximation
using StaticArrays
using JLD2, FileIO
using POMDPs
using Random
using ContinuumWorld
using CMSSPs

# global_logger(SimpleLogger(stderr, Logging.Debug))

rng = MersenneTwister(1234)
const TRIALS = 10

@info "Loading policies"
inhor_fn = "grids-continuum-params1-inhor.jld2"
outhor_fn = "grids-continuum-params1-outhor.jld2"
modal_horizon_policy = load_modal_horizon_policy_localapprox(inhor_fn, outhor_fn)

# Create reqments for solver by loading modal policies and modal MDPs etc
modal_policies = Dict{Int64, ModalFinInfHorPolicy}()
modal_mdps = Dict{Int64, GridsContinuumMDPType}()
for m in CONTINUUM_MODES
    global modal_policies
    global modal_mdps
    modal_policies[m] = ModalFinInfHorPolicy(modal_horizon_policy, nothing)
    modal_mdps[m] = modal_horizon_policy.in_horizon_policy.mdp
end

# Define the params struct and the dmssp problem object
params = modal_horizon_policy.in_horizon_policy.mdp.params
cmssp = create_continuum_cmssp(params, rng)

# Other parameters required by the solver - the replanning period (in terms of timesteps) and the goal modes
deltaT = 5
goal_modes = [CONTINUUM_GOAL_MODE]

for trial = 1:TRIALS

    @info "Generating start state and context"
    start_state = generate_start_state(cmssp, rng)
    set_start_context_set!(cmssp, rng)

    @debug start_state

    @info "Creating solver"
    continuum_solver = GridsContinuumSolverType{typeof(rng)}(params.num_bridge_samples,
                                                modal_policies,
                                                deltaT,
                                                goal_modes,
                                                start_state,
                                                ContinuumDummyBookkeeping(),
                                                rng)

    @info "Now solving"
    solve(continuum_solver, cmssp)
end


text:
using IVPTestSuite  # imports ODE, DASSL and Sundials
const S = Solvers # for convenience

# choose test-case (all of them are in IVPTestSuite.tc_all)
tc_name = [:plei
           :vdpol
           :rober
           :threebody
           :bruss1d
           :chemakzo
           :hires][1]
tc = IVPTestSuite.tc_all[tc_name]

# Pick a solver.  S.allsolvers is the list of all solvers.  Generally
# each family of solvers also has a list.
solver = S.ODEsolvers[ODE.ode4]

###
# make a TestRun which combines a TestCase with a Solver + some extras:
#Example protocol for Fixed step solvers
tsteps = collect(linspace(tc.tspan[1],tc.tspan[2],10^5))
tr = TestRunFixedStep(tc, solver, Dict{Symbol,Any}(), tsteps)
# For adaptive solvers use TestRunAdapt instead. See ex1.jl for such a case.

# run it and display results:
re = IVPTestSuite.run_ode_test_throwerror(tr) # This returns a TestResults instance
re = IVPTestSuite.run_ode_test_throwerror(tr) # Run it twice to get accurate timings
show(re)


text:
# Top level test script for Stan.jl
using StanMCMCChains, Test, Statistics

println("\nRunning tests for StanMCMCChains-j1.2-v5.1.0:\n")

execution_tests = [
  "test_mcmcchains.jl",
  "test_m10.4s.jl",
  "test_m12.6sl.jl"
]

@testset "StanMCMCChains.jl" begin

  for my_test in execution_tests
      println("\n\n  * $(my_test) *\n")
      include(my_test)
  end

end


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

"""
    DeleteRecord()

Deletes a Record from a FeatureGroup. A new record will show up in the OfflineStore when the DeleteRecord API is called. This record will have a value of True in the is_deleted column.

# Required Parameters
- `EventTime`: Timestamp indicating when the deletion event occurred. EventTime can be used to query data at a certain point in time.
- `FeatureGroupName`: The name of the feature group to delete the record from.
- `RecordIdentifierValueAsString`: The value for the RecordIdentifier that uniquely identifies the record, in string format.

"""
delete_record(EventTime, FeatureGroupName, RecordIdentifierValueAsString; aws_config::AWSConfig=global_aws_config()) = sagemaker_featurestore_runtime("DELETE", "/FeatureGroup/$(FeatureGroupName)", Dict{String, Any}("EventTime"=>EventTime, "RecordIdentifierValueAsString"=>RecordIdentifierValueAsString); aws_config=aws_config)
delete_record(EventTime, FeatureGroupName, RecordIdentifierValueAsString, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = sagemaker_featurestore_runtime("DELETE", "/FeatureGroup/$(FeatureGroupName)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("EventTime"=>EventTime, "RecordIdentifierValueAsString"=>RecordIdentifierValueAsString), args)); aws_config=aws_config)

"""
    GetRecord()

Use for OnlineStore serving from a FeatureStore. Only the latest records stored in the OnlineStore can be retrieved. If no Record with RecordIdentifierValue is found, then an empty result is returned.

# Required Parameters
- `FeatureGroupName`: The name of the feature group in which you want to put the records.
- `RecordIdentifierValueAsString`: The value that corresponds to RecordIdentifier type and uniquely identifies the record in the FeatureGroup.

# Optional Parameters
- `FeatureName`: List of names of Features to be retrieved. If not specified, the latest value for all the Features are returned.
"""
get_record(FeatureGroupName, RecordIdentifierValueAsString; aws_config::AWSConfig=global_aws_config()) = sagemaker_featurestore_runtime("GET", "/FeatureGroup/$(FeatureGroupName)", Dict{String, Any}("RecordIdentifierValueAsString"=>RecordIdentifierValueAsString); aws_config=aws_config)
get_record(FeatureGroupName, RecordIdentifierValueAsString, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = sagemaker_featurestore_runtime("GET", "/FeatureGroup/$(FeatureGroupName)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("RecordIdentifierValueAsString"=>RecordIdentifierValueAsString), args)); aws_config=aws_config)

"""
    PutRecord()

Used for data ingestion into the FeatureStore. The PutRecord API writes to both the OnlineStore and OfflineStore. If the record is the latest record for the recordIdentifier, the record is written to both the OnlineStore and OfflineStore. If the record is a historic record, it is written only to the OfflineStore.

# Required Parameters
- `FeatureGroupName`: The name of the feature group that you want to insert the record into.
- `Record`: List of FeatureValues to be inserted. This will be a full over-write. If you only want to update few of the feature values, do the following:   Use GetRecord to retrieve the latest record.   Update the record returned from GetRecord.    Use PutRecord to update feature values.

"""
put_record(FeatureGroupName, Record; aws_config::AWSConfig=global_aws_config()) = sagemaker_featurestore_runtime("PUT", "/FeatureGroup/$(FeatureGroupName)", Dict{String, Any}("Record"=>Record); aws_config=aws_config)
put_record(FeatureGroupName, Record, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = sagemaker_featurestore_runtime("PUT", "/FeatureGroup/$(FeatureGroupName)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Record"=>Record), args)); aws_config=aws_config)


text:
function filter_coeffs(a::Vector{Float64}, b::Vector{Float64}, orders::Vector{Int64}, data::Array{Float64})

    dinv = (b - a)./(orders.-1)
    d = length(a)
    if d == 1
        coeffs = filter_coeffs_1d(dinv,data)
    elseif d == 2
        coeffs = filter_coeffs_2d(dinv,data)
    elseif d == 3
        coeffs = filter_coeffs_3d(dinv,data)
    elseif d == 4
        coeffs = filter_coeffs_4d(dinv,data)
    else
        error()
    end
    return coeffs
end

# function find_coeffs_1d_old(delta_inv, M, data)
#
#     # args: double, int, double[:], double[:
#
#     N = M + 2
#
#     data = data[:]
#
#     if length(data) == N
#         rhs = data
#     else
#         rhs = [0; data; 0]
#     end
#
#     basis = [1.0/6.0, 2.0/3.0, 1.0/6.0]
#
#     vals = repeat( basis, outer=[M])
#     co_x = repeat( collect(2:M+1), inner=[3])
#     co_y = repeat( collect(1:3) , outer=[M]) + co_x .- 2
#
#     db = 4
#     initial = [1.0,-2.0,1.0,0.0]*delta_inv*delta_inv
#     final = [0.0,1.0,-2.0,1.0]*delta_inv*delta_inv
#
#     vals = [initial; vals; final]
#     co_x = [ones(Int,db); co_x; ones(Int,db)*(M+2)]
#     co_y = [1:db; co_y; M+3-db:M+2]
#
#     spmat = sparse(co_x, co_y, vals)
#
#     sol = spmat \ rhs
#
#     return sol
#
# end

function solve_deriv_interp_1d!(bands::Matrix{Float64}, coefs::Vector{Float64})

    M = size(coefs,1) - 2

    # Solve interpolating equations
    # First and last rows are different

    bands[1, 2] /= bands[1, 1]
    bands[1, 3] /= bands[1, 1]
    bands[1, 4] /= bands[1, 1]
    bands[1, 1] = 1.0
    bands[2, 2] -= bands[2, 1] * bands[1, 2]
    bands[2, 3] -= bands[2, 1] * bands[1, 3]
    bands[2, 4] -= bands[2, 1] * bands[1, 4]
    bands[1, 1] = 0.0
    bands[2, 3] /= bands[2, 2]
    bands[2, 4] /= bands[2, 2]
    bands[2, 2] = 1.0

    # Now do rows 2 through M+1
    for row = 3:(M+1)
        bands[row, 2] -= bands[row, 1] * bands[row - 1, 3]
        bands[row, 4] -= bands[row, 1] * bands[row - 1, 4]
        bands[row, 3] /= bands[row, 2]
        bands[row, 4] /= bands[row, 2]
        bands[row, 1] = 0.0
        bands[row, 2] = 1.0
    end
    # Do last row
    bands[M + 2, 2] -= bands[M + 2, 1] * bands[M, 3]
    bands[M + 2, 4] -= bands[M + 2, 1] * bands[M, 4]
    bands[M + 2, 3] -= bands[M + 2, 2] * bands[M + 1, 3]
    bands[M + 2, 4] -= bands[M + 2, 2] * bands[M + 1, 4]
    bands[M + 2, 4] /= bands[M + 2, 3]
    bands[M + 2, 3] = 1.0

    coefs[M + 2] = bands[M + 2, 4]
    # Now back substitute up
    # for row in (M+1):2
    for ii in 1:M
        row = M+2 - ii
        coefs[row] = bands[row, 4] - bands[row, 3] * coefs[row + 1]
    end
    # Finish with first row
    coefs[1] = bands[1, 4] - bands[1, 2] * coefs[2] - bands[1, 3] * coefs[3]

end

function find_coeffs_1d!(delta_inv, M, data, coefs)

    # data is 1d array of length M
    # data is 1d array of length M+2

    basis = [1.0/6.0, 2.0/3.0, 1.0/6.0]
    bands = zeros(M+2, 4)

    # Setup boundary conditions
    abcd_left = zeros(4)
    abcd_right = zeros(4)

    # Left boundary
    abcd_left[1] = 1.0 * delta_inv * delta_inv
    abcd_left[2] = -2.0 * delta_inv * delta_inv
    abcd_left[3] = 1.0 * delta_inv * delta_inv
    abcd_left[4] = 0

    # Right boundary
    abcd_right[1] = 1.0 * delta_inv * delta_inv
    abcd_right[2] = -2.0 * delta_inv * delta_inv
    abcd_right[3] = 1.0 * delta_inv * delta_inv
    abcd_right[4] = 0

    for i in 1:4
        bands[1, i] = abcd_left[i]
        bands[M + 2, i] = abcd_right[i]
    end
    for i in 2:(M+1)
        for j in 1:3
            bands[i, j] = basis[j]
            bands[i, 4] = data[i-1]
        end
    end

    solve_deriv_interp_1d!(bands, coefs)

    return coefs
end

function find_coeffs_1d(delta_inv, M, data)
    coefs = zeros(M+2)
    if length(data) == M
        find_coeffs_1d!(delta_inv, M, data, coefs)
    else
        find_coeffs_1d!(delta_inv, M, data[2:end-1], coefs)
    end
    return coefs
end

function filter_coeffs_1d(dinv, data)

  M = size(data,1)
  N = M+2

  coefs = find_coeffs_1d(dinv[1], M, data)

  return coefs

end

function filter_coeffs_2d(dinv, data)

    Mx = size(data,1)
    My = size(data,2)

    Nx = Mx+2
    Ny = My+2

    coefs = zeros(Nx,Ny)

    # First, solve in the X-direction
    for iy in 1:My
        coefs[:,iy+1] = find_coeffs_1d(dinv[1], Mx, data[:,iy])
    end

    # Now, solve in the Y-direction
    for ix in 1:Nx
        coefs[ix,:] = find_coeffs_1d(dinv[2], My, coefs[ix,:])
    end

    return coefs

end

function filter_coeffs_3d(dinv, data)

    Mx = size(data,1)
    My = size(data,2)
    Mz = size(data,3)

    Nx = Mx+2
    Ny = My+2
    Nz = Mz+2

    coefs = zeros(Nx,Ny,Nz)

    # First, solve in the X-direction
    for iy in 1:My
        for iz in 1:Mz
            coefs[:,iy+1,iz+1] = find_coeffs_1d(dinv[1], Mx, data[:,iy,iz])
        end
    end

    # Now, solve in the Y-direction
    for ix in 1:Nx
        for iz in 1:Mz
            coefs[ix,:,iz+1] = find_coeffs_1d(dinv[2], My, coefs[ix,:,iz+1])
        end
    end

    # Now, solve in the Z-direction
    for ix in 1:Nx
        for iy in 1:Ny
            coefs[ix,iy,:] = find_coeffs_1d(dinv[3], Mz, coefs[ix,iy,:])
        end
    end

    return coefs

end

function filter_coeffs_4d(dinv, data)

    Mx = size(data,1)
    My = size(data,2)
    Mz = size(data,3)
    Mz4 = size(data,4)

    Nx = Mx+2
    Ny = My+2
    Nz = Mz+2
    Nz4 = Mz4+2

    coefs = zeros(Nx,Ny,Nz,Nz4)

    # First, solve in the X-direction
    for iy in 1:My
        for iz in 1:Mz
            for iz4 in 1:Mz4
                coefs[:,iy+1,iz+1,iz4+1] = find_coeffs_1d(dinv[1], Mx, data[:,iy,iz,iz4])
            end
        end
    end

    # Now, solve in the Y-direction
    for ix in 1:Nx
        for iz in 1:Mz
            for iz4 in 1:Mz4
                coefs[ix,:,iz+1,iz4+1] = find_coeffs_1d(dinv[2], My, coefs[ix,:,iz+1,iz4+1])
            end
        end
    end

    # Now, solve in the Z-direction
    for ix in 1:Nx
        for iy in 1:Ny
            for iz4 in 1:Mz4
                coefs[ix,iy,:,iz4+1] = find_coeffs_1d(dinv[3], Mz, coefs[ix,iy,:,iz4+1])
            end
        end
    end
    # Now, solve in the Z4-direction
    for ix in 1:Nx
        for iy in 1:Ny
            for iz in 1:Mz
                coefs[ix,iy,iz,:] = find_coeffs_1d(dinv[4], Mz, coefs[ix,iy,iz,:])
            end
        end
    end

    return coefs

end