text:
"""
Default PowerSimulations Operation Problem Type
"""
struct GenericOpProblem <: DecisionProblem end

"""
    DecisionModel(::Type{M},
    template::ProblemTemplate,
    sys::PSY.System,
    jump_model::Union{Nothing, JuMP.Model}=nothing;
    kwargs...) where {M<:DecisionProblem,
                      T<:PM.AbstractPowerFormulation}

This builds the optimization problem of type M with the specific system and template.

# Arguments

  - `::Type{M} where M<:DecisionProblem`: The abstract operation model type
  - `template::ProblemTemplate`: The model reference made up of transmission, devices,
    branches, and services.
  - `sys::PSY.System`: the system created using Power Systems
  - `jump_model::Union{Nothing, JuMP.Model}`: Enables passing a custom JuMP model. Use with care

# Output

  - `model::DecisionModel`: The operation model containing the model type, built JuMP model, Power
    Systems system.

# Example

```julia
template = ProblemTemplate(CopperPlatePowerModel, devices, branches, services)
OpModel = DecisionModel(MockOperationProblem, template, system)
```

# Accepted Key Words

  - `optimizer`: The optimizer that will be used in the optimization model.
  - `horizon::Int`: Manually specify the length of the forecast Horizon
  - `warm_start::Bool`: True will use the current operation point in the system to initialize variable values. False initializes all variables to zero. Default is true
  - `system_to_file::Bool:`: True to create a copy of the system used in the model. Default true.
  - `export_pwl_vars::Bool`: True to export all the pwl intermediate variables. It can slow down significantly the solve time. Default is false.
  - `allow_fails::Bool`: True to allow the simulation to continue even if the optimization step fails. Use with care, default to false.
  - `optimizer_solve_log_print::Bool`: True to print the optimizer solve log. Default is false.
  - `direct_mode_optimizer::Bool` True to use the solver in direct mode. Creates a [JuMP.direct_model](https://jump.dev/JuMP.jl/dev/reference/models/#JuMP.direct_model). Default is false.
  - `initial_time::Dates.DateTime`: Initial Time for the model solve
  - `time_series_cache_size::Int`: Size in bytes to cache for each time array. Default is 1 MiB. Set to 0 to disable.
"""
mutable struct DecisionModel{M <: DecisionProblem} <: OperationModel
    name::Symbol
    template::ProblemTemplate
    sys::PSY.System
    internal::Union{Nothing, ModelInternal}
    store::DecisionModelStore
    ext::Dict{String, Any}

    function DecisionModel{M}(
        template::ProblemTemplate,
        sys::PSY.System,
        settings::Settings,
        jump_model::Union{Nothing, JuMP.Model}=nothing;
        name=nothing,
    ) where {M <: DecisionProblem}
        if name === nothing
            name = Symbol(typeof(template))
        elseif name isa String
            name = Symbol(name)
        end
        _, _, forecast_count = PSY.get_time_series_counts(sys)
        if forecast_count < 1
            error(
                "The system does not contain forecast data. A DecisionModel can't be built.",
            )
        end
        internal = ModelInternal(
            OptimizationContainer(sys, settings, jump_model, PSY.Deterministic),
        )
        return new{M}(
            name,
            template,
            sys,
            internal,
            DecisionModelStore(),
            Dict{String, Any}(),
        )
    end
end

function DecisionModel{M}(
    template::ProblemTemplate,
    sys::PSY.System,
    jump_model::Union{Nothing, JuMP.Model}=nothing;
    name=nothing,
    optimizer=nothing,
    horizon=UNSET_HORIZON,
    warm_start=true,
    system_to_file=true,
    initialize_model=true,
    initialization_file="",
    deserialize_initial_conditions=false,
    export_pwl_vars=false,
    allow_fails=false,
    optimizer_solve_log_print=false,
    detailed_optimizer_stats=false,
    calculate_conflict=false,
    direct_mode_optimizer=false,
    check_numerical_bounds=true,
    initial_time=UNSET_INI_TIME,
    time_series_cache_size::Int=IS.TIME_SERIES_CACHE_SIZE_BYTES,
) where {M <: DecisionProblem}
    settings = Settings(
        sys;
        horizon=horizon,
        initial_time=initial_time,
        optimizer=optimizer,
        time_series_cache_size=time_series_cache_size,
        warm_start=warm_start,
        system_to_file=system_to_file,
        initialize_model=initialize_model,
        initialization_file=initialization_file,
        deserialize_initial_conditions=deserialize_initial_conditions,
        export_pwl_vars=export_pwl_vars,
        allow_fails=allow_fails,
        calculate_conflict=calculate_conflict,
        optimizer_solve_log_print=optimizer_solve_log_print,
        detailed_optimizer_stats=detailed_optimizer_stats,
        direct_mode_optimizer=direct_mode_optimizer,
        check_numerical_bounds=check_numerical_bounds,
    )
    return DecisionModel{M}(template, sys, settings, jump_model, name=name)
end

"""
    DecisionModel(::Type{M},
    template::ProblemTemplate,
    sys::PSY.System,
    optimizer::MOI.OptimizerWithAttributes,
    jump_model::Union{Nothing, JuMP.Model}=nothing;
    kwargs...) where {M <: DecisionProblem}

This builds the optimization problem of type M with the specific system and template

# Arguments

  - `::Type{M} where M<:DecisionProblem`: The abstract operation model type
  - `template::ProblemTemplate`: The model reference made up of transmission, devices,
    branches, and services.
  - `sys::PSY.System`: the system created using Power Systems
  - `jump_model::Union{Nothing, JuMP.Model}`: Enables passing a custom JuMP model. Use with care

# Output

  - `Stage::DecisionProblem`: The operation model containing the model type, unbuilt JuMP model, Power
    Systems system.

# Example

template = ProblemTemplate(CopperPlatePowerModel, devices, branches, services)
problem = DecisionModel(MyOpProblemType template, system, optimizer)

```
# Accepted Key Words
- `initial_time::Dates.DateTime`: Initial Time for the model solve
- `warm_start::Bool` True will use the current operation point in the system to initialize variable values. False initializes all variables to zero. Default is true
- `export_pwl_vars::Bool` True will write the results of the piece-wise-linear intermediate variables. Slows down the simulation process significantly
- `allow_fails::Bool` True will allow the simulation to continue if the optimizer can't find a solution. Use with care, can lead to unwanted behaviour or results
- `optimizer_solve_log_print::Bool` Uses JuMP.unset_silent() to print the optimizer's log. By default all solvers are set to `MOI.Silent()`
- `name`: name of model, string or symbol; defaults to the type of template converted to a symbol
```
"""
function DecisionModel(
    ::Type{M},
    template::ProblemTemplate,
    sys::PSY.System,
    jump_model::Union{Nothing, JuMP.Model}=nothing;
    kwargs...,
) where {M <: DecisionProblem}
    return DecisionModel{M}(template, sys, jump_model; kwargs...)
end

function DecisionModel(
    template::ProblemTemplate,
    sys::PSY.System,
    jump_model::Union{Nothing, JuMP.Model}=nothing;
    kwargs...,
)
    return DecisionModel{GenericOpProblem}(template, sys, jump_model; kwargs...)
end

"""
    DecisionModel(directory::AbstractString)

Construct an DecisionProblem from a serialized file.

# Arguments

  - `directory::AbstractString`: Directory containing a serialized model
  - `jump_model::Union{Nothing, JuMP.Model}` = nothing: The JuMP model does not get
    serialized. Callers should pass whatever they passed to the original problem.
  - `optimizer::Union{Nothing,MOI.OptimizerWithAttributes}` = nothing: The optimizer does
    not get serialized. Callers should pass whatever they passed to the original problem.
  - `system::Union{Nothing, PSY.System}`: Optionally, the system used for the model.
    If nothing and sys_to_file was set to true when the model was created, the system will
    be deserialized from a file.
"""
function DecisionModel(
    directory::AbstractString,
    optimizer::MOI.OptimizerWithAttributes;
    jump_model::Union{Nothing, JuMP.Model}=nothing,
    system::Union{Nothing, PSY.System}=nothing,
)
    return deserialize_problem(
        DecisionModel,
        directory;
        jump_model=jump_model,
        optimizer=optimizer,
        system=system,
    )
end

# Probably could be more efficient by storing the info in the internal
function get_current_time(model::DecisionModel)
    execution_count = get_internal(model).execution_count
    initial_time = get_initial_time(model)
    interval = get_interval(model.internal.store_parameters)
    return initial_time + interval * execution_count
end

function init_model_store_params!(model::DecisionModel)
    num_executions = get_executions(model)
    horizon = get_horizon(model)
    system = get_system(model)
    interval = PSY.get_forecast_interval(system)
    resolution = PSY.get_time_series_resolution(system)
    base_power = PSY.get_base_power(system)
    sys_uuid = IS.get_uuid(system)
    model.internal.store_parameters = ModelStoreParams(
        num_executions,
        horizon,
        iszero(interval) ? resolution : interval,
        resolution,
        base_power,
        sys_uuid,
        get_metadata(get_optimization_container(model)),
    )
    return
end

function build_pre_step!(model::DecisionModel)
    TimerOutputs.@timeit BUILD_PROBLEMS_TIMER "Build pre-step" begin
        if !is_empty(model)
            @info "OptimizationProblem status not BuildStatus.EMPTY. Resetting"
            reset!(model)
        end
        # Initial time are set here because the information is specified in the
        # Simulation Sequence object and not at the problem creation.

        @info "Initializing Optimization Container For a DecisionModel"
        init_optimization_container!(
            get_optimization_container(model),
            get_network_formulation(get_template(model)),
            get_system(model),
        )
        @info "Initializing ModelStoreParams"
        init_model_store_params!(model)

        @info "Mapping Service Models"
        populate_aggregated_service_model!(get_template(model), get_system(model))
        populate_contributing_devices!(get_template(model), get_system(model))
        add_services_to_device_model!(get_template(model))

        handle_initial_conditions!(model)
        set_status!(model, BuildStatus.IN_PROGRESS)
    end
    return
end

get_horizon(model::DecisionModel) = get_horizon(get_settings(model))

"""
Implementation of build for any DecisionProblem
"""
function build!(
    model::DecisionModel{<:DecisionProblem};
    output_dir::String,
    recorders=[],
    console_level=Logging.Error,
    file_level=Logging.Info,
    disable_timer_outputs=false,
)
    mkpath(output_dir)
    set_output_dir!(model, output_dir)
    set_console_level!(model, console_level)
    set_file_level!(model, file_level)
    TimerOutputs.reset_timer!(BUILD_PROBLEMS_TIMER)
    disable_timer_outputs && TimerOutputs.disable_timer!(BUILD_PROBLEMS_TIMER)
    file_mode = "w"
    add_recorders!(model, recorders)
    register_recorders!(model, file_mode)
    logger = configure_logging(model.internal, file_mode)
    try
        Logging.with_logger(logger) do
            try
                TimerOutputs.@timeit BUILD_PROBLEMS_TIMER "Problem $(get_name(model))" begin
                    build_impl!(model)
                end
                set_status!(model, BuildStatus.BUILT)
                @info "\n$(BUILD_PROBLEMS_TIMER)\n"
            catch e
                set_status!(model, BuildStatus.FAILED)
                bt = catch_backtrace()
                @error "DecisionModel Build Failed" exception = e, bt
            end
        end
    finally
        unregister_recorders!(model)
        close(logger)
    end
    return get_status(model)
end

"""
Default implementation of build method for Operational Problems for models conforming with
DecisionProblem specification. Overload this function to implement a custom build method
"""
function build_model!(model::DecisionModel)
    build_impl!(get_optimization_container(model), get_template(model), get_system(model))
    return
end

function reset!(model::DecisionModel)
    # TODO-PJ: This is needed until we remove the ParameterJuMP dependency
    was_built_for_recurrent_solves = built_for_recurrent_solves(model)
    if was_built_for_recurrent_solves
        set_execution_count!(model, 0)
    end
    model.internal.container = OptimizationContainer(
        get_system(model),
        get_settings(model),
        nothing,
        PSY.Deterministic,
    )
    model.internal.container.built_for_recurrent_solves = was_built_for_recurrent_solves
    model.internal.ic_model_container = nothing
    empty_time_series_cache!(model)
    empty!(get_store(model))
    set_status!(model, BuildStatus.EMPTY)
    return
end

