text:
struct DataFrame <: AbstractFuture
    data::Future
    nrows::Future
    # TODO: Add offset for indexing
    # offset::Future
end

Banyan.convert(::Type{Future}, df::DataFrame) = df.data

# DataFrame creation

# function read_csv(pathname)
#     location = CSVPath(pathname)

#     data = Future()
#     len = Future(location.nrows)

#     sourced(data, location)
#     val(len)

#     pt(data, Block())
#     pt(len, Replicated())

#     @partitioned data begin end

#     FutureDataFrame(data, len)
# end

# TODO: Implement reading/writing for S3FS and HTTP and for CSV, Parquet, Arrow

# compute samples as soon as possible
# delay the computation of sample properties
# transform and select produce new columns though transform keeps previous ones unless they are overwritten
# produced columns either use element-wise functions or they use functions that must be reduced
#
# - columns
#   - which columns are potentially used for grouping
#   - when columns are kept but renamed without changing their divisions
#   - when columns are removed entirely
# - views

# function write_csv(df::FutureDataFrame, pathname)
#     location = CSVPath(pathname)
#     destined(df, location)
#     pt(df, Block())
#     mutated(df)
#     partition_mutated(df, :mut, Block(), CSVPath(pathname))
#     partition(df, Block(), )
#     @partitioned df begin end
# end

function read_csv(path::String; kwargs...)
    df_loc = RemoteSource(path; kwargs...)
    df_loc.src_name == "Remote" || error("$path does not exist")
    df_nrows = Future(df_loc.nrows)
    DataFrame(Future(source=df_loc), df_nrows)
end

read_parquet(p; kwargs...) = read_csv(p; kwargs...)
read_arrow(p; kwargs...) = read_csv(p; kwargs...)

# TODO: For writing functions, if a file is specified, enforce Replicated

function write_csv(df, path; invalidate_source=true, invalidate_sample=true, kwargs...)
    # destined(df, Remote(path, delete_from_cache=true))
    # mutated(df)
    # partitioned_with() do
    #     pt(df, Partitioned(df))
    # end
    # @partitioned df begin end
    # compute(df)
    # sourced(df, Remote(path)) # Allow data to be read from this path if needed in the future
    # destined(df, None())
    partitioned_with() do
        pt(df, Partitioned(df))
    end
    partitioned_computation(
        df,
        destination=RemoteDestination(path; invalidate_source=invalidate_source, invalidate_sample=invalidate_sample, kwargs...),
        new_source=_->RemoteSource(path)
    )
end

write_parquet(A, p; kwargs...) = write_csv(A, p; kwargs...)
write_arrow(A, p; kwargs...) = write_csv(A, p; kwargs...)

function Banyan.write_to_disk(df::DataFrame)
    partitioned_with() do
        pt(df, Partitioned(df))
    end
    partitioned_computation(df, destination=Disk())
end

# TODO: Duplicate above functions for Parquet, Arrow

# DataFrame sample

Banyan.sample_axes(df::DataFrames.DataFrame) = [1]
Banyan.sample_keys(df::DataFrames.DataFrame) = propertynames(df)

# TODO: Make these sample_* functions handle empty data frames

function Banyan.sample_divisions(df::DataFrames.DataFrame, key)
    # There are no divisions for empty data
    if isempty(df)
        return []
    end

    max_ngroups = sample_max_ngroups(df, key)
    ngroups = min(max_ngroups, 512)
    data = sort(map(orderinghash, df[!, key]))
    datalength = length(data)
    grouplength = div(datalength, ngroups)
    # We use `unique` here because if the divisions have duplicates, this could
    # result in different partitions getting the same divisions. The usage of
    # `unique` here is more of a safety precaution. The # of divisions we use
    # is the maximum # of groups.
    # TODO: Ensure that `unique` doesn't change the order
    unique([
        # Each group has elements that are >= start and < end
        (
            data[(i-1)*grouplength + 1],
            data[i == ngroups ? datalength : i*grouplength + 1]
        )
        for i in 1:ngroups
    ])
end

function Banyan.sample_percentile(df::DataFrames.DataFrame, key, minvalue, maxvalue)
    # If the data frame is empty, nothing between `minvalue` and `maxvalue` can
    # exist in `df`. so the percentile is 0.
    if isempty(df) || isnothing(minvalue) || isnothing(maxvalue)
        return 0
    end

    # NOTE: This may cause some problems because the way that data is ultimately split may
    # not allow a really fine division of groups. So in the case of some filtering, the rate
    # of filtering may be 90% but if there are only like 3 groups then maybe it ends up being like
    # 50% and that doesn't get scheduled properly. We can try to catch stuff like maybe by using
    # only 70% of total memory in scheduling or more pointedly by changing this function to
    # call sample_divisions with a reasonable number of divisions and then counting how many
    # divisions the range actually belongs to.

    count(map(o -> begin oh = orderinghash(o); oh >= minvalue && oh <= maxvalue end, df[!, key])) / nrow(df)

    # # minvalue and maxvalue should already be order-preserved hashes
    # # minvalue, maxvalue = orderinghash(minvalue), orderinghash(maxvalue)
    # divisions = sample_divisions(A, key)
    # percentile = 0
    # divpercentile = 1/length(divisions)
    # inminmax = false

    # # Iterate through divisions to compute percentile
    # for (i, (divminvalue, divmaxvalue)) in enumerate(divisions)
    #     # Check if we are between the minvalue and maxvalue
    #     if (i == 1 || minvalue >= divminvalue) && (i == length(divisions) || minvalue < divmaxvalue)
    #         inminmax = true
    #     end

    #     # Add to percentile
    #     if inminmax
    #         percentile += divpercentile
    #     end

    #     # Check if we are no longer between the minvalue and maxvalue
    #     if (i == 1 || maxvalue >= divminvalue) && (i == length(divisions) || maxvalue < divmaxvalue)
    #         inminmax = false
    #     end
    # end

    # percentile
end

Banyan.sample_max_ngroups(df::DataFrames.DataFrame, key) = isempty(df) ? 0 : div(nrow(df), maximum(combine(groupby(df[!, [key,]], key), nrow).nrow))
Banyan.sample_min(df::DataFrames.DataFrame, key) = isempty(df) ? nothing : minimum(map(orderinghash, df[!, key]))
Banyan.sample_max(df::DataFrames.DataFrame, key) = isempty(df) ? nothing : maximum(map(orderinghash, df[!, key]))

# DataFrame properties

DataFrames.nrow(df::DataFrame) = collect(df.nrows)
DataFrames.ncol(df::DataFrame) = ncol(sample(df))
Base.size(df::DataFrame) = (nrow(df), ncol(df))
Base.ndims(df::DataFrame) = 2
Base.names(df::DataFrame, args...) = names(sample(df), args...)
Base.propertynames(df::DataFrame) = propertynames(sample(df))

# DataFrame filtering

# grouping
# - group input on any key
# - group output with same key, divisions, different id
# - group input on any of a set of key (used for column selection), returning key that was used
# - group output with same key, divisions, same id or different if row selection
# - group input on any not in a set of key (used for column selection)
# - group output as unknown, same id or different if row selection
# Grouped(future, job) returns key, gpt, max_ngroups
# GroupedBy(future, key|keys, job) returns key, gpt, max_ngroups
# NotGroupedBy(future, key|keys, job) returns key, gpt, max_ngroups
# GroupedWith(future, key|keys, get_job())
# Balanced(dim=1)
# Distributing() for other cases
# Shuffled() for ID change

# function partition_for_filtering(df, res, res_nrows, args, kwargs)
#     partition(df, Replicated())
#     partition(res, Partitioned(); match=df) # TODO: Make Partitioned the default

#     partition(df, Blocked(;dim=1, balanced=true))
#     partition(df, Blocked(;dim=1, balanced=false))
#     # balanced is required because it determines what constraints are assigned
#     # if any w.r.t. MaxNPartitions, MemoryUsage
#     partition(res, Blocked(;dim=1, id='*', balanced=false))
#     # TODO: Make * in PTs both on client side and in pt_lib_info.json
#     # correspond to randomly generated Int32 IDs

#     partition_later() do job
#         for (key, gpt) in Grouped(df, job)
#             # We can't just constrain with distribution=:grouped because if
#             # this code region is the first in its stage and `df` was already
#             # partitioned with specific divisions in a previous stage, the
#             # scheduler wouldn't know that it can't re-split `df` with a
#             # different max_npartitions. So we must ensure the max_npartitions
#             # constraint is carried over here.
#             # TODO: Maybe introduce the concept of fixed constraints that are
#             # just always applied to certain values. This could be used to
#             # ensure a max_npartitions constraint is always applied
#             partition(df, gpt)
#             partition(
#                 res,
#                 Grouped(;key=key, id='*', balanced=false),
#                 match=df,
#                 on=:divisions
#             )
#         end
#     end

#     partition(res_nrows, Reducing(reducer=+))
#     partition(args, Replicated())
#     partition(kwargs, Replicated())
# end

# function futures_for_filtering(df, args, kwargs)
#     res_nrows = Future(df_nrows)
#     # There may still be skew in the data after filtering so we can't reduce
#     # the known memory usage of the result. We should mark the max_ngroups
#     # properties that were computed for keys as stale but this would mean that
#     # any subsequent GroupWith would result in new quantiles and hence, a
#     # shuffle. If the filtering was actually quite even, we wouldn't want to
#     # shuffle so we need a way to have GroupWith know what the previous
#     # quantiles are.
#     # TODO: Find a way so that re-splitting filtered data results in lower
#     # memory usage; maybe use same memory usage but make max_ngroups stale
#     res = DataFrame(Future(), res_nrows)
#     args = Future(args)
#     kwargs = Future(kwargs)

#     df, res, res_nrows, args, kwargs
# end

function pts_for_filtering(init::AbstractFuture, final::AbstractFuture; with, kwargs...)
    # for (initpt, finalpt) in zip(
    #     with(init; balanced=false, filtered_to=final, kwargs...),
    #     with(final; balanced=false, filtered_from=init, kwargs...),
    # )
    #     # unbalanced -> balanced
    #     pt(init, initpt, match=final, on=["distribution", "key", "divisions", "rev"])
    #     pt(final, Balanced() & Drifted())

    #     # unbalanced -> unbalanced
    #     pt(init, initpt, match=final, on=["distribution", "key", "divisions", "rev"])
    #     pt(final, finalpt & Drifted())

    #     # balanced -> unbalanced
    #     pt(init, Balanced(), match=final, on=["distribution", "key", "divisions", "rev"])
    #     pt(final, finalpt & Drifted())
    # end
    # Initially, we thought that the above was okay (not computing the
    # divisions by calling `with` with `balanced=true`). But then we realized
    # that you might have some blocked data that you then need to call `unique`
    # on and so you need to group it and then you're going to `comput` it right
    # afterwards. In this scenario, you need to have a PT where you compute
    # divisions with a call to `with` where `balanced=true`.
    for (initpt_unbalanced, finalpt_unbalanced, initpt_balanced, finalpt_balanced) in zip(
        # There should be a balanced and unbalanced PT for each possible key
        # that the initial/final data can be grouped on
        # unbalanced
        with(init; balanced=false, filtered_to=final, kwargs...),
        with(final; balanced=false, filtered_from=init, kwargs...),
        # balanced
        with(init; balanced=true, filtered_to=final, kwargs...),
        with(final; balanced=true, filtered_from=init, kwargs...),
    )
        # unbalanced -> balanced
        pt(init, initpt_unbalanced, match=final, on="divisions")
        pt(final, finalpt_balanced & Drifted())

        # unbalanced -> unbalanced
        pt(init, initpt_unbalanced, match=final, on=["distribution", "key", "divisions", "rev"])
        pt(final, finalpt_unbalanced & Drifted())

        # balanced -> unbalanced
        pt(init, initpt_balanced, match=final, on="divisions")
        pt(final, finalpt_unbalanced & Drifted())
    end
end

function DataFrames.dropmissing(df::DataFrame, args...; kwargs...)
    !get(kwargs, :view, false) || throw(ArgumentError("Cannot return view of filtered dataframe"))

    res_nrows = Future()
    res = DataFrame(Future(), res_nrows)
    args = Future(args)
    kwargs = Future(kwargs)

    # partition(df, Replicated())
    # partition(df, Blocked(dim=1, balanced=true))
    # partition(df, Blocked(dim=1, balanced=false))
    # union!(sample(res, :allowedgroupingkeys), sample(df, :allowedgroupingkeys))
    # partition_later() do
    #     union!(sample(df, :allowedgroupingkeys), sample(res, :allowedgroupingkeys))
    #     for key in sample(df, :allowedgroupingkeys)
    #         for balanced in [true, false]
    #             partition(df, Grouped(;key=key, balanced=balanced))
    #         end
    #     end
    # end

    partitioned_using_modules("DataFrames")
    partitioned_using() do
        # We need to maintain these sample properties and hold constraints on
        # memory usage so that we can properly handle data skew
        keep_all_sample_keys(res, df, drifted=true)
        keep_sample_rate(res, df)
    end

    partitioned_with() do
        pts_for_filtering(df, res, with=Distributed)
        pt(res_nrows, Reducing(quote (a, b) -> a .+ b end))
        pt(df, res, res_nrows, args, kwargs, Replicated())
    end

    # partition(res, Partitioned(balanced=false, id="*"), match=df, on=["distribution", "key", "divisions"])
    # partition(res_nrows, Reducing(reducer=+))

    @partitioned df res res_nrows args kwargs begin
        res = DataFrames.dropmissing(df, args...; kwargs...)
        res_nrows = DataFrames.nrow(res)
    end

    res
