text:
using ExtensibleUnions
using Test

function f1 end
ExtensibleUnions._update_all_methods_for_extensiblefunction!(f1)

b = Type{Vector{T} where T}
@test_throws MethodError ExtensibleUnions._update_single_method!(f1, b, Set())
@test_throws MethodError ExtensibleUnions._update_single_method!(f1, b, Set(), nothing => nothing)
@test_throws MethodError ExtensibleUnions._replace_types(b)
@test_throws MethodError ExtensibleUnions._replace_types(b, nothing => nothing)

@test ExtensibleUnions._replace_types(Union{}) == Union{}
@test ExtensibleUnions._replace_types(Union{}, Union{} => Union{Float32, String}) == Union{Float32, String}
@test ExtensibleUnions._replace_types(Union{}, Union{Float32, String} => Union{Float32, Int32, String}) == Union{}
@test ExtensibleUnions._replace_types(Union{Int32}) == Union{Int32}
@test ExtensibleUnions._replace_types(Union{Int32}, Int32 => Float32) == Union{Float32}
@test ExtensibleUnions._replace_types(Union{Int32, Float32}) == Union{Int32, Float32}
@test ExtensibleUnions._replace_types(Union{Int32, Float32}, Union{Int32, Float32} => Union{String, Symbol}) == Union{String, Symbol}
@test ExtensibleUnions._replace_types(Union{Int32, Float32}, Float32 => String) == Union{Int32, String}
@test ExtensibleUnions._replace_types(Union{Int32, Float32, AbstractString}) == Union{Int32, Float32, AbstractString}
@test ExtensibleUnions._replace_types(Union{Int32, Float32, AbstractString}, Float32 => Symbol) == Union{Int32, Symbol, AbstractString}
@test ExtensibleUnions._replace_types(Union{Int32, Float32, AbstractString}, Union{Int32, Float32, AbstractString} => Symbol) == Symbol
@test ExtensibleUnions._replace_types(Union{Int32, Float32, AbstractString, Symbol}) == Union{Int32, Float32, AbstractString, Symbol}
@test ExtensibleUnions._replace_types(Union{Int32, Float32, AbstractString, Symbol}, Float32 => Char) == Union{Int32, Char, AbstractString, Symbol}
@test ExtensibleUnions._replace_types(Union{Int32, Float32, AbstractString, Symbol}, Union{Int32, Float32, AbstractString, Symbol} => Union{Int32, Float32, AbstractString, Symbol, Char}) == Union{Int32, Float32, AbstractString, Symbol, Char}

structlist = Any[]
@genstruct!(structlist, ExtensibleUnion)
@genstruct!(structlist, ExtensibleUnion)
@genstruct!(structlist, ExtensibleUnion)
@genstruct!(structlist)
@genstruct!(structlist)
@genstruct!(structlist)
@genstruct!(structlist)
# funclist = Any[]
# @genfunc!(funclist)
# @genfunc!(funclist)
FooTrait = structlist[1]
BarTrait = structlist[2]
BazTrait = structlist[3]
extensibleunion!(FooTrait)
extensibleunion!(BarTrait)
ExtensibleUnions._update_all_methods_for_extensibleunion!(FooTrait)
ExtensibleUnions._update_all_methods_for_extensibleunion!(BarTrait, nothing=>nothing)
StructA = structlist[4]
StructB = structlist[5]
StructC = structlist[6]
StructD = structlist[7]
# f = funclist[1]
# g = funclist[2]
function f end
function g end
f(x::FooTrait) = "foo"
f(x::BarTrait) = "bar"
extensiblefunction!(f, FooTrait, BarTrait)
extensiblefunction!(f, (FooTrait, BarTrait,))
ExtensibleUnions._update_all_methods_for_extensibleunion!(FooTrait)
ExtensibleUnions._update_all_methods_for_extensibleunion!(BarTrait, nothing=>nothing)
ExtensibleUnions._update_all_methods_for_extensiblefunction!(f)
addtounion!(FooTrait, StructA)
addtounion!(FooTrait, StructB)
addtounion!(FooTrait, StructC)
addtounion!(FooTrait, StructD)
addtounion!(BarTrait, StructA)
addtounion!(BarTrait, StructB)
addtounion!(BarTrait, StructC)
addtounion!(BarTrait, StructD)
ExtensibleUnions._update_all_methods_for_extensibleunion!(FooTrait)
ExtensibleUnions._update_all_methods_for_extensibleunion!(BarTrait, nothing=>nothing)
ExtensibleUnions._update_all_methods_for_extensiblefunction!(f)
g(x::BazTrait) = "baz"
ExtensibleUnions._registry_extensibleunion_to_members[BazTrait] = Set(Any[StructA, StructB, StructC, StructD])
ExtensibleUnions._update_single_method!(g, Tuple{typeof(g), BazTrait}, Set(Any[BazTrait]))


text:
# This file is part of Kpax3. License is MIT.

function merge!(ij::Vector{Int},
                neighbours::Vector{Int},
                S::Int,
                data::Matrix{UInt8},
                priorR::PriorRowPartition,
                priorC::PriorColPartition,
                settings::KSettings,
                support::MCMCSupport,
                state::AminoAcidState)
  # number of clusters after the merge
  k = state.k - 1

  initsupportmerge!(ij, k, data, priorC, support)

  hi = state.R[ij[1]]
  hj = state.R[ij[2]]

  # logarithm of the product of sequential probabilities
  lq = 0.0

  # temporary / support variables
  u = 0
  lcp = zeros(Float64, 2)
  lw = zeros(Float64, 2)
  z = 0.0
  p = 0.0

  vi = state.v[hi] + state.v[hj]
  ni = zeros(Float64, support.m)
  for b in 1:support.m
    ni[b] = state.n1s[hi, b] + state.n1s[hj, b]
  end

  # allocate the neighbours of i and j
  for l in 1:S
    u = neighbours[l]
    lcp[1] = lcp[2] = 0.0

    # compute p(x_{u} | x_{hi,1:(u-1)}) and p(x_{u} | x_{hj,1:(u-1)})
    for b in 1:support.m
      lcp[1] += computeclusteriseqprobs!(data[b, u], b, priorC, support)
      lcp[2] += computeclusterjseqprobs!(data[b, u], b, priorC, support)
    end

    lw[1] = clusterweight(support.vi, priorR)
    lw[2] = clusterweight(support.vj, priorR)

    # (w1 * p1) / (w1 * p1 + w2 * p2) = 1 / (1 + (w2 * p2) / (w1 * p1))
    # => e^(-log(1 + e^(log(w2) + log(p2) - log(w1) - log(p1))))
    z = -log1p(exp(lw[2] + lcp[2] - lw[1] - lcp[1]))
    p = exp(z)

    if state.R[u] == hi
      updateclusteri!(u, data, support)
      lq += z
    else
      updateclusterj!(u, data, support)
      lq += log1p(-p)
    end
  end

  updatelogmargliki!(ni, vi, priorC, support)

  logratiopriorrowmerge!(k, priorR, support)
  logmarglikmerge!(state.cl, state.k, hi, hj, priorC, support)

  ratio = exp(support.lograR + support.logmlikcandidate - support.logmlik + lq)

  if ratio >= 1 || ((ratio > 0) && (rand() <= ratio))
    performmerge!(hi, hj, ni, vi, priorC, settings, support, state)
  end

  nothing
end

function initsupportmerge!(ij::Vector{Int},
                           k::Int,
                           data::Matrix{UInt8},
                           priorC::AminoAcidPriorCol,
                           support::MCMCSupport)
  support.vi = 1
  support.vj = 1

  support.ui[1] = ij[1]
  support.uj[1] = ij[2]

  for b in 1:support.m
    initclusteriweights!(data[b, ij[1]], b, k, priorC, support)
    initclusterjweights!(data[b, ij[2]], b, k, priorC, support)
  end

  nothing
end

function performmerge!(hi::Int,
                       hj::Int,
                       ni::Vector{Float64},
                       vi::Int,
                       priorC::PriorColPartition,
                       settings::KSettings,
                       support::MCMCSupport,
                       state::AminoAcidState)
  for a in 1:state.v[hj]
    state.R[state.unit[hj][a]] = hi
  end

  state.emptycluster[hj] = true

  h = 0
  for a in 1:length(state.emptycluster)
    if !state.emptycluster[a]
      h += 1
      state.cl[h] = a
    end
  end

  state.k = h

  for b in 1:support.m
    state.n1s[hi, b] = ni[b]

    support.lp[1, hi, b] = support.lpi[1, b]
    support.lp[2, hi, b] = support.lpi[2, b]
    support.lp[3, hi, b] = support.lpi[3, b]
    support.lp[4, hi, b] = support.lpi[4, b]
  end

  resize!(state.unit[hi], vi)
  copyto!(state.unit[hi], state.v[hi] + 1, state.unit[hj], 1, state.v[hj])

  state.v[hi] = vi

  state.logpR += support.lograR

  support.logmlik = support.logmlikcandidate

  nothing
end

ext: jl


text:
"""
    mlaplacian(img; dims=coords_spatial(img), r=1)
    mlaplacian(img, se)

Calculate morphological laplacian of the image.

The morphological lapalacian operator is defined as `∇⁺A - ∇⁻A` where `∇⁺A` is the external
gradient `A - erode(A, se)` and `∇⁻A` is the internal gradient `dilate(A, se) - A`. Thus is
`dilate(A, se) + erode(A, se) - 2A`.

`se` is the structuring element that defines the neighborhood of the image. See
[`strel`](@ref) for more details. If `se` is not specified, then it will use the
[`strel_box`](@ref) with an extra keyword `dims` to control the dimensions to filter.

# Examples

```jldoctest; setup = :(using ImageMorphology)
julia> img = falses(7, 7); img[3:5, 3:5] .= true; img[4, 4] = false; img
7×7 BitMatrix:
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0
 0  0  1  1  1  0  0
 0  0  1  0  1  0  0
 0  0  1  1  1  0  0
 0  0  0  0  0  0  0
 0  0  0  0  0  0  0

julia> Int.(mlaplacian(img))
7×7 Matrix{$Int}:
 0  0   0   0   0  0  0
 0  1   1   1   1  1  0
 0  1  -1  -1  -1  1  0
 0  1  -1   1  -1  1  0
 0  1  -1  -1  -1  1  0
 0  1   1   1   1  1  0
 0  0   0   0   0  0  0

julia> Int.(mlaplacian(img, strel_diamond(img))) # use diamond shape SE
7×7 Matrix{$Int}:
 0  0   0   0   0  0  0
 0  0   1   1   1  0  0
 0  1  -1  -1  -1  1  0
 0  1  -1   1  -1  1  0
 0  1  -1  -1  -1  1  0
 0  0   1   1   1  0  0
 0  0   0   0   0  0  0
```

## See also

- [`mlaplacian!`](@ref) is the in-place version of this function.
- [`mgradient`](@ref) for the gradient operator.
- `ImageBase.FiniteDiff` also provides a few finite difference operators, including `fdiff`,
  `fgradient`, etc.
"""
function mlaplacian(img; dims=coords_spatial(img), r=nothing)
    return mlaplacian(img, strel_box(img, dims; r))
end
function mlaplacian(img::AbstractArray{T}, se) where {T}
    out = similar(img, maybe_floattype(T))
    buffer = similar(img, maybe_floattype(T))
    return mlaplacian!(out, img, se, buffer)
end

"""
    mlaplacian!(out, img, buffer; [dims], [r])
    mlaplacian!(out, img, se, buffer)

The in-place version of [`mlaplacian`](@ref) with input image `img` and output image `out`.
The intermediate erosion result is stored in `buffer`.
"""
function mlaplacian!(out, img, buffer; dims=coords_spatial(img), r=nothing)
    return mlaplacian!(out, img, strel_box(img, dims; r), buffer)
end
function mlaplacian!(out, img, se, buffer)
    require_symmetric_strel(se)
    dilate!(out, img, se)
    erode!(buffer, img, se)
    @. out = out + buffer - 2img
    return out
end
mlaplacian!(out::AbstractArray{<:Color3}, img, se, buffer=nothing) = throw(ArgumentError("color image is not supported"))


text:
#################################################################################
#
#   FUNCTIONS FOR ADDING / REMOVING SITES IN LATTICES AND UNITCELLS
#
#   STRUCTURE OF THE FILE:
#   - removing single sites
#       - from unitcells
#       - from lattices
#   - removing multiple sites
#       - from unitcells
#       - from lattices
#   - removing disconnected sites
#       - from unitcells
#       - from lattices
#
#   - adding single sites
#       - to unitcells
#       - to lattices
#   - adding multiple sites
#       - to unitcells
#       - to lattices
#
################################################################################

################################################################################
#
#   REMOVING SINGLE SITES
#       - from unitcells
#       - from lattices
#
################################################################################

# Remove a single site from a unitcell (and return that)
function removeSite!(
        unitcell :: U,
        index    :: Integer
    ) :: S where {D,N,LS,LB, S<:AbstractSite{LS,D}, B<:AbstractBond{LB,N}, U<:AbstractUnitcell{S,B}}

    # store the site to be removed
    s = sites(unitcell)[index]

    # remove the site from the list and pass the new list again to the unitcell
    sites!(unitcell, deleteat!(sites(unitcell), index))

    # get all bonds that are not connected to the removed site
    bond_list = filter(b -> to(b) != index && from(b) != index, bonds(unitcell))
    # relabeling of bond indices
    for b in bond_list
        # set the new indices
        from!(b, from(b)>index ? from(b)-1 : from(b))
        to!(b,   to(b)>index   ? to(b)-1   : to(b))
    end
    # set the bond list in the unitcell object
    bonds!(unitcell, bond_list)
    # return the removed site
    return s
