text:
export rule, marginalrule
export @rule, @marginalrule
export @call_rule, @call_marginalrule

using MacroTools
using .MacroHelpers

import Base: showerror

"""
    rule(fform, on, vconstraint, mnames, messages, qnames, marginals, meta, __node)

This function is used to compute an outbound message for a given node

# Arguments

- `fform`: Functional form of the node in form of a type of the node, e.g. `::Type{ <: NormalMeanVariance }` or `::typeof(+)`
- `on`: Outbound interface's tag for which a message has to be computed, e.g. `::Type{ Val{:out} }` or `::Type{ Val{:μ} }`
- `vconstraint`: Variable constraints for an outbound interface, e.g. `Marginalisation` or `MomentMatching`
- `mnames`: Ordered messages names in form of the Val type, eg. `::Type{ Val{ (:mean, :precision) } }`
- `messages`: Tuple of message of the same length as `mnames` used to compute an outbound message
- `qnames`: Ordered marginal names in form of the Val type, eg. `::Type{ Val{ (:mean, :precision) } }`
- `marginals`: Tuple of marginals of the same length as `qnames` used to compute an outbound message
- `meta`: Extra meta information
- `__node`: Node reference

See also: [`@rule`](@ref), [`marginalrule`](@ref), [`@marginalrule`](@ref)
"""
function rule end

"""
    marginalrule(fform, on, mnames, messages, qnames, marginals, meta, __node)

This function is used to compute a local joint marginal for a given node

# Arguments

- `fform`: Functional form of the node in form of a type of the node, e.g. `::Type{ <: NormalMeanVariance }` or `::typeof(+)`
- `on`: Local joint marginal tag , e.g. `::Type{ Val{ :mean_precision } }` or `::Type{ Val{ :out_mean_precision } }`
- `mnames`: Ordered messages names in form of the Val type, eg. `::Type{ Val{ (:mean, :precision) } }`
- `messages`: Tuple of message of the same length as `mnames` used to compute an outbound message
- `qnames`: Ordered marginal names in form of the Val type, eg. `::Type{ Val{ (:mean, :precision) } }`
- `marginals`: Tuple of marginals of the same length as `qnames` used to compute an outbound message
- `meta`: Extra meta information
- `__node`: Node reference

See also: [`rule`](@ref), [`@rule`](@ref) [`@marginalrule`](@ref)
"""
function marginalrule end

# Macro code

"""
    rule_macro_parse_on_tag(expression)

Do not use this function directly. This function is private and does not belong to the public API.

This function is used to parse an `on` tag for message rules and marginal rules specification.

```
@rule MvNormalMeanPrecision(:out, Marginalisation) (...) = begin
                            ^^^^
                            `on` tag
    ...
end
```

or

```
@rule NormalMixture((:m, k), Marginalisation) (...) = begin
                    ^^^^^^^
                    `on` tag
    ...
end
```

Accepts either a quoted symbol expressions or a (name, index) tuple expression. Returns name expression, index expression and index initialisation expression.

See also: [`@rule`](@ref)
"""
function rule_macro_parse_on_tag(on)
    if @capture(on, :name_)
        # First we check on just quoted symbol expression
        # If captures index exression and index initilisation expression are `nothing`
        return :(Type{ Val{ $(QuoteNode(name)) } }), nothing, nothing
    elseif @capture(on, (:name_, index_Symbol))
        return :(Tuple{ Val{$(QuoteNode(name))}, Int}), index, :($index = on[2])
    else
        error("Error in macro. `on` specification is incorrect: $(on). Must be ither a quoted symbol expression (e.g. `:out` or `:mean`) or tuple expression with quoted symbol and index identifier (e.g. `(:m, k)` or `(:w, k)`)")
    end
end

"""
    rule_macro_parse_fn_args(inputs; specname, prefix, proxy)

Do not use this function directly. This function is private and does not belong to the public API.

This function is used to parse an `arguments` tuple for message rules and marginal rules specification.

```
@rule MvNormalMeanPrecision(:out, Marginalisation) (m_μ::NormalMeanPrecision, m_τ::PointMass) = begin
                                                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
                                                   `arguments` vector
    ...
end
```

Accepts a vector of (name, type) elements, specname, name prefix and proxy type.
Returns parsed names without prefix, proxied types and initialisation code block.

See also: [`@rule`](@ref)
"""
function rule_macro_parse_fn_args(inputs; specname, prefix, proxy)
    # First we filter out only prefixed arguments
    finputs = filter((i) -> startswith(string(first(i)), string(prefix)), inputs)

    lprefix = length(string(prefix))

    # We extract names and types
    names   = map(first, finputs)
    types   = map(last, finputs)

    # Check that all arguments have proper names and not only a single prefix in it
    @assert all((n) -> length(string(n)) > lprefix, names)  || error("Empty $(specname) name found in arguments")

    # Initialisation block is simply a `getdata` call from `specname` for each argument
    init_block = map(enumerate(names)) do (index, iname)
        return :($(iname) = getdata($(specname)[$(index)]))
    end

    # We return names in form of a `Type{ Val{ (:name1, :name2, ...) } }`
    # We return types in form of a `Tuple{ ProxyType{ <: type1 }, ProxyType{ <: type2 } }`
    out_names = isempty(names) ? :Nothing : :(Type{ Val{ $(Expr(:tuple, map(n -> QuoteNode(Symbol(string(n)[(lprefix + 1):end])), names)...)) } })
    out_types = isempty(types) ? :Nothing : :(Tuple{ $(map((t) -> MacroHelpers.proxy_type(proxy, t), types)...) })

    return out_names, out_types, init_block
end

"""
    call_rule_macro_parse_fn_args(inputs; specname, prefix, proxy)

Do not use this function directly. This function is private and does not belong to the public API.

This function is used to parse an `arguments` tuple for message and marginal calling rules specification.

```
@call_rule MvNormalMeanPrecision(:out, Marginalisation) (m_μ = NormalMeanPrecision(...), m_τ = PointMass(...)) = begin
                                                        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
                                                        `arguments` vector
    ...
end
```

Accepts a vector of (name, vale) elements, specname, name prefix and proxy type.
Returns parsed names without prefix and proxied values

See also: [`@rule`](@ref)
"""
function call_rule_macro_parse_fn_args(inputs; specname, prefix, proxy)
    finputs  = filter((i) -> startswith(string(first(i)), string(prefix)), inputs)

    lprefix = length(string(prefix))
    names   = map(first, finputs)
    values  = map(last, finputs)

    @assert all((n) -> length(string(n)) > lprefix, names)  || error("Empty $(specname) name found in arguments")

    # Tuples are special cases
    function apply_proxy(any, proxy)
        if any isa Expr && any.head === :tuple
            output      = Expr(:tuple)
            output.args = map(v -> apply_proxy(v, proxy), any.args)
            return output
        end
        return :($(proxy)($any, false, false))
    end

    names_arg  = isempty(names) ? :nothing : :(Val{ $(Expr(:tuple, map(n -> QuoteNode(Symbol(string(n)[(lprefix + 1):end])), names)...)) })
    values_arg = isempty(names) ? :nothing : :($(map(v -> apply_proxy(v, proxy), values)...), )

    return names_arg, values_arg
end

function rule_function_expression(body::Function, fuppertype, on_type, vconstraint, m_names, m_types, q_names, q_types, metatype, whereargs)
    return quote
        function ReactiveMP.rule(
            fform           :: $(fuppertype),
            on              :: $(on_type),
            vconstraint     :: $(vconstraint),
            messages_names  :: $(m_names),
            messages        :: $(m_types),
            marginals_names :: $(q_names),
            marginals       :: $(q_types),
            meta            :: $(metatype),
            __node
        ) where { $(whereargs...) }
            $(body())
        end
    end
end

function marginalrule_function_expression(body::Function, fuppertype, on_type, m_names, m_types, q_names, q_types, metatype, whereargs)
    return quote
        function ReactiveMP.marginalrule(
            fform           :: $(fuppertype),
            on              :: $(on_type),
            messages_names  :: $(m_names),
            messages        :: $(m_types),
            marginals_names :: $(q_names),
            marginals       :: $(q_types),
            meta            :: $(metatype),
            __node
        ) where { $(whereargs...) }
            $(body())
        end
    end
end

import .MacroHelpers

"""
    Documentation placeholder
"""
macro rule(fform, lambda)
    @capture(fform, fformtype_(on_, vconstraint_, options__)) ||
        error("Error in macro. Functional form specification should in the form of 'fformtype_(on_, vconstraint_, options__)'")

    @capture(lambda, (args_ where { whereargs__ } = body_) | (args_ = body_)) ||
        error("Error in macro. Lambda body specification is incorrect")

    @capture(args, (inputs__, meta::metatype_) | (inputs__, )) ||
        error("Error in macro. Lambda body arguments speicifcation is incorrect")

    fuppertype                       = MacroHelpers.upper_type(fformtype)
    on_type, on_index, on_index_init = rule_macro_parse_on_tag(on)
    whereargs                        = whereargs === nothing ? [] : whereargs
    metatype                         = metatype === nothing ? :Any : metatype

    options = map(options) do option
        @capture(option, name_ = value_) || error("Error in macro. Option specification '$(option)' is incorrect")x
        return (name, value)
    end

    inputs = map(inputs) do input
        @capture(input, iname_::itype_) || error("Error in macro. Input $(input) is incorrect")
        return (iname, itype)
    end

    m_names, m_types, m_init_block = rule_macro_parse_fn_args(inputs, specname = :messages, prefix = :m_, proxy = :(ReactiveMP.Message))
    q_names, q_types, q_init_block = rule_macro_parse_fn_args(inputs, specname = :marginals, prefix = :q_, proxy = :(ReactiveMP.Marginal))

    output = quote
        $(
            rule_function_expression(fuppertype, on_type, vconstraint, m_names, m_types, q_names, q_types, metatype, whereargs) do
                return quote
                    $(on_index_init)
                    $(m_init_block...)
                    $(q_init_block...)
                    $(body)
                end
            end
        )
    end

    foreach(options) do option

        # Symmetrical option
        if first(option) === :symmetrical
            @capture(last(option), [ names__ ]) ||
                error("Invalid symmetrical names specification. Name should be a symbol.")

            @assert length(names) > 1 "Invalid symmetrical names specification. Length of names should be greater than 1."

            names = map(names) do name
                @capture(name, :s_) || error("Invalid symmetrical name specification. Name should be a symbol.")
                return s
            end

            indices = map(names) do name
                return findfirst(input -> isequal(string("m_", name), string(first(input))) || isequal(string("q_", name), string(first(input))), inputs)
            end

            foreach(enumerate(indices)) do (i, index)
                @assert index !== nothing "Name $(names[i]) does not exist in arguments list"
            end

            prefixes = map(index -> string(first(inputs[index]))[1:2], indices)

            @assert length(Set(prefixes)) === 1 "It is not possible to mix symmetric arguments from messages and marginals"

            prefix = first(prefixes)
            swaps  = Iterators.flatten(map(i -> map(j -> (indices[i], indices[j]), i+1:length(indices)), 1:length(indices)))

            is_messages  = isequal(prefix, "m_")
            is_marginals = isequal(prefix, "q_")

            foreach(swaps) do swap

                messages  = is_messages ? MacroHelpers.rearranged_tuple(:messages, length(m_init_block), swap) : :(messages)
                marginals = is_marginals ? MacroHelpers.rearranged_tuple(:marginals, length(m_init_block), swap) : :(marginals)

                swapped_m_types = is_messages ? :(Tuple{$(swapped(m_types.args[2:end], first(swap), last(swap))...)}) : m_types
                swapped_q_types = is_marginals ? :(Tuple{$(swapped(q_types.args[2:end], first(swap), last(swap))...)}) : q_types

                @assert !is_messages || swapped_m_types != m_types "Message types are the same after arguments swap for indices = $(swap)"
                @assert !is_marginals || swapped_q_types != q_types "Marginal types are the same after arguments swap for indices = $(swap)"

                output = quote
                    $output
                    $(
                        rule_function_expression(fuppertype, on_type, vconstraint, m_names, swapped_m_types, q_names, swapped_q_types, metatype, whereargs) do
                            return :(ReactiveMP.rule(fform, on, vconstraint, messages_names, $(messages), marginals_names, $(marginals), meta, __node))
                        end
                    )
                end
            end
        else
            error("Unknown option: $(first(option)) in rule specification")
        end

    end

    return esc(output)
