text:
# --- Types

"""
    Reduced

The type signaling transducible processes to abort.

!!! note
    Call [`reduced`](@ref) function for aborting the transducible
    process since [`reduced`](@ref) makes sure `x` is not doubly
    wrapped.  `Reduced` is meant to be used as `x isa Reduced` for
    checking if the result from [`transduce`](@ref) is due to early
    termination.

See [`reduced`](@ref), [`unreduced`](@ref).

# Examples
```jldoctest
julia> using Transducers

julia> function step_demo(y, x)
           if x > 5
               return reduced(y)
           else
               return y + x
           end
       end;

julia> result = transduce(Map(identity), Completing(step_demo), 0, 1:10)
Reduced(15)

julia> result isa Reduced
true

julia> unreduced(result)
15

julia> result = transduce(Map(identity), Completing(step_demo), 0, 1:4)
10

julia> result isa Reduced
false

julia> unreduced(result)
10
```
"""
struct Reduced{T}
    value::T
end

Base.:(==)(x::Reduced, y::Reduced) = x.value == y.value

Base.show(io::IO, x::Reduced) = _default_show(io, x)

isreduced(::Reduced) = true
isreduced(::Any) = false

"""
    reduced([x = nothing])

Stop transducible process with the final value `x` (default:
`nothing`).  Return `x` as-is if it already is a `reduced` value.

See [`Reduced`](@ref), [`unreduced`](@ref).

$(_thx_clj("ensure-reduced"))

# Examples
```jldoctest
julia> using Transducers

julia> foldl(Enumerate(), "abcdef"; init=0) do y, (i, x)
           if x == 'd'
               return reduced(y)
           end
           return y + i
       end
6

julia> foreach(Enumerate(), "abc") do (i, x)
           println(i, ' ', x)
           if x == 'b'
               return reduced()
           end
       end;
1 a
2 b
```
"""
reduced(x::Reduced) = x
reduced(x) = Reduced(x)
reduced() = reduced(nothing)

"""
    unreduced(x)

Unwrap `x` if it is a `Reduced`; do nothing otherwise.

See [`Reduced`](@ref), [`reduced`](@ref).

$(_thx_clj("unreduced"))
"""
unreduced(x::Reduced) = x.value
unreduced(x) = x

"""
    @return_if_reduced expr

It transforms the given expression to:

```julia
val = expr
val isa Reduced && return val
val
```

See also [`@next`](@ref).

!!! compat "Transducers.jl 0.3"

    In v0.2, the calling convention was `@return_if_reduced
    complete(rf, val)` and it was transformed to `val isa Reduced &&
    return reduced(complete(rf, unreduced(val)))`.  For the rationale
    behind the change, see [this commit
    message](https://github.com/JuliaFolds/Transducers.jl/commit/74f8961fea97b746cb097b27aa5a5761e9bf4dae).

# Examples
```jldoctest; filter = [r"(var\\")?#[0-9]+#val(\\")?", r"#=.*?=#"]
julia> using Transducers: @return_if_reduced

julia> @macroexpand @return_if_reduced f(x)
quote
    #158#val = f(x)
    #= ... =#
    begin
        #158#val isa Transducers.Reduced && return #158#val
        #= ... =#
        #158#val
    end
end
```
"""
macro return_if_reduced(ex)
    code = :(val isa Reduced && return val; val)
    if ex isa Expr && ex.head == :call && length(ex.args) == 3
        val = esc(ex)
        return quote
            if $(esc(ex.args[1])) === complete
                error("""
                Calling `@return_if_reduced complete(rf, val)` is now an error.
                Please use `@return_if_reduced val`.
                """)
            end
            val = $val
            $code
        end
    end
    val = esc(ex)
    :(val = $val; $code)
end

"""
    @next(rf, state, input)

It is expanded to

```julia
result = next(rf, state, input)
result isa Reduced && return result
result
```

This is usually the best way to call `next` as checking for `Reduced`
is required to support early termination.

See also: [`next`](@ref), [`Reduced`](@ref), [`@return_if_reduced`](@ref).
"""
macro next(rf, state, input)
    quote
        result = next($(esc.((rf, state, input))...))
        result isa Reduced && return result
        result
    end
end

abstract type Transducer end
abstract type AbstractFilter <: Transducer end

struct Composition{XO <: Transducer, XI <: Transducer} <: Transducer
    outer::XO
    inner::XI
end

struct IdentityTransducer <: Transducer end

has(xf::Transducer, T::Type) = xf isa T
has(xf::Composition, T::Type) = has(xf.outer, T) || has(xf.inner, T)

Base.broadcastable(xf::Transducer) = Ref(xf)

"""
    Transducer

The abstract type for transducers.

A transducer `xf` can be used as both iterator transformation
`xf(itr)` and reducing function transformation `xf'(rf)`.

See also [`adjoint`](@ref) for `xf'(rf)`.

!!! compat "Transducers.jl 0.4.39"
    The call overload `xf(rf)` requires Transducers.jl 0.4.39 or later.

!!! note
    The call overload `xf(rf)` requires Julia 1.3 or later.
    For older Julia versions, use [`eduction`](@ref).

# Examples
```jldoctest
julia> using Transducers

julia> xs = Map(inv)(2:2:4)
2-element StepRange{Int64, Int64} |>
    Map(inv)

julia> collect(xs)
2-element Vector{Float64}:
 0.5
 0.25

julia> rf = Map(inv)'(+)
Reduction(
    Map(inv),
    BottomRF(
        +))

julia> rf(1, 4)  # +(1, inv(4))
1.25
```
"""
Transducer

"""
    AbstractFilter <: Transducer

The abstract type for filter-like transducers.
"""
AbstractFilter

abstract type AbstractReduction{innertype} <: Function end

if VERSION >= v"1.3"  # post https://github.com/JuliaLang/julia/pull/31916
    @inline (rf::AbstractReduction)(state, input) = next(rf, state, input)
end

InnerType(::Type{<:AbstractReduction{T}}) where T = T

ConstructionBase.constructorof(::Type{T}) where {T <: AbstractReduction} = T

"""
    Transducers.inner(rf::R_)

Return the inner reducing function of `rf`.
"""
inner(rf::AbstractReduction) = rf.inner

"""
    Transducers.xform(rf::R_{X}) -> xf :: X

Return the transducer `xf` associated with `rf`.  Returned transducer
`xf` is "atomic"; i.e., it is not a `Composition` transducer type.
"""
xform(rf::AbstractReduction) = rf.xform

has(rf::AbstractReduction, T::Type{<:Transducer}) = has(Transducer(rf), T)

struct BottomRF{F} <: AbstractReduction{F}
    inner::F
end

ensurerf(rf::AbstractReduction) = rf
ensurerf(f) = BottomRF(f)

# Not calling rf.inner(result, input) etc. directly since it can be
# `Completing` etc.
start(rf::BottomRF, result) = start(inner(rf), result)
@inline next(rf::BottomRF, result, input) = next(inner(rf), result, input)
complete(rf::BottomRF, result) = complete(inner(rf), result)
combine(rf::BottomRF, a, b) = combine(inner(rf), a, b)

Transducer(::BottomRF) = IdentityTransducer()

as(rf::T, ::Type{T}) where T = rf
as(rf, T::Type) = as(inner(rf), T)

# In clojure a reduction function is one of signature
# whatever, input -> whatever
#
# The `Reduction` type corresponds to such a function, but keeps extra information:
# * `xform` and `inner` are a decomposition of the reduction function into
#   a transducer `xform` and an inner reduction function `inner`.
#   `inner` can be either a `Reduction` or a function with arity-2 and arity-1 methods
#
struct Reduction{X <: Transducer, I} <: AbstractReduction{I}
    xform::X
    inner::I

    Reduction{X, I}(xf, inner) where {X, I} = new{X, I}(xf, inner)

    Reduction(xf::X, inner::I) where {X <: Transducer, I} =
        if I <: AbstractReduction
            new{X, I}(xf, inner)
        else
            rf = ensurerf(inner)
            new{X, typeof(rf)}(xf, rf)
        end
end

if VERSION < v"1.3"  # pre https://github.com/JuliaLang/julia/pull/31916
    @inline (rf::Reduction)(state, input) = next(rf, state, input)
end

prependxf(rf::AbstractReduction, xf) = Reduction(xf, rf)
setinner(rf::Reduction, inner) = Reduction(xform(rf), inner)
setxform(rf::Reduction, xform) = Reduction(xform, inner(rf))

Transducer(rf::Reduction) =
    if inner(rf) isa BottomRF
        xform(rf)
    else
        Composition(xform(rf), Transducer(inner(rf)))
    end

# This is a non-ideal definition as it may not return a `Reduction`.
# Making this less non-ideal requires to replace all call/overloads of
# `Reduction` to `AbstractReduction`.
Reduction(::IdentityTransducer, inner) = ensurerf(inner)

Adapt.adapt_structure(to, rf::R) where {R <: Reduction} =
    Reduction(Adapt.adapt(to, xform(rf)), Adapt.adapt(to, inner(rf)))

"""
    Transducers.R_{X}

When defining a transducer type `X`, it is often required to dispatch
on type `rf::R_{X}` (Reducing Function) which bundles the current
transducer `xform(rf)::X` and the inner reducing function
`inner(rf)::R_`.
"""
const R_{X} = Reduction{<:X}

@inline function Reduction(xf_::Composition, f)
    xf = _normalize(xf_)
    # @assert !(xf.outer isa Composition)
    return Reduction(xf.outer, Reduction(xf.inner, f))
end

@inline _normalize(xf) = xf
@inline _normalize(xf::Composition{<:Composition}) = xf.outer |> xf.inner

"""
    f ⨟ g
    g ∘ f
    opcompose(f, g)
    compose(g, f)

Composition of transducers.

!!! compat "Transducers.jl 0.4.39"

    Transducers.jl 0.4.39 or later is required for composing
    transducers with `∘` and other operators and functions derived
    from it.

    Transducers written as `f |> g |> h` in previous versions of
    Transducers.jl can now be written as `f ⨟ g ⨟ h` (in Julia 1.5 or
    later) or `opcompose(f, g, h)`.

!!! note

    "op" in `opcompose` does not stand for _operator_; it stands for
    _opposite_.
"""
@inline Base.:∘(g::Transducer, f::Transducer) = Composition(f, g)
@inline Base.:∘(g::Transducer, f::Composition) = g ∘ f.inner ∘ f.outer
@inline Base.:∘(f::Transducer, ::IdentityTransducer) = f
@inline Base.:∘(::IdentityTransducer, f::Transducer) = f
@inline Base.:∘(f::IdentityTransducer, ::IdentityTransducer) = f
@inline Base.:∘(::IdentityTransducer, f::Composition) = f  # disambiguation

if VERSION >= v"1.3"
    (xf::Transducer)(itr) = eduction(xf, itr)
else
    Base.:|>(itr, xf::Transducer) = eduction(xf, itr)
end

"""
    ReducingFunctionTransform(xf)

The "true" transducer.
"""
struct ReducingFunctionTransform{T <: Transducer} <: Function
    xf::T
end

"""
    xf'

`xf'(rf₁)` is a shortcut for calling `reducingfunction(xf, rf₁)`.

More precisely, adjoint `xf′` of a transducer `xf` is a _reducing
function transform_ `rf₁ -> rf₂`.  That is to say, `xf'` a function
that maps a reducing function `rf₁` to another reducing function
`rf₂`.

# Examples
```jldoctest
julia> using Transducers

julia> y = Map(inv)'(+)(10, 2)
10.5

julia> y == +(10, inv(2))
true
```
"""
Base.adjoint(xf::Transducer) = ReducingFunctionTransform(xf)
Base.adjoint(rxf::ReducingFunctionTransform) = rxf.xf

(f::ReducingFunctionTransform)(rf; kwargs...) =
    reducingfunction(f.xf, rf; kwargs...)

@inline Base.:∘(f::ReducingFunctionTransform, g::ReducingFunctionTransform) = (g' ∘ f')'

"""
    reform(rf, f)

Reset "bottom" reducing function of `rf` to `f`.
"""
reform(rf::Reduction, f) = Reduction(xform(rf), reform(inner(rf), f))
reform(rf::BottomRF, f) = BottomRF(reform(inner(rf), f))
reform(::Any, f) = f

"""
    Transducers.start(rf::R_{X}, state)

This is an optional interface for a transducer.  Default
implementation just calls `start` of the inner reducing function; i.e.,

```julia
start(rf::Reduction, result) = start(inner(rf), result)
```

If the transducer `X` is stateful, it can "bundle" its private state
with `wrap`:

```julia
start(rf::R_{X}, result) = wrap(rf, PRIVATE_STATE, start(inner(rf), result))
```

where `PRIVATE_STATE` is an initial value for the private state that
can be used inside [`next`](@ref) via [`wrapping`](@ref).

See [`Take`](@ref), [`PartitionBy`](@ref), etc. for real-world examples.

Side notes: There is no related API in Clojure's Transducers.
Transducers.jl uses it to implement stateful transducers using "pure"
functions.  The idea is based on a slightly different approach taken
in C++ Transducer library [atria](https://github.com/AbletonAG/atria).
"""
start(rf, init) = initialize(init, rf)
start(rf::Reduction, result) = start(inner(rf), result)
start(rf::R_{AbstractFilter}, result) = start(inner(rf), result)

"""
    Transducers.next(rf::R_{X}, state, input)

This is the only required interface.  It takes the following form
(if `start` is not defined):

```julia
next(rf::R_{X}, result, input) =
    # code calling next(inner(rf), result, possibly_modified_input)
```

When calling `next`, it is almost always a better idea to use the
macro form [`@next`](@ref).  See the details in its documentation.

See [`Map`](@ref), [`Filter`](@ref), [`Cat`](@ref), etc. for
real-world examples.
"""
@inline next(f, result, input) = f(result, input)

# done(rf, result)