end

function Base.filter(f, df::DataFrame; kwargs...)
    !get(kwargs, :view, false) || throw(ArgumentError("Cannot return view of filtered dataframe"))

    f = Future(f)
    res_nrows = Future()
    res = DataFrame(Future(), res_nrows)
    kwargs = Future(kwargs)

    partitioned_using_modules("DataFrames")
    partitioned_using() do
        keep_all_sample_keys(res, df, drifted=true)
        keep_sample_rate(res, df)
    end

    partitioned_with() do
        pts_for_filtering(df, res, with=Distributed)
        pt(res_nrows, Reducing(quote (a, b) -> a .+ b end))
        pt(df, res, res_nrows, f, kwargs, Replicated())
    end

    @partitioned df res res_nrows f kwargs begin
        res = DataFrames.filter(f, df; kwargs...)
        res_nrows = DataFrames.nrow(res)
    end

    res
end

# TODO: Make a `used` field and ensure that splitting/merging functions don't get used if their used are not provided

# DataFrame element-wise

function Missings.allowmissing(df::DataFrame)::DataFrame
    res = Future()

    partitioned_using_modules("DataFrames")
    partitioned_using() do
        keep_all_sample_keys(res, df)
        keep_sample_rate(res, df)
    end

    partitioned_with() do
        pt(df, Distributed(df, scaled_by_same_as=res))
        pt(res, ScaledBySame(as=df), match=df)

        # pt(df, Distributed(df, balanced=true))
        # pt(res, Balanced(), match=df)

        # pt(df, Distributed(df, balanced=false, scaled_by_same_as=res))
        # pt(res, Unbalanced(scaled_by_same_as=df), match=df)

        pt(df, res, Replicated())
    end

    @partitioned df res begin res = DataFrames.allowmissing(df) end

    DataFrame(res, copy(df.nrows))
end

function Missings.disallowmissing(df::DataFrame)::DataFrame
    res = Future()

    partitioned_using_modules("DataFrames")
    partitioned_using() do
        keep_all_sample_keys(res, df)
        keep_sample_rate(res, df)
    end

    partitioned_with() do
        pt(df, Distributed(df, scaled_by_same_as=res))
        pt(res, ScaledBySame(as=df), match=df)

        # pt(df, Distributed(df, balanced=true))
        # pt(res, Balanced(), match=df)

        # pt(df, Distributed(df, balanced=false, scaled_by_same_as=res))
        # pt(res, Unbalanced(scaled_by_same_as=df), match=df)

        pt(df, res, Replicated())
    end

    @partitioned df res begin res = DataFrames.disallowmissing(df) end

    DataFrame(res, copy(df.nrows))
end

function Base.deepcopy(df::DataFrame)::DataFrame
    res = Future()

    partitioned_using_modules("DataFrames")
    partitioned_using() do
        keep_all_sample_keys(res, df)
        keep_sample_rate(res, df)
    end

    partitioned_with() do
        pt(df, Distributed(df, scaled_by_same_as=res))
        pt(res, ScaledBySame(as=df), match=df)

        # pt(df, Distributed(df, balanced=true))
        # pt(res, Balanced(), match=df)

        # pt(df, Distributed(df, balanced=false, scaled_by_same_as=res))
        # pt(res, Unbalanced(scaled_by_same_as=df), match=df)

        pt(df, res, Replicated())
    end

    @partitioned df res begin res = DataFrames.deepcopy(df) end

    DataFrame(res, copy(df.nrows))
end

function Base.copy(df::DataFrame)::DataFrame
    res = Future()

    partitioned_using_modules("DataFrames")
    partitioned_using() do
        keep_all_sample_keys(res, df)
        keep_sample_rate(res, df)
    end

    partitioned_with() do
        pt(df, Distributed(df, scaled_by_same_as=res))
        pt(res, ScaledBySame(as=df), match=df)

        # pt(df, Distributed(df, balanced=true))
        # pt(res, Balanced(), match=df)

        # pt(df, Distributed(df, balanced=false, scaled_by_same_as=res))
        # pt(res, Unbalanced(scaled_by_same_as=df), match=df)

        pt(df, res, Replicated())
    end

    @partitioned df res begin res = DataFrames.copy(df) end

    DataFrame(res, copy(df.nrows))

    # res = DataFrame(Future(), res_nrows)

    # partitioned_using() do
    #     keep_all_sample_keys(res, df)
    #     keep_sample_rate(res, df)
    # end

    # partitioned_with() do
    #     pt(df, Distributed(df))
    #     pt(res, ScaledBySame(as=df), match=df)
    # end

    # partition(df, Replicated())
    # partition(df, Blocked(dim=1, balanced=true))
    # partition(df, Blocked(dim=1, balanced=false))
    # union!(sample(res, :allowedgroupingkeys), sample(df, :allowedgroupingkeys))
    # setsample(res, :keystatistics, sample(df, :keystatistics))
    # partition_later() do
    #     union!(sample(df, :allowedgroupingkeys), sample(res, :allowedgroupingkeys))
    #     setsample(df, :keystatistics, sample(res, :keystatistics))
    #     for key in sample(df, :allowedgroupingkeys)
    #         for balanced in [true, false]
    #             partition(df, Grouped(;key=key, balanced=balanced))
    #         end
    #     end
    # end

    # partition(res, Partitioned(), match=res)
    # mutated(res)

    # @partitioned df res begin
    #     res = copy(df)
    # end

    # DataFrame(res, copy(df.nrows))

    # df_nrows = df.size
    # res_nrows = Future(df_nrows) # This marks df_nrows as mutated if the value can't be copied over.
    # res_data = Future(;same_as=df) # None supports same_as (also copies memory_usage), same_keys_as, same_keys.
    # res = DataFrame(res_data, res_nrows)

    # # If some PAs have assignments for values which aren't assinged in
    # # other PAs, copy the assignment over to those previous PAs

    # partition(df, Replicated())
    # partition(df, Blocked(dim=1))
    # partition_later() do job
    #     # Grouped produces PTs for all different keys that df can be grouped on
    #     # and for both balanced and unbalanced with MaxNPartitions and
    #     # AbsoluteMemoryUsage constraints as needed.
    #     # The by/not_by (for selecting keys by arguments or other information
    #     # from the location of the annotated code) or like/not_like (for
    #     # selecting keys based on other data) options allow for constraining the
    #     # keys that can be used for grouping. A some_key_required option
    #     # indicates that one of the keys must be used for grouping. If set to
    #     # false (the default), it indicates that only a key required by some
    #     # _other_ code region must be used.
    #     for _, gpt in Grouped(df, job)
    #         partition(df, gpt)
    #     end
    # end

    # # We assign either RelativeMemoryUsage or AbsoluteMemoryUsage constraints
    # # or neither if the data is already unbalanced
    # partition(res, Partitioned(); match=df, mutating=true)
    # partition(df_nrows, ReplicatedOrReducing())
    # partition(res_nrows, Partitioned(), match=df_nrows)

    # # We don't mutate res_nrows since whether that needs to be mutated is
    # # determined by the Future constructor that is invoked above

    # @partitioned df df_nrows res res_nrows begin
    #     res = copy(df)
    #     res_nrows = df_nrows
    # end

    # res
end

# DataFrame column manipulation (with filtering)

# function getindex(df::DataFrame, rows, cols)
#     res_data = Future()
#     res_nrows = Future()
#     res = DataFrame(res_data, res_nrows)

#     # TODO: Add in rows and cols

#     partition(df, Replicated())
#     partition(res, Replicated())

#     partition_later() do job
#         selectivity = sample(length(res)) / sample(length(df))
#         partition(df, Distributing(;distribution=:blocked))
#         partition(res, Distributing(;distribution=:blocked, id='*', balanced=false), memory_usage_factor=1/selectivity)
#         # TODO: Figure out exactly how memory usage propagates especially for filters followed by filters

#         for (key, gpt, like) in GroupedLike(df, res)
#             # `gpt` is a PT with distribution=:grouped and key set to something
#             # that df may be grouped by in the future and divisions either
#             # computed using compute_quantiles or left unspecified with
#             # balanced=false
#             #
#             # `GroupedLike` calls `compute_keys` on `res` and checks whether
#             # `df` is grouped by a key that `res` has or not
#             partition(df, gpt, max_npartitions=max_ngroups)
#             if like && rows isa Colon
#                 partition(res, PartitionType(), match=df)
#             elseif like
#                 partition(
#                     res,
#                     Distributing(distribution=:grouped, key=key, id='*'),
#                     match=df,
#                     on=:divisions
#                 )
#             else
#                 partition(res, Distributing(distribution=:unknown, id='*'))
#             end
#         end
#     end

#     # TODO: Use partition for res_nrows

#     @partitioned df res res_nrows begin
#         res = getindex(df)
#         res_nrows = nrow(res)
#     end

#     res
# end

