text:
export rhf

"""
1e integrals; single core version.
"""
function all_1e_ints_singlecore(bfs::BasisSet,mol::Molecule)
    n = length(bfs.bfs)
    S = Array{Float64}(n,n)
    T = Array{Float64}(n,n)
    V = Array{Float64}(n,n)
    for (i,j) in pairs(n)
        a,b = bfs.bfs[i],bfs.bfs[j]
        S[i,j] = S[j,i] = overlap(a,b)
        T[i,j] = T[j,i] = kinetic(a,b)
        V[i,j] = V[j,i] = nuclear_attraction(a,b,mol)
    end
    return S,T,V
end

function all_1e_ints(bfs::BasisSet,mol::Molecule)
    n = length(bfs.bfs)
    S = SharedArray{Float64}(n,n)
    T = SharedArray{Float64}(n,n)
    V = SharedArray{Float64}(n,n)
    @sync @distributed for i = 1:n
        @sync @distributed for j = 1:i
    #for (i,j) in pairs(n)
        a,b = bfs.bfs[i],bfs.bfs[j]
        S[i,j] = S[j,i] = overlap(a,b)
        T[i,j] = T[j,i] = kinetic(a,b)
        V[i,j] = V[j,i] = nuclear_attraction(a,b,mol)
        end
    end
    return S,T,V
end

"""
2e integrals; single core version.
"""
function all_2e_ints_singlecore(bflist,ERI=coulomb_hgp)
    n = length(bflist.bfs)
    totlen = div(n*(n+1)*(n*n+n+2),8)
    ints2e = SharedArray{Float64}(totlen)
    for (i,j,k,l) in iiterator(n)
        ints2e[iindex(i,j,k,l)] = ERI(bflist.bfs[i],bflist.bfs[j],bflist.bfs[k],bflist.bfs[l])
    end
    return ints2e
end

function all_2e_ints(bflist,ERI=coulomb_hgp;verbose::Bool=false) #coloumb
    n = length(bflist.bfs)
    totlen = div(n*(n+1)*(n*n+n+2),8)
    ints2e = SharedArray{Float64}(totlen)
    #for (i,j,k,l) in iiterator(n)
    if verbose println(n^4/4," ints to calculate. Each '.' is ",n^2/2) end
    @sync @distributed for i=1:n # Nb: need sync on both these parallel, to force completion of Ints before SCF
        @sync @distributed for j=1:i
            for (k,l) in pairs(n)
                if triangle(i,j) <= triangle(k,l)
                    ints2e[iindex(i,j,k,l)] = ERI(bflist.bfs[i],bflist.bfs[j],bflist.bfs[k],bflist.bfs[l])
                end
            end
            if verbose print(".") end
        end
    end
    return ints2e
end

"""
 make2JmK

 Coulomb and exchange contribution to the Fock matrix.
 The N^4 component of Hartree Fock.

 TODO: Optimise me :^)
 Nb: Now accepts SharedArray (parallel constructed) Ints
"""
function make2JmK(D::Array{Float64,2},Ints::SharedArray{Float64,1})
    n = size(D,1)
    G = Array{Float64}(undef,n,n) # undef is for 1.0 compat. Unsure what it does.
    D1 = reshape(D,n*n)
    temp = Array{Float64}(undef,n*n)
    for (i,j) in pairs(n)
        kl = 1
        for (k,l) in rpairs(n)
            temp[kl] = 2*Ints[iindex(i,j,k,l)]-Ints[iindex(i,k,j,l)]
            kl += 1
        end
        G[i,j] = G[j,i] = dot(D1,temp)
    end
    return G
end

densitymatrix(U::Array{Float64,2},nocc::Int64) = U[:,1:nocc]*U[:,1:nocc]'

"""
 rhf

Restricted Hartree-Fock function. Will run a self-consistent field (SCF)
calculation on the supplied molecular geometry.

Algorithm (and comments herein) follow Thijssen 2007 2nd Ed pp70-73

"""
function rhf(mol::Molecule,MaxIter::Int64=40; verbose::Bool=false, Econvergence::Float64=1e-6)
    if verbose println("Building basis") end
    bfs = build_basis(mol)
    if verbose println(length(bfs.bfs)," basis functions") end
    # S = Overlap matrix
    # T = one-electron Kinetic Energy
    # V = one-electron potential
    if verbose println("Calculating 1e ints") end
    S,T,V = all_1e_ints(bfs,mol)
    # Ints = two-electron integrals, <pr|g|qs>
    if verbose println("Calculating 2e ints") end
    Ints = all_2e_ints(bfs,verbose=verbose)
    # h = Form one-eletron Hamiltonian
    h = T+V
    # generalised eigenvalue decomposition of one-electron hamiltonian and overlap matrix
    # used as starting guess for self-consistent procedure?
    E,U = eigen(h,S)
    # Enuke = (classical) nuclear repulision
    Enuke = nuclear_repulsion(mol)
    # Define occupied and virtual orbitals
    nclosed,nopen = divrem(nel(mol),2)
    D=densitymatrix(U,nclosed) # initial density matrix
    Eold = Inf
    Energy = 0
    println("Nel=$(nel(mol)) Nclosed=$nclosed")
    if verbose
        println("S=\n$S")
        println("h=\n$h")
        println("T=\n$T")
        println("V=\n$V")
        println("E: $E")
        println("U: $U")
        println("2e ints:\n$Ints")
    end
    println("SCF: Iteration TotalEnergy :=: Enuke + Eone + Etwo")
    for iter in 1:MaxIter
        α=0.4 # mixing of current and prior density matrices
        D = α*D + (1-α)*densitymatrix(U,nclosed) # density matrix; from eigenvalues U
        if verbose
            println("D=\n$D")
        end
        # Coulomb and Exchange contributions to Fock Matrix:
        # This is the N^4, and most time consuming, step of Hartree Fock.
        G = make2JmK(D,Ints)
        # Fock matrix H = one-electron hamiltonian + G matrix
        H = h+G
        # Diagonalise the Fock Matrix ; eig=Generalised eigenvalue decomposition
        E,U = eigen(H,S)
        # expanded versions of trace2 functions
        Eone = sum(D.*h) # one electron contribution
        Etwo = sum(D.*H) # should be = sum of occupied Fock orbitals
        Energy = Enuke + Eone + Etwo #c.f. Thijssen, factor of 2 already taken care of?
        println("HF: $iter  $Energy : $Enuke    $Eone    $Etwo")
        if abs(Energy-Eold)<Econvergence
            break
        end
        Eold  = Energy
    end
    return Energy,E,U
end


text:
# by Lazaro Alonso
# taken from: https://lazarusa.github.io/BeautifulMakie/GeoPlots/submarineCables3D/
using GeoMakie, Downloads
using GeoJSON, GeoInterface
using HDF5, FileIO
using RPRMakie
# data from
# https://github.com/telegeography/www.submarinecablemap.com
urlPoints = "https://raw.githubusercontent.com/telegeography/www.submarinecablemap.com/master/web/public/api/v3/landing-point/landing-point-geo.json"
urlCables = "https://raw.githubusercontent.com/telegeography/www.submarinecablemap.com/master/web/public/api/v3/cable/cable-geo.json"

landPoints = Downloads.download(urlPoints, IOBuffer())
landCables = Downloads.download(urlCables, IOBuffer())

land_geoPoints = GeoJSON.read(seekstart(landPoints))
land_geoCables = GeoJSON.read(seekstart(landCables))

toPoints = GeoMakie.geo2basic(land_geoPoints)
#toLines = GeoMakie.geo2basic(land_geoCables) # this should probably be supported.
feat = GeoInterface.features(land_geoCables)
toLines = GeoInterface.coordinates.(GeoInterface.geometry.(feat))

# broken lines at -180 and 180... they should
# be the same line and be in the same array.

# also the coastlines data
# this data can be download from here:
# https://github.com/lazarusA/BeautifulMakie/tree/main/_assets/data
worldxm = h5open(joinpath(@__DIR__, "world_xm.h5"), "r")
world_110m = read(worldxm["world_110m"])
close(worldxm)
lonw, latw = world_110m["lon"], world_110m["lat"]
# some 3D transformations
function toCartesian(lon, lat; r=1.02, cxyz=(0, 0, 0))
    lat, lon = lat * π / 180, lon * π / 180
    return cxyz[1] + r * cos(lat) * cos(lon), cxyz[2] + r * cos(lat) * sin(lon), cxyz[3] + r * sin(lat)
end

function lonlat3D2(lon, lat; cxyz=(0, 0, 0))
    xyzw = zeros(length(lon), 3)
    for (i, lon) in enumerate(lon)
        x, y, z = toCartesian(lon, lat[i]; cxyz=cxyz)
        xyzw[i, 1] = x
        xyzw[i, 2] = y
        xyzw[i, 3] = z
    end
    return xyzw[:, 1], xyzw[:, 2], xyzw[:, 3]
end

xw, yw, zw = lonlat3D2(lonw, latw)
toPoints3D = [Point3f([toCartesian(point[1], point[2])...]) for point in toPoints]

splitLines3D = []
for i in 1:length(toLines)
    for j in 1:length(toLines[i])
        ptsLines = toLines[i][j]
        tmp3D = []
        for k in 1:length(ptsLines)
            x, y = ptsLines[k]
            x, y, z = toCartesian(x, y)
            push!(tmp3D, [x, y, z])
        end
        push!(splitLines3D, Point3f.(tmp3D))
    end
end

earth_img = load(Downloads.download("https://upload.wikimedia.org/wikipedia/commons/5/56/Blue_Marble_Next_Generation_%2B_topography_%2B_bathymetry.jpg"))
Makie.inline!(true)
# the actual plot !
RPRMakie.activate!(; iterations=100)
scene = with_theme(theme_dark()) do
    fig = Figure(; resolution=(1000, 1000))
    radiance = 30
    lights = [EnvironmentLight(0.5, load(RPR.assetpath("starmap_4k.tif"))),
              PointLight(Vec3f(1, 1, 3), RGBf(radiance, radiance, radiance))]
    ax = LScene(fig[1, 1]; show_axis=false, scenekw=(lights=lights,))
    # lines!(ax, xw, yw, zw, color = :white, linewidth = 0.15)
    n = 1024 ÷ 4 # 2048
    θ = LinRange(0, pi, n)
    φ = LinRange(-pi, pi, 2 * n)
    xe = [cos(φ) * sin(θ) for θ in θ, φ in φ]
    ye = [sin(φ) * sin(θ) for θ in θ, φ in φ]
    ze = [cos(θ) for θ in θ, φ in φ]
    surface!(ax, xe, ye, ze; color=earth_img)
    meshscatter!(toPoints3D; color=1:length(toPoints3D), markersize=0.005, colormap=:plasma)
    colors = Makie.default_palettes.color[]
    c = Iterators.cycle(colors)
    foreach(((l, c),) -> lines!(ax, l; linewidth=2, color=c), zip(splitLines3D, c))
    ax.scene.camera_controls.eyeposition[] = Vec3f(1.5)
    return ax.scene
end

save("submarin_cables.png", scene)

ext: jl


text:
using Modia

HeatPort = Model( T = potential,   # Absolute temperature
                  Q_flow = flow )  # Heat flow into the component

HeatPort = HeatPort | Map(T=Var(start=273.15u"K"),Q_flow=Var(start=0.0u"W"))


text:
immutable RNN <: Model; net; RNN(a...;o...)=new(Net(a...;o...)); end
params(m::RNN)=params(m.net)

function train(m::RNN, data, loss; gclip=0, gcheck=false, maxnorm=nothing, losscnt=nothing, o...)
    reset!(m.net; o...)
    ystack = Any[]
    for item in data
        if item != nothing
            (x,ygold) = item2xy(item)
            ypred = forw(m.net, x...; trn=true, seq=true, o...)
            losscnt != nothing && (losscnt[1] += loss(ypred, ygold); losscnt[2] += 1)
            push!(ystack, copy(ygold))
        else                    # end of sequence
            while !isempty(ystack)
                ygold = pop!(ystack)
                back(m.net, ygold, loss; seq=true, o...)
            end
            gcheck && break
            g = (gclip > 0 || maxnorm!=nothing ? gnorm(m) : 0)
            update!(m; gclip=(g > gclip > 0 ? gclip/g : 0), o...)
            if maxnorm != nothing
                w=wnorm(m)
                w > maxnorm[1] && (maxnorm[1]=w)
                g > maxnorm[2] && (maxnorm[2]=g)
            end
            reset!(m.net; o...)
        end
    end
end

function test(m::RNN, data, loss; o...)
    sumloss = numloss = 0
    reset!(m.net; o...)
    for item in data
        if item != nothing
            (x,ygold) = item2xy(item)
            ypred = forw(m.net, x...; trn=false, seq=true, o...)
            sumloss += loss(ypred, ygold); numloss += 1
        else
            reset!(m.net; o...)
        end
    end
    return sumloss/numloss
end


text:
#= A basic genetic algorithm.
=#
using NKLandscapes

@doc """ function genetic_algorithm( ls::Landscape, pop_size::Int64, maxgens::Int64;
    selection_funct::Function=propsel!, tournsel_k::Int64=2,
    moran_sel_iters::Float64=1.0,   # Fraction of ls.n
    mutation_funct::Function=bwmutate!, mut_prob::Float64=1.0, # Mutation rate is mut_prob/ls.n
    init_genotype::Any=Void, termination_funct::Function=default_termination_funct,
    statistics_funct::Function=default_statistics_funct )

A basic generational genetic algorithm using the functions provided in NKLandscape.jl.
The non-optional arguments are the landscape, which defines the fitness function,
the populaiton size, and the maximum number of generations.
By default, the algorithm uses random initialization, proportional selection,
and bitwise mutation with rate 1/ls.n, and no crossover.  Optional arguments can
specify alternative methods of initialization, selection, and mutation, and termination.
A statistics function is called at the beginning of every generation and at termination.
"""

