text:
# License for this file: MIT (expat)
# Copyright 2021, DLR Institute of System Dynamics and Control

module ModiaPlot_CairoMakie

# Constants
const headingSize = 10

const path = dirname(dirname(@__FILE__))   # Absolute path of package directory
const Version = "0.4.4"
const Date = "2021-11-23"

println("Importing ModiaPlot_CairoMakie Version $Version ($Date) - this takes some time due to CairoMakie import")

import ModiaResult
import Colors
import Measurements
import MonteCarloMeasurements
using  Unitful

using  CairoMakie
include("$(ModiaResult.path)/src/plot.jl")

const showFigureStringInDiagram = false
const callDisplayFunction = false
const reusePossible = true

const Makie_Point2f = isdefined(CairoMakie, :Point2f) ? Point2f : Point2f0
include("$(ModiaResult.path)/src/makie.jl")

showFigure(figureNumber::Int) = nothing

function saveFigure(figureNumber::Int, fileName; kwargs...)::Nothing
    if haskey(figures, figureNumber)
        fig = figures[figureNumber].fig
        fullFileName = joinpath(pwd(), fileName)
        println("... save plot in file: \"$fullFileName\"")
        save(fileName, fig; kwargs...)
    else
        @warn "saveFigure: figure $figureNumber is not defined."
    end
    return nothing
end

function closeFigure(figureNumber::Int)::Nothing
    delete!(figures,figureNumber)
    return nothing
end

"""
    closeAllFigures()

Close all figures.
"""
function closeAllFigures()::Nothing
    if length(figures) > 0
        empty!(figures)
    end
    return nothing
end

end


text:
module __RTS
    function  make_tuple(args...)
        return args
    end

    op_eq = (==)

    op_slt = (<)

    op_minus = (-)

    op_plus = (+)

    function op_external(args...)
        error(join(", ", map(repr, args)))
    end

    function op_write_str(_, s)
        print(s)
    end

    function proj(a, b :: Integer)
        a[b + 1]
    end

    function is_tuple(a)
        a isa Tuple
    end

    function op_str_concat(a, b)
        a = a === nothing ? "" : a
        return a * b
    end

    op_int_str = string

    function assert2(a, b)
        @assert a == b
    end
end

ext: jl


text:
# Solution to Q2c (Implement a CrossEntropy)
import Base

type CrossEntropy{T}
    # Just keep these in memory for later inspection, not really necessary.
    costs::Matrix{T}

    # Buffer of the inputs and targets to fprop, needed during bprop.
    x::Matrix{T}
    t::Matrix{T}

    CrossEntropy() = new(
        zeros(T, 0, 0),
        zeros(T, 0, 0),
        zeros(T, 0, 0),
    )
end

target!{T}(l::CrossEntropy{T}, t::Matrix{T}) = l.t = t

function fprop!{T}(l::CrossEntropy{T}, x::Matrix{T})
    l.x = x

    # This is the plain formula of cross-entropy, no special trick here.
    return l.costs = - sum(l.t .* log.(x), 2)
end

function bprop!{T}(l::CrossEntropy{T}, ∂y::Matrix{T})
    # This comes from Q1c.
    # Note that ∂y is usually a vector of 1/N where N is the mini-batch size,
    # this way we average over the minibatch and get stable performance wrt. N.
    # It could also be used to give different weight to different data!
    return ∂y .* (- l.t ./ l.x)
end

params(::CrossEntropy) = ()
grads(::CrossEntropy) = ()

Base.show(io::IO, l::CrossEntropy) = print(io, "CrossEntropy")


text:
export Default

"""
    Default(<keyword arguments>)
    Default(parameters::SegmentParameters)

`Default` objects `show` a .default FastHenry command.

**Keyword Arguments**

- `w`               -- segment width
- `h`               -- segment height
- `sigma`           -- conductivity
- `rho`             -- resistivity
- `nhinc`, `nwinc`  -- integer number of filaments in height and width
- `rh`, `rw`        -- ratio in the height and width
"""
struct Default <: Element
  wh        :: WH
  sigmarho  :: SigmaRho
end
Default(;w=NaN, h=NaN, sigma=NaN, rho=NaN, nhinc=0, nwinc=0, rh=NaN, rw=NaN) =
    Default(WH(w,h,nhinc,nwinc,rh,rw), SigmaRho(sigma,rho))
Default(sp::SegmentParameters) = Default(sp.wh, sp.sigmarho)


text:
export Settings;

struct Settings
  embeddings_size::Integer
  hidden_size::Integer
  batch_size::Integer
  reg_weight::Float64
  sample_size::Integer

  Settings(;
    embeddings_size::Integer=300,
    hidden_size::Integer=200,
    batch_size::Integer=48,
    reg_weight=1e-8,
    sample_size=3000,
  ) = new(embeddings_size, hidden_size, batch_size, reg_weight, sample_size)
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
# This file is a part of DensityInterface.jl, licensed under the MIT License (MIT).

"""
    DensityInterface

Trait-based interface for mathematical/statistical densities and objects
associated with a density.
"""
module DensityInterface

using InverseFunctions
using Test

include("interface.jl")
include("interface_test.jl")

end # module


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
Methods for working with Iterators.
"""
module Iterators

# small dance to make this work from Base or Intrinsics
import ..@__MODULE__, ..parentmodule
const Base = parentmodule(@__MODULE__)
using .Base:
    @inline, Pair, Pairs, AbstractDict, IndexLinear, IndexCartesian, IndexStyle, AbstractVector, Vector,
    tail, SizeUnknown, HasLength, HasShape, IsInfinite, EltypeUnknown, HasEltype, OneTo,
    @propagate_inbounds, @isdefined, @boundscheck, @inbounds, Generator, AbstractRange,
    LinearIndices, (:), |, +, -, !==, !, <=, <, missing, any, _counttuple

import .Base:
    first, last,
    isempty, length, size, axes, ndims,
    eltype, IteratorSize, IteratorEltype,
    haskey, keys, values, pairs,
    getindex, setindex!, get, iterate,
    popfirst!, isdone, peek

export enumerate, zip, rest, countfrom, take, drop, takewhile, dropwhile, cycle, repeated, product, flatten, partition, flatmap

"""
    Iterators.map(f, iterators...)

Create a lazy mapping.  This is another syntax for writing
`(f(args...) for args in zip(iterators...))`.

!!! compat "Julia 1.6"
    This function requires at least Julia 1.6.

# Examples
```jldoctest
julia> collect(Iterators.map(x -> x^2, 1:3))
3-element Vector{Int64}:
 1
 4
 9
```
"""
map(f, args...) = Base.Generator(f, args...)

tail_if_any(::Tuple{}) = ()
tail_if_any(x::Tuple) = tail(x)

_min_length(a, b, ::IsInfinite, ::IsInfinite) = min(length(a),length(b)) # inherit behaviour, error
_min_length(a, b, A, ::IsInfinite) = length(a)
_min_length(a, b, ::IsInfinite, B) = length(b)
_min_length(a, b, A, B) = min(length(a),length(b))

_diff_length(a, b, A, ::IsInfinite) = 0
_diff_length(a, b, ::IsInfinite, ::IsInfinite) = 0
_diff_length(a, b, ::IsInfinite, B) = length(a) # inherit behaviour, error
_diff_length(a, b, A, B) = max(length(a)-length(b), 0)

and_iteratorsize(isz::T, ::T) where {T} = isz
and_iteratorsize(::HasLength, ::HasShape) = HasLength()
and_iteratorsize(::HasShape, ::HasLength) = HasLength()
and_iteratorsize(a, b) = SizeUnknown()

and_iteratoreltype(iel::T, ::T) where {T} = iel
and_iteratoreltype(a, b) = EltypeUnknown()

## Reverse-order iteration for arrays and other collections.  Collections
## should implement iterate etcetera if possible/practical.
"""
    Iterators.reverse(itr)

Given an iterator `itr`, then `reverse(itr)` is an iterator over the
same collection but in the reverse order.

This iterator is "lazy" in that it does not make a copy of the collection in
order to reverse it; see [`Base.reverse`](@ref) for an eager implementation.

Not all iterator types `T` support reverse-order iteration.  If `T`
doesn't, then iterating over `Iterators.reverse(itr::T)` will throw a [`MethodError`](@ref)
because of the missing [`iterate`](@ref) methods for `Iterators.Reverse{T}`.
(To implement these methods, the original iterator
`itr::T` can be obtained from `r = Iterators.reverse(itr)` by `r.itr`.)

# Examples
```jldoctest
julia> foreach(println, Iterators.reverse(1:5))
5
4
3
2
1
```
"""
reverse(itr) = Reverse(itr)

struct Reverse{T}
    itr::T
end
eltype(::Type{Reverse{T}}) where {T} = eltype(T)
length(r::Reverse) = length(r.itr)
size(r::Reverse) = size(r.itr)
IteratorSize(::Type{Reverse{T}}) where {T} = IteratorSize(T)
IteratorEltype(::Type{Reverse{T}}) where {T} = IteratorEltype(T)
last(r::Reverse) = first(r.itr) # the first shall be last

# reverse-order array iterators: assumes more-specialized Reverse for eachindex
@propagate_inbounds function iterate(A::Reverse{<:AbstractArray}, state=(reverse(eachindex(A.itr)),))
    y = iterate(state...)
    y === nothing && return y
    idx, itrs = y
    (A.itr[idx], (state[1], itrs))
end

# Fallback method of `iterate(::Reverse{T})` which assumes the collection has `getindex(::T) and `reverse(eachindex(::T))`
# don't propagate inbounds for this just in case
function iterate(A::Reverse, state=(reverse(eachindex(A.itr)),))
    y = iterate(state...)
    y === nothing && return y
    idx, itrs = y
    (A.itr[idx], (state[1], itrs))
end

reverse(R::AbstractRange) = Base.reverse(R) # copying ranges is cheap
reverse(G::Generator) = Generator(G.f, reverse(G.iter))
reverse(r::Reverse) = r.itr
reverse(x::Union{Number,AbstractChar}) = x
reverse(p::Pair) = Base.reverse(p) # copying pairs is cheap

iterate(r::Reverse{<:Union{Tuple, NamedTuple}}, i::Int = length(r.itr)) = i < 1 ? nothing : (r.itr[i], i-1)

# enumerate

struct Enumerate{I}
    itr::I
end

"""
    enumerate(iter)

An iterator that yields `(i, x)` where `i` is a counter starting at 1,
and `x` is the `i`th value from the given iterator. It's useful when
you need not only the values `x` over which you are iterating, but
also the number of iterations so far. Note that `i` may not be valid
for indexing `iter`; it's also possible that `x != iter[i]`, if `iter`
has indices that do not start at 1. See the `pairs(IndexLinear(),
iter)` method if you want to ensure that `i` is an index.

# Examples
```jldoctest
julia> a = ["a", "b", "c"];

julia> for (index, value) in enumerate(a)
           println("\$index \$value")
       end
1 a
2 b
3 c
```
"""
enumerate(iter) = Enumerate(iter)

length(e::Enumerate) = length(e.itr)
size(e::Enumerate) = size(e.itr)
@propagate_inbounds function iterate(e::Enumerate, state=(1,))
    i, rest = state[1], tail(state)
    n = iterate(e.itr, rest...)
    n === nothing && return n
    (i, n[1]), (i+1, n[2])
end
last(e::Enumerate) = (length(e.itr), e.itr[end])

eltype(::Type{Enumerate{I}}) where {I} = Tuple{Int, eltype(I)}

IteratorSize(::Type{Enumerate{I}}) where {I} = IteratorSize(I)
IteratorEltype(::Type{Enumerate{I}}) where {I} = IteratorEltype(I)

@inline function iterate(r::Reverse{<:Enumerate})
    ri = reverse(r.itr.itr)
    iterate(r, (length(ri), ri))
end
@inline function iterate(r::Reverse{<:Enumerate}, state)
    i, ri, rest = state[1], state[2], tail(tail(state))
    n = iterate(ri, rest...)
    n === nothing && return n
    (i, n[1]), (i-1, ri, n[2])
end

"""
    pairs(IndexLinear(), A)
    pairs(IndexCartesian(), A)
    pairs(IndexStyle(A), A)

An iterator that accesses each element of the array `A`, returning
`i => x`, where `i` is the index for the element and `x = A[i]`.
Identical to `pairs(A)`, except that the style of index can be selected.
Also similar to `enumerate(A)`, except `i` will be a valid index
for `A`, while `enumerate` always counts from 1 regardless of the indices
of `A`.

Specifying [`IndexLinear()`](@ref) ensures that `i` will be an integer;
specifying [`IndexCartesian()`](@ref) ensures that `i` will be a
[`CartesianIndex`](@ref); specifying `IndexStyle(A)` chooses whichever has
been defined as the native indexing style for array `A`.

Mutation of the bounds of the underlying array will invalidate this iterator.

# Examples
```jldoctest
julia> A = ["a" "d"; "b" "e"; "c" "f"];

julia> for (index, value) in pairs(IndexStyle(A), A)
           println("\$index \$value")
       end
1 a
2 b
3 c
4 d
5 e
6 f

julia> S = view(A, 1:2, :);

julia> for (index, value) in pairs(IndexStyle(S), S)
           println("\$index \$value")
       end
CartesianIndex(1, 1) a
CartesianIndex(2, 1) b
CartesianIndex(1, 2) d
CartesianIndex(2, 2) e
```

