text:
# --- Transducible processes

"""
    reducingfunction(xf, step; simd)
    xf'(step; simd)

Apply transducer `xf` to the reducing function `step` to create a new
reducing function.

!!! compat "Transducers.jl 0.3"

    New in version 0.3.

!!! warning

    Be careful using `reducingfunction` with stateful transducers like
    [`Scan`](@ref) with mutable `init` (e.g., `Scan(push!, [])`).  See
    more in Examples below.

# Arguments
- `xf::Transducer`: A transducer.
- `step`: A callable which accepts 1 and 2 arguments.  If it only
  accepts 2 arguments, wrap it with [`Completing`](@ref) to "add"
  1-argument form (i.e., [`complete`](@ref) protocol).

# Keyword Arguments
- `simd`: `false`, `true`, or `:ivdep`. See [`maybe_usesimd`](@ref).

# Examples
```jldoctest
julia> using Transducers

julia> rf = reducingfunction(Map(x -> x + 1), *);

julia> rf(10, 2) === 10 * (2 + 1)
true
```

## Warning: Be careful when using `reducingfunction` with stateful transducers

Stateful `Transducer`s themselves in Transducers.jl are not inherently
broken with `reducingfunction`.  However, it can produce incorrect
results when combined with mutable states:

```jldoctest reducingfunction; setup = :(using Transducers)
julia> scan_state = [];

julia> rf_bad = opcompose(Scan(push!, scan_state), Cat())'(string);

julia> transduce(rf_bad, "", 1:3)
"112123"
```

The first run works.  However, observe that the vector `scan_state` is
not empty anymore:

```jldoctest reducingfunction
julia> scan_state
3-element Vector{Any}:
 1
 2
 3
```

Thus, the second run produces an incorrect result:

```jldoctest reducingfunction
julia> transduce(rf_bad, "", 1:3)
"123112312123123"
```

One way to solve this issue is to use [`CopyInit`](@ref) or [`OnInit`](@ref).

```jldoctest reducingfunction
julia> scan_state = CopyInit([])
CopyInit(Any[])

julia> rf_good = opcompose(Scan(push!, scan_state), Cat())'(string);

julia> transduce(rf_good, "", 1:3)
"112123"

julia> scan_state
CopyInit(Any[])

julia> transduce(rf_good, "", 1:3)
"112123"
```
"""
@inline reducingfunction(xf::Transducer, step; simd::SIMDFlag = Val(false)) =
    maybe_usesimd(Reduction(xf, step), simd)

# Use `_asmonoid` automatically only when `init` is not specified:
@inline _reducingfunction(
    xf::XF,
    step::Step;
    init = DefaultInit,
    simd::SIMDFlag = Val(false),
    _...,
) where {XF,Step} =
    maybe_usesimd(Reduction(xf, init === DefaultInit ? _asmonoid(step) : step), simd)

"""
    __foldl__(rf, init, reducible::T)

Left fold a `reducible` with reducing function `rf` and initial value
`init`.  This is primary an API for overloading when the reducible
"container" or "context" (e.g., I/O stream) of type `T` can provide a
better reduction mechanism than the default iterator-based one.

For a simple iterable type `MyType`, a valid implementation is:

```julia
function __foldl__(rf, val, itr::MyType)
    for x in itr
        val = @next(rf, val, x)
    end
    return complete(rf, val)
end
```

although in this case default ` __foldl__` can handle `MyType` and
thus there is no need for defining it.  In general, defining
`__foldl__` is useful only when there is a better way to go over items
in `reducible` than `Base.iterate`.

See also: [`@next`](@ref).
"""
__foldl__

# ** Why calling `complete` inside `__foldl__`? **
#
# Initially I was hoping to stop calling `complete` in `__foldl__` so
# that `skipcomplete` wouldn't be necessary.  However, it was
# impossible to make it inferred with some mildly complex cases.  This
# is probably expected since `__foldl__` tends to use something like
# tail-call function-barriers to maximize the chance of union
# splitting.  So, inside `__foldl__`, `complete` has a better chance
# to be inferred to have a concrete input type.  Furthermore, since
# private types are unwrapped inside `complete`, a simpler type would
# be returned from `__foldl__`.  This is probably why
# `complete`-inside-`__foldl__` approach is more compiler-friendly.
# On the other hand, if `complete` is called outside `__foldl__`, the
# compiler has to merge various types of nested objects returned from
# multiple locations.  This easily leads to non-concrete types in the
# inference.

const FOLDL_RECURSION_LIMIT = Val(1)
# const FOLDL_RECURSION_LIMIT = nothing
_dec(::Nothing) = nothing
_dec(::Val{n}) where n = Val(n - 1)

function __foldl__(rf::RF, init::T, coll) where {RF,T}
    ret = iterate(coll)
    ret === nothing && return complete(rf, init)
    x, state = ret
    val = @next(rf, init, x)
    @manual_union_split val isa T begin
        return _foldl_iter(rf, val, coll, state, FOLDL_RECURSION_LIMIT)
    end
end

@inline function _foldl_iter(rf::RF, val::T, iter, state, counter) where {RF,T}
    while true
        ret = iterate(iter, state)
        ret === nothing && break
        x, state = ret
        y = @next(rf, val, x)
        counter === Val(0) || y isa T ||
            return _foldl_iter(rf, y, iter, state, _dec(counter))
        val = y
    end
    return complete(rf, val)
end

@inline __foldl__(rf::RF, init, coll::Tuple) where {RF} =
    complete(rf, @return_if_reduced foldlargs(rf, init, coll...))

@inline function __foldl__(rf::RF, init, arr0::Union{AbstractArray,Broadcasted}) where {RF}
    arr = Broadcast.instantiate(arr0)
    isempty(arr) && return complete(rf, init)
    return _foldl_array(rf, init, arr, _IndexStyle(arr))
end

@inline function _foldl_array(rf::RF, init::T, arr, ::IndexLinear) where {RF,T}
    i = _firstindex(arr)
    acc = @next(rf, init, @inbounds arr[i])
    @manual_union_split acc isa T begin
        if is_prelude(acc)
            return _foldl_linear_rec(rf, acc, arr, i + 1, FOLDL_RECURSION_LIMIT)
        else
            return _foldl_linear_bulk(rf, acc, arr, i + 1)
        end
    end
end

@inline function _foldl_linear_bulk(rf::RF, acc, arr, i0) where {RF}
    @simd_if rf for i in i0:_lastindex(arr)
        acc = @next(rf, acc, @inbounds arr[i])
    end
    return restack(complete(rf, acc))
end

@inline function _foldl_linear_rec(rf::RF, acc::T, arr, i0, counter) where {RF,T}
    for i in i0:_lastindex(arr)
        y = @next(rf, acc, @inbounds arr[i])
        if counter !== Val(0)
            if y isa T
            elseif is_prelude(y)
                return _foldl_linear_rec(rf, y, arr, i + 1, _dec(counter))
            else
                # Otherwise, maybe it could be something like `Union{Float64,Missing}`
                # where Julia's native loop is fast:
                return _foldl_linear_bulk(rf, y, arr, i + 1)
            end
        end
        acc = y
    end
    return complete(rf, acc)
end

@inline function _foldl_array(rf0::RF, init, arr, ::IndexStyle) where {RF,T}
    @inline getvalue(I) = @inbounds arr[I]
    rf = Map(getvalue)'(rf0)
    return __foldl__(rf, init, _CartesianIndices(arr))
end

@inline _getvalues(i) = ()
@inline _getvalues(i, a, rest...) = ((@inbounds a[i]), _getvalues(i, rest...)...)

# TODO: merge this with array implementation
@static if VERSION >= v"1.1-"
    @inline function __foldl__(
            rf, init,
            zs::Iterators.Zip{<:Tuple{Vararg{AbstractArray}}})
        isempty(zs) && return complete(rf, init)
        idxs = eachindex(zs.is...)
        val = @next(rf, init, _getvalues(firstindex(idxs), zs.is...))
        @simd_if rf for i in firstindex(idxs) + 1:lastindex(idxs)
            val = @next(rf, val, _getvalues(i, zs.is...))
        end
        return complete(rf, val)
    end
end

## Convert zip-of-products (ZoP) to product-of-zips (PoZ).

# Arguments for `product` and alike:
unproduct(x::Iterators.ProductIterator) = x.iterators
unproduct(x::CartesianIndices) = x.indices
unproduct(x::AbstractArray) = axes(x)

# After ZoP-to-PoZ transformation, we need to re-shape the value in
# the form that would be fed to the reducing function if it were
# folding `ZoP`:
@inline function _make_zop_getvalues(iterators)
    # Avoid including `iterators` themselves in the closure if
    # possible:
    iterinfo = map(iterators) do itr
        if itr isa CartesianIndices
            Val(CartesianIndices)
        elseif itr isa AbstractArray
            itr
        else
            nothing
        end
    end
    return function (xs)
        map(_unzip((iterinfo, _unzip(xs)))) do (it, x)
            if it === Val(CartesianIndices)
                return CartesianIndex(x)
            elseif it isa AbstractArray
                return @inbounds it[x...]
            else
                return x
            end
        end
    end
end

@static if VERSION >= v"1.1-"
    @inline __foldl__(
        rf,
        init,
        zs::Iterators.Zip{<:Tuple{
            Vararg{Union{AbstractArray,Iterators.ProductIterator}},
        }},
    ) = __foldl__(
        Reduction(Map(_make_zop_getvalues(zs.is)), rf),
        init,
        Iterators.product(map(Base.splat(zip), _unzip(map(unproduct, zs.is)))...),
    )
end

@inline function __foldl__(rf0, init, cartesian::CartesianIndices)
    rf = Map(CartesianIndex)'(rf0)
    val = _foldl_product(rf, init, (), cartesian.indices...)
    val isa Reduced && return val
    return complete(rf, val)
end

@inline function __foldl__(
        rf, init,
        prod::Iterators.ProductIterator{<:Tuple{Any,Any,Vararg{Any}}})
    val = _foldl_product(rf, init, (), prod.iterators...)
    val isa Reduced && return val
    return complete(rf, val)
end

@noinline _foldl_product(rf, val, ::Any) = error("Unreachable")
@inline _foldl_product(rf, val, ::Tuple) = next(rf, val, ())

@inline function _foldl_product(rf, val, outer, iterators...)
    inner_iterators, outer_iterator = _poptail(iterators)
    for input in outer_iterator
        val_ = _foldl_product(rf, val, (input, outer...), inner_iterators...)
        val_ isa Reduced && return val_
        val = val_
    end
    return val
end

@inline function _foldl_product(rf, val, outer, iterator)
    # TODO: Handle the case inner iterators are tuples.  In such case,
    # inner-most non-tuple iterators should use @simd_if.
    @simd_if rf for input in iterator
        val = @next(rf, val, (input, outer...))
    end
    return val
end

function __simple_foldl__(rf, val, itr)
    for x in itr
        val = @next(rf, val, x)
    end
    return complete(rf, val)
end

"""
    simple_transduce(xform, step, init, coll)

Simplified version of [`transduce`](@ref).  For simple transducers Julia
may be able to emit a good code.  This function exists only for
performance tuning.
"""
function simple_transduce(xform, f, init, coll)
    rf = Reduction(xform, f)
    return __simple_foldl__(rf, start(rf, init), coll)
end

"""
    foldl_nocomplete(rf, init, coll)

Call [`__foldl__`](@ref) without calling [`complete`](@ref).
"""
@inline foldl_nocomplete(rf::RF, init, coll) where {RF} =
    __foldl__(skipcomplete(rf), init, coll)

"""
    foldxl(step, xf::Transducer, reducible; init, simd) :: T
    foldxl(step, reducible; init, simd) :: T
    foldl(step, xf::Transducer, reducible; init, simd) :: T
    foldl(step, ed::Eduction; init, simd) :: T
    transduce(xf, step, init, reducible, [executor]; simd) :: Union{T, Reduced{T}}

e**X**tended **l**eft fold.
Compose transducer `xf` with reducing step function `step` and reduce
`itr` using it.

!!! note
    `transduce` differs from `foldxl` as `Reduced{T}` is returned if
    the transducer `xf` or `step` aborts the reduction and `step` is
    _not_ automatically wrapped by [`Completing`](@ref).

This API is modeled after $(_cljref("transduce")).

For parallel versions, see [`foldxt`](@ref) and [`foldxd`](@ref).

See also: [Empty result handling](@ref).

# Arguments
- `xf::Transducer`: A transducer.
- `step`: A callable which accepts 1 and 2 arguments.  If it only
  accepts 2 arguments, wrap it with [`Completing`](@ref) to "add"
  1-argument form (i.e., [`complete`](@ref) protocol).
- `reducible`: A reducible object (array, dictionary, any iterator, etc.).
- `executor`: Specify an executor. See [`SequentialEx`](@ref).
- `init`: An initial value fed to the first argument to reducing step
  function `step`.  This argument can be omitted for well know binary
  operations like `+` or `*`.  Supported binary operations are listed
  in InitialValues.jl documentation.  When `Init` (not the result of
  `Init`, such as `Init(*)`) is given, it is automatically "instantiated"
  as `Init(step)` (where `step` is appropriately unwrapped if `step` is
  a `Completing`).  See [Empty result handling](@ref) in the manual
  for more information.
- `simd`: If `true` or `:ivdep`, enable SIMD using `Base.@simd`.  If
  `:ivdep`, use `@simd ivdep for ... end` variant.  Read Julia manual
  of `Base.@simd` to understand when it is appropriate to use this
  option.  For example, `simd = :ivdep` _must not_ be used with stateful
  transducer like [`Scan`](@ref).  This option has no effect if
  `false` (default).

# Examples
```jldoctest
julia> using Transducers

julia> foldl(Filter(isodd), 1:4, init=0.0) do state, input
           @show state, input
           state + input
       end
(state, input) = (0.0, 1)
(state, input) = (1.0, 3)
4.0

julia> foldxl(+, 1:5 |> Filter(isodd))
9

julia> 1:5 |> Filter(isodd) |> foldxl(+)
9
```

Since [`TeeRF`](@ref) requires the extended fold protocol,
`foldl(TeeRF(min, max), [5, 2, 6, 8, 3])` does not work while it works
with `foldxl`:

```jldoctest; setup = :(using Transducers)
julia> foldxl(TeeRF(min, max), [5, 2, 6, 8, 3])
(2, 8)
```

The unary method of `foldlx` is useful when combined with `|>`:

```jldoctest; setup = :(using Transducers)
julia> (1:5, 4:-1:1) |> Cat() |> Filter(isodd) |> Enumerate() |> foldxl() do a, b
           a[2] < b[2] ? b : a
       end
(3, 5)
```
"""
foldxl
foldxl(rf; kw...) = itr -> foldxl(rf, itr; kw...)
foldxl(rf, itr; kw...) = foldl(rf, extract_transducer(itr)...; kw...)

