text:
@doc doc"""
    ProductManifold{TM<:Tuple, TRanges<:Tuple, TSizes<:Tuple} <: Manifold

Product manifold $M_1 \times M_2 \times \dots \times M_n$ with product geometry.
`TRanges` and `TSizes` statically define the relationship between representation
of the product manifold and representations of point, tangent vectors
and cotangent vectors of respective manifolds.

# Constructor

    ProductManifold(M_1, M_2, ..., M_n)

generates the product manifold $M_1 \times M_2 \times \dots \times M_n$.
Alternatively, the same manifold can be contructed using the `×` operator:
`M_1 × M_2 × M_3`.
"""
struct ProductManifold{TM<:Tuple} <: Manifold
    manifolds::TM
end

ProductManifold(manifolds::Manifold...) = ProductManifold{typeof(manifolds)}(manifolds)

struct ProductMetric <: Metric end

"""
    ProductFVectorDistribution([type::VectorBundleFibers], [x], distrs...)

Generates a random vector at point `x` from vector space (a fiber of a tangent
bundle) of type `type` using the product distribution of given distributions.

Vector space type and `x` can be automatically inferred from distributions `distrs`.
"""
struct ProductFVectorDistribution{
    TSpace<:VectorBundleFibers{<:VectorSpaceType,<:ProductManifold},
    TD<:(NTuple{N,Distribution} where {N}),
    TX,
} <: FVectorDistribution{TSpace,TX}
    type::TSpace
    x::TX
    distributions::TD
end

"""
    ProductPointDistribution(M::ProductManifold, distributions)

Product distribution on manifold `M`, combined from `distributions`.
"""
struct ProductPointDistribution{
    TM<:ProductManifold,
    TD<:(NTuple{N,Distribution} where {N}),
} <: MPointDistribution{TM}
    manifold::TM
    distributions::TD
end

"""
    ProductRetraction(retractions::AbstractRetractionMethod...)

Product retraction of `retractions`. Works on [`ProductManifold`](@ref).
"""
struct ProductRetraction{TR<:Tuple} <: AbstractRetractionMethod
    retractions::TR
end

function ProductRetraction(retractions::AbstractRetractionMethod...)
    return ProductRetraction{typeof(retractions)}(retractions)
end

"""
    InverseProductRetraction(retractions::AbstractInverseRetractionMethod...)

Product inverse retraction of `inverse retractions`. Works on [`ProductManifold`](@ref).
"""
struct InverseProductRetraction{TR<:Tuple} <: AbstractInverseRetractionMethod
    inverse_retractions::TR
end

function InverseProductRetraction(inverse_retractions::AbstractInverseRetractionMethod...)
    return InverseProductRetraction{typeof(inverse_retractions)}(inverse_retractions)
end

"""
    PrecomputedProductOrthonormalBasis(parts::NTuple{N,AbstractPrecomputedOrthonormalBasis} where N, F::AbstractNumbers = ℝ)

A precomputed orthonormal basis of a tangent space of a product manifold.
The tuple `parts` stores bases corresponding to multiplied manifolds.

The type parameter `F` denotes the [`AbstractNumbers`](@ref) that will be used as scalars.
"""
struct PrecomputedProductOrthonormalBasis{
    T<:NTuple{N,AbstractPrecomputedOrthonormalBasis} where {N},
    F,
} <: AbstractPrecomputedOrthonormalBasis{F}
    parts::T
end

function PrecomputedProductOrthonormalBasis(
    parts::NTuple{N,AbstractPrecomputedOrthonormalBasis},
    F::AbstractNumbers = ℝ,
) where {N}
    return PrecomputedProductOrthonormalBasis{typeof(parts),F}(parts)
end

"""
    check_manifold_point(M::ProductManifold, x; kwargs...)

Check whether `x` is a valid point on the [`ProductManifold`](@ref) `M`.

The tolerance for the last test can be set using the `kwargs...`.
"""
function check_manifold_point(M::ProductManifold, x::ProductRepr; kwargs...)
    for t ∈ ziptuples(M.manifolds, submanifold_components(M, x))
        err = check_manifold_point(t...; kwargs...)
        err === nothing || return err
    end
    return nothing
end

function check_manifold_point(M::ProductManifold, x::ProductArray; kwargs...)
    for t ∈ ziptuples(M.manifolds, submanifold_components(M, x))
        err = check_manifold_point(t...; kwargs...)
        err === nothing || return err
    end
    return nothing
end

"""
    check_tangent_vector(M::ProductManifold, x, v; kwargs... )

Check whether `v` is a tangent vector to `x` on the [`ProductManifold`](@ref)
`M`, i.e. atfer [`check_manifold_point`](@ref)`(M, x)`, and all projections to
base manifolds must be respective tangent vectors.

The tolerance for the last test can be set using the `kwargs...`.
"""
function check_tangent_vector(M::ProductManifold, x::ProductRepr, v::ProductRepr; kwargs...)
    perr = check_manifold_point(M, x)
    perr === nothing || return perr
    ts = ziptuples(M.manifolds, submanifold_components(M, x), submanifold_components(M, v))
    for t ∈ ts
        err = check_tangent_vector(t...; kwargs...)
        err === nothing || return err
    end
    return nothing
end
function check_tangent_vector(
    M::ProductManifold,
    x::ProductArray,
    v::ProductArray;
    kwargs...,
)
    perr = check_manifold_point(M, x)
    perr === nothing || return perr
    ts = ziptuples(M.manifolds, submanifold_components(M, x), submanifold_components(M, v))
    for t ∈ ts
        err = check_tangent_vector(t...; kwargs...)
        err === nothing || return err
    end
    return nothing
end

@doc doc"""
    cross(M,N)
    cross(M1, M2, M3,...)

Return the [`ProductManifold`](@ref) For two [`Manifold`](@ref)s `M` and `N`,
where for the case that one of them is a [`ProductManifold`](@ref) itself,
the other is either prepended (if `N` is a product) or appenden (if `M`) is.
If both are product manifold, they are combined into one product manifold,
keeping the order.

For the case that more than one is a product manifold of these is build with the
same approach as above
"""
cross(::Manifold...)
cross(M1::Manifold, M2::Manifold) = ProductManifold(M1, M2)
cross(M1::ProductManifold, M2::Manifold) = ProductManifold(M1.manifolds..., M2)
cross(M1::Manifold, M2::ProductManifold) = ProductManifold(M1, M2.manifolds...)
function cross(M1::ProductManifold, M2::ProductManifold)
    return ProductManifold(M1.manifolds..., M2.manifolds...)
end

function det_local_metric(M::MetricManifold{ProductManifold,ProductMetric}, x::ProductArray)
    dets = map(det_local_metric, M.manifolds, submanifold_components(M, x))
    return prod(dets)
end

@doc doc"""
    distance(M::ProductManifold, x, y)

compute the distance two points `x` and `y` on the [`ProductManifold`](@ref) `M`, which is
the 2-norm of the elementwise distances on the internal manifolds that build `M`.
"""
function distance(M::ProductManifold, x, y)
    return sqrt(sum(
        map(
            distance,
            M.manifolds,
            submanifold_components(M, x),
            submanifold_components(M, y),
        ) .^ 2,
    ))
end

@doc doc"""
    exp(M::ProductManifold, x, v)

compute the exponential map from `x` towards `v` on the [`ProductManifold`](@ref) `M`,
which is the elementwise exponential map on the internal manifolds that build `M`.
"""
exp(::ProductManifold, ::Any...)
function exp(M::ProductManifold, x::ProductRepr, v::ProductRepr)
    return ProductRepr(map(
        exp,
        M.manifolds,
        submanifold_components(M, x),
        submanifold_components(M, v),
    )...)
end

function exp!(M::ProductManifold, y, x, v)
    map(
        exp!,
        M.manifolds,
        submanifold_components(M, y),
        submanifold_components(M, x),
        submanifold_components(M, v),
    )
    return y
end

@doc doc"""
    flat(M::ProductManifold, x, w::FVector{TangentSpaceType})

use the musical isomorphism to transform the tangent vector `w` from the tangent space at
`x` on the [`ProductManifold`](@ref) `M` to a cotangent vector.
This can be done elementwise, so for every entry of `w` (and `x`) sparately
"""
flat(::ProductManifold, ::Any...)

function flat!(M::ProductManifold, v::CoTFVector, x, w::TFVector)
    vfs = map(u -> FVector(CotangentSpace, u), submanifold_components(v))
    wfs = map(u -> FVector(TangentSpace, u), submanifold_components(w))
    map(flat!, M.manifolds, vfs, submanifold_components(M, x), wfs)
    return v
end

function get_basis(M::ProductManifold, x, B::AbstractBasis)
    parts = map(t -> get_basis(t..., B), ziptuples(M.manifolds, submanifold_components(x)))
    return PrecomputedProductOrthonormalBasis(parts)
end
function get_basis(M::ProductManifold, x, B::DiagonalizingOrthonormalBasis)
    vs = map(ziptuples(
        M.manifolds,
        submanifold_components(x),
        submanifold_components(B.v),
    )) do t
        return get_basis(t[1], t[2], DiagonalizingOrthonormalBasis(t[3]))
    end
    return PrecomputedProductOrthonormalBasis(vs)
end
function get_basis(M::ProductManifold, x, B::ArbitraryOrthonormalBasis)
    parts = map(t -> get_basis(t..., B), ziptuples(M.manifolds, submanifold_components(x)))
    return PrecomputedProductOrthonormalBasis(parts)
end

function get_coordinates(M::ProductManifold, x, v, B::PrecomputedProductOrthonormalBasis)
    reps = map(
        get_coordinates,
        M.manifolds,
        submanifold_components(x),
        submanifold_components(v),
        B.parts,
    )
    return vcat(reps...)
end
function get_coordinates(M::ProductManifold, x, v, B::ArbitraryOrthonormalBasis)
    reps = map(
        t -> get_coordinates(t..., B),
        ziptuples(M.manifolds, submanifold_components(x), submanifold_components(v)),
    )
    return vcat(reps...)
end

function get_vector(
    M::ProductManifold{<:NTuple{N,Any}},
    x::ProductRepr,
    v,
    B::PrecomputedProductOrthonormalBasis,
) where {N}
    dims = map(manifold_dimension, M.manifolds)
    dims_acc = accumulate(+, [1, dims...])
    parts = ntuple(N) do i
        get_vector(
            M.manifolds[i],
            submanifold_component(x, i),
            v[dims_acc[i]:dims_acc[i]+dims[i]-1],
            B.parts[i],
        )
    end
    return ProductRepr(parts)
end
function get_vector(
    M::ProductManifold{<:NTuple{N,Any}},
    x::ProductRepr,
    v,
    B::ArbitraryOrthonormalBasis,
) where {N}
    dims = map(manifold_dimension, M.manifolds)
    dims_acc = accumulate(+, [1, dims...])
    parts = ntuple(N) do i
        get_vector(
            M.manifolds[i],
            submanifold_component(x, i),
            v[dims_acc[i]:dims_acc[i]+dims[i]-1],
            B,
        )
    end
    return ProductRepr(parts)
end

function get_vectors(
    M::ProductManifold{<:NTuple{N,Manifold}},
    x::ProductRepr,
    B::PrecomputedProductOrthonormalBasis,
) where {N}
    xparts = submanifold_components(x)
    BVs = map(t -> get_vectors(t...), ziptuples(M.manifolds, xparts, B.parts))
    zero_tvs = map(t -> zero_tangent_vector(t...), ziptuples(M.manifolds, xparts))
    vs = typeof(ProductRepr(zero_tvs...))[]
    for i = 1:N, k = 1:length(BVs[i])
        push!(vs, ProductRepr(zero_tvs[1:i-1]..., BVs[i][k], zero_tvs[i+1:end]...))
    end
    return vs
end

function hat!(M::ProductManifold, v, x, vⁱ)
    dim = manifold_dimension(M)
    @assert length(vⁱ) == dim
    i = one(dim)
    ts = ziptuples(M.manifolds, submanifold_components(M, v), submanifold_components(M, x))
    for t ∈ ts
        dim = manifold_dimension(first(t))
        tvⁱ = @inbounds view(vⁱ, i:(i+dim-1))
        hat!(t..., tvⁱ)
        i += dim
    end
    return v
end

@doc doc"""
    injectivity_radius(M::ProductManifold[, x])

Compute the injectivity radius on the [`ProductManifold`](@ref), which is the
minimum of the factor manifolds.
"""
injectivity_radius(::ProductManifold, ::Any...)
function injectivity_radius(M::ProductManifold, x)
    return min(map(injectivity_radius, M.manifolds, submanifold_components(M, x))...)
end
injectivity_radius(M::ProductManifold) = min(map(injectivity_radius, M.manifolds)...)

@doc doc"""
    inner(M::ProductManifold, x, v, w)

compute the inner product of two tangent vectors `v`, `w` from the tangent space
at `x` on the [`ProductManifold`](@ref) `M`, which is just the sum of the
internal manifolds that build `M`.
"""
function inner(M::ProductManifold, x, v, w)
    subproducts = map(
        inner,
        M.manifolds,
        submanifold_components(M, x),
        submanifold_components(M, v),
        submanifold_components(M, w),
    )
    return sum(subproducts)
end

