text:
# tests parsing in src/traitfns.jl

import Base: ==

f1e  = :(f1{X<:Int,TT; D1{X}, D2{X,TT}}(x::X,y::TT))
f1e_b= :(f1{X<:Int,TT; D1{X}, D2{X,TT}}(x::X,y::TT) = ())

f1e_function = :(function f1{X<:Int,TT; D1{X}, D2{X,TT}}(x::X,y::TT)
                ()
                end)

f1e_p = Traits.ParsedFn(
                  :f1,
                  :(f1{X<:Int,TT}),
                  Any[Expr(:<:, :X, :Int),:TT],
                  Any[:(x::X), :(y::TT)],
                  Any[:(D1{X}), :(D2{X,TT})],
                  :(()))
f1e_pt = Traits.ParsedFn(
                 :f1,
                 :(f1{X1<:Int,X2}),
                 Any[Expr(:<:, :X1, :Int),:X2],
                 Any[:(x1::X1), :(x2::X2)],
                 Any[:(D1{X1}), :(D2{X1,X2})],
                 :(()))
function ==(p::Traits.ParsedFn, q::Traits.ParsedFn)
    out = true
    for n in fieldnames(p)
        if n==:body # tricky to compare...
            continue
        end
        out = out && getfield(p,n)==getfield(q,n)
        if !out
            @show n, getfield(p,n), getfield(q,n)
        end
    end
    out
end

@test Traits.parsetraitfn_head(f1e)==f1e_p
@test Traits.translate_head(Traits.parsetraitfn_head(f1e))==f1e_pt
@test Traits.parsetraitfn(f1e_b)==(f1e_p, f1e_pt)

(a,b) = Traits.parsetraitfn(f1e_function)
a.body =:()
b.body =:()
@test a==f1e_p
@test b==f1e_pt

@test Traits.makefnhead(f1e_p.name, f1e_p.typs, f1e_p.sig).args[1]==f1e_p.fun
@test Traits.makefnhead(f1e_p.name, f1e_p.typs, f1e_p.sig).args[2:end]==f1e_p.sig

@test Traits.makefncall(f1e_p.name, f1e_p.sig)==:(f1(x,y))

@test Traits.get_concrete_type_symb(f1e_p.typs)==Any[:Int, :Any]

@test Traits.make_Type_sig([s.args[2] for s in f1e_p.sig])==Any[:(::Type{X}), :(::Type{TT})]

# next case
f1e2= :(f1{X<:Int; D1{X}, D2{X,X}}(x::X,y::X) = ())
f1e2_pt = Traits.ParsedFn(
                 :f1,
                 :(f1{X1<:Int}),
                 Any[Expr(:<:, :X1, :Int)],
                 Any[:(x1::X1), :(x2::X1)],
                 Any[:(D1{X1}), :(D2{X1,X1})],
                 :(()))
@test Traits.parsetraitfn(f1e2)[2]==f1e2_pt

# next case
f1ee = :(tfd{K, V; D2{K,V}}(x::Vector{K}, y::V) = ())
f1ee_p = Traits.ParsedFn(
                  :tfd,
                  :(tfd{K,V}),
                  Any[:K,:V],
                  Any[:(x::Vector{K}), :(y::V)],
                  Any[:(D2{K,V})],
                  :(()))
f1ee_pt = Traits.ParsedFn(
                  :tfd,
                  :(tfd{X1,X2}),
                  Any[:X1,:X2],
                  Any[:(x1::Vector{X1}), :(x2::X2)],
                  Any[:(D2{X1,X2})],
                  :(()))
@test Traits.parsetraitfn(f1ee)==(f1ee_p, f1ee_pt)

# next case
f1ee2 = :(tfd{K, V; D2{K,V}}(x::Dict{K,Int}, y::V) = ())
f1ee2_p = Traits.ParsedFn(
                  :tfd,
                  :(tfd{K,V}),
                  Any[:K,:V],
                  Any[:(x::Dict{K,Int}), :(y::V)],
                  Any[:(D2{K,V})],
                  :(()))
f1ee2_pt = Traits.ParsedFn(
                  :tfd,
                  :(tfd{X1,X2}),
                  Any[:X1,:X2],
                  Any[:(x1::Dict{X1,Int}), :(x2::X2)],
                  Any[:(D2{X1,X2})],
                  :(()))
@test Traits.parsetraitfn(f1ee2)==(f1ee2_p, f1ee2_pt)

# next case
f1d2 = :(tfd{K, V; D2{K,V}}(x::Dict{K,Int}, y::V, z::Int) = ())
f1d2_p = Traits.ParsedFn(
                  :tfd,
                  :(tfd{K,V}),
                  Any[:K,:V],
                  Any[:(x::Dict{K,Int}), :(y::V), :(z::Int)],
                  Any[:(D2{K,V})],
                  :(()))
f1d2_pt = Traits.ParsedFn(
                  :tfd,
                  :(tfd{X1,X2}),
                  Any[:X1,:X2],
                  Any[:(x1::Dict{X1,Int}), :(x2::X2), :(x3::Int)],
                  Any[:(D2{X1,X2})],
                  :(()))
@test Traits.parsetraitfn(f1d2)==(f1d2_p, f1d2_pt)

# trait function creation
eval(:(@traitfn $f1e_b))
eval(:(@traitfn $f1e_function))
eval(:(@traitfn $f1e2))
# these give warnings for now
eval(:(@traitfn $f1ee))
eval(:(@traitfn $f1ee2))
eval(:(@traitfn $f1d2))

@traitfn yt1{X,Y; Arith{X,Y}}(x::X,y::Y) = x+y
a = yt1(5,6)
@test yt1(5,6)==5+6
@traitfn xt1{X<:Int,Y<:AbstractFloat; Arith{X,Y}}(x::X,y::Y) = x-y
@test_throws MethodError xt1(5,6)

od = ObjectIdDict(); od[5]=8
@traitfn yt1{X,Y; IsAssociative{X}}(x::X,y::Y) = x==y
@test yt1(Dict(5=>7), Dict(5=>8))==false
@traitfn xt1{X<:Dict,Y<:ObjectIdDict; IsAssociative{X}}(x::X,y::Y) = x==y
@test xt1(Dict(5=>7), od)==false
@test_throws MethodError xt1(od, od)

# mixing with normal methods:
xt1(x::Int, y::Int) = 77
@test xt1(5,6)==77

# scope of helper methods when defined in different modules
####
## case 1

    module Mod1
    using Traits
    # to check what happens with exports, etc
    export MTyp1, M1Tr1, tf1, barbar

    type MTyp1
        mt::Int
    end

    @traitdef M1Tr1{X} begin
        barbar(X, Int) -> AbstractString
    end
    @traitimpl M1Tr1{MTyp1} begin
        barbar(a::MTyp1, b::Int) = "MTyp1, barbar"^b
    end

    @traitfn tf1{X, Y<:Int;  M1Tr1{X}}(a::X, b::Y) = barbar(a, b)

    # a normal generic function to be used later
    f1(x::Int) = 5x

    end # module Tst

using Mod1
@test tf1(MTyp1(7), 5)=="MTyp1, barbar"^5

# now try to extend tf1

# 1) add a new type to M1Tr1
type MTyp99
    uu::Int
end
@traitimpl M1Tr1{MTyp99} begin
    barbar(a::MTyp99, b::Int) = "MTyp99, barbar"^b
end
@test tf1(MTyp99(9), 5)=="MTyp99, barbar"^5
@test tf1(Main.MTyp99(9), 5)=="MTyp99, barbar"^5

# 2) add a new method to tf1
@traitdef M1Tr100{X} begin
    foofoo(X) -> AbstractString
end
type MTyp00
    gg::Int
end
@traitimpl M1Tr100{MTyp00} begin
    foofoo(a::MTyp00) = "MTyp00 foofoo"
end

import Mod1.tf1
# @show tf1(Traits._TraitStorage, MTyp1, Int)
# @show tf1(Traits._TraitStorage,Any,Int)
# typs = Any[:(Traits._TraitStorage),:X, Expr(:<:, :X1, :Int)]
# methods(eval(:tf1), Traits. get_concrete_type_Typetuple(typs))
println("  These warnings are ok:")  # Well, I'm not sure whether they are ok.  But at least normal...
@traitfn tf1{X, Y<:Int;  M1Tr100{X}}(a::X, b::Y) = foofoo(a)^b
println("  endof ok-warnings.")
#methods(eval(:tf1), Traits. get_concrete_type_Typetuple(typs))
#@show tf1(Traits._TraitStorage, Any, Int)

@test tf1(MTyp1(7), 5)=="MTyp1, barbar"^5
@test tf1(MTyp00(9), 5)== "MTyp00 foofoo"^5

# @show tf1(Traits._TraitStorage, MTyp1, Int)
# @show tf1(Traits._TraitStorage, Any, Int)

# # error:
@test tf1(MTyp99(9), 5)=="MTyp99, barbar"^5

# case 2: adding a @traitfn to a generic function defined in another module
@traitfn tf1_2{X;  M1Tr100{X}}(a::X) = foofoo(a)^99

############
# error
@traitfn ffgg{T; Eq{T,T}}(x::T,y::T) = 2x==y
@test ffgg(5,6)==false
@test ffgg(6,12)==true

######
# traitmethods in modules
@traitfn ff879{T; Eq{T,T}}(x::T,y::T) = 2x==y
@test length(traitmethods(ff879))==1
module A9374
using Traits
using Base.Test
@traitfn ff879{T; Eq{T,T}}(x::T,y::T) = 2x==y
@test length(traitmethods(ff879))==1
end
@test length(traitmethods(ff879))==1


text:
const sample_typelists = [
    (Int,),
    (Int, Float16),
    (Dict{Any,Any},),
    (Dict{Any,Any}, Int, Integer),
    ()
]

@testset "Encode-decode" begin
    @test Catwalk.encode(Int)() isa Catwalk.TypeListItem
    for typelist in sample_typelists
        @test Catwalk.decode(Catwalk.encode(typelist...)) == typelist
    end
end

@testset "findfirst" begin
    for typelist in sample_typelists
        for l = 1:5
            if length(typelist) >= l
                @test findfirst(typelist[l], Catwalk.encode(typelist...)) == l
            end
        end
    end
end

ext: jl


text:
module MNISTDataset
# package code goes here

using IDXParser

export parseLabels, parseImages

cwd = Base.dirname(Base.source_path())

@doc """
    Parse the file containing the labels for the digits.
    name: Either "train" or "test"
    multi: Whether or not to make the output vector a vector or an nx10 matrix with each column representing one class for the data (i.e. 7 becomes [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0]
    useBetterValues: Instead of using 0.0 and 1.0 for the above values, use 0.1 and 0.9
""" ->
function parseLabels(name; multi=true, useBetterValues=false)
    file   = Base.joinpath(cwd, "..", "data", "$(name)-labels-idx1-ubyte")
    input  = parseIDX(file)

    if multi
        values = useBetterValues ? [0.1, 0.9] : [0.0, 1.0]
        labels = ones(Float64, (length(input), 10)) .* values[1]

        for i = 1:length(input)
            num            = input[i] == 0 ? 10 : input[i]
            labels[i, num] = values[2]
        end
    else
        labels = input
    end

    return labels
end

@doc """
     Parse the file containing the image data for the digits.
     name: Either "train" or "test"
     useDecimal: Transform the value range from 0:255 to 0.0:1.0
     addBias: Add a column to the left-end of the output matrix containing a value for the bias
""" ->
function parseImages(name; useDecimal=true, addBias=false)
    file      = Base.joinpath(cwd, "..", "data", "$(name)-images-idx3-ubyte")
    input     = parseIDX(file)

    numImages = size(input, 1)
    numPixels = prod(size(input)[2:end])
    input_res = reshape(input, (numImages, numPixels))

    if addBias
        bias   = ones(UInt8, (numImages, 1)) .* 255
        input_res = [bias input_res]
    end

    if useDecimal
        images = input_res ./ 255.0
    else
        images = input_res
    end

    return images
end

end # module


text:
module LatinHypercubeSampling

export  randomLHC,
        scaleLHC,
        AudzeEglaisObjective,
        AudzeEglaisObjective!,
        LHCoptim,
        LHCoptim!,
        subLHCoptim,
        subLHCindex,
        Categorical,
        Continuous,
        LHCDimension

using StatsBase
using Random
import Random.randperm

abstract type LHCDimension end

struct Categorical <: LHCDimension
    levels::Int64
    weight::Float64
end

struct Continuous <: LHCDimension
end

Categorical(x) = Categorical(x,0.0)

include("GA.jl")
include("AudzeEglaisObjective.jl")

#make an @threads-equivalent that is a no-op if threading is not requested
macro maybe_threaded(flag, ex)
    esc(quote
        if !$flag
            $ex
        else
            Threads.@threads $ex
        end
    end)
end

function randperm(rng,dim::Continuous,n)
    randperm(rng,n)
end

function randperm(rng,dim::Categorical,n)
    lvls = dim.levels
    out = Array{Int64}(undef,n)

    for i = 1:lvls
        s = round(Int,(i-1)*n/lvls) + 1
        e = round(Int,i*n/lvls)
        out[s:e] .= i
    end
    shuffle!(rng,out)
end

"""
    function scaleLHC(LHC,scale_range)

Scale a Latin Hypercube to the sizes specified in scale_range. Where the range for
each dimension is specified in an vector of tuples.

# Examples
```
julia> plan = randomLHC(5,2)
5×2 Array{Int64,2}:
 5  3
 4  4
 2  2
 1  5
 3  1
julia> scaleLHC(plan,[(-1,1),(10,100)])
5×2 Array{Float64,2}:
  1.0   55.0
  0.5   77.5
 -0.5   32.5
 -1.0  100.0
  0.0   10.0
```
"""
function scaleLHC(LHC::Array{T,2},scale_range::Array{Tuple{U,V},1}) where V where U where T

    scaledLHC = Array{Float64,2}(undef,size(LHC))

    for i in 1:size(LHC,2)
        LHCcol = LHC[:,i]
        old_min, old_max = extrema(LHCcol)
        new_min, new_max = scale_range[i]

        @. scaledLHC[:,i] = (((LHCcol - old_min) * (new_max - new_min)) / (old_max - old_min)) + new_min
    end
    return scaledLHC
end

"""
    function randomLHC(n::Int, d::Int)
Generate a random Latin Hypercube with `d` dimensions and `n` sample points.
"""
function randomLHC(n::Int,d::Int)
    randomLHC(Random.GLOBAL_RNG,n,d)
end

function randomLHC(rng,n::Int,d::Int)

    dims = [Continuous() for i in 1:d]

    return randomLHC(rng,n,dims)

end