function Base.getindex(df::DataFrame, rows=:, cols=:)
    # TODO: Accept either replicated, balanced, grouped by any of columns
    # in res (which is just cols of df if column selector is :), or unknown
    # and always with different ID unless row selector is :
    # and only allow : and copying columns without getting a view
    (rows isa Colon || rows isa BanyanArrays.Vector{Bool}) ||
        throw(ArgumentError("Expected selection of all rows with : or some rows with Vector{Bool}"))
    (cols != !) || throw(ArgumentError("! is not allowed for selecting all columns; use : instead"))

    # TODO: Remove this if not necessary
    if rows isa Colon && cols isa Colon
        return copy(df)
    end

    df_nrows = df.nrows
    return_vector = cols isa Symbol || cols isa String || cols isa Integer
    select_columns = !(cols isa Colon)
    filter_rows = !(rows isa Colon)
    columns = Symbol.(names(sample(df), cols))
    cols = Future(cols)
    rows = rows isa AbstractFuture ? rows : Future(rows)

    res_size =
        if filter_rows
            Future()
        elseif return_vector
            Future(df.nrows, mutation=tuple)
        else
            Future(df.nrows)
        end
    res =
        if return_vector
            BanyanArrays.Vector{eltype(sample(df)[!, collect(cols)])}(Future(), res_size)
        else
            DataFrame(Future(), res_size)
        end

    partitioned_using_modules("DataFrames")
    partitioned_using() do
        keep_all_sample_keys(res, df, drifted=filter_rows)
        keep_sample_rate(res, df)
    end

    partitioned_with() do
        if filter_rows
            for (dfpt_unbalanced, respt_unbalanced, dfpt_balanced, respt_balanced) in zip(
                # unbalanced
                Distributed(df; balanced=false, filtered_to=res),
                Distributed(res; balanced=false, filtered_from=df),
                # balanced
                Distributed(df; balanced=true, filtered_to=res),
                Distributed(res; balanced=true, filtered_from=df),
            )
                # Return Blocked if return_vector or select_columns and grouping by non-selected
                return_blocked = return_vector || (dfpt_balanced.distribution == "grouped" && !(dfpt_balanced.key in columns))

                # unbalanced -> balanced
                pt(df, dfpt_unbalanced, match=(return_blocked ? nothing : res), on=["distribution", "key", "divisions", "rev"])
                pt(res, (return_blocked ? Blocked(along=1) : respt_balanced) & Balanced() & Drifted())

                # unbalanced -> unbalanced
                pt(df, dfpt_unbalanced, match=(return_blocked ? nothing : res), on=["distribution", "key", "divisions", "rev"])
                pt(res, return_blocked ? Blocked(res, along=1, balanced=false, filtered_from=df) : respt_unbalanced & Drifted())
                pt(rows, Blocked(along=1) & ScaledBySame(as=df), match=df, on=["balanced", "id"])

                # balanced -> unbalanced
                pt(df, dfpt_balanced, match=(return_blocked ? nothing : res), on=["distribution", "key", "divisions", "rev"])
                pt(res, return_blocked ? Blocked(res, along=1, balanced=false, filtered_from=df) : respt_unbalanced & Drifted())
                pt(rows, Blocked(along=1) & ScaledBySame(as=df), match=df, on=(dfpt_balanced.distribution == "blocked" ? "balanced" : ["balanced", "id"]))
            end

            # pts_for_filtering(df, res, Blocked)
            # pt(rows, Block(along=1), match=df)

            # # blocked and balanced
            # pt(df, Blocked(along=1) & Balanced())
            # pt(rows, Blocked(along=1) & Balanced())
            # pt(res, Blocked(along=1) & Unbalanced() & Drifted())

            # pt(df, Blocked(df, balanced=false, filtered_to=res))
            # pt(rows, filter_rows ? Blocked() & Unbalanced(scaled_by_same_as=df) : Replicated())

            # for gpt in Grouped(df, filtered_to = filter_rows ? res : nothing)
            #     pt(df, gpt)
            #     # TODO: Handle select_columns
            # end

            pt(res_size, Reducing(return_vector ? quote (a, b) -> Tuple([a[1] + b[1], a[2:end]...]) end : quote (a, b) -> a .+ b end))
        else
            for dpt in Distributed(df, scaled_by_same_as=res)
                pt(df, dpt)
                if return_vector || (dpt.distribution == "grouped" && !(dpt.key in columns))
                    pt(res, Blocked(along=1) & ScaledBySame(as=df), match=df, on=["balanced", "id"])
                else
                    pt(res, ScaledBySame(as=df), match=df)
                end
            end
            pt(res_size, PartitionType(), match=df_nrows)
        end

        # if filter_rows
        #     pt(df, Blocked(df, balanced=true))
        #     pt(rows, Blocked(rows, balanced=true))

        #     pt(df, Blocked(df, balanced=false) | Grouped(df), match=rows, on=["balanced", "id"])
        #     pt(rows, Blocked(along=1))
        # else
        # end
        # pt(rows,  ? Blocked(along=1): Replicated())
        # pt(res)
        # pt(res_size)
        pt(df, res, res_size, rows, cols, Replicated())
        pt(df_nrows, Replicating())
    end

    @partitioned df df_nrows res res_size rows cols begin
        res = df[rows, cols]
        res_size = rows isa Colon ? df_nrows : size(res)
        res_size = res isa Base.Vector ? res_size : first(res_size)
    end

    res

    # by = names(sample(df), cols)
    # onecol = cols isa Symbol || cols isa String || cols isa Integer

    # # TODO: Maybe make @partitioned be called first so that we can have error handling
    # # get triggered first. But make sure that some stuff that needs to happen before
    # # @partitioned happens like creating futures and using mutated(future, new_future)

    # # TODO: Handle case where cols is a single column by returning a Vector

    # # TODO: Compute estimate of memory usage in code regions unless memory is
    # # already specified to be something non-zero

    # df_nrows = df.nrows
    # res_nrows = rows isa Colon ? Future(df_nrows, onecol ? tuple : identity) : Future()
    # res = onceol ? Vector(Future(), res_nrows) : DataFrame(Future(), res_nrows)
    # rows = rows isa Vector ? rows : Future(rows)
    # cols = Future(cols)

    # partition(df, Replicated())
    # partition(rows, Replicated())
    # partition(res, Replicated())

    # for balanced in [true, false]
    #     partition(df, Blocked(;dim=1, balanced=balanced))

    #     if rows isa Colon
    #         partition(rows, Replicated())
    #     elseif balanced
    #         partition(rows, Blocked(;dim=1, balanced=true))
    #     else
    #         partition(rows, Blocked(;dim=1); match=df, on="id")
    #     end

    #     if rows isa Colon
    #         partition(res, Blocked(;dim=1, balanced=balanced), match=df, on=["id"])
    #     else
    #         partition(res, Blocked(;dim=1, balanced=false, id="*"))
    #     end
    # end

    # # Merge sample properties of df to res
    # # TODO: Make sample properies here work with onecol to produce array
    # # TODO: Implement saving/loading functions for dataframe
    # # TODO: Implement array operations, saving/loading for them
    # if !onecol
    #     union!(sample(res, :allowedgroupingkeys), sample(df, :allowedgroupingkeys))
    #     intersect!(sample(res, :allowedgroupingkeys), by)
    #     if rows isa Colon
    #         # Only merge column statistics if the column exists after filtering
    #         for key in by
    #             if key in sample(res, :names)
    #                 setsample(res, :keystatistics, key, sample(df, :keystatistics, key))
    #             end
    #         end
    #     end
    # end

    # partition_later() do
    #     if !onecol
    #         # Merge sample properties of res to df
    #         union!(sample(df, :allowedgroupingkeys), sample(res, :allowedgroupingkeys))
    #     end

    #     # Distributed includes all Blocked and all Grouped on applicable keys
    #     # and both balanced and unbalanced
    #     # TODO: Have annotation code access and set job and sample properties
    #     # while PT constructor code is responsible for automatically
    #     # constructing constraints
    #     # TODO: Implement sample properties
    #     # - copying over from one future to another
    #     # - computing them
    #     # - renaming
    #     # - mutating
    #     for key in sample(df, :allowedgroupingkeys)
    #         for balanced in [true, false]
    #             # PT constructors serve to add in constraints needed to ensure
    #             # the parameters hold true. In other words, they make PTs
    #             # correct by construction
    #             partition(df, Grouped(;key=key, balanced=balanced))

    #             if rows isa Colon
    #                 partition(rows, Replicated())
    #             else
    #                 partition(rows, Blocked(;dim=1); match=df, on="id")
    #             end

    #             if !onecol && key in sample(res, :names)
    #                 if rows isa Colon
    #                     partition(res, Distributed(), match=df, on=["key", "divisions", "balanced", "id"])
    #                 else
    #                     partition(res, Distributed(balanced=false, id="*"), match=df, on=["key", "divisions"])
    #                 end
    #             else
    #                 if rows isa Colon
    #                     partition(res, Blocked(;dim=1), match=df, on=["balanced", "id"])
    #                 else
    #                     # PT constructors should:
    #                     # - Create MaxNPartitions and MemoryUsage constraints if balanced
    #                     # - Create random ID if *
    #                     partition(res, Blocked(;dim=1, balanced=false, id="*"))
    #                 end
    #             end
    #         end
    #     end
    # end

    # partition(cols, Replicated())
    # partition(df_nrows, ReplicatedOrReducing())
    # if rows isa Colon
    #     partition(res_nrows, Reducing(;reducer=(onceol ? .+ : +)))
    # else
    #     partition(res_nrows, ReplicatedOrReducing(match=df))
    # end

    # mutated(res)

    # @partitioned df df_nrows res res_nrows rows cols begin
    #     res = df[rows, cols]
    #     res_nrows = rows isa Colon ? df_nrows : length(res)
    #     res_nrows = res isa Vector ? (rows) : rows
    # end

    # res

    # TODO: Solve same issue with skew for joins and filters
    # TODO: Make the total memory usage of a variable a constraint
    # - grouped -> account for skew
    # If no memory usage is specified, we use whatever the value had before
    # each task has default memory usage which is used if this is the first time
    # the value is being used; otherwise we switch only if the stage has a PA that specifies something
    # TODO: Require a full re-split to change total memory usage of a value
    # don't copy over quantiles or max_npartitions for filtering operations
    # - larger memory usage and same number of partitions
    # - smaller memory usage but smaller number of partitions
    # potential solution:
    # - reuse sample statistics for quantiles and max_ngroups as much as possible
    # - in all operations, handle each possible grouping key by computing quantiles and using appropriate max_ngroups all where skewed=false
    # - in all operations, automatically add in memoryusage constraints using sample statistics if o constraint is provided
    # - in the scheduler, save the memory usage of each value and only update if re-splitting EDIT: actually this isn't needed but instead we
    # should have the option to not default to memory usage from sample so that some PAs that, for example, require skewed=true won't rely on
    # memoryusage constraint or max_npartitions constraints being added in and instead just reuse from the previous partitioning
    # - in PT library, don't split or cast to skewed=true
    # - in all operations that filter or change distribution, return value with skewed=true and a memory usage that is relative to size of sample TODO
    # - in all operations, handle grouping where skewd=true by not setting a new max_ngroups or memoryusage
    # TODO join by using quantiles of either side of join (potentially resulting in skewed result) or the quantiles of the result but with
    # memory usage of inputs adjusted to account for skew in how they are expected to be partitioned to produce a balanced join result
    # TODO: Compute percentile spanned by quantiles associated with min max of join inputs in distribution of join output
    # TODO: Estimate memory usage of skewed output by multiplying outptut memory usage by selectivity for balanced and unknown and for grouped,
    # using percentiles of min amx of keys
    # TODO: Make partition accept parameters for adding in memoryusage constraints if they are needed (no_memory_usage, memory_usage_factor)
    # TODO: Make `evaluate` fill in uninitialized memoryusage constraints, so `partition` will either create MemoryUsage(some initial factor
    # to be multiplied with such as 1) or no constraint at all

    # TODO: MemoryUsage
    # - set memory usage of a variable to what is estimated by the sample if it is balanced
    # - otherwise set it to be relative to some other variable
    # - otherwise abstain and use what it already had in previous stage
    # options
    # - relative to optionally with skew
    # - equal to
    # - sampled
    # - no constraint
    # memory usage of each value has the following
    # - a base value that is either some absolute size or relative to one ore more values
    # - a scaling factor based on skew
    # We can achieve this by the following:
    # - keep track of memory usage of each variable in the backend
    # - for each potential fused PA for a stage, determine memory usage by looking at constraints

    # res_data = cols isa Colon ? Future(;same_keys_as=df, scaled_by_same_as=df) : Future()
    # res_data = Future()

    # TODO: Maybe make some way to mutate the value according to the later
    # approximation of cols
    # TODO: Store length instead of size for DataFrame
    # res_nrows = rows isa Colon ? Future(df.size, ???) : Future()

    # res_nrows = Future()
    # res = DataFrame(res_data, res_nrows)

    # partition(df, Replicating())
    # partition(res, Replicating(), match=df)

    # partition(df, Balanced(dim=1))
    # partition(df, Distributing(distribution=:unknown))
    # if rows isa Colon
    #     partition(res, PartitionType(); match=df, on="id")
    # else:
    #     partition(res, Distributing(distribution=:unknown, id='*'))
    # end

    # None location constructor - determines sample-related stuff

    # find min-max of filtered
    # find quantile percent spanned by min max range
    # take the inverse of the percentage and multiply by mem. usage
    # TODO: On filtering operations, group by quantiles to get
    # - have constraint that specifies adjustment to memory usage
    # - on filtering, look at quantiles to determine adjustment

    # Make compute_divisions have option to check for divisons and only update if
    # the difference is too much
    # TODO: Figure out relationship between filtering and divisions and memory usage

    # partition_later() do job
    #     # TODO: Maybe seperate into PT constructor and constraint constructor for AtMost
    #     for (key, like, gpt, max_ngroups) in GroupedLike(df, res, job)
    #         partition(df, gpt, max_npartitions=max_ngroups)
    #         if like && rows isa Colon
    #             partition(res, PartitionType(), match=df)
    #         elseif like
    #             partition(
    #                 res,
    #                 Distributing(distribution=:grouped, key=key, id='*'),
    #                 match=df,
    #                 on=:divisions
    #             )
    #         else
    #             partition(res, Distributing(distribution=:unknown, id='*'))
    #         end
    #     end
    # end

    # partition(res_nrows, rows isa Colon ? Replicated() : Reduce(dim=1))

    # TODO: Make mut accept multiple arguments
    # TODO: Implement partition

    # mutated(res, res_nrows)
    # @partitioned df res res_nrows rows cols begin
    #     res = getindex(df, rows)
    #     res_nrows = size(res)
    # end
end