@doc doc"""
    inverse_retract(M::ProductManifold, x, y, m::InverseProductRetraction)

Compute the inverse retraction from `x` with respect to `y` on the [`ProductManifold`](@ref)
`M` using an [`InverseProductRetraction`](@ref), which by default encapsulates a inverse
retraction for each manifold of the product. Then this method is performed elementwise,
so the encapsulated inverse retraction methods have to be available per factor.
"""
inverse_retract(::ProductManifold, ::Any, ::Any, ::Any, ::InverseProductRetraction)

function inverse_retract!(M::ProductManifold, v, x, y, method::InverseProductRetraction)
    map(
        inverse_retract!,
        M.manifolds,
        submanifold_components(M, v),
        submanifold_components(M, x),
        submanifold_components(M, y),
        method.inverse_retractions,
    )
    return v
end

is_default_metric(::ProductManifold, ::ProductMetric) = Val(true)

function isapprox(M::ProductManifold, x, y; kwargs...)
    return all(
        t -> isapprox(t...; kwargs...),
        ziptuples(M.manifolds, submanifold_components(M, x), submanifold_components(M, y)),
    )
end
function isapprox(M::ProductManifold, x, v, w; kwargs...)
    return all(
        t -> isapprox(t...; kwargs...),
        ziptuples(
            M.manifolds,
            submanifold_components(M, x),
            submanifold_components(M, v),
            submanifold_components(M, w),
        ),
    )
end

@doc doc"""
    log(M::ProductManifold, x, y)

Compute the logarithmic map from `x` to `y` on the [`ProductManifold`](@ref) `M`,
which can be computed using the logarithmic maps of the manifolds elementwise.
"""
log(::ProductManifold, ::Any...)
function log(M::ProductManifold, x::ProductRepr, y::ProductRepr)
    return ProductRepr(map(
        log,
        M.manifolds,
        submanifold_components(M, x),
        submanifold_components(M, y),
    )...)
end

function log!(M::ProductManifold, v, x, y)
    map(
        log!,
        M.manifolds,
        submanifold_components(M, v),
        submanifold_components(M, x),
        submanifold_components(M, y),
    )
    return v
end

@doc doc"""
    manifold_dimension(M::ProductManifold)

Return the manifold dimension of the [`ProductManifold`](@ref), which is the sum of the
manifold dimensions the product is made of.
"""
manifold_dimension(M::ProductManifold) = mapreduce(manifold_dimension, +, M.manifolds)

@doc doc"""
    norm(M::PowerManifold, x, v)

Compute the norm of `v` from the tangent space of `x` on the [`ProductManifold`](@ref),
i.e. from the element wise norms the 2-norm is computed.
"""
function norm(M::ProductManifold, x, v)
    norms_squared = (
        map(
            norm,
            M.manifolds,
            submanifold_components(M, x),
            submanifold_components(M, v),
        ) .^ 2
    )
    return sqrt(sum(norms_squared))
end

function ProductFVectorDistribution(
    type::VectorBundleFibers{<:VectorSpaceType,<:ProductManifold},
    x::Union{AbstractArray,MPoint,ProductRepr},
    distributions::FVectorDistribution...,
)
    return ProductFVectorDistribution{typeof(type),typeof(distributions),typeof(x)}(
        type,
        x,
        distributions,
    )
end
function ProductFVectorDistribution(
    type::VectorBundleFibers{<:VectorSpaceType,<:ProductManifold},
    distributions::FVectorDistribution...,
)
    x = ProductRepr(map(d -> support(d).x, distributions))
    return ProductFVectorDistribution(type, x, distributions...)
end
function ProductFVectorDistribution(distributions::FVectorDistribution...)
    M = ProductManifold(map(d -> support(d).space.M, distributions)...)
    VS = support(distributions[1]).space.VS
    if !all(d -> support(d).space.VS == VS, distributions)
        error("Not all distributions have support in vector spaces of the same type, which is currently not supported")
    end
    # Probably worth considering sum spaces in the future?
    x = ProductRepr(map(d -> support(d).x, distributions)...)
    return ProductFVectorDistribution(VectorBundleFibers(VS, M), x, distributions...)
end

function ProductPointDistribution(M::ProductManifold, distributions::MPointDistribution...)
    return ProductPointDistribution{typeof(M),typeof(distributions)}(M, distributions)
end
function ProductPointDistribution(distributions::MPointDistribution...)
    M = ProductManifold(map(d -> support(d).manifold, distributions)...)
    return ProductPointDistribution(M, distributions...)
end

function rand(rng::AbstractRNG, d::ProductPointDistribution)
    return ProductRepr(map(d -> rand(rng, d), d.distributions)...)
end
function rand(rng::AbstractRNG, d::ProductFVectorDistribution)
    return ProductRepr(map(d -> rand(rng, d), d.distributions)...)
end

function _rand!(rng::AbstractRNG, d::ProductPointDistribution, x::AbstractArray{<:Number})
    return copyto!(x, rand(rng, d))
end
function _rand!(rng::AbstractRNG, d::ProductPointDistribution, x::ProductRepr)
    map(
        t -> _rand!(rng, t[1], t[2]),
        d.distributions,
        submanifold_components(d.manifold, x),
    )
    return x
end
function _rand!(rng::AbstractRNG, d::ProductFVectorDistribution, v::AbstractArray{<:Number})
    return copyto!(v, rand(rng, d))
end
function _rand!(rng::AbstractRNG, d::ProductFVectorDistribution, v::ProductRepr)
    map(t -> _rand!(rng, t[1], t[2]), d.distributions, submanifold_components(d.space.M, v))
    return v
end

@doc doc"""
    retract(M::ProductManifold, x, v, m::ProductRetraction)

Compute the retraction from `x` with tangent vector `v` on the [`ProductManifold`](@ref) `M`
using an [`ProductRetraction`](@ref), which by default encapsulates retractions of the
base manifolds. Then this method is performed elementwise, so the encapsulated retractions
method has to be one that is available on the manifolds.
"""
retract(::ProductManifold, ::Any...)

function retract!(M::ProductManifold, y, x, v, method::ProductRetraction)
    map(
        retract!,
        M.manifolds,
        submanifold_components(M, y),
        submanifold_components(M, x),
        submanifold_components(M, v),
        method.retractions,
    )
    return y
end

function representation_size(M::ProductManifold)
    return (mapreduce(m -> prod(representation_size(m)), +, M.manifolds),)
end

@doc doc"""
    sharp(M::ProductManifold, x, w::FVector{CotangentSpaceType})

Use the musical isomorphism to transform the cotangent vector `w` from the tangent space at
`x` on the [`ProductManifold`](@ref) `M` to a tangent vector.
This can be done elementwise, so vor every entry of `w` (and `x`) sparately
"""
sharp(::ProductManifold, ::Any...)

function sharp!(M::ProductManifold, v::TFVector, x, w::CoTFVector)
    vfs = map(u -> FVector(TangentSpace, u), submanifold_components(v))
    wfs = map(u -> FVector(CotangentSpace, u), submanifold_components(w))
    map(sharp!, M.manifolds, vfs, submanifold_components(M, x), wfs)
    return v
end

"""
    submanifold(M::ProductManifold, i::Integer)

Extract the `i`th factor of the product manifold `M`.
"""
submanifold(M::ProductManifold, i::Integer) = M.manifolds[i]

"""
    submanifold(M::ProductManifold, i::Val)
    submanifold(M::ProductManifold, i::AbstractVector)

Extract the factor of the product manifold `M` indicated by indices in `i`.
For example, for `i` equal to `Val((1, 3))` the product manifold constructed
from the first and the third factor is returned.

The version with `AbstractVector` is not type-stable, for better preformance use `Val`.
"""
submanifold(M::ProductManifold, i::Val) = ProductManifold(select_from_tuple(M.manifolds, i))
submanifold(M::ProductManifold, i::AbstractVector) = submanifold(M, Val(tuple(i...)))

support(d::ProductPointDistribution) = MPointSupport(d.manifold)
function support(tvd::ProductFVectorDistribution)
    return FVectorSupport(
        tvd.type,
        ProductRepr(map(d -> support(d).x, tvd.distributions)...),
    )
end

function vee!(M::ProductManifold, vⁱ, x, v)
    dim = manifold_dimension(M)
    @assert length(vⁱ) == dim
    i = one(dim)
    ts = ziptuples(M.manifolds, submanifold_components(M, x), submanifold_components(M, v))
    for t ∈ ts
        SM = first(t)
        dim = manifold_dimension(SM)
        tvⁱ = @inbounds view(vⁱ, i:(i+dim-1))
        vee!(SM, tvⁱ, Base.tail(t)...)
        i += dim
    end
    return vⁱ
end

function zero_tangent_vector!(M::ProductManifold, v, x)
    map(
        zero_tangent_vector!,
        M.manifolds,
        submanifold_components(M, v),
        submanifold_components(M, x),
    )
    return v
end