end

macro call_rule(fform, args)
    @capture(fform, fformtype_(on_, vconstraint_)) ||
        error("Error in macro. Functional form specification should in the form of 'fformtype_(on_, vconstraint_)'")

    @capture(args, (inputs__, meta = meta_) | (inputs__, )) ||
        error("Error in macro. Arguments specification is incorrect")

    fuppertype                       = MacroHelpers.upper_type(fformtype)
    fbottomtype                      = MacroHelpers.bottom_type(fformtype)
    on_type, on_index, on_index_init = rule_macro_parse_on_tag(on)

    inputs = map(inputs) do input
        @capture(input, iname_ = ivalue_) || error("Error in macro. Argument $(input) is incorrect")
        return (iname, ivalue)
    end

    m_names_arg, m_values_arg = call_rule_macro_parse_fn_args(inputs, specname = :messages, prefix = :m_, proxy = :(ReactiveMP.Message))
    q_names_arg, q_values_arg = call_rule_macro_parse_fn_args(inputs, specname = :marginals, prefix = :q_, proxy = :(ReactiveMP.Marginal))

    on_arg = MacroHelpers.bottom_type(on_type)

    output = quote
        ReactiveMP.rule($fbottomtype, $on_arg, $(vconstraint)(), $m_names_arg, $m_values_arg, $q_names_arg, $q_values_arg, $meta, nothing)
    end

    return esc(output)
end

"""
    Documentation placeholder
"""
macro test_rules(options, on, test_sequence)

    @capture(options, [ option_entries__ ]) || error("Invalid options specification. Options should be in the form on an array.")

    with_float_conversions = false

    float64_atol  = 1e-12
    float32_atol  = 1e-6
    bigfloat_atol = 1e-12

    foreach(option_entries) do option_entry
        @capture(option_entry, (key_ = value_)) || error("Invalid option entry specification: $(option_entry). Option entry should be in the form of a 'key = value' pair.")
        if key === :with_float_conversions
            if value === :true
                with_float_conversions = true
            elseif value === :false
                with_float_conversions = false
            else
                error("Unknown value $(value) for option $(key). Can be either true or false.")
            end
        elseif key === :atol
            float64_atol  = float(value)
            float32_atol  = float(value)
            bigfloat_atol = float(value)
        elseif key === :float64_atol
            float64_atol = float(value)
        elseif key === :float32_atol
            float32_atol = float(value)
        elseif key === :bigfloat_atol
            bigfloat_atol = float(value)
        else
            error("Unknown option $(key)")
        end
    end

    @capture(test_sequence, [ test_sequence_entries__ ]) || error("Invalid test sequence specification. Test sequence should be in the form of an array.")

    block      = Expr(:block)
    block.args = map(test_sequence_entries) do test_entry
        @capture(test_entry, (input = input_, output = output_)) || error("Invalid test entry specification: $(test_entry). Test entry should be in the form of a named tuple (input = ..., output = ...).")

        test_rule      = Expr(:block)
        test_output_s  = gensym()
        test_rule.args = [
            quote
                begin
                    local $test_output_s = ReactiveMP.@call_rule($on, $input)
                    @test ReactiveMP.custom_isapprox( $test_output_s, $output; atol = $float64_atol)
                    @test ReactiveMP.is_typeof_equal($test_output_s, $output)
                end
            end
        ]

        if with_float_conversions
            @capture(input, (input_entries__, )) || error("Invalid input entries. Input entries should be in the form of a named tuple. ")

            # We filter out indices only for inputs that start with 'm_' or 'q_'
            inputs = map(first, filter(collect(enumerate(input_entries))) do i
                @capture(i[2], (key_ = value_))
                if key !== nothing
                    skey = string(key)
                    return startswith(skey, "m_") || startswith(skey, "q_")
                end
                return false
            end)

            function powerset(x::Vector{T}) where T
                result = Vector{T}[[]]
                for elem in x, j in eachindex(result)
                    push!(result, [result[j] ; elem])
                end
                result
            end

            # Here we create all subsets of a input set, to modify their eltype
            indices_power_set = filter(!isempty, powerset(inputs))

            # We create a modified testset for Float32 inputs
            modified_f32_inputs = map(indices_power_set) do set
                cinput = deepcopy(input)
                for index in set
                    cinput.args[index].args[2] = MacroHelpers.expression_convert_eltype(Float32, cinput.args[index].args[2])
                end
                return (cinput, length(set) === length(inputs))
            end

            for m_f32_input in modified_f32_inputs
                m_f32_output = m_f32_input[2] ? MacroHelpers.expression_convert_eltype(Float32, output) : output
                output_s = gensym()
                push!(test_rule.args, quote
                    begin
                        local $output_s = ReactiveMP.@call_rule($on, $(m_f32_input[1]))
                        @test ReactiveMP.custom_isapprox($output_s, $m_f32_output; atol = $float32_atol)
                        @test ReactiveMP.is_typeof_equal($output_s, $m_f32_output)
                    end
                end)
            end

            # We create a modified testset for BigFloat inputs
            modified_bigf_inputs = map(indices_power_set) do set
                cinput = deepcopy(input)
                for index in set
                    cinput.args[index].args[2] = MacroHelpers.expression_convert_eltype(BigFloat, cinput.args[index].args[2])
                end
                return (cinput, true)
            end

            for m_bigf_input in modified_bigf_inputs
                m_bigf_output = m_bigf_input[2] ? MacroHelpers.expression_convert_eltype(BigFloat, output) : output
                output_s = gensym()
                push!(test_rule.args, quote
                    begin
                        local $output_s = ReactiveMP.@call_rule($on, $(m_bigf_input[1]))
                        @test ReactiveMP.custom_isapprox($output_s, $m_bigf_output; atol = $bigfloat_atol)
                        @test ReactiveMP.is_typeof_equal($output_s, $m_bigf_output)
                    end
                end)
            end
        end

        return test_rule
    end

    return esc(block)
end

macro test_rules(on, test_sequence)
    return :(@test_rules [] $on $test_sequence)
end

"""
    Documentation placeholder
"""
macro marginalrule(fform, lambda)
    @capture(fform, fformtype_(on_)) ||
        error("Error in macro. Functional form specification should in the form of 'fformtype_(on_)'")

    @capture(lambda, (args_ where { whereargs__ } = body_) | (args_ = body_)) ||
        error("Error in macro. Lambda body specification is incorrect")

    @capture(args, (inputs__, meta::metatype_) | (inputs__, )) ||
        error("Error in macro. Lambda body arguments speicifcation is incorrect")

    fuppertype                       = MacroHelpers.upper_type(fformtype)
    on_type, on_index, on_index_init = rule_macro_parse_on_tag(on)
    whereargs                        = whereargs === nothing ? [] : whereargs
    metatype                         = metatype === nothing ? :Any : metatype

    inputs = map(inputs) do input
        @capture(input, iname_::itype_) || error("Error in macro. Input $(input) is incorrect")
        return (iname, itype)
    end

    m_names, m_types, m_init_block = rule_macro_parse_fn_args(inputs, specname = :messages, prefix = :m_, proxy = :(ReactiveMP.Message))
    q_names, q_types, q_init_block = rule_macro_parse_fn_args(inputs, specname = :marginals, prefix = :q_, proxy = :(ReactiveMP.Marginal))

    output = quote
        $(
            marginalrule_function_expression(fuppertype, on_type, m_names, m_types, q_names, q_types, metatype, whereargs) do
                return quote
                    $(on_index_init)
                    $(m_init_block...)
                    $(q_init_block...)
                    $(body)
                end
            end
        )
    end

    return esc(output)
end

macro call_marginalrule(fform, args)
    @capture(fform, fformtype_(on_)) ||
        error("Error in macro. Functional form specification should in the form of 'fformtype_(on_)'")

    @capture(args, (inputs__, meta = meta_) | (inputs__, )) ||
        error("Error in macro. Arguments specification is incorrect")

    fuppertype                       = MacroHelpers.upper_type(fformtype)
    fbottomtype                      = MacroHelpers.bottom_type(fformtype)
    on_type, on_index, on_index_init = rule_macro_parse_on_tag(on)

    inputs = map(inputs) do input
        @capture(input, iname_ = ivalue_) || error("Error in macro. Argument $(input) is incorrect")
        return (iname, ivalue)
    end

    m_names_arg, m_values_arg = call_rule_macro_parse_fn_args(inputs, specname = :messages, prefix = :m_, proxy = :(ReactiveMP.Message))
    q_names_arg, q_values_arg = call_rule_macro_parse_fn_args(inputs, specname = :marginals, prefix = :q_, proxy = :(ReactiveMP.Marginal))

    on_arg = MacroHelpers.bottom_type(on_type)

    output = quote
        ReactiveMP.marginalrule($fbottomtype, $on_arg, $m_names_arg, $m_values_arg, $q_names_arg, $q_values_arg, $meta, nothing)
    end

    return esc(output)
end