end

# Remove a single site from a lattice (and return that)
function removeSite!(
        lattice :: L,
        index   :: Integer
    ) :: SL where {D,N,LS,LB, S<:AbstractSite{LS,D}, B<:AbstractBond{LB,N}, U<:AbstractUnitcell{S,B},
        DL,NL,LSL,LBL,SL<:AbstractSite{LSL,DL}, BL<:AbstractBond{LBL,NL}, L<:AbstractLattice{SL,BL,U}}

    # store the site to be removed
    s = sites(lattice)[index]

    # remove the site from the list and pass the new list again to the lattice
    sites!(lattice, deleteat!(sites(lattice), index))

    # get all bonds that are not connected to the removed site
    bond_list = filter(b -> to(b) != index && from(b) != index, bonds(lattice))
    # relabeling of bond indices
    for b in bond_list
        # set the new indices
        from!(b, from(b)>index ? from(b)-1 : from(b))
        to!(b,   to(b)>index   ? to(b)-1   : to(b))
    end
    # set the bond list in the lattice object
    bonds!(lattice, bond_list)
    # return the removed site
    return s
end

# export the function
export removeSite!

################################################################################
#
#   REMOVING MULTIPLE SITES
#       - from unitcells
#       - from lattices
#
################################################################################

# Remove multiple sites from a unitcell (and return them)
function removeSite!(
        unitcell :: U,
        index    :: Integer...
    ) :: Vector{S} where {D,N,LS,LB, S<:AbstractSite{LS,D}, B<:AbstractBond{LB,N}, U<:AbstractUnitcell{S,B}}

    # return the respective function
    return removeSite!(unitcell, unique(index))
end

# Remove multiple sites from a unitcell (and return them)
function removeSite!(
        unitcell :: U,
        index    :: Vector{<:Integer}
    ) :: Vector{S} where {D,N,LS,LB, S<:AbstractSite{LS,D}, B<:AbstractBond{LB,N}, U<:AbstractUnitcell{S,B}}

    # get the list of sites
    sites_to_return = S[site(unitcell,i) for i in index]

    # get all sites of the unitcell
    site_list = sites(unitcell)
    # build a index list
    indices_new = collect(1:length(site_list))
    index_gets_removed = Bool[false for i in 1:numSites(unitcell)]
    # fill the new index list at positions of removed sites
    for i in index
        indices_new[i] = -1
        index_gets_removed[i] = true
    end
    # fill the new list successively from smallest to largest
    if indices_new[1] == -1
        indices_new[1] = 0
    end
    for i in 2:length(indices_new)
        if indices_new[i] == -1
            indices_new[i] = indices_new[i-1]
        else
            indices_new[i] = indices_new[i-1] + 1
        end
    end
    # get a new site list
    site_list = S[site_list[i] for i in 1:length(site_list) if !index_gets_removed[i]]
    # register new list of sites within the unitcell object
    sites!(unitcell, site_list)

    # get all bonds that are not connected to any removed site
    bond_list = filter(b -> !(index_gets_removed[to(b)] || index_gets_removed[from(b)]), bonds(unitcell))
    # relabeling of bond indices
    for b in bond_list
        # set the new indices
        from!(b, indices_new[from(b)])
        to!(b,   indices_new[to(b)])
    end
    # set the bond list in the unitcell object
    bonds!(unitcell, bond_list)

    # return the removed sites
    return sites_to_return
end

# Remove multiple sites from a lattice (and return them)
function removeSite!(
        lattice :: L,
        index   :: Integer ...
    ) :: Vector{SL} where {D,N,LS,LB, S<:AbstractSite{LS,D}, B<:AbstractBond{LB,N}, U<:AbstractUnitcell{S,B},
        DL,NL,LSL,LBL,SL<:AbstractSite{LSL,DL}, BL<:AbstractBond{LBL,NL}, L<:AbstractLattice{SL,BL,U}}

    # call the respective funtion
    return removeSite!(lattice, unique(index))
end

# Remove multiple sites from a lattice (and return them)
function removeSite!(
        lattice :: L,
        index   :: Vector{<:Integer}
    ) :: Vector{SL} where {D,N,LS,LB, S<:AbstractSite{LS,D}, B<:AbstractBond{LB,N}, U<:AbstractUnitcell{S,B},
        DL,NL,LSL,LBL,SL<:AbstractSite{LSL,DL}, BL<:AbstractBond{LBL,NL}, L<:AbstractLattice{SL,BL,U}}

    # get the list of sites
    sites_to_return = SL[site(lattice,i) for i in index]

    # get all sites of the lattice
    site_list = sites(lattice)
    # build a index list
    indices_new = collect(1:length(site_list))
    index_gets_removed = Bool[false for i in 1:numSites(lattice)]
    # fill the new index list at positions of removed sites
    for i in index
        indices_new[i] = -1
        index_gets_removed[i] = true
    end
    # fill the new list successively from smallest to largest
    if indices_new[1] == -1
        indices_new[1] = 0
    end
    for i in 2:length(indices_new)
        if indices_new[i] == -1
            indices_new[i] = indices_new[i-1]
        else
            indices_new[i] = indices_new[i-1] + 1
        end
    end
    # get a new site list
    site_list = SL[site_list[i] for i in 1:length(site_list) if !index_gets_removed[i]]
    # register new list of sites within the lattice object
    sites!(lattice, site_list)

    # get all bonds that are not connected to any removed site
    bond_list = filter(b -> !(index_gets_removed[to(b)] || index_gets_removed[from(b)]), bonds(lattice))
    # relabeling of bond indices
    for b in bond_list
        # set the new indices
        from!(b, indices_new[from(b)])
        to!(b,   indices_new[to(b)])
    end
    # set the bond list in the lattice object
    bonds!(lattice, bond_list)

    # return the removed sites
    return sites_to_return
end

################################################################################
#
#   REMOVING DISCONNECTED SITES
#       - from unitcells
#       - from lattices
#
################################################################################

# Remove multiple disconnected sites (from site with index) from a unitcell (and return them)
function removeDisconnectedSites!(
        unitcell :: U,
        index    :: Integer = 1
    ) :: Vector{S} where {D,N,LS,LB, S<:AbstractSite{LS,D}, B<:AbstractBond{LB,N}, U<:AbstractUnitcell{S,B}}

    #######################################
    # label all sites by Hoshen Kopelman
    #######################################

    # label variables
    current_label = 2
    site_labels = zeros(Int64, numSites(unitcell)) .- 1
    site_labels[index] = 1

    # the list of real labels
    real_labels = zeros(Int64, numSites(unitcell)) .- 1
    function reallabel(l :: Int64) :: Int64
        if real_labels[l] == -1
            return l
        else
            return reallabel(real_labels[l])
        end
    end

    # the list of bonds (organized)
    bond_list = organizedBondsFrom(unitcell)

    # iterate over all sites
    for s in 1:numSites(unitcell)
        # skip the preset one
        if s == index
            continue
        end
        # check neighbors
        for nb in bond_list[s]
            # check if a label is already set
            if site_labels[to(nb)] > 0
                # set the current label
                if site_labels[s] > 0 && reallabel(site_labels[s])==reallabel(site_labels[to(nb)])
                    # do nothing
                elseif site_labels[s] > 0
                    # make new connection
                    real_labels[max(site_labels[s], site_labels[to(nb)])] = min(site_labels[s], site_labels[to(nb)])
                    # label accordingly
                    site_labels[s] = min(site_labels[s], site_labels[to(nb)])
                else
                    # just use the label
                    site_labels[s] = reallabel(site_labels[to(nb)])
                end
            end
        end
        # if still no label, create new one
        if site_labels[s] < 0
            site_labels[s] = current_label
            current_label = current_label + 1
        end
    end

    # replace all labels by their real counterpart
    site_labels = map(s -> reallabel(s), site_labels)

    # remove all labels that are not 1
    sites_to_remove = [i for i in 1:numSites(unitcell) if site_labels[i]!=1]
    if length(sites_to_remove) == 0
        return S[]
    else
        return removeSite!(unitcell, sites_to_remove)
    end
end

# Remove multiple disconnected sites (from site with index) from a lattice (and return them)
function removeDisconnectedSites!(
        lattice :: L,
        index   :: Integer = 1
    ) :: Vector{SL} where {D,N,LS,LB, S<:AbstractSite{LS,D}, B<:AbstractBond{LB,N}, U<:AbstractUnitcell{S,B},
        DL,NL,LSL,LBL,SL<:AbstractSite{LSL,DL}, BL<:AbstractBond{LBL,NL}, L<:AbstractLattice{SL,BL,U}}

    #######################################
    # label all sites by Hoshen Kopelman
    #######################################

    # label variables
    current_label = 2
    site_labels = zeros(Int64, numSites(lattice)) .- 1
    site_labels[index] = 1

    # the list of real labels
    real_labels = zeros(Int64, numSites(lattice)) .- 1
    function reallabel(l :: Int64) :: Int64
        if real_labels[l] == -1
            return l
        else
            return reallabel(real_labels[l])
        end
    end

    # the list of bonds (organized)
    bond_list = organizedBondsFrom(lattice)

    # iterate over all sites
    for s in 1:numSites(lattice)
        # skip the preset one
        if s == index
            continue
        end
        # check neighbors
        for nb in bond_list[s]
            # check if a label is already set
            if site_labels[to(nb)] > 0
                # set the current label
                if site_labels[s] > 0 && reallabel(site_labels[s])==reallabel(site_labels[to(nb)])
                    # do nothing
                elseif site_labels[s] > 0
                    # make new connection
                    real_labels[max(site_labels[s], site_labels[to(nb)])] = min(site_labels[s], site_labels[to(nb)])
                    # label accordingly
                    site_labels[s] = min(site_labels[s], site_labels[to(nb)])
                else
                    # just use the label
                    site_labels[s] = reallabel(site_labels[to(nb)])
                end
            end
        end
        # if still no label, create new one
        if site_labels[s] < 0
            site_labels[s] = current_label
            current_label = current_label + 1
        end
    end

    # replace all labels by their real counterpart
    site_labels = map(s -> reallabel(s), site_labels)

    # remove all labels that are not 1
    sites_to_remove = [i for i in 1:numSites(lattice) if site_labels[i]!=1]
    if length(sites_to_remove) == 0
        return S[]
    else
        return removeSite!(lattice, sites_to_remove)
    end
end

# export functions
export removeDisconnectedSites!

################################################################################
#
#   ADDING SINGLE SITES
#       - to unitcells
#       - to lattices
#
################################################################################

# Add a single site to a unitcell and return its index in the list
function addSite!(
            unitcell :: U,
            position :: Vector{<:Real},
            label    :: L = getDefaultLabel(L)
        ) :: Int64 where {L,D,B,S<:AbstractSite{L,D},U<:AbstractUnitcell{S,B}}

    # create a new site
    site = newSite(S, position, label)

    # obtain the list of sites in the unitcell
    site_list = sites(unitcell)
    # push to the list of sites
    push!(site_list, site)
    # set the list of sites in the unitcell
    sites!(unitcell, site_list)

    # return the new site
    return length(site_list)
end

# Add a single site to a lattice and return its index in the list
function addSite!(
            lattice  :: LA,
            position :: Vector{<:Real},
            label    :: L = getDefaultLabel(L)
        ) :: Int64 where {L,D,B,S<:AbstractSite{L,D},U,LA<:AbstractLattice{S,B,U}}

    # create a new site
    site = newSite(S, position, label)

    # obtain the list of sites in the lattice
    site_list = sites(lattice)
    # push to the list of sites
    push!(site_list, site)
    # set the list of sites in the lattice
    sites!(lattice, site_list)

    # return the new site
    return length(site_list)
end

# export the function
export addSite!


text:
using ITensors
using ITensorUnicodePlots
using Graphs
using LayeredLayouts
using PastaQ: qft

N = 4
gates = qft(N)

s = siteinds("Qubit", N)

U, s̃ = circuit_network(gates, s)
ψ = MPS(s)
ψ̃ = MPS(s̃)
tn = [ψ..., U..., ψ̃...]

edge_labels = (; tags=true, plevs=true)
layout(g) = layered_layout(solve_positions(Zarate(), g))
@visualize fig tn arrow_show = true edge_labels = edge_labels edge_textsize = 20 layout =
  layout width = 100 height = 50

fig


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

using Markdown
using InteractiveUtils

# ╔═╡ 028bce7d-e731-47e9-9468-84185721f0ab
using Pkg; Pkg.add(["Plots","LinearAlgebra"])

# ╔═╡ c912668e-a13f-11eb-0230-8be9b4bd6d15
using Plots,LinearAlgebra

# ╔═╡ f04b3d61-d68c-4349-ab98-7136acc75447
using CSV,DataFrames

# ╔═╡ aa7fc637-2e9d-4c14-a284-50294ba7cbe5
md"""
# Define all the functions
"""

# ╔═╡ 89e3a6ae-4f5c-4882-81fd-215379d609dd
function ansatz(x,α)
	x = norm(x)
	return exp(-α*x)