See also [`IndexStyle`](@ref), [`axes`](@ref).
"""
pairs(::IndexLinear,    A::AbstractArray) = Pairs(A, LinearIndices(A))
pairs(::IndexCartesian, A::AbstractArray) = Pairs(A, CartesianIndices(axes(A)))

# preserve indexing capabilities for known indexable types
# faster than zip(keys(a), values(a)) for arrays
pairs(tuple::Tuple) = Pairs{Int}(tuple, keys(tuple))
pairs(nt::NamedTuple) = Pairs{Symbol}(nt, keys(nt))
pairs(v::Core.SimpleVector) = Pairs(v, LinearIndices(v))
pairs(A::AbstractArray)  = pairs(IndexCartesian(), A)
pairs(A::AbstractVector) = pairs(IndexLinear(), A)
# pairs(v::Pairs) = v # listed for reference, but already defined from being an AbstractDict

length(v::Pairs) = length(getfield(v, :itr))
axes(v::Pairs) = axes(getfield(v, :itr))
size(v::Pairs) = size(getfield(v, :itr))

@propagate_inbounds function _pairs_elt(p::Pairs{K, V}, idx) where {K, V}
    return Pair{K, V}(idx, getfield(p, :data)[idx])
end

@propagate_inbounds function iterate(p::Pairs{K, V}, state...) where {K, V}
    x = iterate(getfield(p, :itr), state...)
    x === nothing && return x
    idx, next = x
    return (_pairs_elt(p, idx), next)
end

@propagate_inbounds function iterate(r::Reverse{<:Pairs}, state=(reverse(getfield(r.itr, :itr)),))
    x = iterate(state...)
    x === nothing && return x
    idx, next = x
    return (_pairs_elt(r.itr, idx), (state[1], next))
end

@inline isdone(v::Pairs, state...) = isdone(getfield(v, :itr), state...)

IteratorSize(::Type{<:Pairs{<:Any, <:Any, I}}) where {I} = IteratorSize(I)
IteratorSize(::Type{<:Pairs{<:Any, <:Any, <:Base.AbstractUnitRange, <:Tuple}}) = HasLength()

function last(v::Pairs{K, V}) where {K, V}
    idx = last(getfield(v, :itr))
    return Pair{K, V}(idx, v[idx])
end

haskey(v::Pairs, key) = (key in getfield(v, :itr))
keys(v::Pairs) = getfield(v, :itr)
values(v::Pairs) = getfield(v, :data) # TODO: this should be a view of data subset by itr
getindex(v::Pairs, key) = getfield(v, :data)[key]
setindex!(v::Pairs, value, key) = (getfield(v, :data)[key] = value; v)
get(v::Pairs, key, default) = get(getfield(v, :data), key, default)
get(f::Base.Callable, v::Pairs, key) = get(f, getfield(v, :data), key)

# zip

struct Zip{Is<:Tuple}
    is::Is
end

"""
    zip(iters...)

Run multiple iterators at the same time, until any of them is exhausted. The value type of
the `zip` iterator is a tuple of values of its subiterators.

!!! note
    `zip` orders the calls to its subiterators in such a way that stateful iterators will
    not advance when another iterator finishes in the current iteration.

See also: [`enumerate`](@ref), [`splat`](@ref Base.splat).

# Examples
```jldoctest
julia> a = 1:5
1:5

julia> b = ["e","d","b","c","a"]
5-element Vector{String}:
 "e"
 "d"
 "b"
 "c"
 "a"

julia> c = zip(a,b)
zip(1:5, ["e", "d", "b", "c", "a"])

julia> length(c)
5

julia> first(c)
(1, "e")
```
"""
zip(a...) = Zip(a)
function length(z::Zip)
    n = _zip_min_length(z.is)
    n === nothing && throw(ArgumentError("iterator is of undefined length"))
    return n
end
function _zip_min_length(is)
    i = is[1]
    n = _zip_min_length(tail(is))
    if IteratorSize(i) isa IsInfinite
        return n
    else
        return n === nothing ? length(i) : min(n, length(i))
    end
end
_zip_min_length(is::Tuple{}) = nothing
size(z::Zip) = _promote_tuple_shape(Base.map(size, z.is)...)
axes(z::Zip) = _promote_tuple_shape(Base.map(axes, z.is)...)
_promote_tuple_shape((a,)::Tuple{OneTo}, (b,)::Tuple{OneTo}) = (intersect(a, b),)
_promote_tuple_shape((m,)::Tuple{Integer}, (n,)::Tuple{Integer}) = (min(m, n),)
_promote_tuple_shape(a, b) = promote_shape(a, b)
_promote_tuple_shape(a, b...) = _promote_tuple_shape(a, _promote_tuple_shape(b...))
_promote_tuple_shape(a) = a
eltype(::Type{Zip{Is}}) where {Is<:Tuple} = Tuple{ntuple(n -> eltype(fieldtype(Is, n)), _counttuple(Is)::Int)...}
#eltype(::Type{Zip{Tuple{}}}) = Tuple{}
#eltype(::Type{Zip{Tuple{A}}}) where {A} = Tuple{eltype(A)}
#eltype(::Type{Zip{Tuple{A, B}}}) where {A, B} = Tuple{eltype(A), eltype(B)}
@inline isdone(z::Zip) = _zip_any_isdone(z.is, Base.map(_ -> (), z.is))
@inline isdone(z::Zip, ss) = _zip_any_isdone(z.is, Base.map(tuple, ss))
@inline function _zip_any_isdone(is, ss)
    d1 = isdone(is[1], ss[1]...)
    d1 === true && return true
    return d1 | _zip_any_isdone(tail(is), tail(ss))
end
@inline _zip_any_isdone(::Tuple{}, ::Tuple{}) = false

@propagate_inbounds iterate(z::Zip) = _zip_iterate_all(z.is, Base.map(_ -> (), z.is))
@propagate_inbounds iterate(z::Zip, ss) = _zip_iterate_all(z.is, Base.map(tuple, ss))

# This first queries isdone from every iterator. If any gives true, it immediately returns
# nothing. It then iterates all those where isdone returned missing, afterwards all those
# it returned false, again terminating immediately if any iterator is exhausted. Finally,
# the results are interleaved appropriately.
@propagate_inbounds function _zip_iterate_all(is, ss)
    d, ds = _zip_isdone(is, ss)
    d && return nothing
    xs1 = _zip_iterate_some(is, ss, ds, missing)
    xs1 === nothing && return nothing
    xs2 = _zip_iterate_some(is, ss, ds, false)
    xs2 === nothing && return nothing
    return _zip_iterate_interleave(xs1, xs2, ds)
end

@propagate_inbounds function _zip_iterate_some(is, ss, ds::Tuple{T,Vararg{Any}}, f::T) where T
    x = iterate(is[1], ss[1]...)
    x === nothing && return nothing
    y = _zip_iterate_some(tail(is), tail(ss), tail(ds), f)
    y === nothing && return nothing
    return (x, y...)
end
@propagate_inbounds _zip_iterate_some(is, ss, ds::Tuple{Any,Vararg{Any}}, f) =
    _zip_iterate_some(tail(is), tail(ss), tail(ds), f)
_zip_iterate_some(::Tuple{}, ::Tuple{}, ::Tuple{}, ::Any) = ()

function _zip_iterate_interleave(xs1, xs2, ds)
    t = _zip_iterate_interleave(tail(xs1), xs2, tail(ds))
    ((xs1[1][1], t[1]...), (xs1[1][2], t[2]...))
end
function _zip_iterate_interleave(xs1, xs2, ds::Tuple{Bool,Vararg{Any}})
    t = _zip_iterate_interleave(xs1, tail(xs2), tail(ds))
    ((xs2[1][1], t[1]...), (xs2[1][2], t[2]...))
end
_zip_iterate_interleave(::Tuple{}, ::Tuple{}, ::Tuple{}) = ((), ())

function _zip_isdone(is, ss)
    d = isdone(is[1], ss[1]...)
    d´, ds = _zip_isdone(tail(is), tail(ss))
    return (d === true || d´, (d, ds...))
end
_zip_isdone(::Tuple{}, ::Tuple{}) = (false, ())

IteratorSize(::Type{Zip{Is}}) where {Is<:Tuple} = zip_iteratorsize(ntuple(n -> IteratorSize(fieldtype(Is, n)), _counttuple(Is)::Int)...)
IteratorEltype(::Type{Zip{Is}}) where {Is<:Tuple} = zip_iteratoreltype(ntuple(n -> IteratorEltype(fieldtype(Is, n)), _counttuple(Is)::Int)...)

zip_iteratorsize() = IsInfinite()
zip_iteratorsize(I) = I
zip_iteratorsize(a, b) = and_iteratorsize(a,b) # as `and_iteratorsize` but inherit `Union{HasLength,IsInfinite}` of the shorter iterator
zip_iteratorsize(::HasLength, ::IsInfinite) = HasLength()
zip_iteratorsize(::HasShape, ::IsInfinite) = HasLength()
zip_iteratorsize(a::IsInfinite, b) = zip_iteratorsize(b,a)
zip_iteratorsize(a::IsInfinite, b::IsInfinite) = IsInfinite()
zip_iteratorsize(a, b, tail...) = zip_iteratorsize(a, zip_iteratorsize(b, tail...))

zip_iteratoreltype() = HasEltype()
zip_iteratoreltype(a) = a
zip_iteratoreltype(a, tail...) = and_iteratoreltype(a, zip_iteratoreltype(tail...))

reverse(z::Zip) = Zip(Base.map(reverse, z.is)) # n.b. we assume all iterators are the same length
last(z::Zip) = getindex.(z.is, minimum(Base.map(lastindex, z.is)))

# filter

struct Filter{F,I}
    flt::F
    itr::I
end

"""
    Iterators.filter(flt, itr)

Given a predicate function `flt` and an iterable object `itr`, return an
iterable object which upon iteration yields the elements `x` of `itr` that
satisfy `flt(x)`. The order of the original iterator is preserved.

This function is *lazy*; that is, it is guaranteed to return in ``Θ(1)`` time
and use ``Θ(1)`` additional space, and `flt` will not be called by an
invocation of `filter`. Calls to `flt` will be made when iterating over the
returned iterable object. These calls are not cached and repeated calls will be
made when reiterating.

See [`Base.filter`](@ref) for an eager implementation of filtering for arrays.

# Examples
```jldoctest
julia> f = Iterators.filter(isodd, [1, 2, 3, 4, 5])
Base.Iterators.Filter{typeof(isodd), Vector{Int64}}(isodd, [1, 2, 3, 4, 5])

julia> foreach(println, f)
1
3
5

julia> [x for x in [1, 2, 3, 4, 5] if isodd(x)]  # collects a generator over Iterators.filter
3-element Vector{Int64}:
 1
 3
 5
```
"""
filter(flt, itr) = Filter(flt, itr)

function iterate(f::Filter, state...)
    y = iterate(f.itr, state...)
    while y !== nothing
        if f.flt(y[1])
            return y
        end
        y = iterate(f.itr, y[2])
    end
    nothing
end

eltype(::Type{Filter{F,I}}) where {F,I} = eltype(I)
IteratorEltype(::Type{Filter{F,I}}) where {F,I} = IteratorEltype(I)
IteratorSize(::Type{<:Filter}) = SizeUnknown()

reverse(f::Filter) = Filter(f.flt, reverse(f.itr))
last(f::Filter) = first(reverse(f))

# Accumulate -- partial reductions of a function over an iterator

struct Accumulate{F,I,T}
    f::F
    itr::I
    init::T
end

"""
    Iterators.accumulate(f, itr; [init])

Given a 2-argument function `f` and an iterator `itr`, return a new
iterator that successively applies `f` to the previous value and the
next element of `itr`.

This is effectively a lazy version of [`Base.accumulate`](@ref).

!!! compat "Julia 1.5"
    Keyword argument `init` is added in Julia 1.5.

# Examples
```jldoctest
julia> a = Iterators.accumulate(+, [1,2,3,4]);

julia> foreach(println, a)
1
3
6
10

julia> b = Iterators.accumulate(/, (2, 5, 2, 5); init = 100);

julia> collect(b)
4-element Vector{Float64}:
 50.0
 10.0
  5.0
  1.0
```
"""
accumulate(f, itr; init = Base._InitialValue()) = Accumulate(f, itr, init)

function iterate(itr::Accumulate)
    state = iterate(itr.itr)
    if state === nothing
        return nothing
    end
    val = Base.BottomRF(itr.f)(itr.init, state[1])
    return (val, (val, state[2]))
end

function iterate(itr::Accumulate, state)
    nxt = iterate(itr.itr, state[2])
    if nxt === nothing
        return nothing
    end
    val = itr.f(state[1], nxt[1])
    return (val, (val, nxt[2]))
end

length(itr::Accumulate) = length(itr.itr)
size(itr::Accumulate) = size(itr.itr)

IteratorSize(::Type{<:Accumulate{<:Any,I}}) where {I} = IteratorSize(I)
IteratorEltype(::Type{<:Accumulate}) = EltypeUnknown()

# Rest -- iterate starting at the given state

struct Rest{I,S}
    itr::I
    st::S
end

"""
    rest(iter, state)

An iterator that yields the same elements as `iter`, but starting at the given `state`.

See also: [`Iterators.drop`](@ref), [`Iterators.peel`](@ref), [`Base.rest`](@ref).

# Examples
```jldoctest
julia> collect(Iterators.rest([1,2,3,4], 2))
3-element Vector{Int64}:
 2
 3
 4
```
"""
rest(itr,state) = Rest(itr,state)
rest(itr::Rest,state) = Rest(itr.itr,state)
rest(itr) = itr

"""
    peel(iter)

Returns the first element and an iterator over the remaining elements.

If the iterator is empty return `nothing` (like `iterate`).

!!! compat "Julia 1.7"
    Prior versions throw a BoundsError if the iterator is empty.

See also: [`Iterators.drop`](@ref), [`Iterators.take`](@ref).