function Base.setindex!(df::DataFrame, v::Union{BanyanArrays.Vector, BanyanArrays.Matrix, DataFrame}, rows, cols)
    rows isa Colon || throw(ArgumentError("Cannot mutate a subset of rows in place"))

    # selection = names(sample(df), cols)

    res = Future()
    # cols = Future(Symbol.(names(sample(df), cols)))
    cols = Future(cols)

    partitioned_using_modules("DataFrames")
    partitioned_using() do
        keep_all_sample_keys(res, df)
        keep_sample_rate(res, df)
    end

    partitioned_with() do
        for dpt in Distributed(df, scaled_by_same_as=res)
            pt(df, dpt)
            pt(res, ScaledBySame(as=df), match=df)

            # The array that we are inserting into this dataframe must be
            # partitioned with the same ID or it must be perfectly balanced
            # if the original dataframe is also balanced.
            if dpt.distribution == "blocked" && dpt.balanced
                pt(v, Blocked(along=1) & Balanced())
            else
                pt(v, Blocked(along=1), match=df, on=["balanced", "id"])
            end
        end

        pt(df, res, v, cols, Replicated())
    end

    # # union!(sample(res, :allowedgroupingkeys), sample(df, :allowedgroupingkeys))
    # # # Only merge column statistics if the column is not mutated by setindex!
    # # for key in sample(df, :names)
    # #     if !(key in selection)
    # #         setsample(res, :keystatistics, key, sample(df, :keystatistics, key))
    # #     end
    # # end

    # partition(df, Replicated())
    # partition(col, Replicated())

    # partition(df, Blocked(dim=1, balanced=true))
    # partition(v, Blocked(dim=1, balanced=true))

    # partition(df, Blocked(dim=1, balanced=false))
    # partition(res, Partitioned(); match=df)

    # # TODO: Implement partition_later to end the current PA and store
    # # function that can potentially (but may also not) create a new PA

    # partition_later() do
    #     for key in sample(res, :allowedgroupingkeys)
    #         if key in sample(df, :names)
    #             push!(sample(res, :allowedgroupingkeys))
    #         end
    #     end
    #     for key in sample(res, :names)
    #         if !(key in selection)
    #             setsample(df, :keystatistics, key, sample(res, :keystatistics, key))
    #         end
    #     end

    #     for key in sample(df, :allowedgroupingkeys)
    #         for balanced in [true, false]
    #             partition(df, Grouped(;key=key, balanced=balanced))
    #             if key in selection
    #                 partition(res, Blocked(dim=1, balanced=balanced); match=df, on="id")
    #             else
    #                 partition(res, Partitioned(); match=df)
    #             end
    #         end
    #     end
    # end

    # partition(v, Blocked(dim=1, balanced=false), match=df, on="id")
    # partition(cols, Replicated())

    # # res_data = cols isa Colon ? Future()
    # # res_nrows = rows isa Colon ? Future(df.size) : Future()
    # # res = DataFrame(res_data, res_nrows)

    # # mut will replace df with the future and sample and location of res and
    # # record a DestroyRequest on the original future in df. The contents of
    # # res are modified to ensure that when its destructor is called a
    # # DestroyRequest is not called on the new future which is now stored in df.
    # # This could be done by setting the value_id to nothing and having the
    # # finalizer of Future only record a DestroyRequest if the value ID isn't
    # # set to nothing. This should also indicate that a sample should not be
    # # collected for the thing being mutated.
    # # TODO: Implement this version of mut
    # mutated(df, res)

    @partitioned df v cols res begin
        df[:, cols] = v
        res = df
    end

    # partition()  begin job
    #     for gpt, max_ngroups in GroupedBy(df, )
    #         partition()
    #     end
    # end
end

# function Base.getproperty(df::DataFrame, cols)
# end

# function Base.setproperty!(df::DataFrame, cols)
#     # Groups on any of columns of df but resulting df is not grouped on anything that is in the columns
#     # being set
#     # 2 cases
#     # - any of the columns in the input DF - GroupedBy (equivalents for Balanced for arrays)
#     # - any of the columns in the args - GroupedBy
#     # - any of the columns not in the args but in the input DF - NotGroupedBy
#     # - any of the columns in the args but actually use quantile - GroupedWith
#     # 3 cases
#     # - compute atmost and quantiles
#     # - compute atmost
#     # - don't compute anything
#     # define functions for each sample type
#     # - compute_keys
#     # - compute_axes
#     # - compute_max_ngroups
#     # - compute_quantiles
#     # define PT constructors that accept only a Future and columns to group on
#     # Distributed(distribution, key, divisions, id)
#     # - GroupedBy(compute_divisions=true)
#     # - GroupedBy(compute_divisions=true, compute_max_ngroups=true)
#     allcols = hash.(df |> approximate |> eachcol |> keys)
#     partition() begin job
#         for gpt, max_ngroups in GroupedBy(df, )
#             partition()
#         end
#     end
#     # TODO: Accept either replicated, balanced, or same ID unless rows are selected grouped by any of columns
#     # in res (which is just cols of df if column selector is :), or unknown
#     # and always with different ID unless row selector is :
#     # and only allow : and copying columns without getting a view
#     # TODO: Make the Location just produce a sample and then have this sample used
#     # to produce the total memory usage; effectively have the sample translated into
#     # total memory usage when compiled
#     # TODO: Make sure that in generated code on client side, mutated values
#     # have their samples copied before hand

#     # Several PT constructors that may be useful
#     # BlockAndBalanced(keys)
#     # BlockAndUnbalanced(keys)
#     # GroupedAndBalanced(keys) -> Distributing...
#     # GroupedAndUnbalanced(keys) -- key, gpt with constraints

#     # General:
#     # - Replicated
#     # - Blocked by some or all key(s)
#     # - Grouped by some or all key(s)
#     # Considerations:
#     # - copying over column statistics
#     # - partition both data and length
#     # - Balanced/Unbalanced
#     # - MemoryUsage, MaxNPartitions, Match, MatchOn
#     # - partition now or later
#     # - mutating

#     # The important constraints that are produced in PT constructors or in
#     # `partition` are MaxNPartitions and MemoryUsage

#     # MemoryUsage has three options
#     # 1. Specific memory usage that this value is guaranteed to have
#     # 2. Memory usage relative to some other value by some factor (typically based on selectivity, samples, min-max of keys and quantiles)
#     # 3. No constraint if it is required for something to be unbalanced (if no constraints are applied, the most recent size is used in scheduler)
# end

function DataFrames.rename(df::DataFrame, args...; kwargs...)
    res = Future()
    args = Future(args)
    kwargs = Future(kwargs)

    partitioned_using_modules("DataFrames")
    partitioned_using() do
        keep_all_sample_keys_renamed(res, df)
        keep_sample_rate(res, df)
    end

    partitioned_with() do
        # distributed
        for dfpt in Distributed(df, scaled_by_same_as=res)
            pt(dfpt)
            if dfpt.distribution == "grouped"
                groupingkeyindex = indexin(dfpt.key, sample(df, :keys))
                groupingkey = sample(res, :keys)[groupingkeyindex]
                pt(res, Grouped(by=groupingkey) & ScaledBySame(as=df), match=df, on=["balanced", "id", "divisions", "rev"])
            else
                pt(res, ScaledBySame(as=df), match=df)
            end
        end

        # replicated
        pt(df, res, args, kwargs, Replicated())
    end

    @partitioned df res args kwargs begin
        res = DataFrames.rename(df, args...; kwargs...)
    end

    DataFrame(res, copy(df.nrows))

    # TODO: Make partition_delayed calls be processed in reverse
    # TODO: Populate groupingkeys of samples forwards and backwards
    # TODO: Copy over sample properties forwards and backwards

    # res = Future()
    # args = Future(args)
    # kwargs = Future(kwargs)

    # partition(df, Replicated())
    # partition(res, Replicated())

    # for balanced in [true, false]
    #     partition(df, Blocked(;dim=1, balanced=balanced))
    #     partition(res, Blocked(;dim=1, balanced=balanced), match=df, on=["id"])
    # end

    # partition_later() do
    #     allowedgroupingkeys = sample(res, :allowedgroupingkeys)
    #     statistics = sample(res, :keystatistics, key)
    #     for (i, key) in enumerate(sample(res, :names))
    #         prevkey = sample(df, :names)[i]
    #         if key in keys(statistics)
    #             setsample(df, :keystatistics, prevkey, statistics)
    #         end
    #         if key in allowedgroupingkeys
    #             push!(sample(df, :allowedgroupingkeys), prevkey)
    #             for balanced in [true, false]
    #                 partition(df, Grouped(;key=prevkey, balanced=balanced))
    #                 partition(res, Grouped(;key=key), match=df, on=["divisions", "balanced", "id"])
    #             end
    #         end
    #     end
    # end

    # # for k in keys(d)
    # #     gpt, max_ngroups = GroupedBy(df; key=k)
    # #     partition(df, gpt, at_most=max_ngroups)
    # #     partition(df, Grouped(res; key=k), match=df, on=[:divisions, :id])
    # # end
    # # partition() do job
    # #     for (gpt, max_ngroups) in GroupedBy(df, stalekeys)
    # #         partition(df, gpt, at_most=max_ngroups)
    # #         partition(res)
    # #     end
    # # end

    # # res = DataFrame()
    # # partition(df, GroupedBy())
    # # partition(res, GroupedBy(), match=df, on=["distribution", "axis", "divisions"])
    # # res

    # mutated(res)

    # @partitioned df res args kwargs begin
    #     res = rename(df, args...; kwargs...)
    # end

    # allowedgroupingkeys = sample(df, :allowedgroupingkeys)
    # statistics = sample(df, :keystatistics, key)
    # for (i, key) in enumerate(sample(df, :names))
    #     newkey = sample(res, :names)[i]
    #     if key in allowedgroupingkeys
    #         push!(sample(res, :allowedgroupingkeys), newkey)
    #     end
    #     if key in keys(statistics)
    #         setsample(res, :keystatistics, newkey, statistics[key])
    #     end
    # end

    # DataFrame(res, copy(df.nrows))
end

# Make AtMost only accept a value (we can support PT references in the future if needed)
            # TODO: Make scheduler check that the values in AtMost or ScaledBy are actually present to ensure
            # that the constraint can be satisfied for this PT to be used# TODO: Implement select/transform/combine for DataFrame
# TODO: In rename/select/transform/combine figure out what columns to copy over statistics for
# - compute grouping column, owned columns, mutated columns and pass these to PT constructors so that
# they can filter to the relevant subset
# - ensure that we don't reuse statistics for keys that are mutated or assigned to, basically taking into
# account that the key sample properties might be there for keys that aren't actually present in the DF

# function select(df::DataFrame, args...; kwargs...)
#     !get(kwargs, :copycols, false) || throw(ArgumentError("Cannot return view of selected dataframe"))

#     for arg in args
#         source = arg isa Pair
#     end

#     res = Future()

#     # same statistics only if the column is not in the results of a transformation

#     @partitioned df args kwargs begin

#     end

#     DataFrame(res, copy(df.nrows))
# end

# function transform(df::DataFrame, args...; kwargs...)
#     !get(kwargs, :copycols, false) || throw(ArgumentError("Cannot return view of transformed dataframe"))
# end

# function combine(df::DataFrame, args...; kwargs...)
# end

# DataFrame shuffling

function Base.sort(df::DataFrame, cols=:; kwargs...)
    !get(kwargs, :view, false) || throw(ArgumentError("Cannot return view of sorted dataframe"))

    # TODO: Support a UserColOrdering passed into cols
    # # Determine what to sort by and whether to sort in reverse
    # firstcol = first(cols)
    # sortingkey, isreversed = if firstcol isa DataFrames.UserColOrdering
    #     if isempty(firstcols.kwargs)
    #         firstcols.col, get(kwargs, :rev, false)
    #     elseif length(firstcol.kwargs) == 1 && haskey(firstcol.kwargs, :rev)
    #         firstcols.col, get(firstcol.kwargs, :rev, false)
    #     else
    #         throw(ArgumentError("Only rev is supported for ordering"))
    #     end
    # else
    #     first(names(sample(df), firstcol)), get(kwargs, :rev, false)
    # end

    res = Future()
    columns = Symbol.(names(sample(df), cols))
    cols = Future(cols)
    kwargs = Future(kwargs)
    sortingkey = first(columns)
    isreversed = get(collect(kwargs), :rev, false)

    # TODO: Change to_vector(x) to [x;]

    partitioned_using_modules("DataFrames")
    partitioned_using() do
        keep_sample_keys(sortingkey, res, df)
        keep_sample_rate(res, df)
    end

    partitioned_with() do
        # We must construct seperate PTs for balanced=true and balanced=false
        # because these different PTs have different required constraints
        # TODO: Implement reversed in Grouped constructor
        pt(df, Grouped(df, by=sortingkey, rev=isreversed) | Replicated())
        # The Match constraint only applies to PT parameters. So we also need
        # a constraint to ensure that however `df` is grouped (balanced or
        # unbalanced), `res` is scaled by the same factor to account for
        # unbalanced distribution.
        pt(res, ScaledBySame(as=df), match=df)
        pt(df, res, cols, kwargs, Replicated())
    end

    # mutated(res)

    @partitioned df res cols kwargs begin
        res = DataFrames.sort(df, cols; kwargs...)
    end

    # statistics = sample(df, :keystatistics, by)
    # if by in keys(statistics)
    #     setsample(res, :keystatistics, by, statistics[by])
    # end

    DataFrame(res, copy(df.nrows))
end

