text:
# Logic puzzle benchmark from https://github.com/SWI-Prolog/bench
# Where does the zebra live?
# Puzzle solution written by Claude Sammut.

clauses = @julog [
	query(Houses) <<=
		houses(Houses) &
		member(house(red, english, _, _, _), Houses) &
		member(house(_, spanish, dog, _, _), Houses) &
		member(house(green, _, _, coffee, _), Houses) &
		member(house(_, ukrainian, _, tea, _), Houses) &
		right_of(house(green,_,_,_,_), house(ivory,_,_,_,_), Houses) &
		member(house(_, _, snails, _, winstons), Houses) &
		member(house(yellow, _, _, _, kools), Houses) &
		unifies(Houses, [_, _, house(_, _, _, milk, _), _,_]) &
		unifies(Houses, [house(_, norwegian, _, _, _)|_]) &
		next_to(house(_,_,_,_,chesterfields), house(_,_,fox,_,_), Houses) &
		next_to(house(_,_,_,_,kools), house(_,_,horse,_,_), Houses) &
		member(house(_, _, _, orange_juice, lucky_strikes), Houses) &
		member(house(_, japanese, _, _, parliaments), Houses) &
		next_to(house(_,norwegian,_,_,_), house(blue,_,_,_,_), Houses) &
		member(house(_, _, zebra, _, _), Houses) &
		member(house(_, _, _, water, _), Houses),
	houses([
		house(_, _, _, _, _),
		house(_, _, _, _, _),
		house(_, _, _, _, _),
		house(_, _, _, _, _),
		house(_, _, _, _, _)
	]) <<= true,

	right_of(A, B, [B, A | _]) <<= true,
	right_of(A, B, [_ | Y]) <<= right_of(A, B, Y),

	next_to(A, B, [A, B | _]) <<= true,
	next_to(A, B, [B, A | _]) <<= true,
	next_to(A, B, [_ | Y]) <<= next_to(A, B, Y),

	member(X, [X|_]) <<= true,
	member(X, [_|Y]) <<= member(X, Y)
]

println("Zebra puzzle (single run):")
goal = @julog query(Houses)
@time resolve(goal, clauses)


text:
module PolyChaos

using SpecialFunctions, SparseArrays, Distributions
import LinearAlgebra: I, dot, SymTridiagonal, eigen, issymmetric
import FFTW: ifft
import Combinatorics: with_replacement_combinations
import Base: show
import AdaptiveRejectionSampling: RejectionSampler, run_sampler!
import Statistics: mean, std, var
import GaussQuadrature: special_eigenproblem!

include("types.jl")
include("multiIndices.jl")
include("evaluate.jl")
include("quadrature_rules.jl")
include("recurrence_coefficients_monic.jl")
include("stieltjes.jl")
include("densities.jl")
include("auxfuns.jl")
include("scalarproduct.jl")
include("tensor.jl")
include("show.jl")
include("polynomial_chaos.jl")

end

ext: jl


text:
export
    AbstractSN,maxk!,to,@timeit

using TimerOutputs
const to = TimerOutput()
abstract type AbstractSN end

function maxk!(a, k::Int; initialized=false)
    # k:topk
    ix = collect(1:length(a))
    partialsortperm!(ix, a, 1:k, rev=true, initialized=initialized)
    # @views collect(zip(ix[1:k], a[ix[1:k]]))
    @views ix[1:k]
end

# @benchmark

module SN_pq
    using Main

    export
        SearchNeighbors_PQ,
        compute_table,
        compute_scores_

    struct SearchNeighbors_PQ <: AbstractSN
        M::Integer
        Ks::Integer
        D::Integer
        pq_codebook::Array{Float32, 3}
        pq_codes::Array{Int64, 2}
        metric::String
        Ds::Integer
        function SearchNeighbors_PQ(M, Ks, D, pq_codebook, pq_codes, metric)
            Ds = Int64(D/M)
            new(M, Ks, D, pq_codebook, pq_codes, metric, Ds)
        end
    end

    # typeof(x)
    # x = SearchNeighbors_PQ(M=8, Ks=16, D=32, metric = "dot_product")

    function exinit(sp::AbstractSN)
        for i in 1:sp.M
            sp.pq_codes[:,i] .+= sp.Ks*(i-1);
        end
    end

    function compute_table(codebook::Array{Float32,3}, query::Array{<:AbstractFloat,1})
        M, Ks, Ds = size(codebook)
        lookuptable = zeros(M, Ks);
        for i in 1:M
            C = @view codebook[i,:,:];
            q = @view query[ ((i-1)*Ds + 1) : ((i-1)*Ds + Ds)];
            lookuptable[i,:] = C*q;
        end
        return lookuptable
    end

    function compute_scores_(lookuptable::Array{Float64,2}, pq_codes::AbstractArray{Int,2})
        n, M = size(pq_codes);
        scores = zeros(n);
        for i = 1:n
            s_ = 0;
            for j = 1:M
                @inbounds s_ = s_ + lookuptable[j, pq_codes[i,j]];
            end
            scores[i] = s_;
        end
        return scores
    end

    function compute_scores(sp::AbstractSN, query::Array{Float64,1})
        codebook = sp.pq_codebook;
        pq_codes = sp.pq_codes;

        @timeit to "compute_table" lookuptable = compute_table(codebook,query);
        @timeit to "compute_scores_" scores = compute_scores_(lookuptable,pq_codes);

        return scores
    end

    function search_neighbors(sp::AbstractSN, query::Array{Float64,1},topk::Int)
        scores = compute_scores(sp, query)
        @timeit to "maxk!" indices = maxk!(scores, topk);
        return indices
    end

    function get_neighbors(sp::AbstractSN, queries::Array{Float64,2},topk::Int)
        n = size(queries)[1];
        neighbors_matrix = zeros(Int,n,topk)
        # Threads.@threads
        for i in 1:n
            q = queries[i,:];
            neighbors_matrix[i,:] = search_neighbors(sp, q, topk);
        end
        return neighbors_matrix
    end

end

module SN_ivf
    using Main
    export
        coarse_search,
        SearchNeighbors_IVF

    struct SearchNeighbors_IVF <: AbstractSN
        k_v::Integer
        D::Integer
        vq_codebook::Array{<:AbstractFloat, 2}
        vq_codes::Array{<:Integer, 1}
        metric::String

        function SearchNeighbors_IVF(vq_codebook, vq_codes, metric)
            k_v,D = size(vq_codebook)
            new(k_v, D, vq_codebook, vq_codes, metric)
        end
    end

    function coarse_search(si::AbstractSN, query::Vector, num_centroid_tosearch::Integer)
        score = si.vq_codebook*query # (k_v,)
        index = maxk!(score, num_centroid_tosearch)

        score_v = score[index]
        return index, score_v
    end

end

module SN_pqivf

    export SearchNeighbors_PQIVF

    using Main
    using ..SN_pq
    using ..SN_ivf

    struct SearchNeighbors_PQIVF <: AbstractSN
        M::Int8
        Ks::Int8
        D::Int16
        Ds::Int16
        pq_codebook::Array{<:AbstractFloat, 3}
        pq_codes::Array{<:Integer, 2}

        k_v::Int64
        vq_codebook::Array{<:AbstractFloat, 2}
        vq_codes::Array{<:Integer, 1}
        metric::String
        id_index_table::Dict{Int64,Vector{Int64}}

        function SearchNeighbors_PQIVF(M, Ks, pq_codebook, pq_codes, vq_codebook, vq_codes, metric)
            si = SearchNeighbors_IVF(vq_codebook, vq_codes, metric)
            k_v, D, vq_codebook, vq_codes, metric = si.k_v, si.D, si.vq_codebook, si.vq_codes, si.metric

            sp = SearchNeighbors_PQ(M, Ks, D, pq_codebook, pq_codes, metric)
            M,Ks,Ds,pq_codebook,pq_codes = sp.M, sp.Ks, sp.Ds, sp.pq_codebook, sp.pq_codes

            dict = Dict()
            for i in 1:k_v
                dict[i] = []
            end

            n = length(vq_codes)
            for i in 1:n
                codei = vq_codes[i]
                append!(dict[codei],i)
            end

            new(M, Ks, D, Ds, pq_codebook, pq_codes, k_v, vq_codebook, vq_codes, metric, dict)
        end
    end

    function id_indexed(spi::AbstractSN, index::SubArray{Int64, 1}, scores_v::Vector{<:AbstractFloat})
        ids = [] # Array{Int64,1}
        scores_vq = []
        for (i,s) in zip(index,scores_v)
            id_i = spi.id_index_table[i]
            ids = [ids;id_i]
            len = length(id_i)
            scores_vq = [scores_vq;repeat([s],len)]
        end
        ids = convert(Array{Int64,1}, ids)
        scores_vq = convert(Array{Float64,1}, scores_vq)

        return ids, scores_vq
    end

    function compute_scores_ivf(lookuptable::Array{Float64,2}, pq_codes::AbstractArray{Int,2}, ids::Array{Int64,1})
        n0, M = size(pq_codes);
        n = length(ids)
        scores = zeros(n);
        for i = 1:n
            s_ = 0;
            for j = 1:M
                @inbounds s_ = s_ + lookuptable[j, pq_codes[ids[i],j]];
            end
            scores[i] = s_;
        end
        return scores
    end

    function search_neighbors(spi::AbstractSN, q::Vector{<:AbstractFloat}, num_centroid_tosearch::Int, topk::Int)
        @timeit to "coarse_search" index,scores_v = coarse_search(spi, q, num_centroid_tosearch)
        @timeit to "id_indexed" ids, scores_vq = id_indexed(spi, index, scores_v)

        # ix_ = sortperm(ids)
        # ids = ids[ix_]
        # scores_vq = scores_vq[ix_]

        # ids = 1:length(scores_vq)
        # ids = convert(Array{Int64,1}, ids)

        @timeit to "compute_table" lookuptable = compute_table(spi.pq_codebook, q) # (M,Ks)
        @timeit to "ex1" pq_codes_ids = @view spi.pq_codes[ids,:]

        # @timeit to "ex1" pq_codes_ids = spi.pq_codes[ids,:]
        # @timeit to "compute_scores_" scores_pq = compute_scores_(lookuptable, pq_codes_ids)

        @timeit to "compute_scores_ivf" scores_pq = compute_scores_ivf(lookuptable, spi.pq_codes, ids)
        @timeit to "ex2" scores = scores_vq + scores_pq;

        @timeit to "ex3" len_s = length(scores);
        @timeit to "ex4" if topk > len_s
            i_ = maxk!(scores, len_s)
            ids_ = ids[i_]
            ids_ = [ids_;repeat([0],topk - len_s)]
        else
            i_ = maxk!(scores, topk)
            ids_ = ids[i_]
        end
        return ids_
    end

    # @timeit to
    function get_neighbors(spi::AbstractSN, queries::Array{<:AbstractFloat,2}, num_centroid_tosearch::Int,  topk::Int)
        n = size(queries)[1];
        neighbors_matrix = zeros(Int,n,topk)
        # Threads.@threads
        for i in 1:n
            q = queries[i,:];
            # @timeit to "sn"
            neighbors_matrix[i,:] = search_neighbors(spi, q, num_centroid_tosearch, topk);
        end
        return neighbors_matrix
    end

end


text:
# ### stacked bar plot

#md # ![FILE_NAME.png](images/FILE_NAME.png)

using StatsPlots

measles = [38556, 24472, 14556, 18060, 19549, 8122, 28541, 7880, 3283, 4135, 7953, 1884]
mumps = [20178, 23536, 34561, 37395, 36072, 32237, 18597, 9408, 6005, 6268, 8963, 13882]
chickenPox = [37140, 32169, 37533, 39103, 33244, 23269, 16737, 5411, 3435, 6052, 12825, 23332]

## In PyPlot backend, if we use chars like 'A':'L', ticks are displayed with "PyWrap".
ticklabel = string.(collect('A':'L'))
groupedbar([measles mumps chickenPox],
        bar_position = :stack,
        bar_width=0.7,
        xticks=(1:12, ticklabel),
        label=["measles" "mumps" "chickenPox"])


text:
# This deletes arbitrary tokens from files and checks that we can still parse them
# and that iteration functions are still correctly ordered.