# Examples
```jldoctest
julia> (a, rest) = Iterators.peel("abc");

julia> a
'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)

julia> collect(rest)
2-element Vector{Char}:
 'b': ASCII/Unicode U+0062 (category Ll: Letter, lowercase)
 'c': ASCII/Unicode U+0063 (category Ll: Letter, lowercase)
```
"""
function peel(itr)
    y = iterate(itr)
    y === nothing && return y
    val, s = y
    val, rest(itr, s)
end

@propagate_inbounds iterate(i::Rest, st=i.st) = iterate(i.itr, st)
isdone(i::Rest, st...) = isdone(i.itr, st...)

eltype(::Type{<:Rest{I}}) where {I} = eltype(I)
IteratorEltype(::Type{<:Rest{I}}) where {I} = IteratorEltype(I)
rest_iteratorsize(a) = SizeUnknown()
rest_iteratorsize(::IsInfinite) = IsInfinite()
IteratorSize(::Type{<:Rest{I}}) where {I} = rest_iteratorsize(IteratorSize(I))

# Count -- infinite counting

struct Count{T,S}
    start::T
    step::S
end

"""
    countfrom(start=1, step=1)

An iterator that counts forever, starting at `start` and incrementing by `step`.

# Examples
```jldoctest
julia> for v in Iterators.countfrom(5, 2)
           v > 10 && break
           println(v)
       end
5
7
9
```
"""
countfrom(start::T, step::S) where {T,S} = Count{typeof(start+step),S}(start, step)
countfrom(start::Number, step::Number)   = Count(promote(start, step)...)
countfrom(start)                         = Count(start, oneunit(start))
countfrom()                              = Count(1, 1)

eltype(::Type{<:Count{T}}) where {T} = T

iterate(it::Count, state=it.start) = (state, state + it.step)

IteratorSize(::Type{<:Count}) = IsInfinite()

# Take -- iterate through the first n elements

struct Take{I}
    xs::I
    n::Int
    function Take(xs::I, n::Integer) where {I}
        n < 0 && throw(ArgumentError("Take length must be nonnegative"))
        return new{I}(xs, n)
    end
end

"""
    take(iter, n)

An iterator that generates at most the first `n` elements of `iter`.

See also: [`drop`](@ref Iterators.drop), [`peel`](@ref Iterators.peel), [`first`](@ref), [`take!`](@ref).

# Examples
```jldoctest
julia> a = 1:2:11
1:2:11

julia> collect(a)
6-element Vector{Int64}:
  1
  3
  5
  7
  9
 11

julia> collect(Iterators.take(a,3))
3-element Vector{Int64}:
 1
 3
 5
```
"""
take(xs, n::Integer) = Take(xs, Int(n))
take(xs::Take, n::Integer) = Take(xs.xs, min(Int(n), xs.n))

eltype(::Type{Take{I}}) where {I} = eltype(I)
IteratorEltype(::Type{Take{I}}) where {I} = IteratorEltype(I)
take_iteratorsize(a) = HasLength()
take_iteratorsize(::SizeUnknown) = SizeUnknown()
IteratorSize(::Type{Take{I}}) where {I} = take_iteratorsize(IteratorSize(I))
length(t::Take) = _min_length(t.xs, 1:t.n, IteratorSize(t.xs), HasLength())
isdone(t::Take) = isdone(t.xs)
isdone(t::Take, state) = (state[1] <= 0) | isdone(t.xs, tail(state))

@propagate_inbounds function iterate(it::Take, state=(it.n,))
    n, rest = state[1], tail(state)
    n <= 0 && return nothing
    y = iterate(it.xs, rest...)
    y === nothing && return nothing
    return y[1], (n - 1, y[2])
end

# Drop -- iterator through all but the first n elements

struct Drop{I}
    xs::I
    n::Int
    function Drop(xs::I, n::Integer) where {I}
        n < 0 && throw(ArgumentError("Drop length must be nonnegative"))
        return new{I}(xs, n)
    end
end

"""
    drop(iter, n)

An iterator that generates all but the first `n` elements of `iter`.

# Examples
```jldoctest
julia> a = 1:2:11
1:2:11

julia> collect(a)
6-element Vector{Int64}:
  1
  3
  5
  7
  9
 11

julia> collect(Iterators.drop(a,4))
2-element Vector{Int64}:
  9
 11
```
"""
drop(xs, n::Integer) = Drop(xs, Int(n))
drop(xs::Take, n::Integer) = Take(drop(xs.xs, Int(n)), max(0, xs.n - Int(n)))
drop(xs::Drop, n::Integer) = Drop(xs.xs, Int(n) + xs.n)

eltype(::Type{Drop{I}}) where {I} = eltype(I)
IteratorEltype(::Type{Drop{I}}) where {I} = IteratorEltype(I)
drop_iteratorsize(::SizeUnknown) = SizeUnknown()
drop_iteratorsize(::Union{HasShape, HasLength}) = HasLength()
drop_iteratorsize(::IsInfinite) = IsInfinite()
IteratorSize(::Type{Drop{I}}) where {I} = drop_iteratorsize(IteratorSize(I))
length(d::Drop) = _diff_length(d.xs, 1:d.n, IteratorSize(d.xs), HasLength())

function iterate(it::Drop)
    y = iterate(it.xs)
    for i in 1:it.n
        y === nothing && return y
        y = iterate(it.xs, y[2])
    end
    y
end
iterate(it::Drop, state) = iterate(it.xs, state)
isdone(it::Drop, state) = isdone(it.xs, state)

# takewhile

struct TakeWhile{I,P<:Function}
    pred::P
    xs::I
end

"""
    takewhile(pred, iter)

An iterator that generates element from `iter` as long as predicate `pred` is true,
afterwards, drops every element.

!!! compat "Julia 1.4"
    This function requires at least Julia 1.4.

# Examples

```jldoctest
julia> s = collect(1:5)
5-element Vector{Int64}:
 1
 2
 3
 4
 5

julia> collect(Iterators.takewhile(<(3),s))
2-element Vector{Int64}:
 1
 2
```
"""
takewhile(pred,xs) = TakeWhile(pred,xs)

function iterate(ibl::TakeWhile, itr...)
    y = iterate(ibl.xs,itr...)
    y === nothing && return nothing
    ibl.pred(y[1]) || return nothing
    y
end

IteratorSize(::Type{<:TakeWhile}) = SizeUnknown()
eltype(::Type{TakeWhile{I,P}} where P) where {I} = eltype(I)
IteratorEltype(::Type{TakeWhile{I, P}} where P) where {I} = IteratorEltype(I)

# dropwhile

struct DropWhile{I,P<:Function}
    pred::P
    xs::I
end

"""
    dropwhile(pred, iter)

An iterator that drops element from `iter` as long as predicate `pred` is true,
afterwards, returns every element.

!!! compat "Julia 1.4"
    This function requires at least Julia 1.4.

# Examples

```jldoctest
julia> s = collect(1:5)
5-element Vector{Int64}:
 1
 2
 3
 4
 5

julia> collect(Iterators.dropwhile(<(3),s))
3-element Vector{Int64}:
 3
 4
 5
```
"""
dropwhile(pred,itr) = DropWhile(pred,itr)

iterate(ibl::DropWhile,itr) = iterate(ibl.xs, itr)
function iterate(ibl::DropWhile)
    y = iterate(ibl.xs)
    while y !== nothing
        ibl.pred(y[1]) || break
        y = iterate(ibl.xs,y[2])
    end
    y
end

IteratorSize(::Type{<:DropWhile}) = SizeUnknown()
eltype(::Type{DropWhile{I,P}}) where {I,P} = eltype(I)
IteratorEltype(::Type{DropWhile{I,P}}) where {I,P} = IteratorEltype(I)

# Cycle an iterator forever

struct Cycle{I}
    xs::I
end

"""
    cycle(iter)

An iterator that cycles through `iter` forever.
If `iter` is empty, so is `cycle(iter)`.

See also: [`Iterators.repeated`](@ref), [`repeat`](@ref).

# Examples
```jldoctest
julia> for (i, v) in enumerate(Iterators.cycle("hello"))
           print(v)
           i > 10 && break
       end
hellohelloh
```
"""
cycle(xs) = Cycle(xs)

eltype(::Type{Cycle{I}}) where {I} = eltype(I)
IteratorEltype(::Type{Cycle{I}}) where {I} = IteratorEltype(I)
IteratorSize(::Type{Cycle{I}}) where {I} = IsInfinite()

iterate(it::Cycle) = iterate(it.xs)
isdone(it::Cycle) = isdone(it.xs)
isdone(it::Cycle, state) = false
function iterate(it::Cycle, state)
    y = iterate(it.xs, state)
    y === nothing && return iterate(it)
    y
end

reverse(it::Cycle) = Cycle(reverse(it.xs))
last(it::Cycle) = last(it.xs)

# Repeated - repeat an object infinitely many times

struct Repeated{O}
    x::O
end
repeated(x) = Repeated(x)

"""
    repeated(x[, n::Int])

An iterator that generates the value `x` forever. If `n` is specified, generates `x` that
many times (equivalent to `take(repeated(x), n)`).

See also: [`Iterators.cycle`](@ref), [`repeat`](@ref).

# Examples
```jldoctest
julia> a = Iterators.repeated([1 2], 4);

julia> collect(a)
4-element Vector{Matrix{Int64}}:
 [1 2]
 [1 2]
 [1 2]
 [1 2]
```
"""
repeated(x, n::Integer) = take(repeated(x), Int(n))

eltype(::Type{Repeated{O}}) where {O} = O

iterate(it::Repeated, state...) = (it.x, nothing)

IteratorSize(::Type{<:Repeated}) = IsInfinite()
IteratorEltype(::Type{<:Repeated}) = HasEltype()

reverse(it::Union{Repeated,Take{<:Repeated}}) = it
last(it::Union{Repeated,Take{<:Repeated}}) = first(it)

# Product -- cartesian product of iterators
struct ProductIterator{T<:Tuple}
    iterators::T
end

"""
    product(iters...)

Return an iterator over the product of several iterators. Each generated element is
a tuple whose `i`th element comes from the `i`th argument iterator. The first iterator
changes the fastest.

See also: [`zip`](@ref), [`Iterators.flatten`](@ref).

# Examples
```jldoctest
julia> collect(Iterators.product(1:2, 3:5))
2×3 Matrix{Tuple{Int64, Int64}}:
 (1, 3)  (1, 4)  (1, 5)
 (2, 3)  (2, 4)  (2, 5)

julia> ans == [(x,y) for x in 1:2, y in 3:5]  # collects a generator involving Iterators.product
true
```
"""
product(iters...) = ProductIterator(iters)

IteratorSize(::Type{ProductIterator{Tuple{}}}) = HasShape{0}()
IteratorSize(::Type{ProductIterator{T}}) where {T<:Tuple} =
    prod_iteratorsize(ntuple(n -> IteratorSize(fieldtype(T, n)), _counttuple(T)::Int)..., HasShape{0}())

prod_iteratorsize() = HasShape{0}()
prod_iteratorsize(I) = I

prod_iteratorsize(::HasLength, ::HasLength) = HasShape{2}()
prod_iteratorsize(::HasLength, ::HasShape{N}) where {N} = HasShape{N+1}()
prod_iteratorsize(::HasShape{N}, ::HasLength) where {N} = HasShape{N+1}()
prod_iteratorsize(::HasShape{M}, ::HasShape{N}) where {M,N} = HasShape{M+N}()

# products can have an infinite iterator
prod_iteratorsize(::IsInfinite, ::IsInfinite) = IsInfinite()
prod_iteratorsize(a, ::IsInfinite) = IsInfinite()
prod_iteratorsize(::IsInfinite, b) = IsInfinite()
prod_iteratorsize(a, b) = SizeUnknown()
prod_iteratorsize(a, b, tail...) = prod_iteratorsize(a, prod_iteratorsize(b, tail...))

size(P::ProductIterator) = _prod_size(P.iterators)
_prod_size(::Tuple{}) = ()
_prod_size(t::Tuple) = (_prod_size1(t[1], IteratorSize(t[1]))..., _prod_size(tail(t))...)
_prod_size1(a, ::HasShape)  = size(a)
_prod_size1(a, ::HasLength) = (length(a),)
_prod_size1(a, A) =
    throw(ArgumentError("Cannot compute size for object of type $(typeof(a))"))

axes(P::ProductIterator) = _prod_indices(P.iterators)
_prod_indices(::Tuple{}) = ()
_prod_indices(t::Tuple) = (_prod_axes1(t[1], IteratorSize(t[1]))..., _prod_indices(tail(t))...)
_prod_axes1(a, ::HasShape)  = axes(a)
_prod_axes1(a, ::HasLength) = (OneTo(length(a)),)
_prod_axes1(a, A) =
    throw(ArgumentError("Cannot compute indices for object of type $(typeof(a))"))

ndims(p::ProductIterator) = length(axes(p))
length(P::ProductIterator) = prod(size(P))

IteratorEltype(::Type{ProductIterator{Tuple{}}}) = HasEltype()
IteratorEltype(::Type{ProductIterator{Tuple{I}}}) where {I} = IteratorEltype(I)

function IteratorEltype(::Type{ProductIterator{T}}) where {T<:Tuple}
    E = ntuple(n -> IteratorEltype(fieldtype(T, n)), _counttuple(T)::Int)
    any(I -> I == EltypeUnknown(), E) && return EltypeUnknown()
    return E[end]
end

eltype(::Type{ProductIterator{I}}) where {I} = _prod_eltype(I)
_prod_eltype(::Type{Tuple{}}) = Tuple{}
_prod_eltype(::Type{I}) where {I<:Tuple} =
    Tuple{ntuple(n -> eltype(fieldtype(I, n)), _counttuple(I)::Int)...}

iterate(::ProductIterator{Tuple{}}) = (), true
iterate(::ProductIterator{Tuple{}}, state) = nothing

@inline isdone(P::ProductIterator) = any(isdone, P.iterators)
@inline function _pisdone(iters, states)
    iter1 = first(iters)
    done1 = isdone(iter1, first(states)[2]) # check step
    done1 === true || return done1 # false or missing
    done1 = isdone(iter1) # check restart
    done1 === true || return done1 # false or missing
    return _pisdone(tail(iters), tail(states)) # check tail
end
@inline isdone(P::ProductIterator, states) = _pisdone(P.iterators, states)

@inline _piterate() = ()
@inline function _piterate(iter1, rest...)
    next = iterate(iter1)
    next === nothing && return nothing
    restnext = _piterate(rest...)
    restnext === nothing && return nothing
    return (next, restnext...)
end
@inline function iterate(P::ProductIterator)
    isdone(P) === true && return nothing
    next = _piterate(P.iterators...)
    next === nothing && return nothing
    return (Base.map(x -> x[1], next), next)
end

@inline _piterate1(::Tuple{}, ::Tuple{}) = nothing
@inline function _piterate1(iters, states)
    iter1 = first(iters)
    next = iterate(iter1, first(states)[2])
    restnext = tail(states)
    if next === nothing
        isdone(iter1) === true && return nothing
        restnext = _piterate1(tail(iters), restnext)
        restnext === nothing && return nothing
        next = iterate(iter1)
        next === nothing && return nothing
    end
    return (next, restnext...)
end
@inline function iterate(P::ProductIterator, states)
    isdone(P, states) === true && return nothing
    next = _piterate1(P.iterators, states)
    next === nothing && return nothing
    return (Base.map(x -> x[1], next), next)
end

reverse(p::ProductIterator) = ProductIterator(Base.map(reverse, p.iterators))
last(p::ProductIterator) = Base.map(last, p.iterators)

# flatten an iterator of iterators

struct Flatten{I}
    it::I
end

"""
    flatten(iter)

