text:
export directConv, directConv!, directConv2D!, directCrossCorrelation,directCrossCorrelation2D
export BoundaryExtension, ZeroPaddingBE, ConstantBE, PeriodicBE, MirrorBE

# first index
const tilde_i0 = Int(1)

#+BoundaryExtension
#
# Used for tag dispatching, parent of available boundary extensions
#
#!subtypes(BoundaryExtension)
#
abstract type BoundaryExtension end

#+BoundaryExtension
struct ZeroPaddingBE <: BoundaryExtension end
#+BoundaryExtension
struct ConstantBE <: BoundaryExtension end
#+BoundaryExtension
struct PeriodicBE <: BoundaryExtension end
#+BoundaryExtension
struct MirrorBE <: BoundaryExtension end

#+BoundaryExtension,Internal
#
# Range scaling
#
# *Caveat:*
# We do not use Julia =*= operator as it returns a step range:
#!r=6:8
#!-2*r
# What we need is:
#!scale(-2,r)
#
function scale(λ::Int,Ω::UnitRange{Int})
    ifelse(λ>0,
           UnitRange(λ*first(Ω),λ*last(Ω)),
           UnitRange(λ*last(Ω),λ*first(Ω)))
end

#+BoundaryExtension,Internal
#
# In
# $$
# \gamma[k]=\sum\limits_{i\in\Omega^\alpha}\alpha[i]\beta[k+\lambda i],\text{ with }\lambda\in\mathbb{Z}^*
# $$
# the computation of $\gamma[k],\ k\in\Omega^\gamma$ is splitted into two parts:
#  - one part $\Omega^\gamma \cap \Omega^\gamma_1$ *free of boundary effect*,
#  - one part $\Omega^\gamma \setminus \Omega^\gamma_1$ *that requires boundary extension* $\tilde{\beta}=\Phi(\beta,k)$
#
# *Example:*
#!DirectConvolution.compute_Ωγ1(-1:2,-2,1:20)
function compute_Ωγ1(Ωα::UnitRange{Int},
                     λ::Int,
                     Ωβ::UnitRange{Int})

    λΩα = scale(λ,Ωα)

    UnitRange(first(Ωβ)-first(λΩα),
              last(Ωβ)-last(λΩα))
end

#+BoundaryExtension,Internal
#
# Left relative complement
#
# $$
# (A\setminus B)_{\text{Left}}=[  l(A), \min{(u(A),l(B)-1)} ]
# $$
#
# *Example:*
#!DirectConvolution.relativeComplement_left(1:10,-5:5)
#
# $(A\setminus B)=\{6,7,8,9,10\}$ and the left part (elements that are
# $\in A$ but on the left side of $B$) is *empty*.
function relativeComplement_left(A::UnitRange{Int},
                                 B::UnitRange{Int})
    UnitRange(first(A),
              min(last(A),first(B)-1))
end

#+BoundaryExtension,Internal
#
# Left relative complement
#
# $$
# (A\setminus B)_{\text{Right}}=[ \max{(l(A),u(B)+1)}, u(A) ]
# $$
#
# *Example:*
#!DirectConvolution.relativeComplement_right(1:10,-5:5)
#
# $(A\setminus B)=\{6,7,8,9,10\}$ and the right part (elements that are
# $\in A$ but on the right side of $B$) is $\{6,7,8,9,10\}$
function relativeComplement_right(A::UnitRange{Int},
                                  B::UnitRange{Int})
    UnitRange(max(first(A),last(B)+1),
              last(A))
end

#+BoundaryExtension,Internal
#
#!r=-5:10
#!hcat(r,map(x->DirectConvolution.boundaryExtension([1:3;],x,ZeroPaddingBE),r))'
#
function boundaryExtension(β::AbstractArray{T,1},
                           k::Int,
                           ::Type{ZeroPaddingBE}) where {T <: Number}
    kmin = tilde_i0
    kmax = length(β) + kmin - 1

    if (k>=kmin)&&(k<=kmax)
        β[k]
    else
        T(0)
    end
end

#+BoundaryExtension,Internal
#
#!r=-5:10
#!hcat(r,map(x->DirectConvolution.boundaryExtension([1:3;],x,ConstantBE),r))'
#
function boundaryExtension(β::AbstractArray{T,1},
                           k::Int,
                           ::Type{ConstantBE}) where {T <: Number}
    kmin = tilde_i0
    kmax = length(β) + kmin - 1

    if k<kmin
        β[kmin]
    elseif k<=kmax
        β[k]
    else
        β[kmax]
    end
end

#+BoundaryExtension,Internal
#
#!r=-5:10
#!hcat(r,map(x->DirectConvolution.boundaryExtension([1:3;],x,PeriodicBE),r))'
#
function boundaryExtension(β::AbstractArray{T,1},
                           k::Int,
                           ::Type{PeriodicBE}) where {T <: Number}
    kmin = tilde_i0
    kmax = length(β) + kmin - 1

    β[kmin+mod(k-kmin,1+kmax-kmin)]
end

#+BoundaryExtension,Internal
#
#!r=-5:10
#!hcat(r,map(x->DirectConvolution.boundaryExtension([1:3;],x,MirrorBE),r))'
#
function boundaryExtension(β::AbstractArray{T,1},
                           k::Int,
                           ::Type{MirrorBE}) where {T <: Number}
    kmin = tilde_i0
    kmax = length(β) + kmin - 1

    β[kmax-abs(kmax-kmin-mod(k-kmin,2*(kmax-kmin)))]
end

#+Convolution,Internal
function directConv!(tilde_α::AbstractArray{T,1},
                     α_offset::Int,
                     λ::Int,
                     β::AbstractArray{T,1},
                     γ::AbstractArray{T,1},
                     Ωγ::UnitRange{Int},
                     ::Type{LeftBE}=ZeroPaddingBE,
                     ::Type{RightBE}=ZeroPaddingBE;
                     accumulate::Bool=false)::Nothing where {T <: Number,
                                                          LeftBE <: BoundaryExtension,
                                                          RightBE <: BoundaryExtension}
    # Sanity check
    @assert λ!=0
    @assert (first(Ωγ)>=1)&&(last(Ωγ)<=length(γ))

    # Initialization
    Ωα = filter_range(length(tilde_α),α_offset)
    Ωβ = UnitRange(1,length(β))
    tilde_Ωα = 1:length(Ωα)

    if !accumulate
        for k in Ωγ
            γ[k]=0
        end
    end

    rΩγ1=intersect(Ωγ,compute_Ωγ1(Ωα,λ,Ωβ))

    # rΩγ1 part: no boundary effect
    #
    β_offset = λ*(first(Ωα)-tilde_i0)
    for k in rΩγ1
        @simd for i in tilde_Ωα
            @inbounds γ[k]+=tilde_α[i]*β[k+λ*i+β_offset]
        end
    end

    # Left part
    #
    rΩγ1_left = relativeComplement_left(Ωγ,rΩγ1)

    for k in rΩγ1_left
        for i in tilde_Ωα
            γ[k]+=tilde_α[i]*boundaryExtension(β,k+λ*i+β_offset,LeftBE)
        end
    end

    # Right part
    #
    rΩγ1_right = relativeComplement_right(Ωγ,rΩγ1)

    for k in rΩγ1_right
        for i in tilde_Ωα
            γ[k]+=tilde_α[i]*boundaryExtension(β,k+λ*i+β_offset,RightBE)
        end
    end

    nothing
end

# +Convolution L:directConv_details
# Computes a convolution.
#
# Inplace modification of $\gamma[k], k\in\Omega_\gamma$.
# $$
# \gamma[k]=\sum\limits_{i\in\Omega^\alpha}\alpha[i]\beta[k+\lambda i],\text{ with }\lambda\in\mathbb{Z}^*
# $$
# If $k\notin \Omega_\gamma$, $\gamma[k]$ is unmodified.
#
# If *accumulate=false* then an erasing step $\gamma[k]=0,
# k\in\Omega_\gamma$ is performed before computation.
#
# If $\lambda=-1$ you compute a convolution, if $\lambda=+1$ you
# compute a cross-correlation.
#
# *Example:*
# !β=[1:15;];
# !γ=ones(Int,15);
# !α=LinearFilter([0,0,1],0);
# !directConv!(α,1,β,γ,5:10);
# !hcat([1:length(γ);],γ)'
#
function directConv!(α::LinearFilter{T},
                     λ::Int,

                     β::AbstractArray{T,1},

                     γ::AbstractArray{T,1},
                     Ωγ::UnitRange{Int},

                     ::Type{LeftBE}=ZeroPaddingBE,
                     ::Type{RightBE}=ZeroPaddingBE;

                     accumulate::Bool=false)::Nothing where {T <: Number,
                                                          LeftBE <: BoundaryExtension,
                                                          RightBE <: BoundaryExtension}

    directConv!(fcoef(α),
                offset(α),
                λ,

                β,

                γ,
                Ωγ,

                LeftBE,
                RightBE,

                accumulate=accumulate)

    nothing
end

#+Convolution
#
# Computes a convolution.
#
# Convenience function that allocate $\gamma$ and compute all its
# component using [[directConv_details][]]
#
# *Returns:* $\gamma$ a created vector of length identical to the $\beta$ one.
#
# *Example:*
#!β=[1:15;];
#!γ=ones(Int,15);
#!α=LinearFilter([0,0,1],0);
#!γ=directConv(α,1,β);
#!hcat([1:length(γ);],γ)'
#
function directConv(α::LinearFilter{T},

                    λ::Int64,

                    β::AbstractArray{T,1},

                    ::Type{LeftBE}=ZeroPaddingBE,
                    ::Type{RightBE}=ZeroPaddingBE) where {T <: Number,
                                                          LeftBE <: BoundaryExtension,
                                                          RightBE <: BoundaryExtension}
    γ = Array{T,1}(undef,length(β))

    directConv!(α,
                λ,

                β,

                γ,
                UnitRange(1,length(γ)),

                LeftBE,
                RightBE,

                accumulate=false)

    return γ
end

#+Convolution
#
# Computes a convolution.
#
# This is a convenience function where $\lambda=-1$
#
# *Returns:* $\gamma$ a created vector of length identical to the $\beta$ one.
#
function directConv(α::LinearFilter{T},
                    β::AbstractArray{T,1},

                    ::Type{LeftBE}=ZeroPaddingBE,
                    ::Type{RightBE}=ZeroPaddingBE) where {T <: Number,
                                                          LeftBE <: BoundaryExtension,
                                                          RightBE <: BoundaryExtension}

    return directConv(α,-1,β,LeftBE,RightBE)
end

#+Convolution
#
# Computes a cross-correlation
#
# This is a convenience function where $\lambda=+1$
#
# *Returns:* $\gamma$ a created vector of length identical to the $\beta$ one.
#
function directCrossCorrelation(α::LinearFilter{T},
                                β::AbstractArray{T,1},

                                ::Type{LeftBE}=ZeroPaddingBE,
                                ::Type{RightBE}=ZeroPaddingBE) where {T <: Number,
                                                                      LeftBE <: BoundaryExtension,
                                                                      RightBE <: BoundaryExtension}

    return directConv(α,+1,β,LeftBE,RightBE)
end

# 2D

# +Convolution L:directConv2D_inplace
# Computes a 2D (separable) convolution.
#
# For general information about parameters, see [[directConv_details][]]
#
# α_I must be interpreted as filter for *running index I*
#
# CAVEAT: the result overwrites β
#
# TODO: @parallel
function directConv2D!(α_I::LinearFilter{T},
                       λ_I::Int,
                       α_J::LinearFilter{T},
                       λ_J::Int,

                       β::AbstractArray{T,2},

                       min_I_BE::Type{<:BoundaryExtension}=ZeroPaddingBE,
                       max_I_BE::Type{<:BoundaryExtension}=ZeroPaddingBE,
                       min_J_BE::Type{<:BoundaryExtension}=ZeroPaddingBE,
                       max_J_BE::Type{<:BoundaryExtension}=ZeroPaddingBE)::Nothing where {T<:Number}

    γ=similar(β)

    (n,m)=size(β)
    α_I_coef=fcoef(α_I)
    α_I_offset=offset(α_I)
    α_J_coef=fcoef(α_J)
    α_J_offset=offset(α_J)
    Ωγ_I = 1:n
    Ωγ_J = 1:m

    # i running (for filter)
    for j in 1:m

        directConv!(α_I_coef,
                    α_I_offset,
                    λ_I,

                    view(β,:,j),

                    view(γ,:,j),
                    Ωγ_I,

                    min_I_BE,
                    max_I_BE,

                    accumulate=false)
    end

    # j running (for filter)
    for i in 1:n

        directConv!(α_J_coef,
                    α_J_offset,
                    λ_J,

                    view(γ,i,:),

                    view(β,i,:),
                    Ωγ_J,

                    min_J_BE,
                    max_J_BE,

                    accumulate=false)
    end

    nothing
end

# +Convolution
# Computes a 2D cross-correlation
#
# This is a wrapper that calls [[directConv2D_inplace][]]
#
# *Note:* β is not modified, instead the function returns the result.
#
function directCrossCorrelation2D(α_I::LinearFilter{T},
                                  α_J::LinearFilter{T},

                                  β::AbstractArray{T,2},

                                  min_I_BE::Type{<:BoundaryExtension}=ZeroPaddingBE,
                                  max_I_BE::Type{<:BoundaryExtension}=ZeroPaddingBE,
                                  min_J_BE::Type{<:BoundaryExtension}=ZeroPaddingBE,
                                  max_J_BE::Type{<:BoundaryExtension}=ZeroPaddingBE)::Array{T,2} where {T<:Number}

    γ=similar(β)
    γ.=β

    directConv2D!(α_I,+1,α_J,+1,γ,
                  min_I_BE,
                  max_I_BE,
                  min_J_BE,
                  max_J_BE)

    return γ
end


