text:
# Calculates percent identity of two aligned sequences
# No account of the positions with gaps in both sequences in the length

"seq1 and seq2 should have the same len"
function _percentidentity(seq1, seq2, len)
    count = 0
    notcount = 0
    for i in 1:len
        @inbounds aa1 = seq1[i]
        @inbounds aa2 = seq2[i]
        if aa1 == aa2
            if aa1 == GAP || aa1 == XAA
                notcount += 1
            else
                count += 1
            end
        else
            if aa1 == XAA || aa2 == XAA
                notcount += 1
            end
        end
    end
    100.0 * count / (len - notcount)
end

"""
`percentidentity(seq1, seq2)`

Calculates the fraction of identities between two aligned sequences. The identity value is
calculated as the number of identical characters in the i-th position of both sequences
divided by the length of both sequences. Positions with gaps in both sequences doesn't
count to the length of the sequences. Positions with a `XAA` in at least one sequence
aren't counted.
"""
function percentidentity(seq1, seq2)
    len = length(seq1)
    if len != length(seq2)
        throw(ErrorException("""
        Sequences of different length, they aren't aligned or don't come from the same MSA.
        """))
    end
    _percentidentity(seq1, seq2, len)
end

"""
`percentidentity(seq1, seq2, threshold)`

Computes quickly if two aligned sequences have a identity value greater than a given
`threshold` value. Returns a boolean value. Positions with gaps in both sequences
doesn't count to the length of the sequences. Positions with a `XAA` in at least one
sequence aren't counted.
"""
function percentidentity(seq1, seq2, threshold)
    fraction = threshold / 100.0
    len = length(seq1)
    if len != length(seq2)
        throw(ErrorException("""
        Sequences of different length, they aren't aligned or don't come from the same MSA.
        """))
    end
    n = len
    limit_count = n * fraction
    diff = 0
    count = 0
    @inbounds for i in 1:len
        aa1 = seq1[i]
        aa2 = seq2[i]
        if aa1 == XAA || aa2 == XAA
            n -= 1
            continue
        end
        if aa1 == aa2
            if aa1 != GAP
                count += 1
                if count >= limit_count
                    return(true)
                end
            else
                n -= 1
                limit_count = n * fraction
            end
        else
            diff += 1
            if diff > n - limit_count
                return(false)
            end
        end
    end
    (count/n) >= fraction
end

# percentidentity for a MSA
# =========================

function _percentidentity_kernel!(scores, aln::Vector{Vector{Residue}}, nseq, len)
    k = 0
    list = getlist(scores)
    @inbounds for i in 1:(nseq-1)
        a = aln[i]
        for j in (i+1):nseq
            list[k += 1] = _percentidentity(a, aln[j], len)
        end
    end
    scores
end

"""
`percentidentity(msa[, out::Type=Float64])`

Calculates the identity between all the sequences on a MSA. You can indicate the output
element type with the last optional parameter (`Float64` by default). For a MSA with a lot
of sequences, you can use `Float32` or `Flot16` in order to avoid the `OutOfMemoryError()`.
"""
function percentidentity(msa::AbstractMatrix{Residue}, out::Type{T}=Float64) where T
    aln = getresiduesequences(msa)
    nseq = length(aln)
    len = length(aln[1])
    scores = sequencepairsmatrix(msa, T, Val{false}, T(100.0))
    _percentidentity_kernel!(getarray(scores), aln, nseq, len)
    scores
end

# Mean Percent Identity of an MSA
# ===============================

"""
Returns the mean of the percent identity between the sequences of a MSA. If the MSA has 300
sequences or less, the mean is exact. If the MSA has more sequences and the `exact` keyword
is `false` (defualt), 44850 random pairs of sequences are used for the estimation. The
number of samples can be changed using the second argument. Use `exact=true` to perform all
the pairwise comparison (the calculation could be slow).
"""
function meanpercentidentity(msa, nsamples::Int=44850; exact::Bool=false)
    #                 lengthlist(300, false) == 44850
    nseq, ncol = size(msa)
    @assert nsamples > 2 "At least 2 samples are needed."
    nvalues = lengthlist(nseq, Val{false})
    sum = 0.0
    if !exact && nvalues >= nsamples
        samples = Set{Tuple{Int,Int}}()
        sizehint!(samples, nsamples)
        while length(samples) < nsamples
            i = rand(1:(nseq-1))
            j = rand((i+1):nseq)
            if !in((i,j), samples)
                push!(samples, (i,j))
                sum += percentidentity(msa[i,:], msa[j,:])
            end
        end
        return( (sum/nsamples) )
    else # exact and/or few sequences
        for i in 1:(nseq-1)
            @inbounds @simd for j in (i+1):nseq
                sum += percentidentity(msa[i,:], msa[j,:])
            end
        end
        return( (sum/nvalues) )
    end
end

# Similarity percent
# ==================

"""
Calculates the similarity percent between two aligned sequences. The 100% is the length of
the aligned sequences minus the number of columns with gaps in both sequences and the number
of columns with at least one residue outside the alphabet. So, columns with residues outside
the alphabet (other than the specially treated `GAP`) aren't counted to the protein length.
Two residues are considered similar if they below to the same group in a `ReducedAlphabet`.
The `alphabet` (third positional argument) by default is:

```julia
ReducedAlphabet("(AILMV)(NQST)(RHK)(DE)(FWY)CGP")
```

The first group is composed of the non polar residues `(AILMV)`, the second group is composed
of polar residues, the third group are positive residues, the fourth group are negative
residues, the fifth group is composed by the aromatic residues `(FWY)`. `C`, `G` and `P`
are considered unique residues.

**Other residue groups/alphabets:**

**SMS (Sequence Manipulation Suite)** Ident and Sim:

```julia
ReducedAlphabet("(GAVLI)(FYW)(ST)(KRH)(DENQ)P(CM)")
```

*Stothard P (2000) The Sequence Manipulation Suite: JavaScript programs for analyzing and
formatting protein and DNA sequences. Biotechniques 28:1102-1104.*

**Bio3D 2.2** seqidentity:

```julia
ReducedAlphabet("(GA)(MVLI)(FYW)(ST)(KRH)(DE)(NQ)PC")
```

*Grant, B.J. et al. (2006) Bioinformatics 22, 2695--2696.*
"""
function percentsimilarity(seq1::Vector{Residue}, seq2::Vector{Residue},
    alphabet::ResidueAlphabet = ReducedAlphabet("(AILMV)(RHK)(NQST)(DE)(FWY)CGP"))

    len = length(seq1)
    if len != length(seq2)
        throw(ErrorException("""
        Sequences of different lengths, they aren't aligned or don't come from the same MSA.
        """))
    end

    count = 0
    colgap = 0
    colout = 0 # Columns with residues outside the alphabet aren't used
    @inbounds for i in 1:len
        res1 = seq1[i]
        res2 = seq2[i]
        isgap1 = res1 == GAP
        isgap2 = res2 == GAP
        if isgap1 && isgap2
            colgap += 1
            continue
        end
        if  (!isgap1 && !in(res1, alphabet)) ||
            (!isgap2 && !in(res2, alphabet))
            colout += 1
            continue
        end
        if !isgap1 && !isgap2
            count += Int(alphabet[res1] == alphabet[res2])
        end
    end
    alnlen = len - colgap - colout
    (100.0 * count) / alnlen
end

"""
Calculates the similarity percent between all the sequences on a MSA. You can indicate
the output element type with the `out` keyword argument (`Float64` by default). For an
MSA with a lot of sequences, you can use `out=Float32` or `out=Flot16` in order to
avoid the `OutOfMemoryError()`.
"""
function percentsimilarity(msa::AbstractMatrix{Residue}, A...; out::Type=Float64)
    M = getresiduesequences(msa)
    P = sequencepairsmatrix(msa, out, Val{false}, out(100.0))
    @inbounds @iterateupper getarray(P) false begin
        list[k] = percentsimilarity(M[i], M[j], A...)
    end
    P
end


text:
## to make seeds odd because PAML needs odd seeds
function makeOdd!(seeds::Vector{Int64})
    mod = seeds .% 2
    seeds[mod.==0]=seeds[mod.==0] .+ 1
end

## create the control file for PAML
## to simulate aminoacid sequences on a given tree
## note: tree is assumed to be a quartet
## name: name of the output control file (note that the control file
##       needs to be called MCaa.dat, but we will save also as name
##       to avoid overwriting
## tree: tree (quartet, rooted, with bl)
## s: random seed
## L: sequence length
## alpha: rate for gamma
## model: see notebook.md for model options for PAML
## modelfile: we chose one file from the ones provided by PAML with the S matrix
## n= number of leaves
function createCtlFile(name::String, tree::String, s::Integer, L::Integer, alpha::Float64,
                       model::Integer, modelfile::String; n=4::Integer)
    str = """0        * 0: paml format (mc.paml); 1:paup format (mc.nex)
$s       * random number seed (odd number)

$n $L 1   * <# seqs>  <# sites>  <# replicates>

-1         * <tree length, use -1 if tree below has absolute branch lengths>

$tree

$alpha 1        * <alpha; see notes below>  <#categories for discrete gamma>
$model $modelfile * <model> [aa substitution rate file, need only if model=2 or 3]

0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05
0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05

 A R N D C Q E G H I
 L K M F P S T W Y V

// end of file

=============================================================================
Notes for using the option in evolver to simulate amino acid sequences.
Change values of parameters, but do not delete them.  It is o.k. to add
empty lines, but do not break down the same line into two or more lines.

  model = 0 (poisson), 1 (proportional), 2 (empirical), 3 (empirical_F)
  Use 0 for alpha to have the same rate for all sites.
  Use 0 for <#categories for discrete gamma> to use the continuous gamma
  <aa substitution rate file> can be dayhoff.dat, jones.dat, and so on.
  <aa frequencies> have to be in the right order, as indicated.
=================!! Check screen output carefully!! =====================
"""
    outfile = open(name,"w")
    write(outfile, str)
    close(outfile)

    ## we need to name the file MCaa.dat
    outfile2 = open("MCaa.dat","w")
    write(outfile2, str)
    close(outfile2)
end