function DataFrames.innerjoin(dfs::DataFrame...; on, kwargs...)
    length(dfs) >= 2 || throw(ArgumentError("Join requires at least 2 dataframes"))

    # TODO: Make it so that the code region's sampled computation is run first to allow for the function's
    # error handling to kick in first

    # TODO: Change this annotation to allow for grouping on any of the keys we
    # are joining on

    groupingkeys = first(on isa Base.Vector ? on : [on])
    groupingkeys = Symbol.(groupingkeys isa Union{Tuple,Pair} ? [groupingkeys...] : Base.fill(groupingkeys, length(dfs)))

    res_nrows = Future()
    res = DataFrame(Future(), res_nrows)
    on = Future(on)
    kwargs = Future(kwargs)

    partitioned_using_modules("DataFrames")
    partitioned_using() do
        # NOTE: `to_vector` is necessary here (`[on...]` won't cut it) because
        # we allow for a vector of pairs
        # TODO: Determine how to deal with the fact that some groupingkeys can
        # stick around in the case that we are doing a broadcast join
        # TODO: Make it so that the function used here and in groupby/sort
        # simply adds in the grouping key that was used
        keep_sample_keys_named(
            [df => groupingkey for (df, groupingkey) in zip(dfs, groupingkeys)]...,
            res => first(groupingkeys),
            drifted = true,
        )
        # NOTE: We are adjusting the sample rate accordingly, but we still need
        # to note that skew can occur in the selectivity of the join.
        # Therefore, we create ScaleBy constraints just for the
        # selectivity/skew issue - not for the sample rate.
        keep_sample_rate(res, dfs...)
    end

    # TODO: Use something like this for join
    partitioned_with() do
        # unbalanced, ...., unbalanced -> balanced - "partial sort-merge join"
        dfs_with_groupingkeys = [df => groupingkey for (df, groupingkey) in zip(dfs, groupingkeys)]
        for (df, groupingkey) in dfs_with_groupingkeys
            pt(df, Grouped(df, by=groupingkey, balanced=false, filtered_to=(res=>first(groupingkeys))), match=res, on=["divisions", "rev"])
        end
        pt(res, Grouped(res, by=first(groupingkeys), balanced=true, filtered_from=dfs_with_groupingkeys) & Drifted())

        # balanced, unbalanced, ..., unbalanced -> unbalanced
        for i in 1:length(dfs)
            # "partial sort-merge join"
            for (j, (df, groupingkey)) in enumerate(dfs_with_groupingkeys)
                pt(df, Grouped(df, by=groupingkey, balanced=(j==i), filtered_to=(res => first(groupingkeys))), match=dfs[i], on=["divisions", "rev"])
            end
            pt(res, Grouped(res, by=first(groupingkeys), balanced=false, filtered_from=dfs_with_groupingkeys[i]) & Drifted(), match=dfs[i], on=["divisions", "rev"])

            # broadcast join
            pt(dfs[i], Distributed(dfs[i]))
            for (j, df) in enumerate(dfs)
                if j != i
                    pt(df, Replicated())
                end
            end
            pt(res, ScaledBySame(as=dfs[i]), match=dfs[i])

            # TODO: Ensure that constraints are copied backwards properly everywhere
        end

        # TODO: Implement a nested loop join using Cross constraint. To
        # implement this, we may need a new PT constructor thaor some new
        # way of propagating ScaleBy constraints
        # for dpts in IterTools.product([Distributed(dpt, )]...)
        # pt(dfs..., Distributing(), cross=dfs)
        # pg(res, Blocked() & Unbalanced() & Drifted())

        # unbalanced, unbalanced, ... -> unbalanced - "partial sort-merge join"
        for (df, groupingkey) in dfs_with_groupingkeys
            pt(df, Grouped(df, by=groupingkey, balanced=false, filtered_to=(res => first(groupingkeys))), match=res, on=["divisions", "rev"])
        end
        pt(res, Grouped(res, by=first(groupingkeys), balanced=false, filtered_from=dfs_with_groupingkeys) & Drifted())

        # "replicated join"
        pt(res_nrows, Reducing(quote (a, b) -> a .+ b end))
        pt(dfs..., on, kwargs, res, res_nrows, Replicated())

        # TODO: Support nested loop join where multiple are Block and Cross-ed and others are all Replicate
    end

    @partitioned dfs on kwargs res res_nrows begin
        res = DataFrames.innerjoin(dfs...; on=on, kwargs...)
        res_nrows = nrow(res)
    end

    # partition(dfs..., Replicated())
    # partition(df2, Replicated())
    # partition(res, Replicated())
    # partition(res_nrows, Replicated())

    # for balanced in [true, false]
    #     partition(df1, Replicated())
    #     partition(df2, Replicated())
    #     partition(res, Blocked())
    # end

    # # TODO: Support broadcast join
    # # TODO: Implement proper join support where different parties are used for
    # # determining the distribution
    # # TODO: Implement distributed from balanced + unbalanced => unbalanced, unbalanced => unbalanced
    # # TODO: Maybe make MatchOn prefer one argument over the other so that we can take both cases of
    # # different sides of the join having their divisions used

    # on = [k isa Pair ? k : (k => k) for k in (on isa Vector ? on : [on])]
    # on_left = first.(on)
    # union!(sample(res, :allowedgroupingkeys), on_left)
    # union!(sample(df1, :allowedgroupingkeys), on_left)
    # union!(sample(df2, :allowedgroupingkeys), last.(on))
    # # No statistics are copied over because joins are selective and previous
    # # statistics would no longer apply
    # for key in on
    #     for balanced in [true, false]
    #         partition(df1, Grouped(;key=first(key), balanced=balanced))
    #         partition(df2, Grouped(;key=last(key), balanced=balanced))
    #         partition(res, Grouped(;balanced=false, id="*"), match=df1, on=["key", "divisions"])
    #     end
    # end

    # partition(res_nrows, Reducing(;reducer=+))
    # partition(kwargs, Replicated())

    # mutated(res)
    # mutated(res_nrows)

    # @partitioned df1 df2 res res_nrows kwargs begin
    #     res = innerjoin(df1, df2; kwargs...)
    #     res_nrows = nrow(res)
    # end

    res
end

function DataFrames.unique(df::DataFrame, cols=:; kwargs...)
    !get(kwargs, :view, false) || throw(ArgumentError("Returning a view of a Banyan data frame is not yet supported"))
    !(cols isa Pair || (cols isa Function && !(cols isa Colon))) || throw(ArgumentError("Unsupported specification of columns for which to get unique rows"))

    # TOOD: Just reuse select here

    # TODO: Check all usage of first
    res_nrows = Future()
    res = DataFrame(Future(), res_nrows)
    columns = Symbol.(names(sample(df), cols))
    cols = Future(cols)
    kwargs = Future(kwargs)

    partitioned_using_modules("DataFrames")
    partitioned_using() do
        keep_sample_keys(columns, res, df, drifted=true)
        keep_sample_rate(res, df)
    end

    partitioned_with() do
        pts_for_filtering(df, res, with=Grouped, by=columns)
        pt(res_nrows, Reducing(quote (a, b) -> a .+ b end))
        pt(df, res, res_nrows, cols, kwargs, Replicated())
    end

    @partitioned df res res_nrows cols kwargs begin
        res = DataFrames.unique(df, cols; kwargs...)
        res_nrows = DataFrames.nrow(res)
    end

    res

    # res_nrows = Future()
    # res = DataFrame(Future(), res_nrows)
    # cols = Future(cols)

    # partition(df, Replicated())
    # partition(res, Replicated())
    # partition(res_nrows, Replicated())

    # # TODO: Only use at most first 8 for allowedgroupingkeys
    # by = isnothing(cols) ? sample(df, :names) : names(sample(df), cols)
    # union!(sample(df, :allowedgroupingkeys), by)
    # union!(sample(res, :allowedgroupingkeys), by)
    # for key in by
    #     for balanced in [true, false]
    #         partition(df, Grouped(;key=key, balanced=balanced))
    #     end
    # end

    # partition(res, Grouped(;balanced=false, id="*"), match=df, on=["key", "divisions"])
    # partition(res_nrows, Reducing(;reducer=+))
    # partition(cols, Replicated())

    # mutated(res)
    # mutated(res_nrows)

    # @partitioned df res res_nrows cols begin
    #     res = isnothing(cols) ? unique(df) : unique(df, cols)
    #     res_nrows = nrow(res)
    # end

    # res
end

function DataFrames.nonunique(df::DataFrame, cols=:; kwargs...)
    !get(kwargs, :view, false) || throw(ArgumentError("Cannot return view of Banyan dataframe"))
    !(cols isa Pair || (cols isa Function && !(cols isa Colon))) || throw(ArgumentError("Unsupported specification of columns for which to get unique rows"))

    # TOOD: Just reuse select here

    df_nrows = df.nrows
    res_size = Future(df.nrows, mutation=tuple)
    res = BanyanArrays.Vector{Bool}(Future(), res_size)
    columns = Symbol.(names(sample(df), cols))
    cols = Future(cols)
    kwargs = Future(kwargs)

    partitioned_using_modules("DataFrames")
    partitioned_using() do
        keep_sample_rate(res, df)
    end

    partitioned_with() do
        pt(df, Grouped(df, by=columns))
        pt(res, Blocked(along=1) & ScaledBySame(as=df), match=df, on=["balanced", "id"])
        pt(df_nrows, Replicating())
        pt(res_size, PartitionType(), match=df_nrows)
        pt(df, res, df_nrows, res_size, cols, kwargs, Replicated())
    end

    @partitioned df df_nrows res res_size cols kwargs begin
        res = DataFrames.nonunique(df, cols; kwargs...)
        res_size = Tuple(df_nrows)
    end

    res

    # if get(kwargs, :view, false) throw(ArgumentError("Cannot return view of Banyan dataframe")) end

    # res_nrows = Future()
    # res = DataFrame(Future(), res_nrows)
    # cols = Future(cols)

    # partition(df, Replicated())
    # partition(res, Replicated())
    # partition(res_nrows, Replicated())

    # by = isnothing(cols) ? sample(df, :names) : names(sample(df), cols)
    # union!(sample(df, :allowedgroupingkeys), by)
    # union!(sample(res, :allowedgroupingkeys), by)
    # for key in by
    #     for balanced in [true, false]
    #         partition(df, Grouped(;key=key, balanced=balanced))
    #     end
    # end

    # partition(res, Grouped(;balanced=false, id="*"), match=df, on=["key", "divisions"])
    # partition(res_nrows, Reducing(;reducer=+))
    # partition(cols, Replicated())

    # # TODO: Maybe make mut automatic or re-evaluate exactly what it is used for
    # mutated(res)
    # mutated(res_nrows)

    # @partitioned df res res_nrows cols begin
    #     res = isnothing(cols) ? nonunique(df) : nonunique(df, cols)
    #     res_nrows = nrow(res)
    # end

    # res
end

# TODO: Implement SubDataFrame

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
export sync, sync!
function get_sync_t(s::Union{String,DateTime}, t::Array{Int64,1})
  isa(s, DateTime) && return floor(Int64, d2u(s)*sμ)
  if s == "first"
    return minimum(t)
  elseif s == "last"
    return maximum(t)
  else
    return floor(Int64, d2u(DateTime(s))*sμ)
  end
end

# # still used for irregular data
function get_sync_inds(t::AbstractArray{Int64,1}, Ω::Bool, t₀::Int64, t₁::Int64)
  if Ω == true
    return union(findall(t.<t₀), findall(t.>t₁))
  else
    return findall(t.<t₀)
  end
end

function get_del_ranges(xi::Array{Int64, 2}, nx::Int64)
  nw = size(xi, 1)
  x_del = Array{UnitRange, 1}(undef, 0)
  sizehint!(x_del, nw+1)

  # Does first row start at 1?
  if xi[1] != 1
    push!(x_del, UnitRange(1, xi[1]-1))
  end

  # Any gaps in rows 2:n-1
  for i in 2:nw
    if xi[i,1] - xi[i-1,2] > 1
      push!(x_del, UnitRange(xi[i-1,2]+1, xi[i,1]-1))
    end
  end

  # Does last row end at nx?
  if xi[nw,2] < nx
    push!(x_del, UnitRange(xi[nw,2]+1, nx))
  end

  return x_del
end

#=
OUTPUTS
xi  start and end indices of X to keep
W   truncated time windows
=#
function sync_t(t::Array{Int64, 2}, Δ::Int64, t_min::Int64, t_max::Int64)
  (t_max == 0) && (t_max = typemax(Int64))
  xi = x_inds(t)
  W = t_win(t, Δ)
  nw = size(W, 1)
  wi = zeros(Int64, nw, 2)
  wk = Array{Bool, 1}(undef, nw)
  fill!(wk, true)
  for i in 1:nw
    if (W[i,2] < t_min) || (W[i,1] > t_max)
      wk[i] = false
      continue
    end

    # Find last index ≤ t_max
    if W[i,2] > t_max
      k = xi[i,2]
      v = W[i,2]
      while v > t_max
        v -= Δ
        k -= 1
      end
      W[i,2] = v
      xi[i,2] = k
    end

    # Find first index ≥ t_min
    if W[i,1] < t_min
      j = xi[i,1]
      v = W[i,1]
      while v < t_min
        j += 1
        v += Δ
      end
      W[i,1] = v
      xi[i,1] = j
    end
  end

  # eliminated windows
  if minimum(wk) == false
    W = W[wk, :]
    xi = xi[wk, :]
  end
  return xi, W
end
# ==========