function genetic_algorithm( ls::Landscape, pop_size::Int64, maxgens::Int64;
    selection_funct::Function=propsel!,
    tournsel_k::Int64=2,   # number of individuals in a tournament
    moransel_iters::Float64=1.0,   # Fraction of ls.n
    mutation_funct::Function=bwmutate!, mut_prob::Float64=1.0, # Mutation rate is mut_prob/ls.n
    termination_funct::Function=default_termination_funct, init_genotype::Any=Void,
    statistics_funct::Function=default_statistics_funct )
  p::Population
  if typeof(init_genotype) <: Genotype
    p = constant(Population, rg, pop_size)  # Initialize to a single-genotype populaiton
  else
    p = rand(Population, ls, pop_size)
  end
  gen = 1
  while gen <= maxgens
    statistics_funct( gen, p )
    mutation_funct(p, mut_prob/ls.n)
    if selection_funct == tournsel!
      selection_funct(p, tournsel_k)
    elseif selection_funct == moransel!
      selection_funct(p, convert(Int,round(ls.n*moransel_iters)))
    else
      selection_funct(p)
    end
    if termination_funct(p)
      break
    end
    gen += 1
  end
  statistics_funct(gen, p )
  return p
end

@doc """ function default_termination_funct( p::Population )

If this termination function is used, the algorithm runs for maxgens generations.
"""
function default_termination_funct( p::Population )
  return false
end

@doc """ function test_converged( p::Population )

Terminate if the population has converged to a single genotype.
"""
function test_converged( p::Population )
  len = length(p.genotypes)
  i = 1
  while i <= len && p.genotypes[1] == p.genotypes[i]
    i += 1
  end
  return i > len ? true : false
end

@doc """ default_statistics_funct(gen::Int64, p::Population )

Modify to collect statistics about the run.
"""
function default_statistics_funct(gen::Int64, p::Population )
  return false
end

@doc """ simple (gen::Int64, p::Population )

Report the generation number and the maximum fitness at that generation.
"""
function simple_statistics_funct(gen::Int64, p::Population )
  pfits = popfits( p )
  println("Generation: ", gen, "  Max fitness: ", maximum(pfits) )
  #println(p)
  #println(map(fitness,p.genotypes))
end

@doc """ function test_genetic_algorithm(N::Int64, K::Int64 )
  Tests the various options for function genetic_algorithm().
  Note that the use of simple_statistics_funct  does give some printed output.
"""

function test_genetic_algorithm(N::Int64, K::Int64 )
  ls = NKLandscape(N,K)
  popsize = 8
  maxgens = 10
  p0 = genetic_algorithm(ls,popsize,maxgens)
  #println("p0:",p0)

  p1 = genetic_algorithm(ls,popsize,maxgens,statistics_funct=simple_statistics_funct)
  #println("p1:",p1)

  p2 = genetic_algorithm(ls,popsize,maxgens,selection_funct=moransel!,moransel_iters=2.0)
  #println("p2:",p2)

  p3 = genetic_algorithm(ls,popsize,maxgens,selection_funct=tournsel!,tournsel_k=3)
  #println("p3:",p3)

  p4 = genetic_algorithm(ls,popsize,maxgens,mutation_funct=bsmutate!,mut_prob=0.5)
  #println("p4:",p4)

  p5 = genetic_algorithm(ls,popsize,maxgens,termination_funct=test_converged)
  #println("p5:",p5)
end


text:
# This file is a part of SimilaritySearch.jl
export HausdorffDistance, MinHausdorffDistance

"""
    HausdorffDistance(dist::SemiMetric)

Hausdorff distance is defined as the maximum of the minimum between two clouds of points.

```math
Hausdorff(U, V) = \\max{\\max_{u \\in U} nndist(u, V), \\max{v \\in V} nndist(v, U) }
```

where ``nndist(u, V)`` computes the distance of ``u`` to its nearest neighbor in ``V`` using the `dist` metric.
"""
struct HausdorffDistance{Dtype<:SemiMetric} <: SemiMetric
    dist::Dtype
end

function _exhaustive_nndist(dist::SemiMetric, u::T, V) where T
    min_ = typemax(eltype(u))

    @inbounds for j in eachindex(V)
        min_ = min(evaluate(dist, u, V[j]), min_)
    end

    min_
end

function _hausdorff1(dist::SemiMetric, u, v)
    s = 0.0
    @inbounds for i in eachindex(u)
        s = max(s, _exhaustive_nndist(dist, u[i], v))
    end

    s
end

"""
    evaluate(m::HausdorffDistance, u, v)

Computes the Hausdorff distance between two cloud of points.

`u` and `v` are iterables where each object can be measured with the internal distance `dist`
"""
function evaluate(m::HausdorffDistance, u, v)
    if  length(u) == 1 || length(v) == 1
        _hausdorff1(m.dist, u, v)
    else
        max(_hausdorff1(m.dist, u, v), _hausdorff1(m.dist, v, u))
    end
end

"""
    MinHausdorffDistance(dist::SemiMetric)

Similar to [HausdorffDistance](@ref) but using minimum instead of maximum.
"""
struct MinHausdorffDistance{Dtype<:SemiMetric} <: SemiMetric
    dist::Dtype
end

function _minhausdorff1(dist::SemiMetric, u, v)
    s = 0.0
    @inbounds for i in eachindex(u)
        s = min(s, _exhaustive_nndist(dist, u[i], v))
    end

    s
end

"""
    evaluate(m::MinHausdorffDistance, u, v)

Computes a variant of the Hausdorff distance that uses the minimum instead of the maximum. `u` and `v` are iterables where each object can be measured with the internal distance `dist`
"""
function evaluate(m::MinHausdorffDistance, u, v)
    if  length(u) == 1 || length(v) == 1
        _minhausdorff1(m.dist, u, v)
    else
        min(_minhausdorff1(m.dist, u, v), _minhausdorff1(m.dist, v, u))
    end
end

max_stars_repo_name: sadit/SimilaritySearch

max_issues_repo_name: sadit/NNS

max_forks_repo_name: sadit/SimilaritySearch


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

abstract type AbstractClusterRule end

"""
    StaticCluster(clusters::Vector{Float64})

Buffered cuts are sorting according to the supplied weights `clusters`

"""
struct StaticCluster <: AbstractClusterRule
    clusters::Vector{Float64}

    function StaticCluster(clusters::Vector{Float64})
        sum(clusters) ≈ 1.0 || error("Static cluster")
        return new(clusters)
    end
end

function StaticCluster(frequencies::Vector{Int})
    return StaticCluster(frequencies ./ sum(frequencies))
end

function cluster(rule::StaticCluster, cuts::Vector{<:HyperPlane{OptimalityCut,T}}) where T <: AbstractFloat
    clusters::Vector{AggregatedOptimalityCut} = fill(zero(AggregatedOptimalityCut{T}), rule.clusters)
    i = 1
    for (j,n) in enumerate(round.(Int, rule.clusters * length(cuts)))
        for _ = 1:n
            clusters[j] += cuts[i]
            i += 1
            if i == length(cuts) + 1
                return clusters
            end
        end
    end
    return clusters
end

function str(::StaticCluster)
    return "static clustering"
end

struct ClusterByReference <: AbstractClusterRule
    τ::Float64
    distance::Function
end

"""
    ClusterByReference(τ::AbstractFloat; distance::Function = absolute_distance)

Buffered cuts are aggregated if within the tolerance `τ` to a reference cut, according the supplied `distance` function. Behaves as multi-cut otherwise.

The following distance measures are available
- [`absolute_distance`](@ref)
- [`angular_distance`](@ref)
- [`spatioangular_distance`](@ref

"""
function ClusterByReference(τ::AbstractFloat; distance::Function = absolute_distance)
    return ClusterByReference(τ, distance)
end

function cluster(rule::ClusterByReference, cuts::Vector{<:HyperPlane{OptimalityCut,T}}) where T <: AbstractFloat
    clusters::Vector{AggregatedOptimalityCut} = fill(zero(AggregatedOptimalityCut{T}), length(cuts) + 1)
    reference = zero(AggregatedOptimalityCut{T})
    for cut in cuts
        reference += cut
    end
    j = 2
    for cut in cuts
        if rule.distance(cut, reference) <= rule.τ
            # Cluster if distance to reference within tolerance
            clusters[1] += cut
        else
            # Multicut otherwise
            clusters[j] += cut
            j += 1
        end
    end
    return filter(c -> !iszero(c), clusters)
end

function str(::ClusterByReference)
    return "distance to reference based clustering"
end

"""
    Kmedoids(nclusters::Int; distance::Function = absolute_distance)

Buffered cuts are sorted into `nclusters` clusters, using a K-medoids algorithm over a generalized `distance` matrix.

The following distance measures are available
- [`absolute_distance`](@ref)
- [`angular_distance`](@ref)
- [`spatioangular_distance`](@ref

"""
struct Kmedoids <: AbstractClusterRule
    nclusters::Int
    distance::Function
end

function Kmedoids(nclusters::Int; distance::Function = absolute_distance)
    return Kmedoids(nclusters, distance)
end

function cluster(rule::Kmedoids, cuts::Vector{<:HyperPlane{OptimalityCut,T}}) where T <: AbstractFloat
    clusters::Vector{AggregatedOptimalityCut} = fill(zero(AggregatedOptimalityCut{T}), rule.nclusters)
    if length(cuts) < rule.nclusters
        for (i,cut) in enumerate(cuts)
            clusters[i] += cut
        end
    else
        D = [rule.distance(c₁, c₂) for c₁ in cuts, c₂ in cuts]
        try
            kmeds = kmedoids(D, rule.nclusters)
            for (i, j) in enumerate(assignments(kmeds))
                clusters[j] += cuts[i]
            end
        catch
            # Fallback to even clustering
            (ncuts,extra) = divrem(length(cuts), length(clusters))
            extra > 0 && (ncuts += 1)
            i = 1
            for j in eachindex(clusters)
                for _ in 1:ncuts
                    clusters[j] += cuts[i]
                    i += 1
                    i == length(cuts)+1 && return filter(c -> !iszero(c), clusters)
                end
            end
        end
    end
    return filter(c -> !iszero(c), clusters)
end

function str(::Kmedoids)
    return "K-medoids clustering"
end

"""
    Hierarchical(nclusters::Int; distance::Function = absolute_distance, linkage::Symbol = :single)

Buffered cuts are sorted into `nclusters` clusters, using a Hierarchical algorithm, with the given `linkage`, over a generalized `distance` matrix.

The following distance measures are available
- [`absolute_distance`](@ref)
- [`angular_distance`](@ref)
- [`spatioangular_distance`](@ref

"""
struct Hierarchical <: AbstractClusterRule
    nclusters::Int
    distance::Function
    linkage::Symbol
end

function Hierarchical(nclusters::Int; distance::Function = absolute_distance, linkage::Symbol = :single)
    return Hierarchical(nclusters, distance, linkage)
end

function cluster(rule::Hierarchical, cuts::Vector{<:HyperPlane{OptimalityCut,T}}) where T <: AbstractFloat
    clusters::Vector{AggregatedOptimalityCut} = fill(zero(AggregatedOptimalityCut{T}), rule.nclusters)
    D = [rule.distance(c₁, c₂) for c₁ in cuts, c₂ in cuts]
    tree = hclust(D, linkage = rule.linkage)
    for (i, j) in enumerate(Clustering.cutree(tree; k = rule.nclusters))
        clusters[j] += cuts[i]
    end
    return filter(c -> !iszero(c), clusters)
end

function str(::Hierarchical)
    return "hierarchical clustering"
end


text:
using Base.Test
using FlatBuffers

# Store specific byte patterns in these variables for the fuzzer. These
# values are taken verbatim from the C++ function FuzzTest1.
const overflowingInt32Val = read(IOBuffer(UInt8[0x83, 0x33, 0x33, 0x33]), Int32)
const overflowingInt64Val = read(IOBuffer(UInt8[0x84, 0x44, 0x44, 0x44, 0x44, 0x44, 0x44, 0x44]), Int64)