### simulate a quartet with branch lengths uniformly dist
### U(l,u)
## ss: seed
## lba=false by default: long branch attraction case, based on Figure 3 Zou2020:
## two short external branches have lengths b ranging from 0.1 to 1.0,
## the two long branches have lengths a ranging from 2b to 40b,
## and the internal branch has a length c ranging from 0.01b to b
## Root chosen randomly
## the random root was removed to implement the permutation
## strategy in Zou2019
function sampleRootedMetricQuartet(l::Number,u::Number, ss::Integer;
                                   lba=false::Bool, b=0.01::Float64, rab=2.0::Float64, rcb=0.01::Float64)
    Random.seed!(ss)
    quartets = ["((1,2),(3,4));", "((1,3),(2,4));", "((1,4),(2,3));"] ##only thinking of unrooted
    ind = sample(1:3,1)[1]
    q = quartets[ind]
    quartet = lba ? readTopologyLevel1(q) : readTopology(q)
    ## setting branch lenghts
    if lba
        c = rcb*b
        a = rab*b
        ## internal branch:
        intbl = findall([e.istIdentifiable for e in quartet.edge])[1]
        setLength!(quartet.edge[intbl],c)
        edges = [quartet.edge[intbl]]
        ## external edges:
        extbls = setdiff(1:length(quartet.edge),intbl)
        ## first external edge, set BL=b:
        setLength!(quartet.edge[extbls[1]],b)
        push!(edges, quartet.edge[extbls[1]])
        adje = adjacentedges(quartet.edge[extbls[1]])
        ## we find the sister external edge, set BL=a:
        for e in adje
            if e!= quartet.edge[extbls[1]] && !e.istIdentifiable
                e.length = a
                push!(edges,e)
            end
        end
        misse = setdiff(quartet.edge, edges)
        ## second external edge, set BL=b:
        setLength!(misse[1],b)
        misse2 = setdiff(misse, [misse[1]])
        ## sister set BL=a
        misse2[1].length = a
        ## rooting tree to match other simulations:
        rootonedge!(quartet,quartet.edge[intbl])
    else
        for e in quartet.edge
            setLength!(e,rand(Uniform(l,u),1)[1])
        end
    end
    return writeTopology(quartet),ind
end