text:
# ---
# title: 53. Maximum Subarray
# id: problem53
# author: zhwang
# date: 2022-01-16
# difficulty: Easy
# categories: Array, Divide and Conquer, Dynamic Programming
# link: <https://leetcode.com/problems/maximum-subarray/description/>
# hidden: true
# ---
#
# Given an integer array `nums`, find the contiguous subarray (containing at
# least one number) which has the largest sum and return _its sum_.
#
# **Follow up:**  If you have figured out the `O(n)` solution, try coding
# another solution using the **divide and conquer** approach, which is more
# subtle.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [-2,1,-3,4,-1,2,1,-5,4]
#     Output: 6
#     Explanation: [4,-1,2,1] has the largest sum = 6.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1]
#     Output: 1
#
#
# **Example 3:**
#
#
#
#     Input: nums = [0]
#     Output: 0
#
#
# **Example 4:**
#
#
#
#     Input: nums = [-1]
#     Output: -1
#
#
# **Example 5:**
#
#
#
#     Input: nums = [-2147483647]
#     Output: -2147483647
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 2 * 104`
#   * `-231 <= nums[i] <= 231 - 1`
#
#
## @lc code=start
using LeetCode

function maximum_subarray(nums::Vector{Int})::Int
    pre, max_res = 0, nums[1]
    for num in nums
        pre = max(pre + num, num)
        max_res = max(pre, max_res)
    end
    return max_res
end

## @lc code=end

ext: jl


text:
function fft_advection(model, grid, fft_buoy_part)

    ## Grid [spatial & Fourier]
    model = init_grid_k(model)

    ## Initialisation of the spatial fields
    fft_w = sqg_large_uq(model, fft_buoy_part)
    w = real(ifft(fft_w))

    ## Hyperviscosity [order & coefficient]
    model.hv_order = 8
    model.hv_val = compute_hyper_coef(model, w)

    ## Choice of time step : CFL
    model.dt_adv = compute_dt(model, w)

    N_t = model.advection_duration / model.dt_adv

    # Printing some information
    println("Time step: $(model.dt_adv) seconds")
    println("Final time: $(N_t*model.dt_adv/3600/24) days ")

    ## Loop on time
    for t = 1:N_t
        ## solve "Poisson" to get w from buoy
        fft_w = sqg_large_uq(model, grid, fft_buoy_part)
        w = real(ifft(fft_w))
        ## Runge-Kutta 4 scheme
        fft_buoy_part = rk4_fft_advection(model, fft_buoy_part, w)

    end

    tt = floor(Int, Nt * model.dt_adv / (3600 * 24)) # Number of days
    println("$(t*model.dt_adv/(24*3600)) days of advection ")
    day = trunc(Int64, t * model.dt_adv / 24 / 3600)

    return fft_buoy_part

end


text:
struct WiderFactor{S<:Unsigned, T<:Unsigned} <: AbstractFactor{T}
    basefactor::AbstractFactor{S}
end

Base.length(factor::WiderFactor{S, T}) where {S<:Unsigned} where {T<:Unsigned} = length(factor.basefactor)

getlevels(factor::WiderFactor{S, T}) where {S<:Unsigned} where {T<:Unsigned} = getlevels(factor.basefactor)

getname(factor::WiderFactor{S, T}) where {S<:Unsigned} where {T<:Unsigned} = getname(factor.basefactor)

function slice(factor::WiderFactor{S, T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {S<:Unsigned} where {T<:Unsigned}
    f = x -> convert(T, x)
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    mapslice(f, slice(factor.basefactor, fromobs, toobs, slicelength), slicelength, T)
end

function Base.convert(::Type{WiderFactor{S, T}}, x::AbstractFactor{S}) where {S<:Unsigned} where {T<:Unsigned}
    WiderFactor{S, T}(x)
end

function Base.map(factor::WiderFactor, dataframe::AbstractDataFrame)
    map(factor.basefactor, dataframe)
end

function isordinal(factor::WiderFactor)
    isordinal(factor.basefactor)
end

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
using Plots, ThreadTools, Transducers

let s=:SBE5
    rng=(12:1:204).*u"V"
    coil=CoilFunctions.coil[s]
    p=CoilFunctions.power[s]
    f = v -> let
        Ω=CoilFunctions.Ω(s,v)
        turns=optimalcoil(coil,Ω,0.9).turns
        (v=v,turns=turns*24u"V"/v |> upreferred)
    end
    plot(;xlabel="Voltage (V)",ylabel="Relative Turns",title="Relative Turns $s")
    plot!(x->ustrip(x.v),y->y.turns,f.(rng);legend=:none)
end

let s=:SBE9,reqturns=520
    rng=(12:1:204).*u"V"
    coil=CoilFunctions.coil[s]
    p=CoilFunctions.power[s]
    f = v -> let
        Ω=CoilFunctions.Ω(s,v)
        minturns=reqturns*v/24u"V" |> ceil |> Int
        y = CoilFunctions.minimumcoillength(coil,Ω,minturns,0.9)
        (v=v,len=y.coil.len)
    end
    # g = minturns -> let
    #     maximum(x->x.len, f.(rng) )
    # end
    plot(;xlabel="Voltage (V)",ylabel="Relative Turns",title="Relative Turns $s")
    plot(x->ustrip(x.v),y->ustrip(y.len),f.(rng))
end

let s=:SBE7,reqturns=520
    vrng=(12:1:204).*u"V"
    coil=CoilFunctions.coil[s]
    p=CoilFunctions.power[s]
    f = (v,minturns) -> let
        Ω=v^2/p |> u"Ω"
        relative_turns=minturns*v/24u"V" |> ceil |> Int
        y = CoilFunctions.minimumcoillength(coil,Ω,relative_turns,0.9)
        (v=v,len=y.coil.len)
    end
    g = minturns -> let
        (minturns=minturns,len=maximum(x->x.len, f.(vrng,minturns)),)
    end
    # g(500)
    y=tmap(g,400:10:600)
    plot(;xlabel="Turns",ylabel="Minimum Coil Length (mm)",title="Relative Coil Length $s")
    plot!(x->x.minturns,y->ustrip(y.len),y)
end

plt=let
    vrng=(12:1:204).*u"V"
    plt=plot(;xlabel="Turns",ylabel="Minimum Coil Length (mm)",title="Relative Coil Length",
        legend=:topleft,ylims=(15,35))
    for s in (:SBE4,:SBE5,:SBE7,:SBE9)
        g = minturns -> let
            coil=CoilFunctions.coil[s]
            p=CoilFunctions.power[s]
            f = (v,minturns) -> let
                Ω=v^2/p |> u"Ω"
                relative_turns=minturns*v/24u"V" |> ceil |> Int
                y = CoilFunctions.minimumcoillength(coil,Ω,relative_turns,0.9)
                (v=v,len=y.coil.len)
            end
            (minturns=minturns,len=maximum(x->x.len, f.(vrng,minturns)),)
        end
        y=tmap(g,400:10:500)
        plot!(x->x.minturns,y->ustrip(y.len),y;label="$s")
    end
    plt
end

let s=:SBE4
    coil=CoilFunctions.coil[s]
    p=CoilFunctions.power[s]
    vrng = 12u"V":2u"V":208u"V"
    f = v -> let
        Ω = v^2/p |> u"Ω"
        merge((v=v,),CoilFunctions.optimalgeometry(coil,Ω,0.9))
    end
    y=tmap(f,vrng)
    plot(x->ustrip(x.v),y->ustrip(y.len-26u"mm"),y;label="Coil Length Variation")
    plot!(x->ustrip(x.v),y->ustrip(y.insul),y;label="Insulation")
end

let s=:SBE4
    coil=CoilGeometry(CoilFunctions.coil[s];len=24u"mm")
    p=CoilFunctions.power[s]
    vrng = 12u"V":2u"V":208u"V"
    f = v -> let
        Ω = v^2/p |> u"Ω"
        y=CoilFunctions.optimalcoil(coil,Ω,0.9)
        (v=v,insul=(y.coil.od-coil.od)/-2)
    end
    y=tmap(f,vrng)
    plot!(x->ustrip(x.v),y->ustrip(y.insul),y)
end

g=(s) -> let
    p=CoilFunctions.power[s]
    # vrng = 12u"V":1u"V":208u"V"
    vrng = 140u"V":0.1u"V":160u"V"
    rng = 26u"mm" .+ (-3u"mm":0.01u"mm":3u"mm")
    coil=CoilFunctions.coil[s]
    f = (v,len) -> let
        newcoil=CoilGeometry(coil;len=len)
        Ω = v^2/p |> u"Ω"
        y=CoilFunctions.optimalcoil(newcoil,Ω,0.9)
        (v=v,insul=(coil.od-y.coil.od)/2,len=len)
    end
    a=((v,len) for v in vrng, len in rng)
    g = x -> f(x[1],x[2])
    x=tmap(g,a)
    y=minimum(x->x.insul,x;dims=2)
    plot(ustrip.(vrng),ustrip.(y))
end

coil,p  = CoilGeometry(62.5mm,87.5mm,25mm,2mm), 30W
coil,p  = CoilGeometry(82.5mm,110mm,30mm,2mm), 44W
coil,p = CoilGeometry(100mm,140mm,30mm,2mm), 50W
coil,p  = CoilGeometry(140mm,190mm,30mm,2mm), 86W

map([12,24,48,104,180,207].*u"V") do x
    Ω=(x)^2/p |> u"Ω"
    I=(x)/Ω |> u"A"
    opt=CoilFunctions.optimalcoil(coil,Ω,0.9)
    # @info opt
    (V=x,Ω=Ω,gauge=Int(opt.gauge),turns=opt.turns,ampturns=opt.turns*I)
end

a=map((10:1:100).*W) do P
    P
    y=map([24,104]) do x
        Ω=(x*V)^2/P |> u"Ω"
        I=(x*V)/Ω |> u"A"
        opt=CoilFunctions.optimalcoil(coil,Ω)
        (V=x*V,Ω=Ω,gauge=opt.gauge,turns=opt.turns,ampturns=opt.turns*I)
    end
    ustrip(P),y[1].ampturns / y[2].ampturns
end

Plots.plot(x->first(x),x->last(x),a;ylabel="Ratio",xlabel="Power (Watts)")

b=let#
    h=30mm:0.1mm:60mm
    AWG=filter(x->first(x)==32,CoilFunctions.AWG_Chart)
    Ω=(104V)^2/30W |> u"Ω";
    # I=104V /
    a=map(h) do x
        a=CoilGeometry(coil;len=x)
        opt=CoilFunctions.optimalcoil(a,Ω;AWG_Chart=AWG)
    end
    Plots.plot(ustrip.(h),map(x->x.turns,a))

    h=20mm:0.1mm:60mm
    a=map(h) do x
        a=CoilGeometry(coil;len=x)
        opt=CoilFunctions.optimalcoil(a,Ω)
    end
    Plots.plot!(ustrip.(h),map(x->x.turns,a))
end

a = let
    id=100.0u"mm":1u"mm":105u"mm"
    od=140.0u"mm"
    v=24.0u"V"
    p=50.0u"W"
    Ω = v^2/p |> u"Ω"
    rng=(20:0.1:30)*u"mm"
    awg_chart = values(AWG_Chart)
    f = (id,len,gauge) -> begin
        pad=1.5u"mm"
        coil=CoilGeometry(id+2pad,od-2pad,len)
        CoilFunctions.estimatetruefill(coil,gauge,Ω,0.9)
    end
    y=((id,len,gauge) for len in rng for id in id for gauge in awg_chart) |> collect
    # map(x->f(x...),y)
    m = Map(x->f(x...)) |> Filter(x-> Ω/1.1 <= x.Ω <= Ω/0.9)
    tcollect(m,y)

end

plt=let
    s=:SBE5
    coil=CoilFunctions.coil[s]
    p=CoilFunctions.power[s]
    Ω=CoilFunctions.Ω(s,24u"V")
    rng=15u"mm":0.01u"mm":30u"mm"
    f = x -> begin
        let coil=CoilGeometry(coil;len=x)
            optimalcoil(coil,Ω,0.9).coil.od
        end
    end
    y=map(x->f(x),rng)
    plot(ustrip.(rng),ustrip.(y);legend=:none,
        xlabel="Length (mm)",ylabel="Build Diameter (mm)",title="SBE5 Build Diameter")
    plot!(ustrip.(rng[[1,end]]),ustrip.([coil.id,coil.id]))
    plot!(ustrip.(rng[[1,end]]),ustrip.([coil.od,coil.od]))
end

plt=let
    s=:SBE5
    coil=CoilFunctions.coil[s]
    p=CoilFunctions.power[s]
    Ω=CoilFunctions.Ω(s,24u"V")
    rng=15u"mm":0.01u"mm":30u"mm"
    f = (x,gauge) -> begin
        let coil=CoilGeometry(coil;len=x)
            estimatetruefill(coil,gauge,Ω,0.9)
        end
    end
    plt=plot(;legend=:topleft,xlabel="Length (mm)",ylabel="Build Diameter (mm)",title="SBE5 Build Diameter")
    for (key,gauge) in filter(x-> 23 <= first(x) <= 26, AWG_Chart)
        y=map(x->f(x,gauge).coil.od,rng)
        plot!(ustrip.(rng),ustrip.(y),label="$(Int(key)) awg")
    end
    # y=map(x->f(x),rng)
    # plot(ustrip.(rng),ustrip.(y);legend=:none,
    #     xlabel="Length (mm)",ylabel="Build Diameter (mm)",title="SBE5 Build Diameter")
    # plot!(ustrip.(rng[[1,end]]),ustrip.([coil.id,coil.id]))
    # plot!(ustrip.(rng[[1,end]]),ustrip.([coil.od,coil.od]))
    plt
end

max_stars_repo_path: scripts/Voltage Options.jl

max_issues_repo_path: scripts/Voltage Options.jl

max_forks_repo_path: scripts/Voltage Options.jl


text:
export PARSDMM_initialize

"""
subfunctions that initializes and checks, and distributes some quantities required for PARSDMM
"""
function PARSDMM_initialize(
                            x                       ::Vector{TF},
                            l                       ::Union{Vector{Vector{TF}},DistributedArrays.DArray{Array{TF,1},1,Array{Array{TF,1},1}},Array{Any,1}},
                            y                       ::Union{Vector{Vector{TF}},DistributedArrays.DArray{Array{TF,1},1,Array{Array{TF,1},1}},Array{Any,1}},
                            AtA                     ::Union{Vector{Array{TF, 2}}, Vector{SparseMatrixCSC{TF, TI}}, Vector{joAbstractLinearOperator{TF, TF}}, Vector{Union{Array{TF, 2}, SparseMatrixCSC{TF, TI}, joAbstractLinearOperator{TF, TF}}}},
                            TD_OP                   ::Union{Vector{Union{SparseMatrixCSC{TF,TI},JOLI.joAbstractLinearOperator{TF,TF}}},DistributedArrays.DArray{Union{JOLI.joAbstractLinearOperator{TF,TF}, SparseMatrixCSC{TF,TI}},1,Array{Union{JOLI.joAbstractLinearOperator{TF,TF}, SparseMatrixCSC{TF,TI}},1}} },
                            set_Prop,
                            P_sub                   ::Union{Vector{Any},DistributedArrays.DArray{Any,1,Array{Any,1}}},
                            comp_grid,
                            maxit                   ::Integer,
                            rho_ini                 ::Vector{Real},
                            gamma_ini               ::TF,
                            x_min_solver            ::String,
                            rho_update_frequency    ::Integer,
                            adjust_gamma            ::Bool,
                            adjust_rho              ::Bool,
                            adjust_feasibility_rho  ::Bool,
                            m                       ::Vector{TF},
                            parallel                ::Bool,
                            options,
                            zero_ini_guess          ::Bool,
                            feasibility_only=false  ::Bool
                            ) where {TF<:Real,TI<:Integer}

  ind_ref = maxit
  if options.Minkowski == false
    N = length(x)
  elseif options.Minkowski == true && zero_ini_guess==true
    N = length(x)*2
  elseif options.Minkowski == true && zero_ini_guess==false
    @assert length(x)==(2*length(m))
    N = length(x)
  end
  # const N = size(TD_OP[1],2) #this line will give really weird errors later on for some reason..

  # if typeof(AtA[1])==SparseMatrixCSC{TF,TI}
  #   push!(AtA,speye(TF,N));
  # elseif typeof(AtA[1])==Array{TF,2}
  #   push!(AtA,ones(TF,N,1))
  # end
  # push!(TD_OP,speye(TF,N));
  #
  # push!(set_Prop.AtA_offsets,[0])
  # push!(set_Prop.banded,true)
  # push!(set_Prop.AtA_diag,true)
  # push!(set_Prop.dense,false)

  p  = length(TD_OP); #number of terms in in the sum of functions of the projeciton problem
  pp = p-1;
  if feasibility_only==true; pp=p; end;

  rho     = Vector{TF}(undef,p)
  if length(rho_ini)==1
    fill!(rho,rho_ini[1])
  else
    copy!(rho,rho_ini)
  end
  prox = copy(P_sub)
  if feasibility_only==false
    #define prox for all terms in the sum (projectors onto sets)
    #add prox for the data fidelity term 0.5||m-x||_2^2
    m_orig = deepcopy(m) ::Vector{TF}
    prox_data = input -> prox_l2s!(input,rho[end],m_orig)
    push!(prox,prox_data);
  end
  if parallel==true
    prox  = distribute(prox)
    P_sub = distribute(P_sub)
    rho   = distribute(rho)
  end

  if parallel==true && nworkers()<length(TD_OP)
      error("parallel PARSDMM requires at least one JULIA worker per constraint set + 1 (",length(TD_OP),")")
  end

  # detect feasibility of model that needs to be projected
  stop = false
  feasibility_initial = zeros(TF,length(P_sub))
  if options.Minkowski == true
    m = [m ; zeros(TF,length(m)) ]
  end
  if parallel
    feasibility_initial = distribute(feasibility_initial)
    [@sync @spawnat pid m for pid in P_sub.pids]
    [@sync @spawnat pid compute_relative_feasibility(m,feasibility_initial[:L],TD_OP[:L],P_sub[:L]) for pid in P_sub.pids]
    feasibility_initial = convert(Vector{TF},feasibility_initial)

    #using pmap
    #feasibility_initial = pmap((feasibility_initial,TD_OP,P_sub) -> compute_relative_feasibility(m,feasibility_initial,TD_OP,P_sub) , feasibility_initial,TD_OP,P_sub; distributed=true, batch_size=1, on_error=nothing, retry_delays=[], retry_check=nothing)
  else
    for ii = 1:length(P_sub)
      feasibility_initial[ii] = norm(P_sub[ii](TD_OP[ii]*m) .- TD_OP[ii]*m) ./ (norm(TD_OP[ii]*m)+(100*eps(TF)));
    end
  end
  if maximum(feasibility_initial)<options.feas_tol #accept input as feasible and return
      println("input to PARSDMM is feasible, returning")
      stop = true
  end

  # if one of the sets is non-convex, use different lambda and rho update frequency, don't update gamma and set a different fixed gamma
  for ii=1:pp
      if set_Prop.ncvx[ii] == true
          println("non-convex set(s) involved, using special settings")
          rho_update_frequency  = 3;
          adjust_gamma          = false
          gamma_ini             = TF(0.75)
      end
  end

  #allocate arrays of vectors
  Ax_out = zeros(TF,N)

  #if y and l are empty, allocate them and fill with zeros
  if isempty(l)==true
    l=Vector{Vector{TF}}(undef,p)
    for i=1:length(TD_OP); l[i]=zeros(TF,size(TD_OP[i],1)); end
  end
  if isempty(y)==true
    y=Vector{Vector{TF}}(undef,p)
    for i=1:length(TD_OP); y[i]=zeros(TF,size(TD_OP[i],1)); end
  end

  gamma   = Vector{TF}(undef,p);

  #y       = Vector{Vector{TF}}(p);
  y_0     = Vector{Vector{TF}}(undef,p);
  y_old   = Vector{Vector{TF}}(undef,p);

  #l       = Vector{Vector{TF}}(p);
  l_0     = Vector{Vector{TF}}(undef,p);
  l_old   = Vector{Vector{TF}}(undef,p);
  l_hat_0 = Vector{Vector{TF}}(undef,p);
  l_hat   = Vector{Vector{TF}}(undef,p)

  x_0     = Vector{TF}(undef,N);
  x_old   = Vector{TF}(undef,N);
  x_hat   = Vector{Vector{TF}}(undef,p);

  r_dual  = Vector{Vector{TF}}(undef,p);
  r_pri   = Vector{Vector{TF}}(undef,p);

  rhs     = Vector{TF}(undef,N);

  s       = Vector{Vector{TF}}(undef,p);
  s_0     = Vector{Vector{TF}}(undef,p);

  d_l_hat = Vector{Vector{TF}}(undef,p)
  d_H_hat = Vector{Vector{TF}}(undef,p)
  d_l     = Vector{Vector{TF}}(undef,p)
  d_G_hat = Vector{Vector{TF}}(undef,p)

  for ii=1:p #initialize all rho's, gamma's, y's and l's
      gamma[ii]   = gamma_ini;

      if parallel == false
        ly          = size(TD_OP[ii],1)
      else
        ly = 1
      end
      y_0[ii]     = zeros(TF,ly);#copy(y[ii])#zeros(ly);#0.*y[ii];
      y_old[ii]   = zeros(TF,ly);#copy(y[ii])#zeros(ly);#0.*y[ii];

      l_old[ii]   = zeros(TF,ly);#0.*y[ii];
      #l[ii]       = zeros(ly);#0.*y[ii];
      l_0[ii]     = zeros(TF,ly);#0.*y[ii];
      l_hat_0[ii] = zeros(TF,ly);#0.*y[ii];
      l_hat[ii]   = zeros(TF,ly)

      x_hat[ii]   = zeros(TF,ly)
      s_0[ii]     = zeros(TF,ly)
      s[ii]       = zeros(TF,ly)
      r_pri[ii]   = zeros(TF,ly)

      d_l_hat[ii] = zeros(TF,ly)
      d_H_hat[ii] = zeros(TF,ly)
      d_l[ii]     = zeros(TF,ly)
      d_G_hat[ii] = zeros(TF,ly)
  end

  #assemble total transform domain operator in normal equation form (A^t A)
  #note that matrices have already been converted to diagonal storage format if possible for all operators in AtA

  n_sparse_arrays  = 0
  n_JOLI_operators = 0
  n_CDS_matrices   = 0
  for i=1:p; if typeof(AtA[i])==SparseMatrixCSC{TF,TI}; n_sparse_arrays += 1; end; end
  for i=1:p; if typeof(AtA[i]) <: joAbstractLinearOperator; n_JOLI_operators += 1; end; end
  for i=1:p; if typeof(AtA[i])==Array{TF,2}; n_CDS_matrices += 1; end; end

  if n_sparse_arrays == 2 #if all arrays in AtA are sparse matrices
    Q_offsets = []
    Q = rho[1]*AtA[1];
    for i=2:p
        Q .= Q .+ rho[i]*AtA[i]
    end
  elseif n_JOLI_operators >= 1 #at least one operator is a JOLI operator -> wrap all others as a JOLI matrix as well
    Q_offsets = []
    if typeof(AtA[1]) <: joAbstractLinearOperator
      Q = rho[1]*AtA[1]
    else
      Q = rho[1]*joMatrix(AtA[1])
    end
    for i=2:p
      if typeof(AtA[i]) <: joAbstractLinearOperator
        Q = Q + rho[i]*AtA[i]
      else
        Q = Q + rho[i]*joMatrix(AtA[i])
      end
    end
  else #n_CDS_matrices == p should hold #all matrics in AtA are in compressed diagonal storage format
    all_offsets=zeros(TI,999,99)
    for i=1:length(AtA) #find all unique offset
      all_offsets[1:length(set_Prop.AtA_offsets[i]),i]=set_Prop.AtA_offsets[i]
    end
    Q_offsets=convert(Vector{TI},unique(all_offsets))
    Q=zeros(TF,size(AtA[1],1),length(Q_offsets))
    for i=1:length(AtA)
      for j=1:length(set_Prop.AtA_offsets[i])
        #Q_current_col = findin(Q_offsets,set_Prop.AtA_offsets[i][j])
        Q_current_col = findall((in)(set_Prop.AtA_offsets[i][j]),Q_offsets)
        Q[:,Q_current_col] .= Q[:,Q_current_col] .+ rho[i] .* AtA[i][:,j]
      end
    end
  end

  #initialize logs
  log_PARSDMM = log_type_PARSDMM(zeros(maxit,pp),zeros(maxit,p),zeros(maxit,p),zeros(maxit),zeros(maxit),zeros(maxit),
  zeros(maxit),zeros(maxit,p),zeros(maxit,p),zeros(maxit),zeros(maxit),TimerOutput());

  log_PARSDMM.set_feasibility[1,:] = feasibility_initial

  if parallel == true

    gamma = distribute(gamma)

    #y     = distribute(y)
    y_0   = distribute(y_0)
    y_old = distribute(y_old)

    #l       = distribute(l)
    l_0     = distribute(l_0)
    l_old   = distribute(l_old)
    l_hat_0 = distribute(l_hat_0)
    l_hat   = distribute(l_hat)

    x_hat   = distribute(x_hat)
    s_0     = distribute(s_0)
    s       = distribute(s)
    r_pri   = distribute(r_pri)

    d_l_hat = distribute(d_l_hat)
    d_H_hat = distribute(d_H_hat)
    d_l     = distribute(d_l)
    d_G_hat = distribute(d_G_hat)
  end

  #fill distributed vectors with zeros (from 1 entry to N entries, because this is faster than first fill and then dist. Should be able to do in one go ideally)
  if parallel == true
    [@spawnat pid y_0[:L][1]     = zeros(TF,size(TD_OP[:L][1],1)) for pid in y_0.pids]
    [@spawnat pid y_old[:L][1]   = zeros(TF,size(TD_OP[:L][1],1)) for pid in y_0.pids]
    [@spawnat pid l_0[:L][1]     = zeros(TF,size(TD_OP[:L][1],1)) for pid in y_0.pids]
    [@spawnat pid l_old[:L][1]   = zeros(TF,size(TD_OP[:L][1],1)) for pid in y_0.pids]
    [@spawnat pid l_hat_0[:L][1] = zeros(TF,size(TD_OP[:L][1],1)) for pid in y_0.pids]
    [@spawnat pid l_hat[:L][1]   = zeros(TF,size(TD_OP[:L][1],1)) for pid in y_0.pids]
    [@spawnat pid x_hat[:L][1]   = zeros(TF,size(TD_OP[:L][1],1)) for pid in y_0.pids]
    [@spawnat pid s_0[:L][1]     = zeros(TF,size(TD_OP[:L][1],1)) for pid in y_0.pids]
    [@spawnat pid s[:L][1]       = zeros(TF,size(TD_OP[:L][1],1)) for pid in y_0.pids]
    [@spawnat pid r_pri[:L][1]   = zeros(TF,size(TD_OP[:L][1],1)) for pid in y_0.pids]
    [@spawnat pid d_l_hat[:L][1] = zeros(TF,size(TD_OP[:L][1],1)) for pid in y_0.pids]
    [@spawnat pid d_H_hat[:L][1] = zeros(TF,size(TD_OP[:L][1],1)) for pid in y_0.pids]
    [@spawnat pid d_l[:L][1]     = zeros(TF,size(TD_OP[:L][1],1)) for pid in y_0.pids]
    [@spawnat pid d_G_hat[:L][1] = zeros(TF,size(TD_OP[:L][1],1)) for pid in y_0.pids]
  end

  #distribute y and l if they are not already distributed
  if (parallel==true) && (typeof(l)<:DistributedArrays.DArray) == false
    l = distribute(l)
  end
  if parallel==true && (typeof(y)<:DistributedArrays.DArray) == false
    y = distribute(y)
  end
  if parallel==true
    [ @spawnat pid adjust_gamma for pid in y.pids ]
    [ @spawnat pid adjust_rho for pid in y.pids ]
    [ @spawnat pid adjust_feasibility_rho for pid in y.pids ]
  end
  if parallel == true
    [@spawnat pid comp_grid for pid in  y.pids]
  end

  if parallel == true
    set_feas = zeros(TF,length(TD_OP)) #initialize a vector only required for parallel stuff
    set_feas = distribute(set_feas)
  else
    set_feas = []
  end

  if zero_ini_guess == true
    if parallel==false
      for i=1:length(l); fill!(l[i],TF(0.0)); end
      for i=1:length(y); fill!(y[i],TF(0.0)); end
    else
      [@spawnat pid fill!(l[:L][1],TF(0.0)) for pid in l.pids]
      [@spawnat pid fill!(y[:L][1],TF(0.0)) for pid in y.pids]
    end
    fill!(x,TF(0.0))
  end

return ind_ref,N,TD_OP,AtA,p,rho_update_frequency,adjust_gamma,adjust_rho,adjust_feasibility_rho,gamma_ini,rho,gamma,y,y_0,y_old,
       l,l_0,l_old,l_hat_0,x_0,x_old,r_dual,rhs,s,s_0,Q,prox,log_PARSDMM,l_hat,x_hat,r_pri,
       d_l_hat,d_H_hat,d_l,d_G_hat,P_sub,Q_offsets,stop,feasibility_initial,set_feas,Ax_out
end #end parsdmm_initialize

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
using Dates
using Distributed
using JSON
using Test

using Discord
using Discord:
    EMPTY,
    EVENTS_FIRED,
    Bucket,
    Handler,
    Response,
    Limiter,Snowflake,
    allhandlers,
    datetime,
    get_channel_message,
    handlers,
    hasdefault,
    increment,
    insert_or_update!,
    isexpired,
    islimited,
    parse_endpoint,
    process_id,
    readjson,
    should_put,
    snowflake,
    snowflake2datetime,
    worker_id,
    writejson,
    @boilerplate,
    @dict,
    @lower,
    @merge

# A test case which covers most possible field types.
@eval Discord struct Foo
    a::String
    b::DateTime
    c::Snowflake
    d::Vector{String}
    e::Union{Int, Nothing}
    f::Union{Int, Missing}
    g::Union{Vector{String}, Nothing, Missing}
    h::Union{Foo, Missing}
    abcdefghij::Union{Int, Missing}  # 10 characters.
end
@eval Discord @boilerplate Foo :dict :docs :lower :merge

# A simple struct with merge.
@eval Discord struct Bar
    id::Int
    x::Int
    y::Union{Int, Missing}
end
@eval Discord @boilerplate Bar :merge

using Discord: Foo, Bar

# A module with event handlers.
module Handlers

export a, b

using Discord

a(::Client, ::AbstractEvent) = nothing
a(::Client, ::TypingStart) = nothing
b(::Client, ::WebhooksUpdate) = nothing
c(::Client, ::WebhooksUpdate) = nothing

end

@testset "Discord.jl" begin
    c = Client("token")

    @testset "Client" begin
        @testset "Token" begin
            # Tokens should always be prepended by "Bot ".
            c = Client("token")
            @test c.token == "Bot token"
            # Unless it's already there.
            c = Client("Bot token")
            @test c.token == "Bot token"
        end

        @testset "Utils" begin
            # We haven't connected to the gateway yet.
            @test me(c) === nothing

            disable_cache!(c)
            @test !c.use_cache
            enable_cache!(c) do
                @test c.use_cache
            end
            @test !c.use_cache
            enable_cache!(c)
            @test c.use_cache
            disable_cache!(c) do
                @test !c.use_cache
            end
            @test c.use_cache
        end
    end

    @testset "Parsing" begin
        @testset "tryparse" begin
            @test tryparse(c, Int, 123) == (123, nothing)
            @test tryparse(c, Vector{UInt16}, Int[1, 2, 3]) == (UInt16[1, 2, 3], nothing)
            val, e = @test_logs (:error, r"MethodError") tryparse(c, Int, Dict())
            @test val === nothing
            @test e isa MethodError
            @test first(c.state.errors) == Dict()
        end

        @testset "Snowflake" begin
            # Snowflakes usually come in as strings.
            s = snowflake(string(typemax(UInt64)))
            @test isa(s, Snowflake)
            # Snowflake === UInt64.
            @test s == typemax(UInt64)

            # https://discordapp.com/developers/docs/reference#snowflakes-snowflake-id-format-structure-left-to-right
            s = Snowflake(0x06ecefa78e42000c)
            @test snowflake2datetime(s) == DateTime(2018, 10, 9, 1, 55, 31, 1)
            @test worker_id(s) == 0x01
            @test process_id(s) == 0x00
            @test increment(s) == 0x0c
        end

        @testset "DateTime" begin
            # Discord sends dates in some weird, inconsistent ways.
            d = datetime("2018-10-08T05:20:22.782643+00:00")
            @test d == DateTime(2018, 10, 8, 5, 20, 22, 782)
            d = datetime("2018-10-08T05:20:22.782Z")
            @test d == DateTime(2018, 10, 8, 5, 20, 22, 782)

            # But sometimes they also send nice millisecond Unix timestamps.
            d = datetime(1541288588543)
            @test d == DateTime(2018, 11, 3, 23, 43, 08, 543)
        end
    end

    @testset "[read|write]json" begin
        io = IOBuffer()
        val, e = readjson(io)
        @test val === nothing
        @test e == EMPTY

        io = IOBuffer("{bad]")
        val, e = readjson(io)
        @test val === nothing
        @test e !== nothing

        io = IOBuffer("[1,2,3]")
        val, e = readjson(io)
        @test val == [1, 2, 3]
        @test e === nothing

        io = IOBuffer()
        @test writejson(io, @__MODULE__) !== nothing

        io = IOBuffer()
        @test writejson(io, [1, 2, 3]) === nothing
        @test String(take!(io)) == "[1,2,3]"

        io = IOBuffer()
        d = Dict("foo" => 0, "bar" => [1, 2, 3])
        @test writejson(io, d) === nothing
        @test JSON.parse(String(take!(io))) == d
    end

    @testset "Helpers" begin
        ch = DiscordChannel(Dict("id" => "255", "type" => 0, "guild_id" => "1"))
        r = Role(0xff, "", 0, true, 0, 0, false, false)
        u = User(Dict{String, Any}("id" => "255", "username" => "foo"))

        @testset "mention" begin
            @test mention(ch) == "<#255>"
            @test mention(r) == "<@&255>"
            @test mention(u) == "<@255>"
            m = Member(u, "foo", [], now(), true, true)
            @test mention(m) == "<@!255>"
            m = Member(u, nothing, [], now(), true, true)
            @test mention(m) == mention(u)
            m = Member(u, missing, [], now(), true, true)
            @test mention(m) == mention(u)
        end

        @testset "plaintext" begin
            msg = Message(Dict(
                "id" => "1",
                "channel_id" => "1",
                "content" => "<@255> <@!255>",
                "mentions" => [JSON.lower(u)],
            ))
            @test plaintext(msg) == "@foo @foo"
        end
    end

    @testset "REST API" begin
        @testset "Direct endpoint wrapper" begin
            # Direct endpoint wrappers should return a Future.
            f = get_channel_message(c, 123, 456)
            @test f isa Future
            # Since we don't have a valid token, we shouldn't get anything.
            @test fetchval(f) === nothing
            r = fetch(f)
            # But the type should still be sound.
            @test r isa Response{Message}
            @test !r.ok
            @test r.val === nothing
            @test r.http_response !== nothing
            @test r.exception === nothing
        end

        @testset "CRUD" begin
            # This API should behave just like the direct endpoint wrappers.
            f = create(c, Guild; name="foo")
            @test f isa Future
            r = fetch(f)
            @test r isa Response{Guild}
            @test r.val === nothing
            @test !r.ok
            @test r.http_response !== nothing
            @test r.exception === nothing
        end

        @testset "Simultaneous requests" begin
            # We should be able to make a bunch of requests without deadlocking.
            fs = map(i -> retrieve(c, Guild, i), 1:10)
            @test all(f -> f isa Future, fs)
            rs = fetch.(fs)
            @test all(r -> r isa Response{Guild}, rs)
        end

        @testset "Selective caching" begin
            f(c, e) = nothing

            # The regexes should be strict.
            @test all(
                p -> startswith(p, '^') && endswith(p, '$'),
                map(t -> t[1].pattern, vcat(values(EVENTS_FIRED)...)),
            )

            # We don't need to worry about this fast path.
            @eval Discord Base.isopen(::Client) = true

            # GETs should always be cached, because they never fire events.
            @test should_put(c, :GET, "foo")

            # If we have the default handler in place, we shouldn't put.
            add_handler!(c, MessageCreate, f; tag=DEFAULT_HANDLER_TAG)
            @test !should_put(c, :POST, "/channels/1/messages")

            # But if we delete it, then we should.
            delete_handler!(c, MessageCreate, DEFAULT_HANDLER_TAG)
            @test should_put(c, :POST, "/channels/1/messages")

            # For endpoints which fire multiple events, put unless they all have a default.
            add_handler!(c, GuildMemberRemove, f; tag=DEFAULT_HANDLER_TAG)
            add_handler!(c, PresenceUpdate, f; tag=DEFAULT_HANDLER_TAG)
            @test !should_put(c, :DELETE, "/users/@me/guilds/1")
            delete_handler!(c, PresenceUpdate, DEFAULT_HANDLER_TAG)
            @test should_put(c, :DELETE, "/users/@me/guilds/1")
        end
    end

    @testset "Rate limiting" begin
        @testset "Buckets" begin
            l = Limiter()

            # We start with no buckets.
            @test isempty(l.buckets)
            # But when we need one, it gets created for us.
            b = Bucket(l, :GET, "/foo")
            @test collect(keys(l.buckets)) == ["/foo"]
            # An unused bucket isn't limited.
            @test !islimited(l, b)

            # We can treat buckets like a lock.
            lock(b)
            @test b.sem.curr_cnt == 1
            unlock(b)
            @test b.sem.curr_cnt == 0

            # As long as the bucket isn't empty, we aren't limited.
            b.remaining = 1
            # Even if the reset hasn't yet passed.
            b.reset = now(UTC) + Second(1)
            @test !islimited(l, b)

            # But if the bucket is empty, then we're limited.
            b.remaining = 0
            @test islimited(l, b)

            # Once we wait, the bucket is reset.
            wait(l, b)
            b = Bucket(l, :GET, "/foo")
            @test !islimited(l, b)
            @test b.remaining === nothing
            @test b.reset === nothing
        end

        @testset "parse_endpoint" begin
            # The variable parameter doesn't matter.
            @test parse_endpoint("/users/1", :GET) == "/users"

            # Unless it's one of these three.
            @test parse_endpoint("/channels/1", :GET) == "/channels/1"
            @test parse_endpoint("/guilds/1", :GET) == "/guilds/1"
            @test parse_endpoint("/webhooks/1", :GET) == "/webhooks/1"

            # Without a numeric parameter at the end, we get the whole endpoint.
            @test parse_endpoint("/users/@me/channels", :GET) == "/users/@me/channels"

            # Special case 1: Deleting messages.
            @test ==(
                parse_endpoint("/channels/1/messages/1", :DELETE),
                "/channels/1/messages DELETE",
            )
            @test parse_endpoint("/channels/1/messages/1", :GET) == "/channels/1/messages"

            # Special case 2: Invites.
            @test parse_endpoint("/invites/abcdef", :GET) == "/invites"
        end
    end

    @testset "Boilerplate" begin
        local f
        d = Dict(
            "a" => "foo",
            "b" => "2018-10-08T05:20:22.782Z",
            "c" => "1234567890",
            "d" => ["a", "b", "c"],
            "e" => 1,
            "f" => 2,
            "g" => ["a", "b", "c"],
        )
        d["h"] = copy(d)

        @testset "@dict" begin
            f = Foo(d)
            @test f.a == "foo"
            @test f.b == DateTime(2018, 10, 8, 5, 20, 22, 782)
            @test f.c == 1234567890
            @test f.d == ["a", "b", "c"]
            @test f.e == 1
            @test f.f == 2
            @test f.g == ["a", "b", "c"]
            @test f.h isa Foo && f.h.a == "foo"

            # Set e to nothing, f and g to missing.
            d["e"] = nothing
            delete!(d, "f")
            delete!(d, "g")
            f = Foo(d)
            @test f.e === nothing
            @test ismissing(f.f)
            @test ismissing(f.g)

            # Union{T, Nothing, Missing} works too.
            d["g"] = nothing
            f = Foo(d)
            @test f.g === nothing
        end

        @testset "@docs" begin
            docs = string(@doc Foo)
            # Variable names get padded to the longest one.
            @test occursin("$(rpad("a", 10)) :: String", docs)
            # Array{T,1} is replaced with Vector{T}.
            @test occursin("Vector{String}", docs)
            # Union{Missing, Nothing, T} is replaced with Union{T, Missing, Nothing}.
            @test occursin("Union{Vector{String}, Missing, Nothing}", docs)
            # Union{Missing, T} is replaced with Union{T, Missing}.
            @test occursin("Union{Int, Missing}", docs)
            # Union{Nothing, T} is replaced with Union{T, Nothing}.
            @test occursin("Union{Int, Nothing}", docs)
        end

        @testset "@lower" begin
            d["b"] = d["h"]["b"] = round(Int, datetime2unix(f.b))
            d["c"] = d["h"]["c"] = snowflake(f.c)
            d = JSON.lower(f)
            # The result is always a Dict{String, Any}.
            @test d isa Dict{String, Any}
            @test d == Dict(
                "a" => "foo",
                "b" => d["b"],
                "c" => d["c"],
                "d" => ["a", "b", "c"],
                "e" => nothing,
                "g" => nothing,
                "h" => d["h"],
            )
        end

        @testset "@merge" begin
            # Anything except missing values should be taken from f2.
            f2 = Foo("bar", f.b, f.c, f.d, f.e, f.f, f.g, missing, missing)
            f3 = merge(f, f2)
            @test f3.a == f2.a
            @test f3.h == f.h
        end
    end

    @testset "Handlers" begin
        c = Client("token")
        f(c, e) = nothing
        g(c, e) = nothing
        badh(c, e::String) = nothing
        badc(c, e::AbstractEvent) = nothing

        @testset "Adding/deleting regular handlers" begin
            # Deleting handlers without a tag clears all handlers for that type.
            delete_handler!(c, MessageCreate)
            @test !haskey(c.handlers, MessageCreate)

            # Adding handlers without a tag means we can have duplicates.
            add_handler!(c, MessageCreate, f)
            add_handler!(c, MessageCreate, f)
            @test length(get(c.handlers, MessageCreate, Dict())) == 2
            delete_handler!(c, MessageCreate)

            # Using tags prevents duplicates.
            add_handler!(c, MessageCreate, f; tag=:f)
            add_handler!(c, MessageCreate, f; tag=:f)
            @test length(get(c.handlers, MessageCreate, Dict())) == 1
            @test haskey(c.handlers[MessageCreate], :f)

            # With tags, we can delete specific handlers.
            add_handler!(c, MessageCreate, g; tag=:g)
            @test length(get(c.handlers, MessageCreate, Dict())) == 2
            delete_handler!(c, MessageCreate, :g)
            @test length(get(c.handlers, MessageCreate, Dict())) == 1
            @test first(values(c.handlers[MessageCreate])).f == f

            # We can also add handlers from a module.
            add_handler!(c, Handlers)
            @test length(get(c.handlers, AbstractEvent, Dict())) == 1
            @test length(get(c.handlers, TypingStart, Dict())) == 1
            # Only exported functions are considered.
            @test length(get(c.handlers, WebhooksUpdate, Dict())) == 1
            @test first(values(c.handlers[WebhooksUpdate])).f == Handlers.b

            # Adding a module handler with a tag and/or expiry propogates to all handlers.
            empty!(c.handlers) # XXX
            add_handler!(c, Handlers; tag=:h, expiry=Millisecond(50))
            @test all(hs -> all(p -> p.first === :h, hs), values(c.handlers))
            sleep(Millisecond(50))
            @test all(hs -> all(p -> isexpired(p.second), hs), values(c.handlers))

            # We can't add a handler without a valid method.
            @test_throws ArgumentError add_handler!(c, MessageCreate, badh)

            # We can't add a handler that's already expired.
            @test_throws ArgumentError add_handler!(c, Ready, f; expiry=0)
            @test_throws ArgumentError add_handler!(c, Ready, f; expiry=Day(-1))
        end

        @testset "Commands" begin
            delete_handler!(c, MessageCreate)

            # Adding commands adds to the MessageCreate handlers.
            add_command!(c, "!test", f)
            @test length(get(c.handlers, MessageCreate, Dict())) == 1
            # But the handler function is modified.
            @test first(values(c.handlers[MessageCreate])).f != f

            # We can't add a command without a valid method.
            @test_throws ArgumentError add_command!(c, "!test", badc)

            # We can't add a command that's already expired.
            @test_throws ArgumentError add_command!(c, "!test", f; expiry=0)
            @test_throws ArgumentError add_handler!(c, Ready, f; expiry=Day(-1))
        end

        @testset "Handler expiry" begin
            empty!(c.handlers[MessageCreate])

            # By default, handlers don't expire.
            add_handler!(c, MessageCreate, f)
            @test !isexpired(first(values(c.handlers[MessageCreate])))

            add_handler!(c, MessageCreate, f; expiry=Millisecond(100))
            @test count(isexpired, values(c.handlers[MessageCreate])) == 0
            sleep(Millisecond(100))
            @test count(isexpired, values(c.handlers[MessageCreate])) == 1

            # Counting handlers expire when they reach <= 0.
            @test !isexpired(Handler(f, 1))
            @test isexpired(Handler(f, 0))
            @test isexpired(Handler(f, -1))

            # Timed handlers expire when their expiry time is reached.
            @test !isexpired(Handler(f, now() + Day(1)))
            @test isexpired(Handler(f, now() - Day(1)))
        end

        @testset "Handler collection" begin
            delete_handler!(c, Ready)
            add_handler!(c, Ready, f)
            add_handler!(c, Ready, f)
            add_handler!(c, AbstractEvent, f)
            add_handler!(c, FallbackEvent, f)
            @test length(handlers(c, Ready)) == 2

            # No handlers means no handlers.
            empty!(c.handlers)
            @test isempty(allhandlers(c, MessageCreate))
            @test isempty(allhandlers(c, AbstractEvent))
            @test isempty(allhandlers(c, FallbackEvent))

            # Both the specific and catch-all handler should match.
            add_handler!(c, Ready, f)
            add_handler!(c, AbstractEvent, f)
            @test allhandlers(c, Ready) == [
                collect(c.handlers[AbstractEvent]);
                collect(c.handlers[Ready]);
            ]

            # The fallback handler should only match if there are non non-default handlers.
            add_handler!(c, FallbackEvent, f)
            @test allhandlers(c, Ready) == [
                collect(c.handlers[AbstractEvent]);
                collect(c.handlers[Ready]);
            ]
            delete_handler!(c, Ready)
            @test allhandlers(c, Ready) == collect(c.handlers[AbstractEvent])
            add_handler!(c, Ready, f)
            delete_handler!(c, AbstractEvent)
            @test allhandlers(c, Ready) == collect(c.handlers[Ready])
            delete_handler!(c, Ready)
            add_handler!(c, Ready, f; tag=DEFAULT_HANDLER_TAG)
            @test allhandlers(c, Ready) == [
                collect(c.handlers[Ready]);
                collect(c.handlers[FallbackEvent]);
            ]
            delete_handler!(c, Ready)
            @test allhandlers(c, Ready) == collect(c.handlers[FallbackEvent])

            # Expired handlers should be cleaned up.
            first(values(c.handlers[FallbackEvent])).expiry = 0
            allhandlers(c, Ready)  # This triggers the cleanup.
            @test isempty(c.handlers[FallbackEvent])
        end

        @testset "Default handler lookup" begin
            c = Client("token")
            @test hasdefault(c, MessageCreate)
            delete_handler!(c, MessageCreate, DEFAULT_HANDLER_TAG)
            @test !hasdefault(c, MessageCreate)
            add_handler!(c, MessageCreate, f; tag=DEFAULT_HANDLER_TAG)
            @test hasdefault(c, MessageCreate)
        end
    end

    @testset "State" begin
        @testset "insert_or_update!" begin
            d = Dict()
            v = Bar[]

            # Inserting a new entry works fine.
            b = Bar(1, 2, 3)
            insert_or_update!(d, b.id, b)
            @test d[b.id] == b
            # But when we insert a value with the same ID, it's merged.
            b = Bar(1, 3, 3)
            insert_or_update!(d, b.id, b)
            @test d[b.id] == b
            b = Bar(1, 10, missing)
            insert_or_update!(d, b.id, b)
            @test d[b.id] == Bar(b.id, b.x, 3)

            # We can also insert/update without specifying the key, :id is inferred.
            b = Bar(2, 2, 1)
            insert_or_update!(d, b)
            @test d[b.id] == b
            b = Bar(2, 1, missing)
            insert_or_update!(d, b)
            @test d[b.id] == Bar(b.id, b.x, 1)

            # We can also do this on lists.
            b = Bar(1, 2, 3)
            insert_or_update!(v, b.id, b)
            @test length(v) == 1 && first(v) == b
            b = Bar(1, 10, missing)
            insert_or_update!(v, b.id, b)
            @test length(v) == 1 && first(v) == Bar(b.id, b.x, 3)
            b = Bar(1, 2, 2)
            insert_or_update!(v, b)
            @test length(v) == 1 && first(v) == b

            # If we need something fancier to index into a list, we can use a key function.
            b = Bar(1, 4, 5)
            insert_or_update!(v, b.id, b; key=x -> x.y)
            # We inserted a new element because we looked for an element with y == b.id.
            @test length(v) == 2 && last(v) == b
            # If we leave out the insert key, the key function is used on the value too.
            b = Bar(1, 0, 5)
            # We updated the value with y == b.y.
            insert_or_update!(v, b; key=x -> x.y)
            @test length(v) == 2 && last(v) == b

            # The updated value is returned.
            empty!(d)
            empty!(v)
            b = Bar(1, 2, 3)
            @test insert_or_update!(d, b) == b
            @test insert_or_update!(v, b) == b
            b = Bar(1, 3, missing)
            @test insert_or_update!(d, b) == Bar(1, 3, 3)
            @test insert_or_update!(v, b) == Bar(1, 3, 3)
        end
    end
end


text:
#=
Part of this work is derived from astropy/photutils. The relevant derivations
are considered under a BSD 3-clause license. =#

module Aperture

using TypedTables
using Transducers

export photometry,
       Subpixel,
       CircularAperture,
       CircularAnnulus,
       EllipticalAperture,
       EllipticalAnnulus,
       RectangularAperture,
       RectangularAnnulus

"""
    AbstractAperture{T} <: AbstractMatrix{T}

The abstract super-type for Apertures.

Apertures can be thought of as a cutout or stamp of a geometric shape with shading applied. For example, a circular aperture with a diameter of 3 pixels will require a 5x5 pixel grid (when perfectly on-grid) to represent.

```jldoctest ap1
julia> ap = CircularAperture(3, 3, 2.5)
5×5 CircularAperture{Float64} with indices 1:5×1:5:
 0.136857  0.769325  0.983232  0.769325  0.136857
 0.769325  1.0       1.0       1.0       0.769325
 0.983232  1.0       1.0       1.0       0.983232
 0.769325  1.0       1.0       1.0       0.769325
 0.136857  0.769325  0.983232  0.769325  0.136857