# CheckByteLayout verifies the bytes of a Builder in various scenarios.
function CheckByteLayout()
    check = want-> begin
        got = b.bytes[b.head+1:end]
        @test want == got
        return
    end

	# test 1: numbers

    b = FlatBuffers.Builder()
	check(UInt8[])
	FlatBuffers.prepend!(b, true)
	check(UInt8[1])
	FlatBuffers.prepend!(b, Int8(-127))
	check(UInt8[129, 1])
	FlatBuffers.prepend!(b, UInt8(255))
	check(UInt8[255, 129, 1])
	FlatBuffers.prepend!(b, Int16(-32222))
	check(UInt8[0x22, 0x82, 0, 255, 129, 1]) # first pad
	FlatBuffers.prepend!(b, UInt16(0xFEEE))
	check(UInt8[0xEE, 0xFE, 0x22, 0x82, 0, 255, 129, 1]) # no pad this time
	FlatBuffers.prepend!(b, Int32(-53687092))
	check(UInt8[204, 204, 204, 252, 0xEE, 0xFE, 0x22, 0x82, 0, 255, 129, 1])
	FlatBuffers.prepend!(b, UInt32(0x98765432))
	check(UInt8[0x32, 0x54, 0x76, 0x98, 204, 204, 204, 252, 0xEE, 0xFE, 0x22, 0x82, 0, 255, 129, 1])

	# test 1b: numbers 2

	b = FlatBuffers.Builder()
	prepend!(b, 0x1122334455667788)
	check(UInt8[0x88, 0x77, 0x66, 0x55, 0x44, 0x33, 0x22, 0x11])

	# test 2: 1xbyte vector

	b = FlatBuffers.Builder()
	check(UInt8[])
	FlatBuffers.startvector(b, sizeof(Bool), 1, 1)
	check(UInt8[0, 0, 0]) # align to 4bytes
	FlatBuffers.prepend!(b, UInt8(1))
	check(UInt8[1, 0, 0, 0])
	FlatBuffers.endvector(b, 1)
	check(UInt8[1, 0, 0, 0, 1, 0, 0, 0]) # padding

	# test 3: 2xbyte vector

	b = FlatBuffers.Builder()
	FlatBuffers.startvector(b, sizeof(UInt8), 2, 1)
	check(UInt8[0, 0]) # align to 4bytes
	FlatBuffers.prepend!(b, UInt8(1))
	check(UInt8[1, 0, 0])
	FlatBuffers.prepend!(b, UInt8(2))
	check(UInt8[2, 1, 0, 0])
	FlatBuffers.endvector(b, 2)
	check(UInt8[2, 0, 0, 0, 2, 1, 0, 0]) # padding

	# test 3b: 11xbyte vector matches builder size

	b = FlatBuffers.Builder(Any, 12)
	FlatBuffers.startvector(b, sizeof(UInt8), 8, 1)
	start = UInt8[]
	check(start)
    for i = 1:11
		FlatBuffers.prepend!(b, UInt8(i))
		start = append!(UInt8[i], start)
		check(start)
	end
	FlatBuffers.endvector(b, 8)
	check(append!(UInt8[8, 0, 0, 0], start))

	# test 4: 1xuint16 vector

	b = FlatBuffers.Builder()
	FlatBuffers.startvector(b, sizeof(UInt16), 1, 1)
	check(UInt8[0, 0]) # align to 4bytes
	FlatBuffers.prepend!(b, UInt16(1))
	check(UInt8[1, 0, 0, 0])
	FlatBuffers.endvector(b, 1)
	check(UInt8[1, 0, 0, 0, 1, 0, 0, 0]) # padding

	# test 5: 2xuint16 vector

	b = FlatBuffers.Builder()
	FlatBuffers.startvector(b, sizeof(UInt16), 2, 1)
	check(UInt8[]) # align to 4bytes
	FlatBuffers.prepend!(b, UInt16(0xABCD))
	check(UInt8[0xCD, 0xAB])
	FlatBuffers.prepend!(b, UInt16(0xDCBA))
	check(UInt8[0xBA, 0xDC, 0xCD, 0xAB])
	FlatBuffers.endvector(b, 2)
	check(UInt8[2, 0, 0, 0, 0xBA, 0xDC, 0xCD, 0xAB])

	# test 6: CreateString

	b = FlatBuffers.Builder()
	FlatBuffers.createstring(b, "foo")
	check(UInt8[3, 0, 0, 0, 'f', 'o', 'o', 0]) # 0-terminated, no pad
	FlatBuffers.createstring(b, "moop")
	check(UInt8[4, 0, 0, 0, 'm', 'o', 'o', 'p', 0, 0, 0, 0, # 0-terminated, 3-byte pad
		3, 0, 0, 0, 'f', 'o', 'o', 0])

	# test 6b: CreateString unicode

	b = FlatBuffers.Builder()
	# These characters are chinese from blog.golang.org/strings
	# We use escape codes here so that editors without unicode support
	# aren't bothered:
	uni_str = "\u65e5\u672c\u8a9e"
	FlatBuffers.createstring(b, uni_str)
	check(UInt8[9, 0, 0, 0, 230, 151, 165, 230, 156, 172, 232, 170, 158, 0, #  null-terminated, 2-byte pad
		0, 0])

	# test 6c: CreateUInt8String

	b = FlatBuffers.Builder()
	FlatBuffers.createstring(b, "foo")
	check(UInt8[3, 0, 0, 0, 'f', 'o', 'o', 0]) # 0-terminated, no pad
	FlatBuffers.createstring(b, "moop")
	check(UInt8[4, 0, 0, 0, 'm', 'o', 'o', 'p', 0, 0, 0, 0, # 0-terminated, 3-byte pad
		3, 0, 0, 0, 'f', 'o', 'o', 0])

	# test 7: empty vtable
	b = FlatBuffers.Builder()
    FlatBuffers.startobject(b, 0)
	check(UInt8[])
	FlatBuffers.endobject(b)
	check(UInt8[4, 0, 4, 0, 4, 0, 0, 0])

	# test 8: vtable with one true bool
	b = FlatBuffers.Builder()
	check(UInt8[])
    FlatBuffers.startobject(b, 1)
    check(UInt8[])
    FlatBuffers.prependslot!(b, 1, true, false)
    FlatBuffers.endobject(b)
	check(UInt8[
		6, 0, # vtable bytes
		8, 0, # length of object including vtable offset
		7, 0, # start of bool value
		6, 0, 0, 0, # offset for start of vtable (int32)
		0, 0, 0, # padded to 4 bytes
		1, # bool value
	])

	# test 9: vtable with one default bool
	b = FlatBuffers.Builder()
	check(UInt8[])
    FlatBuffers.startobject(b, 1)
    check(UInt8[])
	FlatBuffers.prependslot!(b, 1, false, false)
	FlatBuffers.endobject(b)
	check(UInt8[
		6, 0, # vtable bytes
		4, 0, # end of object from here
		0, 0, # entry 1 is zero
		6, 0, 0, 0, # offset for start of vtable (int32)
	])

	# test 10: vtable with one int16
	b = FlatBuffers.Builder()
    FlatBuffers.startobject(b, 1)
	FlatBuffers.prependslot!(b, 1, Int16(0x789A), Int16(0))
	FlatBuffers.endobject(b)
	check(UInt8[
		6, 0, # vtable bytes
		8, 0, # end of object from here
		6, 0, # offset to value
		6, 0, 0, 0, # offset for start of vtable (int32)
		0, 0, # padding to 4 bytes
		0x9A, 0x78,
	])

	# test 11: vtable with two int16
	b = FlatBuffers.Builder()
    FlatBuffers.startobject(b, 2)
	FlatBuffers.prependslot!(b, 1, Int16(0x3456), Int16(0))
	FlatBuffers.prependslot!(b, 2, Int16(0x789A), Int16(0))
	FlatBuffers.endobject(b)
	check(UInt8[
		8, 0, # vtable bytes
		8, 0, # end of object from here
		6, 0, # offset to value 0
		4, 0, # offset to value 1
		8, 0, 0, 0, # offset for start of vtable (int32)
		0x9A, 0x78, # value 1
		0x56, 0x34, # value 0
	])

	# test 12: vtable with int16 and bool
	b = FlatBuffers.Builder()
    FlatBuffers.startobject(b, 2)
	FlatBuffers.prependslot!(b, 1, Int16(0x3456), Int16(0))
	FlatBuffers.prependslot!(b, 2, true, false)
	FlatBuffers.endobject(b)
	check(UInt8[
		8, 0, # vtable bytes
		8, 0, # end of object from here
		6, 0, # offset to value 0
		5, 0, # offset to value 1
		8, 0, 0, 0, # offset for start of vtable (int32)
		0,          # padding
		1,          # value 1
		0x56, 0x34, # value 0
	])

	# test 12: vtable with empty vector
	b = FlatBuffers.Builder()
	FlatBuffers.startvector(b, sizeof(UInt8), 0, 1)
	vecend = FlatBuffers.endvector(b, 0)
    FlatBuffers.startobject(b, 1)
	FlatBuffers.prependoffsetslot!(b, 1, UInt32(vecend), UInt32(0))
	FlatBuffers.endobject(b)
	check(UInt8[
		6, 0, # vtable bytes
		8, 0,
		4, 0, # offset to vector offset
		6, 0, 0, 0, # offset for start of vtable (int32)
		4, 0, 0, 0,
		0, 0, 0, 0, # length of vector (not in struct)
	])

	# test 12b: vtable with empty vector of byte and some scalars
	b = FlatBuffers.Builder()
	FlatBuffers.startvector(b, sizeof(UInt8), 0, 1)
	vecend = FlatBuffers.endvector(b, 0)
    FlatBuffers.startobject(b, 2)
	FlatBuffers.prependslot!(b, 1, Int16(55), Int16(0))
	FlatBuffers.prependoffsetslot!(b, 2, UInt32(vecend), UInt32(0))
	FlatBuffers.endobject(b)
	check(UInt8[
		8, 0, # vtable bytes
		12, 0,
		10, 0, # offset to value 0
		4, 0, # offset to vector offset
		8, 0, 0, 0, # vtable loc
		8, 0, 0, 0, # value 1
		0, 0, 55, 0, # value 0

		0, 0, 0, 0, # length of vector (not in struct)
	])

	# test 13: vtable with 1 int16 and 2-vector of int16
	b = FlatBuffers.Builder()
	FlatBuffers.startvector(b, sizeof(Int16), 2, 1)
	FlatBuffers.prepend!(b, Int16(0x1234))
	FlatBuffers.prepend!(b, Int16(0x5678))
	vecend = FlatBuffers.endvector(b, 2)
    FlatBuffers.startobject(b, 2)
	FlatBuffers.prependoffsetslot!(b, 2, UInt32(vecend), UInt32(0))
	FlatBuffers.prependslot!(b, 1, Int16(55), Int16(0))
	FlatBuffers.endobject(b)
	check(UInt8[
		8, 0, # vtable bytes
		12, 0, # length of object
		6, 0, # start of value 0 from end of vtable
		8, 0, # start of value 1 from end of buffer
		8, 0, 0, 0, # offset for start of vtable (int32)
		0, 0, # padding
		55, 0, # value 0
		4, 0, 0, 0, # vector position from here
		2, 0, 0, 0, # length of vector (uint32)
		0x78, 0x56, # vector value 1
		0x34, 0x12, # vector value 0
	])

	# test 14: vtable with 1 struct of 1 int8, 1 int16, 1 int32
	b = FlatBuffers.Builder()
    FlatBuffers.startobject(b, 1)
	FlatBuffers.prep!(b, 4+4+4, 0)
	FlatBuffers.prepend!(b, Int8(55))
	FlatBuffers.pad!(b, 3)
	FlatBuffers.prepend!(b, Int16(0x1234))
	FlatBuffers.pad!(b, 2)
	FlatBuffers.prepend!(b, Int32(0x12345678))
	structStart = FlatBuffers.offset(b)
	FlatBuffers.prependstructslot!(b, 1, structStart, 0)
	FlatBuffers.endobject(b)
	check(UInt8[
		6, 0, # vtable bytes
		16, 0, # end of object from here
		4, 0, # start of struct from here
		6, 0, 0, 0, # offset for start of vtable (int32)
		0x78, 0x56, 0x34, 0x12, # value 2
		0, 0, # padding
		0x34, 0x12, # value 1
		0, 0, 0, # padding
		55, # value 0
	])

	# test 15: vtable with 1 vector of 2 struct of 2 int8
	b = FlatBuffers.Builder()
	FlatBuffers.startvector(b, sizeof(Int8)*2, 2, 1)
	FlatBuffers.prepend!(b, Int8(33))
	FlatBuffers.prepend!(b, Int8(44))
	FlatBuffers.prepend!(b, Int8(55))
	FlatBuffers.prepend!(b, Int8(66))
	vecend = FlatBuffers.endvector(b, 2)
    FlatBuffers.startobject(b, 1)
	FlatBuffers.prependoffsetslot!(b, 1, UInt32(vecend), UInt32(0))
	FlatBuffers.endobject(b)
	check(UInt8[
		6, 0, # vtable bytes
		8, 0,
		4, 0, # offset of vector offset
		6, 0, 0, 0, # offset for start of vtable (int32)
		4, 0, 0, 0, # vector start offset

		2, 0, 0, 0, # vector length
		66, # vector value 1,1
		55, # vector value 1,0
		44, # vector value 0,1
		33, # vector value 0,0
	])

	# test 16: table with some elements
	b = FlatBuffers.Builder()
    FlatBuffers.startobject(b, 2)
	FlatBuffers.prependslot!(b, 1, Int8(33), Int8(0))
	FlatBuffers.prependslot!(b, 2, Int16(66), Int16(0))
	off = FlatBuffers.endobject(b)
	FlatBuffers.finish!(b, off) #TODO

	check(UInt8[
		12, 0, 0, 0, # root of table: points to vtable offset

		8, 0, # vtable bytes
		8, 0, # end of object from here
		7, 0, # start of value 0
		4, 0, # start of value 1

		8, 0, 0, 0, # offset for start of vtable (int32)

		66, 0, # value 1
		0,  # padding
		33, # value 0
	])

	# test 17: one unfinished table and one finished table
	b = FlatBuffers.Builder()
    FlatBuffers.startobject(b, 2)
	FlatBuffers.prependslot!(b, 1, Int8(33), Int8(0))
	FlatBuffers.prependslot!(b, 2, Int8(44), Int8(0))
	off = FlatBuffers.endobject(b)
	FlatBuffers.finish!(b, off)

	FlatBuffers.startobject(b, 3)
	FlatBuffers.prependslot!(b, 1, Int8(55), Int8(0))
	FlatBuffers.prependslot!(b, 2, Int8(66), Int8(0))
	FlatBuffers.prependslot!(b, 3, Int8(77), Int8(0))
	off = FlatBuffers.endobject(b)
	FlatBuffers.finish!(b, off)

	check(UInt8[
		16, 0, 0, 0, # root of table: points to object
		0, 0, # padding

		10, 0, # vtable bytes
		8, 0, # size of object
		7, 0, # start of value 0
		6, 0, # start of value 1
		5, 0, # start of value 2
		10, 0, 0, 0, # offset for start of vtable (int32)
		0,  # padding
		77, # value 2
		66, # value 1
		55, # value 0

		12, 0, 0, 0, # root of table: points to object

		8, 0, # vtable bytes
		8, 0, # size of object
		7, 0, # start of value 0
		6, 0, # start of value 1
		8, 0, 0, 0, # offset for start of vtable (int32)
		0, 0, # padding
		44, # value 1
		33, # value 0
	])

	# test 18: a bunch of bools
	b = FlatBuffers.Builder()
    FlatBuffers.startobject(b, 8)
	FlatBuffers.prependslot!(b, 1, true, false)
	FlatBuffers.prependslot!(b, 2, true, false)
	FlatBuffers.prependslot!(b, 3, true, false)
	FlatBuffers.prependslot!(b, 4, true, false)
	FlatBuffers.prependslot!(b, 5, true, false)
	FlatBuffers.prependslot!(b, 6, true, false)
	FlatBuffers.prependslot!(b, 7, true, false)
	FlatBuffers.prependslot!(b, 8, true, false)
	off = FlatBuffers.endobject(b)
	FlatBuffers.finish!(b, off)

	check(UInt8[
		24, 0, 0, 0, # root of table: points to vtable offset

		20, 0, # vtable bytes
		12, 0, # size of object
		11, 0, # start of value 0
		10, 0, # start of value 1
		9, 0, # start of value 2
		8, 0, # start of value 3
		7, 0, # start of value 4
		6, 0, # start of value 5
		5, 0, # start of value 6
		4, 0, # start of value 7
		20, 0, 0, 0, # vtable offset

		1, # value 7
		1, # value 6
		1, # value 5
		1, # value 4
		1, # value 3
		1, # value 2
		1, # value 1
		1, # value 0
	])

	# test 19: three bools
	b = FlatBuffers.Builder()
    FlatBuffers.startobject(b, 3)
	FlatBuffers.prependslot!(b, 1, true, false)
	FlatBuffers.prependslot!(b, 2, true, false)
	FlatBuffers.prependslot!(b, 3, true, false)
	off = FlatBuffers.endobject(b)
	FlatBuffers.finish!(b, off)

	check(UInt8[
		16, 0, 0, 0, # root of table: points to vtable offset

		0, 0, # padding

		10, 0, # vtable bytes
		8, 0, # size of object
		7, 0, # start of value 0
		6, 0, # start of value 1
		5, 0, # start of value 2
		10, 0, 0, 0, # vtable offset from here

		0, # padding
		1, # value 2
		1, # value 1
		1, # value 0
	])

	# test 20: some floats
	b = FlatBuffers.Builder()
    FlatBuffers.startobject(b, 1)
	FlatBuffers.prependslot!(b, 1, Float32(1.0), Float32(0.0))
	off = FlatBuffers.endobject(b)

	check(UInt8[
		6, 0, # vtable bytes
		8, 0, # size of object
		4, 0, # start of value 0
		6, 0, 0, 0, # vtable offset

		0, 0, 128, 63, # value 0
	])