"""
Default solve method for models that conform to the requirements of
DecisionModel{<: DecisionProblem}.

This will call [`build!`](@ref) on the model if it is not already built. It will forward all
keyword arguments to that function.

# Arguments

  - `model::OperationModel = model`: operation model
  - `optimizer::MOI.OptimizerWithAttributes`: The optimizer that is used to solve the model
  - `export_problem_results::Bool`: If true, export ProblemResults DataFrames to CSV files.
  - `serialize::Bool`: If true, serialize the model to a file to allow re-execution later.

# Examples

```julia
results = solve!(OpModel)
results = solve!(OpModel, output_dir="output")
```
"""
function solve!(
    model::DecisionModel{<:DecisionProblem};
    export_problem_results=false,
    console_level=Logging.Error,
    file_level=Logging.Info,
    disable_timer_outputs=false,
    serialize=true,
    kwargs...,
)
    build_if_not_already_built!(
        model;
        console_level=console_level,
        file_level=file_level,
        disable_timer_outputs=disable_timer_outputs,
        kwargs...,
    )
    set_console_level!(model, console_level)
    set_file_level!(model, file_level)
    TimerOutputs.reset_timer!(RUN_OPERATION_MODEL_TIMER)
    disable_timer_outputs && TimerOutputs.disable_timer!(RUN_OPERATION_MODEL_TIMER)
    file_mode = "a"
    register_recorders!(model, file_mode)
    logger = configure_logging(model.internal, file_mode)
    optimizer = get(kwargs, :optimizer, nothing)
    try
        Logging.with_logger(logger) do
            try
                initialize_storage!(
                    get_store(model),
                    get_optimization_container(model),
                    model.internal.store_parameters,
                )
                TimerOutputs.@timeit RUN_OPERATION_MODEL_TIMER "Solve" begin
                    _pre_solve_model_checks(model, optimizer)
                    solve_impl!(model)
                    current_time = get_initial_time(model)
                    write_results!(get_store(model), model, current_time, current_time)
                    write_optimizer_stats!(
                        get_store(model),
                        get_optimizer_stats(model),
                        current_time,
                    )
                end
                if serialize
                    TimerOutputs.@timeit RUN_OPERATION_MODEL_TIMER "Serialize" begin
                        serialize_problem(model, optimizer=optimizer)
                        serialize_optimization_model(model)
                    end
                end
                TimerOutputs.@timeit RUN_OPERATION_MODEL_TIMER "Results processing" begin
                    # TODO: This could be more complicated than it needs to be
                    results = ProblemResults(model)
                    serialize_results(results, get_output_dir(model))
                    export_problem_results && export_results(results)
                end
                @info "\n$(RUN_OPERATION_MODEL_TIMER)\n"
            catch e
                @error "Decision Problem solve failed" exception = (e, catch_backtrace())
                # TODO: Run IIS here if the solve called failed
                set_run_status!(model, RunStatus.FAILED)
            end
        end
    finally
        unregister_recorders!(model)
        close(logger)
    end

    return get_run_status(model)
end

"""
Default solve method for a DecisionModel used inside of a Simulation. Solves problems that conform to the requirements of DecisionModel{<: DecisionProblem}

# Arguments

  - `step::Int`: Simulation Step
  - `model::OperationModel`: operation model
  - `start_time::Dates.DateTime`: Initial Time of the simulation step in Simulation time.
  - `store::SimulationStore`: Simulation output store

# Accepted Key Words

  - `exports`: realtime export of output. Use wisely, it can have negative impacts in the simulation times
"""
function solve!(
    step::Int,
    model::DecisionModel{<:DecisionProblem},
    start_time::Dates.DateTime,
    store::SimulationStore;
    exports=nothing,
)
    # Note, we don't call solve!(decision_model) here because the solve call includes a lot of
    # other logic used when solving the models separate from a simulation
    solve_impl!(model)
    IS.@assert_op get_current_time(model) == start_time
    if get_run_status(model) == RunStatus.SUCCESSFUL
        write_results!(store, model, start_time, start_time; exports=exports)
        write_optimizer_stats!(store, model, start_time)
        advance_execution_count!(model)
    end
    return get_run_status(model)
end

function update_parameters!(
    model::DecisionModel,
    decision_states::DatasetContainer{DataFrameDataset},
)
    for key in keys(get_parameters(model))
        update_parameter_values!(model, key, decision_states)
    end
    if !is_synchronized(model)
        update_objective_function!(get_optimization_container(model))
    end
    return
end


text:
#####
##### Methods for (co)geodesic flow
#####
##### These are temporary, until Manifolds.jl provides the same functionality
#####

export geodesic_flow, geodesic_flow!, cogeodesic_flow, cogeodesic_flow!

@doc doc"""
      geodesic_flow(M, x₀, v₀, t)

Follow the geodesic with initial point $x_0 \in M$ and velocity
$v_0 \in T_{x_0} M$ for time $t$, dragging $v_0$ along.
"""
function geodesic_flow(M, x₀, v₀, t)
      x = similar_result(M, geodesic_flow, x₀, t)
      v = similar_result(M, geodesic_flow, v₀, t)
      geodesic_flow!(M, x, v, x₀, v₀, t)
      return (x, v)
end

@doc doc"""
      geodesic_flow!(M, x, v, x₀, v₀, t)

Follow the geodesic with initial point $x_0 \in M$ and velocity
$v_0 \in T_{x_0} M$ for time $t$, dragging $v_0$ along. The resulting
point and tangent vector are stored in $x$ and $v$, respectively.
"""
function geodesic_flow! end

@traitfn function geodesic_flow!(M::MT,
                                 x,
                                 v,
                                 x₀,
                                 v₀,
                                 t) where {MT;!IsDecoratorManifold{MT}}
      result = DiffResults.DiffResult(x, (v,))
      result = ForwardDiff.derivative!(result, (x,t)->exp!(M,x,x₀,v₀,t), x, t)
      copyto!(v, DiffResults.derivative(result))
      return (x, v)
end

@traitfn function geodesic_flow!(M::MT,
                                 x,
                                 v,
                                 x₀,
                                 v₀,
                                 t) where {MT;IsDecoratorManifold{MT}}
      return geodesic_flow!(base_manifold(M), x, v, x₀, v₀, t)
end

@doc doc"""
      cogeodesic_flow(M, x₀, p₀, t)

Follow the geodesic with initial point $x_0 \in M$ and momentum
$p_0 \in T_{x_0}^* M$ for time $t$, dragging $p_0$ along.
"""
function cogeodesic_flow(M, x₀, p₀, t)
      x = similar_result(M, cogeodesic_flow, x₀)
      p = similar_result(M, cogeodesic_flow, p₀)
      cogeodesic_flow!(M, x, p, x₀, p₀, t)
      return (x, p)
end

@doc doc"""
      cogeodesic_flow!(M, x, p, x₀, p₀, t)

Follow the geodesic with initial point $x_0 \in M$ and momentum
$p_0 \in T_{x_0}^* M$ for time $t$, dragging $p_0$ along. The resulting
point and cotangent vector are stored in $x$ and $p$, respectively.
"""
function cogeodesic_flow! end

@traitfn function cogeodesic_flow!(M::MT,
                                   x,
                                   p,
                                   x₀,
                                   p₀,
                                   t) where {MT;!IsDecoratorManifold{MT}}
      error("cogeodesic_flow! not implemented for Manifold $(typeof(M)), ",
            "input point $(typeof(x₀)), input vector $(typeof(p₀)), time ",
            "$typeof(t), output point $(typeof(x)), and output vector ",
            "$(typeof(p))")
end

@traitfn function cogeodesic_flow!(M::MT,
                                   x,
                                   p,
                                   x₀,
                                   p₀,
                                   t) where {MT;IsDecoratorManifold{MT}}
      return cogeodesic_flow!(base_manifold(M), x, p, x₀, p₀, t)
end

ext: jl


text:
# Compute a complete set of sums using a fast recursive algorithm based on a simple combinatorial relationship between the required sums.
# This code has been partially optimized for quick updates, which means that it is more complicated that it would need to be if we were just updating from scratch every time.

# Throughout these functions:
# - N::Dict{Vector{Int64},<:Integer} is a Dict mapping sorted partition vectors to their volume sums.
# - M::Dict{Vector{Int64},<:Integer} is a Dict mapping sorted partition vectors to their *uncorrected* sums. By *uncorrected*, we mean that these sums must be multipled by a set of combinatorial constants prior to including in the modularity function. Each entry of M differs from the corresponding entry of N by a combinatorial constant.
# - p::Vector{Int64} is a partition vector
# - D::Array{Int64, 1} is the degree sequence of the hypergraph, which does not change.
# - Z::Array{Int64, 1} is the group label vector, which does change.
# - ℓ::Int64 is the number of groups.

# ------------------------------------------------------------------------------
# COMPUTATION OF SUMS FROM SCRATCH
# ------------------------------------------------------------------------------

function correctOvercounting(M::Dict{Vector{Int64},<:Integer}, p::Vector{Int64})
    """
    Internal function: no reason to call this outside of evalSums().
    Computes the second term in the recurrence relation given in the paper.
    M::Dict{Vector{Int64},<:Integer}, a Dict mapping sorted partition vectors to their *uncorrected* volume-sums,
    p::Vector{Int64}, the partition for which we are calculating.
    returns S::bigInt, the second term in the recurrence relation.
    """
    pk = p[end]
    S = 0
    for i = 1:length(p)-1
        p_ = copy(p)[1:(end-1)]
        p_[i] += pk
        S += M[-sort(-p_)]
    end
    return(S)
end

function computeMoments(Z::Vector{<:Integer}, D::Vector{<:Integer}, r::Integer, ℓ::Integer=0)
    """
    Compute the vectors V of group volumes and μ of volume-moments.
    No reason to call this outside of evalSums()
    Z::Vector{Int64}, the vector of cluster labels.
    D::Vector{Int64}, the degree sequence
    r::Int64, the size of the largest hyperedge
    ℓ::Int64, the number of groups (defaults to maximum(Z))
    returns V::Array{Int64, 1} the vector of group volumes and μ::Vector{Int64} of volume-moments.
    NOTE: might want to implement bigInts in μ
    """
    if ℓ==0
        ℓ=maximum(Z)
    end

    V = [sum(D[Z .== j]) for j in 1:ℓ]
    μ = [sum(V.^i) for i = 1:r]
    return(V, μ)
end

function evalConstants(r::Int64; ℓ = r, bigInt::Bool=true)
    """
    Evaluate the combinatorial constants required to convert the "convenient" sums M to the required sums N required for modularity calculations.
    r::The size of the largest hyperedge.
    bigInt:Bool, whether to use bigInts for this computation, recommended unless the instance is VERY small.
    return: C::Dict{Array{Int64, 1}, bigInt} a Dict mapping a partition vector to its associated combinatorial constant.
    """

    C =
        if bigInt C = Dict{Vector{Int64},BigInt}()
        else      C = Dict{Vector{Int64},Int64}()
        end

    for i = 1:r, j = 1:minimum([i,ℓ]), p in Combinatorics.partitions(i, j)
        orderCorrection =
            if bigInt prod([factorial(big(c)) for c in values(StatsBase.countmap(p))])
            else      prod([factorial(c) for c in values(StatsBase.countmap(p))])
            end
        # orderCorrection = 1
        C[p] = Combinatorics.multinomial(big.(p)...) ÷ orderCorrection
    end
    return(C)
end;

function evalSums(Z::Vector{<:Integer}, D::Vector{Int64}, r::Integer; constants::Bool=true, ℓ::Integer=maximum(Z), bigInt::Bool=true)
    """
    Z::Array{Int64, 1} the vector of integer group labels
    D::Array{Int64, 1} the vector of degrees
    r::Int64, the largest hyperedge size to compute
    constants::Bool, whether to postmultiply by the combinatorial constants computed by evalConstants prior to returning the values. Choose true if computing modularity directly using the results, choose false if performing incremental updates a la Louvain
    ℓ::Int64, the number of groups (defaults to maximum(Z))
    bigInt::Bool, whether to use bigInt conversions for D and Z. Recommended unless the instance is VERY small.
    return (V::Array{bigInt,1}, μ::Array{bigInt,1}, M::Dict{Vector{Int64},bigInt}, with V and μ as described in computeMoments(), and M::Dict{Vector{Int64},bigInt} is the array of (optionally uncorrected) volume sums.
    """

    if bigInt
        D = convert(Array{BigInt,1}, D)
    end

    V, μ = computeMoments(Z, D, r, ℓ)

    M = bigInt ? Dict{Vector{Int64}, BigInt}() : Dict{Vector{Int64}, Int64}()

    ix = ℓ == 1 ? length(Z) : ℓ

    for i = 1:r, j = 1:minimum([i,ix]), p in Combinatorics.partitions(i, j)
        M[p] = μ[p[end]]*get(M, p[1:(end-1)], 1) - correctOvercounting(M,p)
    end

    if constants
        C = evalConstants(r; bigInt=bigInt, ℓ = ix)
        for p in keys(M)
            M[p] *= C[p]
        end
    end
    return V, μ, M
end

# Extra methods for evalSums for working with degree dicts and the custom hypergraph class.

function evalSums(Z::Vector{<:Integer}, H::hypergraph, ℓ::Integer=maximum(Z), bigInt::Bool=true, constants::Bool = true)
    r = maximum(keys(H.E))
    return evalSums(Z, H.D, r; constants = constants, ℓ = ℓ, bigInt = bigInt)
end

# ------------------------------------------------------------------------------
# INCREMENTS
# ------------------------------------------------------------------------------

