text:
"""
    @reduce() do (acc₁ [= init₁]; x₁), ..., (accₙ [= initₙ]; xₙ)
        ...
    end
    @reduce(acc₁ ⊗₁= x₁, ..., accₙ ⊗ₙ= xₙ)
    @reduce(acc₁ .⊗₁= x₁, ..., accₙ .⊗ₙ= xₙ)
    @reduce(acc₁ = ⊗₁(init₁, x₁), ..., accₙ = ⊗ₙ(initₙ, xₙ))
    @reduce(acc₁ .= (⊗₁).(init₁, x₁), ..., accₙ = (⊗ₙ).(initₙ, xₙ))

Declare how accumulators are updated in the sequential basecase and
how the resulting accumulators from two basecases are combined.

The arguments `accᵢ` and `xᵢ` must be symbols except for `xᵢ` of the
last three forms in which an expression can be used at `xᵢ`.

In the first form,

```julia
function ((acc₁, acc₂, ..., accₙ), (x₁, x₂, ..., xₙ))
    ...  # body of the `do` block
    return (acc₁, acc₂, ..., accₙ)
end
```

should be an associative function.

In the last three forms, every binary operation `⊗ᵢ` should be an associative
function.

If `initᵢ` is specified, the tuple `(init₁, init₂, ..., initₙ)` should
be the identify of the related associative function.  `accᵢ = initᵢ`
is evaluated for each basecase (each `Task`) in the beginning.

Consider a loop with the following form

```julia
@floop for ...
    # code computing (x₁, x₂, ..., xₙ)
    @reduce() do (acc₁ = init₁; x₁), ..., (accₙ = initₙ; xₙ)
        # code updating (acc₁, acc₂, ..., accₙ) using (x₁, x₂, ..., xₙ)
    end
end
```

This is converted to

```julia
acc₁ = init₁
...
accₙ = initₙ
for ...
    # code computing (x₁, x₂, ..., xₙ)
    # code updating (acc₁, acc₂, ..., accₙ) using (x₁, x₂, ..., xₙ)
end
```

for computing `(acc₁, acc₂, ..., accₙ)` of each basecase.  The
accumulators `accᵢ` of two basecases are combined using "code updating
`(acc₁, acc₂, ..., accₙ)` using `(x₁, x₂, ..., xₙ)`" where `(x₁, x₂,
..., xₙ)` are replaced with `(acc₁, acc₂, ..., accₙ)` of the next
basecase.  Note that "code computing `(x₁, x₂, ..., xₙ)`" is not used
for combining the basecases.

# Examples
```julia
@reduce() do (vmax=-Inf; v), (imax=0; i)
    if isless(vmax, v)
        vmax = v
        imax = i
    end
end

@reduce(s += y, p *= y)

@reduce(xs = append!!(EmptyVector(), x), ys = append!!(EmptyVector(), y))
```
"""
macro reduce(args...)
    # We insert a custom non-`Expr` object `ReduceOpSpec` in the AST
    # so that the argument to `@reduce` can be later analyzed by
    # `floop_parallel` below.  We do this because `floop_parallel` is
    # run after `macroexpand`.  We insert `throw` of `ReduceOpSpec` so
    # that unprocessed `ReduceOpSpec` becomes an error at runtime
    # rather than silently failing.
    :(throw($(ReduceOpSpec(collect(Any, args)))))
end
# TODO: detect free variables in `do` blocks

struct ReduceOpSpec
    args::Vector{Any}
    visible::Vector{Symbol}
end

ReduceOpSpec(args::Vector{Any}) = ReduceOpSpec(args, Symbol[])

"""
    @init begin
        pv₁ = init₁
        ...
        pvₙ = initₙ
    end

Initialize private variables `pvᵢ` with initializer expression `initᵢ` for
each task. This can be used for mutating objects in a data race-free manner.
"""
macro init(ex)
    :(throw($(initspec(ex))))
end

struct InitSpec
    expr::Expr
    lhs::Vector{Union{Symbol,Expr}}
    rhs::Vector{Any}
    visible::Vector{Symbol}
end

InitSpec(expr::Expr, lhs::Vector{Union{Symbol,Expr}}, rhs::Vector{Any}) =
    InitSpec(expr, lhs, rhs, Symbol[])

initspec(@nospecialize x) = invalid_at_init(x)
initspec(ex::Expr) = InitSpec(ex, collect_assignments(ex)...)

invalid_at_init(@nospecialize ex) =
    error("`@init` requires an assignment or a sequence of assignments; got:\n", ex)

# TODO: merge with `assigned_vars`?
collect_assignments(ex) = collect_assignments!(Union{Symbol,Expr}[], [], ex)
collect_assignments!(_, _, @nospecialize(x)) = invalid_at_init(x)
collect_assignments!(lhs, rhs, ::LineNumberNode) = (lhs, rhs)
function collect_assignments!(lhs, rhs, ex::Expr)
    @match ex begin
        Expr(:(=), l::Union{Symbol,Expr}, r) => begin
            push!(lhs, l)
            push!(rhs, r)
        end
        Expr(:block, args...) => begin
            for a in args
                collect_assignments!(lhs, rhs, a)
            end
        end
        # TODO: should we support other things, like side effects, here?
        _ => invalid_at_init(ex)
    end
    return (lhs, rhs)
end

""" Recursively copy `expr::Expr` but not `ReduceOpSpec` or `InitSpec`. """
copyexpr(expr::Expr) = Expr(expr.head, (copyexpr(a) for a in expr.args)...)
copyexpr(@nospecialize x) = x

function analyze_loop_local_variables!(spec::Union{ReduceOpSpec,InitSpec}, scopes)
    @assert isempty(spec.visible)
    append!(spec.visible, (var.name for sc in scopes for var in sc.bounds))
    unique!(spec.visible)
    nothing
end
function analyze_loop_local_variables!(expr::Expr, scopes)
    @match expr begin
        Expr(:scoped, sc, inner) => begin
            push!(scopes, sc)
            analyze_loop_local_variables!(inner, scopes)
            pop!(scopes)
            nothing
        end
        _ => begin
            for a in expr.args
                analyze_loop_local_variables!(a, scopes)
            end
        end
    end
end
analyze_loop_local_variables!(@nospecialize(_), _) = nothing

""" Fill `visible` field of `ReduceOpSpec` and `InitSpec` in-place. """
function fill_loop_local_variables!(expr)
    ex = copyexpr(expr)  # deep copy expr but not ReduceOpSpec or InitSpec
    analyze_loop_local_variables!(solve_from_local!(simplify_ex(ex)), Any[])
    return expr
end
#
# Currently, we overestimate variables visible to `@reduce` and `@init`. For
# example, we treat that `a` is visible to `init` clause of `@reduce` (i.e.,
# `f(a)`) in the following example:
#
#     @reduce y = op(f(a), x)
#     a = 1
#
# This is for (relatively easily) supporting the uses of `let` like this:
#
#     julia> let
#                let
#                    a = 1
#                end
#                a
#            end
#     ERROR: UndefVarError: a not defined
#     ...
#
#     julia> let
#                let
#                    a = 1
#                end
#                @show a
#                a = 0
#                a
#            end
#     a = 1
#     1
#
# The latter is the example where the visibility of `a` at `@show a` is
# "retroactively" changed by the assignment after `@show a`; i.e., we cannot
# use simple single forward-pass algorithm to compute the set of accessible
# variables at given points.

function unpack_kwargs(;
    otherwise = donothing,
    on_expr = otherwise,
    on_init = otherwise,
    kwargs...,
)
    @assert isempty(kwargs)
    return (otherwise, on_expr, on_init)
end

function on_reduce_op_spec(on_spec, ex; kwargs...)
    (otherwise, on_expr, on_init) = unpack_kwargs(; kwargs...)
    @match ex begin
        Expr(:call, throw′, spec::ReduceOpSpec) => on_spec(spec)
        Expr(:call, throw′, spec::InitSpec) => on_init(spec)
        Expr(head, args...) => begin
            new_args = map(args) do x
                on_reduce_op_spec(on_spec, x; kwargs...)
            end
            on_expr(head, new_args...)
        end
        _ => otherwise(ex)
    end
end

on_reduce_op_spec_reconstructing(on_spec, ex; otherwise = identity, on_init = otherwise) =
    on_reduce_op_spec(on_spec, ex; on_expr = Expr, otherwise = otherwise, on_init = on_init)

function floop_parallel(ctx::MacroContext, ex::Expr, simd, executor = nothing)
    if !isexpr(ex, :for, 2)
        error("expected a `for` loop; got:\n", ex)
    end
    iterspec, body = ex.args
    parallel_loop_ex = @match iterspec begin
        Expr(:block, loop_axes...) => begin
            rf_arg, coll = transform_multi_loop(loop_axes)
            as_parallel_loop(ctx, rf_arg, coll, body, simd, executor)
        end
        Expr(:(=), rf_arg, coll) => begin
            as_parallel_loop(ctx, rf_arg, coll, body, simd, executor)
        end
    end
    return parallel_loop_ex
end

function analyze_rf_args(ex::Expr)
    @assert isexpr(ex, :tuple)
    accs = []
    inits = []
    inputs = []
    for arg in ex.args
        @match arg begin
            Expr(:block, acc_init, x) || Expr(:block, acc_init, ::LineNumberNode, x) =>
                begin
                    push!(inputs, x)
                    @match acc_init begin
                        Expr(:(=), a, i) => begin
                            push!(accs, a)
                            push!(inits, i)
                        end
                        a => push!(accs, a)
                    end
                end
            Expr(:tuple, a, x) => begin
                throw(ArgumentError("got `($a, $x)` use `($a; $x)` instead"))
            end
        end
    end
    if !isempty(inits) && length(inits) != length(inputs)
        # TODO: Remove this restriction.  If not all `init`s are
        # specified, use it only for the empty case.
        msg = string(
            "[NOT IMPLEMENTED]",
            " Currently, initial value should be specified for all accumulators",
            " when it is specified for at least one accumulator.",
        )
        throw(ArgumentError(msg))
    end
    if isempty(inits)
        inits = nothing
    end
    return accs, inits, inputs
end

function verify_unique_symbols(all_vars, kind)
    if length(Set(all_vars)) != length(all_vars)
        counts = Dict{Any,Int}()
        for var in all_vars
            counts[var] = get(counts, var, 0) + 1
        end
        dups = sort!([(n, var) for (var, n) in counts if n > 1])
        msg = sprint() do io
            print(io, "Same $kind variable used multiple times.")
            for (n, var) in dups
                println(io)
                print(io, "* `", var, "` used ", n, " times")
            end
        end
        throw(ArgumentError(msg))
    end
end

# To allow something like `@reduce(c += 1)` and `@reduce(c = 0 + 1)`,
# assign the right (second) argument to a temporary variable:
function extract_pre_updates(raw_inputs)
    inputs = []
    pre_updates = []
    for x in raw_inputs
        if x isa Symbol
            push!(inputs, x)
        else
            @gensym input
            push!(pre_updates, :($input = $x))
            push!(inputs, input)
        end
    end
    inputs, pre_updates2 = uniquify_inputs(inputs)
    append!(pre_updates, pre_updates2)
    return (inputs, pre_updates)
end

function dematerialize_dotcall(ex)
    if is_dotcall(ex)
        :($_lazydotcall.($ex))
    else
        ex
    end
end

function _lazydotcall end
struct LazyDotCall{T}
    value::T
end
@inline Broadcast.broadcasted(::typeof(_lazydotcall), x) = LazyDotCall(x)
@inline Broadcast.materialize(x::LazyDotCall) = x.value

function uniquify_inputs(inputs)
    uniquified = empty(inputs)
    pre_updates = Expr[]
    seen = Set{eltype(inputs)}()
    for x in inputs
        if x in seen
            y = gensym(x)
            push!(pre_updates, :($y = $x))
        else
            push!(seen, x)
            y = x
        end
        push!(uniquified, y)
    end
    return uniquified, pre_updates
end

function inject_spec_expr_for_analysis(ex::Expr)
    function on_init(spec::InitSpec)
        quote
            $(spec.expr)
            throw($spec)
        end
    end
    function on_reduce(spec::ReduceOpSpec)
        # If we want to make
        #     @reduce(a = op(a0, x))
        #     @reduce(b = op(f(a), x))
        # work, we need to do something here. But, for the moment, let us
        # pretend that the output of `@reduce` is not visible outside the
        # `@reduce` block.
        :(throw($spec))
    end
    return on_reduce_op_spec_reconstructing(on_reduce, ex; on_init = on_init)
end

