text:
import MathOptInterface
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities

using .ClpCInterface

# Helper functions
_bounds(s::MOI.GreaterThan{Float64}) = (s.lower, Inf)
_bounds(s::MOI.LessThan{Float64}) = (-Inf, s.upper)
_bounds(s::MOI.EqualTo{Float64}) = (s.value, s.value)
_bounds(s::MOI.Interval{Float64}) = (s.lower, s.upper)

# Supported scalar sets
const SCALAR_SETS = Union{
    MOI.GreaterThan{Float64},
    MOI.LessThan{Float64},
    MOI.EqualTo{Float64},
    MOI.Interval{Float64}
}

# Maps Clp's parameters to getter/setter function
const CLP_OPTION_MAP = Dict(
    :PrimalTolerance => (Clp.primal_tolerance, Clp.set_primal_tolerance),
    :DualTolerance => (Clp.dual_tolerance, Clp.set_dual_tolerance),
    :DualObjectiveLimit => (Clp.dual_objective_limit, Clp.set_dual_objective_limit),
    :MaximumIterations => (Clp.maximum_iterations, Clp.set_maximum_iterations),
    :MaximumSeconds => (Clp.maximum_seconds, Clp.set_maximum_seconds),
    :LogLevel => (Clp.log_level, Clp.set_log_level),
    :Scaling => (Clp.scaling_flag, Clp.scaling),
    :Perturbation => (Clp.perturbation, Clp.set_perturbation),
    :Algorithm => (Clp.algorithm, Clp.set_algorithm)
)

const SOLVE_OPTION_MAP = Dict(
   :PresolveType => (Clp.get_presolve_type, Clp.set_presolve_type),
   :SolveType => (Clp.get_solve_type, Clp.set_solve_type),
   :InfeasibleReturn => (Clp.infeasible_return, Clp.set_infeasible_return)
)

mutable struct Optimizer <: MOI.AbstractOptimizer
    # Inner Clp object
    inner::Clp.ClpModel

    # Clp solve options
    solver_options::Clp.ClpSolve
    options_set::Set{Symbol}
    # To handle MOI.OPTIMIZE_NOT_CALLED status
    optimize_called::Bool

    function Optimizer(;kwargs...)
        inner_model = Clp.ClpModel()
        solver_options = Clp.ClpSolve()
        model = new(inner_model, solver_options, Set{Symbol}(), false)

        for (key, value) in kwargs
            MOI.set(model, MOI.RawParameter(String(key)), value)
        end

        return model
    end
end

# ====================
#   empty functions
# ====================
function MOI.is_empty(model::Optimizer)
    # A problem is empty if it has no variable and no linear constraint
    return (Clp.get_num_rows(model.inner) == 0) && (Clp.get_num_cols(model.inner) == 0)
end

function MOI.empty!(model::Optimizer)
    old_model = model.inner
    model.inner = Clp.ClpModel()
    # Copy parameters from old model into new model
    for key in model.options_set
        getter, setter = CLP_OPTION_MAP[key]
        setter(model.inner, getter(old_model))
    end
    model.optimize_called = false
    return
end

MOI.get(::Optimizer, ::MOI.SolverName) = "Clp"

# TODO: improve type-stability of `MOI.RawParameter`-related methods.
MOI.supports(::Optimizer, param::MOI.RawParameter) = true

function MOI.set(model::Optimizer, param::MOI.RawParameter, value)
    key = Symbol(param.name)
    if haskey(CLP_OPTION_MAP, key)
        push!(model.options_set, key)
        CLP_OPTION_MAP[key][2](model.inner, value)
    elseif haskey(SOLVE_OPTION_MAP, key)
        SOLVE_OPTION_MAP[key][2](model.solver_options, value)
    else
        throw(MOI.UnsupportedAttribute(param))
    end
    return
end

function MOI.get(model::Optimizer, param::MOI.RawParameter)
    key = Symbol(param.name)
    if haskey(CLP_OPTION_MAP, key)
        return CLP_OPTION_MAP[key][1](model.inner)
    elseif haskey(SOLVE_OPTION_MAP, key)
        return SOLVE_OPTION_MAP[key][1](model.solver_options)
    else
        throw(MOI.UnsupportedAttribute(param))
    end
end

MOI.supports(::Optimizer, ::MOI.Silent) = true
function MOI.set(model::Optimizer, ::MOI.Silent, value::Bool)
    push!(model.options_set, :LogLevel)
    if value
        Clp.set_log_level(model.inner, 0)
    else
        Clp.set_log_level(model.inner, 1)
    end
    return nothing
end

function MOI.get(model::Optimizer, ::MOI.Silent)
    return Clp.log_level(model.inner) == 0
end

MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true

function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, value)
    push!(model.options_set, :MaximumSeconds)
    if value === nothing
        # De-activate time limit
        Clp.set_maximum_seconds(model.inner, -1.0)
    else
        Clp.set_maximum_seconds(model.inner, value)
    end
    return nothing
end

# Clp behaves weird here
MOI.get(model::Optimizer, ::MOI.TimeLimitSec) = Clp.maximum_seconds(model.inner)

MOI.supports(::Optimizer, ::MOI.NumberOfThreads) = false

# ========================================
#   Supported constraints and objectives
# ========================================
MOI.supports_constraint(::Optimizer,
    ::Type{MOI.ScalarAffineFunction{Float64}}, ::Type{<:SCALAR_SETS}
) = true

MOI.supports_constraint(::Optimizer,
    ::Type{MOI.SingleVariable}, ::Type{<:SCALAR_SETS}
) = true

MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true

MOI.supports(::Optimizer,
    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}
) = true

# =======================
#   `copy_to` function
# =======================

function MOI.copy_to(
    dest::Optimizer,
    src::MOI.ModelLike;
    copy_names=false
)
    # Check that all constraints and objectives
    for (F, S) in MOI.get(src, MOI.ListOfConstraints())
        MOI.supports_constraint(dest, F, S) || throw(MOI.UnsupportedConstraint{F, S}(
            "Clp.Optimizer does not support constraints of type $F-in-$S."
        ))
    end
    fobj_type = MOI.get(src, MOI.ObjectiveFunctionType())
    MOI.supports(dest, MOI.ObjectiveFunction{fobj_type}()) || throw(
        MOI.UnsupportedAttribute(MOI.ObjectiveFunction(fobj))
    )

    # Empty dest
    MOI.empty!(dest)

    mapping = MOIU.IndexMap()

    # First create variables (including bounds)

    # Create variables
    # TODO: it may be more efficient to create the problem all at once
    x_src = MOI.get(src, MOI.ListOfVariableIndices())
    for (j, x) in enumerate(x_src)
        # Variable j corresponds to j-th variable of x_src
        mapping.varmap[x] = MOI.VariableIndex(j)

        # Possible bound combinations:
        #=
            * No bound (default case)
            * Interval
            * EqualTo
            * LessThan
            * GreaterThan
            * LessThan && GreaterThan
        =#
        lb, ub = -Inf, Inf  # Default case: free variable

        idx_RG = MOI.ConstraintIndex{MOI.SingleVariable, MOI.Interval{Float64}}(x.value)
        idx_ET = MOI.ConstraintIndex{MOI.SingleVariable, MOI.EqualTo{Float64}}(x.value)
        idx_LT = MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}}(x.value)
        idx_GT = MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Float64}}(x.value)
        if MOI.is_valid(src, idx_LT)
            # Get the (upper) bound
            s = MOI.get(src, MOI.ConstraintSet(), idx_LT)
            ub = s.upper

            # Track index of upper-bound constraint
            mapping.conmap[idx_LT] = MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}}(j)
        end
        if MOI.is_valid(src, idx_GT)
            # Get the (lower) bound
            s = MOI.get(src, MOI.ConstraintSet(), idx_GT)
            lb = s.lower

            # Track index of upper-bound constraint
            mapping.conmap[idx_GT] = MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Float64}}(j)
        end
        if MOI.is_valid(src, idx_RG)
            # Get upper and lower bound
            s = MOI.get(src, MOI.ConstraintSet(), idx_RG)
            lb, ub = s.lower, s.upper

            # Track index of upper-bound constraint
            mapping.conmap[idx_RG] = MOI.ConstraintIndex{MOI.SingleVariable, MOI.Interval{Float64}}(j)
        end
        if MOI.is_valid(src, idx_ET)
            # Get the bounds
            s = MOI.get(src, MOI.ConstraintSet(), idx_ET)
            lb, ub = s.value, s.value

            # Track index of upper-bound constraint
            mapping.conmap[idx_ET] = MOI.ConstraintIndex{MOI.SingleVariable, MOI.EqualTo{Float64}}(j)
        end

        Clp.add_column(dest.inner, Cint(0), Int32[], Float64[], lb, ub, 0.0)
    end

    # Create linear constraints
    nrows = 0  # number of rows in the problem
    for (F, S) in MOI.get(src, MOI.ListOfConstraints())
        # Check if constraint is a linear constraint
        (F == MOI.ScalarAffineFunction{Float64}) || continue

        for index in MOI.get(src, MOI.ListOfConstraintIndices{F, S}())
            # Retrieve problem data
            # TODO: ensure that `(f, s)` are in canonical form
            f = MOIU.canonical(MOI.get(src, MOI.ConstraintFunction(), index))
            s = MOI.get(src, MOI.ConstraintSet(), index)

            # Extract bounds
            # Constraint writes a'x + a0 - in - set
            # So we offset the bounds by -a0
            lb, ub = _bounds(s) .- f.constant

            # Identify constraint terms
            coeffs = [t.coefficient for t in f.terms]
            cols::Vector{Int} = [
                mapping.varmap[t.variable_index].value
                for t in f.terms
            ]

            # Add constraint
            Clp.add_row(dest.inner, Cint(length(f.terms)), Cint.(cols .- 1), coeffs, lb, ub)

            # Track index
            mapping.conmap[index] = MOI.ConstraintIndex{F, S}(nrows + 1)
            nrows += 1
        end
    end

    # Objective function
    sense = MOI.get(src, MOI.ObjectiveSense())
    if sense == MOI.MIN_SENSE
        Clp.set_obj_sense(dest.inner, 1)   # minimize
    elseif sense == MOI.MAX_SENSE
        Clp.set_obj_sense(dest.inner, -1)  # maximize
    elseif sense == MOI.FEASIBILITY_SENSE
        Clp.set_obj_sense(dest.inner, 0)   # feasibility
    else
        # Should not be reached, but just in case
        error("Unexpected optimization sense $sense")
    end

    fobj = MOI.get(src, MOI.ObjectiveFunction{fobj_type}())
    obj_coeffs = zeros(length(x_src))
    obj_offset = 0.0

    if fobj_type == MOI.ScalarAffineFunction{Float64}
        # Record objective coeffs
        for term in fobj.terms
            x = mapping.varmap[term.variable_index]
            # Here we do += instead of = to accumulate dupplicates
            obj_coeffs[x.value] += term.coefficient
        end
        obj_offset = fobj.constant
    else
        # Should not be reached here, but just in case
        error("Unexpected objective $fobj")
    end

    # Set objective coefficients
    Clp.chg_obj_coefficients(dest.inner, obj_coeffs)
    # Set objective offset
    # Clp seems to negates the objective offset
    Clp.set_objective_offset(dest.inner, -obj_offset)

    return mapping
end

# ===============================
#   Optimize and post-optimize
# ===============================

function MOI.optimize!(model::Optimizer)
    Clp.initial_solve_with_options(model.inner, model.solver_options)
    model.optimize_called = true
    return nothing
end

MOI.get(model::Optimizer, ::MOI.NumberOfVariables) = Clp.get_num_cols(model.inner)

function MOI.get(model::Optimizer, attr::MOI.ObjectiveValue)
    # Check result index
    MOI.check_result_index_bounds(model, attr)

    return Clp.get_obj_value(model.inner)
end

function MOI.get(model::Optimizer, ::MOI.TerminationStatus)

    model.optimize_called || return MOI.OPTIMIZE_NOT_CALLED

    st = Clp.status(model.inner)

    if st == 0
        return MOI.OPTIMAL
    elseif st == 1
        return MOI.INFEASIBLE
    elseif st == 2
        return MOI.DUAL_INFEASIBLE
    elseif st == 3
        # No more granular information that "some limit is reached"
        return MOI.OTHER_LIMIT
    else
        return MOI.OTHER_ERROR
    end
end

function MOI.get(model::Optimizer, ::MOI.RawStatusString)
    model.optimize_called || return "MOI.OPTIMIZE_NOT_CALLED"

    st = Clp.status(model.inner)
    if st == 0
        return "0 - optimal"
    elseif st == 1
        return "1 - primal infeasible"
    elseif st == 2
        return "2 - dual infeasible"
    elseif st == 3
        return "3 - stopped on iterations etc"
    elseif st == 4
        return "4 - stopped due to errors"
    else
        error("Expected integer in [0, 4] but got $st")
    end
end

function MOI.get(model::Optimizer, ::MOI.ResultCount)
    st = MOI.get(model, MOI.TerminationStatus())

    # TODO: check that this is correct
    if st == MOI.OPTIMIZE_NOT_CALLED || st == MOI.OTHER_ERROR
        return 0
    end
    return 1
end