@doc """
sync!(S::GphysData)

Synchronize the start times of all data in S to begin at or after the last
start time in S.

sync!(S[, s=TS, t=TT, pad=false, v=V])

Synchronize all data in S to start no earlier than `TS` and terminate no later
than `TT`, with verbosity level `V`.

By default, a channel with mean `μᵢ = mean(S.x[i])` that begins after `TS` is
prepended with `μᵢ` to begin exactly at `TS`; similarly, if keyword `t` is used,
`μᵢ` is appended so that data ends at `TT`. If `pad=false`, channels that begin
after `TS` or end before `TT` are not extended in either direction.

For regularly-sampled channels, gaps between the specified and true times
are filled with the mean; this isn't possible with irregularly-sampled data.

#### Specifying start time (`s=`)
* s="last": (Default) sync to the last start time of any channel in `S`.
* s="first": sync to the first start time of any channel in `S`.
* A numeric value is treated as an epoch time (`?time` for details).
* A DateTime is treated as a DateTime. (see Dates.DateTime for details.)
* Any string other than "last" or "first" is parsed as a DateTime.

#### Specifying end time (`t=``)
* t="none": (Default) end times are not synchronized.
* t="last": synchronize all channels to end at the last end time in `S`.
* t="first" synchronize to the first end time in `S`.
* numeric, datetime, and non-reserved strings are treated as for `s=`.

See also: `TimeSpec`, `Dates.DateTime`, `parsetimewin`

!!! warning

    `sync!` calls `prune!`; empty channels will be deleted.
""" sync!
function sync!(S::GphysData;
                s::Union{String,DateTime}="last",
                t::Union{String,DateTime}="none",
                pad::Bool=true,
                v::Integer=KW.v,
                )

  # Delete empty traces
  prune!(S)                                         # delete empty channels
  S.n == 0 && return nothing                        # pointless to continue

  do_end = t=="none" ? false : true
  proc_str = string("sync!(S, s = \"", s, "\", t = \"", t, "\")")

  # Do not edit order of operations -------------------------------------------
  start_times = zeros(Int64, S.n)
  if do_end
    end_times = zeros(Int64, S.n)
  end
  fs = S.fs
  irr = falses(S.n)
  z = zero(Int64)

  # (1) Determine start and end times
  for i = 1:S.n
    start_times[i] = starttime(S.t[i], S.fs[i])
    if do_end
      end_times[i] = endtime(S.t[i], S.fs[i])
    end
    if fs[i] == 0.0
      irr[i] = true
    end
  end

  # (2) Determine earliest start and latest end
  t_start = get_sync_t(s, start_times)
  t_end = 0
  t_str = "none"
  if do_end
    t_end = get_sync_t(t, end_times)
    (t_end > t_start) || error("No time overlap with given start & end times!")
    t_str = string(u2d(t_end*μs))
    if v > 0
      @info(@sprintf("Synchronizing %.2f seconds of data\n", (t_end - t_start)*μs))
      if v > 1
        @info(string("t_start = ", u2d(t_start*μs)))
        @info(string("t_end = ", t_str))
      end
    end
  elseif v > 0
    @info(string("Synchronizing to start at ", u2d(t_start*μs)))
  end

  # (3) Synchronization to t_start (and maybe t_end)
  dflag = falses(S.n)
  for i = 1:S.n

    # non-timeseries data
    if fs[i] == 0.0
      t = view(S.t[i], :, 2)
      nt = length(t)
      k = get_sync_inds(t, do_end, t_start, t_end)
      nk = length(k)
      if nk ≥ nt
        dflag[i] = true
        proc_note!(S, i, proc_str, "synchronize, :x unchanged")
        continue
      else
        proc_note!(S, i, proc_str, string("synchronize, deleted ", nk, " samples from :x"))
      end
      deleteat!(S.x[i], k)
      ti = collect(1:nt)
      deleteat!(ti, k)
      S.t[i] = S.t[i][ti,:]
      S.t[i][:,1] .= 1:length(S.x[i])

    # timeseries data
    else
      sync_str = Array{String, 1}(undef, 0)
      desc_str = ""

      # truncate X to values within bounds
      Δ = round(Int64, sμ/fs[i])
      (xi, W) = sync_t(S.t[i], Δ, t_start, t_end)
      if isempty(W)
        dflag[i] = true
        continue
      else
        nx = length(S.x[i])
        x_del = get_del_ranges(xi, nx)
        nr = size(x_del, 1)
        for j in nr:-1:1
          if last(x_del[j]) == nx
            resize!(S.x[i], first(x_del[j])-1)
          else
            deleteat!(S.x[i], x_del[j])
          end
        end

        #= length 0 traces _can_ happen with resampled :x
        where (lx*f_rat) < 0.5; requries short series of
        high-frequency data resampled to too-low fs
        =#
        if length(S.x[i]) == 0
          dflag[i] = true
          continue
        end

        if length(S.x[i]) < nx
          push!(sync_str, string("deleted ", nx-length(S.x[i]), " samples from :x"))
        end

        # prepend points to time series data that begin late
        T = eltype(S.x[i])
        μ = T(mean(S.x[i]))
        ni = div(start_times[i] - t_start, Δ)

        sort_segs!(W)
        if (ni > 0) && (pad == true)
          prepend!(S.x[i], ones(T, ni).*μ)
          W[1] -= ni*Δ

          # logging
          push!(sync_str, string("prepended ", ni, " samples to :x."))
        end

        # append points to time series data that end early
        if do_end
          nj = div(t_end - W[end], Δ)
          if (nj > 0) && (pad == true)
            nx = length(S.x[i])
            resize!(S.x[i], nx+nj)
            V = view(S.x[i], nx+1:nx+nj)
            fill!(V, μ)
            W[end] += nj*Δ

            # logging
            push!(sync_str, string("appended ", nj, " samples to :x."))
          end
        end

        # last step
        S.t[i] = w_time(W, Δ)

        # logging
        if length(sync_str) > 0
          desc_str = string(", ", join(sync_str, ";"))
        end
        proc_note!(S, i, proc_str, desc_str)
      end
    end
  end
  del_flagged!(S, dflag, "length 0 after sync")
  return nothing
end

function sync!(C::SeisChannel;
                pad::Bool=true,
                s::Union{String,DateTime}="last",
                t::Union{String,DateTime}="none",
                v::Integer=KW.v )
  S = SeisData(C)
  sync!(S, pad=pad, s=s, t=t, v=v)
  return nothing
end

@doc (@doc sync!)
function sync(S::GphysData;
                pad::Bool=true,
                s::Union{String,DateTime}="last",
                t::Union{String,DateTime}="none",
                v::Integer=KW.v )

  T = deepcopy(S)
  sync!(T, pad=pad, s=s, t=t, v=v)
  return T
end

function sync(C::SeisChannel;
                pad::Bool=true,
                s::Union{String,DateTime}="last",
                t::Union{String,DateTime}="none",
                v::Integer=KW.v )

  U = deepcopy(C)
  S = SeisData(U)
  sync!(S, pad=pad, s=s, t=t, v=v)
  return S[1]
end

ext: jl


text:
using Test

@testset "Tests on LoadSave" begin

    save_filename = "testsave"
    rm(joinpath(save_filename * ".h5"), force=true) # delete pre-existing file

    param  = ( h₀ = 1, h₁ = 1.5, h₂ = 2, ϵ = 1, μ = 1, N  = 2^6)
    (η,u,v,mesh,par) = CnoidalWaveSerreGreenNaghdi(param)
    param=merge(par,param,( L=par.λ, T = 1, dt = 0.1))
    init = Init(mesh,η,v;label = "SGN cnoidal wave")
    x = mesh.x
    dump(save_filename, x, init)

    model = Airy(param)
    problem = Problem( model, init, param )
    solve!(problem; verbose=false)

    dump(save_filename, problem.data)

    loaded_init = load_init(save_filename)
    @test loaded_init.label == "SGN cnoidal wave"
    @test all( loaded_init.η(x) .≈ init.η(x) )
    @test all( loaded_init.v(x) .≈ init.v(x) )

    loaded_data = load_data(save_filename)
    @test loaded_data == problem.data

    rm(joinpath(save_filename * ".h5"), force=true) # delete pre-existing file

    dump(save_filename, problem)
    loaded_data = load_data(save_filename)
    @test loaded_data == problem.data

    new_problem = Problem( model, init, param )
    @test new_problem !== problem
    load_data!(save_filename, new_problem)
    @test new_problem == problem

end


text:
module DummyPackage

    version() = "v2"

end # module


text:
module io

self.IN = const FileReader('__STDIN__')
self.OUT = const FileWriter('__STDOUT__', True)
self.ERROR = const FileWriter('__STDERR__', True)

class FileInternal {
  field file
  new(fn, rw, a, binary) {
    mode = rw.copy()
    if binary {
      mode.append('b')
    }
    if a {
      mode.append('+')
    }
    file = File__(fn, mode)
    if ~file.success {
      raise Error(concat('Failed to open file. FileInternal(',
                      fn, ',', rw, ',', a, ',', binary,')'))
    }
  }
  method rewind() file.rewind__()
  method gets(n) file.gets__(n)
  method getline() file.getline__()
  method getlines() file.getall__()
  method puts(s) file.puts__(s)
  method close() file.close__()
}

class FileReader {
  field fi
  new(fn, binary=False) {
    fi = const FileInternal(fn, 'r', False, binary)
  }
  method rewind() fi.rewind()
  method gets(n) fi.gets(n)
  method getline() fi.getline()
  method getlines() fi.getlines()
  method getall() fi.getall()
  method close() fi.close()
}

class FileWriter {
  field fi
  new(fn, append=False, binary=False) {
    fi = const FileInternal(fn, 'w', append, binary)
  }
  method rewind() fi.rewind()
  method write(s) fi.puts(s)
  method writeln(s) {
    fi.puts(s)
    fi.puts('\n')
  }
  method close() fi.close()
}

def fprint(f, a) {
  f.write(str(a))
}

def fprintln(f, a) {
  f.writeln(str(a))
}

def print(a) {
  fprint(OUT, a)
}

def println(a) {
  fprintln(OUT, a)
}


text:
# Autogenerated wrapper script for librealsense_jll for x86_64-apple-darwin
export librealsense2