"""
    Documentation placeholder
"""
macro test_marginalrules(options, on, test_sequence)

    @capture(options, [ option_entries__ ]) || error("Invalid options specification. Options should be in the form on an array.")

    with_float_conversions = false

    float64_atol  = 1e-12
    float32_atol  = 1e-6
    bigfloat_atol = 1e-12

    foreach(option_entries) do option_entry
        @capture(option_entry, (key_ = value_)) || error("Invalid option entry specification: $(option_entry). Option entry should be in the form of a 'key = value' pair.")
        if key === :with_float_conversions
            if value === :true
                with_float_conversions = true
            elseif value === :false
                with_float_conversions = false
            else
                error("Unknown value $(value) for option $(key). Can be either true or false.")
            end
        elseif key === :atol
            float64_atol  = float(value)
            float32_atol  = float(value)
            bigfloat_atol = float(value)
        elseif key === :float64_atol
            float64_atol = float(value)
        elseif key === :float32_atol
            float32_atol = float(value)
        elseif key === :bigfloat_atol
            bigfloat_atol = float(value)
        else
            error("Unknown option $(key)")
        end
    end

    @capture(test_sequence, [ test_sequence_entries__ ]) || error("Invalid test sequence specification. Test sequence should be in the form of an array.")

    block      = Expr(:block)
    block.args = map(test_sequence_entries) do test_entry
        @capture(test_entry, (input = input_, output = output_)) || error("Invalid test entry specification: $(test_entry). Test entry should be in the form of a named tuple (input = ..., output = ...).")

        test_rule      = Expr(:block)
        test_output_s  = gensym()
        test_rule.args = [
            quote
                begin
                    local $test_output_s = ReactiveMP.@call_marginalrule($on, $input)
                    @test ReactiveMP.custom_isapprox($test_output_s, $output; atol = $float64_atol)
                    @test ReactiveMP.is_typeof_equal($test_output_s, $output)
                end
            end
        ]

        if with_float_conversions
            @capture(input, (input_entries__, )) || error("Invalid input entries. Input entries should be in the form of a named tuple. ")

            # We filter out indices only for inputs that start with 'm_' or 'q_'
            inputs = map(first, filter(collect(enumerate(input_entries))) do i
                @capture(i[2], (key_ = value_))
                if key !== nothing
                    skey = string(key)
                    return startswith(skey, "m_") || startswith(skey, "q_")
                end
                return false
            end)

            function powerset(x::Vector{T}) where T
                result = Vector{T}[[]]
                for elem in x, j in eachindex(result)
                    push!(result, [result[j] ; elem])
                end
                result
            end

            # Here we create all subsets of a input set, to modify their eltype
            indices_power_set = filter(!isempty, powerset(inputs))

            # We create a modified testset for Float32 inputs
            modified_f32_inputs = map(indices_power_set) do set
                cinput = deepcopy(input)
                for index in set
                    cinput.args[index].args[2] = MacroHelpers.expression_convert_eltype(Float32, cinput.args[index].args[2])
                end
                return (cinput, length(set) === length(inputs))
            end

            for m_f32_input in modified_f32_inputs
                m_f32_output = m_f32_input[2] ? MacroHelpers.expression_convert_eltype(Float32, output) : output
                output_s = gensym()
                push!(test_rule.args, quote
                    begin
                        local $output_s = ReactiveMP.@call_marginalrule($on, $(m_f32_input[1]))
                        @test ReactiveMP.custom_isapprox($output_s, $m_f32_output; atol = $float32_atol)
                        # @test ReactiveMP.is_typeof_equal($output_s, $m_f32_output) # broken
                    end
                end)
            end

            # We create a modified testset for BigFloat inputs
            modified_bigf_inputs = map(indices_power_set) do set
                cinput = deepcopy(input)
                for index in set
                    cinput.args[index].args[2] = MacroHelpers.expression_convert_eltype(BigFloat, cinput.args[index].args[2])
                end
                return (cinput, true)
            end

            for m_bigf_input in modified_bigf_inputs
                m_bigf_output = m_bigf_input[2] ? MacroHelpers.expression_convert_eltype(BigFloat, output) : output
                output_s = gensym()
                push!(test_rule.args, quote
                    begin
                        local $output_s = ReactiveMP.@call_marginalrule($on, $(m_bigf_input[1]))
                        @test ReactiveMP.custom_isapprox($output_s, $m_bigf_output; atol = $bigfloat_atol)
                        # @test ReactiveMP.is_typeof_equal($output_s, $m_bigf_output) # broken
                    end
                end)
            end
        end

        return test_rule
    end

    return esc(block)
end

macro test_marginalrules(on, test_sequence)
    return :(@test_marginalrules [] $on $test_sequence)
end

# Errors

mutable struct NodeErrorStub
    counter :: Int
end

NodeErrorStub() = NodeErrorStub(0)

function interfaceindex(stub::NodeErrorStub, iname::Symbol)
    stub.counter = stub.counter + 1
    return stub.counter
end

function interfaces(stub::NodeErrorStub)
    return fill(nothing, stub.counter)
end

rule_method_error_extract_fform(f::Function) = string("typeof(", f, ")")
rule_method_error_extract_fform(f)           = string(f)

rule_method_error_extract_on(::Type{ Val{ T } })         where T = T
rule_method_error_extract_on(on::Tuple{ Val{ T }, Int }) where T = string("(:", rule_method_error_extract_on(typeof(on[1])), ", k)")

rule_method_error_extract_vconstraint(something) = typeof(something)

rule_method_error_extract_names(::Type{ Val{ T } }) where T = T
rule_method_error_extract_names(::Nothing)                  = ()

rule_method_error_extract_types(t::Tuple)   = map(e -> typeof(getdata(e)), t)
rule_method_error_extract_types(t::Nothing) = ()

rule_method_error_extract_meta(something) = string("meta::", typeof(something))
rule_method_error_extract_meta(::Nothing) = ""

struct RuleMethodError
    fform
    on
    vconstraint
    mnames
    messages
    qnames
    marginals
    meta
    node
end

rule(fform, on, vconstraint, mnames, messages, qnames, marginals, meta, __node) = throw(RuleMethodError(fform, on, vconstraint, mnames, messages, qnames, marginals, meta, __node))

function Base.showerror(io::IO, error::RuleMethodError)

    print(io, "RuleMethodError: no method matching rule for the given arguments")

    try
        node = error.node !== nothing ? error.node : NodeErrorStub()

        spec_fform       = rule_method_error_extract_fform(error.fform)
        spec_on          = rule_method_error_extract_on(error.on)
        spec_vconstraint = rule_method_error_extract_vconstraint(error.vconstraint)

        m_names   = rule_method_error_extract_names(error.mnames)
        m_indices = map(n -> interfaceindex(node, n), m_names)
        q_names   = rule_method_error_extract_names(error.qnames)
        q_indices = map(n -> interfaceindex(node, n), q_names)

        spec_m_names = map(e -> string("m_", e), m_names)
        spec_m_types = rule_method_error_extract_types(error.messages)
        spec_q_names = map(e -> string("q_", e), q_names)
        spec_q_types = rule_method_error_extract_types(error.marginals)

        spec_m = map(m -> string(m[1], "::", m[2]), zip(spec_m_names, spec_m_types))
        spec_q = map(q -> string(q[1], "::", q[2]), zip(spec_q_names, spec_q_types))

        spec = Vector(undef, 2length(interfaces(node)))

        fill!(spec, nothing)

        for (i, j) in enumerate(m_indices)
            spec[2j - 1] = spec_m[i]
        end

        for (i, j) in enumerate(q_indices)
            spec[2j] = spec_q[i]
        end

        filter!(!isnothing, spec)

        arguments_spec = join(spec, ", ")
        meta_spec      = rule_method_error_extract_meta(error.meta)

        possible_fix_definition = """
        @rule $(spec_fform)(:$spec_on, $spec_vconstraint) ($arguments_spec, $meta_spec) = begin
            return ...
        end
        """

        println(io, "\n\nPossible fix, define:\n")
        println(io, possible_fix_definition)
    catch uerror
        println(io, "Unexpected error occured during Base.showerror(::RuleMethodError)")
        println(io, uerror)
        println(io, "rule.fform: ", error.fform)
        println(io, "rule.on: ", error.on)
        println(io, "rule.vconstraint: ", error.vconstraint)
        println(io, "rule.mnames: ", error.mnames)
        println(io, "rule.messages: ", error.messages)
        println(io, "rule.qnames: ", error.qnames)
        println(io, "rule.marginals: ", error.ffmarginalsrm)
        println(io, "rule.meta: ", error.meta)
        println(io, "rule.node: ", error.node)
    end
end

struct MarginalRuleMethodError
    fform
    on
    mnames
    messages
    qnames
    marginals
    meta
    node
end

marginalrule(fform, on, mnames, messages, qnames, marginals, meta, __node) = throw(MarginalRuleMethodError(fform, on, mnames, messages, qnames, marginals, meta, __node))

function Base.showerror(io::IO, error::MarginalRuleMethodError)

    print(io, "MarginalRuleMethodError: no method matching rule for the given arguments")

    try
        node = error.node !== nothing ? error.node : NodeErrorStub()

        spec_fform       = rule_method_error_extract_fform(error.fform)
        spec_on          = rule_method_error_extract_on(error.on)

        m_names   = rule_method_error_extract_names(error.mnames)
        m_indices = map(n -> interfaceindex(node, n), m_names)
        q_names   = rule_method_error_extract_names(error.qnames)
        q_indices = map(n -> interfaceindex(node, n), q_names)

        spec_m_names = map(e -> string("m_", e), m_names)
        spec_m_types = rule_method_error_extract_types(error.messages)
        spec_q_names = map(e -> string("q_", e), q_names)
        spec_q_types = rule_method_error_extract_types(error.marginals)

        spec_m = map(m -> string(m[1], "::", m[2]), zip(spec_m_names, spec_m_types))
        spec_q = map(q -> string(q[1], "::", q[2]), zip(spec_q_names, spec_q_types))

        spec = Vector(undef, 2length(interfaces(node)))

        fill!(spec, nothing)

        for (i, j) in enumerate(m_indices)
            spec[2j - 1] = spec_m[i]
        end

        for (i, j) in enumerate(q_indices)
            spec[2j] = spec_q[i]
        end

        filter!(!isnothing, spec)

        arguments_spec = join(spec, ", ")
        meta_spec      = rule_method_error_extract_meta(error.meta)

        possible_fix_definition = """
        @marginalrule $(spec_fform)(:$spec_on) ($arguments_spec, $meta_spec) = begin
            return ...
        end
        """

        println(io, "\n\nPossible fix, define:\n")
        println(io, possible_fix_definition)
    catch uerror
        println(io, "Unexpected error occured during Base.showerror(::MarginalRuleMethodError)")
        println(io, uerror)
        println(io, "rule.fform: ", error.fform)
        println(io, "rule.on: ", error.on)
        println(io, "rule.mnames: ", error.mnames)
        println(io, "rule.messages: ", error.messages)
        println(io, "rule.qnames: ", error.qnames)
        println(io, "rule.marginals: ", error.ffmarginalsrm)
        println(io, "rule.meta: ", error.meta)
        println(io, "rule.node: ", error.node)
    end
end


text:
using Test
using BenchmarkTools

##

model_pendulum, obj_uncon_pendulum = TrajectoryOptimization.Dynamics.pendulum!

u_min_pendulum = -2
u_max_pendulum = 2
x_min_pendulum = [-20;-20]
x_max_pendulum = [20; 20]

model_dubins, obj_uncon_dubins = TrajectoryOptimization.Dynamics.dubinscar!

# dubins car
u_min_dubins = [-1; -1]
u_max_dubins = [1; 1]
x_min_dubins = [0; -100; -100]
x_max_dubins = [1.0; 100; 100]
obj_con_dubins = ConstrainedObjective(obj_uncon_dubins, u_min=u_min_dubins, u_max=u_max_dubins, x_min=x_min_dubins, x_max=x_max_dubins)

# -Constrained objective
obj_con_pendulum = ConstrainedObjective(obj_uncon_pendulum, u_min=u_min_pendulum, u_max=u_max_pendulum, x_min=x_min_pendulum, x_max=x_max_pendulum)

model = model_dubins
obj = obj_con_dubins
u_max = u_max_dubins
u_min = u_min_dubins

obj = TrajectoryOptimization.update_objective(obj, tf=:min, c=0.0, Q = 1e-3*Diagonal(I,model.n), R = 1e-3*Diagonal(I,model.m), Qf = Diagonal(I,model.n))

# Solver
intergrator_foh = :rk3_foh

dt = 0.005
solver1 = Solver(model,obj,integration=intergrator_foh,N=51)
solver1.opts.minimum_time = true
solver1.opts.infeasible = true
solver1.opts.constrained = true
X0 = line_trajectory(solver1)
U0 = ones(solver1.model.m,solver1.N)
# U0 = [U0; ones(1,solver1.N)]
# u0 = infeasible_controls(solver1,X0)
# U0 = [U0;u0]

solver2 = Solver(model,obj,integration=intergrator_foh,N=51)
solver2.opts.minimum_time = true
solver2.opts.infeasible = true
solver2.opts.constrained = true

results1 = init_results(solver1,X0,U0)
results2 = init_results(solver2,X0,U0)
rollout!(results1,solver1)
rollout!(results2,solver2)
update_jacobians!(results1, solver1)
update_jacobians!(results2, solver2)

println("TEST")
@time v1 = _backwardpass_foh!(results1,solver1)
@time v2 = _backwardpass_foh_speedup!(results2,solver2)

@test isapprox(v1[1:2], v2[1:2])
@test isapprox(to_array(results1.K),to_array(results2.K))
@test isapprox(to_array(results1.b),to_array(results2.b))
@test isapprox(to_array(results1.d),to_array(results2.d))