Given an iterator that yields iterators, return an iterator that yields the
elements of those iterators.
Put differently, the elements of the argument iterator are concatenated.

# Examples
```jldoctest
julia> collect(Iterators.flatten((1:2, 8:9)))
4-element Vector{Int64}:
 1
 2
 8
 9

julia> [(x,y) for x in 0:1 for y in 'a':'c']  # collects generators involving Iterators.flatten
6-element Vector{Tuple{Int64, Char}}:
 (0, 'a')
 (0, 'b')
 (0, 'c')
 (1, 'a')
 (1, 'b')
 (1, 'c')
```
"""
flatten(itr) = Flatten(itr)

eltype(::Type{Flatten{I}}) where {I} = eltype(eltype(I))
eltype(::Type{Flatten{Tuple{}}}) = eltype(Tuple{})
IteratorEltype(::Type{Flatten{I}}) where {I} = _flatteneltype(I, IteratorEltype(I))
IteratorEltype(::Type{Flatten{Tuple{}}}) = IteratorEltype(Tuple{})
_flatteneltype(I, ::HasEltype) = IteratorEltype(eltype(I))
_flatteneltype(I, et) = EltypeUnknown()

flatten_iteratorsize(::Union{HasShape, HasLength}, ::Type{<:NTuple{N,Any}}) where {N} = HasLength()
flatten_iteratorsize(::Union{HasShape, HasLength}, ::Type{<:Tuple}) = SizeUnknown()
flatten_iteratorsize(::Union{HasShape, HasLength}, ::Type{<:Number}) = HasLength()
flatten_iteratorsize(a, b) = SizeUnknown()

_flatten_iteratorsize(sz, ::EltypeUnknown, I) = SizeUnknown()
_flatten_iteratorsize(sz, ::HasEltype, I) = flatten_iteratorsize(sz, eltype(I))
_flatten_iteratorsize(sz, ::HasEltype, ::Type{Tuple{}}) = HasLength()

IteratorSize(::Type{Flatten{I}}) where {I} = _flatten_iteratorsize(IteratorSize(I), IteratorEltype(I), I)

function flatten_length(f, T::Type{<:NTuple{N,Any}}) where {N}
    return N * length(f.it)
end
flatten_length(f, ::Type{<:Number}) = length(f.it)
flatten_length(f, T) = throw(ArgumentError(
    "Iterates of the argument to Flatten are not known to have constant length"))
length(f::Flatten{I}) where {I} = flatten_length(f, eltype(I))
length(f::Flatten{Tuple{}}) = 0

@propagate_inbounds function iterate(f::Flatten, state=())
    if state !== ()
        y = iterate(tail(state)...)
        y !== nothing && return (y[1], (state[1], state[2], y[2]))
    end
    x = (state === () ? iterate(f.it) : iterate(f.it, state[1]))
    x === nothing && return nothing
    y = iterate(x[1])
    while y === nothing
         x = iterate(f.it, x[2])
         x === nothing && return nothing
         y = iterate(x[1])
    end
    return y[1], (x[2], x[1], y[2])
end

reverse(f::Flatten) = Flatten(reverse(itr) for itr in reverse(f.it))
last(f::Flatten) = last(last(f.it))

"""
    Iterators.flatmap(f, iterators...)

Equivalent to `flatten(map(f, iterators...))`.

# Examples
```jldoctest
julia> Iterators.flatmap(n->-n:2:n, 1:3) |> collect
9-element Vector{Int64}:
 -1
  1
 -2
  0
  2
 -3
 -1
  1
  3
```
"""
# flatmap = flatten ∘ map
flatmap(f, c...) = flatten(map(f, c...))

"""
    partition(collection, n)

Iterate over a collection `n` elements at a time.

# Examples
```jldoctest
julia> collect(Iterators.partition([1,2,3,4,5], 2))
3-element Vector{SubArray{Int64, 1, Vector{Int64}, Tuple{UnitRange{Int64}}, true}}:
 [1, 2]
 [3, 4]
 [5]
```
"""
function partition(c, n::Integer)
    n < 1 && throw(ArgumentError("cannot create partitions of length $n"))
    return PartitionIterator(c, Int(n))
end

struct PartitionIterator{T}
    c::T
    n::Int
end
# Partitions are explicitly a linear indexing operation, so reshape to 1-d immediately
PartitionIterator(A::AbstractArray, n::Int) = PartitionIterator(vec(A), n)
PartitionIterator(v::AbstractVector, n::Int) = PartitionIterator{typeof(v)}(v, n)

eltype(::Type{PartitionIterator{T}}) where {T} = Vector{eltype(T)}
# Arrays use a generic `view`-of-a-`vec`, so we cannot exactly predict what we'll get back
eltype(::Type{PartitionIterator{T}}) where {T<:AbstractArray} = AbstractVector{eltype(T)}
# But for some common implementations in Base we know the answer exactly
eltype(::Type{PartitionIterator{T}}) where {T<:Vector} = SubArray{eltype(T), 1, T, Tuple{UnitRange{Int}}, true}

IteratorEltype(::Type{PartitionIterator{T}}) where {T} = IteratorEltype(T)
IteratorEltype(::Type{PartitionIterator{T}}) where {T<:AbstractArray} = EltypeUnknown()
IteratorEltype(::Type{PartitionIterator{T}}) where {T<:Vector} = IteratorEltype(T)

partition_iteratorsize(::HasShape) = HasLength()
partition_iteratorsize(isz) = isz
function IteratorSize(::Type{PartitionIterator{T}}) where {T}
    partition_iteratorsize(IteratorSize(T))
end

function length(itr::PartitionIterator)
    l = length(itr.c)
    return cld(l, itr.n)
end

function iterate(itr::PartitionIterator{<:AbstractRange}, state = firstindex(itr.c))
    state > lastindex(itr.c) && return nothing
    r = min(state + itr.n - 1, lastindex(itr.c))
    return @inbounds itr.c[state:r], r + 1
end

function iterate(itr::PartitionIterator{<:AbstractArray}, state = firstindex(itr.c))
    state > lastindex(itr.c) && return nothing
    r = min(state + itr.n - 1, lastindex(itr.c))
    return @inbounds view(itr.c, state:r), r + 1
end

struct IterationCutShort; end

function iterate(itr::PartitionIterator, state...)
    # This is necessary to remember whether we cut the
    # last element short. In such cases, we do return that
    # element, but not the next one
    state === (IterationCutShort(),) && return nothing
    v = Vector{eltype(itr.c)}(undef, itr.n)
    i = 0
    y = iterate(itr.c, state...)
    while y !== nothing
        i += 1
        v[i] = y[1]
        if i >= itr.n
            break
        end
        y = iterate(itr.c, y[2])
    end
    i === 0 && return nothing
    return resize!(v, i), y === nothing ? IterationCutShort() : y[2]
end

"""
    Stateful(itr)

There are several different ways to think about this iterator wrapper:

1. It provides a mutable wrapper around an iterator and
   its iteration state.
2. It turns an iterator-like abstraction into a `Channel`-like
   abstraction.
3. It's an iterator that mutates to become its own rest iterator
   whenever an item is produced.

`Stateful` provides the regular iterator interface. Like other mutable iterators
(e.g. [`Channel`](@ref)), if iteration is stopped early (e.g. by a [`break`](@ref) in a [`for`](@ref) loop),
iteration can be resumed from the same spot by continuing to iterate over the
same iterator object (in contrast, an immutable iterator would restart from the
beginning).

# Examples
```jldoctest
julia> a = Iterators.Stateful("abcdef");

julia> isempty(a)
false

julia> popfirst!(a)
'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)

julia> collect(Iterators.take(a, 3))
3-element Vector{Char}:
 'b': ASCII/Unicode U+0062 (category Ll: Letter, lowercase)
 'c': ASCII/Unicode U+0063 (category Ll: Letter, lowercase)
 'd': ASCII/Unicode U+0064 (category Ll: Letter, lowercase)

julia> collect(a)
2-element Vector{Char}:
 'e': ASCII/Unicode U+0065 (category Ll: Letter, lowercase)
 'f': ASCII/Unicode U+0066 (category Ll: Letter, lowercase)

julia> Iterators.reset!(a); popfirst!(a)
'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)

julia> Iterators.reset!(a, "hello"); popfirst!(a)
'h': ASCII/Unicode U+0068 (category Ll: Letter, lowercase)
```

```jldoctest
julia> a = Iterators.Stateful([1,1,1,2,3,4]);

julia> for x in a; x == 1 || break; end

julia> peek(a)
3

julia> sum(a) # Sum the remaining elements
7
```
"""
mutable struct Stateful{T, VS}
    itr::T
    # A bit awkward right now, but adapted to the new iteration protocol
    nextvalstate::Union{VS, Nothing}
    taken::Int
    @inline function Stateful{<:Any, Any}(itr::T) where {T}
        new{T, Any}(itr, iterate(itr), 0)
    end
    @inline function Stateful(itr::T) where {T}
        VS = approx_iter_type(T)
        return new{T, VS}(itr, iterate(itr)::VS, 0)
    end
end

function reset!(s::Stateful{T,VS}, itr::T=s.itr) where {T,VS}
    s.itr = itr
    setfield!(s, :nextvalstate, iterate(itr))
    s.taken = 0
    s
end

if Base === Core.Compiler
    approx_iter_type(a::Type) = Any
else
    # Try to find an appropriate type for the (value, state tuple),
    # by doing a recursive unrolling of the iteration protocol up to
    # fixpoint.
    approx_iter_type(itrT::Type) = _approx_iter_type(itrT, Base._return_type(iterate, Tuple{itrT}))
    # Not actually called, just passed to return type to avoid
    # having to typesplit on Nothing
    function doiterate(itr, valstate::Union{Nothing, Tuple{Any, Any}})
        valstate === nothing && return nothing
        val, st = valstate
        return iterate(itr, st)
    end
    function _approx_iter_type(itrT::Type, vstate::Type)
        vstate <: Union{Nothing, Tuple{Any, Any}} || return Any
        vstate <: Union{} && return Union{}
        nextvstate = Base._return_type(doiterate, Tuple{itrT, vstate})
        return (nextvstate <: vstate ? vstate : Any)
    end
end

convert(::Type{Stateful}, itr) = Stateful(itr)

@inline isdone(s::Stateful, st=nothing) = s.nextvalstate === nothing

@inline function popfirst!(s::Stateful)
    vs = s.nextvalstate
    if vs === nothing
        throw(EOFError())
    else
        val, state = vs
        Core.setfield!(s, :nextvalstate, iterate(s.itr, state))
        s.taken += 1
        return val
    end
end

@inline function peek(s::Stateful, sentinel=nothing)
    ns = s.nextvalstate
    return ns !== nothing ? ns[1] : sentinel
end
@inline iterate(s::Stateful, state=nothing) = s.nextvalstate === nothing ? nothing : (popfirst!(s), nothing)
IteratorSize(::Type{Stateful{T,VS}}) where {T,VS} = IteratorSize(T) isa HasShape ? HasLength() : IteratorSize(T)
eltype(::Type{Stateful{T, VS}} where VS) where {T} = eltype(T)
IteratorEltype(::Type{Stateful{T,VS}}) where {T,VS} = IteratorEltype(T)
length(s::Stateful) = length(s.itr) - s.taken

"""
    only(x)

Return the one and only element of collection `x`, or throw an [`ArgumentError`](@ref) if the
collection has zero or multiple elements.

See also [`first`](@ref), [`last`](@ref).

!!! compat "Julia 1.4"
    This method requires at least Julia 1.4.