"""
    Transducers.complete(rf::R_{X}, state)

This is an optional interface for a transducer.  If transducer `X` has
some internal state, this is the last chance to "flush" the result.

See [`PartitionBy`](@ref), etc. for real-world examples.

If `start(rf::R_{X}, state)` is defined, `complete` **must** unwarp
`state` before returning `state` to the outer reducing function.

!!! compat "Transducers.jl 0.3"

    In Transducers.jl 0.2, `complete` had a fallback implementation
    to automatically call `unwrap` when `wrap` is called in `start`.
    Relying on this fallback implementation is now deprecated.
"""
complete(f, result) = f(result)
complete(rf::AbstractReduction, result) =
    # Not using dispatch to avoid ambiguity
    if ownsstate(rf, result)
        Base.depwarn(
            string(
                "`complete` for ", typeof(xform(rf)), " is not defined.",
                " Automatic implementation of `complete` method will be",
                " disabled in the future."
            ),
            :complete)
        complete(inner(rf), unwrap(rf, result)[2])
    else
        complete(inner(rf), result)
    end

"""
    Transducers.combine(rf::R_{X}, state_left, state_right)

This is an optional interface for a transducer.  If transducer `X` is
stateful (i.e., [`wrap`](@ref) is used in [`start`](@ref)), it has to
be able to combine the private states to support fold functions that
require an associative reducing function such as [`foldxt`](@ref).
Typical implementation takes the following form:

```julia
function combine(rf::R_{X}, a, b)
    #   ,---- `ua` and `ub` are the private state of the transducer `X`
    #  /  ,-- `ira` and `irb` are the states of inner reducing functions
    # /  /
    ua, ira = unwrap(rf, a)
    ub, irb = unwrap(rf, b)
    irc = combine(inner(rf), ira, irb)
    uc = # somehow combine private states `ua` and `ub`
    return wrap(rf, uc, irc)
end
```

See [`ScanEmit`](@ref), etc. for real-world examples.
"""
combine(f, a, b) = f(a, b)
combine(rf::Reduction, a, b) =
    # Not using dispatch to avoid ambiguity
    if ownsstate(rf, a)
        # TODO: make sure this branch is compiled out
        error("Stateful transducer ", xform(rf), " does not support `combine`")
    elseif ownsstate(rf, b)
        error("""
        Some thing went wrong in two ways:
        * `combine(rf, a, b)` is called but type of `a` and `b` are different.
        * `xform(rf) = $(xform(rf))` is stateful and does not support `combine`.
        """)
    else
        combine(inner(rf), a, b)
    end

"""
    is_prelude(::T)

Return `true` if it is better to tail-call when the accumulator or the
private state changes its type from `T`.
"""
is_prelude
is_prelude(_) = false
is_prelude(::InitialValues.InitialValue) = true
is_prelude(xs::Tuple) = any(map(is_prelude, xs))
is_prelude(xs::NamedTuple) = is_prelude(Tuple(xs))

privatestate(::T, state, result) where {T <: AbstractReduction} =
    privatestate(T, state, result)

struct PrivateState{T, S, R}
    state::S
    result::R

    # Rename constructor to make sure that it is always constructed
    # through the factory function:
    global privatestate(::Type{T}, state::S, result::R) where {
        T <: AbstractReduction,
        S,
        R,
    } =
        new{T, S, R}(state, result)
end
# TODO: make it a tuple-like so that I can return it as-is

ConstructionBase.constructorof(::Type{<:PrivateState{T}}) where T =
    (state, result) -> privatestate(T, state, result)

@inline psstate(ps) = ps.state
@inline psresult(ps) = ps.result
@inline setpsstate(ps, x) = @set ps.state = x
@inline setpsresult(ps, x) = @set ps.result = x

ownsstate(::Any, ::Any) = false
ownsstate(::R, ::PrivateState{T}) where {R, T} = R === T
# Using `result isa PrivateState{typeof(rf)}` makes it impossible to
# compile Extrema examples in ../examples/tutorial_missings.jl (it
# took more than 10 min).  See also:
# https://github.com/JuliaLang/julia/issues/30125

@inline is_prelude(ps::PrivateState) = is_prelude(psstate(ps)) || is_prelude(psresult(ps))

"""
    unwrap(rf, result)

Unwrap [`wrap`](@ref)ed `result` to a private state and inner result.
Following identity holds:

```julia
unwrap(rf, wrap(rf, state, iresult)) == (state, iresult)
```

This is intended to be used only in [`complete`](@ref).  Inside
[`next`](@ref), use [`wrapping`](@ref).
"""
unwrap(::T, ps::PrivateState{T}) where {T} = (psstate(ps), psresult(ps))

unwrap(::T1, ::PrivateState{T2}) where {T1, T2} =
    error("""
    `unwrap(rf1, ps)` is used for
    typeof(rf1) = $T1
    while `ps` is created by wrap(rf2, ...) where
    typeof(rf2) = $T2
    """)

# TODO: better error message with unmatched `T`

"""
    wrap(rf::R_{X}, state, iresult)

Pack private `state` for reducing function `rf` (or rather the
transducer `X`) with the result `iresult` returned from the inner
reducing function `inner(rf)`.  This packed result is typically passed
to the outer reducing function.

This is intended to be used only in [`start`](@ref).  Inside
[`next`](@ref), use [`wrapping`](@ref).

!!! note "Implementation detail"

    If `iresult` is a [`Reduced`](@ref), `wrap` actually _un_wraps all
    internal state `iresult` recursively.  However, this is an
    implementation detail that should not matter when writing
    transducers.

Consider a reducing step constructed as

    rf = opcompose(xf₁, xf₂, xf₃)'(f)

where each `xfₙ` is a stateful transducer and hence needs a private
state `stateₙ` and this `stateₙ` is constructed in each
`start(::R_{typeof(xfₙ)}, result)`.  Then, calling `start(rf,
result))` is equivalent to

```julia
wrap(rf,
     state₁,                     # private state for xf₁
     wrap(inner(rf),
          state₂,                # private state for xf₂
          wrap(inner(inner(rf)),
               state₃,           # private state for xf₃
               result)))
```

or equivalently

```julia
result₃ = result
result₂ = wrap(inner(inner(rf)), state₃, result₃)
result₁ = wrap(inner(rf),        state₂, result₂)
result₀ = wrap(rf,               state₁, result₁)
```

The inner most step function receives the original `result` as the
first argument while transducible processes such as [`foldl`](@ref)
only sees the outer-most "tree" `result₀` during the reduction.

See [`wrapping`](@ref), [`unwrap`](@ref), and [`start`](@ref).
"""
wrap(rf::T, state, iresult) where {T} = privatestate(rf, state, iresult)
wrap(rf, state, iresult::Reduced) = iresult

"""
    wrapping(f, rf, result)

Function `f` must take two argument `state` and `iresult`, and return
a tuple `(state, iresult)`.  This is intended to be used only in
[`next`](@ref), possibly with a `do` block.

```julia
next(rf::R_{MyTransducer}, result, input) =
    wrapping(rf, result) do my_state, iresult
        # code calling `next(inner(rf), iresult, possibly_modified_input)`
        return my_state, iresult  # possibly modified
    end
```

See [`wrap`](@ref), [`unwrap`](@ref), and [`next`](@ref).
"""
@inline function wrapping(f, rf, result)
    #=
    state0, iresult0 = unwrap(rf, result)
    state1, iresult1 = f(state0, iresult0)
    =#
    # `first`/`last` behaves nicer with type inference for `Union` of `Tuple`s:
    a = unwrap(rf, result)
    state0 = first(a)
    iresult0 = last(a)
    b = f(state0, iresult0)
    state1 = first(b)
    iresult1 = last(b)
    return wrap(rf, state1, iresult1)
end

unwrap_all(ps::PrivateState) = unwrap_all(psresult(ps))
unwrap_all(result) = result
unwrap_all(ps::Reduced) = Reduced(unwrap_all(unreduced(ps)))

# isexpansive(::Any) = true
isexpansive(::Transducer) = true
isexpansive(::AbstractFilter) = false
# isexpansive(rf::Reduction) = isexpansive(xform(rf)) || isexpansive(inner(rf))
isexpansive(xf::Composition) = isexpansive(xf.outer) || isexpansive(xf.inner)
# Should it be a type-level trait?

#=
iscontractive(::Any) = false
iscontractive(::AbstractFilter) = true
iscontractive(rf::Reduction) = iscontractive(xform(rf)) && iscontractive(inner(rf))
=#

# TODO: merge `isexpansive` and `OutputSize`
abstract type OutputSize end
struct SizeStable <: OutputSize end
# struct SizeExpansive <: OutputSize end
struct SizeChanging <: OutputSize end
OutputSize(::Type{<:Transducer}) = SizeChanging()
OutputSize(::Type{Composition{XO,XI}}) where {XO, XI} =
    combine_outputsize(OutputSize(XO), OutputSize(XI))

combine_outputsize(::SizeStable, ::SizeStable) = SizeStable()
combine_outputsize(::OutputSize, ::OutputSize) = SizeChanging()

# outputsize(::XF) where {XF <: Transducers} = OutputSize(XF)

# For `Eduction` (which stores `Reduction` rather than `Transducer`):
outputsize(::Type{Reduction{X,I}}) where {X, I <: Reduction} =
    combine_outputsize(OutputSize(X), outputsize(I))
outputsize(::Type{Reduction{X,I}}) where {X, I} = OutputSize(X)

struct NoComplete <: Transducer end
next(rf::R_{NoComplete}, result, input) = next(inner(rf), result, input)
complete(::R_{NoComplete}, result) = result  # don't call inner complete

"""
    Completing(function)

Wrap a `function` to add a no-op [`complete`](@ref) protocol.  Use it
when passing a `function` without unary method to [`transduce`](@ref)
etc.

$(_thx_clj("completing"))
"""
struct Completing{F} <: _Function  # Note: not a Transducer
    f::F
end

start(rf::Completing, result) = start(rf.f, result)
@inline next(rf::Completing, result, input)  = next(rf.f, result, input)
complete(::Completing, result) = result
combine(rf::Completing, a, b) = combine(rf.f, a, b)

# Apply `Completing` only on the inner-most reducing function so that
# `complete` on transducers are still called.  This is very ugly as it
# does not return a `Completing` object.  But this is required for
# allowing `foldl(reducingfunction(...), ...)` etc.:
Completing(rf::AbstractReduction) = setinner(rf, Completing(inner(rf)))
Completing(rf::BottomRF) = BottomRF(Completing(rf.inner))
Completing(f::Completing) = f

# Currently, `Completing` is recursive while `skipcomplete` is
# non-recursive.  `Completing` only skips `complete` on the inner-most
# (bottom) reducing function (i.e., _not_ including the
# transducers). `skipcomplete` skips `complete` of the outer-most
# reducing function (i.e., including all transducers).

# TODOs for `Completing` and `skipcomplete`:
# 1. Get rid of 1-arg `complete` fallback definition for bottom functions.
# 2. Deprecate `Completing` (after making it a no-op).
skipcomplete(rf::Reduction) = Reduction(NoComplete(), rf)
skipcomplete(f) = Completing(f)

# Since `Completing <: Function`, the default `show` is a bit ugly.
function Base.show(io::IO, rf::Completing)
    @nospecialize
    if rf === Completing(rf.f)
        print(io, Completing, '(', rf.f, ')')
    else
        invoke(show, Tuple{IO,Any}, io, rf)
    end
end

struct SideEffect{F}  # Note: not a Transducer
    f::F
end

# Completing(rf::SideEffect) = rf

start(rf::SideEffect, result) = start(rf.f, result)
complete(::SideEffect, result) = result
@inline next(rf::SideEffect, _, input) = rf.f(input)

"""
    right([l, ]r) -> r

It is simply defined as

```julia
right(l, r) = r
right(r) = r
```

This function is meant to be used as `step` argument for
[`foldl`](@ref) etc. for extracting the last output of the
transducers.

!!! compat "Transducers.jl 0.3"

    Initial value must be manually specified.  In 0.2, it was
    automatically set to `nothing`.

# Examples
```jldoctest
julia> using Transducers

julia> foldl(right, Take(5), 1:10)
5

julia> foldl(right, Drop(5), 1:3; init=0)  # using `init` as the default value
0
```
"""
right(l, r) = r
right(r) = r

InitialValues.@def_monoid right

identityof(::typeof(right), ::Any) = nothing
# This is just a right identity but `right` is useful for left-fold
# context anyway so I guess it's fine.

abstract type Reducible end
abstract type Foldable <: Reducible end
asfoldable(x) = x

"""
    initvalue(initializer::AbstractInitializer) -> init
    initvalue(init) -> init

Materialize the initial value if the input is an `AbstractInitializer`.
Return the input as-is if not.
"""
initvalue(x) = x
_initvalue(rf::Reduction) = initvalue(xform(rf).init)

abstract type AbstractInitializer end

# TODO: Merge `InitOf` to `AbstractInitializer`
# For `Init`, `DefaultInit` and `OptInit`
struct InitOf{IV <: SpecificInitialValue} end
(::InitOf{IV})(::OP) where {IV, OP} = IV{OP}()

# For `Broadcasting`:
Broadcast.broadcastable(f::AbstractInitializer) = Ref(f)
Broadcast.broadcastable(f::InitOf) = Ref(f)

"""
    initialize(initializer, op) -> init
    initialize(init, _) -> init

Return an initial value for `op`.  Throw an error if `initializer`
(e.g., `Init`) creates unknown initial value.

# Examples
```jldoctest
julia> using Transducers
       using Transducers: initialize

julia> initialize(Init, +)
InitialValue(+)

julia> initialize(123, +)
123

julia> unknown_op(x, y) = x + 2y;

julia> initialize(Init, unknown_op)
ERROR: IdentityNotDefinedError: `init = Transducers.Init` is specified but the identity element `InitialValue(op)` is not defined for
    op = unknown_op
[...]
```
"""
initialize(init, op) = init
initialize(f::InitOf, op) = check_init(f(op), f, op)
initialize(init::AbstractInitializer, _) = initvalue(init)

function check_init(init::SpecificInitialValue, f, op)
    InitialValues.isknown(init) || throw(IdentityNotDefinedError(op, f))
    return init
end