@testset "invalid jl file parsing" begin
    function trav(x, f = x->nothing)
        f(x)
        for a in x
            trav(a, f)
        end
    end

    function trav1(x, f = x->nothing)
        f(x)
        if x.args !== nothing
            for a in x
                trav(a, f)
            end
        end
    end

    function check_err_parse(s, n = length(s)÷100)
        check_str(s) # parsing works?
        check_itr_order(s) # iteration produces same text?

        ts = collect(tokenize(s))[1:end-1]
        for _ in 1:n
            length(ts) == 1 && return
            deleteat!(ts, rand(1:length(ts)))
            check_str(untokenize(ts))
        end
    end

    function check_str(s)
        x = try
            CSTParser.parse(s, true)
        catch e
            @info "Couldn't parse:"
            @info s
            rethrow(e)
        end
        try
            trav(x)
        catch e
            @info "Couldn't traverse:"
            @info s
            rethrow(e)
        end
        try
            to_codeobject(x)
        catch e
            @info "Couldn't convert:"
            @info s
            rethrow(e)
        end
        if sizeof(s) != x.fullspan
            @info "sizeof(s) != x.fullspan for :"
            println(">>>>>>>>>>>>>>>>>>")
            println(s)
            println("<<<<<<<<<<<<<<<<<<")
        end
        check_itr_order(s, x)
    end

    function get_segs(x)
        offset = 0
        segs = []
        for i = 1:length(x)
            a = x[i]
            push!(segs, offset .+ (1:a.fullspan))
            offset += a.fullspan
        end
        segs
    end

    function check_itr_order(s, x = CSTParser.parse(s, true))
        length(x) == 0 && return
        segs = get_segs(x)
        s0 = join(String(codeunits(s)[seg]) for seg in segs)
        if s0 == s
            for i = 1:length(x)
                if length(x[i]) > 0
                    seg = segs[i]
                    s2 = String(codeunits(s)[seg])
                    check_itr_order(s2, x[i])
                end
            end
        else
            @info "check_itr_order failed: "
            println(">>>>>>>>>>>>>>>>>>")
            println(s)
            println("<<<<<<<<<<<<<<<<<<")
            error()
        end
    end

    comp(x, y) = x == y
    function comp(x::CSTParser.EXPR, y::CSTParser.EXPR)
        comp(x.head, y.head) &&
            x.span == y.span &&
            x.fullspan == y.fullspan &&
            x.val == y.val &&
            length(x) == length(y) &&
            all(comp(x[i], y[i]) for i = 1:length(x))
    end

    function check_reparse(s0, n = length(s0)÷100)
        for _ in 1:n
            x0 = CSTParser.parse(s0, true)
            CSTParser.has_error(x0) && return
            ts = collect(tokenize(s0))[1:end-1]
            length(ts) < 2 && return
            deleteat!(ts, rand(1:length(ts)))
            s1 = untokenize(ts)
            length(ts) < 3 || isempty(s1) && return
            x1 = CSTParser.parse(s1, true)
            x2 = try
                CSTParser.minimal_reparse(s0, s1, x0, x1)
            catch err
                @info "minimal reparse failed with"
                @info "s0:"
                @info s0
                @info "s1:"
                @info s1
                rethrow(err)
            end
            @test comp(x1, x2) ? true : (@info(string("Comparison failed between s0:\n", s0, "\n\n and s1: \n", s1)); false)
            @test try
                join(String(codeunits(s0)[seg]) for seg in get_segs(x0)) == s0
            catch e
                @info "Couldn't reconstruct original text from EXPR segments for:"
                @info s0
                rethrow(e)
            end
            @test try
                join(String(codeunits(s1)[seg]) for seg in get_segs(x1)) == s1
            catch e
                @info "Couldn't reconstruct original text from EXPR segments for:"
                @info s1
                rethrow(e)
            end
            @test join(String(codeunits(s1)[seg]) for seg in get_segs(x2)) == s1
            s0 = s1
        end
    end

    function check_dir(dir, check)
        for (root, _, files) in walkdir(dir)
            for f in files
                f = joinpath(root, f)
                (!isfile(f) || !endswith(f, ".jl")) && continue
                @info "checking $(nameof(check)) against $f"
                s = String(read(f))
                if isvalid(s) && length(s) >0
                    check(s)
                end
            end
        end
        true
    end

    @test check_dir("..", check_err_parse)
    @test check_dir("..", check_reparse)
end


text:
# addprocs, same fun everywhere!
using Distributed
while nworkers()<4
    addprocs(1)
end
using Pkg
@everywhere using Pkg
@everywhere Pkg.activate(".")
using DEMC

using Random
using LinearAlgebra
using Distributions
using SharedArrays
using ParallelDataTransfer

@everywhere begin
    using DEMC
    using Random
    using LinearAlgebra
    using Distributions
    using SharedArrays
    Random.seed!(31953150+myid())
    using ParallelDataTransfer
end

# set up target distribution: Multivariate Normal
ndim = 10 # Number of dimensions
μ = rand(ndim) # mean of each dimension
A = rand(ndim, ndim)
Σ = A'*A .+ diagm(0 => 2*ones(ndim)) # variance covariance matrix
Σ = Σ./maximum(Σ)/100
distr = MvNormal(μ, Σ)
passobj(myid(), workers(),[:μ, :Σ, :distr])
# log objective function
@everywhere log_obj(mean) = logpdf(distr, mean)

# set up of DEMCz chain
opts = DEMC.demcopt(ndim)
opts.blockindex = [1:ndim] # parameter blocks: here choose all parameters to be updated simultaenously
opts.Nblocks = length(opts.blockindex)
opts.eps_scale = 1e-5*ones(ndim) # scale of random error around DE update
opts.γ = 2. # scale of DE update, 2.38 is the "optimal" number for a normal distribution
opts.N = nworkers() # number of chains
opts.K = 10 # every K steps add current N draws to Z
# Number of iterations in Chain
opts.Ngeneration = 100_000
opts.autostop = :Rhat
opts.autostop_Rhat = 1.075

Z = randn((10*ndim, ndim)) # initial distribution (completely off to make a difficult test case)
mc, Z = DEMC.demcz_sample_par(log_obj, Z,opts; prevrun=nothing, sync_every = 5000, workerpool=CachingPool(workers()))

# drop first half of chain
Ntot = size(mc.chain,3)
Nburn = 30000
if Nburn >= round(Ntot/2)+1
    Nburn = Int(round(Ntot/2))
end
keep = Ntot-Nburn+1:Ntot
Ngen_burned = length(keep)
chain_burned = mc.chain[:,:,keep]
logobj_burned = mc.log_obj[:, keep]
chainflat = DEMC.flatten_chain(chain_burned, opts.N, Ngen_burned, ndim)'
# covariance of estimates
b, Σb = DEMC.mean_cov_chain(chain_burned, opts.N, Ngen_burned, ndim)
println("\n estimates: ", b, "\n dist to true: ", b - μ)

figure_path = pwd()*"/img/normpdf_parallel/"
accept_ratio, Rhat = DEMC.convergence_check(chain_burned, logobj_burned, figure_path; verbose = false)

@test all(Rhat.<1.1)
@test all(accept_ratio.> 0.1)
@test all(accept_ratio.< 0.45)

rmprocs(workers())

# sample = [rand(distr) for s=1:100_000]
# llsample = [logpdf(distr, sample[s]) for s=1:100_000]
# DEMC.histogram(logobj_burned[:], normalize=true)
# DEMC.histogram!(llsample, normalize=true)
# DEMC.savefig(figure_path*"lpdf_histogram.png")


text:
using DynamicalSystems
using Base.Test

println("\nTesting continuous system evolution...")

@testset "Lorenz System" begin

  lo11 = Systems.lorenz()
  lo22 = ContinuousDS(lo11.state, lo11.eom)
  lo33 = Systems.lorenz(big.([0.0, 10.0, 0.0]))
  @testset "Construction" begin
    @test typeof(lo11) <: ContinuousDS
    @test typeof(lo33) <: ContinuousDS
  end
  @testset "Evolve & kwargs" begin
    # test evolve(state):
    st1 = evolve(lo11.state, lo11, 1.0)
    st3 = evolve(lo33.state, lo33, 1.0)
    @test st1 ≈ st3
    # Test evolve(state,keywords):
    st1 = evolve(lo11.state, lo11, 1.0; diff_eq_kwargs=Dict(:abstol=>1e-9, :reltol=>1e-9))
    st3 = evolve(lo33.state, lo33, 1.0, diff_eq_kwargs=Dict(:abstol=>1e-9, :reltol=>1e-9))
    @test st1 ≈ st3
  end

  @testset "Timeseries" begin
    # timeseries pure:
    ts1 = timeseries(lo11, 2.0; mutate = false)
    ts3 = timeseries(lo33, 2.0; mutate = false)
    @test eltype(ts3) == BigFloat
    @test size(ts1) == size(ts3)
    @test ts1[end, :] ≈ ts3[end,:]
    # timeseries with diff_eq_kwargs and dt:
    ts1 = timeseries(lo11, 2.0; dt=0.1,
    diff_eq_kwargs=Dict(:abstol=>1e-9, :reltol=>1e-9), mutate = false)
    ts3 = timeseries(lo33, 2.0; dt=0.1,
    diff_eq_kwargs=Dict(:abstol=>1e-9, :reltol=>1e-9), mutate = false)
    @test size(ts1) == size(ts3)
    @test ts1[end, :] ≈ ts3[end,:]
  end

  @testset "Jacobians" begin
    j1 = jacobian(lo11); j2 = jacobian(lo22); j3 = jacobian(lo33)
    @test eltype(j3) == BigFloat
    @test j1 ≈ j2
    @test j1 ≈ j3
    s1 = evolve(lo11, 1.0)
    s2 = evolve(lo22, 1.0)
    s3 = evolve(lo33, 1.0)
    j1 = lo11.jacob(s1); j2 = lo22.jacob(s2); j3 = lo33.jacob(s3)
    @test eltype(j3) == BigFloat
    @test j1 ≈ j2
    @test j1 ≈ j3
  end

end


text:
using Gadfly

set_default_plot_size(6inch, 3inch)

# make sure Theme.{discrete,continuous}_highlight_color works for Geom.point.

discrete_with_highlights = plot(x=rand(1:3,100),y=rand(100),color=rand(["a","b"],100),
            Geom.point, Stat.x_jitter);
discrete_without_highlights = plot(x=rand(1:3,100),y=rand(100),color=rand(["a","b"],100),
            Geom.point, Stat.x_jitter, Theme(discrete_highlight_color=x->x));
continuous_with_highlights = plot(x=rand(1:3,100),y=rand(100),color=rand(1:2,100),
            Geom.point, Stat.x_jitter);
continuous_without_highlights = plot(x=rand(1:3,100),y=rand(100),color=rand(1:2,100),
            Geom.point, Stat.x_jitter, Theme(continuous_highlight_color=x->x));
gridstack([discrete_with_highlights   discrete_without_highlights
           continuous_with_highlights continuous_without_highlights])


text:
using LinearAlgebra
using Distributions
using ForwardDiff
using Pathfinder
using Test

@testset "multi path pathfinder" begin
    @testset "MvNormal" begin
        n = 10
        nruns = 20
        ndraws = 1000_000
        ndraws_per_run = ndraws ÷ nruns
        Σ = rand_pd_mat(Float64, n)
        μ = randn(n)
        d = MvNormal(μ, Σ)
        logp(x) = logpdf(d, x)
        ∇logp(x) = ForwardDiff.gradient(logp, x)
        x₀s = [rand(Uniform(-2, 2), n) for _ in 1:nruns]
        q, ϕ, component_ids = multipathfinder(
            logp, ∇logp, x₀s, ndraws; ndraws_elbo=100, ndraws_per_run=ndraws_per_run
        )
        @test q isa MixtureModel
        @test ncomponents(q) == nruns
        @test Distributions.component_type(q) <: MvNormal
        @test ϕ isa AbstractMatrix
        @test size(ϕ) == (n, ndraws)
        @test component_ids isa Vector{Int}
        @test length(component_ids) == ndraws
        @test extrema(component_ids) == (1, nruns)
        μ_hat = mean(ϕ; dims=2)
        Σ_hat = cov(ϕ .- μ_hat; dims=2, corrected=false)
        # adapted from the MvNormal tests
        # allow for 15x disagreement in atol, since this method is approximate
        multiplier = 15
        for i in 1:n
            atol = sqrt(Σ[i, i] / ndraws) * 8 * multiplier
            @test μ_hat[i] ≈ μ[i] atol = atol
        end
        for i in 1:n, j in 1:n
            atol = sqrt(Σ[i, i] * Σ[j, j] / ndraws) * 10 * multiplier
            @test isapprox(Σ_hat[i, j], Σ[i, j], atol=atol)
        end
    end
end


text:
#---------------------------------------------------------
# # [Phase Vocoder](@id 02-vocoder)
#---------------------------------------------------------