# Examples
```jldoctest
julia> only(["a"])
"a"

julia> only("a")
'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)

julia> only(())
ERROR: ArgumentError: Tuple contains 0 elements, must contain exactly 1 element
Stacktrace:
[...]

julia> only(('a', 'b'))
ERROR: ArgumentError: Tuple contains 2 elements, must contain exactly 1 element
Stacktrace:
[...]
```
"""
@propagate_inbounds function only(x)
    i = iterate(x)
    @boundscheck if i === nothing
        throw(ArgumentError("Collection is empty, must contain exactly 1 element"))
    end
    (ret, state) = i::NTuple{2,Any}
    @boundscheck if iterate(x, state) !== nothing
        throw(ArgumentError("Collection has multiple elements, must contain exactly 1 element"))
    end
    return ret
end

# Collections of known size
only(x::Ref) = x[]
only(x::Number) = x
only(x::Char) = x
only(x::Tuple{Any}) = x[1]
only(x::Tuple) = throw(
    ArgumentError("Tuple contains $(length(x)) elements, must contain exactly 1 element")
)
only(a::AbstractArray{<:Any, 0}) = @inbounds return a[]
only(x::NamedTuple{<:Any, <:Tuple{Any}}) = first(x)
only(x::NamedTuple) = throw(
    ArgumentError("NamedTuple contains $(length(x)) elements, must contain exactly 1 element")
)

Base.intersect(a::ProductIterator, b::ProductIterator) = ProductIterator(intersect.(a.iterators, b.iterators))

end

max_stars_repo_name: collinwarner/julia

max_issues_repo_name: collinwarner/julia

max_forks_repo_name: collinwarner/julia


text:
using Mimi

@defcomp sealevelrise begin
    regions = Index()

    SLRTHERM = Variable(index=[time])
    THERMEQUIL = Variable(index=[time])
    GSICREMAIN = Variable(index=[time])
    GSICCUM = Variable(index=[time])
    GSICMELTRATE = Variable(index=[time])
    GISREMAIN = Variable(index=[time])
    GISMELTRATE = Variable(index=[time])
    GISEXPONENT = Variable(index=[time])
    GISCUM = Variable(index=[time])
    AISREMAIN = Variable(index=[time])
    AISMELTRATE = Variable(index=[time])
    AISCUM = Variable(index=[time])
    TOTALSLR = Variable(index=[time])

    TATM = Parameter(index=[time])
    thermeq = Parameter()
    therm0 = Parameter()
    thermadj = Parameter()
    gsictotal = Parameter()
    gsicmelt = Parameter()
    gsicexp = Parameter()
    gis0 = Parameter()
    gismelt0 = Parameter()
    gismeltabove = Parameter()
    gismineq = Parameter()
    gisexp = Parameter()
    aismelt0 = Parameter()
    aismeltlow = Parameter()
    aismeltup = Parameter()
    aisratio = Parameter()
    aisinflection = Parameter()
    aisintercept = Parameter()
    aiswais = Parameter()
    aisother = Parameter()

    function run_timestep(p, v, d, t)

    #THERMAL EXPANSION

        #Define function for THERMEQUIL
        v.THERMEQUIL[t] = p.TATM[t] * p.thermeq

            #Define function for SLRTHERM
            if is_first(t)
                v.SLRTHERM[t] = p.therm0 + p.thermadj * (v.THERMEQUIL[t] - p.therm0)
            else
                v.SLRTHERM[t] = v.SLRTHERM[t-1] + p.thermadj * (v.THERMEQUIL[t] - v.SLRTHERM[t-1])
            end

        #GLACIERS AND SMALL ICE CAPS (GSIC)

            #Define function for GSICREMAIN
            if is_first(t)
                v.GSICREMAIN[t] = p.gsictotal
            else
                v.GSICREMAIN[t] = p.gsictotal - v.GSICCUM[t-1]
            end

            #Define function for GSICMELTRATE
            if is_first(t)
                v.GSICMELTRATE[t] = 0.01464  # #Bug in Excel Model?  gsiceq parameter drops out after first period (this is replicated here by just setting period 1 and exluding gsiceq from equation)
            else
                v.GSICMELTRATE[t] = p.gsicmelt * 10 * (v.GSICREMAIN[t] / p.gsictotal)^(p.gsicexp) * p.TATM[t]
            end

            #Define function for GSICCUM
            if is_first(t)
                v.GSICCUM[t] = v.GSICMELTRATE[t]
            else
                v.GSICCUM[t] = v.GSICCUM[t-1] + v.GSICMELTRATE[t]
            end

        #GREENLAND ICE SHEETS (GIS)

            #Define function for GISREMAIN
            if is_first(t)
                v.GISREMAIN[t] = p.gis0
            else
                v.GISREMAIN[t] = v.GISREMAIN[t-1] - v.GISMELTRATE[t-1] / 100
            end

            #Define function for GISMELTRATE
            if is_first(t) || t.t == 2
                v.GISMELTRATE[t] = p.gismelt0
            else
                v.GISMELTRATE[t] = (p.gismeltabove * (p.TATM[t] - p.gismineq) + p.gismelt0) * v.GISEXPONENT[t-1]
            end

            #Define function for GISCUM
            if is_first(t)
                v.GISCUM[t] = p.gismelt0 / 100
            else
                v.GISCUM[t] = v.GISCUM[t-1] + v.GISMELTRATE[t] / 100
            end

            #Define function for GISEXPONENT
            if is_first(t) || t.t == 2
                v.GISEXPONENT[t] = 1.
            else
                v.GISEXPONENT[t] = 1. - (v.GISCUM[t] / p.gis0)^p.gisexp
            end

        #ANTARCTIC ICE SHEET (AIS)

            #Define function for AISMELTRATE
            if t.t <= 11
            v.AISMELTRATE[t] = ifelse(p.TATM[t] < 3., (p.aismeltlow * p.TATM[t] * p.aisratio + p.aisintercept), (p.aisinflection * p.aismeltlow + p.aismeltup * (p.TATM[t] - 3.) + p.aisintercept))
            else
            v.AISMELTRATE[t] = ifelse(p.TATM[t] < 3., (p.aismeltlow * p.TATM[t] * p.aisratio + p.aismelt0), (p.aisinflection * p.aismeltlow + p.aismeltup * (p.TATM[t] - 3.) + p.aismelt0))
            end

            #Define function for AISCUM
            if is_first(t)
                v.AISCUM[t] = v.AISMELTRATE[t] / 100
            else
                v.AISCUM[t] = v.AISCUM[t-1] + v.AISMELTRATE[t] / 100
            end

            #Define function for AISREMAIN
            if is_first(t)
                v.AISREMAIN[t] = p.aiswais + p.aisother
            else
                v.AISREMAIN[t] = v.AISREMAIN[1] - v.AISCUM[t]
            end

        #TOTAL SEA LEVEL RISE AND DAMAGES

            #Define function for TOTALSLR
            v.TOTALSLR[t] = v.SLRTHERM[t] + v.GSICCUM[t] + v.GISCUM[t] + v.AISCUM[t]

    end
end


text:
module QuantumNLDiffEq

import Yao: AbstractBlock, zero_state, expect, dispatch!, dispatch, chain
import Flux: update!, ADAM, mse
import Zygote: Buffer, gradient, ignore

abstract type FeatureMap end
abstract type BoundaryHandling end

struct Product <: FeatureMap end
struct ChebyshevSparse <: FeatureMap end
struct ChebyshevTower <: FeatureMap end
struct Pinned <: BoundaryHandling end
struct Floating <: BoundaryHandling end

include("phi.jl")
include("new_circuit.jl")
include("calculate_diff_evalue.jl")

loss_calc(diff_evalue, evalue, x) = mse(diff_evalue + 8*evalue*(0.1 + tan(8*x)), 0)

function loss(feature_map_circuit::AbstractBlock, cost::AbstractBlock, var::AbstractBlock, M::AbstractVector, N::Int, u_0::Real, theta, pc; mapping::FeatureMap)
	Loss_diff = 0.0
	for x in 1:length(M)
		evalue = expect(cost, zero_state(N)=>new_circuit(feature_map_circuit, var, M[x], theta, N, pc; mapping)) - expect(cost, zero_state(N)=>new_circuit(feature_map_circuit, var, M[1], theta, N, pc; mapping)) + u_0

		diff_evalue = calculate_diff_evalue(feature_map_circuit, cost, var, N, theta, M[x], pc; mapping)

		Loss_diff += loss_calc(diff_evalue, evalue, M[x])
	end
	return real.(Loss_diff/length(M))
end

function train!(feature_map_circuit::AbstractBlock, cost::AbstractBlock, var::AbstractBlock, M::AbstractVector, N::Int, u_0::Real, theta, pc=2; mapping::FeatureMap=ChebyshevSparse(), boundary::BoundaryHandling=Floating(), optimizer=ADAM(0.05), steps=300)
	for _ in 1:steps
		grads = gradient(_theta->loss(feature_map_circuit, cost, var, M, N, u_0, _theta, pc; mapping), theta)[1]
		update!(optimizer, theta, grads)
		dispatch!(var, theta)
	end
end

export loss, train!

end


text:
using ArgParse
using LoggingExtras
using TerminalLoggers

function parse_cmdline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "mcmc"
            help = "Run MCMC to optimize hyperparameters"
            action = :command
        "predict"
            help = "Calculate the posterior of a GP at tdata given data"
            action = :command
        "sample"
            help = "Sample the posterior of a GP"
            action = :command
        "fitplot"
            help = "Diagnostic plots showing the posteriors of different components of the GP"
            action = :command
        "select"
            help = "Output model selection parameters"
            action = :command
        "--verbose", "-v"
            help = "Log level to @info"
            action = :store_true
        "--quiet", "-q"
            help = "Log level to @warning"
            action = :store_true
        "--debug"
            help = "Log level to @debug"
            action = :store_true
        "--log"
            help = "Log to a file as well as stdout"

    end

    @add_arg_table! s["sample"] begin
        "formula"
            help = "GP formula specification"
            required = true
        "--at"
            help = "Range to sample at, eg 'x=-5:0.1:5"
            required = true
        "--plot"
            help = "File to plot to"
        "--output", "-o"
            help = "Table output of GP sample - must end with '.csv' or '.tsv'"
    end

    @add_arg_table! s["mcmc"] begin
        "formula"
            help = "GP formula specification"
            required = true
        "--data", "-i"
            required = true
            help = """
                Table input on which to run inference.
                Must contain columns that correspond to values in 'formula'
                """
        "--infer"
            help = """
                Which model hyperparameter to infer.
                Specify variable names, the hyperparameter(s) will be determined
                based on kernel type (eg length scale for SExp)
                """
            nargs = '+'
            required = true
        "--samples"
            help = """
            Number of samples to take from the chain.
            Default=200
            """
            arg_type=Int
            default=200
        "--output", "-o"
            help = "Table to output sampling chain"
        "--plot"
            help = "File to plot to"
    end

    @add_arg_table! s["select"] begin
        "--formulae"
            help = """
                Compare 2 GP formula specifications, requires '--data' as well.
                Result will be logpdf of formula 2 - logpdf of formula 1.

                A positive value indicates more evidence for formula 2.
                """
            nargs = 2
        "--chains"
            help = """
                Compare 2 sampling chains from 'mcmc' command.
                Result will be the log2 bayes factor.

                A positive value indicates more evidence for chain 1.
                """
            nargs = 2
        "--data", "-i"
            help = """
                Table input on which to run inference.
                Must contain columns that correspond to values in both 'formulae'
                """
        "--plot"
            help = "File to plot to"
    end

    return parse_args(s)
end

# Main
args = parse_cmdline()

function setup_logs!(loglevel, logpath; dryrun=false)
    glog = TerminalLogger(stderr, loglevel)
    if logpath === nothing || dryrun
        global_logger(glog)
    else
        logpath = abspath(expanduser(logpath))
        global_logger(
            TeeLogger(
                MinLevelLogger(FileLogger(logpath), loglevel),
                glog))
    end
end

if args["debug"]
    loglevel = Logging.Debug
elseif args["verbose"]
    loglevel = Logging.Info
elseif args["quiet"]
    loglevel = Logging.Error
else
    loglevel = Logging.Warn
end

setup_logs!(loglevel, args["log"])
@info "Getting started!"

include("cli/GaPLAC_CLI.jl")
using .CLI

args["%COMMAND%"] == "sample" && CLI.Sample.run(args["sample"])
args["%COMMAND%"] == "mcmc" && CLI.MCMC.run(args["mcmc"])
args["%COMMAND%"] == "select" && CLI.Select.run(args["select"])

nothing
# # Redirect logging to a file if necessary
# loglevel = args["debug"] ? Logging.Debug : Logging.Info
# if !isa(args["log"], Nothing)
#     io = open(args["log"], "w+")
#     logger = SimpleLogger(io, loglevel)
#     global_logger(logger)
# else
#     logger = ConsoleLogger(stderr, loglevel; show_limited=false)
#     global_logger(logger)
# end

# # Process commands

# # Clean up
# if !isa(args["log"], Nothing)
#     flush(io)
#     close(io)
# end

max_stars_repo_name: biobakery/GaPLAC

max_issues_repo_name: biobakery/GaPLAC

max_forks_repo_name: biobakery/GaPLAC


text:
using PyPlot
using LinearAlgebra
using Statistics
using Random
using NLopt
using UKrig

##
# data grid and interpolation grid
n₀ = 100 # grid evals
n₁ = 500  # obs points
x₁ = rand(n₁)
y₁ = rand(n₁)

sidex₀ = range(-0.1, 1.1, length=n₀)
sidey₀ = range(-0.1, 1.1, length=n₀)
x₀ = sidex₀ .+ 0 .* sidey₀'
y₀ = 0 .* sidex₀ .+ sidey₀'

##
# set the simulation truth parameters
σεₒ = 0.1
σzₒ = 1.0
ρₒ  = 0.2
νₒ  = 1.25
mₒ  = 2
βₒ  = 2 * (rand(mₒ) .- 1/2)

##
# simulate data
dmat = UKrig.distmat((x₁,y₁), (x₁,y₁))
Σ₁ = (σzₒ^2) .* Mnu.(dmat ./ ρₒ, νₒ)
L₁ = cholesky(Σ₁).L
Z₁ = L₁ * randn(n₁)
Y₁ = Z₁ # .+ fpx.(1:mₒ,x₁')' * βₒ
data = Y₁ + σεₒ .* randn(n₁)

## ===============================
#

musr = 2
ν    = 1.25

loglike_inv_nu1, loglike_chol_nu1, wgrad_loglike_inv_nu1 = loglike_Gnu(
	data, x₁, y₁;
	musr=musr,
	ν=ν, # ≈ differentiability
)

#=
loglike_inv_nu1(0.1, 2.0)
loglike_chol_nu1(0.1, 2.0)
wgrad_loglike_inv_nu1(0.1, 2.0)
ϵ = 1e-7
(loglike_inv_nu1(0.1+ϵ, 2.0) - loglike_inv_nu1(0.1-ϵ, 2.0))/(2ϵ)
(loglike_inv_nu1(0.1, 2.0+ϵ) - loglike_inv_nu1(0.1, 2.0-ϵ))/(2ϵ)
=#

#LN_algm = [:LN_BOBYQA, :LN_COBYLA, :LN_PRAXIS, :LN_NELDERMEAD, :LN_SBPLX]
#LD_algm = [:LD_MMA, :LD_SLSQP, :LD_LBFGS, :LD_TNEWTON]
opt1 = Opt(:LN_BOBYQA, 2) # second arg is the number of variables to optimize
opt2 = Opt(:LN_BOBYQA, 2) # second arg is the number of variables to optimize
opt3 = Opt(:LD_LBFGS, 2) # second arg is the number of variables to optimize
opt1.max_objective = (x,grad) -> loglike_inv_nu1(x[1],x[2])
opt2.max_objective = (x,grad) -> loglike_chol_nu1(x[1],x[2])
opt3.max_objective = function (x,grad)
	ll, grad1, grad2 = wgrad_loglike_inv_nu1(x[1],x[2])
	if length(grad)>0
		grad[1]=grad1
		grad[2]=grad2
	end
	return ll
end
opt1.lower_bounds = [0.0, 0.0]
opt2.lower_bounds = [0.0, 0.0]
opt3.lower_bounds = [0.0, 0.0]
opt1.maxtime = 10 # in seconds
opt2.maxtime = 10
opt3.maxtime = 10
# opt1.maxeval = 50
# opt2.maxeval = 50
# opt3.maxeval = 50
# opt1.ftol_abs = 0.1 * sqrt(length(data)) / sqrt(2)
# opt2.ftol_abs = 0.1 * sqrt(length(data)) / sqrt(2)
# opt3.ftol_abs = 0.1 * sqrt(length(data)) / sqrt(2)
xinit1 = [0.1, 2.0]
xinit2 = [0.1, 2.0]
xinit3 = [0.1, 2.0]
optf1, optx1, ret1 = optimize!(opt1, xinit1)
optf2, optx2, ret2 = optimize!(opt2, xinit2)
optf3, optx3, ret3 = optimize!(opt3, xinit3)
@show sqrt(length(data))/√2

σe = optx1[1]
σg = optx1[2]

rng_σe = σe |> x->range(x - .05 * x, x + .05 * x, length=50)
rng_σg = σg |> x->range(x - .15 * x, x + .15 * x, length=50)
llmat = loglike_chol_nu1.(rng_σe,rng_σg')
pcolor(rng_σe, rng_σg, llmat)

krig1 = krig_Gnu(
	data, x₁, y₁;
	musr=musr,
	ν=ν,
	σe=σe,
	σg=σg,
)

figure()
scatter(x₁, y₁ ,c=data)

figure()
pcolor(x₀, y₀, krig1.(x₀, y₀))


text:
"""
This function add a service to the component without checking if the component and the service are attached to the same system
"""
function add_service_internal!(device::Device, service::Service)
    services = get_services(device)
    for _service in services
        if IS.get_uuid(service) == IS.get_uuid(_service)
            throw(
                ArgumentError(
                    "service $(get_name(service)) is already attached to $(get_name(device))",
                ),
            )
        end
    end

    push!(services, service)
    @debug "Add $service to $(get_name(device))" _group = IS.LOG_GROUP_SYSTEM
end

function add_service_internal!(device::Device, service::AGC)
    if !isa(device, RegulationDevice)
        throw(
            IS.ConflictingInputsError(
                "AGC service can only accept contributing devices of type RegulationDevice",
            ),
        )
    end

    device_bus_area = get_area(get_bus(device))
    service_area = get_area(service)
    if isnothing(device_bus_area) ||
       !(IS.get_uuid(device_bus_area) == IS.get_uuid(service_area))
        throw(
            IS.ConflictingInputsError(
                "Device $(get_name(device)) is not located in the regulation control area",
            ),
        )
    end

    services = get_services(device)
    for _service in services
        if IS.get_uuid(service) == IS.get_uuid(_service)
            throw(
                ArgumentError(
                    "service $(get_name(service)) is already attached to $(get_name(device))",
                ),
            )
        end
    end

    push!(services, service)
    @debug "Add $service to $(get_name(device))" _group = IS.LOG_GROUP_SYSTEM
end

"""
Remove a service from a device.