end

# CheckManualBuild builds a Monster manually.
function CheckManualBuild()
	b = FlatBuffers.Builder()
	str = FlatBuffers.createstring(b, "MyMonster")

	FlatBuffers.startvector(b, 1, 5, 1)
	FlatBuffers.prepend!(b, UInt8(4))
	FlatBuffers.prepend!(b, UInt8(3))
	FlatBuffers.prepend!(b, UInt8(2))
	FlatBuffers.prepend!(b, UInt8(1))
	FlatBuffers.prepend!(b, UInt8(0))
	inv = FlatBuffers.endvector(b, 5)

	FlatBuffers.startobject(b, 13)
	FlatBuffers.prependslot!(b, 2, Int16(20), Int16(100))
	mon2 = FlatBuffers.endobject(b)

	# Test4Vector
	FlatBuffers.startvector(b, 4, 2, 1)

	# Test 0
	FlatBuffers.prep!(b, 2, 4)
	FlatBuffers.pad!(b, 1)
	FlatBuffers.place!(b, Int8(20))
	FlatBuffers.place!(b, Int16(10))

	# Test 1
	FlatBuffers.prep!(b, 2, 4)
	FlatBuffers.pad!(b, 1)
	FlatBuffers.place!(b, Int8(40))
	FlatBuffers.place!(b, Int16(30))

	# end testvector
	test4 = FlatBuffers.endvector(b, 2)

	FlatBuffers.startobject(b, 13)

	# a vec3
	FlatBuffers.prep!(b, 16, 32)
	FlatBuffers.pad!(b, 2)
	FlatBuffers.prep!(b, 2, 4)
	FlatBuffers.pad!(b, 1)
	FlatBuffers.place!(b, UInt8(6))
	FlatBuffers.place!(b, Int16(5))
	FlatBuffers.pad!(b, 1)
	FlatBuffers.place!(b, UInt8(4))
	FlatBuffers.place!(b, Float64(3.0))
	FlatBuffers.pad!(b, 4)
	FlatBuffers.place!(b, Float32(3.0))
	FlatBuffers.place!(b, Float32(2.0))
	FlatBuffers.place!(b, Float32(1.0))
	vec3Loc = FlatBuffers.offset(b)
	# end vec3

	FlatBuffers.prependstructslot!(b, 1, vec3Loc, 0) # vec3. noop
	FlatBuffers.prependslot!(b, 3, Int16(80), Int16(100))     # hp
	FlatBuffers.prependoffsetslot!(b, 4, Int32(str), Int32(0))
	FlatBuffers.prependoffsetslot!(b, 6, Int32(inv), Int32(0)) # inventory
	FlatBuffers.prependslot!(b, 8, UInt8(1), UInt8(0))
	FlatBuffers.prependoffsetslot!(b, 9, Int32(mon2), Int32(0))
	FlatBuffers.prependoffsetslot!(b, 10, Int32(test4), Int32(0))
	mon = FlatBuffers.endobject(b)

	FlatBuffers.finish!(b, mon)

	return b.bytes, b.head
end

# CheckVtableDeduplication verifies that vtables are deduplicated.
function CheckVtableDeduplication()
	b = FlatBuffers.Builder()

	FlatBuffers.startobject(b, 4)
	FlatBuffers.prependslot!(b, 1, UInt8(0), UInt8(0))
	FlatBuffers.prependslot!(b, 2, UInt8(11), UInt8(0))
	FlatBuffers.prependslot!(b, 3, UInt8(22), UInt8(0))
	FlatBuffers.prependslot!(b, 4, Int16(33), Int16(0))
	obj0 = FlatBuffers.endobject(b)

	FlatBuffers.startobject(b, 4)
	FlatBuffers.prependslot!(b, 1, UInt8(0), UInt8(0))
	FlatBuffers.prependslot!(b, 2, UInt8(44), UInt8(0))
	FlatBuffers.prependslot!(b, 3, UInt8(55), UInt8(0))
	FlatBuffers.prependslot!(b, 4, Int16(66), Int16(0))
	obj1 = FlatBuffers.endobject(b)

	FlatBuffers.startobject(b, 4)
	FlatBuffers.prependslot!(b, 1, UInt8(0), UInt8(0))
	FlatBuffers.prependslot!(b, 2, UInt8(77), UInt8(0))
	FlatBuffers.prependslot!(b, 3, UInt8(88), UInt8(0))
	FlatBuffers.prependslot!(b, 4, Int16(99), Int16(0))
	obj2 = FlatBuffers.endobject(b)

	got = b.bytes[b.head+1:end]

	want = UInt8[
		240, 255, 255, 255, # == -12. offset to dedupped vtable.
		99, 0,
		88,
		77,
		248, 255, 255, 255, # == -8. offset to dedupped vtable.
		66, 0,
		55,
		44,
		12, 0, # start of vtable
		8, 0,
		0, 0,
		7, 0,
		6, 0,
		4, 0,
		12, 0, 0, 0, # table0
		33, 0,
		22,
		11
	]

    @test got == want

	table0 = FlatBuffers.Table{Any}(b.bytes, length(b.bytes) - obj0)
	table1 = FlatBuffers.Table{Any}(b.bytes, length(b.bytes) - obj1)
	table2 = FlatBuffers.Table{Any}(b.bytes, length(b.bytes) - obj2)

	function testTable(tab, a, b, c, d)
		# vtable size
        got = FlatBuffers.getoffsetslot(tab, 0, Int16(0))
		@test 12 == got
		# object size
        got = FlatBuffers.getoffsetslot(tab, 2, Int16(0))
		@test 8 == got
		# default value
        got = FlatBuffers.getoffsetslot(tab, 4, Int16(0))
		@test a == got
        got = FlatBuffers.getslot(tab, 6, UInt8(0))
		@test b == got
        val = FlatBuffers.getslot(tab, 8, UInt8(0))
		c != val && throw(ArgumentError("failed 8, 0: $got"))
        got = FlatBuffers.getslot(tab, 10, UInt8(0))
		@test d == got
        return
	end

	testTable(table0, UInt16(0), UInt8(11), UInt8(22), UInt8(33))
	testTable(table1, UInt16(0), UInt8(44), UInt8(55), UInt8(66))
	testTable(table2, UInt16(0), UInt8(77), UInt8(88), UInt8(99))
end

# CheckNotInObjectError verifies that `endobject` fails if not inside an
# object.
function CheckNotInObjectError()
	b = FlatBuffers.Builder()

	@test_throws ArgumentError FlatBuffers.endobject(b)
end

# CheckStringIsNestedError verifies that a string can not be created inside
# another object.
function CheckStringIsNestedError()
	b = FlatBuffers.Builder()
	FlatBuffers.startobject(b, 0)
	@test_throws ArgumentError FlatBuffers.createstring(b, "foo")
end

# CheckByteStringIsNestedError verifies that a bytestring can not be created
# inside another object.
function CheckByteStringIsNestedError()
	b = FlatBuffers.Builder()
	FlatBuffers.startobject(b, 0)
	@test_throws ArgumentError FlatBuffers.createstring(b, "foo")
end

# CheckStructIsNotInlineError verifies that writing a struct in a location
# away from where it is used will cause a panic.
function CheckStructIsNotInlineError()
	b = FlatBuffers.Builder()
	FlatBuffers.startobject(b, 0)
	@test_throws ArgumentError FlatBuffers.prependstructslot!(b, 0, 1, 0)
end

# CheckFinishedBytesError verifies that `FinishedBytes` panics if the table
# is not finished.
function CheckFinishedBytesError()
	b = FlatBuffers.Builder()

	@test_throws ArgumentError FlatBuffers.finishedbytes(b)
end

function CheckCreateByteVector()
	raw = UInt8(0):UInt8(29)

    for size = 1:30
		b1 = FlatBuffers.Builder()
		b2 = FlatBuffers.Builder()
		FlatBuffers.startvector(b1, 1, size, 1)
        for i = size:-1:1
			FlatBuffers.prepend!(b1, raw[i])
		end
		FlatBuffers.endvector(b1, size)
		FlatBuffers.createbytevector(b2, raw[1:size])
        @test b1.bytes == b2.bytes
	end
end

const InitialLCGSeed = 48271
type LCG
    val::UInt32
    LCG() = new(UInt32(InitialLCGSeed))
end
reset!(lcg::LCG) = lcg.val = UInt32(InitialLCGSeed)
function Base.next(lcg)
    n = UInt32((UInt64(lcg.val) * UInt64(279470273)) % UInt64(4294967291))
    lcg.val = n
    return n
end

# Low level stress/fuzz test: serialize/deserialize a variety of
# different kinds of data in different combinations
function checkFuzz(fuzzFields, fuzzObjects, verbose=true)

	# Values we're testing against: chosen to ensure no bits get chopped
	# off anywhere, and also be different from eachother.
	boolVal = true
	int8Val = Int8(-127) # 0x81
	uint8Val = UInt8(0xFF)
	int16Val = Int16(-32222) # 0x8222
	uint16Val = UInt16(0xFEEE)
	int32Val = Int32(overflowingInt32Val)
	uint32Val = UInt32(0xFDDDDDDD)
	int64Val = Int64(overflowingInt64Val)
	uint64Val = UInt64(0xFCCCCCCCCCCCCCCC)
	float32Val = Float32(3.14159)
	float64Val = Float64(3.14159265359)

	testValuesMax = 11 # hardcoded to the number of scalar types

	b = FlatBuffers.Builder()
	l = LCG()

    objects = Vector{Int}(fuzzObjects)

	# Generate fuzzObjects random objects each consisting of
	# fuzzFields fields, each of a random type.
    for i = 1:fuzzObjects
		FlatBuffers.startobject(b, fuzzFields)

        for f = 1:fuzzFields
			choice = next(l) % UInt32(testValuesMax)
            if choice ==  0
				FlatBuffers.prependslot!(b, f, boolVal, false)
			elseif choice ==  1
				FlatBuffers.prependslot!(b, f, int8Val, Int8(0))
			elseif choice ==  2
				FlatBuffers.prependslot!(b, f, uint8Val, UInt8(0))
			elseif choice ==  3
				FlatBuffers.prependslot!(b, f, int16Val, Int16(0))
			elseif choice ==  4
				FlatBuffers.prependslot!(b, f, uint16Val, UInt16(0))
			elseif choice ==  5
				FlatBuffers.prependslot!(b, f, int32Val, Int32(0))
			elseif choice ==  6
				FlatBuffers.prependslot!(b, f, uint32Val, UInt32(0))
			elseif choice ==  7
				FlatBuffers.prependslot!(b, f, int64Val, Int64(0))
			elseif choice ==  8
				FlatBuffers.prependslot!(b, f, uint64Val, UInt64(0))
			elseif choice ==  9
				FlatBuffers.prependslot!(b, f, float32Val, Float32(0))
			elseif choice ==  10
				FlatBuffers.prependslot!(b, f, float64Val, Float64(0))
			end
		end

		off = FlatBuffers.endobject(b)

		# store the offset from the end of the builder buffer,
		# since it will keep growing:
		objects[i] = off
	end

	# Do some bookkeeping to generate stats on fuzzes:
    stats = Dict{String,Int}()
    function check(desc, want, got)
        v = get!(stats, desc, 0)
        stats[desc] = v + 1
        @test want == got
    end

	l = LCG() # Reset.

	# Test that all objects we generated are readable and return the
	# expected values. We generate random objects in the same order
	# so this is deterministic.
    for i = 1:fuzzObjects

		table = FlatBuffers.Table{Any}(b.bytes, length(b.bytes) - objects[i])

        for j = 0:(fuzzFields - 1)
            f = (FlatBuffers.VtableMetadataFields + j) * sizeof(Int16)
			choice = next(l) % UInt32(testValuesMax)

            if choice == 0
				check("bool", boolVal, FlatBuffers.getslot(table, f, false))
			elseif choice == 1
				check("int8", int8Val, FlatBuffers.getslot(table, f, Int8(0)))
			elseif choice == 2
				check("uint8", uint8Val, FlatBuffers.getslot(table, f, UInt8(0)))
			elseif choice == 3
				check("int16", int16Val, FlatBuffers.getslot(table, f, Int16(0)))
			elseif choice == 4
				check("uint16", uint16Val, FlatBuffers.getslot(table, f, UInt16(0)))
			elseif choice == 5
				check("int32", int32Val, FlatBuffers.getslot(table, f, Int32(0)))
			elseif choice == 6
				check("uint32", uint32Val, FlatBuffers.getslot(table, f, UInt32(0)))
			elseif choice == 7
				check("int64", int64Val, FlatBuffers.getslot(table, f, Int64(0)))
			elseif choice == 8
				check("uint64", uint64Val, FlatBuffers.getslot(table, f, UInt64(0)))
			elseif choice == 9
				check("float32", float32Val, FlatBuffers.getslot(table, f, Float32(0)))
			elseif choice == 10
				check("float64", float64Val, FlatBuffers.getslot(table, f, Float64(0)))
			end
		end
	end

	# If enough checks were made, verify that all scalar types were used:
	if fuzzFields*fuzzObjects >= testValuesMax
		if length(stats) != testValuesMax
			throw(ArgumentError("fuzzing failed to test all scalar types"))
		end
	end

	# Print some counts, if needed:
    if verbose
		if fuzzFields == 0 || fuzzObjects == 0
            println("fuzz\tfields: $fuzzFields \tobjects: $fuzzObjects \t[none]\t 0")
		else
			ks = sort!(collect(keys(stats)))
            for k in ks
                println("fuzz\tfields: $fuzzFields \tobjects: $fuzzObjects \t$(k): $(stats[k])")
			end
		end
	end
	return