text:
# load code
include("flow_u1_dm_self.jl")
include("flow_u1_dm_s.jl")
include("flow_u1_dm_t.jl")
include("flow_u1_dm_u.jl")

ext: jl


text:
######################################################################
# Extended copying - customizable behavior of copying structures.
# -----
# Licensed under MIT License
export FieldCopyOverride, xcopy, @xcopy, xcopy_override, @xcopy_override

struct FieldCopyOverride{S} end
struct NoOverride end
const nooverride = NoOverride()

"""`xcopy(tpl, kwargs...)`

 Extended `Base.copy` method which allows setting a field directly rather than copying it from `tpl`.

 See also [`@xcopy`](@ref), [`xcopyargs`](@ref), [`xcopy_override`](@ref), [`xcopy_construct`](@ref).

 # Example

 ```julia
 struct Foo{A}
    a::A
 end

 xcopy(Foo{Int32}(24), a=42) # == Foo{Int32}(42)
 ```
 """
function xcopy(tpl::T; kwargs...) where T
    fnames  = fieldnames(T)
    fvalues = Dict{Symbol, Any}()

    for fname ∈ fnames
        fvalues[fname] = haskey(kwargs, fname) ? kwargs[fname] : xcopy_override(tpl, FieldCopyOverride{fname}())
    end

    unknown = setdiff(keys(kwargs), fnames)
    if !isempty(unknown)
        throw(DomainError(unknown, "Unknown fields"))
    end

    xcopy_construct(tpl, (fvalues[fname] for fname ∈ fnames)...)
end
"""`@xcopy(T)`

 Convenience macro to define `Base.copy(tpl::T; kwargs...) = xcopy(tpl, kwargs...)`

 See [`xcopy`](@ref), [`xcopyargs`](@ref), [`xcopy_override`](@ref)"""
macro xcopy(T)
    esc(:(Base.copy(tpl::$T; kwargs...) = ExtraFun.xcopy(tpl; kwargs...)))
end

"""`xcopyargs(ArgsT::Type, tpl; kwargs...)`

 Like `xcopy`, except `kwargs` must match fields of `ArgsT`. A copy of `tpl` will then be created like such: `typeof(tpl)(argsobj)`
 where `argsobj::ArgsT` is constructed through `tpl` and `kwargs`.

 By default, the algorithm will attempt to retrieve an equally named field from `tpl` to construct `argsobj` with. This
 behavior can be overridden with [`xcopy_override`](@ref).

 See also [`xcopy`](@ref), [`xcopy_override`](@ref), [`@xcopyoverride`](@ref), [`xcopy_construct`](@ref)"""
function xcopyargs(ArgsT::Type, tpl::T; kwargs...) where T
    fnames  = fieldnames(ArgsT)
    fvalues = Dict{Symbol, Any}((fname => xcopy_override(tpl, FieldCopyOverride{fname}()) for fname ∈ fnames))

    for (fname, fvalue) ∈ kwargs
        @assert fname ∈ fnames "Unknown field `$fname`"
        fvalues[fname] = fvalue
    end

    xcopy_construct(tpl, ArgsT((fvalues[fname] for fname ∈ fnames)...))
end
"""`@xcopyargs(T, A)`

 Convenience macro to define `Base.copy(tpl::T; kwargs...) = xcopyargs(A, tpl, kwargs...)`

 See [`xcopy`](@ref), [`xcopyargs`](@ref), [`xcopy_override`](@ref)"""
macro xcopyargs(T, A)
    esc(:(Base.copy(tpl::$T; kwargs...) = ExtraFun.xcopyargs(A, tpl; kwargs...)))
end

"""`xcopy_override(tpl, ::FieldCopyOverride{S}) where S`

 Defines the copied value of field `S` in `tpl`. Defaults to `Base.copy(getproperty(tpl, S))`.

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

 # Example

 ```julia
 struct Foo
    value::Int
 end

 struct Bar
    foo::Foo
 end

 xcopy_override(bar::Bar, ::FieldCopyOverride{:foo}) = Foo(bar.foo.value + 1)
 xcopy(Bar(Foo(24))) # == Bar(Foo(25))
 """
xcopy_override(tpl, ::FieldCopyOverride{S}) where S = (value = getproperty(tpl, S); hassignature(copy, typeof(value)) ? copy(value) : value)
"""`@xcopyoverride(T, S, expr)`

 Convenience macro to define `xcopy_override(tpl::T, ::FieldCopyOverride{S}) = expr`.

 See [`xcopy`](@ref), [`xcopy_override`](@ref)

 # Example

 ```julia
 struct Foo
    value::Int
 end

 @xcopy Foo
 @xcopyoverride Foo :value tpl.value / 2
 ```
 """
macro xcopy_override(T, S, expr)
    if !(S isa QuoteNode && S.value isa Symbol) && !(S isa Symbol)
        error("Invalid property: must be either a symbol or a quote node")
    end
    if S isa Symbol
        S = QuoteNode(S)
    end
    esc(:(ExtraFun.xcopy_override(tpl::$T, ::ExtraFun.FieldCopyOverride{$S}) = $expr))
end

"""`xcopy_construct(tpl, args...; kwargs...)`

 Constructs the actual copy. Useful for running additional initialization code or to add additional constructor arguments.

 # Example

 ```julia
 struct Foo{T<:Real}
    value::T
    deduced::Float64
 end
 Foo(value::Real) = Foo(value, value / 2)

 xcopy_construct(tpl::Foo, args...) = Foo(args[1], args[2]/2)
 xcopy(Foo(24)) # == Foo(24, 12.0)
 ```
 """
xcopy_construct(tpl::T, args...; kwargs...) where T = T(args...; kwargs...)

max_stars_repo_name: Kirusifix/ExtraFun

max_issues_repo_name: Kirusifix/ExtraFun

max_forks_repo_name: Kirusifix/ExtraFun


text:
export invertinvoke, invertir, invertci

"""
`invertinvoke(f, t::Type{<:Tuple}, z, θ)`

Parametric inverse application of method `f(::t...)` to `args` with parameters `θ`
# Inputs:
`f` - function to invert
`t` - Tuple of types which determines method of `f` to invert
`z` - Input to inverse method
`θ` - Parameter values

# Returns
`(a, b, c, ...)` - tuple of values

```
f(x, y, z) = x * y + z
x, y, z = invertinvoke(f, Tuple{Float64, Float64, Float64}, 2.3, rand(3))
@assert f(x, y, z) == 2.3
```
"""
invertinvoke(f, t::Type{T}, z, θ) where {T<:Tuple} =
  choose(f, t, places_from_argtypes(T), ZΘ, z, θ)

"""
`invertir(f, t::Type{T})`

Produce inverse ast for method `f(::T)`
"""
invertir(f, t::Type{T}) where {T<:Tuple} =
  reorientir(typeof(f), t, ZΘ, places_from_argtypes(T))

"Produce inverse ::CodeInfo"
invertci(f, t::Type{T}) where {T<:Tuple} =
  reorient(typeof(f), t, ZΘ, places_from_argtypes(T))

@post length(T.parameters) == length(ret) "Output should be the right length"
@post all(map(isa, ret, T.parameters)) "Each element of ret is correct type"
@post Base.invoke(f, t, ret...) == z "Is valid inverse element"

invertir(f, t::Type{T}) where {T} = reorientir(f, t, places_from_argtypes(T), ZΘ())


text:
module TestConstructors
    using Base.Test
    using DataArrays

    #
    # NA's
    #

    @assert isna(NAtype())
    @assert isna(NA)

    #
    # DataVector's
    #

    dv = DataArray([1, 2, 3], falses(3))
    @assert isequal(dv.data, [1, 2, 3])
    @assert isequal(dv.na, falses(3))
    @assert isequal(dv, DataArray([1, 2, 3], [false, false, false]))
    @assert isequal(dv, DataArray([1, 2, 3]))

    dv = convert(DataArray, trues(3))
    @assert isequal(dv.data, [true, true, true])
    @assert isequal(dv.na, falses(3))
    @assert isequal(dv, convert(DataArray, trues(3)))

    dv = DataArray([1, 2, 3], falses(3))
    @assert isequal(dv, convert(DataArray, 1:3))

    dv = DataArray(Int, 3)
    @assert isequal(eltype(dv), Int)
    @assert isequal(dv.na, trues(3))

    dv = convert(DataArray, zeros(3))
    @assert isequal(dv, convert(DataArray, zeros(3)))

    dv = convert(DataArray, zeros(Int, 3))
    @assert isequal(dv, convert(DataArray, zeros(Int, 3)))

    dv = convert(DataArray, ones(3))
    @assert isequal(dv, convert(DataArray, ones(3)))

    dv = convert(DataArray, ones(Int, 3))
    @assert isequal(dv, convert(DataArray, ones(Int, 3)))

    dv = convert(DataArray, falses(3))
    @assert isequal(dv, convert(DataArray, falses(3)))

    dv = convert(DataArray, trues(3))
    @assert isequal(dv, convert(DataArray, trues(3)))

    #
    # PooledDataArray's
    #

    pdv = PooledDataArray([1, 2, 3], falses(3))
    @assert all(pdv .== [1, 2, 3])
    @assert all(isna(pdv) .== falses(3))

    @assert isequal(pdv, PooledDataArray([1, 2, 3], [false, false, false]))
    @assert isequal(pdv, PooledDataArray([1, 2, 3]))

    pdv = convert(PooledDataArray, trues(3))
    @assert all(pdv .== [true, true, true])
    @assert all(isna(pdv) .== falses(3))
    @assert isequal(pdv, convert(PooledDataArray, trues(3)))

    pdv = PooledDataArray([1, 2, 3], falses(3))
    @assert isequal(pdv, convert(PooledDataArray, 1:3))
    @assert isequal(pdv, convert(PooledDataArray, PooledDataArray([1, 2, 3])))

    pdv = PooledDataArray(Int, 3)
    @assert isequal(eltype(pdv), Int)
    @assert all(isna(pdv) .== trues(3))

    pdv = convert(PooledDataArray, zeros(3))
    @assert isequal(pdv, convert(PooledDataArray, zeros(3)))

    pdv = convert(PooledDataArray, zeros(Int, 3))
    @assert isequal(pdv, convert(PooledDataArray, zeros(Int, 3)))

    pdv = convert(PooledDataArray, ones(3))
    @assert isequal(pdv, convert(PooledDataArray, ones(3)))

    pdv = convert(PooledDataArray, ones(Int, 3))
    @assert isequal(pdv, convert(PooledDataArray, ones(Int, 3)))

    pdv = convert(PooledDataArray, falses(3))
    @assert isequal(pdv, convert(PooledDataArray, falses(3)))

    pdv = convert(PooledDataArray, trues(3))
    @assert isequal(pdv, convert(PooledDataArray, trues(3)))

    #
    # DataMatrix
    #

    dm = DataArray([1 2; 3 4], falses(2, 2))
    @assert isequal(dm.data, [1 2; 3 4])
    @assert isequal(dm.na, falses(2, 2))

    @assert isequal(dm, DataArray([1 2; 3 4], [false false; false false]))
    @assert isequal(dm, DataArray([1 2; 3 4]))

    dm = convert(DataArray, trues(2, 2))
    @assert isequal(dm.data, trues(2, 2))
    @assert isequal(dm.na, falses(2, 2))

    @assert isequal(dm, convert(DataArray, trues(2, 2)))

    dm = DataArray(Int, 2, 2)
    @assert isequal(eltype(dm), Int)
    @assert isequal(dm.na, trues(2, 2))

    convert(DataArray, zeros(2, 2))

    convert(DataArray, zeros(Int, 2, 2))

    convert(DataArray, ones(2, 2))
    convert(DataArray, ones(Int, 2, 2))

    convert(DataArray, falses(2, 2))
    convert(DataArray, trues(2, 2))

    convert(DataArray, eye(3, 2))
    convert(DataArray, eye(2))
    convert(DataArray, diagm(Float64[pi, pi]))
end


text:
include("./Functions.jl")

using Primes
using Functions

answer = 0
for j = 1:9
  for i = 1:(factorial(j))
    k = kthpermutation(i,j)
    if isprime(k)
      if k > answer
        answer = k
      end
    end
  end
end

println(answer)

max_stars_repo_name: vmchale/julia-problem-solving

max_issues_repo_name: vmchale/julia-problem-solving

max_forks_repo_name: vmchale/julia-problem-solving


text:
#===============================================================================

    shared_functions.jl

    functions used by multiple analytic methods

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

import DataStructures.OrderedDict
import JSON.json
import Memento: debug, warn, error, Logger

## used in amplification.jl
## used in shared_functions.jl
"""

    curry(f)(x...)

Currying function.

# Example:
```julia-repl
julia> (1,2,3) |> QpcrAnalysis.curry(map)(x -> 2x)
(2, 4, 6)
```
"""
@inline curry(f) = x -> (y...) -> f(x, y...)

## used in amplification.jl
## used in melting_curve.jl
## used in optical_test_dual_channel.jl
## used in shared_functions.jl
"""
    mold(f)(x...)

Curried `map` function.

# Example:
```julia-repl
julia> (1,2,3) |> QpcrAnalysis.mold(x -> 2x)
(2, 4, 6)
```
"""
mold   = curry(map)         ## mold(f)   = x -> map(f, x)

## used in deconvolute.jl
## used in amplification.jl
## used in melting_curve.jl
## used in shared_functions.jl
"""
    sift(f)(x...)

Curried `filter` function.

# Example:
```julia-repl
julia> 1:10 |> QpcrAnalysis.sift(isodd) |> Tuple
(1, 3, 5, 7, 9)
```
"""
sift   = curry(filter)      ## sift(f)   = x -> filter(f, x)