Throws ArgumentError if the service is not attached to the device.
"""
function remove_service!(device::Device, service::Service)
    if !_remove_service!(device, service)
        throw(
            ArgumentError(
                "service $(get_name(service)) was not attached to $(get_name(device))",
            ),
        )
    end
end

"""
Return true if the service is attached to the device.
"""
function has_service(device::Device, service::Service)
    for _service in get_services(device)
        if IS.get_uuid(_service) == IS.get_uuid(service)
            return true
        end
    end

    return false
end

"""
Return true if a service with type T is attached to the device.
"""
function has_service(device::Device, ::Type{T}) where {T <: Service}
    for _service in get_services(device)
        if isa(_service, T)
            return true
        end
    end

    return false
end

"""
Remove service from device if it is attached.
"""
function _remove_service!(device::Device, service::Service)
    removed = false
    services = get_services(device)

    # The expectation is that there won't be many services in each device, and so
    # a faster lookup method is not needed.
    for (i, _service) in enumerate(services)
        if IS.get_uuid(_service) == IS.get_uuid(service)
            deleteat!(services, i)
            removed = true
            @debug "Removed service $(get_name(service)) from $(get_name(device))" _group =
                IS.LOG_GROUP_SYSTEM
            break
        end
    end

    return removed
end

"""
Remove all services attached to the device.
"""
function clear_services!(device::Device)
    services = get_services(device)
    empty!(services)
end


text:
"""
    abstract type Abstract_BC_QuadSolver_Result end

Quadratic solver result abstraction

"""
abstract type Abstract_BC_QuadSolver_Result end

"""
    converged(::Abstract_BC_QuadSolver_Result)

Return `true` if the solver converged
"""
converged(::Abstract_BC_QuadSolver_Result) = @assert(false,"To implement")

"""
    iteration_count(::Abstract_BC_QuadSolver_Result)

Return the number of consumed iteration
"""
iteration_count(::Abstract_BC_QuadSolver_Result) = @assert(false,"To implement")

"""
    objective_value(::Abstract_BC_QuadSolver_Result)

Returns objective value at the point [`solution`](@ref).
"""
objective_value(r::Abstract_BC_QuadSolver_Result) = @assert(false,"To implement")

"""
    multiplier_τ(::Abstract_BC_QuadSolver_Result)

Returns the multipliers stored in a compact form (see τ definition, TODO)
"""
multiplier_τ(r::Abstract_BC_QuadSolver_Result) = @assert(false,"To implement")

"""
    solution(::Abstract_BC_QuadSolver_Result)