end


text:
#
## Sidef base
#

import Base.-,
       Base.+,
       Base.*,
       Base./,
       Base.//,
       Base.%,
       Base.^,
       Base.$,
       Base.&,
       Base.|,
       Base.<=,
       Base.>=,
       Base.<,
       Base.>,
       #Base.==,
       #Base.!=,
       Base.!,
       Base.abs,
       Base.min,
       Base.map,
       Base.sqrt,
       Base.floor,
       Base.print,
       Base.ifelse,
       Base.string,
       Base.println,
       Base.setindex!,
       Base.getindex;

abstract Sidef_Object
abstract Sidef_Types_Number_Complex

immutable Sidef_Types_Nil_Nil <: Sidef_Object
    value::Void
end

immutable Sidef_Types_Bool_Bool <: Sidef_Object
    value::Bool
end

immutable Sidef_Types_Block_Block <: Sidef_Object
    value::Function
end

immutable Sidef_Types_Number_Number <: Sidef_Object
    value::Number
end

immutable Sidef_Types_String_String <: Sidef_Object
    value::AbstractString
end

immutable Sidef_Types_Range_RangeNumber <: Sidef_Object
    value::Range
end

immutable Sidef_Types_Hash_Hash <: Sidef_Object
    value::Dict{Any,Any}
end

immutable Sidef_Types_Array_Array <: Sidef_Object
    value::Array{Any}
end

const NIL = Sidef_Types_Nil_Nil(Void())

const TRUE = Sidef_Types_Bool_Bool(true)
const FALSE = Sidef_Types_Bool_Bool(false)

const MONE = Sidef_Types_Number_Number(-1)
const ZERO = Sidef_Types_Number_Number(0)
const ONE = Sidef_Types_Number_Number(1)

function interpolate(::Type{Sidef_Types_String_String}, a::Any...)
    str = ""
    for item in a
        str *= string(item)
    end
    Sidef_Types_String_String(str)
end

include("object.jl")
include("bool.jl")
include("nil.jl")
include("number.jl")
include("hash.jl")
include("complex.jl")
include("block.jl")
include("range.jl")
include("string.jl")
include("array.jl")
include("convert.jl")

max_stars_repo_name: trizen/Sidef-Deparse-Julia

max_issues_repo_name: trizen/Sidef-Deparse-Julia

max_forks_repo_name: trizen/Sidef-Deparse-Julia


text:
# Basic Operators
# ---------------

"""
Count how many nucleotides satisfy a condition (i.e. f(seq[i]) -> true).

The first argument should be a function which accepts a nucleotide as its parameter.
"""
function Base.count(f::Function, seq::BioSequence)
    n = 0
    for x in seq
        if f(x)
            n += 1
        end
    end
    return n
end

# Mismatch counting
# -----------------

"""
    mismatches(seq1::BioSequence, seq2::BioSequence[, compatible=false])

Return the number of mismatches between `seq1` and `seq2`.

If `seq1` and `seq2` are of differing lengths, only the first `min(length(seq1),
length(seq2))` nucleotides are compared.  When `compatible` is `true`, sequence
symbols are comapred using `iscompatible`; otherwise using `==`.
"""
function mismatches{A<:Alphabet}(
        seq1::BioSequence{A},
        seq2::BioSequence{A},
        compatible::Bool=false)
    if ((bitsof(A) == 2 || bitsof(A) == 4) && !compatible) ||
        A == DNAAlphabet{2} ||
        A == RNAAlphabet{2}
        return bitparallel_mismatches(seq1, seq2)
    end

    mis = 0
    if compatible
        for (x, y) in zip(seq1, seq2)
            if !iscompatible(x, y)
                mis += 1
            end
        end
    else
        for (x, y) in zip(seq1, seq2)
            if x != y
                mis += 1
            end
        end
    end
    return mis
end

@generated function bitparallel_mismatches{A}(a::BioSequence{A}, b::BioSequence{A})
    n = bitsof(A)
    if n == 2
        bitpar_mismatches = :bitpar_mismatches2
    elseif n == 4
        bitpar_mismatches = :bitpar_mismatches4
    else
        error("n (= $n) ∉ (2, 4)")
    end

    quote
        if length(a) > length(b)
            return bitparallel_mismatches(b, a)
        end
        @assert length(a) ≤ length(b)

        nexta = bitindex(a, 1)
        nextb = bitindex(b, 1)
        stopa = bitindex(a, endof(a) + 1)
        mismatches = 0

        # align reading position of `a.data` so that `offset(nexta) == 0`
        if nexta < stopa && offset(nexta) != 0
            x = a.data[index(nexta)] >> offset(nexta)
            y = b.data[index(nextb)] >> offset(nextb)
            if offset(nextb) > offset(nexta)
                y |= b.data[index(nextb)+1] << (64 - offset(nextb))
            end
            k = 64 - offset(nexta)
            m = mask(k)
            mismatches += $bitpar_mismatches(x & m, y & m)
            nexta += k
            nextb += k
        end
        @assert offset(nexta) == 0

        if offset(nextb) == 0  # data are aligned with each other
            while stopa - nexta ≥ 64
                x = a.data[index(nexta)]
                y = b.data[index(nextb)]
                mismatches += $bitpar_mismatches(x, y)
                nexta += 64
                nextb += 64
            end

            if nexta < stopa
                x = a.data[index(nexta)]
                y = b.data[index(nextb)]
                m = mask(stopa - nexta)
                mismatches += $bitpar_mismatches(x & m, y & m)
            end
        elseif nexta < stopa
            y = b.data[index(nextb)]
            nextb += 64

            while stopa - nexta ≥ 64
                x = a.data[index(nexta)]
                z = b.data[index(nextb)]
                y = y >> offset(nextb) | z << (64 - offset(nextb))
                mismatches += $bitpar_mismatches(x, y)
                y = z
                nexta += 64
                nextb += 64
            end

            if nexta < stopa
                x = a.data[index(nexta)]
                y = y >> offset(nextb)
                if 64 - offset(nextb) < stopa - nexta
                    y |= a.data[index(nextb)] << (64 - offset(nextb))
                end
                m = mask(stopa - nexta)
                mismatches += $bitpar_mismatches(x & m, y & m)
            end
        end

        return mismatches
    end
end

# bit-parallel mismatch count algorithm for 2 and 4-bit encoding
@inline function bitpar_mismatches2(x::UInt64, y::UInt64)
    xyxor = x $ y
    mismatches = UInt64(0)
    mismatches |=  xyxor & 0x5555555555555555
    mismatches |= (xyxor & 0xAAAAAAAAAAAAAAAA) >> 1
    return count_ones(mismatches)
end

@inline function bitpar_mismatches4(x::UInt64, y::UInt64)
    xyxor = x $ y
    mismatches = UInt64(0)
    mismatches |=  xyxor & 0x1111111111111111
    mismatches |= (xyxor & 0x2222222222222222) >> 1
    mismatches |= (xyxor & 0x4444444444444444) >> 2
    mismatches |= (xyxor & 0x8888888888888888) >> 3
    return count_ones(mismatches)
end

# Sequences to Matrix
# -------------------

"""
    seqmatrix{A<:Alphabet}(vseq::Vector{BioSequence{A}}, major::Symbol)

Construct a matrix of nucleotides or amino acids from a vector of `BioSequence`s.

If parameter `major` is set to `:site`, the matrix is created such that one
nucleotide from each sequence is placed in each column i.e. the matrix is laid
out in site-major order.
This means that iteration over one position of many sequences is efficient,
as julia arrays are laid out in column major order.

If the parameter `major` is set to `:seq`, the matrix is created such that each
sequence is placed in one column i.e. the matrix is laid out in sequence-major
order.
This means that iteration across each sequence in turn is efficient,
as julia arrays are laid out in column major order.

# Examples
```julia
julia> seqs = [dna"AAA", dna"TTT", dna"CCC", dna"GGG"]
4-element Array{Bio.Seq.BioSequence{Bio.Seq.DNAAlphabet{4}},1}:
 3nt DNA Sequence:
AAA
 3nt DNA Sequence:
TTT
 3nt DNA Sequence:
CCC
 3nt DNA Sequence:
GGG

julia> seqmatrix(seqs, :site)
4x3 Array{Bio.Seq.DNA,2}:
 DNA_A  DNA_A  DNA_A
 DNA_T  DNA_T  DNA_T
 DNA_C  DNA_C  DNA_C
 DNA_G  DNA_G  DNA_G

 julia> seqmatrix(seqs, :seq)
 3x4 Array{Bio.Seq.DNA,2}:
  DNA_A  DNA_T  DNA_C  DNA_G
  DNA_A  DNA_T  DNA_C  DNA_G
  DNA_A  DNA_T  DNA_C  DNA_G
```
"""
function seqmatrix{A<:Alphabet}(vseq::AbstractVector{BioSequence{A}}, major::Symbol)
    nseqs = length(vseq)
    @assert nseqs > 0 throw(ArgumentError("Vector of BioSequence{$A} is empty."))
    nsites = length(vseq[1])
    @inbounds for i in 2:nseqs
        length(vseq[i]) == nsites || throw(ArgumentError("Sequences in vseq must be of same length"))
    end
    if major == :site
        mat = Matrix{eltype(A)}(nseqs, nsites)
        @inbounds for seq in 1:nseqs, site in 1:nsites
            mat[seq, site] = vseq[seq][site]
        end
        return mat
    elseif major == :seq
        mat = Matrix{eltype(A)}(nsites, nseqs)
        @inbounds for seq in 1:nseqs, site in 1:nsites
            mat[site, seq] = vseq[seq][site]
        end
        return mat
    else
        throw(ArgumentError("major must be :site or :seq"))
    end
end

"""
    seqmatrix{A<:Alphabet,T}(::Type{T}, vseq::Vector{BioSequence{A}}, major::Symbol)

Construct a matrix of `T` from a vector of `BioSequence`s.

If parameter `major` is set to `:site`, the matrix is created such that one
nucleotide from each sequence is placed in each column i.e. the matrix is laid
out in site-major order.
This means that iteration over one position of many sequences is efficient,
as julia arrays are laid out in column major order.

If the parameter `major` is set to `:seq`, the matrix is created such that each
sequence is placed in one column i.e. the matrix is laid out in sequence-major
order.
This means that iteration across each sequence in turn is efficient,
as julia arrays are laid out in column major order.

# Examples
```julia
julia> seqs = [dna"AAA", dna"TTT", dna"CCC", dna"GGG"]
4-element Array{Bio.Seq.BioSequence{Bio.Seq.DNAAlphabet{4}},1}:
 3nt DNA Sequence:
AAA
 3nt DNA Sequence:
TTT
 3nt DNA Sequence:
CCC
 3nt DNA Sequence:
GGG

julia> seqmatrix(seqs, :site, UInt8)
4×3 Array{UInt8,2}:
 0x01  0x01  0x01
 0x08  0x08  0x08
 0x02  0x02  0x02
 0x04  0x04  0x04

julia> seqmatrix(seqs, :seq, UInt8)
3×4 Array{UInt8,2}:
 0x01  0x08  0x02  0x04
 0x01  0x08  0x02  0x04
 0x01  0x08  0x02  0x04
```
"""
function seqmatrix{T,A<:Alphabet}(::Type{T}, vseq::AbstractVector{BioSequence{A}}, major::Symbol)
    nseqs = length(vseq)
    @assert nseqs > 0 throw(ArgumentError("Vector of BioSequence{$A} is empty."))
    nsites = length(vseq[1])
    @inbounds for i in 2:nseqs
        length(vseq[i]) == nsites || throw(ArgumentError("Sequences in vseq must be of same length."))
    end
    if major == :site
        mat = Matrix{T}(nseqs, nsites)
        @inbounds for seq in 1:nseqs, site in 1:nsites
            mat[seq, site] = convert(T, vseq[seq][site])
        end
        return mat
    elseif major == :seq
        mat = Matrix{T}(nsites, nseqs)
        @inbounds for seq in 1:nseqs, site in 1:nsites
            mat[site, seq] = convert(T, vseq[seq][site])
        end
        return mat
    else
        throw(ArgumentError("major must be :site or :seq"))
    end