"""
    OnInit(f)

Call a callable `f` to create an initial value.

See also [`CopyInit`](@ref).

`OnInit` or `CopyInit` must be used whenever using in-place reduction
with [`foldxt`](@ref) etc.

# Examples
```jldoctest OnInit
julia> using Transducers

julia> xf1 = Scan(push!, [])
Scan(push!, Any[])

julia> foldl(right, xf1, 1:3)
3-element Vector{Any}:
 1
 2
 3

julia> xf1
Scan(push!, Any[1, 2, 3])
```

Notice that the array is stored in `xf1` and mutated in-place.  As a
result, second run of `foldl` contains the results from the first
run:

```jldoctest OnInit
julia> foldl(right, xf1, 10:11)
5-element Vector{Any}:
  1
  2
  3
 10
 11
```

This may not be desired.  To avoid this behavior, create an `OnInit`
object which takes a factory function to create a new initial value.

```jldoctest OnInit; filter = r"#+[0-9]+(\\(\\))?"
julia> xf2 = Scan(push!, OnInit(() -> []))
Scan(push!, OnInit(##9#10()))

julia> foldl(right, xf2, 1:3)
3-element Vector{Any}:
 1
 2
 3

julia> foldl(right, xf2, [10.0, 11.0])
2-element Vector{Any}:
 10.0
 11.0
```

Keyword argument `init` for transducible processes also accept an
`OnInit`:

```jldoctest OnInit
julia> foldl(push!, Map(identity), "abc"; init=OnInit(() -> []))
3-element Vector{Any}:
 'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)
 'b': ASCII/Unicode U+0062 (category Ll: Letter, lowercase)
 'c': ASCII/Unicode U+0063 (category Ll: Letter, lowercase)
```

To create a copy of a mutable object, [`CopyInit`](@ref) is easier to
use.

However, more powerful and generic pattern is to use `push!!` from
BangBang.jl and initialize `init` with `Union{}[]` so that it
automatically finds the minimal element type.

```jldoctest OnInit
julia> using BangBang

julia> foldl(push!!, Map(identity), "abc"; init=Union{}[])
3-element Vector{Char}:
 'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)
 'b': ASCII/Unicode U+0062 (category Ll: Letter, lowercase)
 'c': ASCII/Unicode U+0063 (category Ll: Letter, lowercase)
```
"""
struct OnInit{F} <: AbstractInitializer
    f::F
end

initvalue(init::OnInit) = init.f()

Base.show(io::IO, init::OnInit) = _default_show(io, init)

"""
    CopyInit(value)

This is equivalent to `OnInit(() -> deepcopy(value))`.

!!! compat "Transducers.jl 0.3"

    New in version 0.3.

# Examples
```jldoctest
julia> using Transducers

julia> init = CopyInit([]);

julia> foldl(push!, Map(identity), 1:3; init=init)
3-element Vector{Any}:
 1
 2
 3

julia> foldl(push!, Map(identity), 1:3; init=init)  # `init` can be reused
3-element Vector{Any}:
 1
 2
 3
```
"""
struct CopyInit{T} <: AbstractInitializer
    value::T
end

initvalue(init::CopyInit) = deepcopy(init.value)

Base.show(io::IO, init::CopyInit) = _default_show(io, init)

@inline foldlargs(op, x) = x
@inline foldlargs(op, x1, x2, xs...) =
    foldlargs(op,
              @next(op, x1, x2),
              xs...)