# Solution status
function MOI.get(model::Optimizer, attr::MOI.PrimalStatus)
    # Check result index
    attr.N == 1 || return MOI.NO_SOLUTION

    if Clp.is_proven_dual_infeasible(model.inner)
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif Clp.primal_feasible(model.inner)
        return MOI.FEASIBLE_POINT
    else
        return MOI.UNKNOWN_RESULT_STATUS
    end
end

function MOI.get(model::Optimizer, attr::MOI.DualStatus)
    # Check result index
    attr.N == 1 || return MOI.NO_SOLUTION

    if Clp.is_proven_primal_infeasible(model.inner)
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif Clp.dual_feasible(model.inner)
        return MOI.FEASIBLE_POINT
    else
        return MOI.UNKNOWN_RESULT_STATUS
    end
end

# ===================
#   Primal solution
# ===================
function MOI.get(model::Optimizer, attr::MOI.VariablePrimal, x::MOI.VariableIndex)
    # Check result index
    MOI.check_result_index_bounds(model, attr)

    # We assume the variable index is valid,
    # since Clp should be accessed via a CachingOptimizer
    primal_status = MOI.get(model, MOI.PrimalStatus())
    if primal_status == MOI.INFEASIBILITY_CERTIFICATE
        sol = Clp.unbounded_ray(model.inner)
        return sol[x.value]
    elseif primal_status == MOI.FEASIBLE_POINT
        sol = Clp.get_col_solution(model.inner)
        return sol[x.value]
    else
        error("Primal solution not available")
    end
end

function MOI.get(model::Optimizer, attr::MOI.VariablePrimal, xs::Vector{MOI.VariableIndex})
    # Check result index
    MOI.check_result_index_bounds(model, attr)

    col_indices = [idx.value for idx in xs]

    # We assume all variable indices is valid,
    # since Clp should be accessed via a CachingOptimizer
    primal_status = MOI.get(model, MOI.PrimalStatus())
    if primal_status == MOI.INFEASIBILITY_CERTIFICATE
        sol = Clp.unbounded_ray(model.inner)
        return sol[col_indices]
    elseif primal_status == MOI.FEASIBLE_POINT
        sol = Clp.get_col_solution(model.inner)
        return sol[col_indices]
    else
        error("Primal solution not available")
    end
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, S}
) where{S <:SCALAR_SETS}
    # Check result index
    MOI.check_result_index_bounds(model, attr)

    # We assume that the constraint index is valid,
    # since Clp should be accessed via a CachingOptimizer
    # TODO: What happens if model is unbounded / infeasible?

    # Get primal row solution
    sol = Clp.get_row_activity(model.inner)
    return sol[c.value]
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal,
    c::MOI.ConstraintIndex{MOI.SingleVariable, S}
) where{S <: SCALAR_SETS}
    # Check result index
    MOI.check_result_index_bounds(model, attr)

    # TODO: What happens if model is unbounded / infeasible?
    return MOI.get(model, MOI.VariablePrimal(), MOI.VariableIndex(c.value))
end

# =================
#   Dual solution
# =================
# If sense is maximize, we negate all the duals to follow MOI conventions
# Feasibility problems are treated as a minimization
function MOI.get(model::Optimizer, attr::MOI.ConstraintDual,
    c::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, S}
) where{S <: SCALAR_SETS}
    # Check result index
    MOI.check_result_index_bounds(model, attr)

    sense = (Clp.get_obj_sense(model.inner) == -1) ? -1 : 1

    dual_status = MOI.get(model, MOI.DualStatus())
    if dual_status == MOI.FEASIBLE_POINT
        # Return dual solution
        dsol = Clp.get_row_price(model.inner)
        return sense * dsol[c.value]
    elseif dual_status == MOI.INFEASIBILITY_CERTIFICATE
        # Return infeasibility ray
        dsol = Clp.infeasibility_ray(model.inner)
        return sense * dsol[c.value]
    else
        # Throw error if dual solution is not available
        error("Dual solution not available")
    end
end

# Reduced costs
# TODO: what happens if problem is unbounded / infeasible?
function MOI.get(model::Optimizer, attr::MOI.ConstraintDual,
    c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.LessThan{Float64}}
)
    # Check result index
    MOI.check_result_index_bounds(model, attr)

    rc = Clp.get_reduced_cost(model.inner)[c.value]
    sense = (Clp.get_obj_sense(model.inner) == -1) ? -1 : 1

    # Dual should be non-positive
    if sense == 1 && rc <= 0.0
        return rc
    elseif sense == -1 && rc >= 0.0
        return -rc
    else
        return 0.0
    end
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintDual,
    c::MOI.ConstraintIndex{MOI.SingleVariable, MOI.GreaterThan{Float64}}
)
    # Check result index
    MOI.check_result_index_bounds(model, attr)

    rc = Clp.get_reduced_cost(model.inner)[c.value]
    sense = (Clp.get_obj_sense(model.inner) == -1) ? -1 : 1

    # Dual should be non-negative
    if sense == 1 && rc >= 0.0
        return rc
    elseif sense == -1 && rc <= 0.0
        return -rc
    else
        return 0.0
    end
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintDual,
    c::MOI.ConstraintIndex{MOI.SingleVariable, S}
) where{S <: Union{MOI.Interval{Float64}, MOI.EqualTo{Float64}}}
    # Check result index
    MOI.check_result_index_bounds(model, attr)

    sense = (Clp.get_obj_sense(model.inner) == -1) ? -1 : 1
    return sense * Clp.get_reduced_cost(model.inner)[c.value]
end


text:
module DummyStorageBackends

using DataSets

struct DummyBackend
    data
end

function Base.open(f::Function, ::Type{IO}, storage::DummyBackend, path; kws...) where {T}
    @assert isempty(path)
    f(IOBuffer(storage.data))
end

function connect_dummy_backend(f, config, ds)
    storage = DummyBackend(config["data"])
    f(Blob(storage))
end

function __init__()
    DataSets.add_storage_driver("DummyTomlStorage"=>connect_dummy_backend)
end

end

ext: jl


text:
@doc doc"""
bbands{Float64}(x::Vector{Float64}, n::Int64=10, sigma::Float64=2.0)

Bollinger Bands (moving average with standard deviation bands)
""" ->
function bbands{Float64}(x::Vector{Float64}, n::Int64=10, sigma::Float64=2.0; ma::Function=sma, args...)
    @assert n<size(x,1) && n>0 "Argument n is out of bounds."
    out = zeros(size(x,1), 3)  # cols := lower bound, ma, upper bound
    out[:,2] = ma(x, n; args...)
    sd = runsd(x, n, false)
    out[:,1] = out[:,2] - sigma*sd
    out[:,3] = out[:,2] + sigma*sd
    return out
end

@doc doc"""
tr{Float64}(hlc::Array{Float64,2})

True Range
""" ->
function tr{Float64}(hlc::Array{Float64,2})
    if size(hlc,2) != 3
        error("HLC array must have 3 columns.")
    end
    n = size(hlc,1)
    out = zeros(n)
    out[1] = NaN
    @inbounds for i=2:n
        out[i] = max(hlc[i,1]-hlc[i,2], hlc[i,1]-hlc[i-1,3], hlc[i-1,3]-hlc[i,2])
    end
    return out[:,1]
end

@doc doc"""
atr{Float64}(hlc::Array{Float64,2}, n::Int64=14)

Average true range (uses exponential moving average)
""" ->
function atr{Float64}(hlc::Array{Float64,2}, n::Int64=14; ma::Function=ema)
    @assert n<size(hlc,1) && n>0 "Argument n out of bounds."
    return [NaN; ma(tr(hlc)[2:end], n)]
end

@doc doc"""
keltner{Float64}(hlc::Array{Float64,2}, nema::Int64=20, natr::Int64=10, mult::Int64=2)

Keltner Bands
""" ->
function keltner{Float64}(hlc::Array{Float64,2}, nema::Int64=20, natr::Int64=10, mult::Int64=2)
    out = zeros(size(hlc,1), 3)
    out[:,2] = ema(hlc[:,3], nema)
    out[:,1] = out[:,2] - mult*atr(hlc, natr)
    out[:,3] = out[:,2] + mult*atr(hlc, natr)
    return out
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Zellij_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Zellij")
JLLWrappers.@generate_main_file("Zellij", UUID("f51964ae-9c2b-58f9-be08-f952495c30f4"))
end  # module Zellij_jll


text:
# # Tutorial 1

# Crash course in Julia basics:

# Arithmetic, arrays, tuples, strings, dictionaries, functions,
# iteration, random numbers, package loading, plotting

# (40 min)

# ## Julia is a calculator:

1 + 2^3

#-

sqrt(1 + 2^3) # do `sqrt(ans)` in REPL

#-

sin(pi)

# Query a function's document string using `?sin` at the REPL, or in a
# notebook:

@doc sin

# I've forgotten how the arcsin is called. Is it `asin` or `arcsin`? I
# can search all doc-strings containing "sine" to locate the method:

apropos("sine")

# Okay, I see it's called `asin`.

asin(1 + 3*im)

# ## Arrays

# One dimensional vectors:

v = [3, 5, 7]

#-

length(v)

# A "row vector" is a 1 x n array:

row = [3 5 7]

#-

# Multiple row vectors separated by semicolons or new lines define matrices:

A = [3 5 7
     2 4 6
     1 3 5]

#-

size(A)

#-

length(A)

# Accessing elements (Julia indices start at 1 not 0):

A[1, 2]

#-

A[1, 2] == A[2]

#-

A[:, 2] # 2nd column

# Changing elements:

A[1, 1] = 42

#-

inv(A) # inverse

#-

isapprox(inv(A)*v, A\v) # but RHS more efficient

# ## "Variables" in Julia *point* to objects

# Corollary: all passing of function arguments is pass by reference

# Like Python; Unlike R, C or FORTRAN.

#-

w = v

w

#-

v[1] = 42

#-

v

#-

w

#

# ## Tuples

# Similar to vectors but of fixed length and immutable (cannot be changed)

t1 = (1, 2.0, "cat")
typeof(t1)

#-

t1[3]

# Tuples also come in a *named* variety:

t2 = (i = 1, x = 2.0, animal="cat")

#-

t2.x

# ## Strings and relatives

a_string = "the cat"
a_character = 't'
a_symbol = :t

#-

a_string[1] == a_character

# A `Symbol` is string-like but
# [interned](https://en.wikipedia.org/wiki/String_interning). Generally
# use `String` for ordinary textual data, but use `Symbol` for
# language reflection (metaprogramming). For example:

isdefined(Main, :z)

#-

z = 1 + 2im
isdefined(Main, :z)

#-

z.im

#-

fieldnames(typeof(z))

# Symbols are generalized by *expressions*:

ex = :(z == 3)
eval(ex)

# If this is confusing, forget it for now.

# ## Dictionaries

d = Dict('a' => "ant", 'z' => "zebra")

#-

d['a']

#-

d['b'] = "bat"
d

#-

keys(d)

# The expression 'a' => "ant" is itself a stand-alone object:

pair = 'a' => "ant"

first(pair)

# ## Functions

# Three ways to define a generic function:

foo(x) = x^2 # METHOD 1 (inline)
foo(3)

# or

3 |> foo

# or

3 |> x -> x^2 # METHOD 2 (anonymous)

# or

function foo2(x) # METHOD 3 (verbose)
    y = x
    z = y
    w = z
    return w^2
end

foo2(3)

# ## Basic iteration

# Here are three ways to square the integers from 1 to 10.

# METHOD 1 (explicit loop):

squares = [] # or Int[] if performance matters
for x in 1:10
    push!(squares, x^2)
end

squares

# METHOD 2 (comprehension):

[x^2 for x in 1:10]

# METHOD 3 (map):

map(x -> x^2, 1:10)

# METHOD 4 (broadcasting with dot syntax):

(1:10) .^ 2

# ## Random numbers

typeof(2)

#-

rand() # sample a Float64 uniformly from interval [0, 1]

#-

rand(3, 4) # do that 12 times and put in a 3 x 4 array

#-

randn(3, 4) # use normal distribution instead

#-

rand(Int8) # random elment of type Int8

#-

rand(['a', 'b', 'c'], 10) # 10 random elements from a vector

# Some standard libraries are needed to do more, for example:

using Random

#-

randstring(30)

#-

using Statistics  # part of standard library

#-

y = rand(30)
@show mean(y) quantile(y, 0.75);

# (Use the macro @show before stuff you want printed prefixed by
# *what* it is that is being printed.)

# For sampling from more general distributions we need
# Distributions.jl package which is not part of the standard library.

# ## Loading packages

# If not in the REPL:

using Pkg                        # built-in package manager
Pkg.activate("env", shared=true) # create a new pkg env

# Add some packages to your environment (latest compatible versions
# added by default):

Pkg.add("Distributions")
Pkg.add("Plots")

# To load the code for use:

using Distributions
using Plots

N = 1000
samples = rand(Normal(), N);   # equivalent to Julia's built-in `randn(d)`
samples = (samples).^2;        # square element-wise

#-

g = fit(Gamma, samples)

#-

@show mean(g) median(g) pdf(g, 1)

# ## Plotting

f(x) = pdf(g, x)

xs = 0:0.1:4 # floats from 0 to 4 in steps of 0.1
ys = f.(xs)  # apply f element-wise to xs

plot(xs, ys, xrange=(0,4), yrange=(0,0.2))
histogram!(samples , normalize=true, alpha=0.4)