## used in melting_curve.jl
## used in shared_functions.jl
"""
    cast(f)(x...)

Curried `broadcast` function.

# Example:
```julia-repl
julia> Dict(:a=>0:2,:b=>[0,4,8]) |> values |> QpcrAnalysis.cast(mean) |> Tuple
(0.0, 2.5, 5.0)
```
"""
cast   = curry(broadcast)   ## cast(f)   = x -> broadcast(f, x)

## used in calibration.jl
## used in shared_functions.jl
"""
    from(x)(y...)

Curried `range` function.

# Example:
```julia-repl
julia> 10 |> QpcrAnalysis.from(1) |> collect |> Tuple
(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)
```
"""
from   = curry(range)       ## from(b)   = e -> range(b, e)

## used in amplification.jl
## used in allelic_discrimination.jl
"""
    gather(f)(xs...)

Curried `reduce` function.

# Example:
```julia-repl
julia> 1:5 |> QpcrAnalysis.gather(*)
120
```
"""
gather = curry(reduce)      ## gather(f) = x -> reduce(f, x)

## used in amplification.jl
"""
    bless(T)(xs...)

Curried `convert` function.

# Example:
```julia-repl
julia> -1 |> QpcrAnalysis.bless(Complex) |> sqrt
0.0 + 1.0im
```
"""
bless  = curry(convert)     ## bless(T)  = x -> convert(T, x)

## used in calibration.jl
## used in CalibrationInput.jl
"""
    splat(f)(x)

Curried splatting function. Iterate over the collection in the second argument,
passing the elements individually as arguments to the first.

# Examples:
```julia-repl
julia> [3,4] |> splat(//)
3//4

julia> (1,2,3) |> splat(hcat)
1×3 Array{Int64,2}:
 1, 2, 3
```
"""
splat(f) = x -> f(x...)

## used in calibration.jl
"""
    tie(n)(x...)

Curried `cat` function.

# Example:
```julia-repl
julia> reduce(tie(3),[zeros(1,3,1),ones(1,3,1)])
1×3×2 Array{Float64,3}:
[:, :, 1] =
 0.0  0.0  0.0

[:, :, 2] =
 1.0  1.0  1.0
```
"""
tie    = curry(cat)         ## tie(n)    = x -> cat(n, x...)

## used in shared_functions.jl
"""
    tail(f)(x...)

Tail-currying function.

# Example:
```julia-repl
julia> 1:10 |> QpcrAnalysis.tail(filter) |> QpcrAnalysis.tail(map)([iseven,isodd])
2-element Array{Array{Int64,1},1}:
 [2, 4, 6, 8, 10]
 [1, 3, 5, 7, 9]
```
"""
@inline tail(f)  = (x...) -> y -> f(y, x...)

## used in amplification.jl
## used in melting_curve.jl
## thermal_consistency.jl
"""
    field(name ::Symbol)(value)

Tail-curried `getfield` function.

# Example:
```julia-repl
julia> exp(1im * pi) |> QpcrAnalysis.field(:re)
-1.0
```
"""
field = tail(getfield)  ## field(f...) = x -> getfield(x, f...)

## used in amplification.jl
## used in mc_analysis.jl
## used in mc_peak_analysis.jl
## used in shared_functions.jl
"""
    index(inds...)(A)

Tail-curried `getindex` function."

# Example:
```julia-repl
julia> [1:3 4:6] |> QpcrAnalysis.index(1,2)
4
```
"""
index = tail(getindex)  ## index(i...) = x -> getindex(x, i...)

## used in CalibrationData.jl
## used in amp_analysis.jl
## used in mc_analysis.jl
"""
    morph(dims...)(A)

Tail-curried `reshape` function.

# Example:
```julia-repl
julia> [1:3 4:6] |> QpcrAnalysis.reshape(2,3)
2×3 Array{Int64,2}:
 1  3  5
 2  4  6
```
"""
morph = tail(reshape)   ## morph(d...) = x -> reshape(x, d...)

## used in allelic_discrimination.jl
## used in mc_peak_analysis.jl
"""
    furnish(dims...)(x)

Tail-curried `fill` function.

# Example:
```julia-repl
julia> NaN |> QpcrAnalysis.furnish(2,3)
2×3 Array{Int64,2}:
 NaN  NaN  NaN
 NaN  NaN  NaN
```
"""
furnish = tail(fill)    ## furnish(d...) = x -> fill(x, d...)

## used in amplification.jl
## used in allelic_discrimination.jl
## used in melting_curve.jl
"""
    moose(f, op)(itr)

Tail-curried `mapreduce` function.

# Example:
```julia-repl
julia> [1:5;] |> QpcrAnalysis.moose(hcat) do x; x^2 end |> Tuple
(1, 4, 9, 16, 25)
```
"""
@inline moose(f ::Function, op ::Function) = itr -> mapreduce(f, op, itr)

## reporter functions
@inline out(out_format ::OutputFormat) =
    output -> (out_format == json) ? json(output) : output
@inline report(digits ::Integer, x) = round.(x, digits)
const JSON_DIGITS = Int_T(6) ## number of decimal points for floats in JSON output

## used in amplification.jl
## used in mc_peak_analysis.jl
## used in thermal_consistency.jl
"""
    roundoff(digits ::Int_T)(x)

Curried reporter function."

# Example:
```julia-repl
julia> pi |> QpcrAnalysis.roundoff(2)
3.14
```
"""
@inline roundoff(digits ::Int_T) = cast(curry(report))(digits)

@inline get_ordered_keys(dict ::Dict) =
    dict |> keys |> collect |> sort

@inline get_ordered_keys(ordered_dict ::OrderedDict) =
    ordered_dict |> keys |> collect

## used in deconvolution.jl
## used in melting_curve.jl
## used in mc_peak_analysis.jl
## used in optical_test_dual_channel.jl
"""
    sweep(f)(op)(x)

Sweep out a summary function from a data collection `x`.

# Examples:
```julia-repl
julia> (85,100,115) |> QpcrAnalysis.sweep(minimum)(-) |> QpcrAnalysis.sweep(maximum)(/)
(0.0, 0.5, 1.0)

julia> (85,100,115) |> QpcrAnalysis.sweep(mean)(-) |> QpcrAnalysis.sweep(std)(/)
(-1.0, 0.0, 1.0)
```
"""
@inline sweep(summary_func) =
    sweep_func -> x -> broadcast(sweep_func, x, summary_func(x))

## used in enums.jl
"""
    x |> fan(fs)

Apply a collection of functions to a single argument.

# Example:
```julia-repl
julia> [0:10;] |> QpcrAnalysis.fan([minimum, median, maximum]) |> Tuple
(0, 5.0, 10)
"""
fan = fs -> x -> map(f -> f(x), fs)

## used in macros.jl
## used in Field.jl
## used in amp_analysis.jl
## used in mc_peak_analysis.jl
their(f) = x -> map(field(f), x) ## = mold(field(f))

## used in amplification.jl
## used in melting_curve.jl
## used in shared_functions.jl
# thing(x) = !(x === nothing)
thing = (!isequal)(nothing)

## used in amplification.jl
## used in melting_curve.jl
## used in optical_calibration.jl
## used in thermal_consistency.jl
"Handle exceptions by returning an `OrderedDict` containing the error message."
function fail(
    logger      ::Logger,
    err         ::Exception;
    bt          ::Bool = false ## backtrace?
)
	debug(logger, "at fail()")

	if (get(ENV, "JULIA_ENV", nothing) == PRODUCTION_MODE)
		bt = false
	end

    err_msg = sprint(showerror, err)
	try
        warn(logger, err_msg)
        if bt
			st = stacktrace(catch_backtrace())
            stl = collect(enumerate(IndexStyle(st), st))
            warn(logger, "error thrown in " *
                 string(st[1]) * "\nStacktrace:\n" *
                    join(
                        map(stl[1:end]) do tup
                            index, ptr = tup
                            " [$index] $ptr\n"
                        end))
        end ## if bt
    catch e
		warn(logger, "fail catch error: " * sprint(showerror, e))
    end ## try
    OrderedDict(
        :valid => false,
        :error => err_msg)
end ## fail()

## used in amplification.jl
## used in melting_curve.jl
"Finite differencing function. Three methods are implemented: `forward`, `backward`,
and `central`. The default central difference method provides the best approximation
of the derivative for twice-differentiable functions."
function finite_diff(
    X       ::AbstractVector,
    Y       ::AbstractVector; ## X and Y must be of same length
    nu      ::Integer = 1, ## order of derivative
    method  ::FiniteDiffMethod = central
)
    debug(logger, "at finite_diff()")
    dlen = length(X)
    if dlen != length(Y)
        throw(DimensionError, "X and Y must be of same length")
    end ## if
    (dlen == 1) && return zeros(1)
    if (nu == 1)
        (range1, range2) =
            if      (method == central)  tuple(3:dlen+2, 1:dlen)
            elseif  (method == forward)  tuple(3:dlen+2, 1:dlen+1)
            elseif  (method == backward) tuple(2:dlen+1, 1:dlen)
            else
                throw(ArgmentError, "method \"$method\" not recognized")
            end ## if
        (X_p2, Y_p2) = map((X, Y)) do ori
            vcat(
                ori[2] * 2 - ori[1],
                ori,
                ori[dlen-1] * 2 - ori[dlen])
            end ## do ori
        return (Y_p2[range1] .- Y_p2[range2]) ./ (X_p2[range1] .- X_p2[range2])
    end ## nu == 1
    return finite_diff(
        X,
        finite_diff(X, Y; nu = nu - 1, method = method),
        nu = 1;
        method = method)
end ## finite_diff()

## used in standard_curve.jl
"Transform a real number to scientific notation."
function scinot(
    x               ::Real,
    num_sig_digits  ::Integer = 3;
    log_base        ::Integer = 10
)
    isnan(x) && return (NaN_T, NaN_T)
    (x == 0) && return (0, 0)
    _exponent = log(log_base, abs(x)) |> floor
    _mantissa = round(x / log_base ^ _exponent, num_sig_digits)
    return (_mantissa, Int_T(_exponent))
end

## used in calibration.jl
## used in CalibrationData.jl
"Derive the number of wells from the collection of data vectors supplied to the function."
count_wells(fluos ::Array) =
    fluos |> mold(length) |> maximum
count_wells(dict ::Associative) =
    dict |> values |> mold(count_wells ∘ index(FLUORESCENCE_VALUE_KEY)) |> maximum
import Base.length
length(::Void) = 0

## deprecated in favour of Memento.info()
## still used in test functions
"Print function governed by Boolean output flag. Deprecated in favor of `Memento.info`."
function print_v(
    print_func ::Function,
    verbose ::Bool,
    args...;
    kwargs...
)
    if verbose
        print_func(args...; kwargs...)
    end
    return nothing
end

## unused function
## parse AbstractFloat on BBB
# function parse_af{T<:AbstractFloat}( ::Type{T}, strval ::String)
#     str_parts = split(strval, '.')
#     float_parts = map(str_part -> Base.parse(Int32, str_part), str_parts)
#     return float_parts[1] + float_parts[2] / 10^length(str_parts[2])
# end

## unused function
## repeat n times: take the output of a function
## and use it as the input for the same function
# function redo(
#     func ::Function,
#     input,
#     times ::Integer,
#     extra_args...;
#     kwargs...
# )
#     output = input
#     while times > 0
#         output = func(output, extra_args...; kwargs...)
#         times -= 1
#     end
#     return output
# end

## unused functions
# inc_index(i ::Integer, len ::Integer) = (i >= len) ? len : i + 1
# dec_index(i ::Integer) = (i <= 1) ? 1 : i - 1

## unused function
#
## reshape a layered vector into a multi-dimension array
## where outer layer is converted to higher dimension
## and each element has `num_layers_left` layers left
## (e.g. each element is atomic / not an array when `num_layers_lift == 0`,
## a vector of atomic elements when `num_layers_lift == 1`,
## vector of vector of atomic elements when `num_layers_lift == 2`).
# function reshape_lv(
#     layered_vector ::AbstractVector,
#     num_layers_left ::Integer=0
# )
#     md_array = copy(layered_vector) ## safe in case `eltype(layered_vector) <: AbstractArray`
#     while redo(eltype, md_array, num_layers_left + 1) <: AbstractArray
#         md_array = reshape(
#             cat(2, md_array...),
#             length(md_array[1]),
#             size(md_array)...)
#     end
#     return md_array
# end

## legacy function
## deprecated to remove MySql dependency
#
## function: check whether a value different from `calib_info_AIR`
## is passed onto `calib_info`. if not, use `exp_id` to find calibration experiment
## in MySQL database and assumes water "step_id"=2, signal "step_id"=4,
## using FAM to calibrate all the channels.
# function ensure_ci(
#
#     ## remove MySql dependency
#     #
#     # db_conn ::MySQL.MySQLHandle,
#     # calib_info ::Union{Integer,OrderedDict}=calib_info_AIR,
#
#     ## new >>
#     calib_data  ::OrderedDict{String,Any},
#     ## << new
#
#     ## use calibration data from experiment `calib_info_AIR` by default
#     exp_id      ::Integer ##=calib_info_AIR
# )
#     if isa(calib_info, Integer)
#
#         if calib_info == calib_info_AIR
#             calib_id = MySQL.mysql_execute(
#                 db_conn,
#                 "SELECT calibration_id FROM experiments WHERE id=$exp_id"
#             )[1][:calibration_id][1]
#         else
#             calib_id = calib_info
#         end
#
#         step_qry = "SELECT step_id FROM fluorescence_data WHERE experiment_id=$calib_id"
#         step_ids = sort(unique(MySQL.mysql_execute(db_conn, step_qry)[1][:step_id]))
#
#         calib_info = OrderedDict(
#             "water" => OrderedDict(
#                 "calibration_id" => calib_id,
#                 "step_id" => step_ids[1]))
#
#         for i in 2:(length(step_ids))
#             calib_info["channel_$(i-1)"] = OrderedDict(
#                 "calibration_id" => calib_id,
#                 "step_id" => step_ids[i])
#         end ## for
#
#         channel_qry = "SELECT channel FROM fluorescence_data WHERE experiment_id=$calib_id"
#         channels = sort(unique(MySQL.mysql_execute(db_conn, channel_qry)[1][:channel]))
#
#         for channel in channels
#             channel_key = "channel_$channel"
#             if !(channel_key in keys(calib_info))
#                 calib_info[channel_key] = OrderedDict(
#                     "calibration_id" => calib_id,
#                     "step_id" => step_ids[2])
#             end ## if
#         end ## for
#     end ## if isa(calib_info, Integer)
#     return calib_info
# end ## ensure_ci

