text:
abstract type GraphMorphismProblem end
struct IsomorphismProblem <: GraphMorphismProblem end
struct SubGraphIsomorphismProblem <: GraphMorphismProblem end
struct InducedSubGraphIsomorphismProblem <: GraphMorphismProblem end

"""
    IsomorphismAlgorithm

An abstract type used for method dispatch on isomorphism functions.
"""
abstract type IsomorphismAlgorithm end

"""
    could_have_isomorph(g1, g2)

Run quick test to check if `g1 and `g2` could be isomorphic.

If the result is `false`, then `g1` and `g2` are definitely not isomorphic,
but if the result is `true` this is not guaranteed.
"""
function could_have_isomorph(g1::AbstractGraph, g2::AbstractGraph)
    nv(g1) == nv(g2) || return false

    indegs1 = indegree(g1)
    indegs2 = indegree(g2)
    sort!(indegs1)
    sort!(indegs2)
    indegs1 == indegs2 || return false

    if LightGraphs.is_directed(g1) || LightGraphs.is_directed(g2)
        outdegs1 = outdegree(g1)
        outdegs2 = outdegree(g2)
        sort!(outdegs1)
        sort!(outdegs2)
        outdegs1 == outdegs2 || return false
    end

    return true
end

"""
    has_induced_subgraphisomorph(g1, g2, alg::IsomorphismAlgorithm=VF2(); vertex_relation=nothing, edge_relation=nothing)

Return `true` if the graph `g1` contains a vertex induced subgraph that is isomorphic to `g2`.

### Optional Arguments
- `alg`: The algorithm that is used to find the induced subgraph isomorphism. Can be only
    `VF2()` at the moment.
- `vertex_relation`: A binary function that takes a vertex from `g1` and one from `g2`. An
    isomorphism only exists if this function returns `true` for all matched vertices.
- `edge_relation`: A binary function that takes an edge from `g1` and one from `g2`. An
    isomorphism only exists if this function returns `true` for all matched edges.

### Examples
```doctest.jl
julia> has_induced_subgraphisomorph(CompleteGraph(5), CompleteGraph(4))
true
julia> has_induced_subgraphisomorph(CompleteGraph(5), CycleGraph(4))
false

julia> g1 = PathDiGraph(3); color1 = [1, 1, 1]
julia> g2 = PathDiGraph(2); color2 = [1, 2]
julia> color_rel(u, v) = (color1[u] == color2[v])
julia> has_induced_subgraphisomorph(g1, g2)
true
julia> has_induced_subgraphisomorph(g1, g2, vertex_relation=color_rel)
false
```
### See also
[`has_subgraphisomorph`](@ref), [`has_isomorph`](@ref), [`count_induced_subgraphisomorph`](@ref), [`all_induced_subgraphisomorph`](@ref)
"""
function has_induced_subgraphisomorph(g1::AbstractGraph, g2::AbstractGraph, alg::IsomorphismAlgorithm=VF2();
                                 vertex_relation::Union{Nothing, Function}=nothing,
                                 edge_relation::Union{Nothing, Function}=nothing)::Bool

        has_induced_subgraphisomorph(g1, g2, alg; vertex_relation=vertex_relation, edge_relation=edge_relation)

end

"""
    has_subgraphisomorph(g1, g2, alg::IsomorphismAlgorithm=VF2(); vertex_relation=nothing, edge_relation=nothing)

Return `true` if the graph `g1` contains a subgraph that is isomorphic to `g2`.

### Optional Arguments
- `alg`: The algorithm that is used to find the induced subgraph isomorphism. Can be only
    `VF2()` at the moment.
- `vertex_relation`: A binary function that takes a vertex from `g1` and one from `g2`. An
    isomorphism only exists if this function returns `true` for all matched vertices.
- `edge_relation`: A binary function that takes an edge from `g1` and one from `g2`. An
    isomorphism only exists if this function returns `true` for all matched edges.

### Examples
```doctest.jl
julia> has_subgraphisomorph(CompleteGraph(5), CompleteGraph(4))
true
julia> has_subgraphisomorph(CompleteGraph(5), CycleGraph(4))
true

julia> g1 = PathDiGraph(3); color1 = [1, 1, 1]
julia> g2 = PathDiGraph(2); color2 = [1, 2]
julia> color_rel(u, v) = (color1[u] == color2[v])
julia> has_subgraphisomorph(g1, g2)
true
julia> has_subgraphisomorph(g1, g2, vertex_relation=color_rel)
false
```
### See also
[`has_induced_subgraphisomorph`](@ref), [`has_isomorph`](@ref), [`count_subgraphisomorph`](@ref), [`all_subgraphisomorph`](@ref)
"""
function has_subgraphisomorph(g1::AbstractGraph, g2::AbstractGraph, alg::IsomorphismAlgorithm=VF2();
                                 vertex_relation::Union{Nothing, Function}=nothing,
                                 edge_relation::Union{Nothing, Function}=nothing)::Bool
    has_subgraphisomorph(g1, g2, alg;
        vertex_relation=vertex_relation,
        edge_relation=edge_relation)
end

"""
    has_isomorph(g1, g2, alg::IsomorphismAlgorithm=VF2(); vertex_relation=nothing, edge_relation=nothing)

Return `true` if the graph `g1` is isomorphic to `g2`.

### Optional Arguments
- `alg`: The algorithm that is used to find the induced subgraph isomorphism. Can be only
    `VF2()` at the moment.
- `vertex_relation`: A binary function that takes a vertex from `g1` and one from `g2`. An
    isomorphism only exists if this function returns `true` for all matched vertices.
- `edge_relation`: A binary function that takes an edge from `g1` and one from `g2`. An
    isomorphism only exists if this function returns `true` for all matched edges.

### Examples
```doctest.jl
julia> has_isomorph(CompleteGraph(3), CycleGraph(3))
true
julia> has_isomorph(CompleteGraph(4), CycleGraph(4))
false

julia> g1 = PathDiGraph(4); color1 = [1, 2, 1, 1]
julia> g2 = PathDiGraph(4); color2 = [1, 2, 2, 1]
julia> color_rel(u, v) = (color1[u] == color2[v])
julia> has_isomorph(g1, g2)
true
julia> has_isomorph(g1, g2, vertex_relation=color_rel)
false
```
### See also
[`has_induced_subgraphisomorph`](@ref), [`has_subgraphisomorph`](@ref), [`count_subgraphisomorph`](@ref), [`all_subgraphisomorph`](@ref)
"""
function has_isomorph(g1::AbstractGraph, g2::AbstractGraph, alg::IsomorphismAlgorithm=VF2();
                         vertex_relation::Union{Nothing, Function}=nothing,
                         edge_relation::Union{Nothing, Function}=nothing)::Bool
    has_isomorph(g1, g2, alg;
                  vertex_relation=vertex_relation,
                  edge_relation=edge_relation)
end

"""
    count_induced_subgraphisomorph(g1, g2, alg::IsomorphismAlgorithm=VF2(); vertex_relation=nothing, edge_relation=nothing)

Return the number of vertex induced subgraphs of the graph `g1` that are isomorphic to `g2`.

### Optional Arguments
- `vertex_relation`: A binary function that takes a vertex from `g1` and one from `g2`. An
    isomorphism only exists if this function returns `true` for all matched vertices.
- `edge_relation`: A binary function that takes an edge from `g1` and one from `g2`. An
    isomorphism only exists if this function returns `true` for all matched edges.

### Examples
```doctest.jl
julia> count_induced_subgraphisomorph(CompleteGraph(5), CompleteGraph(4))
120
julia> count_induced_subgraphisomorph(CompleteGraph(5), CycleGraph(4))
0

julia> g1 = PathGraph(3); color1 = [1, 1, 2]
julia> g2 = PathGraph(2); color2 = [1, 2]
julia> color_rel(u, v) = (color1[u] == color2[v])
julia> count_induced_subgraphisomorph(g1, g2)
2
julia> count_induced_subgraphisomorph(g1, g2, vertex_relation=color_rel)
1
```
### See also
[`count_subgraphisomorph`](@ref), [`count_isomorph`](@ref), [`has_induced_subgraphisomorph`](@ref), [`all_induced_subgraphisomorph`](@ref)
"""
function count_induced_subgraphisomorph(g1::AbstractGraph, g2::AbstractGraph, alg::IsomorphismAlgorithm=VF2();
                                   vertex_relation::Union{Nothing, Function}=nothing,
                                   edge_relation::Union{Nothing, Function}=nothing)::Int
    count_induced_subgraphisomorph(g1,g2, alg;
                                   vertex_relation=vertex_relation,
                                   edge_relation=edge_relation)
end

"""
    count_subgraphisomorph(g1, g2, alg::IsomorphismAlgorithm=VF2(); vertex_relation=nothing, edge_relation=nothing)

Return the number of subgraphs of the graph `g1` that are isomorphic to `g2`.

### Optional Arguments
- `alg`: The algorithm that is used to find the induced subgraph isomorphism. Can be only
    `VF2()` at the moment.
- `vertex_relation`: A binary function that takes a vertex from `g1` and one from `g2`. An
    isomorphism only exists if this function returns `true` for all matched vertices.
- `edge_relation`: A binary function that takes an edge from `g1` and one from `g2`. An
    isomorphism only exists if this function returns `true` for all matched edges.

### Examples
```doctest.jl
julia> count_subgraphisomorph(CompleteGraph(5), CompleteGraph(4))
120
julia> count_subgraphisomorph(CompleteGraph(5), CycleGraph(4))
120

julia> g1 = CycleDiGraph(3); color1 = [1, 1, 2]
julia> g2 = SimpleDiGraph(2); color2 = [1, 2]
julia> color_rel(u, v) = (color1[u] == color2[v])
julia> count_subgraphisomorph(g1, g2)
6
julia> count_subgraphisomorph(g1, g2, vertex_relation=color_rel)
2
```
### See also
[`count_induced_subgraphisomorph`](@ref), [`count_isomorph`](@ref), [`has_subgraphisomorph`](@ref), [`all_subgraphisomorph`](@ref)
"""
function count_subgraphisomorph(g1::AbstractGraph, g2::AbstractGraph, alg::IsomorphismAlgorithm=VF2();
                                   vertex_relation::Union{Nothing, Function}=nothing,
                                   edge_relation::Union{Nothing, Function}=nothing)::Int

     count_subgraphisomorph(g1::AbstractGraph, g2::AbstractGraph, VF2();
                            vertex_relation=vertex_relation,
                            edge_relation=edge_relation)
end

"""
    count_isomorph(g1, g2, alg::IsomorphismAlgorithm=VF2(); vertex_relation=nothing, edge_relation=nothing)

Return the number of isomorphism from graph `g1` to `g2`.

### Optional Arguments
- `alg`: The algorithm that is used to find the induced subgraph isomorphism. Can be only
    `VF2()` at the moment.
- `vertex_relation`: A binary function that takes a vertex from `g1` and one from `g2`. An
    isomorphism only exists if this function returns `true` for all matched vertices.
- `edge_relation`: A binary function that takes an edge from `g1` and one from `g2`. An
    isomorphism only exists if this function returns `true` for all matched edges.

### Examples
```doctest.jl
julia> count_isomorph(CycleGraph(5), CycleGraph(5))
10
julia> count_isomorph(CompleteGraph(5), CycleGraph(5))
0

julia> g1 = CycleDiGraph(3); color1 = [1, 1, 2]
julia> g2 = CycleDiGraph(3); color2 = [1, 1, 1]
julia> color_rel(u, v) = (color1[u] == color2[v])
julia> count_isomorph(g1, g2)
3
julia> count_isomorph(g1, g2, vertex_relation=color_rel)
0
```
### See also
[`count_induced_subgraphisomorph`](@ref), [`count_subgraphisomorph`](@ref), [`has_isomorph`](@ref), [`all_isomorph`](@ref)
"""
function count_isomorph(g1::AbstractGraph, g2::AbstractGraph, alg::IsomorphismAlgorithm=VF2();
                        vertex_relation::Union{Nothing, Function}=nothing,
                        edge_relation::Union{Nothing, Function}=nothing)::Int
    count_isomorph(g1, g2, alg;
                   vertex_relation=vertex_relation,
                   edge_relation=edge_relation)