function randomLHC(n::Int,dims::T) where T<:AbstractArray{S} where S<:LHCDimension
    randomLHC(Random.GLOBAL_RNG,n,dims)
end

function randomLHC(rng,n::Int,dims::T) where T<:AbstractArray{S} where S<:LHCDimension

    LHC = Array{Int}(undef,n,length(dims))
    for (i, dim) = enumerate(dims)
        LHC[:,i] = randperm(rng,dim,n)
    end

    return LHC

end

"""
    function LHCoptim(n::Int,d::Int,gens;   rng::U=Random.GLOBAL_RNG,
                                            popsize::Int=100,
                                            ntour::Int=2,
                                            ptour=0.8,
                                            dims::Array{T,1}=[Continuous() for i in 1:d],
                                            interSampleWeight::Float64=1.0,
                                            periodic_ae::Bool=false,
                                            ae_power::Union{Int,Float64}=2) where T <: LHCDimension where U <: AbstractRNG
Produce an optimized Latin Hyper Cube with `d` dimensions and `n` sample points.
Optimization is run for `gens` generations. Returns a tuple of the sample plan and
the optimization fitness history.
"""
function LHCoptim(n::Int,d::Int,gens;   rng::U=Random.GLOBAL_RNG,
                                        popsize::Int=100,
                                        ntour::Int=2,
                                        ptour=0.8,
                                        dims::Array{T,1}=[Continuous() for i in 1:d],
                                        interSampleWeight::Float64=1.0,
                                        periodic_ae::Bool=false,
                                        ae_power::Union{Int,Float64}=2,
                                        threading=false) where T <: LHCDimension where U <: AbstractRNG

    #populate first individual
    X = randomLHC(rng,n,d)

    LHCoptim!(X,gens,rng=rng,popsize=popsize,ntour=ntour,ptour=ptour,
                dims=dims,interSampleWeight=interSampleWeight,periodic_ae=periodic_ae,ae_power=ae_power,threading=threading)

end

"""
    function LHCoptim!(X::Array{Int,2},gens;    rng::U=Random.GLOBAL_RNG,
                                                popsize::Int=100,
                                                ntour::Int=2,
                                                ptour::Float64=0.8,
                                                dims::Array{T,1}=[Continuous() for i in 1:size(X,2)],
                                                interSampleWeight::Float64=1.0,
                                                periodic_ae::Bool=false,
                                                ae_power::Union{Int,Float64}=2) where T <: LHCDimension where U <: AbstractRNG
Same as LHCoptim(n::Int,d::Int,gens;popsize::Int=100,ntour::Int=2,ptour=0.8) but using an
existing plan. Useful for continued optimization. Returns a tuple of the sample plan and
the optimization fitness history.
"""
function LHCoptim!(X::Array{Int,2},gens;    rng::U=Random.GLOBAL_RNG,
                                            popsize::Int=100,
                                            ntour::Int=2,
                                            ptour::Float64=0.8,
                                            dims::Array{T,1}=[Continuous() for i in 1:size(X,2)],
                                            interSampleWeight::Float64=1.0,
                                            periodic_ae::Bool=false,
                                            ae_power::Union{Int,Float64}=2,
                                            threading=false) where T <: LHCDimension where U <: AbstractRNG

    #preallocate memory
    n, d = size(X)                             #Num points, num dimensions
    mut_inds = Array{Int}(undef,2)             #Storage of indices to swap in mutation
    tour_inds = Array{Int}(undef,ntour)        #Storage of indices for tournament selection

    #allocate first population
    pop = [randomLHC(rng,n,dims) for i = 1:popsize]
    pop[1] = X

    nextpop = deepcopy(pop)
    fitness = Vector{Float64}(undef, popsize)
    fitnessInds = Vector{Int64}(undef, popsize)
    bestfits = Vector{Float64}(undef, gens+1)

    #crossover for even population size
    popEven::Int64 = iseven(popsize)*-1

    #dynamic mutation rate table
    muts = Array{Int}(undef, gens)
    for l = 1:gens
        muts[l] = round(Int,-n/(gens*0.75)*l+n)
        if muts[l] < 1
            muts[l] = 1
        end
    end

    #evaluate first populations fitness
    @maybe_threaded threading for i = 1:popsize
        fitness[i] = AudzeEglaisObjective(pop[i];
                                            interSampleWeight=interSampleWeight,
                                            dims=dims,
                                            periodic_ae=periodic_ae,
                                            ae_power=ae_power)
    end

    #save the best individual and its fitness
    bestfit, bestind = findmax(fitness)
    nextpop[1] = copy(pop[bestind])
    bestfits[1] = bestfit

    #ensure fixed crossover is only applied to the continuous dimensions
    continuousDims = findall(dims.==Ref(Continuous()))

    #allocate thread-local storage for crossover
    numthreads = Threads.nthreads()
    offsprone = [similar(pop[1][:,1]) for i in 1:numthreads]
    offsprtwo = [similar(pop[1][:,1]) for i in 1:numthreads]

    #ReentrantLock is not threadsafe on Julia < v1.2, so then we need to use a different lock
    #this version check could be removed if support for v1.0 is dropped
    if VERSION < v"1.2"
        randlock = Threads.SpinLock()
    else
        randlock = ReentrantLock()
    end

    #iterate for gens generations
    for k = 1:gens

        #tournament selection
        sortperm!(fitnessInds,fitness)
        for i = 2:popsize
            winnerInd = tournament!(rng,fitnessInds,ntour,tour_inds,ptour)
            nextpop[i] = copy(pop[winnerInd])
        end

        #create children from crossover
        @maybe_threaded threading for i = 2:2:popsize+popEven
            tid = Threads.threadid()
            for j in continuousDims
                r = lock(randlock) do
                    rand(rng)
                end
                if r < 1.0/length(continuousDims)
                    parone = nextpop[i]
                    partwo = nextpop[i+1]

                    fixedcross!(rng,offsprone[tid], offsprtwo[tid], view(parone,:,j), view(partwo,:,j),randlock)
                    nextpop[i][:,j], nextpop[i+1][:,j] = offsprone[tid], offsprtwo[tid]
                end
            end
        end

        #perform inversion and mutation
        for i = 2:popsize
            for j = 1:d
                if rand(rng) < 1.0/d     #probability for inversion
                    nextindividual = nextpop[i]
                    inversion!(rng,view(nextindividual,:,j))
                end
            end

            for j=1:muts[gens]
                mutateLHC!(rng,nextpop[i],mut_inds)
            end
        end

        #evaluate fitness
        @maybe_threaded threading for i = 1:popsize
            fitness[i] = AudzeEglaisObjective(nextpop[i];
                                                interSampleWeight=interSampleWeight,
                                                dims=dims,
                                                periodic_ae=periodic_ae,
                                                ae_power=ae_power)
        end

        #set the first individual to the best and save the fitness
        bestfit, bestind = findmax(fitness)
        nextpop[1] = nextpop[bestind]
        bestfits[k+1] = bestfit
        pop = deepcopy(nextpop)
    end

    return nextpop[1], bestfits

end

"""
    function subLHCoptim(X,n::Int,gens; rng::U=Random.GLOBAL_RNG,
                                        popsize::Int=100,
                                        ntour::Int=2,
                                        ptour::Float64=0.8,
                                        periodic_ae::Bool=false,
                                        ae_power::Union{Int,Float64}=2) where U <: AbstractRNG
Produce an optimized Latin Hyper Cube with `n` sample points from a subset of points in
`X`. Optimization is run for `gens` generations. Returns a tuple of the sample plan and
the optimization fitness history.
"""
function subLHCoptim(X,n::Int,gens; rng::U=Random.GLOBAL_RNG,
                                    popsize::Int=100,
                                    ntour::Int=2,
                                    ptour::Float64=0.8,
                                    periodic_ae::Bool=false,
                                    ae_power::Union{Int,Float64}=2) where U <: AbstractRNG

    #preallocate memory
    nLarge, d = size(X)
    pop = [Array{Int}(undef,n,d) for i = 1:popsize+1]
    nextpop = similar(pop)
    fitness = Vector{Float64}(undef, popsize+1)
    fitnessInds = Vector{Int64}(undef, popsize+1)
    bestfits = Array{Float64}(undef, gens)

    tour_inds = Array{Int}(undef,ntour)     #Storage of indices for tournament selection
    tour_fitinds = Array{Int}(undef,ntour)  #Storage of fitness for tournament selection

    #dynamic mutation rate table
    muts = Array{Int}(undef, gens)
    for l = 1:gens
        muts[l] = round(Int,-n/(gens*0.75)*l+n)
        if muts[l] < 1
            muts[l] = 1
        end
    end

    #populate first population with random subLHC's and evaluate fitness
    for i = 1:popsize+1
        subInds = sample(rng, 1:nLarge, n, replace = false)
        pop[i] = X[subInds,:]
        fitness[i] = AudzeEglaisObjective(pop[i];
                                            periodic_ae=periodic_ae,
                                            ae_power=ae_power,
                                            periodic_n=nLarge)
    end

    #save the best individual and its fitness
    bestfit, bestind = findmax(fitness)
    nextpop[1] = copy(pop[bestind])
    bestfits[1] = bestfit

    #iterate for gens generations
    for k = 1:gens

        #tournament selection
        sortperm!(fitnessInds,fitness)
        for i = 2:popsize+1
            winnerInd = tournament!(rng,fitnessInds,ntour,tour_inds,ptour)
            nextpop[i] = copy(pop[winnerInd])
        end

        #perform mutation
        for i = 2:popsize+1
            for j=1:muts[gens]
                subPoint = sample(rng,1:n) #Randomly chosen point subLHC
                largePoint = sample(rng,1:nLarge) #Randomly chosen point largeLHC

                while nextpop[i][subPoint,:] == X[largePoint,:]
                    largePoint = sample(rng,1:nLarge) #Choose new random point in largeLHC
                end
                nextpop[i][subPoint,:] = X[largePoint,:]
            end
        end

        #evaluate fitness
        for i = 1:popsize+1
            fitness[i] = AudzeEglaisObjective(nextpop[i];
                                                periodic_ae=periodic_ae,
                                                ae_power=ae_power,
                                                periodic_n=nLarge)
        end

        #set the first individual to the best and save the fitness
        bestfit, bestind = findmax(fitness)
        nextpop[1] = nextpop[bestind]
        bestfits[k] = bestfit
        pop = deepcopy(nextpop)
    end

    return pop[1], bestfits

end