"""
    foldl(step, xf::Transducer, reducible; init, simd) :: T
    foldl(step, ed::Eduction; init, simd) :: T

See [`foldxl`](@ref).
"""
foldl

"""
    transduce(xf, step, init, reducible) :: Union{T, Reduced{T}}

See [`foldxl`](@ref).
"""
transduce

const _MAPFOLDL_DEPWARN = (
    "`mapfoldl(::Transducer, rf, itr)` is deprecated. " *
    " Use `foldl(rf, ::Transducer, itr)` if you do not need to call single-argument" *
    " `rf` on `complete`." *
    " Use `foldl(whencomplete(rf, rf), ::Transducer, itr)` to call the" *
    " single-argument method of `rf` on complete."
)

"""
    mapfoldl(xf::Transducer, step, reducible; init, simd)

!!! warning

    $_MAPFOLDL_DEPWARN

Like [`foldl`](@ref) but `step` is _not_ automatically wrapped by
[`Completing`](@ref).

# Examples
```julia
julia> using Transducers

julia> function step_demo(state, input)
           @show state, input
           state + input
       end;

julia> function step_demo(state)
           println("Finishing with state = ", state)
           state
       end;

julia> mapfoldl(Filter(isodd), step_demo, 1:4, init=0.0)
(state, input) = (0.0, 1)
(state, input) = (1.0, 3)
Finishing with state = 4.0
4.0
```
"""
mapfoldl

function transduce(xform::Transducer, f::F, init, coll; kwargs...) where {F}
    rf = _reducingfunction(xform, f; init = init, kwargs...)
    return transduce(rf, init, coll; kwargs...)
end

_unreduced__foldl__(rf, step, coll) = unreduced(__foldl__(rf, step, coll))

# TODO: should it be an internal?
@inline function transduce(
    rf1::RF,
    init,
    coll;
    simd::SIMDFlag = Val(false),
) where {RF<:AbstractReduction}
    # Inlining `transduce` and `__foldl__` were essential for the
    # `restack` below to work.
    rf0, foldable = retransform(rf1, asfoldable(coll))
    rf = maybe_usesimd(rf0, simd)
    state = start(rf, init)
    result = __foldl__(rf, state, foldable)
    if unreduced(result) isa DefaultInitOf
        throw(EmptyResultError(rf0))
        # Should I check if `init` is a `DefaultInit`?
    end
    # At this point, `return result` is the semantically correct thing
    # to do.  What follows are some convoluted instructions to
    # convince the compiler that this function is type-stable (in some
    # cases).  Note that return type would be inference-dependent only
    # if `init` is a `OptInit` type.  In the default case where `init
    # isa DefaultInitOf`, the real code pass is the `throw` above.

    # Unpacking as `ur_result` and re-packing it later somehow helps
    # the compiler to correctly eliminate a possibility in a `Union`.
    ur_result = unreduced(result)
    if ur_result isa InferableInit
        # Using `rf0` instead of `rf` helps the compiler.  Note that
        # this means that we are relying on that enabling SIMD does
        # not change the return type.
        realtype = _nonidtype(Core.Compiler.return_type(
            _unreduced__foldl__,
            typeof((rf0, state, foldable)),
        ))
        if realtype isa Type
            realvalue = convert(realtype, ur_result)
            if result isa Reduced
                return Reduced(realvalue)
            else
                return realvalue
            end
        end
    end
    if result isa Reduced
        return Reduced(ur_result)
    else
        return ur_result
    end
end

Base.mapfoldl(f::F, step::OP, itr::Foldable; kw...) where {F, OP} =
    foldxl(step, Map(f), itr; kw...)

struct Eduction{F, C} <: Foldable
    rf::F
    coll::C
end

Eduction(xform::Transducer, coll) =
    Eduction(Reduction(xform, Completing(push!!)), coll)

Eduction(xform::Transducer, ed::Eduction) =
    Eduction(opcompose(Transducer(ed), xform), ed.coll)

Transducer(ed::Eduction) = Transducer(ed.rf)

transduce(xform::Transducer, f, init, ed::Eduction) =
    transduce(opcompose(Transducer(ed), xform), f, init, ed.coll)

Base.IteratorSize(::Type{Eduction{F,C}}) where{F,C} =
    outputsize(F) isa SizeStable ? Base.IteratorSize(C) : Base.SizeUnknown()

function Base.length(ed::Eduction)
    @argcheck Base.IteratorSize(ed) isa Union{Base.HasLength,Base.HasShape}
    return length(ed.coll)
end

function Base.size(ed::Eduction)
    @argcheck Base.IteratorSize(ed) isa Base.HasShape
    return size(ed.coll)
end

function Base.iterate(ts::Eduction, state = nothing)
    if state === nothing
        cret = iterate(ts.coll)
        cret === nothing && return nothing
        input, cstate = cret

        result = start(ts.rf, Union{}[])

        cdone = false
        rdone = false
        @goto step
        # Assuming the initial `result` can never be `Reduced`.
    else
        cstate, cdone, result, rdone = state
    end
    if !cdone && !rdone
        while isempty(unwrap_all(result))
            cret = iterate(ts.coll, cstate)
            if cret === nothing
                cdone = true
                result = complete(ts.rf, unreduced(result))
                # Stateful transducers may have flushed something.
                # Let's not `return nothing` yet.
                break
            end
            input, cstate = cret

            @label step
            result = next(ts.rf, result, input)
            if isreduced(result)
                rdone = true
                result = unreduced(result)
                break
            end
        end
    end
    buffer = unwrap_all(result)
    isempty(buffer) && return nothing
    y = popfirst!(buffer)
    next_state = (cstate, cdone, result, rdone)
    return (y, next_state)
end

"""
    eduction(xf::Transducer, coll)
    xf(coll)
    coll |> xf

Create a iterable and reducible object.

* [Iterable](https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-iteration-1).
* Reducible; i.e., it can be handled by [`transduce`](@ref) efficiently.

This API is modeled after $(_cljref("eduction")).

!!! note

    Even though `eduction` returns an iterable, it is highly
    recommended to use the `foldl`-based method provided by
    Transducers.jl when the performance is important.

# Examples
```jldoctest
julia> using Transducers

julia> for x in 1:1000 |> Filter(isodd) |> Take(3)  # slow
           @show x
       end
x = 1
x = 3
x = 5

julia> foreach(1:1000 |> Filter(isodd) |> Take(3)) do x  # better
           @show x
       end;
x = 1
x = 3
x = 5
```
"""
eduction(xform, coll) = Eduction(xform, coll)
# Note on `simd` keyword argument: `eduction` ATM does not support
# `simd` argument which could be done in principle.  However, how
# `foldl` and `foreach` with `Eduction` treat `simd` argument must be
# tweaked if that happens.
#
# Note on API:
# Exporting `Eduction` could also work.  But `Base` has, e.g.,
# `skipmissing` so maybe this is better for more uniform API.

"""
    setinput(ed::Eduction, coll)

Set input collection of eduction `ed` to `coll`.

!!! compat "Transducers.jl 0.3"

    Previously, `setinput` combined with `eduction` was a recommended
    way to use transducers in a type stable manner.  As of v0.3, all
    the `foldl`-like functions and `eduction` are type stable for many
    cases.  This workaround is no more necessary.

# Examples
```jldoctest
julia> using Transducers

julia> ed = eduction(Map(x -> 2x), Float64[]);

julia> xs = ones(2, 3);

julia> foldl(+, setinput(ed, xs))
12.0
```
"""
setinput(ed::Eduction, coll) =
    _setinput(eltype(ed.coll), eltype(coll), ed, coll)

_setinput(::Type{T}, ::Type{T}, ed, coll) where T = @set ed.coll = coll
_setinput(::Type, ::Type, ed, coll) = eduction(Transducer(ed), coll)

"""
    append!(xf::Transducer, dest, src) -> dest

This API is modeled after $(_cljref("into")).

!!! warning
    The performance of `append!(dest, src::Eduction)` is poor.
    Use `append!!` instead if two-argument form is preferred.

# Examples
```jldoctest
julia> using Transducers

julia> append!(Drop(2), [-1, -2], 1:5)
5-element Vector{Int64}:
 -1
 -2
  3
  4
  5
```
"""
Base.append!(xf::Transducer, to, from) =
    unreduced(transduce(xf, Completing(push!), to, from))

"""
    BangBang.append!!(xf::Transducer, dest, src) -> dest′
    BangBang.append!!(dest, src::Eduction) -> dest′

Mutate-or-widen version of [`append!`](@ref).

!!! compat "Transducers.jl 0.4.4"

    New in version 0.4.4.

!!! compat "Transducers.jl 0.4.37"

    Performance optimization for `append!!(dest, src::Eduction)`
    requires version 0.4.37.

# Examples
```jldoctest
julia> using Transducers, BangBang

julia> append!!(opcompose(Drop(2), Map(x -> x + 0.0)), [-1, -2], 1:5)
5-element Vector{Float64}:
 -1.0
 -2.0
  3.0
  4.0
  5.0
```
"""
BangBang.append!!
BangBang.append!!(xf::Transducer, to, from) = unreduced(transduce(
    Map(SingletonVector) ∘ xf,
    Completing(append!!),
    to,
    from,
))

function BangBang.__appendto!!__(to, foldable::Foldable)
    xf, from = extract_transducer(foldable)
    return append!!(xf, to, from)
end

"""
    collect(xf::Transducer, itr) :: Vector
    collect(ed::Eduction) :: Vector

Process an iterable `itr` using a transducer `xf` and collect the result
into a `Vector`.

For parallel versions, see [`tcollect`](@ref) and [`dcollect`](@ref).

!!! compat "Transducers.jl 0.4.8"

    `collect` now accepts eductions.

# Examples
```jldoctest
julia> using Transducers

julia> collect(Interpose(missing), 1:3)
5-element Vector{Union{Missing, Int64}}:
 1
  missing
 2
  missing
 3
```
"""
function Base.collect(xf::Transducer, coll)
    result = finish!(unreduced(transduce(
        Map(SingletonVector) ∘ xf,
        wheninit(collector, append!!),
        collector(),
        coll,
    )))
    if result isa Vector{Union{}}
        et = @default_finaltype(xf, coll)
        return et[]
    end
    return result
end
# Base.collect(xf, coll) = append!([], xf, coll)

Base.collect(foldable::Foldable) = collect(extract_transducer(foldable)...)

"""
    copy(xf::Transducer, T, foldable) :: Union{T, Empty{T}}
    copy(xf::Transducer, foldable::T) :: Union{T, Empty{T}}
    copy([T,] eduction::Eduction) :: Union{T, Empty{T}}

Process `foldable` with a transducer `xf` and then create a container of type `T`
filled with the result.  Return
[`BangBang.Empty{T}`](https://juliafolds.github.io/BangBang.jl/dev/#BangBang.NoBang.Empty)
if the transducer does not produce anything.  (This is because there is no
consistent interface to create an empty container given its type and not all
containers support creating an empty container.)

For parallel versions, see [`tcopy`](@ref) and [`dcopy`](@ref).

!!! compat "Transducers.jl 0.4.4"

    New in version 0.4.4.

!!! compat "Transducers.jl 0.4.8"

    `copy` now accepts eductions.

# Examples
```jldoctest
julia> using Transducers
       using BangBang: Empty

julia> copy(Map(x -> x => x^2), Dict, 2:2)
Dict{Int64, Int64} with 1 entry:
  2 => 4

julia> @assert copy(Filter(_ -> false), Set, 1:1) === Empty(Set)

julia> using TypedTables

julia> @assert copy(Map(x -> (a=x, b=x^2)), Table, 1:1) == Table(a=[1], b=[1])

julia> using StructArrays

julia> @assert copy(Map(x -> (a=x, b=x^2)), StructVector, 1:1) == StructVector(a=[1], b=[1])

julia> using DataFrames

julia> @assert copy(
           Map(x -> (A = x.a + 1, B = x.b + 1)),
           DataFrame(a = [1], b = [2]),
       ) == DataFrame(A = [2], B = [3])
```
"""
Base.copy(xf::Transducer, ::Type{T}, foldable) where {T} = append!!(xf, Empty(T), foldable)
Base.copy(xf::Transducer, foldable) = copy(xf, _materializer(foldable), foldable)