end;

# ╔═╡ 1197b58f-7614-4e7c-adeb-96af732b7916
function probability_density(x,α)
	# note that it is not normalised becuase it does not matter!
	return ansatz(x,α)^2
end;

# ╔═╡ db71e022-2374-45bc-98af-3fd5e7256784
function E_local(x,α)
	x = norm(x)
	return -1/x - (α/2)*(α - 2/x)
end;

# ╔═╡ c59f8de5-73db-495f-9a7a-a4b9bf0b58bc
md" # Function for running the metropolis algorithm."

# ╔═╡ 7e5edaff-8597-4f8e-9104-047fff790e99
function metropolis(N,α)
	L = 3/sqrt(2*α)
	x = rand(3)*2*L .- L
	E = 0
	E2 = 0
	Eln_avg = 0
	ln_avg = 0
	rejection_ratio = 0
	#start the loop
	for i in 1:N
		x_t = x + 0.2*(rand(3)*2*L .- L)
		if probability_density(x_t,α) >= probability_density(x,α)
			x = x_t
		else
			rand_number = rand()
			if rand_number < probability_density(x_t,α)/probability_density(x,α)
				x = x_t
			else
				rejection_ratio += 1/N
			end
		end

		E += E_local(x,α)/N
		E2 += E_local(x,α)^2/N
		Eln_avg += -1 * E_local(x,α) * (norm(x)/N)
        ln_avg += -1 * norm(x)/N

	end
	return E,E2,Eln_avg,ln_avg,rejection_ratio
end;

# ╔═╡ 3d1f4ba0-6eec-4d3b-941f-9ad981b1773e
md"## Iterations for finding the α which minimizes the energy expectation value."

# ╔═╡ bf659597-2c99-4af6-8bfd-ba9adbba3a3c
function iterate_alpha()
	α = 0.2
	alpha_iter = 30
	N = 500
	random_walkers = 200
	γ = 0.5
	energy = Vector{Float64}(undef,alpha_iter)
	alpha = Vector{Float64}(undef,alpha_iter)
	variance = Vector{Float64}(undef,alpha_iter)
	for i in 1:alpha_iter
		E = 0
		E2 = 0
		dE_dalpha = 0
		Eln = 0
		ln = 0
		rejection_ratio = 0
		for j in 1:random_walkers
			E, E2, Eln, ln, rejections = metropolis(N,α)
			E += E/random_walkers
			E2 += E2/random_walkers
			Eln += Eln/random_walkers
			ln += ln/random_walkers
			# rejection_ratio += rejections_met/random_walkers
		end

		# Define the next α
		dE_dalpha = 2*(Eln - E*ln)
		α = α - γ*dE_dalpha
		# α = α + 0.05

		# Update the arrays for plotting
		energy[i] = E
		alpha[i] = α
		variance[i] = E2 - E^2

	end
	return energy,alpha,variance
end;

# ╔═╡ d5f7f53c-294a-4c46-a3ed-12a8f3a26724
md"# Running the simulation."

# ╔═╡ cee4be07-8a7c-4ebb-a3a3-16ba865bce64
energy,alpha,variance = iterate_alpha();

# ╔═╡ 575bb7ee-49bc-4d38-b473-086e32684261
md" # Plots"

# ╔═╡ 1411e1d3-86e9-4d7b-af64-133aee70636c
plot(energy,
	label = false,
	title = "Energy expenctation Plot",
	linecolor = :blue,
	ylabel = "<E>",
	xlabel = "timestep"
)

# ╔═╡ 1c094b3f-b76f-4b8c-8aa9-efc8967352fc
plot(variance,
	title = "Variance in energy",
	label = false,
	linecolor = :blue,
	ylabel = "Var(E)",
	xlabel = "timestep"
)

# ╔═╡ f0ad854a-5d32-42f5-81e7-9e850081a64e
plot(alpha,
	title = "Variation Parameter",
	label = false,
	linecolor = :blue,
	ylabel = "α",
	xlabel = "timestep"
)

# ╔═╡ a1e1fc43-87d3-4f1e-9d3a-1549aa115aef
scatter(alpha,
	energy,
	title = "Energy vs α",
	xlabel = "α",
	ylabel = "<E>",
	label = false,
	markercolor = :blue
)

# ╔═╡ Cell order:
# ╠═028bce7d-e731-47e9-9468-84185721f0ab
# ╠═c912668e-a13f-11eb-0230-8be9b4bd6d15
# ╠═f04b3d61-d68c-4349-ab98-7136acc75447
# ╟─aa7fc637-2e9d-4c14-a284-50294ba7cbe5
# ╠═89e3a6ae-4f5c-4882-81fd-215379d609dd
# ╠═1197b58f-7614-4e7c-adeb-96af732b7916
# ╠═db71e022-2374-45bc-98af-3fd5e7256784
# ╟─c59f8de5-73db-495f-9a7a-a4b9bf0b58bc
# ╠═7e5edaff-8597-4f8e-9104-047fff790e99
# ╟─3d1f4ba0-6eec-4d3b-941f-9ad981b1773e
# ╠═bf659597-2c99-4af6-8bfd-ba9adbba3a3c
# ╟─d5f7f53c-294a-4c46-a3ed-12a8f3a26724
# ╠═cee4be07-8a7c-4ebb-a3a3-16ba865bce64
# ╟─575bb7ee-49bc-4d38-b473-086e32684261
# ╠═1411e1d3-86e9-4d7b-af64-133aee70636c
# ╠═1c094b3f-b76f-4b8c-8aa9-efc8967352fc
# ╠═f0ad854a-5d32-42f5-81e7-9e850081a64e
# ╠═a1e1fc43-87d3-4f1e-9d3a-1549aa115aef

max_stars_repo_name: Ved-Mahajan/Variational-Monto-Carlo

max_issues_repo_name: Ved-Mahajan/Variational-Monto-Carlo

max_forks_repo_name: Ved-Mahajan/Variational-Monto-Carlo


text:
using LinearAlgebra
"""
    Compute Lyapunov Exponents and Covariant
    Lyapunov vectors.
    Reference: Ginelli 2013
	Inputs:
		DTu: dxdxm m-length timeseries of jacobian matrices
		du: tangent subspace dimension that the computed CLVs must span.
	Outputs:
		lyap_exps: the first du Lyapunov exponents
		Q: dxduxm m-length timeseries of the first du CLVs

"""
function clvs(DTu::Array{Float64,3},du::Int64)
    d = size(DTu)[1]
    m = size(DTu)[3]
    lyap_exps = zeros(du)
    R = zeros(du,du,m)
    Q = zeros(d,du,m)
    A = qr!(rand(d,du))
	Q[:,:,1] = Array(A.Q)
	R[:,:,1] = A.R
    for i=2:m
        Q[:,:,i] = DTu[:,:,i-1]*Q[:,:,i-1]
        A = qr!(Q[:,:,i])
		Q[:,:,i] = Array(A.Q)
		R[:,:,i] = A.R
        lyap_exps .+= log.(abs.(diag(R[:,:,i])))./m
    end
    C = zeros(du,du,m)
    C[:,:,end] = diagm(ones(du))
    for i=reverse(1:m-1)
        C[:,:,i] = R[:,:,i+1]\C[:,:,i+1]
        [normalize!(view(C,:,j,i)) for j = 1:du]
        Q[:,:,i] = Q[:,:,i]*C[:,:,i]
    end
    return lyap_exps, Q
end

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    DensityRatioValidation(k; [parameters])

Desntity ratio validation where weights are first obtained with density
ratio estimation, and then used in `k`-fold weighted cross-validation.

## Parameters

* `estimator` - Density ratio estimator (default to `LSIF()`)
* `optlib`    - Optimization library (default to `default_optlib(estimator)`)
* `lambda`    - Power of density ratios (default to `1.0`)