end

"""
    all_induced_subgraphisomorph(g1, g2, alg::IsomorphismAlgorithm=VF2(); vertex_relation=nothing, edge_relation=nothing)

Return all isomorphism from vertex induced subgraphs of `g1` to `g2`.
The isomorphisms are returned as an iterator of vectors of tuples, where the i-th vector is
the i-th isomorphism and a tuple (u, v) in this vector means that u ∈ g1 is
mapped to v ∈ g2.

### Optional Arguments
- `alg`: The algorithm that is used to find the induced subgraph isomorphism. Can be only
    `VF2()` at the moment.
- `vertex_relation`: A binary function that takes a vertex from `g1` and one from `g2`. An
    isomorphism only exists if this function returns `true` for all matched vertices.
- `edge_relation`: A binary function that takes an edge from `g1` and one from `g2`. An
    isomorphism only exists if this function returns `true` for all matched edges.

### Examples
```doctest.jl
julia> all_induced_subgraphisomorph(PathGraph(3), SimpleGraph(2)) |> collect
2-element Array{Array{Tuple{Int64,Int64},1},1}:
 [(1, 1), (3, 2)]
 [(3, 1), (1, 2)]

julia> g1 = PathDiGraph(3); color1 = [1, 1, 2]
julia> g2 = PathDiGraph(2); color2 = [1, 2]
julia> color_rel(u, v) = (color1[u] == color2[v])
julia> all_induced_subgraphisomorph(g1, g2) |> collect
2-element Array{Array{Tuple{Int64,Int64},1},1}:
 [(1, 1), (2, 2)]
 [(2, 1), (3, 2)]
julia> all_induced_subgraphisomorph(g1, g2, vertex_relation=color_rel) |> collect
1-element Array{Array{Tuple{Int64,Int64},1},1}:
 [(1, 1), (2, 2)]
```
### See also
[`all_subgraphisomorph`](@ref), [`all_isomorph`](@ref), [`has_induced_subgraphisomorph`](@ref), [`count_induced_subgraphisomorph`](@ref)
"""
function all_induced_subgraphisomorph(g1::AbstractGraph, g2::AbstractGraph, alg::IsomorphismAlgorithm=VF2();
                                      vertex_relation::Union{Nothing, Function}=nothing,
                                      edge_relation::Union{Nothing, Function}=nothing)::Channel{Vector{Tuple{eltype(g1),eltype(g2)}}}
    all_induced_subgraphisomorph(g1, g2, alg;
        vertex_relation=vertex_relation,
        edge_relation=edge_relation)
end

"""
    all_subgraphisomorph(g1, g2, alg::IsomorphismAlgorithm=VF2(); vertex_relation=nothing, edge_relation=nothing)

Return all isomorphism from  subgraphs of `g1` to `g2`.
The isomorphisms are returned as an iterator of vectors of tuples, where the i-th vector is
the i-th isomorphism and a tuple (u, v) in this vector means that u ∈ g1 is
mapped to v ∈ g2.

### Optional Arguments
- `alg`: The algorithm that is used to find the induced subgraph isomorphism. Can be only
    `VF2()` at the moment.
- `vertex_relation`: A binary function that takes a vertex from `g1` and one from `g2`. An
    isomorphism only exists if this function returns `true` for all matched vertices.
- `edge_relation`: A binary function that takes an edge from `g1` and one from `g2`. An
    isomorphism only exists if this function returns `true` for all matched edges.

### Examples
```doctest.jl
julia> all_subgraphisomorph(PathGraph(3), PathGraph(2)) |> collect
4-element Array{Array{Tuple{Int64,Int64},1},1}:
 [(1, 1), (2, 2)]
 [(2, 1), (1, 2)]
 [(2, 1), (3, 2)]
 [(3, 1), (2, 2)]

julia> g1 = PathDiGraph(3); color1 = [1, 1, 2]
julia> g2 = PathDiGraph(2); color2 = [1, 2]
julia> color_rel(u, v) = (color1[u] == color2[v])
julia> all_subgraphisomorph(g1, g2) |> collect
2-element Array{Array{Tuple{Int64,Int64},1},1}:
 [(1, 1), (2, 2)]
 [(2, 1), (3, 2)]
julia> all_subgraphisomorph(g1, g2, vertex_relation=color_rel)
1-element Array{Array{Tuple{Int64,Int64},1},1}:
 [(2, 1), (3, 2)]
```
### See also
[`all_induced_subgraphisomorph`](@ref), [`all_isomorph`](@ref), [`has_subgraphisomorph`](@ref), [`count_subgraphisomorph`](@ref)
"""
function all_subgraphisomorph(g1::AbstractGraph, g2::AbstractGraph, alg::IsomorphismAlgorithm=VF2();
                         vertex_relation::Union{Nothing, Function}=nothing,
                         edge_relation::Union{Nothing, Function}=nothing)::Channel{Vector{Tuple{eltype(g1), eltype(g2)}}}

    all_subgraphisomorph(g1, g2, alg;
        vertex_relation=vertex_relation,
        edge_relation=edge_relation)
end

"""
    all_isomorph(g1, g2, alg::IsomorphismAlgorithm=VF2(); vertex_relation=nothing, edge_relation=nothing)

Return all isomorphism from `g1` to `g2`.
The isomorphisms are returned as an iterator of vectors of tuples, where the i-th vector is
the i-th isomorphism and a tuple (u, v) in this vector means that u ∈ g1 is
mapped to v ∈ g2.

### Optional Arguments
- `alg`: The algorithm that is used to find the induced subgraph isomorphism. Can be only
    `VF2()` at the moment.
- `vertex_relation`: A binary function that takes a vertex from `g1` and one from `g2`. An
    isomorphism only exists if this function returns `true` for all matched vertices.
- `edge_relation`: A binary function that takes an edge from `g1` and one from `g2`. An
    isomorphism only exists if this function returns `true` for all matched edges.

### Examples
```doctest.jl
julia> all_isomorph(StarGraph(4), StarGraph(4)) |> collect
6-element Array{Array{Tuple{Int64,Int64},1},1}:
 [(1, 1), (2, 2), (3, 3), (4, 4)]
 [(1, 1), (2, 2), (4, 3), (3, 4)]
 [(1, 1), (3, 2), (2, 3), (4, 4)]
 [(1, 1), (3, 2), (4, 3), (2, 4)]
 [(1, 1), (4, 2), (2, 3), (3, 4)]
 [(1, 1), (4, 2), (3, 3), (2, 4)]

julia> g1 = CycleDiGraph(3); color1 = [1, 1, 2]
julia> g2 = CycleDiGraph(3); color2 = [2, 1, 1]
julia> color_rel(u, v) = (color1[u] == color2[v])
julia> all_isomorph(g1, g2) |> collect
3-element Array{Array{Tuple{Int64,Int64},1},1}:
 [(1, 1), (2, 2), (3, 3)]
 [(2, 1), (3, 2), (1, 3)]
 [(3, 1), (1, 2), (2, 3)]
julia> all_subgraphisomorph(g1, g2, vertex_relation=color_rel)
1-element Array{Array{Tuple{Int64,Int64},1},1}:
 [(3, 1), (1, 2), (2, 3)]
```
### See also
[`all_induced_subgraphisomorph`](@ref), [`all_subgraphisomorph`](@ref), [`has_isomorph`](@ref), [`count_isomorph`](@ref)
"""
function all_isomorph(g1::AbstractGraph, g2::AbstractGraph, alg::IsomorphismAlgorithm=VF2();
                 vertex_relation::Union{Nothing, Function}=nothing,
                 edge_relation::Union{Nothing, Function}=nothing)::Channel{Vector{Tuple{eltype(g1),eltype(g2)}}}

    all_isomorph(g1, g2, alg;
        vertex_relation=vertex_relation,
        edge_relation=edge_relation)
end


text:
Base.eltype(fg::AFG) = eltype(typeof(fg))
Base.length(fg::AFG) = length(vertices(fg))
nvertices(fg::AFG) = length(fg)
nvertices(::Type{Simplex{n, T}}) where {n, T} = n
nvertices(::Type{HyperCube{N,T}}) where {N,T} = 2^N
nvertices(::Type{HyperRectangle{N,T}}) where {N,T} = 2^N
nvertices(s::Union{AbstractGeometry, Simplex}) = nvertices(typeof(s))

spacedim(s) = length(eltype(s))
spacedim(::HyperCube{N,T}) where {N,T} = N
numtype(s) = eltype(eltype(s))
numtype(::HyperCube{N,T}) where {N,T} = T

Base.push!(fl::AFG, pt) = (push!(vertices(fl), pt); fl)
Base.deleteat!(c::AbstractFlexibleGeometry, i) = (deleteat!(vertices(c), i); c)

Base.copy(fl::FG) where {FG <: AFG} = FG(copy(vertices(fl)))
push(fl::AFG, pt) = push!(copy(fl), pt)

vertices(x::AbstractFlexibleGeometry) = x.vertices
vertices(s::Simplex) = Tuple(s)

standard_cube_vertices(::Type{Val{1}}) = [Vec(0), Vec(1)]
_vcat(v1,v2) = Vec(Tuple(v1)..., Tuple(v2)...)
function _combine_vcat(arr1, arr2)
    T = typeof(_vcat(first(arr1), first(arr2)))
    ret = T[]
    for v in arr1, w in arr2
        push!(ret, _vcat(v,w))
    end
    ret
end

@generated function standard_cube_vertices(::Type{Val{N}}) where N
    @assert N::Int > 0
    quote
        vert_last = standard_cube_vertices($(Val{N-1}))
        vert_1 = standard_cube_vertices(Val{1})
        _combine_vcat(vert_1, vert_last)
    end
end
@generated function vertices_rettype(r::HyperRectangle)
    N = ndims(r)
    T = eltype(r)
    RET = NTuple{(2^N), Vec{N,T}}
    :($RET)
end

function vertices(r::HyperRectangle)
    N = ndims(r)
    ret_type = vertices_rettype(r)
    o = origin(r)
    v = widths(r)
    f(sv) = o + sv .* v
    tuple(map(f, standard_cube_vertices(Val{N}))...)::ret_type
end

vertices(c::HyperCube) = vertices(convert(HyperRectangle, c))

vertexmat(s) = hcat(vertices(s)...)
vertexmatrix(s::AbstractConvexHull) = Matrix(vertexmat(s))::Matrix{numtype(s)}

(::Type{F})(g::Union{AbstractSimplex, AFG, GeometryPrimitive}) where {F <: AFG} = F(vertices(g))
(::Type{F})(v::NTuple) where {F <: AFG} = F(collect(v))
Base.convert(::Type{F}, s::Simplex) where {F <: AFG} = F(s)

(::Type{R})(c::HyperCube) where {R <: HyperRectangle} = R(origin(c), widths(c))
Base.convert(::Type{R}, c::HyperCube) where {R <: HyperRectangle} = R(c)
Base.convert(::Type{F}, c::HyperCube) where {F <: AFG} = F(c)

function Base.isapprox(s1::AbstractConvexHull, s2::AbstractConvexHull;kw...)
    isapprox(vertexmat(s1), vertexmat(s2); kw...)
end

# could add an @symmetric macro, which defines min_euclidean(pt, s) from
# min_euclidean(s, pt) automatically etc.
@inline min_euclidean(pt1::Vec, pt2::Vec) = norm(pt1-pt2)
min_euclidean(pt::Vec, s::Simplex) = √(sqdist(pt, s))
min_euclidean(s::Simplex, pt::Vec) = min_euclidean(pt, s)
min_euclidean(r1::HyperRectangle, r2::Vec) = sqrt(min_euclideansq(r1, r2))
min_euclidean(r1::Vec, r2::HyperRectangle) = sqrt(min_euclideansq(r1, r2))
min_euclidean(r1::HyperRectangle, r2::HyperRectangle) = sqrt(min_euclideansq(r1, r2))
min_euclidean(c1::AbstractConvexHull, pt::Vec) = gjk(c1,pt)
min_euclidean(pt::Vec, c2::AbstractConvexHull) = gjk(pt,c2)
min_euclidean(c1::AbstractConvexHull, c2::AbstractConvexHull) = gjk(c1,c2)