#-

savefig("my_first_plot.png")

# # Exercises

# ## Exercise 1

# Write a function named `total` that adds the elements of its vector input.

# ## Exercise 2

# Generate a 1000 random samples from the standard normal
# distribution. Create a second such sample, and add the two samples
# point-wise.  Compute the (sample) mean and variance of the combined
# samples. In the same
# plot, show a frequency-normalized histogram of the combined samples
# and a plot of the pdf for normal distribution with zero mean and
# variance `2`.

# You can use `std` to compute the standard deviation and `sqrt` to
# compute square roots.

# ## Exercise 3

# The following shows that named tuples share some behaviour with dictionaries:

t = (x = 1, y = "cat", z = 4.5)
keys(t)

#-

t[:y]

#-

# Write a function called `dict` that converts a named tuple to an
# actual dictionary. You can create an empty dictionary using `Dict()`.


text:
module Higher

include("rcd.jl")

end


text:
using Plots
gr()
using DataFrames,CSV

using LsqFit
r2=0

#%% Load the data
include("/home/kevin/Documents/1 - PhD Projects/2 - Confidence Attractor Network/ConfianceAttractorNetwork/src/decision_making.jl")
#include("/home/kevin/Documents/1 - PhD Projects/2 - Confidence Attractor Network/src/reaction_time_analysis_functions.jl")
include("/home/kevin/Documents/1 - PhD Projects/2 - Confidence Attractor Network/ConfianceAttractorNetwork/src/distributions.jl")

function sim_mf_decisionmaking_pattern_fastv2(coh_temp::Array{Float64},Ir::Float64,If::Float64,
  param::Dict{String,Float64},fixed_param::Dict{String,Float64},simu_param::Dict{String,Float64})

  # FI curve parameters - adapted from Wanf fit
  a=fixed_param["a"]
  b=fixed_param["b"]
  d=fixed_param["d"]
  falset= param["false_time"]
  thr = param["Threshold"]

  # time variables of the  simulation

  dt=simu_param["dt"]
  time_wind=simu_param["time_wind"] # mean window of the variables
  slide_wind = simu_param["slide_wind"]
  false_time=param["false_time"]
  falsetIf=param["false_timeIf"]

  #---- Intialise and vectorise variables to be used in loops below ------

  # Final Variables
  res=zeros(length(coh_temp))
  reaction_time=zeros(length(coh_temp))
  diffS=zeros(length(coh_temp))
  S1=zeros(length(coh_temp))
  S2=zeros(length(coh_temp))
  R1=zeros(length(coh_temp))
  R2=zeros(length(coh_temp))

  ###### Simu variables aprameters
  noise_amp=param["noise_amp"]

  mu0=  param[  "mu0"]

  I0E1p=  param[  "I0E1"]
  I0E2p =  param[  "I0E2"]
  JN11 =  param[ "JN11"]
  JN22  = param[ "JN22"]
  JN12 =  param[ "JN12"]
  JN21 =  param[ "JN21"]
  JAext  = param[ "JAext"]
  Tnmda = param[  "Tnmda"]
  gamma = param[  "gamma"]
  Tstim = param[  "Tstim"]
  Tampa = param[  "Tampa"]
  coh=  param[ "coh"]

  false_timeIf=  param[  "false_timeIf"]

  #####

  # Temporary variables

    Isyn1=0.0
  Isyn2=0.0
  Trest=param["T_rest"]
  s1=0.1
  s2=0.1
  phi1=0.0
  phi2=0.0
  I_eta1=noise_amp * randn()
  I_eta2=noise_amp * randn()

  nu1=Float64[2.0]
  nu2=Float64[2.0]
  k=0 # variable for mod 4 in the mean
  th=false
  result=true

  for (j,coh) in enumerate(coh_temp) # loop on trials number

    th=false
    result=true
    t=0

    nu1=Float64[nu1[end]]
    nu2=Float64[nu2[end]]
    k=0
    while th==false # while we do not cross the threshold
      t=t+dt
      k=k+1

      if t<Trest && result==true
        I0E1=I0E1p -Ir*exp(-t/falset)
        I0E2=I0E2p -Ir*exp(-t/falset)

        I_stim1=0
        I_stim2=0

      elseif t<Trest
        I0E1=I0E1p -Ir*exp(-t/falset) -If*exp(-t/falsetIf)
        I0E2=I0E2p -Ir*exp(-t/falset) -If*exp(-t/falsetIf)
        I_stim1=0
        I_stim2=0

      elseif t<Trest+dt/2 && t>Trest-dt/2
        I0E1=I0E1p
        I0E2=I0E2p
        I_stim1 = JAext*mu0*(1+coh/100)
        I_stim2 = JAext*mu0*(1-coh/100)

      else
        I0E1=I0E1p
        I0E2=I0E2p
        I_stim1 = JAext*mu0*(1+coh/100)
        I_stim2 = JAext*mu0*(1-coh/100)
      end

      Isyn1=JN11*s1-JN12*s2+I_stim1 + I_eta1
      Isyn2=JN22*s2-JN21*s1+I_stim2 + I_eta2

      phi1=(a*Isyn1-b)/(1-exp(-d*(a*Isyn1-b)))
      phi2=(a*Isyn2-b)/(1-exp(-d*(a*Isyn2-b)))

      #---- Dynamical equations -------------------------------------------
      # Mean NMDA-mediated synaptic dynamics updating
      s1 = s1 + dt*(-s1/Tnmda + (1-s1)*gamma*nu1[end]/1000);

      s2 = s2 + dt*(-s2/Tnmda + (1-s2)*gamma*nu2[end]/1000);

      # Ornstein-Uhlenbeck generation of noise in pop1 and 2
      I_eta1 = I_eta1 + (dt/Tampa)*(I0E1-I_eta1) + sqrt(dt/Tampa)*noise_amp*randn() ;
      I_eta2 = I_eta2 + (dt/Tampa)*(I0E2-I_eta2) + sqrt(dt/Tampa)*noise_amp*randn() ;

      # To ensure firing rates are always positive. Large noise amplitude
      # may result in unwanted negative values
      if phi1 < 0
        push!(nu1,0)
        phi1 = 0;
      else
        push!(nu1,phi1)
      end;
      if phi2 < 0
        push!(nu2,0)
        phi2 = 0;
      else
        push!(nu2,phi2)
      end;

      if mod(k,4)==0

        # Test the crossing of a threshold
        @views m1= mean(nu1[k-3:k])
        @views m2 = mean(nu2[k-3:k])
        if  t>Trest&& m1>thr
          th=true
          if coh>=0
          @inbounds  res[j]=1
            result=true
          else
            @inbounds  res[j]=-1
            result=false
          end

        elseif t>Trest&& m2>thr
          th = true
          if coh>=0
            @inbounds  res[j]=-1
            result=false
          else
            @inbounds  res[j]=1
            result=true
          end

        end
      end
    end
  @inbounds  reaction_time[j]=t-Trest
  @inbounds    diffS[j]=s2-s1
  @inbounds    R1[j]=nu1[end]
  @inbounds  R2[j]=nu2[end]
  @inbounds    S1[j]=s1
  @inbounds  S2[j]=s2
    th=false
  end

  result=DataFrame()
  result[:RTs]=reaction_time
  result[:NotError]=res
  result[:Distri]=coh_temp
  result[:DiffS]=diffS
  result[:S1]=S1
  result[:S2]=S2
  result[:R1]=R1
  result[:R2]=R2
  result[:LastTrial]=vcat([0],res[1:end-1])
  result[:LastCoh]=vcat([0],coh_temp[1:end-1])
  return result
end

function list_elements(l)
      a = Float64[]
        i=1
        while (length(a) <20)
           if any(a.==l[i])

            else
                push!(a,l[i])
            end
            i=i+1
        end
        return sort(a)
    end

    using Distributions, StatPlots,Colors

colorsT=[:red,:blue,:darkorange,:violet,:green]
colorsT=colormap("Blues",5)
o=4

param0=zeros(40)
paramR0=zeros(40)
paramR1=zeros(40)
param1=zeros(40)
xAcc=zeros(160)
xTr=zeros(160)
f=0
u=1
fig=plot(legend=false)

for Tr in [5000]#[200,500,1500,2500,5000]
    o=o
    c_min=-51
     c_max = 51
     number = 20

       distri = make_distributions(c_min,c_max,number,10000)

     param["T_rest"] = 1000
     Ir=0.035
     If=0.0

      r2= sim_mf_decisionmaking_pattern_fastv2(distri,Ir,If,
        param,fixed_param,simu_param)
        j=0
Accrep=Float64[]
Accalt=Float64[]
for a in list_elements(r2[:Distri])
    j=j+1

for i=2:10000
    if r2[:Distri][i]==a
        if r2[:NotError][i-1]*sign(r2[:Distri][i-1])==r2[:NotError][i]*sign(r2[:Distri][i])

            push!(Accrep,r2[:RTs][i])
        else
            push!(Accalt,r2[:RTs][i])
        end
    end
end

end
histogram!(Accalt,color=:forestgreen,alpha=0.4,bins=100,normed=true)
histogram!(Accrep,color=:darkorange,alpha=0.4,bins=100,normed=true)
print(mean(Accalt)-mean(Accrep))

#savefig(fig,"./LeftRightIndecAcc$Tr.svg")

end
fig
xlims!((0,1000))
savefig(fig,"./Hist0355000.svg")

fig

#### sauvegarde RSI energy package

# /home/kevin/Documents/1 - PhD Projects/1 - Sequential Decision making/3-SequentialEffects/2-Data/testsAltRep
param = Dict(
    "mu0"=>30.0,
    "noise_amp"=>0.02,
    "I0E1"=>0.3255,
    "I0E2"=>0.3255,
    "JN11"=>0.2609,
    "JN22"=>0.2609,
    "JN12"=>0.0497,
    "JN21"=>0.0497,
    "JAext" => 0.00052,
    "Tnmda" => 100.0,
    "Tampa" => 2.0,
    "gamma" => 0.641,
    "Tstim"=>5000.0,
    "T_rest"=>700.0,
    "coh"=>5.0,
    "false_time"=>200.0,
    "false_timeIf"=>500.0,
    "Threshold"=>20.0
    )
    fixed_param=Dict("a"=>270.0,"b"=>108.0,"d"=>0.1540)
  simu_param=Dict("dt"=>0.5,"time_wind"=>4.0,"slide_wind"=>4.0) #4=2/dt
Icdlist=[0.035,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.2,0.3]
#Icdlist=[0.2,0.3]#
Trlist=[500,1000,1500,2000,2500,5000,10000]

for Icd in Icdlist
  for Tr in Trlist

    c_min=-51
     c_max = 51
     number = 20

       distri = make_distributions(c_min,c_max,number,10000)

     param["T_rest"] = Tr
     Ir=Icd
     If=0.0

      r2= sim_mf_decisionmaking_pattern_fastv2(distri,Ir,If,
        param,fixed_param,simu_param)
        CSV.write("/home/kevin/Documents/1 - PhD Projects/1 - Sequential Decision making/3-SequentialEffects/2-Data/testsAltRep/EnergyTest/10000RepAlt-T$Tr-Ir$Icd.csv",r2)

  end
end


text:
import YAML

function get_sim_params(ARGS)

    path = ARGS[1]
    data = YAML.load(open(path))

    key = "grid"
    Nx = data[key]["nx"]
    Ny = data[key]["ny"]

    key = "numerics"
    NPROCX = data[key]["nprocx"]
    NPROCY = data[key]["nprocy"]
    MAX_STEPS = data[key]["max_steps"]
    Dt = data[key]["dt"]
    TOL = data[key]["tol"]

    key = "output"
    filename = data[key]["filename"]
    workdir = dirname(path)
    output_path = joinpath(workdir, basename(filename))

    println()
    println("Loading: ", path)
    println("Output: ", output_path)

    return Int64[Nx, Ny, NPROCX, NPROCY, MAX_STEPS], Float64[Dt, TOL], output_path

end


text:
using SMARTBackendServices
using Test

import JSONWebTokens

include("test_private_key/create.jl")

@testset "SMARTBackendServices.jl" begin
    include("basic.jl")
    include("jwt.jl")
end

include("test_private_key/delete.jl")


text:
using CUDAdrv
using CUDArt
include("CUBLAS.jl-0.0.2/src/CUBLAS.jl")
# using "CUBLAS.jl-0.0.2/src/CUBLAS"
CUDA_enabled = false

#cuda
export BasisTypeCuda, CUDA_store_basis, CUDA_calculate_basis, complexBasis_choice, CUDA_init, FullCorrelation_parallized, new_energy

function CUDA_init()
    CUBLAS_init()
    CUDArt.init(0)
    CUDArt.device(0)
    global md = CUDAdrv.CuModuleFile("$(ENV["THREEPHOTONS_PATH"])/src/cuda_kernel.ptx")
    global calculate_coefficient_matrix_cuda = CUDAdrv.CuFunction(md, "calculate_coefficient_matrix")
    global energy_cuda = CUDAdrv.CuFunction(md, "energy")

    global CUDA_enabled = true
    CUDAdrv.gc()
    println("Initialization of CUDA complete.")
end