end

# Consensus
# ---------

"""
    majorityvote{A<:NucleicAcidAlphabet}(seqs::AbstractVector{BioSequence{A}})

Construct a sequence that is a consensus of a vector of sequences.

The consensus is established by a simple majority vote rule, where amiguous
nucleotides cast an equal vote for each of their possible states.
For each site a winner(s) out of A, T(U), C, or G is determined, in the cases
of ties the ambiguity symbol that unifies all the winners is returned.
E.g if A and T tie, then W is inserted in the consensus. If all A, T, C, and G
tie at a site, then N is inserted in the consensus. Note this means that if a
nucletide e.g. 'C' and a gap '-' draw, the nucleotide will always win over the
gap, even though they tied.

# Examples
```julia
julia> seqs = [dna"CTCGATCGATCC", dna"CTCGAAAAATCA", dna"ATCGAAAAATCG", dna"ATCGGGGGATCG"]

4-element Array{Bio.Seq.BioSequence{Bio.Seq.DNAAlphabet{4}},1}:
 CTCGATCGATCC
 CTCGAAAAATCA
 ATCGAAAAATCG
 ATCGGGGGATCG

julia> majorityvote(seqs)
12nt DNA Sequence:
MTCGAAARATCG
```
"""
function majorityvote{A<:NucleicAcidAlphabet}(seqs::AbstractVector{BioSequence{A}})
    mat = seqmatrix(UInt8, seqs, :site)
    nsites = size(mat, 2)
    nseqs = size(mat, 1)
    result = BioSequence{A}(nsites)
    votes = Array{Int}(16)
    @inbounds for site in 1:nsites
        fill!(votes, 0)
        for seq in 1:nseqs
            nuc = mat[seq, site]
            votes[1] += nuc == 0x00
            votes[2] += (nuc & 0x01) != 0x00
            votes[3] += (nuc & 0x02) != 0x00
            votes[5] += (nuc & 0x04) != 0x00
            votes[9] += (nuc & 0x08) != 0x00
        end
        m = maximum(votes)
        merged = 0x00
        for i in 0x01:0x10
            merged |= ifelse(votes[i] == m, i - 0x01, 0x00)
        end
        result[site] = reinterpret(eltype(A), merged)
    end
    return result
end


text:
abstract type Planet <: CelestialBody end

const PLANET_NAMES = (
    "Mercury",
    "Venus",
    "Earth",
    "Mars",
    "Jupiter",
    "Saturn",
    "Uranus",
    "Neptune",
)

for (i, body) in enumerate(PLANET_NAMES)
    typ = Symbol(body)
    sym = Symbol(lowercase(body))
    id = 100i + 99
    typ_bc = Symbol(body, "Barycenter")
    nam_bc = string(Symbol(body, " Barycenter"))
    sym_bc = Symbol(lowercase(body), "_barycenter")
    add_edge!(BODIES, 0, i)
    add_edge!(BODIES, i, id)
    @eval begin
        struct $typ_bc <: Barycenter end
        const $sym_bc = $typ_bc()
        Base.show(io::IO, ::$typ_bc) = print(io, $nam_bc)
        parent(::$typ_bc) = ssb
        naifid(::$typ_bc) = $i
        from_naifid(::Val{$i}) = $sym_bc
        export $sym_bc, $typ_bc

        struct $typ <: Planet end
        const $sym = $typ()
        parent(::$typ) = $sym_bc
        naifid(::$typ) = $id
        from_naifid(::Val{$id}) = $sym
        export $sym, $typ
    end
end

# FIXME: Move this to a sensible place
export j2
j2(::Earth) = 1.08262668e-3


text:
function get_convergence_checker(obj::AbstractOptimizationProblem)
    return jcall(obj, "getConvergenceChecker", ConvergenceChecker, ())
end

function get_evaluation_counter(obj::AbstractOptimizationProblem)
    return jcall(obj, "getEvaluationCounter", Incrementor, ())
end

function get_iteration_counter(obj::AbstractOptimizationProblem)
    return jcall(obj, "getIterationCounter", Incrementor, ())
end


text:
# Test that planners correctly solve simple problems

@testset "Planners" begin

# Load domains and problems
path = joinpath(dirname(pathof(ParallelPlanners)), "..", "domains", "gridworld")
gridworld = load_domain(joinpath(path, "domain.pddl"))
gw_problem = load_problem(joinpath(path, "problem-1.pddl"))
gw_state = init_state(gw_problem)
gw_spec = Specification(gw_problem)

path = joinpath(dirname(pathof(ParallelPlanners)), "..", "domains", "doors-keys-gems")
doors_keys_gems = load_domain(joinpath(path, "domain.pddl"))
dkg_problem = load_problem(joinpath(path, "problem-1.pddl"))
dkg_state = init_state(dkg_problem)
dkg_spec = Specification(dkg_problem)

path = joinpath(dirname(pathof(ParallelPlanners)), "..", "domains", "blocksworld")
blocksworld = load_domain(joinpath(path, "domain.pddl"))
bw_problem = load_problem(joinpath(path, "problem-0.pddl"))
bw_state = init_state(bw_problem)
bw_spec = Specification(bw_problem)

@testset "Breadth-First Planner" begin

planner = BreadthFirstPlanner()
sol = planner(gridworld, gw_state, gw_spec)
@test is_goal(gw_spec, gridworld, sol.trajectory[end])
@test collect(sol) == @pddl("down", "down", "right", "right", "up", "up")

sol = planner(doors_keys_gems, dkg_state, dkg_spec)
@test is_goal(dkg_spec, doors_keys_gems, sol.trajectory[end])
@test collect(sol) == @pddl("(down)", "(pickup key1)", "(down)",
                            "(unlock key1 right)", "(right)", "(right)",
                            "(up)", "(up)", "(pickup gem1)")

sol = planner(blocksworld, bw_state, bw_spec)
@test is_goal(bw_spec, blocksworld, sol.trajectory[end])
@test collect(sol) == @pddl("(pick-up a)", "(stack a b)",
                            "(pick-up c)", "(stack c a)")

end

@testset "Uniform Cost Planner" begin

clear_available_action_cache!()

planner = UniformCostPlanner()
sol = planner(gridworld, gw_state, gw_spec)
@test is_goal(gw_spec, gridworld, sol.trajectory[end])
@test collect(sol) == @pddl("down", "down", "right", "right", "up", "up")

planner = UniformCostPlanner()
sol = planner(doors_keys_gems, dkg_state, dkg_spec)
@test is_goal(dkg_spec, doors_keys_gems, sol.trajectory[end])
@test collect(sol) == @pddl("(down)", "(pickup key1)", "(down)",
                            "(unlock key1 right)", "(right)", "(right)",
                            "(up)", "(up)", "(pickup gem1)")

planner = UniformCostPlanner()
sol = planner(blocksworld, bw_state, bw_spec)
@test is_goal(bw_spec, blocksworld, sol.trajectory[end])
@test collect(sol) == @pddl("(pick-up a)", "(stack a b)",
                            "(pick-up c)", "(stack c a)")

end

@testset "Greedy Planner" begin

clear_available_action_cache!()

planner = GreedyPlanner(ManhattanHeuristic(@pddl("xpos", "ypos")))
sol = planner(gridworld, gw_state, gw_spec)
@test is_goal(gw_spec, gridworld, sol.trajectory[end])
@test collect(sol) == @pddl("down", "down", "right", "right", "up", "up")

planner = GreedyPlanner(GoalCountHeuristic())
sol = planner(doors_keys_gems, dkg_state, dkg_spec)
@test is_goal(dkg_spec, doors_keys_gems, sol.trajectory[end])
@test collect(sol) == @pddl("(down)", "(pickup key1)", "(down)",
                            "(unlock key1 right)", "(right)", "(right)",
                            "(up)", "(up)", "(pickup gem1)")

planner = GreedyPlanner(HAdd())
sol = planner(blocksworld, bw_state, bw_spec)
@test is_goal(bw_spec, blocksworld, sol.trajectory[end])
@test collect(sol) == @pddl("(pick-up a)", "(stack a b)",
                            "(pick-up c)", "(stack c a)")

end

@testset "A* Planner" begin

clear_available_action_cache!()

planner = AStarPlanner(ManhattanHeuristic(@pddl("xpos", "ypos")))
sol = planner(gridworld, gw_state, gw_spec)
@test is_goal(gw_spec, gridworld, sol.trajectory[end])
@test collect(sol) == @pddl("down", "down", "right", "right", "up", "up")

planner = AStarPlanner(GoalCountHeuristic())
sol = planner(doors_keys_gems, dkg_state, dkg_spec)
@test is_goal(dkg_spec, doors_keys_gems, sol.trajectory[end])
@test collect(sol) == @pddl("(down)", "(pickup key1)", "(down)",
                            "(unlock key1 right)", "(right)", "(right)",
                            "(up)", "(up)", "(pickup gem1)")

planner = AStarPlanner(HAdd())
sol = planner(blocksworld, bw_state, bw_spec)
@test is_goal(bw_spec, blocksworld, sol.trajectory[end])
@test collect(sol) == @pddl("(pick-up a)", "(stack a b)",
                            "(pick-up c)", "(stack c a)")

end

@testset "Backward Planner" begin

clear_relevant_action_cache!()

planner = BackwardPlanner(heuristic=HAddR())
sol = planner(blocksworld, bw_state, bw_spec)
spec = ParallelPlanners.BackwardSearchGoal(bw_spec, bw_state)
@test is_goal(spec, blocksworld, sol.trajectory[1])
@test collect(sol) == @pddl("(pick-up a)", "(stack a b)",
                            "(pick-up c)", "(stack c a)")

end

@testset "Real Time Dynamic Programming" begin

clear_available_action_cache!()
Random.seed!(0)
simulator = StateActionRecorder(100)

heuristic = ManhattanHeuristic(@pddl("xpos", "ypos"))
planner = RTDP(heuristic=heuristic, rollout_noise=1.0, n_rollouts=10)
sol = planner(gridworld, gw_state, gw_spec)
actions, trajectory = simulator(sol, gridworld, gw_state, gw_spec)
@test is_goal(gw_spec, gridworld, trajectory[end])
@test actions == @pddl("down", "down", "right", "right", "up", "up")

heuristic = GoalCountHeuristic()
planner = RTDP(heuristic=heuristic, rollout_noise=1.0, n_rollouts=10)
sol = planner(doors_keys_gems, dkg_state, dkg_spec)
actions, trajectory = simulator(sol, doors_keys_gems, dkg_state, dkg_spec)
@test is_goal(dkg_spec, doors_keys_gems, trajectory[end])
@test actions == @pddl("(down)", "(pickup key1)", "(down)",
                       "(unlock key1 right)", "(right)", "(right)",
                        "(up)", "(up)", "(pickup gem1)")

planner = RTDP(heuristic=HAdd(), rollout_noise=1.0, n_rollouts=10)
sol = planner(blocksworld, bw_state, bw_spec)
actions, trajectory = simulator(sol, blocksworld, bw_state, bw_spec)
@test is_goal(bw_spec, blocksworld, trajectory[end])
@test actions == @pddl("(pick-up a)", "(stack a b)",
                       "(pick-up c)", "(stack c a)")

end

@testset "Monte Carlo Tree Search" begin

clear_available_action_cache!()
Random.seed!(0)
simulator = StateActionRecorder(100)

planner = MCTS(n_rollouts=50)
sol = planner(gridworld, gw_state, GoalReward(gw_problem.goal, 100))
actions, trajectory = simulator(sol, gridworld, gw_state, gw_spec)
@test is_goal(gw_spec, gridworld, trajectory[end])
@test actions == @pddl("down", "down", "right", "right", "up", "up")

planner = MCTS(n_rollouts=50)
sol = planner(doors_keys_gems, dkg_state, GoalReward(dkg_problem.goal, 100))
actions, trajectory = simulator(sol, doors_keys_gems, dkg_state, dkg_spec)
@test is_goal(dkg_spec, doors_keys_gems, trajectory[end])
@test actions == @pddl("(down)", "(pickup key1)", "(down)",
                       "(unlock key1 right)", "(right)", "(right)",
                        "(up)", "(up)", "(pickup gem1)")

planner = MCTS(n_rollouts=50)
sol = planner(blocksworld, bw_state, GoalReward(bw_problem.goal, 100))
actions, trajectory = simulator(sol, blocksworld, bw_state, bw_spec)
@test is_goal(bw_spec, blocksworld, trajectory[end])
@test actions == @pddl("(pick-up a)", "(stack a b)",
                       "(pick-up c)", "(stack c a)")

end

end


text:
module test_04_ConstantSignalsWithUnit

using ModiaResult
using Unitful
using OrderedCollections
ModiaResult.@usingModiaPlot

inertia = [1.1  1.2  1.3;
           2.1  2.2  2.3;
           3.1  3.2  3.3]u"kg*m^2"

result = OrderedDict{String,Any}()

result["time"]     = [0.0, 1.0]*u"s"
result["phi_max"]  = [1.1f0, 1.1f0]*u"rad"
result["i_max"]    = [2, 2]
result["open"]     = [true , true]
result["Inertia"]  = [inertia, inertia]

println("\n... test_04_ConstantSignalsWithUnit.jl:")
ModiaResult.printResultInfo(result)

plot(result, ["phi_max", "i_max", "open", "Inertia[2,2]", "Inertia[1,2:3]", "Inertia"], heading="Constants")

end

max_stars_repo_name: ModiaSim/ModiaResult

max_issues_repo_name: ModiaSim/ModiaResult

max_forks_repo_name: ModiaSim/ModiaResult


text:
using MarketData