function Base.copy(::Type{T}, ed::Foldable) where {T}
    xf, foldable = extract_transducer(ed)
    return copy(xf, T, foldable)
end

function Base.copy(ed::Foldable)
    xf, foldable = extract_transducer(ed)
    return copy(xf, foldable)
end

Base.Set(foldable::Foldable) = copy(Set, foldable)
Base.Dict(foldable::Foldable) = copy(Dict, foldable)

"""
    map!(xf::Transducer, dest, src; simd)

Feed `src` to transducer `xf`, storing the result in `dest`.
Collections `dest` and `src` must have the same shape.  Transducer
`xf` may contain filtering transducers.  If some entries `src` are
skipped, the corresponding entries in `dest` will be unchanged.
Transducer `xf` must not contain any expansive transducers such as
[`MapCat`](@ref).

See also [`copy!`](@ref).

# Examples
```jldoctest
julia> using Transducers

julia> xs = collect(1:5)
       ys = zero(xs)
       map!(Filter(isodd), ys, xs)
5-element Vector{Int64}:
 1
 0
 3
 0
 5

julia> ans === ys
true
```
"""
function Base.map!(xf::Transducer, dest::AbstractArray, src::AbstractArray;
                   simd::SIMDFlag = Val(false))
    _map!(_prepare_map(xf, dest, src, simd)...)
    return dest
end

_map!(rf, coll, dest) = transduce(restack(rf), nothing, coll)

# The idea behind `restack` (previously called `darkritual`):
# Deep-copy `AbstractReduction` so that compiler can treat the all
# reducing function tree nodes as local variables (???).  Aslo, it
# tells compiler that `dest` is a local variable so that it won't
# fetch `dest` via `getproperty` in each iteration.  (This is too much
# magic...  My reasoning of how it works could be completely wrong.
# But at least it should not change the semantics of the function.)
# Probably related: https://github.com/JuliaLang/julia/pull/18632

function _prepare_map(xf, dest, src, simd)
    isexpansive(xf) && error("map! only supports non-expanding transducer")
    # TODO: support Dict
    indices = eachindex(dest, src)

    rf = reducingfunction(
        opcompose(ZipSource(opcompose(GetIndex{true}(src), xf)), SetIndex{true}(dest)),
        (::Vararg) -> nothing,
        simd = simd)

    return rf, indices, dest
end

"""
    copy!(xf::Transducer, dest, src)

Feed `src` to transducer `xf`, storing the result in `dest`.
Collections `dest` and `src` may have the same shape.  Source `src`
must be iterable.  Destination `dest` must implement `empty!` and
`push!`.

See also [`map!`](@ref).

# Examples
```jldoctest
julia> using Transducers

julia> copy!(opcompose(PartitionBy(x -> x ÷ 3), Map(sum)), Int[], 1:10)
4-element Vector{Int64}:
  3
 12
 21
 19
```
"""
Base.copy!(xf::Transducer, dest, src) = append!(xf, empty!(dest), src)

foldxl(step::RF, xform, itr; init = DefaultInit, kw...) where {RF} =
    unreduced(transduce(xform, Completing(step), init, itr; kw...))

Base.foldl(step::RF, xform::Transducer, itr; kw...) where {RF} =
    foldxl(step, xform, itr; kw...)

@inline Base.foldl(step::RF, foldable::Foldable; kw...) where {RF} =
    foldxl(step, foldable; kw...)

"""
    foreach(eff, xf::Transducer, reducible; simd)
    foreach(eff, ed::Eduction; simd)

Feed the results of `xf` processing items in `reducible` into a unary
function `eff`.  This is useful when the primary computation at the
bottom is the side-effect.  It is also equivalent to `foreach(eff,
eduction(xf, coll))`.  Note that

```julia
foreach(eduction(xf, coll)) do x
    ...
end
```

can be more efficient than

```julia
for x in eduction(xf, coll)
    ...
end
```

as the former does not have to translate the transducer protocol to
the iterator protocol.

`foreach` supports all constructs in the native `for` loop as well as
the enhancements [^julia_issue_22891] to `break` with a value (`break
D(x)` below) and append the `else` clause (`E(x)` below).

[^julia_issue_22891]: See also: [break with value + loop else clauses
    (JuliaLang/julia#22891)](https://github.com/JuliaLang/julia/issues/22891)

This native `for` loop

```julia
ans = for x in xs
    A(x)
    B(x) && break
    C(x) && break D(x)
else
    E(x)
end
```

can be written as

```julia
ans = foreach(Map(identity), xs) do x
    A(x)
    B(x) && return reduced()
    C(x) && return reduced(D(x))
    x  # required for passing `x` to `E(x)` below
end |> ifunreduced() do x
    E(x)
end
```

See: [`mapfoldl`](@ref), [`reduced`](@ref), [`ifunreduced`](@ref).

!!! compat "Transducers.jl 0.3"

    `foreach` is changed to return what the `do` block (`eff`
    function) returns as-is in version 0.3.  This was required for
    supporting "for-else" (`|> ifunreduced`).  Previously, it only
    supported break-with-value and always applied `unreduced` before
    it returns.

# Examples
```jldoctest
julia> using Transducers

julia> foreach(eduction(Filter(isodd), 1:4)) do input
           @show input
       end
input = 1
input = 3
3

julia> foreach(Filter(!ismissing), [1, missing, 2, 3]) do input
           @show input
           if iseven(input)
               return reduced()
           end
       end
input = 1
input = 2
Reduced(nothing)
```

It is often useful to append [`|> unreduced`](@ref unreduced) to
unwrap `Reduced` in the final result (note that `|>` here is the
standard function application, not the transducer composition).

```jldoctest; setup = :(using Transducers)
julia> foreach(Filter(!ismissing), [1, missing, 2, 3]) do input
           reduced("got \$input")
       end |> unreduced
"got 1"
```

Combination of break-with-value and for-else is useful for triggering
action after (e.g.) some kind of membership testing failed:

```jldoctest; setup = :(using Transducers)
julia> has2(xs) = foreach(Filter(!ismissing), xs) do input
           input == 2 && reduced(true)
       end |> ifunreduced() do input
           @show input
           false
       end;

julia> has2([1, missing, 2, 3])
true

julia> has2([1, missing])
input = false
false
```

However, note the output `input = false` in the last example.  This is
because how `&&` works in Julia

```jldoctest
julia> false && "otherwise"
false
```

Thus, pure membership testing functions like `has2` above can be
written in a more concise manner:

```jldoctest; setup = :(using Transducers)
julia> simpler_has2(xs) = foreach(Filter(!ismissing), xs) do input
           input == 2 && reduced(true)
       end |> unreduced;

julia> simpler_has2([1, missing, 2, 3])
true

julia> simpler_has2([1, missing])
false
```
"""
Base.foreach(eff::F, xform::Transducer, coll; kwargs...) where {F} =
    transduce(xform, SideEffect(eff), nothing, coll; kwargs...)
function Base.foreach(eff::F, reducible::Reducible; kwargs...) where {F}
    xf, coll = extract_transducer(reducible)
    return transduce(xf, SideEffect(eff), nothing, coll; kwargs...)
end

"""
    ifunreduced(f, [x])

Equivalent to [`unreduced(x)`](@ref unreduced) if `x` is a
[`Reduced`](@ref); otherwise run `f(x)`.  Return a curried version if
`x` is not provided.

See: [`foreach`](@ref).

# Examples
```jldoctest
julia> using Transducers

julia> 1 |> ifunreduced() do x
           println("called with x = ", x)
       end
called with x = 1

julia> reduced(1) |> ifunreduced() do x
           println("called with x = ", x)
       end
1
```

Notice that nothing is printed in the last example.

# Implementation

```julia
ifunreduced(f) = x -> ifunreduced(f, x)
ifunreduced(f, x::Reduced) = unreduced(x)
ifunreduced(f, x) = f(x)
```
"""
ifunreduced(f) = x -> ifunreduced(f, x)
ifunreduced(f, x::Reduced) = unreduced(x)
ifunreduced(f, x) = f(x)

"""
    Channel(xf::Transducer, itr; kwargs...)
    Channel(ed::Eduction; kwargs...)

Pipe items from an iterable `itr` processed by the transducer `xf`
through a channel.  `Channel(xf, itr)` and `Channel(eduction(xf,
itr))` are equivalent.  Note that `itr` itself can be a `Channel`.

Keyword arguments are passed to `Channel(function; kwargs...)`.

# Examples
```jldoctest
julia> using Transducers

julia> ch1 = Channel(Filter(isodd), 1:5);

julia> ch2 = Channel(Map(x -> 2x - 1), ch1);

julia> ed = eduction(Map(x -> 1:x), ch2);

julia> ch3 = Channel(Cat(), ed);

julia> foreach(PartitionBy(isequal(1)), ch3) do input
           @show input
       end;
input = [1, 1]
input = [2, 3, 4, 5]
input = [1]
input = [2, 3, 4, 5, 6, 7, 8, 9]
```
"""
Base.Channel(xform::Transducer, itr; kwargs...) =
    Channel(; kwargs...) do chan
        foreach(x -> put!(chan, x), xform, itr)
        return
    end

Base.Channel(ed::Eduction; kwargs...) =
    Channel(Transducer(ed), ed.coll; kwargs...)

Base.Channel{T}(xform::Transducer, itr, size=0; kwargs...) where {T} =
    _Channel(T, size; kwargs...) do chan
        foreach(x -> put!(chan, x), xform, itr)
        return
    end

Base.Channel{T}(ed::Eduction, size=0; kwargs...) where {T} =
    Channel{T}(Transducer(ed), ed.coll, size; kwargs...)

"""
    AdHocFoldable(foldl, [collection = nothing])

Provide a different way to fold `collection` without creating a
wrapper type.

# Arguments
- `foldl::Function`: a function that implements [`__foldl__`](@ref).
- `collection`: a collection passed to the last argument of
  `foldl`.

# Examples
```jldoctest
julia> using Transducers
       using Transducers: @next, complete
       using ArgCheck

julia> function uppertriangle(A::AbstractMatrix)
           @argcheck !Base.has_offset_axes(A)
           return AdHocFoldable(A) do rf, acc, A
               for j in 1:size(A, 2), i in 1:min(j, size(A, 1))
                   acc = @next(rf, acc, @inbounds A[i, j])
               end
               return complete(rf, acc)
           end
       end;

julia> A = reshape(1:6, (3, 2))
3×2 reshape(::UnitRange{Int64}, 3, 2) with eltype Int64:
 1  4
 2  5
 3  6

julia> collect(Map(identity), uppertriangle(A))
3-element Vector{Int64}:
 1
 4
 5

julia> function circularwindows(xs::AbstractVector, h::Integer)
           @argcheck !Base.has_offset_axes(xs)
           @argcheck h >= 0
           @argcheck 2 * h + 1 <= length(xs)
           return AdHocFoldable(xs) do rf, acc, xs
               buffer = similar(xs, 2 * h + 1)
               @inbounds for i in 1:h
                   buffer[1:h - i + 1] .= @view xs[end - h + i:end]
                   buffer[h - i + 2:end] .= @view xs[1:h + i]
                   acc = @next(rf, acc, buffer)
               end
               for i in h + 1:length(xs) - h
                   acc = @next(rf, acc, @inbounds @view xs[i - h:i + h])
               end
               @inbounds for i in 1:h
                   buffer[1:end - i] .= @view xs[end - 2 * h + i:end]
                   buffer[end - i + 1:end] .= @view xs[1:i]
                   acc = @next(rf, acc, buffer)
               end
               return complete(rf, acc)
           end
       end;

julia> collect(Map(collect), circularwindows(1:9, 2))
9-element Vector{Vector{Int64}}:
 [8, 9, 1, 2, 3]
 [9, 1, 2, 3, 4]
 [1, 2, 3, 4, 5]
 [2, 3, 4, 5, 6]
 [3, 4, 5, 6, 7]
 [4, 5, 6, 7, 8]
 [5, 6, 7, 8, 9]
 [6, 7, 8, 9, 1]
 [7, 8, 9, 1, 2]

julia> expressions(str::AbstractString; kwargs...) =
           AdHocFoldable(str) do rf, val, str
               pos = 1
               while true
                   expr, pos = Meta.parse(str, pos;
                                          raise = false,
                                          depwarn = false,
                                          kwargs...)
                   expr === nothing && break
                   val = @next(rf, val, expr)
               end
               return complete(rf, val)
           end;

julia> collect(Map(identity), expressions(\"\"\"
       x = 1
       y = 2
       \"\"\"))
2-element Vector{Expr}:
 :(x = 1)
 :(y = 2)

julia> counting = AdHocFoldable() do rf, acc, _
           i = 0
           while true
               i += 1
               acc = @next(rf, acc, i)
           end
       end;

julia> foreach(counting) do i
           @show i;
           i == 3 && return reduced()
       end;
i = 1
i = 2
i = 3
```
"""
struct AdHocFoldable{F, C} <: Foldable
    f::F
    coll::C
end

AdHocFoldable(f) = AdHocFoldable(f, nothing)

@inline __foldl__(rf, init, foldable::AdHocFoldable) =
    foldable.f(rf, init, foldable.coll)