"""Datatype for precalculated three photon correlation basis function"""
type BasisTypeCuda <: AbstractBasisType
    d_wignerlist::CudaArray
    d_indices::CudaArray
    d_PAcombos::CudaArray
    d_B::CudaArray
    h_P::Array{Float32,2}
    d_correlation::CudaArray

    basislen::Int64
    N::Int64
    L::Int64
    LMAX::Int64
    lrange::StepRange
    ctr::Dict
    rtc::Dict
    K::Int64
    lambda::Float64
    dq::Float64
end

"""Overrides the original complexBasis function"""
function CUDA_store_basis(basis::BasisType)
    d_wignerlist = CudaArray(convert(Array{Float32}, sdata(basis.wignerlist)))
    d_indices = CudaArray(convert(Array{Int32}, sdata(basis.indices)))
    d_PAcombos = CudaArray(convert(Array{Int32}, sdata(basis.PAcombos)))
    d_B = CudaArray(sdata(basis.B))
    d_correlation = CudaArray(Float32,(2*basis.N^2,round(Int64, basis.K*(basis.K+1)*(basis.K+2)/6)))

    cuda_basis =  BasisTypeCuda(d_wignerlist, d_indices, d_PAcombos, d_B, basis.h_P, d_correlation, basis.basislen, basis.N, basis.L, basis.LMAX, basis.lrange, basis.ctr, basis.rtc, basis.K, basis.lambda, basis.dq)
    println("Initialized CUDA basis.")
    return cuda_basis
end

function CUDA_calculate_basis( L::Int64, LMAX::Int64, N::Int64, K::Int64, lambda::Float64, dq::Float64, forIntensity=true)
    println("Intialize basis for CUDA.")
    basis = calculate_basis(L, LMAX, N, K, lambda, dq, forIntensity)
    if CUDA_enabled
        return CUDA_store_basis(basis)
    else
        return basis
    end
end

function CUDA_calculate_basis( L::Int64, LMAX::Int64, N::Int64, K::Int64, forIntensity=true)
    return CUDA_calculate_basis(L,LMAX,N,K, 0.0, 0.0, forIntensity)
end

global complexBasis_choice = CUDA_calculate_basis