```

This is a useful way of thinking about apertures: if we have some data, we can weight the data with the aperture

```jldoctest ap1
julia> data = fill(2, 5, 5);

julia> idxs = map(intersect, axes(ap), axes(data)) |> CartesianIndices;

julia> weighted_cutout = data[idxs] .* ap[idxs]
5×5 Array{Float64,2}:
 0.273713  1.53865  1.96646  1.53865  0.273713
 1.53865   2.0      2.0      2.0      1.53865
 1.96646   2.0      2.0      2.0      1.96646
 1.53865   2.0      2.0      2.0      1.53865
 0.273713  1.53865  1.96646  1.53865  0.273713
```

Performing aperture photometry is merely summing the weighted cutout shown above.

```jldoctest ap1
julia> flux = sum(weighted_cutout)
39.269908169872416

julia> flux ≈ (π * 2.5^2) * 2 # area of circle times intensity of 2
true
```

What's interesting about the implementation of apertures, though, is they are lazy. This means there is no stored matrix of aperture values; rather, they are calculated on the fly as needed.

```jldoctest ap1
julia> axes(ap)
(1:5, 1:5)

julia> ap[-10, -10] # out-of-bounds, but calculated on the fly
0.0

julia> ap .* ones(5, 7) # broadcasts to eachindex(data), regardless of ap bound
5×7 Array{Float64,2}:
 0.136857  0.769325  0.983232  0.769325  0.136857  0.0  0.0
 0.769325  1.0       1.0       1.0       0.769325  0.0  0.0
 0.983232  1.0       1.0       1.0       0.983232  0.0  0.0
 0.769325  1.0       1.0       1.0       0.769325  0.0  0.0
 0.136857  0.769325  0.983232  0.769325  0.136857  0.0  0.0