ext: jl


text:
for Typ in (:Ldiv, :Rdiv)
    @eval begin
        struct $Typ{StyleA, StyleB, AType, BType}
            A::AType
            B::BType
        end

        @inline $Typ{StyleA, StyleB}(A::AType, B::BType) where {StyleA,StyleB,AType,BType} =
            $Typ{StyleA,StyleB,AType,BType}(A,B)

        @inline $Typ(A::AType, B::BType) where {AType,BType} =
            $Typ{typeof(MemoryLayout(AType)),typeof(MemoryLayout(BType)),AType,BType}(A, B)

        @inline BroadcastStyle(::Type{<:$Typ}) = ApplyBroadcastStyle()
        @inline broadcastable(M::$Typ) = M

        similar(A::$Typ, ::Type{T}, axes) where T = similar(Array{T}, axes)
        similar(A::$Typ, ::Type{T}) where T = similar(A, T, axes(A))
        similar(A::$Typ) = similar(A, eltype(A))

        @inline copy(M::$Typ) = copyto!(similar(M), M)
        @inline materialize(M::$Typ) = copy(instantiate(M))
    end
end

@inline _ldivaxes(::Tuple{}, ::Tuple{}) = ()
@inline _ldivaxes(::Tuple{}, Bax::Tuple) = Bax
@inline _ldivaxes(::Tuple{<:Any}, ::Tuple{<:Any}) = ()
@inline _ldivaxes(::Tuple{<:Any}, Bax::Tuple{<:Any,<:Any}) = (OneTo(1),last(Bax))
@inline _ldivaxes(Aax::Tuple{<:Any,<:Any}, ::Tuple{<:Any}) = (last(Aax),)
@inline _ldivaxes(Aax::Tuple{<:Any,<:Any}, Bax::Tuple{<:Any,<:Any}) = (last(Aax),last(Bax))

@inline ldivaxes(A, B) = _ldivaxes(axes(A), axes(B))

@inline axes(L::Ldiv) = ldivaxes(L.A, L.B)
@inline size(L::Ldiv) = map(length, axes(L))
@inline length(L::Ldiv{<:Any,<:Any,<:Any,<:AbstractVector}) =size(L.A, 2)

@inline size(L::Rdiv) = (size(L.A, 1),size(L.B,1))
@inline axes(L::Rdiv) = (axes(L.A, 1),axes(L.B,1))

@inline ndims(L::Ldiv) = ndims(L.B)
@inline ndims(L::Rdiv) = 2
@inline eltype(M::Ldiv) = promote_type(Base.promote_op(inv, eltype(M.A)), eltype(M.B))
@inline eltype(M::Rdiv) = promote_type(eltype(M.A), Base.promote_op(inv, eltype(M.B)))

check_ldiv_axes(A, B) =
    axes(A,1) == axes(B,1) || throw(DimensionMismatch("First axis of A, $(axes(A,1)), and first axis of B, $(axes(B,1)) must match"))

check_rdiv_axes(A, B) =
    axes(A,2) == axes(B,2) || throw(DimensionMismatch("Second axis of A, $(axes(A,2)), and second axis of B, $(axes(B,2)) must match"))

@inline function instantiate(L::Ldiv)
    check_ldiv_axes(L.A, L.B)
    Ldiv(instantiate(L.A), instantiate(L.B))
end

@inline function instantiate(L::Rdiv)
    check_rdiv_axes(L.A, L.B)
    Rdiv(instantiate(L.A), instantiate(L.B))
end

__ldiv!(::Mat, ::Mat, B) where Mat = error("Overload materialize!(::Ldiv{$(typeof(MemoryLayout(Mat))),$(typeof(MemoryLayout(typeof(B))))})")
__ldiv!(_, F, B) = ldiv!(F, B)
@inline _ldiv!(A, B) = __ldiv!(A, factorize(A), B)
@inline _ldiv!(A::Factorization, B) = ldiv!(A, B)

@inline _ldiv!(dest, A, B) = ldiv!(dest, factorize(A), B)
@inline _ldiv!(dest, A::Factorization, B) = ldiv!(dest, A, B)
@inline _ldiv!(dest, A::Transpose{<:Any,<:Factorization}, B) = ldiv!(dest, A, B)
@inline _ldiv!(dest, A::Adjoint{<:Any,<:Factorization}, B) = ldiv!(dest, A, B)

@inline ldiv(A, B) = materialize(Ldiv(A,B))
@inline rdiv(A, B) = materialize(Rdiv(A,B))