## legacy function
## deprecated to remove MySql dependency
#
# function get_mysql_data_well(
#     well_nums ::AbstractVector, # must be sorted in ascending order
#     qry_2b ::String, # must select "well_num" column
#     db_conn ::MySQL.MySQLHandle,
#     verbose ::Bool,
# )
#
#     well_nums_str = join(well_nums, ',')
#     print_v(println, verbose, "well_nums: $well_nums_str")
#     well_constraint = (well_nums_str == "") ? "" : "AND well_num in ($well_nums_str)"
#     qry = replace(qry_2b, "well_constraint", well_constraint)
#     found_well_namedtuple = MySQL.mysql_execute(db_conn, qry)[1]
#     found_well_nums = sort(unique(found_well_namedtuple[:well_num]))
#     return (found_well_namedtuple, found_well_nums)
# end

## functions
## moved to MySQLforQpcrAnalysis.jl: get_mysql_data_well

## construct DataFrame from dictionary key and value vectors
## `dict_keys` need to be a vector of strings
## to construct DataFrame column indices correctly
# function dictvec2df(
#     dict_keys           ::AbstractVector,
#     dict_values         ::AbstractVector
# )
#     df = DataFrame()
#     for dict_key in dict_keys
#         df[Symbol(dict_key)] = map(index(dict_key), dict_values)
#     end
#     return df
# end

## simple macros
## deprecated in favour of:
## `true  && (Expr)`
## `false || (Expr)`
# macro when(predicate, conditional)
#     return :(if ($predicate)
#         ($conditional)
#     end)
# end
# macro unless(predicate, conditional)
#     return :(if !($predicate)
#         ($conditional)
#     end)
# end

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
### road_utils.jl

import JSON
import Base.Iterators as It
using ResumableFunctions

edge_pairs(G) = It.map(e -> (e.src, e.dst), edges(G))

edge_prop_triples(G) = It.map(e -> (e.src, e.dst, props(G, e.src, e.dst)), edges(G))

"Convert dictionary keys from strings to symbols"
symbolify(d::Dict) = Dict((Symbol(k) => v) for (k, v) ∈ pairs(d))

function read_graph(filename, bbox_ll = (-Inf, -Inf), bbox_ur = (Inf, Inf))::MetaDiGraph
    data = open(filename, "r") do io
		JSON.parse(io)
	end

	G = MetaDiGraph()

    # filter by latitude and longitude
	in_bbox(v_data) = (bbox_ll[1] <= v_data["latitude"] <= bbox_ur[1]) && (bbox_ll[2] <= v_data["longitude"] <= bbox_ur[2])

	for (i, v_data) ∈ enumerate(It.filter(in_bbox, data["nodes"]))
		@assert add_vertex!(G)
		@assert set_props!(G, i, symbolify(v_data))
	end
	set_indexing_prop!(G, :id)

	for e ∈ data["edges"]
		g_src_id = get(G[:id], e["from"], -1)
		g_dst_id = get(G[:id], e["to"], -1)
		if g_src_id < 0 || g_dst_id < 0
			continue
		end
        if g_src_id == g_dst_id
            continue
        end

		@assert add_edge!(G, g_src_id, g_dst_id)
		@assert set_props!(G, g_src_id, g_dst_id, symbolify(e))
	end

	G
end

# Gave up on trying to find the area of a cycle. I hate spherical geometry.

# "Distance between lat/long pair in metres ~ from StackOverflow"
# function geo_distance(lat1, lon1, lat2, lon2)
#     p = π/180
#     a = 0.5 - cos((lat2-lat1)*p)/2 + cos(lat1*p) * cos(lat2*p) * (1-cos((lon2-lon1)*p))/2
#     return 12742 * asin(sqrt(a)) * 1000
# end

# function triangle_area(lat_long_pairs)
#     @assert length(lat_long_pairs) == 3

#     d(u, v) = geo_distance(lat_long_pairs[u]..., lat_long_pairs[v]...)

#     # Heron's formula
#     a, b, c = d(1, 2), d(2, 3), d(3, 1)
#     @assert a >= 0 lat_long_pairs
#     @assert b >= 0 lat_long_pairs
#     @assert c >= 0 lat_long_pairs
#     s = (a + b + c) / 2
#     x = s * (s - a) * (s - b) * (s - c)
#     @assert x >= 0
#     return sqrt(x)
# end

# mean(x) = sum(x) / length(x)

# """
# Crude approximation to the area (m^2) enclosed by a cycle (does not take into account spherical geometry of Earth like `areaquad` in Matlab)

# Note: Assumes vertices not repeated
# """
# function cycle_area(G, cycle)
#     # Replace vertices by their latitude and longitude
#     geo_cycle = collect(It.map(i -> (get_prop(G, i, :latitude), get_prop(G, i, :longitude)), cycle))

#     # Probably not an "actual" midpoint in the case of weird cycles, but good enough
#     midpoint = (mean(first.(geo_cycle)), mean(last.(geo_cycle)))

#     idxmod = i -> mod1(i, length(cycle))

#     # TODO: Formula broken and does not work
#     sum(triangle_area([geo_cycle[idxmod(i)], geo_cycle[idxmod(i + 1)], midpoint]) for i in 1:length(cycle))
# end

"""
Returns a list of lists of directed basic cycles (grouped by their corresponding undirected cycle - so 0/1/2 per group) as tuples (vertex_list, area, length)

Use `It.flatten(basic_cycles(G))` if grouping is irrelevant
"""
function basic_cycles(G)
    # Convert `G` to an undirected graph (vertices *should* stay the same)
    # Can't convert directly because cycle_basis requires G <: AbstractSimpleGraph but MetaGraph, MetaDiGraph <: AbstractMetaGraph
    out = []

    UG = Graph(MetaGraph(G))
    for bs in cycle_basis(UG)
        if length(bs) <= 2
            println("Self-cycle: $(bs); Ignoring...")
            continue
        end

        idx = k -> bs[mod1(k, length(bs))]

        ℓ = (i, j) -> let
            ix, jx = idx(i), idx(j)
            if has_edge(G, ix, jx)
                get_prop(G, ix, jx, :length)
            else
                get_prop(G, jx, ix, :length)
            end
        end

        # area = cycle_area(G, bs)
        area = 0
        circ = sum(ℓ(i, i + 1) for i in 1:length(bs))

        cs = []
        if all(has_edge(G, idx(i), idx(i + 1)) for i in 1:length(bs))
            push!(cs, (bs, area, circ))
        end
        if all(has_edge(G, idx(i), idx(i - 1)) for i in 1:length(bs))
            push!(cs, (reverse(bs), area, circ))
        end
        isempty(cs) || push!(out, cs)
    end
    out
end

"""
Given an iterable of iterables, collect the `n`th element of each into a vector and yield these sequentially.
Necessary because collect(zip(its...)) will cause a StackOverflow due to how Julia does type inference
"""
function zipv(its)
end

@resumable function zipv(its)
    # Given an iterable of iterables, collect the `n`th element of each into a vector and yield these sequentially
    # Necessary because collect(zip(its...)) will cause a StackOverflow because of how Julia does type inference
    its = Iterators.Stateful.(its)
    while true
        nth = []
        for it ∈ its
            isempty(it) && return
            push!(nth, popfirst!(it))
        end
        @yield nth
    end
end

export edge_pairs, edge_prop_triples, symbolify, read_graph, basic_cycles, zipv
# export edge_pairs, edge_prop_triples, symbolify, read_graph, geo_distance, triangle_area, cycle_area, basic_cycles, zipv


text:
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Statistics, DataFrames, CSV
	using Random, Distributions, Plots
	using HypothesisTests
	Random.seed!(0)
end;

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing6.06"

# ╔═╡ 1aa41898-3afa-11eb-2f49-570853ce6435
begin
	mu1, sig1, n1 = 0, 2, 8
	mu2, sig2, n2 = 0, 30, 15
	dist1 = Normal(mu1, sig1)
	dist2 = Normal(mu2, sig2)

	N = 10^4
	tdArray = Array{Tuple{Float64,Float64}}(undef,N)

	df(s1,s2,n1,n2) =
		(s1^2/n1 + s2^2/n2)^2 / ( (s1^2/n1)^2/(n1-1) + (s2^2/n2)^2/(n2-1) )
end;

# ╔═╡ fd06acdc-3b1d-11eb-1d5f-03f5cf0f595a
begin
	for i in 1:N
		x1Data = rand(dist1, n1)
		x2Data = rand(dist2, n2)
		x1Bar,x2Bar = mean(x1Data),mean(x2Data)
		s1,s2 = std(x1Data),std(x2Data)
		tStat = (x1Bar - x2Bar) / sqrt(s1^2/n1 + s2^2/n2)
		tdArray[i] = (tStat , df(s1,s2,n1,n2))
	end
	sort!(tdArray, by = first)

	invVal(v,i) = quantile(TDist(v),i/(N+1))

	xCoords  = Array{Float64}(undef,N)
	yCoords1 = Array{Float64}(undef,N)
	yCoords2 = Array{Float64}(undef,N)

	for i in 1:N
		xCoords[i] = first(tdArray[i])
		yCoords1[i] = invVal(last(tdArray[i]), i)
		yCoords2[i] = invVal(n1+n2-2, i)
	end
end;

# ╔═╡ fd06e9ae-3b1d-11eb-1646-d9c440ed4b98
begin
	scatter(xCoords, yCoords1, c=:blue, label="Calculated v", msw=0)
	scatter!(xCoords, yCoords2, c=:red, label="Fixed v", msw=0)
	plot!([-10,10], [-10,10],
		 c=:black, lw=0.3, xlims=(-8,8), ylims=(-8,8), ratio=:equal, label="",
		 xlabel="Theoretical t-distribution quantiles",
		 ylabel="Simulated t-distribution quantiles", legend=:topleft)
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing6.06"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═1aa41898-3afa-11eb-2f49-570853ce6435
# ╠═fd06acdc-3b1d-11eb-1d5f-03f5cf0f595a
# ╠═fd06e9ae-3b1d-11eb-1646-d9c440ed4b98
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
# # Canvas

## load packages
using PlotPlants
#------------------------------------------------------------------------------

# ## Create canvas
# ### Default 1×1 figure (default)
_fig,_axes = create_canvas("1X1 figure");
_fig
#------------------------------------------------------------------------------

# ### M×N figure
_fig,_axes = create_canvas("M×N figure"; nrow=2, ncol=3);
_fig
#------------------------------------------------------------------------------

# ### Leave some panels empty
## this will omit the fifth panel (2nd row, 2nd column)
_fig,_axes = create_canvas("Omit panels"; nrow=2, ncol=3, axids=[1,2,3,4,6]);
_fig
#------------------------------------------------------------------------------

# ### Change figure size
_fig,_axes = create_canvas("Figure size"; nrow=2, ncol=3, figsize=(6,4));
_fig
#------------------------------------------------------------------------------

# ### Change figure resolution
_fig,_axes = create_canvas("Figure resolution"; nrow=2, ncol=3, dpi=50);
_fig
#------------------------------------------------------------------------------

# ## Panel titles
# ### Default
# Note that by default, `set_titles!` uses latex render, and thus the titles
#     are in a format like `\\textbf{(a)}`. In this way, the title can be
#     cutomized easily if one wants bold label and normal extra label. Also,
#     remember to enable latex render. If `usetex` option is `false`, titles
#     will be bold using the `fontweight` option in `matplotlib`.
_fig,_axes = create_canvas("Panel titles"; nrow=2, ncol=3);
set_titles!(_axes; usetex=false);
_fig
#------------------------------------------------------------------------------

# ### Add extra labels
_fig,_axes = create_canvas("Extra labels"; nrow=2, ncol=3);
set_titles!(_axes; labels=["extra label" for i in 1:6], usetex=false);
_fig
#------------------------------------------------------------------------------

# ### Remove label parenthesis
_fig,_axes = create_canvas("No parenthesis"; nrow=2, ncol=3);
set_titles!(_axes; paren=false, usetex=false);
_fig
#------------------------------------------------------------------------------

# ### Capital labels
_fig,_axes = create_canvas("Capital labels"; nrow=2, ncol=3);
set_titles!(_axes; capital=true, usetex=false);
_fig
#------------------------------------------------------------------------------

# ### Dot after labels
_fig,_axes = create_canvas("Dot after labels"; nrow=2, ncol=3);
set_titles!(_axes; paren=false, dotsep=true, usetex=false);
_fig
#------------------------------------------------------------------------------

# ### Set title fontsize
_fig,_axes = create_canvas("Title font size"; nrow=2, ncol=3);
set_titles!(_axes; fontsize=24, usetex=false);
_fig
#------------------------------------------------------------------------------

# ### Set title location
_fig,_axes = create_canvas("Title location"; nrow=2, ncol=3);
set_titles!(_axes; loc="left", usetex=false);
_fig
#------------------------------------------------------------------------------

# ## Axis labels
# ### Add X-axis labels
_fig,_axes = create_canvas("X axis labels"; nrow=2, ncol=3);
set_xlabels!(_axes, ["X label" for i in 1:6]);
_fig
#------------------------------------------------------------------------------