raw"""
    process_reduce_op_spec(spec) ->
        (; inits, accs, inputs, pre_updates, initializers, updaters)

Process an invocation of `@reduce() do ...` or `@reduce(spec₁, spec₂, ...)`
where each `specᵢ` is one of the following form

    acc ⊗= input
    acc .⊗= input
    acc = ⊗(init, input)

Roughly speaking, the processed results are injected to the loop body as follows
(see `rf_body_with_init` for exactly how this happens):

```
local $(accs...)
for x in xs
    ...

    # <- where `@reduce` is invoked
    $(pre_updates...)       # assign $inputs
    if $THIS_IS_THE_FIRST_ITERATION
        $(initializers...)  # "per-task" lazy initialization occurs here
    else
        $(updaters...)      # the main update logic; e.g., `acc = acc + x`
    end
end
```

To combine `accs`, following expressions are injected to the `combine_function`:

```julia
function $combine_function(all_accs, all_inputs)
    (..., ($(accs...),), ...) = all_accs      # results from left
    (..., ($(inputs...),), ...) = all_inputs  # results from right
    ...  # code from other reductions

    if $RIGHT_IS_EMPTY
        # do nothing; i.e., use the left result as-is
    else
        # Note: no $(pre_updates...) here
        if $LEFT_IS_EMPTY
            $(initializers...)
        else
            $(updaters...)
        end
    end

    ...  # code from other reductions
    return (..., ($(accs...),), ...)  # combined results
end
```
"""
function process_reduce_op_spec(
    spec::ReduceOpSpec,
)::NamedTuple{(:inits, :accs, :inputs, :pre_updates, :initializers, :updaters)}
    opspecs = spec.args::Vector{Any}

    if length(opspecs) == 1 && is_function(opspecs[1])
        # handle: @reduce() do ...
        rf_ex, = opspecs
        # rf_ex = :(((acc1; input1), ..., (accN; inputN)) -> rf_body)
        accs, inits, inputs = analyze_rf_args(rf_ex.args[1])
        rf_body = rf_ex.args[2]
        pre_updates = []
        updaters = [rf_body]
        initializers = default_initializers(inits, accs, inputs, updaters)
        return (;
            inits = inits,
            accs = accs,
            inputs = inputs,
            pre_updates = pre_updates,
            initializers = initializers,
            updaters = updaters,
        )
    end

    if all(is_dot_update, opspecs)  # [^broadcasted_reduction]
        # handle: @reduce(acc₁ .op₁= x₁, ..., accₙ .opₙ= xₙ)
        ops = [Symbol(String(x.head)[2:end-1]) for x in opspecs]
        accs = [x.args[1] for x in opspecs]
        inits = nothing
        (inputs, pre_updates) = extract_pre_updates([x.args[2] for x in opspecs])
        initializers = [:($a = $InitialValue($op)) for (a, op) in zip(accs, ops)]
        true
    elseif all(x -> isexpr(x, :(.=), 2), opspecs)
        if !all(x -> is_dotcall(x.args[2], 2), opspecs)
            msg = sprint() do io
                print(io, "`@reduce lhs .= rhs` syntax requires a binary dot call")
                print(io, " (e.g., `a .+ b` or `f.(a, b)`) on the rhs; got:")
                for (i, x) in pairs(opspecs)
                    if !is_dotcall(x.args[2], 2)
                        println(io)
                        print(io, i, "-th op spec: ", opspecs[i])
                    end
                end
            end
            error(msg)
        end
        # handle: @reduce(acc₁ .= op₁.(init₁, x₁), ..., accₙ .= opₙ.(initₙ, xₙ))
        accs, ops, inits, inputs =
            opspecs |>
            Map() do ex
                acc, rhs = ex.args
                if isexpr(rhs, :call, 3)
                    dotop, init, input = rhs.args
                    str = String(dotop)
                    @assert startswith(str, ".")
                    op = Symbol(str[2:end])
                else
                    @assert rhs.head == :. &&
                            length(rhs.args) == 2 &&
                            isexpr(rhs.args[2], :tuple, 2)
                    op = rhs.args[1]
                    init, input = rhs.args[2].args
                end
                (acc, op, init, input)
            end |>
            foldxl(ProductRF(push!!, push!!, push!!, push!!))
        (inputs, pre_updates) = extract_pre_updates(inputs)

        # Although we can't dematerialize `inits` since it could be hoisted out,
        # we can dematerialize `inits` in the initializers since we know that
        # they are followed by materializing `broadcast_inplace!!` call in the
        # updaters.
        initializers =
            [:($a = $(dematerialize_dotcall(x))) for (a, x) in zip(accs, inits)]
        true
    elseif any(is_dot_update, opspecs) || any(x -> isexpr(x, :(.=), 2), opspecs)
        error(
            "`@reduce` currently requires all syntaxes to be similar",
            " (e.g., don't mix `a .+= x` and `a .= 0 .+ x`); got:\n",
            join(opspecs, "\n"),
        )
    else
        false
    end && begin
        pre_updates = map(dematerialize_dotcall, pre_updates)
        updaters = map(ops, accs, inputs) do op, a, x
            broadcast_inplace!! = GlobalRef(@__MODULE__, :broadcast_inplace!!)
            :($a = $broadcast_inplace!!($op, $a, $x))
            # ^- mutate-or-widen version of `$a .= ($op).($a, $x)`
        end
        append!(initializers, updaters)
        return (;
            inits = inits,
            accs = accs,
            inputs = inputs,
            pre_updates = pre_updates,
            initializers = initializers,
            updaters = updaters,
        )
    end

    if all(is_rebinding_update, opspecs)
        # handle: @reduce(acc₁ op₁= x₁, ..., accₙ opₙ= xₙ)
        ops = [Symbol(String(x.head)[1:end-1]) for x in opspecs]
        accs = [x.args[1] for x in opspecs]
        inits = nothing
        (inputs, pre_updates) = extract_pre_updates([x.args[2] for x in opspecs])
    elseif all(x -> isexpr(x, :(=), 2) && isexpr(x.args[2], :call, 3), opspecs)
        # handle: @reduce(acc₁ = op₁(init₁, x₁), ..., accₙ = opₙ(initₙ, xₙ))
        ops = [x.args[2].args[1] for x in opspecs]
        accs = [x.args[1] for x in opspecs]
        inits = [x.args[2].args[2] for x in opspecs]
        (inputs, pre_updates) = extract_pre_updates([x.args[2].args[3] for x in opspecs])
    else
        error(join(vcat(["unsupported:"], opspecs), "\n"))
    end
    updaters = [:($a = $op($a, $x)) for (op, a, x) in zip(ops, accs, inputs)]
    initializers = default_initializers(inits, accs, inputs, updaters)
    return (;
        inits = inits,
        accs = accs,
        inputs = inputs,
        pre_updates = pre_updates,
        initializers = initializers,
        updaters = updaters,
    )
end
# [^broadcasted_reduction]: Broadcasted reduction of form `@reduce acc .⊕= $rhs`
# is handled by stiching together functional components in BangBang and
# InitialValues.  Thus, FLoops.jl is only responsible for a simple lowering
# which roughly looks like the following:
#
#     acc = InitialValue(⊕)  # aside: actually called in the loop body
#     for x in xs
#         ...  # code before `@reduce`
#
#         input = $rhs  # or `_lazydotcall.($rhs)` if `$rhs` is a dot call
#         acc = broadcast_inplace!!(⊕, acc, input)
#
#         ...  # code after `@reduce`
#     end
#
# where `acc = broadcast_inplace!!(⊕, acc, input)` is a mutate-or-widen version
# of `acc .⊕= input`.  By using `InitialValue(⊕)`, we can uniformly handle the
# reduction by `broadcast_inplace!!`. Furthermore, we can delay the
# initialization/materialization of `acc` until the right moment where size
# information is available through `input`.
#
# In the real code, `acc = InitialValue(⊕)` is called inside the loop body to
# handle the empty iteration consistently (i.e., do not assign variables if
# `@reduce` statement is not evaluated).

function default_initializers(inits, accs, inputs, updaters)
    if inits === nothing
        # e.g., lower `acc += input` to `acc = input`
        initializers = [:($a = $x) for (a, x) in zip(accs, inputs)]
    else
        # e.g., lower `acc = op(init, input)` to `acc = init; acc = op(acc, input)`
        initializers = [:($a = $x) for (a, x) in zip(accs, inits)]
        append!(initializers, updaters)
    end
    return initializers
end

function as_parallel_loop(ctx::MacroContext, rf_arg, coll, body0::Expr, simd, executor)
    body0 = deepcopy(body0)  # To be mutated by `fill_loop_local_variables!`.
    dummy_loop_body = quote  # Dummy expression that simulates the loop body for JuliaVariables.
        $rf_arg = nothing
        $(inject_spec_expr_for_analysis(body0))
    end
    fill_loop_local_variables!(dummy_loop_body)

    accs_symbols = Symbol[]
    inputs_symbols = Symbol[]
    init_exprs = []
    combine_bodies = []
    is_init = Bool[]  # `is_init[i]` is true iff `i`th accumulator is due to `@init`
    all_rf_inits = []
    all_rf_accs = []
    all_rf_inputs = []

    function check_invariance()
        num_state_groups = length(accs_symbols)
        @assert length(inputs_symbols) == num_state_groups
        @assert length(init_exprs) == num_state_groups
        @assert length(is_init) == num_state_groups
        @assert length(all_rf_inits) == num_state_groups
        @assert length(all_rf_accs) == num_state_groups
        @assert length(all_rf_inputs) == num_state_groups
        @assert length(combine_bodies) == num_state_groups

        nums_grouped_states = map(length, all_rf_accs)
        function check_num_states(xs)
            for (x, n) in zip(xs, nums_grouped_states)
                x === nothing && continue
                @assert length(x) == n
            end
        end
        check_num_states(all_rf_inits)
        check_num_states(all_rf_inputs)
    end

    function on_init(spec::InitSpec)
        @gensym grouped_private_states
        push!(accs_symbols, grouped_private_states)
        push!(inputs_symbols, :_)

        accs = spec.lhs
        push!(is_init, true)
        push!(all_rf_inits, nothing)
        push!(all_rf_accs, accs)
        push!(all_rf_inputs, nothing)
        verify_unique_symbols(accs, "private")

        # The corresponding combine function is "keep left" (i.e., do nothing):
        push!(combine_bodies, nothing)

        @gensym init_allocator
        scratch = quote
            $(spec.expr)  # the expression from `@init $initializer`; sets `accs`
            $Base.@inline $init_allocator() = tuple($(spec.rhs...))
            $ScratchSpace($init_allocator, ($(accs...),))
        end

        if isempty(intersect(spec.visible, unbound_rhs(spec.expr)))
            # Hoisting out `@init`, since it is not accessing variables used
            # inside the loop body.
            push!(init_exprs, scratch)
            return quote
                # Just in case it is demoted to `EmptyScratchSpace`:
                $grouped_private_states = $allocate($grouped_private_states)
                ($(accs...),) = $grouped_private_states.value
            end
        else
            push!(init_exprs, _FLoopInit())
            return quote
                if $grouped_private_states isa $_FLoopInit
                    # `$scratch` sets `$accs`.
                    # Assigning to `grouped_private_states` for reusing it next time.
                    $grouped_private_states = $scratch
                else
                    # Just in case it is demoted to `EmptyScratchSpace`:
                    $grouped_private_states = $allocate($grouped_private_states)
                    # After the initialization, just carry it over to the next iteration:
                    ($(accs...),) = $grouped_private_states.value
                end
            end
        end
    end

    body1 = on_reduce_op_spec_reconstructing(body0; on_init = on_init) do spec
        spec = spec::ReduceOpSpec
        @gensym grouped_accs grouped_inputs
        push!(accs_symbols, grouped_accs)
        push!(inputs_symbols, grouped_inputs)
        (inits, accs, inputs, pre_updates, initializers, updaters) =
            process_reduce_op_spec(spec)
        push!(is_init, false)
        push!(all_rf_inits, inits)
        push!(all_rf_accs, accs)
        push!(all_rf_inputs, inputs)
        verify_unique_symbols(accs, "accumulator")
        verify_unique_symbols(inputs, "input")
        # TODO: input symbols just have to be unique within a
        # `@reduce` block.  This restriction (unique across all
        # `@reduce`) can be removed.
        use_oninit =
            inits !== nothing && let visible = spec.visible
                if !isempty(pre_updates)
                    # `pre_updates` defines symbols that may be defined in `inits`
                    visible = mapfoldl(push!, pre_updates; init = copy(visible)) do ex
                        @assert isexpr(ex, :(=), 2)
                        ex.args[1]::Symbol
                    end
                end
                !any(!isempty(intersect(visible, unbound_rhs(ex))) for ex in inits)
            end
        if use_oninit
            # Hoisting out `init` clauses, since it is not accessing variables
            # used inside the loop body.
            push!(init_exprs, :(tuple($(inits...))))
        else
            push!(init_exprs, _FLoopInit())
        end
        function rf_body_with_init(pre_updates = [])
            if use_oninit
                return quote
                    $(pre_updates...)
                    ($(accs...),) = $grouped_accs
                    $(updaters...)
                    $grouped_accs = ($(accs...),)
                end
            end
            quote
                $(pre_updates...)
                if $grouped_accs isa $_FLoopInit
                    # TODO: Check if it is OK to just do `$grouped_accs = $grouped_inputs`
                    # when inside `combine_body`
                    $(initializers...)
                else
                    ($(accs...),) = $grouped_accs
                    $(updaters...)
                end
                $grouped_accs = ($(accs...),)
            end
        end
        combine_body0 = quote
            ($(inputs...),) = $grouped_inputs
            $(rf_body_with_init())
        end
        combine_body = if use_oninit
            combine_body0
        else
            quote
                if $grouped_inputs isa $_FLoopInit
                else
                    $combine_body0
                end
            end
        end
        push!(combine_bodies, combine_body)
        return rf_body_with_init(pre_updates)
    end
    check_invariance()

    body2, info = transform_loop_body(body1, accs_symbols)

    @gensym oninit_function reducing_function combine_function result context_function
    if ctx.module_ === Main
        # Ref: https://github.com/JuliaLang/julia/issues/39895
        oninit_function = Symbol(:__, oninit_function)
        reducing_function = Symbol(:__, reducing_function)
        combine_function = Symbol(:__, combine_function)
    end

    unpackers = map(
        enumerate(zip(is_init, all_rf_accs, all_rf_inits, init_exprs)),
    ) do (i, (nounpack, accs, inits, ex))
        @gensym grouped_accs
        if nounpack
            # This accumulator is from `@init`.
            nothing
        elseif inits === nothing || ex === _FLoopInit()
            quote
                $grouped_accs = $result[$i]
                # Assign to accumulator only if it is updated at least once:
                if $grouped_accs isa $_FLoopInit
                else
                    ($(accs...),) = $grouped_accs
                end
            end
        else
            quote
                $grouped_accs = $result[$i]
                ($(accs...),) = if $grouped_accs isa $_FLoopInit
                    $unreachable_floop()
                    # ($(inits...),)
                else
                    $grouped_accs
                end
            end
        end
    end
    unpack_state = Expr(:block, unpackers...)
    gotos = gotos_for(info.external_labels, unpack_state, result)
    mkdecl(x) = x |> NotA(Nothing) |> Cat() |> Map(a -> :(local $a)) |> collect
    accs_declarations = mkdecl(all_rf_accs)
    inputs_declarations = mkdecl(all_rf_inputs)

    return quote
        $Base.@inline function $oninit_function()
            $(accs_declarations...)
            return tuple($(init_exprs...))
        end
        $Base.@inline function $reducing_function(($(accs_symbols...),), $rf_arg)
            $(accs_declarations...)
            $body2
            return ($(accs_symbols...),)
        end
        $combine_function(_, b::$(Union{Goto,Return})) = b
        function $combine_function(($(accs_symbols...),), ($(inputs_symbols...),))
            $(accs_declarations...)
            $(inputs_declarations...)
            $(combine_bodies...)
            return ($(accs_symbols...),)
        end
        $context_function() = (; ctx = $ctx, id = $(QuoteNode(gensym(:floop_id))))
        $_verify_no_boxes($reducing_function, $context_function)
        $result = $_fold(
            $wheninit(
                $oninit_function,
                $whencombine($combine_function, $reducing_function),
            ),
            $coll,
            $executor,
            $(Val(simd)),
        )
        $result isa $Return && return $result.value
        $(gotos...)
        $unpack_state
        nothing
    end