@benchmark _backwardpass_foh!($results1,$solver1)
@benchmark _backwardpass_foh_speedup!($results2,$solver2)
println("\n")

# Benchmark
# 1. eliminate redundant L calculations
# 189.564 ms (614530 allocations: 35.21 MiB)
# 161.772 ms (576530 allocations: 31.78 MiB)
# 2. distribute dt calculations
# 148.805 ms (614530 allocations: 35.21 MiB)
# 144.791 ms (570530 allocations: 31.57 MiB)
# 3. view into s,S
# 153.813 ms (614530 allocations: 35.21 MiB)
# 134.082 ms (548530 allocations: 30.94 MiB)
# 4. eliminate vec()
# 161.563 ms (614530 allocations: 35.21 MiB)
# 142.320 ms (547522 allocations: 30.92 MiB)
# 5. distribute discrete dynamics calculation
# 156.674 ms (614530 allocations: 35.21 MiB)
# 119.566 ms (491524 allocations: 27.89 MiB)
# 6. don't index into U matrix
# 159.905 ms (614530 allocations: 35.21 MiB)
# 122.204 ms (490524 allocations: 27.80 MiB)
# 7. view into U matrix and calculate transposes only once
# 159.132 ms (614530 allocations: 35.21 MiB)
# 104.974 ms (456524 allocations: 26.30 MiB)
# 8. replace more transposes
# 156.551 ms (614530 allocations: 35.21 MiB)
# 89.473 ms (424524 allocations: 25.34 MiB)
# 9. linear indexing -note, I found that using lin. idx on views was slower
# 160.612 ms (614530 allocations: 35.21 MiB)
# 89.392 ms (418534 allocations: 25.10 MiB)
# 10. more intermediate calculation replacements
# 172.848 ms (614530 allocations: 35.21 MiB)
# 88.955 ms (413533 allocations: 24.56 MiB)

# RK3, RK3 FOH, RK3 FOH speedup
# .0069ms, .0193ms, .0148ms
# 23% speedup
# previously (on Brian's computer) .00413, .0121 -> 2.929x slower
# now 2.14x slower

# 12-3-2018
# 1.
# 466.429 ms (511527 allocations: 38.27 MiB)
# 289.674 ms (338527 allocations: 26.21 MiB)
# 2. Use pre-allocated S, s
# 295.870 ms (511527 allocations: 38.27 MiB)
# 194.656 ms (338523 allocations: 26.21 MiB)
# 3. Using views into all L,l is slow and more allocations
# --- Using minimum_time and infeasible
# 4.
# 22.327 ms (43878 allocations: 4.40 MiB)
# 17.161 ms (35044 allocations: 3.71 MiB)
# 5. - reduced number of calculations (also found small error in Lu)
# 21.847 ms (43878 allocations: 4.40 MiB)
# 18.010 ms (33769 allocations: 3.55 MiB)
# 6. - more calculation reductions
# 22.863 ms (43878 allocations: 4.40 MiB)
# 17.340 ms (33361 allocations: 3.51 MiB)
# 7. -views into preallocated L,l for minimum_time
# 22.817 ms (43878 allocations: 4.40 MiB)
# 16.899 ms (29273 allocations: 3.21 MiB)
# 8. -view into preallocated L,l for infeasible start
# 22.136 ms (43878 allocations: 4.40 MiB)
# 13.648 ms (25710 allocations: 3.01 MiB)
# 9. -removed duplicate assignments
# 22.393 ms (43878 allocations: 4.40 MiB)
# 13.492 ms (24996 allocations: 2.99 MiB)
# 10. -smarter indexing
# 39.727 ms (43878 allocations: 4.40 MiB)
# 19.778 ms (24639 allocations: 2.97 MiB)

model, obj = Dynamics.cartpole_analytical
n,m = model.n, model.m
N = 51
dt = 0.1

obj.x0 = [0;0;0;0]
obj.xf = [0.5;pi;0;0]
obj.tf = 2.0
U0 = ones(m,N)
solver_foh = Solver(model,obj,N=N,opts=opts,integration=:rk3_foh)
solver_zoh = Solver(model,obj,N=N,opts=opts,integration=:rk3)

k = 10
time_per_iter_foh = zeros(k)
time_per_iter_zoh = zeros(k)

for i = 1:k
  res_foh, stat_foh = solve(solver_foh,U0)
  res_zoh, stat_zoh = solve(solver_zoh,U0)
  time_per_iter_foh[i] = stat_foh["runtime"]/stat_foh["iterations"]
  time_per_iter_zoh[i] = stat_zoh["runtime"]/stat_zoh["iterations"]
end
println("Time per iter (foh): $(mean(time_per_iter_foh))")
println("Time per iter (zoh): $(mean(time_per_iter_zoh))")

ext: jl


text:
#=

This file specifies how to pretty-print vSmartMOM_Parameters module types

=#

function Base.show(io::IO,::MIME"text/plain", x::vSmartMOM_Parameters)

    println(io, "\n------------------")
    println(io, "Radiative Transfer")
    println(io, "------------------")
    println(io, "\tSpectral Bands (ν₁:n:νₙ, cm⁻¹):")
    for band in x.spec_bands
        println(io, "\t\t- $(length(band))-length vector from $(round(band[1], digits=2)) to $(round(band[end], digits=2))")
    end
    println(io, "\tSurface BRDFs: ")
    for surface in x.brdf
        println(io, "\t\t- $(surface)")
    end
    println(io, "\tQuadrature Type: $(x.quadrature_type)")
    println(io, "\tPolarization Type: $(x.polarization_type)")
    println(io, "\tmax_m: $(x.max_m)")
    println(io, "\tΔ_angle: $(x.Δ_angle)°")
    println(io, "\tl_trunc: $(x.l_trunc)")
    println(io, "\tDepolarization: $(x.depol)")
    println(io, "\tFloating-point type: $(x.float_type)")
    println(io, "\tArchitecture: $(x.architecture)")

    println(io, "\n--------")
    println(io, "Geometry")
    println(io, "--------")
    println(io, "\tSZA (deg): $(x.sza)")
    println(io, "\tVZA (deg): $(x.vza)")
    println(io, "\tVAZ (deg): $(x.vaz)")
    println(io, "\tObservation Altitude: $(x.obs_alt)")

    println(io, "\n-------------------")
    println(io, "Atmospheric Profile")
    println(io, "-------------------")
    println(io, "\tTemperature: $(length(x.T))-length array")
    println(io, "\tPressure: $(length(x.p))-length array")
    println(io, "\tSpecific humidity: $(length(x.q))-length array")
    if x.profile_reduction_n == -1
        println(io, "\tProfile Reduction: None")
    else
        println(io, "\tProfile Reduction: $(x.profile_reduction_n) layers")
    end

    println(io, "\n----------")
    println(io, "Absorption")
    println(io, "----------")
    if (!isnothing(x.absorption_params))
        println(io, "\tMolecules")
        for molecule_band in x.absorption_params.molecules
            println(io, "\t\t- $(molecule_band)")
        end
        println(io, "\tVMR:")
        for mol in keys(x.absorption_params.vmr)
            println(io, "\t\t- $(mol): $(x.absorption_params.vmr[mol])")
        end
        println(io, "\tBroadening: $(x.absorption_params.broadening_function)")
        println(io, "\tCEF: $(x.absorption_params.CEF)")
        println(io, "\tWing Cutoff: $(x.absorption_params.wing_cutoff) cm⁻¹")
    else
        println(io, "(No Absorption Parameters Specified)")
    end

    println(io, "\n----------")
    println(io, "Scattering")
    println(io, "----------")
    if (!isnothing(x.scattering_params))
        println(io, "\tAerosols:")
        for i in 1:length(x.scattering_params.rt_aerosols)
            println(io, "\t\t- Aerosol #$(i)")
            curr_rt_aerosol = x.scattering_params.rt_aerosols[i]
            println(io, "\t\t  τ_ref: $(curr_rt_aerosol.τ_ref)")
            println(io, "\t\t  Size Distribution: $(curr_rt_aerosol.aerosol.size_distribution)")
            #println(io, "\t\t  σ: $(curr_rt_aerosol.aerosol.size_distribution.σ) μm; geometric standard deviation")
            println(io, "\t\t  nᵣ: $(curr_rt_aerosol.aerosol.nᵣ)")
            println(io, "\t\t  nᵢ: $(curr_rt_aerosol.aerosol.nᵢ)")
            println(io, "\t\t  p₀: $(curr_rt_aerosol.p₀) Pa")
            println(io, "\t\t  σp: $(curr_rt_aerosol.σp) Pa")
        end
        println(io, "\tr_max: $(x.scattering_params.r_max) μm")
        println(io, "\tnquad_radius: $(x.scattering_params.nquad_radius)")
        println(io, "\tλ_ref: $(x.scattering_params.λ_ref) μm")
        print(io, "\tDecomposition Type: $(x.scattering_params.decomp_type)")
    else
        print(io, "(No Scattering Parameters Specified)")
    end

end

# Overload the show method for vSmartMOM_Model
function Base.show(io::IO, x::vSmartMOM_Model)

    println(io, "        params| ::vSmartMOM_Parameters")
    println(io, "aerosol_optics| [i_band][i_aer] ($(length(x.aerosol_optics)) x $(length(x.aerosol_optics[1])))")
    println(io, "greek_rayleigh| ::GreekCoefs")
    println(io, "   quad_points| ::QuadPoints")
    println(io, "         τ_abs| [i_band][i_spec,i_z] ($(length(x.τ_abs)) x (n_spec x $(size(x.τ_abs[1])[2])))")
    println(io, "        τ_rayl| [i_band][i_z] ($(length(x.τ_rayl)) x $(length(x.τ_rayl[1])))")
    println(io, "         τ_aer| [i_band][i_aer,i_z] ($(length(x.τ_aer)) x ($(size(x.τ_aer[1])[1]) x $(size(x.τ_aer[1])[2])))")
    println(io, "      obs_geom| ::ObsGeometry")
    println(io, "       profile| ::AtmosphericProfile")
end


text:
using Test
using {{{PKG}}}

@testset "{{{PKG}}}.jl" begin
    for (root, _, files) in walkdir("../src")
        for file in files
            if endswith(file, "_test.jl")
                include(joinpath(root, file))
            end
        end
    end
end


text:
# #Inflación de media simple
# InflationSimpleMean.jl - Función de inflación de media simple interanual

"""
    InflationSimpleMean <: InflationFunction
Función de inflación para computar la media simple interanual.
"""
struct InflationSimpleMean <: InflationFunction
end

# 2. Extender el método de nombre
measure_name(::InflationSimpleMean) = "Media simple interanual"

# Define cómo opera InflationSimpleMean sobre un objeto de tipo VarCPIBase.
function (inflfn::InflationSimpleMean)(base::VarCPIBase{T}) where T
    #Obtener el indice correspondiente a las variaciones intermensuales
    indmat = capitalize(base.v)
    #Cálculo de la variación interanual
    y2ymat = varinteran(indmat)
    #Cálculo de la media simple de las variaciones interanuales
    y2ysimpmean = mean(y2ymat,dims=2)
    #Utilizar el índice de la media simple de las variaciones intermensuales de los primeros
    #11 meses para completar la serie del índice de media simple de variaciones interanuales
    m2msimpmean = mean(base.v,dims=2)
    simpmeanind = zeros(T,periods(base)+1); simpmeanind[1]=100;
    simpmeanind[2:12] = capitalize(m2msimpmean[1:11])
    for j = 13:periods(base)+1
        simpmeanind[j] = (y2ysimpmean[j-12]/100 + 1)*simpmeanind[j-12]
    end
    #Cálculo de la variación intermensual del índice de media simple
    simpmean = varinterm(simpmeanind[2:end])
    simpmean