function momentIncrements(V::Vector{<:Integer}, μ::Vector{<:Integer}, i::Int64, t::Int64, D::Vector{<:Integer}, Z::Vector{<:Integer})
    """
    Compute the update for the vectors V and μ associated with moving node i from its current group to group t.
    V::Array{Integer, 1}, the vector of group volumes.
    μ::Array{Integer, 1}, the vector of volume-moments.
    i::Int64, the node to move
    t::Int64, the proposed new group for node i
    D::Array{Integer, 1}, the degree vector
    Z::Array{Integer, 1}, the vector of group labels
    returns ΔV, Δμ, the required updates in V and μ
    """

    ΔV = zero(V)
    # update V
    s      = Z[i];
    ΔV[t] += D[i];
    ΔV[s] -= D[i];

    r = length(μ)
    Δμ = [(V[s] + ΔV[s])^j - V[s]^j + (V[t] + ΔV[t])^j - V[t]^j for j in 1:r]

    return ΔV, Δμ
end

function momentIncrements(V::Vector{<:Integer}, μ::Vector{<:Integer}, I_::Vector{<:Integer}, t::Integer, D::Vector{<:Integer}, Z::Vector{<:Integer})
    """
    Compute the update for the vectors V and μ associated with moving all the nodes listed in I from their group to group t.
    THIS CODE ASSUMES THAT THE NODES IN I ARE ALL IN THE SAME GROUP, AND WILL GIVE INCORRECT RESULTS OTHERWISE.
    V::Array{Integer, 1}, the vector of group volumes.
    μ::Array{Integer, 1}, the vector of volume-moments.
    I::Vector{Int64}, the list of multiple nodes to move
    t::Int64, the proposed new group for node i
    D::Array{Integer, 1}, the degree vector
    Z::Array{Integer, 1}, the vector of group labels
    returns ΔV, Δμ, the required updates in V and μ
    """

    ΔV = zero(V)
    # update V
    s  = Z[I_[1]];

    for i ∈ I_
        ΔV[t] += D[i];
        ΔV[s] -= D[i];
    end

    r = length(μ)
    Δμ = [(V[s] + ΔV[s])^j - V[s]^j + (V[t] + ΔV[t])^j - V[t]^j for j in 1:r]

    return ΔV, Δμ
end

function increments(V::Vector{<:Integer}, μ::Vector{<:Integer}, M::Dict{Vector{Int64},<:Integer}, I_::T, t::Integer, D::Vector{<:Integer}, Z::Vector{<:Integer}; bigInt::Bool=true) where T <: Union{<:Integer, Vector{<:Integer}}
    """
    Compute the update for the vectors V and μ, as well as the uncorrected sums M, associated with moving node(s) I from its current group to group t.
    THIS CODE ASSUMES THAT THE NODES IN I ARE ALL IN THE SAME GROUP, AND WILL GIVE INCORRECT RESULTS OTHERWISE.
    V::Array{Integer, 1}, the vector of group volumes.
    μ::Array{Integer, 1}, the vector of volume-moments.
    M::Dict{Array{Int64, 1}, bigInt}, the Dict of uncorrected volume sums.
    i::Int64, the node to move
    t::Int64, the proposed new group for node i
    D::Array{Integer, 1}, the degree vector
    Z::Array{Integer, 1}, the vector of group labels
    returns (ΔV, Δμ, ΔM), the required updates in V, μ, and M
    """

    ΔV, Δμ = momentIncrements(V, μ, I_, t, D, Z)

    # compute increments in M using recursion formula from notes
    if bigInt ΔM = Dict{Vector{Int64}, BigInt}()
    else      ΔM = Dict{Vector{Int64}, Int64}()
    end

    r = maximum([sum(p) for p in keys(M)])

    for i = 1:r, j = 1:i, p in Combinatorics.partitions(i, j)
        ΔM[p] = Δμ[p[end]]*get(M, p[1:(end-1)], 1) + μ[p[end]]*get(ΔM, p[1:(end-1)], 0) + Δμ[p[end]]*get(ΔM, p[1:(end-1)], 0) - correctOvercounting(ΔM,p)
    end
    return(ΔV, Δμ, ΔM)
end

function addIncrements(V::Vector{<:Integer}, μ::Vector{<:Integer}, M::Dict{Vector{Int64},<:Integer}, ΔV::Vector{<:Integer}, Δμ::Vector{<:Integer}, ΔM::Dict{Vector{Int64},<:Integer})
    """
    Add the increments (ΔV, Δμ, ΔM) to (V, μ, M), entrywise in the case of M.
    V::Array{Integer, 1}, the vector of group volumes.
    μ::Array{Integer, 1}, the vector of volume-moments.
    M::Dict{Array{Int64, 1}, bigInt}, the Dict of uncorrected volume sums.
    ΔV, Δμ, ΔM: increments in each of the above quantities returned by increments()
    """
    M̃ = Dict(p => M[p] + ΔM[p] for p in keys(M))
    return(V + ΔV, μ + Δμ, M̃)
end

# ------------------------------------------------------------------------------
# COMPLETE COMPUTATION OF SECOND (VOLUME) TERM IN MODULARITY
# ------------------------------------------------------------------------------

function aggregateSums(M, Ω::IntensityFunction)
    M̂ = Dict()
    for p in keys(M)
        p̂ = Ω.aggregator(p)
        M̂[p̂] = get(M̂, p̂, 0) + M[p]
    end
    return M̂
end

function second_term_eval(H::hypergraph, Z::Vector{<:Integer}, Ω::IntensityFunction; α, ℓ::Int64 = 0, bigInt::Bool=true)
    """
    Naive implementation, computes sums from scratch.
    H::hypergraph
    Z::Array{Int64, 1}, the group label vector.
    ℓ::Int64, number of groups
    Ω: group interation function (e.g., planted partition). Needs to have a mode argument which, when set to value "partition", will cause evaluation on partition vectors rather than label vectors.
    bigInt::Bool, whether to use bigInt conversions. Strongly recommended unless the instance is VERY small.
    """

    obj = 0.0

    if ℓ == 0
        ℓ = maximum(Z)
    end

    V, μ, M = evalSums(Z, H, ℓ, bigInt)
    M̂ = aggregateSums(M, Ω)
    for p̂ in keys(M̂)
        if M̂[p̂] >= 1
            obj += Ω.ω(p̂, α)*M̂[p̂]
        end
    end
    return obj
end

max_stars_repo_name: PhilChodrow/HypergraphModularity

max_issues_repo_name: nveldt/TestPkg

max_forks_repo_name: PhilChodrow/HypergraphModularity


text:
@block KimFung [intermediate] begin

    @cell "simple line plot" [lines, "2d"] begin
        scene = Scene()
        lines!(scene, RNG.rand(10), linewidth=3.0, color=:blue)
        lines!(scene, RNG.rand(10), linewidth=3.0, color=:red)
        lines!(scene, RNG.rand(10), linewidth=3.0, color=:green)
        lines!(scene, RNG.rand(10), linewidth=3.0, color=:purple)
        lines!(scene, RNG.rand(10), linewidth=3.0, color=:orange)
    end

    @cell "multiple functions" [lines, "2d"] begin
        scene = Scene()
        x = range(0, stop = 3pi, step = 0.01)
        lines!(scene, x, sin.(x), color = :black)
        lines!(scene, x, cos.(x), color = :blue)
    end

    @cell "animations" [lines, record, "2d"] begin
        scene = Scene()
        mytime = Node(0.0)
        f(v, t) = sin(v + t)
        g(v, t) = cos(v + t)
        lines!(scene,lift(t -> f.(range(0, stop = 2pi, length = 50), t), mytime),color = :blue)
        lines!(scene,lift(t -> g.(range(0, stop = 2pi, length = 50), t), mytime),color = :orange)
        record(scene, @replace_with_a_path(mp4), range(0, stop = 4pi, length = 100)) do i
            mytime[] = i
        end
    end

    @cell "parametric plots" [poly, "2d"] begin
        x = LinRange(0, 2pi, 100)
        poly(Point2f0.(sin.(x), sin.(2x)), color = "orange", strokecolor = "blue", strokewidth = 4)
    end

    @cell "colors" [lines, "2d"] begin
        lines(1:50, 1:1, linewidth = 20, color = to_colormap(:viridis, 50))
    end

    @cell "attributes" [attributes, lines, color, colormap, linewidth, "2d"] begin
        scene = Scene()
        lines!(scene, RNG.rand(10), color = to_colormap(:viridis, 10), linewidth = 5)
        lines!(scene, RNG.rand(20), color = :red, alpha = 0.5)
    end

    @cell "build plot in pieces" [lines, "2d"] begin
         scene = Scene()
         # initialize the stepper and give it an output destination
         st = Stepper(scene, @replace_with_a_path)
         lines!(scene, RNG.rand(50)/3, color = :purple, linewidth = 5)
         step!(st)
         scatter!(scene, RNG.rand(50), color = :orange, markersize = 1)
         step!(st)
    end

    @cell "barplot" [barplot, "2d"] begin
        barplot(RNG.randn(99))
    end

    @cell "subplots" [lines, scatter, barplot, histogram, vbox, hbox, "2d"] begin
        scene1, scene2, scene3, scene4 = Scene(), Scene(), Scene(), Scene()

        lines!(scene1, RNG.rand(10), color=:red)
        lines!(scene1, RNG.rand(10), color=:blue)
        lines!(scene1, RNG.rand(10), color=:green)

        scatter!(scene2, RNG.rand(10), color=:red)
        scatter!(scene2, RNG.rand(10), color=:blue)
        scatter!(scene2, RNG.rand(10), color=:orange)

        barplot!(scene3, RNG.randn(99))

        v(x::Point2{T}) where T = Point2f0(x[2], 4*x[1])
        streamplot!(scene4, v, -2..2, -2..2)

        vbox(hbox(scene2,scene1),hbox(scene4,scene3))
    end

    @cell "text" [scatter, text, "2d"] begin
        scene = Scene()
        scatter!(scene, RNG.rand(10), color=:red)
        text!(scene,"adding text",textsize = 0.6, position = (5.0, 1.1))
    end

    @cell "contours" [contour, fillrange, "2d"] begin
        x = LinRange(-1, 1, 20)
        y = LinRange(-1, 1, 20)
        z = x .* y'

        vbox(
            contour(x, y, z, levels = 50, linewidth =3),
            contour(x, y, z, levels = 0, linewidth = 0, fillrange = true)
        )
    end

    @cell "3D" [meshscatter, colormap] begin
        x = [2 .* (i/3) .* cos(i) for i in range(0, stop = 4pi, length = 30)]
        y = [2 .* (i/3) .* sin(i) for i in range(0, stop = 4pi, length = 30)]
        z = range(0, stop = 5, length = 30)
        meshscatter(x, y, z, markersize = 0.5, color = to_colormap(:blues, 30))
    end

    @cell "heatmap" [heatmap, "2d"] begin
        x = LinRange(-1, 1, 40)
        y = LinRange(-1, 1, 40)
        z = x .* y'
        heatmap(x, y, z)
    end

    @cell "label rotation, title location" [lines, axisnames, rotation, legend, position, vbox, "2d"] begin
        scene = Scene()
        line = lines!(scene, RNG.rand(10), linewidth=3.0, color=:red)
        scene[Axis][:names, :axisnames] = ("text","")
        scene[Axis][:names, :rotation] = (1.5pi)
        leg = legend([scene[2]], ["line"], position = (1, -1))
        vbox(line, leg)
    end

    @cell "lines with varying colors" [lines, scatter, vbox, colormap, "2d"] begin
        t = range(0, stop=1, length=100)
        θ = (6π) .* t
        x = t .* cos.(θ)
        y = t .* sin.(θ)
        p1 = lines(x,y,color=t,colormap=:colorwheel,linewidth= 8,scale_plot= false)
        p2 = scatter(x,y,color=t,colormap=:colorwheel,linewidth= 8,scale_plot= false)
        vbox(p1, p2)
    end

end


text:
using Configurations
using YAML
using Logging
using Dates

@option struct ConsoleLoggerConfig
    is_expand_stack_trace::Bool = true
end

@option struct RotatingLoggerConfig
    path::String = "logs"
    file_format::String = raw"YYYY-mm-dd.\l\o\g"
end

@option struct DriverLoggerConfig
    console_logger::Union{ConsoleLoggerConfig, Nothing}=ConsoleLoggerConfig()
    rotating_logger::Union{RotatingLoggerConfig, Nothing}=RotatingLoggerConfig()
end

@option struct LokiLoggerConfig
    url::String = "http://127.0.0.1:3100"
end

@option struct LoggingConfig
    # filter
    log_level::String = "Info"
    # transformer
    date_format::String="yyyy-mm-ddTHH:MM:SS.s"
    # sink
    driver_logger::Union{DriverLoggerConfig, Nothing} = DriverLoggerConfig()
    loki_logger::Union{LokiLoggerConfig, Nothing} = nothing
end