end

struct _FLoopInit end

@noinline unreachable_floop() = error("unrechable reached (FLoops.jl bug)")

@inline _fold(rf::RF, coll, ::Nothing, simd) where {RF} =
    _fold(rf, coll, PreferParallel(), simd)
@inline _fold(rf::RF, coll, exc::Executor, simd) where {RF} = unreduced(
    transduce(IdentityTransducer(), rf, DefaultInit, coll, maybe_set_simd(exc, simd)),
)

function Base.showerror(io::IO, opspecs::ReduceOpSpec)
    print(io, "`@reduce(")
    join(io, opspecs.args, ", ")
    print(io, ")` used outside `@floop`")
end

function Base.showerror(io::IO, spec::InitSpec)
    ex = spec.expr
    print(io, "`@init", ex, "` used outside `@floop`")
end


text:
function algorithmB(network::AbstractNetwork{T},
                    trips::Matrix{U},
                    costfn::CostFunction;
                    bushtype=SimpleBush{T},
                    λ=0.01,
                    maxiter=1000,
                    errtol=1e-4) where {T<:Integer, U<:Real}

    # unable to state the type for bushtype as `::Type{<:AbstractBush{T}}`

    zoneids = id.(zones(network))
    nzones = length(zoneids)
    nlinks = numlinks(network)

    # initialize
    bushes = [bushtype(network, r, costfn) for r in zoneids]
    flows = allornothing(network, trips, costfn; basedon=:origin)

    # iteration
    linkcosts, linkcostdrvs = costfn(flows; returnitems=[:costs, :derivs]) # updated as and when flows change
    LARGE_NUM = typemax(eltype(linkcosts))/2 # large constant used below
    nadded, nremoved = 1, 0 # number of links swithing in and out of bushes
    iterno = 0
    outererr = 1.
    while (outererr > errtol) && (iterno <= maxiter) #(nadded + nremoved > 0)
        # equilibrate bushes
        totalshift = zeros(Bool, nlinks, nzones)
        err = 1. # problem solved for the bushes completely
        while err > errtol
            # equilibrate core step
            for bush in bushes
                r = src(bush)
                scosts, sparentvec = acyclic(bush, linkcosts; kind=:short)

                for j in reverse(orderednodes(bush))
                    #lcosts, lparentvec = acyclic(bush, linkcosts - LARGE_NUM * totalshift[:,r]; kind=:long)
                    lcosts, lparentvec = acyclic(bush, linkcosts - LARGE_NUM * iszero.(flows[:,r]); kind=:long)
                    if (j == r) || (lcosts[j] ≈ scosts[j]) || (trips[r,j] == 0.)
                        continue
                    end

                    spath = pathto(sparentvec, j)
                    lpath = pathto(lparentvec, j)

                    max_i = min(length(spath), length(lpath))

                    # find path segments
                    i = 1
                    while (i <= max_i) && (spath[i] == lpath[i])
                        i += 1
                    end
                    sa = la = i

                    i = 0
                    while (i < max_i) && (spath[end-i] == lpath[end-i])
                        i += 1
                    end
                    sb = length(spath)-i
                    lb = length(lpath)-i

                    # link indices
                    ssegids = idx.(Ref(bush), spath[sa:sb])
                    lsegids = idx.(Ref(bush), lpath[la:lb])

                    scost, sdcost = sum(linkcosts[ssegids]), sum(linkcostdrvs[ssegids])
                    lcost, ldcost = sum(linkcosts[lsegids]), sum(linkcostdrvs[lsegids])

                    lflow = minimum(flows[lsegids,r])
                    dx = min(lflow, λ * (lcost - scost)/(ldcost + sdcost))
                    if iszero(lflow)
                        return (r, lsegids, bush, flows, totalshift)
                    end

                    flows[ssegids,r] .+= dx
                    flows[lsegids,r] .-= dx

                    # check any zero flows on longpath
                    for lsi in lsegids
                        if iszero(flows[lsi,r])
                            totalshift[lsi,r] = true
                        end
                    end
                end
            end
            totalshift .&= iszero.(flows)
            tmpd = costfn(flows; returnitems=[:costs, :derivs]) # careful with the order
            linkcosts, linkcostdrvs = tmpd[:costs], tmpd[:derivs]
            numer = 0.
            denom = 0.
            for (i,bush) in enumerate(bushes)
                r = src(bush)
                sflows = allornothing(bush, r, trips[:,r], linkcosts)
                numer += sum(flows[:,i] .* linkcosts)
                denom += sum(sflows .* linkcosts)
            end
            err = numer/denom - 1.
            #println(iterno, "\t", outererr, "\t", err)
        end

        # update bushes
        nadded, nremoved = updatebushes!(bushes, flows, linkcosts)

        outersflows = allornothing(network, trips, linkcosts; basedon=:link)
        outererr = sum(flows .* linkcosts)/sum(outersflows .* linkcosts) - 1.
        iterno += 1
        #println(nadded + nremoved, "\t", iterno, "\t", outererr)
    end
    return (flows, outererr)
end

function updatebushes!(bushes::Vector{<:AbstractBush}, flows::Matrix, linkcosts::Vector)
    nadded, nremoved = 0, 0
    for (iteri,bush) in enumerate(bushes)
        r = src(bush)
        #scosts, sparentvec = acyclic(bush, linkcosts; kind=:short)
        lcosts, lparentvec = acyclic(bush, linkcosts; kind=:long)
        for link in links(net(bush))
            i,j,linkidx = upn(link),dwn(link),idx(net(bush), link)

            if has_link(bush, i, j)
                if iszero(flows[linkidx,iteri]) && !(indegree(bush, j) == 1)
                    rem_link!(bush, i, j)
                    nremoved += 1
                end
            else
                if (lcosts[j] > lcosts[i] + linkcosts[linkidx])
                    add_link!(bush, i, j)
                    nadded += 1
                end
            end
        end
    end
    return (nadded, nremoved)
end

ext: jl


text:
function encode(input::AbstractString)

end

function decode(input::AbstractString)

end


text:
using Documenter, CuArrays

using Pkg
if haskey(ENV, "GITLAB_CI")
    Pkg.add([PackageSpec(name = x; rev = "master")
             for x in ["CUDAapi", "GPUArrays", "CUDAnative", "NNlib", "CUDAdrv"]])
end

makedocs(
    modules = [CuArrays],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "CuArrays.jl",
    pages = [
        "Home" => "index.md",
    ],
    doctest = true
)


text:
"""
	mpvca : example from Bazinga.jl

    Academic example of mathematical program with vanishing constraints (MPVC),
    see Section 9.5.1 [Hoh09].

	Original formulations:

	minimize      4*x1 + 2*x2
	subject to    x1                       >= 0
	              x2                       >= 0
	              (x1 + x2 - 5*sqrt(2))*x1 >= 0
	              (x1 + x2 - 5        )*x2 >= 0
 		          x1 + x2 - 3              >= 0 (ncon == 5)

	Reformulation as a constrained structured problem in the form

	minimize     f(x) + g(x)
	subject to   c(x) in S

	where
	    f(x) = 4*x1 + 2*x2
	    g(x) = IndBox[0,+infty](x)
	    c(x) = [ x1                  ]
	           [ x2                  ]
	           [ x1 + x2 - 5*sqrt(2) ]
	           [ x1 + x2 - 5         ]
		       [ x1 + x2 - 3         ] (ncon == 5)
	    S   = { c | (c1,c3) in SV, (c2,c4) in SV, c5 >= 0 }
	with
	    SV  = { (a,b) | a >= 0, a*b >= 0 }

    References:
    [Hoh09]     Hoheisel, "Mathematical Programs with Vanishing Constraints",
                PhD thesis, University of Würzburg, 2009.
"""

using Bazinga, OptiMo
using Random, LinearAlgebra
using DataFrames, CSV
using Printf, Plots

###################################################################################
# problem definition
###################################################################################
mutable struct MPVCA <: AbstractOptiModel
    meta::OptiModelMeta
end

function MPVCA(; ncon::Int = 4)
    @assert 4 <= ncon <= 5
    name = "mpvca" * "_$ncon"
    meta = OptiModelMeta(2, ncon, x0 = [5.0; 5.5], name = name)
    return MPVCA(meta)
end

# necessary methods:
# obj, grad!: cons!, jprod!, jtprod!, proj!, prox!, objprox!
function OptiMo.obj(prob::MPVCA, x::AbstractVector)
    OptiMo.@lencheck prob.meta.nvar x
    return 4 * x[1] + 2 * x[2]
end

function OptiMo.grad!(prob::MPVCA, x::AbstractVector, dfx::AbstractVector)
    OptiMo.@lencheck prob.meta.nvar x
    dfx .= [4; 2]
    return nothing
end

function OptiMo.cons!(prob::MPVCA, x::AbstractVector, cx::AbstractVector)
    OptiMo.@lencheck prob.meta.nvar x
    OptiMo.@lencheck prob.meta.ncon cx
    cx[1:4] .= [x[1]; x[2]; x[1] + x[2] - 5.0 * sqrt(2.0); x[1] + x[2] - 5.0]
    if prob.meta.ncon > 4
        cx[5] = x[1] + x[2] - 3.0
    end
    return nothing
end