"""
    function subLHCindex(X,Xsub)
Index in the large LHC to get the subLHC.
"""
function subLHCindex(X,Xsub)

    nsub = size(Xsub,1)
    subInds = Array{Int}(undef, nsub)

    for i = 1:nsub
        A = all(Xsub[i,:]' .== X, dims=2)
        subInds[i] = (LinearIndices(A))[findall(A)][1]
    end
    return subInds
end

end # module

max_stars_repo_name: MrUrq/LatinHypercubeSampling

max_issues_repo_name: MrUrq/LatinHypercubeSampling

max_forks_repo_name: MrUrq/LatinHypercubeSampling


text:
###############################################################################
#
# Save the dataset
#
###############################################################################
"""
    save_LUT!(ds::GriddedDataset{FT},
              filename::String) where {FT<:AbstractFloat}
    save_LUT!(data::Array{FT,2},
              filename::String) where {FT<:AbstractFloat}
    save_LUT!(data::Array{FT,3},
              filename::String,
              var_name::String,
              var_attr::Dict{String,String}) where {FT<:AbstractFloat}

Save the dataset, given
- `ds` [`GriddedDataset`](@ref) type struct
- `filename` NC file name to save to
- `data` Temporaty data
- `var_name` Variable name in nc file
- `var_attr` Variable attributes in nc file

Note that `save_LUT!(data, filename)` is designed to use with temporaty data
    only, be cautious when using this function.
"""
function save_LUT!(
            ds::GriddedDataset{FT},
            filename::String
) where {FT<:AbstractFloat}
    save_LUT!(ds.data, filename, ds.var_name, ds.var_attr);

    return nothing
end

function save_LUT!(
            data::Array{FT,2},
            filename::String
) where {FT<:AbstractFloat}
    # attributes
    varatts = Dict("longname" => "Variable" , "units" => "-");

    # create nc file
    save_nc!(filename, "Var", varatts, data);

    return nothing
end

function save_LUT!(
            data::Array{FT,3},
            filename::String,
            var_name::String,
            var_attr::Dict{String,String}
) where {FT<:AbstractFloat}
    save_nc!(filename, var_name, var_attr, data);

    return nothing
end


text:
#######################################
#
# Global field power
#
#######################################

function gfp(x::Array)

    samples, sensors = size(x)

    @info("Computing global field power for $sensors sensors and $samples samples")

    result = zeros(samples, 1)

    for sample = 1:samples
        u = vec(x[sample, :]) .- Statistics.mean(x[sample, :])
        sumsqdif = 0
        for sensor = 1:sensors
            for sensor2 = 1:sensors
                sumsqdif += (u[sensor] - u[sensor2])^2
            end
        end
        result[sample] = sqrt.(sumsqdif / (2 * length(samples)))
    end

    return result
end


text:
using JackPolynomials
using Test

@testset "Jack" begin
  @testset "Jack at x = 0" begin
    @test isapprox(Jack([0.0; 0.0], Int64[], 2.0), 1.0)
    @test isapprox(Jack([0.0; 0.0; 0.0], [3; 2], 2.0), 0.0)
  end
  @testset "Jack = 0 when l(lambda)>l(x)" begin
    @test isapprox(Jack([1.0; 2.0], [3; 2; 1], 2.0), 0.0)
    @test isapprox(Jack([1.0; 2.0], [3; 2; 1], 1.0), 0.0)
  end
  @testset "Jack = 1 when the partition is empty" begin
    @test isapprox(Jack([1.0; 2.0], Int64[], 2.0), 1.0)
    @test isapprox(Jack([1.0; 2.0], [0; 0], 2.0), 1.0)
  end
  @testset "Jack - some values" begin
    @test isapprox(Jack([2.0; 3.0; 4.0; 5.0], [3; 1], 2.5), 111217.0)
    @test isapprox(Jack([2.0; 3.0; 4.0; 5.0], [4], 2.5), 500338.5)
  end
end

@testset "Jack symbolic" begin
  @testset "Jack at x = 0" begin
    @test isapprox(JackPolynomial(2, Int64[], 2.0)(0.0, 0.0), 1.0)
    @test isapprox(JackPolynomial(3, [3; 2], 2.0)(0.0, 0.0, 0.0), 0.0)
  end
  @testset "Jack = 0 when l(lambda)>l(x)" begin
    @test JackPolynomial(2, [3; 2; 1], 2.0) == 0.0
    @test JackPolynomial(2, [3; 2; 1], 1.0) == 0.0
  end
  @testset "Jack = 1 when the partition is empty" begin
    @test JackPolynomial(2, Int64[], 2.0) == 1.0
    @test JackPolynomial(2, [0; 0], 2.0) == 1.0
  end
  @testset "Jack - some values" begin
    @test isapprox(JackPolynomial(4, [3; 1], 2.5)(2.0, 3.0, 4.0, 5.0), 111217.0)
    @test isapprox(JackPolynomial(4, [4], 2.5)(2.0, 3.0, 4.0, 5.0), 500338.5)
  end
end

@testset "Schur" begin
  @testset "Schur expansion of (sum x_i)^n" begin
    x = [3.0; 4.0; 5.0; 6.0]
    s = Schur(x, [4]) + 3*Schur(x, [3,1]) + 2*Schur(x, [2,2]) +
      3*Schur(x, [2,1,1]) + Schur(x, [1,1,1,1])
    @test isapprox(s, sum(x)^4)
  end
  @testset "Schur = 0 when l(lambda)>l(x)" begin
    @test isapprox(Schur([1.0; 2.0], [3; 2; 1]), 0.0)
  end
  @testset "Schur (3,2)" begin
    x = [0.3, 2.0, 5.0]
    expected = x[1]^3*x[2]^2 + x[1]^3*x[3]^2 + x[1]^3*x[2]*x[3] +
      x[1]^2*x[2]^3 + x[1]^2*x[3]^3 + 2*x[1]^2*x[2]*x[3]^2 +
      2*x[1]^2*x[2]^2*x[3] + x[1]*x[2]*x[3]^3 + 2*x[1]*x[2]^2*x[3]^2 +
      x[1]*x[2]^3*x[3] + x[2]^2*x[3]^3 + x[2]^3*x[3]^2
    @test isapprox(Schur(x, [3; 2]), expected)
  end
end

@testset "Zonal" begin
  @testset "Zonal sum to the trace" begin
    x = [3.0; 4.0; 5.0]
    z = Zonal(x, [3]) + Zonal(x, [2,1]) + Zonal(x, [1,1,1])
    @test isapprox(z, sum(x)^3)
  end
  @testset "Zonal = 0 when l(lambda)>l(x)" begin
    @test isapprox(Zonal([1.0; 2.0], [3; 2; 1]), 0.0)
  end
end

@testset "ZonalQ" begin
  @testset "ZonalQ sum to the trace" begin
    x = [3.0; 4.0; 5.0]
    z = ZonalQ(x, [3]) + ZonalQ(x, [2,1]) + ZonalQ(x, [1,1,1])
    @test isapprox(z, sum(x)^3)
  end
end


text:
"""

## Envelope

Holds an array of `Line`s, the upper envelope of those lines, and a vector of `Point`s marking the intersections between lines.

### Fields

* `L      `: Vector of [`Line`]@ref
* `env    `: The upper envelope (i.e pointwise maximum) over `L`, itself a [`Line`]@ref
* `env_set`: `true` if envelope vector was set.
* `isects `: Vector of intersections between `Line`s in `L`
* `removed`: Vector of Points removed from `env` during assembly

"""
mutable struct Envelope{T<:Number}
    L      :: Vector{Line{T}}
    env    :: Line{T}
    env_set :: Bool
    isects :: Vector{Point{T}}
    removed :: Vector{Vector{Point{T}}}
    vbound :: T
    function Envelope(x::T) where {T<:Number}
        this = new{T}()
        this.L = Line{T}[]
        this.env = Line([typemin(T)],[typemin(T)])
        this.env_set = false
        this.isects = Point{T}[]
        this.removed = Vector{Point{T}}[Point{T}[] ]
        this.vbound = zero(T)
        return this
    end
    function Envelope(e::Line{T}) where {T<:Number}
        this = new{T}()
        this.L = Line{T}[]
        this.env = deepcopy(e)
        this.env_set = true
        this.isects = Point{T}[]
        this.removed = Vector{Point{T}}[Point{T}[] ]
        this.vbound = zero(T)
        return this
    end
    function Envelope(l::Vector{Line{T}}) where {T<:Number}
        this = new{T}()
        this.L = deepcopy(l)
        this.env = Line([typemin(T)],[typemin(T)])
        this.env_set = false
        this.isects = Point{T}[]
        this.removed = Vector{Point{T}}[Point{T}[] ]
        this.vbound = zero(T)
        return this
    end
end
function show(io::IO,en::Envelope{T}) where {T<:Number}
    print(io,"$T Envelope\n")
    print(io,"env Line set?: $(en.env_set) \n")
    print(io,"num of `Line`s: $(length(en.L))\n")
    print(io,"num of intersections: $(length(en.isects))\n")
    print(io,"num of pts removed: $(sum(length(i) for i in en.removed))\n")
end
size(e::Envelope) = size(e.L)
eltype(e::Envelope) = eltype(e.L)
bound(e::Envelope) = e.vbound
getx(en::Envelope) = en.env.x
gety(en::Envelope) = en.env.y
gets(en::Envelope) = en.isects
getr(en::Envelope) = en.removed
# getLine(en::Envelope,j::Int) = en.L[j]
# function set!(en::Envelope{T},L::Line{T}) where {T<:Number}
#     en.env = L
# end
# function set!(en::Envelope{T},id::Int,l::Line{T}) where {T<:Number}
#     en.L[id] = l
# end
# Base.setindex!(en::Envelope,l::Line,id::Int) = en.L[id] = l
Base.getindex(en::Envelope,id::Int) = en.L[id]

"""
    splitLine(m::Line)

splits a `Line` object at wrong EGM solution points. Wrong solutions appear in kinked regions.
"""
function splitLine(o::Line{T}) where T<:Number

    # 1) find all jump-backs in x-grid
    # println(o.x)
    ii = o.x[2:end].>o.x[1:end-1]
    info("splitLine: ii = $(find(.!(ii)))")
    info("splitLine: x = $(o.x[find(.!(ii))])")

    # 2) if no backjumps at all, exit
    if all(ii)
        # return as an Envelope
        return Envelope(o)
    else
    # 3) else, identify subsets
        i = 1
        sections = Line{T}[]  # an array of Lines
        new_sections = Line{T}[]  # an array of Lines
        while true
            # println(ii)
            j = findfirst(ii .!= ii[1])  # identifies all indices within kinked region from left to right until the first kink
            # println(j)

            # if no more kinks
            if j==0
                if i > 1
                    # add remaining Line
                    push!(sections,o)
                end
                # then break
                break
            end
            newm,o = splitat(o,j)  # split old Line at j
            push!(sections,newm)
            ii = ii[j:end] # chop off from jump index
            i += 1
        end

        # all the ones with 2 un-sorted x corrdinates are illegal lines from connection two proper ones
        # discard those
        ns = [!issorted(s.x) && length(s.x)==2 for s in sections]
        info("which illegal $ns")

        # 4) get rid of illegal sections on x
        for s in eachindex(sections)
            if !ns[s]
                push!(new_sections,sections[s])
                @assert(issorted(sections[s].x))
            end
        end
        return Envelope(new_sections)
    end
end

function upper_env!(e::Envelope{T}) where T<:Number
    # 5) compute upper envelope of all lines
        # - get all x's from all s and sort into a vector xx
        # - interpolate(extrapolate) all s on xx
        # - how to deal with points at which some Line is infeasible?

    if length(e.L)<2
        error("an upper envelope requires by definition at least 2 lines.")
    end

    # - get all x's from all Lines and sort into a vector xx
    xx = sort(unique(vcat([l.x for l in e.L]...)))
    n = length(xx)

    # - interpolate(extrapolate) all Ls on xx
    # this returns a matrix (length(L),n)
    # i.e. each row is the interpolation
    yy = interp(e.L,xx)

    # find the top line at each point in xx
    val,lin_ind = findmax(yy,1)  # colwise max
    subs = map(x->ind2sub(yy,x),lin_ind)  # get subsript indices of colwise maxima
    r_idx = [i[1] for i in subs] # get row indices only: the row index tells us which Line was optimal at that point.

    # Identify changes in optimal Line
    # switch in top line after index s (indexing global support xx)
    # s tells us after which position in xx we have a change in optimal line
    s = find(r_idx[2:end].!=r_idx[1:end-1])

    info("upper_env: jumps after $s")
    info("upper_env: jumps at $(xx[s])")

    # Assemble Upper Envelope from Line segments
    # ==========================================

    if length(s)==0
        # there is one complete upper envelope already
        # return
        e.env = Line(xx,yy[r_idx[1],:])
        # e.isects = [Point(NaN,NaN)]
        # e.removed = [[Point(NaN,NaN)]]
    else
        # sort out which line is top at which index of xx and compute intersection points in between switches
        # s = 1: there is a switch in top line after the first index
        # s = i: there is a switch in top line after the i-th index

        # return Line
        # The envelope starts with the first line that is on top
        # that line is on top until index s[1] in xx, after which the
        # top line changes.
        env = Line(xx[1:s[1]],yy[r_idx[s[1]],1:s[1]] )
        isec = [Point(NaN,NaN) for i in 1:length(s)]

        for id_s in eachindex(s)   # id_s indexes line segment in resulting envelope: for n intersecting lines, there are n-1 segments

            js = s[id_s]  # value of index: position in xx
            jx  = subs[js][2]   # colindex of element in yy before switch takes place
            jjx = subs[js+1][2] # colindex of element in yy after switch took place

            # switching from Line to Line
            from = r_idx[js]
            to   = r_idx[js+1]
            info("from L number $(r_idx[js])")
            info("to L Number $(r_idx[js+1])")

            # xx coordinates between which the switching occurs
            # remember xx is a vector as long as size(yy,2)
            x_from = xx[jx ]  # only pick col coordinate
            x_to   = xx[jjx]  # only pick col coordinate
            info("x_from = $(xx[jx ])")  # only pick col coordinate
            info("x_to   = $(xx[jjx])")  # only pick col coordinate

            # values of both lines at those corrdinates
            v_from = yy[subs[js]...]
            v_to   = yy[subs[js+1]...]
            info("v_from = $(yy[subs[js]...])")
            info("v_to   = $(yy[subs[js+1]...])")

            # The intersction of both lines ∈ [x_from,x_to]
            # complication: intersection could also be on the egdes of this interval.
            # this is easy to check:
            if isapprox(yy[from,jx] , yy[to,jx])
                # intersection is on lower bound of interval
                isec[id_s] = Point(x_from,v_from,i=id_s)
                # don't add intersection to envelope!
            elseif isapprox(yy[from,jjx] , yy[to,jjx])
                # intersection is on upper bound of interval
                isec[id_s] = Point(x_to,v_to,i=id_s)
                # don't add intersection to envelope!
            else
                # need to to compute intersection
                f_closure(z) = interp(e.L[to],[z])[1] - interp(e.L[from],[z])[1]
                x_x = fzero(f_closure, x_from, x_to)
                v_x = interp(e.L[from],[x_x])[1]

                # record intersection as new point
                isec[id_s] = Point(x_x,v_x,i=id_s,newpoint=true)
                # and add to envelope
                append!(env,x_x,v_x)
                info("added $(Point(x_x,v_x,i=id_s)) to envelope")
                info("x in env? $(in(x_x,env.x))")
            end

            # add next line segment to envelope
            # index range s[id_s]+1:s[id_s+1] is
            #     from current switch (next index after): s[id_s]+1
            #     to last index before next switch: s[id_s+1]
            last_ind = id_s==length(s) ? n : s[id_s+1]
            append!(env,xx[js+1:last_ind],yy[to,js+1:last_ind])
        end  # eachindex(s)

        e.env = env
        e.isects = isec
        info("x = $(getx(e))")
        for l in e.L
            # collect points that were removed from Lines
            info("l.x = $(l.x)")
            info("setdiff(l.x,x) = $(setdiff(getx(e),l.x))")
            ix = map(x->!in(x,getx(e)),l.x)
            iy = map(x->!in(x,gety(e)),l.y)
            jj = find(ix .| iy)
            if length(jj) > 0
                push!(e.removed,[Point(l.x[jx],l.y[jx],i=jx) for jx in jj])
            end
        end
        # say that you have now set an upper envelope on this object
        e.env_set = true
        return nothing
    end
end

max_stars_repo_name: floswald/MultiLine

max_issues_repo_name: floswald/MultiLine

max_forks_repo_name: floswald/MultiLine


text:
## data layout

export DataLayout, dispose,
       byteorder, pointersize, intptr,
       sizeof, storage_size, abi_size,
       abi_alignment, frame_alignment, preferred_alignment,
       element_at, offsetof

# forward definition of DataLayout in src/module.jl

Base.unsafe_convert(::Type{API.LLVMTargetDataRef}, dl::DataLayout) = dl.ref

DataLayout(rep::String) = DataLayout(API.LLVMCreateTargetData(rep))

DataLayout(tm::TargetMachine) = DataLayout(API.LLVMCreateTargetDataLayout(tm))

function DataLayout(f::Core.Function, args...; kwargs...)
    data = DataLayout(args...; kwargs...)
    try
        f(data)
    finally
        dispose(data)
    end
end

dispose(data::DataLayout) = API.LLVMDisposeTargetData(data)

Base.string(data::DataLayout) =
    unsafe_message(API.LLVMCopyStringRepOfTargetData(data))

function Base.show(io::IO, data::DataLayout)
    @printf(io, "DataLayout(%s)", string(data))
end

byteorder(data::DataLayout) = API.LLVMByteOrder(data)

pointersize(data::DataLayout) = API.LLVMPointerSize(data)
pointersize(data::DataLayout, addrspace::Integer) =
    API.LLVMPointerSizeForAS(data, addrspace)

intptr(data::DataLayout; ctx::Context) =
    IntegerType(API.LLVMIntPtrTypeInContext(ctx, data))
intptr(data::DataLayout, addrspace::Integer; ctx::Context) =
    IntegerType(API.LLVMIntPtrTypeForASInContext(ctx, data, addrspace))

Base.sizeof(data::DataLayout, typ::LLVMType) =
    Int(API.LLVMSizeOfTypeInBits(data, typ) / 8)
storage_size(data::DataLayout, typ::LLVMType) = API.LLVMStoreSizeOfType(data, typ)
abi_size(data::DataLayout, typ::LLVMType) = API.LLVMABISizeOfType(data, typ)

abi_alignment(data::DataLayout, typ::LLVMType) =
    API.LLVMABIAlignmentOfType(data, typ)
frame_alignment(data::DataLayout, typ::LLVMType) =
    API.LLVMCallFrameAlignmentOfType(data, typ)
preferred_alignment(data::DataLayout, typ::LLVMType) =
    API.LLVMPreferredAlignmentOfType(data, typ)
preferred_alignment(data::DataLayout, var::GlobalVariable) =
    API.LLVMPreferredAlignmentOfGlobal(data, var)

element_at(data::DataLayout, typ::StructType, offset::Integer) =
    API.LLVMElementAtOffset(data, typ, Culonglong(offset))

offsetof(data::DataLayout, typ::StructType, element::Integer) =
    API.LLVMOffsetOfElement(data, typ, element)


text:
"""
    AbstractFiltration{I,T}

A filtration is used to find the edges in filtration and to create simplices. An
`AbstractFiltration{I,T}`'s simplex type is expected to return simplices of type
`<:AbstractCell{_,T,I}`.

An `AbstractFiltration` constructor must accept the `verbose` keyword argument.

# Interface

* [`nv(::AbstractFiltration)`](@ref)
* [`births(::AbstractFiltration)`](@ref)
* [`vertices(::AbstractFiltration)`](@ref)
* [`edges(::AbstractFiltration)`](@ref)
* [`simplex_type(::Type{AbstractFiltration}, dim)`](@ref)
* [`simplex(::AbstractFiltration, ::Val{dim}, vertices)`](@ref)
* [`unsafe_simplex(::AbstractFiltration, ::Val{dim}, vertices)`](@ref)
* [`unsafe_cofacet`](@ref)`(::AbstractFiltration, simplex, vertices, vertex[, edges])`
* [`threshold(::AbstractFiltration)`](@ref)
* [`columns_to_reduce(::AbstractFiltration, ::Any)`](@ref)
* [`emergent_pairs(::AbstractFiltration)`](@ref)
* [`postprocess_diagram(::AbstractFiltration, ::Any)`](@ref)
"""
abstract type AbstractFiltration{I,T} end

function Base.show(io::IO, flt::AbstractFiltration{I,T}) where {I,T}
    return print(io, nameof(typeof(flt)), "{$I, $T}(nv=$(nv(flt)))")
end

"""
    simplex_type(::Type{<:AbstractFiltration}, D)
    simplex_type(::AbstractFiltration, D)

Return the `D`-dimensional simplex type in the filtration. Only the method for the type
needs to be overloaded.

# Examples

```jldoctest
julia> Ripserer.simplex_type(Rips{Int,Float64}, 1)
Simplex{1, Float64, Int64}

julia> Ripserer.simplex_type(Cubical{2,Float16}, 2)
Cube{2, Float16, 2}

```
"""
simplex_type(flt::AbstractFiltration, dim) = simplex_type(typeof(flt), dim)

vertex_type(flt::AbstractFiltration) = simplex_type(typeof(flt), 0)
edge_type(flt::AbstractFiltration) = simplex_type(typeof(flt), 1)

"""
    nv(::AbstractFiltration)

Return the number of vertices in `filtration`.

# Example

```jldoctest
julia> Ripserer.nv(Rips([1 1; 1 1]))
2

```
"""
nv(::AbstractFiltration)

"""
    edges(::AbstractFiltration)

Get edges (1-simplices) in `filtration`. Edges should be of type
[`simplex_type`](@ref)`(filtration, 1)`.

# Example

```
julia> Ripserer.edges(Rips([0 2 1; 2 0 1; 1 1 0], threshold=2))
3-element Array{Simplex{1,Int64,Int64},1}:
 +Simplex{1}([2, 1], 2)
 +Simplex{1}([3, 1], 1)
 +Simplex{1}([3, 2], 1)

```
"""
edges(::AbstractFiltration)

"""
     simplex(::AbstractFiltration, ::Val{D}, vertices)

Return `D`-simplex constructed from `vertices`. Return `nothing`
if simplex is not in filtration. This function is safe to call with vertices that are out of
order. Default implementation sorts `vertices` and calls [`unsafe_simplex`](@ref).

# Example

```
julia> simplex(Rips([0 2 1; 2 0 1; 1 1 0], threshold=2), Val(1), (1, 2))
1-dimensional Simplex(index=1, birth=2):
  +[2, 1]

```
"""
function simplex(flt::AbstractFiltration, ::Val{D}, vertices) where {D}
    if D == 0 && length(vertices) == 1 && vertices[1] > 0
        # No need to run allunique in this case.
        return unsafe_simplex(flt, Val(0), vertices)
    else
        vs = TupleTools.sort(Tuple(vertices); rev=true)
        if allunique(vs) &&
            all(x -> x > 0, vs) &&
            length(vs) == length(simplex_type(flt, D))
            return unsafe_simplex(flt, Val(D), vs)
        end
    end
    return throw(ArgumentError("invalid vertices $(vertices)"))
end

"""
    unsafe_simplex(::AbstractFiltration, ::Val{D}, vertices)

Return `D`-simplex constructed from `vertices`. Return `nothing`
if simplex is not in filtration. The unsafe in the name implies that it's up to the caller
to ensure vertices are sorted and unique.
"""
function unsafe_simplex(flt, ::Val{D}, vertices) where {D}
    return unsafe_simplex(simplex_type(typeof(flt), D), flt, vertices)
end

"""
    unsafe_cofacet(filtration, simplex, cofacet_vertices, v[, edges])
    unsafe_cofacet(::Type{S}, filtration, simplex, cofacet_vertices, v[, edges])

Return cofacet of `simplex` with vertices equal to `cofacet_vertices`. `v` is the
vertex that was added to construct the cofacet. In the case of sparse rips filtrations, an
additional argument `edges` is used. `edges` is a vector that contains the weights on edges
connecting the new vertex to old vertices. `S` is the simplex type which can be used for
dispatch.

The unsafe in the name implies that it's up to the caller to ensure vertices are sorted and
unique.

Default implementation uses [`unsafe_simplex`](@ref).
"""
@inline @propagate_inbounds function unsafe_cofacet(flt, sx, args...)
    return unsafe_cofacet(simplex_type(typeof(flt), dim(sx) + 1), flt, sx, args...)
end
function unsafe_cofacet(::Type{S}, flt, _, vertices, _, args...) where {S}
    return unsafe_simplex(flt, Val(dim(S)), vertices)::Union{S,Nothing}
end

"""
    vertices(::AbstractFiltration)

Return the vertices in filtration. Defaults to `1:n`. The `shape` and `eltype` of the result
can be anything as long as `result[result[i]] == result[i]` holds.

# Example

```jldoctest
julia> vertices(Rips([0 1 1; 1 0 1; 1 1 0]))
Base.OneTo(3)

julia> vertices(Cubical([0 1 1; 1 0 1; 1 1 0]))
3×3 CartesianIndices{2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}:
 CartesianIndex(1, 1)  CartesianIndex(1, 2)  CartesianIndex(1, 3)
 CartesianIndex(2, 1)  CartesianIndex(2, 2)  CartesianIndex(2, 3)
 CartesianIndex(3, 1)  CartesianIndex(3, 2)  CartesianIndex(3, 3)

```
"""
vertices(flt::AbstractFiltration) = Base.OneTo(nv(flt))

"""
    births(::AbstractFiltration)

Get the birth times of vertices in filtration. Defaults to all births being 0. Must return
array of the same shape as the filtration's `vertices`.

# Examples

```jldoctest
julia> flt = Rips([1 1 2; 1 0 1; 2 1 0]);

julia> Ripserer.births(flt)
3-element view(::Vector{Int64}, 1:4:9) with eltype Int64:
 1
 0
 0

```
"""
births(flt::AbstractFiltration{<:Any,T}) where {T} = zeros(T, size(vertices(flt)))

"""
    threshold(::AbstractFiltration)

Get the threshold of filtration. This is the maximum diameter a simplex in the filtration
can have. Defaults to `Inf`.

# Examples

```jldoctest
julia> threshold(Rips([0 2 1; 2 0 1; 1 1 0]))
1

julia> threshold(Cubical([1 1 2; 3 2 1; 0 0 0]))
3

```
"""
threshold(::AbstractFiltration) = Inf

"""
    adjacency_matrix(::AbstractFiltration)

Return the adjacency matrix. For sparse filtrations, this should return a `SparseMatrixCSC`.

# Examples

```jldoctest
julia> Ripserer.adjacency_matrix(Rips([0 2 1; 2 0 1; 1 1 0]))
3×3 Matrix{Int64}:
 0  2  1
 2  0  1
 1  1  0

julia> Ripserer.adjacency_matrix(Rips([0 10 2; 10 0 1; 2 1 0]; sparse=true))
3×3 SparseArrays.SparseMatrixCSC{Int64, Int64} with 4 stored entries:
 ⋅  ⋅  2
 ⋅  ⋅  1
 2  1  ⋅

julia> Ripserer.adjacency_matrix(Custom([(2, 1) => 1, (5, 1) => 2, (3, 4) => 3]))
5×5 SparseArrays.SparseMatrixCSC{Bool, Int64} with 6 stored entries:
 ⋅  1  ⋅  ⋅  1
 1  ⋅  ⋅  ⋅  ⋅
 ⋅  ⋅  ⋅  1  ⋅
 ⋅  ⋅  1  ⋅  ⋅
 1  ⋅  ⋅  ⋅  ⋅

```
"""
adjacency_matrix(::AbstractFiltration)

"""
    columns_to_reduce(::AbstractFilration, prev_column_itr)

List all columns to reduce in next dimension, possibly computing it from previous
columns. Default implementation uses [`coboundary`](@ref) with the all cofacets parameter
set to `Val(false)`.

# Example

```jldoctest
julia> flt = Rips([0 1 1; 1 0 1; 1 1 0]);

julia> Ripserer.columns_to_reduce(flt, Ripserer.edges(flt)) |> collect
1-element Vector{Simplex{2, Int64, Int64}}:
 +Simplex{2}((3, 2, 1), 1)

```
"""
function columns_to_reduce(flt::AbstractFiltration, prev_column_itr)
    return Iterators.flatten(imap(σ -> coboundary(flt, σ, Val(false)), prev_column_itr))
end

"""
    emergent_pairs(::AbstractFiltration)

Return `true` if the emergent pairs optimization is to be performed. Default to returning
`true`. Should be set to `false` for a filtration type that is unable to produce
(co)boundary simplices in the correct order.
"""
emergent_pairs(::AbstractFiltration) = true

"""
    postprocess_diagram(::AbstractFiltration, diagram)

This function is called on each resulting persistence diagram after all intervals have been
computed. Defaults to not doing anything.
"""
postprocess_diagram(::AbstractFiltration, diagram) = diagram

"""
    find_apparent_pairs(::AbstractFiltration, columns, is_cohomology, verbose)

Find apparent pairs in filtration. Columns are the initial columns that are being reduced.
Should return a tuple of:

* A collection of columns that are left to reduce.
* A collection of apparent pairs `(σ, τ)`. These will be insered into the reduced matrix so
  that `matrix[τ] = [σ]`.

!!! note
    Mutating `columns` and returning it from this function is safe.

!!! note
    This can only be applied to cohomology.

!!! warning
    This is still experimental.
"""
find_apparent_pairs(::AbstractFiltration, columns, _) = columns, ()

function index_overflow_check(
    flt::AbstractFiltration,
    ::Type{F},
    dim_max,
    message="$flt overflows in $(dim_max)D. " *
            "Try using a bigger index type or lower `dim_max`",
) where {F}
    S = simplex_type(flt, dim_max + 1)
    length(S) > nv(flt) && throw(ArgumentError("$S has more than $(nv) vertices."))
    return index_overflow_check(S, F, nv(flt), message)
end

"""
    distance_matrix(::AbstractFiltration)

Return a matrix with distances between vertices of the filtration. These distances are used
to determine edge length when finding shortest represenatative cycle in
[`reconstruct_cycle`](@ref).

Defaults to all distances being 1.

# Examples

```jldoctest
julia> flt = Rips([1 1 2; 1 0 1; 2 1 0]);

julia> Ripserer.distance_matrix(flt)
3×3 Matrix{Int64}:
 1  1  2
 1  0  1
 2  1  0

julia> flt = Custom([(1,2,3,4) => 10.0, (1,2) => 3.0, (1,3) => 4.0, (2,3) => 5.0]);

julia> Ripserer.distance_matrix(flt)
4×4 Ripserer.DefaultDist{Float64}:
 0.0  1.0  1.0  1.0
 1.0  0.0  1.0  1.0
 1.0  1.0  0.0  1.0
 1.0  1.0  1.0  0.0

```
"""
distance_matrix(flt::AbstractFiltration{<:Any,T}) where {T} = DefaultDist{T}(nv(flt))

struct DefaultDist{T} <: AbstractMatrix{T}
    nv::Int
end

Base.size(dd::DefaultDist) = (dd.nv, dd.nv)
Base.getindex(dd::DefaultDist, i, j) = ifelse(i == j, zero(eltype(dd)), one(eltype(dd)))


text:
function isnice(str::AbstractString)
    occursin(r"[aeiou].*[aeiou].*[aeiou]", str) &&
    occursin(r"(.)\1", str) &&
    !occursin(r"ab|cd|pq|xy", str)
end

@assert isnice("ugknbfddgicrmopn")
@assert isnice("aaa")
@assert !isnice("yzbqklnj")
@assert !isnice("haegwjzuvuyypxyu")
@assert !isnice("dvszwmarrgswjxmb")

println("Part 1: ", count(isnice, eachline(filename)))

isnice2(str::AbstractString) = occursin(r"(..).*\1", str) && occursin(r"(.).\1", str)

@assert isnice2("qjhvhtzxzqqjkmpb")
@assert isnice2("xxyxx")
@assert !isnice2("uurcxstgmygtbstg")
@assert !isnice2("ieodomkazucvgmuy")

println("Part 2: ", count(isnice2, eachline("input.txt")))

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
# ---
# title: 1010. Pairs of Songs With Total Durations Divisible by 60
# id: problem1010
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/pairs-of-songs-with-total-durations-divisible-by-60/description/>
# hidden: true
# ---
#
# You are given a list of songs where the ith song has a duration of `time[i]`
# seconds.
#
# Return _the number of pairs of songs for which their total duration in seconds
# is divisible by_ `60`. Formally, we want the number of indices `i`, `j` such
# that `i < j` with `(time[i] + time[j]) % 60 == 0`.
#
#
#
# **Example 1:**
#
#
#
#     Input: time = [30,20,150,100,40]
#     Output: 3
#     Explanation: Three pairs have a total duration divisible by 60:
#     (time[0] = 30, time[2] = 150): total duration 180
#     (time[1] = 20, time[3] = 100): total duration 120
#     (time[1] = 20, time[4] = 40): total duration 60
#
#
# **Example 2:**
#
#
#
#     Input: time = [60,60,60]
#     Output: 3
#     Explanation: All three pairs have a total duration of 120, which is divisible by 60.
#
#
#
#
# **Constraints:**
#
#   * `1 <= time.length <= 6 * 104`
#   * `1 <= time[i] <= 500`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using OptimizationDynamics
const iLQR = OptimizationDynamics.IterativeLQR
using LinearAlgebra
using Random

# ## visualization
vis = Visualizer()
render(vis);

# ## mode
MODE = :translate
MODE = :rotate

# ## gradient bundle
GB = false

# ## state-space model
h = 0.1
T = 26

im_dyn = ImplicitDynamics(planarpush, h, eval(r_pp_func), eval(rz_pp_func), eval(rθ_pp_func);
    r_tol=1.0e-8, κ_eval_tol=1.0e-4, κ_grad_tol=1.0e-2, nc=1, nb=9, info=(GB ? GradientBundle(planarpush, N=50, ϵ=1.0e-4) : nothing))

nx = 2 * planarpush.nq
nu = planarpush.nu

# ## iLQR model
ilqr_dyn = iLQR.Dynamics((d, x, u, w) -> f(d, im_dyn, x, u, w),
	(dx, x, u, w) -> GB ? fx_gb(dx, im_dyn, x, u, w) : fx(dx, im_dyn, x, u, w),
	(du, x, u, w) -> GB ? fu_gb(du, im_dyn, x, u, w) : fu(du, im_dyn, x, u, w),
	nx, nx, nu)

model = [ilqr_dyn for t = 1:T-1];

# ## initial conditions and goal
r_dim = 0.1
if MODE == :translate
	q0 = [0.0, 0.0, 0.0, -r_dim - 1.0e-8, 0.0]
	q1 = [0.0, 0.0, 0.0, -r_dim - 1.0e-8, 0.0]
	x_goal = 1.0
	y_goal = 0.0
	θ_goal = 0.0 * π
	qT = [x_goal, y_goal, θ_goal, x_goal - r_dim, y_goal - r_dim]
	xT = [qT; qT]
elseif MODE == :rotate
	q0 = [0.0, 0.0, 0.0, -r_dim - 1.0e-8, -0.01]
	q1 = [0.0, 0.0, 0.0, -r_dim - 1.0e-8, -0.01]
	x1 = [q1; q1]
	x_goal = 0.5
	y_goal = 0.5
	θ_goal = 0.5 * π
	qT = [x_goal, y_goal, θ_goal, x_goal-r_dim, y_goal-r_dim]
	xT = [qT; qT]
end

# ## objective
function objt(x, u, w)
	J = 0.0

	q1 = x[1:planarpush.nq]
	q2 = x[planarpush.nq .+ (1:planarpush.nq)]
	v1 = (q2 - q1) ./ h

	J += 0.5 * transpose(v1) * Diagonal([1.0, 1.0, 1.0, 0.1, 0.1]) * v1
	J += 0.5 * transpose(x - xT) * Diagonal([1.0, 1.0, 1.0, 0.1, 0.1, 1.0, 1.0, 1.0, 0.1, 0.1]) * (x - xT)
	J += 0.5 * (MODE == :translate ? 1.0e-1 : 1.0e-2) * transpose(u) * u

	return J
end

function objT(x, u, w)
	J = 0.0

	q1 = x[1:planarpush.nq]
	q2 = x[planarpush.nq .+ (1:planarpush.nq)]
	v1 = (q2 - q1) ./ h

	J += 0.5 * transpose(v1) * Diagonal([1.0, 1.0, 1.0, 0.1, 0.1]) * v1
	J += 0.5 * transpose(x - xT) * Diagonal([1.0, 1.0, 1.0, 0.1, 0.1, 1.0, 1.0, 1.0, 0.1, 0.1]) * (x - xT)

	return J
end

ct = iLQR.Cost(objt, nx, nu)
cT = iLQR.Cost(objT, nx, 0)
obj = [[ct for t = 1:T-1]..., cT];

# ## constraints
ul = [-5.0; -5.0]
uu = [5.0; 5.0]

function stage_con(x, u, w)
    [
     ul - u; # control limit (lower)
     u - uu; # control limit (upper)
    ]
end

function terminal_con(x, u, w)
    [
     (x - xT)[collect([(1:3)..., (6:8)...])]; # goal
    ]
end

cont = iLQR.Constraint(stage_con, nx, nu, idx_ineq=collect(1:(2 * nu)))
conT = iLQR.Constraint(terminal_con, nx, 0)
cons = [[cont for t = 1:T-1]..., conT];

# ## rollout
x1 = [q0; q1]
ū = MODE == :translate ? [t < 5 ? [1.0; 0.0] : [0.0; 0.0] for t = 1:T-1] : [t < 5 ? [1.0; 0.0] : t < 10 ? [0.5; 0.0] : [0.0; 0.0] for t = 1:T-1]
x̄ = iLQR.rollout(model, x1, ū)
q̄ = state_to_configuration(x̄)
visualize!(vis, planarpush, q̄, Δt=h);

# ## solver
solver = iLQR.solver(model, obj, cons,
	opts=iLQR.Options(
		linesearch = :armijo,
		α_min=1.0e-5,
		obj_tol=1.0e-3,
		grad_tol=1.0e-3,
		max_iter=10,
		max_al_iter=10,
		con_tol=0.005,
		ρ_init=1.0,
		ρ_scale=10.0,
		verbose=false))
iLQR.initialize_controls!(solver, ū)
iLQR.initialize_states!(solver, x̄);

# ## solve
iLQR.reset!(solver.s_data)
@time iLQR.solve!(solver);

@show iLQR.eval_obj(solver.m_data.obj.costs, solver.m_data.x, solver.m_data.u, solver.m_data.w)
@show solver.s_data.iter[1]
@show norm(terminal_con(solver.m_data.x[T], zeros(0), zeros(0)), Inf)
@show solver.s_data.obj[1] # augmented Lagrangian cost

# ## solution
x_sol, u_sol = iLQR.get_trajectory(solver)
q_sol = state_to_configuration(x_sol)
visualize!(vis, planarpush, q_sol, Δt=h);

# ## benchmark
solver.options.verbose = false
@benchmark iLQR.solve!($solver, x̄, ū) setup=(x̄=deepcopy(x̄), ū=deepcopy(ū));


text:
using Distributions
using Counters

"""
`binom_rv(n,p)` generates a single random value according to
the binomial distribution `B(n,p)`.
"""
binom_rv(n::Int,p::Real=0.5) = rand(Binomial(n,p))

"""
`binomial_counts(n,p,reps)` generates `reps` binomial random values `B(n,p)`
returning a `Counter` that reports how many times each value was observed.
"""
function binomial_counts(n::Int, p::Real, reps::Int)
  c = Counter{Int}()
  for k=1:reps
    x = binom_rv(n,p)
    c[x] += 1
  end
  return c
end

"""
`parallel_binomial_counts(n,p,reps,rounds)` generates `reps*rounds`
binomial random values `B(n,p)` returning a `Counter` that reports how
many times each value was observed. In this regard, this is the same
as `binomial_counts`. However, this is done by multiple processors by
making `rounds` calls to `binomial_counts(n,p,reps)` and combining
the results.
"""
function parallel_binomial_counts(n::Int, p::Real, reps::Int, rounds::Int)
  counts = @parallel (+) for k=1:rounds
    binomial_counts(n,p,reps)
  end
  return counts
end


text:
function push_msg!(a::Agent,msg::OFMessage)
    put!(a.queue,msg)
end

"""
It records messages pending to process in a given tick
"""
function push_pending!(a::AbstractAgent,msg::OFMessage)
    push!(a.pending,msg)
end

function send_msg!(receiver::Int64,msg::OFMessage,model)
    # log_info(model.ticks,receiver,"==> receiver queue: $(length(getindex(model,receiver).queue.data)) ===> max $(getindex(model,receiver).queue.sz_max)")
    ag = getindex(model,receiver)
    #TODO implement links and get delay of link in ticks
    queue = typeof(ag) == SimNE ? ag.queue : ag.queue
    # if model.ticks >= 30
        # log_info(model.ticks,msg.dpid,"Sent to $receiver msg: $msg ==> max: $(ag.queue.sz_max) current: $(length(ag.queue.data))")
    # end
    put!(queue,msg)
end

"""
    Function that makes a control agent trigger a flow install in the sne
"""
function install_flow!(a::Agent,path::Array{Int64,1},model::ABM,msg::OFMessage=nothing)
    # find which ones of path I am controlling

    es = get_controlled_assets(a.id,model)
    lpath = isempty(path) ? es : path
    eois = intersect(es,lpath)
    lpath = [ v for v in lpath]

    # log_info(model.ticks,a.id,"{$(get_controlled_assets(a.id,model))} install_flow! => path: $path -- es: $es -- eois: $eois - msg: -> $msg")

    #check that all nodes of path are up node are up if centralised
    is_install = model.ctrl_model == CENTRALISED ? length(collect(eois)) != length(path) ? false : true : true

    if is_install
        for e in eois
            i = length(lpath) > 1 ? first(indexin(e,lpath)) : 1
            sne = getindex(model,e)
            i_prev = i > 1 ? i - 1 : i

            ports = get_port_edge_list(sne)

            # log_info("[$(model.ticks)]{$(a.id)}($(sne.id)) - ports: $(ports) - i: $i - i_prev: $i_prev - e: $e -- lpath : $lpath")
            #Regardless of where the traffic comes
            r_src = "*" #string("h",first(lpath))
            r_dst = last(lpath)
            in_port = 0
            if i == 1
                in_port = msg.in_port
                #TODO of_msg remove from pending
            else
                in_port = first([ first(p) for p in ports if parse(Int,p[2][2:end]) == lpath[i_prev]])
            end
            out_port = 0

            if i < length(lpath)
                next_port = filter(p->parse(Int,p[2][2:end]) == lpath[i+1],ports)
                out_port = isempty(next_port) ? -1 : first(first(next_port))
            else
                #make sure it reaches the host regardles of the incoming port
                in_port = "*"
            end

            if out_port >= 0
                flow = Flow(  sne.id
                        ,MRule(string(in_port)
                        ,string(r_src)
                        ,string(r_dst))
                        ,[out_port]
                        ,OFS_Output)
                qid = msg.id
                install_msg = OFMessage(next_ofmid!(model), model.ticks,e,1,OFPR_ADD_FLOW,(flow=flow,qid=qid))
                send_msg!(e,install_msg,model)
            else
                log_info(model.ticks,a.id,"ERROR: Unable to find port")
            end

            next_sne = out_port > 0 ? lpath[i+1] : lpath[i]
            #  log_info(model.ticks,a.id,"--> ($(sne.id),$(next_sne))")
        end
    else
        log_info(model.ticks,a.id,"Install flow, path discarded: $path")
    end
end

function process_msg!(a::Agent,msg::OFMessage,model)

    @match msg.reason begin
        Ofp_Protocol(1) =>
                        begin
                            # log_info(model.ticks,a.id,"in_pkt -> $msg ===> prv_queries:  $(a.previous_queries)")
                            in_packet_handler(a,msg,model)
                        end
        Ofp_Protocol(2) =>
                            begin
                                # sneid_print = first(get_controlled_assets(a.id,model))
                                # sne_print = getindex(model,sneid_print)

                                # log_info(model.ticks,a.id,"BEFORE port_delete -> $msg ==>  $(get_state(sne_print).flow_table) ===> all ports: $(get_port_edge_list(sne_print))")

                                port_delete_handler(a,msg,model)

                                # log_info(model.ticks,a.id,"AFTER port_delete -> $msg ==>  $(get_state(sne_print).flow_table) ===> all ports: $(get_port_edge_list(sne_print))")
                            end
        Ofp_Protocol(5) => begin
                                node_rejoin_handler(a,msg,model)
                           end

        _ => begin
            log_info("[$(model.ticks)]($(a.id)) -> match default!")
            end
    end

    state = get_state(a)
    state.in_of_msg+=1.0
    set_state!(a,state)
end

"""
msg: SimNE.id, in_port, DPacket
"""
function in_packet_handler(a::Agent,msg::OFMessage,model)
    path::Array{Int64,1} = []
    found = false

    if msg.dpid != dst
        path = do_query!(msg,a,model)
        found = isempty(path) ? false : true
    else
        found = true
    end

    if found
        install_flow!(a,path,model,msg)
        get_state(a).active_paths[(first(path),last(path))] = path
    end
end

"""
    Get query exclusions for a given message
    @Deprecated
"""
function get_exclusions(msg::OFMessage,model)
    exc = Int64[]

    # Avoid returning original source node
    if msg.dpid != msg.data.src
        push!(exc,msg.data.src)
    end

    # If packet is not originated in a host connected to this sne
    if msg.in_port != 0
        #find out previous sne node
        ports = get_port_edge_list(getindex(model,msg.dpid))
        prv_id = parse(Int64,filter(x->x[1]==msg.in_port,ports)[1][2][2:end])
        # if previous node has not been already pushed to exc
        if msg.data.src != prv_id
            push!(exc,prv_id)
        end
    end

    return exc
end

function pending_pkt_handler(a::Agent,model)
    # if model.ticks in 80:1:90 && a.id == 10
    # end
    new_pending = []
    if !isempty(a.pending)
        # log_info(model.ticks,a.id,25,"BEFORE pending: $(length(a.pending))")
        for msgt in a.pending
            # log_info("[$(model.ticks)]($(a.id)) pending_msgt: $msgt")
            if msgt[3]
                put!(a.queue,msgt[2]) #msgt[2]: msg
            else
                push!(new_pending,msgt)
            end
         end
         a.pending = new_pending
        #  log_info(model.ticks,a.id,25,"AFTER pending: $(length(a.pending))")
    end
end

# """
#     Handle control actions when a control agent is dropped
# """
# function do_drop!(a::Agent,model)
    # TODO Check if this is useful. For time being, assuming
    # Control agent remains up but just forward messages.
#     lg = a.params[:ctl_graph]
#     lv = to_local_vertex(lg,a.id)
#     nbs = [ lg[nb,:aid] for nb in neighbors(lg,lv) ]

#     #Move connections to neighbors
#     for nb in nbs
#         nbs_nb = setdiff(nbs,[nb])
#         body = Dict(new_nbs=>nbs_nb,dp_aid=>a.id)
#         nb_msg = AGMessage(next_amid!(model),model.ticks,a.id,nb,NEW_NB,body)
#         send_msg!(nb,nb_msg_model)
#     end

#     #remove edges from "actual" control graph
#     soft_remove_vertex!(model.ctl_graph,a.id)

#     #TODO it can't set agent down, needs to check if any controlled asset is up.
#     set_down!(a)
# end

"""
Set controller agent down
"""
function set_down!(a::Agent)
    state = get_state(a)
    state.up = false
    set_state!(a,state)
end

"""
Set controller agent up
"""
function set_up!(a::Agent)
    state = get_state(a)
    state.up = true
    set_state!(a,state)
end

"""
    It processes OF msg sent by controlled NE  to
    remove a given port from its graph
    This is done by the agent controlling the neighbour of a dropped node
"""
function port_delete_handler(a::Agent,msg::OFMessage,model)
    # init_agent!(a,model)
    remove_drop_sne!(a,msg.data,model.ticks)
    # if model.ctrl_model != GraphModel(1)

    #     for ctl_p in a.ctl_paths
    #         log_info(model.ticks,a.id," ctl_p: $ctl_p")
    #         propagate_drop_sne!(a,msg.data,ctl_p,model)
    #     end

    # # else
    #     # lv = to_local_vertex(a.params[:ntw_graph],msg.data)
    #     # a.params[:ntw_graph] = soft_remove_vertex!(a.params[:ntw_graph],lv)

    # end
#flows involving this NE should have been deleted at NE
end

"""
It handles msg from sne when a neighbor sne has rejoined the underlying network.
"""
function node_rejoin_handler(a,msg,model)
   #update local graph
    src = first(msg.data)
    dst = last(msg.data)
    lvc_s = to_local_vertex(a.ntw_graph,src)
    lvc_d = to_local_vertex(a.ntw_graph,dst)
    lvb_s = to_local_vertex(a.base_ntw_graph,src)
    lvb_d = to_local_vertex(a.base_ntw_graph,dst)

    if lvb_s > 0 && lvb_d > 0
        if !has_edge(a.base_ntw_graph,lvb_s,lvb_d)
            add_edge!(a.base_ntw_graph,lvb_s,lvb_d)
        end
        if !has_edge(a.base_ntw_graph,lvb_d,lvb_s)
            add_edge!(a.base_ntw_graph,lvb_d,lvb_s)
        end
    end
    if lvc_s > 0 && lvc_d > 0
        if !has_edge(a.ntw_graph,lvc_s,lvc_d)
            add_edge!(a.ntw_graph,lvc_s,lvc_d)
        end
        if !has_edge(a.ntw_graph,lvc_d,lvc_s)
            add_edge!(a.ntw_graph,lvc_d,lvc_s)
        end
    end

    #trigger update paths
# log_info(model.ticks,a.id,"Existing paths: $(a.paths)")
# log_info(model.ticks,a.id,"rejoin msg received: $(msg)")
end

"""
Check confindence of a path
@Deprecated (for time being)
"""
function do_confidence_check!(a,model)
#TODO
end

function calculate_metrics_step!(a::Agent,model::ABM)
end


text:
module Structs

export Query

struct Query
    name::String
    data::Int
end

end


text:
# This collection of functionalities should probably be moved to a better home.

"""
    BoxCox

See "Learning non-Gaussian Time Series using the Box-Cox Gaussian Process".
"""
struct BoxCox{Tλ<:Real} <: Bijector{0}
    λ::Tλ
end

(b::BoxCox)(x) = (sign(x) * abs(x)^b.λ - 1) / b.λ
Bijectors.inv(b::BoxCox) = InvBoxCox(b.λ)
Bijectors.logabsdetjac(b::BoxCox, x) = (b.λ - 1) *  log(abs(x))

"""
    InvBoxCox

See "Learning non-Gaussian Time Series using the Box-Cox Gaussian Process".
"""
struct InvBoxCox{Tλ<:Real} <: Bijector{0}
    λ::Tλ
end

(b::InvBoxCox)(y) = sign(b.λ * y + 1) * abs(b.λ * y + 1)^(1 / b.λ)
Bijectors.inv(b::InvBoxCox) = BoxCox(b.λ)
Bijectors.logabsdetjac(b::InvBoxCox, y) = Bijectors.logabsdetjac(inv(b), inv(b)(y))

affine(a, b) = Bijectors.Shift(a) ∘ Bijectors.Scale(b)


text:
using FlightGNC
using UnPack, ComponentArrays, Transducers
using LinearAlgebra, DifferentialEquations  # for callbacks
using Plots

function main(p_M_0::Vector, v_M_0::Vector, p_T_0::Vector, v_T_0::Vector, s_guidance::BPNG; A_M_max = Inf)
    @unpack dim, σ_M_lim = s_guidance
    # Environments
    pursuer = PointMassMissile(dim, A_M_max)
    evador = PointMassMissile(dim, Inf)
    env = PursuerEvadorMissile(pursuer, evador)

    # Initial condition
    x0_pursuer = State(pursuer)(p_M_0, v_M_0)
    x0_evador = State(evador)(p_T_0, v_T_0)
    x0 = State(env)(x0_pursuer, x0_evador)

    # Simulation parameters
    Δt      = 0.01
    t_sim_f = 300

    # callbacks
    function condition_stop(σ_M_lim)
        return function (u, t, integrator)
            (p_M, v_M, p_T, v_T) = (u.pursuer.p, u.pursuer.v, u.evador.p , u.evador.v)
            r = norm(p_T-p_M)
            ṙ = dot(p_T-p_M, v_T-v_M) / r
            σ_M = acos( clamp( dot(normalize(p_T-p_M), normalize(v_M)), -1, 1 ) )

            (r < 2.5)  || (r < 4 && ṙ >= 0) # || (σ_M > σ_M_lim)
        end
    end
    affect!(integrator) = terminate!(integrator)  # See DiffEq.jl documentation
    cb_stop    = DiscreteCallback(condition_stop(σ_M_lim), affect!)
    cb = CallbackSet(cb_stop)  # useful for multiple callbacks

    # Execute Simulation
    # prob: DE problem, df: DataFrame
    simulator = Simulator(x0,
                        apply_inputs(Dynamics!(env, s_guidance);
                                    u_pursuer = BPNG_cmd(s_guidance),
                                    u_evador = (x, params, t) -> zeros(evador.dim));
                        Problem = :ODE,
                        solver = Tsit5(),
                        tf = t_sim_f)

    # Non-interactive simulation: solve approach (automatically reinitialised)
    @time df = solve(simulator;
                    savestep = Δt,
                    callback = cb,
                    reltol = 1e-8,
                    abstol = 1e-8)

    # @time prob, df = FSimBase.sim(
    #                      x0,  # initial condition
    #                      apply_inputs(Dynamics!(env, s_guidance);
    #                                   u_pursuer = BPNG_cmd(s_guidance),
    #                                   u_evador = (x, params, t) -> zeros(evador.dim));  # dynamics!; apply_inputs is exported from FS and is so useful for systems with inputs
    #                      tf = t_sim_f,
    #                      savestep = Δt,  # savestep is NOT simulation step
    #                      solver = Tsit5(),
    #                      callback = cb,
    #                      reltol = 1e-8,
    #                      abstol = 1e-8
    #                     )  # sim is exported from FS

	return df
end

function sim_plot(p_M_0::Vector, v_M_0::Vector, p_T_0::Vector, v_T_0::Vector, s_BPNG::BPNG, A_M_max::Number)
    # Execute simulation
    df = main(p_M_0, v_M_0, p_T_0, v_T_0, s_BPNG; A_M_max = A_M_max)

    ts = df.time
    p_Ms = df.sol |> Map(datum -> datum.pursuer.p) |> collect
    p_Ts = df.sol |> Map(datum -> datum.evador.p)  |> collect

    p_Ms = hcat(p_Ms...)'
    p_Ts = hcat(p_Ts...)'

    # Plotting
    legend_string = ["Missile" "Target"]

    if dim == 2
        f = fig_print([p_Ms[:,1] p_Ts[:,1]], [p_Ms[:,2] p_Ts[:,2]], "", legend_string, "x [m]", "y [m]"; ar_val = :equal, save_file = 0)
        # plot!(f, xlims=(0, 6E3), ylims=(-3E3, 3E3))

        # f = fig_print(ts, [p_Ms[:,1] p_Ts[:,2]], "", ["x" "y"], "t [s]", "pos [m]"; save_file = 0)

    elseif dim == 3
        f_2D= fig_print([p_Ms[:,1] p_Ts[:,1]], [p_Ms[:,2] p_Ts[:,2]], "", legend_string, "x [m]", "y [m]"; ar_val = :equal, save_file = 0)
        # display(f_2D)

        # f_3D = plot([p_Ms[:,1] p_Ts[:,1]], [p_Ms[:,2] p_Ts[:,2]], [p_Ms[:,3] p_Ts[:,3]], label = legend_string,
        #         xlabel = "x [m]", ylabel = "y [m]", zlabel = "z [m]", aspect_ratio = :equal,
        #         camera =(45,30)
        #         )
        # f_3D = equal_AR_3D([p_Ms[:,1] p_Ts[:,1]], [p_Ms[:,2] p_Ts[:,2]], [p_Ms[:,3] p_Ts[:,3]], [], f_3D)

        f_3D = equal_AR_3D([p_Ms[:,1] p_Ts[:,1]], [p_Ms[:,2] p_Ts[:,2]], [p_Ms[:,3] p_Ts[:,3]], legend_string)
        plot!(f_3D, xlabel = "x [m]", ylabel = "y [m]", zlabel = "z [m]", camera = (45,30))

        f = plot(f_2D, f_3D, layout = (2,1), size = (600, 1200))
    end

    f = fig_print([], [], "Trajectory", [], [], f; ar_val = :equal) 	# To just save the result
    # display(f)

    return df, f
end

# -----------------------------------------------------------------------------------------
# Simulation Input Default Values: Initial condition / Design parameters
dim     = 3

p_M_0   = zeros(dim)

V_M_0   = 300
γ_M_0 	= deg2rad(45)
χ_M_0 	= deg2rad(90)

p_T_0   = [10E3; 0; 0] #[10E3; 5E3; 5E3]
v_T_0   = zeros(dim) # 100*[-1; 0; 0]

γ_f_d   = deg2rad(-90)
χ_f_d 	= deg2rad(180)

if dim == 2
    γ_M_0   = 0
    γ_f_d   = 0
end

v_M_0   = V_M_0*[cos(γ_M_0)*sin(χ_M_0); cos(γ_M_0)*cos(χ_M_0); sin(γ_M_0)]
v̂_f_d   = [cos(γ_f_d)*sin(χ_f_d); cos(γ_f_d)*cos(χ_f_d); sin(γ_f_d)]

if dim == 2
    (p_M_0, v_M_0, p_T_0, v_T_0) = (p_M_0, v_M_0, p_T_0, v_T_0) .|> x -> x[1:2]
end

σ_M_lim = deg2rad(90)
A_M_max = 300

N       = 3
s_Bias  = ComponentArray(α = 1, δ = 0.01, n = 1, r_ref = 10E3, k = 3, m = 0, k̂_d = [1; 0; 0], i_Ω_μ = 0, i_σ_M_lim = 1)

s_BPNG  = BPNG(N, dim, σ_M_lim, v̂_f_d, Bias_IACG_StationaryTarget, s_Bias)
# Bias options: Bias_zero, Bias_IACG_StationaryTarget, Bias_IACG_StationaryTarget_2D


text:
#
# TODO: should ITensorMap be a special version of
# an ITensorNetwork with input and output indices specified?
#

# T is however the nodes are indexed
# TODO: how to deal with 2D, multiple networks, etc.?
# struct IndexSetNetwork{T}
#   # Use Vector{SortedVector{Pair{T, IndexSet}}}
#   data::Vector{Vector{Pair{T, IndexSet}}}
# end

# Make MPS, PEPS, etc. wrappers around `ITensorNetwork`,
# which could be subclasses of `AbstractITensorNetwork`
#
# Also `InfiniteITensorNetwork`, `AbstractInfiniteITensorNetwork`:
#
# struct InfiniteMPS <: AbstractInfiniteITensorNetwork
#   tensornetwork::InfiniteITensorNetwork
# end
#
# struct ITensorNetwork <: AbstractITensorNetwork
#   itensors::Vector{ITensor}
#   indnetwork::IndexSetNetwork # Adjacency list of IndexSets
#   input_inds::IndexSet
#   output_inds::IndexSet
# end
#
# TODO: how to deal with networks of ITensorNetwork,
# for example a network of MPS and MPO?
# ITensorNetworkNetwork that is a tree of ITensorNetwork?
#

abstract type AbstractITensorMap end

input_inds(T::AbstractITensorMap) = T.input_inds
output_inds(T::AbstractITensorMap) = T.output_inds
(T::AbstractITensorMap)(v::ITensor) = replaceinds(T * v, output_inds(T) => input_inds(T))

# convert from Tuple to Vector
tuple_to_vector(t::Tuple) = collect(t)
tuple_to_vector(v::Vector) = v

# Represents the action of applying the
# vector of ITensors to a starting state and then mapping
# them back (from output_inds to input_inds)
# TODO: rename ITensorNetworkMap?
# TODO: maybe parametrize the type to allow storing just 1 ITensor?
# TODO: contraction order optimization!
struct ITensorMap <: AbstractITensorMap
  itensors::Vector{ITensor}
  scalar::Number
  input_inds::Vector{Index}
  output_inds::Vector{Index}
  function ITensorMap(itensors::Vector{ITensor}, scalar, input_inds, output_inds)
    inds_in = tuple_to_vector(input_inds)
    inds_out = tuple_to_vector(output_inds)
    #inds_eltype = promote_type(eltype(input_inds), eltype(output_inds))
    #return new{inds_eltype}(itensors, inds_in, inds_out)
    return new(itensors, scalar, inds_in, inds_out)
  end
end

function ITensorMap(itensors::Vector{ITensor}, input_inds, output_inds)
  return ITensorMap(itensors, true, input_inds, output_inds)
end

function (M1::ITensorMap * M2::ITensorMap)
  # TODO: check the directions are correct
  @assert output_inds(M2) == input_inds(M1)
  return ITensorMap(
    vcat(M1.itensors, M2.itensors), M1.scalar * M2.scalar, input_inds(M2), output_inds(M1)
  )
end

function default_input_inds(itensors::Vector{ITensor})
  return filter(i -> plev(i) == 0, noncommoninds(itensors...))
end

function ITensorMap(
  itensors::Vector{ITensor};
  input_inds=default_input_inds(itensors),
  output_inds=dag(input_inds'),
)
  return ITensorMap(itensors, input_inds, output_inds)
end

Base.iterate(T::ITensorMap, args...) = iterate(T.itensors, args...)

function Base.transpose(T::ITensorMap)
  return ITensorMap(reverse(T.itensors), output_inds(T), input_inds(T))
end

# This is actually a Hermitian conjugation, not priming
function Base.adjoint(T::ITensorMap)
  return ITensorMap(reverse(dag.(T.itensors)), dag(output_inds(T)), dag(input_inds(T)))
end

# TODO: make a default constructor that searches for pairs of primed and unprimed indices
# TODO: this would be useful for ITensor matrix factorizations!
# function ITensorMap(itensors::Vector{ITensor}, pair_match::Pair = 0 => 1)
#   # pair_match could be:
#   # pair_match = 0 => 2
#   # pair_match = "Input" => "Output"
#   # pair_match = ("Input", 0) => ("Output", 1)
#   external_inds = unioninds(siteinds(itensors)...)
#   input_match = first(pair_match)
#   output_match = first(pair_match)
#   # TODO: preprocess pair_match to be of the form `(tags, plev)`
#   ind_match(i::Index, match) = hastags(i, match[1]) && hasplev(i, match[2])
#   input_inds = filter(i -> ind_match(i, input_match), external_inds)
#   output_inds = dag(replaceprime(replacetags(input_inds, first.(pair_match)), last.(pair_match)))
#   @assert = hassameinds(output_inds, filter(i -> ind_match(i, output_match), external_inds))
#   return ITensorMap(itensors, input_inds, output_inds)
# end

function set_scalar(T::ITensorMap, scalar::Number)
  return ITensorMap(T.itensors, scalar, input_inds(T), output_inds(T))
end

# Lazily scale by a scalar
(T::ITensorMap * c::Number) = set_scalar(T, T.scalar * c)
(c::Number * T::ITensorMap) = set_scalar(T, c * T.scalar)

-(T::ITensorMap) = set_scalar(T, -T.scalar)

# TODO: assert isempty(uniqeinds(v, T))?
# Apply the operator as T|v⟩
(T::ITensorMap * v::ITensor) = T.scalar * contract(pushfirst!(copy(T.itensors), v)) #*(v, T...)
# Apply the operator as ⟨v̅|T (simple left multiplication, without conjugation)
# This applies the ITensorMap tensors in reverse (maybe this is not always the best contraction
# ordering)
(v::ITensor * T::ITensorMap) = T.scalar * contract(pushfirst!(reverse(T.itensors, v))) #*(v, reverse(T)...)

# TODO: implement Base.iterate(::Base.Iterators.Reverse{MPS})

# TODO: use something like:
# neighbors(ψ, ϕ, (1, 1))
# neighbors(ψ, ϕ, (2, 1))
# neighbors(ψ, ϕ, (1, N))
# neighbors(ψ, ϕ, (2, N))
# Transfer matrix made from two MPS: T|v⟩ -> |w⟩
function ITensorMap(ψ::MPS, ϕ::MPS; input_inds=nothing, output_inds=nothing)
  N = length(ψ)
  @assert length(ϕ) == N
  itensors::Vector{ITensor} = reverse(collect(Iterators.flatten(Iterators.zip(ψ, ϕ))))
  if isnothing(input_inds)
    input_inds = unioninds(
      uniqueinds(ψ[N], ψ[N - 1], ϕ[N]), uniqueinds(ϕ[N], ϕ[N - 1], ψ[N])
    )
  end
  if isnothing(output_inds)
    output_inds = unioninds(uniqueinds(ψ[1], ψ[2], ϕ[1]), uniqueinds(ϕ[1], ϕ[2], ψ[1]))
  end
  return ITensorMap(itensors, input_inds, output_inds)
end

# Represents a sum of ITensor maps
struct ITensorMapSum <: AbstractITensorMap
  itensormaps::Vector{ITensorMap}
  input_inds::Vector{Index}
  output_inds::Vector{Index}
  function ITensorMapSum(itensormaps::Vector{ITensorMap})
    # TODO: check that all input_inds and output_inds are the same
    return new(itensormaps, input_inds(first(itensormaps)), output_inds(first(itensormaps)))
  end
end

(M1::ITensorMap + M2::ITensorMap) = ITensorMapSum([M1, M2])
(M1::ITensorMapSum + M2::ITensorMap) = ITensorMapSum(push!(copy(M1.itensormaps), M2))
(M1::ITensorMap + M2::ITensorMapSum) = M2 + M1
(M1::ITensorMap - M2::ITensorMap) = M1 + (-M2)
(M1::ITensorMapSum - M2::ITensorMap) = M1 + (-M2)
(M1::ITensorMap - M2::ITensorMapSum) = M1 + (-M2)
(M::ITensorMapSum * v::ITensor) = sum([m * v for m in M.itensormaps])

(M::ITensorMapSum * c::Number) = ITensorMapSum([m * c for m in M.itensormaps])
(c::Number * M::ITensorMapSum) = M * c
-(M::ITensorMapSum) = -1 * M

# Required by Arpack.eigs
Base.size(T::AbstractITensorMap) = (dim(output_inds(T)), dim(input_inds(T)))
LinearAlgebra.issymmetric(::AbstractITensorMap) = false
# TODO: promote the element types of all of the ITensors
# in the ITensorMap
Base.eltype(T::AbstractITensorMap) = eltype(T.itensors[1])
function LinearAlgebra.mul!(y, A::AbstractITensorMap, x)
  xt = itensor(x, dag(input_inds(A)); tol=1e-15)
  # XXX: something wrong with in-place ITensor
  # contraction, not overwriting input data
  yt = A(xt)
  yt = permute(yt, input_inds(A))
  y .= vec(array(yt))
  return y
end


text:
module TestMean
    using StreamStats
    using Distributions
    using Base.Test

    # Means of uniform draws
    for n in rand(1:1_000_000, 100)
        xs = rand(n)
        stat = StreamStats.Mean()
        for x in xs
            update!(stat, x)
        end
        online_m = mean(stat)
        online_ms = state(stat)
        online_n = nobs(stat)
        batch_m = mean(xs)
        @test_approx_eq(online_m, batch_m)
        @test_approx_eq(online_ms, batch_m)
        @test online_n == n
    end

    # Means of non-uniform draws
    ds = [Normal(0, 1), Normal(1e5, 1e8), LogNormal(10, 10), Cauchy(0.0, 1e-8)]
    for d in ds
        for n in rand(1:1_000_000, 10)
            xs = rand(d, n)
            stat = StreamStats.Mean()
            for x in xs
                update!(stat, x)
            end
            online_m = mean(stat)
            online_n = nobs(stat)
            batch_m = mean(xs)
            @test abs(online_m - batch_m) < 1e-8 * abs(batch_m)
            @test online_n == n
        end
    end
end


text:
using KendallTau
using Test

X = Float64[1 0; 2 1; 3 0; 4 1; 5 10]
Y = Float64[5 5 6; 3 4 1; 4 0 4; 2 6 1; 5 7 10]
Xm = [1 0; missing 1; 2 1; 3 0; 4 1; 5 10]
Ym = [5 5 6; 1 2 3; 3 4 1; 4 0 4; 2 6 1; 5 7 10]

x1 = X[:,1]
x2 = X[:,2]
y = Y[:,1]

#= corspearman

@test corspearman(x1, y) ≈ -0.102597835208515
@test corspearman(x2, y) ≈ -0.081110710565381

@test corspearman(X, y) ≈ [-0.102597835208515, -0.081110710565381]
@test corspearman(y, X) ≈ [-0.102597835208515 -0.081110710565381]

c11 = corspearman(x1, x1)
c12 = corspearman(x1, x2)
c22 = corspearman(x2, x2)
@test c11 ≈ 1.0
@test c22 ≈ 1.0
@test corspearman(X, X) ≈ [c11 c12; c12 c22]
@test corspearman(X)    ≈ [c11 c12; c12 c22]

@test corspearman(X, Y) ==
     [corspearman(X[:,i], Y[:,j]) for i in axes(X, 2), j in axes(Y, 2)] =#

# corkendall

# Check error, handling of NaN, Inf etc
@test_throws DimensionMismatch("Vectors must have same length") corkendall([1,2,3,4], [1,2,3])
@test isnan(corkendall([1,2], [3,NaN]))
@test isnan(corkendall([1,1,1], [1,2,3]))
@test corkendall([-Inf,-0.0,Inf], [1,2,3]) == 1.0

# Test, with exact equality, some known results.
# RealVector, RealVector
@test corkendall(x1, y) == -1 / sqrt(90)
@test corkendall(x2, y) == -1 / sqrt(72)
# RealMatrix, RealVector
@test corkendall(X, y)  == [-1 / sqrt(90), -1 / sqrt(72)]
# RealVector, RealMatrix
@test corkendall(y, X)  == [-1 / sqrt(90) -1 / sqrt(72)]

# n = 78_000 tests for overflow errors on 32 bit
# Testing for overflow errors on 64bit would require n be too large for practicality
# This also tests merge_sort! since n is (much) greater than SMALL_THRESHOLD.
n = 78_000
# Test with many repeats
@test corkendall(repeat(x1, n), repeat(y, n)) ≈ -1 / sqrt(90)
@test corkendall(repeat(x2, n), repeat(y, n)) ≈ -1 / sqrt(72)
@test corkendall(repeat(X, n), repeat(y, n))  ≈ [-1 / sqrt(90), -1 / sqrt(72)]
@test corkendall(repeat(y, n), repeat(X, n))  ≈ [-1 / sqrt(90) -1 / sqrt(72)]
@test corkendall(repeat([0,1,1,0], n), repeat([1,0,1,0], n)) == 0.0

# Test with no repeats, note testing for exact equality
@test corkendall(collect(1:n), collect(1:n))          == 1.0
@test corkendall(collect(1:n), reverse(collect(1:n))) == -1.0

# All elements identical should yield NaN
@test isnan(corkendall(repeat([1], n), collect(1:n)))

# Test handling of missings
@test corkendall(Xm, Xm, skipmissing=:pairwise) == corkendall(X, X)
@test corkendall(Xm, Xm, skipmissing=:complete) == corkendall(X, X)
@test corkendall(Xm, Ym, skipmissing=:complete) == corkendall(X, Y)
@test corkendall(Xm, Ym, skipmissing=:pairwise) ≈  [-1/√90 0.4 1/√90;
                                                  -2/√154 7/√165 -1/√154]

# Test not-correct values of skipmissing
@test_throws ArgumentError corkendall(Xm)
@test_throws ArgumentError corkendall(X, skipmissing=:foo)

c11 = corkendall(x1, x1)
c12 = corkendall(x1, x2)
c22 = corkendall(x2, x2)

# RealMatrix, RealMatrix
@test corkendall(X, X) ≈ [c11 c12; c12 c22]
# RealMatrix
@test corkendall(X)    ≈ [c11 c12; c12 c22]

@test c11 == 1.0
@test c22 == 1.0
@test c12 == 3 / sqrt(20)

# Finished testing for overflow, so redefine n for speedier tests
n = 100

@test corkendall(repeat(X, n), repeat(X, n)) ≈ [c11 c12; c12 c22]
@test corkendall(repeat(X, n))               ≈ [c11 c12; c12 c22]

# All eight three-element permutations
z = [1  1  1;
     1  1  2;
     1  2  2;
     1  2  2;
     1  2  1;
     2  1  2;
     1  1  2;
     2  2  2]

@test corkendall(z)         == [1 0 1 / 3; 0 1 0; 1 / 3 0 1]
@test corkendall(z, z)      == [1 0 1 / 3; 0 1 0; 1 / 3 0 1]
@test corkendall(z[:,1], z) == [1 0 1 / 3]
@test corkendall(z, z[:,1]) == [1; 0; 1 / 3]

z = float(z)
@test corkendall(z)         == [1 0 1 / 3; 0 1 0; 1 / 3 0 1]
@test corkendall(z, z)      == [1 0 1 / 3; 0 1 0; 1 / 3 0 1]
@test corkendall(z[:,1], z) == [1 0 1 / 3]
@test corkendall(z, z[:,1]) == [1; 0; 1 / 3]

w = repeat(z, n)
@test corkendall(w)         == [1 0 1 / 3; 0 1 0; 1 / 3 0 1]
@test corkendall(w, w)      == [1 0 1 / 3; 0 1 0; 1 / 3 0 1]
@test corkendall(w[:,1], w) == [1 0 1 / 3]
@test corkendall(w, w[:,1]) == [1; 0; 1 / 3]

KendallTau.midpoint(1,10)        == 5
KendallTau.midpoint(1,widen(10)) == 5

# NaN handling

Xnan = copy(X)
Xnan[1,1] = NaN
Ynan = copy(Y)
Ynan[2,1] = NaN

# for f in (corspearman, corkendall)
for f in (corkendall,)
    @test isnan(f([1.0, NaN, 2.0], [2.0, 1.0, 3.4]))
    @test all(isnan, f([1.0, NaN], [1 2; 3 4]))
    @test all(isnan, f([1 2; 3 4], [1.0, NaN]))
    @test isequal(f([1 NaN; NaN 4]), [1 NaN; NaN 1])
    @test all(isnan, f([1 NaN; NaN 4], [1 NaN; NaN 4]))
    @test all(isnan, f([1 NaN; NaN 4], [NaN 1; NaN 4]))

    @test isequal(f(Xnan, Ynan),
                   [f(Xnan[:,i], Ynan[:,j]) for i in axes(Xnan, 2), j in axes(Ynan, 2)])
    @test isequal(f(Xnan),
                   [i == j ? 1.0 : f(Xnan[:,i], Xnan[:,j])
                    for i in axes(Xnan, 2), j in axes(Xnan, 2)])
     for k in 1:2
          @test isequal(f(Xnan[:,k], Ynan),
                        [f(Xnan[:,k], Ynan[:,j]) for i in 1:1, j in axes(Ynan, 2)])
          # TODO: fix corkendall (PR#659)
          if f === corkendall
               @test isequal(f(Xnan, Ynan[:,k]),
                             [f(Xnan[:,i], Ynan[:,k]) for i in axes(Xnan, 2)])
          else
               @test isequal(f(Xnan, Ynan[:,k]),
                             [f(Xnan[:,i], Ynan[:,k]) for i in axes(Xnan, 2), j in 1:1])
          end
     end

     # Wrong dimensions
     @test_throws DimensionMismatch f([1], [1, 2])
     @test_throws DimensionMismatch f([1], [1 2; 3 4])
     @test_throws DimensionMismatch f([1 2; 3 4], [1])
     @test_throws ArgumentError f([1 2; 3 4:4 6], [1 2; 3 4])

end


text:
function do_exp_plot(id1, id2)

    tot_p = nothing
    for (iJR, Data) in [
        (ChK.iJR904, ChK.KayserData),
        (ChF.iJR904, ChF.FolsomData),
        (ChN.iJR904, ChN.NanchenData),
        # (ChH.iJR904, ChH.HeerdenData),
    ]

        src = nameof(Data)
        # unit1 = Data.unit(id1)
        # unit2 = Data.unit(id2)

        DAT = ChE.load_DAT(src)
        EXPS = DAT[:EXPS]

        local_p = plot(;
            xlabel = string(id1),
            ylabel = string(id2)
        )
        tot_p = isnothing(tot_p) ? deepcopy(local_p) : tot_p

        xs = abs.([Data.val(id1, exp, 0.0) for exp in EXPS])
        sidxs = sortperm(xs)
        xs = xs[sidxs]
        ys = abs.([Data.val(id2, exp, 0.0) for exp in EXPS][sidxs])

        marker = (8, source_markers[Data])
        for p in [local_p, tot_p]
            scatter!(p, xs, log.(ys);
                label = "", marker,
                alpha = 0.8, color = :black
            )
            # plot!(p, xs, ys;
            #     label = "", lw = 2, ls = :dash,
            #     alpha = 0.8, color = :black
            # )
        end

        plot!(local_p; title = string(src))
        mysavefig(local_p, "experimental_corr"; id1, id2, src)
    end

    mysavefig(tot_p, "experimental_corr"; id1, id2)
end

## ---------------------------------------------------------------------------
let
    # ids = [:sGLC, :sAC, :D, :X, :uGLC, :uAC]
    ids = [:D, :X]
    for i in 1:length(ids)
        for j in (i + 1):length(ids)
            do_exp_plot(ids[i], ids[j])
        end
    end
end

## ---------------------------------------------------------------------------
let
    fontsize = 13
    commom_params = (;
        # dpi = 1000,
        thickness_scaling = 1.3,
        xguidefontsize = fontsize, yguidefontsize = fontsize,
        ylim = [0.0, Inf]
    )

    tot_p = nothing
    for (iJR, Data) in [
        (ChK.iJR904, ChK.KayserData),
        (ChF.iJR904, ChF.FolsomData),
        (ChN.iJR904, ChN.NanchenData),
        # (ChH.iJR904, ChH.HeerdenData),
    ]

        src = nameof(Data)

        DAT = ChE.load_DAT(src)
        EXPS = DAT[:EXPS]
        # EXPS = Data.EXPS

        local_p = plot(;
            xlabel = "D (1/ h)",
            ylabel = "sGLC / cGLC",
            commom_params...
        )
        tot_p = isnothing(tot_p) ? deepcopy(local_p) : tot_p

        Ds = abs.([Data.val(:D, exp, 0.0) for exp in EXPS])
        sidxs = sortperm(Ds)
        Ds = Ds[sidxs]
        sGLC = abs.([Data.val(:sGLC, exp, 0.0) for exp in EXPS][sidxs])
        cGLC = abs.([Data.val(:cGLC, exp, 0.0) for exp in EXPS][sidxs])

        marker = (8, source_markers[Data])
        for p in [local_p, tot_p]

            # sGLC
            scatter!(p, Ds, sGLC ./ cGLC;
                label = "", marker,
                alpha = 0.8, color = :black
            )
            plot!(p, Ds, sGLC ./ cGLC;
                label = "", lw = 2, ls = :dash,
                alpha = 0.8, color = :black
            )
        end

        plot!(local_p; title = string(src))
        mysavefig(local_p, "exp_GLC_vs_D"; src)
    end

    mysavefig(tot_p, "exp_GLC_vs_D")
end


text:
using SolveDSGE

filename = "model5a.txt"
path5 = joinpath(@__DIR__,filename)
process_model(path5)

processed_filename = "model5a_processed.txt"
processed_path5 =  joinpath(@__DIR__,processed_filename)

dsge5 = retrieve_processed_model(processed_path5)

x5 = [0.05, 0.05, -0.05, 0.3, 0.3, 0.3, 0.005, 0.005]

tol = 1e-8
maxiters = 1000

ss5 = compute_steady_state(dsge5,x5,tol,maxiters)
ss5 = ss5.zero

N   = PerturbationScheme(ss5,1.0,"first")
NN  = PerturbationScheme(ss5,1.0,"second")
NNN = PerturbationScheme(ss5,1.0,"third")

soln_fo5  = solve_model(dsge5,N)
soln_so5  = solve_model(dsge5,NN)
soln_to5  = solve_model(dsge5,NNN)

P = ChebyshevSchemeStoch(ss5,chebyshev_nodes,[11,11,11,11],9,4,[0.03 0.13 0.01 0.34;-0.03 -0.13 -0.021 0.31],tol,1e-5,maxiters)
PP = ChebyshevSchemeStoch(ss5,chebyshev_nodes,[21,21,21,21],9,4,[0.03 0.13 0.01 0.34;-0.03 -0.13 -0.021 0.31],tol,1e-5,maxiters)
PPP = ChebyshevSchemeStoch(ss5,chebyshev_nodes,[31,31,31,31],9,5,[0.03 0.13 0.01 0.34;-0.03 -0.13 -0.021 0.31],tol,1e-5,maxiters)

soln_nl5a = solve_model(dsge5,P)
soln_nl5b = solve_model(dsge5,soln_to5,PP)
soln_nl5c = solve_model(dsge5,soln_nl5b,PPP)

L = SmolyakSchemeStoch(ss5,chebyshev_gauss_lobatto,9,3,[0.03 0.13 0.01 0.34;-0.03 -0.13 -0.021 0.31],tol,1e-5,maxiters)
LL = SmolyakSchemeStoch(ss5,chebyshev_gauss_lobatto,9,4,[0.03 0.13 0.01 0.34;-0.03 -0.13 -0.021 0.31],tol,1e-5,maxiters)
LLL = SmolyakSchemeStoch(ss5,chebyshev_gauss_lobatto,9,5,[0.03 0.13 0.01 0.34;-0.03 -0.13 -0.021 0.31],tol,1e-5,maxiters)

soln_nl5d = solve_model(dsge5,L)
soln_nl5e = solve_model(dsge5,soln_to5,LL)
soln_nl5f = solve_model(dsge5,soln_nl5e,LLL)

M = PiecewiseLinearSchemeStoch(ss5,[11,11,11,11],9,[0.03 0.13 0.01 0.34;-0.03 -0.13 -0.021 0.31],1e-3,1e-5,maxiters)
MM = PiecewiseLinearSchemeStoch(ss5,[21,21,21,21],9,[0.03 0.13 0.01 0.34;-0.03 -0.13 -0.021 0.31],tol,1e-5,maxiters)

soln_nl5g = solve_model(dsge5,M)
soln_nl5h = solve_model(dsge5,soln_to5,MM)
soln_nl5i = solve_model(dsge5,soln_nl5b,MM)

simulated_data5 = simulate(soln_nl5c,ss5[1:4],100)
pos_imps, neg_imps = impulses(soln_nl5c,50,[1],10000)


text:
module DiffEqSensitivity

using DiffEqBase, ForwardDiff, Tracker, FiniteDiff, Statistics
using DiffEqCallbacks, QuadGK, RecursiveArrayTools, LinearAlgebra
using RecursiveArrayTools, QuasiMonteCarlo, Adapt
using Parameters: @unpack, @with_kw
using FFTW, Distributions, Requires
using StochasticDiffEq
using SharedArrays
import DiffEqNoiseProcess
import ZygoteRules, Zygote, ReverseDiff

abstract type SensitivityFunction end
abstract type GSAMethod end
abstract type TransformedFunction end

include("require.jl")
include("local_sensitivity/sensitivity_algorithms.jl")
include("local_sensitivity/derivative_wrappers.jl")
include("local_sensitivity/sensitivity_interface.jl")
include("local_sensitivity/forward_sensitivity.jl")
include("local_sensitivity/adjoint_common.jl")
include("local_sensitivity/backsolve_adjoint.jl")
include("local_sensitivity/interpolating_adjoint.jl")
include("local_sensitivity/quadrature_adjoint.jl")
include("local_sensitivity/concrete_solve.jl")
include("local_sensitivity/second_order.jl")
include("local_sensitivity/steadystate_adjoint.jl")
include("local_sensitivity/sde_tools.jl")
include("global_sensitivity/morris_sensitivity.jl")
include("global_sensitivity/sobol_sensitivity.jl")
include("global_sensitivity/regression_sensitivity.jl")
include("global_sensitivity/DGSM_sensitivity.jl")
include("global_sensitivity/eFAST_sensitivity.jl")

export extract_local_sensitivities

export ODEForwardSensitivityFunction, ODEForwardSensitivityProblem, SensitivityFunction,
       ODEAdjointSensitivityProblem, ODEAdjointProblem, AdjointSensitivityIntegrand,
       SDEAdjointProblem,
       adjoint_sensitivities, adjoint_sensitivities_u0, Sobol, Morris, gsa,
       SensitivityAlg, RegressionGSA, DGSM, eFAST

export BacksolveAdjoint, QuadratureAdjoint, InterpolatingAdjoint,
       TrackerAdjoint, ZygoteAdjoint, ReverseDiffAdjoint,
       ForwardSensitivity, ForwardDiffSensitivity,
       ForwardDiffOverAdjoint,
       SteadyStateAdjoint

export second_order_sensitivities, second_order_sensitivity_product

export TrackerVJP, ZygoteVJP, ReverseDiffVJP

export StochasticTransformedFunction
end # module


text:
# set current limit
# SI
export psp_setilimit

"""
    psp_setilimit(io, current_limit)

Set the current limit.
"""
function psp_setilimit(io_psp::IO,ilim::Real)
  @assert ilim<10 "current limit must be < 10A"
  ilimstring = @sprintf("%01.2f",ilim)
  set_current_limit_string = "SI "*ilimstring*"\r"
  write(io_psp, set_current_limit_string)
  return nothing
end


text:
#   This file is part of Reduce.jl. It is licensed under the MIT license
#   Copyright (C) 2017 Michael Reed

const switchbas = [
    :complex,
    :div,
    :lcm,
    :gcd,
]

const switches = [
    :factor,
    :ifactor,
    :expandlog,
    :combinelog,
    :precise,
    :combineexpt,
    :rounded,
    :evallhseq,
    :horner,
    :trigform,
    :fullroots,
    :cramer,
    :multiplicities,
    :allbranch,
    :arbvars,
    :solvesingular,
    :varopt,
    :intstr,
    :allfac,
    :rat,
    :revpri,
    :fort,
    :ratarg,
    :ezgcd,
    :heugcd,
    :bezout,
    :modular,
    :rational,
    :roundbf,
    :adjprec,
    :roundall,
    :balancedmod,
    :nocommutedf,
    :commutedf,
    :simpnoncomdf,
    :expanddf,
    :allowdfint,
    :dfint,
    :trint,
    :trintsubst,
    :failhard,
    :nolnr,
    :nointsubst,
]

push!(VERSION < v"1.0-" ? switchbas : switches, :expand)

const switchtex = [
    :nat,
    :latex
]

const onswitch = Symbol.("@",[switchbas;switches])
const offswitch = Symbol.("@off_",[switchbas;switches])

Expr(:toplevel,[:(import Base: $i) for i ∈ switchbas]...) |> eval
:(export $([switchbas;switches;switchtex]...)) |> eval
:(export $([onswitch;offswitch;Symbol.("@",switchtex)]...)) |> eval

for fun in [switchbas;switches;switchtex]
    @eval begin
        $(parsegen(fun,:switch))
        $fun(s::Bool) = (rcall("$(s ? "on" : "off") "*$(string(fun))); s)
    end
end

for fun in [switchbas;switches]
    @eval begin
        $(unfoldgen(fun,:switch))
        function $fun(expr::String;be=0)
            convert(String, $fun(RExpr(expr);be=be))
        end
    end
end

for fun in switchtex
    @eval begin
        $(unfoldgen(fun,:switch))
        function $fun(expr::String;be=0)
            convert(String, $fun(RExpr(expr);be=be))
        end
        macro $fun(expr)
            $fun(expr)
        end
    end
end

for fun in [switchbas;switches]
    @eval begin
        macro $fun(expr)
            (r,on,off) = macroshift(expr)
            push!(on,$(string(fun)))
            Expr(:quote,rcall(r;on=on,off=off))
        end
        macro $(Symbol("off_",fun))(expr)
            (r,on,off) = macroshift(expr)
            push!(off,$(string(fun)))
            Expr(:quote,rcall(r;on=on,off=off))
        end
    end
end

function macroshift(r)
    if typeof(r) == Expr && r.head == :macrocall
        (expr,on,off) = macroshift(r.args[3])
        if r.args[1] ∈ onswitch
            (expr,push!(on,string(r.args[1])[2:end]),off)
        elseif r.args[1] ∈ offswitch
            (expr,on,push!(off,string(r.args[1])[6:end]))
        else
            (eval(r),String[],String[])
        end
    else
        ((typeof(r)==Expr && r.head==:quote) ? r.args[1] : r,String[],String[])
    end
end

@doc """
    multiplicities(::Bool)

If you want the multiplicities explicitly displayed, the switch `multiplicities` can be turned on. For example
```Julia
julia> Algebra.on(:multiplicities); Algebra.solve(:(x^2==2x-1),:x)
```
yields the result
```Julia
(:(x = 1), :(x = 1))
```
""" Reduce.multiplicities


text:
using Documenter, ObjectFile

makedocs(modules = [ObjectFile],
         sitename = "ObjectFile.jl",
         pages = [
            "Home" => "index.md"
         ]
)

deploydocs(
    repo = "github.com/staticfloat/ObjectFile.jl.git",
    target = "build"
)


text:
@testset "Numpy" begin

    @testset "load" begin
        @test Pickle.npyload("./test_pkl/test-np.pkl") == pyload("./test_pkl/test-np.pkl")
        @test Pickle.npyload("./test_pkl/test-np-end.pkl") == pyload("./test_pkl/test-np-end.pkl")
        @test Pickle.npyload("./test_pkl/test-np-multi.pkl") == pyload("./test_pkl/test-np-multi.pkl")

        jlnpstrings = Pickle.npyload("./test_pkl/test-np-string.pkl")
        npstrings = pyload("./test_pkl/test-np-string.pkl")
        @test jlnpstrings["little-endian"] == convert(Array{String}, collect(npstrings["little-endian"]))
        @test jlnpstrings["big-endian"] == convert(Array{String}, collect(npstrings["big-endian"]))

        jlscalars = Pickle.npyload("./test_pkl/test-np-scalar.pkl")
        # [np.uint32(4), np.int64(6), np.float16(3), np.float64(5), np.unicode_("abcα😀")]
        npscalars = collect(pyload("./test_pkl/test-np-scalar.pkl"))
        @test jlscalars[1] == tryparse(UInt32, npscalars[1].__repr__())
        @test jlscalars[2] == tryparse(Int64, npscalars[2].__repr__())
        @test jlscalars[3] == tryparse(Float16, npscalars[3].__repr__())
        @test jlscalars[4] == Float64(npscalars[4])
        @test jlscalars[5] == String(npscalars[5])

    end

end