## converts the file of sequences into one hot encoded tensor
## e.g. if file has 4 sequences of aminoacid (20 letters) of length L
## this function should return a tensor of 4x20xL where each letter is represented
## by its one-hot code: l1->000...01
## BUT in Zou2019 they then vectorize the tensor into matrix (e.g. 80xL),
## so we will return the matrix directly
## we need the length (L) to initialize the matrix
## note: it assumes phylip format so that the first row is "n L"
function convert2onehot(name::String, L::Integer)
    mc = readlines(name)
    mc = mc[mc .!= ""]
    mat = zeros(L,1) ##output matrix
    ## assumes phylip format so that first row is "n L"
    ## we do i=2 first to have the same alll vector
    i=2
    mm = split(mc[i])
    @show mm[1]
    seq = join(mm[2:end])
    alll = unique(seq)  ##e.g. the 20 aminoacids
    m = zeros(1,length(alll)) ## output matrix
    for s in seq
        m = [m;Flux.onehot(s,alll)']
    end
    m = m[2:end,:]
    mat = hcat(mat,m)

    for i in 3:length(mc)
        mm = split(mc[i])
        @show mm[1]
        seq = join(mm[2:end])
        m = zeros(1,length(alll)) ## output matrix
        for s in seq
            m = [m;Flux.onehot(s,alll)']
        end
        m = m[2:end,:]
        mat = hcat(mat,m)
    end
    mat = mat[:,2:end]
    return mat
end

### simulate a quintet with branch lengths uniformly dist
### U(l,u)
## ss: seed
### Rooted chosen randomly
## Note that we do not simulate the tree with coalescent (as rcoal)
## to make it easier to keep the "index" of the tree
function sampleRootedMetricQuintet(l::Number,u::Number, ss::Integer)
    Random.seed!(ss)
    perm = validPermutations()
    quintet = readTopology("(((t1,t2),(t3,t4)),t5);")
    ind = sample(1:length(perm),1)[1]
    renameTips!(quintet, string.(perm[ind]))
    ## choose root randomly (or leave as is=balanced tree)
    r = rand(Uniform(0,1),1)[1]
    if(r<0.2)
        rootatnode!(quintet,"1")
    elseif(r<0.4)
        rootatnode!(quintet,"2")
    elseif(r<0.6)
        rootatnode!(quintet,"3")
    elseif(r<0.8)
        rootatnode!(quintet,"4")
    end
    ## setting branch lenghts
    for e in quintet.edge
        setLength!(e,rand(Uniform(l,u),1)[1])
    end
    return writeTopology(quintet),ind
end

## rename one tip in HybridNetworks object
function renameTip!(net::HybridNetwork, oldtip::String, newtip::String)
    ind = findall([n.name for n in net.node] .== oldtip)
    length(ind) > 0 || error("$oldtip not found in network")
    length(ind) < 2 || error("$oldtip found in multiple leaves in network")
    net.node[ind[1]].name = newtip
end

## rename tips in HybridNetworks object
## oldtip: vector with old tips to change in the order of newtip vector
function renameTips!(net::HybridNetwork, oldtip::Vector{String}, newtip::Vector{String})
    for i in 1:length(oldtip)
        renameTip!(net, oldtip[i], newtip[i])
    end
end

## rename tips in HybridNetworks object
## in this version, you do not give the vector of old names
## you simply change the tip names in the order that the names appear in net
renameTips!(net::HybridNetwork, newtip::Vector{String}) = renameTips!(net, net.names, newtip)

## manually created list of valid permutations
## this only works for quintets
function validPermutations()
    perm = [
    [1,2,3,4,5],
[1,2,3,5,4],
[1,2,4,5,3],
[1,3,2,4,5],
[1,3,2,5,4],
[1,3,4,5,2],
[1,4,2,3,5],
[1,4,2,5,3],
[1,4,3,5,2],
[1,5,2,3,4],
[1,5,2,4,3],
[1,5,3,4,2],
[2,3,4,5,1],
[2,4,3,5,1],
        [2,5,3,4,1]]
    return perm
end

## from master PhyloNetworks but can't update the package now
function adjacentedges(centeredge::PhyloNetworks.Edge)
    n = centeredge.node
    length(n) == 2 || error("center edge is connected to $(length(n)) nodes")
    @inbounds edges = copy(n[1].edge) # shallow copy, to avoid modifying the first node
    @inbounds for ei in n[2].edge
        ei === centeredge && continue # don't add the center edge again
        # a second edge between nodes n[1] and n[2] would appear twice
        push!(edges, ei)
    end
    return edges
end

## function to read the simulated sequences from PAML
## and write without taxon names for the NN
## sequences saved to outfile: 4 x L
## append=true means that the file already exists and we append the
## sequences to the existing file
function writeSequence2File(name::String, L::Integer, outfile::String ; append=true::Bool)
    if append
        f = open(outfile, "a")
    else
        f = open(outfile, "w")
    end
    mc = readlines(name)
    mc = mc[mc .!= ""]
    ## assumes phylip format so that first row is "n L"
    ## we do i=2 first to have the same alll vector
    for i in 2:length(mc)
        mm = split(mc[i])
        @show mm[1]
        seq = join(mm[2:end])
        write(f,seq)
        write(f,"\n")
    end
    close(f)
end

ext: jl

max_stars_repo_name: Forgotten/nn-phylogenetics

max_issues_repo_name: Forgotten/nn-phylogenetics

max_forks_repo_name: Forgotten/nn-phylogenetics


text:
"""
TODO...
"""
function gibbs_pt(init::Vector{T}, args::Vector{Dict{Symbol, Any}},
                  update::Function;
                  monitors::Vector{Vector{Symbol}}=deepcopy(monitor_default),
                  thins::Vector{Int}=deepcopy(thin_default),
                  nmcmc::Int64=1000, nburn::Int=0,
                  printFreq::Int=0,
                  save_all_states::Bool=false,
                  printlnAfterMsg::Bool=true) where T

  @assert printFreq >= -1
  if printFreq == 0
    numPrints = 10
    printFreq = Int(ceil((nburn + nmcmc) / numPrints))
  end

  states = deepcopy(init)
  num_chains = length(states)

  # Checking number of monitors.
  numMonitors = length(monitors)
  println("Number of monitors: $(numMonitors)"); flush(stdout)
  @assert numMonitors == length(thins)

  # Check monitor
  if numMonitors == 0
    println("Using default monitor."); flush(stdout)
    fnames = [fname for fname in fieldnames(typeof(init))]
    append!(monitors, [fnames])
    append!(thins, 1)
    numMonitors = 1
  end

  # Number of Samples for each Monitor
  numSamps = [div(nmcmc, thins[i]) for i in 1:numMonitors]

  if printFreq > 0
    println("Preallocating memory..."); flush(stdout)
  end

  # Create object to return
  if save_all_states
    @time out = [fill([deepcopyFields(s, monitors[i]) for s in states],
                      numSamps[i]) for i in 1:numMonitors]
  else
    @time out = [fill(deepcopyFields(states[1], monitors[i]), numSamps[i])
                 for i in 1:numMonitors]
  end

  function printMsg(i::Int)
    if (printFreq > 0) && (i % printFreq == 0) && (i > 1)
      if :ll in keys(args[1])
        loglikeMsg = "-- loglike: $(last(args[1][:ll]))"
      else
        loglikeMsg = ""
      end

      print("$(showtime()) -- $(i)/$(nburn+nmcmc) $loglikeMsg")

      if printlnAfterMsg
        println()
      end

      flush(stdout)
    end
  end

  # burn in
  for i in 1:nburn
    states, args = update(states, args, i)
    printMsg(i)
  end

  counters = zeros(Int, numMonitors)
  # Gibbs loop
  for i in 1:nmcmc
    states, args = update(states, args, i + nburn)
    printMsg(i + nburn)

    for j in 1:numMonitors
      if i % thins[j] == 0
        if save_all_states
          substate = [deepcopyFields(s, monitors[j]) for s in states]
        else
          substate = deepcopyFields(states[1], monitors[j])
        end
        counters[j] += 1
        out[j][counters[j]] = substate
      end
    end
  end

  return (out, states, args)
end


text:
#################### Approximate Bayesian Computation ####################

#################### Types and Constructors ####################

mutable struct ABCTune
  datakeys::Vector{Symbol}
  Tsim::Vector{Vector{Float64}}
  epsilon::Vector{Float64}
  epsilonprime::Vector{Float64}

  function ABCTune()
    new(
      Vector{Symbol}(),
      Vector{Vector{Float64}}(),
      Vector{Float64}()
    )
  end
end

#################### Sampler Constructor ####################
"""
              ABC(params::ElementOrVector{Symbol},

              scale::ElementOrVector{T}, summary::Function,

              epsilon::Real; kernel::KernelDensityType=SymUniform,

              dist::Function=(Tsim, Tobs) -> sqrt(sum(abs2, Tsim - Tobs)),

              proposal::SymDistributionType=Normal, maxdraw::Integer=1,

              nsim::Integer=1, decay::Real=1.0, randeps::Bool=false,

              args...) where {T<:Real}

Construct a `Sampler` object for ABC sampling. Parameters are assumed to be continuous, but may be constrained or unconstrained.

Returns a `Sampler{ABCTune}` type object.

* `params`: stochastic node(s) to be updated with the sampler.  Constrained parameters are mapped to unconstrained space according to transformations defined by the Stochastic `unlist()` function.

* `scale` : scaling value or vector of the same length as the combined elements of nodes `params` for the `proposal` distribution.  Values are relative to the unconstrained parameter space, where candidate draws are generated.

* `summary` : function that takes a vector of observed or simulated data and returns a summary statistic or vector of statistics.

* `epsilon` : target tolerance for determining how similar observed and simulated data summary statistics need to be in order to accept a candidate draw.

* `kernel` : weighting kernel density of type `Biweight`, `Cosine`, `Epanechnikov`, `Normal`, `SymTriangularDist`, `SymUniform`, or `Triweight` to use in measuring similarity between observed and simulated data summary statistics.  Specified `epsilon` determines the standard deviation of Normal kernels and widths of the others.

* `dist` : positive function for the kernel density to compute distance between vectors of observed (`Tobs`) and simulated (`Tsim`) data summary statistics (default: Euclidean distance).

* `proposal` : symmetric distribution of type `Biweight`, `Cosine`, `Epanechnikov`, `Normal`, `SymTriangularDist`, `SymUniform`, or `Triweight` to be centered around current parameter values and used to generate proposal draws.  Specified `scale` determines the standard deviations of Normal proposals and widths of the others.

* `maxdraw` : maximum number of unaccepted candidates to draw in each call of the sampler.  Draws are generated until one is accepted or the maximum is reached.  Larger values increase acceptance rates at the expense of longer runtimes.

* `nsim` : number of data sets to simulate in deciding whether to accept a candidate draw.  Larger values lead to closer approximations of the target distribution at the expense of longer runtimes.

* `decay` : if `0 < decay <= 1`, the rate at which internal tolerances are monotonically decreased from the initial distance between observed and simulated summary statistics toward the maximum of each subsequent distance and `epsilon`; if `decay = 0`, internal tolerances are fixed at `epsilon`.

* `randeps` : whether to perturb internal tolerances by random exponential variates.

* `args...` : additional keyword arguments to be passed to the `dist` function.

"""
function ABC(params::ElementOrVector{Symbol},
              scale::ElementOrVector{T}, summary::Function,
              epsilon::Real; kernel::KernelDensityType=SymUniform,
              dist::Function=(Tsim, Tobs) -> sqrt(sum(abs2, Tsim - Tobs)),
              proposal::SymDistributionType=Normal, maxdraw::Integer=1,
              nsim::Integer=1, decay::Real=1.0, randeps::Bool=false,
              args...) where {T<:Real}
  0 <= decay <= 1 || throw(ArgumentError("decay is not in [0, 1]"))

  params = asvec(params)
  kernelpdf = (epsilon, d) -> pdf(kernel(0.0, epsilon), d)
  local obsdata, simdata, summarizenodes

  samplerfx = function(model::Model, block::Integer)
    tune = gettune(model, block)

    ## current parameter and density values
    theta0 = unlist(model, block, true)
    logprior0 = logpdf(model, params, true)
    pi_epsilon0 = 0.0
    pi_error0 = 1.0

    ## initialize tuning parameters
    local Tobs
    if model.iter == 1
      ## data node symbols
      targets = keys(model, :target, params)
      stochastics = keys(model, :stochastic)
      tune.datakeys = intersect(setdiff(targets, params), stochastics)

      ## local functions to get and summarize data nodes
      obsdata = key -> unlist(model[key])
      simdata = key -> unlist(model[key], rand(model[key]))
      summarizenodes = length(tune.datakeys) > 1 ?
        data -> vcat(map(key -> summary(data(key)), tune.datakeys)...) :
        data -> asvec(summary(data(tune.datakeys[1])))

      ## observed data summary statistics
      Tobs = summarizenodes(obsdata)

      tune.Tsim = Array{Vector{Float64}}(undef, nsim)
      tune.epsilon = Array{Float64}(undef, nsim)
      tune.epsilonprime = Array{Float64}(undef, nsim)
      for i in 1:nsim
        ## simulated data summary statistics for current parameter values
        tune.Tsim[i] = summarizenodes(simdata)
        d = dist(tune.Tsim[i], Tobs; args...)

        ## starting tolerance
        tune.epsilon[i] = decay > 0 ? max(epsilon, d) : epsilon
        if randeps
          dexp = Exponential(tune.epsilon[i])
          tune.epsilonprime[i] = rand(dexp)
          pi_error0 = pdf(dexp, tune.epsilonprime[i])
        else
          tune.epsilonprime[i] = tune.epsilon[i]
        end

        ## kernel density evaluation
        pi_epsilon0 += kernelpdf(tune.epsilonprime[i], d) * pi_error0
      end
    else
      Tobs = summarizenodes(obsdata)
      for i in 1:nsim
        d = dist(tune.Tsim[i], Tobs; args...)
        if randeps
          dexp = Exponential(tune.epsilon[i])
          pi_error0 = pdf(dexp, tune.epsilonprime[i])
        end
        pi_epsilon0 += kernelpdf(tune.epsilonprime[i], d) * pi_error0
      end
    end

    Tsim1 = similar(tune.Tsim)
    epsilon1 = similar(tune.epsilon)
    epsilonprime1 = similar(tune.epsilonprime)

    for k in 1:maxdraw
      ## candidate draw and prior density value
      theta1 = theta0 + scale .* rand(proposal(0.0, 1.0), length(theta0))
      logprior1 = mapreduce(keyval -> logpdf(model[keyval[1]], keyval[2]), +,
                            relist(model, theta1, block, true))
      if logprior1 == -Inf
        continue
      end
      relist!(model, theta1, block, true)

      ## tolerances and kernel density
      pi_epsilon1 = 0.0
      pi_error1 = 1.0
      for i in 1:nsim
        ## simulated data summary statistics for candidate draw
        Tsim1[i] = summarizenodes(simdata)
        d = dist(Tsim1[i], Tobs; args...)

        ## monotonically decrease tolerance to target
        epsilon1[i] = (1 - decay) * tune.epsilon[i] +
                      decay * max(epsilon, min(d, tune.epsilon[i]))
        if randeps
          dexp = Exponential(epsilon1[i])
          epsilonprime1[i] = rand(dexp)
          pi_error1 = pdf(dexp, epsilonprime1[i])
        else
          epsilonprime1[i] = epsilon1[i]
        end

        ## kernel density evaluation
        pi_epsilon1 += kernelpdf(epsilonprime1[i], d) * pi_error1
      end

      ## accept/reject the candidate draw
      if rand() < pi_epsilon1 / pi_epsilon0 * exp(logprior1 - logprior0)
        theta0 = theta1
        tune.Tsim = Tsim1
        tune.epsilon = epsilon1
        tune.epsilonprime = epsilonprime1
        break
      end
    end

    relist(model, theta0, block, true)
  end

  Sampler(params, samplerfx, ABCTune())
end


text:
using Plots
using Random
using DataFrames
using AnnotationFusion

Random.seed!(4)

n = 100
scores = rand(1:9, n)

annotations = DataFrame(items = 1:n, A = scores, B = round.(Int, clamp.(1.2 .* scores, 1, 9)))

μ = fuse(annotations, :items, Mean())

te = fuse(annotations, :items, TE())

copeland = fuse(annotations, :items, Copeland())

plot(annotations.A, label="Annotator A", xlabel="Item index", ylabel="Score", title="Scores")
plot!(annotations.B, label="Annotator B")
plot!(μ, label="Mean")
plot!(te, label="Triplet Embedding")
plot!(copeland, label="Copeland")


text:
"""

```
readPDB(filename; selection)
```

Reads a PDB file and stores the data in a vector of type `Atom`.

### Example

```julia-repl
julia> protein = readPDB("../test/structure.pdb")
   Array{Atoms,1} with 62026 atoms with fields:
   index name resname chain   resnum  residue        x        y        z  beta occup model segname index_pdb
       1    N     ALA     A        1        1   -9.229  -14.861   -5.481  0.00  1.00     1    PROT         1
       2  HT1     ALA     A        1        1  -10.048  -15.427   -5.569  0.00  0.00     1    PROT         2
       3  HT2     ALA     A        1        1   -9.488  -13.913   -5.295  0.00  0.00     1    PROT         3
                                                       ⋮
   62024  OH2    TIP3     C     9339    19638   13.485   -4.534  -34.438  0.00  1.00     1    WAT2     62024
   62025   H1    TIP3     C     9339    19638   13.218   -3.647  -34.453  0.00  1.00     1    WAT2     62025
   62026   H2    TIP3     C     9339    19638   12.618   -4.977  -34.303  0.00  1.00     1    WAT2     62026

julia> protein = readPDB("../test/structure.pdb","resname = ALA")
   Array{Atoms,1} with 72 atoms with fields:
   index name resname chain   resnum  residue        x        y        z  beta occup model segname index_pdb
       1    N     ALA     A        1        1   -9.229  -14.861   -5.481  0.00  1.00     1    PROT         1
       2  HT1     ALA     A        1        1  -10.048  -15.427   -5.569  0.00  0.00     1    PROT         2
       3  HT2     ALA     A        1        1   -9.488  -13.913   -5.295  0.00  0.00     1    PROT         3
                                                       ⋮
    1338  HB3     ALA     A       95       95   11.464   -3.278   -4.953  0.00  1.00     1    PROT      1338
    1339    C     ALA     A       95       95   14.815   -3.057   -5.633  0.00  1.00     1    PROT      1339
    1340    O     ALA     A       95       95   14.862   -2.204   -6.518  0.00  1.00     1    PROT      1340
```

"""
function readPDB(file::String, selection::String)
  query = parse_query(selection)
  return readPDB(file, only = atom -> apply_query(query,atom) )
end

function readPDB(file::String; only = all)

  # Check if structure is in mmCIF format
  mmCIF, mmCIF_fields = check_mmCIF(file)

  # Read file
  pdbfile = open(file,"r")
  natoms = 0
  index = 0
  imodel = 1
  iresidue = 1
  atoms = Atom[]
  local lastatom
  for line in eachline(pdbfile)
    if occursin("END",line)
      imodel = imodel + 1
    end
    atom = read_atom(line, mmCIF = mmCIF, mmCIF_fields = mmCIF_fields)
    if atom != nothing
      index = index + 1
      atom.index = index
      atom.model = imodel
      if index > 1
        if ! same_residue(atom,lastatom)
          iresidue += 1
        end
      end
      atom.residue = iresidue
      if only(atom)
        natoms = natoms + 1
        push!(atoms,atom)
      end
      lastatom = atom
    end
  end
  close(pdbfile)
  if natoms == 0
    error(" Could not find any atom in PDB file matching the selection. ")
  end

  return atoms
end

function Base.show(io::IO, atoms::AbstractVector{Atom})
  println(" Structure file with ", length(atoms), " atoms. ")
end


text:
# length(pdp::VkPhysicalDeviceProperties) = 1
# iterate(pdp::VkPhysicalDeviceProperties) = (pdp, nothing)


text:
using BasicBSpline
using Plots
gr()

##
p = 2
k = Knots(1:8)
P = BSplineSpace(p, k)
plot([t -> bsplinebasis₊₀(P, i, t) for i in 1:dim(P)], 1, 8, ylims = (0, 1))
savefig("docs/src/img/bsplinebasisplot.png")

##
p = 2
k = Knots(1:8)
P = BSplineSpace(p, k)
plot([t -> bsplinebasis′₊₀(P, i, t) for i in 1:dim(P)], 1, 8)
savefig("docs/src/img/bsplinebasisderivativeplot.png")

##
p = 2
k = Knots(1:8)
P = BSplineSpace(p, k)
plot(t -> sum(bsplinebasis₊₀(P, i, t) for i in 1:dim(P)), 1, 8, ylims = (0, 1))
savefig("docs/src/img/sumofbsplineplot.png")

##
p = 2
k = Knots(1:8) + p * Knots([1, 8])
P = BSplineSpace(p, k)
plot(t -> sum(bsplinebasis₊₀(P, i, t) for i in 1:dim(P)), 1, 8, ylims = (0, 1))
savefig("docs/src/img/sumofbsplineplot2.png")

##
p = 2
k = Knots(1:8) + p * Knots([1, 8])
P = BSplineSpace(p, k)
plot(t -> sum(bsplinebasis(P, i, t) for i in 1:dim(P)), 1, 8, ylims = (0, 1))
savefig("docs/src/img/sumofbsplineplot3.png")

##
P1 = BSplineSpace(1, Knots([1, 3, 5, 8]))
P2 = BSplineSpace(1, Knots([1, 3, 5, 6, 8, 9]))
P3 = BSplineSpace(2, Knots([1, 1, 3, 3, 5, 5, 8, 8]))
plot([t -> bsplinebasis₊₀(P1, i, t) for i in 1:dim(P1)], 1, 9, ylims = (0, 1))
savefig("docs/src/img/bsplineplotP1.png")
plot([t -> bsplinebasis₊₀(P2, i, t) for i in 1:dim(P2)], 1, 9, ylims = (0, 1))
savefig("docs/src/img/bsplineplotP2.png")
plot([t -> bsplinebasis₊₀(P3, i, t) for i in 1:dim(P3)], 1, 9, ylims = (0, 1))
savefig("docs/src/img/bsplineplotP3.png")

plot(
    plot([t -> bsplinebasis₊₀(P1, i, t) for i in 1:dim(P1)], 1, 9, ylims = (0, 1), legend = false),
    plot([t -> bsplinebasis₊₀(P2, i, t) for i in 1:dim(P2)], 1, 9, ylims = (0, 1), legend = false),
    plot([t -> bsplinebasis₊₀(P3, i, t) for i in 1:dim(P3)], 1, 9, ylims = (0, 1), legend = false),
    layout = (3, 1),
    link = :x,
)
savefig("docs/src/img/subbsplineplot.png")

##
A12 = changebasis(P1, P2)
A13 = changebasis(P1, P3)

plot(
    plot([t -> bsplinebasis₊₀(P1, i, t) for i in 1:dim(P1)], 1, 9, ylims = (0, 1), legend = false),
    plot([t -> sum(A12[i, j] * bsplinebasis₊₀(P2, j, t) for j in 1:dim(P2)) for i in 1:dim(P1)], 1, 9, ylims = (0, 1), legend = false),
    plot([t -> sum(A13[i, j] * bsplinebasis₊₀(P3, j, t) for j in 1:dim(P3)) for i in 1:dim(P1)], 1, 9, ylims = (0, 1), legend = false),
    layout = (3, 1),
    link = :x,
)
savefig("docs/src/img/subbsplineplot2.png")

##
p = 3
k = Knots([0.00, 1.50, 2.50, 5.50, 8.00, 9.00, 9.50, 10.0])
P0 = FastBSplineSpace(0, k)
P1 = FastBSplineSpace(1, k)
P2 = FastBSplineSpace(2, k)
P3 = FastBSplineSpace(3, k)
plot(
    plot([t -> bsplinebasis(P0, i, t) for i in 1:dim(P0)], 0, 10, ylims = (0, 1), legend = false),
    plot([t -> bsplinebasis(P1, i, t) for i in 1:dim(P1)], 0, 10, ylims = (0, 1), legend = false),
    plot([t -> bsplinebasis(P2, i, t) for i in 1:dim(P2)], 0, 10, ylims = (0, 1), legend = false),
    plot([t -> bsplinebasis(P3, i, t) for i in 1:dim(P3)], 0, 10, ylims = (0, 1), legend = false),
    layout = (2, 2),
)
savefig("docs/src/img/cover.png")


text:
export fitness_evaluate, select_indexes

using Random

"""
    null_evaluate(i::CGPInd, j::CGPInd)

Default evaluation function for two CGP individuals setting minimum fitness.
"""
function null_evaluate(i::CGPInd, j::CGPInd)
    return -Inf
end

"""
	fitness_evaluate(e::NSGA2Evo, fitness::Function)

Evaluation function for the NSGA-II algorithm.
"""
function fitness_evaluate(e::NSGA2Evo{T}, fitness::Function) where T
	@sync @inbounds for i in eachindex(e.population)
		Threads.@spawn begin
			e.population[i].fitness .= fitness(e.population[i], e.gen, e.atari_games[i])
		end
    end
end

function fitness_evaluate_two_atari_runs(e::NSGA2Evo{T}, fitness::Function) where T
	n = length(e.population)
	o1, f1 = zeros(n), zeros(Int64, n)
	o2, f2 = zeros(n), zeros(Int64, n)
	@sync @inbounds for i in eachindex(e.population)
		Threads.@spawn begin
			o1[i], f1[i] = fitness(e.population[i], e.gen, e.atari_games[i][1])
		end
		Threads.@spawn begin
			o2[i], f2[i] = fitness(e.population[i], 25000+e.gen, e.atari_games[i][2])
		end
    end
	@inbounds for i in eachindex(e.population)
		e.population[i].fitness .= [o1[i], o2[i]]
		e.population[i].reached_frames = f1[i] + f2[i]
	end
end

"""
    fitness_evaluate(e::CGPEvolution; fitness::Function=null_evaluate)

Sets the fitness of each individual to the Array of values returned by fitness.
Multithreaded option enabled in this version.
"""
function fitness_evaluate(e::CGPEvolution, fitness::Function)
    @sync for i in eachindex(e.population)
        Threads.@spawn begin
            e.population[i].fitness .= fitness(e.population[i])
        end
    end
end

function fitness_evaluate_ij(
	e::DualCGPEvolution,
	f::Array{Float64, 2},
	i::Int64,
	j::Int64,
	fitness::Function=null_evaluate
)
	enco_i = IPCGPInd(e.encoder_config, e.encoder_population[i].chromosome)
	cont_j = CGPInd(e.controller_config, e.controller_population[j].chromosome)
	f[i,j] = fitness(enco_i, cont_j)[1] # Currently, only pick 1st fitness dimension
end

"""
    fitness_evaluate(e::DualCGPEvolution, fitness::Function=null_evaluate)

Sets the fitness of each individual to the maximum value of the fitness matrix
in the dual CGP evolution framework.
"""
function fitness_evaluate(e::DualCGPEvolution, fitness::Function=null_evaluate)
    n_enco = e.encoder_config.n_population
    n_cont = e.controller_config.n_population
    fitness_matrix = zeros(n_enco, n_cont)

	# FORMER method with @sync
    #=
	@sync for i in 1:n_enco
        for j in 1:n_cont
            encoder_i = IPCGPInd(e.encoder_config, e.encoder_population[i].chromosome)
            controller_j = CGPInd(e.controller_config, e.controller_population[j].chromosome)
            Threads.@spawn begin
                fitness_matrix[i, j] = fitness(encoder_i, controller_j)[1] # Currently, only pick 1st fitness dimension
            end
        end
    end
	=#

	# NEW method with Threads.@threads
	indexes = [(i, j) for i in 1:n_enco for j in 1:n_cont]
    Threads.@threads for l in 1:(n_enco*n_cont)
        i, j = indexes[l]
		fitness_evaluate_ij(e, fitness_matrix, i, j, fitness)
    end

    for i in eachindex(e.encoder_population)
        e.encoder_population[i].fitness[1] = maximum(fitness_matrix[i,:])
    end
    for j in eachindex(e.controller_population)
        e.controller_population[j].fitness[1] = maximum(fitness_matrix[:,j])
    end
end

"""
	select_indexes(e::DualCGPGAEvo)

Select the indexes of the pairs to be evaluated for this generation.
	- 1. Select the elite pairs
	- 2. Select at least one pair per row/col
	- 3. Select random pairs until e.n_eval pairs have been selected
"""
function select_indexes(e::DualCGPGAEvo)
	indexes = Vector{Tuple{Int64, Int64}}()
	nrows = size(e.elites_matrix, 1)
	ncols = size(e.elites_matrix, 2)
	nmax = max(nrows, ncols)
	# 1. Select elites and maybe other pairs set in eval_matrix (e.g. mutant)
	for i in 1:nrows
		for j in 1:ncols
			if e.eval_matrix[i, j]
				push!(indexes, (i, j))
			end
		end
	end
	# 2. Select at least one pair per row/col
	shuffledrows = shuffle(collect(1:nrows))
	shuffledcols = shuffle(collect(1:ncols))
	candidates = Vector{Tuple{Int64, Int64}}()
	for i in 1:nmax
		i_modrow = i-nrows*divrem(i-1,nrows)[1]
		i_modcol = i-ncols*divrem(i-1,ncols)[1]
		c = (shuffledrows[i_modrow], shuffledcols[i_modcol])
		push!(candidates, c)
	end
	eliterows = [i[1] for i in indexes]
	elitecols = [i[2] for i in indexes]
	for c in candidates # Only push candidate if the row or col is not evaluated
		if c[1] ∉ eliterows || c[2] ∉ elitecols
			push!(indexes, c)
		end
	end
	# 3. Select additional random evaluations
	while length(indexes) < e.n_eval
		c = (rand(1:nrows), rand(1:ncols))
		if c ∉ indexes
			push!(indexes, c)
		end
	end
	indexes
end

function fitness_evaluate_ij!(
	e::DualCGPGAEvo,
	i::Int64,
	j::Int64,
	fitness::Function
)
	enco_i = IPCGPInd(e.encoder_config, e.encoder_sympop[i].chromosome)
	cont_j = CGPInd(e.controller_config, e.controller_sympop[j].chromosome)
	e.fitness_matrix[i, j] = fitness(enco_i, cont_j, e.gen)[1]
end

function set_ind_fitnesses!(e::DualCGPGAEvo)
	for i in eachindex(e.encoder_sympop)
        e.encoder_sympop[i].fitness[1] = maximum(e.fitness_matrix[i,:])
    end
    for j in eachindex(e.controller_sympop)
        e.controller_sympop[j].fitness[1] = maximum(e.fitness_matrix[:,j])
    end
end

function set_elites!(e::DualCGPGAEvo)
	e.elites_matrix = falses(size(e.elites_matrix)...) # zero elites matrix
	ci = CartesianIndices(size(e.fitness_matrix))
    # p = sortperm(vec(e.fitness_matrix))[end-e.n_elite+1:end]
	p = partialsortperm(vec(e.fitness_matrix), 1:e.n_elite; rev=true)
    elite_indexes = ci[p]
	for i in elite_indexes
		e.elites_matrix[i] = true
	end
end

"""
	fitness_evaluate(e::DualCGPGAEvo, fitness::Function=null_evaluate)

GA sparse fitness evaluation method.
"""
function fitness_evaluate(e::DualCGPGAEvo, fitness::Function=null_evaluate)
	# 1. Select indexes of individuals to evaluate
	indexes = select_indexes(e)
	# 2. Evaluate those individuals
	Threads.@threads for l in eachindex(indexes)
        i, j = indexes[l]
		fitness_evaluate_ij!(e, i, j, fitness)
    end
	# 3. Set individual's fitnesses
	set_ind_fitnesses!(e)
	# 4. Set elites
	set_elites!(e)
end


text:
module StaticLint

include("exception_types.jl")

using SymbolServer, CSTParser

using CSTParser: EXPR, isidentifier, setparent!, valof, headof, hastrivia, parentof, isoperator, ispunctuation, to_codeobject
# CST utils
using CSTParser: is_getfield, isassignment, isdeclaration, isbracketed, iskwarg, iscall, iscurly, isunarycall, isunarysyntax, isbinarycall, isbinarysyntax, issplat, defines_function, is_getfield_w_quotenode, iswhere, iskeyword, isstringliteral, isparameters, isnonstdid, istuple
using SymbolServer: VarRef

const noname = EXPR(:noname, nothing, nothing, 0, 0, nothing, nothing, nothing)

include("coretypes.jl")
include("bindings.jl")
include("scope.jl")
include("subtypes.jl")
include("methodmatching.jl")

mutable struct Meta
    binding::Union{Nothing,Binding}
    scope::Union{Nothing,Scope}
    ref::Union{Nothing,Binding,SymbolServer.SymStore}
    error
end
Meta() = Meta(nothing, nothing, nothing, nothing)

function Base.show(io::IO, m::Meta)
    m.binding !== nothing && show(io, m.binding)
    m.ref !== nothing && printstyled(io, " * ", color = :red)
    m.scope !== nothing && printstyled(io, " new scope", color = :green)
    m.error !== nothing && printstyled(io, " lint ", color = :red)
end
hasmeta(x::EXPR) = x.meta isa Meta
hasbinding(m::Meta) = m.binding isa Binding
hasref(m::Meta) = m.ref !== nothing
hasscope(m::Meta) = m.scope isa Scope
scopeof(m::Meta) = m.scope
bindingof(m::Meta) = m.binding

"""
    ExternalEnv

Holds a representation of an environment cached by SymbolServer.
"""
mutable struct ExternalEnv
    symbols::SymbolServer.EnvStore
    extended_methods::Dict{SymbolServer.VarRef,Vector{SymbolServer.VarRef}}
    project_deps::Vector{Symbol}
end

abstract type State end
mutable struct Toplevel{T} <: State
    file::T
    included_files::Vector{String}
    scope::Scope
    in_modified_expr::Bool
    modified_exprs::Union{Nothing,Vector{EXPR}}
    delayed::Vector{EXPR}
    resolveonly::Vector{EXPR}
    env::ExternalEnv
    server
    flags::Int
end

Toplevel(file, included_files, scope, in_modified_expr, modified_exprs, delayed, resolveonly, env, server) =
    Toplevel(file, included_files, scope, in_modified_expr, modified_exprs, delayed, resolveonly, env, server, 0)

function (state::Toplevel)(x::EXPR)
    resolve_import(x, state)
    mark_bindings!(x, state)
    add_binding(x, state)
    mark_globals(x, state)
    handle_macro(x, state)
    s0 = scopes(x, state)
    resolve_ref(x, state)
    followinclude(x, state)

    old_in_modified_expr = state.in_modified_expr
    if state.modified_exprs !== nothing && x in state.modified_exprs
        state.in_modified_expr = true
    end
    if CSTParser.defines_function(x) || CSTParser.defines_macro(x) || headof(x) === :export
        if state.in_modified_expr
            push!(state.delayed, x)
        else
            push!(state.resolveonly, x)
        end
    else
        old = flag!(state, x)
        traverse(x, state)
        state.flags = old
    end

    state.in_modified_expr = old_in_modified_expr
    state.scope != s0 && (state.scope = s0)
    return state.scope
end

mutable struct Delayed <: State
    scope::Scope
    env::ExternalEnv
    server
    flags::Int
end

Delayed(scope, env, server) = Delayed(scope, env, server, 0)

function (state::Delayed)(x::EXPR)
    mark_bindings!(x, state)
    add_binding(x, state)
    mark_globals(x, state)
    handle_macro(x, state)
    s0 = scopes(x, state)
    resolve_ref(x, state)

    old = flag!(state, x)
    traverse(x, state)
    state.flags = old
    if state.scope != s0
        for b in values(state.scope.names)
            infer_type_by_use(b, state.env)
            check_unused_binding(b, state.scope)
        end
        state.scope = s0
    end
    return state.scope
end

mutable struct ResolveOnly <: State
    scope::Scope
    env::ExternalEnv
    server
end

function (state::ResolveOnly)(x::EXPR)
    if hasscope(x)
        s0 = state.scope
        state.scope = scopeof(x)
    else
        s0 = state.scope
    end
    resolve_ref(x, state)

    traverse(x, state)
    if state.scope != s0
        state.scope = s0
    end
    return state.scope
end

# feature flags that can disable or enable functionality further down in the CST
const NO_NEW_BINDINGS = 0x1

function flag!(state, x::EXPR)
    old = state.flags
    if CSTParser.ismacrocall(x) && (valof(x.args[1]) == "@." || valof(x.args[1]) == "@__dot__")
        state.flags |= NO_NEW_BINDINGS
    end
    return old
end

"""
    semantic_pass(file, modified_expr=nothing)

Performs a semantic pass across a project from the entry point `file`. A first pass traverses the top-level scope after which secondary passes handle delayed scopes (e.g. functions). These secondary passes can be, optionally, very light and only seek to resovle references (e.g. link symbols to bindings). This can be done by supplying a list of expressions on which the full secondary pass should be made (`modified_expr`), all others will receive the light-touch version.
"""
function semantic_pass(file, modified_expr = nothing)
    server = file.server
    env = getenv(file, server)
    setscope!(getcst(file), Scope(nothing, getcst(file), Dict(), Dict{Symbol,Any}(:Base => env.symbols[:Base], :Core => env.symbols[:Core]), nothing))
    state = Toplevel(file, [getpath(file)], scopeof(getcst(file)), modified_expr === nothing, modified_expr, EXPR[], EXPR[], env, server)
    state(getcst(file))
    for x in state.delayed
        if hasscope(x)
            traverse(x, Delayed(scopeof(x), env, server))
            for (k, b) in scopeof(x).names
                infer_type_by_use(b, env)
                check_unused_binding(b, scopeof(x))
            end
        else
            traverse(x, Delayed(retrieve_delayed_scope(x), env, server))
        end
    end
    if state.resolveonly !== nothing
        for x in state.resolveonly
            if hasscope(x)
                traverse(x, ResolveOnly(scopeof(x), env, server))
            else
                traverse(x, ResolveOnly(retrieve_delayed_scope(x), env, server))
            end
        end
    end
end

"""
    traverse(x, state)

Iterates across the child nodes of an EXPR in execution order (rather than
storage order) calling `state` on each node.
"""
function traverse(x::EXPR, state)
    if (isassignment(x) && !(CSTParser.is_func_call(x.args[1]) || CSTParser.iscurly(x.args[1]))) || CSTParser.isdeclaration(x)
        state(x.args[2])
        state(x.args[1])
    elseif CSTParser.iswhere(x)
        for i = 2:length(x.args)
            state(x.args[i])
        end
        state(x.args[1])
    elseif headof(x) === :generator || headof(x) === :filter
        @inbounds for i = 2:length(x.args)
            state(x.args[i])
        end
        state(x.args[1])
    elseif headof(x) === :call && length(x.args) > 1 && headof(x.args[2]) === :parameters
        state(x.args[1])
        @inbounds for i = 3:length(x.args)
            state(x.args[i])
        end
        state(x.args[2])
    elseif x.args !== nothing && length(x.args) > 0
        @inbounds for i = 1:length(x.args)
            state(x.args[i])
        end
    end
end

"""
    followinclude(x, state)

Checks whether the arguments of a call to `include` can be resolved to a path.
If successful it checks whether a file with that path is loaded on the server
or a file exists on the disc that can be loaded.
If this is successful it traverses the code associated with the loaded file.

"""
function followinclude(x, state::State)
    if CSTParser.iscall(x) && length(x.args) > 0 && isidentifier(x.args[1]) && valofid(x.args[1]) == "include"

        init_path = path = get_path(x, state)
        if isempty(path)
        elseif isabspath(path)
            if hasfile(state.server, path)
            elseif canloadfile(state.server, path)
                loadfile(state.server, path)
            else
                path = ""
            end
        elseif !isempty(getpath(state.file)) && isabspath(joinpath(dirname(getpath(state.file)), path))
            # Relative path from current
            if hasfile(state.server, joinpath(dirname(getpath(state.file)), path))
                path = joinpath(dirname(getpath(state.file)), path)
            elseif canloadfile(state.server, joinpath(dirname(getpath(state.file)), path))
                path = joinpath(dirname(getpath(state.file)), path)
                loadfile(state.server, path)
            else
                path = ""
            end
        elseif !isempty((basepath = _is_in_basedir(getpath(state.file)); basepath))
            # Special handling for include method used within Base
            path = joinpath(basepath, path)
            if hasfile(state.server, path)
                # skip
            elseif canloadfile(state.server, path)
                loadfile(state.server, path)
            else
                path = ""
            end
        else
            path = ""
        end
        if hasfile(state.server, path)
            if path in state.included_files
                seterror!(x, IncludeLoop)
                return
            end
            oldfile = state.file
            state.file = getfile(state.server, path)
            push!(state.included_files, getpath(state.file))
            setroot(state.file, getroot(oldfile))
            setscope!(getcst(state.file), nothing)
            state(getcst(state.file))
            state.file = oldfile
            pop!(state.included_files)
        elseif !is_in_fexpr(x, CSTParser.defines_function) && !isempty(init_path)
            seterror!(x, MissingFile)
        end
    end
end

"""
    get_path(x::EXPR)

Usually called on the argument to `include` calls, and attempts to determine
the path of the file to be included. Has limited support for `joinpath` calls.
"""
function get_path(x::EXPR, state)
    if CSTParser.iscall(x) && length(x.args) == 2
        parg = x.args[2]

        if CSTParser.isstringliteral(parg)
            if occursin("\0", valof(parg))
                seterror!(parg, IncludePathContainsNULL)
                return ""
            end
            path = CSTParser.str_value(parg)
            path = normpath(path)
            Base.containsnul(path) && throw(SLInvalidPath("Couldn't convert '$x' into a valid path. Got '$path'"))
            return path
        elseif CSTParser.ismacrocall(parg) && valof(parg.args[1]) == "@raw_str" && CSTParser.isstringliteral(parg.args[3])
            if occursin("\0", valof(parg.args[3]))
                seterror!(parg.args[3], IncludePathContainsNULL)
                return ""
            end
            path = normpath(CSTParser.str_value(parg.args[3]))
            Base.containsnul(path) && throw(SLInvalidPath("Couldn't convert '$x' into a valid path. Got '$path'"))
            return path
        elseif CSTParser.iscall(parg) && isidentifier(parg.args[1]) && valofid(parg.args[1]) == "joinpath"
            path_elements = String[]

            for i = 2:length(parg.args)
                arg = parg[i]
                if _is_macrocall_to_BaseDIR(arg) # Assumes @__DIR__ points to Base macro.
                    push!(path_elements, dirname(getpath(state.file)))
                elseif CSTParser.isstringliteral(arg)
                    if occursin("\0", valof(arg))
                        seterror!(arg, IncludePathContainsNULL)
                        return ""
                    end
                    push!(path_elements, string(valof(arg)))
                else
                    return ""
                end
            end
            isempty(path_elements) && return ""

            path = normpath(joinpath(path_elements...))
            Base.containsnul(path) && throw(SLInvalidPath("Couldn't convert '$x' into a valid path. Got '$path'"))
            return path
        end
    end
    return ""
end

include("server.jl")
include("imports.jl")
include("references.jl")
include("macros.jl")
include("linting/checks.jl")
include("type_inf.jl")
include("utils.jl")
include("interface.jl")
end


text:
## Step size rules

import Base: step, vecnorm

export StepSizeRule,
	ConstantStepSize, BacktrackingStepSize,
	HopefulStepSize, DecreasingStepSize,
	step

## In general

abstract StepSizeRule
stop(ssr::StepSizeRule) = true # whether your should ever stop before maxiters

## More specifically,

type ConstantStepSize<:StepSizeRule
	stepsize
end
step(s::ConstantStepSize, args...; kwargs...) = s.stepsize
ConstantStepSize() = ConstantStepSize(1)

type DecreasingStepSize<:StepSizeRule
	initial_stepsize::AbstractFloat
	iteration_counter::Int
end
function step(s::DecreasingStepSize, args...; kwargs...)
	s.iteration_counter += 1
	return s.initial_stepsize/s.iteration_counter
end
DecreasingStepSize() = DecreasingStepSize(1.0)
DecreasingStepSize(initial_stepsize) = DecreasingStepSize(initial_stepsize, 0)
stop(s::DecreasingStepSize) = false

# this stepsize starts out every iteration at initial_stepsize and backtracks until it finds sufficient decrease
type BacktrackingStepSize<:StepSizeRule
	initial_stepsize::Float64
	decrease_by::Float64
	suff_decrease::Float64
end
function step(s::BacktrackingStepSize, objective::Function, x0, grad_x0;
	objval = objective(x0),
	normgradsq = vecnorm(grad_x0)^2)
	stepsize = s.initial_stepsize
	while objective(x0 - stepsize*grad_x0) >= objval - s.suff_decrease*stepsize*normgradsq
		stepsize *= s.decrease_by
	end
	return stepsize
end
BacktrackingStepSize() = BacktrackingStepSize(1.0, .9, .5)

# This stepsize gets bigger if we keep getting lucky and smaller otherwise
type HopefulStepSize<:StepSizeRule
	initial_stepsize::Float64
	max_stepsize::Float64
	decrease_by::Float64
	increase_by::Float64
	suff_decrease::Float64
end
function step(s::HopefulStepSize, objective::Function, x0, grad_x0;
	objval = objective(x0),
	normgradsq = vecnorm(grad_x0)^2)
	stepsize = s.initial_stepsize
	while objective(x0 - stepsize*grad_x0) > objval - s.suff_decrease*stepsize*normgradsq
		stepsize *= s.decrease_by
	end
	# hope and change!
	if stepsize < s.initial_stepsize
		s.initial_stepsize = stepsize
	else
		s.initial_stepsize *= s.increase_by
	end
	return stepsize
end
HopefulStepSize() = HopefulStepSize(1.0, Inf, .8, 1.5, .1)
HopefulStepSize(initial_stepsize) = HopefulStepSize(initial_stepsize, Inf, .8, 1.5, .1)

max_stars_repo_name: madeleineudell/FirstOrderOptimization


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.1
#     language: julia
#     name: julia-1.6
# ---

# %%
ci = CartesianIndex(2, 3, 4)

# %%
dump(ci)

# %%
Tuple(ci)

# %%
collect(Tuple(ci))

# %%

max_stars_repo_path: 0009/convert CartesianIndex to Tuple and Vector.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb04dp(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    NMAX = 10
    LDA  = NMAX
    LDC = NMAX
    LDDE = NMAX
    LDVW = NMAX
    A = Array{Float64,2}(undef, LDA,NMAX)
    DE = Array{Float64,2}(undef, LDDE,NMAX+1)
    C = Array{Float64,2}(undef, LDC,NMAX)
    VW = Array{Float64,2}(undef, LDVW,NMAX+1)
    LSCALE = Array{Float64,1}(undef, NMAX)
    RSCALE = Array{Float64,1}(undef, NMAX)
    DWORK = Array{Float64,1}(undef, 8*NMAX)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    N = parse(BlasInt, vs[1])
    JOB =  vs[2][1]
    THRESH = parse(Float64, replace(vs[3],'D'=>'E'))
    if ( N<=0 || N>NMAX )
        @error "Illegal N=$N"
    end

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       A[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (N,N+1)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       DE[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       C[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (N,N+1)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       VW[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
    close(f)
# interp call 1

    ILO, INFO, IWARN = SLICOT.mb04dp!(JOB, N, THRESH, A, DE, C, VW, LSCALE, RSCALE, DWORK)
    @test INFO == 0
    INFO == 0 || return
    println(io, "ILO = $ILO")

# interp output 1
    println(io, "A:")
    _nc = N
    _nr = N
    show(io, "text/plain", A[1:_nr,1:_nc])
    println(io)

# interp output 2
    println(io, "DE:")
    _nc = N+1
    _nr = N
    show(io, "text/plain", DE[1:_nr,1:_nc])
    println(io)

# interp output 3
    println(io, "C:")
    _nc = N
    _nr = N
    show(io, "text/plain", C[1:_nr,1:_nc])
    println(io)

# interp output 4
    println(io, "VW:")
    _nc = N+1
    _nr = N
    show(io, "text/plain", VW[1:_nr,1:_nc])
    println(io)

# interp output 5
    println(io, "LSCALE:")
    _nr = N
    show(io, "text/plain", LSCALE[1:_nr])
    println(io)

# interp output 6
    println(io, "RSCALE:")
    _nr = N
    show(io, "text/plain", RSCALE[1:_nr])
    println(io)

if ( LSAME( JOB, 'S' ) || LSAME( JOB, 'B' ) )
if ( !( THRESH==-2 || THRESH==-4 ) )
    println(io, "initial norms: ", DWORK[1:2])
    println(io, "initial norms: ", DWORK[3:4])
    println(io, "final threshold: ", DWORK[5])
else
    println(io, "IWARN = $IWARN")
end # if
end # if

end # run_mb04dp()


text:
"""
	gmtsimplify(cmd0::String="", arg1=[], kwargs...)

Line reduction using the Douglas-Peucker algorithm.

Full option list at [`gmtsimplify`](http://gmt.soest.hawaii.edu/doc/latest/gmtsimplify.html)

Parameters
----------

- **T** : **tol** : **tolerance** : -- Number or Str --    Flags = tolerance[unit]

    Specifies the maximum mismatch tolerance in the user units. If the data is not Cartesian then append the distance unit.
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/gmtsimplify.html#t)
- $(GMT.opt_V)
- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_b)
- $(GMT.opt_d)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_o)
- $(GMT.opt_swap_xy)
"""
function gmtsimplify(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && return monolitic("gmtsimplify", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)

	cmd = parse_V_params("", d)
	cmd, = parse_b(cmd, d)
	cmd, = parse_d(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_g(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, = parse_i(cmd, d)
	cmd, = parse_o(cmd, d)
	cmd, = parse_swap_xy(cmd, d)

	cmd = add_opt(cmd, 'T', d, [:T :tol :tolerance])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	return common_grd(d, cmd, got_fname, 1, "gmtsimplify", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
gmtsimplify(arg1=[], cmd0::String=""; kw...) = gmtsimplify(cmd0, arg1; kw...)


text:
module GaussianMixtureTest

using Test
using ForneyLab
using ForneyLab: outboundType, isApplicable
using ForneyLab: VBGaussianMixtureM, VBGaussianMixtureW, VBGaussianMixtureZBer, VBGaussianMixtureZCat, VBGaussianMixtureOut

#-------------
# Update rules
#-------------

@testset "VBGaussianMixtureM" begin
    @test VBGaussianMixtureM <: NaiveVariationalRule{GaussianMixture}
    @test outboundType(VBGaussianMixtureM) == Message{GaussianMeanPrecision}
    @test !isApplicable(VBGaussianMixtureM, [Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution])
    @test !isApplicable(VBGaussianMixtureM, [ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution])
    @test !isApplicable(VBGaussianMixtureM, [ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution])
    @test isApplicable(VBGaussianMixtureM, [ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution])
    @test isApplicable(VBGaussianMixtureM, [ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution])
    @test isApplicable(VBGaussianMixtureM, [ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution])

    @test ruleVBGaussianMixtureM(ProbabilityDistribution(Univariate, GaussianMeanVariance, m=8.5, v=0.5), ProbabilityDistribution(Univariate, Bernoulli, p=0.2), nothing, ProbabilityDistribution(Univariate, Gamma, a=1.0, b=2.0), ProbabilityDistribution(Univariate, GaussianMeanVariance, m=0.0, v=1.0), ProbabilityDistribution(Univariate, Gamma, a=2.0, b=1.0)) == Message(Univariate, GaussianMeanPrecision, m=8.5, w=0.1)
    @test ruleVBGaussianMixtureM(ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[8.5], v=mat(0.5)), ProbabilityDistribution(Univariate, Bernoulli, p=0.2), nothing, ProbabilityDistribution(MatrixVariate, Wishart, nu=2.0, v=mat(0.25)), ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[0.0], v=mat(1.0)), ProbabilityDistribution(MatrixVariate, Wishart, nu=4.0, v=mat(0.5))) == Message(Multivariate, GaussianMeanPrecision, m=[8.5], w=mat(0.1))
    @test ruleVBGaussianMixtureM(ProbabilityDistribution(Univariate, GaussianMeanVariance, m=8.5, v=0.5), ProbabilityDistribution(Univariate, Bernoulli, p=0.2), ProbabilityDistribution(Univariate, GaussianMeanVariance, m=0.0, v=1.0), ProbabilityDistribution(Univariate, Gamma, a=1.0, b=2.0), nothing, ProbabilityDistribution(Univariate, Gamma, a=2.0, b=1.0)) == Message(Univariate, GaussianMeanPrecision, m=8.5, w=1.6)
    @test ruleVBGaussianMixtureM(ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[8.5], v=mat(0.5)), ProbabilityDistribution(Univariate, Bernoulli, p=0.2), ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[0.0], v=mat(1.0)), ProbabilityDistribution(MatrixVariate, Wishart, nu=2.0, v=mat(0.25)), nothing, ProbabilityDistribution(MatrixVariate, Wishart, nu=4.0, v=mat(0.5))) == Message(Multivariate, GaussianMeanPrecision, m=[8.5], w=mat(1.6))
    @test ruleVBGaussianMixtureM(ProbabilityDistribution(Univariate, GaussianMeanVariance, m=8.5, v=0.5), ProbabilityDistribution(Univariate, Categorical, p=[0.2, 0.8]), nothing, ProbabilityDistribution(Univariate, Gamma, a=1.0, b=2.0), ProbabilityDistribution(Univariate, GaussianMeanVariance, m=0.0, v=1.0), ProbabilityDistribution(Univariate, Gamma, a=2.0, b=1.0)) == Message(Univariate, GaussianMeanPrecision, m=8.5, w=0.1)
    @test ruleVBGaussianMixtureM(ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[8.5], v=mat(0.5)), ProbabilityDistribution(Univariate, Categorical, p=[0.2, 0.8]), nothing, ProbabilityDistribution(MatrixVariate, Wishart, nu=2.0, v=mat(0.25)), ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[0.0], v=mat(1.0)), ProbabilityDistribution(MatrixVariate, Wishart, nu=4.0, v=mat(0.5))) == Message(Multivariate, GaussianMeanPrecision, m=[8.5], w=mat(0.1))
end

@testset "VBGaussianMixtureW" begin
    @test VBGaussianMixtureW <: NaiveVariationalRule{GaussianMixture}
    @test outboundType(VBGaussianMixtureW) == Message{Union{Gamma, Wishart}}
    @test !isApplicable(VBGaussianMixtureW, [Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution])
    @test !isApplicable(VBGaussianMixtureW, [ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution])
    @test !isApplicable(VBGaussianMixtureW, [ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution])
    @test isApplicable(VBGaussianMixtureW, [ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution])
    @test isApplicable(VBGaussianMixtureW, [ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing])
    @test isApplicable(VBGaussianMixtureW, [ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, Nothing])

    @test ruleVBGaussianMixtureW(ProbabilityDistribution(Univariate, GaussianMeanVariance, m=8.5, v=0.5), ProbabilityDistribution(Univariate, Bernoulli, p=0.2), ProbabilityDistribution(Univariate, GaussianMeanVariance, m=5.0, v=2.0), nothing, ProbabilityDistribution(Univariate, GaussianMeanVariance, m=10.0, v=3.0), ProbabilityDistribution(Univariate, Gamma, a=2.0, b=1.0)) == Message(Univariate, Gamma, a=1.1, b=1.475)
    @test ruleVBGaussianMixtureW(ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[8.5], v=mat(0.5)), ProbabilityDistribution(Univariate, Bernoulli, p=0.2), ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[5.0], v=mat(2.0)), nothing, ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[10.0], v=mat(3.0)), ProbabilityDistribution(MatrixVariate, Wishart, nu=4.0, v=mat(0.5))) == Message(MatrixVariate, Wishart, nu=2.2, v=mat(0.33898305084745767))
    @test ruleVBGaussianMixtureW(ProbabilityDistribution(Univariate, GaussianMeanVariance, m=8.5, v=0.5), ProbabilityDistribution(Univariate, Bernoulli, p=0.2), ProbabilityDistribution(Univariate, GaussianMeanVariance, m=5.0, v=2.0), ProbabilityDistribution(Univariate, Gamma, a=1.0, b=2.0), ProbabilityDistribution(Univariate, GaussianMeanVariance, m=10.0, v=3.0), nothing) == Message(Univariate, Gamma, a=1.4, b=2.3000000000000003)
    @test ruleVBGaussianMixtureW(ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[8.5], v=mat(0.5)), ProbabilityDistribution(Univariate, Bernoulli, p=0.2), ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[5.0], v=mat(2.0)), ProbabilityDistribution(MatrixVariate, Wishart, nu=2.0, v=mat(0.25)), ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[10.0], v=mat(3.0)), nothing) == Message(MatrixVariate, Wishart, nu=2.8, v=mat(0.2173913043478261))
    @test ruleVBGaussianMixtureW(ProbabilityDistribution(Univariate, GaussianMeanVariance, m=8.5, v=0.5), ProbabilityDistribution(Univariate, Categorical, p=[0.2, 0.8]), ProbabilityDistribution(Univariate, GaussianMeanVariance, m=5.0, v=2.0), nothing, ProbabilityDistribution(Univariate, GaussianMeanVariance, m=10.0, v=3.0), ProbabilityDistribution(Univariate, Gamma, a=2.0, b=1.0)) == Message(Univariate, Gamma, a=1.1, b=1.475)
    @test ruleVBGaussianMixtureW(ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[8.5], v=mat(0.5)), ProbabilityDistribution(Univariate, Categorical, p=[0.2, 0.8]), ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[5.0], v=mat(2.0)), nothing, ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[10.0], v=mat(3.0)), ProbabilityDistribution(MatrixVariate, Wishart, nu=4.0, v=mat(0.5))) == Message(MatrixVariate, Wishart, nu=2.2, v=mat(0.33898305084745767))