```

This allows extremely efficient computation of aperture photometry from small to medium sized apertures.

```julia
julia> using BenchmarkTools

julia> @btime sum(idx -> \$ap[idx] * \$data[idx], \$idxs)
  1.097 μs (0 allocations: 0 bytes)
39.26990816987243
```

This is essentially the full implementation of [`photometry`](@ref), save for the packing of additional information into a tabular form.
"""
abstract type AbstractAperture{T}  <: AbstractMatrix{T} end

"""
    bounds(::AbstractAperture)

Return the (`xlow`, `xhigh`, `ylow`, `yhigh`) bounds for a given Aperture.
"""
bounds(::AbstractAperture)

center(ap::AbstractAperture) = ap.y, ap.x # greedy

"""
    size(::AbstractAperture)

Return (`ny`, `nx`) of the aperture.
"""
function Base.size(ap::AbstractAperture)
    xmin, xmax, ymin, ymax = bounds(ap)
    return ymax - ymin + 1, xmax - xmin + 1
end

function Base.axes(ap::AbstractAperture)
    xmin, xmax, ymin, ymax = bounds(ap)
    return ymin:ymax, xmin:xmax
end

"""
    Subpixel(ap, N=1) <: AbstractAperture

Use a subpixel quadrature approximation for pixel shading instead of exact geometric methods.

For any pixel laying on the border of `ap`, this alters the shading algorithm by breaking the border pixel up into `(N, N)` subpixels. The shading value is the fraction of these subpixels within the geometric border of `ap`.

Using a subpixel shading method is sometimes faster than exact methods at the cost of accuracy. For [`CircularAperture`](@ref) the subpixel method is only faster than the exact method for `N` ~ 7. for [`EllipticalAperture`](@ref) the cutoff is `N` ~ 12, and for [`RectangularAperture`](@ref) the cutoff is `N` ~ 20.

# Examples

```jldoctest
julia> ap = CircularAperture(3, 3, 2.5)
5×5 CircularAperture{Float64} with indices 1:5×1:5:
 0.136857  0.769325  0.983232  0.769325  0.136857
 0.769325  1.0       1.0       1.0       0.769325
 0.983232  1.0       1.0       1.0       0.983232
 0.769325  1.0       1.0       1.0       0.769325
 0.136857  0.769325  0.983232  0.769325  0.136857