Please see [DensityRatioEstimation.jl]
(https://github.com/JuliaEarth/DensityRatioEstimation.jl)
for a list of supported estimators.

## References

* TODO
"""
struct DensityRatioValidation{T,E,O} <: AbstractErrorEstimator
  k::Int
  lambda::T
  dre::E
  optlib::O
  loss::Dict{Symbol,SupervisedLoss}
end

function DensityRatioValidation(k::Int; lambda=1.0, loss=Dict(),
                                estimator=LSIF(), optlib=default_optlib(estimator))
  @assert k > 0 "number of folds must be positive"
  @assert 0 ≤ lambda ≤ 1 "lambda must lie in [0,1]"
  T = typeof(lambda)
  E = typeof(estimator)
  O = typeof(optlib)
  DensityRatioValidation{T,E,O}(k, lambda, estimator, optlib, loss)
end

function error(solver::AbstractLearningSolver,
               problem::LearningProblem,
               eestimator::DensityRatioValidation)
  # weight samples based on the features of target data
  tdata = targetdata(problem)
  vars = collect(features(task(problem)))
  weighter = DensityRatioWeighter(tdata, vars,
                                  estimator=eestimator.dre,
                                  optlib=eestimator.optlib)
  wcv = WeightedCrossValidation(weighter,
                                eestimator.k, shuffle=true,
                                lambda=eestimator.lambda,
                                loss=eestimator.loss)
  error(solver, problem, wcv)
end


text:
using Plots
import Distributions
using Random
using Turing
using Statistics
using StatsPlots
using FillArrays

dataerr = 1.0
d = Cauchy(0,dataerr)

x = 0:10
y = 10.0 .* exp.(-x*0.1) .+ rand(d,length(x))

scatter(x,y,label=false)
#plot!(x_fit,y_fit,label=false)

@model function linear_fit(x,y)
    σ ~ Uniform(0.0,20.0)
    b ~ Uniform(-100.0,100.0)
    m ~ Uniform(-100.0,100.0)
    for i in 1:length(x)
        y[i] ~ Cauchy(b+m*x[i],σ)
    end
end

model_fit = linear_fit(x,y)
chain = sample(model_fit,NUTS(0.65),1000)
plot(chain)

mean(vec(chain[:b]))
std(vec(chain[:b]))

mean(vec(chain[:m]))
std(vec(chain[:m]))

@model function nonlinear_fit(x,y,f)
    σ ~ Uniform(0.0,20.0)
    b ~ Uniform(0.0,100.0)
    m ~ Uniform(0.0,1.0)
    for i in 1:length(x)
        y[i] ~ Cauchy(f(x[i],b,m),σ)
    end
end

function f(x,b,m)
    b*exp(-x*m)
end

model_nonlinear = nonlinear_fit(x,y,f)
chain2 = sample(model_nonlinear,NUTS(0.65),1000)

plot(chain2)


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    DomainView(domain, inds)

Return a view of `domain` at `inds`.
"""
struct DomainView{𝒟,I} <: AbstractDomain
  domain::𝒟
  inds::I
end

nelms(dv::DomainView) = length(dv.inds)
ncoords(dv::DomainView) = ncoords(dv.domain)
coordtype(dv::DomainView) = coordtype(dv.domain)
coordinates!(buff::AbstractVector, dv::DomainView, ind::Int) =
  coordinates!(buff, dv.domain, dv.inds[ind])

Base.collect(dv::DomainView) = PointSet(coordinates(dv))

# ------------
# IO methods
# ------------
function Base.show(io::IO, dv::DomainView)
  N = ncoords(dv)
  T = coordtype(dv)
  npts = nelms(dv)
  print(io, "$npts DomainView{$T,$N}")
end

function Base.show(io::IO, ::MIME"text/plain", dv::DomainView)
  println(io, dv)
  Base.print_array(io, coordinates(dv))
end


text:
function check_component( cmpnnt_root, is_subsys = false )  #=

checks if a component, described by the Industrial System Description Format (iSDF), is grammatcally correct or not.

A component is a functioning unit defined in the iSDF. A component usually consists of a physical part embedded with its control part but can consist only of a physical part or a control part. One or more components can be used to define another component.

The first argument "cmpnnt_root" gives the root element of the component. If the second argument "is_subsys" is unity, the checking is applied to a subsystem instead of a component.

This function requires that the Julia package EzXML has been loaded prior to execution.

Author: Taku Noda
Started on Oct. 31, 2019
Update on June 29, 2021

Copyright 2021 Central Research Institute of Electric Power Industry

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.

=#

path = nodepath( cmpnnt_root )
tag = cmpnnt_root.name
if !is_subsys
  if tag != "component"
    isdf_error( "In a component definition, the tag of the root XML element must be <component>. The tag <$tag>, which is not valid, has been found.", path )
  end
  if !haskey( cmpnnt_root, "name" )
    isdf_error( "The <component> ... </component> section does not have a \"name\" attribute. A <component> ... </component> section must have a \"name\" attribute so as to identify the component by its name.", path )
  end
  name = cmpnnt_root[ "name" ]
  if isempty( name )
    isdf_error( "The \"name\" attribute of the <component> ... </component> section is empty. It cannot be empty.", path )
  end
  if !haskey( cmpnnt_root, "type" )
    isdf_error( "The component \"$name\" does not have a \"type\" attribute. A <component> ... </component> section must have a \"type\" attribute so as to specify the type of the component.", path )
  end
  type = cmpnnt_root[ "type" ]
  if isempty( type )
    isdf_error( "The \"type\" attribute of the component \"$name\" is empty. It cannot be empty.", path )
  end
else  # if is_subsys,
  if tag != "subsystem"
    isdf_error( "In a subsystem definition, the tag of the root XML element must be <subsystem>. The tag <$tag>, which is not valid, has been found.", path )
  end
  if !haskey( cmpnnt_root, "name" )
    isdf_error( "The <subsystem> ... </subsystem> section does not have a \"name\" attribute. A <subsystem> ... </subsystem> section must have a \"name\" attribute so as to identify the subsystem by its name.", path )
  end
  name = cmpnnt_root[ "name" ]
  if isempty( name )
    isdf_error( "The \"name\" attribute of the <subsystem> ... </subsystem> section is empty. It cannot be empty.", path )
  end
  if !haskey( cmpnnt_root, "type" )
    isdf_error( "The subsystem \"$name\" does not have a \"type\" attribute. A <subsystem> ... </subsystem> section must have a \"type\" attribute so as to specify the type of the subsystem.", path )
  end
  type = cmpnnt_root[ "type" ]
  if isempty( type )
    isdf_error( "The \"type\" attribute of the subsystem \"$name\" is empty. It cannot be empty.", path )
  end
end

num_calcs = 0
num_branches = 0
num_ctrls = 0
num_uses = 0
for elm in eachelement( cmpnnt_root )
  path = nodepath( elm )
  tag = elm.name
  if tag == "calc"  # checking validity of calc sections:
    num_calcs = num_calcs + 1
    if num_calcs > 1
      if !is_susys
        isdf_error( "More than one <calc> ... </calc> section has been found. A component definition can have one or no <calc> ... </calc> section.", path )
      else
        isdf_error( "More than one <calc> ... </calc> section has been found. A subsystem definition can have one or no <calc> ... </calc> section.", path )
      end
    end
  elseif tag == "branch"  # checking validity of branches:
    num_branches = num_branches + 1
    ord_branch = ordinal_string( num_branches )
    if !haskey( elm, "name" )
      isdf_error( "The $ord_branch <branch> ... </branch> section does not have a \"name\" attribute. A <branch> ... </branch> section must have a \"name\" attribute so as to identify the branch by its name.", path )
    end
    name = elm[ "name" ]
    if isempty( name )
      isdf_error( "The \"name\" attribute of the $ord_branch <branch> ... </branch> section is empty. It cannot be empty.", path )
    end
    if !haskey( elm, "type" )
      isdf_error( "The branch \"$name\" does not have a \"type\" attribute. A <branch> ... </branch> section must have a \"type\" attribute so as to specify the type of the branch.", path )
    end
    type = elm[ "type" ]
    if isempty( type )
      isdf_error( "The \"type\" attribute of the branch \"$name\" is empty. It cannot be empty.", path )
    end
    # checking if the given type is an existing type or not:
    exists = false
    for type1 in BranchTypes
      if type == type1
        exists = true
        break
      end
    end
    if exists == false
      isdf_error( "Undefined branch type \"$type\" has been found for the branch \"$name\"." )
    end
    if ( type == "R" )||( type == "L" )||( type == "G" )||( type == "C" )
      if !haskey( elm, "nph" )  # if single-phase,
        if !haskey( elm, "val" )
          isdf_error( "The branch \"$name\" does not have a \"val\" attribute. A single-phase R, L, G or C branch must have a \"val\" attribute so as to specify the value of the branch. Since no \"nph\" attribute is given, it is assumed that this branch is single-phase.", path )
        end
        val = elm[ "val" ]
        if isempty( val )
          isdf_error( "The \"val\" attribute of the branch \"$name\" is empty. It cannot be empty.", path )
        end
        if !haskey( elm, "lhs" )
          isdf_error( "The branch \"$name\" does not have a \"lhs\" attribute. A single-phase R, L, G or C branch must have a \"lhs\" attribute so as to specify the node to which the left-hand-side pin is connected. Since no \"nph\" attribute is given, it is assumed that this branch is single-phase.", path )
        end
        lhs = elm[ "lhs" ]
        if isempty( lhs )
          isdf_error( "The \"lhs\" attribute of the branch \"$name\" is empty. It cannot be empty.", path )
        end
        # eliminating unnecessary white spaces:
        elm[ "lhs" ] = replace( lhs, " " => "" )
        if !haskey( elm, "rhs" )
          isdf_error( "The branch \"$name\" does not have a \"rhs\" attribute. A single-phase R, L, G or C branch must have a \"rhs\" attribute so as to specify the node to which the right-hand-side pin is connected. Since no \"nph\" attribute is given, it is assumed that this branch is single-phase.", path )
        end
        rhs = elm[ "rhs" ]
        if isempty( rhs )
          isdf_error( "The \"rhs\" attribute of the branch \"$name\" is empty. It cannot be empty.", path )
        end
        # eliminating unnecessary white spaces:
        elm[ "rhs" ] = replace( rhs, " " => "" )
      else  # if mutiphase,
        nph = elm[ "nph" ]
        if isempty( nph )
          isdf_error( "The \"nph\" attribute of the branch \"$name\" is empty. It cannot be empty.", path )
        end
        if parse( Int, nph ) < 2
          omega_error( "For the branch \"$name\", the number of phases is given as the \"nph\" attribute, and this means that the branch is multiphase. However, it is set to a value less than two. The \"nph\" attribute of a multiphase branch must be set to an integer greater than or equal to two." )
        end
        # checking the <val> section:
        hastag = false
        val = ""
        for elm1 in eachelement( elm )
          path1 = nodepath( elm1 )
          if elm1.name == "val"
            if hastag
              isdf_error( "More than one <val> ... </val> section have been found in the definition of the branch \"$name\". Multiple definition is not allowed.", path1 )
            end
            val = elm1.content
            if isempty( val )
              isdf_error( "The <val> ... </val> section of the branch \"name\" is empty.", path1 )
            end
            mat = Meta.parse( '['*replace( val, r"\t|\n" => "" )*']' )
            if ( size( mat, 1 ) != nph )||( size( mat, 2 ) != nph )
              isdf_error( "The size of the resistance, inductance, conductance or capacitance matrix for the branch \"$name\" is not $nph by $nph. Its size must be $nph by $nph.", path )
            end
            hastag = true
          end
        end
        if !hastag
          isdf_error( "The branch \"$name\" does not have a <val> ... </val> section. A multiphase R, L, G or C branch must have a <val> ... </val> section so as to specify its matrix.", path )
        end
        # checking the <lhs> section:
        hastag = false
        lhs = ""
        for elm1 in eachelement( elm )
          path1 = nodepath( elm1 )
          if elm1.name == "lhs"
            if hastag
              isdf_error( "More than one <lhs> ... </lhs> section have been found in the definition of the branch \"$name\". Multiple definition is not allowed.", path1 )
            end
            lhs = elm1.content
            if isempty( lhs )
              isdf_error( "The <lhs> ... </lhs> section of the branch \"name\" is empty.", path1 )
            end
            # eliminating unnecessary white spaces:
            elm1.content = replace( lhs, " " => "" )
            hastag = true
          end
        end
        if !hastag
          isdf_error( "The branch \"$name\" does not have a <lhs> ... </lhs> section. A multiphase R, L, G or C branch must have a <lhs> ... </lhs> section so as to specify the nodes to which the left-hand-side pins are connected.", path )
        end
        # checking the <rhs> section:
        hastag = false
        rhs = ""
        for elm1 in eachelement( elm )
          path1 = nodepath( elm1 )
          if elm1.name == "rhs"
            if hastag
              isdf_error( "More than one <rhs> ... </rhs> section have been found in the definition of the branch \"$name\". Multiple definition is not allowed.", path1 )
            end
            rhs = elm1.content
            if isempty( rhs )
              isdf_error( "The <rhs> ... </rhs> section of the branch \"name\" is empty.", path1 )
            end
            # eliminating unnecessary white spaces:
            elm1.content = replace( rhs, " " => "" )
            hastag = true
          end
        end
        if !hastag
          isdf_error( "The branch \"$name\" does not have a <rhs> ... </rhs> section. A multiphase R, L, G or C branch must have a <rhs> ... </rhs> section so as to specify the nodes to which the right-hand-side pins are connected.", path )
        end
      end
    else  # if NOT R, L, G or C,
      # common part for checking lhs, rhs and nph:
      if !haskey( elm, "lhs" )
        isdf_error( "The branch \"$name\" does not have a \"lhs\" attribute. A branch whose type is \"$type\" must have a \"lhs\" attribute so as to specify the node to which the left-hand-side pin is connected.", path )
      end
      lhs = elm[ "lhs" ]
      if isempty( lhs )
        isdf_error( "The \"lhs\" attribute of the branch \"$name\" is empty. It cannot be empty.", path )
      end
      # eliminating unnecessary white spaces:
      elm[ "lhs" ] = replace( lhs, " " => "" )
      if !haskey( elm, "rhs" )
        isdf_error( "The branch \"$name\" does not have a \"rhs\" attribute. A branch whose type is \"$type\" must have a \"rhs\" attribute so as to specify the node to which the right-hand-side pin is connected.", path )
      end
      rhs = elm[ "rhs" ]
      if isempty( rhs )
        isdf_error( "The \"rhs\" attribute of the branch \"$name\" is empty. It cannot be empty.", path )
      end
      # eliminating unnecessary white spaces:
      elm[ "rhs" ] = replace( rhs, " " => "" )
      if haskey( elm, "nph" )
        nph = parse( Int, elm[ "nph" ] )
        if( nph != 1 )
          isdf_error( "A value being not unity is set to the \"nph\" attribute of the branch \"$name\" whose type is \"$type\". Currently, the number of phases is restricted to unity for this type. Note that omitting the \"nph\" attribute implies that it is unity.", path )
        end
      end
      # checking NR, NL, NC and Sw:
      if ( type == "NR" )||( type == "NL" )||( type == "NC" )||( type == "Sw" )
        # checking if the content exists:
        if isempty( elm.content )
          isdf_error( "The branch \"$name\" does not have a point list to represent its nonlinear characteristic or a list of closing and opening time points. A NR (nonlinear resistor), NL (nonlinear inductor) or NC (nonlinear capacitor) must have the point list as the content of the <branch> ... </branch> section. A Sw (time-controlled switch) must have the time points as the content.", path )
        end
        if type == "Sw"
          val = Meta.parse( '('*replace( elm.content, r"\t|\n" => "" )*')' )
          if size( val, 1 ) != 1
            isdf_error( "The list of closing and opening time points for the Sw branch \"$name\" is not given as a single line. It must be given as a single-line sequence of numbers separated by commas.", path )
          end
        else  # if NR, NL or NC
          val = Meta.parse( '['*replace( elm.content, r"\t|\n" => "" )*']' )
          if size( val, 2 ) != 2
            isdf_error( "The point list to represent the nonlinear characteristic of the NR, NL or NC branch \"$name\" is not given as a two-column matrix. It must be given as a two-column matrix.", path )
          end
        end
      elseif ( type == "E" )||( type == "J" )
        # checking if the equation exists:
        if isempty( elm.content )
          isdf_error( "The branch \"$name\" does not have an  equation to determine its output. An E (voltage source) or J (current source) branch must have a Julia equation to determine its voltage or current with respect to time t as the content of the <branch> ... </branch> section.", path )
        end
      elseif ( type == "CR" )||( type == "CE" )||( type == "CJ" )
        if !haskey( elm, "ctrl" )
          isdf_error( "The branch \"$name\" does not have a \"ctrl\" attribute. A CR (controlled resistor), CE (controlled voltage source) or CJ (controlled current source) branch must have a \"ctrl\" attribute so as to specify the name of a control block which gives the resistance, voltage or current value.", path )
        end
        ctrl = elm[ "ctrl" ]
        if isempty( ctrl )
          isdf_error( "The \"ctrl\" attribute of the branch \"$name\" is empty. It cannot be empty.", path )
        end
        if !haskey( elm, "inival" )
          isdf_error( "The branch \"$name\" does not have an \"inival\" attribute. A CR (controlled resistor), CE (controlled voltage source) or CJ (controlled current source) branch must have an \"inival\" attribute so as to specify the initial value of the control quantity which is used until the first control signal is received.", path )
        end
        inival = elm[ "inival" ]
        if isempty( inval )
          isdf_error( "The \"inival\" attribute of the branch \"$name\" is empty. It cannot be empty.", path )
        end
      elseif type == "CSw"
        if !haskey( elm, "ctrl" )
          isdf_error( "The branch \"$name\" does not have a \"ctrl\" attribute. A CSw (controlled switch) branch must have a \"ctrl\" attribute so as to specify the name of a control block which gives the switch state.", path )
        end
        ctrl = elm[ "ctrl" ]
        if isempty( ctrl )
          isdf_error( "The \"ctrl\" attribute of the branch \"$name\" is empty. It cannot be empty.", path )
        end
        if !haskey( elm, "inistate" )
          isdf_error( "The branch \"$name\" does not have an \"inistate\" attribute. A CSw (controlled switch) branch must have an \"inistate\" attribute so as to specify the initial state of the switch which is used until the first control signal is received.", path )
        end
        inistate = elm[ "inistate" ]
        if isempty( instate )
          isdf_error( "The \"inistate\" attribute of the branch \"$name\" is empty. It cannot be empty.", path )
        end
      end
    end
  elseif tag == "ctrl"  # checking validity of control blocks:
    num_ctrls = num_ctrls + 1
    ord_ctrl = ordinal_string( num_ctrls )
    if !haskey( elm, "name" )
      isdf_error( "The $ord_ctrl <ctrl> ... </ctrl> section does not have a \"name\" attribute. A <ctrl> ... </ctrl> section must have a \"name\" attribute so as to identify the control block by its name.", path )
    end
    name = elm[ "name" ]
    if isempty( name )
      isdf_error( "The \"name\" attribute of the $ord_ctrl <ctrl> ... </ctrl> section is empty. It cannot be empty.", path )
    end
    # checking the <input> section:
    hastag = false
    input = ""
    for elm1 in eachelement( elm )
      path1 = nodepath( elm1 )
      if elm1.name == "input"
        if hastag
          isdf_error( "More than one <input> ... </input> section have been found in the definition of the control block \"$name\". Multiple definition is not allowed.", path1 )
        end
        input = elm1.content
        if isempty( input )
          isdf_error( "The <input> ... </input> section of the control block \"name\" is empty.", path1 )
        end
        hastag = true
      end
    end
    if !hastag
      isdf_error( "The control block \"$name\" does not have a <input> ... </input> section. A control block must have an <input> ... </input> section so as to specify the list of input variables to this control block.", path )
    end
    # checking the <Julia> section:
    hastag = false
    Julia = ""
    for elm1 in eachelement( elm )
      path1 = nodepath( elm1 )
      if elm1.name == "Julia"
        if hastag
          isdf_error( "More than one <Julia> ... </Julia> section have been found in the definition of the control block \"$name\". Multiple definition is not allowed.", path1 )
        end
        Julia = elm1.content
        if isempty( Julia )
          isdf_error( "The <Julia> ... </Julia> section of the control block \"name\" is empty.", path1 )
        end
        hastag = true
      end
    end
    if !hastag
      isdf_error( "The control block \"$name\" does not have a <Julia> ... </Julia> section. A control block must have a <Julia> ... </Julia> section so as to specify its algorithm as a Julia expression.", path )
    end
  elseif tag == "use"
    num_uses = num_uses + 1
    ord_use = ordinal_string( num_uses )
    if !haskey( elm, "component" )
      isdf_error( "The $ord_use <use> ... </use> section does not have a \"component\" attribute. A <use> ... </use> section must have a \"component\" attribute so as to identify the component to be used as a subcomponent in this <use> ... </use> section.", path )
    end
    component = elm[ "component" ]
    if isempty( component )
      isdf_error( "The \"component\" attribute of the $ord_use <use> ... </use> section is empty. It cannot be empty.", path )
    end
    if !haskey( elm, "as" )
      isdf_error( "The $ord_use <use> ... </use> section does not have an \"as\" attribute. A <use> ... </use> section must have an \"as\" attribute so as to specify the name of the instance to be generated by this <use> ... </use> section.", path )
    end
    as = elm[ "as" ]
    if isempty( as )
      isdf_error( "The \"as\" attribute of the $ord_use <use> ... </use> section is empty. It cannot be empty.", path )
    end
  end
end

return nothing

end

max_stars_repo_name: takunoda/omegasim

max_issues_repo_name: takunoda/omegasim

max_forks_repo_name: takunoda/omegasim


text:
#MAIN ENTRY POINT
#final columns for the decompressed file

#TODO: Check the portfolio sorts and the factor returns
#(replicate on 2015/12/31)

#merges crsp and compustat
# note refresh merge will be automatically run if any of the other two are run
function formfactors(;
  refreshcomp::Bool = true,
  refreshcrsp::Bool = true,
  refreshmerge::Bool = true,
  refreshfactors::Bool = true,
  datapath::String=DATA_PATH,
  factname::String = FACT_NAME,
  injlsstream::Function = IN_JLS_STREAM,
  outjlsstream::Function = OUT_JLS_STREAM,
  validatemerged::Bool = true,
  parallel::Bool = false,
  outpath::String = OUT_PATH,
  writecsv::Bool = true)::Nothing

  local univ::DataFrame
  local fact::DataFrame

  #NOTE: Need to move this to the entry point if the entry point changes
  PARALLEL[] = parallel

  #refreshfactors must be true if any of the others are true
  refreshfactors = (refreshfactors || refreshmerge || refreshcrsp ||
    refreshcomp || validatemerged)

  #makes a serialization object so we don't have to keep parsing the CSV
  if refreshfactors || (!isfile("$datapath\\$factname.jls.lz4"))
    univ = makeuniv(refreshmerge=refreshmerge,  refreshcomp=refreshcomp,
      refreshcrsp=refreshcrsp, validatemerged=validatemerged) do univ::DataFrame
        formportfolios!(univ)
    end

    fact = formfactors(univ)
    fact = mergeoutside!(fact)
    fact = boundfactors(fact)

    #NOTE: may want to push the portfolio formation into univ once it is done
    #if so, the below will be just the factor portfolio
    outjlsstream("$datapath\\$factname.jls.lz4") do s
      serialize(s, fact)
    end

  else
    injlsstream("$datapath\\$factname.jls.lz4") do s
      fact = deserialize(s) #NOTE: DEBUG ONLY
    end
  end

  writecsv && (fact |> CSV.write("$outpath\\$factname.csv"))

  return nothing
end

#assigns the group for a portfolio
function assigngroups!(vals::AbstractVector{T},
  quantiles::AbstractVector{MFloat64},
  groups::AbstractVector{MInt},
  breakpoints::NTuple) where T

  #compute the quantiles
  qdist = StatsBase.ecdf((val::T->val).(vals))
  @inbounds @simd for i ∈ 1:length(vals)
    f::Float64 = vals[i]
    quantiles[i] = qdist(f)
  end

  #compute the groups
  groups .= (q::Float64-> findfirst(b::Float64->
      q≤b, breakpoints)).(quantiles)

  return nothing
end

function formportfolios!(univ::DataFrame;
    pcs::PortfolioConstructions = PortfolioConstructions(),
    weighton::Symbol = WEIGHT_ON, returnfield::Symbol=RETURN_FIELD)

  local sunivs::GroupedDataFrame

  factorfields::Vector{Symbol} = Fvals(pcs)
  for f ∈ factorfields #might be able to do this earlier for performance
    univ = univ[(!ismissing).(univ[!,f]),:]
  end

  N::Int = size(univ,1)
  #first pre-allocate the requisite fields
  for pc ∈ pcs
    univ[!, pc.Fquantile] = Vector{MFloat64}(undef, N)
    univ[!, pc.Fgroup] = Vector{MInt}(undef, N)
    univ[!, pc.Fweight] = Vector{MFloat64}(undef, N)
  end

  sunivs = groupby(univ, :date)
  for i ∈ 1:length(sunivs)
    suniv::SubDataFrame = sunivs[i]
    for pc ∈ pcs
      assigngroups!(suniv[!,pc.Fval],
        suniv[!, pc.Fquantile], suniv[!, pc.Fgroup], pc.breakpoints)
    end
  end

  println("size: ", N)
  #group field
  groupfields = Fgroups(pcs)
  univ.returndollars = Vector{MFloat64}(undef, N)
  for suniv ∈ groupby(univ, [:date; groupfields])
    for pc ∈ pcs
      suniv.returndollars .= suniv[!,weighton] .* suniv[!, returnfield]
    end
  end

  return univ
end

#this actually computes the factor returns
function formfactors(univ::DataFrame;
    pcs::PortfolioConstructions = PortfolioConstructions(),
    weighton::Symbol = WEIGHT_ON)

  fact::DataFrame = DataFrame(date=unique(univ.date))
  sort!(fact, :date)
  Nfact::Int = size(fact,1)
  fact.retid = collect(1:Nfact)

  #pre-allocate the factor return fields
  for pc ∈ pcs
     fact[!,pc.Fval] = zeros(Float64, Nfact)
     fact[!,pc.Fnval] = zeros(Int, Nfact)
  end

  #helper function to create the portfolio name from specific sort buckets
  @inline makeportname(v::Union{AbstractVector,NTuple})::Symbol = Symbol(
    "port_", join((string).(v), "_"))

  @inline makenportname(v::Union{AbstractVector,NTuple})::Symbol = (
        Symbol("N_", makeportname(v)))
  @inline makenportname(portname::Symbol) = Symbol("N_", portname)

  #now make the portfolio sorts
  sorts::Vector{NTuple{pcs.N, Int}} = computesorts(pcs)
  portnames::Vector{Symbol} = (t::NTuple->makeportname(t)).(sorts)
  for s ∈ portnames
    fact[!,s] = Vector{MFloat64}(undef, Nfact)
    fact[!,makenportname(s)] = Vector{MInt}(undef, Nfact)
  end

  #index the rows by date
  rfacts::Dict = Dict(r.date => r for r ∈ eachrow(fact))

  #now compute the returns
  groupfields::Vector{Symbol} = Fgroups(pcs)
  for suniv ∈ groupby(univ, :date)
    rfact = rfacts[suniv.date[1]]

    #split each date by triple sort
    for ssuniv ∈ groupby(suniv, groupfields)

      local portname::Symbol = makeportname((f->ssuniv[1,f]).(groupfields))
      local nportname::Symbol = makenportname(portname)

      #this is the overall portfolio return
      portret::Float64 = sum(ssuniv.returndollars) / sum(ssuniv[!,weighton])
      rfact[portname] = portret
      rfact[nportname] = size(ssuniv,1)
      #println("portret: $portret")

    end
  end

  #since we have the returns for each of the sorts on all dates, forming
  #the factors is just a matter of weighting the returns

  #iterate through each of the sort portfolos
  for coord ∈ sorts
    local portname::Symbol = makeportname(coord)
    local nportname::Symbol = makenportname(portname)

    #apply the appropriate weight on the sort portfolio to each factor
    for i ∈ 1:length(coord)
      pc::PortfolioConstruction = pcs[i]
      groupweight::Float64 = pc.groupweights[coord[i]]

      #compute the factor contribution
      fact[!, pc.Fval] .+= groupweight .* fact[!, portname]
      fact[!, pc.Fnval] .+= (fact[!, nportname]) .* (groupweight>0)

    end
  end

  return fact
end

boundfactors(fact::DataFrame;
    mindate::Date = MIN_DATE,
    maxdate::Date = MAX_DATE) = fact[(d->(d>mindate) && (d<maxdate)).(fact.date),:]

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
"""
    earthmover(a, b, C)

Compute the optimal transport map
```math
argmin_{P ∈ U(a, b)} ⟨C, P⟩
```
and the dual potentials using the network simplex method.

# References

Peyré, G., & Cuturi, M.. (2018). Computational Optimal Transport. [arXiv:1803.00567](https://arxiv.org/abs/1803.00567).
"""
function earthmover(
    a::AbstractVector{<:Real},
    b::AbstractVector{<:Real},
    C::AbstractMatrix{<:Real}
)
    P = SparseArrays.spzeros(Base.promote_eltype(a, b), length(a), length(b))
    earthmover!(P, a, b, C)
end

function earthmover!(P, a, b, C)
    # find a feasible initial transport map
    initialmap!(P, a, b)

    # construct the initial graph
    n = length(a)
    m = length(b)
    graph = LightGraphs.SimpleGraph(n + m)
    for ij in findall(!iszero, P)
        i, j = Tuple(ij)

        # convert to vertices
        u = i
        v = j + n

        # add new edge
        LightGraphs.add_edge!(graph, u, v)
    end

    # compute the dual potentials
    T = Base.promote_eltype(P, C)
    f = Vector{T}(undef, n)
    g = Vector{T}(undef, m)
    dualpotentials!!(f, g, C, graph)

    while true
        # check if the dual potentials are feasible
        res = findfirstinfeasible(f, g, C)
        res === nothing && break

        # obtain vertices of the edge that will be added
        i, j = res
        u = i
        v = j + n

        # check if these vertices are already connected
        path = LightGraphs.a_star(graph, v, u)

        if !isempty(path)
            # determine the largest possible increase and which edge to remove
            e = path[1]
            θmin = P[LightGraphs.dst(e), LightGraphs.src(e) - n]
            kmin = 1
            for k in 3:2:length(path)
                e = path[k]
                θ = P[LightGraphs.dst(e), LightGraphs.src(e) - n]
                if θ < θmin
                    θmin = θ
                    kmin = k
                end
            end

            # update the transport map
            P[i, j] += θmin
            for k in 1:length(path)
                e = path[k]
                if isodd(k)
                    # try to avoid floating point errors
                    if k == kmin
                        P[LightGraphs.dst(e), LightGraphs.src(e) - n] = 0
                    else
                        P[LightGraphs.dst(e), LightGraphs.src(e) - n] -= θmin
                    end
                else
                    P[LightGraphs.src(e), LightGraphs.dst(e) - n] += θmin
                end
            end

            # remove the replaced edge
            LightGraphs.rem_edge!(graph, path[kmin])
        end

        # add the new edge
        LightGraphs.add_edge!(graph, u, v)

        # recompute the dual potentials
        dualpotentials!!(f, g, C, graph)
    end

    P, f, g
end

function findfirstinfeasible(f, g, C)
    n = length(f)
    m = length(g)

    @inbounds for j in 1:m, i in 1:n
        f[i] + g[j] > C[i, j] && return (i, j)
    end

    nothing
end

function dualpotentials!!(f, g, C, graph)
    n = length(f)
    m = length(g)

    visited = falses(n+m)
    S = Int[]
    @inbounds for v in LightGraphs.vertices(graph)
        # do not revisit processed vertices
        visited[v] && continue

        # handle first vertex in a new tree
        if v > n
            g[v - n] = 0
        else
            f[v] = 0
        end
        push!(S, v)
        visited[v] = true

        # handle all other vertices in this tree by DFS
        while !isempty(S)
            s = S[end]
            u = 0
            for t in LightGraphs.outneighbors(graph, s)
                if !visited[t]
                    u = t
                    if t > n
                        g[t - n] = C[s, t - n] - f[s]
                    else
                        f[t] = C[t, s - n] - g[s - n]
                    end
                    break
                end
            end
            if u == 0
                pop!(S)
            else
                visited[u] = true
                push!(S, u)
            end
        end
    end

    nothing
end

function initialmap!(P, a, b)
    n = length(a)
    m = length(b)
    size(P) == (n, m) ||
        throw(DimensionMismatch("dimensions of histograms and transport map do not match"))

    # initialize row and column sums
    T = eltype(P)
    r = T(a[1])
    c = T(b[1])

    i = j = 1
    while i ≤ n && j ≤ m
        # update initial map
        t = min(r, c)
        P[i, j] = t

        # update r and c
        r -= t
        c -= t

        if iszero(r)
            i += 1
            if i ≤ n
                r = T(a[i])
            end
        end

        if iszero(c)
            j += 1
            if j ≤ m
                c = T(b[j])
            end
        end
    end

    P
end


text:
##################################################################
# Abstract types

"""
$(TYPEDEF)
Apex type for grid components.
"""
abstract type AbstractGridComponent <: AbstractExtendableGridApexType end

"""
$(TYPEDEF)
1D Array of floating point data
"""
abstract type AbstractGridFloatArray1D <: AbstractGridComponent end

"""
$(TYPEDEF)
2D Array of floating point data
"""
abstract type AbstractGridFloatArray2D <: AbstractGridComponent end

"""
$(TYPEDEF)
1D Array of interger data
"""
abstract type AbstractGridIntegerArray1D <: AbstractGridComponent end

"""
$(TYPEDEF)
2D Array of integer data
"""
abstract type AbstractGridIntegerArray2D <: AbstractGridComponent end

"""
$(TYPEDEF)
Integer number
"""
abstract type AbstractGridIntegerConstant <: AbstractGridComponent end

"""
$(TYPEDEF)
Floating point number
"""
abstract type AbstractGridFloatConstant <: AbstractGridComponent end

"""
$(TYPEDEF)

Any kind of adjacency between grid components
"""
abstract type AbstractGridAdjacency <: AbstractGridComponent end

"""
$(TYPEDEF)

Array of element geometry information.
"""
abstract type AbstractElementGeometries <: AbstractGridComponent end

"""
$(TYPEDEF)

Array of element region number information.
"""
abstract type AbstractElementRegions <: AbstractGridComponent end

##################################################################
# Component key types

"""
$(TYPEDEF)

Node coordinates
"""
abstract type Coordinates <: AbstractGridFloatArray2D end

"""
$(TYPEDEF)

Adjacency describing nodes per grid cell
"""
abstract type CellNodes <: AbstractGridAdjacency end

"""
$(TYPEDEF)

Adjacency describing nodes per grid boundary face
"""
abstract type BFaceNodes <: AbstractGridAdjacency end

"""
$(TYPEDEF)

Description of cell geometries
"""
abstract type CellGeometries <: AbstractElementGeometries end

"""
Description of boundary face geometries

$(TYPEDEF)
"""
abstract type BFaceGeometries <: AbstractElementGeometries end

"""
$(TYPEDEF)

Cell region number per cell
"""
abstract type CellRegions <: AbstractElementRegions end

"""
$(TYPEDEF)

Boundary region number per boundary face
"""
abstract type BFaceRegions <: AbstractElementRegions end

"""
$(TYPEDEF)

Number of cell regions
"""
abstract type NumCellRegions <: AbstractGridIntegerConstant end

"""
$(TYPEDEF)

Number of boundary face regions
"""
abstract type NumBFaceRegions <: AbstractGridIntegerConstant end

"""
$(TYPEDEF)

Coordinate system
"""
abstract type CoordinateSystem <: AbstractGridComponent end

"""
$(TYPEDEF)

Boundary edge region number per boundary edge
"""
abstract type BEdgeRegions <: AbstractElementRegions end

"""
$(TYPEDEF)

Number of boundary edge regions
"""
abstract type NumBEdgeRegions <: AbstractGridIntegerConstant end

############################################################
# Grid type

"""
$(TYPEDEF)
Grid type wrapping Dict
"""
mutable struct ExtendableGrid{Tc,Ti}
    components::Dict{Type{<:AbstractGridComponent},Any}
    ExtendableGrid{Tc,Ti}() where{Tc,Ti} =new(Dict{Type{<:AbstractGridComponent},Any}())
end

#############################################################
#
"""
````
const ElementInfo{T}=Union{Vector{T},VectorOfConstants{T}}`
````

Union type for element information arrays. If all elements have
the same information, it can be stored in an economical form
as a [`VectorOfConstants`](@ref).
"""
const ElementInfo{T}=Union{Vector{T},VectorOfConstants{T}}

############################################################
# Generic component methods

"""
$(TYPEDSIGNATURES)
Default veryform method.

"veryform"  means "verify and/or transform"  and is called to check
and possibly transform components to be added to the grid via `setindex!`.

The default method just passes data through.
"""
veryform(grid::ExtendableGrid,v,::Type{<:AbstractGridComponent})=v

"""
$(TYPEDSIGNATURES)

Set new grid component
"""
Base.setindex!(grid::ExtendableGrid,v,T::Type{<:AbstractGridComponent})= grid.components[T]=veryform(grid,v,T)

"""
$(TYPEDSIGNATURES)

To be called by getindex. This triggers lazy creation of
non-existing gridcomponents
"""
Base.get!(grid::ExtendableGrid,T::Type{<:AbstractGridComponent})= get!( ()->veryform(grid,instantiate(grid,T),T),  grid.components,T)

"""
$(TYPEDSIGNATURES)

"Hook" for methods instantiating lazy components.
"""
function instantiate end

"""
````
Base.getindex(grid::ExtendableGrid,T::Type{<:AbstractGridComponent})
````

Generic method for obtaining grid component.

This method is mutating in the sense that non-existing grid components
are created on demand.

Due to the fact that components are stored as Any the return
value triggers type instability.
"""
Base.getindex(grid::ExtendableGrid,T::Type{<:AbstractGridComponent})=get!(grid,T)

############################################################
# Verfication of inserted data

"""
````
veryform(grid::ExtendableGrid{Tc,Ti},v,T::Type{<:AbstractGridAdjacency}) where{Tc,Ti}
````

Check proper type of adjacencies upon insertion
"""
veryform(grid::ExtendableGrid{Tc,Ti},v,T::Type{<:AbstractGridAdjacency}) where{Tc,Ti}= typeof(v)<:Adjacency{Ti} ? v : throw("Type mismatch")

############################################################
# Type stable component access

"""
````
Base.getindex(grid::ExtendableGrid{Tc,Ti},T::Type{<:AbstractGridAdjacency})::Adjacency{Ti} where{Tc,Ti}
````

Type stable return of adjacency component
"""
function Base.getindex(grid::ExtendableGrid{Tc,Ti},T::Type{<:AbstractGridAdjacency})::Adjacency{Ti} where{Tc,Ti}
    get!(grid,T)
end

"""
````
Base.getindex(grid::ExtendableGrid{Tc,Ti},T::Type{<:AbstractGridFloatArray2D})::Array{Tc,2}
````

Type stable c method to obtain 2D array from grid
"""
function Base.getindex(grid::ExtendableGrid{Tc,Ti},T::Type{<:AbstractGridFloatArray2D})::Array{Tc,2} where{Tc,Ti}
    get!(grid,T)
end

"""
````
Base.getindex(grid::ExtendableGrid{Tc,Ti},T::Type{<:AbstractGridFloatArray1D})::Array{Tc,1} where{Tc,Ti}
````

Type stable method to obtain 1D array from grid
"""
function Base.getindex(grid::ExtendableGrid{Tc,Ti},T::Type{<:AbstractGridFloatArray1D})::Array{Tc,1} where{Tc,Ti}
    get!(grid,T)
end

"""
````
Base.getindex(grid::ExtendableGrid{Tc,Ti},T::Type{<:AbstractElementGeometries})::ElementInfo{DataType}
````

Type stable method to obtain  element type
"""
function Base.getindex(grid::ExtendableGrid{Tc,Ti},T::Type{<:AbstractElementGeometries})::ElementInfo{DataType} where{Tc,Ti}
    get!(grid,T)
end

"""
````
Base.getindex(grid::ExtendableGrid{Tc,Ti},T::Type{<:AbstractElementRegions})::ElementInfo{Ti} where{Tc,Ti}
````

Type stable method to obtain element region number
"""
function Base.getindex(grid::ExtendableGrid{Tc,Ti},T::Type{<:AbstractElementRegions})::ElementInfo{Ti} where{Tc,Ti}
    get!(grid,T)
end

"""
````
Base.getindex(grid::ExtendableGrid{Tc,Ti},T::Type{<:AbstractGridIntegerConstant})::Ti where{Tc,Ti}
````

Type stable method to obtain  integer constant stored on elements
"""
function Base.getindex(grid::ExtendableGrid{Tc,Ti},T::Type{<:AbstractGridIntegerConstant})::Ti where{Tc,Ti}
    get!(grid,T)
end

"""
````
Base.getindex(grid::ExtendableGrid{Tc,Ti},T::Type{<:AbstractGridFloatConstant})::Tc where{Tc,Ti}
````

Type stable  method to obtain float constant stored on elements
"""
function Base.getindex(grid::ExtendableGrid{Tc,Ti},T::Type{<:AbstractGridFloatConstant})::Tc where{Tc,Ti}
    get!(grid,T)
end

############################################################
# Instantiation methods

"""
$(TYPEDSIGNATURES)

Instantiate number of cell regions
"""
instantiate(grid, ::Type{NumCellRegions})=maximum(grid[CellRegions])

"""
$(TYPEDSIGNATURES)

Instantiate number of bface regions
"""
instantiate(grid, ::Type{NumBFaceRegions})=maximum(grid[BFaceRegions])

"""
$(TYPEDSIGNATURES)

Instantiate number of boundary edge regions
"""
instantiate(grid, ::Type{NumBEdgeRegions})=maximum(grid[BEdgeRegions])

function prepare_bedgeregions!(grid::ExtendableGrid)
    bedges             = grid[BEdgeNodes]
    bedgeregions       = zeros(Int32, num_bedges(grid))
    grid[BEdgeRegions] = bedgeregions
end

instantiate(grid, ::Type{BEdgeRegions})=prepare_bedgeregions!(grid)

#############################################################
# General methods
"""
$(TYPEDSIGNATURES)

Print the hierarchy of grid component key types (subtypes of [`AbstractGridComponent`](@ref).
This includes additionally user defined subptypes.
"""
gridcomponents()=AbstractTrees.print_tree(AbstractGridComponent,5,indicate_truncation=false)

"""
$(TYPEDSIGNATURES)

Keys in grid
"""
Base.keys(g::ExtendableGrid)=Base.keys(g.components)

"""
$(TYPEDSIGNATURES)

Check if key is in grid
"""
Base.haskey(g::ExtendableGrid,k) = Base.haskey(g.components,k)

"""
$(SIGNATURES)

Space dimension of grid
"""
dim_space(grid::ExtendableGrid)= size(grid[Coordinates],1)

"""
$(SIGNATURES)

Grid dimension dimension of grid (larges element dimension)
"""
dim_grid(grid::ExtendableGrid)=  dim_element(grid[CellGeometries][1])

"""
$(SIGNATURES)

Number of nodes in grid
"""
num_nodes(grid::ExtendableGrid)::Int= size(grid[Coordinates],2)

"""
$(TYPEDSIGNATURES)

Number of cells in grid
"""
num_cells(grid::ExtendableGrid)= num_sources(grid[CellNodes])

"""
$(TYPEDSIGNATURES)

Number of boundary faces in grid.
"""
num_bfaces(grid::ExtendableGrid)= haskey(grid,BFaceNodes) ? num_sources(grid[BFaceNodes]) : 0

"""
$(TYPEDSIGNATURES)

Number of boundary edges in grid.
"""
num_bedges(grid::ExtendableGrid)= haskey(grid,BEdgeNodes) ? num_sources(grid[BEdgeNodes]) : 0

"""
$(TYPEDSIGNATURES)

Maximum  cell  region number
"""
num_cellregions(grid::ExtendableGrid)=grid[NumCellRegions]

"""
$(TYPEDSIGNATURES)

Maximum  boundary face region numbers
"""
num_bfaceregions(grid::ExtendableGrid)=grid[NumBFaceRegions]

"""
$(TYPEDSIGNATURES)

Maximum boundary edge region numbers
"""
num_bedgeregions(grid::ExtendableGrid)=grid[NumBEdgeRegions]

"""
$(SIGNATURES)
Type of coordinates in grid
"""
coord_type(grid::ExtendableGrid{Tc, Ti}) where {Tc,Ti}=Tc

"""
$(SIGNATURES)

Type of indices
"""
index_type(grid::ExtendableGrid{Tc, Ti}) where {Tc,Ti}=Ti

"""
$(TYPEDSIGNATURES)

Map a function onto node coordinates of grid
"""
function Base.map(f::Function, grid::ExtendableGrid)
    coord=grid[Coordinates]
    dim=dim_space(grid)
    if dim==1
        @views Base.map(f,coord[1,:])
    elseif dim==2
        @views Base.map(f,coord[1,:], coord[2,:])
    else
        @views Base.map(f,coord[1,:], coord[2,:], coord[3,:])
    end
end

function Base.show(io::IO,grid::ExtendableGrid)
    if num_edges(grid)>0
        str=@sprintf("%s;\ndim: %d nodes: %d cells: %d bfaces: %d, edges: %d\n",
                     typeof(grid),dim_space(grid),num_nodes(grid), num_cells(grid), num_bfaces(grid), num_edges(grid))
    else
        str=@sprintf("%s;\ndim: %d nodes: %d cells: %d bfaces: %d\n",
                     typeof(grid),dim_space(grid),num_nodes(grid), num_cells(grid), num_bfaces(grid))
    end
    println(io,str)
end

"""
$(SIGNATURES)

Recursively check seeming equality of two grids. Seemingly means
that long arrays are only compared via random samples
"""
function seemingly_equal(grid1::ExtendableGrid, grid2::ExtendableGrid)
    for key in keys(grid1)
        if !haskey(grid2,key)
            return false
        end
        if !seemingly_equal(grid1[key],grid2[key])
            return false
        end
    end
    return true
end

"""
$(SIGNATURES)

Check for seeming equality of two arrays by random sample.
"""
function seemingly_equal(array1::AbstractArray, array2::AbstractArray)
    if size(array1)!=size(array2)
        return false
    end
    l=length(array1)
    ntests=Float64(min(l,50))
    p=min(0.5,ntests/l)
    for i in randsubseq( 1:l,p)
        if !seemingly_equal(array1[i],array2[i])
            return false
        end
    end
    true
end

seemingly_equal(a1::VariableTargetAdjacency, a2::VariableTargetAdjacency)=seemingly_equal(a1.colentries,a2.colentries)&& seemingly_equal(a1.colstart, a2.colstart)
seemingly_equal(x1::Type, x2::Type)=(x1==x2)
seemingly_equal(x1::Number, x2::Number)=(x1≈x2)
seemingly_equal(x1::Any, x2::Any)=(x1==x2)

Base.extrema(grid::ExtendableGrid)=Base.extrema(grid[Coordinates],dims=2)

function bbox(grid)
    e=extrema(grid)
    map(a->a[1],e),map(a->a[2],e)
end


text:
__precompile__(true)

module MultivariateAnomalies
#######################
  #import MultivariateStats
  import Combinatorics
  using LinearAlgebra
  using Distances
  using LIBSVM
  using Base.Cartesian
  using Statistics
  using StatsBase

  # DistDensity
  export
    dist_matrix,
    dist_matrix!,
    init_dist_matrix,
    knn_dists,
    init_knn_dists,
    knn_dists!,
    kernel_matrix,
    kernel_matrix!,
# Detection Algorithms
    REC,
    init_REC,
    REC!,
    KDE,
    init_KDE,
    KDE!,
    init_T2,
    T2,
    T2!,
    KNN_Gamma,
    init_KNN_Gamma,
    KNN_Gamma!,
    KNN_Delta,
    KNN_Delta!,
    init_KNN_Delta,
    UNIV,
    UNIV!,
    init_UNIV,
    SVDD_train,
    SVDD_predict,
    KNFST_train,
    KNFST_predict,
    KNFST_predict!,
    init_KNFST,
    Dist2Centers,
# FeatureExtraction
    sMSC,
    TDE,
    mw_VAR,
    mw_VAR!,
    mw_AVG,
    mw_AVG!,
    mw_COR,
    EWMA,
    EWMA!,
    get_MedianCycles,
    get_MedianCycle,
    get_MedianCycle!,
    init_MedianCycle,
    mapMovingWindow,
# AUC
    auc,
# Scores
    get_quantile_scores,
    get_quantile_scores!,
    compute_ensemble,
# high level functions
    getParameters,
    detectAnomalies,
    init_detectAnomalies,
    detectAnomalies!,
# online algorithms
    Euclidean_distance!,
    Mahalanobis_distance!,
    SigmaOnline!,
    KDEonline!,
    KNNonline!,
    REConline!

# Distance and Density Estimation
include("DistDensity.jl")
# Multivariate Anomaly Detection Algorithms
include("DetectionAlgorithms.jl")
# AUC computations
include("AUC.jl")
# Feature Extraction techniques
include("FeatureExtraction.jl")
# post processing for Anomaly Scores
include("Scores.jl")
# high level functions
include("HighLevelFunctions.jl")
# online algorithms
include("OnlineAlgorithms.jl")

#######################
end


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

# Transformation
# --------------

# Logarithmic transformtion
struct LogT
    pseudocount::Float32
end
function transform!(t::LogT, X::Matrix{Float32})
    @fastmath @inbounds for i in 1:length(X)
        X[i] = log(X[i] + t.pseudocount)
    end
    return X
end

# Freeman-Tukey transformation
struct FTT end
function transform!(::FTT, X::Matrix{Float32})
    @fastmath @inbounds for i in 1:length(X)
        X[i] = sqrt(X[i]) + sqrt(X[i]+1.0f0)
    end
    return X
end

# Dimensionality Reduction
# ------------------------

mutable struct PCA
    dims::Int
    randomize::Bool
    proj::Matrix{Float32}
    PCA(dims::Integer; randomize::Bool=true) = new(dims, randomize)
end

# NOTE: This assumes features are already centered (if needed).
function fit!(pca::PCA, X::AbstractMatrix)
    if pca.randomize
        pca.proj, = rsvd(X, pca.dims)
    else
        pca.proj, = tsvd(X, pca.dims)
    end
    return pca
end
reducedims(pca::PCA, X::AbstractMatrix) = pca.proj'X

# Preprocessor
# ------------

# permutation, log transformation, scaling, projection, etc.
struct Preprocessor
    n::Int
    featurenames::Vector{String}
    sums::Vector{Float32}
    soss::Vector{Float32}
    transformer::Union{LogT,FTT}
    dimreducer::PCA
    normalize::Bool
    standardize::Bool
    scalefactor::Float32
end

function Preprocessor(
        Y::ExpressionMatrix,
        transformer,
        dimreducer,
        normalize::Bool,
        standardize::Bool,
        scalefactor::Float32)
    @assert Y.data isa Matrix{Float32}
    X::Matrix{Float32} = copy(Y.data)
    # normalize cell counts
    if normalize
        X .*= scalefactor ./ sum(X, dims=1)
    end
    # apply variance-stabilizing transformation
    transform!(transformer, X)
    # standardize features
    m, n = size(X)
    sums = vec(sum(X, dims=2))
    soss = vec(sum(X.^2, dims=2))
    μ = sums ./ n
    X .-= μ
    if standardize
        σ = sqrt.((soss .- 2 .* μ .* sums) ./ n .+ μ.^2)
        if any(σ .== 0)
            throw(ArgumentError("found $(sum(σ .== 0)) features with no variance; filter out these features"))
        end
        X .*= inv.(σ)
    end
    # fit dimreducer
    fit!(dimreducer, X)
    return Preprocessor(
        n,
        Y.featurenames,
        sums,
        soss,
        transformer,
        dimreducer,
        normalize,
        standardize,
        scalefactor,
    )
end

function preprocess(proc::Preprocessor, Y::ExpressionMatrix, inferstats::Symbol)
    perm = zeros(Int, length(proc.featurenames))
    for (i, name) in enumerate(proc.featurenames)
        perm[i] = get(Y.featuremap, name, 0)
    end
    X = permuterows(perm, Y.data)::Matrix{Float32}  # should be inferable
    if proc.normalize
        X .*= proc.scalefactor ./ sum(X, dims=1)
    end
    transform!(proc.transformer, X)
    m, n = size(X)
    @assert length(proc.sums) == length(proc.soss) == m
    if inferstats == :query
        μ = vec(mean(X, dims=2))
        σ = vec(std(X, dims=2))
        @inbounds for i in 1:m
            # insert the std of the database if no variability in the query cells
            if σ[i] == 0
                σ[i] = sqrt((proc.soss[i] - 2 * μ[i] * proc.sums[i]) / proc.n + μ[i]^2)
            end
        end
    elseif inferstats == :database
        μ = proc.sums ./ proc.n
        σ = sqrt.((proc.soss .- 2 .* μ .* proc.sums) ./ proc.n .+ μ.^2)
    else
        @assert inferstats == :both
        sums = vec(sum(X, dims=2))
        soss = vec(sum(X.^2, dims=2))
        μ = zeros(Float32, m)
        σ = zeros(Float32, m)
        @inbounds for i in 1:m
            μ[i], σ[i] = mean_and_std(
                proc.sums[i], proc.soss[i], proc.n,
                sums[i], soss[i], n,
            )
        end
    end
    if proc.standardize
        invstd = inv.(σ)
        @inbounds for j in 1:n, i in 1:m
            X[i,j] = (X[i,j] - μ[i]) * invstd[i]
        end
    else
        X .-= μ
    end
    return reducedims(proc.dimreducer, X)
end

# sos: sum of squares
@inline function mean_and_std(
        sum1::Float32, sos1::Float32, n1::Integer,
        sum2::Float32, sos2::Float32, n2::Integer,
    )
    sum = sum1 + sum2
    d = inv(n1 + n2)
    μ = sum * d
    σ = sqrt((sos1 + sos2 - 2 * μ * sum) * d + μ^2)
    return μ, σ
end

indims(p::Preprocessor) = length(p.featurenames)
outdims(p::Preprocessor) = p.dimreducer.dims

function permuterows(perm::Vector{Int}, Y::AbstractMatrix)
    m = length(perm)
    n = size(Y, 2)
    Y_permuted = Matrix{Float32}(undef, m, n)
    @inbounds for j in 1:n, i in 1:m
        Y_permuted[i,j] = perm[i] > 0 ? Y[perm[i],j] : 0
    end
    return Y_permuted
end


text:
using Distributed
using JLD2
using ProgressMeter

@time @everywhere include("../BioClim/src/required.jl")
cd("/Users/kiristern/Documents/GitHub/SDM")
## Get & prepare data
@time @everywhere begin
    # Load data from CSV files
    pred_prey = CSV.read("data/pred_prey/vulpes_vulpes_america.csv", header=true)
     # Prepare data (select columns, arrange values)
    df = prepare_gbif_data(pred_prey)
    # Separate species
    taxa_occ = [df[df.species .== u,:] for u in unique(df.species)]
    #define predator coordinate range using min/max ranges
    long_predator = extrema(taxa_occ[1].longitude)
    lat_predator = extrema(taxa_occ[1].latitude)
    #define coordinate range for all species using min/max ranges
    long_all = extrema(df.longitude)
    lat_all = extrema(df.latitude)
end

#take worldclim values for all coordinate range
@time wc_vars = pmap(x -> worldclim(x)[long_all, lat_all], 1:19)
#make predictions for all species
@time predictions = @showprogress pmap(x -> species_bclim(x, wc_vars), taxa_occ)

#new predictions of species only within predator coordinate range
new_predictions = [p[long_predator, lat_predator] for p in predictions]
##another way to write loop
#new_predictions = [predictions[i][long_predator, lat_predator] for i in 1:4]

# convert year "string" to "integer"
#df.longitude = parse.(Float64, df.longitude)

# Get the worldclim data by their layer number
cd("/Users/kiristern/Documents/GitHub/BioClim/")

# #prediction pred
# @info "Extract and crop bioclim variables"
# @time wc_vars = [clip(worldclim(i),  taxa_occ[1]) for i in 1:19];
# # Make the prediction for each layer
# @info "Predictions for each layer"
# @time predictions = [bioclim(wc_vars[i], taxa_occ[1]) for i in 1:length(wc_vars)];
# # Make the final prediction by taking the minimum
# @info "Minimum-consensus aggregation"
# @time prediction1 = reduce(minimum, predictions);
# # Get the threshold for NaN given a percentile
# @info "Threshold estimation"
# @time threshold = first(quantile(prediction1[taxa_occ[1]], [0.05]))
# @info "5% threshold:\t$(round(threshold; digits=3))"
# # Filter the predictions based on the threshold
# @info "Final prediction filtering"
# @time for i in eachindex(prediction1.grid)
#     prediction1.grid[i] < threshold && (prediction1.grid[i] = NaN)
# end
#
# #prediction prey1
# @info "Extract and crop bioclim variables"
# @time wc_vars = [clip(worldclim(i), taxa_occ[2]) for i in 1:19];
# # Make the prediction for each layer
# @info "Predictions for each layer"
# @time predictions = [bioclim(wc_vars[i], taxa_occ[2]) for i in 1:length(wc_vars)];
# # Make the final prediction by taking the minimum
# @info "Minimum-consensus aggregation"
# @time prediction2 = reduce(minimum, predictions);
# # Get the threshold for NaN given a percentile
# @info "Threshold estimation"
# @time threshold = first(quantile(prediction2[taxa_occ[2]], [0.05]))
# @info "5% threshold:\t$(round(threshold; digits=3))"
# # Filter the predictions based on the threshold
# @info "Final prediction filtering"
# @time for i in eachindex(prediction2.grid)
#     prediction2.grid[i] < threshold && (prediction2.grid[i] = NaN)
# end
#
# @info "Extract and crop bioclim variables"
# @time wc_vars = [clip(worldclim(i), taxa_occ[3]) for i in 1:19];
# # Make the prediction for each layer
# @info "Predictions for each layer"
# @time predictions = [bioclim(wc_vars[i], taxa_occ[3]) for i in 1:length(wc_vars)];
# # Make the final prediction by taking the minimum
# @info "Minimum-consensus aggregation"
# @time prediction3 = reduce(minimum, predictions);
# # Get the threshold for NaN given a percentile
# @info "Threshold estimation"
# @time threshold = first(quantile(prediction3[taxa_occ[3]], [0.05]))
# @info "5% threshold:\t$(round(threshold; digits=3))"
# # Filter the predictions based on the threshold
# @info "Final prediction filtering"
# @time for i in eachindex(prediction3.grid)
#     prediction3.grid[i] < threshold && (prediction3.grid[i] = NaN)
# end
#
# @info "Extract and crop bioclim variables"
# @time wc_vars = [clip(worldclim(i), taxa_occ[4]) for i in 1:19];
# # Make the prediction for each layer
# @info "Predictions for each layer"
# @time predictions = [bioclim(wc_vars[i], taxa_occ[4]) for i in 1:length(wc_vars)];
# # Make the final prediction by taking the minimum
# @info "Minimum-consensus aggregation"
# @time prediction4 = reduce(minimum, predictions);
# # Get the threshold for NaN given a percentile
# @info "Threshold estimation"
# @time threshold = first(quantile(prediction4[taxa_occ[4]], [0.05]))
# @info "5% threshold:\t$(round(threshold; digits=3))"
# # Filter the predictions based on the threshold
# @info "Final prediction filtering"
# @time for i in eachindex(prediction4.grid)
#     prediction4.grid[i] < threshold && (prediction4.grid[i] = NaN)
# end

# @info "Extract and crop bioclim variables"
# @time wc_vars = [clip(worldclim(i), taxa_occ[5]) for i in 1:19];
# # Make the prediction for each layer
# @info "Predictions for each layer"
# @time predictions = [bioclim(wc_vars[i], taxa_occ[5]) for i in 1:length(wc_vars)];
# # Make the final prediction by taking the minimum
# @info "Minimum-consensus aggregation"
# @time prediction5 = reduce(minimum, predictions);
# # Get the threshold for NaN given a percentile
# @info "Threshold estimation"
# @time threshold = first(quantile(prediction5[taxa_occ[5]], [0.05]))
# @info "5% threshold:\t$(round(threshold; digits=3))"
# # Filter the predictions based on the threshold
# @info "Final prediction filtering"
# @time for i in eachindex(prediction5.grid)
#     prediction5.grid[i] < threshold && (prediction5.grid[i] = NaN)
# end
#
# @info "Extract and crop bioclim variables"
# @time wc_vars = [clip(worldclim(i), taxa_occ[6]) for i in 1:19];
# # Make the prediction for each layer
# @info "Predictions for each layer"
# @time predictions = [bioclim(wc_vars[i], taxa_occ[6]) for i in 1:length(wc_vars)];
# # Make the final prediction by taking the minimum
# @info "Minimum-consensus aggregation"
# @time prediction6 = reduce(minimum, predictions);
# # Get the threshold for NaN given a percentile
# @info "Threshold estimation"
# @time threshold = first(quantile(prediction6[taxa_occ[6]], [0.05]))
# @info "5% threshold:\t$(round(threshold; digits=3))"
# # Filter the predictions based on the threshold
# @info "Final prediction filtering"
# @time for i in eachindex(prediction6.grid)
#     prediction6.grid[i] < threshold && (prediction6.grid[i] = NaN)
# end
#
# @info "Extract and crop bioclim variables"
# @time wc_vars = [clip(worldclim(i), taxa_occ[7]) for i in 1:19];
# # Make the prediction for each layer
# @info "Predictions for each layer"
# @time predictions = [bioclim(wc_vars[i], taxa_occ[7]) for i in 1:length(wc_vars)];
# # Make the final prediction by taking the minimum
# @info "Minimum-consensus aggregation"
# @time prediction7 = reduce(minimum, predictions);
# # Get the threshold for NaN given a percentile
# @info "Threshold estimation"
# @time threshold = first(quantile(prediction7[taxa_occ[7]], [0.05]))
# @info "5% threshold:\t$(round(threshold; digits=3))"
# # Filter the predictions based on the threshold
# @info "Final prediction filtering"
# @time for i in eachindex(prediction7.grid)
#     prediction7.grid[i] < threshold && (prediction7.grid[i] = NaN)
# end
#
# @info "Extract and crop bioclim variables"
# @time wc_vars = [clip(worldclim(i), taxa_occ[8]) for i in 1:19];
# # Make the prediction for each layer
# @info "Predictions for each layer"
# @time predictions = [bioclim(wc_vars[i], taxa_occ[8]) for i in 1:length(wc_vars)];
# # Make the final prediction by taking the minimum
# @info "Minimum-consensus aggregation"
# @time prediction8 = reduce(minimum, predictions);
# # Get the threshold for NaN given a percentile
# @info "Threshold estimation"
# @time threshold = first(quantile(prediction8[taxa_occ[8]], [0.05]))
# @info "5% threshold:\t$(round(threshold; digits=3))"
# # Filter the predictions based on the threshold
# @info "Final prediction filtering"
# @time for i in eachindex(prediction8.grid)
#     prediction8.grid[i] < threshold && (prediction8.grid[i] = NaN)
# end

# @info "Extract and crop bioclim variables"
# @time wc_vars = [clip(worldclim(i), taxa_occ[9]) for i in 1:19];
# # Make the prediction for each layer
# @info "Predictions for each layer"
# @time predictions = [bioclim(wc_vars[i], taxa_occ[9]) for i in 1:length(wc_vars)];
# # Make the final prediction by taking the minimum
# @info "Minimum-consensus aggregation"
# @time prediction9 = reduce(minimum, predictions);
# # Get the threshold for NaN given a percentile
# @info "Threshold estimation"
# @time threshold = first(quantile(prediction9[taxa_occ[9]], [0.05]))
# @info "5% threshold:\t$(round(threshold; digits=3))"
# # Filter the predictions based on the threshold
# @info "Final prediction filtering"
# @time for i in eachindex(prediction9.grid)
#     prediction9.grid[i] < threshold && (prediction9.grid[i] = NaN)
# end
#
# @info "Extract and crop bioclim variables"
# @time wc_vars = [clip(worldclim(i), taxa_occ[10]) for i in 1:19];
# # Make the prediction for each layer
# @info "Predictions for each layer"
# @time predictions = [bioclim(wc_vars[i], taxa_occ[10]) for i in 1:length(wc_vars)];
# # Make the final prediction by taking the minimum
# @info "Minimum-consensus aggregation"
# @time prediction10 = reduce(minimum, predictions);
# # Get the threshold for NaN given a percentile
# @info "Threshold estimation"
# @time threshold = first(quantile(prediction10[taxa_occ[10]], [0.05]))
# @info "5% threshold:\t$(round(threshold; digits=3))"
# # Filter the predictions based on the threshold
# @info "Final prediction filtering"
# @time for i in eachindex(prediction10.grid)
#     prediction10.grid[i] < threshold && (prediction10.grid[i] = NaN)
# end
#
# @info "Extract and crop bioclim variables"
# @time wc_vars = [clip(worldclim(i), taxa_occ[11]) for i in 1:19];
# # Make the prediction for each layer
# @info "Predictions for each layer"
# @time predictions = [bioclim(wc_vars[i], taxa_occ[11]) for i in 1:length(wc_vars)];
# # Make the final prediction by taking the minimum
# @info "Minimum-consensus aggregation"
# @time prediction11 = reduce(minimum, predictions);
# # Get the threshold for NaN given a percentile
# @info "Threshold estimation"
# @time threshold = first(quantile(prediction11[taxa_occ[11]], [0.05]))
# @info "5% threshold:\t$(round(threshold; digits=3))"
# # Filter the predictions based on the threshold
# @info "Final prediction filtering"
# @time for i in eachindex(prediction11.grid)
#     prediction11.grid[i] < threshold && (prediction11.grid[i] = NaN)
# end
#
# @info "Extract and crop bioclim variables"
# @time wc_vars = [clip(worldclim(i), taxa_occ[12]) for i in 1:19];
# # Make the prediction for each layer
# @info "Predictions for each layer"
# @time predictions = [bioclim(wc_vars[i], taxa_occ[12]) for i in 1:length(wc_vars)];
# # Make the final prediction by taking the minimum
# @info "Minimum-consensus aggregation"
# @time prediction12 = reduce(minimum, predictions);
# # Get the threshold for NaN given a percentile
# @info "Threshold estimation"
# @time threshold = first(quantile(prediction12[taxa_occ[12]], [0.05]))
# @info "5% threshold:\t$(round(threshold; digits=3))"
# # Filter the predictions based on the threshold
# @info "Final prediction filtering"
# @time for i in eachindex(prediction12.grid)
#     prediction12.grid[i] < threshold && (prediction12.grid[i] = NaN)
# end
#
# @info "Extract and crop bioclim variables"
# @time wc_vars = [clip(worldclim(i), taxa_occ[13]) for i in 1:19];
# # Make the prediction for each layer
# @info "Predictions for each layer"
# @time predictions = [bioclim(wc_vars[i], taxa_occ[13]) for i in 1:length(wc_vars)];
# # Make the final prediction by taking the minimum
# @info "Minimum-consensus aggregation"
# @time prediction13 = reduce(minimum, predictions);
# # Get the threshold for NaN given a percentile
# @info "Threshold estimation"
# @time threshold = first(quantile(prediction13[taxa_occ[13]], [0.05]))
# @info "5% threshold:\t$(round(threshold; digits=3))"
# # Filter the predictions based on the threshold
# @info "Final prediction filtering"
# @time for i in eachindex(prediction13.grid)
#     prediction13.grid[i] < threshold && (prediction13.grid[i] = NaN)
# end
#
# @info "Extract and crop bioclim variables"
# @time wc_vars = [clip(worldclim(i), taxa_occ[14]) for i in 1:19];
# # Make the prediction for each layer
# @info "Predictions for each layer"
# @time predictions = [bioclim(wc_vars[i], taxa_occ[14]) for i in 1:length(wc_vars)];
# # Make the final prediction by taking the minimum
# @info "Minimum-consensus aggregation"
# @time prediction14 = reduce(minimum, predictions);
# # Get the threshold for NaN given a percentile
# @info "Threshold estimation"
# @time threshold = first(quantile(prediction14[taxa_occ[14]], [0.05]))
# @info "5% threshold:\t$(round(threshold; digits=3))"
# # Filter the predictions based on the threshold
# @info "Final prediction filtering"
# @time for i in eachindex(prediction14.grid)
#     prediction14.grid[i] < threshold && (prediction14.grid[i] = NaN)
# end
#
# @info "Extract and crop bioclim variables"
# @time wc_vars = [clip(worldclim(i), taxa_occ[15]) for i in 1:19];
# # Make the prediction for each layer
# @info "Predictions for each layer"
# @time predictions = [bioclim(wc_vars[i], taxa_occ[15]) for i in 1:length(wc_vars)];
# # Make the final prediction by taking the minimum
# @info "Minimum-consensus aggregation"
# @time prediction15 = reduce(minimum, predictions);
# # Get the threshold for NaN given a percentile
# @info "Threshold estimation"
# @time threshold = first(quantile(prediction8[taxa_occ[15]], [0.05]))
# @info "5% threshold:\t$(round(threshold; digits=3))"
# # Filter the predictions based on the threshold
# @info "Final prediction filtering"
# @time for i in eachindex(prediction15.grid)
#     prediction15.grid[i] < threshold && (prediction15.grid[i] = NaN)
# end
#
# @info "Extract and crop bioclim variables"
# @time wc_vars = [clip(worldclim(i), taxa_occ[16]) for i in 1:19];
# # Make the prediction for each layer
# @info "Predictions for each layer"
# @time predictions = [bioclim(wc_vars[i], taxa_occ[16]) for i in 1:length(wc_vars)];
# # Make the final prediction by taking the minimum
# @info "Minimum-consensus aggregation"
# @time prediction16 = reduce(minimum, predictions);
# # Get the threshold for NaN given a percentile
# @info "Threshold estimation"
# @time threshold = first(quantile(prediction16[taxa_occ[16]], [0.05]))
# @info "5% threshold:\t$(round(threshold; digits=3))"
# # Filter the predictions based on the threshold
# @info "Final prediction filtering"
# @time for i in eachindex(prediction16.grid)
#     prediction16.grid[i] < threshold && (prediction16.grid[i] = NaN)
# end

# # predictions = [prediction1, prediction2, prediction3, prediction4]
# # prediction5, prediction6, prediction7, prediction8]
#                # prediction9, prediction10, prediction11,
#                # prediction12, prediction13, prediction14, prediction15, prediction16]
#
# ## Export predictions
# @save "../SDM/data/predictions_vulpes_vulpes_am.jld2" predictions
#
# # Test import
# @load "data/predictions_vulpes_vulpes.jld2" predictions

max_stars_repo_name: kiristern/SDM

max_issues_repo_name: kiristern/SDM

max_forks_repo_name: kiristern/SDM