function OptiMo.jprod!(
    prob::MPVCA,
    x::AbstractVector,
    v::AbstractVector,
    Jv::AbstractVector,
)
    OptiMo.@lencheck prob.meta.nvar x v
    OptiMo.@lencheck prob.meta.ncon Jv
    v12 = v[1] + v[2]
    if prob.meta.ncon > 4
        Jv .= [v[1]; v[2]; v12; v12; v12]
    else
        Jv .= [v[1]; v[2]; v12; v12]
    end
    return nothing
end

function OptiMo.jtprod!(
    prob::MPVCA,
    x::AbstractVector,
    v::AbstractVector,
    Jtv::AbstractVector,
)
    OptiMo.@lencheck prob.meta.nvar x Jtv
    OptiMo.@lencheck prob.meta.ncon v
    if prob.meta.ncon > 4
        v345 = v[3] + v[4] + v[5]
    else
        v345 = v[3] + v[4]
    end
    Jtv .= [v[1] + v345; v[2] + v345]
    return nothing
end

function OptiMo.prox!(prob::MPVCA, x::AbstractVector, a::Real, z::AbstractVector)
    OptiMo.@lencheck prob.meta.nvar x z
    z .= max.(0.0, x)
    return nothing
end

function OptiMo.objprox!(prob::MPVCA, x::AbstractVector, a::Real, z::AbstractVector)
    OptiMo.@lencheck prob.meta.nvar x z
    OptiMo.prox!(prob, x, a, z)
    return 0.0
end

function OptiMo.proj!(prob::MPVCA, cx::AbstractVector, px::AbstractVector)
    OptiMo.@lencheck prob.meta.ncon cx px
    # vanishing constraint
    px[1] = max(0.0, cx[1])
    px[3] = cx[3]
    if (px[1] + px[3] < 0)
        px[1] = 0.0
    else
        px[3] = max(px[3], 0.0)
    end
    # vanishing constraint
    px[2] = max(0.0, cx[2])
    px[4] = cx[4]
    if (px[2] + px[4] < 0)
        px[2] = 0.0
    else
        px[4] = max(px[4], 0.0)
    end
    # inequality (ncon == 5)
    if prob.meta.ncon > 4
        px[5] = max(0.0, cx[5])
    end
    return nothing
end

###############################################################################
###############################################################################
###############################################################################
#=function MPVCA()
    ncon = 4
    name = "mpvca" * "_5g"
    meta = OptiModelMeta(2, ncon, x0 = [5.0; 5.0], name = name)
    return MPVCA(meta)
end
function OptiMo.prox!(prob::MPVCA, x::AbstractVector, a::Real, z::AbstractVector)
    OptiMo.@lencheck prob.meta.nvar x z
    c = 3.0
    @assert(c ≥ 0.0)
    if x[1] ≥ c || x[2] ≥ c
        z .= max.(0.0, x)
    elseif x[1] + x[2] ≥ c
        z .= max.(0.0, x)
    elseif x[2] - x[1] ≥ c
        z .= [0.0; c]
    elseif x[1] - x[2] ≥ c
        z .= [c; 0.0]
    else
        z[1] = min(c, (c+x[1]-x[2])/2.0)
        z[2] = min(c, (c-x[1]+x[2])/2.0)
    end
    return nothing
end
=#

###############################################################################
###############################################################################
###############################################################################
foldername = "/home/alberto/Documents/Bazinga.jl/demo/data/"

# problem build
#problem = MPVCA(ncon = 4)
problem = MPVCA(ncon = 5)
#problem = MPVCA()

# solver build
solver =
    Bazinga.ALPX(max_iter = 10, max_sub_iter = 1000, verbose = false, subsolver = :zerofpr)

R = eltype(problem.meta.x0)
nvar = problem.meta.nvar

xmin = -5.0
xmax = 20.0

data = DataFrame()

xgrid = [(i, j) for i = xmin:0.5:xmax, j = xmin:0.5:xmax];
xgrid = xgrid[:];
ntests = length(xgrid)

for i = 1:ntests
    x0 = [xgrid[i][1]; xgrid[i][2]]

    out = solver(problem, x0 = x0)

    @printf "."
    if mod(i, 50) == 0
        @printf "\n"
    end

    push!(data, (id = i, xi_1 = x0[1], xi_2 = x0[2], xf_1 = out.x[1], xf_2 = out.x[2], iter=out.iterations, time=out.time))

end
@printf "\n"

filename = problem.meta.name * "_grid"
CSV.write(foldername * filename * ".csv", data, header=false)

################################################################################
tolx = 1e-3

global c_00 = 0
global c_05 = 0
global c_11 = 0
global c_un = 0

pyplot()

feasset = Shape([(0.0,xmax),(0.0,5*sqrt(2)),(5*sqrt(2),0.0),(xmax,0.0),(xmax,xmax),(0.0,xmax)])
hplt = plot(feasset, color=plot_color(:grey,0.4), linewidth=0, legend = false)
xlims!(xmin,xmax)
ylims!(xmin,xmax)

for i in 1:ntests
    xi = [data[i,2]; data[i,3]]
    xf = [data[i,4]; data[i,5]]

    if norm(xf,Inf) <= tolx
        global c_00 += 1
        scatter!(hplt, [xi[1]], [xi[2]], color=:blue, marker=:circle, markerstrokewidth=0, legend = false)
    elseif norm(xf - [0.0;5.0],Inf) <= tolx
        global c_05 += 1
        scatter!(hplt, [xi[1]], [xi[2]], color=:red, marker=:diamond, markerstrokewidth=0, legend = false)
    elseif xf[1] >= 0 && xf[2] >= 0 && xf[1]+xf[2] <= 3
        global c_11 += 1
        scatter!(hplt, [xi[1]], [xi[2]], color=:green, marker=:star5, markerstrokewidth=0, legend = false)
    else
        global c_un += 1
    end
end

savefig(foldername * filename * ".pdf")


text:
include("SumSpace.jl")
include("ArraySpace.jl")
include("ProductSpaceOperators.jl")
include("BlockOperators.jl")
include("SubSpace.jl")

⊕(A::Space,B::Space)=domainscompatible(A,B)?SumSpace(A,B):PiecewiseSpace(A,B)
⊕(f::Fun,g::Fun)=Fun(interlace(coefficients(f),coefficients(g)),space(f)⊕space(g))

# Conversion from Vector to Tuple
# If a Vector fun has different spaces in each component we represenent
#  it by a Tuple fun, so this allows conversion.

function coefficients(f::Vector,a::VectorSpace,b::TupleSpace)
    A=a.space
    n=length(a)
    @assert n==length(b.spaces)
    ret=copy(f)
    for k=1:n
        ret[k:n:end]=coefficients(ret[k:n:end],A,b.spaces[k])
    end
    ret
end

#split the cfs into component spaces
function coefficients(cfs::Vector,A::SumSpace,B::SumSpace)
    if spacescompatible(A,B)
        cfs
    else
        mapreduce(f->Fun(f,B),+,vec(Fun(cfs,A))).coefficients
    end
end
function coefficients(cfs::Vector,A::PiecewiseSpace,B::PiecewiseSpace)
    if spacescompatible(A,B)
        cfs
    else
        mapreduce(f->Fun(f,B),+,pieces(Fun(cfs,A))).coefficients
    end
end

# add constant to each piece
function coefficients(cfsin::Vector,A::ConstantSpace,B::PiecewiseSpace)
    m=length(B.spaces)

    #assume first that  each space needs 1 coefficient to rep constant
    ret=zeros(eltype(cfsin),m)

    for k=1:m
        cfs=coefficients(cfsin,A,B.spaces[k])
        mm=m*(length(cfs)-1)+k  # max ret, should normally be k
        ret=pad!(ret,max(length(ret),mm))
        ret[k:m:mm]=cfs
    end

    ret
end

for TYP in (:SumSpace,:PiecewiseSpace)
    @eval begin
        # we need to be able to call this to avoid confusion
        function sumspacecoefficients(cfsin::Vector,A::Space,B::$TYP)
            m=length(B.spaces)

            for k=1:m
                if isconvertible(A,B[k])
                    cfs=coefficients(cfsin,A,B[k])
                    ret=zeros(eltype(cfs),k-1)  # all spaces have dimension at least 1,
                                                # so there are k-1 zero coefficients
                                                # corresponding to all spaces less than
                                                # k

                    j=k    # represents the row of ret to be added
                    row=1  # represents the row of cfs to be added

                    @assert length(cfs) ≤ dimension(B[k])
                    dk=length(cfs)

                    # we only allow at most dimension
                    while row ≤ dk
                        push!(ret,cfs[row])  # this sets ret[j] to cfs[row]
                        for λ = k+1:m
                            if dimension(B[λ]) ≥ row # loop through the rest of the spaces
                                push!(ret,0)
                                j+=1
                            end
                        end
                        row += 1   # move on to the next row
                        for λ = 1:k-1
                            if dimension(B[λ]) ≥ row # loop through the previous spaces
                                push!(ret,0)
                                j+=1
                            end
                        end
                        j+=1  # we always increment by 1 for the current space
                    end

                    return ret
                end
            end

            defaultcoefficients(cfsin,A,B)
        end
        coefficients(cfsin::Vector,A::Space,B::$TYP)=sumspacecoefficients(cfsin,A,B)
    end
end

## LowRank Constructors

# convert a vector of functionals and an operator to a LowRnakPertOperator
# the rangespace is a DirectSumSpace specified by ST of the input rangespaces
# the default is a  TupleSpace, but support is there for PiecewiseSpace
# for bcs
for TYP in (:PiecewiseSpace,:TupleSpace)
    @eval function LowRankPertOperator{OT<:Operator}(A::Vector{OT},::Type{$TYP})
        A=promotedomainspace(A)
        for k=1:length(A)-1
            @assert isafunctional(A[k])
        end
        @assert isbanded(A[end])
        L=LowRankOperator(A[1:end-1],$TYP)
        # add zero functionals to shift down
        BB=[fill(ZeroOperator(domainspace(BB),ConstantSpace()),length(A)-1);A[end]]
        S=InterlaceOperator(BB,domainspace(BB),$TYP(map(rangespace,A)))
        L+S
    end
end

LowRankPertOperator{OT<:Operator}(A::Vector{OT})=LowRankPertOperator(A,TupleSpace)

function LowRankOperator{FT<:Operator}(Bin::Vector{FT},::Type{PiecewiseSpace})
    B=promotedomainspace(Bin)
    rsp=PiecewiseSpace(map(rangespace,B))
    LowRankOperator(
        Fun{typeof(rsp),Float64}[Fun([zeros(k-1);1],rsp) for k=1:length(B)],
        B)
end

function LowRankOperator{FT<:Operator}(Bin::Vector{FT},::Type{TupleSpace})
    B=promotedomainspace(Bin)
    rsp=TupleSpace(tuple(map(rangespace,B)...,ZeroSpace()))  #TODO: Why the hack?
    LowRankOperator(
        Fun{typeof(rsp),Float64}[Fun([zeros(k-1);1],rsp) for k=1:length(B)],
        B)
end

LowRankOperator{FT<:Operator}(Bin::Vector{FT})=LowRankOperator(Bin,TupleSpace)


text:
module Rainflow
import Base.show

export sort_peaks, find_boundary_vals, count_cycles, sum_cycles

""" This function sort out points where the slope is changing sign"""
function sort_peaks(signal::AbstractArray{Float64,1}, dt=collect(1.:length(signal)))
    slope = diff(signal)
    # Determines if the point is local extremum
    is_extremum = vcat(true, (slope[1:end-1].*slope[2:end]).<=0., true)
    return signal[is_extremum] , dt[is_extremum]
end

struct Cycle  # This is the information stored for each cycle found
    count::Float64
    range::Float64
    mean::Float64
    Rvalue::Float64   # value
    v_s::Float64   # value start
    t_s::Float64   # time start
    v_e::Float64   # value end
    t_e::Float64   # time end
end

show(io::IO,x::Cycle) = print(io, "Cycle: count=", x.count, ", range=",x.range, ", mean=",x.mean, ", R=", x.Rvalue)

function cycle(count::Float64, v_s::Float64, t_s::Float64, v_e::Float64, t_e::Float64)
    Cycle(count, abs(v_s-v_e), (v_s+v_e)/2, min(v_s,v_e)/max(v_s,v_e), v_s, t_s, v_e, t_e)
end