julia> sub_ap = Subpixel(ap, 5)
5×5 Subpixel{Float64,CircularAperture{Float64}} with indices 1:5×1:5:
 0.12  0.76  1.0  0.76  0.12
 0.76  1.0   1.0  1.0   0.76
 1.0   1.0   1.0  1.0   1.0
 0.76  1.0   1.0  1.0   0.76
 0.12  0.76  1.0  0.76  0.12
```

!!! note
    `photutils` offers a `center` shading method which is equivalent to using the `Subpixel` method with 1 subpixel. To avoid unneccessary namespace cluttering, we simply instruct users to use `Subpixel(ap)` instead.
"""
struct Subpixel{T,AP<:AbstractAperture{T}} <: AbstractAperture{T}
    ap::AP
    N::Int
end

function Base.show(io::IO, c::Subpixel)
    print(io, "Subpixel(")
    print(io, c.ap)
    print(io, ", $(c.N))")
end

center(sp_ap::Subpixel) = center(sp_ap.ap)
bounds(sp_ap::Subpixel) = bounds(sp_ap.ap)
overlap(sp_ap::Subpixel, i, j) = overlap(sp_ap.ap, i, j)

Subpixel(ap::AbstractAperture) = Subpixel(ap, 1)

@enum OverlapFlag Inside Outside Partial

function Base.getindex(ap::AbstractAperture, i::Int, j::Int)
    flag = overlap(ap, i, j)
    flag === Outside && return 0.0
    flag === Inside && return 1.0
    cy, cx = center(ap)

    return partial(ap, j - cx, i - cy)
end

# This bypasses checking aperture axes for broadcasting
Broadcast.combine_axes(ap::AbstractAperture, arrs...) = Broadcast.combine_axes(arrs...)
Broadcast.combine_axes(arr, ap::AbstractAperture) = axes(arr)

###########

"""
    photometry(::AbstractAperture, data::AbstractMatrix, [error])
    photometry(::AbstractVector{<:AbstractAperture}, data::AbstractMatrix, [error])

Perform aperture photometry on `data` given aperture(s). If `error` (the pixel-wise standard deviation) is provided, will calculate sum error. If a list of apertures is provided the output will be a `TypedTables.Table`, otherwise a `NamedTuple`.

!!! tip
    This code is automatically multi-threaded. To take advantage of this please make sure `JULIA_NUM_THREADS` is set before starting your runtime.
"""
function photometry(ap::AbstractAperture, data::AbstractMatrix, error)
    cy, cx = center(ap)
    meta = (xcenter = cx, ycenter = cy)
    idxs = map(intersect, axes(ap), axes(data), axes(error))
    any(isempty, idxs) && return (meta..., aperture_sum = 0.0, aperture_sum_err = NaN)

    aperture_sum = sum(CartesianIndices(idxs) |> Map(idx -> ap[idx] * data[idx]))
    aperture_sum_var = sum(CartesianIndices(idxs) |> Map(idx -> ap[idx] * error[idx]^2))
    aperture_sum_err = sqrt(aperture_sum_var)

    return (meta..., aperture_sum = aperture_sum, aperture_sum_err = aperture_sum_err)
end

function photometry(ap::AbstractAperture, data::AbstractMatrix)
    cy, cx = center(ap)
    meta = (xcenter = cx, ycenter = cy)
    idxs = map(intersect, axes(ap), axes(data))
    any(isempty, idxs) && return (meta..., aperture_sum = 0.0)
    aperture_sum = sum(CartesianIndices(idxs) |> Map(idx -> ap[idx] * data[idx]))
    return (meta..., aperture_sum = aperture_sum)
end

function photometry(aps::AbstractVector{<:AbstractAperture}, data::AbstractMatrix, error)
    rows = tcollect(aps |> Map(ap -> photometry(ap, data, error)))
    return Table(rows)
end

function photometry(aps::AbstractVector{<:AbstractAperture}, data::AbstractMatrix)
    rows = tcollect(aps |> Map(ap -> photometry(ap, data)))
    return Table(rows)
end

include("circular.jl")
include("elliptical.jl")
include("rectangle.jl")
include("overlap.jl")
include("plotting.jl")

end


text:
# Assumes Toolbox already exists in EChart, which is true by composite type definition
function toolbox!(ec::EChart; 	show::Bool = true,
								#mark::Bool = true,
								restore::Bool = true,
								dataView::Bool = true,
								saveAsImage::Bool = true,
								chartTypes::AbstractVector = []
								)

	chartlookup = Dict(	"line" => "Line",
						"bar" => "Bar",
						"stack" => "Stack",
						"tiled" => "Tiled",
						"force" => "Force",
						"chord" => "Chord",
						"pie" => "Pie",
						"funnel" => "Funnel"
						)

	ec.toolbox.show = show

	restore ? ec.toolbox.feature["restore"] = Dict("show" => true, "title" => "Restore") : ec.toolbox.feature["restore"] = Dict("show" => false)
	dataView ? ec.toolbox.feature["dataView"] = Dict("show" => true, "title" => "Data View", "lang" => ["Data View", "Cancel", "Refresh"]) : ec.toolbox.feature["dataView"] = Dict("show" => false)
	saveAsImage ? ec.toolbox.feature["saveAsImage"] = Dict("show" => true, "title" => "Save As PNG") : ec.toolbox.feature["saveAsImage"] = Dict("show" => false)
	chartTypes != [] ? ec.toolbox.feature["magicType"] = Dict("show" => true, "type" => chartTypes, "title" => chartlookup) : nothing

	return ec

end


text:
using .Samsara, Conception
using UUIDs, Random

""" EuclideanNode is a (linked-list type) node with one one node in
each cartinal direction: N, S, E, W
 ~> Implies some Euclidean properties
"""
abstract type EuclideanNode <: Conception.Time.TemporalType end

_all_known_CartinalNodes = Samsara.EuclideanNode[]

""" LinkedCardinalNode is a node in a linked list, representing the states of a simulation.
Each node has a number of possible transitions, of which can be effectuated by events.
This is perfect for emulating a system for IV-learning.

Each link is of type Conception.TemporalType (like SAT):
    can be activated,
    have activation time,
    and be a part of a MuExS.
"""
mutable struct LinkedCardinalNode <: EuclideanNode
    _node::Conception.AbstractConcept
    _node_E::Union{Conception.Conception.AbstractConcept, Nothing}
    _node_W::Union{Conception.Conception.AbstractConcept, Nothing}

    function LinkedCardinalNode(id=nothing;
                        node_to_E =nothing,
                        node_to_W =nothing,
                        connected_T=nothing,
                        in_MuEx::Union{AbstractMuExS, Nothing}=nothing)
        # Id
        the_SAT = SAT(id, inMuExS=in_MuEx, connected_with=connected_T)

        LinkedCardinalNode(the_SAT, node_to_E=node_to_E, node_to_W=node_to_W)

        #= # Init
        the_node = new(the_SAT, node_to_E, node_to_W)

        #Conception.add_element_to_MuExS!(in_MuEx, the_node)

        push!(_all_known_CartinalNodes, the_node)
        return the_node
        =#
    end
    function LinkedCardinalNode(sat_arg::SAT;
                        node_to_E =nothing,
                        node_to_W =nothing)
        the_node = new(sat_arg, node_to_E, node_to_W)
        push!(_all_known_CartinalNodes, the_node)
        the_node
    end
end

function Conception.get_id(it::LinkedCardinalNode)
    get_id(it._node)
end

Conception.is_active(it::LinkedCardinalNode) = is_active(it._node)
Conception.exists_in_MuExS(caseNode::LinkedCardinalNode, state_set::Conception.MuExS) = Conception.exists_in_MuExS(caseNode._node, state_set)

" member of muex, forwarded to node "
member_of_MuExS(it::LinkedCardinalNode) = it._node.member_of_MuExS

function Conception.activate!(it::LinkedCardinalNode)
    activate!(it._node)
end
function Conception.set_connected_concept(it::LinkedCardinalNode, with::Conception.AbstractConcept)
    Conception.set_connected_concept(it._node, with)
end

""" show(IO, LinkeNode)
Print LinkeNode with node_to_W and node_to_E, with '|' if absent.
"""
function Base.show(io::IO, arg::LinkedCardinalNode)
    if isnothing(arg._node_E)
        text_nE = "|"
    else
        text_nE = string(get_id(arg._node_E))
    end
    if isnothing(arg._node_W)
        text_nW = "|"
    else
        text_nW = string(get_id(arg._node_W))
    end
    print(io, "[ "*text_nE*"→|"*string(get_id(arg._node))*"|→"*text_nW*" ]")
end

""" _set_node_to_W!(nodeA, nodeB)
Set note to the west of nodeA to become nodeB. Note that nodeB can be Nothing
"""
function _set_node_to_W!(nodeA::Conception.AbstractConcept,
                         nodeB::Union{Conception.AbstractConcept, Nothing})
    nodeA._node_W = nodeB
end

""" _set_node_to_E!(nodeA, nodeB)
Set note to the west of nodeA to become nodeB. Note that nodeB can be Nothing
"""
function _set_node_to_E!(nodeA::Conception.AbstractConcept,
                         nodeB::Union{Conception.AbstractConcept, Nothing})
    nodeA._node_E = nodeB
end

" Find link that contains SAT "
function find_link_with_SAT(node::SAT)
    for it in _all_known_CartinalNodes
        if it._node == node
            return it
        end
    end
    return nothing
end

""" east_of(node)
return the node that lies to the east of node::LinkedCardinalNode)
"""
function east_of(node::LinkedCardinalNode)
    return node._node_E
end
""" east_of(node) implemented for SAT (the content) """
function east_of(node::SAT)
    east_of(find_link_with_SAT(node))
end

""" node(it::LinkedCardinalNode)
return the actual SAT of a LinkedCardinalNode
"""
function the_node(it::LinkedCardinalNode)
    it._node
end

""" west_of(node)
return the node that lies to the west of node::LinkedCardinalNode)
"""
function west_of(node::LinkedCardinalNode)
    return node._node_W
end
""" west_of(node) implemented for SAT (the content) """
function west_of(node::SAT)
    west_of(find_link_with_SAT(node))
end

function linked_list_factory(N::Int; in_MuExS =nothing)
    retList = LinkedCardinalNode[]
    # first item      ( N > 0 )
    previous_node = LinkedCardinalNode("n1", in_MuEx=in_MuExS)
    push!(retList, previous_node)
    N == 1 && return retList
    # .. then the rest ( N > 1 )
    #  split is made for the sake of registering _node_E
    for i in 2:N
        the_node = LinkedCardinalNode(SAT("n"*string(i), inMuExS=in_MuExS))
        #the_node = LinkedCardinalNode("n"*string(i), in_MuEx=in_MuExS)
        push!(retList, the_node)
        _set_node_to_W!(the_node, previous_node)
        _set_node_to_E!(previous_node, the_node)
        previous_node = the_node
    end
    retList
end

max_stars_repo_name: leikanger/Samsara

max_issues_repo_name: leikanger/Samsara

max_forks_repo_name: leikanger/Samsara


text:
using Pkg; Pkg.activate(pwd())
using DEMC
using Random
using LinearAlgebra
using Distributions
Random.seed!(319531501)

# generate data
nobs = 1000
npar = 25
Σ = ones(npar, npar)
for i = 1:npar
    for iprime = 1:i
        if i != iprime
            Σ[i, iprime] = 0.25
            Σ[iprime, i] = 0.25
        end
    end
end

DataDistribution = MvNormal(zeros(npar), Σ)
X = zeros(nobs, npar+1)
X[:,1] .= 1.
for iobs = 1:nobs
    X[iobs, 2:end] .= rand(DataDistribution)
end

β = 1 .+ rand(npar+1)*3
y = X*β + randn(nobs)

# log objective function - for correct std errors need to use optimal weights
log_obj(b) = -0.5 * sum((y.-X*b).^2)

# set up of DEMCz chain
ndim = length(β)
opts = DEMC.demcopt(ndim)
opts.blockindex = [1:ndim] # parameter blocks: here choose all parameters to be updated simultaenously
opts.Nblocks = length(opts.blockindex)
opts.eps_scale = 1e-5*ones(ndim) # scale of random error around DE update
opts.γ = 2. # scale of DE update, 2.38 is the "optimal" number for a normal distribution
opts.N = 5 # number of chains
opts.K = 10 # every K steps add current N draws to Z
# Number of iterations in Chain
opts.Ngeneration = 100000
opts.autostop = :Rhat
opts.autostop_every = 2000
opts.autostop_Rhat = 1.1
opts.print_step = 1000
opts.verbose = false

Z = randn((10*ndim, ndim))

# run chain
mc, Z = DEMC.demcz_sample(log_obj, Z, opts)

# drop first half of chain
N, Npar, Ntot = size(mc.chain)
keep = Int(Ntot-opts.autostop_every)+1:Ntot
Ngen_burned = length(keep)
chain_burned = mc.chain[:,:,keep]
chainflat = DEMC.flatten_chain(chain_burned, N, Ngen_burned, Npar)'

# covariance of estimates
b, Σb = DEMC.mean_cov_chain(chain_burned)
# correct std errors for incorrect weighting matrix
# OLS
bols = (X'*X)\(X'*y)
ui = y .- X*bols
s2hat = (1/(nobs-npar)) * sum(ui.^2)
Σols = inv(X'X)*s2hat

# difference simulated to OLS
println("\n MCMC estimates: ", b, "\n dist to OLS: ", b .- bols)

accept_ratio, Rhat = DEMC.convergence_check(mc.chain, mc.log_obj, "none"; verbose = false, parnames = [])


text: """ Add from private repo folder once "copyright" situation is resolved"""

max_stars_repo_name: mforets/InvariantSets

max_issues_repo_name: mforets/InvariantSets


text:
# Out-of-place C = A*B
@inline function gemm(A::AbstractMatrix, B::AbstractMatrix)
    m, n, p = size(A,1), size(B,2), size(A,2)

    C = Matrix{eltype(A)}(undef, m, n)

    @turbo for i in 1:m, k in 1:n
        aux = zero(eltype(A[1]))
        for j in 1:p
            aux += A[i,j]*B[j,k]
        end
        C[i,k] = aux
    end
    return C

end

# threaded version of vectorized gemm
@inline function gemmt(A::AbstractMatrix, B::AbstractMatrix)
    m, n, p = size(A,1), size(B,2), size(A,2)

    C = Matrix{eltype(A)}(undef, m, n)

    @tturbo for i in 1:m, k in 1:n
        aux = zero(eltype(A[1]))
        for j in 1:p
            aux += A[i,j]*B[j,k]
        end
        C[i,k] = aux
    end
    return C

end

# In-place C = A*B
@inline function gemm!(C::AbstractArray{T},A::AbstractArray{T},B::AbstractArray{T}) where T

    m, n, p = size(A,1), size(B,2), size(A,2)
    @turbo for i in 1:m, k in 1:n
        aux = zero(eltype(A[1]))
        for j in 1:p
            aux += A[i,j]*B[j,k]
        end
        C[i,k] = aux
    end
end

@inline function gemmt!(C::AbstractArray{T},A::AbstractArray{T},B::AbstractArray{T}) where T

    m, n, p = size(A,1), size(B,2), size(A,2)
    @tturbo for i in 1:m, k in 1:n
        aux = zero(eltype(A[1]))
        for j in 1:p
            aux += A[i,j]*B[j,k]
        end
        C[i,k] = aux
    end
end

# x + p*y -> y (vec+scalar*vec)
function xpy!(y::Vector{T}, x::Vector{T}, p::T) where T
    @tturbo for i in eachindex(y)
        y[i] = x[i] + p * y[i]
    end
end

# x + p*y -> y (vec+scalar*vec)
# xpy(y::CuArray, x::CuArray, p::Number) = x .+ p .* y

# out-of-place Sparse CSC matrix times dense vector
function spmv(A::AbstractSparseMatrix, x::DenseVector)
    out = zeros(eltype(x), A.m)
    Aj = rowvals(A)
    nzA = nonzeros(A)
    @inbounds for col in 1:A.n
        xj = x[col]
        @fastmath for j in nzrange(A, col)
            out[Aj[j]] += nzA[j]*xj
        end
    end
    out
end

function _spmvt!(out, A, Aj, xj, nzA, col)
    @fastmath for j in nzrange(A, col)
        out[Aj[j]] += nzA[j]*xj
    end
end

# in-of-place Sparse CSC matrix times dense vector
function spmv!(out::DenseVector, A::AbstractSparseMatrix, x::DenseVector)
    fill!(out, zero(eltype(out)))
    Aj = rowvals(A)
    nzA = nonzeros(A)
    @inbounds for col in 1:A.n
        xj = x[col]
        @fastmath for j in nzrange(A, col)
            out[Aj[j]] += nzA[j]*xj
        end
    end
    out
end

mynorm(r::Vector) = sqrt(mydot(r,r))

# Dot product
function mydot(a::AbstractArray{T}, b::AbstractArray{T}) where {T}
    n = zero(T)
    @turbo for i in eachindex(a)
        n += a[i]*b[i]
    end
    n
end

# Threaded dot product
function mydott(a::AbstractArray{T}, b::AbstractArray{T}) where {T}
    n = zero(T)
    @tturbo for i in eachindex(a)
        n += a[i]*b[i]
    end
    n
end

max_stars_repo_name: albert-de-montserrat/Persephone

max_issues_repo_name: albert-de-montserrat/Persephone

max_forks_repo_name: albert-de-montserrat/Persephone


text:
struct S a::T end
struct T end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
function download()
  # File Names
  DataDir = joinpath(dirname(@__FILE__),"..","data")
  CompressFile = joinpath(DataDir,"bindata.tar.gz")
  CifarOrigDir = joinpath(DataDir,"cifar-10-batches-bin")
  CifarFinalDir = joinpath(DataDir,"bin")

  # Download Compressed Files
  download("http://cbcl.mit.edu/projects/cbcl/software-datasets/faces.tar.gz",CompressFile)

  # # Decompress Dataset
  cd(DataDir)
  run(`tar -zxvf $CompressFile`)
  run(`tar -zxvf face.test.tar.gz`)
  run(`tar -zxvf face.train.tar.gz`)
end


text:
for T in (:CLCL, :CLOP, :OPCL, :OPOP)
   @eval begin
       function (sqrt){Q<:Qualia,C<:Clay}(a::Flex{$T,Q,C})
           lo = (sqrt)(a.lo, RoundDown)
           hi = (sqrt)(a.hi, RoundUp)
           Flex{$T,Q,C}(lo,hi)
       end
   end
end


text:
export parse_newick!, nw_parse_children, nw_parse_name, read_newick, read_tree
export fasta2tree!, seq2num

let n::Int64=0
	global increment_n() = (n+=1)
	global reset_n() = (n=0)
end

"""
	read_tree(nw_file::String; NodeDataType=EvoData)

Read Newick file `nw_file` and create a `Tree{NodeDataType}` object from it.
`NodeDataType` must be a subtype of `TreeNodeData`, and must have a *callable default outer constructor*. In other words, the call `NodeDataType()` must exist and return a valid instance of `NodeDataType`. This defaults to `EvoData`.
"""
function read_tree(nw_file::String; NodeDataType=EvoData)
	# println("Checking Tree")
	@time tree = node2tree(read_newick(nw_file; NodeDataType=NodeDataType))
	check_tree(tree)
	return tree
end

"""
	read_newick(nw_file::String; NodeDataType=EvoData)

Read Newick file `nw_file` and create a graph of `TreeNode{NodeDataType}` objects in the process. Return the root of said graph. `node2tree` or `read_tree` must be called to obtain a `Tree{NodeDataType}` object.
`NodeDataType` must be a subtype of `TreeNodeData`, and must have a *callable default outer constructor*. In other words, the call `NodeDataType()` must exist and return a valid instance of `NodeDataType`. This defaults to `EvoData`.
"""
function read_newick(nw_file::String; NodeDataType=EvoData)
	@assert NodeDataType <: TreeNodeData
	f = open(nw_file)
	nw = readlines(f)
	close(f)
	if length(nw) > 1
		error("File $nw_file has more than one line.")
	elseif length(nw) == 0
		error("File $nw_file is empty")
	end
	nw = nw[1]
	if nw[end] != ';'
		error("File $nw_file does not end with ';'")
	end
	nw = nw[1:end-1]

	reset_n()
	root = parse_newick(nw, NodeDataType=NodeDataType)
	return root
end
function parse_newick(nw::String; NodeDataType=EvoData)
	root = TreeNode(NodeDataType())
	parse_newick!(nw, root, NodeDataType)
	root.isroot = true # Rooting the tree with outer-most node of the newick string
	return root
end

"""
	parse_newick!(nw::String, root::TreeNode)

Parse the tree contained in Newick string `nw`, rooting it at `root`.
"""
function parse_newick!(nw::String, root::TreeNode, NodeDataType)

	# Setting isroot to false. Special case of the root is handled in main calling function
	root.isroot = false
	# Getting label of the node, after last parenthesis
	parts = map(x->String(x), split(nw, ")"))
	lab, tau = nw_parse_name(String(parts[end]))
	if lab == ""
		lab = "NODE_$(increment_n())"
	end
	root.label, root.data.tau = (NodeDataType == LBIData && ismissing(tau) ? (lab,0.) : (lab,tau))

	if length(parts) == 1 # Is a leaf. Subtree is empty
		root.isleaf = true
	else # Has children
		root.isleaf = false
		if parts[1][1] != '('
			println(parts[1][1])
			error("Parenthesis mismatch.")
		else
			parts[1] = parts[1][2:end] # Removing first bracket
		end
		children = join(parts[1:end-1], ")") # String containing children, now delimited with ','
		l_children = nw_parse_children(children) # List of children (array of strings)

		for sc in l_children
			nc = TreeNode(NodeDataType())
			parse_newick!(sc, nc, NodeDataType) # Will set everything right for subtree corresponding to nc
			nc.anc = root
			push!(root.child, nc)
		end
	end

end

"""
	nw_parse_children(s::String)

Idea from http://stackoverflow.com/a/26809037
Split a string of children in newick format to an array of strings.

## Example
`"A,(B,C),D"` --> `["A","(B,C)","D"]`
"""
function nw_parse_children(s::String)
	parcount = 0
	l_children = []
	current = ""
	cstart = 1
	cend = 1
	for (i,c) in enumerate("$(s),")
		if c == ',' && parcount == 0
			cend = i-1
			push!(l_children, s[cstart:cend])
			cstart = i+1
		else
			if c == '('
				parcount +=1
			elseif c == ')'
				parcount -=1
			end
		end
	end
	return l_children
end

"""
	nw_parse_name(s::String)

Parse Newick string of child into name and time to ancestor. Default value for missing time is `missing`.
"""
function nw_parse_name(s::String)
	temp = split(s, ":")
	if occursin(':', s) # Node has a time
		if length(temp) == 2 # Node also has a name, return both
			tau = (tau = tryparse(Float64,temp[2]); typeof(tau)==Nothing ? missing : tau) # Dealing with unparsable times
			return string(temp[1]), tau
		else # Return empty name
			tau = (tau = tryparse(Float64,temp[1]); typeof(tau)==Nothing ? missing : tau) # Dealing with unparsable times
			return "", tau
		end
	else # Node does not have a time, return string as name
		return s, missing
	end
end

"""
	fasta2tree!(tree::Tree{EvoData}, fastafile::String ; seqtype=:nucleotide)

Read `fastafile` and stores sequences in nodes of `tree` with corresponding label.
Implemented sequence types are `:nucleotide` and `:binary`.
"""
function fasta2tree!(tree::Tree{EvoData}, fastafile::String ; seqtype=:nucleotide)
	for (name,seq) in FastaReader(fastafile)
		if haskey(tree.lnodes, name)
			# storeseq!(tree.lnodes[name], seq, seqtype)
			tree.lnodes[name].data.sequence = [c for c in seq]
		end
	end
end

"""
	storeseq!(node::TreeNode{EvoData}, seq, seqtype)

Store sequence `seq` at node `node`.
"""
function storeseq!(node::TreeNode{EvoData}, seq, seqtype)
	if seqtype==:nucleotide
		numseq = seq2num(seq, seqtype)
		node.data.sequence = numseq
		node.data.q = 5
	elseif seqtype==:binary
		if isa(seq, String)
			node.data.sequence = map(x->parse(Int64, x), collect(seq))
		elseif isa(seq, Array{Int64})
			node.data.sequence = seq
		else
			error("Unrecognized format for sequence `seq`.")
		end
		node.data.q = 2
	else
		println("Accepted sequence types are `:binary` and `:nucleotide`.")
		error("Unknown seq-type.")
	end
end

"""
"""
function seq2num(seq::String, seqtype)
	if seqtype == :nucleotide
		mapping = "ACGT-NWSMKRYBDHV"
		numseq = zeros(Int64, length(seq))
		for (i,c) in enumerate(seq)
			num = findall(x->x==c, mapping)
			if isempty(num)
				@warn "String $seq could not be matched to nucleotides -- position ($i, $c)"
			else
				numseq[i] = findall(x->x==c, mapping)[1]
				# numseq[i] = findall(x->x==c, mapping)[1]<5 ? findall(x->x==c, mapping)[1] : 5
			end
		end
	end
	return numseq
end

"""
"""
function seq2num(c::Char, seqtype)
	if seqtype == :nucleotide
		mapping = "ACGT-NWSMKRYBDHV"
		num = findall(x->x==c, mapping)
		if isempty(num)
			@warn "String $seq could not be matched to nucleotides -- position ($i, $c)"
		else
			numseq = findall(x->x==c, mapping)[1]
		end
	end
	return numseq
end

max_stars_repo_name: PierreBarrat/FluPredictibility

max_issues_repo_name: PierreBarrat/FluPredictibility

max_forks_repo_name: PierreBarrat/FluPredictibility


text:
"""
    solve!(solution::AbstractRungeKuttaSolution, problem::AbstractInitialValueProblem, solver::AbstractRungeKuttaSolver) :: AbstractRungeKuttaSolution

returns the [`AbstractRungeKuttaSolution`](@ref) of an [`AbstractInitialValueProblem`](@ref).
"""
function NSDEBase.solve!(solution::AbstractRungeKuttaSolution, problem::AbstractInitialValueProblem, solver::AbstractRungeKuttaSolver)
    cache = RungeKuttaCache(problem, solver)
    @↓ u0, (t0, tN) ← tspan = problem
    @↓ u, t = solution
    @↓ n = cache
    u0_T = eltype(u0)
    u0_L = length(u0)
    N = length(t)
    # WHILE instead of FOR loop -> adaptive methods
    while n < N && t[n] < tN
        step!(solution, problem, solver, cache)
        adaptive_step!(solution, solver, cache)
        @↓ n = cache
        # Stop if update is too small -> adaptive methods
        if n > 1 && t[n] ≈ t[n-1]
            break
        end
        # Increase size of output vectors -> adaptive methods
        # if n == N && t[n] < tN
        #     N += 1
        #     push!(u, Vector{u0_T}(undef, u0_L))
        #     resize!(t, N)
        # end
    end
    # Decrease size of output vectors -> adaptive methods
    if n < N
        resize!(u, n)
        resize!(t, n)
    end
    return solution
end

"""
    solve(problem::AbstractInitialValueProblem, solver::AbstractRungeKuttaSolver; savestages::Bool = false) :: AbstractRungeKuttaSolution

returns the [`AbstractRungeKuttaSolution`](@ref) of an [`AbstractInitialValueProblem`](@ref). `savestages` flags when to save all stages into `solution.k`.
"""
function NSDEBase.solve(problem::AbstractInitialValueProblem, solver::AbstractRungeKuttaSolver; savestages::Bool = false)
    solution = RungeKuttaSolution(problem, solver; savestages=savestages)
    NSDEBase.solve!(solution, problem, solver)
    return solution
end