#=
This page illustrates the phase vocoder feature of the Julia package
[`Sound`](https://github.com/JeffFessler/Sound.jl).

This page was generated from a single Julia file:
[02-vocoder.jl](@__REPO_ROOT_URL__/02-vocoder.jl).
=#

#md # In any such Julia documentation,
#md # you can access the source code
#md # using the "Edit on GitHub" link in the top right.

#md # The corresponding notebook can be viewed in
#md # [nbviewer](http://nbviewer.jupyter.org/) here:
#md # [`02-vocoder.ipynb`](@__NBVIEWER_ROOT_URL__/02-vocoder.ipynb),
#md # and opened in [binder](https://mybinder.org/) here:
#md # [`02-vocoder.ipynb`](@__BINDER_ROOT_URL__/02-vocoder.ipynb).

# ### Setup

# Packages needed here.

using Sound: phase_vocoder, soundsc, hann
using DSP: spectrogram
using MIRTjim: jim, prompt
using InteractiveUtils: versioninfo

# The following line is helpful when running this file as a script;
# this way it will prompt user to hit a key after each figure is displayed.

isinteractive() ? jim(:prompt, true) : prompt(:draw);

# ### Overview

#=
Here we illustrate applying the phase vocoder
to stretch time (by the default factor of 2).
=#

S = 8192*2^2
x1 = cos.(2π*400*(1:2S)/S) .* hann(2S)
x2 = cos.(2π*300*(1:S)/S).^3 .* hann(S) # nonlinearity for harmonics
x = [x1; x2] # 2-note song
y = phase_vocoder(x, S)
length(x), length(y)

# spectrograms
sx = spectrogram(x / maximum(abs, x); fs=S)
sy = spectrogram(y / maximum(abs, y); fs=S)
dB = x -> max(log10.(x), -10)
fun = (s) -> jim(s.time, s.freq, dB.(s.power)', aspect_ratio=:auto,
    ylims = (0,1800), color=:viridis, clim = (-10,0),
    yticks = 0:300:1800, yflip=false,
    xlabel="time", ylabel="freq",
    prompt = false,
)
jim(fun(sx), fun(sy))

# listen to before/after:
soundsc([x; y], S)

# ### Reproducibility

# This page was generated with the following version of Julia:

io = IOBuffer(); versioninfo(io); split(String(take!(io)), '\n')

# And with the following package versions

import Pkg; Pkg.status()


text:
# Here's simple addition
1 + 2 # 3

# Operators lik `+` are just functions in Julia
+ # + (generic function with 171 methods)

# So, you can do this (where 1 and 2 are arguments to the function `+`)
+(1, 2) # 3

# Or even this (Does this remind you of Lisp?)
+(1, 2, 3, 4) # 10

# Let's define a function
function add(a, b)
    # Julia function returns the value of the last expression, 'return' is only optional
    +(a, b)
end

add # add (generic function with 1 method)
add(1, 2) # 3

# With a nested function
function add2(a)
    function nested(b)
        +(a, b)
    end
end

add2 # add2 (generic function with 1 method)
add2(1)(2) # 3

# Anonymous function
add3 = ((a, b) -> +(a, b))

add3 # (anonymous function)
add3(1, 2) # 3

# Anonymous function with multi-line body
add4 = begin
    (a, b) ->
    +(a, b)
end

add4 # (anonymous function)
add4(1, 2) # 3

# With "varargs"
function add5(a, b, c...)
    +(a, b, sum(c))
end

add5 # (anonymous function)
add5(1, 2, 3, 4) # 10

# With optional arguments
function add6(a, b, c=0)
    +(a, b, c)
end

add6 # (anonymous function)
add6(1, 2, 3) # 6
add6(1, 2, c=3) # ERROR: ArgumentError: function add6 does not accept keyword arguments

# With keyword arguments
function add7(a, b; c=0)
    +(a, b, c)
end

add7(1, 2, c=3) # 6


text:
using ConservationLawsParticles
using RecursiveArrayTools, DifferentialEquations, Plots

# model
V(t, x) = 0.0 #2 + 0.5sin(2x) + 0.3sin(2sqrt(2)*x) + 0.2cos(2sqrt(7)*x)
W(t, x) = 0.0 #x^2/2 - abs(x)
W′(t, x) = 0.0 #x - sign(x)
mob(ρ) = 0.0 # max((ρ-1)^2 * (2ρ+1), 0)
smodel = DiffusiveSampledModel((V,), ((W′,),), (mob,), (Diffusion(1),))
# imodel = DiffusiveIntegratedModel((V,), ((W,),), (mob,))

# initial condition
x0 = ArrayPartition(gaussian_particles(3, 201))

# time span
tspan = (0., 15.)

# ODE system for the particles
dprob = ODEProblem(velocities_diff!, x0, tspan, smodel)

abstol = reltol = 2e-8

@time dsol = solve(dprob, BS5(); abstol=abstol, reltol=reltol);

# plot the particle trajectories
plot(legend=false, title="Trajectories", xlabel="time", ylabel="position")
plot!(dsol, vars=1:5:201; color=:green)
savefig("diffusion1.png")

# plot an animation of the density
anim = @animate for t in range(0.0, 7.0; step=1/24)
    plot(title="Density", xlabel="position", ylabel="density",
        xlims=(-10,10), ylims=(0,1/sqrt(pi)))
    densityplot!(dsol(t); color=:green, label="particles")
    plot!(x -> exp(-x^2/(1+4t))/sqrt(pi*(1+4t)); color=:orange, label="analytic")
end
gif(anim, "diffusion1.gif")


text:
module XCB

const xcb = :libxcb

include("c_types.jl")
include("c_interface.jl")

end # module


text:
module Ext
using ..Mosek

immutable VarByIndex index :: Int32 end
immutable VarByName name :: String end
immutable BarvarByIndex index :: Int64 end
immutable BarvarByName name :: String end
immutable ConByIndex index :: Int32 end
immutable ConByName  name  :: String end
immutable ConeByIndex index :: Int32 end
immutable ConeByName  name  :: String end

immutable Variable
    t::Mosek.Task
    index :: Int32
end

immutable SemidefiniteVariable
    t::Mosek.Task
    index :: Int32
end

immutable Constraint
    t::Mosek.Task
    index :: Int32
end

immutable ConeConstraint
    t::Mosek.Task
    index :: Int32
end

immutable Obj
end

immutable Objective
    t::Mosek.Task
end

immutable Symmat
    index :: Int
end

Base.getindex(t::Mosek.Task, index :: VarByIndex) = Variable(t,index.index)
function Base.getindex(t::Mosek.Task, name :: VarByName)
    ok,index = getvarnameindex(t,name.name)
    if ok > 0
        getindex(t,VarByIndex(index))
    else
        throw(KeyError("'$name' not found"))
    end
end

Base.getindex(t::Mosek.Task, index :: BarvarByIndex) = SemidefiniteVariable(t,index.index)
function Base.getindex(t::Mosek.Task, name :: BarvarByName)
    ok,index = getbarvarnameindex(t,name.name)
    if ok > 0
        getindex(t,BarvarByIndex(index))
    else
        throw(KeyError("'$name' not found"))
    end
end

Base.getindex(t::Mosek.Task, index :: ConByIndex) = Constraint(t,index.index)
function Base.getindex(t::Mosek.Task, name :: ConByName)
    ok,index = getconnameindex(t,name.name)
    if ok > 0
        getindex(t,ConByIndex(index))
    else
        throw(KeyError("'$name' not found"))
    end
end

Base.getindex(t::Mosek.Task, index :: ConeByIndex) = ConeConstraint(t,index.index)
function Base.getindex(t::Mosek.Task, name :: ConeByName)
    ok,index = getconenameindex(t,name.name)
    if ok > 0
        getindex(t,ConeByIndex(index))
    else
        throw(KeyError("'$name' not found"))
    end
end

Base.getindex(t::Mosek.Task, index :: Obj) = Objective(t)
function Base.getindex(t::Mosek.Task, index :: Symmat)
    dim,nz,tp = getsymmatinfo(t,index.index)
    subi,subj,valij = getsparsesymmat(t,index.index)
    for i in 1:nz
        if subi[i] != subj[i]
            push!(subi,subj[i])
            push!(subj,subi[i])
            push!(valij,valij[i])
        end
    end
    sparse(subi,subj,valij,dim,dim)
end

function Base.show(f::IO, var :: Variable)
    bk,bl,bu = getvarbound(var.t,var.index)
    name = getvarname(var.t,var.index)

    if bk == MSK_BK_FX
        print(f,"Variable('$name' = $bl)")
    else
        blstr =
            if     bk == MSK_BK_LO || bk == MSK_BK_RA "[$bl"
            else   "]-Inf"
            end
        bustr =
            if     bk == MSK_BK_UP || bk == MSK_BK_RA "$bu]"
            else   "+Inf["
            end

        print(f,"Variable('$name' ∈ $blstr;$bustr)")
    end
end

function Base.show(f::IO, var :: SemidefiniteVariable)
    if var.index > 0 && var.index < getnumbarvar(var.t)
        dim,ns,tp = getsymmatinfo(var.t,var.index)
        name = mkbarvarname(var.t,var.index)

        print(f,"SemidefiniteVariable('$name' ∈ S($dim))")
    else
        error("Invalid reference")
    end
end

function Base.show(f::IO, cone :: ConeConstraint)
    name = getconename(cone.t,cone.index)
    ct,cp,nummem,submem = getcone(cone.t,cone.index)

    dom =
        if ct == MSK_CT_QUAD "Q"
        elseif ct == MSK_CT_RQUAD "Q_r"
        else "?"
        end

    varnames = join(map(
        j -> let varname = getvarname(cone.t,j)
               if length(varname) == 0
                 "#x$j"
               else
                 Mosek.escapename(varname)
               end
             end,
        submem),",")

    print(f,"ConeByIndex('$name': ($varnames) ∈ $dom)")
end

mkvarname(t,j) =
    let n = getbarvarname(t,j)
        if length(n) == 0
            "#x$j"
        else
            Mosek.escapename(n)
        end
    end

mkbarvarname(t,j) =
    let n = getbarvarname(t,j)
        if length(n) == 0
            "#X̂$j"
        else
            Mosek.escapename(n)
        end
    end

function Base.show(f::IO, con :: Objective)
    t = con.t
    name = getobjname(t)

    print(f,"Objective('$name': ")

    for j in 1:getnumvar(t)
        cj = getcj(t,j)
        if cj < 0 || cj > 0
            varname = Mosek.escapename(getvarname(t,j))
            if length(varname) == 0
                varname = "#x$j"
            end
            print(f,"$(Mosek.fmtcof(cj)) $varname ")
        end
    end

    if getnumbarvar(t) > 0
        _,barcidx = getbarcsparsity(t)

        for idx in barcidx
            (barcj,num,sub,w) = getbarcidx(t,idx)

            barvarname = mkbarvarname(t,barcj)

            if num == 1
                print(f,"$(Mosek.fmtcof(w[1])) #MX$(sub[1]) ⋅ $barvarname) ")
            elseif num > 1
                print(f,"(")
                printf(f,"$(Mosek.fmtcof(w[1])) #MX$(sub[1])")
                print(f,") ⋅ $barvarname ")
            end
        end
    end

    numqobjnz = getnumqobjnz(t)
    if numqobjnz > 0
        (nnz,qcsubi,qcsubj,qcval) = getqobj(t)
        for k in 1:nnz
            if qcsubi[k] == qcsubj[k]
                ni = mkvarname(t,qcsubi[k])
                print(f,"$(Mosek.fmtcof(qcval[k])) $ni ^ 2")
            else
                ni = mkvarname(t,qcsubi[k])
                nj = mkvarname(t,qcsubj[k])
                print(f,"$(Mosek.fmtcof(qcval[k])) $ni ⋅ $nj")
            end
        end
    end

    print(f,")")
end

function Base.show(f::IO, con :: Constraint)
    t = con.t
    i = con.index
    bk,bl,bu = getconbound(t,i)
    name = getconname(t,i)

    print(f,"Constraint('$name': ")
    if bk == MSK_BK_FX
        print(f,"$bl = ")
    else
        blstr =
            if     bk == MSK_BK_LO || bk == MSK_BK_RA "[$bl"
            else   "]-Inf"
            end
        bustr =
            if     bk == MSK_BK_UP || bk == MSK_BK_RA "$bu]"
            else   "+Inf["
            end

        print(f,"$blstr;$bustr ∋ ")
    end

    (nzi,subi,vali) = getarow(t,Int32(i))

    for k in 1:nzi
        varname = Mosek.escapename(getvarname(t,subi[k]))
        if length(varname) == 0
            varname = "#x$(subi[k])"
        end
        print(f,"$(Mosek.fmtcof(vali[k])) $varname ")
    end

    numqconnz = getnumqconknz(t,i)
    if numqconnz > 0
        (nnz,qcsubi,qcsubj,qcval) = getqconk(t,i)
        for k in 1:nnz
            if qcsubi[k] == qcsubj[k]
                ni = mkvarname(t,qcsubi[k])
                print(f,"$(Mosek.fmtcof(qcval[k])) $ni ^ 2")
            else
                ni = mkvarname(t,qcsubi[k])
                nj = mkvarname(t,qcsubj[k])
                print(f,"$(Mosek.fmtcof(qcval[k])) $ni ⋅ $nj")
            end
        end
    end

    _,baraidx = getbarasparsity(t)
    if length(baraidx) > 0
        # count
        numbarnzi = count(idx -> getbaraidxij(t,idx)[1] == i, baraidx)
        idxs = Vector{Int}(numbarnzi)
        k = 0
        for idx in baraidx
            (barai,baraj) = getbaraidxij(t,idx)
            if barai == i
                k += 1
                idxs[k] = idx
            end
        end

        for idx in idxs
            (barai,baraj,num,sub,w) = getbaraidx(t,idx)
            barvarname = mkbarvarname(t,baraj)

            if num == 1
                print(f,"$(Mosek.fmtcof(w[1])) #MX$(sub[1]) ⋅ $barvarname ")
            elseif num > 1
                print(f,"(")
                printf(f,"$(Mosek.fmtcof(w[1])) #MX$(sub[1])")
                print(f,") ⋅ $barvarname ")
            end
        end
    end

    print(f,")")
end

function findvars(t::Mosek.Task, name::String)
    res = VarByIndex[]
    for j in 1:getnumvar(t)
        if getvarname(t,j) == name
            push!(res,VarByIndex(j))
        end
    end
    res
end

function findvars(t::Mosek.Task, regex::Regex)
    res = VarByIndex[]
    for j in 1:getnumvar(t)
        if ismatch(regex,getvarname(t,j))
            push!(res,VarByIndex(j))
        end
    end
    res
end

function findcons(t::Mosek.Task, name::String)
    res = ConByIndex[]
    for j in 1:getnumcon(t)
        if getconname(t,j) == name
            push!(res,ConByIndex(j))
        end
    end
    res
end

function findcons(t::Mosek.Task, regex::Regex)
    res = ConByIndex[]
    for j in 1:getnumcon(t)
        if ismatch(regex,getconname(t,j))
            push!(res,ConByIndex(j))
        end
    end
    res
end

function findcones(t::Mosek.Task, name::String)
    res = ConeByIndex[]
    for j in 1:getnumcone(t)
        if getconename(t,j) == name
            push!(res,ConeByIndex(j))
        end
    end
    res
end

function findcones(t::Mosek.Task, regex::Regex)
    res = ConeByIndex[]
    for j in 1:getnumcone(t)
        if ismatch(regex,getconename(t,j))
            push!(res,ConeByIndex(j))
        end
    end
    res
end

Var{I <: Integer}(index::I) = VarByIndex(convert(Int32,index))
Var(name::String) = VarByName(name)

Barvar{I <: Integer}(index::I) = BarvarByIndex(convert(Int32,index))
Barvar(name::String) = BarvarByName(name)

Con{I <: Integer}(index::I) = ConByIndex(convert(Int32,index))
Con(name::String) = ConByName(name)

Cone{I <: Integer}(index::I) = ConeByIndex(convert(Int32,index))
Cone(name::String) = ConeByName(name)

export Var,Barvar,Con,Cone,Obj,Symmat,Barvar,findvars,findcons,findcones

end


text:
module MakieCore

using Observables
using GeometryBasics
using Colors
using StaticArrays
using Random
using LinearAlgebra

include("utils.jl")
include("quaternions.jl")
include("types.jl")
include("conversion.jl")
include("scatter.jl")
include("backend.jl")

end


text:
module Bump10MSun

using Cosmology
using HDF5
using StatsBase
using StatsFuns
using Turing
using Unitful
using UnitfulAstro
using Trapz

include("model.jl")
export mlow, mhigh, zref
export make_log_dN, broken_pl_model
export make_dNdm1, make_dNdm2, make_dNdq

include("weights.jl")
export li_prior_wt, md_sfr, md_sfr_zwt, read_selection
export resample_selection

include("utils.jl")
export with_seed, append_generated_quantities

end # module


text:
__precompile__(true)

module VarianceCovarianceEstimators

using StatsBase
import StatsBase: vcov

## Abstract and Structs
abstract type AbstractVCE end
abstract type AbstractHCVCE <: AbstractVCE end
abstract type AbstractCRVCE <: AbstractVCE end
struct OLS <: AbstractVCE end
struct HC1 <: AbstractHCVCE end
struct HC2 <: AbstractHCVCE end
struct HC3 <: AbstractHCVCE end
struct CRVE1 <: AbstractCRVCE end
struct CRVE2 <: AbstractCRVCE end
struct CRVE3 <: AbstractCRVCE end
VCESymHCMapper = Dict([(:OLS, OLS()),(:HC1, HC1()),(:HC2, HC2()),(:HC3, HC3())])
VCESymCRMapper = Dict([(:HC1, CRVE1()), (:HC2, CRVE2()), (:HC3, CRVE3())])

## Methods packages should implement
function bread(obj::RegressionModel)
	mm = modelmatrix(obj)
	output = inv(factorize(Hermitian(mm.'mm)))
	return output
end
clusters(obj::RegressionModel) = Vector{Vector{Int64}}()

## API
function vcov(model::RegressionModel, estimator::Symbol)
    groups = clusters(model)
    if isempty(groups)
        @assert haskey(VCESymHCMapper, estimator) (@sprintf "%s is not a valid variance covariance estimator. Valid estimators are: OLS, HC1, HC2, HC3." estimator)
        output = vcov(model, get(VCESymHCMapper, estimator, OLS()))
    else
        @assert haskey(VCESymCRMapper, estimator) (@sprintf "%s is not a valid cluster robust variance covariance estimator. Valid estimators are: HC1, HC2, HC3." estimator)
        output = vcov(model, get(VCESymCRMapper, estimator, HC1()), groups)
    end
    return output
end

## Core Methods
function vcov(model::RegressionModel, estimator::OLS)
	û = residuals(model)
	output = û.'û / dof_residual(model) * bread(model)
end
function vcov(model::RegressionModel, estimator::AbstractHCVCE)
    mm = modelmatrix(model)
    Bread = bread(model)
    û = residuals(model)
    output = Bread * meat(estimator, mm, û) * Bread
    if isa(estimator, HC1)
        output .*= nobs(model) / dof_residual(model)
    end
    return output
end
function vcov(model::RegressionModel,
			  estimator::AbstractCRVCE,
			  Clusters::Vector{Vector{Vector{T}}}) where T <: Integer
    mm = modelmatrix(model)
    Bread = bread(model)
    û = residuals(model)
    R = eye(length(û))
    for dimension ∈ Clusters
        for level ∈ dimension
            for comparison ∈ dimension
                values = intersect(level, comparison)
                R[values, values] .= one(Float64)
            end
        end
    end
    if isa(estimator, CRVE2) | isa(estimator, CRVE3)
        h = Dict{Vector{Int64},Vector{Float64}}()
        output = zeros(length(û))
        for each ∈ 1:size(R, 1)
            group = find(R[:,each])
            output[each] = 1 - first(get!(h, group, hatvalues(mm[group,:]))[findfirst(equalto(each), group)])
        end
        if isa(estimator, CRVE2)
            û = û ./ sqrt.(output)
        elseif isa(estimator, CRVE3)
            û = û ./ output
        end
    end
    output = Bread * mm.' * (û * û.' .* R) * mm * Bread
    # g = length(unique(map(col -> find(R[:,col]), 1:size(R, 2))))
	g = minimum(length.(Clusters))
    rdf = dof_residual(model)
    output .*= g / (g - 1) * (nobs(model) - 1) / reduce(-, size(mm))
    return output, min(rdf, g - 1)
end

## Helper

gram(obj::AbstractMatrix) = obj'obj
function hatvalues(X::AbstractMatrix; Γ::AbstractMatrix = zeros(size(X, 2), size(X, 2)))
	sum((X * inv(factorize(Hermitian(X'X + Γ)))) .* X, 2)
end

meat(estimator::HC1, X::AbstractMatrix, û::AbstractVector) = gram(abs.(û) .* X)
meat(estimator::HC2, X::AbstractMatrix, û::AbstractVector) = gram(abs.(û) ./ sqrt.(broadcast(-, 1., hatvalues(X))) .* X)
meat(estimator::HC3, X::AbstractMatrix, û::AbstractVector) = gram(abs.(û) ./ broadcast(-, 1., hatvalues(X)) .* X)

end


text:
facts("dockerfile creation") do
end


text:
"""
    circumcenter_radius2(pts)

Calculate circumcenter and circumradius squared from points.

Based on https://github.com/hirani/pydec/blob/master/pydec/math/circumcenter.py and
http://www.ics.uci.edu/~eppstein/junkyard/circumcenter.html
"""
# TODO: try this with SMatrix?
function circumcenter_radius2(pts)
    if length(pts) == 1
        return pts[1], 0.0
    elseif length(pts) == 2
        return 0.5 * sum(pts), 0.25 * sum(abs2, pts[1] - pts[2])
    else
        P = to_matrix(pts)
        n = length(pts)

        A = ones(n + 1, n + 1)
        A[end, end] = 0.0
        A[1:(end - 1), 1:(end - 1)] .= 2 * P' * P

        b = vec(sum(abs2, P; dims=1))
        push!(b, 1.0)

        fact = lu(A; check=false)
        if issuccess(fact)
            x = fact \ b
            bary_coords = x[1:(end - 1)]

            center = P * bary_coords
            radius = P[:, 1] - center
            return center, sum(abs2, radius)
        else
            return fill(NaN, length(pts[1])), Inf
        end
    end
end

# TODO: could be faster, takes a long time to compile
function _build_dims!(dicts, triangulation, points, ::Val{D}, verbose) where {D}
    if verbose
        progbar = Progress(size(triangulation, 2); desc="Collecting $D-simplcies... ")
    end
    for face in eachcol(triangulation)
        for σ in IterTools.subsets(face, Val(D + 1))
            @assert issorted(σ; rev=true)
            σ_idx = index(σ)
            if !haskey(dicts[D + 1], σ_idx)
                _, σ_r2 = circumcenter_radius2(points[SVector(σ)])
                if isfinite(σ_r2)
                    dicts[D + 1][σ_idx] = σ_r2
                else
                    continue
                end
            end
            # Propagate birth time to facets.
            σ_r2 = dicts[D + 1][σ_idx]
            for i in 1:(D + 1)
                τ = TupleTools.deleteat(σ, i)
                τ_idx = index(τ)
                if haskey(dicts[D], τ_idx)
                    dicts[D][τ_idx] = min(σ_r2, dicts[D][τ_idx])
                elseif length(τ) > 1
                    τ_c, τ_r2 = circumcenter_radius2(points[SVector(τ)])
                    if sum(abs2, τ_c - points[σ[i]]) < τ_r2
                        dicts[D][τ_idx] = σ_r2
                    end
                end
            end
        end
        verbose && next!(progbar)
    end
end

function _fix_dim!(dicts, threshold, ::Val{D}, verbose) where {D}
    for (idx, birth) in dicts[D + 1]
        σ = Tuple(_vertices(idx, Val(D + 1)))
        σ_idx = index(σ)
        if D > 1
            for i in 1:(D + 1)
                τ = TupleTools.deleteat(σ, i)
                τ_idx = index(τ)
                dicts[D][τ_idx] = min(dicts[D][τ_idx], birth)
            end
        end
        corrected_birth = 2 * √dicts[D + 1][σ_idx]
        dicts[D + 1][σ_idx] = corrected_birth
    end
end

"""
    alpha_simplices(points[, verbose])

Collect all simplices and their birth times in alpha filtration.

Based on https://github.com/scikit-tda/cechmate/blob/master/cechmate/filtrations/alpha.py
"""
function alpha_simplices(points, threshold, verbose, ::Type{I}) where {I}
    @prog_print verbose "Building triangulation... "
    triangulation = I.(delaunay(to_matrix(points)))
    sort!.(eachcol(triangulation), rev=true, alg=InsertionSort)
    @prog_println verbose "done."

    largest_face = tuple(maximum(eachcol(triangulation))...)
    index_overflow_check(largest_face)

    dim = length(points[1])
    dicts = [Dict{I,Float64}() for _ in 0:dim]

    # Build the filtration
    for d in dim:-1:1
        _build_dims!(dicts, triangulation, points, Val(d), verbose)
    end
    for i in 1:length(points)
        dicts[1][i] = 0.0
    end
    if verbose
        progbar = Progress(dim; desc="Fixing birth times...     ")
    end
    # Make sure all simplices are born after their facets and sqrt the birth times.
    for d in dim:-1:1
        _fix_dim!(dicts, threshold, Val(d), verbose)
        verbose && next!(progbar)
    end

    return dicts
end

"""
    Alpha{I, P<:SVector} <: AbstractFiltration{I, Float64}

`Alpha` filtrations are filtrations of the Delaunay complex.

They have much fewer simplices than `Rips`, so they are efficient even with large datasets,
as long as their dimensionality is low.  What "low" means depends on the data, but this is
definitely a good choice for 3D or lower. For high dimensional data, filtration construction
may take a long time.

!!! note
    Unlike most implementations, this one uses circumdiameters instead of
    circumradii. This makes the scale of the results comparable to `Rips`. If you need
    radius based values, divide your data or the resulting interval endpoints by 2.

!!! warning
    This filtration uses [MiniQhull.jl](https://github.com/gridap/MiniQhull.jl). Please see
    the installation instructions if constructions cause errors. MiniQhull currently has
    problems running on Windows. See [this
    issue](https://github.com/gridap/MiniQhull.jl/issues/5) for more info.

# Constructors

* `Alpha(points; threshold, verbose)`: `points` should be a vector of `Tuple`s, `SVector`s
  or similar.
* `Alpha{I}(args...)`: `I` sets the size of integer used to represent simplices. Try using
  `I=Int128` if construction complains about overflow.

# Reference

Edelsbrunner, H. (1993, July). The union of balls and its dual shape. [In Proceedings of the
ninth annual symposium on Computational geometry
(pp. 218-231)](https://dl.acm.org/doi/abs/10.1145/160985.161139).

# Example

```jldoctest
julia> data = [(sin(t), cos(t), (t - π)^2) for t in range(0, 2π, length=101)[1:end-1]];

julia> alpha = Alpha(data)
Alpha{Int64, Float64}(nv=100)

julia> rips = Rips(data)
Rips{Int64, Float64}(nv=100, sparse=false)

julia> length(Ripserer.edges(alpha))
197

julia> length(Ripserer.edges(rips))
3613

julia> sort(ripserer(alpha)[2], by=persistence)[end]
[0.375, 2.01) with:
 birth_simplex: Simplex{1, Float64, Int64}
 death_simplex: Simplex{2, Float64, Int64}

julia> sort(ripserer(rips)[2], by=persistence)[end]
[0.375, 2.01) with:
 birth_simplex: Simplex{1, Float64, Int64}
 death_simplex: Simplex{2, Float64, Int64}
```
"""
struct Alpha{I,P<:SVector} <: AbstractCustomFiltration{I,Float64}
    dicts::Vector{Dict{I,Float64}}
    adj::SparseMatrixCSC{Bool,Int}
    threshold::Float64
    points::Vector{P}
end
function Alpha{I}(points; threshold=nothing, verbose=false) where {I}
    pts = SVector.(points)
    threshold = isnothing(threshold) ? 2radius(pts) : threshold
    dicts = alpha_simplices(pts, threshold, verbose, I)
    adj = _adjacency_matrix(dicts)
    return Alpha{I,eltype(pts)}(dicts, adj, threshold, pts)
end
function Alpha(points; kwargs...)
    return Alpha{Int}(points; kwargs...)
end

adjacency_matrix(alpha::Alpha) = alpha.adj
simplex_dicts(alpha::Alpha) = alpha.dicts
threshold(alpha::Alpha) = alpha.threshold

struct AlphaDist{A} <: AbstractMatrix{Float64}
    alpha::A
end

Base.size(ad::AlphaDist) = (nv(ad.alpha), nv(ad.alpha))
function Base.getindex(ad::AlphaDist, i::Integer, j::Integer)
    return Euclidean()(ad.alpha.points[i], ad.alpha.points[j])
end

distance_matrix(alpha::Alpha) = AlphaDist(alpha)


text:
mutable struct BNCache
  mean
  ivar
end

BNCache() = BNCache(nothing, nothing)

@inline _wsize(y) = (fill(1, ndims(y)-2)..., size(y)[end-1], 1)

# NOTE: CuDNN supports only 4D and 5D Tensors for BatchNorm Operations
# so reshape a 2D Tensor into 4D
batchnorm(g::DenseCuArray{T}, b::DenseCuArray{T}, x::DenseCuArray{T,2},
          running_mean::DenseCuArray{T}, running_var::DenseCuArray{T}, momentum;
          cache = nothing, alpha = T(1), beta = T(0),
          eps = T(1e-5), training = true) where T<:Union{Float32, Float64} =
  dropdims(batchnorm(g, b, reshape(x, 1, 1, size(x, 1), size(x, 2)), running_mean, running_var, momentum,
            cache = cache, alpha = alpha, beta = beta, eps = eps, training = training), dims = (1, 2))

function batchnorm(g::DenseCuArray{T}, b::DenseCuArray{T}, x::Union{DenseCuArray{T,4},DenseCuArray{T,5}},
                   running_mean::DenseCuArray{T}, running_var::DenseCuArray{T}, momentum;
                   cache = nothing, alpha = T(1), beta = T(0),
                   eps = T(1e-5), training = true) where T<:Union{Float32, Float64}
  cudnnBNForward!(similar(x), g, b, x, running_mean, running_var, momentum, cache = cache,
      alpha = alpha, beta = beta, eps = eps, training = training)
end

function cudnnBNForward!(y::DenseCuArray{T}, g::DenseCuArray{T}, b::DenseCuArray{T}, x::DenseCuArray{T},
                        running_mean::DenseCuArray{T}, running_var::DenseCuArray{T},
                        momentum; cache = nothing,
                        alpha = T(1), beta = T(0),
                        eps = T(1e-5), training = true) where T<:Union{Float32, Float64}
  dims = _wsize(x)
  if eps < CUDNN_BN_MIN_EPSILON
    # warn("eps ",eps," is too small for CuDNN so eps has been assigned the value ", CUDNN_BN_MIN_EPSILON)
    eps = CUDNN_BN_MIN_EPSILON
  end
  xd = cudnnTensorDescriptor(x)
  yd = cudnnTensorDescriptor(y)
  gd = cudnnTensorDescriptor(CUDNN_TENSOR_NCHW, cudnnDataType(T), Cint(length(dims)), dim4(dims,Val(CUDNN_TENSOR_NCHW)))

  if training

    if cache !== nothing
      mean = zeros(CuArray{T}, dims...)
      ivar = ones(CuArray{T}, dims...)
    else
      mean = CU_NULL
      ivar = CU_NULL
    end

    cudnnBatchNormalizationForwardTraining(handle(), CUDNN_BATCHNORM_SPATIAL, scalingParameter(T, alpha), scalingParameter(T, beta), xd, x, yd, y, gd, g, b, momentum, running_mean, running_var, eps, mean, ivar)

    if cache !== nothing
      cache.mean = mean
      cache.ivar = ivar
    end
  else
    cudnnBatchNormalizationForwardInference(handle(), CUDNN_BATCHNORM_SPATIAL, scalingParameter(T, alpha), scalingParameter(T, beta), xd, x, yd, y, gd, g, b, running_mean, running_var, eps)
  end
  return y
end

function ∇batchnorm(g::DenseCuArray{T}, b::DenseCuArray{T}, x::DenseCuArray{T, 2}, dy::DenseCuArray{T, 2},
           running_mean::DenseCuArray{T}, running_var::DenseCuArray{T}, momentum;
           cache = nothing, eps = T(1e-5), alpha = T(1),
           beta = T(0), training = true) where T<:Union{Float32, Float64}
  dg, db, dx = ∇batchnorm(g, b, reshape(x, 1, 1, size(x, 1), size(x, 2)), reshape(dy, 1, 1, size(dy, 1),
                          size(dy, 2)), running_mean, running_var, momentum, cache = cache, eps = eps,
                          alpha = alpha, beta = beta, training = training)
  (dg, db, dropdims(dx, dims = (1, 2)))
end

function ∇batchnorm(g::DenseCuArray{T}, b::DenseCuArray{T}, x::DenseCuArray{T}, dy::DenseCuArray{T},
                    running_mean::DenseCuArray{T}, running_var::DenseCuArray{T}, momentum;
                    cache = nothing, eps = T(1e-5), alpha = T(1),
                    beta = T(0), training = true) where T<:Union{Float32, Float64}
  dg = similar(g)
  db = similar(b)
  dx = similar(x)
  cudnnBNBackward!(dg, g, db, dx, x, dy, running_mean, running_var, T(momentum),
    training = training, cache = cache, eps = eps, alpha = alpha, beta = beta)
  (dg, db, dx)
end

function cudnnBNBackward!(dg::DenseCuArray{T}, g::DenseCuArray{T}, db::DenseCuArray{T},
                          dx::DenseCuArray{T}, x::DenseCuArray{T}, dy::DenseCuArray{T},
                          running_mean::DenseCuArray{T}, running_var::DenseCuArray{T},
                          momentum; cache = nothing, eps = T(1e-5),
                          alpha = T(1), beta = T(0),
                          dalpha = T(1), dbeta = T(0), training = true) where T<:Union{Float32, Float64}
  if training
    xd = cudnnTensorDescriptor(x)
    dyd = cudnnTensorDescriptor(dy)
    dxd = cudnnTensorDescriptor(dx)
    gd = cudnnTensorDescriptor(CUDNN_TENSOR_NCHW, cudnnDataType(T), Cint(length(_wsize(x))), dim4(_wsize(x),Val(CUDNN_TENSOR_NCHW)))
    if cache !== nothing
      mean, ivar = cache.mean, cache.ivar
      info("mean and ivar are fetched from the cache")
    else
      mean, ivar = CU_NULL, CU_NULL
    end

    if eps < CUDNN_BN_MIN_EPSILON
      eps = CUDNN_BN_MIN_EPSILON
    end

    cudnnBatchNormalizationBackward(handle(), CUDNN_BATCHNORM_SPATIAL, scalingParameter(T, alpha), scalingParameter(T, beta), scalingParameter(T, dalpha), scalingParameter(T, dbeta), xd, x, dyd, dy, dxd, dx, gd, g, dg, db, eps, mean, ivar)
  else
    ivar = 1 ./ sqrt.(reshape(running_var, _wsize(x)) .+ eps)
    dx .= dy .* reshape(g, _wsize(x)) .* ivar
    rdims = ((1:ndims(x)-2)..., ndims(x))
    dg .= vec(sum(dy .* (x .- reshape(running_mean, _wsize(x))) .* ivar, dims=rdims))
    db .= vec(sum(dy, dims=rdims))
  end
end


text:
###
  Name     : c9_32_iris.jl
  Book     : Hands-on Data Science with Anaconda )
  Publisher: Packt Publishing Ltd.
  Author   : Yuxing Yan and James Yan
  Date     : 4/6/2018
  email    : yany@canisius.edu
             paulyxy@hotmail.com
###

using Gadfly
using RDatasets
using Clustering
iris = dataset("datasets", "iris")
head(iris)
features=permutedims(convert(Array, iris[:,1:4]),[2, 1])
result=kmeans(features,3)
nameX="PetalLength"
nameY="PetalWidth"
assignments=result.assignments
plot(iris, x=nameX,y=nameY,color=assignments,Geom.point)


text:
using Pkg; for p in ("Knet","Colors","Images"); haskey(Pkg.installed(),p) || Pkg.add(p); end
using Knet, Colors, Images, Statistics
include(Pkg.dir("Knet","data","mnist.jl")) #MNIST data loader functions
global atype = gpu() >= 0 ? KnetArray{Float32} : Array{Float32}

#A generic MLP function with customizable activation functions
function mlp(w,x;p=0.0,activation=elu,outputactivation=sigm)
    for i=1:2:length(w)
        x = w[i]*dropout(mat(x),p) .+ w[i+1] # mat() used for flattening images to a vector.
        i<length(w)-1 && (x = activation.(x))
    end
    return outputactivation.(x) #output layer
end

global const 𝜀=Float32(1e-8)
D(w,x;p=0.0) = mlp(w,x;p=p)
G(w,z;p=0.0) = mlp(w,z;p=p)
𝑱d(𝗪d,x,Gz) = -mean(log.(D(𝗪d,x) .+ 𝜀)+log.((1+𝜀) .- D(𝗪d,Gz)))/2
𝑱g(𝗪g, 𝗪d, z) = -mean(log.(D(𝗪d,G(𝗪g,z)) .+ 𝜀))
𝒩(input, batch) = atype(randn(Float32, input, batch))  #SampleNoise

∇d  = grad(𝑱d) # Discriminator gradient
∇g  = grad(𝑱g) # Generator gradient

function initweights(hidden,input, output)
    𝗪 = Any[];
    x = input
    for h in [hidden... output]
        push!(𝗪, atype(xavier(h,x)), atype(zeros(h, 1))) #FC Layers weights and bias
        x = h
    end
    return 𝗪  #return model params
end

function generate_and_save(𝗪,number,𝞗;fldr="generations/")
    Gz = G(𝗪[1], 𝒩(𝞗[:ginp], number)) .> 0.5
    Gz = permutedims(reshape(Gz,(28,28,number)), (2,1,3))
    [save(fldr*string(i)*".png",Gray.(Gz[:,:,i])) for i=1:number]
end

#(if) train ? it updates model parameters : (else) it print losses
function runmodel(𝗪, data, 𝞗; dtst=nothing, optim=nothing, train=false, saveinterval=20)
    gloss = dloss = total=0.0;
    B = 𝞗[:batchsize]
    for i=1:(train ? 𝞗[:epochs] : 1)
        for (x,_) in data
            total+=2B
            Gz = G(𝗪[1], 𝒩(𝞗[:ginp], B)) #Generate Fake Images
            train ? update!(𝗪[2], ∇d(𝗪[2],x,Gz), optim[2]) : (dloss += 2B*𝑱d(𝗪[2], x, Gz))

            z=𝒩(𝞗[:ginp],2B) #Sample z from Noise
            train ? update!(𝗪[1], ∇g(𝗪[1], 𝗪[2], z), optim[1]) : (gloss += 2B*𝑱g(𝗪[1],𝗪[2],z))
        end
        train ? runmodel(𝗪, dtst, 𝞗; train=false) : println((gloss/total, dloss/total))
        i % saveinterval == 0 && generate_and_save(𝗪, 100, 𝞗)  # save 10 images
    end
end

function main()
    𝞗=Dict(:batchsize=>32,:epochs=>80,:ginp=>256,:genh=>[512],:disch=>[512],:optim=>Adam,:lr=>0.0002);
    xtrn,ytrn,xtst,ytst = mnist()
    global dtrn = minibatch(xtrn, ytrn, 𝞗[:batchsize]; xtype=atype)
    global dtst = minibatch(xtst, ytst, 𝞗[:batchsize]; xtype=atype)
    𝗪  = (𝗪g,𝗪d)  = initweights(𝞗[:genh], 𝞗[:ginp], 784), initweights(𝞗[:disch], 784, 1)
    𝚶  = (𝚶pg,𝚶pd) = optimizers(𝗪g, 𝞗[:optim]; lr=𝞗[:lr]), optimizers(𝗪d,𝞗[:optim]; lr=𝞗[:lr])
    generate_and_show(𝗪,100,𝞗)
    runmodel(𝗪, dtst, 𝞗; optim=𝚶, train=false) # initial losses
    runmodel(𝗪, dtrn, 𝞗; optim=𝚶, train=true, dtst=dtst) # training
    𝗪,𝚶,𝞗,(dtrn,dtst)    # return weights,optimizers,options and dataset
end
main() #enjoy!


text:
export ContinuousEuclid

"""
$(TYPEDEF)

Most basic continuous scalar variable in a `::DFG.AbstractDFG` object.

DevNotes
- TODO Consolidate with ContinuousEuclid{1}
"""
@defVariable ContinuousScalar 1 (:Euclid,)

"""
$(TYPEDEF)

Continuous variable of dimension `.dims` on manifold `.manifolds`.
"""
struct ContinuousMultivariate{T1 <: Tuple} <: InferenceVariable
  dims::Int
  manifolds::T1
  # ContinuousMultivariate{T}() where {T} = new()
  # ContinuousMultivariate{T}(x::Int; manifolds::T=(:Euclid,)) where {T <: Tuple} = new(x, manifolds)
end

function ContinuousMultivariate(x::Int;
                                manifolds::T1=(:Euclid,)  )  where {T1 <: Tuple}
  #
  maniT = length(manifolds) < x ? ([manifolds[1] for i in 1:x]...,) : manifolds
  ContinuousMultivariate{typeof(maniT)}(x, maniT)
end

"""
    ContinuousEuclid{N}
Continuous Euclidean variable of dimension `N`.
"""
struct ContinuousEuclid{N} <: InferenceVariable end

ContinuousEuclid(x::Int) = ContinuousEuclid{x}()

getDimension(::ContinuousEuclid{N}) where N = N::Int
getManifolds(::ContinuousEuclid{N}) where N = ntuple(i -> :Euclid, N)


text:
module EmissionExpressions

using Symbolics

export parameters, spectrum_variable, f_actual_experimental_matlab_fitting, f_experimental_paper_AQT
export f_CUI_background_paper

parameters = @variables γ_QD γ_SP ω_QD ω_SP Ω_R
spectrum_variable = Symbolics.variable(:ω)

ex_actual_experimental_matlab_fitting() = let (a,b,c,d,G) = parameters, x_ev = spectrum_variable
    A = im*(x_ev-c + (c-d));
    B = (a+b)/4 - im*(c-d)/2 - im*(x_ev-c);
#    G = sqrt((e/2)^2 + ((c-d)/2)^2 - ((a-b)/4)^2);
    a/π * abs2((b/2-A)/(B^2+G^2))
end

f_expr_actual_experimental_matlab_fitting() = build_function(ex_actual_experimental_matlab_fitting(),γ_QD,γ_SP,ω_QD,ω_SP,Ω_R,spectrum_variable)

"""
    f_actual_experimental_matlab_fitting() -> f(γ_QD,γ_SP,ω_QD,ω_SP,Ω_R,ω)

Returns a function representing the fit code in https://github.com/mollyamay/strong-coupling-modeling/blob/b95b4dde3cbabd8a8d1dfe4266dbfb9a2ab0683c/f.m

In this expression, we assume that `G` corresponds to Ω_R (see source code)
Also, we removed a scaling parameter and an offset parameter used in the original fit code.
This means that we are free to offset and scale as necessary.
"""
f_actual_experimental_matlab_fitting() = eval(f_expr_actual_experimental_matlab_fitting())

ex_experimental_paper_AQT() = let (γ_QD,γ_SP,ω_QD,ω_SP,Ω_R) = parameters, ω = spectrum_variable
    δ = ω_SP - ω_QD;
    γ_QD/(2π) * abs2( (γ_SP/2 - im*(ω - ω_SP)) / (((γ_SP + γ_QD)/4 + im*(δ)/2 - im*(ω-ω_QD))^2 + Ω_R^2) )
end

f_expr_experimental_paper_AQT() = build_function(ex_experimental_paper_AQT(),γ_QD,γ_SP,ω_QD,ω_SP,Ω_R,spectrum_variable)

"""
    f_experimental_paper_AQT() -> f(γ_QD,γ_SP,ω_QD,ω_SP,Ω_R,ω)
Returns function representing equation 12 of:
May, M. A., Fialkow, D., Wu, T., Park, K., Leng, H., Kropp, J. A., Gougousi, T., Lalanne, P., Pelton, M., & Raschke, M. B. (2020). Nano‐Cavity QED with Tunable Nano‐Tip Interaction. Advanced Quantum Technologies, 3(2), 1900087. https://doi.org/10.1002/qute.201900087
"""
f_experimental_paper_AQT() = eval(f_expr_experimental_paper_AQT())

#in the paper γ = γ_QD/2, κ =γ_SP/2 are defined as half of the lifetimes
ex_CUI_background_paper() = let γ = γ_QD/2, κ = γ_SP/2, ω_0 = ω_QD, ω_c = ω_SP, g = Ω_R, ω = spectrum_variable
    Δ = ω_0 -ω_c;
    K = κ + γ;
    Γ = κ - γ;
    γ/π * abs2( (κ - im*(ω - ω_c))/ (  (K/2 - im*(2ω - ω_0 - ω_c)/2)^2 +g^2    )   )
end

f_expr_CUI_background_paper() = build_function(ex_CUI_background_paper(),γ_QD,γ_SP,ω_QD,ω_SP,Ω_R,spectrum_variable)

"""
    f_CUI_background_paper() -> f(γ_QD,γ_SP,ω_QD,ω_SP,Ω_R,ω)

Returns function representing equation 14 of:
Cui, G., & Raymer, M. G. (2006). Emission spectra and quantum efficiency of single-photon sources in the cavity-QED strong-coupling regime. Physical Review A, 73(5), 053807. https://doi.org/10.1103/PhysRevA.73.053807

"""
f_CUI_background_paper() = eval(f_expr_CUI_background_paper())

end


text:
"""
    te_cond(source::AbstractArray{<:Real, 1},
        response::AbstractArray{<:Real, 1},
        cond::AbstractArray{<:Real, 1},
        k::Int, l::Int, m::Int, n::Int;
        η = 1, τ = 1,
        estimator = VisitationFrequency(),
        n_subdivs = 1,
        b = 2)

# Conditional TE with default discretization scheme(s)

Calculate transfer entropy from `source` to `response` conditioned on
`cond` using the provided
`estimator` on a rectangular discretization of a `k + l + m + n`-dimensional
delay embedding of the input data, using an embedding delay of `τ` across
all embedding components. `η` is the prediction lag.

## Arguments
- **`source`**: The source data series.
- **`target`**: The target data series.
- **`cond`**: The conditional data series.
- **`k`**: The dimension of the ``T_{f}`` component of the embedding.
- **`l`**: The dimension of the ``T_{pp}`` component of the embedding.
- **`m`**: The dimension of the ``S_{pp}`` component of the embedding.
- **`n`**: The dimension of the ``C_{pp}`` component of the embedding.

## Keyword arguments

- **`τ`**: The embedding lag. Default is `τ = 1`.
- **`η`**: The prediction lag. Default is `η = 1`.
- **`estimator`**: The transfer entropy estimator to use. The default
    is `VisitationFrequency()`.
- **`n_subdivs`**: The number of different partitions of varying coarseness
    to compute TE over, as described below. Default is `n_subdivs = 2`.
    (this way, TE is computed over two separate partitions). Unless `n_subdivs = 0`,
    make sure that `n_subdivs` is the same across analyses if they
    are to be compared. T TE
- **`b`**: Base of the logarithm. The default (`b = 2`) gives the TE in bits.

## More about the embedding

To compute transfer entropy, we need an appropriate delay embedding
of `source` (``S``), `target` (``T``) and `cond` (``C``). For convenience, define

```math
\\begin{align}
T_f^{(k)} &= \\{(T(t+\\eta_k), \\ldots, T(t+\\eta_2), T(t+\\eta_1)) \\} \\\\
T_{pp}^{(l)} &= \\{ (T(t), T(t-\\tau_1), T(t-\\tau_2), \\ldots, T(t - \\tau_{l - 1})) \\} \\\\
S_{pp}^{(m)} &= \\{ (S(t), S(t-\\tau_1), S(t-\\tau_2), \\ldots, S(t-\\tau_{m - 1})) \\} \\\\
C_{pp}^{(n)} &= \\{ (C(t), C(t-\\tau_1), C(t-\\tau_2), \\ldots, C(t-\\tau_{n - 1})) \\}
\\end{align}
```

where ``T_f`` denotes the `k`-dimensional set of vectors furnishing the future states of ``T``,
``T_{pp}`` denotes the `l`-dimensional set of vectors furnishing the past and present states of ``T``,
,``S_{pp}`` denotes the `m`-dimensional set of vectors furnishing the past and present of ``S``,
and ``C_{pp}`` denotes the `n`-dimensional set of vectors furnishing the past and present of ``C``.
``\\eta`` is the prediction lag. This convenience function uses ``\\tau_1`` = `τ`,
``\\tau_2`` = `2*τ`, ``\\tau_3`` = `3*τ`, and so on.

Combined, we get the generalised embedding ``\\mathbb{E} = (T_f^{(k)}, T_{pp}^{(l)}, S_{pp}^{(m)}, C_{pp}^{(n)})``,
which is discretized as described below.

## More about discretization

To compute TE, we coarse-grain the `k+l+m+n`-dimensional generalised embedding
``\\mathbb{E} = (T_f^{(k)}, T_{pp}^{(l)}, S_{pp}^{(m)}, C_{pp}^{(n)})`` into hyperrectangular boxes.
The magnitude of the TE may be biased by the particular choice of binning scheme,
so we compute TE across a number of different box sizes, determined as follows.

Let ``L`` be the number of observations in ``S`` (and ``T`` and ``C``). The coarsest box size is
determined by subdiving the i-th coordinate axis into
``N = ceiling(L^\\frac{1}{k + l + m + n + 1})`` intervals of equal lengths,
resulting in a box size of ``|max(dim_{i}) - min(dim_{i})|/N``. The next box size
is given by ``|max(dim_{i}) - min(dim_{i})|/(N+1)``, then
``|max(dim_{i}) - min(dim_{i})|/(N+2)``, and so on, until the finest box
size which is given by ``|max(dim_{i}) - min(dim_{i})|/(N+N_{subdivs})``.

## Transfer entropy computation

TE is then computed as

```math
\\begin{align}
TE_{S \\rightarrow T} = \\int_{\\mathbb{E}} P(T_f, T_{pp}, S_{pp}, C_{pp}) \\log_{b}{\\left(\\frac{P(T_f | T_{pp}, S_{pp}, C_{pp})}{P(T_f | T_{pp}, C_{pp})}\\right)}
\\end{align}
```
using the provided `estimator` (default = `VisitationFrequency`) for
each of the discretizations. A vector of the TE estimates for each discretization
is returned.
"""
function te_cond(source::AbstractArray{<:Real, 1},
        response::AbstractArray{<:Real, 1},
        cond::AbstractArray{<:Real, 1},
        k::Int, l::Int, m::Int, n::Int;
        η = 1, τ = 1,
        estimator::TransferEntropyEstimator = VisitationFrequency(),
        n_subdivs = 2,
        b = 2)

    k + l + m + n >= 4 || throw(ArgumentError("`dim = k + l + m + n` must be 4 or higher for conditional TE"))

    pts, vars = te_embed(source, response, cond, k, l, m, n, η = η, τ = τ)

    # Determine appropriate binnings from time series length (roughly according to
    # Krakovska et al. (2018)'s recommendations)
    L = length(source)

    # one subdivision coarser than Krakovska as the coarest partition
    nbins_coarsest = ceil(Int, L^(1/(k + l + m + n + 1))) - 1

    # with the default n_subdivs = 2, we cover both the Krakovska recommendation,
    # as well as one finer bin size and one coarser bin size.
    nbins_finest = nbins_coarsest + n_subdivs
    binning_scheme = map(n-> RectangularBinning(n), nbins_coarsest:nbins_finest)

    # Compute TE over different partitions
    # ====================================
    tes = map(binscheme -> transferentropy(pts, vars, binscheme, estimator, b = b), binning_scheme)

    return tes
end

"""
    te_cond(source::AbstractArray{<:Real, 1},
        response::AbstractArray{<:Real, 1},
        cond::AbstractArray{<:Real, 1},
        k::Int, l::Int, m::Int, n::Int,
        binning_scheme::Vector{RectangularBinning};
        η = 1, τ = 1,
        estimator = VisitationFrequency(),
        b = 2)

# Conditional TE with default discretization scheme(s

Calculate transfer entropy from `source` to `response` conditioned on
`cond` using the provided
`estimator` on a rectangular discretization of a `k + l + m + n`-dimensional
delay embedding of the input data, using an embedding delay of `τ` across
all embedding components. `η` is the prediction lag.

## Arguments
- **`source`**: The source data series.
- **`target`**: The target data series.
- **`cond`**: The conditional data series.
- **`k`**: The dimension of the ``T_{f}`` component of the embedding.
- **`l`**: The dimension of the ``T_{pp}`` component of the embedding.
- **`m`**: The dimension of the ``S_{pp}`` component of the embedding.
- **`n`**: The dimension of the ``C_{pp}`` component of the embedding.
- **`binning_scheme`**: The binning scheme(s) used to construct the partitions
    over which TE is computed. Must be either one or several instances of
    `RectangularBinning`s (provided as a vector). TE is computed for each
    of the resulting partitions.

## Keyword arguments

- **`τ`**: The embedding lag. Default is `τ = 1`.
- **`η`**: The prediction lag. Default is `η = 1`.
- **`estimator`**: The transfer entropy estimator to use. The default
    is `VisitationFrequency()`.
- **`b`**: Base of the logarithm. The default (`b = 2`) gives the TE in bits.

## More about the embedding

To compute transfer entropy, we need an appropriate delay embedding
of `source` (``S``), `target` (``T``) and `cond` (``C``). For convenience, define

```math
\\begin{align}
T_f^{(k)} &= \\{(T(t+\\eta_k), \\ldots, T(t+\\eta_2), T(t+\\eta_1)) \\} \\\\
T_{pp}^{(l)} &= \\{ (T(t), T(t-\\tau_1), T(t-\\tau_2), \\ldots, T(t - \\tau_{l - 1})) \\} \\\\
S_{pp}^{(m)} &= \\{ (S(t), S(t-\\tau_1), S(t-\\tau_2), \\ldots, S(t-\\tau_{m - 1})) \\} \\\\
C_{pp}^{(n)} &= \\{ (C(t), C(t-\\tau_1), C(t-\\tau_2), \\ldots, C(t-\\tau_{n - 1})) \\}
\\end{align}
```

where ``T_f`` denotes the `k`-dimensional set of vectors furnishing the future states of ``T``,
``T_{pp}`` denotes the `l`-dimensional set of vectors furnishing the past and present states of ``T``,
,``S_{pp}`` denotes the `m`-dimensional set of vectors furnishing the past and present of ``S``,
and ``C_{pp}`` denotes the `n`-dimensional set of vectors furnishing the past and present of ``C``.
``\\eta`` is the prediction lag. This convenience function uses ``\\tau_1`` = `τ`,
``\\tau_2`` = `2*τ`, ``\\tau_3`` = `3*τ`, and so on.

Combined, we get the generalised embedding ``\\mathbb{E} = (T_f^{(k)}, T_{pp}^{(l)}, S_{pp}^{(m)}, C_{pp}^{(n)})``,
which is discretized as described below.

## More about discretization

The discretization scheme must be either a single `RectangularBinning` instance, or a vector of
`RectangularBinning` instances. Run `?RectangularBinning` after loading `CausalityTools` for
details.

## Transfer entropy computation

TE is then computed as

```math
\\begin{align}
TE_{S \\rightarrow T} = \\int_{\\mathbb{E}} P(T_f, T_{pp}, S_{pp}, C_{pp}) \\log_{b}{\\left(\\frac{P(T_f | T_{pp}, S_{pp}, C_{pp})}{P(T_f | T_{pp}, C_{pp})}\\right)}
\\end{align}
```
using the provided `estimator` (default = `VisitationFrequency`) for
each of the discretizations. A vector of the TE estimates for each discretization
is returned.
"""
function te_cond(source::AbstractArray{<:Real, 1},
        response::AbstractArray{<:Real, 1},
        cond::AbstractArray{<:Real, 1},
        k::Int, l::Int, m::Int, n::Int,
        binning_scheme::Vector{RectangularBinning};
        η = 1, τ = 1,
        estimator = VisitationFrequency(),
        b = 2)

    k + l + m + n >= 4 || throw(ArgumentError("`dim = k + l + m + n` must be 4 or higher for conditional TE"))

    pts, vars = te_embed(source, response, cond, k, l, m, n, η = η, τ = τ)

    # Compute TE over the partitions constructed from the provided binning schemes
    # ====================================
    tes = [transferentropy(pts, vars, bs, estimator, b = b) for bs in binning_scheme]

    return tes
end

function te_cond(source::AbstractArray{<:Real, 1},
        response::AbstractArray{<:Real, 1},
        cond::AbstractArray{<:Real, 1},
        k::Int, l::Int, m::Int, n::Int,
        binning_scheme::RectangularBinning;
        η = 1, τ = 1,
        estimator = VisitationFrequency(),
        b = 2)

    k + l + m + n >= 4 || throw(ArgumentError("`dim = k + l + m + n` must be 4 or higher for conditional TE"))

    pts, vars = te_embed(source, response, cond, k, l, m, n, η = η, τ = τ)

    # Compute TE over the partitions constructed from the provided binning schemes
    # ===========================================================================
    te = transferentropy(pts, vars, binning_scheme, estimator, b = b)

    return te
end

export te_cond


text:
using WebIO
using JSExpr

abstract type AbstractGameWindow end

# general gui for event based games like go, chess etc
# use interface.js and interface.css as well
mutable struct GameWindow <: AbstractGameWindow
 w::WebIO.Scope
 action::WebIO.Observable
 board::WebIO.Observable
 message::WebIO.Observable
 config
end

function GameWindow(imports; actionInit=Dict(), boardInit=Dict(), messageInit="", config=Dict())
 w = Scope(imports=imports)
 action = Observable(w, "action", actionInit)
 board = Observable(w, "board", boardInit)
 message = Observable(w, "message", messageInit)

 GameWindow(w, action, board, message, config)
end

# setup
WebIO.onimport(g::AbstractGameWindow, f = @js () -> begin
 game.config = $(g.config)

 function action(a)
  $(g.action)[] = a
 end
 game.action = action

 __init__()

end) = WebIO.onimport(g.w, f)

# for rendering
(g::GameWindow)(h) = (g.w)(h)
WebIO.render(g::AbstractGameWindow) =
 WebIO.render(g.w)

# input(g) do x -> ... end
# ( to get user input after game.action is invoked )
oninput(f, g::AbstractGameWindow) =
 on(f, g.action)

# to update screen
function update!(g::AbstractGameWindow, x)
 g.board[] = x
end

# to set message
function message(g::AbstractGameWindow, x)
 g.message[] = x
end

# to setup js subscribers and the script to call init
function setup(g::AbstractGameWindow)
 onjs(g.board, @js x -> begin
  game.update(x)
 end)

 onjs(g.message, @js x-> begin
  game.showMsg(x)
 end)

 onimport(g)
end


text:
###############################################################################
# Trace View/Accessor Functions
###############################################################################
export @tracevals, @traceitems, @tracevalsdic, @tracesdf, @tracesdic,
@plotvals, @savetraces, @loadtraces, setsnaps

"""
`setsnaps(traces::Vector{TraceItem})`

Set the global RickTracy.happysnaps used by default in traceviews to traces.
Useful when loading traces from file
"""
function setsnaps(traces::Vector{TraceItem})
    global happysnaps = traces
end

"""
`tracevals(query, snaps=happysnaps)`
Get the values of the `snaps` that match `query`
"""
tracevals(query=Dict{Symbol, Any}(), snaps=happysnaps) = begin
    pluck(traceitems(query, snaps), :val)
end

"""
`traceitems(query, snaps=happysnaps)`
Get all TraceItems from `snaps` that match the `query`
"""
traceitems(query=Dict{Symbol, Any}(), snaps=happysnaps) = begin
    filterquery(query, snaps)
end

"""
`tracevalsdic(query, snaps=happysnaps)`
Get a Dict{String, Vector{Any}} mapping expressions => a vector of the values
they held at your tracepoints that match `query`
"""
tracevalsdic(query=Dict{Symbol, Any}(), snaps=happysnaps) = begin
    res = DefaultDict(String, Vector{Any}, Vector{Any})
    for si in traceitems(query, snaps)
        push!(res[si.exprstr], si.val)
    end
    res
end

"""
`tracesdic(traces=happysnaps)`
Get a Dict{Symbol, Vector{Any}} mapping attributes of your TraceItems
(:location, :exprstr, :val, :ts) => the vector of values of that attribute for
traces that match `query`
"""
tracesdic(query=Dict{Symbol, Any}(), snaps=happysnaps) = begin
    fulldic = Dict{Symbol, Vector{Any}}()
    traces = traceitems(query, snaps)
    foreach(fieldnames(TraceItem)) do attr
        fulldic[attr] = pluck(traces, attr)
    end
    fulldic
end

"""
`tracesdf(query, snaps=happysnaps)`
Get a DataFrame from your traces that match `query`
"""
tracesdf(query=Dict{Symbol, Any}(), snaps=happysnaps) = begin
    @eval import DataFrames
    DataFrames.DataFrame(tracesdic(query, snaps))
end

"""
Plot the values your expr string took
"""
plotvals(query=Dict{Symbol, Any}(), snaps=happysnaps) = begin
    dicsnaps = tracevalsdic(query, snaps)
    #for some reason transposing a vector of strings throws a depwarn
    #so we use hcat with '...' instead of collect(keys(dicsnaps))'
    @eval import Plots
    Plots.plot(collect(values(dicsnaps)),
        label=reshape(collect(keys(dicsnaps)), 1, length(dicsnaps)))
end

function parse_query_args(exprs; kwargspec=trace_kwargspec, no_defaults=true, kwargs...)
    query, kwextras, exprs, arginfo =
        kwargparse(kwargspec, exprs; no_defaults=no_defaults, kwargs...)
    !isempty(exprs) && (query[:exprstr] = map(string, exprs))
    query, kwextras, exprs, arginfo
end

getquery(exprs) = parse_query_args(exprs)[1]

function getpath_varname_query(exprs; path="traces.jld", varname="traces")
    kwargspec = filter(trace_kwargspec) do k,v; k != :_and_the_rest end
    query, kwextras, exprs, arginfo =
        parse_query_args(exprs; kwargspec=kwargspec)
    path    = get(kwextras, :path,    path)
    varname = get(kwextras, :varname, varname)
    path, varname, query
end

"""
`@tracevals [loc=some_location] [expr1] [expr2] [expr3] ...`
returns a vector of all the values that variables/expressions
took at all tracepoints, optionally limited to a particular tracepoint location or expression(s).
"""
macro tracevals(exprs...)
    query = getquery(exprs)
    :(RickTracy.tracevals($query))
end

"""
`@traceitems [loc=any_location] [expr1] [expr2] [expr3] ...`
returns a Vector of all the raw `TraceItem ` snaps, optionally
limited to a particular tracepoint location or expression(s).
"""
macro traceitems(exprs...)
    query = getquery(exprs)
    :(RickTracy.traceitems($query))
end

"""
`@tracevalsdic [loc=a_location] [expr1] [expr2] [expr3] ...`: returns a Dict mapping expressions=>values the variable/expression took, optionally
limited to a particular tracepoint location or expression(s).
"""
macro tracevalsdic(exprs...)
    query = getquery(exprs)
    :(RickTracy.tracevalsdic($query))
end

"""
`@tracesdf [loc=any_location] [expr1] [expr2] [expr3] ...`: returns a DataFrame with fields `:location`, `:exprstr`, `:val`, `:lcount`, `:ts`, with each row holding a snap of one expression, optionally limit to a particular tracepoint location or expression(s).
"""
macro tracesdf(exprs...)
    query = getquery(exprs)
    :(RickTracy.tracesdf($query))
end

"""
`@tracesdic [loc=any_location] [expr1] [expr2] [expr3] ...`: returns a Dict{Symbol, Vector{Any}} with keys `:location`, `:exprstr`, `:val`, `:lcount`, `:ts`, with values being a vector of the value for that field for all snaps, optionally limited to a particular tracepoint location or expression(s).
"""
macro tracesdic(exprs...)
    query = getquery(exprs)
    :(RickTracy.tracesdic($query))
end

"""
`@plotvals [loc=any_location] [expr1] [expr2] [expr3] ...`: Make a plot of all values that `@tracevals` would return for the same arguments
N.b. will break if any values are non-numeric and probably in lots of other
cases too. You'll probably want to use the optional filter by the
specified location, and the given expression(s).
"""
macro plotvals(exprs...)
    query = getquery(exprs)
    :(RickTracy.plotvals($query))
end

"""
`@savetraces [path="$(pwd())/traces.jld"] [loc=any_location] [expr1] [expr2] [expr3] ...`
Saves the raw `TraceItem` snaps to path, optionally
limited to a particular tracepoint location or expression(s).
"""
macro savetraces(exprs...)
    path, varname, query = getpath_varname_query(exprs)
    :(RickTracy.savetraces(; traces=RickTracy.traceitems($query),
        path=$path, varname=$varname))
end

macro loadtraces(exprs...)
    path, varname, query = getpath_varname_query(exprs)
    :(RickTracy.loadtraces(; path=$path, varname=$varname))
end

max_stars_repo_name: JobJob/RickTracy


text:
using GuessworkQuantumSideInfo: ket, dm, BB84_states, iid_copies, randdm, guesswork,
                                guesswork_upper_bound, guesswork_MISDP
using StableRNGs, UnPack
using Mosek, Gurobi, Pajarito, Cbc, SCS

if !@isdefined(GRB_ENV)
    const GRB_ENV = Gurobi.Env()
end

function Y_states(::Type{T}) where {T}
    return [begin
                θⱼ = T(j) * (T(2) * π) / T(3)
                dm(cos(θⱼ) * ket(T, 1, 2) + sin(θⱼ) * ket(T, 2, 2))
            end
            for j in 1:3]
end

function two_random_qubits()
    rng = StableRNG(746)
    return [randdm(rng, Float64, 2) for _ in 1:2]
end

function three_random_qubits()
    rng = StableRNG(514)
    return [randdm(rng, Float64, 2) for _ in 1:3]
end

function two_random_qutrits()
    rng = StableRNG(948)
    return [randdm(rng, Float64, 3) for _ in 1:2]
end

function three_random_qutrits()
    rng = StableRNG(219)
    return [randdm(rng, Float64, 3) for _ in 1:3]
end

function iid_problem(states, n)
    J = length(states)
    p = ones(T, J^n) / J^n
    ρBs = iid_copies(states, n)
    (p=p, ρBs = ρBs)
end

T = Float64
problems = []
for n in 1:2
    append!(problems, [
                (iid_problem(two_random_qubits(), n)..., numeric_type=T, problem="2qubits($n)"),
                (iid_problem(three_random_qubits(), n)..., numeric_type=T, problem="3qubits($n)"),
                (iid_problem(Y_states(T), n)..., numeric_type=T, problem="Y($n)"),
                (iid_problem(BB84_states(T), n)..., numeric_type=T, problem="BB84($n)"),
                (iid_problem(two_random_qutrits(), n)..., numeric_type=T, problem="2qutrits($n)"),
                (iid_problem(three_random_qutrits(), n)..., numeric_type=T, problem="3qutrits($n)"),
                ])
end

algos = []

for (misdp_solver, settings) in [(PajaritoSolver(cont_solver=MosekSolver(LOG=0),
                     mip_solver=GurobiSolver(GRB_ENV, OutputFlag=0, IntFeasTol=1e-9, FeasibilityTol=1e-8, MIPGap=0),
                     mip_solver_drives=false, rel_gap = 1e-5,
                     log_level=0), "Pajarito(Mosek, Gurobi, MSD=false)"),
                     (PajaritoSolver(cont_solver=MosekSolver(LOG=0),
                     mip_solver=GurobiSolver(GRB_ENV, OutputFlag=0, IntFeasTol=1e-9, FeasibilityTol=1e-8, MIPGap=1e-5),
                     mip_solver_drives=true, rel_gap = 1e-5,
                     log_level=0), "Pajarito(Mosek, Gurobi, MSD=true)"),
     (PajaritoSolver(cont_solver=SCSSolver(verbose=false, eps=1e-6), mip_solver=CbcSolver(loglevel=0, integerT=1e-8),
                     mip_solver_drives=false, rel_gap = 1e-5,
                     log_level=0), "Pajarito(SCS, Cbc, MSD=false)")]
    push!(algos,
        (; f = (prob -> @timed guesswork_MISDP(prob.p, prob.ρBs, size(first(prob.ρBs),1); solver=misdp_solver, verbose=false).optval), algo = "MISDP_dB",
                                                                                settings))
    push!(algos,
          (; f = (prob -> @timed guesswork_MISDP(prob.p, prob.ρBs, size(first(prob.ρBs), 1)^2; solver=misdp_solver, verbose=false).optval), algo = "MISDP",
                                                                                settings))

    push!(algos,
        (; f = (prob -> @timed guesswork_MISDP(prob.p, prob.ρBs, size(first(prob.ρBs), 1);
        solver=misdp_solver, verbose=false).optval), algo = "MISDP (dB^2)", settings))
end

for (sdp_solver, settings) in
    [((() -> MosekSolver(LOG=0)), "Mosek"), ((() -> SCSSolver(verbose=0, eps=1e-6)), "SCS")]
    push!(algos,
          (; f = (prob -> @timed guesswork(prob.p, prob.ρBs; solver=sdp_solver(), verbose=false).optval), algo = "SDP", settings))

    push!(algos,
          (; f = (prob -> @timed guesswork(prob.p, prob.ρBs; dual=true, solver=sdp_solver(), verbose=false).optval), algo = "dual_SDP", settings))

    for max_time in (20, 60, 60*4)
        push!(algos,
          (; f = (prob -> @timed guesswork_upper_bound(prob.p, prob.ρBs; make_solver=sdp_solver, verbose=false, max_time).optval), algo = "guesswork_upper_bound(max_time=$max_time)", settings))
    end
end


text:
# Stream management

export
    CuStream, CuDefaultStream, priority, priority_range, synchronize

mutable struct CuStream
    handle::CUstream
    ctx::CuContext
end

Base.unsafe_convert(::Type{CUstream}, s::CuStream) = s.handle

Base.:(==)(a::CuStream, b::CuStream) = a.handle == b.handle
Base.hash(s::CuStream, h::UInt) = hash(s.handle, h)

@enum_without_prefix CUstream_flags_enum CU_

"""
    CuStream(; flags=STREAM_DEFAULT, priority=nothing)

Create a CUDA stream.
"""
function CuStream(; flags::CUstream_flags=STREAM_DEFAULT,
                    priority::Union{Nothing,Integer}=nothing)
    handle_ref = Ref{CUstream}()
    if priority === nothing
        cuStreamCreate(handle_ref, flags)
    else
        priority in priority_range() || throw(ArgumentError("Priority is out of range"))
        cuStreamCreateWithPriority(handle_ref, flags, priority)
    end

    ctx = CuCurrentContext()
    obj = CuStream(handle_ref[], ctx)
    finalizer(unsafe_destroy!, obj)
    return obj
end

function unsafe_destroy!(s::CuStream)
    if isvalid(s.ctx)
        cuStreamDestroy(s)
    end
end

"""
    CuDefaultStream()

Return the default stream.
"""
@inline CuDefaultStream() = CuStream(convert(CUstream, C_NULL), CuContext(C_NULL))

"""
    synchronize(s::CuStream)

Wait until a stream's tasks are completed.
"""
synchronize(s::CuStream) = cuStreamSynchronize(s)

"""
    query(s::CuStream)

Return `false` if a stream is busy (has task running or queued)
and `true` if that stream is free.
"""
function query(s::CuStream)
    res = unsafe_cuStreamQuery(s)
    if res == ERROR_NOT_READY
        return false
    elseif res == SUCCESS
        return true
    else
        throw_api_error(res)
    end
end

"""
    priority_range()

Return the valid range of stream priorities as a `StepRange` (with step size  1). The lower
bound of the range denotes the least priority (typically 0), with the upper bound
representing the greatest possible priority (typically -1).
"""
function priority_range()
    least_ref = Ref{Cint}()
    greatest_ref = Ref{Cint}()
    cuCtxGetStreamPriorityRange(least_ref, greatest_ref)
    step = least_ref[] < greatest_ref[] ? 1 : -1
    return least_ref[]:Cint(step):greatest_ref[]
end

"""
    priority_range(s::CuStream)

Return the priority of a stream `s`.
"""
function priority(s::CuStream)
    priority_ref = Ref{Cint}()
    cuStreamGetPriority(s, priority_ref)
    return priority_ref[]
end


text:
"""
Original publication:
A. Ghasemi, H. R. Rabiee, M. Fadaee, M. T. Manzuri,
and M. H. Rohban. Active learning from positive and
unlabeled data. In 2011 IEEE 11th International
Conference on Data Mining Workshops, pages
244–250, Dec 2011.
"""
struct MinimumMarginPQs <: DataBasedPQs
    p_x::Array{Float64,1}
    bw_method::Union{String, U} where U <: Real
    p_inlier::Float64
    function MinimumMarginPQs(x::Array{T, 2}; bw_method="scott", p_inlier=nothing) where T <: Real
        ((p_inlier === nothing) || !(0 <= p_inlier <= 1)) && throw(ArgumentError("Invalid inlier probability $(p_inlier)."))
        return new(multi_kde(x, bw_method)(x), bw_method, p_inlier)
    end
end

function qs_score(qs::MinimumMarginPQs, x::Array{T, 2}, labels::Dict{Symbol, Array{Int, 1}})::Array{Float64, 1} where T <: Real
    haskey(labels, :Lin) || throw(MissingLabelTypeException(:Lin))
    p_x_inlier = multi_kde(x[:, labels[:Lin]], qs.bw_method)(x)
    return -abs.((2 * p_x_inlier * qs.p_inlier - qs.p_x) ./ qs.p_x)
end


text:
function solve()
    n = parse(Int, readline())
    if n%1000 == 0
        0
    else
        1000 - n % 1000
    end
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
using Test, NaiveNASlib, NaiveNASlib.Advanced, NaiveNASlib.Extend

include("testutil.jl")

# import Documenter makes include("tests/runtests.jl") workflows cumbersome as NaiveNASlib does no depend on it
# here is my workaround:
# ]activate test
# dev NaiveNASlib
# include("tests/runtests.jl")
# Delete test/Manifest.toml afterwards or else ]test from main project might not work

@testset "NaiveNASlib.jl" begin

	@info "Testing computation"
	include("vertex.jl")
	include("compgraph.jl")

	include("prettyprint.jl")

	@info "Testing mutation"

	include("mutation/vertex.jl")
	include("mutation/graph.jl")

	@info "Testing size mutation"
	include("mutation/size.jl")

	@info "Testing index mutation"
	include("mutation/select.jl")

	@info "Testing structural mutation"
	include("mutation/structure.jl")

	@info "Testing api"
	include("api/vertex.jl")

	@info "Testing doc examples"
	include("examples.jl")

	if Int !== Int32
		# Don't test documentation unless 64-bit os since some example print numerical types
		import Documenter
		Documenter.doctest(NaiveNASlib)
	end
end