""" Count the cycles from the data """
function count_cycles(peaks::Array{Float64,1},t::Array{Float64,1})
    list = copy(peaks) # Makes a copy because there is going to be sorted in the vectors
    time = copy(t)
    currentindex = 1
    nextindex = 2
    cycles = Cycle[]
    sizehint!(cycles, length(list)) # This reduces the memory consumption a bit
    @inbounds begin
    while length(list) > (currentindex+1)
            currentvalue = abs(list[currentindex+1]-list[currentindex])
            nextvalue = abs(list[nextindex+1]-list[nextindex])
        if nextvalue > currentvalue
            if currentindex == 1 # This case counts a half and cycle deletes the poit that is counted
                push!(cycles,cycle(0.5 ,list[currentindex], time[currentindex], list[currentindex+1],time[currentindex+1]))
                popfirst!(list) # Removes the first entrance in ext and time
                popfirst!(time)
            else # This case counts one cycle and deletes the point that is counted
                push!(cycles,cycle(1. ,list[currentindex], time[currentindex], list[currentindex+1], time[currentindex+1]))
                deleteat!(list, currentindex:(currentindex+1)) # Removes the i and i+1 entrance in ext and time
                deleteat!(time, currentindex:(currentindex+1))
            end
            currentindex = 1
            nextindex = 2
        else
            currentindex += 1
            nextindex += 1
        end
    end
    for currentindex=1:length(list)-1 # This counts the rest of the points that have not been counted as a half cycle

        push!(cycles,cycle(0.5 ,list[currentindex], time[currentindex], list[currentindex+1],time[currentindex+1]))
    end
    end
    return cycles
end

mutable struct Cycles_bounds #
    min_mean::Float64
    max_mean::Float64
    max_range::Float64
    min_R::Float64
    max_R::Float64
end

show(io::IO,x::Cycles_bounds) = print(io, "Cycles_bounds : min mean value=", x.min_mean, ", max mean value=", x.max_mean, ", max range=",x.max_range, ", min R=", x.min_R, ", max R=",x.max_R)

""" Find the minimum and maximum mean value and maximum range from a vector of cycles"""
function find_boundary_vals(cycles::Array{Cycle,1})
    bounds = Cycles_bounds(Inf, -Inf, -Inf, Inf, -Inf)
    for cycle in cycles
        cycle.mean > bounds.max_mean && setfield!(bounds, :max_mean, cycle.mean)
        cycle.mean < bounds.min_mean && setfield!(bounds, :min_mean, cycle.mean)
        cycle.range > bounds.max_range && setfield!(bounds, :max_range, cycle.range)
        cycle.Rvalue < bounds.min_R && setfield!(bounds, :min_R, cycle.Rvalue)
        cycle.Rvalue > bounds.max_R && setfield!(bounds, :max_R, cycle.Rvalue)
    end
    return bounds
end

""" Returns the range index the value is belonging in """
function find_range(interval::Array{T,1},value) where {T <: Real}
    for i=1:length(interval)-1
        if interval[i] <= value < interval[i+1]
            return i
        end
    end
    error("The value where not in range")
end

Interval{T} = Union{Array{T,1}, StepRangeLen{T}}

""" Sums the cycle count given intervals of range_intervals and mean_intervals. The range_intervals and mean_intervals is given in fraction of range size"""
function sum_cycles(cycles::Array{Cycle,1}, range_intervals::Interval{T}, mean_intervals::Interval{T}) where {T <: Real}
    bounds = find_boundary_vals(cycles)
    bins = zeros(length(range_intervals)-1, length(mean_intervals)-1)
    range_in = (range_intervals*bounds.max_range)/100
    mean_in = (mean_intervals*(bounds.max_mean-bounds.min_mean))/100
    mean_in = mean_in .+ bounds.min_mean
    issorted(mean_intervals) || error("The array needs to be sorted in raising order")
    issorted(range_intervals) || error("The array needs to be sorted in raising order")
    nr_digits = 14  # The rounding is performed due to numerical noise in the floats when comparing
    mean_i = collect(mean_in)
    range_i = collect(range_in)
    #ensure the cycles is in the intervals by adding a small error to the end values of the interal.
    error_m = (bounds.max_mean-bounds.min_mean)*1e-14
    mean_i[end]+=error_m
    mean_i[1]-=error_m
    error_r = bounds.max_range*1e-14
    range_i[end]+=error_r
    range_i[1]-=error_r
    #show(mean_intervals)
    for cycle in cycles
        i = find_range(range_i, cycle.range)
        j = find_range(mean_i, cycle.mean)
        bins[i,j] += cycle.count
    end
    return bins
end

function sum_cycles(cycles::Array{Cycle,1}, nr_ranges::Int=10, nr_means::Int=1)
    range_intervals = range(0,stop=100,length=nr_ranges+1)
    mean_intervals = range(0,stop=100,length=nr_means+1)
    sum_cycles(cycles, range_intervals, mean_intervals)
end

end # module


text:
"""
    function addParameters_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

Generate the variables of the model and declare them.
"""
function addParameters_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

    #Parameter declare###########################################################################

    if length(abm.declaredSymbols["Local"])>0
        append!(p.declareVar.args,
            (quote
                localV = Array(Float64.([com.local_;Base.zeros(nMax-N,$(length(abm.declaredSymbols["Local"])))]))
            end).args
        )

        push!(p.args,:localV)

        if !isempty(p.update["Local"])
            append!(p.declareVar.args,
                (quote
                    localVCopy = zeros(Float64,size(localV)[1],$(length(keys(p.update["Local"]))))
                end).args
            )

            push!(p.args,:localVCopy)
        end
    end

    if length(abm.declaredSymbols["Identity"])>0
        append!(p.declareVar.args,
            (quote
                identityV = Array([Int.(com.identity_);Base.zeros(Int,nMax-N,$(length(abm.declaredSymbols["Identity"])))])
            end).args
        )

        push!(p.args,:identityV)

        if !isempty(p.update["Identity"])
            append!(p.declareVar.args,
                (quote
                    identityVCopy = zeros(Int,size(identityV)[1],$(length(keys(p.update["Identity"]))))
                end).args
            )

            push!(p.args,:identityVCopy)
        end
    end

    if length(abm.declaredSymbols["Global"])>0
        append!(p.declareVar.args,
            (quote
                globalV = Array(Float64.(com.global_))
            end).args
        )

        push!(p.args,:globalV)
        if !isempty(p.update["Global"]) && !isempty([ v for v in keys(p.update["Global"]) if v in abm.declaredSymbols["Global"] ])
            append!(p.declareVar.args,
                (quote
                    globalVCopy = zeros(Float64,$(length(p.update["Global"])))
                end).args
            )

            push!(p.args,:globalVCopy)
        end
    end

    if length(abm.declaredSymbols["GlobalArray"])>0
        for (j,i) in enumerate(abm.declaredSymbols["GlobalArray"])
            append!(p.declareVar.args,
            (quote
                $(Meta.parse(string(i))) = Array(copy(com.globalArray_[$j]))
            end).args
            )

            push!(p.args,Meta.parse(string(i)))
        end
        if !isempty(p.update["GlobalArray"])
            for (j,i) in enumerate(keys(p.update["GlobalArray"]))
                append!(p.declareVar.args,
                (quote
                    $(Meta.parse(string(i,GLOBALARRAYCOPY))) = Array(copy(com.globalArray_[$j]))
                end).args
                )

                push!(p.args,Meta.parse(string(i,GLOBALARRAYCOPY)))
            end
        end
    end

    if length(abm.declaredSymbols["Medium"])>0
        if abm.dims >= 1
            if space.medium[1].minBoundaryType == "Periodic"
                push!(p.declareVar.args, :(Nx_ = size(com.medium_)[1]+2))
                push!(p.declareVar.args, :(mediumV = cat(zeros(1,size(com.medium_)[2:end]...),com.medium_,zeros(1,size(com.medium_)[2:end]...),dims=1)))
            else
                push!(p.declareVar.args, :(Nx_ = size(com.medium_)[1]))
                push!(p.declareVar.args, :(mediumV = copy(com.medium_)))
            end

            push!(p.declareVar.args, :(dxₘ_ = ($(space.box[1].max)-$(space.box[1].min))/$(space.medium[1].N)))

            push!(p.args,:Nx_)
            push!(p.args,:dxₘ_)
        end
        if abm.dims >= 2
            if space.medium[1].minBoundaryType == "Periodic"
                push!(p.declareVar.args, :(Ny_ = size(com.medium_)[2]+2))
                push!(p.declareVar.args, :(mediumV = cat(zeros(Nx_,1,size(com.medium_)[3:end]...),mediumV,zeros(Nx_,1,size(com.medium_)[3:end]...),dims=2)))
            else
                push!(p.declareVar.args, :(Ny_ = size(com.medium_)[2]))
            end

            push!(p.declareVar.args, :(dyₘ_ = ($(space.box[2].max)-$(space.box[2].min))/$(space.medium[2].N)))

            push!(p.args,:Ny_)
            push!(p.args,:dyₘ_)
        end
        if abm.dims >= 3
            if space.medium[1].minBoundaryType == "Periodic"
                push!(p.declareVar.args, :(Nz_ = size(com.medium_)[3]+2))
                push!(p.declareVar.args, :(mediumV = cat(zeros(Nx_,Ny_,1,size(com.medium_)[4:end]...),mediumV,zeros(Nx_,Ny_,1,size(com.medium_)[4:end]...),dims=3)))
            else
                push!(p.declareVar.args, :(Nz_ = size(com.medium_)[3]))
            end

            push!(p.declareVar.args, :(dzₘ_ = ($(space.box[3].max)-$(space.box[3].min))/$(space.medium[3].N)))

            push!(p.args,:Nz_)
            push!(p.args,:dzₘ_)
        end
        push!(p.declareVar.args, :(mediumV = Array(mediumV)))

        push!(p.args,:mediumV)

        if "UpdateMedium" in keys(abm.declaredUpdates)
            push!(p.declareVar.args, :(mediumVCopy = copy(mediumV)))
            push!(p.args,:mediumVCopy)
        end

        #Add identifiers of medium position
        push!(p.declareVar.args,
            :(idMediumV = zeros(Int,N,$(p.agent.dims)))
        )
        push!(p.args,:idMediumV)

        if p.agent.dims == 1
            code = simpleFirstLoopWrapInFunction_(platform,:idMediumCompute!,
                quote
                    idMediumV[ic1_,1] = floor(Int,(localV[ic1_,1]-$(p.space.box[1].min))/($(p.space.box[1].max)-$(p.space.box[1].min))*Nx_)+1
                end
            )
        elseif p.agent.dims == 2
            code = simpleFirstLoopWrapInFunction_(platform,:idMediumCompute!,
                quote
                    idMediumV[ic1_,1] = floor(Int,(localV[ic1_,1]-$(p.space.box[1].min))/($(p.space.box[1].max)-$(p.space.box[1].min))*Nx_)+1
                    idMediumV[ic1_,2] = floor(Int,(localV[ic1_,2]-$(p.space.box[2].min))/($(p.space.box[2].max)-$(p.space.box[2].min))*Ny_)+1
                end
            )
        elseif p.agent.dims == 3
            code = simpleFirstLoopWrapInFunction_(platform,:idMediumCompute!,
                quote
                    idMediumV[ic1_,1] = floor(Int,(localV[ic1_,1]-$(p.space.box[1].min))/($(p.space.box[1].max)-$(p.space.box[1].min))*Nx_)+1
                    idMediumV[ic1_,2] = floor(Int,(localV[ic1_,2]-$(p.space.box[2].min))/($(p.space.box[2].max)-$(p.space.box[2].min))*Ny_)+1
                    idMediumV[ic1_,3] = floor(Int,(localV[ic1_,3]-$(p.space.box[3].min))/($(p.space.box[3].max)-$(p.space.box[3].min))*Nz_)+1
                end
            )
        end
        push!(p.declareF.args, code)
        push!(p.execInit.args, :(@platformAdapt idMediumCompute!(ARGS)))
        push!(p.execInloop.args, :(@platformAdapt idMediumCompute!(ARGS)))

    end

    return nothing
end