@inline mapfoldlargs(f, op, xs...) = foldlargs(Map(f)'(op), xs...)

"""
    @default_finaltype(xf::Transducer, coll)

Infer the type of the object that would be fed into the second
argument `input` of the bottom reducing function `rf(acc, input)`.

See: `Base.@default_eltype`
"""
macro default_finaltype(xf, coll)
    quote
        # `eltype(coll)` does not work when `coll` is an iterator;
        # that's why we are using `Base.@default_eltype`:
        intype = Base.@default_eltype($(esc(coll)))
        Core.Compiler.return_type(
            _getoutput,
            Tuple{typeof($(esc(xf))), intype},
        ) |> _real_state_type
    end
end

# # How `@default_finaltype` works
#
# When `xf` has a filter-like transducer, the inference would say that
# the output of `next(reducingfunction(xf, right), ..., x)` is
# `Union{Nothing, THE_DESIRED_TYPE}` because `nothing` is the default
# "identity" of `right`.  So, we set the initial value of `acc` as a
# dummy singleton `_FakeState()` and strip it off from the inferred
# `Union{_FakeState, THE_DESIRED_TYPE}` to get `THE_DESIRED_TYPE`.
#
# Note that using `_FakeState` is important for supporting cases like
# `Union{Nothing, Int}` being the correct output type of `xf`.

struct _FakeState end

function _getoutput(xf, x)
    rf = reducingfunction(xf, right)
    return unreduced(complete(rf, next(rf, start(rf, _FakeState()), x)))
end

_real_state_type(T) = T
_real_state_type(::Type{Union{T, _FakeState}}) where {T} = @isdefined(T) ? T : Any

"""
    Init
    Init(op) :: InitialValues.InitialValue

The canonical initializer; i.e., a singleton placeholder usable for
`init` argument of `foldl` for binary functions with known initial
values.

When `init = Init` is passed to `foldl` etc., `Init(op)` is called for
the bottom reducing function `op` during the [`start`](@ref) phase.
`Init(op)` returns
[`InitialValue(op)`](https://juliafolds.github.io/InitialValues.jl/dev/#InitialValues.InitialValue)
which acts as the canonical initial value of `op`.

# Examples
```jldoctest Init
julia> using Transducers

julia> foldl(+, 1:3 |> Filter(isodd); init = Init)
4

julia> foldl(+, 2:2:4 |> Filter(isodd); init = Init)
InitialValue(+)
```

# Extended help

Note that `op` passed to `foldl` etc. must be known to
InitialValues.jl:

```jldoctest Init
julia> unknown_op(a, b) = a + b;

julia> foldl(unknown_op, 2:2:4 |> Filter(isodd); init = Init)
ERROR: IdentityNotDefinedError: `init = Transducers.Init` is specified but the identity element `InitialValue(op)` is not defined for
    op = unknown_op
[...]
```

`InitialValues.asmonoid` can be used to wrap a binary function to add
("adjoin") the identity value to its domain:

```jldoctest Init
julia> using InitialValues: asmonoid

julia> foldl(asmonoid(unknown_op), 2:2:4 |> Filter(isodd); init = Init)
InitialValue(::InitialValues.AdjoinIdentity{typeof(unknown_op)})
```

When [`start(rf, Init)`](@ref start) is called with a composite
reducing function `rf`, `Init(rf₀)` is called for the bottom reducing
function `rf₀` of `rf`:

```jldoctest Init
julia> rf = Take(3)'(+);  # `+` is the bottom reducing function

julia> acc = Transducers.start(rf, Init);

julia> Transducers.unwrap(rf, acc)
(3, InitialValue(+))
```
"""
Init
@def_singleton Init = InitOf{InitialValues.InitialValueOf}()

"""
    DefaultInit(op)

`DefaultInit` is like [`Init`](@ref) but **strictly** internal
to Transducers.jl.  It is used for checking if the bottom reducing
function is never called.
"""
DefaultInit
struct DefaultInitOf{OP} <: SpecificInitialValue{OP} end
@def_singleton DefaultInit = InitOf{DefaultInitOf}()

struct OptInitOf{OP} <: SpecificInitialValue{OP} end
@def_singleton OptInit = InitOf{OptInitOf}()
# It seems that compiler can infer more when passing around a
# `Function` than a `Type` (since a `Function` is a singleton?).
# That's why `OptInit` is defined as a factory.

InferableInit{OP} = Union{DefaultInitOf{OP}, OptInitOf{OP}}

_nonidtype(::Any) = nothing
_nonidtype(::Type{Union{S, T}}) where {T, S <: InferableInit} = T

# Defining `_asmonoid` internally as the resulting reducing function
# is not really a monoid (only the bottom reducing function becomes a
# monoid).
# TODO: better name for `_asmonoid`
@inline _asmonoid(rf) = asmonoid(rf)
@inline _asmonoid(rf::Reduction) = Reduction(xform(rf), _asmonoid(inner(rf)))
@inline _asmonoid(rf::BottomRF) = BottomRF(_asmonoid(inner(rf)))
@inline _asmonoid(rf::Completing) = Completing(_asmonoid(rf.f))

struct EmptyResultError <: Exception
    rf
end

function Base.showerror(io::IO, e::EmptyResultError)
    println(
        io,
        "EmptyResultError: ",
        "Reducing function `", _realbottomrf(e.rf), "` is never called. ")
    print(
        io,
        "The input collection is empty or the items are all filtered out ",
        "by some transducer(s). ",
        "It is recommended to specify `init` to avoid this kind of errors.")
    # TODO: improve error message
end

_realbottomrf(op) = op
_realbottomrf(rf::AbstractReduction) = _realbottomrf(as(rf, BottomRF).inner)
_realbottomrf(rf::Completing) = rf.f

check_empty(rf, result) =
    if unreduced(result) isa DefaultInitOf
        throw(EmptyResultError(rf))
    else
        result
    end

struct IdentityNotDefinedError <: Exception
    op
    idfactory
end

function Base.showerror(io::IO, e::IdentityNotDefinedError)
    Init = e.idfactory
    op = e.op
    print(io, "IdentityNotDefinedError: ")
    print(io, strip("""
    `init = $Init` is specified but the identity element `InitialValue(op)` is not defined for
        op = $op
    Note that `op` must be a well known binary operations like `+` or `*`.
    See InitialValues.jl documentation for more information.
    """))
end


text:
# Setup:

const REPO = "https://github.com/ablaom/MachineLearningInJulia2020"
using Pkg
uuid = Pkg.TOML.parsefile("Project.toml")["uuid"]
uuid == "4764ce03-6504-4302-ab9a-b32cdba420f3" ||
    error("It appears this file is not in the same directory as other "*
          "files it needs (and in particular, the intended "*
          "Project.toml file). "*
          "A complete tutorial can be obtained by cloning the"*
          "MachineLearningInJulia repository from $REPO. ")
Pkg.activate(DIR)
Pkg.instantiate()
using CategoricalArrays
import MLJLinearModels
import DataFrames
import CSV
import DecisionTree
import MLJFlux
import Plots
import MLJLinearModels
import MultivariateStats
using MLJ

ext: jl


text:
# See https://gist.github.com/singularitti/e9e04c501ddfe40ba58917a754707b2e
const INTEGER = raw"([-+]?[0-9]+)"
const FIXED_POINT_REAL = raw"([-+]?[0-9]*\.[0-9]+|[0-9]+\.?[0-9]*)"
const GENERAL_REAL = raw"([-+]?(?:[0-9]*\.[0-9]+|[0-9]+\.?[0-9]*)(?:[eE][-+]?[0-9]+)?)"
const EQUAL_SIGN = raw"\s*=\s*"

# This format is from https://github.com/QEF/q-e/blob/4132a64/Modules/environment.f90#L215-L224.
const PARALLEL_INFO = r"(?<kind>(?:Parallel version [^,]*|Serial version))(?:, running on\s*(?<num>[0-9]+) processors)?"
const READING_INPUT_FROM = r"(?:Reading input from \s*(.*|standard input))"
const PWSCF_VERSION = r"Program PWSCF v\.(?<version>[0-9]\.[0-9]+\.?[0-9]?)"
# This format is from https://github.com/QEF/q-e/blob/4132a64/PW/src/summary.f90#L374-L375.
const FFT_DIMENSIONS = Regex("Dense  grid:\\s*$INTEGER\\s+G-vectors\\s+FFT dimensions: \\(\\s*$INTEGER,\\s*$INTEGER,\\s*$INTEGER\\)")
# The following format is from https://github.com/QEF/q-e/blob/7357cdb/PW/src/summary.f90#L100-L119.
const SUMMARY_BLOCK = r"(bravais-lattice index\X+?)\s*celldm"  # Match between "bravais-lattice index" & the 1st of the "celldm"s, `+?` means un-greedy matching (required)
# 'bravais-lattice index     = ',I12
const BRAVAIS_LATTICE_INDEX = Regex("bravais-lattice index$EQUAL_SIGN$INTEGER")
# 'lattice parameter (alat)  = ',F12.4,'  a.u.'
const LATTICE_PARAMETER = Regex("lattice parameter \\(alat\\)$EQUAL_SIGN$FIXED_POINT_REAL")
# 'unit-cell volume          = ',F12.4,' (a.u.)^3'
const UNIT_CELL_VOLUME = Regex("unit-cell volume$EQUAL_SIGN$FIXED_POINT_REAL")
# 'number of atoms/cell      = ',I12
const NUMBER_OF_ATOMS_PER_CELL = Regex("number of atoms\\/cell$EQUAL_SIGN$INTEGER")
# 'number of atomic types    = ',I12
const NUMBER_OF_ATOMIC_TYPES = Regex("number of atomic types$EQUAL_SIGN$INTEGER")
# 'number of electrons       = ',F12.2,' (up:',f7.2,', down:',f7.2,')'
const NUMBER_OF_ELECTRONS = Regex("number of electrons$EQUAL_SIGN$FIXED_POINT_REAL" *
                                  "(?:\\(up:\\s*$FIXED_POINT_REAL, down:\\s*$FIXED_POINT_REAL\\))?")
# 'number of Kohn-Sham states= ',I12
const NUMBER_OF_KOHN_SHAM_STATES = Regex("number of Kohn-Sham states$EQUAL_SIGN$INTEGER")
# 'kinetic-energy cutoff     = ',F12.4,'  Ry'
const KINETIC_ENERGY_CUTOFF = Regex("kinetic-energy cutoff$EQUAL_SIGN$FIXED_POINT_REAL\\s+Ry")
# 'charge density cutoff     = ',F12.4,'  Ry'
const CHARGE_DENSITY_CUTOFF = Regex("charge density cutoff$EQUAL_SIGN$FIXED_POINT_REAL\\s+Ry")
# 'cutoff for Fock operator  = ',F12.4,'  Ry'
const CUTOFF_FOR_FOCK_OPERATOR = Regex("cutoff for Fock operator$EQUAL_SIGN$FIXED_POINT_REAL\\s+Ry")
# 'convergence threshold     = ',1PE12.1
const CONVERGENCE_THRESHOLD = Regex("convergence threshold$EQUAL_SIGN$GENERAL_REAL")
# 'mixing beta               = ',0PF12.4
const MIXING_BETA = Regex("mixing beta$EQUAL_SIGN$FIXED_POINT_REAL")
# 'number of iterations used = ',I12,2X,A,' mixing'
const NUMBER_OF_ITERATIONS_USED = Regex("number of iterations used$EQUAL_SIGN$INTEGER\\s+([-+\\w]+)\\s+mixing")
const EXCHANGE_CORRELATION = r"Exchange-correlation\s*=\s*(.*)"
# "nstep                     = ",I12
const NSTEP = Regex("nstep$EQUAL_SIGN$INTEGER")
# The following format is from https://github.com/QEF/q-e/blob/4132a64/Modules/fft_base.f90#L70-L91.
const FFT_BASE_INFO = r"""\s*(?<head>Parallelization info|G-vector sticks info)
\s*--------------------
\s*sticks:   dense  smooth     PW     G-vecs:    dense   smooth      PW
(?<body>(?:\s*Min.*)?(?:\s*Max.*)?\s*Sum.*)"""m
const SYM_OPS = r"\h*(No symmetry found|(?<n>\d+)\s*Sym\. Ops\..*found)"
# The following format is from https://github.com/QEF/q-e/blob/4132a64/PW/src/summary.f90#L341-L381.
const K_POINTS_BLOCK = r"""
number of k points=\s*(?<nk>[0-9]+)\h*(?<metainfo>.*)
\s*(?:cart\. coord\. in units 2pi\/alat\s*(?<cart>\X+?)^\s*$|Number of k-points >= 100: set verbosity='high' to print them\.)
\s*(?:cryst\. coord\.\s*(?<cryst>\X+?)
\s*Dense  grid)?"""m
# The following format is from https://github.com/QEF/q-e/blob/4132a64/PW/src/summary.f90#L353-L354.
# '(8x,"k(",i5,") = (",3f12.7,"), wk =",f12.7)'
const K_POINTS_ITEM = Regex("k\\(.*\\) = \\(\\s*$FIXED_POINT_REAL\\s*$FIXED_POINT_REAL\\s*$FIXED_POINT_REAL\\s*\\), wk =\\s*$FIXED_POINT_REAL")
# The following format is from https://github.com/QEF/q-e/blob/4132a64/PW/src/output_tau.f90#L47-L60.
const CELL_PARAMETERS_BLOCK = r"""
CELL_PARAMETERS \h+
\( (?<option>\w+) =? \s* (?<alat>[-+]?[0-9]*\.[0-9]{8})? \) \h*  # Match `alat`: `F12.8`
(?<data>
    (?: \s*
        (?:
            [-+]?[0-9]*\.[0-9]+ \s*  # Match element
        ){3}  # I need exactly 3 elements per vector
    ){3}  # I need exactly 3 vectors
)
"""x
const CELL_PARAMETERS_ITEM = r"""
\s*
([-+]?[0-9]*\.[0-9]+) \s*  # x
([-+]?[0-9]*\.[0-9]+) \s*  # y
([-+]?[0-9]*\.[0-9]+) \s*  # z
"""x
# The following format is from https://github.com/QEF/q-e/blob/4132a64/PW/src/output_tau.f90#L64-L109.
const ATOMIC_POSITIONS_BLOCK = r"""
ATOMIC_POSITIONS \h*                   # Atomic positions start with that string
\( (?<option>\w+) \)                   # Option of the card
(?<data>
    (?:
        \s*
        [A-Za-z]+[A-Za-z0-9]{0,2} \s+  # Atom spec
        (?:
            [-+]?[0-9]*\.[0-9]+ \s*  # Match element
        ){3}                           # I need exactly 3 floats per vector.
        (?:
            [-+]?[0-9]+ \s*
        ){0,3}                     # I need exactly 3 integers in `if_pos`, if there is any.
    )+
)
"""x
const ATOMIC_POSITIONS_ITEM = r"""
\s*
([A-Za-z]+[A-Za-z0-9]{0,2}) \s+  # Atom spec
([-+]?[0-9]*\.[0-9]+) \s*  # x
([-+]?[0-9]*\.[0-9]+) \s*  # y
([-+]?[0-9]*\.[0-9]+) \s*  # z
([-+]?[0-9]+)? \s*            # if_pos(1)
([-+]?[0-9]+)? \s*            # if_pos(2)
([-+]?[0-9]+)? \s*            # if_pos(3)
"""x
const FINAL_COORDINATES_BLOCK = r"""
Begin final coordinates
(\X+?)
End final coordinates
"""
const STRESS_BLOCK = r"""
^[ \t]*
total\s+stress\s*\(Ry\/bohr\*\*3\)\s+
\(kbar\)\s+P=\s*([\-|\+]? (?: [0-9]*[\.][0-9]+ | [0-9]+[\.]?[0-9]*)
    ([E|e|d|D][+|-]?[0-9]+)?)
\R
(
(?:
\s*
(?:
[\-|\+]? (?: [0-9]*[\.][0-9]+ | [0-9]+[\.]?[0-9]*)
    ([E|e|d|D][+|-]?[0-9]+)?[ \t]*
){6}
){3}
)
"""mx
const SELF_CONSISTENT_CALCULATION_BLOCK = r"(Self-consistent Calculation\X+?End of self-consistent calculation)"
const ITERATION_BLOCK = r"(?<=iteration #)(.*?)(?=iteration #|End of self-consistent calculation)"s
# This format is from https://github.com/QEF/q-e/blob/4132a64/PW/src/electrons.f90#L920-L921.
# '     iteration #',I3,'     ecut=', F9.2,' Ry',5X,'beta=',F5.2
const ITERATION_HEAD = Regex("\\s*$INTEGER\\s+ecut=\\s*$FIXED_POINT_REAL\\s+Ry\\s+beta=\\s*$FIXED_POINT_REAL")
# These formats are from https://github.com/QEF/q-e/blob/4132a64/PW/src/c_bands.f90#L129-L130
# and https://github.com/QEF/q-e/blob/4132a64/PW/src/c_bands.f90#L65-L73.
const C_BANDS = Regex(
    """
    (?<diag>Davidson diagonalization.*|CG style diagonalization|PPCG style diagonalization)
    \\h*ethr =\\h*$GENERAL_REAL,  avg # of iterations =\\h*$FIXED_POINT_REAL""",
    "m",
)
# This format is from https://github.com/QEF/q-e/blob/4132a64/PW/src/electrons.f90#L917-L918.
# '     total cpu time spent up to now is ',F10.1,' secs'
const TOTAL_CPU_TIME = Regex("total cpu time spent up to now is\\s*$FIXED_POINT_REAL\\s* secs")
const KS_ENERGIES_BLOCK = r"""
(Number\s+of\s+k-points\s+>=\s+100:\s+set\s+verbosity='high'\s+to\s+print\s+the\s+bands\.
|
\s*(------\s+SPIN\s+UP\s+------------|------\s+SPIN\s+DOWN\s+----------)?
(?:
\s*(?:k\s+=(?:\s*[-+]?[0-9]*\.[0-9]+|[0-9]+\.?[0-9]*){3}\s*\(\s*(?:[-+]?[0-9]+)\s*PWs\)\s+bands\s+\(ev\):
|
k\s+=(?:\s*[-+]?[0-9]*\.[0-9]+|[0-9]+\.?[0-9]*){3}\s+band\s+energies\s+\(ev\):
)
(?:(?:\s*[-+]?[0-9]*\.[0-9]+|[0-9]+\.?[0-9]*){1,8})+
)+
)"""mx
const KS_ENERGIES_BANDS = r"""
\s*k\s+=\s*(?<k>.*)\s*\(\s*(?<PWs>[-+]?[0-9]+)\s*PWs\)\s+bands\s+\(ev\):
(?<band>(?:\s*[-+]?(?:[0-9]*\.[0-9]+|[0-9]+\.?[0-9]*)(?:[eE][-+]?[0-9]+)?)+)"""
const KS_ENERGIES_BAND_ENERGIES = r"""
\s*k\s+=\s*(?<k>.*)\s*band\s+energies\s+\(ev\):
(?<band>(?:\s*[-+]?(?:[0-9]*\.[0-9]+|[0-9]+\.?[0-9]*)(?:[eE][-+]?[0-9]+)?)+)"""
# This format is from https://github.com/QEF/q-e/blob/4132a64/PW/src/electrons.f90#L1257-L1261.
const UNCONVERGED_ELECTRONS_ENERGY = Regex(
    """
    ^[^!]\\s+total energy\\s+=\\s*$FIXED_POINT_REAL\\s+Ry(?:\\s*Harris-Foulkes estimate\\s+=\\s*$FIXED_POINT_REAL\\s+Ry)?
    \\s*estimated scf accuracy\\s+<\\s*$GENERAL_REAL\\s+Ry""",
    "m",
)
const CONVERGED_ELECTRONS_ENERGY = Regex(
    """
    ^!\\h+total energy\\s+=\\s*$FIXED_POINT_REAL\\s+Ry(?:\\s*Harris-Foulkes estimate\\s+=\\s*$FIXED_POINT_REAL\\s+Ry)?(?:\\s*total all-electron energy\\s+=\\s*$FIXED_POINT_REAL\\s+Ry)?
    \\s*estimated scf accuracy\\s+<\\s*$GENERAL_REAL\\s+Ry
    \\s*(?<ae>total all-electron energy =.*Ry)?\\s*(?<decomp>The total energy is the sum of the following terms:
    \\s*one-electron contribution =.*Ry
    \\s*hartree contribution      =.*Ry
    \\s*xc contribution           =.*Ry
    \\s* ewald contribution        =.*Ry)?\\s*(?<one>one-center paw contrib.*Ry
    \\s*-> PAW hartree energy AE =.*Ry
    \\s*-> PAW hartree energy PS =.*Ry
    \\s*-> PAW xc energy AE      =.*Ry
    \\s*-> PAW xc energy PS      =.*Ry
    \\s*-> total E_H with PAW    =.*Ry
    \\s*-> total E_XC with PAW   =.*Ry)?\\s*(?<smearing>smearing contrib.*Ry)?""",
    "m",
)
const TIME_BLOCK = r"(init_run\X+?This run was terminated on:.*)"
# This format is from https://github.com/QEF/q-e/blob/4132a64/PW/src/print_clock_pw.f90#L29-L33.
const SUMMARY_TIME_BLOCK = r"""
(?<head>)
(?<body>
init_run\s+:.*
\s*electrons\s+:.*
\s*(?:update_pot\s+.*)?  # This does not always exist.
\s*(?:forces\s+:.*)?     # This does not always exist.
\s*(?:stress\s+:.*)?     # This does not always exist.
)
"""mx
const TIME_ITEM = Regex("\\s*([\\w0-9:]+)\\s+:\\s*$(FIXED_POINT_REAL)s\\sCPU\\s*$(FIXED_POINT_REAL)s\\sWALL\\s\\(\\s*$INTEGER\\scalls\\)")
# This format is from https://github.com/QEF/q-e/blob/4132a64/PW/src/print_clock_pw.f90#L35-L36.
const INIT_RUN_TIME_BLOCK = r"Called by (?<head>init_run):(?<body>\X+?)^\s*$"m
# This format is from https://github.com/QEF/q-e/blob/4132a64/PW/src/print_clock_pw.f90#L53-L54.
const ELECTRONS_TIME_BLOCK = r"Called by (?<head>electrons):(?<body>\X+?)^\s*$"m
# This format is from https://github.com/QEF/q-e/blob/4132a64/PW/src/print_clock_pw.f90#L78-L79.
const C_BANDS_TIME_BLOCK = r"Called by (?<head>c_bands):(?<body>\X+?)^\s*$"m
const SUM_BAND_TIME_BLOCK = r"Called by (?<head>sum_band):(?<body>\X+?)^\s*$"m
const EGTERG_TIME_BLOCK = r"Called by (?<head>\*egterg):(?<body>\X+?)^\s*$"m
const H_PSI_TIME_BLOCK = r"Called by (?<head>h_psi):(?<body>\X+?)^\s*$"m
const GENERAL_ROUTINES_TIME_BLOCK = r"(?<head>General routines)(?<body>\X+?)^\s*$"m
const PARALLEL_ROUTINES_TIME_BLOCK = r"(?<head>Parallel routines)(?<body>\X+?)^\s*$"m
const TERMINATED_DATE = r"This run was terminated on:(.+)"  # TODO: Date
const JOB_DONE = r"JOB DONE\."
# These formats are from https://github.com/QEF/q-e/blob/4132a64/UtilXlib/error_handler.f90#L48-L68.
const ERROR_BLOCK = r"%{78}(?<body>\X+?)\s*%{78}"
const ERROR_IN_ROUTINE = r"Error in routine\s+(.*)\s+\((.*)\):"


text:
function triangularity(g::GEQDSKFile)
    Rmin, Rmax = extrema(filter(!iszero,g.rbbbs))
    Rgeo = (Rmin + Rmax)/2
    a = (Rmax - Rmin)/2
    Rupper = g.rbbbs[argmax(g.zbbbs)]
    Rlower = g.rbbbs[argmin(g.zbbbs)]
    delta_upper = (Rgeo - Rupper)/a
    delta_lower = (Rgeo - Rlower)/a

    return delta_lower, delta_upper
end

function ellipticity(g::GEQDSKFile)
    Rmin, Rmax = extrema(filter(!iszero,g.rbbbs))
    a = (Rmax - Rmin)/2
    Zmin, Zmax = extrema(g.zbbbs)

    return (Zmax - Zmin)/(2a)
end

function elongation(g::GEQDSKFile)
    ellipticity(g)
end

function major_radius(g::GEQDSKFile)
    0.5*(+(extrema(filter(!iszero,g.rbbbs))...))
end

function minor_radius(g::GEQDSKFile)
    -0.5*(-(extrema(filter(!iszero,g.rbbbs))...))
end

function aspect_ratio(g::GEQDSKFile)
    major_radius(g)/minor_radius(g)
end


text:
module DrawData

include("./dataset.jl")

export Dataset

end # module DrawData


text:
module RelayGLM

export Poisson, Binomial, Gaussian, Identity, Exponential, Logistic
export NullPrior, RidgePrior, SmoothingPrior, LassoPrior
export PredictorSet, Predictor, GLM, DefaultBasis, CosineBasis
export wasrelayed, set_parameter!, get_basis, cross_validate, get_predictor, get_basis
export get_coef, get_error, get_lambda, metric, ncol
export shuffle, JSDivergence, ROCArea, PRArea, BinomialLikelihood, RelativeBinomialLikelihood, RRI

include("./RelayUtils.jl")
using .RelayUtils

include("./GLMMetrics.jl")
using .GLMMetrics

include("./interface.jl")

import .GLMFit.GLMModels

include("./RelayISI.jl")

# ============================================================================ #
function predict(::Type{A}, ret::Vector{Float64}, xf::Vector{<:Real}, b::Real) where {A <: GLMModels.ActivationFunction}
    idx_max = round(Int, ret[end] / 0.001)
    retidx = ts2idx.(ret, 0.001)
    r, _ = GLMModels.activate(A,
        vec(
            spikefilter(retidx, reverse(xf, dims=1), (1, idx_max))
        ) .+ b,
        0.001
    )
    return r
end
# ---------------------------------------------------------------------------- #
function predict(::Type{A}, ret::Vector{Float64}, lgn::Vector{Float64}, xf::Vector{<:Real}, hf::Vector{<:Real}, b::Real) where {A <: GLMModels.ActivationFunction}
    idx_max = round(Int, max(ret[end], lgn[end]) / 0.001)

    retidx = ts2idx.(ret, 0.001)
    lgnidx = ts2idx.(lgn, 0.001)

    r, _ = GLMModels.activate(A,
        vec(
            spikefilter(retidx, reverse(xf, dims=1), (1, idx_max)) .+
            spikefilter(lgnidx, reverse(hf, dims=1), (1, idx_max))
        ) .+ b,
        0.001
    )
    return r
end
# ============================================================================ #
end


text:
"Get Beta Distribution parameters using standard tabulated values θₗ,s, Eq. 2.9 and 2.10 in Bonan"
function βparameters(θₗ,s)
    α = (1 - (s^2+θₗ^2) / (θₗ*π/2) ) / ( ((s^2+θₗ^2) / θₗ^2) - 1 )
    β = ( (π/2)/θₗ -1 ) * α
    return α,β
end

"Eq 46 in Shultis and Myneni"
function H(μ::FT,μₗ::FT) where FT
    # cot(θ) * cot(θₗ)
    r = cot(acos(μ)) * cot(acos(μₗ));
    if r > 1
        return μ * μₗ
    elseif r < -1
        return FT(0)
    else
        # Eq 47
        ϕ = acos(-r);
        return FT(1/π) * (μ*μₗ*ϕ + sqrt(FT(1)-μ^2) * sqrt(1-μₗ^2) * sin(ϕ))
    end
end

"Eq 45 in Shultis and Myneni, fixed grid in μ for both μ and μ' "
function compute_Ψ(μ::Array{FT,1}, μₗ::FT) where FT
    H⁺ = H.(μ,μₗ)
    H⁻ = H.(-μ,μₗ)
    Ψ⁺ = H⁺ * H⁺' + H⁻ * H⁻'
    Ψ⁻ = H⁺ * H⁻' + H⁻ * H⁺'
    return Ψ⁺, Ψ⁻
end

"Eq 45 in Shultis and Myneni, fixed grid in μ for both μ and μ' "
function compute_Ψ(μ::Array{FT,1},μꜛ::Array{FT,1}, μₗ::FT) where FT
    Ψ⁺ = H.(μ,μₗ) * H.(μꜛ,μₗ)'  + H.(-μ,μₗ) * H.(-μꜛ,μₗ)'
    Ψ⁻ = H.(μ,μₗ) * H.(-μꜛ,μₗ)' + H.(-μ,μₗ) * H.(+μꜛ,μₗ)'
    return Ψ⁺, Ψ⁻
end

" Eq 37 in Shultis and Myneni "
function scattering_model_lambertian(Ωⁱⁿ::dirVector, Ωᵒᵘᵗ::dirVector, Ωᴸ::dirVector, r, t)
    α = (Ωⁱⁿ ⋅ Ωᴸ) * (Ωᵒᵘᵗ ⋅ Ωᴸ)
    if a<0
        return r * abs(Ωᵒᵘᵗ ⋅ Ωᴸ) / π
    else
        return t * abs(Ωᵒᵘᵗ ⋅ Ωᴸ) / π
    end
end

#See Eq 18, Canopy RT book
function getSpecularΩ(Ωⁱⁿ::dirVector{FT}, Ωᵒᵘᵗ::dirVector{FT}) where FT
    sa = Ωⁱⁿ ⋅ Ωᵒᵘᵗ
    sa > 1 ? sa = FT(1) : nothing

    α = acos(sa)/2
    # Relative azimuth angle:
    ϕ = (Ωᵒᵘᵗ.ϕ - Ωⁱⁿ.ϕ)
    # Leaf polar angle:
    θₗ = acos((cos(Ωⁱⁿ.θ) + cos(Ωᵒᵘᵗ.θ))/2cos(α));

    t = (sin(Ωⁱⁿ.θ) + sin(Ωᵒᵘᵗ.θ) * cos(ϕ))/(2cos(α)*sin(θₗ))
    if isnan(t)
        t = FT(1.0)
    end
    # Leaf azimuth angle:
    ϕₗ = acos(max(FT(-1),min(t,FT(1))))
    ϕ > π   ? ϕₗ = 2π-ϕₗ-Ωⁱⁿ.ϕ : ϕₗ = ϕₗ+Ωⁱⁿ.ϕ
    #θₗ > π/2 ? θₗ = θₗ - π/2  : nothing
    out = dirVector(θₗ,ϕₗ)
    #ϕ > π   ? out = dirVector(θₗ,2π-ϕₗ-Ωⁱⁿ.ϕ) : out = dirVector(θₗ,ϕₗ+Ωⁱⁿ.ϕ)
    #@show (Ωⁱⁿ ⋅ out) , (Ωᵒᵘᵗ ⋅ out)
    #@show abs(((Ωⁱⁿ ⋅ out) - (Ωᵒᵘᵗ ⋅ out)) / (Ωᵒᵘᵗ ⋅ out)) * 100
    #@show θₗ , ϕₗ
    @assert abs(((Ωⁱⁿ ⋅ out) - (Ωᵒᵘᵗ ⋅ out)) / (Ωᵒᵘᵗ ⋅ out)) * 100 < 1 "Leaf angle wrong, $Ωⁱⁿ, $Ωᵒᵘᵗ, $out"
    return out

end

function getSpecularΩ(Ωⁱⁿ::dirVector_μ{FT}, Ωᵒᵘᵗ::dirVector_μ{FT}) where FT
    sa = Ωⁱⁿ ⋅ Ωᵒᵘᵗ
    #@show sa
    sa > 1 ? sa = FT(1) : nothing
    # Scattering angle
    α = acos(sa)/FT(2)

    # Relative azimuth angle:
    ϕ = (Ωᵒᵘᵗ.ϕ - Ωⁱⁿ.ϕ)
    # Leaf polar angle:
    μₗ = ((Ωⁱⁿ.μ + Ωᵒᵘᵗ.μ)/2cos(α));

    t = (sqrt(1-Ωⁱⁿ.μ^2) + sqrt(1-Ωᵒᵘᵗ.μ^2) * cos(ϕ))/(2cos(α)*sqrt(1-μₗ^2))
    #@show t, sqrt(1-Ωⁱⁿ.μ^2) + sqrt(1-Ωᵒᵘᵗ.μ)
    if isnan(t)
        t = FT(1.0)
    end
    # Leaf azimuth angle:
    ϕₗ = acos(max(FT(-1),min(t,FT(1))))
    ϕ > π   ? ϕₗ = 2π-ϕₗ-Ωⁱⁿ.ϕ : ϕₗ = ϕₗ+Ωⁱⁿ.ϕ
    out = dirVector_μ(μₗ,ϕₗ)
    @assert Ωⁱⁿ ⋅ out ≈  Ωᵒᵘᵗ ⋅ out "Leaf angle wrong in specular reflection, $Ωⁱⁿ, $Ωᵒᵘᵗ, $out"
    return out, α
end


text:
"""
    calculate_measure!(measure::ConditionalExpectedReturn{Inf,N,S,R}, w)

Returns worst case return in the convex hull uncertainty set, defined by the following dual problem:

```math
    \\max_{\\theta} \\quad  \\theta \\\\
    s.t.  \\quad \\theta \\leq r_s ' w \\quad \\forall s = 1:\\mathcal{S} \\\\
```

Further information:
  - Fernandes, B., Street, A., ValladA˜ £o, D., e Fernandes, C. (2016). An adaptive robust portfolio optimization model with loss constraints based on data-driven polyhedral uncertainty sets. European Journal of Operational Research, 255(3):961 – 970. ISSN 0377-2217. URL.

"""
function calculate_measure!(measure::ConditionalExpectedReturn{Inf,N,D,R}, w) where {D<:ContinuousMultivariateSampleable,N,R}
    model = owner_model(w)
    s = ambiguityset(measure)
    samples = rand(s, sample_size(measure))

    # auxilary variables
    θ = @variable(model)
    # convex hull
    @constraint(model, sum(samples' * w , dims=2) .>= θ)

    return θ
end

"""conditional_expectation = -cvar = -expected_shortfall"""
function calculate_measure!(measure::ConditionalExpectedReturn{α,N,D,EstimatedCase}, w) where {D<:ContinuousMultivariateSampleable,α,N}
    model = owner_model(w)
    s = ambiguityset(measure)

    # parameters
    numS = sample_size(measure)
    samples = rand(s, numS)
    sample_probability = fill(1/numS,numS)

    # dual variables
    @variable(model, z)
    @variable(model, y[i=1:numS] >= 0)

    @constraints(
        model,
        begin
            ys[s=1:numS], y[s] >= z - sum(dot(samples[:, s], w))
        end
    )
    return z - dot(sample_probability,y) / (1.0 - α)
end

max_stars_repo_name: andrewrosemberg/PortfolioOpt

max_issues_repo_name: andrewrosemberg/PortfolioOpt

max_forks_repo_name: andrewrosemberg/PortfolioOpt


text:
# TODO: NEED VARIABLE SUBTYPES...

is_convex(::typeof(VARIABLE), x::Interval{T}) where T = true
is_concave(::typeof(VARIABLE), x::Interval{T}) where T = true

is_loglogconvex(::typeof(VARIABLE), x::Interval{T}) where T = true
is_loglogconcave(::typeof(VARIABLE), x::Interval{T}) where T = true

is_increasing(::typeof(VARIABLE), x::Interval{T}) where T = true
is_decreasing(::typeof(VARIABLE), x::Interval{T}) where T = false

is_positive(::typeof(VARIABLE), x::Interval{T}) where T = x >= zero(T)
is_negative(::typeof(VARIABLE), x::Interval{T}) where T = x <= zero(T)
is_locked(::typeof(VARIABLE), x::Interval{T}) where T = one(T) ∉ x

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#=
const CR = JuMP.ConstraintRef

function classify_sip!(mSIP::SIPModel, c::AbstractConstraint)
    error("Function classify_sip! not defined for a constraint of type = $(typeof(c))")
end

function update_status(status::SIPCons, term_status::SIPCons)
    if status == SIPNOTSET
        return term_status
    elseif (status == DECISION) && (term_status == DECISION)
        return DECISION
    elseif (status == UNCERTAIN) && (term_status == UNCERTAIN)
        return UNCERTAIN
    end
    return SEMIINFINITE
end

classify_sip(mSIP::SIPModel, v::VariableRef) = mSIP.p[v] ? UNCERTAIN : DECISION
function classify_sip(mSIP::SIPModel, expr::GenericAffExpr{T,VariableRef}) where T
    status = SIPNOTSET
    for term in linear_terms(expr)
        status = update_status(status, classify_sip(mSIP, term[1]))
        if status == SEMIINFINITE
            break
        end
    end
    return status
end
function classify_sip(mSIP::SIPModel, expr::GenericQuadExpr{T,VariableRef}) where T
    status = SIPNOTSET
    for term in quad_terms(expr)
        status = update_status(status, classify_sip(mSIP, term[2]))
        status = update_status(status, classify_sip(mSIP, term[3]))
        if status == SEMIINFINITE
            break
        end
    end
    for term in linear_terms(expr)
        status = update_status(status, classify_sip(mSIP, term[1]))
        if status == SEMIINFINITE
            break
        end
    end
    return status
end

for typ in (:(Vector{GenericQuadExpr{T,VariableRef}}),
            :(Vector{GenericAffExpr{T,VariableRef}}),
            :(Vector{VariableRef}))
    @eval function classify_sip(mSIP::SIPModel, expr::Vector{$typ{T,VariableRef}}) where T
        status = SIPNOTSET
        for ex in expr
            status = update_status(status, classify_sip(mSIP, ex))
            if status == SEMIINFINITE
                break
            end
        end
        return status
    end
end

function classify_sip!(mSIP::SIPModel, cr::ScalarConstraint{F,S}) where {F<:JuMP.AbstractJuMPScalar, S<:JuMP.AbstractScalarSet}
    mSIP.constraint_type[cr] = classify_sip(mSIP, constraint_object(cr).func)
    nothing
end

function classify_sip!(mSIP::SIPModel, cr::CR{M,C,S}) where {M<:JuMP.AbstractModel, C ,S<:JuMP.AbstractShape}
    mSIP.constraint_type[cr] = classify_sip(constraint_object(cr))
    nothing
end

function classify_sip!(mSIP::SIPModel)
    constraint_types = list_of_constraint_types(mSIP.m)
    for (func_typ, set_typ) in constraint_types
        constraint_refs = all_constraints(mSIP.m, func_typ, set_typ)
        for cr in constraint_refs
            classify_sip!(mSIP, cr)
        end
    end
    # TODO: ADD CLASSIFICATION OF NONLINEAR EQUATIONS
    return nothing
end
=#


text:
# Below are ideas for potential types at the _high_ level DSL. These would be used to
# manipulate ModelingToolkit IR -> generate functions for consumption by the DiffEq backend

abstract type Geometry end
abstract type Point <: Geometry end
abstract type Sphere <: Geometry end
abstract type Cylinder <: Geometry end

# Dynamics types describe inherent non-linearities/dynamics
abstract type Dynamics end
abstract type HodgkinHuxley <: Dynamics end
abstract type Izhikevich <: Dynamics end
abstract type MorrisLecar <: Dynamics end

const SectionID = UInt64

struct Compartment{G <: Geometry, D <: Dynamics, S <: Synapse}
    id::CompartmentID
    dims
end

abstract type AbstractNeuron end

abstract type SingleCompartment <: AbstractNeuron end

struct MultiCompartment <: AbstractNeuron
    compartments::IdDict{CompartmentID, Compartment}
    morphology::SimpleGraph # e.g. from LightGraphs.jl or some acyclic graph type
    synapses::Vector{CompartmentID}
end

abstract type ArtificialNN <: AbstractNeuron end

chnl = @channel begin
   # name, describe, and register a custom channel model
end

dyn = @dynamics begin
    # name and describe a set of membrane properties/non-linear dynamics
    # potentially combining different custom channel models, etc
end

neuron = @neuron begin
   # compose channel models, dynamics, morphology and topology into the description of a cell
end

max_stars_repo_name: wsphillips/FutureNeuralSimulation

max_issues_repo_name: wsphillips/FutureNeuralSimulation

max_forks_repo_name: wsphillips/FutureNeuralSimulation


text:
# Possible syntax: container.[indices]
"""
getindices(container, indices)

Return an indexable container with indices `keys(indices)` and values `container[i]` for
`i ∈ values(indices)`. This generalizes scalar `getindex(container, index)` for multiple
indices, for dictionaries, tuples, named tuples, etc.

# Examples

```julia
julia> d = Dict(:a => "Alice", :b => "Bob", :c => "Charlie")
Dict{Symbol,String} with 3 entries:
  :a => "Alice"
  :b => "Bob"
  :c => "Charlie"

julia> getindices(d, [:a, :c])
2-element Array{String,1}:
 "Alice"
 "Charlie"

julia> getindices(d, (:a, :c))
("Alice", "Charlie")

julia> getindices(d, Dict("Wife" => :a, "Husband" => :c))
Dict{String,String} with 2 entries:
  "Wife"    => "Alice"
  "Husband" => "Charlie"
```
"""
function getindices(container, indices)
    # TODO figure out how to do @inbounds
    return map(i -> container[i], indices)
end
getindices(container, ::Colon) = getindices(container, keys(container))
getindices(t::Tuple, ::Colon) = t
getindices(nt::NamedTuple, ::Colon) = nt

function getindices(container, indices::AbstractDict)
    out = empty(indices, keytype(indices), _valtype(container))
    for (i, v) in indices
        out[i] = container[v]
    end
    return out
end
function getindices(d::AbstractDict, ::Colon)
    return copy(d)
end

# Make a 0-D array instead of scalar
Base.@propagate_inbounds function getindices(a::AbstractArray, is::Int...)
    out = similar(a, eltype(a), ())
    out[] = a[is...]
    return out
end

# Ambiguities
Base.@propagate_inbounds getindices(a::AbstractArray, is::Union{Int,Colon,AbstractArray}...) = a[is...]
getindices(a::AbstractVector, ::Colon) = a[:]


text:
module MapGen

include("SAILx/src/maze.jl")

function random_map(h, w)
    maze, available = generate_maze(h, w; numdel=2)
    navimap = generate_navi_map(maze, "RANDOMMAZE"; itemcountprobs=[0.0 0.0 0.05 0.05 0.1 0.1 0.1 0.1 0.1 0.2 0.2], iprob=0.6)
end

end

max_stars_repo_name: ozanarkancan/navi-app

max_issues_repo_name: ozanarkancan/navi-app

max_forks_repo_name: ozanarkancan/navi-app


text:
"""
    DefaultManifold <: AbstractManifold

This default manifold illustrates the main features of the interface and provides a skeleton
to build one's own manifold. It is a simplified/shortened variant of `Euclidean` from
`Manifolds.jl`.

This manifold further illustrates how to type your manifold points and tangent vectors. Note
that the interface does not require this, but it might be handy in debugging and educative
situations to verify correctness of involved variabes.
"""
struct DefaultManifold{T<:Tuple,𝔽} <: AbstractManifold{𝔽} end
function DefaultManifold(n::Vararg{Int,N}; field = ℝ) where {N}
    return DefaultManifold{Tuple{n...},field}()
end

distance(::DefaultManifold, p, q) = norm(p - q)

embed!(::DefaultManifold, q, p) = copyto!(q, p)

embed!(::DefaultManifold, Y, p, X) = copyto!(Y, X)

exp!(::DefaultManifold, q, p, X) = (q .= p .+ X)

function get_basis_orthonormal(::DefaultManifold, p, N::RealNumbers)
    return CachedBasis(
        DefaultOrthonormalBasis(N),
        [_euclidean_basis_vector(p, i) for i in eachindex(p)],
    )
end
function get_basis_orthogonal(::DefaultManifold, p, N::RealNumbers)
    return CachedBasis(
        DefaultOrthogonalBasis(N),
        [_euclidean_basis_vector(p, i) for i in eachindex(p)],
    )
end
function get_basis_default(::DefaultManifold, p, N::RealNumbers)
    return CachedBasis(
        DefaultBasis(N),
        [_euclidean_basis_vector(p, i) for i in eachindex(p)],
    )
end
function get_basis_diagonalizing(M::DefaultManifold, p, B::DiagonalizingOrthonormalBasis)
    vecs = get_vectors(M, p, get_basis(M, p, DefaultOrthonormalBasis()))
    eigenvalues = zeros(real(eltype(p)), manifold_dimension(M))
    return CachedBasis(B, DiagonalizingBasisData(B.frame_direction, eigenvalues, vecs))
end

# Complex manifold, real basis -> coefficients c are complesx -> reshape
# Real manifold, real basis -> reshape
function get_coordinates_orthonormal!(M::DefaultManifold, c, p, X, ::RealNumbers)
    return copyto!(c, reshape(X, number_of_coordinates(M, ℝ)))
end
function get_coordinates_diagonalizing!(
    M::DefaultManifold,
    c,
    p,
    X,
    ::DiagonalizingOrthonormalBasis{ℝ},
)
    return copyto!(c, reshape(X, number_of_coordinates(M, ℝ)))
end
function get_coordinates_orthonormal!(::DefaultManifold, c, p, X, ::ComplexNumbers)
    m = length(X)
    return copyto!(c, [reshape(real(X), m); reshape(imag(X), m)])
end
function get_vector_orthonormal!(M::DefaultManifold, Y, p, c, ::RealNumbers)
    return copyto!(Y, reshape(c, representation_size(M)))
end
function get_vector_diagonalizing!(
    M::DefaultManifold,
    Y,
    p,
    c,
    ::DiagonalizingOrthonormalBasis{ℝ},
)
    return copyto!(Y, reshape(c, representation_size(M)))
end
function get_vector_orthonormal!(
    M::DefaultManifold{T,ℂ},
    Y,
    p,
    c,
    ::ComplexNumbers,
) where {T}
    n = div(length(c), 2)
    return copyto!(Y, reshape(c[1:n] + c[(n + 1):(2n)] * 1im, representation_size(M)))
end

injectivity_radius(::DefaultManifold) = Inf

@inline inner(::DefaultManifold, p, X, Y) = dot(X, Y)

log!(::DefaultManifold, Y, p, q) = (Y .= q .- p)

@generated function manifold_dimension(::DefaultManifold{T,𝔽}) where {T,𝔽}
    return length(T.parameters) == 0 ? 1 : *(T.parameters...) * real_dimension(𝔽)
end

number_system(::DefaultManifold{T,𝔽}) where {T,𝔽} = 𝔽

norm(::DefaultManifold, p, X) = norm(X)

project!(::DefaultManifold, q, p) = copyto!(q, p)
project!(::DefaultManifold, Y, p, X) = copyto!(Y, X)

@generated representation_size(::DefaultManifold{T}) where {T} = tuple(T.parameters...)

function Base.show(io::IO, ::DefaultManifold{N,𝔽}) where {N,𝔽}
    return print(io, "DefaultManifold($(join(N.parameters, ", ")); field = $(𝔽))")
end

function parallel_transport_to!(::DefaultManifold, Y, p, X, q)
    return copyto!(Y, X)
end

zero_vector(::DefaultManifold, p) = zero(p)

zero_vector!(::DefaultManifold, Y, p) = fill!(Y, 0)


text:
connector_comp_name(i::Int) = Symbol("ConnectorComp$i")

# Return the datatype to use for instance variables/parameters
function _instance_datatype(md::ModelDef, def::DatumDef)
    dtype = def.datatype == Number ? number_type(md) : def.datatype
    dims = dimensions(def)
    num_dims = dim_count(def)

    if num_dims == 0
        T = ScalarModelParameter{dtype}

    elseif dims[1] != :time
        T = Array{dtype, num_dims}

    else
        if isuniform(md)
            first, stepsize = first_and_step(md)
            T = TimestepArray{FixedTimestep{first, stepsize}, Union{dtype, Missing}, num_dims}
        else
            times = time_labels(md)
            T = TimestepArray{VariableTimestep{(times...,)}, Union{dtype, Missing}, num_dims}
        end
    end

    # @info "_instance_datatype returning $T"
    return T
end

# Create the Ref or Array that will hold the value(s) for a Parameter or Variable
function _instantiate_datum(md::ModelDef, def::DatumDef)
    dtype = _instance_datatype(md, def)
    dims = dimensions(def)
    num_dims = length(dims)

    # Scalar datum
    if num_dims == 0
        value = dtype(0)

    # Array datum, with :time dimension
    elseif dims[1] == :time

        if num_dims == 1
            value = dtype(dim_count(md, :time))
        else
            counts = dim_counts(md, Vector{Symbol}(dims))
            value = dtype <: AbstractArray ? dtype(undef, counts...) : dtype(counts...)
        end

    # Array datum, without :time dimension
    else
        # TBD: Handle unnamed indices properly
        counts = dim_counts(md, Vector{Symbol}(dims))
        value = dtype <: AbstractArray ? dtype(undef, counts...) : dtype(counts...)
    end

    return value
end

"""
    _instantiate_component_vars(md::ModelDef, comp_def::ComponentDef)

Instantiate a component `comp_def` in the model `md` and its variables (but not its parameters).
Return the resulting ComponentInstance.
"""
function _instantiate_component_vars(md::ModelDef, comp_def::ComponentDef)
    comp_name = name(comp_def)
    var_defs = variables(comp_def)

    names  = ([name(vdef) for vdef in var_defs]...,)
    types  = Tuple{[_instance_datatype(md, vdef) for vdef in var_defs]...}
    values = [_instantiate_datum(md, def) for def in var_defs]

    return ComponentInstanceVariables(names, types, values)
end

# Save a reference to the model's dimension dictionary to make it
# available in calls to run_timestep.
function save_dim_dict_reference(mi::ModelInstance)
    dim_dict = dim_value_dict(mi.md)

    for ci in values(mi.components)
        ci.dim_dict = dim_dict
    end

    return nothing
end

function build(m::Model)
    # Reference a copy in the ModelInstance to avoid changes underfoot
    m.mi = build(copy(m.md))
    return nothing
end

function build(md::ModelDef)

    add_connector_comps(md)

    # check if all parameters are set
    not_set = unconnected_params(md)
    if ! isempty(not_set)
        params = join(not_set, " ")
        msg = "Cannot build model; the following parameters are not set: $params"
        error(msg)
    end

    var_dict = Dict{Symbol, Any}()                 # collect all var defs and
    par_dict = Dict{Symbol, Dict{Symbol, Any}}()   # store par values as we go

    comp_defs = compdefs(md)
    for comp_def in comp_defs
        comp_name = name(comp_def)
        var_dict[comp_name] = _instantiate_component_vars(md, comp_def)
        par_dict[comp_name] = Dict()  # param value keyed by param name
    end

    # Iterate over connections to create parameters, referencing storage in vars
    for ipc in internal_param_conns(md)
        comp_name = ipc.src_comp_name

        vars = var_dict[comp_name]
        var_value_obj = get_property_obj(vars, ipc.src_var_name)

        par_values = par_dict[ipc.dst_comp_name]
        par_values[ipc.dst_par_name] = var_value_obj
    end

    for ext in external_param_conns(md)
        comp_name = ext.comp_name
        param = external_param(md, ext.external_param)
        par_values = par_dict[comp_name]
        par_values[ext.param_name] = param isa ScalarModelParameter ? param : value(param)
    end

    # Make the external parameter connections for the hidden ConnectorComps.
    # Connect each :input2 to its associated backup value.
    for (i, backup) in enumerate(md.backups)
        comp_name = connector_comp_name(i)
        param = external_param(md, backup)

        par_values = par_dict[comp_name]
        par_values[:input2] = param isa ScalarModelParameter ? param : value(param)
    end

    mi = ModelInstance(md)

    # instantiate parameters
    for comp_def in comp_defs
        comp_name = name(comp_def)

        vars = var_dict[comp_name]

        par_values = par_dict[comp_name]
        pnames = Tuple(parameter_names(comp_def))
        pvals  = [par_values[pname] for pname in pnames]
        ptypes = Tuple{map(typeof, pvals)...}
        pars = ComponentInstanceParameters(pnames, ptypes, pvals)

        first = first_period(md, comp_def)
        last = last_period(md, comp_def)

        ci = ComponentInstance{typeof(vars), typeof(pars)}(comp_def, vars, pars, first, last, comp_name)
        add_comp!(mi, ci)
    end

    save_dim_dict_reference(mi)
    return mi
end

"""
    create_marginal_model(base::Model, delta::Float64=1.0)

Create a `MarginalModel` where `base` is the baseline model and `delta` is the
difference used to create the `marginal` model.  Return the resulting `MarginaModel`
which shares the internal `ModelDef` between the `base` and `marginal`.
"""
function create_marginal_model(base::Model, delta::Float64=1.0)
    # Make sure the base has a ModelInstance before we copy since this
    # copies the ModelDef to avoid being affected by later changes.
    if base.mi === nothing
        build(base)
    end

    # Create a marginal model, which shares the internal ModelDef between base and marginal
    mm = MarginalModel(base, delta)
end

function Base.run(mm::MarginalModel, ntimesteps::Int=typemax(Int))
    run(mm.base, ntimesteps=ntimesteps)
    run(mm.marginal, ntimesteps=ntimesteps)
end


text:
module XLSX

import Dates
import Printf.@printf
import ZipFile
import EzXML

# https://github.com/fhs/ZipFile.jl/issues/39
if !hasmethod(Base.bytesavailable, Tuple{ZipFile.ReadableFile})
    Base.bytesavailable(f::ZipFile.ReadableFile) = f.uncompressedsize - f._pos
end

const SPREADSHEET_NAMESPACE_XPATH_ARG = [ "xpath" => "http://schemas.openxmlformats.org/spreadsheetml/2006/main" ]

include("types.jl")
include("cellref.jl")
include("sst.jl")
include("stream.jl")
include("table.jl")
include("relationship.jl")
include("read.jl")
include("workbook.jl")
include("worksheet.jl")
include("cell.jl")
include("styles.jl")
include("write.jl")

end # module XLSX


text:
#####
##### `reshape`
#####

function rrule(::typeof(reshape), A::AbstractArray, dims::Tuple{Vararg{Int}})
    A_dims = size(A)
    function reshape_pullback(Ȳ)
        return (NO_FIELDS, reshape(Ȳ, A_dims), DoesNotExist())
    end
    return reshape(A, dims), reshape_pullback
end

function rrule(::typeof(reshape), A::AbstractArray, dims::Int...)
    A_dims = size(A)
    function reshape_pullback(Ȳ)
        ∂A = reshape(Ȳ, A_dims)
        ∂dims = broadcast(_ -> DoesNotExist(), dims)
        return (NO_FIELDS, ∂A, ∂dims...)
    end
    return reshape(A, dims...), reshape_pullback
end

#####
##### `hcat` (🐈)
#####

function rrule(::typeof(hcat), A::AbstractArray, Bs::AbstractArray...)
    function hcat_pullback(Ȳ)
        Xs = (A, Bs...)
        ntuple(length(Bs) + 2) do full_i
            full_i == 1 && return NO_FIELDS

            i = full_i - 1
            l = mapreduce(j->size(Xs[j], 2), Base.add_sum, 1:i-1; init=0)
            u = l + size(Xs[i], 2)
            dim = u > l + 1 ? (l+1:u) : u
            # NOTE: The copy here is defensive, since `selectdim` returns a view which we can
            # materialize with `copy`
            copy(selectdim(Ȳ, 2, dim))
        end
    end
    return hcat(A, Bs...), hcat_pullback
end

function rrule(::typeof(reduce), ::typeof(hcat), As::AbstractVector{<:AbstractVecOrMat})
    function reduce_hcat_pullback(ΔY)
        sizes = size.(As, 2)
        cumsizes = cumsum(sizes)
        ∂As = map(cumsizes, sizes) do post, diff
            pre = post - diff + 1
            return ΔY[:, pre:post]
        end
        return (NO_FIELDS, DoesNotExist(), ∂As)
    end
    return reduce(hcat, As), reduce_hcat_pullback
end

#####
##### `vcat`
#####

function rrule(::typeof(vcat), A::AbstractArray, Bs::AbstractArray...)
    function vcat_pullback(Ȳ)
        n = size(A, 1)
        ∂A = copy(selectdim(Ȳ, 1, 1:n))
        ∂Bs = ntuple(length(Bs)) do i
            l = n + mapreduce(j->size(Bs[j], 1), Base.add_sum, 1:i-1; init=0)
            u = l + size(Bs[i], 1)
            copy(selectdim(Ȳ, 1, l+1:u))
        end
        return (NO_FIELDS, ∂A, ∂Bs...)
    end
    return vcat(A, Bs...), vcat_pullback
end

function rrule(::typeof(reduce), ::typeof(vcat), As::AbstractVector{<:AbstractVecOrMat})
    function reduce_vcat_pullback(ΔY)
        sizes = size.(As, 1)
        cumsizes = cumsum(sizes)
        ∂As = map(cumsizes, sizes) do post, diff
            pre = post - diff + 1
            return ΔY[pre:post, :]
        end
        return (NO_FIELDS, DoesNotExist(), ∂As)
    end
    return reduce(vcat, As), reduce_vcat_pullback
end

#####
##### `fill`
#####

function rrule(::typeof(fill), value::Any, dims::Tuple{Vararg{Int}})
    function fill_pullback(Ȳ)
        return (NO_FIELDS, sum(Ȳ), DoesNotExist())
    end
    return fill(value, dims), fill_pullback
end

function rrule(::typeof(fill), value::Any, dims::Int...)
    function fill_pullback(Ȳ)
        return (NO_FIELDS, sum(Ȳ), ntuple(_->DoesNotExist(), length(dims))...)
    end
    return fill(value, dims), fill_pullback
end


text:
using Pkg; Pkg.activate(@__DIR__)
using Comrade
using Distributions
using Pathfinder
using AdvancedHMC
using Plots
using Zygote
using Parameters
using GalacticOptim
using Optim
using NLopt

# load eht-imaging we use this to load eht data
load_ehtim()
# To download the data visit https://doi.org/10.25739/g85n-f134
obs = ehtim.obsdata.load_uvfits(joinpath(@__DIR__, "SR1_M87_2017_096_hi_hops_netcal_StokesI.uvfits"))
obs.add_scans()
# kill 0-baselines since we don't care about
# large scale flux and make scan-average data
obs = obs.flag_uvdist(uv_min=0.1e9).avg_coherent(0.0, scan_avg=true).add_fractional_noise(0.02)
# extract log closure amplitudes and closure phases
vis = extract_vis(obs)
lklhd = RadioLikelihood(vis)

# Build the Model. Here we we a struct to hold some caches
# This will be useful to hold precomputed caches

struct Model{C,G}
    cache::C
    gcache::G
    fov::Float64
    npix::Int
end

function (model::Model)(θ)
    ;c, f, gamp, gphase = θ
    # Construct the image model
    cmat = f.*reshape(c, model.npix, model.npix)
    img = IntensityMap(cmat, model.fov, model.fov, BSplinePulse{3}())
    m = modelimage(img, cache)
    # Now corrupt the model with Gains
    g = gamp.*cis.(gphase)
    Comrade.GainModel(gcache, g, m)
end

# First we define the station gain priors
distamp = (AA = LogNormal(0.0, 0.1),
         AP = LogNormal(0.0, 0.1),
         LM = LogNormal(0.0, 0.5),
         AZ = LogNormal(0.0, 0.1),
         JC = LogNormal(0.0, 0.1),
         PV = LogNormal(0.0, 0.1),
         SM = LogNormal(0.0, 0.1))
# We will use AA as a reference array
distphase = (AA = VonMises(0.0, 10000.0),
             AP = VonMises(0.0, 1/π^2),
             LM = VonMises(0.0, 1/π^2),
             AZ = VonMises(0.0, 1/π^2),
             JC = VonMises(0.0, 1/π^2),
             PV = VonMises(0.0, 1/π^2),
             SM = VonMises(0.0, 1/π^2))

npix = 8
fov = μas2rad(80.0)
prior = (
          c = Dirichlet(npix*npix, 1.0),
          f = Uniform(0.2, 0.9),
          gamp = Comrade.GainPrior(distamp, scantable(vis)),
          gphase = Comrade.GainPrior(distphase, scantable(vis))
        )

# Construct cache's and design matrices
gcache = Comrade.GainCache(scantable(vis))
cache = create_cache(Comrade.DFTAlg(vis), IntensityMap(rand(npix,npix), μas2rad(fov), μas2rad(fov), BSplinePulse{3}()))

mms = Model(cache, gcache, fov, npix)

post = Posterior(lklhd, prior, mms)
tpost = asflat(post)
# We will use HMC to sample the posterior.
# First to reduce burn in we use pathfinder
ndim = dimension(tpost)
f, tr = GalacticOptim.OptimizationFunction(post, GalacticOptim.AutoReverseDiff())
x0 = inverse(tr, rand(post.prior))
prob = GalacticOptim.OptimizationProblem(f, x0, nothing, lb=fill(-15.0, dimension(tr)), ub=fill(15.0, dimension(tr)))
sol, xopt = solve(prob, BBO_adaptive_de_rand_1_bin_radiuslimited(), tr; maxiters=1000_000)

prob = GalacticOptim.OptimizationProblem(f, sol.u, nothing)
sol, xopt = solve(prob, LBFGS(), tr)

residual(mms(xopt), vis)

q, phi, _ = multipathfinder(post, 100)
# now we sample using hmc
metric = DiagEuclideanMetric(ndim)
chain, stats = sample(post, HMC(;metric, autodiff=AD.ReverseDiffBackend()), 2000; nadapts=1000, init_params=xopt)


text:
#!/usr/bin/env julia

using Base.Test
using DataFrames
import SloanDigitalSkySurvey: SDSS, PSF, WCSUtils

include("test_sdss.jl")
include("test_wcs.jl")


text:
##############################################
######### Define all processing ##############
##############################################

#This is an abstract function
"""
It will generally be called as
       new_timed_events = process_event(time, state, event)
It will generate 0 or more new timed events based on the current event
"""
function process_event end

function process_event(time::Float64, state::State, es_event::EndSimEvent)::Vector{TimedEvent}
    println("Ending simulation at time $time.")
    return []
end

function process_event(time::Float64, state::State, ls_event::LogStateEvent)::Vector{TimedEvent}
    println("Logging state at time $time.")
    println(state)
    return []
end;

"""
Process an arrival event adding a new job to the system

Returns all events created in the process
"""
function process_event(time::Float64, state::State, arrival_event::ExternalArrivalEvent)::Vector{TimedEvent}
    new_timed_events = TimedEvent[]

    q = get_entry_queue(state)
    if !is_full(q, state)
        new_job(q, time, state)
        #if this is the only job on the server engage service
        num_in_queue(q, state) == 1 && push!(new_timed_events,
                                    TimedEvent(EndOfServiceAtQueueEvent(q), time + next_service_time(state, q)))
    else
        #Finds new queue using overflow matrix
        next_q = get_next_queue(q, state, :overflow)
        if next_q > 0
            transit_time = time + travel_time(state)
            new_transit(time, transit_time, state)
            push!(new_timed_events, TimedEvent(InTransitEvent(next_q), transit_time))
        else
            failed_arrival(time, state)
        end
    end

    #prepare next arrival
    push!(new_timed_events, TimedEvent(ExternalArrivalEvent(), time + next_arrival_time(state)))
    return new_timed_events
end

"""
Process an end of service event by sending the next job in the specified queue either to transit to its next destination
or out of the sysytem.

Returns all events created in the process
"""
function process_event(time::Float64, state::State, eos_event::EndOfServiceAtQueueEvent)::Vector{TimedEvent}
    q = eos_event.q
    new_timed_events = TimedEvent[]

    #if another customer in the queue then start a new service
    if num_in_queue(q, state) > 1
        push!(new_timed_events, TimedEvent(EndOfServiceAtQueueEvent(q), time + next_service_time(state, q)))
    end

    next_q = get_next_queue(q, state, :routing)
    if next_q > 0
        transit_time = time + travel_time(state)
        queue_to_transit(q, time, transit_time, state)
        push!(new_timed_events, TimedEvent(InTransitEvent(next_q), transit_time))
    else
        pop_queue(q, time, state)
    end
    return new_timed_events
end

"""
Process a transit event by removing the next job from transit and attemping to add it to its destination queue

Returns all events created in the process
"""
function process_event(time::Float64, state::State, transit_event::InTransitEvent)::Vector{TimedEvent}
    new_timed_events = TimedEvent[]

    q = transit_event.q
    if !is_full(q, state)
        transit_to_queue(q, time, state)
        #if this is the only job on the server engage service
        num_in_queue(q, state) == 1 && push!(new_timed_events,
                                    TimedEvent(EndOfServiceAtQueueEvent(q), time + next_service_time(state, q)))
    else
        #Finds new queue using overflow matrix
        next_q = get_next_queue(q, state, :overflow)
        if next_q > 0
            transit_time = time + travel_time(state)
            update_transit(time, transit_time, state)
            push!(new_timed_events, TimedEvent(InTransitEvent(next_q), transit_time))
        else
            pop_transit(time, state)
        end
    end
    return new_timed_events
end

""" simulate(scenario, [job_tracking, max_time, log_times, call_back])

Runs a simulation of the given scenario until the max_time is reached.
The simulation can be run with, no job tracking (:none) as default, only tracking job entry and exit times (:times) or
tracking full job history (:full).

An overview of the system state will be printed at the provided log_times.
The call_back function will be called each time an event occurs.
"""
function simulate(scenario::NetworkParameters;
                    job_tracking::Symbol = :none,
                    max_time::Real = 10.0,
                    log_times::Vector{Float64} = Float64[],
                    call_back = (time,state) -> clear_left(state))
    # Initializes the the state depeding on simulation mode
    if job_tracking == :none
        state = NetworkState(scenario)
    elseif job_tracking == :times
        state = TrackedState(scenario)
    elseif job_tracking == :full
        state = FullTrackedState(scenario)
    else
        error("$mode is not a valid simulation mode")
    end

    init_timed_event = TimedEvent(ExternalArrivalEvent(), 0.0)

    #The event queue
    priority_queue = BinaryMinHeap{TimedEvent}()

    #Put the standard events in the queue
    push!(priority_queue, init_timed_event)
    push!(priority_queue, TimedEvent(EndSimEvent(), max_time))
    for lt in log_times
        push!(priority_queue,TimedEvent(LogStateEvent(), lt))
    end

    time = 0.0
    call_back(time, state)

    #The main discrete event simulation loop
    while true
        #Get the next event
        timed_event = pop!(priority_queue)

        #advance the time
        time = timed_event.time

        #Act on the event
        new_timed_events = process_event(time, state, timed_event.event)

        #if the event was an end of simulation then stop
        isa(timed_event.event, EndSimEvent) && break

        #The event may spawn 0 or more events which we put in the priority queue
        for new_event in new_timed_events
            push!(priority_queue, new_event)
        end

        call_back(time, state)
    end
end


text:
include("transforms.jl")
include("kernel.jl")
include("model.jl")
include("prediction.jl")
include("statistics.jl")
include("validation.jl")
using Plots
using Random
using Printf

##SYNTHETIC DATA TEST

## TESTING CONFIDENCE INTERVALS FOR SYNTHETIC DATASET
## HOW DOES PROXIMITY TO KNOWN DATA POINTS AFFECT WIDTH OF CONFIDENCE INTERVAL?
## HOW TO THETA AND LAMBDA EFFECT THE WIDTH OF THE CONFIDENCE INTERVAL?

function plt(f, a, b, numpts=500)
    h = (b-a)/(numpts-1)
    y = zeros(1,numpts)
    x = collect(a:h:b)
    y = map(f, x)
    display(plot(x, y))
end

#Random.seed!(1234);
Random.seed!(4312)
h = 50
smaster=collect(0:1/h:1-1/h)
s0 = [4]
K1 = K(smaster, smaster, 120, rbf)
zo = sample(K1) #zero mean data with fixed covariance structure
mean = (minimum(zo) < 0 ? abs(minimum(zo)) + 12 : 12.1);
zo = zo .+  mean; # sample positive data with fixed mean
zmaster = map(x -> invBoxCox(x, 2), zo)
pλ = x -> 1 #let lambda be uniformly distributed (value of constant is immaterial, because of renormalization)
pθ = x -> 1 #let theta be uniformly distributed
#display(plot(repeat(smaster, 1, 2), hcat(zo, z), seriestype = :scatter, layout=2, title = "boxcox"))
ind_master = shuffle(collect(1:h))
ind = ind_master[1:15]
s = (collect(0:1/h:1-1/h))[ind]; z  = zmaster[ind]
display(plot(smaster, zmaster, seriestype = :scatter))
display(plot!(s, z, seriestype = :scatter))
#savefig("masters.png")
#println(s)
range_theta = [50 250]
range_lambda = [-3 6]
X = ones(length(s), 1) #assume constant mean b/c of lack of covariate information
X0 = ones(1, 1)

totest = [30 31 32 33 34 35 36] #selected testing location indices (of smaster)
#totest = [11 12 13 14 15]
results = Array{Float64}(undef, length(totest)) #preallocate space for medians and interval widths
intwidths = Array{Float64}(undef, length(totest))

if false #sanity check + plot data
    pdf = model(X, X0, s, [smaster[totest[3]]], boxCox, boxCoxPrime, pθ, pλ, z, range_theta, range_lambda)
    cdf = z0 ->  int1D(pdf, 0, z0, "3")
    @time begin
        un = cdf(23)
    end
    cdfn = x -> cdf(x)/un #normalized cdf
    pdfn = x -> pdf(x)/un #normalized pdf
    plt(pdfn, 0, 23)
    savefig("results//synthetic//data3.jpg")
    med = bisection(x -> cdfn(x)-0.5, 1e-3, 25, 1e-3, 10)
    println(med)
end

if false #test prediction at multiple locations
    for i=1:length(totest) #single-point predictions
        println(i)
        @time begin
        pdf = model(X, X0, s, [smaster[totest[i]]], boxCox, boxCoxPrime, pθ, pλ, z, range_theta, range_lambda)
        cdf = z0 ->  int1D(pdf, 0, z0, "3")
        @time begin
        un = cdf(23)
        end
        cdfn = x -> cdf(x)/un #normalized cdf
        pdfn = x -> pdf(x)/un #normalized pdf
        @time begin
        med = bisection(x -> cdfn(x)-0.5, 1e-3, 25, 1e-3, 10)
        #med, val = gridsearch(pdfn, 0, 30, 20) #MAP estimation
        @printf("ith pt med: %f \n", med)
        end
        @time begin
        intwidths[i] = confidence(pdfn, med) #compute width of credible interval corresponding to 95% confidence level
        @printf("ith pt intwidth: %f \n", intwidths[i])
        end
        results[i] = med
        end
    end
display(plot(smaster, zmaster, seriestype = :scatter, label = "latent process", title = "\\theta [50, 150], \\lambda [1, 3], actual (\\theta = 120, \\lambda = 5)"));
display(plot!(s, z, seriestype = :scatter, label = "training"));
display(plot!(smaster[totest'], results, seriestype = :scatter, label = "predicted median"));
#display(plot!(test_x, test_y, seriestype = :scatter));
display(plot!(repeat(smaster[totest]', 2), vcat(results-intwidths, results+intwidths), seriestype = :scatter, color=:red, label = "95% confidence"));
#display(plot!(test_x, test_y))
savefig("synthetic7.png")
end

if true #cross validation
    ss, Xs, Ys = cross_validate(X, s, boxCox, boxCoxPrime, pθ, pλ, z, range_theta, range_lambda, 1000, 2, 20)
    display(plot(Xs, Ys, layout = length(z), legend=nothing, xtickfont = font(10, "Courier"), lw=0.5))
    display(plot!(z', layout = length(z), legend=nothing, seriestype = :vline, xtickfont = font(10, "Courier"), lw=0.5))
    savefig("results//synthetic//synthetic_delete_one_cross_validate4.jpg")
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using SpecialFunctions
using PRAS # From added registry

function map_log_gamma(v::Vector{N}) where {N <: Number}
  return map(loggamma, v)
end


text:
function createPlotWindow{F<:(@compat AbstractFloat)}(X::Vector{F}, Y::Vector{F};
                                            width::Int = 40,
                                            height::Int = 15,
                                            margin::Int = 3,
                                            padding::Int = 1,
                                            title::(@compat AbstractString) = "",
                                            border::Symbol = :solid,
                                            labels::Bool = true,
                                            xlim::Vector = [0.,0.],
                                            ylim::Vector = [0.,0.])
  length(xlim) == length(ylim) == 2 || throw(ArgumentError("xlim and ylim must only be vectors of length 2"))
  margin >= 0 || throw(ArgumentError("Margin must be greater than or equal to 0"))
  length(X) == length(Y) || throw(DimensionMismatch("X and Y must be the same length"))
  width = max(width, 5)
  height = max(height, 2)
  minX = float(minimum(xlim))
  maxX = float(maximum(xlim))
  if minX == 0. && maxX == 0.
    minX = minimum(X)
    maxX = maximum(X)
  end
  diffX = maxX - minX
  if diffX == 0
    t = minX
    minX = t / 2
    maxX = t + t / 2
  end
  minX, maxX = xlim == [0.,0.] ? plottingRange(minX, maxX) : plottingRangeNarrow(minX, maxX)
  minY = float(minimum(ylim))
  maxY = float(maximum(ylim))
  if minY == 0. && maxY == 0.
    minY = minimum(Y)
    maxY = maximum(Y)
  end
  diffY = maxY - minY
  if diffY == 0
    t = minY
    minY = t / 2
    maxY = t + t / 2
  end
  minY, maxY = ylim == [0.,0.] ? plottingRange(minY, maxY) : plottingRangeNarrow(minY, maxY)
  plotOriginX = minX; plotWidth = maxX - plotOriginX
  plotOriginY = minY; plotHeight = maxY - plotOriginY

  canvas = BrailleCanvas(width, height,
                         plotOriginX = plotOriginX, plotOriginY = plotOriginY,
                         plotWidth = plotWidth, plotHeight = plotHeight)
  newPlot = Plot(canvas, title=title, margin=margin,
                 padding=padding, border=border, showLabels=labels)

  minXString=string(isinteger(minX) ? safeRound(minX) : minX)
  maxXString=string(isinteger(maxX) ? safeRound(maxX) : maxX)
  minYString=string(isinteger(minY) ? safeRound(minY) : minY)
  maxYString=string(isinteger(maxY) ? safeRound(maxY) : maxY)
  annotate!(newPlot, :l, 1, maxYString)
  annotate!(newPlot, :l, height, minYString)
  annotate!(newPlot, :bl, minXString)
  annotate!(newPlot, :br, maxXString)
  if minY < 0 < maxY
    for i in linspace(minX, maxX, width*2)
      setPoint!(newPlot, i, 0., :white)
    end
  end
  if minX < 0 < maxX
    for i in linspace(minY, maxY, height*4)
      setPoint!(newPlot, 0., i, :white)
    end
  end
  newPlot
end

function scatterplot!{T<:Canvas,F<:Real,R<:Real}(plot::Plot{T}, X::AbstractVector{F}, Y::AbstractVector{R};
                                                 color::Symbol = :auto,
                                                 name::(@compat AbstractString) = "",
                                                 args...)
  X = convert(Vector{(@compat AbstractFloat)},X)
  Y = convert(Vector{(@compat AbstractFloat)},Y)
  color = color == :auto ? nextColor!(plot) : color
  name == "" || autoAnnotate!(plot, :r, name, color)
  setPoint!(plot, X, Y, color)
end

function scatterplot{F<:Real,R<:Real}(X::AbstractVector{F}, Y::AbstractVector{R};
                                      color::Symbol=:auto,
                                      name::(@compat AbstractString) = "",
                                      args...)
  X = convert(Vector{(@compat AbstractFloat)},X)
  Y = convert(Vector{(@compat AbstractFloat)},Y)
  newPlot = createPlotWindow(X, Y; args...)
  color = color == :auto ? nextColor!(newPlot) : color
  name == "" || autoAnnotate!(newPlot, :r, name, color)
  setPoint!(newPlot, X, Y, color)
end

function lineplot!{T<:Canvas,F<:Real,R<:Real}(plot::Plot{T}, X::AbstractVector{F}, Y::AbstractVector{R};
                                              color::Symbol=:auto,
                                              name::(@compat AbstractString) = "",
                                              args...)
  X = convert(Vector{(@compat AbstractFloat)},X)
  Y = convert(Vector{(@compat AbstractFloat)},Y)
  color = color == :auto ? nextColor!(plot) : color
  name == "" || autoAnnotate!(plot, :r, name, color)
  drawLine!(plot, X, Y, color)
end

function lineplot{F<:Real,R<:Real}(X::AbstractVector{F}, Y::AbstractVector{R};
                                   color::Symbol=:auto,
                                   name::(@compat AbstractString) = "",
                                   args...)
  X = convert(Vector{(@compat AbstractFloat)},X)
  Y = convert(Vector{(@compat AbstractFloat)},Y)
  newPlot = createPlotWindow(X, Y; args...)
  color = color == :auto ? nextColor!(newPlot) : color
  name == "" || autoAnnotate!(newPlot, :r, name, color)
  drawLine!(newPlot, X, Y, color)
end

function lineplot!{T<:Canvas}(plot::Plot{T}, intercept::Real, slope::Real; args...)
  xmin = plot.graphics.plotOriginX
  xmax = plot.graphics.plotOriginX + plot.graphics.plotWidth
  ymin = plot.graphics.plotOriginY
  ymax = plot.graphics.plotOriginY + plot.graphics.plotHeight
  lineplot!(plot, [xmin, xmax], [intercept + xmin*slope, intercept + xmax*slope]; args...)
end

function lineplot!{T<:Canvas}(plot::Plot{T}, Y::Function, x::Range; args...)
  X = collect(x)
  lineplot!(plot, Y, X; args...)
end

function lineplot(Y::Function, x::Range; args...)
  X = collect(x)
  lineplot(Y, X; args...)
end

function lineplot!{T<:Canvas,R<:Real}(plot::Plot{T}, Y::Function, X::AbstractVector{R};
                                      name::(@compat AbstractString) = "",
                                      args...)
  y = convert(Vector{Float64}, [Y(i) for i in X])
  name = name == "" ? string(Y, "(x)") : name
  lineplot!(plot, X, y; name = name, args...)
end

function lineplot{R<:Real}(Y::Function, X::Vector{R};
                           name::(@compat AbstractString) = "",
                           args...)
  y = convert(Vector{Float64}, [Y(i) for i in X])
  name = name == "" ? string(Y, "(x)") : name
  newPlot = lineplot(X, y; name = name, args...)
  xlabel!(newPlot, "x")
  ylabel!(newPlot, "f(x)")
end

function lineplot!{T<:Canvas}(plot::Plot{T}, Y::Function, startx::Real, endx::Real; args...)
  diff = abs(endx - startx)
  X = collect(startx:(diff/(3*ncols(plot.graphics))):endx)
  lineplot!(plot, Y, X; args...)
end

function lineplot(Y::Function, startx::Real, endx::Real; width::Int = 40, args...)
  diff = abs(endx - startx)
  X = collect(startx:(diff/(3*width)):endx)
  lineplot(Y, X; width = width, args...)
end

function lineplot(Y::AbstractVector{Function}, startx::Real, endx::Real; args...)
  n = length(Y)
  @assert n > 0
  newPlot = lineplot(Y[1], startx, endx; args...)
  for i = 2:n
    lineplot!(newPlot, Y[i], startx, endx; args...)
  end
  newPlot
end

function lineplot{F<:Real,R<:Real}(X::Range{F}, Y::Range{R}; args...)
  lineplot(collect(X), collect(Y); args...)
end

function lineplot{F<:Real}(X::Range, Y::AbstractVector{F}; args...)
  lineplot(collect(X), Y; args...)
end

function lineplot{F<:Real}(X::AbstractVector{F}, Y::Range; args...)
  lineplot(X, collect(Y); args...)
end

function lineplot!{T<:Canvas,F<:Real,R<:Real}(plot::Plot{T}, X::Range{F}, Y::Range{R}; args...)
  lineplot!(plot, collect(X), collect(Y); args...)
end

function lineplot!{T<:Canvas,F<:Real}(plot::Plot{T}, X::Range, Y::AbstractVector{F}; args...)
  lineplot!(plot, collect(X), Y; args...)
end

function lineplot!{T<:Canvas,F<:Real}(plot::Plot{T}, X::AbstractVector{F}, Y::Range; args...)
  lineplot!(plot, X, collect(Y); args...)
end

function computeStairLines{F<:Real,R<:Real}(X::AbstractVector{F}, Y::AbstractVector{R}, style::Symbol)
  if style == :post
    xVec = zeros(length(X) * 2 - 1)
    yVec = zeros(length(X) * 2 - 1)
    xVec[1] = X[1]
    yVec[1] = Y[1]
    o = 0
    for i = 2:(length(X))
      xVec[i + o] = X[i]
      xVec[i + o + 1] = X[i]
      yVec[i + o] = Y[i-1]
      yVec[i + o + 1] = Y[i]
      o += 1
    end
    xVec, yVec
  elseif style == :pre
    xVec = zeros(length(X) * 2 - 1)
    yVec = zeros(length(X) * 2 - 1)
    xVec[1] = X[1]
    yVec[1] = Y[1]
    o = 0
    for i = 2:(length(X))
      xVec[i + o] = X[i-1]
      xVec[i + o + 1] = X[i]
      yVec[i + o] = Y[i]
      yVec[i + o + 1] = Y[i]
      o += 1
    end
    xVec, yVec
  end
end

function stairs!{T<:Canvas,F<:Real,R<:Real}(plot::Plot{T}, X::AbstractVector{F}, Y::AbstractVector{R};
                                            style::Symbol = :post,
                                            args...)
  xVec, yVec = computeStairLines(X, Y, style)
  lineplot!(plot, xVec, yVec; args...)
end

function stairs{F<:Real,R<:Real}(X::Vector{F}, Y::Vector{R};
                                 style::Symbol = :post,
                                 args...)
  xVec, yVec = computeStairLines(X, Y, style)
  lineplot(xVec, yVec; args...)
end