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, AbstractDict, IndexLinear, IndexCartesian, IndexStyle, AbstractVector, Vector,
    tail, tuple_type_head, tuple_type_tail, tuple_type_cons, SizeUnknown, HasLength, HasShape,
    IsInfinite, EltypeUnknown, HasEltype, OneTo, @propagate_inbounds, Generator, AbstractRange,
    linearindices, (:), |, +, -, !==, !, <=, <

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

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

_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 start/next/done 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 [`start`](@ref), [`next`](@ref), and [`done`](@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
first(r::Reverse) = last(r.itr) # and the last shall be first

# reverse-order array iterators: assumes more-specialized Reverse for eachindex
@inline start(A::Reverse{<:AbstractArray}) = (itr = reverse(eachindex(A.itr)); (itr, start(itr)))
@propagate_inbounds next(A::Reverse{<:AbstractArray}, i) = ((idx, s) = next(i[1], i[2]); (A.itr[idx], (i[1], s)))
@propagate_inbounds done(A::Reverse{<:AbstractArray}, i) = done(i[1], i[2])

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

start(r::Reverse{<:Tuple}) = length(r.itr)
done(r::Reverse{<:Tuple}, i::Int) = i < 1
next(r::Reverse{<:Tuple}, i::Int) = (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 `enumerate(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)
@inline start(e::Enumerate) = (1, start(e.itr))
@propagate_inbounds function next(e::Enumerate, state)
    n = next(e.itr,state[2])
    (state[1],n[1]), (state[1]+1,n[2])
end
@inline done(e::Enumerate, state) = done(e.itr, state[2])

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 start(r::Reverse{<:Enumerate})
    ri = reverse(r.itr.itr)
    return (length(ri), ri, start(ri))
end
@inline function next(r::Reverse{<:Enumerate}, state)
    n = next(state[2],state[3])
    (state[1],n[1]), (state[1]-1,state[2],n[2])
end
@inline done(r::Reverse{<:Enumerate}, state) = state[1] < 1

"""
    Iterators.Pairs(values, keys) <: AbstractDict{eltype(keys), eltype(values)}

Transforms an indexable container into an Dictionary-view of the same data.
Modifying the key-space of the underlying data may invalidate this object.
"""
struct Pairs{K, V, I, A} <: AbstractDict{K, V}
    data::A
    itr::I
    Pairs(data::A, itr::I) where {A, I} = new{eltype(I), eltype(A), I, A}(data, itr)
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()` ensures that `i` will be an integer;
specifying `IndexCartesian()` ensures that `i` will be a
`CartesianIndex`; 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(A::AbstractArray)  = pairs(IndexCartesian(), A)
pairs(A::AbstractVector) = pairs(IndexLinear(), A)
pairs(tuple::Tuple) = Pairs(tuple, keys(tuple))
pairs(nt::NamedTuple) = Pairs(nt, keys(nt))
# pairs(v::Pairs) = v # listed for reference, but already defined from being an AbstractDict

length(v::Pairs)  = length(v.itr)
axes(v::Pairs) = axes(v.itr)
size(v::Pairs)    = size(v.itr)
@inline start(v::Pairs) = start(v.itr)
@propagate_inbounds function next(v::Pairs{K, V}, state) where {K, V}
    indx, n = next(v.itr, state)
    item = v.data[indx]
    return (Pair{K, V}(indx, item), n)
end
@inline done(v::Pairs, state) = done(v.itr, state)

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

reverse(v::Pairs) = Pairs(v.data, reverse(v.itr))

haskey(v::Pairs, key) = (key in v.itr)
keys(v::Pairs) = v.itr
values(v::Pairs) = v.data
getindex(v::Pairs, key) = v.data[key]
setindex!(v::Pairs, value, key) = (v.data[key] = value; v)
get(v::Pairs, key, default) = get(v.data, key, default)
get(f::Base.Callable, collection::Pairs, key) = get(f, v.data, key)

# zip

abstract type AbstractZipIterator end

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()

struct Zip1{I} <: AbstractZipIterator
    a::I
end
zip(a) = Zip1(a)
length(z::Zip1) = length(z.a)
size(z::Zip1) = size(z.a)
axes(z::Zip1) = axes(z.a)
eltype(::Type{Zip1{I}}) where {I} = Tuple{eltype(I)}
@inline start(z::Zip1) = start(z.a)
@propagate_inbounds function next(z::Zip1, st)
    n = next(z.a,st)
    return ((n[1],), n[2])
end
@inline done(z::Zip1, st) = done(z.a,st)

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

struct Zip2{I1, I2} <: AbstractZipIterator
    a::I1
    b::I2
end
zip(a, b) = Zip2(a, b)
length(z::Zip2) = _min_length(z.a, z.b, IteratorSize(z.a), IteratorSize(z.b))
size(z::Zip2) = promote_shape(size(z.a), size(z.b))
axes(z::Zip2) = promote_shape(axes(z.a), axes(z.b))
eltype(::Type{Zip2{I1,I2}}) where {I1,I2} = Tuple{eltype(I1), eltype(I2)}
@inline start(z::Zip2) = (start(z.a), start(z.b))
@propagate_inbounds function next(z::Zip2, st)
    n1 = next(z.a,st[1])
    n2 = next(z.b,st[2])
    return ((n1[1], n2[1]), (n1[2], n2[2]))
end
@inline done(z::Zip2, st) = done(z.a,st[1]) | done(z.b,st[2])

IteratorSize(::Type{Zip2{I1,I2}}) where {I1,I2} = zip_iteratorsize(IteratorSize(I1),IteratorSize(I2))
IteratorEltype(::Type{Zip2{I1,I2}}) where {I1,I2} = and_iteratoreltype(IteratorEltype(I1),IteratorEltype(I2))

struct Zip{I, Z<:AbstractZipIterator} <: AbstractZipIterator
    a::I
    z::Z
end

"""
    zip(iters...)

For a set of iterable objects, return an iterable of tuples, where the `i`th tuple contains
the `i`th component of each input iterable.

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

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

julia> c = zip(a,b)
Base.Iterators.Zip2{UnitRange{Int64},Array{String,1}}(1:5, ["e", "d", "b", "c", "a"])

julia> length(c)
5

julia> first(c)
(1, "e")
```
"""
zip(a, b, c...) = Zip(a, zip(b, c...))
length(z::Zip) = _min_length(z.a, z.z, IteratorSize(z.a), IteratorSize(z.z))
size(z::Zip) = promote_shape(size(z.a), size(z.z))
axes(z::Zip) = promote_shape(axes(z.a), axes(z.z))
eltype(::Type{Zip{I,Z}}) where {I,Z} = tuple_type_cons(eltype(I), eltype(Z))
@inline start(z::Zip) = tuple(start(z.a), start(z.z))
@propagate_inbounds function next(z::Zip, st)
    n1 = next(z.a, st[1])
    n2 = next(z.z, st[2])
    (tuple(n1[1], n2[1]...), (n1[2], n2[2]))
end
@inline done(z::Zip, st) = done(z.a,st[1]) | done(z.z,st[2])

IteratorSize(::Type{Zip{I1,I2}}) where {I1,I2} = zip_iteratorsize(IteratorSize(I1),IteratorSize(I2))
IteratorEltype(::Type{Zip{I1,I2}}) where {I1,I2} = and_iteratoreltype(IteratorEltype(I1),IteratorEltype(I2))

reverse(z::Zip1) = Zip1(reverse(z.a))
reverse(z::Zip2) = Zip2(reverse(z.a), reverse(z.b))
reverse(z::Zip) = Zip(reverse(z.a), reverse(z.z))

# 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),Array{Int64,1}}(isodd, [1, 2, 3, 4, 5])

julia> foreach(println, f)
1
3
5
```
"""
filter(flt, itr) = Filter(flt, itr)

start(f::Filter) = start_filter(f.flt, f.itr)
function start_filter(pred, itr)
    s = start(itr)
    while !done(itr,s)
        v,t = next(itr,s)
        if pred(v)
            return (false, v, t)
        end
        s=t
    end
    (true,)
end

next(f::Filter, s) = advance_filter(f.flt, f.itr, s)
function advance_filter(pred, itr, st)
    _, v, s = st
    while !done(itr,s)
        w,t = next(itr,s)
        if pred(w)
            return v, (false, w, t)
        end
        s=t
    end
    v, (true, v, s)
end

done(f::Filter, s) = s[1]

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))

# 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`.

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

"""
    peel(iter)

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

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

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

julia> collect(rest)
2-element Array{Char,1}:
 'b'
 'c'
```
"""
function peel(itr)
    s = start(itr)
    done(itr, s) && throw(BoundsError())
    val, s = next(itr, s)
    val, rest(itr, s)
end

start(i::Rest) = i.st
@propagate_inbounds next(i::Rest, st) = next(i.itr, st)
done(i::Rest, st) = done(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{S<:Number}
    start::S
    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::Number, step::Number) = Count(promote(start, step)...)
countfrom(start::Number)               = Count(start, oneunit(start))
countfrom()                            = Count(1, 1)

eltype(::Type{Count{S}}) where {S} = S

start(it::Count) = it.start
next(it::Count, state) = (state, state + it.step)
done(it::Count, state) = false

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`.

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

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

julia> collect(Iterators.take(a,3))
3-element Array{Int64,1}:
 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())

start(it::Take) = (it.n, start(it.xs))

@propagate_inbounds function next(it::Take, state)
    n, xs_state = state
    v, xs_state = next(it.xs, xs_state)
    return v, (n - 1, xs_state)
end

function done(it::Take, state)
    n, xs_state = state
    return n <= 0 || done(it.xs, xs_state)
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 Array{Int64,1}:
  1
  3
  5
  7
  9
 11

julia> collect(Iterators.drop(a,4))
2-element Array{Int64,1}:
  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 start(it::Drop)
    xs_state = start(it.xs)
    for i in 1:it.n
        if done(it.xs, xs_state)
            break
        end

        _, xs_state = next(it.xs, xs_state)
    end
    xs_state
end

@propagate_inbounds next(it::Drop, state) = next(it.xs, state)
done(it::Drop, state) = done(it.xs, state)

# Cycle an iterator forever

struct Cycle{I}
    xs::I
end

"""
    cycle(iter)

An iterator that cycles through `iter` forever.

# 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()

function start(it::Cycle)
    s = start(it.xs)
    return s, done(it.xs, s)
end

function next(it::Cycle, state)
    s, d = state
    if done(it.xs, s)
        s = start(it.xs)
    end
    v, s = next(it.xs, s)
    return v, (s, false)
end

done(it::Cycle, state) = state[2]

reverse(it::Cycle) = Cycle(reverse(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)`).

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

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

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

start(it::Repeated) = nothing
next(it::Repeated, state) = (it.x, nothing)
done(it::Repeated, state) = false

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

reverse(it::Union{Repeated,Take{<:Repeated}}) = 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.

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

IteratorSize(::Type{ProductIterator{Tuple{}}}) = HasShape{0}()
IteratorSize(::Type{ProductIterator{T}}) where {T<:Tuple} =
    prod_iteratorsize( IteratorSize(tuple_type_head(T)), IteratorSize(ProductIterator{tuple_type_tail(T)}) )

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()

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_indices1(t[1], IteratorSize(t[1]))..., _prod_indices(tail(t))...)
_prod_indices1(a, ::HasShape)  = axes(a)
_prod_indices1(a, ::HasLength) = (OneTo(length(a)),)
_prod_indices1(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))
_length(p::ProductIterator) = prod(map(unsafe_length, axes(p)))