text:
function rps()
	print("Welcome to Rock, paper, scissors! Go ahead and type your pick.\n
		r(ock), p(aper), or s(cissors)\n
		Enter 'q' to quit.\n>")
	comp_score = 0
	user_score = 0
	options = ["r","p","s"]
	new_pick() = options[rand(1:3)]
	i_win(m,y) = ((m == "r" && y == "s")|(m == "s" && y == "p")|(m == "p" && y == "r"))
	while true
		input = string(readline(STDIN)[1])
		input == "q" && break
		!ismatch(r"^[rps]",input) && begin print("Invalid guess: Please enter 'r', 'p', or 's'\n"); continue end
		answer = new_pick()
		if input == answer
			print("\nTie!\nScore still: \nyou $user_score\nme  $comp_score\n>")
			continue
		else
			i_win(answer,input) ? (comp_score += 1) : (user_score += 1)
			print(i_win(answer,input) ? "\nSorry you lose!\n" : "\nYou win!","Score: \nyou $user_score\nme  $comp_score\n>")
		end
	end
end

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module Gabor
using DSP

function pad_zeros(v, n)
    d = max(n - length(v),0)
    pre,post = fld(d,2),cld(d,2)
    [zeros(pre);v;zeros(post)]
end

function gabor{T}(t, f::Vector{T}, σ = 1.0, tol = 1e-5,
                  wf = n -> hanning(n))
    Δt = sqrt(-2σ^2*log(tol))
    Δi = round(Int, 2Δt/(t[2]-t[1]))

    w = wf(2Δi+1)

    nt = length(t)
    G = zeros(Complex{real(T)}, (2Δi+1, nt))
    for i = 1:nt
        sel = max(1,i-Δi):min(nt,i+Δi)
        G[:,i] = fftshift(fft(pad_zeros(exp(-(t[sel]-t[i]).^2/2(σ^2)).*f[sel], 2Δi+1).*w))
    end
    fftshift(collect(fftfreq(2Δi+1,1.0/(t[2]-t[1])))), G/√(2π)
end

export gabor

end


text:
"""
Convenience code for plotting result of Scalpels

For algorithm information, see Collier-Cameron, Ford, Shahaf et al. 2020

Author: Eric Ford
Date:   September 2020
"""

function make_plots_scalpels(rvs::AbstractVector{T1}, ccfs::AbstractArray{T2,2}
                ; σ_rvs::AbstractVector{T3} = ones(length(rvs)),
                max_num_basis::Integer = 3,
                v_grid::AbstractRange = 1:size(ccfs,1),
                times::AbstractVector{T4} = collect(1:length(rvs)),
                output_path::String = "", save_fig::Bool = false ) where
                { T1<:Real, T2<:Real, T3<:Real, T4<:Real  }

    @assert length(rvs) == length(σ_rvs)
    @assert length(v_grid) == size(ccfs,1)
    @assert length(rvs) == size(ccfs,2)
    @assert 1 <= max_num_basis < length(rvs)

    mean_rv = mean(rvs, weights(1.0 ./ σ_rvs.^2 ))
    rvs_centered = rvs .- mean_rv

    ccfs_minus_mean = ccfs .- mean(ccfs,dims=2)
    zvals =ccfs_minus_mean
    colorscale = cgrad(:balance)

    acfs = autocor(ccfs,0:size(ccfs,1)-1)
    acfs_minus_mean = acfs .- mean(acfs,dims=2)
    Δv_grid = convert(Float64,v_grid.step).*(0:size(acfs,1)-1)
    zvals = acfs_minus_mean
    heatmap(Δv_grid, 1:size(acfs,2),zvals',c=colorscale, clims=(-maximum(abs.(zvals)),maximum(abs.(zvals))) )
    title!("ACF(δv)-<ACF(δv)>")
    xlabel!("Δv (m/s)")
    ylabel!("Obs #")
    if save_fig
        savefig(output_path * "acf_heatmap.png")
    end

    svd_ccfs = svd(ccfs_minus_mean')
    plot(v_grid,svd_ccfs.V[:,1:max_num_basis])
    title!("CCF(v) basis functions")
    xlabel!("v (m/s)")
    if save_fig
        savefig(output_path * "ccf_basis.png")
    end

    svd_acfs = svd(acfs_minus_mean')
    plot(Δv_grid,svd_acfs.V[:,1:max_num_basis])
    title!("ACF(v) basis functions")
    xlabel!("Δv (m/s)")
    if save_fig
        savefig(output_path  * "acf_basis.png")
    end

    alpha = svd_acfs.U'*rvs_centered
    idx = sortperm(abs.(alpha),rev=true)
    local rvs_clean, plt
    for num_basis in 1:max_num_basis
        U_keep = view(svd_acfs.U,:,idx[1:num_basis])
        Δrv_shape = U_keep*U_keep'*rvs_centered
        rvs_clean = rvs .- Δrv_shape
        plt = scatter(times,rvs, label="RVs orig")
        scatter!(plt, times,rvs_clean, label="RVs cleaned", legend=:bottomright)
        xlabel!("Time (d)")
        ylabel!("RV (m/s)")
        title!("Cleaning RVs via Scalpels (" * string(num_basis) * " basis vectors)")
        if save_fig
            savefig( output_path * "scalpels_rvs_" * string(num_basis) * ".png")
        end
    end
    return plt
end


text:
export support,
    domaintype,
    codomaintype,
    isdiscrete,
    iscontinuous,
    isnormalized,
    weight,
    weightfunction,
    points,
    weights,
    discrete_weight,
    AbstractLebesgueMeasure,
    LebesgueMeasure,
    UnitLebesgueMeasure,
    DomainLebesgueMeasure,
    LegendreMeasure,
    JacobiMeasure,
    ChebyshevTMeasure,
    ChebyshevMeasure,
    ChebyshevUMeasure,
    LaguerreMeasure,
    HermiteMeasure,
    GaussianMeasure,
    DiracMeasure,
    point,
    lebesguemeasure

"Supertype of measures. The support of an `AbstractMeasure{T}` is a `Domain{T}`."
abstract type AbstractMeasure{T} end

domaintype(μ::AbstractMeasure) = domaintype(typeof(μ))
domaintype(::Type{<:AbstractMeasure{T}}) where {T} = T

"What is the codomain type of the measure?"
codomaintype(μ::AbstractMeasure) = codomaintype(typeof(μ))
codomaintype(::Type{<:AbstractMeasure{T}}) where {T} = prectype(T)

prectype(::Type{<:AbstractMeasure{T}}) where {T} = prectype(T)
numtype(::Type{<:AbstractMeasure{T}}) where {T} = numtype(T)

"Is the measure normalized?"
isnormalized(μ::AbstractMeasure) = false

convert(::Type{AbstractMeasure{T}}, μ::AbstractMeasure{T}) where {T} = μ
convert(::Type{AbstractMeasure{T}}, μ::AbstractMeasure{S}) where {S,T} = similar(μ, T)

"""
A `Measure{T}` is a continuous measure that is defined in terms of a
weightfunction: `dμ = w(x) dx`.
"""
abstract type Measure{T} <: AbstractMeasure{T} end

"""
A `DiscreteMeasure` is defined in terms of a discrete set of points and an
associated weight vector.
"""
abstract type DiscreteMeasure{T} <: AbstractMeasure{T} end

"Is the measure discrete?"
isdiscrete(μ::Measure) = false
isdiscrete(μ::DiscreteMeasure) = true

"Is the measure continuous?"
iscontinuous(μ::Measure) = true
iscontinuous(μ::DiscreteMeasure) = false

## Support for continuous measures

"Return the support of the measure"
support(μ::Measure{T}) where {T} = FullSpace{T}()

"Evaluate the weight function associated with the measure."
function weight(μ::AbstractMeasure{T}, x::S) where {S,T}
    U = promote_type(S,T)
    weight(convert(AbstractMeasure{U}, μ), convert(U, x))
end

# If the argument x is a vector: ensure the element types match
function weight(μ::AbstractMeasure{Vector{S}}, x::AbstractVector{T}) where {S,T}
    U = promote_type(S,T)
    weight(convert(Measure{Vector{U}}, μ), convert(AbstractVector{U}, x))
end
# If the measure expects SVector, convert x to SVector (we can assume matching eltype)
weight(μ::AbstractMeasure{SVector{N,S}}, x::AbstractVector{T}) where {N,S,T} =
    weight(μ, convert(SVector{N,S}, x))

# Accept matching types, and matching vectors
weight(μ::AbstractMeasure{T}, x::T) where {T} = weight1(μ, x)
weight(μ::AbstractMeasure{Vector{T}}, x::AbstractVector{T}) where {T} = weight1(μ, x)
weight(μ::AbstractMeasure{SVector{N,T}}, x::SVector{N,T}) where {N,T} = weight1(μ, x)

# Check for support, then invoke unsafe_weight
weight1(μ::AbstractMeasure, x) =
    x ∈ support(μ) ? unsafe_weight(μ, x) : zero(codomaintype(μ))

weightfunction(m::AbstractMeasure) = x->weight(m, x)
unsafe_weightfunction(m::AbstractMeasure) = x->unsafe_weight(m, x)

## Support for discrete measures

points(μ::DiscreteMeasure) = μ.points
weights(μ::DiscreteMeasure) = μ.weights

support(μ::DiscreteMeasure) = points(μ)

function discrete_weight(μ::DiscreteMeasure, i)
    # Perform a bounds check and invoke unsafe_discrete_weight,
    # so that concrete measures may implement e.g. an on-the-fly formula for
    # the weights without bounds checking
    @boundscheck checkbounds(μ, i)
    unsafe_discrete_weight(μ, i)
end
checkbounds(μ::DiscreteMeasure, i) = checkbounds(points(μ), i)

function unsafe_discrete_weight(μ::DiscreteMeasure, i)
    @inbounds weights(μ)[i]
end

isnormalized(μ::DiscreteMeasure) = sum(weights(μ)) ≈ 1

length(μ::DiscreteMeasure) = length(points(μ))

## Basic measures

"Supertype of Lebesgue measures"
abstract type AbstractLebesgueMeasure{T} <: Measure{T} end

unsafe_weight(μ::AbstractLebesgueMeasure, x) = one(codomaintype(μ))

islebesguemeasure(m::AbstractMeasure) = false
islebesguemeasure(m::AbstractLebesgueMeasure) = true

"The Lebesgue measure on the space `FullSpace{T}`."
struct LebesgueMeasure{T} <: AbstractLebesgueMeasure{T}
end

LebesgueMeasure() = LebesgueMeasure{Float64}()
similar(μ::LebesgueMeasure, ::Type{T}) where {T} = LebesgueMeasure{T}()

"The Lebesgue measure on the unit interval `[0,1]`."
struct UnitLebesgueMeasure{T} <: AbstractLebesgueMeasure{T}
end

UnitLebesgueMeasure() = UnitLebesgueMeasure{Float64}()
similar(μ::UnitLebesgueMeasure, ::Type{T}) where {T <: Real} = LebesgueMeasure{T}()
support(μ::UnitLebesgueMeasure{T}) where {T} = UnitInterval{T}()

isnormalized(μ::UnitLebesgueMeasure) = true

"Lebesgue measure supported on a general domain."
struct DomainLebesgueMeasure{T} <: AbstractLebesgueMeasure{T}
    domain  ::  Domain{T}
end

similar(μ::DomainLebesgueMeasure, ::Type{T}) where {T} = DomainLebesgueMeasure{T}(μ.domain)
support(m::DomainLebesgueMeasure) = m.domain

#################
# Applications
#################

"A point measure"
struct DiracMeasure{T} <: DiscreteMeasure{T}
    point   ::  T
end

similar(μ::DiracMeasure, ::Type{T}) where {T} = DiracMeasure{T}(μ.point)

point(μ::DiracMeasure) = μ.point
support(μ::DiracMeasure) = Point(μ.point)
isnormalized(μ::DiracMeasure) = true
unsafe_weight(μ::DiracMeasure, x) = convert(codomaintype(μ), Inf)

## Some widely used measures associated with orthogonal polynomials follow

"The Lebesgue measure on the interval `[-1,1]`."
struct LegendreMeasure{T} <: AbstractLebesgueMeasure{T}
end
LegendreMeasure() = LegendreMeasure{Float64}()

similar(μ::LegendreMeasure, ::Type{T}) where {T <: Real} = LegendreMeasure{T}()
support(μ::LegendreMeasure{T}) where {T} = ChebyshevInterval{T}()

"The Jacobi measure on the interval `[-1,1]`."
struct JacobiMeasure{T} <: Measure{T}
    α   ::  T
    β   ::  T

    JacobiMeasure{T}(α = zero(T), β = zero(T)) where {T} = new(α, β)
end
JacobiMeasure() = JacobiMeasure{Float64}()
JacobiMeasure(α, β) = JacobiMeasure(promote(α, β)...)
JacobiMeasure(α::T, β::T) where {T<:AbstractFloat} = JacobiMeasure{T}(α, β)
JacobiMeasure(α::N, β::N) where {N<:Number} = JacobiMeasure(float(α), float(β))

similar(μ::JacobiMeasure, ::Type{T}) where {T <: Real} = JacobiMeasure{T}(μ.α, μ.β)
support(μ::JacobiMeasure{T}) where {T} = ChebyshevInterval{T}()
unsafe_weight(μ::JacobiMeasure, x) = (1+x)^μ.α * (1-x)^μ.β

jacobi_α(μ::JacobiMeasure) = μ.α
jacobi_β(μ::JacobiMeasure) = μ.β

"""
The `Chebyshev` or `ChebyshevT` measure is the measure on `[-1,1]` with the
Chebyshev weight `w(x) = 1/√(1-x^2)`.
"""
struct ChebyshevTMeasure{T} <: Measure{T}
end
ChebyshevTMeasure() = ChebyshevTMeasure{Float64}()

const ChebyshevMeasure = ChebyshevTMeasure

chebyshev_weight_firstkind(x) = 1/sqrt(1-x^2)

similar(μ::ChebyshevTMeasure, ::Type{T}) where {T <: Real} = ChebyshevTMeasure{T}()
support(μ::ChebyshevTMeasure{T}) where {T} = ChebyshevInterval{T}()
unsafe_weight(μ::ChebyshevTMeasure, x) = chebyshev_weight_firstkind(x)

"""
The ChebyshevU measure is the measure on `[-1,1]` with the Chebyshev weight
of the second kind `w(x) = √(1-x^2).`
"""
struct ChebyshevUMeasure{T} <: Measure{T}
end
ChebyshevUMeasure() = ChebyshevUMeasure{Float64}()

chebyshev_weight_secondkind(x) = sqrt(1-x^2)

similar(μ::ChebyshevUMeasure, ::Type{T}) where {T <: Real} = ChebyshevUMeasure{T}()
support(μ::ChebyshevUMeasure{T}) where {T} = ChebyshevInterval{T}()
unsafe_weight(μ::ChebyshevUMeasure, x) = chebyshev_weight_secondkind(x)

convert(::Type{JacobiMeasure}, μ::LegendreMeasure{T}) where {T} =
    JacobiMeasure{T}(0, 0)
convert(::Type{JacobiMeasure}, μ::ChebyshevTMeasure{T}) where {T} =
    JacobiMeasure{T}(-one(T)/2, -one(T)/2)
convert(::Type{JacobiMeasure}, μ::ChebyshevUMeasure{T}) where {T} =
    JacobiMeasure{T}(one(T)/2, one(T)/2)

function convert(::Type{ChebyshevTMeasure}, μ::JacobiMeasure{T}) where {T}
    (jacobi_α(μ) ≈ -one(T)/2 && jacobi_β(μ) ≈ -one(T)/2) || throw(InexactError(:convert, ChebyshevTMeasure, μ))
    ChebyshevTMeasure{T}()
end

function convert(::Type{ChebyshevUMeasure}, μ::JacobiMeasure{T}) where {T}
    (jacobi_α(μ) ≈ one(T)/2 && jacobi_β(μ) ≈ one(T)/2) || throw(InexactError(:convert, ChebyshevUMeasure, μ))
    ChebyshevUMeasure{T}()
end

function convert(::Type{LegendreMeasure}, μ::JacobiMeasure{T}) where {T}
    (μ.α ≈ 0 && μ.β ≈ 0) || throw(InexactError(:convert, LegendreMeasure, μ))
    LegendreMeasure{T}()
end

jacobi_α(μ::LegendreMeasure{T}) where {T} = zero(T)
jacobi_β(μ::LegendreMeasure{T}) where {T} = zero(T)
jacobi_α(μ::ChebyshevTMeasure{T}) where {T} = -one(T)/2
jacobi_β(μ::ChebyshevTMeasure{T}) where {T} = -one(T)/2
jacobi_α(μ::ChebyshevUMeasure{T}) where {T} = one(T)/2
jacobi_β(μ::ChebyshevUMeasure{T}) where {T} = one(T)/2

"The generalised Laguerre measure on the halfline `[0,∞)`."
struct LaguerreMeasure{T} <: Measure{T}
    α   ::  T

    LaguerreMeasure{T}(α = zero(T)) where {T} = new(α)
end
LaguerreMeasure() = LaguerreMeasure{Float64}()
LaguerreMeasure(α::T) where {T<:AbstractFloat} = LaguerreMeasure{T}(α)
LaguerreMeasure(α) = LaguerreMeasure(float(α))

similar(μ::LaguerreMeasure, ::Type{T}) where {T <: Real} = LaguerreMeasure{T}(μ.α)
support(μ::LaguerreMeasure{T}) where {T} = HalfLine{T}()
isnormalized(m::LaguerreMeasure) = m.α == 0
unsafe_weight(μ::LaguerreMeasure, x) = exp(-x) * x^μ.α

laguerre_α(μ::LaguerreMeasure) = μ.α

"The Hermite measure with weight exp(-x^2) on the real line."
struct HermiteMeasure{T} <: Measure{T}
end
HermiteMeasure() = HermiteMeasure{Float64}()

hermite_weight(x) = exp(-x^2)

similar(μ::HermiteMeasure, ::Type{T}) where {T <: Real} = HermiteMeasure{T}()
unsafe_weight(μ::HermiteMeasure, x) = hermite_weight(x)

"The Gaussian measure with weight exp(-|x|^2/2)."
struct GaussianMeasure{T} <: Measure{T}
end
GaussianMeasure() = GaussianMeasure{Float64}()

similar(μ::GaussianMeasure, ::Type{T}) where {T} = GaussianMeasure{T}()
isnormalized(μ::GaussianMeasure) = true
unsafe_weight(μ::GaussianMeasure, x) = 1/(2*convert(prectype(μ), pi))^(length(x)/2) * exp(-norm(x)^2)

"The Lebesgue measure associated with the given domain"
lebesguemeasure(domain::UnitInterval{T}) where {T} = UnitLebesgueMeasure{T}()
lebesguemeasure(domain::ChebyshevInterval{T}) where {T} = LegendreMeasure{T}()
lebesguemeasure(domain::FullSpace{T}) where {T} = LebesgueMeasure{T}()
lebesguemeasure(domain::Domain{T}) where {T} = DomainLebesgueMeasure{T}(domain)


text:
# tests for DG functionality

# input file to modify, same for all test functions
const test_dg_inputfile = "input_vals_channel_dg.jl"
#const test_dg_moddict = Dict{String, Any}("Flux_name" => "RoeFlux", "use_DG" => true, "new_fname" => "input_vals_channel_dg")

"""
  This functino tests tests calculating fluxes as used for the DG face
  integrals.
"""
function test_dg_flux(mesh, sbp, eqn, opts)
  @testset "----- Testing DG flux -----" begin

    # test the Roe Flux
    uL = [1.0, 2.0, 3.0, 7.0]
    uR = copy(uL)
    flux_roe = zeros(4)
    flux_euler = zeros(4)
    func = EulerEquationMod.FluxDict["RoeFlux"]
    for i=1:mesh.numInterfaces
      iface = mesh.interfaces[i]
      for j=1:mesh.sbpface.numnodes
        eqn.aux_vars_bndry[1, j, i] = EulerEquationMod.calcPressure(eqn.params, uL)
        aux_vars = eqn.aux_vars_face[:, j, i]

        nrm_scaled = sview(mesh.nrm_face, :, j, i)
        EulerEquationMod.calcEulerFlux(eqn.params, uL, aux_vars, nrm_scaled, flux_euler)

        func(eqn.params, uL, uR, aux_vars, nrm_scaled, flux_roe)

        @test isapprox( flux_roe, flux_euler) atol=1e-13
      end
    end

    # now test calcFaceFlux
    fill!(eqn.flux_face, 0.0)
    for i=1:mesh.numInterfaces
      for j=1:mesh.sbpface.numnodes
        eqn.q_face[:, 1, j, i] = uL
        eqn.q_face[:, 2, j, i] = uL
        eqn.aux_vars_face[1, j, i] = EulerEquationMod.calcPressure(eqn.params, uL)
      end
    end

    EulerEquationMod.calcFaceFlux(mesh, sbp, eqn, func, mesh.interfaces, eqn.flux_face)
    for i=1:mesh.numInterfaces
      iface = mesh.interfaces[i]
      for j=1:mesh.sbpface.numnodes
        aux_vars = eqn.aux_vars_face[:, j, i]
        nrm = mesh.sbpface.normal[:, iface.faceL]

        nrm_scaled = sview(mesh.nrm_face, :, j, i)
        EulerEquationMod.calcEulerFlux(eqn.params, uL, aux_vars, nrm_scaled, flux_euler)

        @test isapprox( eqn.flux_face[:, j, i], flux_euler) atol=1e-13
      end
    end

  end  # end facts block

  return nothing
end  # end function

#test_dg_flux(mesh, sbp, eqn, opts)
add_func2!(EulerTests, test_dg_flux, test_dg_inputfile, [TAG_FLUX, TAG_SHORTTEST])

"""
  This function tests DG boundary integrals, including interpolation to
  interfaces.
"""
function test_dg_boundary(mesh, sbp, eqn, opts)

  @testset "----- Testing DG Boundary -----" begin

    EulerEquationMod.ICRho1E2U3(mesh, sbp, eqn, opts, eqn.q_vec)
    EulerEquationMod.interpolateBoundary(mesh, sbp, eqn, opts, eqn.q, eqn.q_bndry, eqn.aux_vars_bndry)
    mesh.bndry_funcs[1:end] = EulerEquationMod.BCDict["Rho1E2U3BC"](mesh, eqn)

    # check that the interpolation worked
    for i=1:mesh.numBoundaryFaces
      for j=1:mesh.sbpface.numnodes
        @test isapprox( eqn.q_bndry[:, j, i], [1.0, 0.35355, 0.35355, 2.0]) atol=1e-13
      end
    end

    uL = eqn.q_bndry[:, 1, 1]
    flux_euler = zeros(4)
    EulerEquationMod.getBCFluxes(mesh, sbp, eqn, opts)

    for i=1:mesh.numBoundaryFaces
      bndry_i = mesh.bndryfaces[i]
      for j=1:mesh.sbpface.numnodes
        eqn.aux_vars_bndry[1, j, i] = EulerEquationMod.calcPressure(eqn.params, eqn.q_bndry[:, j, i])
        aux_vars = eqn.aux_vars_bndry[:, j, i]
        nrm_scaled = sview(mesh.nrm_bndry, :, j, i)

        EulerEquationMod.calcEulerFlux(eqn.params, uL, aux_vars, nrm_scaled, flux_euler)

        @test isapprox( eqn.bndryflux[:, j, i], flux_euler) atol=1e-13
      end
    end

  end  # end facts block

  return nothing
end  # end function

#test_dg_boundary(mesh, sbp, eqn, opts)
add_func2!(EulerTests, test_dg_boundary, test_dg_inputfile, [TAG_BC, TAG_SHORTTEST])

"""
  This functions tests that a uniform flow gives zero residual
"""
function test_dg_uniform(mesh, sbp, eqn, opts)

  @testset "----- Testing Uniform Channel -----" begin

    physicsRhs(mesh, sbp, eqn, opts, eqn.res_vec, (evalResidual,))

    for i=1:mesh.numDof
      @test isapprox( eqn.res_vec[i], 0.0) atol=1e-13
    end

  end  # end facts block

  return nothing
end  # end function

#test_dg_uniform(mesh, sbp, eqn, opts)
add_func2!(EulerTests, test_dg_uniform, test_dg_inputfile, [TAG_SHORTTEST])


text:
struct CivecmI2 <: AbstractCivecm
    endogenous::Matrix{Float64}
    exogenous::Matrix{Float64}
    lags::Int64
    rankI1::Int64
    rankI2::Int64
    α::Matrix{Float64}
    ρδ::Matrix{Float64}
    Hρδ::Matrix{Float64}
    hρδ::Vector{Float64}
    τ::Matrix{Float64}
    Hτ::Matrix{Float64}
    hτ::Vector{Float64}
    τ⊥::Matrix{Float64}
    ζt::Matrix{Float64}
    llConvCrit::Float64
    maxiter::Int64
    convCount::Ref{Int64}
    method::String
    verbose::Bool
    Z0::Matrix{Float64}
    Z1::Matrix{Float64}
    Z2::Matrix{Float64}
    Z3::Matrix{Float64}
    R0::Matrix{Float64}
    R1::Matrix{Float64}
    R2::Matrix{Float64}
end

function civecmI2(
    endogenous::Matrix{Float64};
    exogenous::Matrix{Float64}=Matrix{Float64}(undef, size(endogenous, 1), 0),
    lags::Int64=2,
    rankI1::Int64=size(endogenous, 2),
    rankI2::Int64=0
    )

    ss, p = size(endogenous)
    iT = ss - lags
    pexo = size(exogenous, 2)
    p1 = p + pexo
    obj = CivecmI2(endogenous,
                    exogenous,
                    lags,
                    rankI1,
                    rankI2,
                    Matrix{Float64}(undef, p , rankI1),
                    Matrix{Float64}(undef, p1, rankI1),
                    kron(Matrix{Float64}(I, rankI1, rankI1), [zeros(rankI1 + rankI2); ones(p1 - rankI1 - rankI2)]),
                    vec([Matrix{Float64}(I, rankI1 + rankI2, rankI1); zeros(p1 - rankI1 - rankI2, rankI1)]),
                    Matrix{Float64}(undef, p1, rankI1 + rankI2),
                    Matrix{Float64}(undef, p1*(rankI1 + rankI2), p1*(rankI1 + rankI2)),
                    zeros(p1*(rankI1 + rankI2)),
                    Matrix{Float64}(undef, p1, p1 - rankI1 - rankI2),
                    Matrix{Float64}(undef, rankI1 + rankI2, p),
                    1.0e-8,
                    5000,
                    Ref(0),
                    "ParuoloRahbek",
                    false,
                    Matrix{Float64}(undef, iT, p),
                    Matrix{Float64}(undef, iT, p1),
                    Matrix{Float64}(undef, iT, p1),
                    Matrix{Float64}(undef, iT, p*(lags - 2) + pexo*(lags - 1)),
                    Matrix{Float64}(undef, iT, p),
                    Matrix{Float64}(undef, iT, p1),
                    Matrix{Float64}(undef, iT, p1))
    auxilliaryMatrices!(obj)
    # estimate!(obj)
    return obj
end

function auxilliaryMatrices!(obj::CivecmI2)
    iT, p = size(obj.Z0)
    pexo = size(obj.exogenous, 2)
    for j = 1:p
        for i = 1:iT
            obj.Z0[i,j] = obj.endogenous[i+obj.lags,j] - 2obj.endogenous[i+obj.lags-1,j] + obj.endogenous[i+obj.lags-2,j]
            obj.Z1[i,j] = obj.endogenous[i+obj.lags-1,j] - obj.endogenous[i+obj.lags-2,j]
            obj.Z2[i,j] = obj.endogenous[i+obj.lags-1,j]
        end
    end
    for k = 1:obj.lags - 2
        for j = 1:p
            for i = 1:iT
                obj.Z3[i,p*(k-1)+j] = obj.endogenous[i+obj.lags-k,j] - 2obj.endogenous[i+obj.lags-k-1,j] + obj.endogenous[i+obj.lags-k-2,j]
            end
        end
    end
    for j = 1:pexo
        for i = 1:iT
            obj.Z1[i,p+j] = obj.exogenous[i+obj.lags-1,j] - obj.exogenous[i+obj.lags-2,j]
            obj.Z2[i,p+j] = obj.exogenous[i+obj.lags-1,j]
            obj.Z3[i,p*(obj.lags-2)+j] = obj.exogenous[i+obj.lags,j] - 2obj.exogenous[i+obj.lags-1,j] + obj.exogenous[i+obj.lags-2,j]
        end
    end
    for k = 1:obj.lags - 2
        for j = 1:pexo
            for i = 1:iT
                obj.Z3[i,p*(obj.lags-2)+pexo*k+j] = obj.exogenous[i+obj.lags-k,j] - 2obj.exogenous[i+obj.lags-k-1,j] + obj.exogenous[i+obj.lags-k-2,j]
            end
        end
    end
    if size(obj.Z3, 2) > 0 && norm(obj.Z3) > size(obj.Z3, 1)*eps()
        obj.R0[:] = mreg(obj.Z0, obj.Z3)[2]
        obj.R1[:] = mreg(obj.Z1, obj.Z3)[2]
        obj.R2[:] = mreg(obj.Z2, obj.Z3)[2]
    else
        obj.R0[:] = obj.Z0
        obj.R1[:] = obj.Z1
        obj.R2[:] = obj.Z2
    end
    return obj
end

copy(obj::CivecmI2) = CivecmI2(copy(obj.endogenous),
                               copy(obj.exogenous),
                               obj.lags,
                               obj.rankI1,
                               obj.rankI2,
                               copy(obj.α),
                               copy(obj.ρδ),
                               copy(obj.Hρδ),
                               copy(obj.hρδ),
                               copy(obj.τ),
                               copy(obj.Hτ),
                               copy(obj.hτ),
                               copy(obj.τ⊥),
                               copy(obj.ζt),
                               obj.llConvCrit,
                               obj.maxiter,
                               obj.convCount,
                               obj.method,
                               obj.verbose,
                               copy(obj.Z0),
                               copy(obj.Z1),
                               copy(obj.Z2),
                               copy(obj.Z3),
                               copy(obj.R0),
                               copy(obj.R1),
                               copy(obj.R2))

function setrank(obj::CivecmI2, rankI1::Int64, rankI2::Int64)
    if rankI1 + rankI2 > size(obj.endogenous, 2)
        error("Illegal choice of rank")
    else
        p   = size(obj.endogenous, 2)
        p1  = p + size(obj.exogenous, 2)
        α   = Matrix{Float64}(undef, p , rankI1)
        ρδ  = Matrix{Float64}(undef, p1, rankI1)
        Hρδ = Matrix{Float64}(I, p1*rankI1, p1*rankI1)
        hρδ = zeros(p1*rankI1)
        τ   = Matrix{Float64}(undef, p1, rankI1 + rankI2)
        Hτ  = Matrix{Float64}(I, p1*(rankI1 + rankI2), p1*(rankI1 + rankI2))
        hτ  = zeros(p1*(rankI1 + rankI2))
        τ⊥  = Matrix{Float64}(undef, p1, p1 - rankI1 - rankI2)
        ζt  = Matrix{Float64}(undef, rankI1 + rankI2, p)

        newobj = CivecmI2(obj.endogenous,
                          obj.exogenous,
                          obj.lags,
                          rankI1,
                          rankI2,
                          α,
                          ρδ,
                          Hρδ,
                          hρδ,
                          τ,
                          Hτ,
                          hτ,
                          τ⊥,
                          ζt,
                          obj.llConvCrit,
                          obj.maxiter,
                          obj.convCount,
                          obj.method,
                          obj.verbose,
                          obj.Z0,
                          obj.Z1,
                          obj.Z2,
                          obj.Z3,
                          obj.R0,
                          obj.R1,
                          obj.R2)
    end
    return estimate!(newobj)
end

function estimate!(obj::CivecmI2)
    if obj.method == "ParuoloRahbek"
        # if obj.rankI1 == 0 || obj.rankI2 == size(obj.R0,2) - obj.rankI1 return estimate2step(obj) end
        return estimateτSwitch!(obj)
    end
    error("No method named %obj.method")
end

function estimateτSwitch!(obj::CivecmI2)
    # Timer
    tt = time()

    # Dimentions
    p = size(obj.R0, 2)
    iT, p1 = size(obj.R2)
    rs = obj.rankI1 + obj.rankI2

    # Moment matrices
    S10 = obj.R1'obj.R0/iT
    S11 = obj.R1'obj.R1/iT
    S20 = obj.R2'obj.R0/iT
    S21 = obj.R2'obj.R1/iT
    S22 = obj.R2'obj.R2/iT

    # Conditioning estimate (to summarize error magnitude during calculations)
    condS = cond([obj.R2 obj.R1] |> t -> t't)

    # Memory allocation
    Rτ      = Matrix{Float64}(undef, iT, p1)
    R1τ     = Matrix{Float64}(undef, iT, rs)
    workX   = Matrix{Float64}(undef, rs, p1)
    mX      = Matrix{Float64}(undef, iT, p1)
    workY   = Matrix{Float64}(undef, rs, p)
    mY      = Matrix{Float64}(undef, iT, p)
    α⊥      = Matrix{Float64}(undef, p , p - obj.rankI1)
    workRRR = Vector{Float64}(undef, obj.rankI1)
    ρ       = view(obj.ρδ, 1:rs, 1:obj.rankI1)
    ρort    = Matrix{Float64}(undef, rs, rs - obj.rankI1)
    δ       = view(obj.ρδ, rs+1:p1, 1:obj.rankI1)
    φ_ρδ    = Matrix{Float64}(undef, size(obj.Hρδ, 2), obj.rankI1)
    φ_τ     = Vector{Float64}(undef, size(obj.Hτ, 2))
    res     = Matrix{Float64}(undef, iT, p)
    Ω       = Matrix{Float64}(I, p, p)
    A       = Matrix{Float64}(undef, rs, rs)
    B       = S22
    C       = Matrix{Float64}(undef, rs, rs)
    D       = S11
    E       = Vector{Float64}(undef, p1*rs)
    ABCD    = Matrix{Float64}(undef, p1*rs, p1*rs)

    # Choose initial values from two step estimation procedure
    estimate2step!(obj)
    # obj.τ[:] = obj.Hτ*randn(size(obj.Hτ, 2)) + obj.hτ
    # obj.τ[:] = obj.Hτ*(obj.Hτ\(vec(obj.τ) - obj.hτ)) + obj.hτ # but choose a value within the restrited parameter space by projecting onto it

    # Algorithm
    ll  = -floatmax()
    ll0 = -floatmax()
    j = 1
    for j = 1:obj.maxiter
        if obj.verbose
            time() - tt > 1 && println("\nIteration:", j)
        end
        obj.τ⊥[:] = nullspace(obj.τ')[:,1:p1 - obj.rankI1 - obj.rankI2]
        Rτ[:,1:rs] = obj.R2*obj.τ
        Rτ[:,rs+1:end] = obj.R1*obj.τ⊥
        R1τ[:] = obj.R1*obj.τ
        workX[:], mX[:] = mreg(Rτ, R1τ)
        workY[:], mY[:] = mreg(obj.R0, R1τ)
        if j == 1
            # Initiate parameters
            obj.α[:], workRRR[:], obj.ρδ[:] = rrr(mY, mX, obj.rankI1)
            # obj.ρδ[:] = obj.Hρδ*φ_ρδ
            rmul!(obj.α, Diagonal(workRRR))
            obj.ζt[:], res[:] = mreg(obj.R0 - Rτ*obj.ρδ*obj.α', R1τ)
            Ω = res'res/iT
            if obj.verbose
                # if time() - tt > 1
                    # println("\nτ:\n", obj.τ)
                    println("ll:", loglikelihood(obj))
                # end
            end
        else
            switch!(mY, mX, obj.ρδ, obj.α, Ω, obj.Hρδ, obj.hρδ, maxiter = obj.maxiter, xtol = obj.llConvCrit)
            obj.ζt .= R1τ\(obj.R0 - Rτ*obj.ρδ*obj.α')
        end
        # ll = loglikelihood(obj)
        if obj.verbose
            if time() - tt > 1
                # println("\nτ:\n", obj.τ)
                println("Right after switching given τ\nll:", ll)
            end
        end
        if ll - ll0 < -obj.llConvCrit
            println("Old likelihood: $(ll0)\nNew likelihood: $(ll)\nIteration: $(j)")
            error("Likelihood cannot decrease")
        elseif abs(ll - ll0) < obj.llConvCrit && j > 1 # Use abs to avoid spurious stops due to noise
            obj.verbose && @printf("Convergence in %d iterations.\n", j - 1)
            obj.convCount[] = j
            break
        end
        if isnan(ll)
            @warn "nans in loglikehood. Aborting!"
            obj.convCount[] = obj.maxiter
            break
        end
        ll0 = ll
        α⊥ = nullspace(obj.α')[:,1:p-obj.rankI1]
        # A = ρ*obj.α'*(Ω\obj.α)*ρ'
        # B[:] = S22
        κ = obj.ζt*α⊥
        # C = κ*(cholesky!(α⊥'Ω*α⊥)\κ')
        # D[:] = S11
        Ωα = Ω\obj.α
        ψ = obj.τ⊥*δ + obj.τ*obj.ζt*(Ωα/(obj.α'*Ωα))
        # E[:] = S20*Ωα*ρ' - S21*ψ*obj.α'Ωα*ρ' + S10*α⊥*(cholesky!(α⊥'Ω*α⊥)\κ')
        # ABCD = kron(A,B) + kron(C,D)
        # φ_τ[:] = qr!(obj.Hτ'ABCD*obj.Hτ, Val(true))\(obj.Hτ'*(E - ABCD*obj.hτ))
        sqrtΩ = sqrt(Ω)
        tmpX = kron(sqrtΩ\obj.α*ρ', obj.R2) + kron(pinv(sqrtΩ*α⊥)'κ', obj.R1)
        φ_τ[:] = (tmpX*obj.Hτ)\(vec((obj.R0-obj.R1*ψ*obj.α')/sqrtΩ)-tmpX*obj.hτ)
        # φ_τ[:] = (obj.Hτ'tmpX'tmpX*obj.Hτ)\(obj.Hτ'tmpX'*(vec((obj.R0-obj.R1*ψ*obj.α')/sqrtΩ)-tmpX*obj.hτ))
        obj.τ[:] = obj.Hτ*φ_τ + obj.hτ

        myres = obj.R0 - obj.R2*obj.τ*ρ*obj.α' - obj.R1*(ψ*obj.α' + obj.τ*κ*((α⊥'Ω*α⊥)\α⊥'Ω))
        ll = loglikelihood(myres)
        if obj.verbose
            if time() - tt > 1
                # println("\nτ:\n", obj.τ)
                println("Rigth after estimation of τ\nll:", ll)
                tt = time()
            end
        end
        if ll - ll0 < -obj.llConvCrit
            println("Old likelihood: $(ll0)\nNew likelihood: $(ll)\nIteration: $(j)")
            error("Likelihood cannot decrease")
        elseif abs(ll - ll0) < obj.llConvCrit # Use abs to avoid spurious stops due to noise
            obj.verbose && @printf("Convergence in %d iterations.\n", j - 1)
            obj.convCount[] = j
            break
        end
        if isnan(ll)
            @warn "nans in loglikehood. Aborting!"
            obj.convCount[] = obj.maxiter
            break
        end
        ll0 = ll
    end
    return obj
end

function estimate2step!(obj)
    _, Res0 = mreg(obj.R0, obj.R1)
    _, Res1 = mreg(obj.R2, obj.R1)
    obj.α[:], vals, β = rrr(Res0, Res1, obj.rankI1)
    rmul!(obj.α, Diagonal(vals))
    Γt = obj.R1\(obj.R0 - obj.R2*β*obj.α')
    β⊥ = nullspace(β')
    ξ, vals, η = rrr((obj.R0 - obj.R1*β*bar(β)'Γt)*nullspace(obj.α'), obj.R1*β⊥, obj.rankI2)
    ξ *= Diagonal(vals)
    obj.τ[:,1:obj.rankI1] = β
    obj.τ[:,obj.rankI1 + 1:end] = β⊥*η
    obj.τ⊥[:] = β⊥*nullspace(η')
    obj.ρδ[1:obj.rankI1 + obj.rankI2,:] = Matrix{Float64}(I, obj.rankI1 + obj.rankI2, obj.rankI1)
    obj.ρδ[obj.rankI1 + obj.rankI2 + 1:end,:] = bar(β⊥*nullspace(η'))'Γt*bar(obj.α)
    obj.ζt[1:obj.rankI1,:] = bar(β)'Γt
    obj.ζt[obj.rankI1 + 1:end,:] = bar(β⊥*η)'Γt
    return obj
end

function ranktest(obj::CivecmI2)
    ip       = size(obj.endogenous, 2)
    r, s     = obj.rankI1, obj.rankI2
    ll0      = loglikelihood(setrank(obj, ip, 0))
    tmpTrace = zeros(ip, ip + 1)
    for i = 0:ip - 1
        for j = 0:ip - i
            obj.verbose && println("r=$(i), s=$(j)")
            tmpTrace[i + 1, i + j + 1] = 2 * (ll0 - loglikelihood(setrank(obj, i, j)))
        end
    end
    setrank(obj, r, s)
    tmpTrace
end

function ranktest(rng::AbstractRNG, obj::CivecmI2, reps::Int64)
    vals  = ranktest(obj)
    pvals = ranktestPvaluesSimluateAsymp(rng, obj, vals, reps)
    return (vals, pvals)
end
ranktest(obj::CivecmI2, reps::Int64) = ranktest(Random.default_rng(), obj, reps)

function ranktestPvaluesSimluateAsymp(rng::AbstractRNG, obj::CivecmI2, testvalues::Matrix, reps::Int64)
    (iT, ip) = size(obj.endogenous)
    pvals = zeros(ip, ip + 1)
    rankdist = zeros(reps)

    # Handling the progress bar
    prgr = Progress(ip*(ip + 1)>>1, 0.5, "Simulating I(2) rank test...")
    for i = 0:ip - 1
        for j = 0:ip - i
            next!(prgr)
            for k = 1:reps
                rankdist[k] = I2TraceSimulate(randn(rng, iT, ip - i), j, obj.exogenous)
            end
            pvals[i + 1, i + j + 1] = mean(rankdist .> testvalues[i + 1, i + j + 1])
        end
    end
    return pvals
end
function ranktestPvaluesBootstrap(obj::CivecmI2, testvalues::Matrix, reps::Int64)
    iT, p    = size(obj.R0)
    r, s     = obj.rankI1, obj.rankI2
    bootobj  = copy(obj)
    objres   = Matrix{Float64}(undef, iT, p)
    workres  = Matrix{Float64}(undef, iT, p)
    mm       = Vector{Float64}(undef, p)
    bootbool = BitArray(reps)
    pvals    = zeros(p, p+1)
    for i = 0:p-1
        for j = 0:p-i
            objvar = convert(VAR, setrank(obj, i, j))
            objres[:] = residuals(obj)
            mm[:] = mean(objres, 1)
            tmpval = testvalues[i+1,i+j+1]
            for k = 1:reps
                bn = randn(iT)
                for l = 1:p
                    for m = 1:iT
                        workres[m,l] = (objres[m,l] - mm[l])*bn[m] + mm[l]
                    end
                end
                bootobj.endogenous[:] = simulate(objvar, workres)
                auxilliaryMatrices(bootobj)
                bootbool[k] = 2*(loglikelihood(setrank(bootobj, p, 0)) - loglikelihood(setrank(bootobj, i, j))) > tmpval
            end
            pvals[i+1,j+i+1] = mean(bootbool)
        end
    end
    setrank(obj, r, s)
    return pvals
end

function residuals(obj::CivecmI2)
    res = obj.R0 - obj.R2*obj.τ*ρ(obj)*obj.α' - obj.R1*obj.τ⊥*δ(obj)*obj.α' - obj.R1*obj.τ*obj.ζt
    return res
end

function show(io::IO, ::MIME"text/plain", obj::CivecmI2)
    println(io, "β':")
    show(io, MIME"text/plain"(), copy(β(obj)'))

    println(io, "\n\nτ⊥δ:")
    show(IOContext(io, :compact=>true), MIME"text/plain"(), obj.τ⊥*δ(obj)) #'

    println(io, "\n\nτ':")
    show(io, MIME"text/plain"(), copy(τ(obj)'))
end

# Coefficients
α(obj::CivecmI2) = fit.α
β(obj::CivecmI2) = τ(obj)*ρ(obj)
τ(obj::CivecmI2) = obj.τ
# τ(obj::CivecmI2) = copy(obj.τ)
ρ(obj::CivecmI2) = obj.ρδ[1:obj.rankI1+obj.rankI2,:]
δ(obj::CivecmI2) = obj.ρδ[obj.rankI1+obj.rankI2+1:end,:]


text:
using FITSIO
using Statistics
using DelimitedFiles

# read in SOHO data with long timescale variations filtered out
function read_soho(path)
    red = FITS("$path/SPM_red_intensity_series.fits")
    blue = FITS("$path/SPM_blue_intensity_series.fits")
    green = FITS("$path/SPM_green_intensity_series.fits")

    tstart = read_header(blue[1])["JDSTART"]
    tend = read_header(blue[1])["JDEND"]

    red = read(red[1])
    blue = read(blue[1])
    green = read(green[1])

    time = (collect(LinRange(tstart, tend, length(blue))) .- tstart)
    return time, red, green, blue
end

function write_transit(t, r, g, b, path)
    open(path, "w") do io
        writedlm(io, [t r g b])
    end
end

function DefineCeleriteDist()

    extensions = quote

        ## Load needed packages and import methods to be extended
        using Distributions
        import Distributions: length, insupport, _logpdf

        ## Type declaration
        mutable struct CeleriteDist <: ContinuousMultivariateDistribution
            t::Array{Real, 1}
            c::Array{Real, 1}
            log_s::Real
            log_q::Real
            log_w::Real
            log_wn::Array{Real, 1}
            mu::Array{Real, 1}
            gp::celerite.Celerite

            function CeleriteDist(t, log_s, log_q, log_w; log_wn=[-12], c=[1], mu=zeros(Base.length(t)))
                gp = celerite.Celerite(celerite.SHOTerm(log_s, log_q, log_w), broadcast(*, c, c'))
                new(t, c, log_s, log_q, log_w, log_wn, mu, gp)
            end
        end

        ## The following method functions must be implemented

        ## Dimension of the distribution
        length(d::CeleriteDist) = length(d.t)*length(d.c)

        ## Logical indicating whether x is in the support
        function insupport(d::CeleriteDist, x::AbstractVector{T}) where {T<:Real}
            length(d) == length(x) && all(isfinite.(x))
        end

        function _logpdf(d::CeleriteDist, x::AbstractVector{T}) where {T<:Real}
            m = Base.length(d.c)
            n = length(d.t)
            wn_vec = zeros(n*m)
            for i in 1:m
                wn_vec[i:m:end] .= exp(d.log_wn[i])
            end
            try
                celerite.compute!(d.gp, d.t, wn_vec)
                z = zeros(length(d.t)*m)
                for i in 1:m
                    z[i:m:end] .= x[i:m:end] .- d.mu
                end
                return celerite.log_likelihood(d.gp, z)
            catch
                return -Inf
            end
        end

        function _rand!(d::CeleriteDist, x::AbstractVector{T}) where {T<:Real}
            m = Base.length(d.c)
            n = length(d.t)
            wn_vec = zeros(n*m)
            for i in 1:m
                wn_vec[i:m:end] .= exp(d.log_wn[i])
            end
            celerite.compute!(d.gp, d.t, wn_vec)
            u = randn(Base.length(d.t)*Base.length(d.c))
            x[:] = celerite.simulate_gp(d.gp, u)
            for i in 1:m
                x[i:m:end] .+= d.mu
            end
        end

        mean(d::CeleriteDist) = convert(Vector, d.mu)

    end
    return extensions
end


text:
include("header.jl")

#=
# Submerged surfaces
In this notebook, we will learn how to calculate force and moment on
flat surfaces immersed in static fluids
=#

# ### Set up the module
using MAE103
#-
using Plots

#=
# Example
A circular gate of 3 m diameter is located in an inclined wall (60 degrees
from horizontal) of a large reservoir filled with water. A horizontal shaft runs through
the center of the gate. The shaft is at a depth of 12 m. The other side
of the gate is exposed to air at atmospheric pressure.
* What is the force on the gate?
* Where does this force effectively act?
* What is the moment applied by the water on the gate shaft?
=#
#-
#=
Here, we are working with a simple, symmetric shape (a circle), and we
know that the centroid of a circle is at its center. Since the shaft runs
through the center, we know that the centroid is at a depth of 12 m.
=#
hc = Depth(12)
d = Diameter(3)
Θ = 60*π/180

#=
The force is easy to compute, since it is just the specific weight times
area times depth.

$$F_R = \gamma h_c A$$

(The atmospheric pressure contributes nothing in this problem, since it
acts the same on both sides.)
=#
A = Area(Circle,d/2)
FR = Force(SpecificWeight(Water)*A*hc)

#=
To find the center of pressure, we use the formula

$$y_R = \dfrac{I_{xc}}{y_c A} + y_c$$

This the location along the $y$ axis, which runs from the free surface downward, parallel
to the plate. In the formula, $I_{xc}$ is the second area moment about the $x$ axis running
through the centroid of the plate, and $y_c$ is the $y$ coordinate of the centroid.

But $y_c$ is simply $h_c/\sin\theta$:
=#
yc = Length(hc/sin(Θ))

#=
and $I_{xc}$ we can compute with a handy function:
=#
Ixc = SecondAreaMomentX(Circle,d/2)

#=
Now put them together to get the center of pressure
=#
yR = Length(Ixc/(yc*A) + yc)

#=
How far is this from the centroid?
=#
yR - yc

#=
Remember, it is lower than the centroid, because there are higher pressures
acting below the centroid than above it.

To determine the moment $M_R$ about the shaft, we simply need to remember that the
force is acting at the center of pressure, so the moment is equal to the force times the distance
of this center of pressure from the shaft (i.e., from the centroid)
=#
MR = Moment((yR - yc)*FR)

#=
#### A follow-up question
What happens to these values as the plate gets deeper?

Obviously force gets larger. What about $y_R$ and $M_R$? Let's set the depth
at 20 m, and recalculate:
=#
hc = Depth(20)
FR = Force(SpecificWeight(Water)*A*hc)
yc = Length(hc/sin(Θ))
yR = Length(Ixc/(yc*A) + yc)
yR - yc
#=
so it gets a little closer to the centroid as the plate gets lower. What about
the moment?
=#
MR = Moment((yR - yc)*FR)
#=
It's exactly the same! In fact, **the moment about the centroid does not
depend on how deep the plate is**.
=#

#=
#### Another follow-up question
What happens if the plate becomes more vertical? We'll set $h_c$ to 12 m again
and make $Θ = 90^\circ$. Only $y_c$ changes:
=#
hc = Depth(12)
Θ = 90*π/180
FR = Force(SpecificWeight(Water)*A*hc)
#=
Note that **force is not affected** by a change in angle, since only depth
of the centroid matters.
=#
yc = Length(hc/sin(Θ))
yR = Length(Ixc/(yc*A) + yc)
yR - yc

#=
So the center of pressure is a little further now from the centroid than
when it was at 60 degrees.

The moment:
=#
MR = Moment((yR - yc)*FR)

#=
So **moment gets larger as the surface becomes more vertical**. This is because the
center of pressure has moved away from the centroid. This is intuitive if you
think about the other extreme: a completely horizontal surface. In that situation,
the pressure acts uniformly across the plate, so there is **zero net moment**
on a horizontal surface.
=#


text:
using Documenter, Surrogates
makedocs(
    sitename="Surrogates.jl",
    pages = [
    "index.md"
    "User guide" => [
        "Samples" => "samples.md",
        "Surrogates" => "surrogate.md",
        "Optimization" => "optimizations.md"
        ]
    "Tutorials" => [
        "Basics" => "tutorials.md",
        "Custom Kriging with Stheno" => "stheno.md"
        ]
    "Contributing" => "contributing.md"
    ]
)

deploydocs(
   repo = "github.com/JuliaDiffEq/Surrogates.jl.git",
)


text:
# want to test and benchmark easily all the time evolution methods
# will do this for both StaticArrays and non StaticArrays, (will assume some 5 spin system to make it sufficiently big)
using BenchmarkTools

# Static Stuff!

SSxI2 = SMatrix{4,4}(kron(SSx, I(2)))
SSyI2 = SMatrix{4,4}(kron(SSy, I(2)))
I4 = SMatrix{4,4}(I(4) .+ 0.0im)

const Hsys_static = SMatrix{4,4}(kron(SSx, SSz))
const HCtrl_static2 = [SSxI2, SSyI2]

# Normal Array (small)

SSxI2_mutable = Array(SSxI2)
SSyI2_mutable = Array(SSyI2)
I4_mutable = Array(I4)

const Hsys_mutable = Array(SMatrix{4,4}(kron(SSx, SSz)))
const HCtrl_mutable2 = Array.([SSxI2, SSyI2])

# general settings
n_slices = 500
n_pulses = 2
input_arr = rand(2, 500)

T = 2.0
dt = T / n_slices

outmut = [similar(Hsys_mutable) for i = 1:n_slices]
outstat = [similar(Hsys_static) for i = 1:n_slices]

# benchmark static
BenchmarkTools.@benchmark pw_evolve(
    $Hsys_static,
    $HCtrl_static2,
    $input_arr,
    $n_pulses,
    $dt,
    $n_slices,
    $I4,
)

# benchmark mutable
BenchmarkTools.@benchmark pw_evolve(
    $Hsys_mutable,
    $HCtrl_mutable2,
    $input_arr,
    $n_pulses,
    $dt,
    $n_slices,
    $I4_mutable,
)

# benchmark static _T
BenchmarkTools.@benchmark pw_evolve_T(
    $Hsys_static,
    $HCtrl_static2,
    $input_arr,
    $n_pulses,
    $dt,
    $n_slices,
    $I4,
)

# benchmark mutable
BenchmarkTools.@benchmark pw_evolve_T(
    $Hsys_mutable,
    $HCtrl_mutable2,
    $input_arr,
    $n_pulses,
    $dt,
    $n_slices,
    $I4_mutable,
)

# benchmark static

# benchmark mutable
BenchmarkTools.@benchmark QuOptimalControl.pw_prop_save!(
    $Hsys_mutable,
    $HCtrl_mutable2,
    $input_arr,
    $n_pulses,
    $n_slices,
    $dt,
    $outmut,
)

BenchmarkTools.@benchmark QuOptimalControl.pw_prop_save!(
    $Hsys_mutable,
    $HCtrl_mutable2,
    $input_arr,
    $n_pulses,
    $n_slices,
    $dt,
    $outmut,
)

BenchmarkTools.@benchmark QuOptimalControl.pw_prop_save!(
    $Hsys_mutable,
    $HCtrl_mutable2,
    $input_arr,
    $n_pulses,
    $n_slices,
    $dt,
    $outmut,
)

BenchmarkTools.@benchmark QuOptimalControl.pw_prop_save!(
    $Hsys_mutable,
    $HCtrl_mutable2,
    $input_arr,
    $n_pulses,
    $n_slices,
    $dt,
    $outmut,
)

BenchmarkTools.@benchmark QuOptimalControl.pw_prop_save!(
    $Hsys_mutable,
    $HCtrl_mutable2,
    $input_arr,
    $n_pulses,
    $n_slices,
    $dt,
    $outmut,
)

BenchmarkTools.@benchmark QuOptimalControl.pw_prop_save!(
    $Hsys_static,
    $HCtrl_static2,
    $input_arr,
    $n_pulses,
    $n_slices,
    $dt,
    $outstat,
)

QuOptimalControl.pw_prop_save!(
    $Hsys_static,
    $HCtrl_static2,
    $input_arr,
    $n_pulses,
    $n_slices,
    $dt,
    $outstat,
)

# now we test some of the stuff that actually gets used
BenchmarkTools.@benchmark pw_evolve_save(
    $Hsys_static,
    $HCtrl_static2,
    $input_arr,
    $n_pulses,
    $dt,
    $n_slices,
)

BenchmarkTools.@benchmark pw_evolve_save(
    $Hsys_mutable,
    $HCtrl_mutable2,
    $input_arr,
    $n_pulses,
    $dt,
    $n_slices,
)

BenchmarkTools.@benchmark pw_evolve_save(
    $Hsys_mutable,
    $HCtrl_mutable2,
    $input_arr,
    $n_pulses,
    $dt,
    $n_slices,
)


text:
include("PointsAndVectors/PointsAndVectors.jl")
include("LinePrimitives/LinePrimitives.jl")
include("SurfacePrimitives/SurfacePrimitives.jl")
include("VolumePrimitives/VolumePrimitives.jl")
include("Meshing.jl")
include("CSG/CSG.jl")


text:
module Eigen

using LinearAlgebra

abstract type IterativeStoppingCriteria
end

mutable struct ExactlyNIterations <: IterativeStoppingCriteria
    n::Integer
end

mutable struct Residual <: IterativeStoppingCriteria
    value
end

mutable struct CompositeStoppingCriteria <: IterativeStoppingCriteria
    criteria
end

mutable struct EigenEstimates
    eigenvalue::Number
    eigenvector::AbstractVector

    EigenEstimates(Matrix::AbstractMatrix, Guess::AbstractVector) = new(rayleighquotient(Matrix, Guess), Guess)
    EigenEstimates(val::Number, vec::AbstractVector) = new(val, vec)
end

function normalizestrategy(x::AbstractVector)
    return LinearAlgebra.normalize!(x)
end

function normalizestrategy(x::AbstractVector{<:Integer})
    return LinearAlgebra.normalize(x)
end

function normalizestrategy(estimates::EigenEstimates)
    estimates.eigenvector = normalizestrategy(estimates.eigenvector)
    return estimates
end

function stoppingcriteria(TargetResidual::Real)
    Residual(TargetResidual)
end

function stoppingcriteria(TargetResidual::Real, MaxIterations::Integer)
    CompositeStoppingCriteria([ExactlyNIterations(MaxIterations), Residual(TargetResidual)])
end

mutable struct EigenResult
    eigenvalue::Number
    eigenvector::AbstractVector
    isverified::Bool

    EigenResult(Estimates::EigenEstimates, isverified::Bool) = new(Estimates.eigenvalue, Estimates.eigenvector, isverified)
end

function shouldstop(N::ExactlyNIterations, IterationsExecuted::Integer, Matrix::AbstractMatrix, CurrentValues)
    return IterationsExecuted >= N.n, false
end

function shouldstop(R::Residual, IterationsExecuted::Integer, Matrix::AbstractMatrix, V::AbstractVector)
    estimate = rayleighquotient(Matrix, V)
    diff = Matrix * V - estimate * V
    return LinearAlgebra.norm(diff) <= R.value, true
end

function shouldstop(R::Residual, IterationsExecuted::Integer, Matrix::AbstractMatrix, CurrentValues::EigenEstimates)
    diff = Matrix * CurrentValues.eigenvector - CurrentValues.eigenvalue * CurrentValues.eigenvector
    return LinearAlgebra.norm(diff) <= R.value, true
end

function shouldstop(Composite::CompositeStoppingCriteria, IterationsExecuted::Integer, Matrix::AbstractMatrix, CurrentValues)
    for c in Composite.criteria
        stop, foundresult = shouldstop(c, IterationsExecuted, Matrix, CurrentValues)
        if stop
            return stop, foundresult
        end
    end

    return false, false
end

function rayleighquotient(Matrix::AbstractMatrix, X::AbstractVector)
    LinearAlgebra.transpose(X) * Matrix * X
end

function iterationmethod(Matrix::AbstractMatrix,
    StartingValues,
    IterationAction::Function,
    StoppingCriteria::IterativeStoppingCriteria,
    overflowStrategy,
    makeUpperHessenberg = true)

    if (makeUpperHessenberg)
        Matrix = LinearAlgebra.UpperHessenberg(LinearAlgebra.hessenberg(Matrix).H)
    end

    iterationmethodinternal(Matrix,
        StartingValues,
        IterationAction,
        StoppingCriteria,
        overflowStrategy)
end

function iterationmethod(Matrix::LinearAlgebra.UpperHessenberg,
    StartingValues,
    IterationAction::Function,
    StoppingCriteria::IterativeStoppingCriteria,
    overflowStrategy)

    iterationmethodinternal(Matrix,
        StartingValues,
        IterationAction,
        StoppingCriteria,
        overflowStrategy)

end

function iterationmethodinternal(Matrix::AbstractMatrix,
    StartingValues,
    IterationAction::Function,
    StoppingCriteria::IterativeStoppingCriteria,
    overflowStrategy)

    current = StartingValues
    iteration = 0

    stop, foundresult = shouldstop(StoppingCriteria, iteration, Matrix, current)
    while !stop
        current = IterationAction(Matrix, current, iteration)
        current = overflowStrategy(current)
        iteration += 1
        stop, foundresult = shouldstop(StoppingCriteria, iteration, Matrix, current)
    end

    return iterationresult(Matrix, current, foundresult)
end

function iterationresult(Matrix::AbstractMatrix, V::AbstractVector, FoundResult::Bool)
    EigenResult(EigenEstimates(Matrix, V), FoundResult)
end

function iterationresult(Matrix::AbstractMatrix, Estimates::EigenEstimates, FoundResult::Bool)
    EigenResult(Estimates, FoundResult)
end

function powermethod(Matrix::AbstractMatrix,
    Guess::AbstractVector,
    StoppingCriteria::IterativeStoppingCriteria,
    overflowstrategy = normalizestrategy)

    iterationmethod(Matrix,
        Guess,
        (m,c,i) -> m * c,
        StoppingCriteria,
        overflowstrategy)
end

function inverseiteration(Matrix::AbstractMatrix,
    Shift,
    Guess::AbstractVector,
    StoppingCriteria::IterativeStoppingCriteria,
    overflowstrategy = normalizestrategy)

    iterationmethod(Matrix - LinearAlgebra.UniformScaling(Shift),
        LinearAlgebra.normalize(Guess),
        (m, c, i) -> m \ c,
        StoppingCriteria,
        overflowstrategy)
end

function rayleighiteration(Matrix::AbstractMatrix,
    Guess::AbstractVector,
    StoppingCriteria::IterativeStoppingCriteria,
    overflowstrategy = normalizestrategy)

    iterationmethod(Matrix,
        EigenEstimates(Matrix, Guess),
        function(m,c,i)
            vecestimate = (m - LinearAlgebra.UniformScaling(c.eigenvalue)) \ c.eigenvector
            valestimate = LinearAlgebra.transpose(vecestimate) * m * vecestimate / (LinearAlgebra.transpose(vecestimate) * vecestimate)
            EigenEstimates(valestimate, vecestimate)
        end,
        StoppingCriteria,
        overflowstrategy)
end

mutable struct ArnoldiResult
    H::AbstractMatrix
    numColumns::Integer
    wasStopped::Bool
end

function arnoldi(Matrix::AbstractMatrix,
    Guess::AbstractVector,
    NumColumns::Integer)

    rows = size(Matrix, 1)
    columns = size(Matrix, 2)

    isComplete = (NumColumns == columns)

    Q = zeros(rows, NumColumns)
    H = zeros(rows, NumColumns)

    q = normalize(Guess)
    Q[1:rows, 1] = q

    usedColumns = 0
    wasStopped = false
    for i = 1:NumColumns
        v = Matrix * Q[1:rows, i]

        numJs = i
        if (i == columns)
            numJs = numJs - 1
        end

        for j = 1:numJs
            H[j, i] = transpose(Q[1:rows, j]) * v
            v = v - H[j, i] * Q[1:rows, j]
        end

        H[numJs + 1, i] = norm(v)

        usedColumns = usedColumns + 1
        if (H[numJs + 1, i] < 0.0000001)
            wasStopped = true
            break
        end

        if (i < NumColumns)
            Q[1:rows, i + 1] = v / H[numJs + 1, i]
        end
    end

    return ArnoldiResult(H, usedColumns, wasStopped)
end

function QRAlgorithm(Matrix::AbstractMatrix)
    iterations = 10
    for i = 1:iterations
        factors = LinearAlgebra.qr(Matrix)
        Matrix = factors.R * factors.Q
    end

    Matrix
end

function QRAlgorithmWithShifts(Matrix::AbstractMatrix,
    shifts)

    numShifts = length(shifts)
    numRows = size(Matrix, 1)

    for i = 1:numShifts
        shiftMatrix = shifts[i] * Matrix(I, numRows, numRows)

        factors = LinearAlgebra.qr(Matrix) - shiftMatrix
        Matrix = factors.R * factors.Q + shiftMatrix
    end

    Matrix
end

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 ErrorDetails <: SwaggerModel
    code::Any # spec type: Union{ Nothing, String } # spec name: code
    target::Any # spec type: Union{ Nothing, String } # spec name: target
    message::Any # spec type: Union{ Nothing, String } # spec name: message

    function ErrorDetails(;code=nothing, target=nothing, message=nothing)
        o = new()
        validate_property(ErrorDetails, Symbol("code"), code)
        setfield!(o, Symbol("code"), code)
        validate_property(ErrorDetails, Symbol("target"), target)
        setfield!(o, Symbol("target"), target)
        validate_property(ErrorDetails, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        o
    end
end # type ErrorDetails

const _property_map_ErrorDetails = Dict{Symbol,Symbol}(Symbol("code")=>Symbol("code"), Symbol("target")=>Symbol("target"), Symbol("message")=>Symbol("message"))
const _property_types_ErrorDetails = Dict{Symbol,String}(Symbol("code")=>"String", Symbol("target")=>"String", Symbol("message")=>"String")
Base.propertynames(::Type{ ErrorDetails }) = collect(keys(_property_map_ErrorDetails))
Swagger.property_type(::Type{ ErrorDetails }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ErrorDetails[name]))}
Swagger.field_name(::Type{ ErrorDetails }, property_name::Symbol) =  _property_map_ErrorDetails[property_name]

function check_required(o::ErrorDetails)
    true
end

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