# ### Change X-axis label fontsize
_fig,_axes = create_canvas("X axis label fontsize"; nrow=2, ncol=3);
set_xlabels!(_axes, ["X label" for i in 1:6]; fontsize=12);
_fig
#------------------------------------------------------------------------------

# ### Add Y-axis labels
_fig,_axes = create_canvas("Y axis labels"; nrow=2, ncol=3);
set_ylabels!(_axes, ["Y label" for i in 1:6]);
_fig
#------------------------------------------------------------------------------

# ### Change Y-axis label fontsize
_fig,_axes = create_canvas("Y axis label fontsize"; nrow=2, ncol=3);
set_ylabels!(_axes, ["Y label" for i in 1:6]; fontsize=12);
_fig
#------------------------------------------------------------------------------

# ### Add both X and Y axis labels
_fig,_axes = create_canvas("X and Y axis"; nrow=2, ncol=3);
set_xylabels!(_axes, ["X label" for i in 1:6],
                     ["Y label" for i in 1:6]; fontsize=12);
_fig
#------------------------------------------------------------------------------

# ## Axis limits
# ### Set X-axis limits
_fig,_axes = create_canvas("X limits"; nrow=2, ncol=3);
set_xlims!(_axes, [[-1,1] for i in 1:6]);
_fig
#------------------------------------------------------------------------------

# ### Set Y-axis limits
_fig,_axes = create_canvas("Y limits"; nrow=2, ncol=3);
set_ylims!(_axes, [[-1,1] for i in 1:6]);
_fig
#------------------------------------------------------------------------------

# ### Set both X and Y axis limits
_fig,_axes = create_canvas("X and Y limits"; nrow=2, ncol=3);
set_xylims!(_axes, [[-1,1] for i in 1:6], [[-1,1] for i in 1:6]);
_fig
#------------------------------------------------------------------------------

# ## Axis ticks and ticklabels
# ### Set X-axis ticks and ticklabels
_fig,_axes = create_canvas("X ticks"; nrow=2, ncol=3);
set_xticks!(_axes, [[-1,0,1] for i in 1:6]);
set_xticklabels!(_axes, [["A","B","C"] for i in 1:6]);
_fig
#------------------------------------------------------------------------------

# ### Set Y-axis ticks and ticklabels
_fig,_axes = create_canvas("Y ticks"; nrow=2, ncol=3);
set_yticks!(_axes, [[-1,0,1] for i in 1:6]);
set_yticklabels!(_axes, [["D","E","F"] for i in 1:6]);
_fig
#------------------------------------------------------------------------------

# ### Set Both X and Y axis ticks and ticklabels
_fig,_axes = create_canvas("X and Y ticks"; nrow=2, ncol=3);
set_xyticks!(_axes, [[-1,0,1] for i in 1:6], [[-1,0,1] for i in 1:6]);
set_xyticklabels!(_axes, [["A","B","C"] for i in 1:6],
                         [["D","E","F"] for i in 1:6]);
_fig
#------------------------------------------------------------------------------

# ## Set tight layout
_fig,_axes = create_canvas("Tight layout"; nrow=2, ncol=3);
set_titles!(_axes; loc="left", usetex=false);
set_xylabels!(_axes, ["X label" for i in 1:6],
                     ["Y label" for i in 1:6]; fontsize=12);
set_xyticks!(_axes, [[-1,0,1] for i in 1:6], [[-1,0,1] for i in 1:6]);
set_xyticklabels!(_axes, [["A","B","C"] for i in 1:6],
                         [["D","E","F"] for i in 1:6]);
_fig.set_tight_layout(true);
_fig
#------------------------------------------------------------------------------


text:
using EggAnalysis2
using Base.Test

# write your own tests here
# @test 1 == 2

include("eegtests.jl")


text:
# db.jl

using Serialization

# Database item
struct DatabaseItem{T}
   time::Float64
   agent::String
   action::String
   key::String
   item::T
end # DatabaseItem

function read_from_file(file::String)
    io = open(file, "r")

    r = []
    while !eof(io)
        push!(r, deserialize(io))
    end

    close(io)

    return r
end

function add_to_file(file::String, data::Array{T, 1} where T <: Any)
    io = open(file, "a+")

    [serialize(io, r) for r in data]

    close(io)
end


text:
"""
    AbstractΔSizeStrategy

Abstract base type for strategies for how to change the size.

Only used as a transition until JuMP approach has been fully verified.
"""
abstract type AbstractΔSizeStrategy end

"""
    OnlyFor <: AbstractΔSizeStrategy

Change size only for the provided vertex.
"""
struct OnlyFor <: AbstractΔSizeStrategy end

"""
    AbstractJuMPΔSizeStrategy <: AbstractΔSizeStrategy

Abstract type for strategies to change or align the sizes of vertices using JuMP.
"""
abstract type AbstractJuMPΔSizeStrategy <: AbstractΔSizeStrategy end

"""
    ΔSizeFailError <: AbstractJuMPΔSizeStrategy
    ΔSizeFailError(msg::String)

Throws an `ErrorException` with message `msg`.
"""
struct ΔSizeFailError <: AbstractJuMPΔSizeStrategy
    msg::String
end

"""
    ΔSizeFailNoOp <: AbstractJuMPΔSizeStrategy
    ΔSizeFailNoOp()

Does not perform any action.
"""
struct ΔSizeFailNoOp <: AbstractJuMPΔSizeStrategy end

"""
    LogΔSizeExec <: AbstractJuMPΔSizeStrategy
    LogΔSizeExec(msg::String)
    LogΔSizeExec(level::Logging.LogLevel, msg::String)
    LogΔSizeExec(level::Logging.LogLevel, msg::String, andthen::AbstractJuMPΔSizeStrategy)

Logs `msg` at log level `level`, then executes `AbstractJuMPΔSizeStrategy andthen`.
"""
struct LogΔSizeExec <: AbstractJuMPΔSizeStrategy
    level::LogLevel
    msg::String
    andthen::AbstractJuMPΔSizeStrategy
end
LogΔSizeExec(msg::String) = LogΔSizeExec(Logging.Info, msg)
LogΔSizeExec(level::LogLevel, msg::String) = LogΔSizeExec(level, msg, ΔSizeFailNoOp())

"""
    DefaultJuMPΔSizeStrategy <: AbstractJuMPΔSizeStrategy

Default strategy intended to be used when adding some extra constraints or objectives to a model on top of the default.
"""
struct DefaultJuMPΔSizeStrategy <: AbstractJuMPΔSizeStrategy end

struct Exact end
struct Relaxed end

"""
    ΔNout{T} <: AbstractJuMPΔSizeStrategy
    ΔNout{T}(vertex::AbstractVertex, Δ::Integer)
    ΔNoutExact(vertex::AbstractVertex, Δ::Integer, fallback::AbstractJuMPΔSizeStrategy)
    ΔNoutRelaxed(vertex::AbstractVertex, Δ::Integer, fallback::AbstractJuMPΔSizeStrategy)

Strategy for changing nout of `vertex` by `Δ`, i.e new size is `nout(vertex) + Δ`.

If `T == Exact`, size change will be added as a constraint to the model which means that the operation will fail if it is not possible to change `nout(vertex)` by exactly `Δ`. If the operation fails, it will be retried with the `fallback` strategy (default `ΔNoutRelaxed`).

If `T == Relaxed`, size change will be added as an objective to the model which means that `nout(vertex)` might not change by exactly `Δ`. In addition, a constraint that `nout(vertex)` must change is also added.

If the operation fails, it will be retried with the `fallback` strategy (default `ΔNout{Relaxed}` if `T==Exact` and `ΔSizeFailError` if `T==Relaxed`).
"""
struct ΔNout{T} <: AbstractJuMPΔSizeStrategy
    vertex::AbstractVertex
    Δ::Integer
    fallback::AbstractJuMPΔSizeStrategy
end
ΔNoutExact(v::AbstractVertex, Δ::Integer) = ΔNout{Exact}(v, Δ, LogΔSizeExec(Logging.Warn, "Could not change nout of $v by $(Δ)! Relaxing constraints...", ΔNoutRelaxed(v, Δ)))
ΔNoutRelaxed(v::AbstractVertex, Δ::Integer) = ΔNout{Relaxed}(v, Δ, ΔSizeFailError("Could not change nout of $v by $(Δ)!!"))
fallback(s::ΔNout) = s.fallback

"""
    ΔNin{T} <: AbstractJuMPΔSizeStrategy
    ΔNin{T}(vertex::AbstractVertex, Δs::Vector{Maybe{Int}}, fallback::AbstractJuMPΔSizeStrategy)
    ΔNinExact(vertex::AbstractVertex, Δs::Vector{Maybe{Int}})
    ΔNinRelaxed(vertex::AbstractVertex, Δs::Vector{Maybe{Int}})

Strategy for changing nin of `vertex` by `Δs`, i.e new size is `nin(vertex) .+ Δs`. Note that `Δs` must have the same number of elements as `nin(vertex)`.

Use `missing` to indicate "no change required" as 0 will be interpreted as "must not change".

If `T == Exact`, size change will be added as a constraint to the model which means that the operation will fail if it is not possible to change `nin(vertex)` by exactly `Δs`.

If `T == Relaxed`, size change will be added as an objective to the model which means that `nin(vertex)` might not change by exactly `Δs`. In addition, a constraint that `nin(vertex)` must change is also added.

If the operation fails, it will be retried with the `fallback` strategy (default `ΔNin{Relaxed}` if `T==Exact` and `ΔSizeFailError` if `T==Relaxed`).
"""
struct ΔNin{T} <: AbstractJuMPΔSizeStrategy
    vertices::Vector{<:AbstractVertex}
    Δs::Vector{Int}
    fallback::AbstractJuMPΔSizeStrategy
    function ΔNin{T}(v, Δs, fallback) where T
        @assert size(Δs) == size(inputs(v)) "Must supply same number of Δs as v has inputs! Got $Δs for $v."
        inds = .!ismissing.(Δs)
        new(inputs(v)[inds], Δs[inds], fallback)
    end
end
ΔNinExact(v::AbstractVertex, Δs::Vector{<:Maybe{Int}}) = ΔNin{Exact}(v, Δs, LogΔSizeExec(Logging.Warn, "Could not change nin of $v by $(join(Δs, ", "))! Relaxing constraints...", ΔNinRelaxed(v, Δs)))
ΔNinExact(v::AbstractVertex, Δ::Integer) = ΔNin{Exact}(v, [Δ])
ΔNinRelaxed(v::AbstractVertex, Δs::Vector{<:Maybe{Int}}) = ΔNin{Relaxed}(v, Δs, ΔSizeFailError("Could not change nin of $vertex by $(join(Δs, ", "))!!"))
ΔNinRelaxed(v::AbstractVertex, Δ::Integer) = ΔNin{Relaxed}(v, [Δ])
fallback(s::ΔNin) = s.fallback

"""
    AlignNinToNout <: AbstractJuMPΔSizeStrategy
    AlignNinToNout(vstrat=DefaultJuMPΔSizeStrategy())

Adds variables and constraints for `nin(vi) == nout.(inputs(vi))`.

If it fails, the operation will be retried with the `fallback` strategy (default `ΔSizeFailError`).
"""
struct AlignNinToNout <: AbstractJuMPΔSizeStrategy
    nindict::Dict{AbstractVertex, Vector{JuMP.VariableRef}}
    vstrat::AbstractJuMPΔSizeStrategy
    fallback::AbstractJuMPΔSizeStrategy
end
AlignNinToNout(;vstrat=DefaultJuMPΔSizeStrategy(), fallback=ΔSizeFailError("Failed to align Nin to Nout!!")) = AlignNinToNout(vstrat, fallback)
AlignNinToNout(vstrat, fallback) = AlignNinToNout(Dict{AbstractVertex, JuMP.VariableRef}(), vstrat, fallback)
fallback(s::AlignNinToNout) = s.fallback

"""
    AlignNinToNoutVertices <: AbstractJuMPΔSizeStrategy
    AlignNinToNoutVertices(vin, vout, inds::Integer...;vstrat=AlignNinToNout(), fallback=ΔSizeFailError())
    AlignNinToNoutVertices(vin, vout, inds::AbstractArray{<:Integer},vstrat=AlignNinToNout(), fallback=ΔSizeFailError())

Same as [`AlignNinToNout`](@ref) with an additional constraint that `nin(s.vin)[s.ininds] == nout(s.vout)` where `s` is a `AlignNinToNoutVertices`.

Useful in the context of removing vertices and/or edges.

If it fails, the operation will be retried with the `fallback` strategy (default `ΔSizeFailError`).
"""
struct AlignNinToNoutVertices <: AbstractJuMPΔSizeStrategy
    vin::AbstractVertex
    vout::AbstractVertex
    ininds::AbstractArray{<:Integer}
    vstrat::AlignNinToNout
    fallback::AbstractJuMPΔSizeStrategy
end
AlignNinToNoutVertices(vin, vout, inds::Integer...;vstrat=AlignNinToNout(), fallback=failToAlign(vin, vout)) = AlignNinToNoutVertices(vin, vout, collect(inds), vstrat, fallback)
AlignNinToNoutVertices(vin, vout, inds::AbstractArray{<:Integer}, vstrat::AlignNinToNout, fallback::AbstractJuMPΔSizeStrategy=failToAlign(vin,vout)) = AlignNinToNoutVertices(vin, vout, inds, vstrat, fallback)
AlignNinToNoutVertices(vin, vout, inds::AbstractArray{<:Integer}, innerstrat::AbstractJuMPΔSizeStrategy, fallback=failToAlign(vin, vout)) = AlignNinToNoutVertices(vin, vout, inds, AlignNinToNout(vstrat=innerstrat), fallback)
fallback(s::AlignNinToNoutVertices) = s.fallback

failToAlign(vin, vout) = ΔSizeFailError("Could not align nout of $vin to nin of $(vout)!!")