IteratorEltype(::Type{ProductIterator{Tuple{}}}) = HasEltype()
IteratorEltype(::Type{ProductIterator{Tuple{I}}}) where {I} = IteratorEltype(I)
function IteratorEltype(::Type{ProductIterator{T}}) where {T<:Tuple}
    I = tuple_type_head(T)
    P = ProductIterator{tuple_type_tail(T)}
    IteratorEltype(I) == EltypeUnknown() ? EltypeUnknown() : IteratorEltype(P)
end

eltype(::Type{<:ProductIterator{I}}) where {I} = _prod_eltype(I)
_prod_eltype(::Type{Tuple{}}) = Tuple{}
_prod_eltype(::Type{I}) where {I<:Tuple} =
    Base.tuple_type_cons(eltype(tuple_type_head(I)),_prod_eltype(tuple_type_tail(I)))

start(::ProductIterator{Tuple{}}) = false
next(::ProductIterator{Tuple{}}, state) = (), true
done(::ProductIterator{Tuple{}}, state) = state

function start(P::ProductIterator)
    iterators = P.iterators
    iter1 = first(iterators)
    state1 = start(iter1)
    d, states, nvalues = _prod_start(tail(iterators))
    d |= done(iter1, state1)
    return (d, (state1, states...), nvalues)
end
function next(P::ProductIterator, state)
    iterators = P.iterators
    d, states, nvalues = state
    iter1 = first(iterators)
    value1, state1 = next(iter1, states[1])
    tailstates = tail(states)
    values = (value1, map(get, nvalues)...) # safe if not done(P, state)
    if done(iter1, state1)
        d, tailstates, nvalues = _prod_next(tail(iterators), tailstates, nvalues)
        if !d # only restart iter1 if not completely done
            state1 = start(iter1)
        end
    end
    return values, (d, (state1, tailstates...), nvalues)
end
done(P::ProductIterator, state) = state[1]

struct MaybeValue{T}
    x::T
    MaybeValue{T}() where {T} = new{T}()
    MaybeValue{T}(x::T) where {T} = new{T}(x)
end

get(v::MaybeValue) = v.x

_prod_start(iterators::Tuple{}) = false, (), ()
function _prod_start(iterators)
    iter1 = first(iterators)
    state1 = start(iter1)
    d, tailstates, tailnvalues = _prod_start(tail(iterators))
    if done(iter1, state1)
        d = true
        nvalue1 = MaybeValue{eltype(iter1)}()
    else
        value1, state1 = next(iter1, state1)
        nvalue1 = MaybeValue{eltype(iter1)}(value1)
    end
    return (d, (state1, tailstates...), (nvalue1, tailnvalues...))
end

_prod_next(iterators::Tuple{}, states, nvalues) = true, (), ()
function _prod_next(iterators, states, nvalues)
    iter1 = first(iterators)
    state1 = first(states)
    if !done(iter1, state1)
        value1, state1 = next(iter1, state1)
        nvalue1 = MaybeValue{eltype(iter1)}(value1)
        return false, (state1, tail(states)...), (nvalue1, tail(nvalues)...)
    else
        d, tailstates, tailnvalues = _prod_next(tail(iterators), tail(states), tail(nvalues))
        if d # all iterators are done
            nvalue1 = MaybeValue{eltype(iter1)}()
        else
            value1, state1 = next(iter1, start(iter1)) # iter cannot be done immediately
            nvalue1 = MaybeValue{eltype(iter1)}(value1)
        end
        return d, (state1, tailstates...), (nvalue1, tailnvalues...)
    end
end