end

@testset "VBGaussianMixtureZBer" begin
    @test VBGaussianMixtureZBer <: NaiveVariationalRule{GaussianMixture}
    @test outboundType(VBGaussianMixtureZBer) == Message{Bernoulli}
    @test !isApplicable(VBGaussianMixtureZBer, [ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution])
    @test isApplicable(VBGaussianMixtureZBer, [ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution])

    @test ruleVBGaussianMixtureZBer(ProbabilityDistribution(Univariate, GaussianMeanVariance, m=8.5, v=0.5), nothing, ProbabilityDistribution(Univariate, GaussianMeanVariance, m=5.0, v=2.0), ProbabilityDistribution(Univariate, Gamma, a=1.0, b=2.0), ProbabilityDistribution(Univariate, GaussianMeanVariance, m=10.0, v=3.0), ProbabilityDistribution(Univariate, Gamma, a=2.0, b=1.0)) == Message(Univariate, Bernoulli, p=0.7713458788198754)
    @test ruleVBGaussianMixtureZBer(ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[8.5], v=mat(0.5)), nothing, ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[5.0], v=mat(2.0)), ProbabilityDistribution(MatrixVariate, Wishart, nu=2.0, v=mat(0.25)), ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[10.0], v=mat(3.0)), ProbabilityDistribution(MatrixVariate, Wishart, nu=4.0, v=mat(0.5))) == Message(Univariate, Bernoulli, p=0.7713458788198754)