end


text:
using EMIRT
using Watershed
using Agglomeration
using Process
using ImageView

#include(joinpath(dirname(@__FILE__), "../plugins/show.jl"))
include(joinpath(Pkg.dir(), "EMIRT/plugins/show.jl"))

# define
faffs = Dict{Symbol, AbstractString}(
        :jnet       => "/usr/people/jingpeng/seungmount/research/kisuklee/Sharing/Jingpeng/blend_test/JNet/chunk_33405_8905_229.aff.h5",
        :ensemble   => "/usr/people/jingpeng/seungmount/research/kisuklee/Sharing/Jingpeng/blend_test/Ensemble/chunk_33405_8905_229.aff.h5",
        :multiscale => "/usr/people/jingpeng/seungmount/research/kisuklee/Sharing/Jingpeng/blend_test/MS/chunk_33405_8905_229.aff.h5",
        :msf        => "/usr/people/jingpeng/seungmount/research/kisuklee/Sharing/Jingpeng/blend_test/MSF/chunk_33405_8905_229.aff.h5" )

flbl = "/usr/people/jingpeng/seungmount/Omni/TracerTasks/ZFishEnsembleValidation/zfish_chunk_33405_8905_229.omni.seg.ben.h5"
lbl  = readseg(flbl)

ecs = ScoreCurves()
for (name,faff) in faffs
  aff = readaff(faff)
  seg = atomicseg(aff; is_threshold_relative=true)
  segmentPairs, segmentPairAffinities = Process.forward(aff, seg)
  sgm = EMIRT.SegMST(seg, segmentPairs, segmentPairAffinities)

  errorcurve = sgm2ec(sgm, lbl, 0:0.2:1)
  append!(ecs, errorcurve; tag=name)
end

save("/tmp/ecs.zfish.h5", ecs);
# transform to dataframes
df = ecs2df(ecs)

# %% plot
plot(df, x="thd", y="rf", Geom.line,
               Guide.xlabel("threshold"),
               Guide.ylabel("rand f score"))

plot(df, x="thd", y="re", Geom.line,
              Guide.xlabel("threshold"),
              Guide.ylabel("rand error"))

# %%
# run(`julia overlay_img_seg.jl $fimg $flbl`)


text:
function bitmask_type(spec::SpecBitmask)
    @match spec.width begin
        8 => :UInt8
        16 => :UInt16
        32 => :UInt32
        64 => :UInt64
        _ => error("Failed to generate bitmask type with width $(spec.width) from $spec")
    end
end

function wrap(spec::SpecBitmask)
    name = bitmask_flag_type(spec.name)
    Dict(
        :category => :enum,
        :macro => Symbol("@bitmask_flag"),
        :values => map(x -> :($(remove_vk_prefix(x.name)) = $(2^(x.position))), spec.bits),
        :decl => :($name::$(bitmask_type(spec))),
    )
end


text:
using OperatorFlux
using Test
using FFTW
using Flux
using LinearAlgebra
using ChainRulesCore
using ChainRulesTestUtils

@testset "ChainRulesCore Extensions" begin
    # FFTW.r2r
    M = 4
    trafo = ChebyshevTransform(modes = (4,))
    c = rand(3, M, 5)
    test_rrule(
        FFTW.r2r,
        c,
        FFTW.REDFT00,
        2:2,
        check_thunked_output_tangent = false,
    )

    # truncate modes
    M = 4
    trafo = FourierTransform(modes = (2, 2))
    c = rand(3, M, M, 5)
    test_rrule(
        OperatorFlux.truncate_modes,
        trafo,
        c,
        check_thunked_output_tangent = false,
    )

    trafo = ChebyshevTransform(modes = (2, 2))
    test_rrule(
        OperatorFlux.truncate_modes,
        trafo,
        c,
        check_thunked_output_tangent = false,
    )

    for mode in [(5,), (5, 6), (5, 7, 9)]
        z = rand(3, mode..., 5)
        trafos = (LegendreTransform(modes = mode),)
        for trafo in trafos
            test_rrule(
                OperatorFlux.truncate_modes,
                trafo,
                z,
                check_thunked_output_tangent = false,
            )
        end
    end

    # pad modes
    c = rand(3, 2, 8, 5)
    trafo = FourierTransform(modes = (2, 4))
    test_rrule(
        OperatorFlux.pad_modes,
        trafo,
        c,
        (4, 8),
        check_thunked_output_tangent = false,
    )

    trafo = ChebyshevTransform(modes = (2, 4))
    c = rand(3, 2, 4, 5)
    test_rrule(
        OperatorFlux.pad_modes,
        trafo,
        c,
        (2, 4),
        check_thunked_output_tangent = false,
    )

    for mode in [(5,), (5, 6), (5, 7, 9)]
        z = rand(3, mode..., 5)
        trafos = (LegendreTransform(modes = mode),)
        for trafo in trafos
            test_rrule(
                OperatorFlux.pad_modes,
                trafo,
                z,
                mode,
                check_thunked_output_tangent = false,
            )
        end
    end
end


text:
struct DensePathway{I,J,S} <: Pathway
    pre::I      # presynaptic units
    syns::S     # synaptic state, of length pre×post
    post::J     # postsynaptic units
end

"""
Apply `f` on every connection in the pathway where `select` is true
for the presynaptic element.
"""
function Base.broadcast(f, path::DensePathway, select=(_)->true)
    @inbounds for i in eachindex(path.pre)
        if select(Elem(path.pre,i))
            @simd for j in eachindex(path.post)
                f(Elem(path.pre,i), Elem(path.syns,j,i), Elem(path.post,j))
            end
        end
    end
    nothing
end


text:
@testset "Linear algebra" begin

    @testset "SVector as a (mathematical) vector space" begin
        c = 2
        v1 = @SVector [2,4,6,8]
        v2 = @SVector [4,3,2,1]

        @test v1 + c === @SVector [4,6,8,10]
        @test v1 - c === @SVector [0,2,4,6]
        @test v1 * c === @SVector [4,8,12,16]
        @test v1 / c === @SVector [1.0,2.0,3.0,4.0]

        @test v1 + v2 === @SVector [6, 7, 8, 9]
        @test v1 - v2 === @SVector [-2, 1, 4, 7]
    end

    @testset "eye()" begin
        @test eye(SMatrix{2,2,Int}) === @SMatrix [1 0; 0 1]
        @test eye(SMatrix{2,2}) === @SMatrix [1.0 0.0; 0.0 1.0]
        @test eye(SMatrix{2}) === @SMatrix [1.0 0.0; 0.0 1.0]

        @test eye(MMatrix{2,2,Int})::MMatrix == @MMatrix [1 0; 0 1]
        @test eye(MMatrix{2,2})::MMatrix == @MMatrix [1.0 0.0; 0.0 1.0]
        @test eye(MMatrix{2})::MMatrix == @MMatrix [1.0 0.0; 0.0 1.0]
    end

    @testset "cross()" begin
        @test cross(SVector(1,2,3), SVector(4,5,6)) === SVector(-3, 6, -3)
    end

    @testset "transpose() and conj()" begin
        @test conj(SVector(1+im, 2+im)) === SVector(1-im, 2-im)

        @test @SVector([1, 2, 3]).' === @SMatrix([1 2 3])
        @test @SMatrix([1 2; 0 3]).' === @SMatrix([1 0; 2 3])
        @test @SMatrix([1 2 3; 4 5 6]).' === @SMatrix([1 4; 2 5; 3 6])

        @test @SVector([1, 2, 3])' === @SMatrix([1 2 3])
        @test @SMatrix([1 2; 0 3])' === @SMatrix([1 0; 2 3])
        @test @SMatrix([1 2 3; 4 5 6])' === @SMatrix([1 4; 2 5; 3 6])
    end

    @testset "vcat() and hcat()" begin
        @test vcat(SVector(1,2,3), SVector(4,5,6)) === SVector(1,2,3,4,5,6)
        @test hcat(SVector(1,2,3), SVector(4,5,6)) === @SMatrix [1 4; 2 5; 3 6]

        @test vcat(@SMatrix([1;2;3]), SVector(4,5,6)) === @SMatrix([1;2;3;4;5;6])
        @test vcat(SVector(1,2,3), @SMatrix([4;5;6])) === @SMatrix([1;2;3;4;5;6])
        @test hcat(@SMatrix([1;2;3]), SVector(4,5,6)) === @SMatrix [1 4; 2 5; 3 6]
        @test hcat(SVector(1,2,3), @SMatrix([4;5;6])) === @SMatrix [1 4; 2 5; 3 6]

        @test vcat(@SMatrix([1;2;3]), @SMatrix([4;5;6])) === @SMatrix([1;2;3;4;5;6])
        @test hcat(@SMatrix([1;2;3]), @SMatrix([4;5;6])) === @SMatrix [1 4; 2 5; 3 6]
    end
end


text:
module ThreadedSparseArrays

export ThreadedSparseMatrixCSC

using LinearAlgebra
import LinearAlgebra: mul!
using SparseArrays
import SparseArrays: getcolptr, AbstractSparseMatrixCSC
const AdjOrTransDenseMatrix = if VERSION < v"1.6.0-rc2"
    SparseArrays.AdjOrTransStridedOrTriangularMatrix
else
    SparseArrays.AdjOrTransDenseMatrix
end

# * Threading utilities
struct RangeIterator
    k::Int
    d::Int
    r::Int
end

"""
    RangeIterator(n::Int,k::Int)

Returns an iterator splitting the range `1:n` into `min(k,n)` parts of (almost) equal size.
"""
RangeIterator(n::Int, k::Int) = RangeIterator(min(n,k),divrem(n,k)...)
Base.length(it::RangeIterator) = it.k
endpos(it::RangeIterator, i::Int) = i*it.d+min(i,it.r)
Base.iterate(it::RangeIterator, i::Int=1) = i>it.k ? nothing : (endpos(it,i-1)+1:endpos(it,i), i+1)

# * ThreadedSparseMatrixCSC

"""
    ThreadedSparseMatrixCSC(A)

Thin container around `A::SparseMatrixCSC` that will enable certain
threaded multiplications of `A` with dense matrices.
"""
struct ThreadedSparseMatrixCSC{Tv,Ti,At} <: AbstractSparseMatrixCSC{Tv,Ti}
    A::At
    ThreadedSparseMatrixCSC(A::At) where {Tv,Ti,At<:AbstractSparseMatrixCSC{Tv,Ti}} =
        new{Tv,Ti,At}(A)
end

Base.size(A::ThreadedSparseMatrixCSC, args...) = size(A.A, args...)

for f in [:rowvals, :nonzeros, :getcolptr]
    @eval SparseArrays.$(f)(A::ThreadedSparseMatrixCSC) = SparseArrays.$(f)(A.A)
end

function mul!(C::StridedVecOrMat, A::ThreadedSparseMatrixCSC, B::Union{StridedVector,AdjOrTransDenseMatrix}, α::Number, β::Number)
    size(A, 2) == size(B, 1) || throw(DimensionMismatch())
    size(A, 1) == size(C, 1) || throw(DimensionMismatch())
    size(B, 2) == size(C, 2) || throw(DimensionMismatch())
    nzv = nonzeros(A)
    rv = rowvals(A)
    if β != 1
        β != 0 ? rmul!(C, β) : fill!(C, zero(eltype(C)))
    end
    @sync for r in RangeIterator(size(C,2), Threads.nthreads())
        Threads.@spawn for k in r
            @inbounds for col = 1:size(A, 2)
                αxj = B[col,k] * α
                for j = getcolptr(A)[col]:(getcolptr(A)[col + 1] - 1)
                    C[rv[j], k] += nzv[j]*αxj
                end
            end
        end
    end
    C
end