Returns the founded solution
"""
solution(r::Abstract_BC_QuadSolver_Result) = @assert(false,"To implement")


text:
#
# Set operations
# These are custom implementations of the set operations in
# abstractset.jl. However, they do not convert input collections
# to Sets. Therefore, they have higher complexity (O(N²) instead
# of O(N) for N elements) but they are faster for small sets.
# In addition, they assume that elements of the input collections
# are already unique, i.e. that the inputs are "set-like".
# Therefore, you should call a function like `unique` before inputting
# if you're not sure the collections themselves will be unique.
#

# A version of Base.setdiff that scales quadratically in the number of elements
# and assumes the elements of each input set are already unique.
_setdiff(s) = Base.copymutable(s)
_setdiff(s, itrs...) = _setdiff!(Base.copymutable(s), itrs...)
function _setdiff!(s, itrs...)
  for x in itrs
    _setdiff!(s, x)
  end
  return s
end
function _setdiff!(s, itr)
  isempty(s) && return s
  for x in itr
    n = findfirst(==(x), s)
    !isnothing(n) && deleteat!(s, n)
  end
  return s
end

# A version of Base.intersect that scales quadratically in the number of elements
# and assumes the elements of each input set are already unique.
_intersect(s) = Base.copymutable(s)
_intersect(s, itr, itrs...) = _intersect!(_intersect(s, itr), itrs...)
# XXX: Base has `s` and `itr` swapped in the definition, which one is correct?
# Is this special case needed, or is `filter!` sufficient?
_intersect(s, itr) = Base.mapfilter(in(itr), push!, s, Base.emptymutable(s, eltype(s)))
function _intersect!(s, itrs...)
  for x in itrs
    _intersect!(s, x)
  end
  return s
end
_intersect!(s, s2) = filter!(in(s2), s)

# A version of Base.symdiff that scales quadratically in the number of elements
# and assumes the elements of each input set are already unique.
_symdiff(s) = Base.copymutable(s)
function _symdiff(s, itrs...)
  return _symdiff!(Base.emptymutable(s, Base.promote_eltype(s, itrs...)), s, itrs...)
end
function _symdiff!(s, itrs...)
  for x in itrs
    _symdiff!(s, x)
  end
  return s
end
function _symdiff!(s, itr)
  if isempty(s)
    append!(s, itr)
    return s
  end
  for x in itr
    n = findfirst(==(x), s)
    !isnothing(n) ? deleteat!(s, n) : push!(s, x)
  end
  return s
end

# A version of Base.union that scales quadratically in the number of elements
# and assumes the elements of each input set are already unique.
_union(s) = Base.copymutable(s)
function _union(s, sets...)
  return _union!(Base.emptymutable(s, Base.promote_eltype(s, sets...)), s, sets...)
end
function _union!(s, sets...)
  for x in sets
    _union!(s, x)
  end
  return s
end
function _union!(s, itr)
  Base.haslength(itr) && sizehint!(s, length(s) + Int(length(itr))::Int)
  for x in itr
    x ∉ s && push!(s, x)
  end
  return s
end


text:
using Turing
using Distributions

const p = 0.5

#@sf geom() = normal(0, 1) > 0. ? 1+geom() : 1
@sf geom(p) = bool(bernoulli(p)) ? 1+geom(p) : 1

@sf test_geom() = begin
    a = geom(0.5)
    normal(a, 3; condition=15)
    a
end

as = float(sample(test_geom, 10^6))

ms = mean(as)
vs = sqrt(var(as))

#Analytic
n = 30
xs = 1:n
prior = p.^xs
likelihood = pdf(Normal(15, 3), xs)
posterior = prior .* likelihood
norm_posterior = posterior/sum(posterior)
ma = sum(norm_posterior.*xs)
va = sqrt(sum(norm_posterior.*xs.^2) - ma^2)

println((ms, ma))
println((vs, va))


text:
module HTsi

using Interpolations
using RemoteFiles
using OptionalData
using StaticArrays
using Statistics
using Dates
using Printf

const path_si_root = joinpath(homedir(), "RD", "SPI")

include("space_indices.jl")

greet() = print("Hello World! HTsi")

end # module


text:
using Arrhenius
using Sundials
using LinearAlgebra
using OrdinaryDiffEq
using ForwardDiff
using ForwardDiff: jacobian, jacobian!
using DiffEqSensitivity
using YAML
using BSON: @save, @load
using Plots, Printf, Random
using Statistics
using ProgressBars
using Base.Threads
using BandedMatrices
using SparseArrays
using Dierckx # for interpolation

ENV["GKSwstype"] = "100"

# load input.yaml
runtime = YAML.load_file("./input.yaml")
expr_name = runtime["expr_name"]
is_restart = runtime["is_restart"]

# load $expr_name/config.yaml
conf = YAML.load_file("$expr_name/config.yaml")
mech = conf["mech"]
fuel = conf["fuel"]
fuel2air = Float64(conf["fuel2air"])
oxygen = conf["oxygen"]
inert = conf["inert"]

dT = Float64(conf["dTign"])
dTabort = Float64(conf["dTabort"])
method = conf["method"]

n_plot = Int64(conf["n_plot"])

if is_restart
    println("Continue to run $expr_name ...\n")
else
    println("Runing $expr_name ...\n")
end

exp_path = string(expr_name)
fig_path = string(exp_path, "/figs")
ckpt_path = string(exp_path, "/checkpoint")

if !is_restart
    if ispath(fig_path)
        rm(fig_path, recursive=true)
    end
    if ispath(ckpt_path)
        rm(ckpt_path, recursive=true)
    end
end

if ispath(fig_path) == false
    mkdir(fig_path)
    mkdir(string(fig_path, "/conditions"))
end

if ispath(ckpt_path) == false
    mkdir(ckpt_path)
end

max_stars_repo_name: DENG-MIT/ArrheniusActiveSubspace

max_issues_repo_name: DENG-MIT/ArrheniusActiveSubspace

max_forks_repo_name: DENG-MIT/ArrheniusActiveSubspace


text:
# dummy allocator that passes through any requests, calling into the GC if that fails.

pool_init() = return

function pool_alloc(dev, sz)
    block = nothing
    for phase in 1:3
        if phase == 2
            @pool_timeit "$phase.0 gc (incremental)" GC.gc(false)
        elseif phase == 3
            @pool_timeit "$phase.0 gc (full)" GC.gc(true)
        end

        @pool_timeit "$phase.1 alloc" begin
            block = actual_alloc(dev, sz)
        end
        block === nothing || break
    end

    return block
end

function pool_free(dev, block)
    actual_free(dev, block)
    return
end

pool_reclaim(dev, target_bytes::Int=typemax(Int)) = return 0

cached_memory(dev=device()) = 0


text:
using Tar

dir = "./mort.soa.org_raw"

out_dir = mkpath("./mort.soa.org")

# create tarballs from the
for (root, dirs, files) in walkdir(dir)
    name = splitpath(root)[end]
    println("Directories in $root,$dir,$name")
    Tar.create(root,"$(joinpath(out_dir,name)).tar.gz")
end

max_stars_repo_name: JuliaActuary/Artifacts

max_issues_repo_name: JuliaActuary/Artifacts

max_forks_repo_name: JuliaActuary/Artifacts


text:
# Two separate comparisons

# How does a given visibility vary with inclination? (fixed nu)
# How does a given visibility vary with nu? (fixed incl)

using DiskJockey.constants
using DiskJockey.model
using DiskJockey.visibilities
using DiskJockey.image
using HDF5

import PyPlot.plt
using LaTeXStrings

cmap = plt[:get_cmap]("plasma")

import YAML
config = YAML.load(open("config.yaml"))

grd = config["grid"]
grid = Grid(grd)

# Load the parameters file from the configuration setting, and make this an object.
pars = convert_dict(config["parameters"], config["model"])
vel = pars.vel # [km/s]
npix = config["npix"] # number of pixels
species = config["species"]
transition = config["transition"]
lam0 = lam0s[species*transition]
model = config["model"]

nincl = 10
nchan = 10

# Instead of reading in wavelengths from data file, use those we've specified
# Spanning the central 3 chanels of DM Tau, but w/ high resolution.
freqs = linspace(2.3052281389490067e11, 2.3052365978614038e11, nchan)
# Convert from Hz to wavelengths in μm
lams = cc ./freqs * 1e4 # [μm]

# The list of inclinations we want to create
incls = linspace(140., 158., nincl)

# Make a single plot
# Choose a visibility at random
# And plot it's value (real and imag) as we go from inc = 140 to 160

# chan_index = 6 # Choose just the first channel, for now (fixed nu)

dvarr = DataVis(config["data_file"], true)
nvis = length(dvarr[1].VV)

fig, ax = plt[:subplots](nrows=2, figsize=(6,6), sharex=true)

# Load one model into memory
mv = DataVis("model_140.hdf5", true)

mvs = Array(Complex128, (nchan, nvis))

for i=1:nchan
    #Load the visibilities corresponding to each channel
    mvs[i, :] = mv[i].VV
end

for vis_index in rand(1:nvis, 50)
    vis = mvs[:, vis_index]

    ax[1][:plot](real(vis)/mean(real(vis)), lw=0.2)
    ax[2][:plot](imag(vis)/mean(imag(vis)), lw=0.2)
end

ax[1][:set_ylabel]("real")
ax[2][:set_ylabel]("imag")
ax[2][:set_xlabel]("channel")

fig[:savefig]("visibilities_chan.png")

max_stars_repo_name: iancze/DiskEmulator

max_issues_repo_name: iancze/DiskEmulator

max_forks_repo_name: iancze/DiskEmulator


text:
function backproject_transform(pixels, K_inv, R, t)
    points = Vector{SVector{3, Float64}}(undef, length(pixels))
    for (i, px) in enumerate(pixels)
        p = K_inv * SVector{3, Float64}(px..., 1)
        points[i] = R * p + t
    end
    points
end

@testset "P3P" begin
    n_points = 100

    fcmin, fcmax = 1e-3, 100
    K = SMatrix{3, 3, Float64}(
        rand() * (fcmax - fcmin) + fcmin, 0, 0,
        0, rand() * (fcmax - fcmin) + fcmin, 0,
        rand() * (fcmax - fcmin) + fcmin, rand() * (fcmax - fcmin) + fcmin, 1,
    )
    K_inv = K |> inv

    θ = π / 4
    R = RotXYZ(rand() * θ, rand() * θ, rand() * θ)
    t = SVector{3, Float64}(rand(), rand() * 2, rand() * 3)
    P_target = RecoverPose.get_transformation(R', -R' * t)

    max_res = 2 * min(K[1, 3], K[2, 3])
    pmin = SVector{2}(1, 1)
    pmax = SVector{2}(max_res, max_res)
    δ = pmax - pmin

    pixels = [floor.(rand(SVector{2, Float64}) .* δ .+ pmin) for i in 1:n_points]
    pdn_pixels = pre_divide_normalize(pixels, K)
    points = backproject_transform(pixels, K_inv, R, t)

    models = p3p(points[1:3], pdn_pixels[1:3], K)
    n_inliers, (KP, inliers, error) = p3p_select_model(models, points, pixels; threshold=0.1)
    @test n_inliers == n_points
    @test sum(inliers) == n_inliers
    @test all(isapprox.(K_inv * KP, P_target; atol=1e-2))

    models = p3p(points[1:3], pixels[1:3], K)
    n_inliers, (KP, inliers, error) = p3p_select_model(models, points, pixels)
    @test n_inliers == n_points
    @test sum(inliers) == n_inliers
    @test all(isapprox.(K_inv * KP, P_target; atol=1e-2))

    # Test P3P RANSAC.

    n_inliers, (KP, inliers, error) = p3p_ransac(points, pixels, pdn_pixels, K)
    @test n_inliers == n_points
    @test sum(inliers) == n_inliers
    @test all(isapprox.(K_inv * KP, P_target; atol=1e-2))

    n_inliers, (KP, inliers, error) = p3p_ransac(points, pixels, K)
    @test n_inliers == n_points
    @test sum(inliers) == n_inliers
    @test all(isapprox.(K_inv * KP, P_target; atol=1e-2))
end

@testset "P3P RANSAC with noise and outliers." begin
    n_points = 100
    noise_scale = 1e-3
    rot_atol = 0.5
    t_atol = 2

    fcmin, fcmax = 1e-3, 100
    K = SMatrix{3, 3, Float64}(
        rand() * (fcmax - fcmin) + fcmin, 0, 0,
        0, rand() * (fcmax - fcmin) + fcmin, 0,
        rand() * (fcmax - fcmin) + fcmin, rand() * (fcmax - fcmin) + fcmin, 1,
    )
    K_inv = K |> inv

    θ = π / 8
    R = RotXYZ(rand() * θ, rand() * θ, rand() * θ)
    t = SVector{3, Float64}(rand(), rand() * 2, rand() * 3)
    P_target = RecoverPose.get_transformation(R', -R' * t)

    max_res = 2 * min(K[1, 3], K[2, 3])
    pmin = SVector{2}(1, 1)
    pmax = SVector{2}(max_res, max_res)
    δ = pmax - pmin

    pixels = [floor.(rand(SVector{2, Float64}) .* δ .+ pmin) for i in 1:n_points]
    pdn_pixels = pre_divide_normalize(pixels, K)
    points = backproject_transform(pixels, K_inv, R, t)
    points = [p .+ rand(SVector{3, Float64}) .* noise_scale for p in points]

    # Add outliers.
    outlier_scale = 5
    points[1] = points[1] .+ rand() * outlier_scale
    points[5] = points[5] .+ rand() * outlier_scale
    points[end] = points[end] .+ rand() * outlier_scale
    points[end - 5] = points[end - 5] .+ rand() * outlier_scale
    points[n_points ÷ 2] = points[n_points ÷ 2] .+ rand() * outlier_scale

    n_inliers, (KP, inliers, error) = p3p_ransac(points, pixels, pdn_pixels, K)
    P = K_inv * KP
    @test n_inliers ≥ n_points - 10
    @test !inliers[1] && !inliers[5] && !inliers[end] && !inliers[end - 5] &&
        !inliers[n_points ÷ 2]

    @test all(isapprox.(P[1:3, 1:3], P_target[1:3, 1:3]; atol=rot_atol))
    @test all(isapprox.(P[1:3, 4], P_target[1:3, 4]; atol=t_atol))

    n_inliers, (KP, inliers, error) = p3p_ransac(points, pixels, K)
    P = K_inv * KP
    @test n_inliers ≥ n_points - 10
    @test !inliers[1] && !inliers[5] && !inliers[end] && !inliers[end - 5] &&
        !inliers[n_points ÷ 2]

    @test all(isapprox.(P[1:3, 1:3], P_target[1:3, 1:3]; atol=rot_atol))
    @test all(isapprox.(P[1:3, 4], P_target[1:3, 4]; atol=t_atol))
end


text:
# ---
# title: 1093. Statistics from a Large Sample
# id: problem1093
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math, Two Pointers
# link: <https://leetcode.com/problems/statistics-from-a-large-sample/description/>
# hidden: true
# ---
#
# We sampled integers between `0` and `255`, and stored the results in an array
# `count`:  `count[k]` is the number of integers we sampled equal to `k`.
#
# Return the minimum, maximum, mean, median, and mode of the sample
# respectively, as an array of **floating point numbers**.  The mode is
# guaranteed to be unique.
#
# _(Recall that the median of a sample is:_
#
#   * _The middle element, if the elements of the sample were sorted and the number of elements is odd;_
#   * _The average of the middle two elements, if the elements of the sample were sorted and the number of elements is even.)_
#
#
#
# **Example 1:**
#
#
#
#     Input: count = [0,1,3,4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
#     Output: [1.00000,3.00000,2.37500,2.50000,3.00000]
#
#
# **Example 2:**
#
#
#
#     Input: count = [0,4,3,2,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
#     Output: [1.00000,4.00000,2.18182,2.00000,1.00000]
#
#
#
#
# **Constraints:**
#
#   1. `count.length == 256`
#   2. `1 <= sum(count) <= 10^9`
#   3. The mode of the sample that count represents is unique.
#   4. Answers within `10^-5` of the true value will be accepted as correct.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module JuliaParserHooks
using JuliaParser
using JuliaParser.Lexer: ¬
using JuliaParser: Lexer

# mute override warnings
if ccall(:jl_generating_output, Cint, ()) == 0
    ORIG_STDERR = STDERR
    redirect_stderr()
end

immutable REPLDiagnostic
    fname::AbstractString
    text::AbstractString
    d::Any
end

extract_diag(d) = isa(d, Main.JuliaParser.Diagnostics.Incomplete) ? d.d : d

function Base.showerror(io::IO, d::REPLDiagnostic, bt)
    print_with_color(:white,io,"")
    JuliaParser.Diagnostics.display_diagnostic(io, d.text, extract_diag(d.d); filename = d.fname)
end

function Base.showerror(io::IO, d::REPLDiagnostic)
    print_with_color(:white,io,"")
    JuliaParser.Diagnostics.display_diagnostic(io, d.text, extract_diag(d.d); filename = d.fname)
end

function Base.REPL.display_error(io::IO, d::REPLDiagnostic, bt)
    print_with_color(:white,io,"")
    JuliaParser.Diagnostics.display_diagnostic(io, d.text, extract_diag(d.d); filename = d.fname)
end

function Core.include(fname::String)
    io = open(fname)
    file = readstring(io)
    ts = Lexer.TokenStream{Lexer.SourceLocToken}(file)
    ts.filename = fname
    local result = nothing
    while !Lexer.eof(ts)
        ast = try
            Parser.parse(ts)
        catch e
            !isa(e, Main.JuliaParser.Diagnostics.AbstractDiagnostic) && rethrow(e)
            rethrow(REPLDiagnostic(fname, file, e))
        end
        result = ccall(:jl_toplevel_eval, Any, (Any,), ¬ast)
    end
    result
end

function Base.include_string(code, filename = "string")
    ts = Lexer.TokenStream{Lexer.SourceLocToken}(code)
    ts.filename = filename
    local result = nothing
    while !Lexer.eof(ts)
        ast = try
            Parser.parse(ts)
        catch e
            !isa(e, Main.JuliaParser.Diagnostics.AbstractDiagnostic) && rethrow(e)
            rethrow(REPLDiagnostic(fname, file, e))
        end
        result = ccall(:jl_toplevel_eval, Any, (Any,), ¬ast)
    end
    result
end

# Needs to be replaced due to # 265
function Base.include_from_node1(_path::AbstractString)
    path, prev = Base._include_dependency(_path)
    tls = task_local_storage()
    tls[:SOURCE_PATH] = path
    local result
    try
        if myid()==1
            # sleep a bit to process file requests from other nodes
            nprocs()>1 && sleep(0.005)
            result = Core.include(path)
            nprocs()>1 && sleep(0.005)
        else
            result = include_string(remotecall_fetch(readstring, 1, path), path)
        end
    finally
        if prev === nothing
            delete!(tls, :SOURCE_PATH)
        else
            tls[:SOURCE_PATH] = prev
        end
    end
    result
end

function Base.parse(str::AbstractString, pos::Int; greedy::Bool=true, raise::Bool=true)
    io = IOBuffer(str)
    seek(io, pos-1)
    ts = Main.JuliaParser.Lexer.TokenStream{Main.JuliaParser.Lexer.SourceLocToken}(io)
    local result = nothing
    ast = try
        Parser.parse(ts; production = greedy ? Parser.parse_stmts : Parser.parse_atom)
    catch e
        !isa(e, Main.JuliaParser.Diagnostics.AbstractDiagnostic) && rethrow(e)
        e = REPLDiagnostic("REPL", str, e)
        raise && rethrow(e)
        return Expr(isa(e, Main.JuliaParser.Diagnostics.Incomplete) ? :incomplete :
            :error, e), position(io) + 1
    end
    ¬ast, position(io) + 1
end

function Base.parse_input_line(code::String; filename::String="none")
    ts = Lexer.TokenStream{Main.JuliaParser.Lexer.SourceLocToken}(code)
    ts.filename = filename
    local result = nothing
    ast = try
        Parser.parse(ts)
    catch e
        !isa(e, Main.JuliaParser.Diagnostics.AbstractDiagnostic) && rethrow(e)
        return Expr(isa(e, Main.JuliaParser.Diagnostics.Incomplete) ? :incomplete :
            :error, REPLDiagnostic("REPL", code, e))
    end
    ¬ast
end

function Base.incomplete_tag(e::Expr)
    isa(e.args[1], REPLDiagnostic) || return :none
    isa(e.args[1].d, Main.JuliaParser.Diagnostics.Incomplete) || return :none
    e.args[1].d.tag
end

if ccall(:jl_generating_output, Cint, ()) == 0
    REDIRECTED_STDERR = STDERR
    redirect_stderr(ORIG_STDERR)
end

end


text:
###
### mappedarray
###
function MappedArrays.mappedarray(f, data::NamedDimsArray{L}) where {L}
    return NamedDimsArray{L}(mappedarray(f, parent(data)))
end

function MappedArrays.mappedarray(::Type{T}, data::NamedDimsArray{L}) where {T,L}
    return NamedDimsArray{L}(mappedarray(T, parent(data)))
end

function MappedArrays.mappedarray(f, data::NamedDimsArray...)
    dn = _unify_names(map(dimnames, data))
    return NamedDimsArray{dn}(mappedarray(f, map(parent, data)...))
end

function MappedArrays.mappedarray(::Type{T}, data::NamedDimsArray...) where T
    dn = _unify_names(map(dimnames, data))
    return NamedDimsArray{dn}(mappedarray(T, map(parent, data)...))
end

function MappedArrays.mappedarray(f, finv::Function, data::NamedDimsArray{L}) where {L}
    return NamedDimsArray{L}(mappedarray(f, parent(data)))
end

function MappedArrays.mappedarray(f, finv::Function, data::NamedDimsArray...)
    dn = _unify_names(map(dimnames, data))
    return NamedDimsArray{dn}(mappedarray(f, finv, map(parent, data)...))
end

function MappedArrays.mappedarray(::Type{T}, finv::Function, data::NamedDimsArray...) where T
    dn = _unify_names(map(dimnames, data))
    return NamedDimsArray{dn}(mappedarray(T, finv, map(parent, data)...))
end

function MappedArrays.mappedarray(f, ::Type{Finv}, data::NamedDimsArray...) where Finv
    dn = _unify_names(map(dimnames, data))
    return NamedDimsArray{dn}(mappedarray(f, Finv, map(parent, data)...))
end

function MappedArrays.mappedarray(::Type{T}, ::Type{Finv}, data::NamedDimsArray...) where {T,Finv}
    dn = _unify_names(map(dimnames, data))
    return NamedDimsArray{dn}(mappedarray(T, Finv, map(parent, data)...))
end

_unify_names(x::Tuple{Any}) = first(x)
function _unify_names(x::Tuple{Any,Vararg{Any}})
    return NamedDims.unify_names_longest(first(x), _unify_names(tail(x)))
end

# TODO this should be defined in NamedDims
@inline function NamedDims.dimnames(::Type{Base.PermutedDimsArray{T,N,permin,permout,A}}) where {T,N,permin,permout,A}
    dn = dimnames(A)
    return map(i -> getfield(dn, i), permin)
end
@inline function NamedDims.dimnames(::Type{Base.ReinterpretArray{T,N,S,A}}) where {T,N,S,A}
    return dimnames(A)
end
NamedDims.dimnames(::Type{ReadonlyMappedArray{T,N,A,F}}) where {T,N,A,F} = dimnames(A)
NamedDims.dimnames(::Type{MappedArray{T,N,A,F,Finv}}) where {T,N,A,F,Finv} = dimnames(A)
function NamedDims.dimnames(::Type{ReadonlyMultiMappedArray{T,N,AAs,F}}) where {T,N,AAs,F}
    return _multi_array_dimnames(AAs, ntuple(_ -> :_, Val(N)))
end
function NamedDims.dimnames(::Type{MultiMappedArray{T,N,AAs,F,Finv}}) where {T,N,AAs,F,Finv}
    return _multi_array_dimnames(AAs, ntuple(_ -> :_, Val(N)))
end

@inline function _multi_array_dimnames(::Type{T}, dnames::Tuple{Vararg{Symbol}}) where {T}
    for T_i in T.parameters
        dnames = NamedDims.unify_names_longest(dnames, dimnames(T_i))
    end
    return dnames
end

#=_multi_array_dimnames(::Tuple{T}) where {T} = dimnames(T)
_multi_array_dimnames(::Tuple{T1,T2}) where {T1,T2} = NamedDims.unify_names_longest(dimnames(T1), dimnames(T2))
@inline function _multi_array_dimnames(x::Tuple{T1,T2,Vararg{Any}}) where {T1,T2}
    return NamedDims.unify_names_longest(dimnames(T1), _multi_array_dimnames(tail(x)))
end
=#

ArrayInterface.parent_type(::Type{<:NamedDimsArray{L,T,N,A}}) where {L,T,N,A} = A

has_dimnames(x) = has_dimnames(typeof(x))
has_dimnames(::Type{T}) where {T<:NamedDimsArray} = true
function has_dimnames(::Type{T}) where {T}
    if parent_type(T) <: T
        return false
    else
        return has_dimnames(parent_type(T))
    end
end
function Metadata.attach_metadata(data::NamedDimsArray, m::Metadata.METADATA_TYPES)
    return NamedDimsArray{dimnames(data)}(attach_metadata(parent(data), m))
end

"""
    NamedAxisArray(parent::AbstractArray; kwargs...) = NamedAxisArray(parent, kwargs)
    NamedAxisArray(parent::AbstractArray, axes::NamedTuple{L,AbstractAxes})