end

@testset "VBGaussianMixtureZCat" begin
    @test VBGaussianMixtureZCat <: NaiveVariationalRule{GaussianMixture}
    @test outboundType(VBGaussianMixtureZCat) == Message{Categorical}
    @test isApplicable(VBGaussianMixtureZCat, [ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution])
    @test !isApplicable(VBGaussianMixtureZCat, [ProbabilityDistribution, Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution])

    @test ruleVBGaussianMixtureZCat(ProbabilityDistribution(Univariate, GaussianMeanVariance, m=8.5, v=0.5), nothing, ProbabilityDistribution(Univariate, GaussianMeanVariance, m=5.0, v=2.0), ProbabilityDistribution(Univariate, Gamma, a=1.0, b=2.0), ProbabilityDistribution(Univariate, GaussianMeanVariance, m=10.0, v=3.0), ProbabilityDistribution(Univariate, Gamma, a=2.0, b=1.0)) == Message(Univariate, Categorical, p=[0.7713458788198754, 0.22865412118012463])
    @test ruleVBGaussianMixtureZCat(ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[8.5], v=mat(0.5)), nothing, ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[5.0], v=mat(2.0)), ProbabilityDistribution(MatrixVariate, Wishart, nu=2.0, v=mat(0.25)), ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[10.0], v=mat(3.0)), ProbabilityDistribution(MatrixVariate, Wishart, nu=4.0, v=mat(0.5))) == Message(Univariate, Categorical, p=[0.7713458788198754, 0.22865412118012463])