function mul!(C::StridedVecOrMat, adjA::Adjoint{<:Any,<:ThreadedSparseMatrixCSC}, B::AdjOrTransDenseMatrix, α::Number, β::Number)
    A = adjA.parent
    size(A, 2) == size(C, 1) || throw(DimensionMismatch())
    size(A, 1) == size(B, 1) || throw(DimensionMismatch())
    size(B, 2) == size(C, 2) || throw(DimensionMismatch())
    colptrA = getcolptr(A)
    nzv = nonzeros(A)
    rv = rowvals(A)
    if β != 1
        β != 0 ? rmul!(C, β) : fill!(C, zero(eltype(C)))
    end
    @sync for r in RangeIterator(size(C,2), Threads.nthreads())
        Threads.@spawn for k in r
            @inbounds for col = 1:size(A, 2)
                tmp = zero(eltype(C))
                for j = getcolptr(A)[col]:(getcolptr(A)[col + 1] - 1)
                    tmp += adjoint(nzv[j])*B[rv[j],k]
                end
                C[col,k] += tmp * α
            end
        end
    end
    C
end
function mul!(C::StridedVecOrMat, adjA::Adjoint{<:Any,<:ThreadedSparseMatrixCSC}, B::StridedVector, α::Number, β::Number)
    A = adjA.parent
    size(A, 2) == size(C, 1) || throw(DimensionMismatch())
    size(A, 1) == size(B, 1) || throw(DimensionMismatch())
    size(B, 2) == size(C, 2) || throw(DimensionMismatch())
    @assert size(B,2)==1
    colptrA = getcolptr(A)
    nzv = nonzeros(A)
    rv = rowvals(A)
    if β != 1
        β != 0 ? rmul!(C, β) : fill!(C, zero(eltype(C)))
    end
    @sync for r in RangeIterator(size(A,2), Threads.nthreads())
        Threads.@spawn @inbounds for col = r
            tmp = zero(eltype(C))
            for j = getcolptr(A)[col]:(getcolptr(A)[col + 1] - 1)
                tmp += adjoint(nzv[j])*B[rv[j]]
            end
            C[col] += tmp * α
        end
    end
    C
end

function mul!(C::StridedVecOrMat, transA::Transpose{<:Any,<:ThreadedSparseMatrixCSC}, B::AdjOrTransDenseMatrix, α::Number, β::Number)
    A = transA.parent
    size(A, 2) == size(C, 1) || throw(DimensionMismatch())
    size(A, 1) == size(B, 1) || throw(DimensionMismatch())
    size(B, 2) == size(C, 2) || throw(DimensionMismatch())
    nzv = nonzeros(A)
    rv = rowvals(A)
    if β != 1
        β != 0 ? rmul!(C, β) : fill!(C, zero(eltype(C)))
    end
    @sync for r in RangeIterator(size(C,2), Threads.nthreads())
        Threads.@spawn for k in r
            @inbounds for col = 1:size(A, 2)
                tmp = zero(eltype(C))
                for j = getcolptr(A)[col]:(getcolptr(A)[col + 1] - 1)
                    tmp += transpose(nzv[j])*B[rv[j],k]
                end
                C[col,k] += tmp * α
            end
        end
    end
    C
end
function mul!(C::StridedVecOrMat, transA::Transpose{<:Any,<:ThreadedSparseMatrixCSC}, B::StridedVector, α::Number, β::Number)
    A = transA.parent
    size(A, 2) == size(C, 1) || throw(DimensionMismatch())
    size(A, 1) == size(B, 1) || throw(DimensionMismatch())
    size(B, 2) == size(C, 2) || throw(DimensionMismatch())
    @assert size(B,2)==1
    nzv = nonzeros(A)
    rv = rowvals(A)
    if β != 1
        β != 0 ? rmul!(C, β) : fill!(C, zero(eltype(C)))
    end
    @sync for r in RangeIterator(size(A,2), Threads.nthreads())
        Threads.@spawn @inbounds for col = r
            tmp = zero(eltype(C))
            for j = getcolptr(A)[col]:(getcolptr(A)[col + 1] - 1)
                tmp += transpose(nzv[j])*B[rv[j]]
            end
            C[col] += tmp * α
        end
    end
    C
end

function mul!(C::StridedVecOrMat, X::AdjOrTransDenseMatrix, A::ThreadedSparseMatrixCSC, α::Number, β::Number)
    mX, nX = size(X)
    nX == size(A, 1) || throw(DimensionMismatch())
    mX == size(C, 1) || throw(DimensionMismatch())
    size(A, 2) == size(C, 2) || throw(DimensionMismatch())
    rv = rowvals(A)
    nzv = nonzeros(A)
    if β != 1
        β != 0 ? rmul!(C, β) : fill!(C, zero(eltype(C)))
    end
    @sync for r in RangeIterator(size(A,2), Threads.nthreads())
        Threads.@spawn for col in r
            @inbounds for k=getcolptr(A)[col]:(getcolptr(A)[col+1]-1)
                j = rv[k]
                αv = nzv[k]*α
                for multivec_row=1:mX
                    C[multivec_row, col] += X[multivec_row, j] * αv
                end
            end
        end
    end
    C
end

end # module


text:
module GraphSearch

export graph_search!, AbstractTreeSearchAlgorithm

include("tree_search_node.jl")
include("tree_search_algo.jl")

end


text:
## core string functions ##

length(s::String) = error("you must implement length(",typeof(s),")")
next(s::String, i::Int) = error("you must implement next(",typeof(s),",Int)")
next(s::DirectIndexString, i::Int) = (s[i],i+1)
next(s::String, i::Integer) = next(s,int(i))

## generic supplied functions ##

start(s::String) = 1
done(s::String,i) = (i > length(s))
isempty(s::String) = done(s,start(s))
ref(s::String, i::Int) = next(s,i)[1]
ref(s::String, i::Integer) = s[int(i)]
ref(s::String, x::Real) = s[iround(x)]
ref{T<:Integer}(s::String, r::Range1{T}) = s[int(first(r)):int(last(r))]

symbol(s::String) = symbol(cstring(s))
string(s::String) = s

print(s::String) = for c=s; print(c); end
print(x...) = for i=x; print(i); end
println(args...) = print(args..., '\n')

show(s::String) = print_quoted(s)

(*)(s::String...) = strcat(s...)
(^)(s::String, r::Integer) = repeat(s,r)

size(s::String) = (length(s),)
size(s::String, d::Integer) = d==1 ? length(s) :
    error("in size: dimension ",d," out of range")

strlen(s::DirectIndexString) = length(s)
function strlen(s::String)
    i = start(s)
    if done(s,i)
        return 0
    end
    n = 1
    while true
        c, j = next(s,i)
        if done(s,j)
            return n
        end
        n += 1
        i = j
    end
end

isvalid(s::DirectIndexString, i::Integer) = (start(s) <= i <= length(s))
function isvalid(s::String, i::Integer)
    try
        next(s,i)
        true
    catch
        false
    end
end

prevind(s::DirectIndexString, i::Integer) = i-1
thisind(s::DirectIndexString, i::Integer) = i
nextind(s::DirectIndexString, i::Integer) = i+1

prevind(s::String, i::Integer) = thisind(s,thisind(s,i)-1)

function thisind(s::String, i::Integer)
    for j = i:-1:1
        if isvalid(s,j)
            return j
        end
    end
    return 0 # out of range
end

function nextind(s::String, i::Integer)
    for j = i+1:length(s)
        if isvalid(s,j)
            return j
        end
    end
    length(s)+1 # out of range
end

ind2chr(s::DirectIndexString, i::Integer) = i
chr2ind(s::DirectIndexString, i::Integer) = i

function ind2chr(s::String, i::Integer)
    s[i] # throws error if invalid
    j = 1
    k = start(s)
    while true
        c, l = next(s,k)
        if i <= k
            return j
        end
        j += 1
        k = l
    end
end

function chr2ind(s::String, i::Integer)
    if i < 1
        return i
    end
    j = 1
    k = start(s)
    while true
        c, l = next(s,k)
        if i == j
            return k
        end
        j += 1
        k = l
    end
end

function strchr(s::String, c::Char, i::Integer)
    i = nextind(s,i)
    while !done(s,i)
        d, j = next(s,i)
        if c == d
            return i
        end
        i = j
    end
    return 0
end
strchr(s::String, c::Char) = strchr(s, c, start(s))
contains(s::String, c::Char) = (strchr(s,c)!=0)

function chars(s::String)
    cx = Array(Char,strlen(s))
    i = 0
    for c in s
        cx[i += 1] = c
    end
    return cx
end

function cmp(a::String, b::String)
    i = start(a)
    j = start(b)
    while !done(a,i) && !done(b,i)
        c, i = next(a,i)
        d, j = next(b,j)
        if c != d
            return c < d ? -1 : +1
        end
    end
    done(a,i) && !done(b,j) ? -1 :
    !done(a,i) && done(b,j) ? +1 : 0
end

isequal(a::String, b::String) = cmp(a,b) == 0
isless(a::String, b::String)  = cmp(a,b) <  0

# faster comparisons for byte strings

cmp(a::ByteString, b::ByteString)     = lexcmp(a.data, b.data)
isequal(a::ByteString, b::ByteString) = length(a)==length(b) && cmp(a,b)==0

## character column width function ##

charwidth(c::Char) = max(0,int(ccall(:wcwidth, Int32, (Char,), c)))
strwidth(s::String) = (w=0; for c in s; w += charwidth(c); end; w)
strwidth(s::ByteString) = ccall(:u8_strwidth, Int, (Ptr{Uint8},), s.data)
# TODO: implement and use u8_strnwidth that takes a length argument

## generic string uses only length and next ##

type GenericString <: String
    string::String
end

length(s::GenericString) = length(s.string)
next(s::GenericString, i::Int) = next(s.string, i)

## plain old character arrays ##

type CharString <: String
    chars::Array{Char,1}

    CharString(a::Array{Char,1}) = new(a)
    CharString(c::Char...) = new([ c[i] | i=1:length(c) ])
end
CharString(x...) = CharString(map(char,x)...)

next(s::CharString, i::Int) = (s.chars[i], i+1)
length(s::CharString) = length(s.chars)
strlen(s::CharString) = length(s)

string(c::Char) = CharString(c)
string(c::Char, x::Char...) = CharString(c, x...)

## substrings reference original strings ##

type SubString <: String
    string::String
    offset::Int
    length::Int

    SubString(s::String, i::Int, j::Int) = new(s, i-1, j-i+1)
    SubString(s::SubString, i::Int, j::Int) =
        new(s.string, i-1+s.offset, j-i+1)
end
SubString(s::String, i::Integer, j::Integer) = SubString(s, int(i), int(j))

function next(s::SubString, i::Int)
    if i < 1 || i > s.length
        error("string index out of bounds")
    end
    c, i = next(s.string, i+s.offset)
    c, i-s.offset
end

length(s::SubString) = s.length
# TODO: strlen(s::SubString) = ??
# default implementation will work but it's slow
# can this be delegated efficiently somehow?
# that may require additional string interfaces

function ref(s::String, r::Range1{Int})
    if first(r) < 1 || length(s) < last(r)
        error("in substring slice: index out of range")
    end
    SubString(s, first(r), last(r))
end

## efficient representation of repeated strings ##

type RepString <: String
    string::String
    repeat::Integer
end

length(s::RepString) = length(s.string)*s.repeat
strlen(s::RepString) = strlen(s.string)*s.repeat

function next(s::RepString, i::Int)
    if i < 1 || i > length(s)
        error("string index out of bounds")
    end
    j = mod1(i,length(s.string))
    c, k = next(s.string, j)
    c, k-j+i
end

function repeat(s::String, r::Integer)
    r <  0 ? error("can't repeat a string ",r," times") :
    r == 0 ? "" :
    r == 1 ? s  :
    RepString(s,r)