@inline materialize!(M::Ldiv) = _ldiv!(M.A, M.B)
@inline materialize!(M::Rdiv) = materialize!(Lmul(M.B', M.A'))'
@inline copyto!(dest::AbstractArray, M::Rdiv) = copyto!(dest', Ldiv(M.B', M.A'))'

if VERSION ≥ v"1.1-pre"
    @inline copyto!(dest::AbstractArray, M::Ldiv) = _ldiv!(dest, M.A, M.B)
else
    @inline copyto!(dest::AbstractArray, M::Ldiv) = _ldiv!(dest, M.A, copy(M.B))
end

const MatLdivVec{styleA, styleB, T, V} = Ldiv{styleA, styleB, <:AbstractMatrix{T}, <:AbstractVector{V}}
const MatLdivMat{styleA, styleB, T, V} = Ldiv{styleA, styleB, <:AbstractMatrix{T}, <:AbstractMatrix{V}}
const BlasMatLdivVec{styleA, styleB, T<:BlasFloat} = MatLdivVec{styleA, styleB, T, T}
const BlasMatLdivMat{styleA, styleB, T<:BlasFloat} = MatLdivMat{styleA, styleB, T, T}

const MatRdivMat{styleA, styleB, T, V} = Rdiv{styleA, styleB, <:AbstractMatrix{T}, <:AbstractMatrix{V}}
const BlasMatRdivMat{styleA, styleB, T<:BlasFloat} = MatRdivMat{styleA, styleB, T, T}

# function materialize!(L::BlasMatLdivVec{<:AbstractColumnMajor,<:AbstractColumnMajor})

# end

macro _layoutldiv(Typ)
    esc(quote
        LinearAlgebra.ldiv!(A::$Typ, x::AbstractVector) = ArrayLayouts.materialize!(ArrayLayouts.Ldiv(A,x))
        LinearAlgebra.ldiv!(A::$Typ, x::AbstractMatrix) = ArrayLayouts.materialize!(ArrayLayouts.Ldiv(A,x))
        LinearAlgebra.ldiv!(A::$Typ, x::StridedVector) = ArrayLayouts.materialize!(ArrayLayouts.Ldiv(A,x))
        LinearAlgebra.ldiv!(A::$Typ, x::StridedMatrix) = ArrayLayouts.materialize!(ArrayLayouts.Ldiv(A,x))

        Base.:\(A::$Typ, x::AbstractVector) = ArrayLayouts.materialize(ArrayLayouts.Ldiv(A,x))
        Base.:\(A::$Typ, x::AbstractMatrix) = ArrayLayouts.materialize(ArrayLayouts.Ldiv(A,x))

        Base.:\(x::AbstractMatrix, A::$Typ) = ArrayLayouts.materialize(ArrayLayouts.Ldiv(x,A))
        Base.:\(x::Diagonal, A::$Typ) = ArrayLayouts.materialize(ArrayLayouts.Ldiv(x,A))

        Base.:\(x::$Typ, A::$Typ) = ArrayLayouts.materialize(ArrayLayouts.Ldiv(x,A))

        Base.:/(A::$Typ, x::AbstractVector) = ArrayLayouts.materialize(ArrayLayouts.Rdiv(A,x))
        Base.:/(A::$Typ, x::AbstractMatrix) = ArrayLayouts.materialize(ArrayLayouts.Rdiv(A,x))

        Base.:/(x::AbstractMatrix, A::$Typ) = ArrayLayouts.materialize(ArrayLayouts.Rdiv(x,A))
        Base.:/(x::Diagonal, A::$Typ) = ArrayLayouts.materialize(ArrayLayouts.Rdiv(x,A))

        Base.:/(x::$Typ, A::$Typ) = ArrayLayouts.materialize(ArrayLayouts.Rdiv(x,A))
    end)
end

macro layoutldiv(Typ)
    esc(quote
        ArrayLayouts.@_layoutldiv $Typ
        ArrayLayouts.@_layoutldiv UpperTriangular{T, <:$Typ{T}} where T
        ArrayLayouts.@_layoutldiv UnitUpperTriangular{T, <:$Typ{T}} where T
        ArrayLayouts.@_layoutldiv LowerTriangular{T, <:$Typ{T}} where T
        ArrayLayouts.@_layoutldiv UnitLowerTriangular{T, <:$Typ{T}} where T

        ArrayLayouts.@_layoutldiv UpperTriangular{T, <:SubArray{T,2,<:$Typ{T}}} where T
        ArrayLayouts.@_layoutldiv UnitUpperTriangular{T, <:SubArray{T,2,<:$Typ{T}}} where T
        ArrayLayouts.@_layoutldiv LowerTriangular{T, <:SubArray{T,2,<:$Typ{T}}} where T
        ArrayLayouts.@_layoutldiv UnitLowerTriangular{T, <:SubArray{T,2,<:$Typ{T}}} where T
    end)
end


text:
#
# This file is part of the YAActL.jl Julia package, MIT license
#
# Paul Bayer, 2020
#

using YAActL, Test

struct Incr <: Message
    x::Int
end
struct Decr <: Message
    x::Int
end

t = Ref{Task}()
lp = LinkParams(taskref=t)

a = [0]
inca(msg::Incr) = a[1] += msg.x
inca(msg::Decr) = ( become(deca); send!(self(), msg) )
deca(msg::Decr) = a[1] -= msg.x
deca(msg::Incr) = ( become(inca); send!(self(), msg))
fail(msg::Incr) = undefined

A = Actor(lp, inca)
@test A isa Link
@test t[].state == :runnable
act = YAActL.diag!(A)
@test act.self == A
@test act.bhv.f == inca

send!(A, Incr(10))
sleep(0.1)
@test act.res == 10
@test a[1] == 10

send!(A, Decr(5))
send!(A, Decr(5))
sleep(0.2)
@test a[1] == 0

@test !istaskfailed(A)
@test info(A) == :runnable

exit!(A)
sleep(0.1)
@test !isopen(A.chn)

B = Actor(fail)
send!(B, Incr(0))
sleep(0.1)
@test istaskfailed(B)
s = info(B)
@test s.exception isa UndefVarError


text:
struct Level{TA, TP, TR}
    A::TA
    P::TP
    R::TR
end

struct MultiLevel{S, Pre, Post, TA, TP, TR, TW}
    levels::Vector{Level{TA, TP, TR}}
    final_A::TA
    coarse_solver::S
    presmoother::Pre
    postsmoother::Post
    workspace::TW
end

struct MultiLevelWorkspace{TX, bs}
    coarse_xs::Vector{TX}
    coarse_bs::Vector{TX}
    res_vecs::Vector{TX}
end
function MultiLevelWorkspace(::Type{Val{bs}}, ::Type{T}) where {bs, T<:Number}
    if bs === 1
        TX = Vector{T}
    else
        TX = Matrix{T}
    end
    MultiLevelWorkspace{TX, bs}(TX[], TX[], TX[])
end
Base.eltype(w::MultiLevelWorkspace{TX}) where TX = eltype(TX)
blocksize(w::MultiLevelWorkspace{TX, bs}) where {TX, bs} = bs

function residual!(m::MultiLevelWorkspace{TX, bs}, n) where {TX, bs}
    if bs === 1
        push!(m.res_vecs, TX(undef, n))
    else
        push!(m.res_vecs, TX(undef, n, bs))
    end
end
function coarse_x!(m::MultiLevelWorkspace{TX, bs}, n) where {TX, bs}
    if bs === 1
        push!(m.coarse_xs, TX(undef, n))
    else
        push!(m.coarse_xs, TX(undef, n, bs))
    end
end
function coarse_b!(m::MultiLevelWorkspace{TX, bs}, n) where {TX, bs}
    if bs === 1
        push!(m.coarse_bs, TX(undef, n))
    else
        push!(m.coarse_bs, TX(undef, n, bs))
    end
end

abstract type CoarseSolver end
struct Pinv{T} <: CoarseSolver
    pinvA::Matrix{T}
    Pinv{T}(A) where T = new{T}(pinv(Matrix(A)))
end
Pinv(A) = Pinv{eltype(A)}(A)

(p::Pinv)(x, b) = mul!(x, p.pinvA, b)

Base.length(ml::MultiLevel) = length(ml.levels) + 1

function Base.show(io::IO, ml::MultiLevel)
    op = operator_complexity(ml)
    g = grid_complexity(ml)
    c = ml.coarse_solver
    total_nnz = nnz(ml.final_A)
    if !isempty(ml.levels)
        total_nnz += sum(nnz(level.A) for level in ml.levels)
    end
    lstr = ""
    if !isempty(ml.levels)
    for (i, level) in enumerate(ml.levels)
        lstr = lstr *
            @sprintf "   %2d   %10d   %10d [%5.2f%%]\n" i size(level.A, 1) nnz(level.A) (100 * nnz(level.A) / total_nnz)
    end
    end
    lstr = lstr *
        @sprintf "   %2d   %10d   %10d [%5.2f%%]" length(ml.levels) + 1 size(ml.final_A, 1) nnz(ml.final_A) (100 * nnz(ml.final_A) / total_nnz)

    opround = round(op, digits = 3)
    ground = round(op, digits = 3)

    str = """
    Multilevel Solver
    -----------------
    Operator Complexity: $opround
    Grid Complexity: $ground
    No. of Levels: $(length(ml))
    Coarse Solver: $c
    Level     Unknowns     NonZeros
    -----     --------     --------
    $lstr
    """
    print(io, str)
end

function operator_complexity(ml::MultiLevel)
    if !isempty(ml.levels)
        (sum(nnz(level.A) for level in ml.levels) +
                nnz(ml.final_A)) / nnz(ml.levels[1].A)
    else
        1.
    end
end

function grid_complexity(ml::MultiLevel)
    if !isempty(ml.levels)
        (sum(size(level.A, 1) for level in ml.levels) +
                size(ml.final_A, 1)) / size(ml.levels[1].A, 1)
    else
        1.
    end
end

abstract type Cycle end
struct V <: Cycle
end

struct W <: Cycle
end

struct F <: Cycle
end

"""
    solve(ml::MultiLevel, b::AbstractArray, cycle, kwargs...)

Execute multigrid cycling.

Arguments
=========
* ml::MultiLevel - the multigrid hierarchy
* b::Vector - the right hand side
* cycle -  multigird cycle to execute at each iteration. Defaults to AMG.V()

Keyword Arguments
=================
* reltol::Float64 - relative tolerance criteria for convergence, the absolute tolerance will be `reltol * norm(b)`
* abstol::Float64 - absolute tolerance criteria for convergence
* maxiter::Int64 - maximum number of iterations to execute
* verbose::Bool - display residual at each iteration
* log::Bool - return vector of residuals along with solution

"""
function solve(ml::MultiLevel, b::AbstractArray, args...; kwargs...)
    n = length(ml) == 1 ? size(ml.final_A, 1) : size(ml.levels[1].A, 1)
    V = promote_type(eltype(ml.workspace), eltype(b))
    x = zeros(V, size(b))
    return solve!(x, ml, b, args...; kwargs...)
end
function solve!(x, ml::MultiLevel, b::AbstractArray{T},
                                    cycle::Cycle = V();
                                    maxiter::Int = 100,
                                    abstol::Real = zero(real(eltype(b))),
                                    reltol::Real = sqrt(eps(real(eltype(b)))),
                                    verbose::Bool = false,
                                    log::Bool = false,
                                    calculate_residual = true) where {T}

    A = length(ml) == 1 ? ml.final_A : ml.levels[1].A
    V = promote_type(eltype(A), eltype(b))
    log && (residuals = Vector{V}())
    normres = normb = norm(b)
    if normb != 0
        abstol = max(reltol * normb, abstol)
    end
    log && push!(residuals, normb)

    res = ml.workspace.res_vecs[1]
    itr = lvl = 1
    while itr <= maxiter && (!calculate_residual || normres > abstol)
        if length(ml) == 1
            ml.coarse_solver(x, b)
        else
            __solve!(x, ml, cycle, b, lvl)
        end
        if calculate_residual
            mul!(res, A, x)
            reshape(res, size(b)) .= b .- reshape(res, size(b))
            normres = norm(res)
            log && push!(residuals, normres)
        end
        itr += 1
    end

    # @show residuals
    log ? (x, residuals) : x
end

function __solve_next!(x, ml, cycle::V, b, lvl)
    __solve!(x, ml, cycle, b, lvl)
end

function __solve_next!(x, ml, cycle::W, b, lvl)
    __solve!(x, ml, cycle, b, lvl)
    __solve!(x, ml, cycle, b, lvl)
end

function __solve_next!(x, ml, cycle::F, b, lvl)
    __solve!(x, ml, cycle, b, lvl)
    __solve!(x, ml, V(), b, lvl)
end

function __solve!(x, ml, cycle::Cycle, b, lvl)
    A = ml.levels[lvl].A
    ml.presmoother(A, x, b)

    res = ml.workspace.res_vecs[lvl]
    mul!(res, A, x)
    reshape(res, size(b)) .= b .- reshape(res, size(b))

    coarse_b = ml.workspace.coarse_bs[lvl]
    mul!(coarse_b, ml.levels[lvl].R, res)

    coarse_x = ml.workspace.coarse_xs[lvl]
    coarse_x .= 0
    if lvl == length(ml.levels)
        ml.coarse_solver(coarse_x, coarse_b)
    else
        coarse_x = __solve_next!(coarse_x, ml, cycle, coarse_b, lvl + 1)
    end

    mul!(res, ml.levels[lvl].P, coarse_x)
    x .+= res

    ml.postsmoother(A, x, b)

    x
end


text:
using Distributions
using Flux.CUDA

"""
	Implements all needed parts for constructing f-AnoGAN model
		(https://www.sciencedirect.com/science/article/pii/S1361841518302640).
"""

struct fAnoGAN
	prior
	generator
	discriminator
	encoder
end

Flux.@functor fAnoGAN

function Flux.trainable(model::fAnoGAN)
	(generator = model.generator, discriminator = model.discriminator, encoder = model.encoder)
end

function izi(m::fAnoGAN, x)
	z = m.encoder(x)
	x̂ = m.generator(z)
	f = m.discriminator[1:end-1]
	fx = f(x)
	fx̂ = f(x̂)
	return x̂, fx, fx̂
end

"""
	function izif_loss(m::fAnoGAN, x, κ)
loss L_{izi_f} = 1/n * ||x-D(x)||^2 + κ / n_g ||f(x) - f(G(x))||^2
"""
function izif_loss(m::fAnoGAN, x, κ)
	x̂, fx, fx̂ = izi(m, x)
	loss = Flux.mse(x,x̂) + κ .* Flux.mse(fx,fx̂)
	return loss
end

"""
	function anomaly_score(m::fAnoGAN, x, κ=1f0)
Anomaly score computed according to formula from fAnoGAN paper.
Version izi_f (image-z-image)
		A(x) = 1/n * ||x-D(x)||^2 + κ / n_g ||f(x) - f(G(x))||^2
"""
function anomaly_score(m::fAnoGAN, x, κ=1f0)
	x̂, fx, fx̂ = izi(m, x)
	xx̂= vec(Flux.sum((x̂ .- x).^2, dims=[1,2,3]))
	fxfx̂ = vec(Flux.sum((fx̂ .- fx).^2, dims=1))
	return xx̂ .+ κ .* fxfx̂
end

"""
	function anomaly_score_gpu(m::fAnoGAN, real_input, κ = 1.0; batch_size=64, to_testmode::Bool=true)
Function which compute anomaly score on GPU.
"""
function anomaly_score_gpu(m::fAnoGAN, real_input, κ = 1f0; batch_size=64, to_testmode::Bool=true)
	real_input = Flux.Data.DataLoader(real_input, batchsize=batch_size)
	(to_testmode == true) ? Flux.testmode!(m) : nothing
	m = m |> gpu
	scores = Array{Float32}([])
	for x in real_input
		x = x |> gpu
		score = anomaly_score(m, x, κ)
		scores = cat(scores, score |> cpu, dims=1)
	end
	(to_testmode == true) ? Flux.testmode!(m, false) : nothing
	return scores
end

# loss functions for Wassersten GAN
wdloss(gen, dis, x, z) = -(Flux.mean(dis(x)) - Flux.mean(dis(gen(z))))
wdloss(model::fAnoGAN, x) = wdloss(model.generator, model.discriminator, x, rand(model.prior, size(x, ndims(x))))
wdloss(model::fAnoGAN, x, z) = wdloss(model.generator, model.discriminator, x, z)

wgloss(gen, dis, z) = - Flux.mean(dis(gen(z)))
wgloss(model::fAnoGAN, z) = wgloss(model.generator, model.discriminator, z)

function fanogan_constructor(
		;idim=(2,2,1),
		zdim::Int=1,
		activation="relu",
		hdim=1024,
		kernelsizes=(1,1),
		channels=(1,1),
		scalings=(1,1),
		init_seed=nothing,
		batchnorm=false,
		kwargs...
	)

	(init_seed !== nothing) ? Random.seed!(init_seed) : nothing

	generator = conv_decoder(idim, zdim, reverse(kernelsizes), reverse(channels),
		reverse(scalings); activation=activation, batchnorm=batchnorm)

	discriminator = Chain(conv_encoder(idim, 1, kernelsizes, channels, scalings,
			activation=activation, batchnorm=batchnorm)...)

	encoder_ =  Chain(conv_encoder(idim, zdim, kernelsizes, channels, scalings,
	activation=activation, batchnorm=batchnorm)..., x->tanh.(x))

	(init_seed !== nothing) ? Random.seed!() : nothing

	model = fAnoGAN(MvNormal(zdim, 1f0), generator, discriminator, encoder_)
end

"""
	function StatsBase.fit!(model::fAnoGAN, data::Tuple, params)

General function for fitting of fAnoGAN.
	model::fAnoGAN
	data::Tupe   		... ((x_train, y_train), (x_valid, y_valid), (x_test, y_test))
	params::NamedTuple  ... all needed parameters for training
		kappa  		 		... scaling coeficent for izi_f loss function
		weight_clip  		... weight clipping parameter (https://arxiv.org/abs/1701.07875)
		lr_gan		 		... learning rate for WGAN
		lr_enc		 		... learning rate for encoder (image-z-iamge_f)
		batch_size::Int 	... batch size
		check_every::Int  	... number of iterations between EarlyStopping evaluation check
		patience::Int		... patience before stoping training
		iters::Int     		... max number of iteration
		mtt::Int			... max train time (if training is too slow)
		n_critic::Int   	... number of discriminator updater per one generator update
		usegpu::Bool    	... to use GPU or not

	Example:
		params = (kappa = 1f0, weight_clip=0.1, lr_gan = 0.001, lr_enc = 0.001, batch_size=128,
			patience=10, check_every=30, iters=10000, mtt=82800, n_critic=1, usegpu=true)

"""
function StatsBase.fit!(model::fAnoGAN, data::Tuple, params::NamedTuple)

	model = params.usegpu ? model |>gpu : model |>cpu
	best_model = deepcopy(model)
	history = MVHistory()
	train_loader, val_loader = prepare_dataloaders(data, params)
	val_batches = length(val_loader)
	"""
		WGAN training
	"""

	opt_G = RMSProp(params.lr_gan)
	opt_D = RMSProp(params.lr_gan)
	opt_E = ADAM(params.lr_enc)
	ps_E = Flux.params(model.encoder)
	ps_D = Flux.params(model.discriminator) #model.gan.discriminator
	ps_G = Flux.params(model.generator) # model.gan.generator
	# no early stopping here <– bahavior of loss function will not allow it
	progress = Progress(length(train_loader))
	for (iter, X) in enumerate(train_loader)
		X = params.usegpu ? gpu(getobs(X)) : getobs(X)
		last_loss_D = 0
		for n = 1:params.n_critic
			z = rand(model.prior, size(X, ndims(X)))
			z = params.usegpu ? z|>gpu : z
			loss_D, back_D = Flux.pullback(ps_D) do
				wdloss(model, X, z)
			end
			grad_D = back_D(1f0)
			Flux.Optimise.update!(opt_D, ps_D, grad_D)
			clip_weights!(ps_D, params.weight_clip)
			last_loss_D = loss_D
		end
		z = rand(model.prior, size(X, ndims(X)))
		z = params.usegpu ? z|>gpu : z
		loss_G, back_G = Flux.pullback(ps_G) do
				wgloss(model, z)
		end
		grad_G = back_G(1f0)
		Flux.Optimise.update!(opt_G, ps_G, grad_G)
		push!(history, :loss_G, iter, loss_G)
		push!(history, :loss_D, iter, last_loss_D)
		next!(progress; showvalues=[
			(:iters, "$(iter)/$(params.iters)"),
			(:loss_G, loss_G),
			(:loss_D, last_loss_D)
			])
		if mod(iter, params.check_every) == 0
			val_loss_G = 0
			val_loss_D = 0
			Flux.testmode!(model)
			for X_val in val_loader
				X_val = params.usegpu ? X_val |> gpu : X_val
				z = rand(model.prior, size(X_val, ndims(X_val)))
				z = params.usegpu ? z|>gpu : z
				val_loss_G += wgloss(model, z)
				val_loss_D += wdloss(model, X_val, z)
			end
			Flux.testmode!(model, false)
			push!(history, :validation_D, iter, val_loss_D/val_batches)
			push!(history, :validation_G, iter, val_loss_G/val_batches)
		end
	end

	"""
		Encoder training
	"""
	# possible memory problems
	CUDA.reclaim() # somehow GPU will run out of memory. This helps to clear cached memory and continue with it

	best_model = deepcopy(model)
	# early stopping
	best_val_loss = Inf
	patience_ = deepcopy(params.patience)
	progress = Progress(length(train_loader))

	start_time = time()
	for (iter, X) in enumerate(train_loader)
		X = params.usegpu ? gpu(getobs(X)) : getobs(X)
		loss_E, back_E = Flux.pullback(ps_E) do
			izif_loss(model, X, params.kappa)
		end
		grad_E = back_E(1f0)
		Flux.Optimise.update!(opt_E, ps_E, grad_E)

		push!(history, :loss_E, iter, loss_E)
		next!(progress; showvalues=[
			(:iters, "$(iter)/$(params.iters)"),
			(:loss_izif, loss_E)
			])

		if mod(iter, params.check_every) == 0
			val_loss_E = 0
			Flux.testmode!(model)
			for X_val in val_loader
				X_val = params.usegpu ? X_val |> gpu : X_val
				val_loss_E += izif_loss(model, X_val, params.kappa)
			end
			Flux.testmode!(model, false)
			push!(history, :validation_E, iter, val_loss_E/val_batches)
			if val_loss_E < best_val_loss
				best_val_loss = val_loss_E
				patience_ = params.patience
				best_model = deepcopy(model)
			else
				patience_ -= 1
				if patience_ == 0
					@info "Stopped training of Autoencoder after $(iter) iterations"
					break
				end
			end
		end
		if time() - start_time > params.mtt # stop early if time is running out
			best_model = deepcopy(model)
			@info "Stopped training after $(i) iterations, $((time() - start_time)/3600) hours."
			break
		end
	end
	return best_model, history, sum(map(p->length(p), Flux.params(model)))
end


text:
using StatsBase
using Base.Test
import Base: maxabs
import StatsBase: norepeat, randi

srand(1234)

#### randi

n = 10^5

x = [randi(10) for i = 1:n]
@test isa(x, Vector{Int})
@test extrema(x) == (1, 10)
@test_approx_eq_eps proportions(x, 1:10) fill(0.1, 10) 5.0e-3

x = [randi(3, 12) for i = 1:n]
@test isa(x, Vector{Int})
@test extrema(x) == (3, 12)
@test_approx_eq_eps proportions(x, 3:12) fill(0.1, 10) 5.0e-3

#### sample with replacement

function check_sample_wrep(a::AbstractArray, vrgn, ptol::Real; ordered::Bool=false)
    vmin, vmax = vrgn
    (amin, amax) = extrema(a)
    @test vmin <= amin <= amax <= vmax
    n = vmax - vmin + 1
    p0 = fill(1/n, n)
    if ordered
        @test issorted(a)
        if ptol > 0
            @test_approx_eq_eps proportions(a, vmin:vmax) p0 ptol
        end
    else
        @test !issorted(a)
        ncols = size(a,2)
        if ncols == 1
            @test_approx_eq_eps proportions(a, vmin:vmax) p0 ptol
        else
            for j = 1:ncols
                aj = view(a, :, j)
                @test_approx_eq_eps proportions(aj, vmin:vmax) p0 ptol
            end
        end
    end
end

import StatsBase: direct_sample!

a = direct_sample!(1:10, zeros(Int, n, 3))
check_sample_wrep(a, (1, 10), 5.0e-3; ordered=false)

a = direct_sample!(3:12, zeros(Int, n, 3))
check_sample_wrep(a, (3, 12), 5.0e-3; ordered=false)

a = direct_sample!([11:20;], zeros(Int, n, 3))
check_sample_wrep(a, (11, 20), 5.0e-3; ordered=false)

a = sample(3:12, n)
check_sample_wrep(a, (3, 12), 5.0e-3; ordered=false)

a = sample(3:12, n; ordered=true)
check_sample_wrep(a, (3, 12), 5.0e-3; ordered=true)

a = sample(3:12, 10; ordered=true)
check_sample_wrep(a, (3, 12), 0; ordered=true)

#### sampling pairs

srand(1);

@test samplepair(2)  === (1, 2)
@test samplepair(10) === (7, 3)

@test samplepair([3, 4, 2, 6, 8]) === (4, 3)
@test samplepair([1, 2])          === (1, 2)

#### sample without replacement

function check_sample_norep(a::AbstractArray, vrgn, ptol::Real; ordered::Bool=false)
    # each column of a for one run

    vmin, vmax = vrgn
    (amin, amax) = extrema(a)
    @test vmin <= amin <= amax <= vmax
    n = vmax - vmin + 1

    for j = 1:size(a,2)
        aj = view(a,:,j)
        @assert norepeat(aj)
        if ordered
            @assert issorted(aj)
        end
    end

    if ptol > 0
        p0 = fill(1/n, n)
        if ordered
            @test_approx_eq_eps proportions(a, vmin:vmax) p0 ptol
        else
            b = transpose(a)
            for j = 1:size(b,2)
                bj = view(b,:,j)
                @test_approx_eq_eps proportions(bj, vmin:vmax) p0 ptol
            end
        end
    end
end

import StatsBase: knuths_sample!, fisher_yates_sample!, self_avoid_sample!
import StatsBase: seqsample_a!, seqsample_c!

a = zeros(Int, 5, n)
for j = 1:size(a,2)
    knuths_sample!(3:12, view(a,:,j))
end
check_sample_norep(a, (3, 12), 5.0e-3; ordered=false)

a = zeros(Int, 5, n)
for j = 1:size(a,2)
    fisher_yates_sample!(3:12, view(a,:,j))
end
check_sample_norep(a, (3, 12), 5.0e-3; ordered=false)

a = zeros(Int, 5, n)
for j = 1:size(a,2)
    self_avoid_sample!(3:12, view(a,:,j))
end
check_sample_norep(a, (3, 12), 5.0e-3; ordered=false)

a = zeros(Int, 5, n)
for j = 1:size(a,2)
    seqsample_a!(3:12, view(a,:,j))
end
check_sample_norep(a, (3, 12), 5.0e-3; ordered=true)

a = zeros(Int, 5, n)
for j = 1:size(a,2)
    seqsample_c!(3:12, view(a,:,j))
end
check_sample_norep(a, (3, 12), 5.0e-3; ordered=true)

a = sample(3:12, 5; replace=false)
check_sample_norep(a, (3, 12), 0; ordered=false)

a = sample(3:12, 5; replace=false, ordered=true)
check_sample_norep(a, (3, 12), 0; ordered=true)


text:
# We create one trace per derived function, which tracks all the salsa computations called
# by the current derived function, so that they can be recorded as input dependencies. We
# create one instance per derived function to limit contention in multithreaded code, and
# protect insertions with a lock to allow concurrency within the body of the function.
#
# TODO: For best parallel performance, we might consider using Thread-local storage (not to
# be confused with Task-local storage) to collect dependencies, so that we can remove the
# lock entirely: By storing one vector per OS-thread, and aggregating at the end. The
# potential downside is losing some ordering information if part of the function is serial
# and part is parallel.
#
# NOTE: This is a mutable struct to allow us to modify the call_stack. This is fine because
# it's not going to be isbits anyway due to all of the Vectors, Sets, and Locks, it
# contains, and we are going to keep these in the trace pool (see explanation below).
mutable struct TraceOfDependencyKeys
    # Performance Optimization: De-duplicating Derived Function Traces:
    #   We want to maintain the order in which the dependencies were encountered (to
    #   preserve correctness when checking dependencies for invalidation), but we also
    #   don't want to track exact duplicate dependencies, so we use the seen_deps to filter.
    # TODO: We may want to consider using a linked list instead of a vector to prevent
    # array growth events from causing contention when spawning parallel derived functions.
    ordered_deps::Vector{DependencyKey}
    seen_deps::Set{DependencyKey}

    # We lock around modifications to the list of dependencies, since a derived function may
    # spawn multiple threads, which may themselves call derived functions. This doesn't add
    # too much contention, because it's only blocking around the spawning of the immediate
    # children of node.
    lock::Base.ReentrantLock

    # The trace of Salsa functions that represents the current computation in progress.
    # Since derived functions may spawn threads, this represents a pure stack trace back to
    # the root call into Salsa, which potentially traverses across threads.  It's an
    # immutable structure (linked list) to make it thread-safe.
    call_stack::Union{Nothing,SalsaStackFrame}

    # This is set to false during recursive "still_valid" checks to avoid unused dependency
    # tracking, to save time and allocs.
    should_trace::Bool

    # We always create a new, empty TraceOfDependencyKeys for each derived function, since
    # we're only tracing the immediate dependencies of that function.
    function TraceOfDependencyKeys()
        # Pre-allocate all traces to be non-empty, to minimize allocations at runtime.
        # These traces are all constructed once ahead of time in the per-thread trace pools,
        # so this initialization is only done once during Module __init__().
        #
        # TODO: Tune this. Too big wastes RAM (though, it's fixed cost up front).
        #       Current size, 30, adds about 1MiB, which seems not bad.
        N = 30
        return new(
            sizehint!(Vector{DependencyKey}(), N),
            sizehint!(Set{DependencyKey}(), N),
            Base.ReentrantLock(),
            nothing,
            true,
        )
    end
end

function push_key!(trace::TraceOfDependencyKeys, depkey)
    # (Don't need to lock around this since it can never be modified in parallel.)
    if !trace.should_trace
        return
    end
    @lock trace.lock begin
        # Performance Optimization: De-duplicating Derived Function Traces
        if depkey ∉ trace.seen_deps
            push!(trace.ordered_deps, depkey)
            push!(trace.seen_deps, depkey)
        end
    end
    return nothing
end

function collect_trace(trace::TraceOfDependencyKeys)
    # NOTE: We do not need to lock around this access, because the user's function is
    # guaranteed to have finished by this point, and it is _not supported_ to register
    # more derived functions on the current trace after that. We should have exclusive
    # access to the trace at this point, and it's about to be handed back to the trace pool!

    # Derived functions that spawn tasks currently *must* wait until all tasks have
    # completed before returning. Calling new derived functions from spawned tasks after
    # returning is an error! (But unfortunately we have no way to detect/prevent it.)

    return copy(trace.ordered_deps)
end


text:
using Test
using FlightMechanics
using FlightMechanics.Models

# TEST check_state_aerostate_env_coherence
# STATE
pos = PositionEarth(0.0, 0.0, 0.0)
att = Attitude(0.0, 0.0, 0.0)
vel = [100.0, 0.0, 0.0]
ang_vel = [0.0, 0.0, 0.0]
accel = [0.0, 0.0, 0.0]
ang_accel = [0.0, 0.0, 0.0]
state = State(pos, att, vel, ang_vel, accel, ang_accel)

# ENVIRONMENT
atm = AtmosphereISA(pos)
wind = ConstantWind(0.0, 10.0, 0.0)
grav = EarthConstantGravity()
env = Environment(atm, wind, grav)

# AEROSTATE
tas = 110
alpha, beta = 0.0, 0.0
aerost = AeroState(tas, alpha, beta, 0.0)

@test isequal(check_state_aerostate_env_coherence(state, aerost, env), true)

# TEST generate_state_aerostate for the same case
state, aerost = generate_state_aerostate(pos, att, tas, alpha, beta, env,
                                         ang_vel, accel, ang_accel)
@test isequal(check_state_aerostate_env_coherence(state, aerost, env), true)
@test isapprox(get_horizon_velocity(state), [100, 0.0, 0.0])

# Opposite direction
# AEROSTATE
tas = 90
alpha, beta = 0.0, 0.0
aerost_exp = AeroState(tas, alpha, beta, 0.0)

att = Attitude(pi, 0.0, 0.0)
# TEST generate_state_aerostate
state, aerost = generate_state_aerostate(pos, att, tas, alpha, beta, env, ang_vel, accel,
                                ang_accel)
@test isequal(check_state_aerostate_env_coherence(state, aerost, env), true)
@test isapprox(get_horizon_velocity(state), [-100, 0.0, 0.0])

# Opposite direction check for fail
# AEROSTATE
tas = 90
alpha, beta = 0.0, 0.0
aerost_exp = AeroState(tas, alpha, beta, 0.0)

att = Attitude(pi, 0.0, 0.0)
# TEST generate_state_aerostate
state, aerost = generate_state_aerostate(pos, att, tas, alpha, beta, env, ang_vel, accel,
                                ang_accel)
state_wrong = State(get_position(state),
                    get_attitude(state),
                    get_body_velocity(state) .+ 0.5,
                    get_body_ang_velocity(state),
                    get_body_accel(state),
                    get_body_ang_accel(state)
                    )
@test isequal(check_state_aerostate_env_coherence(state_wrong, aerost, env),
              false)


text:
"""
    phase(loc::String, type::DataType, digit::Int)
Takes a String of numbers returns a typed locus appropriate for PopGen.jl as used in the
`genepop` and `csv` file parsers. Use `type` to specify output type (`Int8` or `Int16`),
and `digit` to specify the number of digits/characters used per allele in a locus.

## Examples
`ph_locus = phase("128114", Int16, 3)`

`map(i -> phase(i, Int16, 3), ["112131", "211112", "001003", "516500"])`

`[phase(i, Int8, 2) for i in ["0101", "0103", "0202", "0103"]]`
"""
function phase(loc::String, type::DataType, digit::Int)
    phased = [join(k) for k in Iterators.partition(loc, digit)]
    typed = parse.(type, phased) |> sort!
    tupled = Tuple(typed)
    iszero(tupled |> collect) && return missing
    return tupled
end

### Variant Call Format parsing ###

"""
    bcf(infile::String)
Load a BCF file into memory as a PopObj object. Population and [optional]
location information need to be provided separately.
- `infile` : path to BCF file
"""
function bcf(infile::String)
    bcf_file = BCF.Reader(open(infile, "r"))

    # get sample names from header
    sample_names = header(bcf_file).sampleID

    # fill in pop/lat/long with missing
    population = fill("missing", length(sample_names))
    loc_xy = Vector{Union{Missing, Float32}}()
    append!(loc_xy, fill(missing, length(sample_names)))

    # get loci names
    locinames = Vector{String}()

    ## array of genotypes
    locus_array = Vector{Vector{Union{Missing, Tuple{Vararg}}}}()

    # get genotypes
    for record in bcf_file
        # fix locus names to be syntax-safe
        chr_safe = replace(BCF.chrom(record), r"\.|\-|\=|\/" => "_")
        chr_safer = replace(chr_safe, "|" => "_")
        pos = BCF.pos(record) |> string
        push!(locinames, chr_safer*"_"*pos)

        # get the genotype information
        geno_raw = [split(i, ('/', '|')) for i in BCF.genotype(record, :, "GT")] |> sort

        # change missing data "." to "-1"
        geno_corr_miss = map(i -> replace(i, "." => "-1"), geno_raw)

        # convert everything to an integer
        geno_int = map(i -> parse.(Int8, i), geno_corr_miss)

        # add 1 to shift genos so 0 is 1 and -1 is 0 etc.
        geno_shift = map(i -> i .+ Int8(1), geno_int)
        geno_final = [replace(i, 0 => missing) for i in geno_shift]
        geno_tuple = [Tuple(i) for i in geno_final]
        push!(locus_array, geno_tuple)
    end
    close(vcf_file)

    # intelligently scan for ploidy
    ploidy = Vector{Int8}()
    @inbounds for i in 1:length(sample_names)
        @inbounds for j in 1:length(locinames)
            locus_array[j][i] === missing && continue   # if missing, go to next locus
            ploid = length(locus_array[j][i])   # if not, get the # of alleles
            push!(ploidy, ploid)    # push that to the ploidy vector
            break   # break out of the loop and begin next sample
        end
    end

    @info "\n$(abspath(infile))
$(length(sample_names)) samples detected
population info must be added <---
$(length(locinames)) loci detected"

    # create loci dataframe
    loci_df = DataFrame([j => k for (j,k) in zip(Symbol.(locinames), locus_array)])

    # create samples df
    samples_df = DataFrame(
        name = sample_names,
        population = population,
        ploidy = ploidy,
        latitude = loc_xy,
        longitude = loc_xy
    )
    PopObj(samples_df, loci_df)
end

### delimited parsing ###

"""
    delimited(infile::String; delim::Union{Char,String,Regex}, digits::Int64 = 3)
Load a delimited-type file into memory as a PopObj object. *There should be no empty cells
in your file*
### Arguments
- `infile` : path to file
- `delim` : delimiter characters. By default auto-parses comma, space, or tab (use ONLY one!)
- `digits` : number of digits denoting each allele (default = 3)
- `marker` : "snp" (default) or "msat"

### File formatting:
- First row is column names in this order:
    1. name
    2. population
    3. longitude
    4. latitude
    5. locus_1_name
    6. locus_2_name
    7. etc...

#### Formatting location data
If location data is missing for a sample (which is ok!), make sure the value is written
as *0*, otherwise there will be transcription errors!

## Example
`lizardsCA = Read.delimited("CA_lizards.csv", digits = 3);`

### Formatting example:
name,population,long,lat,Locus1,Locus2,Locus3   \n
sierra_01,mountain,11.11,-22.22,001001,002002,001001   \n
sierra_02,mountain,11.12,-22.21,001001,001001,001002   \n
snbarb_03,coast,0,0,001001,001001,001002 \n
snbarb_02,coast,11.14,-22.24,001001,001001,001001 \n
snbarb_03,coast,11.15,0,001002,001001,001001 \n
"""
function delimited(
    infile::String;
    delim::Union{Char,String,Regex} = r"\,|\t|\s",
    digits::Int = 3,
    marker = "snp",
)

    # instantiate empty vectors
    indnames = Vector{String}()
    popid = Vector{Union{String,Int32}}()
    ploidy = Vector{Int8}()
    locx = Vector{Union{Missing, Float32}}()
    locy = Vector{Union{Missing, Float32}}()
    sample_geno_array = Array{Vector{Union{Missing, Tuple}},1}()

    # establish types for genotypes
    if lowercase(marker) == "snp"
        geno_type = Int8
    else
        geno_type = Int16
    end

    file = open(infile, "r")

    # get header information
    col_names = map(join, split(readline(file), delim))
    replace!(col_names, "." => "_")
    locinames = col_names[5:end]

    @inbounds for ln in eachline(file)
        tmp = map(join, split(ln, delim))
        # just in case -9 = missing
        replace!(tmp, "-9" => "0"^digits)
        push!(indnames, popfirst!(tmp))
        push!(popid, popfirst!(tmp))
        push!(locx, parse(Float32, popfirst!(tmp)))
        push!(locy, parse(Float32, popfirst!(tmp)))

        # phase the loci
        phasedloci = map(i -> phase(i, geno_type, digits), tmp)
        push!(sample_geno_array, phasedloci)

        # get the ploidy via # of alleles in the first non-missing locus
        sample_ploidy = (skipmissing(phasedloci) |> collect)[1]  |> length
        push!(ploidy, sample_ploidy)

    end
    close(file)

    @info "\n$(abspath(infile))
$(length(indnames)) samples detected
$(length(unique(popid))) populations detected
$(length(locinames)) loci detected"

    # convert the entire thing into a matrix
    geno_matrix = hcat(sample_geno_array...)

    # convert is back into an array of arrays, transposed across loci
    # instead of samples
    loci_array = collect.(eachrow(geno_matrix))

    # replace 0.0 with missing in location
    replace!(locx, 0.0 => missing)
    replace!(locy, 0.0 => missing)

    # genotype DataFrame
    loci_df = DataFrame([j => k for (j,k) in zip(Symbol.(locinames), loci_array)])

    samples_df = DataFrame(
        name = indnames,
        population = popid,
        ploidy = ploidy,
        longitude = locx,
        latitude = locy,
    )

    return PopObj(samples_df, loci_df)
end

const csv = delimited

### GenePop parsing ###
"""
    genepop(infile::String; digits::Int = 3, popsep::String = "POP", marker::String = "snp")
Load a Genepop format file into memory as a PopObj object.
### Arguments
- `infile` : path to Genepop file
- `digits` : number of digits denoting each allele
- `popsep` : word that separates populations in `infile` (default: "POP")
- `marker` : "snp" (default) or "msat"
### File must follow standard Genepop formatting:
- First line is a comment (and skipped)
- Loci are listed after first line as one-per-line without commas or in single comma-separated row
- A line with a particular keyword (default "POP") must delimit populations
- Sample name is followed by a comma and space or tab (", ")
- File is tab or space delimted

## Example

 `waspsNY = genepop("wasp_hive.gen", digits = 3, popsep = "POP")`

### Genepop file example:
Wasp populations in New York \n
Locus1\n
Locus2\n
Locus3\n
POP\n
Oneida_01,  250230 564568 110100\n
Oneida_02,  252238 568558 100120\n
Oneida_03,  254230 564558 090100\n
POP\n
Newcomb_01,  254230 564558 080100\n
Newcomb_02,  000230 564558 090080\n
Newcomb_03,  254230 000000 090100\n
Newcomb_04,  254230 564000 090120\n
"""
function genepop(
    infile::String;
    digits::Int = 3,
    popsep::String = "POP",
    marker::String = "snp",
)
    if lowercase(marker) == "snp"
        geno_type = Int8
    else
        geno_type = Int16
    end

    gpop = split(open(readlines, infile)[2:end], popsep)

    @info "\n$(abspath(infile))
$(sum(length.(gpop[2:end]))) samples detected
$(length(gpop)-1) populations detected using \"$popsep\" seperator
$(length(gpop[1])) loci detected"

    if length(gpop[1]) == 1     # loci horizontally stacked.
        locinames = strip.(split(gpop[1] |> join, ",") |> Vector{String})
        replace!(locinames, "." => "_")
    else                        # loci vertically stacked
        locinames = replace(gpop[1], "." => "_")
    end

    sample_meta_names = [:name, :population, :ploidy, :latitude, :longitude]
    sample_meta_array = [
        Vector{String}(),
        Vector{Union{String,Int}}(),
        Vector{Int8}(),
        Vector{Union{Missing,Float32}}(),
        Vector{Union{Missing, Float32}}()
    ]

    sample_geno_array = Vector{Vector{Union{Missing, Tuple}}}()

    for i = 2:length(gpop)
        for j in gpop[i]
            samplerow = map(join, split(strip(j), r"\,\t|\,\s|\s|\t"))

            # just in case missing genotypes are coded as -9 (msats)
            replace!(samplerow, "-9" => "0"^digits)

            # sample name and population
            push!(sample_meta_array[1], popfirst!(samplerow))
            push!(sample_meta_array[2], i-1)

            # phase the loci into tuples
            converted_row = map(k -> phase(k, geno_type, digits), samplerow)
            push!(sample_geno_array, converted_row)

            # get the ploidy via # of alleles in the first non-missing locus
            push!(sample_meta_array[3], length((skipmissing(converted_row) |> collect)[1]))

            # long and lat
            push!(sample_meta_array[4], missing) ; push!(sample_meta_array[5], missing)
        end
    end

    # convert the entire thing into a matrix
    geno_matrix = hcat(sample_geno_array...)

    # convert back into an array of arrays of loci instead of samples
    loci_array = collect.(eachrow(geno_matrix))

    # create the two dataframes necessary for a PopObj
    samples = DataFrame([j => k for (j,k) in zip(sample_meta_names, sample_meta_array)])
    loci = DataFrame([j => k for (j,k) in zip(Symbol.(locinames), loci_array)])

    return PopObj(samples, loci)
end

#= Alternative line-by-line genepop reader

# Alternative line-by-line reader
## NOT EXPORTED

function gpop2(
    infile::String;
    digits::Int = 3,
    popsep::String = "POP",
    marker::String = "snp",
)

    if lowercase(marker) == "snp"
        geno_type = Int8
    else
        geno_type = Int16
    end

    sample_meta_names = [:name, :population, :ploidy, :latitude, :longitude]
    sample_meta_array = [
        Vector{String}(),
        Vector{Union{String,Int}}(),
        Vector{Int8}(),
        Vector{Union{Missing,Float32}}(),
        Vector{Union{Missing, Float32}}()
    ]

    locinames = Vector{String}()
    sample_geno_array = Vector{Vector{Union{Missing, Tuple}}}()
    popcount = 1

    gpop = open(infile, "r")

    # skip first line
    readline(gpop)

    # test for loci horizontal/vertical formatting
    locus_name = readline(gpop)
    if occursin(",", locus_name) == true
        append!(locinames, strip.(split(locus_name |> join, ",")))
        replace!(locinames, "." => "_")
    else
        while locus_name != popsep
            push!(locinames, locus_name)
            locus_name = readline(gpop)
        end
    end

    for ln in readlines(gpop)
        if ln == popsep
            popcount += 1
            continue
        else
            samplerow = map(join, split(strip(ln), r"\,\t|\,\s|\s|\t"))

            # just in case missing genotypes are coded as -9 (msats)
            replace!(samplerow, "-9" => "0"^digits)

            # sample name and population
            push!(sample_meta_array[1], popfirst!(samplerow))
            push!(sample_meta_array[2], popcount)

            # phase the loci into tuples
            converted_row = map(k -> phase(k, geno_type, digits), samplerow)
            push!(sample_geno_array, converted_row)

            # get the ploidy via # of alleles in the first non-missing locus
            push!(sample_meta_array[3], length((skipmissing(converted_row) |> collect)[1]))

            # long and lat
            push!(sample_meta_array[4], missing) ; push!(sample_meta_array[5], missing)
        end
    end

    @info "\n$(abspath(infile))
$(length(sample_meta_array[1])) samples detected
$(popcount) populations detected using \"$popsep\" seperator
$(length(locinames)) loci detected"
    # convert the entire thing into a matrix
    geno_matrix = hcat(sample_geno_array...)

    # convert back into an array of arrays of loci instead of samples
    loci_array = collect.(eachrow(geno_matrix))

    # create the two dataframes necessary for a PopObj
    samples = DataFrame([j => k for (j,k) in zip(sample_meta_names, sample_meta_array)])
    loci = DataFrame([j => k for (j,k) in zip(Symbol.(locinames), loci_array)])

    return PopObj(samples, loci)
end
=#

### VCF parsing ###

"""
    vcf(infile::String)
Load a VCF file into memory as a PopObj object. Population and [optional]
location information need to be provided separately.
- `infile` : path to VCF file
"""
function vcf(infile::String)
    vcf_file = VCF.Reader(open(infile, "r"))

    # get sample names from header
    sample_names = header(vcf_file).sampleID

    # fill in pop/lat/long with missing
    population = fill("missing", length(sample_names))
    loc_xy = Vector{Union{Missing, Float32}}()
    append!(loc_xy, fill(missing, length(sample_names)))

    # get loci names
    locinames = Vector{String}()

    ## array of genotypes
    locus_array = Vector{Vector{Union{Missing, Tuple{Vararg}}}}()

    # get genotypes
    for record in vcf_file
        # fix locus names to be syntax-safe
        chr_safe = replace(VCF.chrom(record), r"\.|\-|\=|\/" => "_")
        chr_safer = replace(chr_safe, "|" => "_")
        pos = VCF.pos(record) |> string
        push!(locinames, chr_safer*"_"*pos)

        # get the genotype information
        geno_raw = [split(i, ('/', '|')) for i in VCF.genotype(record, :, "GT")] |> sort

        # change missing data "." to "-1"
        geno_corr_miss = map(i -> replace(i, "." => "-1"), geno_raw)

        # convert everything to an integer
        geno_int = map(i -> parse.(Int8, i), geno_corr_miss)

        # add 1 to shift genos so 0 is 1 and -1 is 0 etc.
        geno_shift = map(i -> i .+ Int8(1), geno_int)
        geno_final = [replace(i, 0 => missing) for i in geno_shift]
        geno_tuple = [Tuple(i) for i in geno_final]
        push!(locus_array, geno_tuple)
    end
    close(vcf_file)

    # intelligently scan for ploidy
    ploidy = Vector{Int8}()
    @inbounds for i in 1:length(sample_names)
        @inbounds for j in 1:length(locinames)
            locus_array[j][i] === missing && continue   # if missing, go to next locus
            ploid = length(locus_array[j][i])   # if not, get the # of alleles
            push!(ploidy, ploid)    # push that to the ploidy vector
            break   # break out of the loop and begin next sample
        end
    end

    @info "\n$(abspath(infile))
$(length(sample_names)) samples detected
population info must be added <---
$(length(locinames)) loci detected"

    # create loci dataframe
    loci_df = DataFrame([j => k for (j,k) in zip(Symbol.(locinames), locus_array)])

    # create samples df
    samples_df = DataFrame(
        name = sample_names,
        population = population,
        ploidy = ploidy,
        latitude = loc_xy,
        longitude = loc_xy
    )
    PopObj(samples_df, loci_df)
end

"""
    PopGen.read(infile::String; kwargs...)
Wraps `csv()`, `genepop()`, `bcf()`, and `vcf()` to read a file in as a `PopObj`. File type is
inferred from the file extension:
- delimited: .csv | .tsv | .txt
- genepop: .gen | .genepop
- variant call format: .vcf
This function uses the same keyword arguments (and defaults) as the file importing
functions it wraps; please see their respective docstrings in the Julia help console.
(e.g. `?genepop`) for specific usage details.

## Examples
`read("cavernous_assfish.gen", markers = "msat", digits = 3)`

`read("bos_tauros.csv")`

`read("juglans_nigra.vcf")`
"""
function Base.read(infile::String; kwargs...)
    kwargs = Dict(kwargs)
    ext = split(infile, ".")[end]
    if ext == "gen" || ext == "genepop"
        if haskey(kwargs, :digits)
            digits = kwargs[:digits]
        else
            digits = 3
        end
        if haskey(kwargs, :popsep)
            popsep = kwargs[:popsep]
        else
            popsep = "POP"
        end
        if haskey(kwargs, :marker)
            marker = kwargs[:marker]
        else
            marker = "snp"
        end
        return genepop(infile, digits = digits, popsep = popsep, marker = marker)

    elseif ext == "csv" || ext == "txt" || ext == "tsv"
        if haskey(kwargs, :delim)
            delim = kwargs[:delim]
        else
            delim = r"\,|\s|\t"
        end
        if haskey(kwargs, :digits)
            digits = kwargs[:digits]
        else
            digits = 3
        end
        if haskey(kwargs, :marker)
            marker = kwargs[:marker]
        else
            marker = "snp"
        end
        return csv(infile, delim = delim, digits = digits, marker = marker)

    elseif ext == "vcf" || ext == "bcf"
            ext == "vcf" && return vcf(infile)
            ext == "bcf" && return bcf(infile)
    else
        @error "file type not recognized by filename extension \n delimited: .csv | .tsv | .txt \n genepop: .gen | .genepop \n variant call formant: .bcf | .vcf"
    end
end


text:
module Morph
# Partial backwards syntax comparability for morphological operations:
#
#  - strlen
#
#       SE = strel(shape, parameters)
#       SE = strel('arbitrary', NHOOD)
#       SE = strel('arbitrary', NHOOD, HEIGHT)
#       SE = strel('ball', R, H, N)
#       SE = strel('diamond', R)
#       SE = strel('disk', R, N)
#       SE = strel('line', LEN, DEG)
#       SE = strel('octagon', R)
#       SE = strel('pair', OFFSET)
#       SE = strel('periodicline', P, V)
#       SE = strel('rectangle', MN)
#       SE = strel('square', W)
#
#  - imopen
#
#       IM2 = imopen(IM,SE)
#       IM2 = imopen(IM,NHOOD)
#
#  - imdilate
#  - imerode
#
#
# Copyright © 2014-2015 Vardan Andriasyan, Yauhen Yakimovich, Artur Yakimovich.
#
#  MIT license.
#
# Permission is hereby granted, free of charge, to any person
# obtaining a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including without
# limitation the rights to use, copy, modify, merge, publish, distribute,
# sublicense, and/or sell copies of the Software, and to permit persons to whom
# the Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

export strel

function strel(shape::AbstractArray{Bool,2}, args...)
    # Construct a proper type for StructElement.
    return shape
end
function strel(shape::String, args...)
    if shape == "disk"
        assert(length(args) in (1,2))
        radius = args[1];
        nhood = (length(args) == 2) ? args[2] : 0;
        assert(isa(radius, Number) && radius > 0)
        if radius == 1
            assert(isa(nhood, Number) && nhood in (0,2,4,8))
            return bool([0 1 0;1 1 1; 0 1 0]);
        else
            error("This case is not implemented yet!")
        end
    else
        error("This case is not implemented yet!")
    end
end

end #End of Morph


text:
using JLD2
using DataFrames
using PyPlot
using IcanProj

"""
Compute monthly summary statistics.
"""
function monthly_summary(dict_wsh)

    dict_monthly = Dict()

    for (key_wsh, value_wsh) in dict_wsh

        dict_res = Dict()

        for (key_res, value_res) in value_wsh

            # Crop first time step containing weired data

            df = value_res[2:end, :]

            df[:year] = Dates.year.(df[:time])

            df[:month] = Dates.month.(df[:time])

            delete!(df, :time)

            # Variables for summation

            df_tmp = df[[:prec, :evap, :runoff, :baseflow, :year, :month]]

            df_tmp[:total_runoff] = df_tmp[:runoff] + df_tmp[:baseflow]

            df_tmp = df_tmp |> groupby([:year, :month]) |> sum

            delete!(df_tmp, :year)

            df_sum = df_tmp |> groupby(:month) |> mean

            # Variables for averaging

            df_tmp = df[[:swe, :snow_depth, :month]]

            df_mean = df_tmp |> groupby(:month) |> mean

            # Add to dictionary

            df_all = [df_mean df_sum]

            delete!(df_all, :month_1)

            dict_res[key_res] = [df_mean df_sum]

        end

        dict_monthly[key_wsh] = dict_res

    end

    return dict_monthly

end

"""
Plot monthly variabels.
"""
function plot_monthly(dict_monthly, wsh_name, var_name, ylabel_name, plot_title, file_name)

    fig = figure(figsize=(15, 10))

    subplots_adjust(hspace = 0.0)

    i = 1

    for (key_wsh, value_wsh) in dict_monthly

        for (key_res, df) in value_wsh

            ax = fig[:add_subplot]("42$(i)")

            ax[:text](0.05, 0.8, "$(wsh_name[key_wsh])", verticalalignment = "bottom", horizontalalignment = "left", transform=ax[:transAxes])

            plot(df[:month], df[var_name], label = key_res)

            xlabel("Month")

            ylabel(ylabel_name)

        end

        i += 1

        if i == 6
            legend()
        end

    end

    suptitle(plot_title)

    file_name = Pkg.dir("IcanProj", "plots", "monthly_$(file_name).png")

    rm(file_name, force=true)
    savefig(file_name, dpi = 300)

    close(fig)

end

"""
Compute yearly average evapotranspiration.
"""
function yearly_average_evapotranspiration(dict_monthly)

    df_evap = DataFrame()

    df_evap[:Resolution] = collect(keys(dict_monthly["224.1"]))

    for (key_wsh, dict_wsh) in dict_monthly

        evap_sum = []

        for (key_res, df_res) in dict_wsh

            push!(evap_sum, sum(df_res[:evap_sum_mean]))

        end

        df_evap[Symbol(wsh_name[key_wsh])] = evap_sum

    end

    return df_evap

end

# Load all necessary results

opt = get_options()

path = opt["eval_folder"]

dict_wsh = get_summary_tables(path, opt)

dict_monthly = monthly_summary(dict_wsh)

@load Pkg.dir("IcanProj", "data", "wsh_info.jld2") wsh_info

wsh_name = get_wsh_name(wsh_info)

# Plot monthly averaged states and fluxes

plot_monthly(dict_monthly, wsh_name, :swe_mean, "SWE (mm)", "Snow water equivalent", "swe")

plot_monthly(dict_monthly, wsh_name, :snow_depth_mean, "HS (cm)", "Snow depth", "hs")

plot_monthly(dict_monthly, wsh_name, :evap_sum_mean, "Evap (mm/month)", "Evapotranspiration", "evap")

plot_monthly(dict_monthly, wsh_name, :total_runoff_sum_mean, "Runoff (mm/month)", "Runoff", "runoff")

# Dataframe with yearly mean evapotranspiration

df_evap = yearly_average_evapotranspiration(dict_monthly)

# Plot time series

function plot_time_series(dict_wsh, wsh_key, var_name, ylabel_name)

    dict_res = dict_wsh[wsh_key]

    for (key_res, df_res) in dict_res

        #plot(df_res[:time], df_res[:runoff]+df_res[:baseflow], label = key_res)
        plot(df_res[:time], df_res[var_name], label = key_res)

    end

    legend()
    title(wsh_name[wsh_key])
    ylabel(ylabel_name)

end

wsh_key = "22.4"
var_name = :evap
ylabel_name = "Evapotranspiration (mm/3h)"
plot_time_series(dict_wsh, wsh_key, var_name, ylabel_name)


text:
function acqr_kalman(y,m::acqr;s::Bool=true)

	# parameters
	a = m.a
	c = m.c
	q = m.q
	r = m.r
	m1 = m.m1
	P1 = m.P1

    if s
        # stationary
        (xtt,Ptt,xtt1,Ptt1,et,St,Kt,loglik) = acqr_kfilter_s(y,a,c,q,r,m1)
		(xtN,PtN,Pt1tN) = acqr_ksmoother_s(a,xtt,Ptt,xtt1,Ptt1)
	else
		(xtt,Ptt,xtt1,Ptt1,et,St,Kt,loglik) = acqr_kfilter(y,a,c,q,r,m1,P1)
		(xtN,PtN,Pt1tN) = acqr_ksmoother(a,xtt,Ptt,xtt1,Ptt1)
	end

	salida = acqrKal(y,m,xtt1,Ptt1,xtt,Ptt,xtN,PtN,s)

	return salida
end


text:
"""
    symmetric_A(A; ignore_asymmetry=false, chol=false)

Check if `A` is approximately symmetric and then return a symmetrized version
or, optionally, the Cholesky factorization

# Examples
```jldoctest
julia> A = [1. 1; 0 1];

julia> GPLinearODEMaker.symmetric_A(A; ignore_asymmetry=true)
2×2 LinearAlgebra.Symmetric{Float64,Array{Float64,2}}:
 1.0  0.5
 0.5  1.0
```
"""
function symmetric_A(A::Union{Matrix{T},Symmetric}; ignore_asymmetry::Bool=false, chol::Bool=false) where {T<:Real}

    # an arbitrary threshold that is meant to catch numerical errors
    thres = max(1e-6 * maximum(abs.(A)), 1e-8)

    if size(A, 1) == size(A, 2)
        max_dif = maximum(abs.(A - transpose(A)))

        if max_dif == zero(max_dif)
            A = Symmetric(A)

        elseif (max_dif < thres) || ignore_asymmetry
            # return the symmetrized version of the matrix
            A = symmetrize_A(A)
        else
            @warn "Array dimensions match, but the max dif ($max_dif) is greater than the threshold ($thres)"
            chol = false
        end
    else
        if !ignore_asymmetry
            @warn "Array dimensions do not match. The matrix can't be symmetric"
        end
        chol = false
    end

    return chol ? ridge_chol(A) : A

end

"""
    symmetrize_A(A)

Symmetrize `A` (i.e. add its transpose and divide by 2)

# Examples
```jldoctest
julia> A = [1. 1; 0 1];

julia> GPLinearODEMaker.symmetrize_A(A)
2×2 LinearAlgebra.Symmetric{Float64,Array{Float64,2}}:
 1.0  0.5
 0.5  1.0
```
"""
symmetrize_A(A::Union{Matrix{T},Symmetric}) where {T<:Real} =
    Symmetric((A + transpose(A)) / 2)
symmetrize_A(A::Symmetric) where {T<:Real} = A

"""
    ridge_chol(A)

Perform a Cholesky factorization on `A`, adding a small ridge when necessary

# Examples
```jldoctest
julia> A = [4. 2;2 10];

julia> GPLinearODEMaker.ridge_chol(A)
LinearAlgebra.Cholesky{Float64,Array{Float64,2}}
U factor:
2×2 LinearAlgebra.UpperTriangular{Float64,Array{Float64,2}}:
 2.0  1.0
  ⋅   3.0
```
"""
function ridge_chol(A::Union{Matrix{T},Symmetric}) where {T<:Real}

    # only add a small ridge (based on the smallest eigenvalue) if necessary
    try
        return cholesky(A)
    catch
        smallest_eigen = IterativeSolvers.lobpcg(A, false, 1).λ[1]
        ridge = 1.10 * abs(smallest_eigen)
        @warn "added a ridge"
        println("ridge size:          10^$(log10(ridge))")
        println("max value of array:  10^$(log10(maximum(abs.(A))))")
        return cholesky(A + UniformScaling(ridge))
    end

end
ridge_chol(A::Cholesky) where {T<:Real} = A

"""
    reconstruct_array(non_zero_entries, template_A)

Fill a matrix with `non_zero_entries` at the locations of the non-zero entries
in `template_A`

# Examples
```jldoctest
julia> template_A = [4. 0;5 0];

julia> non_zero_entries = [2., 3];

julia> GPLinearODEMaker.reconstruct_array(non_zero_entries, template_A)
2×2 Array{Float64,2}:
 2.0  0.0
 3.0  0.0
```
"""
function reconstruct_array(non_zero_entries, template_A::Matrix{T}) where {T<:Real}
    @assert length(findall(!iszero, template_A))==length(non_zero_entries)
    new_A = zeros(size(template_A))
    new_A[findall(!iszero, template_A)] = non_zero_entries
    return new_A
end

"""
    assert_positive(vars...)

# Examples
```jldoctest
julia> GPLinearODEMaker.assert_positive([1,2,3,4])

julia> GPLinearODEMaker.assert_positive([-1,2,3,4])
ERROR: AssertionError: passed a negative/0 variable that needs to be positive
[...]
```
"""
function assert_positive(vars...)
    for i in vars
        @assert all(i .> 0) "passed a negative/0 variable that needs to be positive"
    end
end

"""
    ndims(A::Cholesky)

Extends ndims to return 2 if passed a Cholesky object
"""
ndims(A::Cholesky) where {T<:Real} = 2

"""
    sendto(workers; args...)

Send the `args` to the `workers` specified by thier number IDs.
Stolen shamelessly from ParallelDataTransfer.jl

# Examples
```julia-repl
julia> sendto([1, 2], x=100, y=rand(2, 3));

julia> z = randn(10, 10); sendto(workers(), z=z);
```
"""
function sendto(workers::Union{T,Vector{T}}; args...) where {T<:Integer}
    for worker in workers
        for (var_name, var_value) in args
            @spawnat(worker, Core.eval(Main, Expr(:(=), var_name, var_value)))
        end
    end
end

"""
    auto_addprocs(;add_procs=0)

Adds either as many workers as there are CPU threads minus 2 if none are active
or the amount specified by `add_procs`
"""
function auto_addprocs(;add_procs::Integer=0)
    # only add as any processors as possible if we are on a consumer chip
    if (add_procs==0) && (nworkers()==1)
        add_procs = length(Sys.cpu_info()) - 2
    end
    addprocs(add_procs)
    println("added $add_procs workers")
end

"""
    powers_of_negative_one(power)

Finds -1 ^ power without calling ^

# Examples
```jldoctest
julia> GPLinearODEMaker.powers_of_negative_one(0) == GPLinearODEMaker.powers_of_negative_one(2) == 1
true
julia> GPLinearODEMaker.powers_of_negative_one(1) == GPLinearODEMaker.powers_of_negative_one(3) == -1
true
```
"""
powers_of_negative_one(power::Integer) = iseven(power) ? 1 : -1

"""
    remove_zeros(V)

Return the a version of the passed vector after removing all zero entries.
Could possibly be replaced with a view-based version.

# Examples
```jldoctest
julia> V = [1, 2, 0, 3, 0, 4, 5];

julia> GPLinearODEMaker.remove_zeros(V)
5-element Array{Int64,1}:
 1
 2
 3
 4
 5
```
"""
remove_zeros(V::Vector{T} where T<:Real) = V[findall(!iszero, V)]

"""
    log_laplace_approximation(H, g, logh; λ=1)

Compute the logarithm of the Laplace approxmation for the integral of a function
of the following form
∫ exp(-λ g(y)) h(y) dy ≈ exp(-λ g(y*)) h(y*) (2π/λ)^(d/2) |H(y*)|^(-1/2)
where y* is the value of y at the global mode and H is the (Hessian) matrix of
second order partial derivatives of g(y) (see slide 10 of
http://www.stats.ox.ac.uk/~steffen/teaching/bs2HT9/laplace.pdf). When used to
calculate evidences, one can set λ = 1, g(y) = -log-likelihood,
h(y) = model prior, and H(y) = the Fisher information matrix (FIM) or the
Hessian matrix of the negative log-likelihood. Could possiblly improve with
methods from Ruli et al. 2016 (https://arxiv.org/pdf/1502.06440.pdf)?
"""
function log_laplace_approximation(
    H::Union{Symmetric,Matrix{T}},
    g::Real,
    logh::Real;
    λ::Real=1
    ) where {T<:Real}

    @assert size(H, 1) == size(H, 2)
    n = size(H, 1)

    return logh - λ * g + 0.5 * (n * log(2 * π / λ) - logdet(H))

end

max_stars_repo_name: christiangil/GPLinearODEMaker

max_issues_repo_name: christiangil/GPLinearODEMaker

max_forks_repo_name: christiangil/GPLinearODEMaker


text:
using JuMP

using LinearAlgebra

function formulate_nearest_good_box_qp(
    Q::Matrix,
    c::Vector,
    y_max::Float64,
    x_target::Vector{Float64},
)
    n = length(c)
    @assert n == size(Q, 1) == size(Q, 2)
    model = JuMP.Model()
    JuMP.@variable(model, 0 ≤ x[1:n] ≤ 1)
    JuMP.@variable(model, y <= y_max)
    JuMP.@constraint(model, y ≥ -0.5dot(x, Q * x) - dot(c, x))

    JuMP.@variable(m, t[1:n] >= 0)
    for i = 1:n
        JuMP.@constraint(m, t[i] >= x[i] - x_target[i])
        JuMP.@constraint(m, t[i] >= -x[i] + x_target[i])
    end
    JuMP.@objective(model, Min, sum(t))
    return JuMP.backend(model)
end