"""Calculates the full three photon correlation"""
function FullCorrelation_parallized(intensity::SphericalHarmonicsVolume, basis::BasisTypeCuda, minimal::Bool=true, normalize::Bool=false, return_raw::Bool=false)
    # println("CUDA correlation.")
    #Prepare all arrays on GPU
    klength = Integer(basis.K*(basis.K+1)*(basis.K+2)/6)
    numcoeff = num_coeff(basis.LMAX)
    d_coeff = CudaArray(Complex{Float32}[intensity.coeff[k][i] for k = 1:basis.K, i=1:numcoeff]')
    d_PA = CudaArray(basis.h_P)

    threadsperblock = 1024
    blockspergrid = ceil(Int32, Base.size(basis.d_PAcombos)[2] / threadsperblock)
    CUDAdrv.cudacall(calculate_coefficient_matrix_cuda, blockspergrid, threadsperblock,
    (Ptr{Cfloat}, Cint, Ptr{Cfloat}, Ptr{Cint}, Cint, Ptr{Cint}, Cint, Ptr{Cfloat}, Cint),
    d_coeff, numcoeff, basis.d_wignerlist, basis.d_indices, Base.size(basis.d_indices)[2], basis.d_PAcombos, Base.size(basis.d_PAcombos)[2], d_PA, klength)

    gemm!('N','N',Float32(1.0), basis.d_B, d_PA, Float32(0.0), basis.d_correlation)

    #Wait for computation to finish and transfer result to host
    res = to_host(basis.d_correlation)

    #Enforce GC to avoid crashes
    CUDArt.free(d_PA)
    CUDArt.free(d_coeff)
    # CUDAdrv.gc()

    if return_raw return res end

    #Translate matrix into 5dimensional correlation format
    t = zeros(Float64, basis.N, 2*basis.N, basis.K, basis.K, basis.K)
    i = 1
    for k1=1:basis.K for k2=1:k1 for k3=1:k2
        t[:,:, k3, k2, k1] = reshape(res[:, i], basis.N, 2*basis.N)
        i += 1
    end end end
    return normalize ? t / sum(abs, t) : t
end

# """Calculate the energy on the GPU to avoid data transfer of large correlation Matrix
# However, this function is not much faster than the original.
# Bottlenck is still the large matrix operation."""
# function new_energy(intensity::SphericalHarmonicsVolume, basis::BasisTypeCuda, c3ref::C3, K3_range::UnitRange{Int64}, measure::String="Bayes", negativity_factor::Float64=0.0)
#     # println("CUDA correlation.")
#     #Prepare all arrays on GPU
#     klength = Integer(basis.K*(basis.K+1)*(basis.K+2)/6)
#     numcoeff = num_coeff(basis.LMAX)
#     d_coeff = CudaArray(Complex{Float32}[intensity.coeff[k][i] for k = 1:basis.K, i=1:numcoeff]')
#     d_PA = CudaArray(basis.h_P)
#
#     threadsperblock = 1024
#     blockspergrid = ceil(Int32, Base.size(basis.d_PAcombos)[2] / threadsperblock)
#     CUDAdrv.cudacall(calculate_coefficient_matrix_cuda, blockspergrid, threadsperblock, (d_coeff, numcoeff, basis.d_wignerlist, basis.d_indices, Base.size(basis.d_indices)[2], basis.d_PAcombos, Base.size(basis.d_PAcombos)[2], d_PA, klength))
#
#     gemm!('N','N',Float32(1.0), basis.d_B, d_PA, Float32(0.0), basis.d_correlation)
#
#     Kcombos = [(k1,k2,k3) for k1 in K3_range for k2=minimum(K3_range):k1 for k3=minimum(K3_range):k2]
#     Kcombos_length = length(Kcombos)
#     Kcombos_matrix = Int32[Kcombos[i][j] for i = 1:length(Kcombos), j=1:3]
#
#     d_Kcombos = CudaArray(Kcombos_matrix')
#     d_result = CudaArray(Float32, Kcombos_length)
#     c3ref_compressed = Array(Float32, 2*basis.N^2, Kcombos_length)
#
#     for i=1:Kcombos_length
#         k1,k2,k3 = Kcombos[i]
#         c3ref_compressed[:,i] = reshape(c3ref[:,:, k3,k2,k1], 2*basis.N^2)
#     end
#     d_c3ref = CudaArray(c3ref_compressed')
#
#     #  __global__ void energy(const float *c3, const float *c3ref, const int* Kcombos, const int combolength, const int slab_size, float *result)
#     threadsperblock = 1024
#     blockspergrid = ceil(Int32, Kcombos_length/threadsperblock)
#
#     launch(energy_cuda, blockspergrid, threadsperblock, (basis.d_correlation, d_c3ref, d_Kcombos, Kcombos_length, 2*basis.N^2, d_result))
#
#     return sum(abs, to_host(d_result))
# end


text:
export clusterize,
       KMeansClusterizer,
       YuEigenvectorRotation,
       EigenvectorClusterizer,
       EigenvectorClusteringResult

using Clustering

import Clustering.assignments
import Clustering.ClusteringResult
struct EigenvectorClusteringResult{T<:Integer} <: ClusteringResult
    assignments::Vector{T}
end
function assignments(r::EigenvectorClusteringResult)
    return r.assignments
end
abstract type EigenvectorClusterizer end
"""
```julia
struct KMeansClusterizer <: EigenvectorClusterizer
    k::Integer
    init::Symbol
end
```
"""
struct KMeansClusterizer <: EigenvectorClusterizer
    k::Integer
    init::Symbol
    maxiter::Integer
end
function KMeansClusterizer(k::Integer)
    return KMeansClusterizer(k, :kmpp, 50000)
end
function clusterize(t::KMeansClusterizer, E)
    model = kmeans(Matrix(E'), t.k, init =t.init, maxiter=t.maxiter)

    return EigenvectorClusteringResult(assignments(model))
end

"""
Multiclass Spectral Clustering
"""
struct YuEigenvectorRotation <: EigenvectorClusterizer
    maxIter::Integer
end
function YuEigenvectorRotation()
    return YuEigenvectorRotation(500)
end
function clusterize(cfg::YuEigenvectorRotation, X_star_hat::Matrix)
     (N,k) = size(X_star_hat)
     X_star_hat = spdiagm(0=>1 ./ sqrt.(vec(mapslices(norm, X_star_hat, dims=[2]))))*X_star_hat
     hasConverged = false
     R_star = zeros(k,k)
     R_star[:, 1] = [X_star_hat[rand(1:N), i] for i = 1:k ]
     c = zeros(N)
     for j=2:k
        c = c + abs.(X_star_hat*R_star[:,j-1])
        i = findmin(c)[2]
        R_star[:, j] = X_star_hat[i, :]'
     end
     lastObjectiveValue = Inf
     nIter = 0
     ncut_value = 0
     X_star = nothing
     while !hasConverged
        nIter = nIter+ 1
        X_hat = X_star_hat*R_star
        #non maximum supression
        labels = vec([I[2] for I in findmax(X_hat, dims=2)[2]])
        X_star = zeros(size(X_star_hat))
        for (i, l) = enumerate(labels)
            X_star[i, l] = l
        end
        F = svd(X_star' * X_star_hat, full=true)
        U, S, Vh = (F.U, F.S, F.Vt)
        ncutValue = sum(S)
        if ((abs(ncutValue - lastObjectiveValue) < eps()) || (nIter > cfg.maxIter))
            hasConverged = true
        else
            lastObjectiveValue = ncutValue
            R_star = Vh'*U'
        end
    end
    labels = vec([I[2] for I in findmax(X_star, dims=2)[2]])
    return EigenvectorClusteringResult(labels)
end
"""
```julia
function clusterize{T<:EigenvectorEmbedder, C<:EigenvectorClusterizer}(cfg::T, clus::C, X)
```

Given a set of patterns `X` generates an eigenvector space according to `T<:EigenvectorEmbeddder` and then clusterize the eigenvectors using the algorithm defined
by `C<:EigenvectorClusterize`.

"""
function clusterize(cfg::T, clus::C, X, params...) where T<:EigenvectorEmbedder where C<:EigenvectorClusterizer
    E = embedding(cfg,X, params...)
    return clusterize(clus, E)
end


text:
# Model interconnections

@doc """
`series(s1::LTISystem, s2::LTISystem)`

Connect systems in series, equivalent to `s2*s1`
""" ->
series(s1::LTISystem, s2::LTISystem) = s2*s1

@doc """
`series(s1::LTISystem, s2::LTISystem)`

Connect systems in parallel, equivalent to `s2+s1`
""" ->
parallel(s1::LTISystem, s2::LTISystem) = s1 + s2

append() = LTISystem[]
@doc """
`append(systems::StateSpace...), append(systems::TransferFunction...)`

Append systems in block diagonal form
""" ->
function append(systems::StateSpace...)
    Ts = systems[1].Ts
    if !all([s.Ts == Ts for s in systems])
        error("Sampling time mismatch")
    end
    A = blkdiag([s.A for s in systems]...)
    B = blkdiag([s.B for s in systems]...)
    C = blkdiag([s.C for s in systems]...)
    D = blkdiag([s.D for s in systems]...)
    states = vcat([s.statenames for s in systems]...)
    outputs = vcat([s.outputnames for s in systems]...)
    inputs = vcat([s.inputnames for s in systems]...)
    return StateSpace(A, B, C, D, Ts, states, inputs, outputs)
end

function append(systems::TransferFunction...)
    Ts = systems[1].Ts
    if !all([s.Ts == Ts for s in systems])
        error("Sampling time mismatch")
    end
    mat = blkdiag([s.matrix for s in systems]...)
    inputs = vcat([s.inputnames for s in systems]...)
    outputs = vcat([s.outputnames for s in systems]...)
    return TransferFunction(mat, Ts, inputs, outputs)
end

append(systems::LTISystem...) = append(promote(systems...)...)

#This is needed until julia removes deprecated vect()
function Base.vect(sys::Union{LTISystem, Real}...)
    T = Base.promote_typeof(sys...)
    copy!(Array(T,length(sys)), sys)
end

function Base.vcat(systems::StateSpace...)
    # Perform checks
    nu = systems[1].nu
    if !all([s.nu == nu for s in systems])
        error("All systems must have same input dimension")
    end
    Ts = systems[1].Ts
    if !all([s.Ts == Ts for s in systems])
        error("Sampling time mismatch")
    end
    A = blkdiag([s.A for s in systems]...)
    B = vcat([s.B for s in systems]...)
    C = blkdiag([s.C for s in systems]...)
    D = vcat([s.D for s in systems]...)
    states = vcat([s.statenames for s in systems]...)
    outputs = vcat([s.outputnames for s in systems]...)
    inputs = systems[1].inputnames
    if !all([s.inputnames == inputs for s in systems])
        inputs = fill(String(""),size(inputs, 1))
    end
    return StateSpace(A, B, C, D, Ts, states, inputs, outputs)
end

function Base.vcat(systems::TransferFunction...)
    # Perform checks
    nu = systems[1].nu
    if !all([s.nu == nu for s in systems])
        error("All systems must have same input dimension")
    end
    Ts = systems[1].Ts
    if !all([s.Ts == Ts for s in systems])
        error("Sampling time mismatch")
    end
    mat = vcat([s.matrix for s in systems]...)
    outputs = vcat([s.outputnames for s in systems]...)
    inputs = systems[1].inputnames
    if !all([s.inputnames == inputs for s in systems])
        inputs = fill(String(""),size(inputs, 1))
    end
    return TransferFunction(mat, Ts, inputs, outputs)
end

Base.vcat(systems::LTISystem...) = vcat(promote(systems...)...)

function Base.vcat{T<:Real}(systems::Union{VecOrMat{T},T,TransferFunction}...)
    if promote_type(map(e->typeof(e),systems)...) <: TransferFunction
        vcat(map(e->convert(TransferFunction,e),systems)...)
    else
        cat(1,systems...)
    end
end

function Base.hcat(systems::StateSpace...)
    # Perform checks
    ny = systems[1].ny
    if !all([s.ny == ny for s in systems])
        error("All systems must have same output dimension")
    end
    Ts = systems[1].Ts
    if !all([s.Ts == Ts for s in systems])
        error("Sampling time mismatch")
    end
    A = blkdiag([s.A for s in systems]...)
    B = blkdiag([s.B for s in systems]...)
    C = hcat([s.C for s in systems]...)
    D = hcat([s.D for s in systems]...)
    states = vcat([s.statenames for s in systems]...)
    inputs = vcat([s.inputnames for s in systems]...)
    outputs = systems[1].outputnames
    if !all([s.outputnames == outputs for s in systems])
        outputs = fill(String(""),size(outputs, 1))
    end
    return StateSpace(A, B, C, D, Ts, states, inputs, outputs)
end

function Base.hcat(systems::TransferFunction...)
    # Perform checks
    ny = systems[1].ny
    if !all([s.ny == ny for s in systems])
        error("All systems must have same output dimension")
    end
    Ts = systems[1].Ts
    if !all([s.Ts == Ts for s in systems])
        error("Sampling time mismatch")
    end
    mat = hcat([s.matrix for s in systems]...)
    inputs = vcat([s.inputnames for s in systems]...)
    outputs = systems[1].outputnames
    if !all([s.outputnames == outputs for s in systems])
        outputs = fill(String(""),size(outputs, 1))
    end
    return TransferFunction(mat, Ts, inputs, outputs)
end

Base.hcat(systems::LTISystem...) = hcat(promote(systems...)...)

function Base.hcat{T<:Real}(systems::Union{T,VecOrMat{T},TransferFunction}...)
    if promote_type(map(e->typeof(e),systems)...) <: TransferFunction
        hcat(map(e->convert(TransferFunction,e),systems)...)
    else
        cat(2,systems...)
    end
end

# Empty definition to get rid of warning
Base.blkdiag() = []
function Base.blkdiag(mats::Matrix...)
    rows = Int[size(m, 1) for m in mats]
    cols = Int[size(m, 2) for m in mats]
    T = eltype(mats[1])
    for ind=1:length(mats)
        T = promote_type(T, eltype(mats[ind]))
    end
    res = zeros(T, sum(rows), sum(cols))
    m = 1
    n = 1
    for ind=1:length(mats)
        mat = mats[ind]
        i = rows[ind]
        j = cols[ind]
        res[m:m + i - 1, n:n + j - 1] = mat
        m += i
        n += j
    end
    return res
end


text:
_shared_vector_warning = """
!!! warning
    The vector passed to the inner reducing function is valid only
    during its _immediate_ reduction step.  It must be reduced
    immediately _or_ copied.
"""

_shared_notes_unfold = """
Note that input is ignored.  To use the input in the downstream
reduction steps, use [`TeeZip`](@ref) or [`Zip`](@ref).
"""

# https://clojure.github.io/clojure/clojure.core-api.html#clojure.core/map
# https://clojuredocs.org/clojure.core/map
"""
    Map(f)

Apply unary function `f` to each input and pass the result to the
inner reducing step.

$(_thx_clj("map"))

# Examples
```jldoctest
julia> using Transducers

julia> collect(Map(x -> 2x), 1:3)
3-element Array{Int64,1}:
 2
 4
 6
```
"""
struct Map{F} <: Transducer
    f::F
end

outtype(xf::Map, intype) = Union{Base.return_types(xf.f, (intype,))...}
next(rf::R_{Map}, result, input) = next(rf.inner, result, rf.xform.f(input))

# https://clojure.github.io/clojure/clojure.core-api.html#clojure.core/replace
# https://clojuredocs.org/clojure.core/replace
"""
    Replace(dict)

Replace each input with the value in the dictionary `dict` if it
matches with a key.  Otherwise output the input as-is.

$(_thx_clj("replace"))

# Examples
```jldoctest
julia> using Transducers

julia> collect(Replace(Dict('a' => 'A')), "abc")
3-element Array{Char,1}:
 'A'
 'b'
 'c'
```
"""
struct Replace{D} <: Transducer
    d::D  # dictionary-like object
end

outtype(xf::Replace, intype) = Union{intype, valtype(xf.d)}
next(rf::R_{Replace}, result, input) =
    next(rf.inner, result, get(rf.xform.d, input, input))

# https://clojure.github.io/clojure/clojure.core-api.html#clojure.core/cat
# https://clojuredocs.org/clojure.core/cat
"""
    Cat()

Concatenate/flatten nested iterators.

$(_thx_clj("cat"))

# Examples
```jldoctest
julia> using Transducers

julia> collect(Cat(), [[1, 2], [3], [4, 5]]) == 1:5
true
```
"""
struct Cat <: Transducer
end

outtype(::Cat, intype) = ieltype(intype)

# Inner transducer has to be started once the input is known.  That's
# why `start` for `Cat` has to bail out immediately; i.e., it's not a
# bug that `start(rf.inner, iresult)` is not called here:
start(rf::R_{Cat}, result) = wrap(rf, Unseen(), result)

next(rf::R_{Cat}, result, input) =
    wrapping(rf, result) do istate0, iresult
        if istate0 isa Unseen
            istate1 = start(rf.inner, iresult)
        else
            istate1 = istate0
        end
        istate2 = __foldl__(rf.inner, istate1, input, nocomplete)
        istate2, istate2
    end

# https://clojure.github.io/clojure/clojure.core-api.html#clojure.core/mapcat
# https://clojuredocs.org/clojure.core/mapcat
"""
    MapCat(f)

Concatenate output of `f` which is expected to return an iterable.

$(_thx_clj("mapcat"))

# Examples
```jldoctest
julia> using Transducers

julia> collect(MapCat(x -> 1:x), 1:3)
6-element Array{Int64,1}:
 1
 1
 2
 1
 2
 3
```
"""
const MapCat = Composition{<:Map, <:Cat}

MapCat(f) = Map(f) |> Cat()

# https://clojure.github.io/clojure/clojure.core-api.html#clojure.core/filter
# https://clojuredocs.org/clojure.core/filter
"""
    Filter(pred)

Skip items for which `pred` is evaluated to `false`.

$(_thx_clj("filter"))

# Examples
```jldoctest
julia> using Transducers

julia> collect(Filter(iseven), 1:3)
1-element Array{Int64,1}:
 2
```
"""
struct Filter{P} <: AbstractFilter
    pred::P
end

next(rf::R_{Filter}, result, input) =
    rf.xform.pred(input) ? next(rf.inner, result, input) : result

# https://clojure.github.io/clojure/clojure.core-api.html#clojure.core/take
# https://clojuredocs.org/clojure.core/take
"""
    Take(n)

Take `n` items from the input sequence.

$(_thx_clj("take"))

# Examples
```jldoctest
julia> using Transducers

julia> collect(Take(2), 1:10)
2-element Array{Int64,1}:
 1
 2
```
"""
struct Take <: AbstractFilter
    n::Int
end

start(rf::R_{Take}, result) = wrap(rf, rf.xform.n, result)

next(rf::R_{Take}, result, input) =
    wrapping(rf, result) do n, iresult
        if n > 0
            iresult = next(rf.inner, iresult, input)
            n -= 1
        end
        if n <= 0
            iresult = ensure_reduced(iresult)
        end
        return n, iresult
    end

# https://clojure.github.io/clojure/clojure.core-api.html#clojure.core/take-while
# https://clojuredocs.org/clojure.core/take-while
"""
    TakeWhile(pred)

Take items while `pred` returns `true`.  Abort the transducible
process when `pred` returns `false` for the first time.

$(_thx_clj("take-while"))

# Examples
```jldoctest
julia> using Transducers

julia> collect(TakeWhile(x -> x < 3), 1:10)
2-element Array{Int64,1}:
 1
 2
```
"""
struct TakeWhile{P} <: AbstractFilter
    pred::P
end

next(rf::R_{TakeWhile}, result, input) =
    rf.xform.pred(input) ? next(rf.inner, result, input) : ensure_reduced(result)

# https://clojure.github.io/clojure/clojure.core-api.html#clojure.core/take-nth
# https://clojuredocs.org/clojure.core/take-nth
"""
    TakeNth(n)

Output every `n` item to the inner reducing step.

$(_thx_clj("take-nth"))

# Examples
```jldoctest
julia> using Transducers

julia> collect(TakeNth(3), 1:9)
3-element Array{Int64,1}:
 1
 4
 7
```
"""
struct TakeNth <: AbstractFilter
    n::Int
end

start(rf::R_{TakeNth}, result) = wrap(rf, rf.xform.n, result)

next(rf::R_{TakeNth}, result, input) =
    wrapping(rf, result) do c, iresult
        if c == rf.xform.n
            iresult = next(rf.inner, iresult, input)
            c = 1
        else
            c += 1
        end
        return c, iresult
    end

# https://clojure.github.io/clojure/clojure.core-api.html#clojure.core/drop
# https://clojuredocs.org/clojure.core/drop
"""
    Drop(n)

Drop first `n` items.

$(_thx_clj("drop"))

# Examples
```jldoctest
julia> using Transducers

julia> collect(Drop(3), 1:5)
2-element Array{Int64,1}:
 4
 5
```
"""
struct Drop <: AbstractFilter
    n::Int
end

start(rf::R_{Drop}, result) = wrap(rf, 0, result)

next(rf::R_{Drop}, result, input) =
    wrapping(rf, result) do c, iresult
        if c >= rf.xform.n
            c, next(rf.inner, iresult, input)
        else
            c += 1
            c, iresult
        end
    end

# https://clojure.github.io/clojure/clojure.core-api.html#clojure.core/drop-last
# https://clojuredocs.org/clojure.core/drop-last
"""
    DropLast(n)

Drop last `n` items.

$(_thx_clj("drop-last"))

# Examples
```jldoctest
julia> using Transducers

julia> collect(DropLast(2), 1:5)
3-element Array{Int64,1}:
 1
 2
 3
```
"""
struct DropLast <: AbstractFilter
    n::Int

    function DropLast(n)
        @assert n > 0
        return new(n)
    end
end

function start(rf::R_{DropLast}, result)
    n = rf.xform.n + 1
    return wrap(rf, (-n, Vector{InType(rf)}(undef, n)), result)
end

next(rf::R_{DropLast}, result, input) =
    wrapping(rf, result) do (c, buffer), iresult
        c += 1
        n = length(buffer)
        if c <= 0
            buffer[c + n] = input
            if c == 0
                (c, buffer), next(rf.inner, iresult, buffer[1])
            else
                (c, buffer), iresult
            end
        elseif c >= n
            buffer[c] = input
            (0, buffer), next(rf.inner, iresult, buffer[1])
        else
            buffer[c] = input
            (c, buffer), next(rf.inner, iresult, buffer[c + 1])
        end
    end

# https://clojure.github.io/clojure/clojure.core-api.html#clojure.core/drop-while
# https://clojuredocs.org/clojure.core/drop-while
"""
    DropWhile(pred)

Drop items while `pred` returns `true` consecutively.  It becomes a
no-op after `pred` returns a `false`.

$(_thx_clj("drop-while"))

# Examples
```jldoctest
julia> using Transducers

julia> collect(DropWhile(x -> x < 3), [1:5; 1:2])
5-element Array{Int64,1}:
 3
 4
 5
 1
 2
```
"""
struct DropWhile{F} <: AbstractFilter
    pred::F
end

start(rf::R_{DropWhile}, result) = wrap(rf, true, result)

next(rf::R_{DropWhile}, result, input) =
    wrapping(rf, result) do dropping, iresult
        if dropping
            dropping = rf.xform.pred(input)
            dropping && return (dropping, iresult)
        end
        dropping, next(rf.inner, iresult, input)
    end

# https://clojure.github.io/clojure/clojure.core-api.html#clojure.core/partition-all
# https://clojuredocs.org/clojure.core/partition-all
"""
    Window(size, step = size, flush = false)
    Window(size; step = size, flush = false)

Sliding window of width `size` and interval `step`.

$_shared_vector_warning

$(_thx_clj("partition-all"))

# Examples
```jldoctest
julia> using Transducers

julia> collect(Window(3) |> Map(copy), 1:8)
2-element Array{Array{Int64,1},1}:
 [1, 2, 3]
 [4, 5, 6]

julia> collect(Window(3; flush=true) |> Map(copy), 1:8)
3-element Array{Array{Int64,1},1}:
 [1, 2, 3]
 [4, 5, 6]
 [7, 8]

julia> collect(Window(3; step=1) |> Map(copy), 1:8)
6-element Array{Array{Int64,1},1}:
 [1, 2, 3]
 [2, 3, 4]
 [3, 4, 5]
 [4, 5, 6]
 [5, 6, 7]
 [6, 7, 8]
```
"""
struct Window <: Transducer
    size::Int
    step::Int
    flush::Bool
end

Window(size, step; flush = false) = Window(size; step = step, flush = flush)
Window(size; step = size, flush = false) = Window(size, step, flush)

outtype(::Window, intype) = DenseSubVector{intype}

function start(rf::R_{Window}, result)
    buf = Vector{InType(rf)}()
    sizehint!(buf, rf.xform.size)
    return wrap(rf, (0, 0, buf), result)
end

function next(rf::R_{Window}, result, input)
    wrapping(rf, result) do (i, s, buf), iresult
        _window_next(rf, i, s, buf, iresult, input)
    end
end

function _window_next(rf, i, s, buf, iresult, input)
    i += 1
    @assert 1 <= i <= rf.xform.size
    len = length(buf)
    if s == 0 && len < rf.xform.size
        push!(buf, input)
        if i == rf.xform.size
            # This is the first time `length(buf) == rf.xform.size` is
            # true.
            iinput = (@view buf[1:end]) :: DenseSubVector
            # Do not use `@view buf[:]` to make the indices
            # `UnitRange` instead of `OneTo`.
            @assert i == rf.xform.size
            i = 0
            s = rf.xform.step
            iresult = next(rf.inner, iresult, iinput)
            # Throw away buf[1] which is just consumed by the downstream:
            deleteat!(buf, 1)
        end
        return (i, s, buf), iresult
    elseif len < 2 * rf.xform.size - 1
        push!(buf, input)
    else
        # TODO: optimize for size < step case
        buf[i + rf.xform.size - 1] = input
    end
    if i == rf.xform.size
        # Wrapping the window.  Next window will be i=1.
        # Copy tail to head.
        copyto!(buf,
                1,                  # destination offset
                buf,
                rf.xform.size + 1,  # source offset
                rf.xform.size - 1)  # size
    end
    s -= 1
    @assert 0 <= s < rf.xform.step
    if s == 0
        s = rf.xform.step
        iinput = @view buf[i:i + rf.xform.size - 1] # unsafe_view? @inbounds?
        iinput :: DenseSubVector
        iresult = next(rf.inner, iresult, iinput)
    end
    if i == rf.xform.size
        i = 0
    end
    return (i, s, buf), iresult
end

function complete(rf::R_{Window}, result)
    (i, s, buf), iresult = unwrap(rf, result)
    if rf.xform.flush && s != rf.xform.step
        iinput = @view buf[i:i + rf.xform.size - 1] # unsafe_view? @inbounds?
        iinput = @view iinput[s + 1:end]
        iinput :: DenseSubVector
        iresult = next(rf.inner, iresult, iinput)
    end
    return complete(rf.inner, iresult)
end

# TODO: implement SVector version of Window

# https://clojure.github.io/clojure/clojure.core-api.html#clojure.core/partition-by
# https://clojuredocs.org/clojure.core/partition-by
"""
    PartitionBy(f)

Group input sequence into chunks in which `f` returns a same value
consecutively.

$_shared_vector_warning

$(_thx_clj("partition-by"))

# Examples
```jldoctest
julia> using Transducers

julia> collect(PartitionBy(x -> (x + 1) ÷ 3) |> Map(copy), 1:9)
4-element Array{Array{Int64,1},1}:
 [1]
 [2, 3, 4]
 [5, 6, 7]
 [8, 9]
```
"""
struct PartitionBy{F} <: Transducer
    f::F
end

struct Unseen end

outtype(::PartitionBy, intype) = Vector{intype}

function start(rf::R_{PartitionBy}, result)
    iinput = InType(rf)[]
    return wrap(rf, (iinput, Unseen()), start(rf.inner, result))
end

@inline function next(rf::R_{PartitionBy}, result, input)
    wrapping(rf, result) do (iinput, pval), iresult
        val = rf.xform.f(input)
        if pval isa Unseen || val == pval
            push!(iinput, input)
        else
            iresult = next(rf.inner, iresult, iinput)
            empty!(iinput)
            isreduced(iresult) || push!(iinput, input)
        end
        return (iinput, val), iresult
    end
end

function complete(rf::R_{PartitionBy}, ps)
    (iinput, _), iresult = unwrap(rf, ps)
    if !isempty(iinput)
        iresult = next(rf.inner, iresult, iinput)
    end
    return complete(rf.inner, iresult)
end

# https://clojure.github.io/clojure/clojure.core-api.html#clojure.core/keep
# https://clojuredocs.org/clojure.core/keep
"""
    Keep(f)

Pass non-`nothing` output of `f` to the inner reducing step.

$(_thx_clj("keep"))

# Examples
```jldoctest
julia> using Transducers

julia> xf = Keep() do x
           if x < 3
               x + 1
           end
       end;

julia> collect(xf, 1:5)
2-element Array{Int64,1}:
 2
 3
```
"""
struct Keep{F} <: Transducer
    f::F
end

outtype(xf::Keep, intype) =
    Core.Compiler.typesubtract(Union{Base.return_types(xf.f, (intype,))...},
                               Nothing)

function next(rf::R_{Keep}, result, input)
    iinput = rf.xform.f(input)
    return iinput === nothing ? result : next(rf.inner, result, iinput)
end

# https://clojure.github.io/clojure/clojure.core-api.html#clojure.core/distinct
# https://clojuredocs.org/clojure.core/distinct
"""
    Distinct()

Pass only unseen item to the inner reducing step.

$(_thx_clj("distinct"))

# Examples
```jldoctest
julia> using Transducers

julia> collect(Distinct(), [1, 1, 2, 1, 3, 3, 2])
3-element Array{Int64,1}:
 1
 2
 3
```
"""
struct Distinct <: AbstractFilter
end

function start(rf::R_{Distinct}, result)
    seen = Set(InType(rf)[])
    return wrap(rf, seen, start(rf.inner, result))
end

function next(rf::R_{Distinct}, result, input)
    wrapping(rf, result) do seen, iresult
        if input in seen
            return seen, iresult
        else
            push!(seen, input)
            return seen, next(rf.inner, iresult, input)
        end
    end
end

# https://clojure.github.io/clojure/clojure.core-api.html#clojure.core/interpose
# https://clojuredocs.org/clojure.core/interpose
"""
    Interpose(sep)

Interleave input items with a `sep`.

$(_thx_clj("interpose"))

# Examples
```jldoctest
julia> using Transducers

julia> collect(Interpose(missing), 1:3)
5-element Array{Union{Missing, Int64},1}:
 1
  missing
 2
  missing
 3
```
"""
struct Interpose{T} <: Transducer
    sep::T
end

outtype(xf::Interpose, intype) = Union{typeof(xf.sep), intype}

start(rf::R_{Interpose}, result) = wrap(rf, Val(true), start(rf.inner, result))

next(rf::R_{Interpose}, result, input) =
    wrapping(rf, result) do isfirst, iresult
        if isfirst isa Val{false}
            iresult = next(rf.inner, iresult, rf.xform.sep)
        end
        return Val(false), next(rf.inner, iresult, input)
    end

# https://clojure.github.io/clojure/clojure.core-api.html#clojure.core/dedupe
# https://clojuredocs.org/clojure.core/dedupe
"""
    Dedupe()

De-duplicate _consecutive_ items.

$(_thx_clj("dedupe"))

# Examples
```jldoctest
julia> using Transducers

julia> collect(Dedupe(), [1, 1, 2, 1, 3, 3, 2])
5-element Array{Int64,1}:
 1
 2
 1
 3
 2
```
"""
struct Dedupe <: AbstractFilter
end

start(rf::R_{Dedupe}, result) = wrap(rf, Unseen(), start(rf.inner, result))

next(rf::R_{Dedupe}, result, input) =
    wrapping(rf, result) do prev, iresult
        if prev isa Unseen || prev != input
            return input, next(rf.inner, iresult, input)
        else
            return input, iresult
        end
    end

"""
    Scan(f, [init])

Accumulate input with binary function `f` and pass the accumulated
result so far to the inner reduction step.

The inner reducing step receives the sequence `y₁, y₂, y₃, ..., yₙ, ...`
when the sequence `x₁, x₂, x₃, ..., xₙ, ...` is fed to `Scan(f)`.

    y₁ = f(x₁, init)
    y₂ = f(x₂, y₁)
    y₃ = f(x₃, y₂)
    ...
    yₙ = f(xₙ, yₙ₋₁)

This is a generalized version of the
[_prefix sum_](https://en.wikipedia.org/wiki/Prefix_sum) also known as
_cumulative sum_, _inclusive scan_, or _scan_.

Note that the associativity of `f` is not required when the transducer
is used in a process that gurantee an order, such as [`mapfoldl`](@ref).

Unless `f` is a function with known identity element such as `+`, `*`,
`min`, `max`, and `append!`, the initial state `init` must be
provided.

# Examples
```jldoctest
julia> using Transducers

julia> collect(Scan(*), 1:3)
3-element Array{Int64,1}:
 1
 2
 6

julia> collect(Map(x -> x + im) |> Scan(*), 1:3)
3-element Array{Complex{Int64},1}:
 1 + 1im
 1 + 3im
 0 + 10im

julia> collect(Scan(*, 10), 1:3)
3-element Array{Int64,1}:
 10
 20
 60
```
"""
struct Scan{F, T} <: Transducer
    f::F
    init::T
end

Scan(f) = Scan(f, nothing)

_lefttype(xf::Scan{<:Any, Nothing}, intype) = typeof(identityof(xf.f, intype))
_lefttype(xf::Scan{<:Any, T}, _) where {T} = T

outtype(xf::Scan, intype) =
    Union{Base.return_types(xf.f, (_lefttype(xf, intype), intype))...}

function start(rf::R_{Scan}, result)
    if rf.xform.init === nothing
        init = identityof(rf.xform.f, InType(rf))
    else
        init = rf.xform.init :: InType(rf)
    end
    return wrap(rf, init, start(rf.inner, result))
end

function next(rf::R_{Scan}, result, input)
    wrapping(rf, result) do acc, iresult
        acc = rf.xform.f(acc, input)
        return acc, next(rf.inner, iresult, acc)
    end
end

"""
    Iterated(f, init[, T::Type])

Generate a sequence `init`, `f(init)`, `f(f(init))`, `f(f(f(init)))`,
and so on.

$(_shared_notes_unfold)

Use the third argument `T` to specify the output type of `f`.

See also:
[`IterTools.iterated`](https://juliacollections.github.io/IterTools.jl/latest/#IterTools.iterated)

# Examples
```jldoctest
julia> using Transducers

julia> collect(Iterated(x -> 2x, 1), 1:5)
5-element Array{Int64,1}:
  1
  2
  4
  8
 16

julia> collect(TeeZip(Iterated(x -> 2x, 1)), 1:5)
5-element Array{Tuple{Int64,Int64},1}:
 (1, 1)
 (2, 2)
 (3, 4)
 (4, 8)
 (5, 16)
```
"""
struct Iterated{F, T} <: Transducer
    f::F
    init::T

    Iterated(f::F, init, T::Type) where F = new{F, T}(f, init)
end

Iterated(f, init::T) where T = Iterated(f, init, _type_fixedpoint(f, T))

function _type_fixedpoint(f, X, limit = 10)
    for _ in 1:limit
        Y = Union{Base.return_types(f, (X,))...}
        X === Y && return Y
        X = Y
    end
    return Any
end

outtype(xf::Iterated{<:Any, T}, ::Any) where T = T
start(rf::R_{Iterated}, result) = wrap(rf, rf.xform.init, result)
next(rf::R_{Iterated}, result, ::Any) =
    wrapping(rf, result) do istate, iresult
        return rf.xform.f(istate), next(rf.inner, iresult, istate)
    end

"""
    Count(start = 1, step)

Generate a sequence `start`, `start + step`, `start + step + step`,
and so on.

$(_shared_notes_unfold)

`start` defaults 1 and `step` defaults to `oneunit(start)`.

See also:
[`Iterators.countfrom`](https://docs.julialang.org/en/v1/base/iterators/).

# Examples
```jldoctest
julia> using Transducers

julia> collect(TeeZip(Count()), -3:-1)
3-element Array{Tuple{Int64,Int64},1}:
 (-3, 1)
 (-2, 2)
 (-1, 3)
```
"""
struct Count{T} <: Transducer
    start::T
    step::T
end

Count(start = 1) = Count(start, oneunit(start))

outtype(xf::Count{T}, ::Any) where T = T
start(rf::R_{Count}, result) = wrap(rf, rf.xform.start, result)
next(rf::R_{Count}, result, ::Any) =
    wrapping(rf, result) do istate, iresult
        return istate + rf.xform.step, next(rf.inner, iresult, istate)
    end

"""
    TeeZip(xform::Transducer)

Branch input into two "flows", inject one into `xform` and then merge
the output of `xform` with the original input.

To illustrate how it works, consider the following usage

```
xf0 |> TeeZip(xf1) |> xf2
```

where `xf0`, `xf1`, and `xf2` are some transducers.  Schematically,
the output `yn` from `xfn` flows as follows:

```text
xf0      xf1                       xf2
---- y0 ------ y1 ---.-- (y0, y1) ----->
      |              |
       `-------------'
    "Tee"          "Zip"
```

# Examples
```jldoctest
julia> using Transducers

julia> collect(TeeZip(Filter(isodd) |> Map(x -> x + 1)), 1:5)
3-element Array{Tuple{Int64,Int64},1}:
 (1, 2)
 (3, 4)
 (5, 6)
```
"""
struct TeeZip{T} <: Transducer
    xform::T
end
# The idea is to insert an object `Joiner` to the bottom of
# `Reduction` dynamically which calls to the rest of inner reductions
# after the value is zipped (joined).

struct Joiner{F, T}  # Note: not a Transducer
    inner::F  # original inner reduction
    value::T  # original input

    Joiner(inner::F, value::T) where {F,T} = new{F,T}(inner, value)
    Joiner{F,T}(inner) where {F,T} = new{F,T}(inner)
end

next(rf::Joiner, result, input) = next(rf.inner, result, (rf.value, input))
complete(rf::Joiner, result) = complete(rf.inner, result)

outtype(xf::TeeZip, intype) = Tuple{intype, outtype(xf.xform, intype)}

@inline rejoin(rf::R_{TeeZip},
               f = Joiner{typeof(rf.inner), InType(rf)}(rf.inner)) =
    Reduction(rf.xform.xform, f, InType(rf))

start(rf::R_{TeeZip}, result) = start(rejoin(rf), result)
next(rf::R_{TeeZip}, result, input) =
    next(rejoin(rf, Joiner(rf.inner, input)), result, input)
complete(rf::R_{TeeZip}, result) = complete(rejoin(rf), result)

# add joint
Base.adjoint(xf::Transducer) = TeeZip(xf)

"""
    Zip(xforms...)

Zip outputs of transducers `xforms` in a tuple and pass it to the
inner reduction step.

!!! warning
    Head transducers drive tail transducers.  Be careful when using it
    with transducers other than [`Map`](@ref), especially the
    contractive ones like [`PartitionBy`](@ref) and the expansive ones
    like [`MapCat`](@ref).

# Examples
```jldoctest
julia> using Transducers

julia> collect(Zip(Map(identity), Map(x -> 10x), Map(x -> 100x)), 1:3)
3-element Array{Tuple{Int64,Int64,Int64},1}:
 (1, 10, 100)
 (2, 20, 200)
 (3, 30, 300)
```
"""
Zip(xforms...) =
    Map(_zip_init) |> _Zip(xforms...) |> Map(last)
# TODO: add `lower(xf)` mechanism so that constructing Zip does not
# immidiately create a complex composite transducer.

_Zip() = IdentityTransducer()
_Zip(xf1, xforms...) =
    TeeZip(Map(first) |> xf1) |> Map(_zip_between) |> _Zip(xforms...)

_zip_init(y0) = (y0, ())
_zip_between(((y0, ys), yn)) = (y0, (ys..., yn))

"""
    GetIndex(array)
    GetIndex{inbounds}(array)

Transform an integer input `i` to `array[i]`.

# Examples
```jldoctest
julia> using Transducers

julia> collect(GetIndex(1:10), [2, 3, 4])
3-element Array{Int64,1}:
 2
 3
 4

julia> collect(GetIndex{true}(1:10), [2, 3, 4])
3-element Array{Int64,1}:
 2
 3
 4
```
"""
struct GetIndex{inbounds, A} <: Transducer
    array::A
end

GetIndex{inbounds}(array::A) where {inbounds, A} = GetIndex{inbounds, A}(array)
GetIndex(array) = GetIndex{false}(array)

outtype(xf::GetIndex, ::Type{<:Integer}) = eltype(xf.array)
outtype(::GetIndex, T) =
    error("Unexpected non-integer input type for GetIndex:\n", T)

next(rf::R_{GetIndex{true}}, result, input) =
    next(rf.inner, result, @inbounds rf.xform.array[input])
next(rf::R_{GetIndex{false}}, result, input) =
    next(rf.inner, result, rf.xform.array[input])

"""
    SetIndex(array)
    SetIndex{inbounds}(array)

Perform `array[i] = v` for each input pair `(i, v)`.

# Examples
```jldoctest
julia> using Transducers

julia> ys = zeros(3);

julia> mapfoldl(SetIndex(ys), first ∘ tuple, nothing, [(1, 11.1), (3, 33.3)])

julia> ys
3-element Array{Float64,1}:
 11.1
  0.0
 33.3
```
"""
struct SetIndex{inbounds, A} <: Transducer
    array::A
end

SetIndex{inbounds}(array::A) where {inbounds, A} = SetIndex{inbounds, A}(array)
SetIndex(array) = SetIndex{false}(array)

outtype(xf::SetIndex, ::Type{<:Integer}) = eltype(xf.array)
outtype(::SetIndex, T) =
    error("Unexpected non-integer input type for SetIndex:\n", T)

next(rf::R_{SetIndex{true}}, result, input::NTuple{2, Any}) =
    next(rf.inner, result, (@inbounds rf.xform.array[input[1]] = input[2];))
next(rf::R_{SetIndex{false}}, result, input::NTuple{2, Any}) =
    next(rf.inner, result, (rf.xform.array[input[1]] = input[2];))
# Index is `input[1]` due to `TeeZip`'s definition.  Is it better to
# flip, to be compatible with `Base.setindex!`?

"""
    Merge(iterator)

Merge the output from `iterator` to the stream processed by the inner
reduction step.

# Examples
```jldoctest
julia> using Transducers

julia> collect(Merge(Iterators.cycle("hello")), 1:8)
8-element Array{Tuple{Int64,Char},1}:
 (1, 'h')
 (2, 'e')
 (3, 'l')
 (4, 'l')
 (5, 'o')
 (6, 'h')
 (7, 'e')
 (8, 'l')
```
"""
struct Merge{T} <: Transducer
    iterator::T
end

outtype(xf::Merge, intype) = Tuple{intype, ieltype(xf.iterator)}
start(rf::R_{Merge}, result) = wrap(rf, iterate(rf.xform.iterator), result)
next(rf::R_{Merge}, result, input) =
    wrapping(rf, result) do istate, iresult
        istate === nothing && return istate, ensure_reduced(iresult)
        y, s = istate
        iresult2 = next(rf.inner, iresult, (input, y))
        return iterate(rf.xform.iterator, s), iresult2
    end


text:
function Base.map(f, a::StaticArrayLite)
    T = Core.Compiler.return_type(f, Tuple{eltype(a)})
    out = similar(a, T, size(a))
    map!(f, out, a)
    return freeze(out)
end

function Base.map(f, a::StaticArrayLite, b::StaticArrayLite)
    T = Core.Compiler.return_type(f, Tuple{eltype(a), eltype(b)})
    out = similar(a, T, size(a))
    map!(f, out, a, b)
    return freeze(out)
end

function Base.map!(f, out::StaticArrayLite)
    for i in 1:length(out)
        @inbounds out[i] = f()
    end
end

function Base.map!(f, out::StaticArrayLite, a::StaticArrayLite)
    for i in 1:length(out)
        @inbounds out[i] = f(a[i])
    end
end

function Base.map!(f, out::StaticArrayLite, a::StaticArrayLite, b::StaticArrayLite)
    for i in 1:length(out)
        @inbounds out[i] = f(a[i], b[i])
    end
end

function Base.mapreduce(f, op, a::StaticArrayLite; init)
    out = init
    for i in 1:length(a)
        out = op(f(@inbounds(a[i])), out)
    end
    return out
end

function Base.mapreduce(f, op, a::StaticArrayLite, b::StaticArrayLite; init)
    out = init
    for i in 1:length(a)
        out = op(f(@inbounds(a[i]),@inbounds(b[i])), out)
    end
    return out
end


text:
abstract SmoothingKernel
immutable IdentityKernel <: SmoothingKernel
end

immutable TruncatedKernel <: SmoothingKernel
    ξ::Real
    St::Real
    smoother::Function
    κ₁::Real
    κ₂::Real
end

immutable BartlettKernel <: SmoothingKernel
    ξ::Real
    St::Real
    smoother::Function
    κ₁::Real
    κ₂::Real
end

function TruncatedKernel(ξ::Real)
    function smoother(G::Array{Float64, 2})
        T, M = size(G)
        nG   = zeros(T, M)
        for m=1:M
            for t=1:T
			             low = max((t-T), -ξ)
			             high = min(t-1, ξ)
					           for s = low:high
                    @inbounds nG[t, m] += G[t-s, m]
                end
            end
        end
        return(nG/(2.0*ξ+1.0))
    end
    TruncatedKernel(ξ, (2.0*ξ+1.0)/2.0, smoother, 1.0, 1.0)
end

function BartlettKernel(ξ::Real)
    function smoother(G::Array{Float64, 2})
        T, M = size(G)
        nG   = zeros(T, M)
        St   = (2.0*ξ+1.0)/2.0
        for m=1:M
            for t=1:T
			             low = max((t-T), -ξ)
			             high = min(t-1, ξ)
					           for s = low:high
                    κ = 1.0-s/St
                    @inbounds nG[t, m] += κ*G[t-s, m]
                end
            end
        end
        return(nG/(2*ξ+1))
    end
    BartlettKernel(ξ, (2.0*ξ+1.0)/2.0, smoother, 1.0, 2.0/3.0)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct CreateOrUpdateTrustedIdProviderProperties <: SwaggerModel
    idProvider::Any # spec type: Union{ Nothing, String } # spec name: idProvider

    function CreateOrUpdateTrustedIdProviderProperties(;idProvider=nothing)
        o = new()
        validate_property(CreateOrUpdateTrustedIdProviderProperties, Symbol("idProvider"), idProvider)
        setfield!(o, Symbol("idProvider"), idProvider)
        o
    end
end # type CreateOrUpdateTrustedIdProviderProperties

const _property_map_CreateOrUpdateTrustedIdProviderProperties = Dict{Symbol,Symbol}(Symbol("idProvider")=>Symbol("idProvider"))
const _property_types_CreateOrUpdateTrustedIdProviderProperties = Dict{Symbol,String}(Symbol("idProvider")=>"String")
Base.propertynames(::Type{ CreateOrUpdateTrustedIdProviderProperties }) = collect(keys(_property_map_CreateOrUpdateTrustedIdProviderProperties))
Swagger.property_type(::Type{ CreateOrUpdateTrustedIdProviderProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_CreateOrUpdateTrustedIdProviderProperties[name]))}
Swagger.field_name(::Type{ CreateOrUpdateTrustedIdProviderProperties }, property_name::Symbol) =  _property_map_CreateOrUpdateTrustedIdProviderProperties[property_name]

function check_required(o::CreateOrUpdateTrustedIdProviderProperties)
    (getproperty(o, Symbol("idProvider")) === nothing) && (return false)
    true
end

function validate_property(::Type{ CreateOrUpdateTrustedIdProviderProperties }, name::Symbol, val)
end


text:
#---------------------------------------------------------
# # [A dust model](@id dust-model)
#---------------------------------------------------------

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/4_dustmodel.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/4_dustmodel.ipynb)

# In this tutorial we will simulate particles of dust that graviationally settle.
# 2D maps of dust deposition are available from within the AIBECS.
# To use these, we will convert them into sources in the top layer of the model grid.
# Once "born", depending on the settling velocity, these dust particles can spread horizontally until they eventually reach the sea floor and are buried in the sediments.

# The governing equation of the 3D concentration field of dust, denoted $x_\mathsf{dust}$, is:
#
# $$\left[\frac{\partial}{\partial t} + \nabla \cdot (\boldsymbol{u} + \boldsymbol{w} + \mathbf{K}\nabla)\right] x_\mathsf{dust} = s_\mathsf{dust}$$
#
# where $\nabla \cdot \left[ \boldsymbol{u} - \mathbf{K} \nabla \right]$ and $\nabla \cdot \boldsymbol{w}$ represent the ocean circulation transport and the vertical settling of dust particles, respectively.
# ([Tracer transport operators are described in the documentation](@ref tracer-transport-operators).)
# The only source term, $s_\mathsf{dust}$, is the dust deposition.

# In AIBECS, we must recast this equation in the generic form
#
# $$\left[\frac{\partial}{\partial t} + \mathbf{T}(\boldsymbol{p})\right] \boldsymbol{x} = \boldsymbol{G}(\boldsymbol{x},\boldsymbol{p}).$$

# We start by telling Julia we want to use the AIBECS and the OCIM0.1 transport matrix for the ocean circulation.

using AIBECS
grd, T_OCIM = OCIM0.load()

# The transport of dust is the sum of the ocean circulation and a settling-transport:

function T_dust(p)
    @unpack fsedremin = p
    return T_OCIM + transportoperator(grd, z -> w(z,p), fsedremin=fsedremin)
end

# For the settling transport of dust particles, we have used the `transportoperator` function
# for which we need to define the settling velocity `w(z,p)` as a function of depth `z` and of the parameters `p`.

#md # !!! note
#md #     We have imposed a `fsedremin` remineralization at the sediment with its keyword argument.
#md #     (Default is `fsedremin=1.0`.)
#md #     This turns the settling transport into a sink that removes all the dust reaching the sea floor.

#nb # > *Note*:
#nb # > We have imposed a `fsedremin` remineralization at the sediment with its keyword argument.
#nb # > (Default is `fsedremin=1.0`.)
#nb # > This turns the settling transport into a sink that removes all the dust reaching the sea floor.

# Following the assumption that $w(z) = w_0 + w' z$ increases linearly with depth, we write it as

function w(z,p)
    @unpack w₀, w′ = p
    return @. w₀ + w′ * z
end
z = depthvec(grd)

# The parameters `fsedremin`, `w₀`, and `w′` will be defined shortly.

# ##### Dust depositon

# The AIBECS allows you to use climatological dust deposition fields from the Mahowald lab.

s_A_2D = AeolianSources.load()

# First, we must regrid these 2D maps onto the surface layer of `grd`, the 3D grid of the ocean circulation.
# To do that, we interpolate the 2D map of dust deposition onto the ocean grid's lat and lon:

s_dust_2D_monthly = s_A_2D[:Dust] # kg m⁻² s⁻¹
s_dust_2D_annual = permutedims(dropdims(sum(s_dust_2D_monthly, dims=3), dims=3), (2,1)) / 12
s_dust_2D = regrid(s_dust_2D_annual, s_A_2D[:lat], s_A_2D[:lon], grd)

# (This required a bit of work to reshape `s_dust_2D_annual` before using `regrid`.)

# Then we create a blank 3D array of zeros, of which we paint the top layer:

s_dust_3D = zeros(size(grd)...)
s_dust_3D[:,:,1] .= s_dust_2D

# we convert it to a SI volumetric unit (i.e., in g m⁻³)
s_dust_3D = ustrip.(upreferred.(s_dust_3D * u"kg/m^2/s" ./ grd.δz_3D))

# Finally, we vectorize it via

s_dust = vectorize(s_dust_3D, grd) # which is the same as `s_dust_3D[iswet(grd)]`

# We then write the generic $\boldsymbol{G}$ function, which is simply the constant dust source:

G_dust(x, p) = s_dust

# ##### Parameters

# We specify some initial values for the parameters and also include units.

import AIBECS: @units, units
import AIBECS: @initial_value, initial_value
@initial_value @units struct DustModelParameters{U} <: AbstractParameters{U}
    w₀::U        | 0.1 | u"km/yr"
    w′::U        | 0.1 | u"km/yr/km"
    fsedremin::U | 5.0 | u"percent"
end

# Finally, thanks to the initial values we provided, we can instantiate the parameter vector succintly as

p = DustModelParameters()

# We generate the state function `F` and its Jacobian `∇ₓF`,

F, ∇ₓF = state_function_and_Jacobian(T_dust, G_dust)

# generate the steady-state problem,

nb = sum(iswet(grd))
x = ones(nb) # initial guess
prob = SteadyStateProblem(F, ∇ₓF, x, p)

# and solve it

s = solve(prob, CTKAlg()).u

# Let's now run some vizualizations using the Plots.jl recipes.

using Plots

# Let's have a look at a map of the dust per m² integrated over the whole water column:

plotverticalintegral(s * u"g/m^3", grd, color=cgrad(:turbid, rev=true, scale=:exp))

# compared to the original source

plt1 = plotverticalintegral(s_dust * u"g/m^3/s" .|> u"mg/m^3/yr", xlabel="", ylabel="", grd, color=cgrad(:starrynight, scale=:exp), title="Surface flux")

# Huh? This is odd isn't it? What is happening? It seems that dust settles slowly enough to actually explore the ocean basins a little bit.
# This is of course dependent on our choice of model parameters.

# Let's look at what is exported below 500 m and compare:

plt2 = plothorizontalslice(s * u"g/m^3" .* w(z,p) * u"m/s" .|> u"mg/yr/m^2", grd, depth=500u"m", color=cgrad(:starrynight, scale=:exp), xlabel="", ylabel="", title="Flux at 500m")
plot(plt1, plt2, layout=(2,1))

# We can also take a look at the global mean profile (the horizontal average)

plothorizontalmean(s * u"g/m^3" .|> u"mg/m^3", grd)

# As you may see, it also seems that with this model, dust (here with a 1% fraction allowed to stay in the system when it reaches the seafloor)
# accumulates a little bit in the deepest boxes of the model grid.

# Don't hesitate to play around with the parameters and see how things change!
# For example, impose a larger fraction that stays in the bottom or change other parameters via

p = DustModelParameters(w₀=0.1u"km/yr", w′=0.1u"km/yr/km", fsedremin=95.0u"percent")
prob = SteadyStateProblem(F, ∇ₓF, x, p)
s = solve(prob, CTKAlg()).u
plotverticalintegral(s * u"g/m^3", grd, color=cgrad(:turbid, rev=true, scale=:exp))

# And watch how all the dust accumulates in the deepest holes of the ocean basins!


text:
function FieldSpacecraftState(arg0::Field, arg1::SpacecraftState)
    return FieldSpacecraftState((Field, SpacecraftState), arg0, arg1)
end

function FieldSpacecraftState(arg0::FieldAbsolutePVCoordinates)
    return FieldSpacecraftState((FieldAbsolutePVCoordinates,), arg0)
end

function FieldSpacecraftState(arg0::FieldAbsolutePVCoordinates, arg1::FieldAttitude)
    return FieldSpacecraftState((FieldAbsolutePVCoordinates, FieldAttitude), arg0, arg1)
end

function FieldSpacecraftState(arg0::FieldAbsolutePVCoordinates, arg1::FieldAttitude, arg2::Map)
    return FieldSpacecraftState((FieldAbsolutePVCoordinates, FieldAttitude, Map), arg0, arg1, arg2)
end

function FieldSpacecraftState(arg0::FieldAbsolutePVCoordinates, arg1::FieldAttitude, arg2::RealFieldElement)
    return FieldSpacecraftState((FieldAbsolutePVCoordinates, FieldAttitude, RealFieldElement), arg0, arg1, arg2)
end

function FieldSpacecraftState(arg0::FieldAbsolutePVCoordinates, arg1::FieldAttitude, arg2::RealFieldElement, arg3::Map)
    return FieldSpacecraftState((FieldAbsolutePVCoordinates, FieldAttitude, RealFieldElement, Map), arg0, arg1, arg2, arg3)
end

function FieldSpacecraftState(arg0::FieldAbsolutePVCoordinates, arg1::Map)
    return FieldSpacecraftState((FieldAbsolutePVCoordinates, Map), arg0, arg1)
end

function FieldSpacecraftState(arg0::FieldAbsolutePVCoordinates, arg1::RealFieldElement)
    return FieldSpacecraftState((FieldAbsolutePVCoordinates, RealFieldElement), arg0, arg1)
end

function FieldSpacecraftState(arg0::FieldAbsolutePVCoordinates, arg1::RealFieldElement, arg2::Map)
    return FieldSpacecraftState((FieldAbsolutePVCoordinates, RealFieldElement, Map), arg0, arg1, arg2)
end

function FieldSpacecraftState(arg0::FieldOrbit)
    return FieldSpacecraftState((FieldOrbit,), arg0)
end

function FieldSpacecraftState(arg0::FieldOrbit, arg1::FieldAttitude)
    return FieldSpacecraftState((FieldOrbit, FieldAttitude), arg0, arg1)
end

function FieldSpacecraftState(arg0::FieldOrbit, arg1::FieldAttitude, arg2::Map)
    return FieldSpacecraftState((FieldOrbit, FieldAttitude, Map), arg0, arg1, arg2)
end

function FieldSpacecraftState(arg0::FieldOrbit, arg1::FieldAttitude, arg2::RealFieldElement)
    return FieldSpacecraftState((FieldOrbit, FieldAttitude, RealFieldElement), arg0, arg1, arg2)
end

function FieldSpacecraftState(arg0::FieldOrbit, arg1::FieldAttitude, arg2::RealFieldElement, arg3::Map)
    return FieldSpacecraftState((FieldOrbit, FieldAttitude, RealFieldElement, Map), arg0, arg1, arg2, arg3)
end

function FieldSpacecraftState(arg0::FieldOrbit, arg1::Map)
    return FieldSpacecraftState((FieldOrbit, Map), arg0, arg1)
end

function FieldSpacecraftState(arg0::FieldOrbit, arg1::RealFieldElement)
    return FieldSpacecraftState((FieldOrbit, RealFieldElement), arg0, arg1)
end

function FieldSpacecraftState(arg0::FieldOrbit, arg1::RealFieldElement, arg2::Map)
    return FieldSpacecraftState((FieldOrbit, RealFieldElement, Map), arg0, arg1, arg2)
end

function add_additional_state(obj::FieldSpacecraftState, arg0::JString, arg1::Vector{RealFieldElement})
    return jcall(obj, "addAdditionalState", FieldSpacecraftState, (JString, Vector{RealFieldElement}), arg0, arg1)
end

function ensure_compatible_additional_states(obj::FieldSpacecraftState, arg0::FieldSpacecraftState)
    return jcall(obj, "ensureCompatibleAdditionalStates", void, (FieldSpacecraftState,), arg0)
end

function get_a(obj::FieldSpacecraftState)
    return jcall(obj, "getA", RealFieldElement, ())
end

function get_abs_pva(obj::FieldSpacecraftState)
    return jcall(obj, "getAbsPVA", FieldAbsolutePVCoordinates, ())
end

function get_additional_state(obj::FieldSpacecraftState, arg0::JString)
    return jcall(obj, "getAdditionalState", Vector{RealFieldElement}, (JString,), arg0)
end

function get_additional_states(obj::FieldSpacecraftState)
    return jcall(obj, "getAdditionalStates", Map, ())
end

function get_attitude(obj::FieldSpacecraftState)
    return jcall(obj, "getAttitude", FieldAttitude, ())
end

function get_date(obj::FieldSpacecraftState)
    return jcall(obj, "getDate", FieldAbsoluteDate, ())
end

function get_e(obj::FieldSpacecraftState)
    return jcall(obj, "getE", RealFieldElement, ())
end

function get_equinoctial_ex(obj::FieldSpacecraftState)
    return jcall(obj, "getEquinoctialEx", RealFieldElement, ())
end

function get_equinoctial_ey(obj::FieldSpacecraftState)
    return jcall(obj, "getEquinoctialEy", RealFieldElement, ())
end

function get_frame(obj::FieldSpacecraftState)
    return jcall(obj, "getFrame", Frame, ())
end

function get_hx(obj::FieldSpacecraftState)
    return jcall(obj, "getHx", RealFieldElement, ())
end

function get_hy(obj::FieldSpacecraftState)
    return jcall(obj, "getHy", RealFieldElement, ())
end

function get_i(obj::FieldSpacecraftState)
    return jcall(obj, "getI", RealFieldElement, ())
end

function get_keplerian_mean_motion(obj::FieldSpacecraftState)
    return jcall(obj, "getKeplerianMeanMotion", RealFieldElement, ())
end

function get_keplerian_period(obj::FieldSpacecraftState)
    return jcall(obj, "getKeplerianPeriod", RealFieldElement, ())
end

function get_le(obj::FieldSpacecraftState)
    return jcall(obj, "getLE", RealFieldElement, ())
end

function get_lm(obj::FieldSpacecraftState)
    return jcall(obj, "getLM", RealFieldElement, ())
end

function get_lv(obj::FieldSpacecraftState)
    return jcall(obj, "getLv", RealFieldElement, ())
end

function get_mass(obj::FieldSpacecraftState)
    return jcall(obj, "getMass", RealFieldElement, ())
end

function get_mu(obj::FieldSpacecraftState)
    return jcall(obj, "getMu", RealFieldElement, ())
end

function get_orbit(obj::FieldSpacecraftState)
    return jcall(obj, "getOrbit", FieldOrbit, ())
end

function get_pv_coordinates(obj::FieldSpacecraftState)
    return jcall(obj, "getPVCoordinates", TimeStampedFieldPVCoordinates, ())
end

function get_pv_coordinates(obj::FieldSpacecraftState, arg0::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedFieldPVCoordinates, (Frame,), arg0)
end

function has_additional_state(obj::FieldSpacecraftState, arg0::JString)
    return jcall(obj, "hasAdditionalState", jboolean, (JString,), arg0)
end

function interpolate(obj::FieldSpacecraftState, arg0::FieldAbsoluteDate, arg1::Stream)
    return jcall(obj, "interpolate", FieldSpacecraftState, (FieldAbsoluteDate, Stream), arg0, arg1)
end

function interpolate(obj::FieldTimeInterpolable, arg0::FieldAbsoluteDate, arg1::Collection)
    return jcall(obj, "interpolate", FieldTimeInterpolable, (FieldAbsoluteDate, Collection), arg0, arg1)
end

function is_orbit_defined(obj::FieldSpacecraftState)
    return jcall(obj, "isOrbitDefined", jboolean, ())
end

function shifted_by(obj::FieldSpacecraftState, arg0::RealFieldElement)
    return jcall(obj, "shiftedBy", FieldSpacecraftState, (RealFieldElement,), arg0)
end

function shifted_by(obj::FieldSpacecraftState, arg0::jdouble)
    return jcall(obj, "shiftedBy", FieldSpacecraftState, (jdouble,), arg0)
end

function to_spacecraft_state(obj::FieldSpacecraftState)
    return jcall(obj, "toSpacecraftState", SpacecraftState, ())
end

function to_string(obj::FieldSpacecraftState)
    return jcall(obj, "toString", JString, ())
end

function to_transform(obj::FieldSpacecraftState)
    return jcall(obj, "toTransform", FieldTransform, ())
end