"""
    function addUpdate_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

Generate the functions to update all the modified values.
"""
function addUpdate_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

    create = false

    for i in keys(abm.declaredUpdates)
        if !(i in ["UpdateLocalInteraction","UpdateInteraction","EventDivision"]) && !emptyquote_(abm.declaredUpdates[i])
            create = true
        end
    end

    if create
        gen = quote end #General update
        up = quote end
        for i in keys(abm.declaredSymbols)
            if  i == "GlobalArray"
                for j in keys(p.update[i])
                    n = Meta.parse(string(j,GLOBALARRAYCOPY))
                    push!(gen.args,:($j.=$n))
                end
            elseif i in ["Local","Identity"]
                var = Meta.parse(string(lowercase(i),"V"))
                varCopy = Meta.parse(string(lowercase(i),"VCopy"))
                for j in keys(p.update[i])
                    pos = findfirst(abm.declaredSymbols[i] .== j)
                    posCopy = p.update[i][j]
                    push!(up.args,:($var[ic1_,$pos]=$varCopy[ic1_,$posCopy]))
                end
            elseif i == "Global"
                var = Meta.parse(string(lowercase(i),"V"))
                varCopy = Meta.parse(string(lowercase(i),"VCopy"))

                upb = quote end #Assign updates of the global to the first thread
                for j in keys(p.update[i])
                    pos = findfirst(abm.declaredSymbols[i] .== j)
                    posCopy = p.update[i][j]
                    push!(upb.args,:($var[$pos]=$varCopy[$posCopy]))
                end

                push!(up.args,
                    :(begin
                        if ic1_ == 1
                            $(upb)
                        end
                    end))
            elseif i == "Medium"
                if "UpdateMedium" in keys(abm.declaredUpdates)
                    push!(gen.args,:(mediumV.=mediumVCopy))
                end
            end

        end

        #Construct global and local update function
        f = simpleFirstLoopWrapInFunction_(platform,:updateLocGlob_!,up)
        push!(p.declareF.args,f)
        #Construct general update function
        push!(p.declareF.args,
            :(begin
                function update_!(ARGS_)
                    $gen
                    @platformAdapt updateLocGlob_!(ARGS_)
                    return
                end
            end)
            )

        #Add to execution cue
        push!(p.execInloop.args,:(update_!(ARGS_)))
    end

    return nothing
end

"""
    function addCopyInitialisation_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

Generate the functions to update all the modified values.
"""
function addCopyInitialisation_!(p::Program_,abm::Agent,space::SimulationSpace,platform::String)

    create = false

    for i in keys(abm.declaredUpdates)
        if !(i in ["UpdateLocalInteraction","UpdateInteraction","EventDivision"]) && !emptyquote_(abm.declaredUpdates[i])
            create = true
        end
    end

    if create
        gen = quote end #General update
        up = quote end
        for i in keys(abm.declaredSymbols)
            if  i == "GlobalArray"
                for j in keys(p.update[i])
                    n = Meta.parse(string(j,GLOBALARRAYCOPY))
                    push!(gen.args,:($n.=$j))
                end
            elseif i in ["Local","Identity"]
                var = Meta.parse(string(lowercase(i),"V"))
                varCopy = Meta.parse(string(lowercase(i),"VCopy"))
                for j in keys(p.update[i])
                    pos = findfirst(abm.declaredSymbols[i] .== j)
                    posCopy = p.update[i][j]
                    push!(up.args,:($varCopy[ic1_,$posCopy]=$var[ic1_,$pos]))
                end
            elseif i == "Global"
                var = Meta.parse(string(lowercase(i),"V"))
                varCopy = Meta.parse(string(lowercase(i),"VCopy"))

                upb = quote end #Assign updates of the global to the first thread
                for j in keys(p.update[i])
                    pos = findfirst(abm.declaredSymbols[i] .== j)
                    posCopy = p.update[i][j]
                    push!(upb.args,:($varCopy[$posCopy]=$var[$pos]))
                end

                push!(up.args,
                    :(begin
                        if ic1_ == 1
                            $(upb)
                        end
                    end))
            end

        end

        #Construct global and local update function
        f = simpleFirstLoopWrapInFunction_(platform,:initialiseLocGlobCopy_!,up)
        push!(p.declareF.args,f)
        #Construct general update function
        push!(p.declareF.args,
            :(begin
                function initialiseCopy_!(ARGS_)
                    $gen
                    @platformAdapt initialiseLocGlobCopy_!(ARGS_)
                    return
                end
            end)
            )

        #Add to execution cue
        push!(p.execInit.args,:(initialiseCopy_!(ARGS_)))
    end

    return nothing
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
module OpenQASM

using RBNF

include("types.jl")
include("parse.jl")
include("tools.jl")

using .Types: print_qasm
using .Tools: cmp_ast

"""
    parse(qasm::String)

Parse a piece of QASM program at top-level to AST.
"""
function parse(src::String)
    ast, ctx = RBNF.runparser(Parse.mainprogram, RBNF.runlexer(Parse.QASMLang, src))
    ctx.tokens.current > ctx.tokens.length || throw(Meta.ParseError("invalid syntax in QASM program"))
    return ast
end

"""
    parse(qasm::String)

Parse a piece of QASM gate program.
"""
function parse_gate(src::String)
    ast, ctx = RBNF.runparser(Parse.gate, RBNF.runlexer(Parse.QASMLang, src))
    ctx.tokens.current > ctx.tokens.length || throw(Meta.ParseError("invalid syntax in QASM program"))
    return ast
end

end


text:
using Fn
using MacroTools

using Test

@testset "Fn.jl" begin
    @testset "no _ creates 0-argument function" begin
        @test @capture(Fn.fn(:(exp(5 * √π))), function () exp(5 * √π) end)
    end

    @testset "_ creates single argument function" begin
        @test @capture(Fn.fn(:(5 + _)), function (x_) 5 + x_ end)
    end

    @testset "numbered underscores introduce multiple arguments" begin
        @test @capture(Fn.fn(:(sin(_1) * cos(_2))), function (x_, y_) sin(x_) * cos(y_) end)
    end

    @testset "numbered arguments are ordered by number, not by occurence" begin
        @test @capture(Fn.fn(:(sin(_2) * cos(_1))), function (x_, y_) sin(y_) * cos(x_) end)
    end

    @testset "missing numbers introduce unused extra arguments" begin
        expr = Fn.fn(:(sin(_4) * cos(_1)))
        @test @capture(expr, function (w_, x_, y_, z_) sin(z_) * cos(w_) end)
    end

    @testset "leading zeros are not accepted" begin
        # They would make everything more complicated
        expr = Fn.fn(:(sin(_01) * cos(_02)))
        @test @capture(expr, function () sin(x_) * cos(y_) end)
        # Putting _01 and _02 into the @capture expression confuses MacroTools
        @test x == :_01 && y == :_02
    end

    @testset "mixing _ and _n raises an error" begin
        @test_throws ErrorException Fn.fn(:(_ + _1))
    end

    @testset "@fn does not replace _ in nested @fn" begin
        fn = Fn.fn(:(map(@fn(_ + 3), _)))

        matched = @capture(fn, function (x_) map(@fn(y_ + 3), x_) end)
        @test matched
        @test y === :_
    end

    @testset "@fn escapes function body" begin
        h(x) = x^2
        f = @fn(h(_))
        @test f(5) == 25
    end
end


text:
"""
    call_limit_stop_check(status, information, options)

Limit the number of function evaluations, i.e.,
return `status.f_calls >= options.f_calls_limit`.
"""
function call_limit_stop_check(status, information, options)
    cond =  status.f_calls >= options.f_calls_limit
    options.debug && cond && @info("Stopped since call_limit was met.")
    cond
end

"""
    iteration_stop_check(status, information, options)
Used to limit the number of iterations.
"""
function iteration_stop_check(status, information, options)
    cond =  status.iteration >= options.iterations
    options.debug && cond && @info("Stopped since iteration limit was met.")
    cond
end

"""
    time_stop_check(status, information, options)
Used to limit the time (in seconds), i.e., `status.overall_time >= options.time_limit`.
"""
function time_stop_check(status, information, options)
    cond =  status.overall_time >= options.time_limit
    options.debug && cond && @info("Stopped since time limit was met.")
    cond
end

"""
    accuracy_stop_check(status, information, options)

If the optimum is provided, then check if the accuracy is met via
`abs(status.best_sol.f - information.f_optimum) < options.f_tol`.
"""
function accuracy_stop_check(status::State{xf_indiv}, information, options)
    cond =  !isnan(information.f_optimum) && abs(status.best_sol.f - information.f_optimum) < options.f_tol
    options.debug && cond && @info("Stopped since accuracy was met.")
    cond
end

function accuracy_stop_check(status::State{xfgh_indiv}, information, options)
    isnan(information.f_optimum) && (return false)

    vio = violationsSum(status.best_sol.g, status.best_sol.h, ε=options.h_tol)
    cond =  vio ≈ 0.0 && abs(status.best_sol.f - information.f_optimum) < options.f_tol

    options.debug && cond && @info("Stopped since accuracy and feasibility was met.")
    cond
end

function accuracy_stop_check(status::State{xFgh_indiv}, information, options)
    return false
end

"""
    var_stop_check(status, information, options)

Check if the variance is close to zero in objective space.
"""
function var_stop_check(status, information, options)
    if typeof(status.best_sol.f) <: Array
        return false
    end
    cond =  var(map( s->s.f, status.population )) ≈ 0.0
    options.debug && cond && @info("Stopped since varF was met.")
    cond
end

function diversity_stop_check(status, information, options)
    typeof(status.best_sol.f) <: Array && (return false)

    cond =  var(map( s->s.f, status.population )) ≈ 0.0
    options.debug && cond && @info("Stopped since varF was met.")
    cond
end

function stop_check(status, information, options)
    accuracy_stop_check(status, information, options) ||
    var_stop_check(status, information, options)
end

"""
    diff_check(status, information, options; d = options.f_tol, p = 0.5)

Check the difference between best and worst objective function values in current
population (where at least %p of solution are feasible). Return `true` when such difference
is `<= d`, otherwise return `false`.

> Ref. Zielinski, K., & Laur, R. (n.d.). Stopping Criteria for Differential Evolution in
> Constrained Single-Objective Optimization. Studies in Computational Intelligence,
> 111–138. doi:10.1007/978-3-540-68830-3_4 (https://doi.org/10.1007/978-3-540-68830-3_4)
"""
function diff_check(status, information, options; d = options.f_tol, p = 0.3)
    mask = is_feasible.(status.population)
    p_feasibles = mean(mask)

    # not enough feasible samples?
    p_feasibles < p && (return false)

    fmin = minimum(s -> fval(s), status.population[mask])
    fmax = maximum(s -> fval(s), status.population[mask])

    fmax - fmin <= d
end


text:
#=
Author:= Shreeda Bhat
License:= MIT
UTF-8 =#

################## Data Types in Julia ###################
#=

Julia's basic data types are Int64,Float64,Char and String

=#
a = 1
println(typeof(a))#Int32     Depends on the machine
b = 1+2+3
println(typeof(b))#Int32
c = -2
println(typeof(c))#Int32
d = -100+2
println(typeof(d))#Int32
e = eps(0.0)
println(typeof(e))#Float64
f = 'α'
println(typeof(f))# Greek alphabet, in unicode '' - single quotes is considered as char "" -> String
g = "α"
println(typeof(g))# String

h = "h"
println(typeof(h))#char
i = 'H'
println(typeof(i))#String
j = 1+2+3.0
println(typeof(j))#Float64

#################Let's play a bit with julia##########################
#= You can convert unicode to symbols like alpha beta  + -

how to write unicode in julia

press \ with unicode and press tab =#

# println(/pm )

println("±")

# just like tex and latex

######################Boolean type############################
t = true
println(typeof(t))## in julia it's 8 bits

######################let's play more around julia##############

println(2 == 2.0)#true


text:
module Outliernicers

using Random
using Dates
using DataFrames
using Statistics
using StatsBase: iqr, quantile, sample

using ..ValDateFilters
using ..AbsTypes
using ..Utils
import ..AbsTypes: fit, fit!, transform, transform!

export fit, fit!, transform, transform!
export Outliernicer

"""
    Outliernicer(Dict(
       :dateinterval => Dates.Hour(1),
       :nnsize => 1,
       :missdirection => :symmetric,
       :scale => 1.25
    ))

Detects outliers below or above (median-scale*iqr,median+scale*iqr)
and calls DateValNNer to replace them with nearest neighbors.

Example:

    fname = joinpath(dirname(pathof(TSML)),"../data/testdata.csv")
    csvfilter = CSVDateValReader(Dict(:filename=>fname,:dateformat=>"dd/mm/yyyy HH:MM"))
    valgator = DateValgator(Dict(:dateinterval=>Dates.Hour(1)))
    valnner = DateValNNer(Dict(:dateinterval=>Dates.Hour(1)))
    stfier = Statifier(Dict(:processmissing=>true))
    mono = Monotonicer(Dict())
    outliernicer = Outliernicer(Dict(:dateinterval=>Dates.Hour(1)))

    mpipeline = @pipeline csvfilter |> valgator |> mono |> valnner |> outliernicer |> stfier
    results = fit_transform!(mpipeline)

Implements: `fit!`, `transform!`
"""
mutable struct Outliernicer <: Transformer
   name::String
   model::Dict{Symbol,Any}

   function Outliernicer(args=Dict())
      default_args = Dict{Symbol,Any}(
         :name => "outlrncr",
         :dateinterval => Dates.Hour(1),
         :nnsize => 1,
         :missdirection => :symmetric,
         :scale => 1.25,
         :strict => false
      )
      cargs=nested_dict_merge(default_args,args)
      cargs[:name] = cargs[:name]*"_"*randstring(3)
      new(cargs[:name],cargs)
   end
end