Type alias for `NamedDimsArray` whose parent array is an `AxisArray`. If key word
arguments are provided then each key word becomes the name of a dimension and its
assigned value is sent to the corresponding axis when constructing the underlying
`AxisArray`.

## Examples
```jldoctest
julia> using AxisIndices

julia> A = NamedAxisArray{(:x, :y, :z)}(reshape(1:24, 2, 3, 4), ["a", "b"], ["one", "two", "three"], 2:5)
2×3×4 NamedDimsArray(AxisArray(reshape(::UnitRange{Int64}, 2, 3, 4)
  • axes:
     x = ["a", "b"]
     y = ["one", "two", "three"]
     z = 2:5
))
[:, :, 2] =
       "one"   "two"   "three"
  "a"  1       3       5
  "b"  2       4       6

[:, :, 3] =
       "one"   "two"    "three"
  "a"  7        9       11
  "b"  8       10       12

[:, :, 4] =
       "one"    "two"    "three"
  "a"  13       15       17
  "b"  14       16       18

[:, :, 5] =
       "one"    "two"    "three"
  "a"  19       21       23
  "b"  20       22       24

julia> B = A["a", :, :]
3×4 NamedDimsArray(AxisArray(::Array{Int64,2}
  • axes:
     y = ["one", "two", "three"]
     z = 2:5
))
           2  3   4   5
  "one"    1   7  13  19
  "two"    3   9  15  21
  "three"  5  11  17  23

julia> C = B["one",:]
4-element NamedDimsArray(AxisArray(::Array{Int64,1}
  • axes:
     z = 2:5
))
     1
  2   1
  3   7
  4  13
  5  19

```
"""
const NamedAxisArray{L,T,N,P,AI} = NamedDimsArray{L,T,N,AxisArray{T,N,P,AI}}

NamedAxisArray{L}(x::AxisArray) where {L} = NamedDimsArray{L}(x)

function NamedAxisArray{L}(x::AbstractArray, axs::Tuple) where {L}
    return NamedAxisArray{L}(AxisArray(x, axs))
end

function NamedAxisArray{L}(x::AbstractArray, args::AbstractVector...) where {L}
    return NamedAxisArray{L}(x, args)
end

function NamedAxisArray(x::AbstractArray, axs::NamedTuple{L}) where {L}
    return NamedAxisArray{L}(x, values(axs))
end

function NamedAxisArray{L,T,N}(init::ArrayInitializer, axs::Tuple) where {L,T,N}
    return NamedAxisArray{L}(AxisArray{T,N}(init, axs))
end

function NamedAxisArray{L,T,N}(init::ArrayInitializer, args::AbstractVector...) where {L,T,N}
    return NamedAxisArray{L,T,N}(init, args)
end

function NamedAxisArray{L,T}(init::ArrayInitializer, axs::Tuple) where {L,T}
    return NamedAxisArray{L}(AxisArray{T}(init, axs))
end

function NamedAxisArray{L,T}(init::ArrayInitializer, args::AbstractVector...) where {L,T}
    return NamedAxisArray{L,T}(init, args)
end

NamedAxisArray(x::AbstractArray; kwargs...) = NamedAxisArray(x, kwargs.data)

for f in (:getindex, :view, :dotview)
    @eval begin
        @propagate_inbounds function Base.$f(A::NamedAxisArray; named_inds...)
            inds = NamedDims.order_named_inds(A; named_inds...)
            return Base.$f(A, inds...)
        end

        @propagate_inbounds function Base.$f(a::NamedAxisArray, raw_inds...)
            inds = to_indices(parent(a), raw_inds)  # checkbounds happens within to_indices
            data = @inbounds(Base.$f(parent(a), inds...))
            data isa AbstractArray || return data # Case of scalar output
            L = NamedDims.remaining_dimnames_from_indexing(dimnames(a), inds)
            if L === ()
                # Cases that merge dimensions down to vector like `mat[mat .> 0]`,
                # and also zero-dimensional `view(mat, 1,1)`
                return data
            else
                return NamedDimsArray{L}(data)
            end
        end
    end
end

Base.getindex(A::NamedAxisArray, ::Ellipsis) = A

const ReinterpretNamedAxisArray{T,N,S,L,A<:NamedAxisArray{L,S,N}} = ReinterpretArray{T,N,S,A}

function Base.axes(A::ReinterpretNamedAxisArray{T,N,S}) where {T,N,S}
    paxs = axes(parent(A))
    axis_1 = first(paxs)
    len = div(length(axis_1) * sizeof(S), sizeof(T))
    return tuple(StaticRanges.resize_last(axis_1, len), tail(paxs)...)
end

function Base.BroadcastStyle(a::AxisArrayStyle{A}, b::NamedDims.NamedDimsStyle{B}) where {A,B}
    return NamedDims.NamedDimsStyle(a, B())
end
function Base.BroadcastStyle(a::NamedDims.NamedDimsStyle{M}, b::AxisArrayStyle{B}) where {B,M}
    return NamedDims.NamedDimsStyle(M(), b)
end

function Base.show(io::IO, ::MIME"text/plain", X::NamedDimsArray{<:Any,<:Any,<:Any,<:Union{<:AxisArray,Metadata.MetaArray}})
    if isempty(X) && (get(io, :compact, false) || X isa Vector)
        return show(io, X)
    end
    # 0) show summary before setting :compact
    summary(io, X)
    isempty(X) && return
    Base.show_circular(io, X) && return

    # 1) compute new IOContext
    if !haskey(io, :compact) && length(axes(X, 2)) > 1
        io = IOContext(io, :compact => true)
    end
    if get(io, :limit, false) && eltype(X) === Method
        # override usual show method for Vector{Method}: don't abbreviate long lists
        io = IOContext(io, :limit => false)
    end

    if get(io, :limit, false) && displaysize(io)[1]-4 <= 0
        return print(io, " …")
    else
        println(io)
    end

    # 2) update typeinfo
    #
    # it must come after printing the summary, which can exploit :typeinfo itself
    # (e.g. views)
    # we assume this function is always called from top-level, i.e. that it's not nested
    # within another "show" method; hence we always print the summary, without
    # checking for current :typeinfo (this could be changed in the future)
    io = IOContext(io, :typeinfo => eltype(X))

    # 2) show actual content
    recur_io = IOContext(io, :SHOWN_SET => X)
    Base.print_array(recur_io, parent(X))
end
function Base.summary(io::IO, x::NamedDimsArray{<:Any,<:Any,<:Any,<:Union{<:AxisArray,Metadata.MetaArray}})
    return Base.showarg(io, x, true)
end

function Base.showarg(io::IO, x::NamedDimsArray{<:Any,<:Any,<:Any,<:AxisArray}, toplevel)
    if toplevel
        print(io, Base.dims2string(length.(axes(x))), " ")
    end

    print(io, "NamedDimsArray(")
    print(io, "AxisArray(")
    Base.showarg(io, parent(parent(x)), false)
    print(io, "\n")
    print_axes_summary(io, NamedTuple{dimnames(x)}(axes(x)))
    print(io, "\n))")
end

function Base.showarg(io::IO, x::NamedDimsArray{<:Any,<:Any,<:Any,<:Metadata.MetaArray}, toplevel)
    if toplevel
        print(io, Base.dims2string(length.(axes(x))), " ")
    end
    print(io, "NamedDimsArray(")
    print(io, "attach_metadata(")

    if parent_type(parent_type(x)) <: AxisArray
        print(io, "AxisArray(")
        Base.showarg(io, parent(parent(parent(x))), false)
        print(io, "\n")
        print_axes_summary(io, NamedTuple{dimnames(x)}(axes(x)))
        print(io, "\n)")
    else
        Base.showarg(io, parent(x), false)
    end
    print(io, ", ", Metadata.showarg_metadata(x))
    println(io)
    Metadata.metadata_summary(io, x)
    print(io, "\n)") # closing bracket on attach_metadata
    print(io, ")")   # closing bracket on NamedDimsArray
end