facts("Visualizations") do

    context("Raw data visualizations") do
        plot(cl)
        plot(ohlc)
        histogram(cl)
        density(cl)
    end

    context("Serial correlation visualizations") do
        directedscatter(cl[1:5])
        @fact_throws directedscatter(cl)
        directedscatter(cl, naive=true)
        directedscatter(cl, naive=true, nlags=2)

        laghistogram(cl[1:5])
        @fact_throws laghistogram(cl)
        laghistogram(cl, naive=true)
        laghistogram(cl, naive=true, nlags=2)

        @fact_throws acf(cl[1:5], lags=1:10)
        acf(cl[1:5], lags=1:3)
        @fact_throws acf(cl)
        acf(cl, naive=true)
        acf(cl, naive=true, lags=1:5)

        @fact_throws pacf(cl[1:5], lags=1:10)
        pacf(cl[1:5], lags=1:2)
        @fact_throws pacf(cl)
        pacf(cl, naive=true)
        pacf(cl, naive=true, lags=1:10)

        @fact_throws ccf(ohlc["Open", "Close"][1:5])
        ccf(ohlc["Open", "Close"][1:5], lags=-1:1)
        @fact_throws ccf(ohlc["Open", "Close"])
        ccf(ohlc["Open", "Close"], naive=true)
        ccf(ohlc["Open", "Close"], naive=true, lags=-5:5)

    end

    context("Spectral analysis") do
        cumulative_periodogram(cl[1:5])
        @fact_throws cumulative_periodogram(cl)
        cumulative_periodogram(cl, naive=true)
    end

end


text:
using PyCall

@pyimport optunity
include("nogc.jl")

map_args = args -> map(e -> e[1], args)

function wrap_call(f::Function,all_args)
	wrapper_(;kw...) = f(Dict(kw))

	function wrapper(;kw...)
		for e in kw
			ex = e[1]; @eval $ex = $e[2]
		end
		args = intersect(all_args,map_args(kw))
		f(map(e -> (@eval $e), args)...)
	end

	try
		methods(f).defs.sig <: Tuple{Dict} ? wrapper_ : wrapper
	catch ArgumentError
		wrapper
	end
end

function minimize(f::Function;kws...)
	results = @nogc optunity.minimize(wrap_call(f,map_args(kws));kws...)
	results[1], results[2]
end

function maximize(f::Function;kws...)
	results = @nogc optunity.maximize(wrap_call(f,map_args(kws));kws...)
	results[1], results[2]
end

max_stars_repo_name: xrounder/optunity

max_issues_repo_name: xrounder/optunity

max_forks_repo_name: xrounder/optunity


text:
using HealthMLBase
using Test

@testset "HealthMLBase.jl" begin
end


text:
"""
REWRITE:
back(r::Net,ygold,loss) computes the gradients of weights and
activations according to the given loss function and gold output.
ygold represents an individual item minibatch that may or may not be
an element of a sequence.  The seq keyword argument determines which:
initback sets incr=true for par if seq, back pops from stack if seq.
The loss gradient of the output, ygrad, is computed using
loss(ypred,ygold,ygrad).  ypred is retrieved from r.out[N] where N is
the index of the last op.  ygrad is written to r.dif[N].  If r.op[N]
has multiple outputs (toincr[N]), r.dif[N] is incremented.  If the
optional loss argument is not provided, ygold is used as the loss
gradient.  If ygold=nothing means the loss gradient from the output is
taken to be 0.  Gradient computation proceeds backwards from N..1.

"""
function back(f::Net, ygold=nothing, loss=copyloss; seq=false, getdx=false, o...)
    getdx = getdxbool(getdx, ninputs(f)); dx=Any[]
    initback(f, ygold, loss, getdx, seq)
    gotreturn = false
    gotinputs = 0
    for n = length(f):-1:1
        skipped(f,n,seq) && continue
        y = reg(f,n)
        getp(y,:grad) || continue # :grad means we need the gradient of the output of this operation, does not mean go back on this operation. but no :grad means no need to go back.
        yout = (!seq ? y.out : stack_output(f,n))
        if !gotreturn
            gotreturn = true
            if ygold == nothing # represents zero gradient
                getp(y,:incr) || (y.dif = nothing)
            else
                !isreturn(y) && Base.warn_once("ygold specified when there is no return")
                if getp(y,:incr)
                    loss(yout, ygold, y.tmp; o...) # loss needs o... for e.g. mask
                    y.dif = axpy!(1,y.tmp,y.dif0)
                else
                    y.dif = loss(yout, ygold, y.dif0; o...)
                end
            end
        elseif isreturn(y)
            error("Got return in non-final instruction")
        end
        # println(:back, (findfirst(regs(f), y), typeof(y.op), y.argv, Knet.vecnorm0(y.dif)))
        xx = (!seq ? inputregs(f,y) : stack_inputregs(f,n))
        if isempty(xx)
            # nothing to do
        elseif y.dif == nothing
            for x in xx
                getp(x,:grad) && !getp(x,:incr) && (x.dif = nothing)
            end
        else
            xdif = difs(xx)
            xout = get1(!seq ? outs(xx) : stack_inputs(f,n))
            back(y.op, y.dif, xdif...; x=xout, y=yout, o...)
            for x in xx
                x.dif = (getp(x,:incr) ? axpy!(1, x.tmp, x.dif0) :
                         getp(x,:grad) ? x.dif0 : nothing)
            end
            if getp(y,:incr) && !isa(y.op, Par)
                fillsync!(y.dif,0)
            end
        end
        if isa(y.op,Input)
            getdx[end-gotinputs] && unshift!(dx, y.dif)
            gotinputs += 1
        end
    end
    seq && (f.sp -= length(f))  # this way f.stack[f.sp] gives the output of n'th op next call
    return get1(dx)
end

sback(f::Net, ygold=nothing, loss=copyloss; o...)=back(f,ygold,loss;o...,seq=true)

# For FNNs we can tell whether an op was skipped using the :forw flag.
# For RNNs the :forw flag may change through the iterations, so we check the stack.
skipped(f::Net,n::Int,seq::Bool)=(seq ? !stack_forw(f,n) : !getp(reg(f,n),:forw))

# We assume f.reg[end].out is saved in f.stack[f.sp].out
# In the calls below n is the network index
stack_entry(f::Net,n::Int)=f.stack[f.sp - length(f) + n]
stack_output(f::Net,n::Int)=stack_entry(f,n).out
stack_argv(f::Net,n::Int)=stack_entry(f,n).argv
stack_forw(f::Net,n::Int)=stack_entry(f,n).forw

# For inputs the situation is a bit more complicated.
# We want to use their values prior to n'th op.
# And we want to use the indices from the argv on stack.

stack_inputregs(f::Net,n::Int)=f.reg[stack_argv(f,n)]

function stack_inputs(f::Net,n::Int)
    r = reg(f,n)
    back_reads_x(r.op) || return nothing
    a = stack_argv(f,n)
    x = cell(length(a))
    @inbounds for i=1:length(a)
        stack_index = f.sp - length(f) + a[i]
        a[i] >= n && (stack_index -= length(f)) # use prior value of read-before-write register
        if stack_index <= 0
            x[i] = nothing
        else
            s = f.stack[stack_index]
            s.forw || Base.warn_once("Stack input has forw=false.")
            x[i] = s.out
        end
    end
    return x
end

# this is used when no loss fn specified, in which case we assume ygold is actually ygrad
copyloss(ypred,ygold,ygrad;o...)=(ygrad===ygold ? ygrad : copysync!(ygrad,ygold))

# turn various forms of getdx into boolean vector
function getdxbool(getdx, n)
    (isa(getdx, Vector{Bool}) && length(getdx)==n ? getdx :
     isa(getdx, Bool) ? fill(getdx, n) :
     isa(getdx, Vector{Int}) ? (tmp=falses(n);tmp[getdx]=true;tmp) :
     isa(getdx, Int) ? (tmp=falses(n);tmp[getdx]=true;tmp) :
     error("getdx=$getdx ninputs=$(n)"))
end

# This is extremely slow:
#get1(x)=(!isempty(methods(length, (typeof(x),))) && length(x)==1?x[1]:x)
get1(x)=(x==nothing?x:length(x)==1?x[1]:x)

### DEAD CODE:

# function sback_delete(f::Net, ygold=nothing, loss=copyloss; getdx=false, o...)
#     getdx = getdxbool(getdx, ninputs(f))
#     initback(f, ygold, loss, getdx) # TODO: rethink lastforw==lastback in a seq context
#     gotreturn = false

#     for n = length(f):-1:1
#         s = pop!(f)
#         s == nothing && continue
#         (y, xsave, ysave) = s
#         @assert y === reg(f,n)
#         getp(y,:grad) || continue # :grad means we need the gradient of the output of this operation, does not mean go back on this operation. but no :grad means no need to go back.

#         if !gotreturn
#             gotreturn = true
#             if ygold == nothing # represents zero gradient
#                 getp(y,:incr) || (y.dif = nothing)
#             else
#                 !isreturn(y) && Base.warn_once("ygold specified when there is no return")
#                 ysave == nothing && error("return value was not saved")
#                 if getp(y,:incr)
#                     loss(ysave, ygold, y.tmp; o...) # loss needs o... for e.g. mask
#                     y.dif = axpy!(1,y.tmp,y.dif0)
#                 else
#                     y.dif = loss(ysave, ygold, y.dif0; o...)
#                 end
#             end
#         elseif isreturn(y)
#             error("Got return in non-final instruction")
#         end

#         xx = inputregs(f,y)
#         if y.dif == nothing
#             for x in xx
#                 getp(x,:grad) && !getp(x,:incr) && (x.dif = nothing)
#             end
#         else
#             xdif = inputdifs(f,y)  # map too slow? map(x->(!getp(x,:grad) ? nothing : getp(x,:incr) ? x.tmp : x.dif0), xx)
#             back(y.op, y.dif, xdif...; x=get1(xsave), y=ysave, o...)
#             for x in xx
#                 x.dif = (getp(x,:incr) ? axpy!(1, x.tmp, x.dif0) :
#                          getp(x,:grad) ? x.dif0 : nothing)
#             end
#             if getp(y,:incr) && !isa(y.op, Par)
#                 # what if y.op=Arr?  then it will have no inputs, thus :grad=:incr=false
#                 # where does Par.dif get zeroed out? at reset!
#                 fillsync!(y.dif,0)
#             end
#         end
#     end

#     if any(getdx)
#         dx = Any[]; nx = 0
#         for p in regs(f)
#             isa(p.op,Input) && getdx[nx+=1] && push!(dx, p.dif)
#         end
#         return get1(dx)
#     end
# end

# # back(r::Net,dy::Vector) for a sequence
# function back(r::Net, dy::Vector, dx...; a...)
#     dxi = map(x->(x==nothing ? x : x[end]), dx)
#     initback(r, dy[end], dxi...; seq=true, a...)
#     for i=length(dy):-1:1
#         dxi = map(x->(x==nothing ? x : x[i]), dx)
#         back(r, dy[i], dxi...; seq=true, a...)
#     end
# end

# DONE: truncated bptt
# - go forward k1 steps, run back for k2, update, recover state
# - if k1==k2 we just need the keepstate option to forw
# - if k1>k2 the stack won't be cleared
# - if k1<k2 the stack will be overdrawn

    # for i = ninputs(r):-1:1
    #     n = i+N
    #     r.tosave[n] && pop(r,n)                                    # ; r.tosave[n] && dbg(r,:out,n)
    #     dx == nothing || copysync!(dx[i], r.dif[n])
    # end

    # if dx != ()
    #     lastinput = 0
    #     for n = 1:N
    #         isa(r.op[n], Input) || continue
    #         copysync!(dx[lastinput += 1], r.dif[n])
    #     end
    # end
# function stack_output(f::Net,n::Int)
#     o = f.stack[f.sp - length(f) + n]
#     p = f.reg[n]
#     (o!=nothing) || !(isreturn(p) || back_reads_y(p.op)) || error("sp=$(f.sp) n=$n N=$(length(f)) op=$(typeof(f.reg[n].op))")
#     return o
# end


text:
"""
    YASStyle()

Formatting style based on [YASGuide](https://github.com/jrevels/YASGuide)
and [JuliaFormatter#198](https://github.com/domluna/JuliaFormatter.jl/issues/198).

Configurable options with different defaults to [`DefaultStyle`](@ref) are:
- `always_for_in` = true
- `whitespace_ops_in_indices` = true
- `remove_extra_newlines` = true
- `import_to_using` = true
- `pipe_to_function_call` = true
- `short_to_long_function_def` = true
- `always_use_return` = true
- `whitespace_in_kwargs` = false
"""
struct YASStyle <: AbstractStyle end
@inline getstyle(s::YASStyle) = s

function options(style::YASStyle)
    return (;
        always_for_in = true,
        whitespace_ops_in_indices = true,
        remove_extra_newlines = true,
        import_to_using = true,
        pipe_to_function_call = true,
        short_to_long_function_def = true,
        always_use_return = true,
        whitespace_in_kwargs = false,
    )
end

function nestable(::YASStyle, cst::CSTParser.EXPR)
    (CSTParser.defines_function(cst) || nest_assignment(cst)) && return false
    (cst[2].kind === Tokens.PAIR_ARROW || cst[2].kind === Tokens.ANON_FUNC) && return false
    return true
end