using libusb_jll
JLLWrappers.@generate_wrapper_header("librealsense")
JLLWrappers.@declare_library_product(librealsense2, "@rpath/librealsense2.2.41.dylib")
function __init__()
    JLLWrappers.@generate_init_header(libusb_jll)
    JLLWrappers.@init_library_product(
        librealsense2,
        "lib/librealsense2.2.41.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@field SpatialOrder = (1, 2) # default value for SpatialOrder (xy)

"""
Do a custom convert for SpatialOrder. Partial is left untyped, since this should apply
to all SpatialOrders in all parent composables.
If you need to overwrite behaviour if SpatialOrder is part of another composable, you just need to type parent
"""
function convertfor(Partial, ::Type{SpatialOrder}, data, value)
    @needs data: image = ImageData
    N = ndims(image)
    if isa(value, Tuple) &&
        if eltype(value) == Int || length(value) != N
            throw(UsageError(SpatialOrder, value))
        end
        return value
    end
    str = if isa(value, Symbol)
        string(value)
    elseif isa(value, String)
        value
    else
        throw(UsageError(SpatialOrder, value))
    end
    if length(str) != N
        throw(UsageError(SpatialOrder, value))
    end
    ntuple(length(str)) do i
        idx = findfirst(('x', 'y', 'z', 't'), str[i])
        if idx == 0
            throw(UsageError(SpatialOrder, value))
        end
        idx
    end
end

"""
Determines what the order of the dimensions is.
Can be a symbol or string like: xyz, yx, etc,
or a tuple with the indices of the dimensions, exactly like how you would pass it
to permutedims.
"""
SpatialOrder

@field ImageData
@field Ranges

"""
Ranges indicate, on what an otherwise dimensionless visualization should be mapped.
E.g. use Ranges to indicate that an image should be mapped to a certain range.
"""
Ranges

function default(Partial, ::Type{Ranges}, data)
    @needs data: image = ImageData
    s = get(data, SpatialOrder) # if SpatialOrder in x, gets that, if not gets default(x, SpatialOrder)
    (0:size(image, s[1]), 0:size(image, s[2]))
end

@composed type Image
    ImageData
    Ranges
    SpatialOrder::NTuple{2, Int}
end


text:
function inline_getproperty_annotation(name, type, sname)
    (name == "sType" || is_enum(type) || is_bitmask(type) || is_base_type(type) || is_base_type(converted_type(fieldtype_transform(name, type, sname))) || is_optional_parameter(name, sname)) && return "noproperty"
    is_handle(type) && return "handle"
    is_vulkan_type(type) && return "vulkan_type"
    "unknown"
end

function annotate_pass(statements, pass)
    body_statements = getproperty.(statements, :body)
    Statement.(rstrip.(body_statements, '\n') .* "    # $(typeof(pass).name)" .* map(x -> endswith(x, "\n") ? "\n" : "", body_statements), getproperty.(statements, :assigned_id))
end


text:
function SpinLab(;frame=true)
## ASSETS
path = @__DIR__
assets = AssetRegistry.register(dirname(path*"/ui/assets/"))
scripts = AssetRegistry.register(dirname(path*"/ui/scripts/"))
css = AssetRegistry.register(dirname(path*"/ui/css/"))

logo = joinpath(assets, "Logo.png")
loading = joinpath(assets, "Loading.gif")
popper = joinpath(scripts, "popper.min.js")
bsjs = joinpath(scripts, "bootstrap.min.js")
bscss = joinpath(css,"bootstrap.min.css")
jquery = joinpath(scripts,"jquery-3.4.1.slim.min.js")
# KaTeX
katexrender = joinpath(scripts, "auto-render.min.js")
katexjs = joinpath(scripts,"katex.min.js")
katexcss = joinpath(css,"katex.min.css")
# User defined
customcss = joinpath(css,"custom.css")
customjs = joinpath(scripts,"custom.js")
customjs2 = joinpath(scripts,"custom2.js")
# Custom icons
icons = joinpath(css,"icons.css")
# Others
imphantom = joinpath(assets, "phantom.png")
imscanner = joinpath(assets, "scanner.png")
impulses = joinpath(assets, "pulses.png")
## WINDOW
global w = Blink.Window(Dict(
    "title"=>"SpinLab",
    # "autoHideMenuBar"=>false,
    "frame"=>frame, #removes title bar
    "node-integration" => true,
    :icon=>path*"/ui/assets/Logo_icon.png"
    ),async=false);
## LOADING BAR
loadbar = """
"""
## NAV BAR
navbar = open(f->read(f, String), path*"/ui/html/navbar.html")
navbar = replace(navbar, "LOGO"=>logo)
## CONTENT
index = open(f->read(f, String), path*"/ui/html/index.html")
## FOOTER
footer = open(f->read(f, String), path*"/ui/html/footer.html")
## CSS
loadcss!(w, bscss)
loadcss!(w, customcss)
# KATEX
loadcss!(w, katexcss)
loadjs!(w, katexjs)
loadjs!(w, katexrender)
# JQUERY, BOOSTRAP JS
# mathjax = "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.7/MathJax.js?config=TeX-AMS-MML_SVG"
# loadjs!(w, mathjax)
loadjs!(w, customjs)    #must be before jquery
loadjs!(w, jquery)
loadjs!(w, customjs2)   #must be after jquery
loadjs!(w, popper)
loadjs!(w, bsjs)        #after jquery
# LOAD IMAGES AND ICONS
loadcss!(w, icons)
index = replace(index, "PHANTOM"=>imphantom)
index = replace(index, "SCANNER"=>imscanner)
index = replace(index, "PULSES"=>impulses)
## MENU FUNCTIONS
handle(w, "index") do args...
    @js_ w (@var loading = $loadbar; document.getElementById("content").innerHTML=loading)
    body!(w,*(navbar,index,footer))
end
handle(w, "pulses") do args...
    @js_ w (@var loading = $loadbar; document.getElementById("content").innerHTML=loading)
    include(path*"/ui/PulsesGUI.jl")
end
handle(w, "phantom") do args...
    @js_ w (@var loading = $loadbar; document.getElementById("content").innerHTML=loading)
    include(path*"/ui/PhantomGUI.jl")
end
handle(w, "sim") do args...
    @js_ w (@var loading = $loadbar; document.getElementById("content").innerHTML=loading)
    include(path*"/ui/SimulatorGUI.jl")
end
handle(w, "recon") do args...
    @js_ w (@var loading = $loadbar; document.getElementById("content").innerHTML=loading)
    include(path*"/ui/ReconGUI.jl")
end
handle(w, "simulate") do args...
    @info "Running simulation..."
    @js_ w (@var loading = $loadbar; document.getElementById("simulate").innerHTML=loading)
    Δt = 4e-6 #<- simulate param
    t = collect(0:Δt:MRIsim.dur(seq))
    Nphant, Nt = prod(size(phantom)), length(t)
    N_parts = floor(Int, Nphant*Nt/2.7e6)
    println("Dividing simulation in Nblocks=$N_parts")
    S = @time MRIsim.run_sim_time_iter(phantom,seq,t;N_parts)
    global signal = S ./prod(size(phantom)) #Acquired data
    #Recon, will be replaced by call to MRIReco.jl
    S = nothing #remove aux signal S
    Nx = Ny = 99 #hardcoded by now
    global kdata = reshape(signal,(Nx,Ny)) #Turning into kspace image
    global kdata[:,2:2:Ny] = kdata[Nx:-1:1,2:2:Ny] #Flip in freq-dir for the EPI
    global kdata = convert(Array{Complex{Float64},2},kdata)

    @js_ w document.getElementById("simulate").innerHTML="""<button type="button" onclick='Blink.msg("simulate", 1)' class="btn btn-primary btn-lg btn-block">Run simulation!</button>"""
    @js_ w Blink.msg("recon", 0)
end
handle(w, "close") do args...
    global phantom = nothing
    global seq = nothing
    global scanner = nothing

    close(w)
end
## PRINTING INFO
#PHANTOM init
@info "Loading Phantom (default)"
global phantom = brain_phantom2D(;axis="coronal")
println("Phantom object \"$(phantom.name)\" successfully loaded!")
#SEQ init
@info "Loading Sequence (default) "
B1 = 6e-6; durRF = π/2/(2π*γ*B1) #90-degree hard excitation pulse
EX = PulseDesigner.RF_hard(B1, durRF)
Gmax = 60e-3
EPI,_,_,_ = PulseDesigner.EPI_base(40/100, 99, 4e-6, Gmax)
TE = 25e-3
d = delay(TE-dur(EPI)/2-dur(EX))
DELAY = Sequence([d;d])
global seq = EX + DELAY + EPI
println("EPI successfully loaded! (TE = $(TE*1e3) ms)")
#Init
global scanner = []
global signal = 0
global kdata = [0.0im 0.; 0. 0.]
#GPUs
if has_cuda()
    @info "Loading GPUs"
    print_gpus()
end
#Update GUI's home
body!(w,*(navbar,index,footer),async=false)
nothing
end


text:
#!/usr/bin/julia

function uniform_random(a, b, seed)
    local t

    seed = 2045 * seed + 1
    #println("seed1:", seed)
    seed = seed - ( seed / 1048576) * 1048576
    #println("seed2:", seed)
    t = seed/1048576.0
    t = a + (b-a) * t

    return t
end

seed=13579

for i=1:10

    #println("seed:",seed)

    x=uniform_random(0, 1.0, seed)
    print(x,'\n')
end


text:
"""
Computes Jacobi polynomial and derivatives of Jacobi polynomials.

This function computes the jacobi polynomial of degree `n` with weights
`a` and `b` at point x (\$P_n^{a,b}(x)\$). There are several variants of the function
with default values:

 * `jacobi(x, n)` (`a` and `b` are both zero)
 * `jacobi(x, n, a)` (`b` is zero)
 * `jacobi(x, n, a, b) where `x` is an array
 * `jacobi!(x, n, a, b, y)` modifying array version of the function

The derivative of Jacobi polynomials (\$\\frac{dP_n^{a,b}(x)}{dx}\$)are computed with functions that add a `d` in front of its name:

 * `djacobi`
 * `djacobi!`

### Examples
```julia
using Jacobi

x = 0.3
a = 0.2
b = 0.1
m = 5
y = jacobi(x, m, a, b)
x1 = linspace(-1,1,21)
y1 = jacobi(x1, m, a, b)
jacobi!(x1, m, a, b, y1)

dy = djacobi(x, m, a, b)
dy1 = djacobi(x1, m, a, b)
djacobi!(x1, m, a, b, dy1)

```

"""
function jacobi(x, n, a, b)
    ox = one(x)
    zx = zero(x)
    if n==0
        return ox
    elseif n==1
        return ox/2 * (a - b + (a + b + 2)*x)
    end

    p0 = ox
    p1 = ox/2 * (a - b + (a + b + 2)*x)
    p2 = zx;

    for i = 1:(n-1)
	a1 = 2*(i+1)*(i+a+b+1)*(2*i+a+b);
	a2 = (2*i+a+b+1)*(a*a-b*b);
	a3 = (2*i+a+b)*(2*i+a+b+1)*(2*i+a+b+2);
	a4 = 2*(i+a)*(i+b)*(2*i+a+b+2);
	p2 = ox/a1*( (a2 + a3*x)*p1 - a4*p0);

        p0 = p1
        p1 = p2
    end

    return p2
end
jacobi(x, n) = jacobi(x, n, zero(x), zero(x))
jacobi(x, n, a) = jacobi(x, n, a, zero(x))

djacobi(x, n, a, b) =  one(x)/2 * (a + b + n + 1) * jacobi(x, n-1, a+1, b+1)

djacobi(x, n) = djacobi(x, n, zero(x), zero(x))
djacobi(x, n, a) = djacobi(x, n, a, zero(x))

eps1(::Type{T}) where {T<:AbstractFloat} = eps(T)
eps1(::Type{Complex{T}}) where {T<:AbstractFloat} = eps(T)

"""
Compute the zeros of Jacobi polynomials

This function computes the zeros of Jacobi polynomials:

 * \$P_m^{a,b}(x) = 0 \$

The `jacobi_zeros!` is the modifying version and the memory where the zeros
will be stored are preallocated. The non-modifying version, `jacobi_zeros`
allocates a the memory and calls the modifying version. The function `legendre_zeros`
compute the zeros of Legendre polynomials (`a = b = 0`)

### Examples
```
using Jacobi
z = jacobi_zeros(7, 0.3, 0.2)
```
"""
function jacobi_zeros!(m, alpha, beta, x::AbstractArray{T}) where {T<:Number}

    o = one(T)
    z = zero(T)

    a = convert(T,alpha)
    b = convert(T,beta)

    MAXITER = 500
    EPS::T = 100 * eps1(T)
    local i; local k; local iter=0

    for k = 1:m
        # Initial guess.
        r = -cos( (2k-o)/(2*m) * pi)
        if (k > 1)
            r = (r + x[k-1]) / 2
        end
        iter = 0
        while(true)
            s = z
            for i = 1:(k-1)
                s += o/(r - x[i])
            end

            poly = jacobi(r, m, a, b)
            delta = -poly / (djacobi(r, m, a, b) - poly*s)

            r += delta
            iter += 1

            if iter > MAXITER
                throw("Program did not converge")
            end

            if abs(delta) < abs(EPS)
                break
            end
        end
        x[k] = r
    end

    return x

end

function jacobi_zeros(m, a, b, ::Type{T}=Float64) where {T<:Number}
    jacobi_zeros!(m, a, b, zeros(T,m))
end

jacobi_zeros(m) = jacobi_zeros(m, 0.0, 0.0)
jacobi_zeros(m, a) = jacobi_zeros(m, a, zero(a))

@doc (@doc jacobi) djacobi
@doc (@doc jacobi) djacobi!
@doc (@doc jacobi) jacobi!
@doc (@doc jacobi_zeros!) jacobi_zeros


text:
module PlotRecipes

import PyPlot
PyPlot.svg(true)
using Plots
using Statistics
#using Plotly
using PlotlyJS
using Clustering
using SparseArrays

using RealNeuralNetworks.Neurons
using RealNeuralNetworks.Neurons.Segments
using RealNeuralNetworks.Neurons.Segments.Synapses

function plot_synapse_distributions( cellList;
        synapseDistribution = map(get_synapse_to_soma_path_length_lists, cellList),
        saveName = joinpath(homedir(), "synapse_distribution.svg"),
        prerange=(0,50), postrange=(0,100), bins=30)
    fig = PyPlot.figure()
    for (i,cellId) in enumerate(cellList)
        @show i
        fig[:add_subplot](4,2,(i-1)*2+1)
        PyPlot.plt[:hist](synapseDistribution[i][1]./1000; range=prerange, bins=bins)
        PyPlot.plt[:xlabel]("distrance from presynapse to soma ($cellId)")

        fig[:add_subplot](4,2,(i-1)*2+2)
        PyPlot.plt[:hist](synapseDistribution[i][2]./1000; range=postrange, bins=bins)
        PyPlot.plt[:xlabel]("distrance from postsynapse to soma ($cellId)")

    end
    PyPlot.savefig(saveName)
end

function reorder_distance_matrix(distanceMatrix::Array, clust::Hclust)
	# rearrange distance matrix according to order
	reorderedDistanceMatrix = zeros(distanceMatrix)
	for (i,o1) in enumerate(clust.order)
		for (j, o2) in enumerate(clust.order)
			reorderedDistanceMatrix[i,j] = distanceMatrix[o1,o2]
		end
    end
	reorderedDistanceMatrix
end

function plot_connectivity_matrix(orderedConnMatrix::SparseMatrixCSC{T};
                                    title="connectivity matrix",
                                    max_size::Real=10) where T
    X,Y,V = SparseArrays.findnz(orderedConnMatrix)
    # transform to float and normalize to max_size
    V = Float32.(V)
    V ./= maximum(V) / Float32(max_size)

    trace = PlotlyJS.scatter(; x=X, y=Y,
        mode="markers",
        marker=attr(size=V, line_width=0, opacity=0.4))
    data = [trace,]

    layout = Layout(; title=title,
        width=800,
        height=800,
        xaxis = attr(title="presynapse", range=(1, size(orderedConnMatrix,1))),
        yaxis= attr(title="postsynapse", range=(1, size(orderedConnMatrix,2)))
    )
    return PlotlyJS.plot(data, layout)
end

function plot(neuron::Neuron; nodeStep::Integer=1, semantic::Bool=true,
                showgrid::Bool=true, showSynapse::Bool=true)
    traces = PlotlyJS.GenericTrace[]
    # plot soma
    root = Neurons.get_root_node(neuron)
    root_trace = PlotlyJS.scatter3d(;x=[root[1]], y=[root[2]], z=[root[3]],
                                    mode="markers", marker_size=7)
    push!(traces, root_trace)

    for (i,segment) in enumerate(neuron)
        nodeList = Neurons.get_segment_node_list(neuron, i)
        x = map(n->n[1], nodeList[1:nodeStep:end]) |> Vector{Float32}
        y = map(n->n[2], nodeList[1:nodeStep:end]) |> Vector{Float32}
        z = map(n->n[3], nodeList[1:nodeStep:end]) |> Vector{Float32}
        if semantic
            class = Neurons.Segments.get_class(segment)
            if class == Neurons.Segments.AXON_CLASS
                color = "rgb(255,0,0)"
            elseif class == Neurons.Segments.DENDRITE_CLASS
                color = "rgb(0,0,255)"
            else
                color = "rgb(127,127,127)"
            end
            segmentTrace = PlotlyJS.scatter3d(;x=x,y=y,z=z, mode="lines",
                                                line_color=color)
        else
            segmentTrace = PlotlyJS.scatter3d(;x=x,y=y,z=z, mode="lines")
        end
        push!(traces, segmentTrace)

        if showSynapse
            preSynapseList = Neurons.Segments.concatenate_pre_synapses(segment)
            postSynapseList = Neurons.Segments.concatenate_post_synapses(segment)
            preTrace = PlotlyJS.scatter3d(;
                                x=map(c->Synapses.get_psd_coordinate(c)[1], preSynapseList),
                                y=map(c->Synapses.get_psd_coordinate(c)[2], preSynapseList),
                                z=map(c->Synapses.get_psd_coordinate(c)[3], preSynapseList),
                                mode="markers", marker_color="rgb(179,66,244)", marker_size=2)
            postTrace = PlotlyJS.scatter3d(;
                                x=map(c->Synapses.get_psd_coordinate(c)[1], postSynapseList),
                                y=map(c->Synapses.get_psd_coordinate(c)[2], postSynapseList),
                                z=map(c->Synapses.get_psd_coordinate(c)[3], postSynapseList),
                                mode="markers", marker_color="rgb(66,223,244)", marker_size=2)
            push!(traces, preTrace)
            push!(traces, postTrace)
        end
    end
    layout = PlotlyJS.Layout(; showlegend=false, showticklabels=false,
                            scene=attr(
                                xaxis=attr(showgrid=showgrid, showticklabels=false),
                                yaxis=attr(showgrid=showgrid, showticklabels=false),
                                zaxis=attr(showgrid=showgrid, showticklabels=false))
    )
    PlotlyJS.plot(traces, layout)
end

function plot_v2(neuron::Neuron; nodeStep::Integer=10)
    #using PyPlot
    #PyPlot.pygui(true)

    # plot soma
    root = Neurons.get_root_node(neuron)
    PyPlot.scatter3D([root[1]], [root[2]], [root[3]])

    for segment in neuron
        nodeList = Neurons.Segments.get_node_list(segment)
        x = map(n->n[1], nodeList[1:nodeStep:end])
        y = map(n->n[2], nodeList[1:nodeStep:end])
        z = map(n->n[3], nodeList[1:nodeStep:end])
        PyPlot.plot(x,y,z)
    end
end

function plot_v1(neuron::Neuron; nodeStep::Integer=10)
    segmentList = neuron.segmentList
    plotly()
    for branch in segmentList
        nodeList = Neurons.Segments.get_node_list(branch)
        x = map(n->n[1], nodeList[1:nodeStep:end])
        y = map(n->n[2], nodeList[1:nodeStep:end])
        z = map(n->n[3], nodeList[1:nodeStep:end])
        plot!(x,y,z) #color=rand(Colors.RGB))
    end

    root = Neurons.get_root_node(neuron)
    plot!([root[1]], [root[2]], [root[3]], m=(2, :circle), leg=false)
end

"""
    plot_maximum_intensity_projection(vol::Array; colorMap="gray")
color map could be "jet" if you like colorful display.
"""
function plot_maximum_intensity_projection(vol::Array{T,3}; colorMap="gray") where T
    fig = PyPlot.figure()
    fig[:add_subplot](2,2,1)
    xy = maximum(vol, 3)[:,:,1]
    PyPlot.imshow(xy, colorMap)
    # colorbar()
    fig[:add_subplot](2,2,2)
    xz = maximum(vol, 2)[:,1,:]
    PyPlot.imshow(xz, colorMap)
    fig[:add_subplot](2,2,3)
    yz = maximum(vol, 1)[1,:,:] |> rotl90
    PyPlot.imshow(yz, colorMap)
    PyPlot.colorbar()
end

function indexmap(x::Vector)
    ret = Dict()
    for (i,v) in enumerate(x)
        ret[v] = i
    end
    ret
end

function hclustplot(hc::Hclust, useheight::Bool)
    o = indexmap(hc.order)
    n = [x for x in 1:length(o)]

    pos = treepositions(hc, useheight)

    xs = []
    ys = []
    for i in 1: size(hc.merges, 1)
        x1 = pos[hc.merges[i,1]][1]
        x2 = pos[hc.merges[i,2]][1]
        append!(xs, [x1,x1,x2,x2])

        y1 = pos[hc.merges[i,1]][2]
        y2 = pos[hc.merges[i,2]][2]
        useheight ? h = hc.heights[i] : h = 1
        newy = maximum([y1,y2]) + h
        append!(ys, [y1,newy,newy,y2])
    end
    return (reshape(xs, 4, size(hc.merges, 1)), reshape(ys, 4, size(hc.merges, 1)))
end

function treepositions(hc::Hclust, useheight::Bool)
    order = indexmap(hc.order)
    positions = Dict{}()
    for (k,v) in order
        positions[-k] = (v, 0)
    end
    for i in 1:size(hc.merges,1)
        xpos = mean([positions[hc.merges[i,1]][1], positions[hc.merges[i,2]][1]])
        if hc.merges[i,1] < 0 && hc.merges[i,2] < 0
            useheight ? ypos = hc.heights[i] : ypos = 1
        else
            useheight ? h = hc.heights[i] : h = 1
            ypos = maximum([positions[hc.merges[i,1]][2], positions[hc.merges[i,2]][2]]) + h
        end

        positions[i] = (xpos, ypos)
    end
    return positions
end

function plot(clust::Hclust)
	Plots.gr()
	#Plots.plotlyjs()
	Plots.plot(hclustplot(clust, true), seriestype=:path, color=:black,
    grid=false, legend=false) #,  xticks=classificationIdList[clust.order])
end

end # module


text:
# Solve the Landau-Zener problem i ψ' = H(t) ψ, with H(t) = [t 1;1 -t]

using Test
using StaticArrays, LinearAlgebra
using OrdinaryDiffEq, DiffEqBase

H(t) = -im*(@SMatrix [t 1;1 -t])

fun(ψ,p,t) = H(t)*ψ
fun_inplace(dψ,ψ,p,t) = (dψ .= H(t)*ψ)

T = 0.1
tspan = (0,T)
explicit = [Midpoint(),RK4(),DP5(),Tsit5(),Vern7()]
implicit_autodiff = [ImplicitEuler(),Trapezoid(),Kvaerno3(),Rosenbrock23()]
implicit_noautodiff = [ImplicitEuler(autodiff=false),Trapezoid(autodiff=false),Kvaerno3(autodiff=false),Rosenbrock23(autodiff=false)]

@testset "Complex Tests on Explicit Methods" for alg in explicit
  for f in (fun, fun_inplace)
    ψ0 = [1.0+0.0im; 0.0]
    prob = ODEProblem(f,ψ0,(-T,T))
    sol = solve(prob,alg)
    @test norm(sol(T)) ≈ 1 atol=1e-2
  end
  ψ0 = @SArray [1.0+0.0im; 0.0]
  prob = ODEProblem(fun,ψ0,(-T,T))
  sol = solve(prob,alg)
  @test norm(sol(T)) ≈ 1 atol=1e-2
end

@testset "Complex Tests on Implicit Autodiff Methods" for alg in implicit_autodiff
  @test_broken begin
    for f in (fun, fun_inplace)
      ψ0 = [1.0+0.0im; 0.0]
      prob = ODEProblem(f,ψ0,(-T,T))
      sol = solve(prob,alg)
      @test norm(sol(T)) ≈ 1 atol=1e-2
    end
    ψ0 = @SArray [1.0+0.0im; 0.0]
    prob = ODEProblem(fun,ψ0,(-T,T))
    sol = solve(prob,alg)
    @test norm(sol(T)) ≈ 1 atol=1e-2
  end
end

@testset "Complex Tests on Implicit Finite Diff Methods" for alg in implicit_noautodiff
    ψ0 = [1.0+0.0im; 0.0]
    prob = ODEProblem(fun_inplace,ψ0,(-T,T))
    sol = solve(prob,alg)
    @test norm(sol(T)) ≈ 1 atol=1e-2
end

@testset "Complex Tests on Implicit Finite Diff Out-of-place Methods" begin
  for alg in implicit_noautodiff
    ψ0 = [1.0+0.0im; 0.0]
    prob = ODEProblem(fun,ψ0,(-T,T))
    sol = solve(prob,alg)
    @test norm(sol(T)) ≈ 1 atol=1e-2
  end
end

@testset "Complex Tests on Implicit Finite Diff Out-of-place Methods SArray" for alg in implicit_noautodiff
  @test_broken begin
    for alg in implicit_noautodiff
      ψ0 = @SArray [1.0+0.0im; 0.0]
      prob = ODEProblem(fun,ψ0,(-T,T))
      sol = solve(prob,alg)
      @test norm(sol(T)) ≈ 1 atol=1e-2
    end
  end
end


text:
module ShaleDrillingLikelihood_Flow_Test

# using Revise
using ShaleDrillingLikelihood

using Test
using StatsFuns
using Random
using BenchmarkTools
using Base.Threads
# using Profile
# using ProfileView
using InteractiveUtils
using Base.Iterators

using Calculus

using ShaleDrillingLikelihood: TestDrillModel,
    flow,
    dflow,
    full_payoff,
    dflow!,
    flowdψ,
    SimulationDraw,
    reward,
    idx_drill_ψ, theta_drill_ψ,
    idx_drill_x, theta_drill_x,
    idx_drill_z, theta_drill_z,
    idx_drill_ρ, theta_drill_ρ,
    idx_drill_d, theta_drill_d,
    ObservationDrill,
    check_flow_grad

@testset "Drilling Model Flow Payoffs" begin

    @testset "TestDrillModel" begin
        model = TestDrillModel()
        coefnames(model)
        theta = [1.5, 2.0, -3.0, -0.6, 0.5]
        @test theta_drill_ψ(model, theta) == 1.5
        @test theta_drill_x(model, theta) == 2.0
        @test theta_drill_z(model, theta) == -3.0
        @test theta_drill_d(model, theta) == -0.6
        @test theta_drill_ρ(model, theta) == 0.5

        ichars = [tuple(x) for x in -1.0:0.5:1.0]
        zchars = [tuple(z) for z in -1.0:0.5:1.0]
        y = collect(0:2)
        x = collect(1:4)

        obstypes = product(0:2, product(ichars, zchars, y, x))

        for (d, obstype) in obstypes
            obs = ObservationDrill(model, obstype...)
            check_flow_grad(reward(model), d, obs, theta)
        end
    end

end # Drilling Model Flow Payoffs

end # module


text:
using LeastSquaresSVM
using Documenter

# * Build the complete documentation
makedocs(;
    modules=[LeastSquaresSVM],
    authors="Edwin Bedolla",
    repo="https://github.com/edwinb-ai/LeastSquaresSVM/blob/{commit}{path}#L{line}",
    sitename="LeastSquaresSVM",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://edwinb-ai.github.io/LeastSquaresSVM/stable/",
        assets=String[]
    ),
    pages=[
        "Home" => "index.md",
        "Examples" => [
            "example1.md",
            "example2.md",
            "example3.md",
            "example4.md"
        ],
        "Reference" => "reference.md",
    ]
)

deploydocs(;
    repo="github.com/edwinb-ai/LeastSquaresSVM.git",
    devbranch="main",
    push_preview=true
)

max_stars_repo_name: gmagannaDevelop/Elysivm

max_issues_repo_name: gmagannaDevelop/Elysivm

max_forks_repo_name: gmagannaDevelop/Elysivm