reverse(p::ProductIterator) = ProductIterator(map(reverse, 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 Array{Int64,1}:
 1
 2
 8
 9
```
"""
flatten(itr) = Flatten(itr)

eltype(::Type{Flatten{I}}) where {I} = eltype(eltype(I))
IteratorEltype(::Type{Flatten{I}}) where {I} = _flatteneltype(I, IteratorEltype(I))
_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))

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

function flatten_length(f, T::Type{<:NTuple{N,Any}}) where {N}
    fieldcount(T)*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))

function start(f::Flatten)
    local inner, s2
    s = start(f.it)
    d = done(f.it, s)
    # this is a simple way to make this function type stable
    d && throw(ArgumentError("argument to Flatten must contain at least one iterator"))
    while !d
        inner, s = next(f.it, s)
        s2 = start(inner)
        !done(inner, s2) && break
        d = done(f.it, s)
    end
    return s, inner, s2
end

@propagate_inbounds function next(f::Flatten, state)
    s, inner, s2 = state
    val, s2 = next(inner, s2)
    while done(inner, s2) && !done(f.it, s)
        inner, s = next(f.it, s)
        s2 = start(inner)
    end
    return val, (s, inner, s2)
end

@inline function done(f::Flatten, state)
    s, inner, s2 = state
    return done(f.it, s) && done(inner, s2)
end

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

"""
    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 Array{Array{Int64,1},1}:
 [1, 2]
 [3, 4]
 [5]
```
"""
partition(c::T, n::Integer) where {T} = PartitionIterator{T}(c, Int(n))

struct PartitionIterator{T}
    c::T
    n::Int
end

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

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

function length(itr::PartitionIterator)
    l = length(itr.c)
    return div(l, itr.n) + ((mod(l, itr.n) > 0) ? 1 : 0)
end

start(itr::PartitionIterator) = start(itr.c)

done(itr::PartitionIterator, state) = done(itr.c, state)

function next(itr::PartitionIterator{<:Vector}, state)
    l = state
    r = min(state + itr.n-1, length(itr.c))
    return view(itr.c, l:r), r + 1
end

function next(itr::PartitionIterator, state)
    v = Vector{eltype(itr.c)}(undef, itr.n)
    i = 0
    while !done(itr.c, state) && i < itr.n
        i += 1
        v[i], state = next(itr.c, state)
    end
    return resize!(v, i), state
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`), if iteration is stopped early (e.g. by a `break` in a `for` 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).

# Example:
```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 Array{Char,1}:
 'b'
 'c'
 'd'

julia> collect(a)
2-element Array{Char,1}:
 'e'
 'f'
```

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

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

julia> Base.peek(a)
3

# Sum the remaining elements
julia> sum(a)
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(itr::T) where {T}
        state = start(itr)
        VS = fixpoint_iter_type(T, Union{}, typeof(state))
        new{T, VS}(itr, done(itr, state) ? nothing : next(itr, state)::VS, 0)
    end
end

function reset!(s::Stateful{T,VS}, itr::T) where {T,VS}
    s.itr = itr
    state = start(itr)
    if done(itr, state)
        s.nextvalstate = nothing
    else
        s.nextvalstate = next(itr, state)::VS
    end
    s.taken = 0
    s
end

if Base === Core.Compiler
    fixpoint_iter_type(a, b, c) = 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.
    function fixpoint_iter_type(itrT::Type, valT::Type, stateT::Type)
        nextvalstate = Base._return_type(next, Tuple{itrT, stateT})
        nextvalstate <: Tuple{Any, Any} || return Any
        nextvalstate = Tuple{
            typejoin(valT, fieldtype(nextvalstate, 1)),
            typejoin(stateT, fieldtype(nextvalstate, 2))}
        return (Tuple{valT, stateT} == nextvalstate ? nextvalstate :
            fixpoint_iter_type(itrT,
                fieldtype(nextvalstate, 1),
                fieldtype(nextvalstate, 2)))
    end
end

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

@inline isempty(s::Stateful) = s.nextvalstate === nothing

@inline function popfirst!(s::Stateful)
    vs = s.nextvalstate
    if vs === nothing
        throw(EOFError())
    else
        val, state = vs
        # Until the optimizer can handle setproperty! better here, use explicit setfield!
        setfield!(s, :nextvalstate, done(s.itr, state) ? nothing : next(s.itr, state))
        s.taken += 1
        return val
    end
end

@inline peek(s::Stateful, sentinel=nothing) = s.nextvalstate !== nothing ? s.nextvalstate[1] : sentinel
@inline start(s::Stateful) = nothing
@inline next(s::Stateful, state) = popfirst!(s), nothing
@inline done(s::Stateful, state) = isempty(s)
IteratorSize(::Type{Stateful{VS,T}} where VS) where {T} =
    isa(IteratorSize(T), SizeUnknown) ? SizeUnknown() : HasLength()
eltype(::Type{Stateful{VS, T}} where VS) where {T} = eltype(T)
IteratorEltype(::Type{Stateful{VS,T}} where VS) where {T} = IteratorEltype(T)
length(s::Stateful) = length(s.itr) - s.taken

end


text:
# Wrapper for whatever linear programming we need, in case we don't want to
# use JuMP (e.g., for latency reasons).

module LP

import JuMP
using JuMP: @variable, @objective, @constraint
using MathOptInterface: OPTIMAL
using Clp

mutable struct Model
    LP          :: JuMP.Model
    vars
    base
    shape
    maxval
end
function Model(shape)

    P = JuMP.Model(JuMP.with_optimizer(Clp.Optimizer, LogLevel = 0))
    Y = shape
    n = length(Y)

    @variable(P, x[i=1:n, j=1:Y[i]] >= 0)

    @variable(P, z >= 0)
    @constraint(P, z == sum(x[i, Y[i]] for i=1:n))

    @objective(P, Min, sum(x))

    @constraint(P, [i=1:n, j=1:Y[i]-1], x[i,j] + 1 ≤ x[i,j+1])

    Model(P, x, [], Y, z)

end

set_base!(m::Model, idx) = (m.base = idx)
get_base(m::Model) = m.base

optimize!(m::Model) = JuMP.optimize!(m.LP)

function maybe(func!::Function, m::Model)
    cons = func!(m)
    optimize!(m)
    status = JuMP.termination_status(m.LP)
    for con in cons
        JuMP.delete(m.LP, con)
    end
    return status == OPTIMAL
end

maybe_lt(m::Model, lhs, rhs) = maybe(m->lt!(m, lhs, rhs), m)
maybe_gt(m::Model, lhs, rhs) = maybe(m->gt!(m, lhs, rhs), m)

vars(m::Model, idx) = [m.vars[i, j] for (i,j) in zip(m.base, idx)]

lhs_minus_rhs(m::Model, lhs, rhs) = sum(vars(m, lhs)) - sum(vars(m, rhs))

gt!(m::Model, lhs, rhs) = [@constraint(m.LP, lhs_minus_rhs(m, lhs, rhs) ≥ 1)]
lt!(m::Model, lhs, rhs) = gt!(m, rhs, lhs)
eq!(m::Model, lhs, rhs) = [@constraint(m.LP, lhs_minus_rhs(m, lhs, rhs) == 0)]

value(m::Model, i, j) = JuMP.value(m.vars[i, j])

fix(m::Model, i, j, v) = JuMP.fix(m.vars[i, j], v; force = true)
fixtotal(m::Model, v) = JuMP.fix(m.maxval, v; force = true)

end

ext: jl


text:
using Sundials, Plots
pyplot(size=(800,800))

old = pwd()
ProjDir = dirname(@__FILE__)
cd(ProjDir) #do

  ## f routine. Compute function f(t,y).

  function f(t, y, ydot)
      ydot[1] = -0.04*y[1] + 1.0e4*y[2]*y[3]
      ydot[3] = 3.0e7*y[2]*y[2]
      ydot[2] = -ydot[1] - ydot[3]
  end

  #t = [0.0, 4 * logspace(-1., 7., 9)]
  t = [0.01, 0.1, 1.0, 10.0, 100.0, 1000.0, 10000.0, 100000.0, 1.0e6, 1.0e7]
  res = Sundials.cvode(f, [1.0, 0.1, 0.0], t)
  res |> display
  println()

  p1 = plot(t, res[:,1], color=:darkblue,
    xlabel="x", xscale=:log10, ylabel="y[1]",
    title="Example 7.5 (from Sundials.jl)"
    )
  p2 = plot(t, res[:,2], color=:darkred,
    xlabel="x", xscale=:log10, ylabel="y[2]"
    )
  p3 = plot(t, res[:,3], color=:darkgreen,
    xlabel="x", xscale=:log10, ylabel="y[3]"
    )

  plot(p1, p2, p3, layout=(3, 1))
  savefig("Sundials02.png")
  gui()

  #end


text:
# Helper function to detect if we're currently compiling
currently_compiling() = ccall(:jl_generating_output, Cint, ()) != 0

# Helper function to get the UUID of a module, throwing an error if it can't.
function get_uuid(m::Module)
    uuid = Base.PkgId(m).uuid
    if uuid === nothing
        throw(ArgumentError("Module $(m) does not correspond to a loaded package!"))
    end
    return uuid
end

function find_first_project_with_uuid(uuid::UUID)
    # Find first element in `Base.load_path()` that contains this UUID
    # This code should look similar to the search in `Base.get_preferences()`
    for env in Base.load_path()
        project_toml = Base.env_project_file(env)
        if !isa(project_toml, String)
            continue
        end

        # Check to see if this project has a name mapping
        pkg_name = Base.get_uuid_name(project_toml, uuid)
        if pkg_name !== nothing
            return (project_toml, pkg_name)
        end
    end
    return (nothing, nothing)
end


text:
using AtomicData
using Test
using Unitful
import PhysicalConstants.CODATA2018: h, k_B, c_0

@testset "Read utils" begin
    MgII = get_atomic_stage("Mg", "II")
    Mg2 = get_atomic_stage("Mg", 2)
    Mg2rm = get_atomic_stage("Mg", AtomicData.RomanNumeral(2))
    @test MgII.g == Mg2.g == Mg2rm.g
    @test MgII.χ == Mg2.χ == Mg2rm.χ
    # Reading all neutral stages to check for parsing errors
    function read_all_neutral()
        elements = ["Ag", "Al", "Ar", "As", "Au", "B", "Ba", "Be", "Bi", "Br", "C", "Ca",
                    "Cd", "Ce", "Cl", "Co", "Cr", "Cs", "Cu", "Dy", "Er", "Eu", "F", "Fe",
                    "Ga", "Gd", "Ge", "H", "He", "Hf", "Hg", "Ho", "I", "In", "Ir", "K",
                    "Kr", "La", "Li", "Lu", "Mg", "Mn", "Mo", "N", "Na", "Nb", "Nd", "Ne",
                    "Ni", "O", "Os", "P", "Pb", "Pd", "Pr", "Pt", "Rb", "Re", "Rh", "Ru",
                    "S", "Sb", "Sc", "Se", "Si", "Sm", "Sn", "Sr", "Ta", "Tb", "Te", "Th",
                    "Ti", "Tl", "Tm", "U", "V", "W", "Xe", "Y", "Yb", "Zn", "Zr"]
        for e in elements
            get_atomic_stage(e, 1)
        end
        return true
    end
    @test read_all_neutral()
    @test_throws ArgumentError AtomicData.NIST_wavenumber_to_energy("NOT_NUMBER_1")
    @test_throws ErrorException AtomicData.NIST_wavenumber_to_energy([1, 2])
    @test_throws ErrorException read_NIST("H", 3)
    @test_throws ErrorException get_atomic_stage("Mm", 2)
    @test_throws ErrorException get_atomic_stage("Mg", "II"; source="MY_SOURCE")
end

@testset "Partition function" begin
    temp = 1u"K" / ustrip(k_B)
    @test partition_function([1, 1, 1], [1, 1, 1]u"J", temp) ≈ exp(-1) * 3
    @test partition_function([1, 1, 1], [1, 1, 1]u"J", 0u"K") ≈ 0
    @test partition_function([0, 0, 0], [1, 1, 1]u"J", temp) ≈ 0
    @test partition_function([1, 1, 1], [0, 0, 0]u"J", temp) ≈ 3
    SchI = AtomicStage("Sch", 1, [1, 1, 1], [0, 0, 0]u"J", 0u"J")
    @test partition_function(SchI, temp) ≈ 3
    AlI = get_atomic_stage("Al", "I")
    temp = [5000, 10000]u"K"
    itp = partition_function_interpolator(AlI, temp)
    @test itp == partition_function_interpolator("Al", "I", temp)
    @test itp(5000u"K") == partition_function(AlI, 5000u"K")
    @test itp(10000u"K") == partition_function(AlI, 10000u"K")
    @test_throws ErrorException partition_function_interpolator("Al", "I", temp;
                                                                source="MY_SOURCE")
end

@testset "Roman Numerals" begin
    @test convert(Bool, AtomicData.RomanNumeral(10)) == true
    @test convert(Float64, AtomicData.RomanNumeral("II")) == 2.0
    @test promote(AtomicData.RomanNumeral("V"), 1) == (5, 1)
    @test length(AtomicData.RomanNumeral("X")) == 1
    @test hash(AtomicData.RomanNumeral("L")) == xor(hash("L"), hash(50))
end


text:
using Gaston

# plot example
t = 0:0.0001:.15
carrier = cos(2pi*t*200)
modulator = 0.7+0.5*cos(2pi*t*15)
am = carrier.*modulator
plot(t,am,"color","blue","legend","AM DSB-SC","linewidth",1.5,
    t,modulator,"color","black","legend","Envelope",
    t,-modulator,"color","black","title","AM DSB-SC example",
    "xlabel","Time (s)","ylabel","Amplitude",
    "box","horizontal top left")
set_filename("plotex.pdf")
set_print_linewidth(2)
set_print_fontsize(12)
set_print_size("6.5in,3.9in")
printfigure()

# histogram example
y = sqrt( randn(1000).^2 + randn(1000).^2 )
histogram(y,"bins",25,"norm",1,"color","blue","title",
    "Rayleigh density (25 bins)")
set_filename("histoex.pdf")
printfigure()

# imagesc example
z = zeros(10,10,3)
z[:,:,1] = 255*rand(10,10)
z[:,:,2] = 128*rand(10,10)+40
z[:,:,3] = 64*rand(10,10)+190
imagesc(z,"title","imagesc demo")
set_filename("imagescex.pdf")
printfigure()

# surf example
run(`sleep 0.5`)
gnuplot_send("set view 67,25")
surf(-3:.1:3,-3:0.1:3,(x,y)->cos(x)*sin(y),"plotstyle",
    "points","xlabel","coord 1","ylabel","coord 2","zlabel","coord 3",
    "title","surf demo","color","blue")
set_filename("surfex.pdf")
printfigure()


text:
using InterpretableModels: getfactor, parambounds, getintercept, featurename
using InterpretableModels: setfeaturenames!, modelasstring, setparams!, hasparams
using InterpretableModels: minlag, actualminlag, maxlag, actualmaxlag, params, horizon
using InterpretableModels: predict, predict!, fit!, isvalidparams, windowspans
using DataFrames

@testset "MovingAverageLITRM" begin

m = MovingAverageLITRM(2, 4;
        minlag = 3, maxlag = 8,
        hasintercept = true, mincoef = -10.0, maxcoef = 97.5)

@test !isclassificationmodel(m)
@test  isregressionmodel(m)
@test  istimeseriesmodel(m)
@test !hasparams(m)
@test minlag(m) == 3
@test actualminlag(m) == 3
@test horizon(m) == 3
@test maxlag(m) == 8
@test actualmaxlag(m) == 8

pserr1 = [1.0, 2.24, 5.5, 2.1, 4.1, 1.2, 8.0, 0.2, 0.6]
@test !isvalidparams(m, pserr1) # Has a window length < 1.0 (0.2)

pserr2 = [1.0, 2.24, 5.5, 0.1, 4.1, 1.2, 8.0, 1.2, 0.6]
@test !isvalidparams(m, pserr2)

ps = [1.0, 2.24, 5.5, 2.1, 4.1, 1.2, 8.0, 3.2, 0.6]
@test windowspans(m, ps) == [(5, 2, 6), (8, 1, 8)]
@test isvalidparams(m, ps)
setparams!(m, ps)
@test hasparams(m)
@test minlag(m) == 3
@test actualminlag(m) == 5 # Since 5.5 < 8.0
@test actualmaxlag(m) == 8 # Since 8.0 > 5.5
@test horizon(m) == 5
@test params(m) == ps

c, f, ws, wlen, we = getfactor(m, 1)
@test c == 1.0
@test isa(f, Int)
@test f == 2
@test isa(ws, Int)
@test ws == 5
@test isa(wlen, Int)
@test wlen == 2
@test isa(we, Int)
@test we == 6

c, f, ws, wlen, we = getfactor(m, 2)
@test c == 4.1
@test isa(f, Int)
@test f == 1
@test isa(ws, Int)
@test ws == 8
@test isa(wlen, Int)
@test wlen == 1
@test isa(we, Int)
@test we == 8

@test getintercept(m) == 0.6

@testset "parambounds" begin
    bs = parambounds(m)
    @test length(bs) == (4*2+1) # 4 values per factor and 1 for intercept
    for fi in 1:2
        o = (fi-1)*4
        @test bs[1+o] == (-10.0, 97.5)
        @test bs[2+o][1] == 1.0  # feature 1 is lowest possible feature number
        @test isapprox(bs[2+o][2], 1.0+4.0, atol=1e-4)  # Should almost be 1+nfeatures
        @test bs[3+o][1] == 3.0                         # minlag
        @test isapprox(bs[3+o][2], 1.0+8.0, atol=1e-4)  # maxlag
        @test bs[4+o][1] == 1.0                         # minlag
        @test isapprox(bs[4+o][2], 2.0+8.0-3.0, atol=1e-4)
    end
    @test bs[end] == (-10.0, 97.5)
end # @testset "parambounds" begin

@testset "featurenames and setfeaturenames!" begin
    @test featurename(m, 1) == "feature_1"
    @test featurename(m, 2) == "feature_2"
    setfeaturenames!(m, ["arn", "be"])
    @test featurename(m, 1) == "arn"
    @test featurename(m, 2) == "be"
end # @testset "featurenames and setfeaturenames!" begin

@testset "modelasstring" begin
    setparams!(m, [1.2, 1.6, 5.8, 2.1, 6.987, 2.3, 7.3, 1.4, 0.762])
    @test modelasstring(m) == "0.76 + 1.2 * MAvg(arn[T-5], 2) + 6.99 * be[T-7]"
    df = DataFrame(MyFactor1 = Int[], AnotherFac = Float64[])
    @test modelasstring(m, df) == "0.76 + 1.2 * MAvg(MyFactor1[T-5], 2) + 6.99 * AnotherFac[T-7]"
end

@testset "predict, predict! and fit!" begin
    m = MovingAverageLITRM(2, 3;
            minlag = 1, maxlag = 3,
            hasintercept = true, mincoef = -10.0, maxcoef = 97.5)

    df = DataFrame(
        FacA = Int[    1,    2,   3,    4,    5],
        Fac2 = Float64[3.14, 5.4, 9.8, 10.0, 25.0],
        Fac3 = Float64[4.00, 6.0, 8.0, 10.0, 12.0])

    X = Array{Float64}(df)

    trueparams = [-2.5, 1.0, 1.0, 2.0, 4.0, 3.0, 2.0, 1.0, 1.34]
    setparams!(m, trueparams)
    # Since we have a lag range from 1-2 we can only get three predictions
    # from this model, i.e. for rows 3, 4 and 5.
    # True expression: y = 1.34 .+ MAVG(X[T-1, 1], 2) .* (-2.5) .+ X[T-2, 3] .* 4.0
    y = [Inf, Inf,
            1.34 + (-2.5)*((X[3-1, 1] + X[3-2, 1])/2) + 4.0*X[3-2, 3],
            1.34 + (-2.5)*((X[4-1, 1] + X[4-2, 1])/2) + 4.0*X[4-2, 3],
            1.34 + (-2.5)*((X[5-1, 1] + X[5-2, 1])/2) + 4.0*X[5-2, 3]]
    yhat = predict(m, X, trueparams)
    @test length(yhat) == 3
    @test yhat[1] == y[3]
    @test yhat[2] == y[4]
    @test yhat[3] == y[5]

    yhatpreds = zeros(Float64, 3)
    predict!(m, yhatpreds, X, trueparams)
    @test yhatpreds[1] == y[3]
    @test yhatpreds[2] == y[4]
    @test yhatpreds[3] == y[5]

    # If we don't have the true params we can fit the model to data to recover them.
    dummyparams = [3.5, 1.5, 1.5, 2.1, 2.1, 2.0, 1.7, 1.6, 1.34]
    setparams!(m, dummyparams) # Set a random param vector to ensure a better one is set after fitting
    fit!(m, X, y; MaxTime = 5.0)
    @test length(params(m)) == (1+2*4)
    @test params(m) != dummyparams

    # Since the problem is underconstrained we cannot expect to find the "true"
    # params since there are many that works...
    yhat = predict(m, X)
    @test isapprox(yhat[end-1], y[4], atol=1e-1)
    @test isapprox(yhat[end], y[5], atol=1e-1)
end

end # @testset "MovingAverageLITRM" begin


text:
Refpath="/tmp/benchC.txt"
Tarpath="/tmp/benchP.txt"
run(`bash -c "cmp --silent $Refpath $Tarpath || echo \"files are different\""`)

using DelimitedFiles
using Plots

img = readdlm(Tarpath, ',', Int)
M,N = img |> size
p=heatmap(1:N,1:M,img,aspect_ratio=1, size=(400,400))
savefig(p, joinpath("output", "mandP.svg")) #hide


text:
#=
cudatest:
- Julia version:
- Author: jerrymei
- Date: 2018-05-29
=#
using CUDAdrv
using CUDAnative
import NVTX

NVTX.stop() # Use `--profile-from-start=off`
ITERATIONS = 5

include("conv_2d_cuda.jl")

img=CuArray(rand(Float32,64,144,144,64))
fil_trans=CuArray(rand(Float32,64,3,3,64))
dest=CuArray(rand(Float32,64,144,144,64))
shmem = sizeof(Float32) * ((4 * 4) + 32 * (12 * 12))

# CUDAnative.@device_code_typed @cuda (10,10,shmem) fil_trans_kernel(fil,fil_trans)

function benchmark_fil_trans()
    NVTX.@range "kernel" begin
        @cuda blocks=(2,18 * 18,64) threads=(32,4,4) shmem=shmem maxregs=64 img_trans_kernel(img,fil_trans,dest,18)
#         @cuda blocks=(2,18 * 18,64) threads=(32,4,4) shmem=shmem maxregs=64 img_trans_kernel(img,fil_trans,dest,144,144,64)
#         @cuda blocks=(72 * 72,2,64) threads=32 shmem=shmem img_trans_kernel(img,fil_trans,dest,144,144,64,64)
#         @cuda blocks=(496,64,64) threads=(32,4,4) shmem=shmem img_trans_kernel(img,fil_trans,dest,(24,24),(144,144),64)
        CUDAdrv.synchronize()
    end
end

#warmup
NVTX.mark("Warmup benchmark")
benchmark_fil_trans()

NVTX.@activate CUDAdrv.@profile begin
    for i in 1:ITERATIONS
        benchmark_fil_trans()
    end
end

max_stars_repo_name: JerryLingjieMei/WinogradConvolution

max_issues_repo_name: JerryLingjieMei/WinogradConvolution

max_forks_repo_name: JerryLingjieMei/WinogradConvolution


text:
# check sparse matrix construction
@test isequal(full(sparse(complex(ones(5,5),ones(5,5)))), complex(ones(5,5),ones(5,5)))

# check matrix operations
se33 = speye(3)
do33 = ones(3)
@test isequal(se33 * se33, se33)

# check sparse binary op
@test all(full(se33 + convert(SparseMatrixCSC{Float32,Int32}, se33)) == 2*eye(3))
@test all(full(se33 * convert(SparseMatrixCSC{Float32,Int32}, se33)) == eye(3))

# check horiz concatenation
@test all([se33 se33] == sparse([1, 2, 3, 1, 2, 3], [1, 2, 3, 4, 5, 6], ones(6)))

# check vert concatenation
@test all([se33; se33] == sparse([1, 4, 2, 5, 3, 6], [1, 1, 2, 2, 3, 3], ones(6)))

# check h+v concatenation
se44 = speye(4)
sz42 = spzeros(4, 2)
sz41 = spzeros(4, 1)
sz34 = spzeros(3, 4)
se77 = speye(7)
@test all([se44 sz42 sz41; sz34 se33] == se77)

# check blkdiag concatenation
@test all(blkdiag(se33, se33) == sparse([1, 2, 3, 4, 5, 6], [1, 2, 3, 4, 5, 6], ones(6)))

# check concatenation promotion
sz41_f32 = spzeros(Float32, 4, 1)
se33_i32 = speye(Int32, 3, 3)
@test all([se44 sz42 sz41_f32; sz34 se33_i32] == se77)

# check mixed sparse-dense concatenation
sz33 = spzeros(3, 3)
de33 = eye(3)
@test  all([se33 de33; sz33 se33] == full([se33 se33; sz33 se33 ]))

# check splicing + concatenation on
# random instances, with nested vcat
# also side-checks sparse ref
for i = 1 : 10
    a = sprand(5, 4, 0.5)
    @test all([a[1:2,1:2] a[1:2,3:4]; a[3:5,1] [a[3:4,2:4]; a[5,2:4]]] == a)
end

# sparse ref
a116 = reshape(1:16, 4, 4)
s116 = sparse(a116)
p = [4, 1, 2, 3, 2]
@test full(s116[p,:]) == a116[p,:]
@test full(s116[:,p]) == a116[:,p]
@test full(s116[p,p]) == a116[p,p]

# sparse assign
p = [4, 1, 3]
a116[p, p] = -1
s116[p, p] = -1
@test a116 == s116

p = [2, 1, 4]
a116[p, p] = reshape(1:9, 3, 3)
s116[p, p] = reshape(1:9, 3, 3)
@test a116 == s116

# matrix-vector multiplication (non-square)
for i = 1:5
    a = sprand(10, 5, 0.5)
    b = rand(5)
    @test maximum(abs(a*b - full(a)*b)) < 100*eps()
end

# complex matrix-vector multiplication and left-division
for i = 1:5
    a = speye(5) + 0.1*sprandn(5, 5, 0.2)
    b = randn(5,3) + im*randn(5,3)
    @test (maximum(abs(a*b - full(a)*b)) < 100*eps())
    @test (maximum(abs(a'b - full(a)'b)) < 100*eps())
    @test (maximum(abs(a\b - full(a)\b)) < 1000*eps())
    @test (maximum(abs(a'\b - full(a')\b)) < 1000*eps())
    @test (maximum(abs(a.'\b - full(a.')\b)) < 1000*eps())

    a = speye(5) + 0.1*sprandn(5, 5, 0.2) + 0.1*im*sprandn(5, 5, 0.2)
    b = randn(5,3)
    @test (maximum(abs(a*b - full(a)*b)) < 100*eps())
    @test (maximum(abs(a'b - full(a)'b)) < 100*eps())
    @test (maximum(abs(a\b - full(a)\b)) < 1000*eps())
    @test (maximum(abs(a'\b - full(a')\b)) < 1000*eps())
    @test (maximum(abs(a.'\b - full(a.')\b)) < 1000*eps())

    a = speye(5) + tril(0.1*sprandn(5, 5, 0.2))
    b = randn(5,3) + im*randn(5,3)
    @test (maximum(abs(a*b - full(a)*b)) < 100*eps())
    @test (maximum(abs(a'b - full(a)'b)) < 100*eps())
    @test (maximum(abs(a\b - full(a)\b)) < 1000*eps())
    @test (maximum(abs(a'\b - full(a')\b)) < 1000*eps())
    @test (maximum(abs(a.'\b - full(a.')\b)) < 1000*eps())

    a = speye(5) + tril(0.1*sprandn(5, 5, 0.2) + 0.1*im*sprandn(5, 5, 0.2))
    b = randn(5,3)
    @test (maximum(abs(a*b - full(a)*b)) < 100*eps())
    @test (maximum(abs(a'b - full(a)'b)) < 100*eps())
    @test (maximum(abs(a\b - full(a)\b)) < 1000*eps())
    @test (maximum(abs(a'\b - full(a')\b)) < 1000*eps())
    @test (maximum(abs(a.'\b - full(a.')\b)) < 1000*eps())

    a = speye(5) + triu(0.1*sprandn(5, 5, 0.2))
    b = randn(5,3) + im*randn(5,3)
    @test (maximum(abs(a*b - full(a)*b)) < 100*eps())
    @test (maximum(abs(a'b - full(a)'b)) < 100*eps())
    @test (maximum(abs(a\b - full(a)\b)) < 1000*eps())
    @test (maximum(abs(a'\b - full(a')\b)) < 1000*eps())
    @test (maximum(abs(a.'\b - full(a.')\b)) < 1000*eps())

    a = speye(5) + triu(0.1*sprandn(5, 5, 0.2) + 0.1*im*sprandn(5, 5, 0.2))
    b = randn(5,3)
    @test (maximum(abs(a*b - full(a)*b)) < 100*eps())
    @test (maximum(abs(a'b - full(a)'b)) < 100*eps())
    @test (maximum(abs(a\b - full(a)\b)) < 1000*eps())
    @test (maximum(abs(a'\b - full(a')\b)) < 1000*eps())
    @test (maximum(abs(a.'\b - full(a.')\b)) < 1000*eps())

    a = speye(5) + triu(0.1*sprandn(5, 5, 0.2))
    b = randn(5,3) + im*randn(5,3)
    @test (maximum(abs(a*b - full(a)*b)) < 100*eps())
    @test (maximum(abs(a'b - full(a)'b)) < 100*eps())
    @test (maximum(abs(a\b - full(a)\b)) < 1000*eps())
    @test (maximum(abs(a'\b - full(a')\b)) < 1000*eps())
    @test (maximum(abs(a.'\b - full(a.')\b)) < 1000*eps())

    a = spdiagm(randn(5)) + im*spdiagm(randn(5))
    b = randn(5,3)
    @test (maximum(abs(a*b - full(a)*b)) < 100*eps())
    @test (maximum(abs(a'b - full(a)'b)) < 100*eps())
    @test (maximum(abs(a\b - full(a)\b)) < 1000*eps())
    @test (maximum(abs(a'\b - full(a')\b)) < 1000*eps())
    @test (maximum(abs(a.'\b - full(a.')\b)) < 1000*eps())

    b = randn(5,3) + im*randn(5,3)
    @test (maximum(abs(a*b - full(a)*b)) < 100*eps())
    @test (maximum(abs(a'b - full(a)'b)) < 100*eps())
    @test (maximum(abs(a\b - full(a)\b)) < 1000*eps())
    @test (maximum(abs(a'\b - full(a')\b)) < 1000*eps())
    @test (maximum(abs(a.'\b - full(a.')\b)) < 1000*eps())
end

# matrix multiplication and kron
for i = 1:5
    a = sprand(10, 5, 0.7)
    b = sprand(5, 10, 0.3)
    @test maximum(abs(a*b - full(a)*full(b))) < 100*eps()
    @test full(kron(a,b)) == kron(full(a), full(b))
end

# reductions
@test sum(se33)[1] == 3.0
@test sum(se33, 1) == [1.0 1.0 1.0]
@test sum(se33, 2) == [1.0 1.0 1.0]'
@test prod(se33)[1] == 0.0
@test prod(se33, 1) == [0.0 0.0 0.0]
@test prod(se33, 2) == [0.0 0.0 0.0]'

# spdiagm
@test full(spdiagm((ones(2), ones(2)), (0, -1), 3, 3)) ==
                       [1.0  0.0  0.0; 1.0  1.0  0.0;  0.0  1.0  0.0]

# elimination tree
## upper triangle of the pattern test matrix from Figure 4.2 of
## "Direct Methods for Sparse Linear Systems" by Tim Davis, SIAM, 2006
rowval = int32([1,2,2,3,4,5,1,4,6,1,7,2,5,8,6,9,3,4,6,8,10,3,5,7,8,10,11])
colval = int32([1,2,3,3,4,5,6,6,6,7,7,8,8,8,9,9,10,10,10,10,10,11,11,11,11,11,11])
A = sparse(rowval, colval, ones(length(rowval)))
P,post = Base.LinAlg.etree(A, true)
@test P == int32([6,3,8,6,8,7,9,10,10,11,0])
@test post == int32([2,3,5,8,1,4,6,7,9,10,11])

# issue #4986, reinterpret
sfe22 = speye(Float64, 2)
mfe22 = eye(Float64, 2)
@test reinterpret(Int64, sfe22) == reinterpret(Int64, mfe22)

# issue #5190
@test_throws DimensionMismatch sparsevec([3,5,7],[0.1,0.0,3.2],4)

# issue #5169
@test nnz(sparse([1,1],[1,2],[0.0,-0.0])) == 0

# issue #5386
K,J,V = findnz(SparseMatrixCSC(2,1,[1,3],[1,2],[1.0,0.0]))
@test length(K) == length(J) == length(V) == 1

# issue #5437
@test nnz(sparse([1,2,3],[1,2,3],[0.0,1.0,2.0])) == 2

# issue #5824
@test sprand(4,5,0.5).^0 == sparse(ones(4,5))

# issue #5985
@test sprandbool(4, 5, 0.0) == sparse(zeros(Bool, 4, 5))
@test sprandbool(4, 5, 1.00) == sparse(ones(Bool, 4, 5))
sprb45 = sprandbool(4, 5, 0.5)
@test length(sprb45) == 20
@test 4 <= sum(sprb45)[1] <= 16

# issue #5853, sparse diff
for i=1:2, a={[1 2 3], [1 2 3]', speye(3)}
    @test all(diff(sparse(a),i) == diff(a,i))
end

# test for "access to undefined error" types that initially allocate elements as #undef
@test all(sparse(1:2, 1:2, Number[1,2])^2 == sparse(1:2, 1:2, [1,4]))
sd1 = diff(sparse([1,1,1], [1,2,3], Number[1,2,3]), 1)

# issue #6036
P = spzeros(Float64, 3, 3)
for i = 1:3
    P[i,i] = i
end

@test minimum(P) === 0.0
@test maximum(P) === 3.0
@test minimum(-P) === -3.0
@test maximum(-P) === 0.0

@test maximum(P, (1,)) == [1.0 2.0 3.0]
@test maximum(P, (2,)) == reshape([1.0,2.0,3.0],3,1)
@test maximum(P, (1,2)) == reshape([3.0],1,1)

@test maximum(sparse(-ones(3,3))) == -1
@test minimum(sparse(ones(3,3))) == 1

# Unary functions
a = sprand(5,15, 0.5)
afull = full(a)
for op in (:sin, :cos, :tan, :iceil, :ifloor, :ceil, :floor, :abs, :abs2)
    @eval begin
        @test ($op)(afull) == full($(op)(a))
    end
end

# getindex tests
ni = 23
nj = 32
a116 = reshape(1:(ni*nj), ni, nj)
s116 = sparse(a116)

ad116 = diagm(diag(a116))
sd116 = sparse(ad116)

for (aa116, ss116) in [(a116, s116), (ad116, sd116)]
    ij=11; i=3; j=2
    @test ss116[ij] == aa116[ij]
    @test ss116[(i,j)] == aa116[i,j]
    @test ss116[i,j] == aa116[i,j]
    @test ss116[i-1,j] == aa116[i-1,j]
    ss116[i,j] = 0
    @test ss116[i,j] == 0
    ss116 = sparse(aa116)

    # range indexing
    @test full(ss116[i,:]) == aa116[i,:]
    @test full(ss116[:,j]) == aa116[:,j]'' # sparse matrices/vectors always have ndims==2:
    @test full(ss116[i,1:2:end]) == aa116[i,1:2:end]
    @test full(ss116[1:2:end,j]) == aa116[1:2:end,j]''
    @test full(ss116[i,end:-2:1]) == aa116[i,end:-2:1]
    @test full(ss116[end:-2:1,j]) == aa116[end:-2:1,j]''
    # float-range indexing is not supported

    # sorted vector indexing
    @test full(ss116[i,[3:2:end-3]]) == aa116[i,[3:2:end-3]]
    @test full(ss116[[3:2:end-3],j]) == aa116[[3:2:end-3],j]''
    @test full(ss116[i,[end-3:-2:1]]) == aa116[i,[end-3:-2:1]]
    @test full(ss116[[end-3:-2:1],j]) == aa116[[end-3:-2:1],j]''

    # unsorted vector indexing with repetition
    p = [4, 1, 2, 3, 2, 6]
    @test full(ss116[p,:]) == aa116[p,:]
    @test full(ss116[:,p]) == aa116[:,p]
    @test full(ss116[p,p]) == aa116[p,p]

    # bool indexing
    li = randbool(size(aa116,1))
    lj = randbool(size(aa116,2))
    @test full(ss116[li,j]) == aa116[li,j]''
    @test full(ss116[li,:]) == aa116[li,:]
    @test full(ss116[i,lj]) == aa116[i,lj]
    @test full(ss116[:,lj]) == aa116[:,lj]
    @test full(ss116[li,lj]) == aa116[li,lj]
end

# workaround issue #7197: comment out let-block
#let S = SparseMatrixCSC(3, 3, Uint8[1,1,1,1], Uint8[], Int64[])
S1290 = SparseMatrixCSC(3, 3, Uint8[1,1,1,1], Uint8[], Int64[])
    S1290[1,1] = 1
    S1290[5] = 2
    S1290[end] = 3
    @test S1290[end] == (S1290[1] + S1290[2,2])
    @test 6 == sum(diag(S1290))
    @test (full(S1290)[[3,1],1])'' == full(S1290[[3,1],1])
# end

# setindex tests
let a = spzeros(Int, 10, 10)
    @test countnz(a) == 0
    a[1,:] = 1
    @test countnz(a) == 10
    @test a[1,:] == sparse(ones(Int,1,10))
    a[:,2] = 2
    @test countnz(a) == 19
    @test a[:,2] == 2*sparse(ones(Int,10,1))

    a[1,:] = 1:10
    @test a[1,:] == sparse([1:10]')
    a[:,2] = 1:10
    @test a[:,2] == sparse([1:10])
end

let A = spzeros(Int, 10, 20)
    A[1:5,1:10] = 10
    A[1:5,1:10] = 10
    @test countnz(A) == 50
    @test A[1:5,1:10] == 10 * ones(Int, 5, 10)
    A[6:10,11:20] = 0
    @test countnz(A) == 50
    A[6:10,11:20] = 20
    @test countnz(A) == 100
    @test A[6:10,11:20] == 20 * ones(Int, 5, 10)
    A[4:8,8:16] = 15
    @test countnz(A) == 121
    @test A[4:8,8:16] == 15 * ones(Int, 5, 9)
end

let ASZ = 1000, TSZ = 800
    A = sprand(ASZ, 2*ASZ, 0.0001)
    B = copy(A)
    nA = countnz(A)
    x = A[1:TSZ, 1:(2*TSZ)]
    nx = countnz(x)
    A[1:TSZ, 1:(2*TSZ)] = 0
    nB = countnz(A)
    @test nB == (nA - nx)
    A[1:TSZ, 1:(2*TSZ)] = x
    @test countnz(A) == nA
    @test A == B
    A[1:TSZ, 1:(2*TSZ)] = 10
    @test countnz(A) == nB + 2*TSZ*TSZ
    A[1:TSZ, 1:(2*TSZ)] = x
    @test countnz(A) == nA
    @test A == B
end

let A = speye(Int, 5), I=[1:10], X=reshape([trues(10), falses(15)],5,5)
    @test A[I] == A[X] == reshape([1,0,0,0,0,0,1,0,0,0], 10, 1)
    A[I] = [1:10]
    @test A[I] == A[X] == reshape([1:10], 10, 1)
end

let S = sprand(50, 30, 0.5, x->int(rand(x)*100)), I = sprandbool(50, 30, 0.2)
    FS = full(S)
    FI = full(I)
    @test sparse(FS[FI]) == S[I] == S[FI]
    @test sum(S[FI]) + sum(S[!FI]) == sum(S)

    sumS1 = sum(S)
    sumFI = sum(S[FI])
    S[FI] = 0
    @test sum(S[FI]) == 0
    sumS2 = sum(S)
    @test (sum(S) + sumFI) == sumS1

    S[FI] = 10
    @test sum(S) == sumS2 + 10*sum(FI)
    S[FI] = 0
    @test sum(S) == sumS2

    S[FI] = [1:sum(FI)]
    @test sum(S) == sumS2 + sum(1:sum(FI))
end

let S = sprand(50, 30, 0.5, x->int(rand(x)*100))
    N = length(S) >> 2
    I = randperm(N) .* 4
    J = randperm(N)
    sumS1 = sum(S)
    sumS2 = sum(S[I])
    S[I] = 0
    @test sum(S) == (sumS1 - sumS2)
    S[I] = J
    @test sum(S) == (sumS1 - sumS2 + sum(J))
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
include("basic.jl")
include("beta_saving.jl")
include("alpha_sampling.jl")
include("tensor.jl")
include("custom_rd.jl")
include("rel_feat.jl")
include("lambda_sampling.jl")
include("solver.jl")
include("sparse_csr.jl")
include("sparsebin_csr.jl")
include("sbm.jl")
include("parallel_latent_basic.jl")

include("macau_hmc.jl")


text:
# In-place update of the initial state distribution.
function update_a!(a::AbstractVector, α::AbstractMatrix, β::AbstractMatrix)
    @argcheck size(α, 1) == size(β, 1)
    @argcheck size(α, 2) == size(β, 2) == size(a, 1)

    K = length(a)
    c = 0.0

    for i in OneTo(K)
        a[i] = α[1, i] * β[1, i]
        c += a[i]
    end

    for i in OneTo(K)
        a[i] /= c
    end
end

# In-place update of the transition matrix.
function update_A!(
    A::AbstractMatrix,
    ξ::AbstractArray,
    α::AbstractMatrix,
    β::AbstractMatrix,
    LL::AbstractMatrix,
)
    @argcheck size(α, 1) == size(β, 1) == size(LL, 1) == size(ξ, 1)
    @argcheck size(α, 2) ==
              size(β, 2) ==
              size(LL, 2) ==
              size(A, 1) ==
              size(A, 2) ==
              size(ξ, 2) ==
              size(ξ, 3)

    T, K = size(LL)

    @inbounds for t in OneTo(T - 1)
        m = vec_maximum(view(LL, t + 1, :))
        c = 0.0

        for i in OneTo(K), j in OneTo(K)
            ξ[t, i, j] = α[t, i] * A[i, j] * exp(LL[t+1, j] - m) * β[t+1, j]
            c += ξ[t, i, j]
        end

        for i in OneTo(K), j in OneTo(K)
            ξ[t, i, j] /= c
        end
    end

    fill!(A, 0.0)

    @inbounds for i in OneTo(K)
        c = 0.0

        for j in OneTo(K)
            for t in OneTo(T - 1)
                A[i, j] += ξ[t, i, j]
            end
            c += A[i, j]
        end

        for j in OneTo(K)
            A[i, j] /= c
        end
    end
end

# In-place update of the observations distributions.
function update_B!(B::AbstractVector, γ::AbstractMatrix, observations, estimator)
    @argcheck size(γ, 1) == size(observations, 1)
    @argcheck size(γ, 2) == size(B, 1)
    K = length(B)
    for i in OneTo(K)
        if sum(γ[:, i]) > 0
            B[i] = estimator(typeof(B[i]), permutedims(observations), γ[:, i])
        end
    end
end

function fit_mle!(
    hmm::AbstractHMM,
    observations;
    display = :none,
    maxiter = 100,
    tol = 1e-3,
    robust = false,
    estimator = fit_mle
)
    @argcheck display in [:none, :iter, :final]
    @argcheck maxiter >= 0

    T, K = size(observations, 1), size(hmm, 1)
    history = EMHistory(false, 0, [])

    # Allocate memory for in-place updates
    c = zeros(T)
    α = zeros(T, K)
    β = zeros(T, K)
    γ = zeros(T, K)
    ξ = zeros(T, K, K)
    LL = zeros(T, K)

    loglikelihoods!(LL, hmm, observations)
    robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))

    forwardlog!(α, c, hmm.a, hmm.A, LL)
    backwardlog!(β, c, hmm.a, hmm.A, LL)
    posteriors!(γ, α, β)

    logtot = sum(c)
    (display == :iter) && println("Iteration 0: logtot = $logtot")

    for it = 1:maxiter
        update_a!(hmm.a, α, β)
        update_A!(hmm.A, ξ, α, β, LL)
        update_B!(hmm.B, γ, observations, estimator)

        # Ensure the "connected-ness" of the states,
        # this prevents case where there is no transitions
        # between two extremely likely observations.
        robust && (hmm.A .+= eps())

        @check isprobvec(hmm.a)
        @check istransmat(hmm.A)

        loglikelihoods!(LL, hmm, observations)
        robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))

        forwardlog!(α, c, hmm.a, hmm.A, LL)
        backwardlog!(β, c, hmm.a, hmm.A, LL)
        posteriors!(γ, α, β)

        logtotp = sum(c)
        (display == :iter) && println("Iteration $it: logtot = $logtotp")

        push!(history.logtots, logtotp)
        history.iterations += 1

        if abs(logtotp - logtot) < tol
            (display in [:iter, :final]) &&
                println("EM converged in $it iterations, logtot = $logtotp")
            history.converged = true
            break
        end

        logtot = logtotp
    end

    if !history.converged
        if display in [:iter, :final]
            println("EM has not converged after $(history.iterations) iterations, logtot = $logtot")
        end
    end

    history
end

mutable struct EMHistory
    converged::Bool
    iterations::Int
    logtots::Vector{Float64}
end

"""
    fit_mle(hmm, observations; ...) -> AbstractHMM

Estimate the HMM parameters using the EM (Baum-Welch) algorithm, with `hmm` as the initial state.

**Keyword Arguments**
- `display::Symbol = :none`: when to display convergence logs, can be set to `:iter` or `:final`.
- `init::Symbol = :none`: if set to `:kmeans` the HMM parameters will be initialized using a K-means clustering.
- `maxiter::Integer = 100`: maximum number of iterations to perform.
- `tol::Integer = 1e-3`: stop the algorithm when the improvement in the log-likelihood is less than `tol`.

**Output**
- `<:AbstractHMM`: a copy of the original HMM with the updated parameters.
"""
function fit_mle(hmm::AbstractHMM, observations; init = :none, kwargs...)
    hmm = copy(hmm)

    if init == :kmeans
        kmeans_init!(hmm, observations, display = get(kwargs, :display, :none))
    end

    history = fit_mle!(hmm, observations; kwargs...)
    hmm, history
end


text:
"""
    InverseMultiquadric <: Kernel

TODO
"""
struct InverseMultiquadric <: Kernel
    # TODO
end


text:
module PostgresORM

  # Exposed types of PostgresORM
  export IEntity, IEntity, Modification

  # Exposed functions of PostgresORM
  export get_orm, create_entity!,create_in_bulk_using_copy,retrieve_entity,
         retrieve_one_entity,update_entity!,update_vector_property!,
         delete_entity,delete_entity_alike,execute_plain_query,
         execute_query_and_handle_result

  using Dates, UUIDs

  function get_orm end

  include("./enums/CRUDType.jl")

  # Add the types used by the package
  # NOTE: We make it available at the root of the module because they are used
  #         by the calling libraries
  include("./model/abstract_types.jl")
  include("./model/Modification.jl")

  module ModificationORM
    using ..PostgresORM
    # using ..PostgresORMUtil
    # using ..Controller
    include("./ModificationORM.jl")
  end

  # This is a failed attempt to put the API in the module.
  # It works, but despite Revise, the web server does not take into account the
  #   changes. Therefore we prefer to have it put in the main module so that we
  #   can easily reload it.
  # module WebApi
  #   using ..PostgresORM
  #   export web_api
  #   include("./web-api-definition.jl")
  # end

  greet() = return ("Hello World from PostgresORM!")

  module PostgresORMUtil
      using ..PostgresORM
      using DataFrames, Dates, TimeZones

      export opendbconn, closedbconn, dict2namedtuple, namedtuple2dict, tovector,
            get_nonmissing_typeof_uniontype, dataframe2vector_of_namedtuples,
            dataframerow2namedtuple, getdictvalues, getpropertiesvalues,
            setpropertiesvalues!, remove_spaces_and_split, diff_dict,
            string2enum, string2vector_of_enums,
            int2enum, enum2int, dictstringkeys2symbol, dictnothingvalues2missing,
            getproperties_asdict, string2zoneddatetime

      # Implementations of 'pluralize' are kept in one module per language
      module Pluralize
        module ENG
          include("./PostgresORMUtil/Pluralize/ENG.jl")
        end
        module FRA
          include("./PostgresORMUtil/Pluralize/FRA.jl")
        end
      end

      include("./PostgresORMUtil/pluralize.jl")
      include("./PostgresORMUtil/utils.jl")

  end # module PostgresORMUtil

  # Provides functions to get information about the database structures
  module SchemaInfo
      include("./SchemaInfo/SchemaInfo-def.jl") # This is only the definition of the
                                                 #   module. See below for the actual
                                                 #   implementation.
  end #module SchemaInfo

  module Controller

    using ..CRUDType, ..PostgresORM
    using ..PostgresORMUtil
    # using .ModificationORM # no need (because ModificationORM is a children module ?)
    using Tables, DataFrames, Query, LibPQ, Dates, UUIDs, TickTock, TimeZones,
          JSON
    using IterTools:imap

    include("./Controller/coreORM.utils.part1.jl")
    include("./Controller/coreORM.create.jl")
    include("./Controller/coreORM.retrieve.jl")
    include("./Controller/coreORM.update.jl")
    include("./Controller/coreORM.delete.jl")
    include("./Controller/coreORM.utils.part2.jl")

    # module ModificationORM
    #   using ..PostgresORM
    #   using ..PostgresORMUtil
    #   using ..Controller
    #   include("./Controller/ModificationORM.jl")
    # end

  end # module Controller

  module Tool
    using LibPQ, StringCases
    using ..PostgresORM
    using ..Controller
    using ..CRUDType
    using ..PostgresORMUtil, ..SchemaInfo
    # using .ModificationORM # no need (because ModificationORM is a children module ?)
    using Tables, DataFrames, Query, LibPQ, Dates, TimeZones, UUIDs, TickTock
    include("./Tool/Tool.jl")
    include("./Tool/Tool-typescript.jl")

  end # module Tool

  # Implementation of the SchemaInfo module
  include("./SchemaInfo/SchemaInfo-imp.jl")

  include("./exposed-functions-from-submodules.jl")

end # ENDIF module PostgresORM


text:
using Plots

function khayyam_triangle(lentingle)
    output = []
    for i in 1:lentingle
        push!(output, [])
        for j in 1:i
            if j == 1
                append!(output[i], 1)
            elseif j == i
                append!(output[i], 1)
            else j != 1 && j != i
                add = output[i - 1][j] + output[i - 1][j - 1]
                append!(output[i], add)
            end
        end
    end
    return output
end

function SPT_KhPT(order)
    order = 2^order
    triangle = khayyam_triangle(order)
    x = []
    y = []
    for i in 1:order
        for j in 1:i
            if triangle[i][j]%2 != 0
                push!(x, j - i/2)
                push!(y, -i)
            end
        end
    end
    return x, y
end

msize = [7,1,0.1]
olist = [5, 7, 10]
for i in 1:3
    x,y = SPT_KhPT(olist[i])
    scatter(x,y,markersize = msize[i], legend = false, border=:none, dpi=300)
    savefig("C:\\Users\\Yaghoub\\Documents\\GitHub\\ComputationalPhysics-Fall2021\\ProblemSet1\\Figs\\Q5\\SPTKP-O$(i).png")
end


text:
module QuadraticModelsGurobi

export gurobi

using Gurobi
using QuadraticModels, SolverCore
using LinearAlgebra, SparseArrays

const gurobi_statuses = Dict(1 => :unknown,
                             2 => :acceptable,
                             3 => :infeasible,
                             4 => :infeasible,
                             5 => :unbounded,
                             6 => :exception,
                             7 => :max_iter,
                             8 => :exception,
                             9 => :max_time,
                             10 => :exception,
                             11 => :user,
                             12 => :exception,
                             13 => :exception,
                             14 => :exception,
                             15 => :exception)

function sparse_csr(I, J, V, m=maximum(I), n=maximum(J))
    csrrowptr = zeros(Int, m+1)
    # Compute the CSR form's row counts and store them shifted forward by one in csrrowptr
    coolen = length(I)
    min(length(J), length(V)) >= coolen || throw(ArgumentError("J and V need length >= length(I) = $coolen"))
    @inbounds for k in 1:coolen
        Ik = I[k]
        if 1 > Ik || m < Ik
            throw(ArgumentError("row indices I[k] must satisfy 1 <= I[k] <= m"))
        end
        csrrowptr[Ik+1] += 1
    end

    # Compute the CSR form's rowptrs and store them shifted forward by one in csrrowptr
    countsum = 1
    csrrowptr[1] = 1
    @inbounds for i in 2:(m+1)
        overwritten = csrrowptr[i]
        csrrowptr[i] = countsum
        countsum += overwritten
    end

    # Counting-sort the column and nonzero values from J and V into csrcolval and csrnzval
    # Tracking write positions in csrrowptr corrects the row pointers
    csrcolval = zeros(Int, length(I))
    csrnzval = zeros(length(I))
    @inbounds for k in 1:coolen
        Ik, Jk = I[k], J[k]
        if 1 > Jk || n < Jk
            throw(ArgumentError("column indices J[k] must satisfy 1 <= J[k] <= n"))
        end
        csrk = csrrowptr[Ik+1]
        csrrowptr[Ik+1] = csrk + 1
        csrcolval[csrk] = Jk
        csrnzval[csrk] = V[k]
    end
    csrrowptr = csrrowptr[1:end-1]
    return csrrowptr, csrcolval, csrnzval
end

function gurobi(QM::QuadraticModel; method=2, kwargs...)
    env = Gurobi.Env()
    # -1=automatic, 0=primal simplex, 1=dual simplex, 2=barrier,
    # 3=concurrent, 4=deterministic concurrent, 5=deterministic concurrent simplex.
    # default to barrier
    GRBsetintparam(env, "Method", method)
    # use kwargs change to presolve, scaling and crossover mode
    # example: gurobi(QM, presolve=0) (see gurobi doc for other options)
    for (k, v) in kwargs
        if k==:presolve
            GRBsetintparam(env, "Presolve", v) # 0 = no presolve
        elseif k==:scaling
            GRBsetintparam(env, "ScaleFlag", v) # 0 = no scaling
        elseif k==:crossover
            GRBsetintparam(env, "Crossover", v) # 0 = no crossover
        elseif k==:display
            GRBsetintparam(env, "OutputFlag", v) # 0 = no display
        end
    end

    model = Ref{Ptr{Cvoid}}()
    GRBnewmodel(env, model, "", QM.meta.nvar, QM.data.c, QM.meta.lvar, QM.meta.uvar, C_NULL, C_NULL)
    GRBsetdblattr(model.x, "ObjCon", QM.data.c0)
    if QM.meta.nnzh > 0
        Hvals = zeros(eltype(QM.data.Hvals), length(QM.data.Hvals))
        for i=1:length(QM.data.Hvals)
            if QM.data.Hrows[i] == QM.data.Hcols[i]
                Hvals[i] = QM.data.Hvals[i] / 2
            else
                Hvals[i] = QM.data.Hvals[i]
            end
        end
        GRBaddqpterms(model.x, length(QM.data.Hcols), convert(Array{Cint,1}, QM.data.Hrows.-1),
                      convert(Array{Cint,1}, QM.data.Hcols.-1), Hvals)
    end

    Acsrrowptr, Acsrcolval, Acsrnzval = sparse_csr(QM.data.Arows,QM.data.Acols,
                                                   QM.data.Avals, QM.meta.ncon,
                                                   QM.meta.nvar)
    GRBaddrangeconstrs(model.x, QM.meta.ncon, length(Acsrcolval), convert(Array{Cint,1}, Acsrrowptr.-1),
                       convert(Array{Cint,1}, Acsrcolval.-1), Acsrnzval, QM.meta.lcon, QM.meta.ucon, C_NULL)

    GRBoptimize(model.x)

    x = zeros(QM.meta.nvar)
    GRBgetdblattrarray(model.x, "X", 0, QM.meta.nvar, x)
    y = zeros(QM.meta.ncon)
    GRBgetdblattrarray(model.x, "Pi", 0, QM.meta.ncon, y)
    s = zeros(QM.meta.nvar)
    GRBgetdblattrarray(model.x, "RC", 0, QM.meta.nvar, s)
    status = Ref{Cint}()
    GRBgetintattr(model.x, "Status", status)
    baritcnt = Ref{Cint}()
    GRBgetintattr(model.x, "BarIterCount", baritcnt)
    objval = Ref{Float64}()
    GRBgetdblattr(model.x, "ObjVal", objval)
    p_feas = Ref{Float64}()
    GRBgetdblattr(model.x, "ConstrResidual", p_feas)
    d_feas = Ref{Float64}()
    GRBgetdblattr(model.x, "DualResidual", d_feas)
    elapsed_time = Ref{Float64}()
    GRBgetdblattr(model.x, "Runtime", elapsed_time)
    stats = GenericExecutionStats(get(gurobi_statuses, status[], :unknown),
                                  QM, solution = x,
                                  objective = objval[],
                                  iter = Int64(baritcnt[]),
                                  primal_feas = p_feas[],
                                  dual_feas = d_feas[],
                                  multipliers = y,
                                  elapsed_time = elapsed_time[])
    return stats
end

end


text:
# algorithm on page 504 of https://people.maths.ox.ac.uk/trefethen/barycentric.pdf
function barycentric_lagrange_weights(x)
    n = length(x) - 1
    w = similar(x)
    w[1] = 1.0
    for j in 1:n
        for k in 0:j-1
            w[k+1] = (x[k+1] - x[j+1]) * w[k+1]
        end
        w[j+1] = prod(0:j-1) do k
            x[j+1] - x[k+1]
        end
    end
    for j in 0:n
        w[j+1] = 1/w[j+1]
    end
    return w
end

# equation 4.2 of https://people.maths.ox.ac.uk/trefethen/barycentric.pdf
function barycentric_lagrange_matrix(xs,xt,w=barycentric_lagrange_weights(xs))
    ns = length(xs)
    nt = length(xt)
    d = similar(xt)
    for i in 1:nt
        d[i] = sum(zip(xs,w)) do (x,w)
            w /(xt[i]-x)
        end
    end
    A= [w[j]/(xt[i]-xs[j]) / d[i] for i in 1:nt, j in 1:ns]
    return A
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
using LyapunovExponents

num_rec = 20
sampling_interval = 5
t_renorm = 1
t_clv = t_renorm * num_rec * sampling_interval
prob = CLVProblem(LyapunovExponents.lorenz_63().prob,
                  t_renorm = t_renorm,
                  t_forward_tran = 4000,
                  t_backward_tran = 4000,
                  t_clv = t_clv)
solver = init(prob)

x_history = [Vector{Float64}(3) for _ in 1:num_rec]
G_history = [Matrix{Float64}(3, 3) for _ in 1:num_rec]
C_history = [Matrix{Float64}(3, 3) for _ in 1:num_rec]

forward = @time forward_dynamics!(solver; progress=1)
@time for (i, G) in indexed_forward_dynamics!(forward)
    k, r = divrem(i, sampling_interval)
    j = k + 1
    if r == 1 && j <= num_rec
        x_history[j] .= phase_state(forward)
        G_history[j] .= G
    end
end

@time for (i, C) in @time indexed_backward_dynamics!(solver; progress=1)
    k, r = divrem(i, sampling_interval)
    j = k + 1
    if r == 1
        C_history[j] .= C
    end
end

CLV_history = [G * C for (G, C) in zip(G_history, C_history)]

using DifferentialEquations
sol = solve(ODEProblem(
    prob.phase_prob.f,
    x_history[1],
    (0.0, t_clv * 3),
    prob.phase_prob.p
))

using Plots
plt = plot(sol, vars=(2, 3), linewidth=0.5, linealpha=0.5, label="")
vec_scale = 3
for (n, (x, V)) in enumerate(zip(x_history, CLV_history))
    for i in 1:3
        plot!(plt,
              [x[2], x[2] + vec_scale * V[2, i]],
              [x[3], x[3] + vec_scale * V[3, i]],
              color = i + 1,
              arrow = 0.4,
              label = (n == 1 ? "CLV$i" : ""))
    end
end
plt


text:
@testset "TFIDF" begin

    doc1 = "a a a sample text text"
    doc2 = "another example example text text"
    doc3 = ""
    doc4 = "another another text text text text"

    # TODO: this should work!
    # crps = Corpus(map(StringDocument, [doc1 doc2 doc3 doc4]))

    crps = Corpus(Any[StringDocument(doc1), StringDocument(doc2), StringDocument(doc3), StringDocument(doc4)])

    update_lexicon!(crps)
    m = DocumentTermMatrix(crps)

    # Terms are in alphabetical ordering
    correctweights =[0.5  0.0  0.0  1/6  1/3
                     0.0  0.2  0.4  0.0  0.4
                     0.0  0.0  0.0  0.0  0.0
                     0.0  1/3  0.0  0.0  2/3]

    myweights = tf(m)
    @test myweights == correctweights

    myweights = tf(dtm(m))
    @test myweights ≈ sparse(correctweights)
    @test typeof(myweights) <: SparseMatrixCSC

    myweights = tf(dtm(m, :dense))
    @test isnan(sum(myweights)) == 0
    @test myweights ≈ correctweights
    @test typeof(myweights) <: Matrix

    myweights = float(dtm(m));
    tf!(myweights)
    @test myweights ≈ correctweights
    @test typeof(myweights) <: SparseMatrixCSC

    myweights = float(dtm(m, :dense));
    tf!(myweights)
    @test myweights ≈ correctweights
    @test typeof(myweights) <: Matrix

    # Terms are in alphabetical ordering
    correctweights = [0.6931471805599453 0.0 0.0 0.23104906018664842 0.09589402415059362
	              0.0 0.13862943611198905 0.5545177444479562 0.0 0.11507282898071235
	              0.0 0.0 0.0 0.0 0.0
	              0.0 0.23104906018664842  0.0 0.0 0.19178804830118723]

    myweights = tf_idf(m)
    @test myweights ≈ correctweights

    myweights = tf_idf(dtm(m))
    @test myweights ≈ correctweights
    @test typeof(myweights) <: SparseMatrixCSC

    myweights = tf_idf(dtm(m, :dense))
    @test isnan(sum(myweights)) == 0
    @test myweights ≈ correctweights
    @test typeof(myweights) <: Matrix

    myweights = float(dtm(m));
    tf_idf!(myweights)
    @test myweights ≈ correctweights
    @test typeof(myweights) <: SparseMatrixCSC

    myweights = float(dtm(m, :dense));
    tf_idf!(myweights)
    @test myweights ≈ correctweights
    @test typeof(myweights) <: Matrix
end


text:
# Sensitivity analysis for c-posterior of autoregressive model of variable order.
module Run
using PyPlot
draw() = (pause(0.001); get_current_fig_manager()[:window][:raise_]())
logsumexp(x) = (m=maximum(x); (m==-Inf ? -Inf : log(sum(exp(x-m)))+m))

# Settings
n = 10^4
alphas = [50,100,500,1200,2000] # robustification param (max sample size)
K = 25 # maximum value of k to use

# Data generation parameters
s0 = 1.0 # std dev of noise
sm = 0.5 # scale of misspecification
a0 = [0.25,0.25,-0.25,0.25] # autoregression coefficients
k0 = length(a0)
f_noise(t) = randn()*s0 + sm*sin(t) # Time-varying noise

# Model parameters
s = s0
sa = 1.0 # std dev of prior on a's

# Code for autoregressive model
include("AR.jl")

figure(1,figsize=(8,3.2)); clf(); hold(true)
subplots_adjust(bottom=0.2)

srand(0) # Reset RNG
x = generate(n,a0,f_noise) # Generate data
for (i_a,alpha) in enumerate(alphas)
    zeta = (1/n)/(1/n + 1/alpha) # robustification param

    # Compute marginal likelihood for each k
    log_m = zeros(K+1)
    for k = 0:K
        log_m[k+1] = log_marginal(x,k,s,sa,zeta)
    end

    # Compute posterior on k, using an improper uniform prior on k
    pk = exp(log_m - logsumexp(log_m))

    # Plot
    colors = "bgyrm"
    shapes = "ds^v*"
    lab = "\$\\alpha = " * (alpha<Inf ? repr(int(alpha)) : "\\infty") * "\$"
    plot(0:K,pk,"$(colors[i_a])$(shapes[i_a])-",label=lab,markersize=8,linewidth=2)
end

title("\$\\mathrm{Sensitivity\\,analysis}\$",fontsize=17)
ylabel("\$\\Pi(k|\\mathrm{data})\$",fontsize=15)
xlabel("\$k \\,\\,\\mathrm{(\\#\\,of\\,coefficients)}\$",fontsize=16)
xlim(0,12)
xticks(0:12)
ylim(0,1)
legend(numpoints=1)
draw()
savefig("AR-sensitivity.png",dpi=150)

end # module

max_stars_repo_name: jwmi/CoarsenedPosterior

max_issues_repo_name: jwmi/CoarsenedPosterior

max_forks_repo_name: jwmi/CoarsenedPosterior


text:
module Blobs

export Blob

using ..Elements
import ..Motions: induce_velocity, mutually_induce_velocity!, self_induce_velocity!, advect,
                  allocate_velocity

#== Type definition ==#

"""
    Blob <: Elements.Element

An immutable structure representing a regularized point source/vortex

## Fields
- `z`: position
- `S`: strength/circulation
- `δ`: regularization radius

## Constructors

"""
struct Blob{S <: Number,R <: Real} <: Element
    z::Complex{R}
    S::Union{R,Complex{R}}
    δ::Float64
    Blob{S}(z::Complex{R}, s::T, δ)           where {S <: Complex, T <: Real, R <: Real} = (U = promote_type(R,T); new{Complex{U},U}(z, im*convert(U,s),δ))
    Blob{S}(z::Complex{R}, s::Complex{T}, δ)  where {S <: Complex, T <: Real, R <: Real} = (U = promote_type(R,T); new{Complex{U},U}(z, convert(Complex{U},s),δ))
    Blob{S}(z::Complex{R}, s::T, δ)           where {S <: Real, T <: Real, R <: Real}    = (U = promote_type(R,T); new{U,U}(z, convert(U,s),δ))
end
# struct Blob{T <: Number, R <: Number} <: Element
#     z::Complex{R}
#     S::T
#     δ::Float64
#     Blob{T,R}(z, s::Real, δ) where {T <: Complex, R <: Number} = new(z, im*s, δ)
#     Blob{T,R}(z, s, δ)       where {T <: Complex, R <: Number} = new(z, s, δ)
#     Blob{T,R}(z, s::Real, δ) where {T <: Real, R <: Number}    = new(z, s, δ)
# end

Elements.kind(::Blob) = Singleton
Elements.kind(::Type{Blob{T,R}}) where {T,R} = Singleton

Elements.property_type(::Blob{T,R}) where {T,R} = R
Elements.property_type(::Type{Blob{T,R}}) where {T,R} = R

#Elements.promote_property_type(::Blob{T,R}) where {T,R} = promote_type(T,R)
#Elements.promote_property_type(::Type{Blob{T,R}}) where {T,R} = promote_type(T,R)

#(b::Blob{T,R})(; z = b.z, S = b.S, δ = b.δ) where {T,R} = Blob{T,R}(z, S, δ)

#== Methods to be extended ==#

Elements.position(b::Blob) = b.z
Elements.blobradius(b::Blob) = b.δ
Elements.streamfunction(z::Complex{T}, b::Blob) where {T} = real(-0.5b.S*log(z - b.z)/π)

Elements.complexpotential(z::Complex{T}, b::Blob) where {T} = -0.5im*b.S*log(z - b.z)/π

blob_kernel(z, δ) = 0.5im*z/(π*(abs2(z) + δ^2))

# ensures that the velocity is of same type as position
allocate_velocity(v::Vector{Blob{T,R}}) where {T,R} = zeros(Complex{R},length(v))

function induce_velocity(z::Complex{T}, b::Blob, t) where {T}
    b.S'*blob_kernel(z - b.z, b.δ)
end

function induce_velocity(target::Blob, source::Blob, t)
    δ = √(0.5(target.δ^2 + source.δ^2))
    source.S'*blob_kernel(target.z - source.z, δ)
end

function mutually_induce_velocity!(ws₁, ws₂,
                                   blobs₁::Vector{Blob{T₁}},
                                   blobs₂::Vector{Blob{T₂}}, t) where {T₁, T₂}
    for (s, source) in enumerate(blobs₁)
        for (t, target) in enumerate(blobs₂)
            δ = √(0.5(target.δ^2 + source.δ^2))
            K = blob_kernel(target.z - source.z, δ)
            ws₂[t] += source.S'*K
            ws₁[s] -= target.S'*K
        end
    end
    nothing
end

function self_induce_velocity!(ws, blobs::Vector{Blob{T,R}}, t) where {T,R}
    N = length(blobs)

    for s in 1:N, t in s+1:N
        δ = √(0.5(blobs[t].δ^2 + blobs[s].δ^2))
        K = blob_kernel(blobs[t].z -blobs[s].z, δ)
        ws[t] += blobs[s].S'*K
        ws[s] -= blobs[t].S'*K
    end
    ws
end

function advect(b::Blob{T}, w::ComplexF64, Δt::Float64) where {T}
    Blob{T}(b.z + w*Δt, b.S, b.δ)
end
#
#function Base.show(io::IO, b::Blob)
#    print(io, "Vortex Blob: z = $(round(b.z, digits=3)), Γ = $(round(b.Γ, digits=3)), δ = $(round(b.δ, digits=3))")
#end

end