#TODO: Ugh, this is too many abstraction layers for too little benefit. Refactor so
# all MutationVertex has state?
nin(v::InputSizeVertex) = []
nin(v::MutationVertex) = nin(v, op(v))
nin(v::AbstractVertex, op::MutationState) = nin(op)
nin(v::AbstractVertex, op::MutationOp) = nin(trait(v), v)
nin(t::DecoratingTrait, v::AbstractVertex) = nin(base(t), v)

# SizeTransparent does not need mutation state to keep track of sizes
nin(::SizeTransparent, v::AbstractVertex) = nout.(inputs(v))

nout(v::InputSizeVertex) = v.size
nout(v::MutationVertex) = nout(v, op(v))
nout(v::AbstractVertex, op::MutationState) = nout(op)
nout(v::AbstractVertex, op::MutationOp) = nout(trait(v), v)
nout(t::DecoratingTrait, v::AbstractVertex) = nout(base(t), v)

# SizeTransparent does not need mutation state to keep track of sizes
nout(t::SizeInvariant, v::AbstractVertex) = nin(v)[1]
nout(t::SizeStack, v::AbstractVertex) = sum(nin(v))

nout_org(v::AbstractVertex) = nout_org(trait(v), v)
nout_org(t::DecoratingTrait, v) = nout_org(base(t), v)
nout_org(::MutationSizeTrait, v::MutationVertex) = nout_org(op(v))
nout_org(::Immutable, v) = nout(v)

nin_org(v::InputSizeVertex) = []
nin_org(v::AbstractVertex) = nin_org(trait(v), v)
nin_org(t::DecoratingTrait, v) = nin_org(base(t), v)
nin_org(::MutationSizeTrait, v::MutationVertex) = nin_org(op(v))
nin_org(::Immutable, v) = nin(v)

"""
    minΔnoutfactor(v::AbstractVertex)

Returns the smallest `k` so that allowed changes to `nout` of `v` as well as `nin` of its outputs are `k * n` where `n` is an integer.
Returns `missing` if it is not possible to change `nout`.
"""
minΔnoutfactor(v::AbstractVertex) = missing
minΔnoutfactor(v::MutationVertex) = minΔnoutfactor(trait(v), v)
minΔnoutfactor(t::DecoratingTrait, v::AbstractVertex) = minΔnoutfactor(base(t), v)
minΔnoutfactor(::MutationTrait, v::AbstractVertex) = lcmsafe(vcat(minΔninfactor_only_for.(outputs(v))..., minΔnoutfactor_only_for(v)))
minΔnoutfactor(t::SizeTransparent, v::AbstractVertex) = minΔninfactor(t, v)

"""
    minΔnoutfactor(v::AbstractVertex, [s=VisitState{Int}()])

Returns the smallest `k` so that allowed changes to `nin` of `v` as well as `nout` of its inputs are `k * n` where `n` is an integer.
Returns `missing` if it is not possible to change `nin`.
"""
minΔninfactor(v::AbstractVertex) = lcmsafe(vcat(minΔnoutfactor_only_for.(inputs(v))..., minΔninfactor_only_for(v)))
minΔninfactor(v::MutationVertex) = minΔninfactor(trait(v), v)
minΔninfactor(t::DecoratingTrait, v::AbstractVertex) = minΔninfactor(base(t), v)
minΔninfactor(::MutationTrait, v::AbstractVertex) = lcmsafe(vcat(minΔnoutfactor_only_for.(inputs(v))..., minΔninfactor_only_for(v)))
minΔninfactor(::SizeTransparent, v::AbstractVertex) = lcmsafe([minΔnoutfactor_only_for(v), minΔninfactor_only_for(v)])

"""
    minΔnoutfactor_only_for(v::AbstractVertex, [s=VisitState{Int}()])

Returns the smallest `k` so that allowed changes to `nout` of `v` are `k * n` where `n` is an integer.
Returns `missing` if it is not possible to change `nout`.
"""
minΔnoutfactor_only_for(v::AbstractVertex, s=[]) = outvisit(v, s) ? 1 : minΔnoutfactor_only_for(base(v),s)
minΔnoutfactor_only_for(v::MutationVertex, s=[]) = outvisit(v, s) ? 1 : minΔnoutfactor_only_for(trait(v), v, s)
minΔnoutfactor_only_for(v::InputVertex, s=[]) = missing
minΔnoutfactor_only_for(v::CompVertex, s=[]) = minΔnoutfactor(v.computation)
minΔnoutfactor(f::Function) = 1 # TODO: Move to test as this does not make alot of sense
minΔnoutfactor_only_for(t::DecoratingTrait, v::AbstractVertex, s) = minΔnoutfactor_only_for(base(t), v, s)
minΔnoutfactor_only_for(::Immutable, v::AbstractVertex, s) = missing
minΔnoutfactor_only_for(::SizeAbsorb, v::AbstractVertex, s) = minΔnoutfactor_only_for(base(v),s)

minΔnoutfactor_only_for(::SizeInvariant, v::AbstractVertex, s) = lcmsafe(vcat(map(vi -> minΔnoutfactor_only_for(vi,s), inputs(v)), map(vo -> minΔninfactor_only_for(vo,s), outputs(v))))
function minΔnoutfactor_only_for(::SizeStack, v::AbstractVertex, s)
    absorbing = findterminating(v, inputs, outputs)

    # This is not strictly the minimum as using only one of the factors would work as well
    # However, this would create a bias as the same factor would be used all the time
    # Life is really hard sometimes :(

    # Count thingy is for duplicate outputs. Must be counted twice as it is impossible
    # to only change one of them, right?
    factors = [count(x->x==va,absorbing) * minΔnoutfactor_only_for(va,s) for va in unique(absorbing)]
    return lcmsafe(factors)
end

"""
    minΔninfactor_only_for(v::AbstractVertex)

Returns the smallest `k` so that allowed changes to `nin` of `v` are `k * n` where `n` is an integer.
Returns `missing` if it is not possible to change `nin`.
"""
minΔninfactor_only_for(v::AbstractVertex, s=[]) = invisit(v, s) ? 1 : minΔninfactor_only_for(base(v),s)
minΔninfactor_only_for(v::MutationVertex, s=[]) = invisit(v, s) ? 1 : minΔninfactor_only_for(trait(v), v, s)
minΔninfactor_only_for(v::InputVertex,s=[]) = missing
minΔninfactor_only_for(v::CompVertex,s=[]) = minΔninfactor(v.computation)
minΔninfactor(f::Function) = 1 # TODO: Move to test as this does not make alot of sense
minΔninfactor_only_for(t::DecoratingTrait, v::AbstractVertex, s) = minΔninfactor_only_for(base(t), v, s)
minΔninfactor_only_for(::Immutable, v::AbstractVertex, s) = missing
minΔninfactor_only_for(::SizeAbsorb, v::AbstractVertex, s) = minΔninfactor_only_for(base(v), s)
minΔninfactor_only_for(t::SizeInvariant, v::AbstractVertex, s) = minΔnoutfactor_only_for(t, v, s)
function minΔninfactor_only_for(::SizeStack, v::AbstractVertex, s)
    absorbing = findterminating(v, outputs, inputs)
    # This is not strictly the minimum as using only one of the factors would work as well
    # However, this would create a bias as the same factor would be used all the time
    # Life is really hard sometimes :(

    # Count thingy is for duplicate inputs. Must be counted twice as it is impossible
    # to only change one of them, right?
    factors = [count(x->x==va,absorbing) * minΔnoutfactor_only_for(va, s) for va in unique(absorbing)]
    return lcmsafe(factors)
end

#lcm which also checks for missing and arrays of undefined type
function lcmsafe(x)
    isempty(x) && return 1
    return any(ismissing.(x)) ? missing : lcm(Integer.(x))
end

invisit(v, visited) = hasvisited(v, Input(), visited)
outvisit(v, visited) = hasvisited(v, Output(), visited)

function hasvisited(v, d::Direction, visited)
    (v, d) in visited && return true
    push!(visited, (v, d))
    return false
end

"""
    findterminating(v::AbstractVertex, direction::Function, other::Function= v -> [], visited = [])

Return an array of all vertices which terminate size changes (i.e does not propagate them) seen through the given direction (typically inputs or outputs). A vertex will be present once for each unique path through which its seen.

The `other` direction may be specified and will be traversed if a SizeInvariant vertex is encountered.

Will return the given vertex if it is terminating.

# Examples
```julia-repl

julia> v1 = inputvertex("v1", 3);

julia> v2 = inputvertex("v2", 3);

julia> v3 = conc(v1,v2,v1,dims=1);

julia> name.(findterminating(v1, outputs, inputs))
1-element Array{String,1}:
 "v1"

julia> name.(findterminating(v3, outputs, inputs))
0-element Array{Any,1}

julia> name.(findterminating(v3, inputs, outputs))
3-element Array{String,1}:
 "v1"
 "v2"
 "v1"

 julia> v5 = v3 + inputvertex("v4", 9);

 julia>  # Note, + creates a SizeInvariant vertex and this causes its inputs to be seen through the output direction

 julia> name.(findterminating(v3, outputs, inputs))
 1-element Array{String,1}:
  "v4"
```
"""
function findterminating(v::AbstractVertex, direction::Function, other::Function=v->[], visited = Set{AbstractVertex}())
    v in visited && return []
    push!(visited, v)
    res = findterminating(trait(v), v, direction, other, visited)
    delete!(visited, v)
    return res
 end
findterminating(t::DecoratingTrait, v, d::Function, o::Function, visited) = findterminating(base(t), v, d, o, visited)
findterminating(::SizeAbsorb, v, d::Function, o::Function, visited) = [v]
findterminating(::Immutable, v, d::Function, o::Function, visited) = [v]

findterminating(::SizeStack, v, d::Function, o::Function, visited) = collectterminating(v, d, o, visited)
findterminating(::SizeInvariant, v, d::Function, o::Function, visited) = vcat(collectterminating(v, d, o, visited), collectterminating(v, o, d, visited))
collectterminating(v, d::Function, o::Function, visited) = mapfoldl(vf -> findterminating(vf, d, o, visited), vcat, d(v), init=[])

# Just another name for Δsize with the corresponding direction
Δnout(v::AbstractVertex, Δ::Integer) = Δsize(Output(), v, Δ)
Δnin(v::AbstractVertex, Δs::Maybe{<:Integer}...) = Δsize(Input(), v, Δs...)

"""
    Δsize(d::Direction, v::AbstractVertex, Δ...)

Change size of `v` by `Δ` in direction `d`.
"""
Δsize(d::Input, v, Δs::Maybe{T}...) where T <:Integer = Δsize(ΔNinExact(v, collect(Maybe{T}, Δs)), all_in_Δsize_graph(v, d))
Δsize(d::Output, v, Δ::T) where T <: Integer = Δsize(ΔNoutExact(v, Δ), all_in_Δsize_graph(v, d))

"""
    Δsize(s::AbstractΔSizeStrategy, vertices::AbstractArray{<:AbstractVertex})

Calculate new sizes for (potentially) all provided `vertices` using the strategy `s` and apply all changes.
"""
function Δsize(s::AbstractΔSizeStrategy, vertices::AbstractVector{<:AbstractVertex})
    execute, nins, nouts = newsizes(s, vertices)
    if execute
        Δsize(nins, nouts, vertices)
    end
end

"""
    Δsize(nouts::AbstractVector{<:Integer}, vertices::AbstractVector{<:AbstractVertex})

Set output size of `vertices[i]` to `nouts[i]` for all `i` in `1:length(vertices)`.
Set input size of all keys `vi` in `nins` to `nins[vi]`.
"""
function Δsize(nins::Dict, nouts::AbstractVector{<:Integer}, vertices::AbstractVector{<:AbstractVertex})
    Δnouts = nouts .- nout.(vertices)

    for (i, vi) in enumerate(vertices)
        ninΔs = get(() -> nin(vi), nins, vi) .- nin(vi)
        Δnin(OnlyFor(), vi, ninΔs...)
        Δnout(OnlyFor(), vi, Δnouts[i])
    end

    for (i, vi) in enumerate(vertices)
        ninΔs = get(() -> nin(vi), nins, vi) .- nin(vi)
        after_Δnin(vi, ninΔs...)
        after_Δnout(vi, Δnouts[i])
    end
end

function Δnin(s::OnlyFor, v, Δs::Maybe{<:Integer}...)
    any(skipmissing(Δs) .!= 0) || return
    Δnin(s, trait(v), v, Δs)
end
Δnin(s::AbstractΔSizeStrategy, t::DecoratingTrait, v, Δs) = Δnin(s, base(t), v, Δs)
function Δnin(s::AbstractΔSizeStrategy, t::SizeChangeLogger, v, Δs)

    @logmsg t.level "Change nin of $(infostr(t, v)) by $(join(compressed_string.(Δs), ", "))"
    Δnin(s, base(t), v, Δs)
end
Δnin(::OnlyFor, ::MutationSizeTrait, v, Δs) = Δnin(op(v), Δs...)

function Δnout(s::OnlyFor, v, Δ::Integer)
    Δ == 0 && return
    Δnout(s, trait(v), v, Δ)
end
Δnout(s::AbstractΔSizeStrategy, t::DecoratingTrait, v, Δ) = Δnout(s, base(t), v, Δ)
function Δnout(s::AbstractΔSizeStrategy, t::SizeChangeLogger, v, Δ)
    @logmsg t.level "Change nout of $(infostr(t, v)) by $(compressed_string(Δ))"
    Δnout(s, base(t), v, Δ)
end
Δnout(::OnlyFor, ::MutationSizeTrait, v, Δ) = Δnout(op(v), Δ)

after_Δnin(v, Δs...) = after_Δnin(trait(v), v, Δs)
after_Δnin(t::DecoratingTrait, v, Δs) = after_Δnin(base(t), v, Δs)
after_Δnin(t::SizeChangeValidation, v, Δs) = validate_Δnin(v, Δs, () -> after_Δnin(base(t), v, Δs))
function after_Δnin(t, v, Δs) end