end

@testset "VBGaussianMixtureOut" begin
    @test VBGaussianMixtureOut <: NaiveVariationalRule{GaussianMixture}
    @test outboundType(VBGaussianMixtureOut) == Message{GaussianWeightedMeanPrecision}
    @test isApplicable(VBGaussianMixtureOut, [Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution])
    @test isApplicable(VBGaussianMixtureOut, [Nothing, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution, ProbabilityDistribution])

    @test ruleVBGaussianMixtureOut(nothing, ProbabilityDistribution(Univariate, Bernoulli, p=0.2), ProbabilityDistribution(Univariate, GaussianMeanVariance, m=5.0, v=2.0), ProbabilityDistribution(Univariate, Gamma, a=1.0, b=2.0), ProbabilityDistribution(Univariate, GaussianMeanVariance, m=10.0, v=3.0), ProbabilityDistribution(Univariate, Gamma, a=2.0, b=1.0)) == Message(Univariate, GaussianWeightedMeanPrecision, xi=16.5, w=1.7)
    @test ruleVBGaussianMixtureOut(nothing, ProbabilityDistribution(Univariate, Bernoulli, p=0.2), ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[5.0], v=mat(2.0)), ProbabilityDistribution(MatrixVariate, Wishart, nu=2.0, v=mat(0.25)), ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[10.0], v=mat(3.0)), ProbabilityDistribution(MatrixVariate, Wishart, nu=4.0, v=mat(0.5))) == Message(Multivariate, GaussianWeightedMeanPrecision, xi=[16.5], w=mat(1.7))
    @test ruleVBGaussianMixtureOut(nothing, ProbabilityDistribution(Univariate, Categorical, p=[0.2, 0.8]), ProbabilityDistribution(Univariate, GaussianMeanVariance, m=5.0, v=2.0), ProbabilityDistribution(Univariate, Gamma, a=1.0, b=2.0), ProbabilityDistribution(Univariate, GaussianMeanVariance, m=10.0, v=3.0), ProbabilityDistribution(Univariate, Gamma, a=2.0, b=1.0)) == Message(Univariate, GaussianWeightedMeanPrecision, xi=16.5, w=1.7)
    @test ruleVBGaussianMixtureOut(nothing, ProbabilityDistribution(Univariate, Categorical, p=[0.2, 0.8]), ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[5.0], v=mat(2.0)), ProbabilityDistribution(MatrixVariate, Wishart, nu=2.0, v=mat(0.25)), ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[10.0], v=mat(3.0)), ProbabilityDistribution(MatrixVariate, Wishart, nu=4.0, v=mat(0.5))) == Message(Multivariate, GaussianWeightedMeanPrecision, xi=[16.5], w=mat(1.7))