end

## reversed strings without data movement ##

type RevString <: String
    string::String
end

length(s::RevString) = length(s.string)
strlen(s::RevString) = strlen(s.string)

start(s::RevString) = (n=length(s); n-thisind(s.string,n)+1)
function next(s::RevString, i::Int)
    n = length(s); j = n-i+1
    (s.string[j], n-thisind(s.string,j-1)+1)
end

reverse(s::String) = RevString(s)
reverse(s::RevString) = s.string

## ropes for efficient concatenation, etc. ##

# Idea: instead of this standard binary tree structure,
# how about we keep an array of substrings, with an
# offset array. We can do binary search on the offset
# array so we get O(log(n)) indexing time still, but we
# can compute the offsets lazily and avoid all the
# futzing around while the string is being constructed.

type RopeString <: String
    head::String
    tail::String
    depth::Int32
    length::Int

    RopeString(h::RopeString, t::RopeString) =
        depth(h.tail) + depth(t) < depth(h.head) ?
            RopeString(h.head, RopeString(h.tail, t)) :
            new(h, t, max(h.depth,t.depth)+1, length(h)+length(t))

    RopeString(h::RopeString, t::String) =
        depth(h.tail) < depth(h.head) ?
            RopeString(h.head, RopeString(h.tail, t)) :
            new(h, t, h.depth+1, length(h)+length(t))

    RopeString(h::String, t::RopeString) =
        depth(t.head) < depth(t.tail) ?
            RopeString(RopeString(h, t.head), t.tail) :
            new(h, t, t.depth+1, length(h)+length(t))

    RopeString(h::String, t::String) =
        new(h, t, 1, length(h)+length(t))
end

depth(s::String) = 0
depth(s::RopeString) = s.depth

function next(s::RopeString, i::Int)
    if i <= length(s.head)
        return next(s.head, i)
    else
        c, j = next(s.tail, i-length(s.head))
        return c, j+length(s.head)
    end
end

length(s::RopeString) = s.length
strlen(s::RopeString) = strlen(s.head) + strlen(s.tail)

strcat() = ""
strcat(s::String) = s
strcat(x...) = strcat(map(string,x)...)
strcat(s::String, t::String...) =
    (t = strcat(t...); isempty(s) ? t : isempty(t) ? s : RopeString(s, t))

print(s::RopeString) = print(s.head, s.tail)

## transformed strings ##

type TransformedString <: String
    transform::Function
    string::String
end

length(s::TransformedString) = length(s.string)
strlen(s::TransformedString) = strlen(s.string)

function next(s::TransformedString, i::Int)
    c, j = next(s.string,i)
    c = s.transform(c, i)
    return c, j
end

## uppercase and lowercase transformations ##

uppercase(c::Char) = ccall(:towupper, Char, (Char,), c)
lowercase(c::Char) = ccall(:towlower, Char, (Char,), c)

uppercase(s::String) = TransformedString((c,i)->uppercase(c), s)
lowercase(s::String) = TransformedString((c,i)->lowercase(c), s)

ucfirst(s::String) = TransformedString((c,i)->i==1 ? uppercase(c) : c, s)
lcfirst(s::String) = TransformedString((c,i)->i==1 ? lowercase(c) : c, s)

const uc = uppercase
const lc = lowercase

## string map ##

function map(f::Function, s::String)
    out = memio(length(s))
    for c in s
        write(out, f(c)::Char)
    end
    takebuf_string(out)
end

## conversion of general objects to strings ##

string(x) = print_to_string(show, x)
cstring(x...) = print_to_string(print, x...)

function cstring(p::Ptr{Uint8})
    p == C_NULL ? error("cannot convert NULL to string") :
    ccall(:jl_cstr_to_string, Any, (Ptr{Uint8},), p)::ByteString
end

## string promotion rules ##

promote_rule(::Type{UTF8String} , ::Type{ASCIIString}) = UTF8String
promote_rule(::Type{UTF8String} , ::Type{CharString} ) = UTF8String
promote_rule(::Type{ASCIIString}, ::Type{CharString} ) = UTF8String

## printing literal quoted string data ##

# TODO: this is really the inverse of print_unbackslashed

function print_quoted_literal(s::String)
    print('"')
    for c = s; c == '"' ? print("\\\"") : print(c); end
    print('"')
end

## string escaping & unescaping ##

escape_nul(s::String, i::Int) =
    !done(s,i) && '0' <= next(s,i)[1] <= '7' ? L"\x00" : L"\0"

is_hex_digit(c::Char) = '0'<=c<='9' || 'a'<=c<='f' || 'A'<=c<='F'
need_full_hex(s::String, i::Int) = !done(s,i) && is_hex_digit(next(s,i)[1])

function print_escaped(s::String, esc::String)
    i = start(s)
    while !done(s,i)
        c, j = next(s,i)
        c == '\0'       ? print(escape_nul(s,j)) :
        c == '\e'       ? print(L"\e") :
        c == '\\'       ? print("\\\\") :
        contains(esc,c) ? print('\\', c) :
        iswprint(c)     ? print(c) :
        7 <= c <= 13    ? print('\\', "abtnvfr"[c-6]) :
        c <= '\x7f'     ? print(L"\x", hex(c, 2)) :
        c <= '\uffff'   ? print(L"\u", hex(c, need_full_hex(s,j) ? 4 : 2)) :
                          print(L"\U", hex(c, need_full_hex(s,j) ? 8 : 4))
        i = j
    end
end

escape_string(s::String) = print_to_string(length(s), print_escaped, s, "\"")
print_quoted(s::String) = (print('"'); print_escaped(s, "\"\$"); print('"'))
#"  # work around syntax highlighting problem
quote_string(s::String) = print_to_string(length(s)+2, print_quoted, s)

# bare minimum unescaping function unescapes only given characters

function print_unescaped_chars(s::String, esc::String)
    if !contains(esc,'\\')
        esc = strcat("\\", esc)
    end
    i = start(s)
    while !done(s,i)
        c, i = next(s,i)
        if c == '\\' && !done(s,i) && contains(esc,s[i])
            c, i = next(s,i)
        end
        print(c)
    end
end

unescape_chars(s::String, esc::String) =
    print_to_string(length(s), print_unescaped_chars, s, esc)

# general unescaping of traditional C and Unicode escape sequences

function print_unescaped(s::String)
    i = start(s)
    while !done(s,i)
        c, i = next(s,i)
        if !done(s,i) && c == '\\'
            c, i = next(s,i)
            if c == 'x' || c == 'u' || c == 'U'
                n = k = 0
                m = c == 'x' ? 2 :
                    c == 'u' ? 4 : 8
                while (k+=1) <= m && !done(s,i)
                    c, j = next(s,i)
                    n = '0' <= c <= '9' ? n<<4 + c-'0' :
                        'a' <= c <= 'f' ? n<<4 + c-'a'+10 :
                        'A' <= c <= 'F' ? n<<4 + c-'A'+10 : break
                    i = j
                end
                if k == 1
                    error("\\x used with no following hex digits")
                end
                if m == 2 # \x escape sequence
                    write(uint8(n))
                else
                    print(char(n))
                end
            elseif '0' <= c <= '7'
                k = 1
                n = c-'0'
                while (k+=1) <= 3 && !done(s,i)
                    c, j = next(s,i)
                    n = '0' <= c <= '7' ? n<<3 + c-'0' : break
                    i = j
                end
                if n > 255
                    error("octal escape sequence out of range")
                end
                write(uint8(n))
            else
                print(c == 'a' ? '\a' :
                      c == 'b' ? '\b' :
                      c == 't' ? '\t' :
                      c == 'n' ? '\n' :
                      c == 'v' ? '\v' :
                      c == 'f' ? '\f' :
                      c == 'r' ? '\r' :
                      c == 'e' ? '\e' : c)
            end
        else
            print(c)
        end
    end
end

unescape_string(s::String) = print_to_string(length(s), print_unescaped, s)

## checking UTF-8 & ACSII validity ##

byte_string_classify(s::ByteString) =
    ccall(:u8_isvalid, Int32, (Ptr{Uint8}, Int), s.data, length(s))
    # 0: neither valid ASCII nor UTF-8
    # 1: valid ASCII
    # 2: valid UTF-8

is_valid_ascii(s::ByteString) = byte_string_classify(s) == 1
is_valid_utf8 (s::ByteString) = byte_string_classify(s) != 0

check_ascii(s::ByteString) = is_valid_ascii(s) ? s : error("invalid ASCII sequence")
check_utf8 (s::ByteString) = is_valid_utf8(s)  ? s : error("invalid UTF-8 sequence")

## string interpolation parsing ##

function _jl_interp_parse(s::String, unescape::Function, printer::Function)
    sx = {}
    i = j = start(s)
    while !done(s,j)
        c, k = next(s,j)
        if c == '$'
            if !isempty(s[i:j-1])
                push(sx, unescape(s[i:j-1]))
            end
            ex, j = parseatom(s,k)
            push(sx, ex)
            i = j
        elseif c == '\\' && !done(s,k)
            if s[k] == '$'
                if !isempty(s[i:j-1])
                    push(sx, unescape(s[i:j-1]))
                end
                i = k
            end
            c, j = next(s,k)
        else
            j = k
        end
    end
    if !isempty(s[i:])
        push(sx, unescape(s[i:j-1]))
    end
    length(sx) == 1 && isa(sx[1],ByteString) ? sx[1] :
        expr(:call, :print_to_string, printer, sx...)
end

_jl_interp_parse(s::String, u::Function) = _jl_interp_parse(s, u, print)
_jl_interp_parse(s::String) = _jl_interp_parse(s, x->check_utf8(unescape_string(x)))

function _jl_interp_parse_bytes(s::String)
    writer(x...) = for w=x; write(w); end
    _jl_interp_parse(s, unescape_string, writer)
end

## core string macros ##

macro   str(s); _jl_interp_parse(s); end
macro S_str(s); _jl_interp_parse(s); end
macro I_str(s); _jl_interp_parse(s, x->unescape_chars(x,"\"")); end
macro E_str(s); check_utf8(unescape_string(s)); end
macro B_str(s); _jl_interp_parse_bytes(s); end
macro b_str(s); ex = _jl_interp_parse_bytes(s); :(($ex).data); end

## shell-like command parsing ##

function _jl_shell_parse(s::String, interp::Bool)

    in_single_quotes = false
    in_double_quotes = false

    args = {}
    arg = {}
    i = start(s)
    j = i

    function update_arg(x)
        if !isa(x,String) || !isempty(x)
            push(arg, x)
        end
    end
    function append_arg()
        if isempty(arg); arg = {"",}; end
        push(args, arg)
        arg = {}
    end

    while !done(s,j)
        c, k = next(s,j)
        if !in_single_quotes && !in_double_quotes && iswspace(c)
            update_arg(s[i:j-1])
            append_arg()
            j = k
            while !done(s,j)
                c, k = next(s,j)
                if !iswspace(c)
                    i = j
                    break
                end
                j = k
            end
        elseif interp && !in_single_quotes && c == '$'
            update_arg(s[i:j-1]); i = k; j = k
            if done(s,k)
                error("\$ right before end of command")
            end
            if iswspace(s[k])
                error("space not allowed right after \$")
            end
            ex, j = parseatom(s,j)
            update_arg(ex); i = j
        else
            if !in_double_quotes && c == '\''
                in_single_quotes = !in_single_quotes
                update_arg(s[i:j-1]); i = k
            elseif !in_single_quotes && c == '"'
                in_double_quotes = !in_double_quotes
                update_arg(s[i:j-1]); i = k
            elseif c == '\\'
                if in_double_quotes
                    if done(s,k)
                        error("unterminated double quote")
                    end
                    if s[k] == '"' || s[k] == '$'
                        update_arg(s[i:j-1]); i = k
                        c, k = next(s,k)
                    end
                elseif !in_single_quotes
                    if done(s,k)
                        error("dangling backslash")
                    end
                    update_arg(s[i:j-1]); i = k
                    c, k = next(s,k)
                end
            end
            j = k
        end
    end

    if in_single_quotes; error("unterminated single quote"); end
    if in_double_quotes; error("unterminated double quote"); end

    update_arg(s[i:])
    append_arg()

    if !interp
        return args
    end

    # construct an expression
    exprs = {}
    for arg in args
        push(exprs, expr(:tuple, arg))
    end
    expr(:tuple,exprs)