after_Δnout(v, Δ) = after_Δnout(trait(v), v, Δ)
after_Δnout(t::DecoratingTrait, v, Δ) = after_Δnout(base(t), v, Δ)
after_Δnout(t::SizeChangeValidation, v, Δ) = validate_Δnout(v, Δ, () -> after_Δnout(base(t), v, Δ))
function after_Δnout(t, v, Δ) end

sizeΔ(Δ::Integer) = Δ
sizeΔ(Δ::AbstractArray) = length(Δ)
function validate_Δnin(v::AbstractVertex, Δ, Δfun, validvisit = true)

    # Yeah, this is checking more than one thing. Cba to have three different structs and methods for validation
    length(Δ) == length(inputs(v)) || throw(ArgumentError("Length of Δ must be equal to number of inputs for $(v)! length(Δ) = $(length(Δ)), length(inputs(v)) = $(length(inputs(v)))"))

    if validvisit
        # TODO base(v) makes this a bit weaker than I would have wanted. Right now it is only because testcases use smaller factors to trigger SizeStack to do unusual stuff
        Δninfactor = minΔninfactor_only_for(base(v))
        any(Δi -> sizeΔ(Δi) % Δninfactor != 0, skipmissing(Δ)) && throw(ArgumentError("Nin change of $Δ to $v is not an integer multiple of $(Δninfactor)!"))
    end

    Δfun()

    if validvisit
        nout.(inputs(v)) == nin(v) || throw(ArgumentError("Nin change of $Δ to $v did not result in expected size! Expected: $(nout.(inputs(v))), actual: $(nin(v))"))
    end
end

function validate_Δnout(v::AbstractVertex, Δ, Δfun, validvisit=true)

    if validvisit
        # TODO base(v) makes this a bit weaker than I would have wanted. Right now it is only because testcases use smaller factors to trigger SizeStack to do unusual stuff
        Δnoutfactor = minΔnoutfactor_only_for(base(v))
        sizeΔ(Δ) % Δnoutfactor != 0 && throw(ArgumentError("Nout change of $Δ to $v is not an integer multiple of $(Δnoutfactor)!"))
    end

    Δfun()

    if validvisit
        nin_of_outputs = unique(mapreduce(vi -> nin(vi)[inputs(vi) .== v], vcat, outputs(v), init=nout(v)))

        nin_of_outputs == [nout(v)] || throw(ArgumentError("Nout change of $Δ to $v resulted in size mismatch! Nin of outputs: $nin_of_outputs, nout of this: $([nout(v)])"))
    end
end

newsizes(s::ΔSizeFailError, vertices::AbstractVector{<:AbstractVertex}) = error(s.msg)
newsizes(s::ΔSizeFailNoOp, vertices::AbstractVector{<:AbstractVertex}) = false, Dict(vertices .=> nin.(vertices)), nout.(vertices)
function newsizes(s::LogΔSizeExec, vertices::AbstractVector{<:AbstractVertex})
    @logmsg s.level s.msg
    return newsizes(s.andthen, vertices)
end

"""
    newsizes(s::AbstractΔSizeStrategy, vertices::AbstractArray{<:AbstractVertex})

Return a vector of new outputs sizes for and a `Dict` of new input sizes for all provided `vertices` using the strategy `s`.

Result vector is index aligned with `vertices`.
Result `Dict` has a vector of input sizes for each element of `vertices` which has an input (i.e everything except input vertices).
"""
function newsizes(s::AbstractJuMPΔSizeStrategy, vertices::AbstractVector{<:AbstractVertex})

    model = sizemodel(s, vertices)

    noutvars = @variable(model, noutvars[i=1:length(vertices)], Int)

    noutdict = Dict(zip(vertices, noutvars))
    for v in vertices
        vertexconstraints!(v, s, (model=model, noutdict=noutdict))
    end

    sizeobjective!(s, model, noutvars, vertices)

    JuMP.optimize!(model)

    if accept(s, model)
        return true, ninsAndNouts(s, vertices, noutvars)...
    end
    return newsizes(fallback(s), vertices)
end

"""
    sizemodel(s::AbstractJuMPΔSizeStrategy, vertices)

Return a `JuMP.Model` for executing strategy `s` on `vertices`.
"""
sizemodel(s::AbstractJuMPΔSizeStrategy, vertices) = JuMP.Model(JuMP.optimizer_with_attributes(Cbc.Optimizer, "loglevel"=>0))

# Just a shortcut for broadcasting on dicts
getall(d::Dict, ks, deffun=() -> missing) = get.(deffun, [d], ks)

# First we dispatch on trait in order to filter out immutable vertices
vertexconstraints!(v::AbstractVertex, s, data) = vertexconstraints!(trait(v), v, s, data)
vertexconstraints!(t::DecoratingTrait, v, s, data) = vertexconstraints!(base(t), v, s,data)
function vertexconstraints!(::Immutable, v, s, data)
    @constraint(data.model, data.noutdict[v] == nout(v))
    @constraint(data.model, getall(data.noutdict, inputs(v)) .== nin(v))
    vertexconstraints!(s, v, data)
end
vertexconstraints!(::MutationTrait, v, s, data) = vertexconstraints!(s, v, data)

# This must be applied to immutable vertices as well
function vertexconstraints!(v::AbstractVertex, s::AlignNinToNout, data)
    vertexconstraints!(v, s.vstrat, data)
    # Code below secretly assumes vo is in data.noutdict (ninarr will be left with undef entries otherwise).
    for vo in filter(vo -> vo in keys(data.noutdict), outputs(v))
        ninvar = @variable(data.model, integer=true)
        @constraint(data.model, data.noutdict[v] == ninvar)

        ninarr = get!(() -> Vector{JuMP.VariableRef}(undef, length(inputs(vo))), s.nindict, vo)
        ninarr[inputs(vo) .== v] .= ninvar
    end
end

function vertexconstraints!(v::AbstractVertex, s::AlignNinToNoutVertices, data)
    # Any s.ininds which are mapped to vertices which are not yet added in s.vstrat.nindict[s.vout] will cause undefined reference below
    # There is a check to wait for them to be added, but if they are not in data.noutdict they will never be added, so we need to check for that too. Not 100% this can even happen, so I'm just waiting for this in itself to trigger a bug. Sigh...
    neededinds = filter(i -> i != nothing, indexin(keys(data.noutdict), inputs(s.vout)))

    condition() = s.vout in keys(s.vstrat.nindict) && all(i -> isassigned(s.vstrat.nindict[s.vout], i), neededinds)

    # Just to make sure we only do this once without having to store hasadded as a field in AlignNinToNoutVertices:
    #  - Check the condition before vertexconstraints! -> Only when condition changes from false to true do we add as we assume this is the first time the condition becomes fulfilled
    hasadded = condition()
    vertexconstraints!(v, s.vstrat, data)

    if !hasadded && condition()
        @constraint(data.model, data.noutdict[s.vin] .== s.vstrat.nindict[s.vout][s.ininds])
    end
end

"""
    vertexconstraints!(s::AbstractJuMPΔSizeStrategy, v, data)

Add constraints for `AbstractVertex v` using strategy `s`.

Extra info like the model and variables is provided in `data`.
"""
function vertexconstraints!(s::AbstractJuMPΔSizeStrategy, v, data)
    ninconstraint!(s, v, data)
    compconstraint!(s, v, (data..., vertex=v))
    sizeconstraint!(s, v, data)
end

"""
    sizeconstraint!(s::AbstractJuMPΔSizeStrategy, v, data)

Add size constraints for `AbstractVertex v` using strategy `s`.

Extra info like the model and variables is provided in `data`.
"""
sizeconstraint!(s::AbstractJuMPΔSizeStrategy, v, data) = @constraint(data.model, data.noutdict[v] >= 1)

function sizeconstraint!(s::ΔNout{Exact}, v, data)
    if v == s.vertex
        @constraint(data.model, data.noutdict[v] == nout(v) + s.Δ)
    else
        sizeconstraint!(DefaultJuMPΔSizeStrategy(), v, data)
    end
end

function sizeconstraint!(s::ΔNin{Exact}, v, data)
    if v in s.vertices
        Δ = s.Δs[s.vertices .== v][1]
        @constraint(data.model, data.noutdict[v] == nout(v) + Δ)
    else
        sizeconstraint!(DefaultJuMPΔSizeStrategy(), v, data)
    end
end

"""
    ninconstraint!(s, v, data)

Add input size constraints for `AbstractVertex v` using strategy `s`.

Extra info like the model and variables is provided in `data`.
"""
ninconstraint!(s, v, data) = ninconstraint!(s, trait(v), v, data)
ninconstraint!(s, t::DecoratingTrait, v, data) = ninconstraint!(s, base(t), v, data)
function ninconstraint!(s, ::MutationTrait, v, data) end
ninconstraint!(s, ::SizeStack, v, data) = @constraint(data.model, sum(getall(data.noutdict, inputs(v))) == data.noutdict[v])
ninconstraint!(s, ::SizeInvariant, v, data) = @constraint(data.model, getall(data.noutdict, unique(inputs(v))) .== data.noutdict[v])

"""
    ninconstraint!(s, v, data)

Add constraints on the computation (e.g. neural network layer) for `AbstractVertex v` using strategy `s`.

Extra info like the model and variables is provided in `data`.
"""
compconstraint!(s, v::AbstractVertex, data) = compconstraint!(s, base(v), data)
compconstraint!(s, v::CompVertex, data) = compconstraint!(s, v.computation, data)
function compconstraint!(s, v::InputVertex, data) end
function compconstraint!(s, f, data) end

"""
    sizeobjective!(s::AbstractJuMPΔSizeStrategy, model, noutvars, sizetargets)

Add the objective for `noutvars` using strategy `s`.
"""
sizeobjective!(s::AbstractJuMPΔSizeStrategy, model, noutvars, vertices) = @objective(model, Min, objective!(s, model, noutvars, vertices))

function objective!(s, model, noutvars, vertices)
    sizetargets = nout.(vertices)
    # L1 norm prevents change in vertices which does not need to change.
    # Max norm tries to spread out the change so no single vertex takes most of the change.
    return norm!(SumNorm(0.1 => L1NormLinear(), 0.8 => MaxNormLinear()), model, @expression(model, objective[i=1:length(noutvars)], noutvars[i] - sizetargets[i]), sizetargets)
end

objective!(s::ΔNout{Relaxed}, model, noutvars, vertices) = noutrelax!(model, [s.vertex], [s.Δ], noutvars, vertices)
objective!(s::ΔNin{Relaxed}, model, noutvars, vertices) = noutrelax!(model, s.vertices, s.Δs, noutvars, vertices)

function noutrelax!(model, vs, Δs, noutvars, vertices)
    inds = mapreduce(v -> vertices .== v, (i1,i2) -> i1 .| i2, vs)
    def_obj = objective!(DefaultJuMPΔSizeStrategy(), model, noutvars[.!inds], vertices[.!inds])
    sizetarget = nout.(vs) + Δs
    Δnout_obj = norm!(L1NormLinear(), model, @expression(model, noutvars[inds] .- sizetarget))
    # Force it to change as s.Δ might be too small
    # Trick from http://lpsolve.sourceforge.net/5.1/absolute.htm
    Δnout_const = @expression(model, noutvars[inds] - nout.(vs))
    B = @variable(model, [1:length(vs)], binary=true)
    M = 1e5
    ϵ = 1e-2 # abs(Δnout_const) must be larger than this
    @constraint(model, Δnout_const .+ M .* B .>= ϵ)
    @constraint(model, Δnout_const .+ M .* B .<= M .- ϵ)

    return @expression(model, def_obj + 1e6*sum(Δnout_obj))
end

"""
    accept(::AbstractJuMPΔSizeStrategy, model::JuMP.Model)

Return true of the solution for `model` is accepted using strategy `s`.
"""
accept(::AbstractJuMPΔSizeStrategy, model::JuMP.Model) = JuMP.termination_status(model) != MOI.INFEASIBLE && JuMP.primal_status(model) == MOI.FEASIBLE_POINT

function ninsAndNouts(::AbstractJuMPΔSizeStrategy, vs, noutvars)
    nouts = round.(Int, JuMP.value.(noutvars))
    mapnout(i::Integer) = nouts[i]
    mapnout(i::Nothing) = missing

    nins = Dict(vs .=> map(vi -> mapnout.(indexin(inputs(vi), vs)), vs))
    return nins, nouts
end

function ninsAndNouts(s::AlignNinToNout, vs, noutvars)
    nouts = round.(Int, JuMP.value.(noutvars))
    nins = Dict(key => round.(Int, JuMP.value.(value)) for (key, value) in s.nindict)
    return nins,nouts
end

ninsAndNouts(s::AlignNinToNoutVertices, vs, noutvars) = ninsAndNouts(s.vstrat, vs, noutvars)

# TODO: Remove since only used for debugging. If only it wasn't so bloody cumbersome to just list the constraints in a JuMP model....
nconstraints(model) = mapreduce(tt -> JuMP.num_constraints.(model,tt...), +,  filter(tt -> tt != (JuMP.VariableRef, MOI.Integer), JuMP.list_of_constraint_types(model)), init=0)

# TODO: Remove since only used for debugging. If only it wasn't so bloody cumbersome to just list the constraints in a JuMP model....
function list_constraints(model)
    for tt in filter(tt -> tt != (JuMP.VariableRef, MOI.Integer), JuMP.list_of_constraint_types(model))
        display(JuMP.all_constraints(model, tt...))
    end
end


text:
using Tau
using Base.Test

@test tau == 2*pi
@test float32(tau) == 2*float32(pi)
@test float64(float32(tau)) == float64(2*float32(pi))
@test big(tau) == 2(big(pi))
@test isa(tau, MathConst)