end

@testset "averageEnergy" begin
    # Univariate
    marg_out = ProbabilityDistribution(Univariate, GaussianMeanVariance, m=0.0, v=1.0)
    marg_switch = ProbabilityDistribution(Univariate, Bernoulli, p=0.2)
    marg_mean_1 = ProbabilityDistribution(Univariate, GaussianMeanVariance, m=1.0, v=2.0)
    marg_prec_1 = ProbabilityDistribution(Univariate, Gamma, a=2.0, b=3.0)
    marg_mean_2 = ProbabilityDistribution(Univariate, GaussianMeanVariance, m=3.0, v=4.0)
    marg_prec_2 = ProbabilityDistribution(Univariate, Gamma, a=4.0, b=5.0)
    ref_val = 0.2*averageEnergy(GaussianMeanPrecision, marg_out, marg_mean_1, marg_prec_1) +
              0.8*averageEnergy(GaussianMeanPrecision, marg_out, marg_mean_2, marg_prec_2)

    @test averageEnergy(GaussianMixture, marg_out, marg_switch, marg_mean_1, marg_prec_1, marg_mean_2, marg_prec_2) == ref_val

    marg_switch = ProbabilityDistribution(Univariate, Categorical, p=[0.2, 0.8])
    @test averageEnergy(GaussianMixture, marg_out, marg_switch, marg_mean_1, marg_prec_1, marg_mean_2, marg_prec_2) == ref_val

    # Multivariate
    marg_out = ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[0.0], v=mat(1.0))
    marg_switch = ProbabilityDistribution(Univariate, Bernoulli, p=0.2)
    marg_mean_1 = ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[1.0], v=mat(2.0))
    marg_prec_1 = ProbabilityDistribution(MatrixVariate, Wishart, nu=4.0, v=mat(1/6))
    marg_mean_2 = ProbabilityDistribution(Multivariate, GaussianMeanVariance, m=[3.0], v=mat(4.0))
    marg_prec_2 = ProbabilityDistribution(MatrixVariate, Wishart, nu=8.0, v=mat(0.1))
    ref_val = 0.2*averageEnergy(GaussianMeanPrecision, marg_out, marg_mean_1, marg_prec_1) +
              0.8*averageEnergy(GaussianMeanPrecision, marg_out, marg_mean_2, marg_prec_2)

    @test averageEnergy(GaussianMixture, marg_out, marg_switch, marg_mean_1, marg_prec_1, marg_mean_2, marg_prec_2) == ref_val

    marg_switch = ProbabilityDistribution(Univariate, Categorical, p=[0.2, 0.8])
    @test averageEnergy(GaussianMixture, marg_out, marg_switch, marg_mean_1, marg_prec_1, marg_mean_2, marg_prec_2) == ref_val