@option struct Config
    banner::Bool = Base.JLOptions().banner != 0
    color::Bool = Base.have_color
    logging::LoggingConfig = LoggingConfig()
end


text:
include("train.jl")

using UnicodePlots

# get the data

# characters: [^0-9a-zA-Z&:,./()[]_-] THIS IS OLD
raw = []
open("case_names/data.txt") do f

    line = 0
    while !eof(f)
        push!(raw, readline(f))
    end # while

end # do

#println(raw[1:20])

#chars = "ABCDEFGHIJKLMNOPQRSTUVWXYZ -[]()0123456789"
chars = "ABCDEFGHIJKLMNOPQRSTUVWXYZ -"
raw = uppercase.(raw)
binary_split = [[[0 #==1/(length(chars)*2)==# for _ in 1:(length(chars)+1)] for i = 1:length(x)+1] for x in raw]
for s = 1:length(binary_split)

    for c = 1:length(raw[s])
        binary_split[s][c][findfirst(raw[s][c], chars)] = 1
    end # for
    binary_split[s][end][end] = 1

end # for

input_chars = 3
# split the cases to serve as inputs for the RNN.
DMatrix = Vector{Float64}[]
Y = Vector{Float64}[]

for s in binary_split

    for c = 1:length(s) -1

        push!(DMatrix, [])
        push!(Y, s[c+1])

        if c < input_chars
            for i = 1:(input_chars - c)
                append!( DMatrix[end], zeros(length(chars)) )
            end # for
        end # if
        if c != 0
            for i = (max(c - input_chars + 1, 1)):c
                append!( DMatrix[end], s[c][1:end-1] )
            end # for
        end # if

    end # for

end # for

input_dim = (length(chars)) * input_chars
hidden_dim = (length(chars)+1) * 2
output_dim = (length(chars)+1)
dims = [hidden_dim for i in 1:3]
prepend!(dims, input_dim)

append!(dims, output_dim)

println(string("Dimensions: ", dims))

nn_results = train_network(dims, reverse(DMatrix)[1:500], reverse(Y)[1:500], chars, epochs=50, η = 0.001)

max_stars_repo_name: Squalm/NeuralNetworks

max_issues_repo_name: Squalm/NeuralNetworks

max_forks_repo_name: Squalm/NeuralNetworks


text:
import Pkg; Pkg.add("BinaryBuilder")
using BinaryBuilder

# Collection of sources required to build OpenBLAS
name = "OpenBLAS"
version = v"0.3.13"
sources = [
    "https://github.com/xianyi/OpenBLAS/archive/v0.3.13.tar.gz" =>
    "d2b11c47774b9216660e76e2fc67e87079f26fa1",
    "./bundled",
]

# Bash recipe for building across all platforms
script = raw"""
# We always want threading
flags=(USE_THREAD=1 GEMM_MULTITHREADING_THRESHOLD=50 NO_AFFINITY=1)

# We are cross-compiling
flags+=(CROSS=1 "HOSTCC=$CC_FOR_BUILD" PREFIX=/ "CROSS_SUFFIX=${target}-")

# We need to use our basic objconv, not a prefixed one:
flags+=(OBJCONV=objconv)

# Set BINARY=64 on x86_64 platforms (but not AArch64 or powerpc64le)
if [[ ${target} == x86_64-* ]]; then
    flags+=(BINARY=64)
fi

# On Intel architectures, engage DYNAMIC_ARCH
if [[ ${proc_family} == intel ]]; then
    flags+=(DYNAMIC_ARCH=1)
# Otherwise, engage a specific target
elif [[ ${target} == aarch64-* ]]; then
    flags+=(TARGET=ARMV8)
elif [[ ${target} == arm-* ]]; then
    flags+=(TARGET=ARMV7)
elif [[ ${target} == powerpc64le-* ]]; then
    flags+=(TARGET=POWER8)
fi

# Enter the fun zone
cd ${WORKSPACE}/srcdir/OpenBLAS-*/

# Patch so that our LDFLAGS make it all the way through
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/osx_exports_ldflags.patch"

# Build the library
make "${flags[@]}" -j${nproc}

# Install the library
make "${flags[@]}" "PREFIX=$prefix" install

# Force the library to be named the same as in Julia-land.
# Move things around, fix symlinks, and update install names/SONAMEs.
ls -la ${prefix}/lib
for f in ${prefix}/lib/libopenblas*p-r0*; do
    name=${LIBPREFIX}.0.${f#*.}

    # Move this file to a julia-compatible name
    mv -v ${f} ${prefix}/lib/${name}

    # If there were links that are now broken, fix 'em up
    for l in $(find ${prefix}/lib -xtype l); do
        if [[ $(basename $(readlink ${l})) == $(basename ${f}) ]]; then
            ln -vsf ${name} ${l}
        fi
    done

    # If this file was a .so or .dylib, set its SONAME/install name
    if [[ ${f} == *.so.* ]] || [[ ${f} == *.dylib ]]; then
        if [[ ${target} == *linux* ]] || [[ ${target} == *freebsd* ]]; then
            patchelf --set-soname ${name} ${prefix}/lib/${name}
        elif [[ ${target} == *apple* ]]; then
            install_name_tool -id ${name} ${prefix}/lib/${name}
        fi
    fi
done
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.
platforms = supported_platforms()

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: scienceasdf/OpenBLASBuilder

max_issues_repo_name: scienceasdf/OpenBLASBuilder

max_forks_repo_name: scienceasdf/OpenBLASBuilder


text:
export top,
    skip,
    mergex,
    mergeleftx,
    postfix

using Tables

"""
Take top `n` elements from an iterator
"""
top(n) = rows -> Iterators.take(rows, n)

Tables.istable(t::Iterators.Take{I}) where {I} = true
Tables.rowaccess(t::Iterators.Take{I}) where {I} = true
Tables.schema(t::Iterators.Take{I}) where {I} = Tables.schema(t.xs) #Tables.Schema(eltype(I))

"""
skip `n` elements from an iterator
"""
skip(n) = rows -> Iterators.drop(rows, n)

Tables.istable(t::Iterators.Drop{I}) where {I} = true
Tables.rowaccess(t::Iterators.Drop{I}) where {I} = true
Tables.schema(t::Iterators.Drop{I}) where {I} = Tables.Schema(eltype(I))

"""
append symbol `s` to `symbols` postfixing it with `postfix` until unique
"""
# should we call it push_uniq!  ?
function mergesymbol(symbols::Vector{Symbol}, s::Symbol; postfix = "x")
    if s in symbols
	mergesymbol(symbols, Symbol(string(s)*postfix))
    else
	push!(symbols, s)
    end
    symbols
end

# merge_uniq better name?
@generated function mergex(x::NamedTuple{lN,lT}, y::NamedTuple{rN,rT}) where {lN,lT,rN,rT}
    N = Symbol[lN...]
    for s in rN
	N = mergesymbol(N, s)
    end
    NT = Tuple(N)
    T = Tuple{lT.parameters..., rT.parameters...}
    :(NamedTuple{$NT,$T}((x..., y...)))
end

@generated function postfix(x::NamedTuple{N,T}, ::Type{Val{K}}) where {N,T,K}
    n = length(N)
    newN = Vector{Symbol}(undef, n)
    for i = 1:n
	newN[i] = Symbol(string(N[i])*string(K))
    end
    :(NamedTuple{$(Tuple(newN)),$T}((x...,)))
end

postfix(name::String) = rows -> (postfix(x, Val{Symbol(name)}) for x in rows)

# for leftjoin

# like mergex, but convert coltypes of right stream y to Union{Missing, coltype}
# merge_maybe better name?
@generated function mergeleftx(x::NamedTuple{lN,lT}, y::NamedTuple{rN,rT}) where {lN,lT,rN,rT}
    N = Symbol[lN...]
    for s in rN
	N = mergesymbol(N, s)
    end
    NT = Tuple(N)
    T = Tuple{lT.parameters..., (Union{Missing, t} for t in rT.parameters)...}
    :(NamedTuple{$NT,$T}((x..., y...)))
end

# like mergex, but insert all columns from right stream as missing
@generated function mergeleftx(x::NamedTuple{lN,lT}, ::Type{NamedTuple{rN,rT}}) where {lN,lT,rN,rT}
    N = Symbol[lN...]
    for s in rN
	N = mergesymbol(N, s)
    end
    NT = Tuple(N)
    T = Tuple{lT.parameters..., (Union{Missing, t} for t in rT.parameters)...}
    y = repeat([missing], length(rT.parameters))
    :(NamedTuple{$NT,$T}((x..., $(y...))))

end

# do not export this
# should be called as GorJulia.pkgpath
pkgpath(paths...) = joinpath(@__DIR__, "..", paths...)


text:
module ApproxFunFourier
using Base, LinearAlgebra, Reexport, AbstractFFTs, FFTW, InfiniteArrays, FastTransforms, IntervalSets,
            DomainSets

@reexport using ApproxFunBase

import AbstractFFTs: Plan, fft, ifft
import FFTW: plan_r2r!, fftwNumber, REDFT10, REDFT01, REDFT00, RODFT00, R2HC, HC2R,
                r2r!, r2r,  plan_fft, plan_ifft, plan_ifft!, plan_fft!

import ApproxFunBase: normalize!, flipsign, FiniteRange, Fun, MatrixFun, UnsetSpace, VFun, RowVector,
                UnivariateSpace, AmbiguousSpace, SumSpace, SubSpace, WeightSpace, NoSpace, Space,
                HeavisideSpace, PointSpace,
                IntervalOrSegment, RaggedMatrix, AlmostBandedMatrix,
                AnyDomain, ZeroSpace, ArraySpace, TrivialInterlacer, BlockInterlacer,
                AbstractTransformPlan, TransformPlan, ITransformPlan,
                ConcreteConversion, ConcreteMultiplication, ConcreteDerivative, ConcreteIntegral, CalculusOperator,
                ConcreteVolterra, Volterra, VolterraWrapper,
                MultiplicationWrapper, ConversionWrapper, DerivativeWrapper, Evaluation, EvaluationWrapper,
                Conversion, defaultConversion, defaultcoefficients, default_Fun, Multiplication, Derivative, Integral, bandwidths,
                ConcreteEvaluation, ConcreteDefiniteLineIntegral, ConcreteDefiniteIntegral, ConcreteIntegral,
                DefiniteLineIntegral, DefiniteIntegral, ConcreteDefiniteIntegral, ConcreteDefiniteLineIntegral, IntegralWrapper,
                ReverseOrientation, ReverseOrientationWrapper, ReverseWrapper, Reverse, NegateEven,
                Dirichlet, ConcreteDirichlet, DirichletWrapper,
                TridiagonalOperator, SubOperator, Space, @containsconstants, spacescompatible,
                hasfasttransform, canonicalspace, domain, setdomain, prectype, domainscompatible,
                plan_transform, plan_itransform, plan_transform!, plan_itransform!, transform, itransform, hasfasttransform,
                CanonicalTransformPlan, ICanonicalTransformPlan,
                Integral,
                domainspace, rangespace,
                union_rule, conversion_rule, maxspace_rule, conversion_type, maxspace, hasconversion, points,
                rdirichlet, ldirichlet, lneumann, rneumann, ivp, bvp,
                linesum, differentiate, integrate, linebilinearform, bilinearform,
                UnsetNumber, coefficienttimes, subspace_coefficients, sumspacecoefficients, specialfunctionnormalizationpoint,
                Segment, IntervalOrSegmentDomain, PiecewiseSegment, isambiguous, Vec, eps, isperiodic,
                arclength, complexlength,
                invfromcanonicalD, fromcanonical, tocanonical, fromcanonicalD, tocanonicalD, canonicaldomain, setcanonicaldomain, mappoint,
                reverseorientation, checkpoints, evaluate, mul_coefficients, coefficients, coefficientmatrix, isconvertible,
                clenshaw, ClenshawPlan, sineshaw,
                toeplitz_getindex, toeplitz_axpy!, sym_toeplitz_axpy!, hankel_axpy!, ToeplitzOperator, SymToeplitzOperator, hankel_getindex,
                SpaceOperator, ZeroOperator, InterlaceOperator,
                interlace!, reverseeven!, negateeven!, cfstype, pad!, alternatesign!, mobius,
                extremal_args, hesseneigvals, chebyshev_clenshaw, recA, recB, recC, roots,splitatroots,
                chebmult_getindex, intpow, alternatingsum,
                domaintype, diagindshift, rangetype, weight, isapproxinteger, default_Dirichlet, scal!, dotu,
                components, promoterangespace, promotedomainspace, choosedomainspace,
                block, blockstart, blockstop, blocklengths, isblockbanded, pointscompatible,
                AbstractProductSpace, MultivariateFun, BivariateSpace,
                @wrapperstructure, @wrapperspaces, @wrapper, @calculus_operator, resizedata!, slnorm, affine_setdiff,
                complexroots, EmptyDomain, RectDomain

import DomainSets: Domain, indomain, UnionDomain, ProductDomain, FullSpace, Point, elements, DifferenceDomain,
            Interval, ChebyshevInterval, boundary, ∂, rightendpoint, leftendpoint,
            dimension

import Base: values, convert, getindex, setindex!, *, +, -, ==, <, <=, >, |, !, !=, eltype, iterate,
                >=, /, ^, \, ∪, transpose, size, tail, broadcast, broadcast!, copyto!, copy, to_index, (:),
                similar, map, vcat, hcat, hvcat, show, summary, stride, sum, cumsum, sign, imag, conj, inv,
                complex, reverse, exp, sqrt, abs, abs2, sign, issubset, values, in, first, last, rand, intersect, setdiff,
                isless, union, angle, join, isnan, isapprox, isempty, sort, merge, promote_rule,
                minimum, maximum, extrema, argmax, argmin, findmax, findmin, isfinite,
                zeros, zero, one, promote_rule, repeat, length, resize!, isinf,
                getproperty, findfirst, unsafe_getindex, fld, cld, div, real, imag,
                @_inline_meta, eachindex, firstindex, lastindex, keys, isreal, OneTo,
                Array, Vector, Matrix, view, ones, @propagate_inbounds, print_array,
                split

import LinearAlgebra: BlasInt, BlasFloat, norm, ldiv!, mul!, det, eigvals, dot, cross,
                qr, qr!, rank, isdiag, istril, istriu, issymmetric, ishermitian,
                Tridiagonal, diagm, diagm_container, factorize, nullspace,
                Hermitian, Symmetric, adjoint, transpose, char_uplo

import InfiniteArrays: InfRanges, AbstractInfUnitRange, OneToInf

import FastTransforms: ChebyshevTransformPlan, IChebyshevTransformPlan, plan_chebyshevtransform,
                        plan_chebyshevtransform!, plan_ichebyshevtransform, plan_ichebyshevtransform!

export Fourier, Taylor, Hardy, CosSpace, SinSpace, Laurent, PeriodicDomain

include("utils.jl")
include("Domains/Domains.jl")

for T in (:CosSpace,:SinSpace)
    @eval begin
        struct $T{D<:PeriodicDomain,R} <: Space{D,R}
            domain::D
            $T{D,R}(d::PeriodicDomain) where {D,R} = new(convert(D,d))
            $T{D,R}(d::D) where{D,R} = new(d)
        end
        $T(d::PeriodicDomain) = $T{typeof(d),real(prectype(d))}(d)
        $T(d) = $T(convert(PeriodicDomain, d))
        $T() = $T(PeriodicSegment())
        spacescompatible(a::$T,b::$T) = domainscompatible(a,b)
        hasfasttransform(::$T) = true
        canonicalspace(S::$T) = Fourier(domain(S))
        setdomain(S::$T,d::PeriodicDomain) = $T(d)
    end
end

"""
`CosSpace()` is the space spanned by `[1,cos θ,cos 2θ,...]`
"""
CosSpace()

"""
`SinSpace()` is the space spanned by `[sin θ,sin 2θ,...]`
"""
SinSpace()

# s == true means analytic inside, taylor series
# s == false means anlytic outside and decaying at infinity
"""
`Hardy{false}()` is the space spanned by `[1/z,1/z^2,...]`.
`Hardy{true}()` is the space spanned by `[1,z,z^2,...]`.
"""
struct Hardy{s,D<:PeriodicDomain,R} <: Space{D,R}
    domain::D
    Hardy{s,D,R}(d) where {s,D,R} = new{s,D,R}(d)
    Hardy{s,D,R}() where {s,D,R}  = new{s,D,R}(D())
end

# The <: PeriodicDomain is crucial for matching Basecall overrides
"""
`Taylor()` is the space spanned by `[1,z,z^2,...]`.
This is a type alias for `Hardy{true}`.
"""
const Taylor{D<:PeriodicDomain,R} = Hardy{true,D,R}

@containsconstants CosSpace
@containsconstants Taylor

Base.promote_rule(::Type{Fun{S,V}},::Type{T}) where {T<:Number,S<:Union{Hardy{true},CosSpace},V} =
    Fun{S,promote_type(V,T)}
Base.promote_rule(::Type{Fun{S}},::Type{T}) where {T<:Number,S<:Union{Hardy{true},CosSpace}} =
    Fun{S,T}

(H::Type{Hardy{s}})(d::PeriodicDomain) where {s} = Hardy{s,typeof(d),complex(prectype(d))}(d)
(H::Type{Hardy{s}})() where {s} = Hardy{s}(Circle())

canonicalspace(S::Hardy) = S
setdomain(S::Hardy{s},d::PeriodicDomain) where {s} = Hardy{s}(d)

spacescompatible(a::Hardy{s},b::Hardy{s}) where {s} = domainscompatible(a,b)
hasfasttransform(::Hardy) = true

for (Typ,Plfft!,Plfft,Pltr!,Pltr) in ((:TransformPlan,:plan_fft!,:plan_fft,:plan_transform!,:plan_transform),
                           (:ITransformPlan,:plan_ifft!,:plan_ifft,:plan_itransform!,:plan_itransform))
    @eval begin
        $Pltr!(sp::Hardy,x::AbstractVector{T}) where {T<:Complex} = $Typ(sp,$Plfft!(x),Val{true})
        $Pltr!(::Hardy,x::AbstractVector{T}) where {T<:Real} =
            error("In place variants not possible with real data.")

        $Pltr(sp::Hardy,x::AbstractVector{T}) where {T<:Complex} = $Typ(sp,$Pltr!(sp,x),Val{false})
        function $Pltr(sp::Hardy,x::AbstractVector{T}) where T
            plan = $Pltr(sp,Array{Complex{T}}(undef,length(x))) # we can reuse vector in itransform
            $Typ{T,typeof(sp),false,typeof(plan)}(sp,plan)
        end

        *(P::$Typ{T,HS,false},vals::AbstractVector{T}) where {T<:Complex,HS<:Hardy} = P.plan*copy(vals)
        *(P::$Typ{T,HS,false},vals::AbstractVector{T}) where {T,HS<:Hardy} = P.plan*Vector{Complex{T}}(vals)
    end
end

*(P::TransformPlan{T,Hardy{true,DD,RR},true},vals::AbstractVector{T}) where {T,DD,RR} =
    lmul!(one(T)/length(vals),P.plan*vals)
*(P::ITransformPlan{T,Hardy{true,DD,RR},true},cfs::AbstractVector{T}) where {T,DD,RR} =
    lmul!(length(cfs),P.plan*cfs)
*(P::TransformPlan{T,Hardy{false,DD,RR},true},vals::AbstractVector{T}) where {T,DD,RR} =
    lmul!(one(T)/length(vals),reverse!(P.plan*vals))
*(P::ITransformPlan{T,Hardy{false,DD,RR},true},cfs::AbstractVector{T}) where {T,DD,RR} =
    lmul!(length(cfs),P.plan*reverse!(cfs))

transform(sp::Hardy,vals::AbstractVector,plan) = plan*vals
itransform(sp::Hardy,vals::AbstractVector,plan) = plan*vals

evaluate(f::AbstractVector,S::Taylor{D,R},z) where {D<:PeriodicDomain,R} = horner(f,fromcanonical(Circle(),tocanonical(S,z)))
function evaluate(f::AbstractVector,S::Taylor{D,R},z) where {D<:Circle,R}
    z=mappoint(S,Circle(),z)
    d=domain(S)
    horner(f,z)
end

function evaluate(f::AbstractVector,S::Hardy{false,D},z) where D<:PeriodicDomain
    z=mappoint(S,Circle(),z)
    z=1/z
    z*horner(f,z)
end
function evaluate(f::AbstractVector,S::Hardy{false,D},z) where D<:Circle
    z=mappoint(S,Circle(),z)
    z=1/z
    z*horner(f,z)
end

##TODO: fast routine

function horner(c::AbstractVector,kr::AbstractRange{Int},x)
    T = promote_type(eltype(c),eltype(x))
    if isempty(c)
        return zero(x)
    end

    ret = zero(T)
    @inbounds for k in reverse(kr)
        ret = muladd(x,ret,c[k])
    end

    ret
end

function horner(c::AbstractVector,kr::AbstractRange{Int},x::AbstractVector)
    n,T = length(x),promote_type(eltype(c),eltype(x))
    if isempty(c)
        return zero(x)
    end

    ret = zeros(T,n)
    @inbounds for k in reverse(kr)
        ck = c[k]
        @simd for i = 1:n
            ret[i] = muladd(x[i],ret[i],ck)
        end
    end

    ret
end

horner(c::AbstractVector,x) = horner(c,1:length(c),x)
horner(c::AbstractVector,x::AbstractArray) = horner(c,1:length(c),x)
horner(c::AbstractVector,kr::AbstractRange{Int},x::AbstractArray) = reshape(horner(c,kr,vec(x)),size(x))

## Cos and Sin space

points(sp::CosSpace,n) = points(domain(sp),2n-2)[1:n]

plan_transform(::CosSpace,x::AbstractVector) = plan_chebyshevtransform(x, Val(2))
plan_itransform(::CosSpace,x::AbstractVector) = plan_ichebyshevtransform(x, Val(2))
transform(::CosSpace,vals,plan) = plan*vals
itransform(::CosSpace,cfs,plan) = plan*cfs

clenshaw(::CosSpace, c::AbstractVector, x) = chebyshev_clenshaw(c, x)

clenshaw(sp::CosSpace, c::AbstractVector, x::AbstractArray) =
    clenshaw(c,x,ClenshawPlan(promote_type(eltype(c),eltype(x)),sp,length(c),length(x)))

clenshaw(c::AbstractVector, x::Vector, plan::ClenshawPlan{<:CosSpace}) = chebyshev_clenshaw(c,x,plan)

function evaluate(f::AbstractVector,S::CosSpace,t)
    if t ∈ domain(S)
        clenshaw(S,f,cos(tocanonical(S,t)))
    else
        zero(cfstype(f))
    end
end

points(sp::SinSpace,n)=points(domain(sp),2n+2)[2:n+1]

for (Typ,Pltr!,Pltr) in ((:TransformPlan,:plan_transform!,:plan_transform),
                         (:ITransformPlan,:plan_itransform!,:plan_itransform))
    @eval begin
        $Pltr!(sp::SinSpace{DD},x::AbstractVector{T}) where {DD,T<:fftwNumber} =
            $Typ(sp,plan_r2r!(x,RODFT00),Val{true})
        $Pltr(sp::SinSpace{DD},x::AbstractVector{T}) where {DD,T<:fftwNumber} =
            $Typ(sp,$Pltr!(sp,x),Val{false})
        $Pltr!(sp::SinSpace{DD},x::AbstractVector{T}) where {DD,T} =
            error("transform for SinSpace only implemented for fftwNumbers")
        $Pltr(sp::SinSpace{DD},x::AbstractVector{T}) where {DD,T} =
            error("transform for SinSpace only implemented for fftwNumbers")

        *(P::$Typ{T,SinSpace{D,R},false},vals::AbstractVector{T}) where {T,D,R} = P.plan*copy(vals)
    end
end

*(P::TransformPlan{T,SinSpace{DD,RR},true},vals::AbstractVector{T}) where {T,DD,RR} =
    lmul!(one(T)/(length(vals)+1),P.plan*vals)
*(P::ITransformPlan{T,SinSpace{DD,RR},true},cfs::AbstractVector{T}) where {T,DD,RR} =
    lmul!(one(T)/2,P.plan*cfs)

transform(sp::SinSpace,vals::AbstractVector,plan) = plan*vals
itransform(sp::SinSpace,vals::AbstractVector,plan) = plan*vals

evaluate(f::AbstractVector,S::SinSpace,t) = sineshaw(f,tocanonical(S,t))

## Laurent space
"""
`Laurent()` is the space spanned by the complex exponentials
```
    1,exp(-im*θ),exp(im*θ),exp(-2im*θ),…
```
See also `Fourier`.
"""
const Laurent{DD,RR} = SumSpace{Tuple{Hardy{true,DD,RR},Hardy{false,DD,RR}},DD,RR}

##FFT That interlaces coefficients

plan_transform!(sp::Laurent{DD,RR},x::AbstractVector{T}) where {DD,RR,T<:Complex} =
    TransformPlan(sp,plan_fft!(x),Val{true})
plan_itransform!(sp::Laurent{DD,RR},x::AbstractVector{T}) where {DD,RR,T<:Complex} =
    ITransformPlan(sp,plan_ifft!(x),Val{true})

plan_transform!(sp::Laurent{DD,RR},x::AbstractVector{T}) where {DD,RR,T} =
    error("In place variants not possible with real data.")
plan_itransform!(sp::Laurent{DD,RR},x::AbstractVector{T}) where {DD,RR,T} =
    error("In place variants not possible with real data.")

plan_transform(sp::Laurent{DD,RR},x::AbstractVector{T}) where {T<:Complex,DD,RR} =
    TransformPlan(sp,plan_transform!(sp,x),Val{false})
plan_itransform(sp::Laurent{DD,RR},x::AbstractVector{T}) where {T<:Complex,DD,RR} =
    ITransformPlan(sp,plan_itransform!(sp,x),Val{false})

function plan_transform(sp::Laurent{DD,RR},x::AbstractVector{T}) where {T,DD,RR}
    plan = plan_transform(sp,Array{Complex{T}}(undef, length(x))) # we can reuse vector in itransform
    TransformPlan{T,typeof(sp),false,typeof(plan)}(sp,plan)
end
function plan_itransform(sp::Laurent{DD,RR},x::AbstractVector{T}) where {T,DD,RR}
    plan = plan_itransform(sp,Array{Complex{T}}(undef, length(x))) # we can reuse vector in itransform
    ITransformPlan{T,typeof(sp),false,typeof(plan)}(sp,plan)
end

function *(P::TransformPlan{T,Laurent{DD,RR},true},vals::AbstractVector{T}) where {T,DD,RR}
    n = length(vals)
    vals = lmul!(inv(convert(T,n)),P.plan*vals)
    reverseeven!(interlace!(vals,1))
end

function *(P::ITransformPlan{T,Laurent{DD,RR},true},cfs::AbstractVector{T}) where {T,DD,RR}
    n = length(cfs)
    reverseeven!(cfs)
    cfs[:]=[cfs[1:2:end];cfs[2:2:end]]  # TODO: deinterlace!
    lmul!(n,cfs)
    P.plan*cfs
end

*(P::TransformPlan{T,Laurent{DD,RR},false},vals::AbstractVector{T}) where {T<:Complex,DD,RR} =
    P.plan*copy(vals)
*(P::TransformPlan{T,Laurent{DD,RR},false},vals::AbstractVector{T}) where {T,DD,RR} =
    P.plan*Vector{Complex{T}}(vals)
*(P::ITransformPlan{T,Laurent{DD,RR},false},vals::AbstractVector{T}) where {T<:Complex,DD,RR} =
    P.plan*copy(vals)
*(P::ITransformPlan{T,Laurent{DD,RR},false},vals::AbstractVector{T}) where {T,DD,RR} =
    P.plan*Vector{Complex{T}}(vals)

transform(::Laurent{DD,RR},vals,plan) where {DD,RR} = plan*vals
itransform(::Laurent{DD,RR},cfs,plan) where {DD,RR} = plan*cfs

transform(sp::Laurent{DD,RR},vals::AbstractVector) where {DD,RR} = plan_transform(sp,vals)*vals
itransform(sp::Laurent{DD,RR},cfs::AbstractVector) where {DD,RR} = plan_itransform(sp,cfs)*cfs

function evaluate(f::AbstractVector,S::Laurent{DD,RR},z) where {DD,RR}
    z = mappoint(domain(S),Circle(),z)
    invz = 1/z
    horner(f,1:2:length(f),z) + horner(f,2:2:length(f),invz).*invz
end

function Base.conj(f::Fun{Laurent{DD,RR}}) where {DD,RR}
    ncoefficients(f) == 0 && return f

    cfs = Array{cfstype(f)}(undef, iseven(ncoefficients(f)) ? ncoefficients(f)+1 : ncoefficients(f))
    cfs[1] = conj(f.coefficients[1])
    cfs[ncoefficients(f)] = 0
    for k=2:2:ncoefficients(f)-1
        cfs[k] = conj(f.coefficients[k+1])
    end
    for k=3:2:ncoefficients(f)+1
        cfs[k] = conj(f.coefficients[k-1])
    end
    Fun(space(f),cfs)
end

## Fourier space

"""
`Fourier()` is the space spanned by the trigonemtric polynomials
```
    1,sin(θ),cos(θ),sin(2θ),cos(2θ),…
```
See also `Laurent`.
"""
const Fourier{DD,RR} = SumSpace{Tuple{CosSpace{DD,RR},SinSpace{DD,RR}},DD,RR}

Laurent(d::PeriodicDomain) = Laurent{typeof(d),complex(prectype(d))}(d)
Fourier(d::PeriodicDomain) = Fourier{typeof(d),real(prectype(d))}(d)

for Typ in (:Laurent,:Fourier)
    @eval begin
        $Typ() = $Typ(PeriodicSegment())
        $Typ(d) = $Typ(convert(PeriodicDomain, d))

        hasfasttransform(::$Typ{D,R}) where {D,R} = true
    end
end

Laurent(S::Fourier{DD,RR}) where {DD,RR} = Laurent(domain(S))
Fourier(S::Laurent{DD,RR}) where {DD,RR} = Fourier(domain(S))

for T in (:CosSpace,:SinSpace)
    @eval begin
        # override default as canonicalspace must be implemented
        maxspace(::$T,::Fourier{D,R}) where {D,R} = NoSpace()
        maxspace(::Fourier{D,R},::$T) where {D,R} = NoSpace()
    end
end

points(sp::Fourier{D,R},n) where {D,R}=points(domain(sp),n)

struct IFourierTransformPlan{T,SP,PL} <: AbstractTransformPlan{T}
    space::SP
    plan::PL
    work::Vector{T}
end

plan_transform!(sp::Fourier{D,R},x::AbstractVector{T}) where {T<:fftwNumber,D,R} =
    TransformPlan(sp,plan_r2r!(x, R2HC),Val{true})
plan_itransform!(sp::Fourier{D,R},x::AbstractVector{T}) where {T<:fftwNumber,D,R} =
    IFourierTransformPlan(sp, plan_r2r!(x, HC2R), copy(x)) # copy(x) is workspace data

for (Typ,Pltr!,Pltr) in ((:TransformPlan,:plan_transform!,:plan_transform),
                         (:ITransformPlan,:plan_itransform!,:plan_itransform))
    @eval begin
        $Pltr(sp::Fourier{DD,RR},x::AbstractVector{T}) where {T<:fftwNumber,DD,RR} =
            $Typ(sp,$Pltr!(sp,x),Val{false})
        $Pltr!(sp::Fourier{DD,RR},x::AbstractVector{T}) where {T,DD,RR} =
            error("transform for Fourier only implemented for fftwNumbers")
        $Pltr(sp::Fourier{DD,RR},x::AbstractVector{T}) where {T,DD,RR} =
            error("transform for Fourier only implemented for fftwNumbers")

        *(P::$Typ{T,Fourier{DD,RR},false},vals::AbstractVector{T}) where {T,DD,RR} = P.plan*copy(vals)
        *(P::$Typ{T,Fourier{DD,RR},false},vals::AbstractVector{<:Complex}) where {T,DD,RR} =
                P.plan*real(vals) .+ im.*(P.plan*imag(vals))
    end
end

mul!(cfs::AbstractVector{T}, P::TransformPlan{T,Fourier{DD,RR},true}, vals::AbstractVector{T}) where {T,DD,RR} =
    P*copyto!(cfs, vals)

mul!(cfs::AbstractVector{T}, P::TransformPlan{T,Fourier{DD,RR},false}, vals::AbstractVector{T}) where {T,DD,RR} =
    P.plan*copyto!(cfs, vals)

function *(P::TransformPlan{T,Fourier{DD,RR},true},vals::AbstractVector{T}) where {T,DD,RR}
    n = length(vals)
    cfs = lmul!(convert(T,2)/n,P.plan*vals)
    cfs[1] /= 2
    if iseven(n)
        cfs[n÷2+1] /= 2
    end

    negateeven!(reverseeven!(interlace!(cfs,1)))
end

mul!(vals::AbstractVector{T}, P::IFourierTransformPlan{T,Fourier{DD,RR}}, cfs::AbstractVector{T}) where {T,DD,RR} =
    P*copyto!(vals, cfs)

mul!(vals::AbstractVector{T}, P::ITransformPlan{T,Fourier{DD,RR},false},cfs::AbstractVector{T}) where {T,DD,RR} =
    mul!(vals, P.plan, cfs)

function *(P::IFourierTransformPlan{T,Fourier{DD,RR}},cfs::AbstractVector{T}) where {T,DD,RR}
    n = length(cfs)
    reverseeven!(negateeven!(cfs))

    # allocation free version of
    #    cfs[:] = [view(cfs,1:2:n); view(cfs, 2:2:n)]
    P.work .= cfs
    j = 1
    for k=1:2:n
        cfs[j] = P.work[k]
        j+=1
    end
    for k=2:2:n
        cfs[j] = P.work[k]
        j+=1
    end

    if iseven(n)
        cfs[n÷2+1] *= 2
    end
    cfs[1] *= 2
    P.plan*lmul!(inv(convert(T,2)),cfs)
end

transform(sp::Fourier{DD,RR},vals::AbstractVector,plan) where {DD,RR} = plan*vals
itransform(sp::Fourier{DD,RR},cfs::AbstractVector,plan) where {DD,RR} = plan*cfs

transform(sp::Fourier{DD,RR},vals::AbstractVector) where {DD,RR} = plan_transform(sp,vals)*vals
itransform(sp::Fourier{DD,RR},cfs::AbstractVector) where {DD,RR} = plan_itransform(sp,cfs)*cfs

canonicalspace(S::Laurent{DD,RR}) where {DD<:PeriodicSegment,RR} = Fourier(domain(S))
canonicalspace(S::Fourier{DD,RR}) where {DD<:Circle,RR} = Laurent(domain(S))
canonicalspace(S::Laurent{DD,RR}) where {DD<:PeriodicLine,RR} = S

## Ones and zeros

for sp in (:Fourier,:CosSpace,:Laurent,:Taylor)
    @eval begin
        one(::Type{T},S::$sp{DD,RR}) where {T<:Number,DD,RR} = Fun(S,fill(one(T),1))
        one(S::$sp{DD,RR}) where {DD,RR} = Fun(S,fill(1.0,1))
    end
end

function Fun(::typeof(identity), S::Taylor{DD,RR}) where {DD<:Circle,RR}
    d=domain(S)
    if d.orientation
        Fun(S,[d.center,d.radius])
    else
        error("Cannot create identity on $S")
    end
end

Fun(::typeof(identity), S::Fourier{DD,RR}) where {DD<:Circle,RR} =
    Fun(Fun(identity, Laurent(domain(S))),S)

reverseorientation(f::Fun{Fourier{DD,RR}}) where {DD,RR} =
    Fun(Fourier(reverseorientation(domain(f))),alternatesign!(copy(f.coefficients)))
function reverseorientation(f::Fun{Laurent{DD,RR}}) where {DD,RR}
    # exp(im*k*x) -> exp(-im*k*x), or equivalentaly z -> 1/z
    n=ncoefficients(f)
    ret=Array{cfstype(f)}(undef, iseven(n) ? n+1 : n)  # since z -> 1/z we get one more coefficient
    ret[1]=f.coefficients[1]
    for k=2:2:length(ret)-1
        ret[k+1]=f.coefficients[k]
    end
    for k=2:2:n-1
        ret[k]=f.coefficients[k+1]
    end
    iseven(n) && (ret[n] = 0)

    Fun(Laurent(reverseorientation(domain(f))),ret)
end

include("calculus.jl")
include("specialfunctions.jl")
include("FourierOperators.jl")
include("LaurentOperators.jl")
include("LaurentDirichlet.jl")
include("roots.jl")

Fun(::typeof(identity), d::Circle) = Fun(Laurent(d),[d.center,0.,d.radius])

Space(d::PeriodicDomain) = Fourier(d)
Space(d::Circle) = Laurent(d)

## Evaluation

Evaluation(d::PeriodicDomain,x::Number,n...) = Evaluation(Laurent(d),complex(x),n...)

## Definite Integral

DefiniteIntegral(d::PeriodicDomain) = DefiniteIntegral(Laurent(d))
DefiniteLineIntegral(d::PeriodicDomain) = DefiniteLineIntegral(Laurent(d))

## Toeplitz
union_rule(A::Space{<:PeriodicSegment}, B::Space{<:IntervalOrSegment}) =
    union(Space(Interval(domain(A))), B)

    ## Derivative

function invfromcanonicalD(S::Laurent{PeriodicLine{false}})
    d=domain(S)
    @assert d.center==0  && d.L==1.0
    a=Fun(Laurent(),[1.,.5,.5])
end

function invfromcanonicalD(S::LaurentDirichlet{PeriodicLine{false}})
    d=domain(S)
    @assert d.center==0  && d.L==1.0
    a=Fun(Laurent(),[1.,.5,.5])
end

Space(d::PeriodicCurve{S}) where {S<:Fourier} = Fourier(d)
Space(d::PeriodicCurve{S}) where {S<:Laurent} = Laurent(d)

function Fun(T::ToeplitzOperator)
    if length(T.nonnegative)==1
       Fun(Taylor(),[T.nonnegative;T.negative])
     elseif length(T.negative)==0
         Fun(Hardy{false}(),T.nonnegative)
     else
         Fun(Laurent(Circle()),interlace(T.nonnegative,T.negative))
     end
 end

end #module


text:
# Variables in Gurobi model

# variable kinds

const GRB_CONTINUOUS = convert(Cchar, 'C')
const GRB_BINARY     = convert(Cchar, 'B')
const GRB_INTEGER    = convert(Cchar, 'I')

# add_var!
function add_var!(model::Model, numnz::Integer, vind::Vector, vval::Vector{Float64}, c::Float64, lb::Float64, ub::Float64, vtype::Cchar)
    if checkvalue(c, GRB_INFINITY)
        _objwarning(c)
    end
    if checkvalue(lb, GRB_BOUNDMAX) || checkvalue(ub, GRB_BOUNDMAX)
        _boundwarning(lb, ub)
    end
    ret = @grb_ccall(addvar, Cint, (
        Ptr{Cvoid},    # model
        Cint,         # numnz
        Ptr{Cint},    # vind
        Ptr{Float64}, # vval
        Float64,      # obj
        Float64,      # lb
        Float64,      # ub
        UInt8,        # vtype
        Ptr{UInt8}    # name
        ),
        model, numnz, ivec(vind.-1), vval, c, lb, ub, vtype, C_NULL)

    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function add_var!(model::Model, vtype::Cchar, c::Float64, lb::Float64, ub::Float64)
    if checkvalue(c, GRB_INFINITY)
        _objwarning(c)
    end
    if checkvalue(lb, GRB_BOUNDMAX) || checkvalue(ub, GRB_BOUNDMAX)
        _boundwarning(lb, ub)
    end
    ret = @grb_ccall(addvar, Cint, (
        Ptr{Cvoid},    # model
        Cint,         # numnz
        Ptr{Cint},    # vind
        Ptr{Float64}, # vval
        Float64,      # obj
        Float64,      # lb
        Float64,      # ub
        UInt8,        # vtype
        Ptr{UInt8}    # name
        ),
        model, 0, C_NULL, C_NULL, c, lb, ub, vtype, C_NULL)

    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

add_var!(model::Model, vtype::GChars, c::Real, lb::Real, ub::Real) = add_var!(model, cchar(vtype), Float64(c), Float64(lb), Float64(ub))
add_var!(model::Model, vtype::GChars, c::Real) = add_var!(model, vtype, c, -Inf, Inf)

add_cvar!(model::Model, c::Real, lb::Real, ub::Real) = add_var!(model, GRB_CONTINUOUS, c, lb, ub)
add_cvar!(model::Model, c::Real) = add_cvar!(model, c, -Inf, Inf)

add_bvar!(model::Model, c::Real) = add_var!(model, GRB_BINARY, c, 0., 1.)

add_ivar!(model::Model, c::Real, lb::Real, ub::Real) = add_var!(model, GRB_INTEGER, c, lb, ub)
add_ivar!(model::Model, c::Real) = add_ivar!(model, c, -Inf, Inf)

# add_vars!

function add_vars!(model::Model, vtypes::CVec, c::FVec, lb::FVec, ub::FVec)
    if checkvalue(c, GRB_INFINITY)
        _objwarning(c)
    end
    if checkvalue(lb, GRB_BOUNDMAX) || checkvalue(ub, GRB_BOUNDMAX)
        _boundwarning(lb, ub)
    end
    # check dimensions
    n = length(vtypes)
    _chklen(c, n)
    _chklen(lb, n)
    _chklen(ub, n)

    # main call
    ret = @grb_ccall(addvars, Cint, (
        Ptr{Cvoid},  # model
        Cint,       # numvars
        Cint,       # numnz
        Ptr{Cint},  # vbeg
        Ptr{Cint},  # vind
        Ptr{Float64}, # vval
        Ptr{Float64}, # obj
        Ptr{Float64}, # lb
        Ptr{Float64}, # ub
        Ptr{Cchar},   # vtypes
        Ptr{Ptr{UInt8}}, # varnames
        ),
        model, n, 0, C_NULL, C_NULL, C_NULL, c, lb, ub, vtypes, C_NULL)

    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function add_vars!(model::Model, vtypes::GCharOrVec, c::Vector, lb::Bounds, ub::Bounds)
    n = length(c)
    add_vars!(model, cvecx(vtypes, n), fvec(c), fvecx(lb, n), fvecx(ub, n))
end

add_cvars!(model::Model, c::Vector, lb::Bounds, ub::Bounds) = add_vars!(model, GRB_CONTINUOUS, c, lb, ub)
add_cvars!(model::Model, c::Vector) = add_cvars!(model, c, -Inf, Inf)

add_bvars!(model::Model, c::Vector) = add_vars!(model, GRB_BINARY, c, 0, 1)

add_ivars!(model::Model, c::Vector, lb::Bounds, ub::Bounds) = add_vars!(model, GRB_INTEGER, c, lb, ub)
add_ivars!(model::Model, c::Vector) = add_ivars!(model, GRB_INTEGER, c, -Inf, Inf)

del_vars!(model::Model, idx::T) where {T<:Real} = del_vars!(model, Cint[idx])
del_vars!(model::Model, idx::Vector{T}) where {T<:Real} = del_vars!(model, convert(Vector{Cint},idx))
function del_vars!(model::Model, idx::Vector{Cint})
    numdel = length(idx)
    ret = @grb_ccall(delvars, Cint, (
                     Ptr{Cvoid},
                     Cint,
                     Ptr{Cint}),
                     model, convert(Cint,numdel), idx.-Cint(1))
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
end

function get_vars(model::Model, start::Integer, len::Integer)
# http://www.gurobi.com/documentation/8.1/refman/c_grbgetvars.html

    # NR_CONTRAINTS
    #--------------
    m = num_constrs(model)
    # NR_VARIABLES
    #-------------
    n = num_vars(model)

    # INPUT VALIDATION
    #-----------------
    @assert start > 0               "Indexing in Julia starts from 1."
    @assert start <= n              string("Index out of bounds: There are only ", n, " variables attached to this model.")
    @assert len > 0                 "At least one variable must be selected; len > 0."
    @assert start * len <= m * n    "Maximal amount of possible non-zero elements surpassed."

    # FUNCTION CALLS
    #---------------
    numnzP = Ref{Cint}()
    vbeg = Array{Cint}(undef, len)

    ret = @grb_ccall(getvars, Cint, (
                        Ptr{Cvoid},     #Model
                        Ptr{Cint},      #numnzP
                        Ptr{Cint},      #vbeg
                        Ptr{Cint},      #vind
                        Ptr{Cdouble},   #vval
                        Cint,           #start
                        Cint),          #len
                        model, numnzP, vbeg, C_NULL, C_NULL, Cint(start - 1), Cint(len))
    if ret != 0
        throw(GurobiError(model.env, ret))
    end

    nnz = numnzP[]
    vind = Array{Cint}(undef, nnz)
    vval = Array{Cdouble}(undef, nnz)

    ret = @grb_ccall(getvars, Cint, (
                        Ptr{Cvoid},     #Model
                        Ptr{Cint},      #numnzP
                        Ptr{Cint},      #vbeg
                        Ptr{Cint},      #vind
                        Ptr{Cdouble},   #vval
                        Cint,           #start
                        Cint),          #len
                        model, numnzP, vbeg, vind, vval, Cint(start - 1), Cint(len))

    if ret != 0
        throw(GurobiError(model.env, ret))
    end

    # ADJUSTING INDICES TO JULIA'S INDEXING.
    #---------------------------------------
    for i in 1:size(vbeg, 1)
        vbeg[i] += 1
    end

    for i in 1:size(vind, 1)
        vind[i] += 1
    end

    # SPARSE ARRAY
    #-------------
    push!(vbeg, nnz)
    I = Array{Int64}(undef, nnz)
    J = Array{Int64}(undef, nnz)
    V = Array{Float64}(undef, nnz)
    for i in 1:length(vbeg) - 1
        for j in vbeg[i]:vbeg[i + 1]
            I[j] = vind[j]
            J[j] = i + start - 1
            V[j] = vval[j]
        end
    end

    pop!(vbeg)

    # return vbeg, vind, vval
    return SparseArrays.sparse(I, J, V, m, n)
end


text:
using Test
using Koala
using DataFrames
using Statistics

# loss functions:
y = [1, 2, 3, 4]
yhat = y .+ 1
@test isapprox(rms(y, yhat), 1.0)
@test isapprox(rmsl(y, yhat),
               sqrt((log(1/2)^2 + log(2/3)^2 + log(3/4)^2 + log(4/5)^2)/4))
@test isapprox(rmslp1(y, yhat),
               sqrt((log(2/3)^2 + log(3/4)^2 + log(4/5)^2 + log(5/6)^2)/4))
@test isapprox(rmsp(y, yhat), sqrt((1 + 1/4 + 1/9 + 1/16)/4))

# for when ROC is added as working dependency:
# y = ["n", "p", "n", "p", "n", "p"]
# yhat = [0.1, 0.2, 0.3, 0.6, 0.7, 0.8]
# @test auc("p")(y, yhat) ≈ 2/3

# helpers:
dict = Dict{String,Int}()
dict["c"] = 3
dict["b"] = 2
dict["a"] = 1
@test keys_ordered_by_values(dict) == ["a", "b", "c"]
bootstrap_resample_of_mean(randn(100))
v = [1,2,5,4,missing,7,missing]
@test hasmissing(v)
@test ismissingtype(eltype(v))
@test Koala.principaltype(eltype(v)) == Int
@test countmissing(v) == 2

v[5] = 42; v[7] = 42
@test eltype(purify(v)) == Int

X, y = load_boston();
train, test = partition(eachindex(y), 0.8); # 80:20 split
describe(X)
# get_meta(X)

transformer = Koala.FeatureSelector(features=[:Indus, :Chas])
transformerM = Machine(transformer, X)
@test transform(transformerM, X) == X[[:Indus, :Chas]]
[transformer, transformer]
[transformerM, transformerM]

rgs = ConstantRegressor()
[rgs, rgs]
mach = Machine(rgs, X, y, train, transformer_X=transformer)
show(stdout, mach, 3)
show(mach, 2)
[mach, mach]

fit!(mach, train)
score = err(mach, test)
println("score = $score")
fit!(mach) # fit again without recomputing cache
@test 9.5 < score && 9.6 > score

@test score == err(mach, test, raw=true)
@test predict(mach, X, test) ==
    predict(mach.model, mach.predictor, mach.Xt[test,:], false, false)

learning_curve(mach, train, test, [2, 4, 8, 1000], raw=false)
learning_curve(mach, train, test, [2000, 3000], restart=false)
cv(mach, vcat(test, train))

u,v = @curve r range(0, stop=10, length=50) (r^2 + 1)
u,v = @pcurve r range(0, stop=10, length=50) (r^2 + 1)
u,v,w = @curve r range(0, stop=10, length=5) s range(0,stop=5,length=4) r*s^2

# test split_seen_unseen
v = ['a', 'b', 'b', 'c',
     'a', 'd', 'a', 'b', 'e']
trainrows = 1:4
testrows = 5:9
@test split_seen_unseen(v, trainrows, testrows) == ([5, 7, 8], [6, 9])

w = ["log", "house", "house", "house",
     "brick", "house", "log", "log", "log"]
df = DataFrame(v=v, w=w)
@test split_seen_unseen(df, trainrows, testrows) == ([7, 8], [5, 6, 9])

# time tests:
v=rand(UInt32, 10^7)
trainrows, testrows = partition(eachindex(v), 0.9)
@time map(length, split_seen_unseen(v, trainrows, testrows))

# test drop_unseen capability
X[4] = map(Char, X[4])
X[9] = map(Char, X[9])
X[10] = map(Char, X[10])
train = 1:length(y) - 6
test = length(y) - 5 : length(y)
model = ConstantRegressor()
mach = Machine(model, X, y, train, drop_unseen=true)
fit!(mach, train)
err(mach, test)
learning_curve(mach, train, test, [2, 4, 8, 1000])
learning_curve(mach, train, test, [2000, 3000], restart=false)

t = Koala.RowsTransformer(10)
tM  = Machine(t, [2, 4, 5, 9])
@test transform(tM, 3:7) == [2, 3, 4]
df = DataFrame(letters=['a','b','c','d','e','f','g','h','i','j'])
df_shrunk = df[[1, 3, 6, 7, 8], :]
@test df_shrunk[transform(tM, 3:7), :letters] == ['c', 'f', 'g']

rgs = ConstantRegressor()
train, test = partition(eachindex(y), 0.8); # 80:20 split
mach = Machine(rgs, X, y, train, drop_unseen=true)
fit!(mach, train)
score = err(mach, test)
unseen = mach.rows_with_unseen
seen = filter(test) do i
    !(i in unseen)
end
@test mean(y[train]) ≈ mach.predictor
yhat = mean(y[train])*ones(length(y))
@test rms(yhat, y, seen) ≈ score
println("score = $score")
fit!(mach) # fit again without recomputing cache
@test 9 < score && 10 > score
@test score == err(mach, test, raw=true)
learning_curve(mach, train, test, [2, 4, 8, 1000], raw=false)
learning_curve(mach, train, test, [2000, 3000], restart=false)
cv(mach, vcat(test, train))

# test compete function:
@test compete(randn(1000) .+ 5, randn(1000))[1] == '1'

load_iris()
load_ames()


text:
using Dates
using Test
using HTTP
using HTTP: Headers, URI
using HTTP.AWS4AuthRequest: sign_aws4!

# Based on https://docs.aws.amazon.com/general/latest/gr/signature-v4-test-suite.html
# and https://docs.aws.amazon.com/AmazonS3/latest/API/sig-v4-header-based-auth.html

function test_sign!(method, headers, params, body=""; opts...)
    sign_aws4!(method,
               URI("https://example.amazonaws.com/" * params),
               headers,
               Vector{UInt8}(body);
               timestamp=DateTime(2015, 8, 30, 12, 36),
               aws_service="service",
               aws_region="us-east-1",
               # NOTE: These are the example credentials as specified in the AWS docs,
               # they are not real
               aws_access_key_id="AKIDEXAMPLE",
               aws_secret_access_key="wJalrXUtnFEMI/K7MDENG+bPxRfiCYEXAMPLEKEY",
               include_md5=false,
               include_sha256=false,
               opts...)
    headers
end

function test_auth_string(headers, sig, key="AKIDEXAMPLE", date="20150830", service="service")
    d = [
        "AWS4-HMAC-SHA256 Credential" => "$key/$date/us-east-1/$service/aws4_request",
        "SignedHeaders" => headers,
        "Signature" => sig,
    ]
    join(map(p->join(p, '='), d), ", ")
end

header_keys(headers) = sort!(map(first, headers))

const required_headers = ["Authorization", "host", "x-amz-date"]

@testset "AWS Signature Version 4" begin
    # The signature for requests with no headers where the path ends up as simply /
    slash_only_sig = "5fa00fa31553b73ebf1942676e86291e8372ff2a2260956d9b8aae1d763fbf31"
    noheaders = [
        ("get-vanilla", "", slash_only_sig),
        ("get-vanilla-empty-query-key", "?Param1=value1", "a67d582fa61cc504c4bae71f336f98b97f1ea3c7a6bfe1b6e45aec72011b9aeb"),
        ("get-utf8", "ሴ", "8318018e0b0f223aa2bbf98705b62bb787dc9c0e678f255a891fd03141be5d85"),
        ("get-relative", "example/..", slash_only_sig),
        ("get-relative-relative", "example1/example2/../..", slash_only_sig),
        ("get-slash", "/", slash_only_sig),
        ("get-slash-dot-slash", "./", slash_only_sig),
        ("get-slashes", "example/", "9a624bd73a37c9a373b5312afbebe7a714a789de108f0bdfe846570885f57e84"),
        ("get-slash-pointless-dot", "./example", "ef75d96142cf21edca26f06005da7988e4f8dc83a165a80865db7089db637ec5"),
        ("get-space", "example space/", "652487583200325589f1fba4c7e578f72c47cb61beeca81406b39ddec1366741"),
        ("post-vanilla", "", "5da7c1a2acd57cee7505fc6676e4e544621c30862966e37dddb68e92efbe5d6b"),
        ("post-vanilla-empty-query-value", "?Param1=value1", "28038455d6de14eafc1f9222cf5aa6f1a96197d7deb8263271d420d138af7f11"),
    ]
    @testset "$name" for (name, p, sig) in noheaders
        m = startswith(name, "get") ? "GET" : "POST"
        headers = test_sign!(m, Headers([]), p)
        @test header_keys(headers) == required_headers
        d = Dict(headers)
        @test d["x-amz-date"] == "20150830T123600Z"
        @test d["host"] == "example.amazonaws.com"
        @test d["Authorization"] == test_auth_string("host;x-amz-date", sig)
    end

    yesheaders = [
        ("get-header-key-duplicate", "", "",
         Headers(["My-Header1" => "value2",
                  "My-Header1" => "value2",
                  "My-Header1" => "value1"]),
         "host;my-header1;x-amz-date",
         "c9d5ea9f3f72853aea855b47ea873832890dbdd183b4468f858259531a5138ea"),
        ("get-header-value-multiline", "", "",
         Headers(["My-Header1" => "value1\n  value2\n    value3"]),
         "host;my-header1;x-amz-date",
         "ba17b383a53190154eb5fa66a1b836cc297cc0a3d70a5d00705980573d8ff790"),
        ("get-header-value-order", "", "",
         Headers(["My-Header1" => "value4",
                  "My-Header1" => "value1",
                  "My-Header1" => "value3",
                  "My-Header1" => "value2"]),
         "host;my-header1;x-amz-date",
         "08c7e5a9acfcfeb3ab6b2185e75ce8b1deb5e634ec47601a50643f830c755c01"),
        ("get-header-value-trim", "", "",
         Headers(["My-Header1" => " value1",
                  "My-Header2" => " \"a   b   c\""]),
         "host;my-header1;my-header2;x-amz-date",
         "acc3ed3afb60bb290fc8d2dd0098b9911fcaa05412b367055dee359757a9c736"),
        ("post-header-key-sort", "", "",
         Headers(["My-Header1" => "value1"]),
         "host;my-header1;x-amz-date",
         "c5410059b04c1ee005303aed430f6e6645f61f4dc9e1461ec8f8916fdf18852c"),
        ("post-header-value-case", "", "",
         Headers(["My-Header1" => "VALUE1"]),
         "host;my-header1;x-amz-date",
         "cdbc9802e29d2942e5e10b5bccfdd67c5f22c7c4e8ae67b53629efa58b974b7d"),
        ("post-x-www-form-urlencoded", "", "Param1=value1",
         Headers(["Content-Type" => "application/x-www-form-urlencoded",
                  "Content-Length" => "13"]),
         "content-type;host;x-amz-date",
         "ff11897932ad3f4e8b18135d722051e5ac45fc38421b1da7b9d196a0fe09473a"),
        ("post-x-www-form-urlencoded-parameters", "", "Param1=value1",
         Headers(["Content-Type" => "application/x-www-form-urlencoded; charset=utf8",
                  "Content-Length" => "13"]),
         "content-type;host;x-amz-date",
         "1a72ec8f64bd914b0e42e42607c7fbce7fb2c7465f63e3092b3b0d39fa77a6fe"),
    ]
    @testset "$name" for (name, p, body, h, sh, sig) in yesheaders
        hh = sort(map(first, h))
        m = startswith(name, "get") ? "GET" : "POST"
        test_sign!(m, h, p, body)
        @test header_keys(h) == sort(vcat(required_headers, hh))
        d = Dict(h) # collapses duplicates but we don't care here
        @test d["x-amz-date"] == "20150830T123600Z"
        @test d["host"] == "example.amazonaws.com"
        @test d["Authorization"] == test_auth_string(sh, sig)
    end

    @testset "AWS Security Token Service" begin
        # Not a real security token, provided by AWS as an example
        token = string("AQoDYXdzEPT//////////wEXAMPLEtc764bNrC9SAPBSM22wDOk4x4HIZ8j4FZTwd",
                       "QWLWsKWHGBuFqwAeMicRXmxfpSPfIeoIYRqTflfKD8YUuwthAx7mSEI/qkPpKPi/k",
                       "McGdQrmGdeehM4IC1NtBmUpp2wUE8phUZampKsburEDy0KPkyQDYwT7WZ0wq5VSXD",
                       "vp75YU9HFvlRd8Tx6q6fE8YQcHNVXAkiY9q6d+xo0rKwT38xVqr7ZD0u0iPPkUL64",
                       "lIZbqBAz+scqKmlzm8FDrypNC9Yjc8fPOLn9FX9KSYvKTr4rvx3iSIlTJabIQwj2I",
                       "CCR/oLxBA==")
        @testset "Token included in signature" begin
            sh = "host;x-amz-date;x-amz-security-token"
            sig = "85d96828115b5dc0cfc3bd16ad9e210dd772bbebba041836c64533a82be05ead"
            h = test_sign!("POST", Headers([]), "", aws_session_token=token)
            d = Dict(h)
            @test d["Authorization"] == test_auth_string(sh, sig)
            @test haskey(d, "x-amz-security-token")
        end
        @testset "Token not included in signature" begin
            sh = "host;x-amz-date"
            sig = "5da7c1a2acd57cee7505fc6676e4e544621c30862966e37dddb68e92efbe5d6b"
            h = test_sign!("POST", Headers([]), "", aws_session_token=token, token_in_signature=false)
            d = Dict(h)
            @test d["Authorization"] == test_auth_string(sh, sig)
            @test haskey(d, "x-amz-security-token")
        end
    end

    @testset "AWS Simple Storage Service" begin
        s3url = "https://examplebucket.s3.amazonaws.com"
        opts = (timestamp=DateTime(2013, 5, 24),
                aws_service="s3",
                aws_region="us-east-1",
                # NOTE: These are the example credentials as specified in the AWS docs,
                # they are not real
                aws_access_key_id="AKIAIOSFODNN7EXAMPLE",
                aws_secret_access_key="wJalrXUtnFEMI/K7MDENG/bPxRfiCYEXAMPLEKEY",
                include_md5=false)

        @testset "GET Object" begin
            sh = "host;range;x-amz-content-sha256;x-amz-date"
            sig = "f0e8bdb87c964420e857bd35b5d6ed310bd44f0170aba48dd91039c6036bdb41"
            h = Headers(["Range" => "bytes=0-9"])
            sign_aws4!("GET", URI(s3url * "/test.txt"), h, UInt8[]; opts...)
            d = Dict(h)
            @test d["Authorization"] == test_auth_string(sh, sig, opts.aws_access_key_id, "20130524", "s3")
            @test haskey(d, "x-amz-content-sha256") # required for S3 requests
        end

        @testset "PUT Object" begin
            sh = "date;host;x-amz-content-sha256;x-amz-date;x-amz-storage-class"
            sig = "98ad721746da40c64f1a55b78f14c238d841ea1380cd77a1b5971af0ece108bd"
            h = Headers(["Date" => "Fri, 24 May 2013 00:00:00 GMT",
                         "x-amz-storage-class" => "REDUCED_REDUNDANCY"])
            sign_aws4!("PUT", URI(s3url * "/test\$file.text"), h, UInt8[];
                       # Override the SHA-256 of the request body, since the actual body is not provided
                       # for this example in the documentation, only the SHA
                       body_sha256=hex2bytes("44ce7dd67c959e0d3524ffac1771dfbba87d2b6b4b4e99e42034a8b803f8b072"),
                       opts...)
            d = Dict(h)
            @test d["Authorization"] == test_auth_string(sh, sig, opts.aws_access_key_id, "20130524", "s3")
            @test haskey(d, "x-amz-content-sha256")
        end

        @testset "GET Bucket Lifecycle" begin
            sh = "host;x-amz-content-sha256;x-amz-date"
            sig = "fea454ca298b7da1c68078a5d1bdbfbbe0d65c699e0f91ac7a200a0136783543"
            h = Headers([])
            sign_aws4!("GET", URI(s3url * "/?lifecycle"), h, UInt8[]; opts...)
            d = Dict(h)
            @test d["Authorization"] == test_auth_string(sh, sig, opts.aws_access_key_id, "20130524", "s3")
            @test haskey(d, "x-amz-content-sha256")
        end

        @testset "GET Bucket (List Objects)" begin
            sh = "host;x-amz-content-sha256;x-amz-date"
            sig = "34b48302e7b5fa45bde8084f4b7868a86f0a534bc59db6670ed5711ef69dc6f7"
            h = Headers([])
            sign_aws4!("GET", URI(s3url * "/?max-keys=2&prefix=J"), h, UInt8[]; opts...)
            d = Dict(h)
            @test d["Authorization"] == test_auth_string(sh, sig, opts.aws_access_key_id, "20130524", "s3")
            @test haskey(d, "x-amz-content-sha256")
        end
    end

    @testset "HTTP.request with AWS authentication" begin
        resp = HTTP.request("GET",
                            "https://httpbin.org/headers";
                            aws_authorization=true,
                            timestamp=DateTime(2015, 8, 30, 12, 36),
                            aws_service="service",
                            aws_region="us-east-1",
                            # NOTE: These are the example credentials as specified in the AWS docs,
                            # they are not real
                            aws_access_key_id="AKIDEXAMPLE",
                            aws_secret_access_key="wJalrXUtnFEMI/K7MDENG+bPxRfiCYEXAMPLEKEY",
                            include_md5=false,
                            include_sha256=false)
        @test resp.status == 200
        headers = JSON.parse(String(resp.body))["headers"]
        @test headers["Host"] == "httpbin.org"
        @test headers["X-Amz-Date"] == "20150830T123600Z"
        auth = "AWS4-HMAC-SHA256 " *
               "Credential=AKIDEXAMPLE/20150830/us-east-1/service/aws4_request, " *
               "SignedHeaders=host;x-amz-date, " *
               "Signature=5d7c2f752f5597c1a2ad2f65008974307a7b25bd7cc4f79146b06435e91fa95e"
        @test headers["Authorization"] == auth
    end
end