"""
    fit!(st::Outliernicer, features::T, labels::Vector=[])

Check that `features` are two-colum data.
"""
function fit!(st::Outliernicer, features::DataFrame, labels::Vector=[])::Nothing
   ncol(features) == 2 || throw(ArgumentError("dataframe must have 2 columns: Date, Val"))
   return nothing
end

function fit(st::Outliernicer, features::DataFrame, labels::Vector=[])::Outliernicer
   fit!(st,features,labels)
   return deepcopy(st)
end

"""
    transform!(st::Outliernicer, features::T) where {T<:Union{Vector,Matrix,DataFrame}}

Locate outliers based on IQR factor and calls DateValNNer to replace them with nearest neighbors.
"""
function transform!(st::Outliernicer, features::DataFrame)::DataFrame
   features != DataFrame() || return DataFrame()
   ncol(features) == 2 || throw(ArgumentError("dataframe must have 2 columns: Date, Val"))
   sum(names(features) .== ("Date","Value"))  == 2 || throw(ArgumentError("wrong column names"))
   mfeatures=deepcopy(features)
   rvals = mfeatures.Value
   # compute the outlier range
   # setup to store both missing and numbers
   mvals = Array{Union{Missing,eltype(rvals)},1}(missing,length(rvals))
   mvals .= rvals
   crvals = skipmissing(rvals) # stat of non-missing
   miqr = iqr(crvals)
   med = median(crvals) # median
   scale = st.model[:scale]
   lower=med - scale*miqr; upper=med + scale*miqr
   missindx = findall(x -> !ismissing(x) && (x > upper || x < lower),rvals)
   mvals[missindx] .= missing
   mfeatures.Value = mvals
   # use ValNNer to replace missings
   valnner = DateValNNer(st.model)
   fit!(valnner,mfeatures)
   resdf = transform!(valnner,mfeatures)
   resdf.Value = collect(skipmissing(resdf.Value))
   resdf
end

function transform(st::Outliernicer, features::DataFrame)::DataFrame
   return transform!(st,features)
end

end


text:
module TestTypes

using Test
using CategoricalDistributions
using CategoricalArrays
using StableRNGs
using ScientificTypesBase
import Random

# coverage of constructor testing is expanded in the other test files

@testset "constructors" begin
    @test_throws ArgumentError UnivariateFinite(Dict('f'=>0.7, 'q'=>0.2))
    @test_throws ArgumentError UnivariateFinite(Dict('f'=>0.7, 'q'=>0.2),
                                                pool=missing)
    v = categorical(collect("asqfasqffqsaaaa"), ordered=true)
    f = v[4]
    @test_logs((:warn, r"Ignoring"),
               UnivariateFinite(Dict('f'=>0.7, 'q'=>0.3),
                                pool=f, ordered=true))

    @test_logs((:warn, r"No "),
               UnivariateFinite(['f', 'q', 's'],  [0.7, 0.2, 0.1]))

end

@testset "array constructors" begin

    rng = StableRNG(111)
    n   = 10
    c   = 3

    probs  = rand(rng, n)
    supp = ["class1", "class2"]

    UnivariateFinite(supp, probs, pool=missing, augment=true);

    # dimension mismatches:
    badprobs = rand(rng, 40, 3)
    @test_throws(CategoricalDistributions.err_dim(supp, badprobs),
                 UnivariateFinite(supp, badprobs, pool=missing))

    # dimension mismatch, augmented case:
    probs2 = [0.1 0.5 0.1;
              0.3 0.2 0.1]
    supp2 = ["no", "yes", "maybe"]
    @test_throws(CategoricalDistributions.err_dim_augmented(supp2, probs2),
                 UnivariateFinite(supp2, probs2, augment=true, pool=missing))

    # not augmentable:
    @test_throws(CategoricalDistributions.ERR_AUG,
                 UnivariateFinite(["no", "yes", "maybe"],
                                  [0.6 0.5;   # sum exceeding one!
                                   0.3 0.2],
                                  augment=true,
                                  pool=missing))

    # autosupport:
    u = UnivariateFinite(probs, pool=missing, augment=true);
    probs = probs ./ sum(probs)
    u = UnivariateFinite(probs, pool=missing);
    @test u isa UnivariateFinite
    probs = rand(rng, 10, 2)
    probs = probs ./ sum(probs, dims=2)
    u = UnivariateFinite(probs, pool=missing);
    @test u.scitype == Multiclass{2}
    u = UnivariateFinite(probs, pool=missing, augment=true);
    @test u.scitype == Multiclass{3}

    probs = [-1,0,1]
    @test_throws(CategoricalDistributions.ERR_01,
                 UnivariateFinite(probs, pool=missing, augment=true))

    v = categorical(1:3)
    @test_logs((:warn, r"Ignoring"),
               UnivariateFinite(v[1:2], rand(rng, 3),
                                augment=true, pool=missing));
    @test_logs((:warn, r"Ignoring"),
               UnivariateFinite(v[1:2], rand(rng, 3),
                                augment=true, ordered=true));

    # using `UnivariateFiniteArray` as a constructor just falls back
    # to `UnivariateFinite` constructor:
    probs  = rand(rng, n)
    d1 = UnivariateFiniteArray(supp, probs, pool=missing, augment=true)
    d2 = UnivariateFinite(supp, probs, pool=missing, augment=true)
    @test d1.prob_given_ref == d2.prob_given_ref
end

@testset "display" begin
    d = UnivariateFinite(["x", "y"], [0.3, 0.7], pool=missing)
    v = UnivariateFinite(["x", "y"], rand(3), augment=true, pool=missing)
    io = IOBuffer()
    print(io, d, v, [d, ], [v, ])
end

end

true


text:
input_file = "input_file.txt"

function parse_line(line)
    # lines look like this: 8-12 f: fmffhjmfkxfrvcfrbr
    line_chunks = split(line, ":")

    password = strip(line_chunks[2])

    # this part: 8-12 f
    left_chunk = strip(line_chunks[1])

    tmp1 = split(left_chunk, "-")
    min_count = parse(Int32, tmp1[1])

    tmp2 = split(tmp1[2], " ")
    max_count = parse(Int32, tmp2[1])
    letter = only(tmp2[2]) # converts this from substring to character

    min_count, max_count, letter, password
end

function validate_password1(min_count, max_count, letter, password)
    string_array = collect(password)

    letter_count = 0
    for c in string_array
        equality = c == letter
        # println("$c $(typeof(c)) $letter $(typeof(letter)) $equality")
        if equality
            letter_count += 1
        end
    end
    # println(letter_count)
    valid_password = letter_count >= min_count && letter_count <= max_count
end

function validate_password2(pos1, pos2, letter, password)
    string_array = collect(password)
    valid_password = xor(string_array[pos1] == letter, string_array[pos2] == letter)
end

# This should possibly be an enum but they seem clunky and un-idiomatic in Julia?
password_mode = 2
open(input_file) do file
    valid_passwords = 0
    for line in eachline(file)
        min_count, max_count, letter, password = parse_line(line)
        if password_mode == 1
            if validate_password1(min_count, max_count, letter, password)
                valid_passwords += 1
            end
        elseif password_mode == 2
            if validate_password1(min_count, max_count, letter, password)
                valid_passwords += 1
            end
        end
    end
    valid_passwords
end


text:
#File: fMean.jl
#Author:
#Date: 29-June-2020

#MEAN
function fMean(x,f)
#-------------------------------------------------------------------
#x:     returns vector
#f:     reporting frequency e.g. 12 (monthly)
#mn:    annualised mean
#-------------------------------------------------------------------

    fMean = mean(x).*f;
    mn = fMean;
    return mn
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

module InterfaceTests

using JuliaFEM.Testing

function test_foo()
    @test 1+2 == 3
end

end


text:
using ForwardDiff

"Parameters for right-hand side function of formal Lagrangian Runge-Kutta methods."
mutable struct ParametersFLRK{DT,TT,VT,D,S} <: Parameters{DT,TT}
    v::VT
    Δt::TT

    a::Matrix{TT}
    â::Matrix{TT}
    c::Vector{TT}

    t::TT

    q::Vector{DT}
    p::Vector{DT}
end

function ParametersFLRK(DT, D, v::VT, Δt::TT, tab) where {TT,VT}
    ParametersFLRK{DT,TT,VT,D,tab.s}(v, Δt, tab.a, tab.â, tab.c, 0, zeros(DT,D), zeros(DT,D))
end

struct NonlinearFunctionCacheFLRK{DT}
    Q::Vector{Vector{DT}}
    V::Vector{Vector{DT}}
    Y::Vector{Vector{DT}}

    function NonlinearFunctionCacheFLRK{DT}(d, s) where {DT}

        # create internal stage vectors
        Q = create_internal_stage_vector(DT, d, s)
        V = create_internal_stage_vector(DT, d, s)
        Y = create_internal_stage_vector(DT, d, s)

        new(Q, V, Y)
    end
end

function compute_stages!(x::Vector{ST}, Q::Vector{Vector{ST}}, V::Vector{Vector{ST}}, Y::Vector{Vector{ST}},
                                    params::ParametersFLRK{DT,TT,VT,D,S}) where {ST,DT,TT,VT,D,S}
    local tᵢ::TT

    @assert S == length(Q) == length(V) == length(Y)

    for i in 1:S
        @assert D == length(Q[i]) == length(V[i]) == length(Y[i])
        tᵢ = params.t + params.Δt * params.c[i]

        # copy x to Y
        for k in 1:D
            Y[i][k] = x[D*(i-1)+k]
        end

        # compute Q = q + Δt Y
        Q[i] .= params.q .+ params.Δt .* Y[i]

        # compute V = v(Q)
        params.v(tᵢ, Q[i], V[i])
    end
end

"Compute stages of formal Lagrangian Runge-Kutta methods."
@generated function function_stages!(x::Vector{ST}, b::Vector{ST}, params::ParametersFLRK{DT,TT,VT,D,S}) where {ST,DT,TT,VT,D,S}

    cache = NonlinearFunctionCacheFLRK{ST}(D, S)

    quote
        compute_stages!(x, $cache.Q, $cache.V, $cache.Y, params)

        local y1::ST
        local y2::ST

        # compute b = - (Y-AV)
        for i in 1:S
            for k in 1:D
                y1 = 0
                y2 = 0
                for j in 1:S
                    y1 += params.a[i,j] * $cache.V[j][k]
                    y2 += params.â[i,j] * $cache.V[j][k]
                end
                b[D*(i-1)+k] = - $cache.Y[i][k] + (y1 + y2)
            end
        end
    end
end

"Formal Lagrangian Runge-Kutta integrator."
struct IntegratorFLRK{DT, TT, AT, FT, GT, VT, ΩT, dHT, SPT, ST, IT <: InitialGuessODE{DT,TT,VT}, N} <: DeterministicIntegrator{DT,TT}
    equation::VODE{DT,TT,AT,FT,GT,VT,ΩT,dHT,N}
    tableau::TableauFIRK{TT}
    Δt::TT

    params::SPT
    solver::ST
    iguess::IT

    q::Vector{Vector{TwicePrecision{DT}}}
    p::Vector{Vector{TwicePrecision{DT}}}

    v::Vector{DT}
    y::Vector{DT}

    Q::Vector{Vector{DT}}
    V::Vector{Vector{DT}}
    P::Vector{Vector{DT}}
    F::Vector{Vector{DT}}
    G::Vector{Vector{DT}}
    ϑ::Vector{Vector{DT}}
    Y::Vector{Vector{DT}}
    Z::Vector{Vector{DT}}

    J::Vector{Matrix{DT}}
    A::Array{DT,2}
end

function IntegratorFLRK(equation::VODE{DT,TT,AT,FT,GT,VT,ΩT,dHT,N}, tableau::TableauFIRK{TT}, Δt::TT) where {DT,TT,AT,FT,GT,VT,ΩT,dHT,N}
    D = equation.d
    M = equation.n
    S = tableau.q.s

    # create solution vector for internal stages / nonlinear solver
    x = zeros(DT, D*S)

    # create solution vectors
    q = create_solution_vector(DT, D, M)
    p = create_solution_vector(DT, D, M)

    # create velocity and update vector
    v = zeros(DT,D)
    y = zeros(DT,D)

    # create internal stage vectors
    Q = create_internal_stage_vector(DT, D, S)
    V = create_internal_stage_vector(DT, D, S)
    P = create_internal_stage_vector(DT, D, S)
    F = create_internal_stage_vector(DT, D, S)
    G = create_internal_stage_vector(DT, D, S)
    ϑ = create_internal_stage_vector(DT, D, S)
    Y = create_internal_stage_vector(DT, D, S)
    Z = create_internal_stage_vector(DT, D, S)

    J = Array{Matrix{DT}}(undef, S)
    # J = Vector{Matrix{DT}}(S)
    for i in 1:S
        J[i] = zeros(DT,D,D)
    end

    A = zeros(DT,D*S,D*S)

    # create params
    params = ParametersFLRK(DT, D, equation.v, Δt, tableau.q)

    # create rhs function for nonlinear solver
    function_stages = (x,b) -> function_stages!(x, b, params)

    # create solver
    solver = get_config(:nls_solver)(x, function_stages)

    # create initial guess
    iguess = InitialGuessODE(get_config(:ig_interpolation), equation, Δt)

    # create integrator
    IntegratorFLRK{DT, TT, AT, FT, GT, VT, ΩT, dHT, typeof(params), typeof(solver), typeof(iguess), N}(
                                        equation, tableau, Δt, params, solver, iguess,
                                        q, p, v, y, Q, V, P, F, G, ϑ, Y, Z, J, A)