# degree-based trig functions
for T = (Float32,Float64)
    for x = -2:0.1:2
        @test_approx_eq_eps sintau(convert(T,x))::T convert(T,sin(tau*x)) eps(tau*convert(T,x))
        @test_approx_eq_eps costau(convert(T,x))::T convert(T,cos(tau*x)) eps(tau*convert(T,x))
    end
    for x = 0.0:1.0:4.0
        @test sintau(convert(T,x)) === zero(T)
        @test sintau(-convert(T,x)) === -zero(T)
    end
end

# check type stability
for T = (Float32,Float64,BigFloat)
    for f = (sintau,costau)
        @test Base.return_types(f,(T,)) == [T]
    end
end

# Taken from /julia/test/mod2pi.jl
@test_throws ErrorException modtau(int64(2)^60-1)

@test_approx_eq modtau(10)          mod(10, tau)
@test_approx_eq modtau(-10)         mod(-10, tau)
@test_approx_eq modtau(355)         3.1416227979431572
@test_approx_eq modtau(int32(355))  3.1416227979431572
@test_approx_eq modtau(355.0)       3.1416227979431572
@test_approx_eq modtau(355.0f0)     3.1416228f0
@test modtau(int64(2)^60) == modtau(2.0 ^ 60)


text:
#===============================================================================

    McOutput.jl

    containers for output from melting curve analysis

    Author: Tom Price
    Date:   July 2019

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

import StaticArrays: SVector, SArray

abstract type McOutput end

struct McLongOutput <: McOutput
    wells                       ::SVector{W,Symbol} where {W}
    channels                    ::SVector{C,Int_T} where {C}
    raw_data                    ::Array{<: Real,3}
    background_subtracted_data  ::Array{<: Real,3}
    k_deconv                    ::DeconvolutionMatrices
    deconvoluted_data           ::Array{Float_T,3}
    norm_data                   ::SArray{S,<: Real,3} where {S}
    norm_wells                  ::SVector{V,Symbol} where {V}
    calibrated_data             ::Array{Float_T,3}
    peak_output                 ::Array{McPeakOutput,2} ## dim1 is well and dim2 is channel
    # tf_bychwl                   ::OrderedDict{Int_T,Vector{OrderedDict{String,Vector{Float_T}}}}
end

## the following type not necessary because `output_dict` in mc_analysis()
## can be constructed directly from `mc_matrix`:
# abstract type McShortOutput <: McOutput end


text:
using Pkg
using InteractiveUtils
using CUDA
using Oceananigans.Architectures

function versioninfo_with_gpu()
    s = sprint(versioninfo)
    @hascuda begin
        gpu_name = CUDA.device() |> CUDA.name
        s = s * "  GPU: $gpu_name\n"
    end
    return s
end

function oceananigans_versioninfo()
    project = Pkg.project()

    # If Oceananigans is listed as a dependency in a Project.toml
    # (or in the base Julia env)
    if "Oceananigans" in keys(project.dependencies)
        uuid = project.dependencies["Oceananigans"]
        pkg_info = Pkg.dependencies()[uuid]
        s = "Oceananigans v$(pkg_info.version)"
        s *= isnothing(pkg_info.git_revision) ? "" : "#$(pkg_info.git_revision)"
        return s
    end

    # If we're using the Oceananigans development environment,
    # i.e. running from the git repository. Really we should not
    # use untagged versions for real science. It's not as reproducible.
    if "Oceananigans" == project.name
        return "Oceananigans v$(project.version) (DEVELOPMENT BRANCH)"
    end

    @warn "Could not determine Oceananigans version info."

    return ""
end


text:
function Base.show(io::IO, ir::YaoIR)
    indent = get(io, :indent, 0)
    printstyled(io, tab^indent, "circuit", color = :light_blue, bold = true)
    print(io, tab^indent, " ")
    print_head(io, ir.name)
    print_args(io, ir.args)

    if !isempty(ir.whereparams)
        print_where(io, ir.whereparams)
    end
    println(io)
    print(io, ir.body)
end

function print_head(io::IO, name)
    name isa Expr && name.head === :(::) && print(io, "(")
    print_name(io, name)
    name isa Expr && name.head === :(::) && print(io, ")")
end

print_name(io::IO, name::Symbol) = printstyled(io, name)

function print_name(io::IO, name::Expr)
    if name.head === :(::)
        print(io, name.args[1], "::")
        printstyled(io, name.args[2]; color = :green, bold = true)
    elseif name.head === :(kw)
        print(io, name.args[1], "=", name.args[2])
    else # not sure this is the case, but print it anyway
        print(io, name)
    end
end

function print_args(io::IO, args::Vector{Any})
    print(io, "(")
    for (k, each) in enumerate(args)
        print_name(io, each)
        if k != lastindex(args)
            print(io, ", ")
        end
    end
    print(io, ")")
end

function print_where(io::IO, whereparams::Vector{Any})
    printstyled(io, " where ", color = :light_blue)
    print(io, "{")
    for (k, each) in enumerate(whereparams)
        print(io, each)
        if k != lastindex(whereparams)
            print(io, ", ")
        end
    end
    print(io, "}")
end

const printers = Dict{Symbol,Any}()
# modified show(io::IO, b::Block) from IRTools: print.jl

function print_quantum(io, ex)
    printstyled(io, ex.args[1]; color = :light_blue, bold = true)
    print(io, "(")
    #         location          gate
    print(io, ex.args[2])
    for each in ex.args[3:end]
        print(io, ", ")
        print(io, each)
    end
    print(io, ")")
end

printers[:gate] = print_quantum
printers[:ctrl] = print_quantum
printers[:measure] = print_quantum
printers[:register] = function print_register_meta(io, ex)
    if ex.args[2] === :new
        printstyled(io, "%new%"; color = :light_blue, bold = true)
        print(io, "(", ex.args[3], ")")
    elseif ex.args[2] === :prev
        printstyled(io, "%prev%"; color = :light_blue, bold = true)
    else
        print(io, ex)
    end
end


text:
using Test, AlgoritmiStruttureDati.Laboratorio3
a = [1, 2, 3, 4, 5]

@testset "sumOfMult.jl" begin
    @test sumOfMult(a) == 40
end


text:
"""
SnoopCompile allows you to collect and analyze data on actions of Julia's compiler.

The capabilities depend on your version of Julia; in general, the capabilities that
require more recent Julia versions are also the most powerful and useful. When possible,
you should prefer them above the more limited tools available on earlier versions.

## All Julia versions

- `@snoopc`: record Julia's code generation
- `SnoopCompile.read`: parse data collected from `@snoopc`
- `SnoopCompile.parcel`: split precompile statements into modules/packages
- `SnoopCompile.write`: write module-specific precompile files (")

## At least Julia 1.2

- `@snoopi`: record entrances to Julia's type-inference (`parcel` and `write` work on these data, too)

## At least Julia 1.6

### Invalidations

- `@snoopr`: record invalidations
- `uinvalidated`: collect unique method invalidations from `@snoopr`
- `invalidation_trees`: organize invalidation data into trees
- `filtermod`: select trees that invalidate methods in particular modules
- `findcaller`: find a path through invalidation trees reaching a particular method
- `ascend`: interactive analysis of an invalidation tree

### LLVM

- `@snoopl`: record data about the actions of LLVM, the library used to generate native code
- `read_snoopl`: parse data collected by `@snoopl`

### "Deep" data on inference

- `@snoopi_deep`: record more extensive data about type-inference (`parcel` and `write` work on these data, too)
- `flamegraph`: prepare a visualization from `@snoopi_deep`
- `flatten`: reduce the tree format recorded by `@snoopi_deep` to list format
- `accumulate_by_source`: aggregate list items by their source
- `inference_triggers`: extract data on the triggers of inference
- `callerinstance`, `callingframe`, `skiphigherorder`, and `InferenceTrigger`: manipulate stack frames from `inference_triggers`
- `ascend`: interactive analysis of an inference-triggering call chain
- `runtime_inferencetime`: profile-guided deoptimization
"""
module SnoopCompile

using SnoopCompileCore
export @snoopc
# More exports are defined below in the conditional loading sections

using Core: MethodInstance, CodeInfo
using InteractiveUtils
using Serialization
using Printf
using OrderedCollections
import YAML  # For @snoopl
using Requires

# Parcel Regex
const anonrex = r"#{1,2}\d+#{1,2}\d+"         # detect anonymous functions
const kwrex = r"^#kw##(.*)$|^#([^#]*)##kw$"   # detect keyword-supplying functions
const kwbodyrex = r"^##(\w[^#]*)#\d+"         # detect keyword body methods
const genrex = r"^##s\d+#\d+$"                # detect generators for @generated functions
const innerrex = r"^#[^#]+#\d+"               # detect inner functions

# This is for SnoopCompile's own directives. You don't want to call this from packages because then
# SnoopCompile becomes a dependency of your package. Instead, make sure that `writewarnpcfail` is set to `true`
# in `SnoopCompile.write` and a copy of this macro will be placed at the top
# of your precompile files.
macro warnpcfail(ex::Expr)
    modl = __module__
    file = __source__.file === nothing ? "?" : String(__source__.file)
    line = __source__.line
    quote
        $(esc(ex)) || @warn """precompile directive
     $($(Expr(:quote, ex)))
 failed. Please report an issue in $($modl) (after checking for duplicates) or remove this directive.""" _file=$file _line=$line
    end
end

# Parcel
include("parcel_snoopc.jl")

if isdefined(SnoopCompileCore, Symbol("@snoopi"))
    include("parcel_snoopi.jl")
    export @snoopi
end

if isdefined(SnoopCompileCore, Symbol("@snoopi_deep"))
    include("parcel_snoopi_deep.jl")
    include("deep_demos.jl")
    export @snoopi_deep, exclusive, inclusive, flamegraph, flatten, accumulate_by_source, collect_for, runtime_inferencetime, staleinstances
    export InferenceTrigger, inference_triggers, callerinstance, callingframe, skiphigherorder, trigger_tree, suggest, isignorable
    export report_callee, report_caller, report_callees
end

if isdefined(SnoopCompileCore, Symbol("@snoopl"))
    export @snoopl
end
# To support reading of results on an older Julia version, this isn't conditional
include("parcel_snoopl.jl")
export read_snoopl

if isdefined(SnoopCompileCore, Symbol("@snoopr"))
    include("invalidations.jl")
    export @snoopr, uinvalidated, invalidation_trees, filtermod, findcaller, ascend
end

if isdefined(SnoopCompileCore, Symbol("@snoopr")) && isdefined(SnoopCompileCore, Symbol("@snoopi_deep"))
    include("invalidation_and_inference.jl")
    export precompile_blockers
end

# Write
include("write.jl")

function __init__()
    if isdefined(SnoopCompile, :runtime_inferencetime)
        @require PyPlot = "d330b81b-6aea-500a-939a-2ce795aea3ee" include("visualizations.jl")
    end
end

end # module


text:
using ClimateModels
using MITgcmTools
using Pluto
using PlutoUI
import Plots
import CairoMakie

#MC=MITgcm_config(configuration="global_with_exf")
#setup(MC)
# build(MC,"--allow-skip")
# launch(MC)

tmp=ModelConfig(model=ClimateModels.RandomWalker)
setup(tmp)
launch(tmp)

max_stars_repo_name: crbinz/GlobalOceanNotebooks

max_issues_repo_name: crbinz/GlobalOceanNotebooks

max_forks_repo_name: crbinz/GlobalOceanNotebooks


text:
export CURANDError

struct CURANDError <: Exception
    code::curandStatus_t
end

Base.convert(::Type{curandStatus_t}, err::CURANDError) = err.code

Base.showerror(io::IO, err::CURANDError) =
    print(io, "CURANDError: ", description(err), " (code $(reinterpret(Int32, err.code)), $(name(err)))")

name(err::CURANDError) = string(err.code)

## COV_EXCL_START
function description(err)
    if err.code == CURAND_STATUS_SUCCESS
        "generator was created successfully"
    elseif err.code == CURAND_STATUS_VERSION_MISMATCH
        "header file and linked library version do not match"
    elseif err.code == CURAND_STATUS_NOT_INITIALIZED
        "generator not initialized"
    elseif err.code == CURAND_STATUS_ALLOCATION_FAILED
        "memory allocation failed"
    elseif err.code == CURAND_STATUS_TYPE_ERROR
        "generator is wrong type"
    elseif err.code == CURAND_STATUS_OUT_OF_RANGE
        "argument out of range"
    elseif err.code == CURAND_STATUS_LENGTH_NOT_MULTIPLE
        "length requested is not a multple of dimension"
    elseif err.code == CURAND_STATUS_DOUBLE_PRECISION_REQUIRED
        "GPU does not have double precision required by MRG32k3a"
    elseif err.code == CURAND_STATUS_LAUNCH_FAILURE
        "kernel launch failure"
    elseif err.code == CURAND_STATUS_PREEXISTING_FAILURE
        "preexisting failure on library entry"
    elseif err.code == CURAND_STATUS_INITIALIZATION_FAILED
        "initialization of CUDA failed"
    elseif err.code == CURAND_STATUS_ARCH_MISMATCH
        "architecture mismatch, GPU does not support requested feature"
    elseif err.code == CURAND_STATUS_INTERNAL_ERROR
        "internal library error"
    else
        "no description for this error"
    end
end
## COV_EXCL_STOP

## API call wrapper

# outlined functionality to avoid GC frame allocation
@noinline function throw_api_error(res)
    throw(CURANDError(res))
end

function initialize_api()
    CUDAnative.prepare_cuda_call()
end

macro check(ex)
    quote
        res = @retry_reclaim CURAND_STATUS_ALLOCATION_FAILED $(esc(ex))
        if res != CURAND_STATUS_SUCCESS
            throw_api_error(res)
        end

        return
    end
end