function p_import(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = FST(cst, nspaces(s))
    add_node!(t, pretty(ys, cst[1], s), s)
    add_node!(t, Whitespace(1), s)

    for (i, a) in enumerate(cst.args[2:end])
        if CSTParser.is_comma(a)
            add_node!(t, pretty(ys, a, s), s, join_lines = true)
            add_node!(t, Placeholder(1), s)
        elseif CSTParser.is_colon(a)
            add_node!(t, pretty(ys, a, s), s, join_lines = true)
            add_node!(t, Whitespace(1), s)
        else
            add_node!(t, pretty(ys, a, s), s, join_lines = true)
        end
    end
    t
end
@inline p_using(ys::YASStyle, cst::CSTParser.EXPR, s::State) = p_import(ys, cst, s)
@inline p_export(ys::YASStyle, cst::CSTParser.EXPR, s::State) = p_import(ys, cst, s)

function p_curly(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = FST(cst, nspaces(s))
    for (i, a) in enumerate(cst)
        n = pretty(ys, a, s)
        if CSTParser.is_comma(a) && i == length(cst) - 1
            continue
        elseif CSTParser.is_comma(a) && i < length(cst) && !is_punc(cst[i+1])
            add_node!(t, n, s, join_lines = true)
            add_node!(t, Placeholder(0), s)
        else
            add_node!(t, n, s, join_lines = true)
        end
    end
    t
end
@inline p_braces(ys::YASStyle, cst::CSTParser.EXPR, s::State) = p_curly(ys, cst, s)

function p_tupleh(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = FST(cst, nspaces(s))
    for (i, a) in enumerate(cst)
        if CSTParser.is_comma(a) && i + 1 == length(cst)
            n = pretty(ys, a, s)
            if n_args(cst) == 1
                add_node!(t, n, s, join_lines = true)
            elseif !is_closer(cst[i+1])
                add_node!(t, n, s, join_lines = true)
                add_node!(t, Placeholder(1), s)
            end
        elseif CSTParser.is_comma(a) && i < length(cst) && !is_punc(cst[i+1])
            add_node!(t, pretty(ys, a, s), s, join_lines = true)
            add_node!(t, Placeholder(1), s)
        elseif a.typ === CSTParser.BinaryOpCall && a[2].kind === Tokens.EQ
            add_node!(t, p_kw(ys, a, s), s, join_lines = true)
        else
            add_node!(t, pretty(ys, a, s), s, join_lines = true)
        end
    end
    t
end

function p_tupleh(ys::YASStyle, nodes::Vector{CSTParser.EXPR}, s::State)
    t = FST(CSTParser.TupleH, nspaces(s))
    for (i, a) in enumerate(nodes)
        n = pretty(ys, a, s)
        if CSTParser.is_comma(a) && i + 1 == length(nodes)
            if n_args(nodes) == 1
                add_node!(t, n, s, join_lines = true)
            elseif !is_closer(nodes[i+1])
                add_node!(t, n, s, join_lines = true)
                add_node!(t, Placeholder(1), s)
            end
        elseif CSTParser.is_comma(a) && i < length(nodes) && !is_punc(nodes[i+1])
            add_node!(t, n, s, join_lines = true)
            add_node!(t, Placeholder(1), s)
        else
            add_node!(t, n, s, join_lines = true)
        end
    end
    t
end

# InvisBrackets
function p_invisbrackets(
    ys::YASStyle,
    cst::CSTParser.EXPR,
    s::State;
    nonest = false,
    nospace = false,
)
    t = p_invisbrackets(DefaultStyle(ys), cst, s, nonest = nonest, nospace = nospace)
    for (i, n) in enumerate(t.nodes)
        if n.typ === PLACEHOLDER
            t[i] = Whitespace(length(n))
        end
    end
    t
end

function p_call(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = FST(cst, nspaces(s))
    for (i, a) in enumerate(cst)
        n = pretty(ys, a, s)
        if CSTParser.is_comma(a) && i + 1 == length(cst)
            continue
        elseif CSTParser.is_comma(a) && i < length(cst) && !is_punc(cst[i+1])
            add_node!(t, n, s, join_lines = true)
            add_node!(t, Placeholder(1), s)
        else
            add_node!(t, n, s, join_lines = true)
        end
    end
    t
end
@inline p_vect(ys::YASStyle, cst::CSTParser.EXPR, s::State) = p_call(ys, cst, s)

function p_ref(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = FST(cst, nspaces(s))
    nospace = !s.opts.whitespace_ops_in_indices

    for (i, a) in enumerate(cst)
        if CSTParser.is_comma(a) && i < length(cst) && !is_punc(cst[i+1])
            add_node!(t, pretty(ys, a, s), s, join_lines = true)
            add_node!(t, Placeholder(1), s)
        elseif a.typ === CSTParser.BinaryOpCall ||
               a.typ === CSTParser.InvisBrackets ||
               a.typ === CSTParser.ChainOpCall ||
               a.typ === CSTParser.Comparison

            n = pretty(ys, a, s, nonest = true, nospace = nospace)
            add_node!(t, n, s, join_lines = true)
        else
            add_node!(t, pretty(ys, a, s), s, join_lines = true)
        end
    end
    t
end

function p_comprehension(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = FST(cst, nspaces(s))

    if is_block(cst[2]) || (cst[2].typ === CSTParser.Generator && is_block(cst[2][1]))
        t.nest_behavior = AlwaysNest
    end

    add_node!(t, pretty(ys, cst[1], s), s, join_lines = true)
    add_node!(t, pretty(ys, cst[2], s), s, join_lines = true)
    add_node!(t, pretty(ys, cst[3], s), s, join_lines = true)
    t
end

function p_typedcomprehension(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = FST(cst, nspaces(s))

    if is_block(cst[3]) || (cst[3].typ === CSTParser.Generator && is_block(cst[3][1]))
        t.nest_behavior = AlwaysNest
    end

    add_node!(t, pretty(ys, cst[1], s), s, join_lines = true)
    add_node!(t, pretty(ys, cst[2], s), s, join_lines = true)
    add_node!(t, pretty(ys, cst[3], s), s, join_lines = true)
    add_node!(t, pretty(ys, cst[4], s), s, join_lines = true)
    t
end

function p_macrocall(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    fst = p_macrocall(DefaultStyle(ys), cst, s)
    fst.typ == CSTParser.MacroCall || return fst
    is_closer(cst.args[end]) || return fst

    # remove initial and last placeholders
    # @call(PLACEHOLDER, args..., PLACEHOLDER) -> @call(args...)
    idx = findfirst(n -> n.typ === PLACEHOLDER, fst.nodes)
    idx !== nothing && (fst[idx] = Whitespace(0))
    idx = findlast(n -> n.typ === PLACEHOLDER, fst.nodes)
    idx !== nothing && (fst[idx] = Whitespace(0))

    return fst
end

function p_parameters(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = FST(cst, nspaces(s))
    for (i, a) in enumerate(cst)
        n = pretty(ys, a, s)
        if i == length(cst) && CSTParser.is_comma(a)
            # do nothing
        elseif CSTParser.is_comma(a) && i < length(cst) && !is_punc(cst[i+1])
            add_node!(t, n, s, join_lines = true)
            add_node!(t, Placeholder(1), s)
        else
            add_node!(t, n, s, join_lines = true)
        end
    end
    t
end

function p_whereopcall(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = FST(cst, nspaces(s))
    add_node!(t, pretty(ys, cst[1], s), s)
    add_node!(t, Whitespace(1), s)
    add_node!(t, pretty(ys, cst[2], s), s, join_lines = true)
    add_node!(t, Whitespace(1), s)

    add_braces =
        !CSTParser.is_lbrace(cst[3]) &&
        cst.parent.typ !== CSTParser.Curly &&
        cst[3].typ !== CSTParser.Curly &&
        cst[3].typ !== CSTParser.BracesCat

    if add_braces
        brace = FST(CSTParser.PUNCTUATION, -1, t.endline, t.endline, "{")
        add_node!(t, brace, s, join_lines = true)
    end

    for (i, a) in enumerate(cst.args[3:end])
        n =
            a.typ === CSTParser.BinaryOpCall ? pretty(ys, a, s, nospace = true) :
            pretty(ys, a, s)

        if CSTParser.is_comma(a) && i + 2 == length(cst) - 1
            continue
        elseif CSTParser.is_comma(a) && i + 2 < length(cst) && !is_punc(cst[i+3])
            add_node!(t, n, s, join_lines = true)
            add_node!(t, Placeholder(0), s)
        else
            add_node!(t, n, s, join_lines = true)
        end
    end

    if add_braces
        brace = FST(CSTParser.PUNCTUATION, -1, t.endline, t.endline, "}")
        add_node!(t, brace, s, join_lines = true)
    end
    t
end

function p_generator(ys::YASStyle, cst::CSTParser.EXPR, s::State)
    t = FST(cst, nspaces(s))

    for (i, a) in enumerate(cst)
        n = pretty(ys, a, s)
        if a.typ === CSTParser.KEYWORD
            incomp = parent_is(
                a,
                is_iterable,
                ignore = n -> is_gen(n) || n.typ === CSTParser.InvisBrackets,
            )

            if is_block(cst[i-1])
                add_node!(t, Newline(), s)
            elseif incomp
                add_node!(t, Placeholder(1), s)
            else
                add_node!(t, Whitespace(1), s)
            end
            add_node!(t, n, s, join_lines = true)
            add_node!(t, Whitespace(1), s)

            if a.kind === Tokens.FOR
                for j = i+1:length(cst)
                    eq_to_in_normalization!(cst[j], s.opts.always_for_in)
                end
            end

            if !is_gen(cst.args[i+1])
                tup = p_tupleh(ys, cst.args[i+1:length(cst)], s)
                add_node!(t, tup, s, join_lines = true)
                return t
            end
        elseif CSTParser.is_comma(a) && i < length(cst) && !is_punc(cst[i+1])
            add_node!(t, n, s, join_lines = true)
            add_node!(t, Placeholder(1), s)
        else
            add_node!(t, n, s, join_lines = true)
        end
    end

    t
end
@inline p_filter(ys::YASStyle, cst::CSTParser.EXPR, s::State) = p_generator(ys, cst, s)
@inline p_flatten(ys::YASStyle, cst::CSTParser.EXPR, s::State) = p_generator(ys, cst, s)


text:
using BenchmarkTools
include("../src/MFCC.jl")
using .MFCC

# @benchmark MFCC.mfcc(x) setup=(x=randn(10_000))
# 800 us

@benchmark MFCC.powspec(x) setup=(x=randn(10_000))
# 400 us

# @benchmark MFCC.audspec(x) setup=(x=MFCC.powspec(randn(10_000)))


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

using Markdown
using InteractiveUtils

# ╔═╡ 66923a50-0fa8-11eb-1547-cf8e289818b0
using Revise

# ╔═╡ 8336bf6e-0fa8-11eb-1f48-37cfe47a6d0f
using Pkg

# ╔═╡ 93125894-0fa8-11eb-0ff6-211feff97889
using PolaronMobility

# ╔═╡ f0529460-0fa8-11eb-1bec-1feb30645d42
using Plots

# ╔═╡ 8ae04a5a-0fa8-11eb-2821-8dc46ee73269
Pkg.develop("PolaronMobility")

# ╔═╡ 386fd0a4-0faa-11eb-14e3-a3cdf3f3cb5d
α=5.0

# ╔═╡ a86c04ae-0faa-11eb-3565-7bd3108906b2
feynmanvw(α) # Variational solution, athermal (original) action

# ╔═╡ 3e91bf88-0faa-11eb-33db-8b7fd24cf1eb
βred=5

# ╔═╡ b3e10abe-0faa-11eb-308e-91c4c0886b8b
feynmanvw(α, βred) # Variational solution, thermal (Osaka) action

# ╔═╡ 4306f2f4-0faa-11eb-1d32-dd8ade1d930a
ω=2.25E12

# ╔═╡ 4b26ec8c-0faa-11eb-2af1-e3f5e39fb358
meff=0.12

# ╔═╡ 34a7c84e-0faa-11eb-08b6-83f4ab86b689
nurange=0:0.1:22

# ╔═╡ ad10a50c-0fa8-11eb-30a9-1b1bb4fc26b7
X=PolaronMobility.ImX(nurange,feynmanvw(α, βred)...,βred,α,ω,meff*PolaronMobility.MassElectron)

# ╔═╡ f384207c-0fa8-11eb-29b6-9b35b33aa526
plot(X.nu, X.μ, label="Mobility")

# ╔═╡ a3cb113e-0fa9-11eb-1093-d7f12e1cf14b
plot(X.nu, X.ImX, label="ImX")

# ╔═╡ Cell order:
# ╠═66923a50-0fa8-11eb-1547-cf8e289818b0
# ╠═8336bf6e-0fa8-11eb-1f48-37cfe47a6d0f
# ╠═8ae04a5a-0fa8-11eb-2821-8dc46ee73269
# ╠═93125894-0fa8-11eb-0ff6-211feff97889
# ╠═386fd0a4-0faa-11eb-14e3-a3cdf3f3cb5d
# ╠═a86c04ae-0faa-11eb-3565-7bd3108906b2
# ╠═3e91bf88-0faa-11eb-33db-8b7fd24cf1eb
# ╠═b3e10abe-0faa-11eb-308e-91c4c0886b8b
# ╠═4306f2f4-0faa-11eb-1d32-dd8ade1d930a
# ╠═4b26ec8c-0faa-11eb-2af1-e3f5e39fb358
# ╠═34a7c84e-0faa-11eb-08b6-83f4ab86b689
# ╠═ad10a50c-0fa8-11eb-30a9-1b1bb4fc26b7
# ╠═f0529460-0fa8-11eb-1bec-1feb30645d42
# ╠═f384207c-0fa8-11eb-29b6-9b35b33aa526
# ╠═a3cb113e-0fa9-11eb-1093-d7f12e1cf14b


text:
using Documenter, SpaceExParser

DocMeta.setdocmeta!(SpaceExParser, :DocTestSetup,
                    :(using SpaceExParser); recursive=true)

makedocs(
    sitename = "SpaceExParser.jl",
    modules = [SpaceExParser],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        assets = ["assets/aligned.css"]),
    strict = true,
    pages = [
        "Home" => "index.md",
        "Examples" => Any["Introduction" => "examples/examples.md", "Bouncing ball" => "examples/bball.md"],
        "Library" => Any["Types" => "lib/types.md", "Methods" => "lib/methods.md"],
        "About" => "about.md"
    ]
)

deploydocs(
    repo = "github.com/JuliaReach/SpaceExParser.jl.git",
    push_preview = true
)