end

function initialize!(int::IntegratorFLRK{DT,TT}, sol::SolutionPODE, m::Int) where {DT,TT}
    @assert m ≥ 1
    @assert m ≤ sol.ni

    # copy initial conditions from solution
    get_initial_conditions!(sol, int.q[m], int.p[m], m)

    # initialise initial guess
    initialize!(int.iguess, m, sol.t[0], int.q[m])
end

function initial_guess!(int::IntegratorFLRK, m::Int)
    # compute initial guess for internal stages
    for i in 1:int.tableau.q.s
        evaluate!(int.iguess, m, int.y, int.v, int.tableau.q.c[i])
        for k in 1:int.equation.d
            int.V[i][k] = int.v[k]
        end
    end
    for i in 1:int.tableau.q.s
        for k in 1:int.equation.d
            int.solver.x[int.equation.d*(i-1)+k] = 0
            for j in 1:int.tableau.q.s
                int.solver.x[int.equation.d*(i-1)+k] += int.tableau.q.a[i,j] * int.V[j][k]
            end
        end
    end
end

"Integrate ODE with fully implicit Runge-Kutta integrator."
function integrate_step!(int::IntegratorFLRK{DT,TT}, sol::SolutionPODE{DT,TT}, m::Int, n::Int) where {DT,TT}
    @assert m ≥ 1
    @assert m ≤ sol.ni

    @assert n ≥ 1
    @assert n ≤ sol.ntime

    # set time for nonlinear solver
    int.params.t  = sol.t[0] + (n-1)*int.Δt
    int.params.q .= int.q[m]
    int.params.p .= int.p[m]

    # compute initial guess
    initial_guess!(int, m)

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute vector field at internal stages
    compute_stages!(int.solver.x, int.Q, int.V, int.Y, int.params)

    # compute final update for q
    update_solution!(int.q[m], int.V, int.tableau.q.b, int.tableau.q.b̂, int.Δt)

    # create temporary arrays
    tV = zeros(DT, int.equation.d)
    δP = zeros(DT, int.equation.d*int.tableau.q.s)

    # compute ϑ = α(Q), V(Q) = int.equation.v(t, Q, V)
    # and f_0(Q, V(Q)) = int.equation.f(t, Q, V, F)
    for i in 1:int.tableau.q.s
        tᵢ = int.params.t + int.Δt * int.tableau.q.c[i]
        int.equation.α(tᵢ, int.Q[i], int.ϑ[i])
        int.equation.v(tᵢ, int.Q[i], int.V[i])
        int.equation.g(tᵢ, int.Q[i], int.V[i], int.F[i])
    end

    # compute Jacobian of v via ForwardDiff
    for i in 1:int.tableau.q.s
        tᵢ = int.params.t + int.Δt * int.tableau.q.c[i]
        v_rev! = (v,q) -> int.equation.v(tᵢ,q,v)
        ForwardDiff.jacobian!(int.J[i], v_rev!, tV, int.Q[i])
    end

    # contract J with ϑ and add to G
    for l in 1:int.tableau.q.s
        for i in 1:int.equation.d
            int.G[l][i] = 0
            for j in 1:int.equation.d
                int.G[l][i] += int.ϑ[l][j] * int.J[l][j,i]
            end
        end
    end

    # solve linear system AP=δP for P

    # compute δP
    for l in 1:int.tableau.q.s
        for i in 1:int.equation.d
            # set δP = p
            δP[(l-1)*int.equation.d+i] = int.params.p[i]
            # add A(F+G) to δP
            for k in 1:int.tableau.q.s
                δP[(l-1)*int.equation.d+i] += int.Δt * int.tableau.q.a[l,k] * (int.F[k][i] + int.G[k][i])
            end
        end
    end

    # construct A = identity(sd×sd) + A ⊗ J
    for k in 1:int.tableau.q.s
        for l in 1:int.tableau.q.s
            for i in 1:int.equation.d
                for j in 1:int.equation.d
                    int.A[(k-1)*int.equation.d+i, (l-1)*int.equation.d+j] = int.Δt * int.tableau.q.a[k,l] * int.J[l][j,i]
                end
            end
        end
    end
    for i in 1:int.equation.d*int.tableau.q.s
        int.A[i,i] += one(DT)
    end

    # solve AP = δP and copy result to int.P
    lu = LUSolver(int.A, δP)
    factorize!(lu)
    solve!(lu)
    tP = reshape(lu.x, (int.equation.d, int.tableau.q.s))
    for l in 1:int.tableau.q.s
        int.P[l] .= tP[:,l]
    end

    # contract J with P and subtract from G, so that G = (ϑ-P)J
    for l in 1:int.tableau.q.s
        for i in 1:int.equation.d
            int.G[l][i] = 0
            for j in 1:int.equation.d
                int.G[l][i] += (int.ϑ[l][j] - int.P[l][j]) * int.J[l][j,i]
            end
        end
    end

    # println(int.ϑ)
    # println(int.P)
    # println(int.ϑ .- int.P)
    # println()

    # compute final update for p
    update_solution!(int.p[m], int.F, int.tableau.q.b, int.tableau.q.b̂, int.Δt)
    update_solution!(int.p[m], int.G, int.tableau.q.b, int.tableau.q.b̂, int.Δt)

    # copy solution to initial guess
    update!(int.iguess, m, sol.t[0] + n*int.Δt, int.q[m])
end


text:
# Checks for various fairness criteria.

"""
    check_complete(A)

Check that the allocation is complete, or effective, in the sense that each item
has been allocated to at least one agent.
"""
function check_complete(A)
    for g in items(A)
        length(owners(A, g)) >= 1 || return false
    end
    return true
end

"""
    check_partition(A)

Check that the allocation is a partition, i.e., that each item has been
allocated to exactly one agent.
"""
function check_partition(A)
    for g in items(A)
        length(owners(A, g)) == 1 || return false
    end
    return true
end

# General check for EF-like properties (EF_ = EF, EF1, EFX, ...), where the
# value of the other agent's bundle is somehow modified, e.g., by removing
# some item. This modified value is provided by the function value_:
function check_ef_(V, A, value_)

    N = agents(V)

    for i in N, j in N

        i !== j || continue

        if value(V, i, A) < value_(V, i, bundle(A, j))
            return false
        end

    end

    return true

end

"""
    check_ef(V, A)

Check whether the allocation `A` is *envy-free* for the profile `V`, i.e., if no
agent strictly prefers another agent's bundle.
"""
check_ef(V, A) = check_ef_(V, A, value)

"""
    check_ef1(V, A)

Check whether the allocation `A` is *envy-free up to one item* for the profile
`V`, i.e., if no agent strictly prefers another agent's bundle, given that an
appropriate (e.g., the most valuable) item is removed.
"""
check_ef1(V, A) = check_ef_(V, A, value_1)

"""
    check_efx(V, A)

Check whether the allocation `A` is *envy-free up to any item* for the profile
`V`, i.e., if no agent strictly prefers another agent's bundle, given that an
appropriate (e.g., the least valuable) item is removed.
"""
check_efx(V, A) = check_ef_(V, A, value_x)

"""
    check(V, A, C)

Check that the allocation `A` obeys the `Constraint` `C`, given the profile `V`.
"""
function check end

"""
    check(V, A, C::Counts)

Check whether the allocation `A` respects the cardinality constraints `C`.
"""
function check(V, A, C::Counts)

    for c in C

        counts = zeros(Int, na(A))

        for g in c, i in owners(A, g)
            counts[i] += 1
        end

        maximum(counts) <= c.threshold || return false

    end

    return true

end

"""
    check(V, A, C::Conflicts)

Check whether the allocation `A` respects the item conflicts `C`.
"""
function check(V, A, C::Conflicts)

    G = graph(C)

    for e in edges(G)
        isempty(owners(A, src(e)) ∩ owners(A, dst(e))) || return false
    end

    return true

end


text:
include("libavutil_h.jl")

#include("audioconvert.jl")
#include("avutil.jl")
#include("dict.jl")
#include("fifo.jl")
#include("file.jl")
#include("imgutils.jl")
include("log.jl")
include("mem.jl")
include("opt.jl")
include("pixdesc.jl")
#include("rational.jl")
#include("samplefmt.jl")

### Manually define av_pix_fmt_desc_get to match later versions
export av_pix_fmt_desc_get

const p_av_pix_fmt_descriptors = cglobal((:av_pix_fmt_descriptors, :libavutil), AVPixFmtDescriptor)
av_pix_fmt_desc_get(n::Integer) = p_av_pix_fmt_descriptors + n*sizeof(AVPixFmtDescriptor)


text:
__precompile__()

"""
    FunManifolds

Main module for `FunManifolds.jl` -- a Julia package for functional differential geometry.
"""
module FunManifolds

using LinearAlgebra
using Interpolations
using Manifolds
using ManifoldsBase
using Markdown: @doc_str
using QuadGK
using StaticArrays

import Base: +, -, *, isapprox

import ManifoldsBase:
    allocate,
    allocate_result,
    allocate_result_type,
    allocation_promotion_function,
    array_value,
    base_manifold,
    check_point,
    check_point__transparent,
    check_vector,
    decorated_manifold,
    decorator_transparent_dispatch,
    default_decorator_dispatch,
    distance,
    embed,
    embed!,
    exp,
    exp!,
    exp!__intransparent,
    geodesic,
    # TODO: uncomment the import if `flat!` goes to ManifoldsBase
    # flat!__intransparent,
    get_basis,
    get_coordinates,
    get_coordinates!,
    get_embedding,
    get_vector,
    get_vector!,
    get_vectors,
    injectivity_radius,
    inner,
    inner__intransparent,
    is_point,
    is_vector,
    inverse_retract,
    inverse_retract!,
    log,
    log!,
    manifold_dimension,
    mid_point,
    mid_point!,
    norm,
    number_eltype,
    number_of_coordinates,
    project,
    project!,
    representation_size,
    retract,
    retract!,
    shortest_geodesic,
    vector_transport_direction,
    vector_transport_direction!,
    vector_transport_to,
    vector_transport_to!,
    zero_vector,
    zero_vector!

import Manifolds:
    apply,
    apply!,
    compose,
    compose!,
    get_iterator,
    identity_element,
    identity_element!,
    inv,
    inv!,
    optimal_alignment,
    optimal_alignment!,
    zero_vector

using Manifolds: _read, _write, AbstractRiemannianDiffBackend

mutable struct GeneralParams
    quad_rel_tol::Union{Real,Nothing}
    quad_abs_tol::Union{Real,Nothing}
end

const global PARAMS = GeneralParams(nothing, nothing)

function rtoldefault(M::AbstractManifold, x1, x2)
    return 1.e-8
end

function atoldefault(M::AbstractManifold, x1, x2)
    return 0.0
end

function concretize_tols(M::AbstractManifold, x1, x2; reltol=nothing, abstol=nothing)
    rtol = if reltol === nothing
        rtoldefault(M, x1, x2)
    else
        reltol
    end

    atol = if abstol === nothing
        atoldefault(M, x1, x2)
    else
        abstol
    end

    return (rtol, atol)
end

struct ProjectedDifferenceBackend{TDT<:Union{Number,Nothing}} <:
       AbstractRiemannianDiffBackend
    dt::TDT
end

ProjectedDifferenceBackend() = ProjectedDifferenceBackend{Float64}(1e-7)

include("interpolation.jl")

include("DiscretizedCurves.jl")
include("FunctionCurve.jl")
include("CurveWarping.jl")
include("CurveWarpingSRSF.jl")
include("curve_warping_alignment.jl")
include("functional_transformations.jl")

export CurveWarpingAction,
    CurveWarpingGroup,
    CurveWarpingSpace,
    CurveWarpingSRSFAction,
    CurveWarpingSRSFGroup,
    CurveWarpingSRSFSpace,
    DiscretizedCurves,
    FunctionCurveSpace,
    reverse_srsf,
    reverse_srvf,
    srsf,
    srvf,
    tsrvf,
    transport_srvf,
    transport_srvf!,
    UniformDiscretizedCurves

end #module