end

end #module


text:
push!(LOAD_PATH, "../src/")
using Documenter, PanGraph

makedocs(
    sitename = "PanGraph.jl",
    authors = "Nicholas Noll, Marco Molari, Richard Neher",
    modules = [PanGraph],
    pages = [
        "Home" => "index.md",
        "Tutorial" => [
            "tutorials/tutorial_1.md",
            "tutorials/tutorial_2.md",
            "tutorials/tutorial_3.md",
            "tutorials/tutorial_4.md",
        ],
        "Library" => [
            "lib/align.md",
            "lib/block.md",
            "lib/edge.md",
            "lib/graph.md",
            "lib/gfa.md",
            "lib/mash.md",
            "lib/minimap.md",
            "lib/node.md",
            "lib/path.md",
            "lib/simulate.md",
            "lib/utility.md",
        ],
        "Command Line" => [
            "cli/build.md",
            "cli/export.md",
            "cli/generate.md",
            "cli/marginalize.md",
            "cli/polish.md",
        ]
    ]
)

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
"""
    luhn(str)

Return true if `str` encodes a valid number in the luhn formula.

A valid string will:

 - Contain only 0-9 and space
 - Contain at least two numbers
 - After a transformation, the sum of the transformed digits will be evenly divisible by 10

The transformation is:

 - Every second digit, counting the rightmost digit as first, is doubled
 - If a doubled digit is greater than 9, subtract 9 from it.
"""
function luhn(str)
    acc = 0
    len = 0
    for chr in Iterators.reverse(str)
        if isdigit(chr)
            len += 1
            x = parse(Int, chr)
            acc += iseven(len) ? (2x < 10 ? 2x : 2x-9) : x
        elseif chr != ' '
            return false
        end
    end
    len > 1 ? acc % 10 == 0 : false
end