end
_jl_shell_parse(s::String) = _jl_shell_parse(s,true)

function shell_split(s::String)
    parsed = _jl_shell_parse(s,false)
    args = String[]
    for arg in parsed
       push(args, strcat(arg...))
    end
    args
end

function print_shell_word(word::String)
    if isempty(word)
        print("''")
    end
    has_single = false
    has_special = false
    for c in word
        if iswspace(c) || c=='\\' || c=='\'' || c=='"' || c=='$'
            has_special = true
            if c == '\''
                has_single = true
            end
        end
    end
    if !has_special
        print(word)
    elseif !has_single
        print('\'', word, '\'')
    else
        print('"')
        for c in word
            if c == '"' || c == '$'
                print('\\')
            end
            print(c)
        end
        print('"')
    end
end

function print_shell_escaped(cmd::String, args::String...)
    print_shell_word(cmd)
    for arg in args
        print(' ')
        print_shell_word(arg)
    end
end

shell_escape(cmd::String, args::String...) =
    print_to_string(print_shell_escaped, cmd, args...)

## interface to parser ##

function parse(s::String, pos, greedy)
    # returns (expr, end_pos). expr is () in case of parse error.
    ex, pos = ccall(:jl_parse_string, Any,
                    (Ptr{Uint8}, Int32, Int32),
                    cstring(s), pos-1, greedy ? 1:0)
    if isa(ex,Expr) && is(ex.head,:error)
        throw(ParseError(ex.args[1]))
    end
    if ex == (); throw(ParseError("end of input")); end
    ex, pos+1 # C is zero-based, Julia is 1-based
end

parse(s::String)          = parse(s, 1, true)
parse(s::String, pos)     = parse(s, pos, true)
parseatom(s::String)      = parse(s, 1, false)
parseatom(s::String, pos) = parse(s, pos, false)

## miscellaneous string functions ##

function lpad(s::String, n::Integer, p::String)
    m = n - strlen(s)
    if m <= 0; return s; end
    l = strlen(p)
    if l==1
        return p^m * s
    end
    q = div(m,l)
    r = m - q*l
    cstring(p^q*p[1:chr2ind(p,r)]*s)
end

function rpad(s::String, n::Integer, p::String)
    m = n - strlen(s)
    if m <= 0; return s; end
    l = strlen(p)
    if l==1
        return s * p^m
    end
    q = div(m,l)
    r = m - q*l
    cstring(s*p^q*p[1:chr2ind(p,r)])
end

lpad(s, n::Integer, p) = lpad(string(s), n, string(p))
rpad(s, n::Integer, p) = rpad(string(s), n, string(p))

lpad(s, n::Integer) = lpad(string(s), n, " ")
rpad(s, n::Integer) = rpad(string(s), n, " ")

function split(s::String, delims, include_empty::Bool)
    i = 1
    strs = String[]
    len = length(s)
    while true
        tokstart = tokend = i
        while !done(s,i)
            (c,i) = next(s,i)
            if contains(delims, c)
                break
            end
            tokend = i
        end
        tok = s[tokstart:(tokend-1)]
        if include_empty || !isempty(tok)
            push(strs, tok)
        end
        if !((i <= len) || (i==len+1 && tokend!=i))
            break
        end
    end
    strs
end

split(s::String) = split(s, (' ','\t','\n','\v','\f','\r'), false)
split(s::String, x) = split(s, x, true)
split(s::String, x::Char, incl::Bool) = split(s, (x,), incl)

function print_joined(strings, delim, last)
    i = start(strings)
    if done(strings,i)
        return
    end
    str, i = next(strings,i)
    print(str)
    while !done(strings,i)
        str, i = next(strings,i)
        print(done(strings,i) ? last : delim)
        print(str)
    end
end

function print_joined(strings, delim)
    i = start(strings)
    while !done(strings,i)
        str, i = next(strings,i)
        print(str)
        if !done(strings,i)
            print(delim)
        end
    end
end
print_joined(strings) = print_joined(strings, "")

join(args...) = print_to_string(print_joined, args...)

chop(s::String) = s[1:thisind(s,length(s))-1]
chomp(s::String) = (i=thisind(s,length(s)); s[i]=='\n' ? s[1:i-1] : s)
chomp(s::ByteString) = s.data[end]==0x0a ? s[1:end-1] : s

function lstrip(s::String)
    i = start(s)
    while !done(s,i)
        c, j = next(s,i)
        if !iswspace(c)
            return s[i:end]
        end
        i = j
    end
    ""
end

function rstrip(s::String)
    r = reverse(s)
    i = start(r)
    while !done(r,i)
        c, j = next(r,i)
        if !iswspace(c)
            return s[1:end-i+1]
        end
        i = j
    end
    ""
end

strip(s::String) = lstrip(rstrip(s))

## string to integer functions ##

function parse_int{T<:Integer}(::Type{T}, s::String, base::Integer)
    if !(2 <= base <= 36); error("invalid base: ",base); end
    i = start(s)
    if done(s,i)
        error("premature end of integer (in ",show_to_string(s),")")
    end
    c,i = next(s,i)
    sgn = one(T)
    if T <: Signed && c == '-'
        sgn = -sgn
        if done(s,i)
            error("premature end of integer (in ",show_to_string(s),")")
        end
        c,i = next(s,i)
    end
    base = convert(T,base)
    n::T = 0
    while true
        d = '0' <= c <= '9' ? c-'0' :
            'A' <= c <= 'Z' ? c-'A'+10 :
            'a' <= c <= 'z' ? c-'a'+10 : typemax(Int)
        if d >= base
            error(show_to_string(c)," is not a valid digit (in ",show_to_string(s),")")
        end
        # TODO: overflow detection?
        n = n*base + d
        if done(s,i)
            break
        end
        c,i = next(s,i)
    end
    return flipsign(n,sgn)
end

parse_int(s::String, base::Integer) = parse_int(Int,s,base)
parse_int(T::Type, s::String)       = parse_int(T,s,10)
parse_int(s::String)                = parse_int(Int,s,10)

parse_bin(T::Type, s::String) = parse_int(T,s,2)
parse_oct(T::Type, s::String) = parse_int(T,s,8)
parse_hex(T::Type, s::String) = parse_int(T,s,16)

parse_bin(s::String) = parse_int(Int,s,2)
parse_oct(s::String) = parse_int(Int,s,8)
parse_hex(s::String) = parse_int(Int,s,16)

integer (s::String) = int(s)
unsigned(s::String) = uint(s)
int     (s::String) = parse_int(Int,s)
uint    (s::String) = parse_int(Uint,s)
int8    (s::String) = parse_int(Int8,s)
uint8   (s::String) = parse_int(Uint8,s)
int16   (s::String) = parse_int(Int16,s)
uint16  (s::String) = parse_int(Uint16,s)
int32   (s::String) = parse_int(Int32,s)
uint32  (s::String) = parse_int(Uint32,s)
int64   (s::String) = parse_int(Int64,s)
uint64  (s::String) = parse_int(Uint64,s)

## integer to string functions ##

const _jl_dig_syms = "0123456789abcdefghijklmnopqrstuvwxyz".data

function int2str(n::Union(Int64,Uint64), b::Integer, l::Int)
    if b < 2 || b > 36; error("int2str: invalid base ", b); end
    neg = n < 0
    n = unsigned(abs(n))
    b = convert(typeof(n), b)
    ndig = ndigits(n, b)
    sz = max(convert(Int, ndig), l) + neg
    data = Array(Uint8, sz)
    i = sz
    if ispow2(b)
        digmask = b-1
        shift = trailing_zeros(b)
        while i > neg
            ch = n & digmask
            data[i] = _jl_dig_syms[int(ch)+1]
            n >>= shift
            i -= 1
        end
    else
        while i > neg
            ch = n % b
            data[i] = _jl_dig_syms[int(ch)+1]
            n = div(n,b)
            i -= 1
        end
    end
    if neg
        data[1] = '-'
    end
    ASCIIString(data)
end
int2str(n::Integer, b::Integer)         = int2str(n, b, 0)
int2str(n::Integer, b::Integer, l::Int) = int2str(int64(n), b, l)

string(x::Signed) = dec(int64(x))
cstring(x::Signed) = dec(int64(x))

## string to float functions ##

function float64_isvalid(s::String, out::Array{Float64,1})
    s = cstring(s)
    return (ccall(:jl_strtod, Int32, (Ptr{Uint8},Ptr{Float64}), s, out)==0)
end

function float32_isvalid(s::String, out::Array{Float32,1})
    s = cstring(s)
    return (ccall(:jl_strtof, Int32, (Ptr{Uint8},Ptr{Float32}), s, out)==0)
end

begin
    local tmp::Array{Float64,1} = Array(Float64,1)
    local tmpf::Array{Float32,1} = Array(Float32,1)
    global float64, float32
    function float64(s::String)
        if !float64_isvalid(s, tmp)
            throw(ArgumentError("float64(String): invalid number format"))
        end
        return tmp[1]
    end

    function float32(s::String)
        if !float32_isvalid(s, tmpf)
            throw(ArgumentError("float32(String): invalid number format"))
        end
        return tmpf[1]
    end
end

float(x::String) = float64(x)
parse_float(x::String) = float64(x)
parse_float(::Type{Float64}, x::String) = float64(x)
parse_float(::Type{Float32}, x::String) = float32(x)

# copying a byte string (generally not needed due to "immutability")

strcpy{T<:ByteString}(s::T) = T(copy(s.data))

# lexicographically compare byte arrays (used by Latin-1 and UTF-8)

function lexcmp(a::Array{Uint8,1}, b::Array{Uint8,1})
    c = ccall(:memcmp, Int32, (Ptr{Uint8}, Ptr{Uint8}, Uint),
              a, b, min(length(a),length(b)))
    c < 0 ? -1 : c > 0 ? +1 : cmp(length(a),length(b))
end

# find the index of the first occurrence of a byte value in a byte array

function memchr(a::Array{Uint8,1}, b::Integer)
    p = pointer(a)
    q = ccall(:memchr, Ptr{Uint8}, (Ptr{Uint8}, Int32, Uint), p, b, length(a))
    q == C_NULL ? 0 : q - p + 1
end

# concatenate byte arrays into a single array

memcat() = Array(Uint8,0)
memcat(a::Array{Uint8,1}) = copy(a)

function memcat(arrays::Array{Uint8,1}...)
    n = 0
    for a in arrays
        n += length(a)
    end
    arr = Array(Uint8, n)
    ptr = pointer(arr)
    offset = 0
    for a in arrays
        ccall(:memcpy, Ptr{Uint8}, (Ptr{Uint8}, Ptr{Uint8}, Uint),
              ptr+offset, a, length(a))
        offset += length(a)
    end
    return arr
end

# concatenate the data fields of byte strings

memcat(s::ByteString) = memcat(s.data)
memcat(sx::ByteString...) = memcat(map(s->s.data, sx)...)

max_stars_repo_name: btashton/pygments

max_issues_repo_name: btashton/pygments

max_forks_repo_name: btashton/pygments