text:
#=
    check_blocks(ϕ, Xhat0, U, overapproximate, n, b, N, blocks, prop)

Property checking of a given number of two-dimensional blocks of an affine
system with nondeterministic inputs.

The variants have the following structure:

INPUT:

- `ϕ` -- sparse matrix of a discrete affine system
- `Xhat0` -- initial set as a cartesian product over 2d blocks
- `U` -- input set of undeterministic inputs
- `overapproximate_inputs` -- function for overapproximation of inputs
- `n` -- ambient dimension
- `N` -- number of sets computed
- `blocks` -- the block indices to be computed
- `partition` -- the partition into blocks
- `prop` -- property to be checked

OUTPUT:

The first time index where the property is violated, and 0 if the property is satisfied.
=#

# helper functions
@inline proj(bi::UnitRange{Int}, n::Int) =
         sparse(1:length(bi), bi, ones(length(bi)), length(bi), n)
@inline proj(bi::Int, n::Int) = sparse([1], [bi], ones(1), 1, n)
@inline row(ϕpowerk::AbstractMatrix, bi::UnitRange{Int}) = ϕpowerk[bi, :]
@inline row(ϕpowerk::AbstractMatrix, bi::Int) = ϕpowerk[[bi], :]
@inline row(ϕpowerk::SparseMatrixExp, bi::UnitRange{Int}) = get_rows(ϕpowerk, bi)
@inline row(ϕpowerk::SparseMatrixExp, bi::Int) = Matrix(get_row(ϕpowerk, bi))
@inline block(ϕpowerk_πbi::AbstractMatrix, bj::UnitRange{Int}) = ϕpowerk_πbi[:, bj]
@inline block(ϕpowerk_πbi::AbstractMatrix, bj::Int) = ϕpowerk_πbi[:, [bj]]

# sparse
function check_blocks(ϕ::SparseMatrixCSC{NUM, Int},
                       Xhat0::Vector{<:LazySet{NUM}},
                       U::Union{ConstantInput, Void},
                       overapproximate_inputs::Function,
                       n::Int,
                       N::Int,
                       blocks::AbstractVector{Int},
                       partition::AbstractVector{<:Union{AbstractVector{Int}, Int}},
                       eager_checking::Bool,
                       prop::Property
                       )::Int where {NUM}
    violation_index = 0
    if !check_property(CartesianProductArray(Xhat0[blocks]), prop)
        if eager_checking
            return 1
        end
        violation_index = 1
    elseif N == 1
        return violation_index
    end

    b = length(blocks)
    Xhatk = Vector{LazySet{NUM}}(b)
    ϕpowerk = copy(ϕ)

    if U != nothing
        Whatk = Vector{LazySet{NUM}}(b)
        inputs = next_set(U)
        @inbounds for i in 1:b
            bi = partition[blocks[i]]
            Whatk[i] = overapproximate_inputs(1, blocks[i], proj(bi, n) * inputs)
        end
    end

    k = 2
    p = Progress(N, 1, "Computing successors ")
    @inbounds while true
        update!(p, k)
        for i in 1:b
            bi = partition[blocks[i]]
            Xhatk_bi = ZeroSet(length(bi))
            for (j, bj) in enumerate(partition)
                block = ϕpowerk[bi, bj]
                if findfirst(block) != 0
                    Xhatk_bi = Xhatk_bi + block * Xhat0[j]
                end
            end
            Xhatk[i] = (U == nothing ? Xhatk_bi : Xhatk_bi + Whatk[i])
        end

        if !check_property(CartesianProductArray(Xhatk), prop)
            if eager_checking
                return k
            elseif violation_index == 0
                violation_index = k
            end
        end
        if k == N
            break
        end

        if U != nothing
            for i in 1:b
                bi = partition[blocks[i]]
                Whatk[i] = overapproximate_inputs(k, blocks[i],
                    Whatk[i] + row(ϕpowerk, bi) * inputs)
            end
        end

        ϕpowerk = ϕpowerk * ϕ
        k += 1
    end

    return violation_index
end

# dense
function check_blocks(ϕ::AbstractMatrix{NUM},
                       Xhat0::Vector{<:LazySet{NUM}},
                       U::Union{ConstantInput, Void},
                       overapproximate_inputs::Function,
                       n::Int,
                       N::Int,
                       blocks::AbstractVector{Int},
                       partition::AbstractVector{<:Union{AbstractVector{Int}, Int}},
                       eager_checking::Bool,
                       prop::Property
                       )::Int where {NUM}
    violation_index = 0
    if !check_property(CartesianProductArray(Xhat0[blocks]), prop)
        if eager_checking
            return 1
        end
        violation_index = 1
    elseif N == 1
        return violation_index
    end

    b = length(blocks)
    Xhatk = Vector{LazySet{NUM}}(b)
    ϕpowerk = copy(ϕ)
    ϕpowerk_cache = similar(ϕ)

    if U != nothing
        Whatk = Vector{LazySet{NUM}}(b)
        inputs = next_set(U)
        @inbounds for i in 1:b
            bi = partition[blocks[i]]
            Whatk[i] = overapproximate_inputs(1, blocks[i], proj(bi, n) * inputs)
        end
    end

    arr_length = (U == nothing) ? length(partition) : length(partition) + 1
    k = 2
    p = Progress(N, 1, "Computing successors ")
    @inbounds while true
        update!(p, k)
        for i in 1:b
            bi = partition[blocks[i]]
            arr = Vector{LazySet{NUM}}(arr_length)
            for (j, bj) in enumerate(partition)
                arr[j] = ϕpowerk[bi, bj] * Xhat0[j]
            end
            if U != nothing
                arr[arr_length] = Whatk[i]
            end
            Xhatk[i] = MinkowskiSumArray(arr)
        end

        if !check_property(CartesianProductArray(Xhatk), prop)
            if eager_checking
                return k
            elseif violation_index == 0
                violation_index = k
            end
        end
        if k == N
            break
        end

        if U != nothing
            for i in 1:b
                bi = partition[blocks[i]]
                Whatk[i] = overapproximate_inputs(k, blocks[i],
                    Whatk[i] + row(ϕpowerk, bi) * inputs)
            end
        end

        A_mul_B!(ϕpowerk_cache, ϕpowerk, ϕ)
        copy!(ϕpowerk, ϕpowerk_cache)
        k += 1
    end

    return violation_index
end

# lazy_expm sparse
function check_blocks(ϕ::SparseMatrixExp{NUM},
                       assume_sparse::Val{true},
                       Xhat0::Vector{<:LazySet{NUM}},
                       U::Union{ConstantInput, Void},
                       overapproximate_inputs::Function,
                       n::Int,
                       N::Int,
                       blocks::AbstractVector{Int},
                       partition::AbstractVector{<:Union{AbstractVector{Int}, Int}},
                       eager_checking::Bool,
                       prop::Property
                       )::Int where {NUM}
    violation_index = 0
    if !check_property(CartesianProductArray(Xhat0[blocks]), prop)
        if eager_checking
            return 1
        end
        violation_index = 1
    elseif N == 1
        return violation_index
    end

    b = length(blocks)
    Xhatk = Vector{LazySet{NUM}}(b)
    ϕpowerk = SparseMatrixExp(copy(ϕ.M))

    if U != nothing
        Whatk = Vector{LazySet{NUM}}(b)
        inputs = next_set(U)
        @inbounds for i in 1:b
            bi = partition[blocks[i]]
            Whatk[i] = overapproximate_inputs(1, blocks[i], proj(bi, n) * inputs)
        end
    end

    k = 2
    p = Progress(N, 1, "Computing successors ")
    @inbounds while true
        update!(p, k)
        for i in 1:b
            bi = partition[blocks[i]]
            ϕpowerk_πbi = row(ϕpowerk, bi)
            Xhatk_bi = ZeroSet(length(bi))
            for (j, bj) in enumerate(partition)
                πbi = block(ϕpowerk_πbi, bj)
                if findfirst(πbi) != 0
                    Xhatk_bi = Xhatk_bi + πbi * Xhat0[j]
                end
            end
            Xhatk[i] = (U == nothing ? Xhatk_bi : Xhatk_bi + Whatk[i])
            if U != nothing
                Whatk[i] = overapproximate_inputs(k, blocks[i],
                    Whatk[i] + ϕpowerk_πbi * inputs)
            end
        end

        if !check_property(CartesianProductArray(Xhatk), prop)
            if eager_checking
                return k
            elseif violation_index == 0
                violation_index = k
            end
        end
        if k == N
            break
        end

        ϕpowerk.M .= ϕpowerk.M + ϕ.M
        k += 1
    end

    return violation_index
end

# lazy_expm dense
function check_blocks(ϕ::SparseMatrixExp{NUM},
                       assume_sparse::Val{false},
                       Xhat0::Vector{<:LazySet{NUM}},
                       U::Union{ConstantInput, Void},
                       overapproximate_inputs::Function,
                       n::Int,
                       N::Int,
                       blocks::AbstractVector{Int},
                       partition::AbstractVector{<:Union{AbstractVector{Int}, Int}},
                       eager_checking::Bool,
                       prop::Property
                       )::Int where {NUM}
    violation_index = 0
    if !check_property(CartesianProductArray(Xhat0[blocks]), prop)
        if eager_checking
            return 1
        end
        violation_index = 1
    elseif N == 1
        return violation_index
    end

    b = length(blocks)
    Xhatk = Vector{LazySet{NUM}}(b)
    ϕpowerk = SparseMatrixExp(copy(ϕ.M))

    if U != nothing
        Whatk = Vector{LazySet{NUM}}(b)
        inputs = next_set(U)
        @inbounds for i in 1:b
            bi = partition[blocks[i]]
            Whatk[i] = overapproximate_inputs(1, blocks[i], proj(bi, n) * inputs)
        end
    end

    arr_length = (U == nothing) ? length(partition) : length(partition) + 1
    k = 2
    p = Progress(N, 1, "Computing successors ")
    @inbounds while true
        update!(p, k)
        for i in 1:b
            bi = partition[blocks[i]]
            arr = Vector{LazySet{NUM}}(arr_length)
            ϕpowerk_πbi = row(ϕpowerk, bi)
            for (j, bj) in enumerate(partition)
                arr[j] = block(ϕpowerk_πbi, bj) * Xhat0[j]
            end
            if U != nothing
                arr[arr_length] = Whatk[i]
            end
            Xhatk[i] = MinkowskiSumArray(arr)
            if U != nothing
                Whatk[i] = overapproximate_inputs(k, blocks[i],
                    Whatk[i] + ϕpowerk_πbi * inputs)
            end
        end

        if !check_property(CartesianProductArray(Xhatk), prop)
            if eager_checking
                return k
            elseif violation_index == 0
                violation_index = k
            end
        end
        if k == N
            break
        end

        ϕpowerk.M .= ϕpowerk.M + ϕ.M
        k += 1
    end

    return violation_index
end


text:
#demo_psf_r.jl: psf read tests
#-------------------------------------------------------------------------------

using CMDimCircuits
CMDimCircuits.@using_CData()

try
import LibPSFC #To access sample data ONLY!!!
catch e
	msg = "Demo requires sampled data stored in the \"LibPSFC\" package."
	msg *= "\nTo execute demo, install with the following:"
	msg *= "\n\n    using Pkg; Pkg.add(PackageSpec(path=\"git://github.com/ma-laforge/LibPSFC.jl\"))"

	@error msg
	rethrow(e)
end

#Get a demo display:
include(CMDimCircuits.demoplotcfgscript)
#Normally use something like:
#CMDimData.@includepkg EasyPlotInspect

#==Constants
===============================================================================#
vvst = cons(:a, labels = set(yaxis="Amplitude (V)", xaxis="Time (s)"))
color1 = cons(:a, line = set(color=:red))
color2 = cons(:a, line = set(color=:blue))
color3 = cons(:a, line = set(color=:green))

#==Input data
===============================================================================#
testfile = "timeSweep"
filepath = joinpath(LibPSFC.rootpath, "core/data", testfile)

@info("Reading $filepath...")
EDAData.open_psf(filepath) do data
	global inp, outp0, outp1
	@show data
#	@show names(data) #Too many
	inp = read(data, "INP")
	outp0 = read(data, "OUTP<0>")
	outp1 = read(data, "OUTP<1>")
end

#==Generate plot
===============================================================================#
plot1 = push!(cons(:plot, vvst, title="VCO"),
	cons(:wfrm, inp, color1, label="Vinp"),
	cons(:wfrm, outp0, color2, label="Voutp<0>"),
	cons(:wfrm, outp1, color3, label="Voutp<1>"),
)

pcoll = push!(cons(:plot_collection, title="EDAData Tests: psf Format"), plot1)
	pcoll.ncolumns = 1

#==Display results as a plot
===============================================================================#
EasyPlot.displaygui(pcoll)

#==Return pcoll to user (call evalfile(...))
===============================================================================#
pcoll

ext: jl


text:
#=
# 1D Poisson, Dirichlet bc
# CG, Linear element
# Simplest test possible
=#
if !@isdefined(Femshop)
    include("../Femshop.jl");
    using .Femshop
end
init_femshop("poisson2d");

# Try making an optional log
@useLog("poisson2dlog")

# Set up the configuration (order doesn't matter)
@domain(2, SQUARE, UNIFORM_GRID)    # dimension, geometry, decomposition
@solver(CG)                         # DG, CG, etc.
@functionSpace(LEGENDRE, 1)         # function, order (or use testFunction and trialFunction)
@nodes(LOBATTO)                     # elemental node arrangement

# Specify the problem
@mesh(QUADMESH, [20,10], 1, [0,2,0,1])                   # .msh file or generate our own

@variable(u)                        # same as @variable(u, SCALAR)

@testSymbol(v)                    # sets the symbol for a test function

@boundary(u, 1, DIRICHLET, 0)

# Write the weak form
@coefficient(f, "-2*pi*pi*sin(pi*x)*sin(pi*y)")
@weakForm(u, "-dot(grad(u), grad(v)) - f*v")

solve(u);

# exact solution is sin(pi*x)*sin(pi*y)
# check error
maxerr = 0;
exact(x,y) = sin(pi*x)*sin(pi*y);

for i=1:size(Femshop.grid_data.allnodes,2)
    x = Femshop.grid_data.allnodes[1,i];
    y = Femshop.grid_data.allnodes[2,i];
    err = abs(u.values[i] - exact(x,y));
    global maxerr;
    maxerr = max(err,maxerr);
end
println("max error = "*string(maxerr));

# solution is stored in the variable's "values"
# using Plots
# pyplot();
# display(plot(Femshop.grid_data.allnodes[1,:], Femshop.grid_data.allnodes[2,:], u.values[:], st=:surface))

# check
log_dump_config();
log_dump_prob();

@finalize()


text:
# by lazarusA # HIDE
#using GLMakie, Random # HIDE
using CairoMakie, RDatasets, Colors, ColorSchemes
CairoMakie.activate!() #HIDE
let
    cars = dataset("datasets", "mtcars")
    byCat = cars.Cyl
    categ = unique(byCat)
    fig = Figure(resolution = (700, 850), font = "sans")
    ax1 = Axis(fig, xlabel = "MPG", ylabel = "density", xgridstyle=:dash,
            ygridstyle=:dash, rightspinevisible = false, topspinevisible = false)
    ax2 = Axis(fig, xlabel = "MPG", ylabel = "density")
    colors1 = ColorScheme(range(colorant"#65ADC2", colorant"#E84646",
                    length=length(categ)))
    colors2 = ColorScheme(range(colorant"orange", colorant"dodgerblue",
                    length=length(categ)))
    legends1 = []
    legends2 = []
    for (i,c) in enumerate(categ)
        indc = findall(x->x == c, byCat)
        line1 = density!(ax1, cars.MPG[indc], color = (colors1[i],0.5),
                    strokewidth=1.25, strokecolor = colors1[i])
        line2 = density!(ax2, cars.MPG[indc], color = (colors2[i],0.5),
                    strokewidth=1.25, strokecolor = colors2[i])
        push!(legends1, line1)
        push!(legends2, line2)
    end
    leg1 = Legend(fig, legends1, string.(categ), "Cyl", orientation = :horizontal,
    tellheight = true,tellwidth = false,framevisible = false,titleposition = :left)
    leg2 = Legend(fig, legends2, string.(categ), "Cyl")

    fig[1,1] = leg1
    fig[2,1] = ax1
    fig[3,1] = ax2
    fig[3,2] = leg2
    fig
    #save("FigAirqualityBox.png", fig, px_per_unit = 2) # HIDE
    save(joinpath(@__DIR__, "output", "DensitiesDataSet.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie", "RDatasets", "Colors", "ColorSchemes"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using MathOptInterface
const MOI = MathOptInterface

@testset "NaiveAllDiff (two binary vars)" begin
    optimizer = SCIP.Optimizer(display_verblevel=0)
    atol, rtol = 1e-6, 1e-6

    # add two binary variables: x, y
    x, y = MOI.add_variables(optimizer, 2)
    MOI.add_constraint(optimizer, MOI.SingleVariable(x), MOI.ZeroOne())
    MOI.add_constraint(optimizer, MOI.SingleVariable(y), MOI.ZeroOne())

    # maximize 2x + y
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0], [x, y]), 0.0))
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    # add constraint handler with constraint all-diff(x, y)
    alldiffch = NaiveAllDiff.NADCH(optimizer)
    SCIP.include_conshdlr(optimizer, alldiffch; needs_constraints=true)

    alldiffcons = NaiveAllDiff.NADCons([x, y])
    cr = SCIP.add_constraint(optimizer, alldiffch, alldiffcons)

    # solve problem and query result
    MOI.optimize!(optimizer)
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 2.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), x) ≈ 1.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), y) ≈ 0.0 atol=atol rtol=rtol
end

@testset "NaiveAllDiff (3 bin.vars, 2 pairwise conss)" begin
    optimizer = SCIP.Optimizer(display_verblevel=0)
    atol, rtol = 1e-6, 1e-6

    # add three binary variables
    x, y, z = MOI.add_variables(optimizer, 3)
    MOI.add_constraint(optimizer, MOI.SingleVariable(x), MOI.ZeroOne())
    MOI.add_constraint(optimizer, MOI.SingleVariable(y), MOI.ZeroOne())
    MOI.add_constraint(optimizer, MOI.SingleVariable(z), MOI.ZeroOne())

    # maximize 2x + 3y + 2z
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 3.0, 2.0],
                                                           [x, y, z]), 0.0))
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    # add constraint handler with constraints all-diff(x, y), all-diff(y, x)
    alldiffch = NaiveAllDiff.NADCH(optimizer)
    SCIP.include_conshdlr(optimizer, alldiffch; needs_constraints=true)

    SCIP.add_constraint(optimizer, alldiffch, NaiveAllDiff.NADCons([x, y]))
    SCIP.add_constraint(optimizer, alldiffch, NaiveAllDiff.NADCons([y, z]))

    # solve problem and query result
    MOI.optimize!(optimizer)
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 4.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), x) ≈ 1.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), y) ≈ 0.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), z) ≈ 1.0 atol=atol rtol=rtol
end

@testset "NaiveAllDiff (3 bin.vars, 3 pairwise conss)" begin
    optimizer = SCIP.Optimizer(display_verblevel=0)
    atol, rtol = 1e-6, 1e-6

    # add three binary variables
    x, y, z = MOI.add_variables(optimizer, 3)
    MOI.add_constraint(optimizer, MOI.SingleVariable(x), MOI.ZeroOne())
    MOI.add_constraint(optimizer, MOI.SingleVariable(y), MOI.ZeroOne())
    MOI.add_constraint(optimizer, MOI.SingleVariable(z), MOI.ZeroOne())

    # maximize 2x + 3y + 2z
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 3.0, 2.0],
                                                           [x, y, z]), 0.0))
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    # add constraint handler with constraints all-diff(x, y, z)
    alldiffch = NaiveAllDiff.NADCH(optimizer)
    SCIP.include_conshdlr(optimizer, alldiffch; needs_constraints=true)

    SCIP.add_constraint(optimizer, alldiffch, NaiveAllDiff.NADCons([x, y]))
    SCIP.add_constraint(optimizer, alldiffch, NaiveAllDiff.NADCons([x, z]))
    SCIP.add_constraint(optimizer, alldiffch, NaiveAllDiff.NADCons([y, z]))

    # solve problem and query result
    MOI.optimize!(optimizer)
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.INFEASIBLE
    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.NO_SOLUTION

    @test MOI.get(optimizer, MOI.NodeCount()) >= 8  # complete enumeration
end

@testset "NaiveAllDiff (3 bin.vars, 1 cons with all)" begin
    optimizer = SCIP.Optimizer(display_verblevel=0)
    atol, rtol = 1e-6, 1e-6

    # add three binary variables
    x, y, z = MOI.add_variables(optimizer, 3)
    MOI.add_constraint(optimizer, MOI.SingleVariable(x), MOI.ZeroOne())
    MOI.add_constraint(optimizer, MOI.SingleVariable(y), MOI.ZeroOne())
    MOI.add_constraint(optimizer, MOI.SingleVariable(z), MOI.ZeroOne())

    # maximize 2x + 3y + 2z
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 3.0, 2.0],
                                                           [x, y, z]), 0.0))
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    # add constraint handler with constraints all-diff(x, y, z)
    alldiffch = NaiveAllDiff.NADCH(optimizer)
    SCIP.include_conshdlr(optimizer, alldiffch; needs_constraints=true)

    SCIP.add_constraint(optimizer, alldiffch, NaiveAllDiff.NADCons([x, y, z]))

    # solve problem and query result
    MOI.optimize!(optimizer)
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.INFEASIBLE
    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.NO_SOLUTION

    @test MOI.get(optimizer, MOI.NodeCount()) >= 8  # complete enumeration
end

@testset "NaiveAllDiff (3 int.vars, 1 cons with all)" begin
    optimizer = SCIP.Optimizer(display_verblevel=0)
    atol, rtol = 1e-6, 1e-6

    # add three integer variables, in {0, 1, 2}
    x, y, z = MOI.add_variables(optimizer, 3)
    for v in [x, y, z]
        MOI.add_constraint(optimizer, MOI.SingleVariable(v), MOI.Integer())
        MOI.add_constraint(optimizer, MOI.SingleVariable(v), MOI.Interval(0.0, 2.0))
    end

    # maximize 2x + y
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0],
                                                           [x, y]), 0.0))
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    # add constraint handler with constraints all-diff(x, y, z)
    alldiffch = NaiveAllDiff.NADCH(optimizer)
    SCIP.include_conshdlr(optimizer, alldiffch; needs_constraints=true)

    SCIP.add_constraint(optimizer, alldiffch, NaiveAllDiff.NADCons([x, y, z]))

    # solve problem and query result
    MOI.optimize!(optimizer)
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

    @test MOI.get(optimizer, MOI.ObjectiveValue()) ≈ 5.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), x) ≈ 2.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), y) ≈ 1.0 atol=atol rtol=rtol
    @test MOI.get(optimizer, MOI.VariablePrimal(), z) ≈ 0.0 atol=atol rtol=rtol
end

@testset "NoGoodCounter (2 binary vars)" begin
    optimizer = SCIP.Optimizer(display_verblevel=0,
                               misc_allowdualreds=SCIP.FALSE)
    atol, rtol = 1e-6, 1e-6

    # add binary variables
    x, y = MOI.add_variables(optimizer, 2)
    MOI.add_constraint(optimizer, MOI.SingleVariable(x), MOI.ZeroOne())
    MOI.add_constraint(optimizer, MOI.SingleVariable(y), MOI.ZeroOne())

    # maximize 2x + y
    MOI.set(optimizer, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.([2.0, 1.0],
                                                           [x, y]), 0.0))
    MOI.set(optimizer, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    # add constraint handler with constraints
    counter = NoGoodCounter.Counter(optimizer, [x, y])
    SCIP.include_conshdlr(optimizer, counter; needs_constraints=false)

    # solve problem and query result
    MOI.optimize!(optimizer)
    @test MOI.get(optimizer, MOI.TerminationStatus()) == MOI.INFEASIBLE
    @test MOI.get(optimizer, MOI.PrimalStatus()) == MOI.NO_SOLUTION

    @test length(counter.solutions) == 4
end


text:
solution(conv::ConvergenceInfo) = _solution(conv.solver, conv.info, conv.options)
_solution(solver, info, options) =  info.solution


text:
export BrickletHallEffectEdgeCountConfig
struct BrickletHallEffectEdgeCountConfig
    edge_type::Integer
    debounce::Integer
end

export BrickletHallEffectEdgeInterrupt
struct BrickletHallEffectEdgeInterrupt
    count::Integer
    value::Bool
end

export BrickletHallEffectIdentity
struct BrickletHallEffectIdentity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletHallEffect
"""
Detects presence of magnetic field
"""
mutable struct BrickletHallEffect <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletHallEffect(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_hall_effect")
        deviceInternal = package.BrickletHallEffect(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_value
"""
    $(SIGNATURES)

Returns *true* if a magnetic field of 3.5 millitesla or greater is detected.
"""
function get_value(device::BrickletHallEffect)
    return device.deviceInternal.get_value()
end

export get_edge_count
"""
    $(SIGNATURES)

Returns the current value of the edge counter. You can configure
edge type (rising, falling, both) that is counted with
:func:`Set Edge Count Config`.

If you set the reset counter to *true*, the count is set back to 0
directly after it is read.
"""
function get_edge_count(device::BrickletHallEffect, reset_counter)
    return device.deviceInternal.get_edge_count(reset_counter)
end

export set_edge_count_config
"""
    $(SIGNATURES)

The edge type parameter configures if rising edges, falling edges or
both are counted. Possible edge types are:

* 0 = rising
* 1 = falling
* 2 = both

A magnetic field of 3.5 millitesla or greater causes a falling edge and a
magnetic field of 2.5 millitesla or smaller causes a rising edge.

If a magnet comes near the Bricklet the signal goes low (falling edge), if
a magnet is removed from the vicinity the signal goes high (rising edge).

Configuring an edge counter resets its value to 0.

If you don't know what any of this means, just leave it at default. The
default configuration is very likely OK for you.
"""
function set_edge_count_config(device::BrickletHallEffect, edge_type, debounce)
    device.deviceInternal.set_edge_count_config(edge_type, debounce)
end

export get_edge_count_config
"""
    $(SIGNATURES)

Returns the edge type and debounce time as set by :func:`Set Edge Count Config`.
"""
function get_edge_count_config(device::BrickletHallEffect)
    return device.deviceInternal.get_edge_count_config()
end

export set_edge_interrupt
"""
    $(SIGNATURES)

Sets the number of edges until an interrupt is invoked.

If *edges* is set to n, an interrupt is invoked for every n-th detected edge.

If *edges* is set to 0, the interrupt is disabled.
"""
function set_edge_interrupt(device::BrickletHallEffect, edges)
    device.deviceInternal.set_edge_interrupt(edges)
end

export get_edge_interrupt
"""
    $(SIGNATURES)

Returns the edges as set by :func:`Set Edge Interrupt`.
"""
function get_edge_interrupt(device::BrickletHallEffect)
    return device.deviceInternal.get_edge_interrupt()
end

export set_edge_count_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Edge Count` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Edge Count` callback is only triggered if the edge count has changed
since the last triggering.
"""
function set_edge_count_callback_period(device::BrickletHallEffect, period)
    device.deviceInternal.set_edge_count_callback_period(period)
end

export get_edge_count_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Edge Count Callback Period`.
"""
function get_edge_count_callback_period(device::BrickletHallEffect)
    return device.deviceInternal.get_edge_count_callback_period()
end

export edge_interrupt
"""
    $(SIGNATURES)

This callback is triggered every n-th count, as configured with
:func:`Set Edge Interrupt`. The parameters are the
current count and the current value (see :func:`Get Value` and
:func:`Get Edge Count`).
"""
function edge_interrupt(device::BrickletHallEffect)
    return device.deviceInternal.edge_interrupt()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletHallEffect)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletHallEffect, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
using DataStructures

"""
This implements an incorrect algorithm which is the most common
intuitive choice. Each trajectory's distribution looks correct,
but the master equation isn't correct.
"""

struct NaiveSampler
    firing_queue::MutableBinaryHeap{NRTransition,Base.Order.ForwardOrdering}
    # This maps from transition to entry in the firing queue.
    transition_entry::Dict{Any,Int}
    disabled::Set{Any}
    init::Bool
end

"""
Construct a NaiveSampler.
This doesn't require a clock that remembers integrated hazard.
This sampler is inappropriate if any transition is either
re-enabled after being disabled or has its distribution
modified while enabled. It's OK if each transition fires
only once.

BUT I can't create a model where this sampler's output
varies from First Reaction. If anyone can show me when
this fails to work, I'd be grateful. Even the sis.jl example works.
"""
function NaiveSampler()
    heap=mutable_binary_minheap(NRTransition)
    state=Dict{Any,Int}()
    NaiveSampler(heap, state, Set{Any}(), true)
end

# Finds the next one without removing it from the queue.
function Next(propagator::NaiveSampler, system, rng)
    if propagator.init
        Hazards(system, rng) do clock, now, updated, rng2
            NaiveObserve(propagator, clock, now, updated, rng2)
        end
        propagator.init=false
    end

    NotFound=NRTransition(nothing, Inf)
    if !isempty(propagator.firing_queue)
        least=top(propagator.firing_queue)
    else
        least=NotFound
    end
    @debug("SampleSemiMarkov.next queue length ",
            length(propagator.firing_queue), " least ", least)
    (least.time, least.key)
end

"""
Returns an observer of intensities to decide what to
do when they change.
"""
function Observer(propagator::NaiveSampler)
    function nrobserve(clock, time, updated, rng)
        NaiveObserve(propagator, clock, time, updated, rng)
    end
end

function NaiveObserve(propagator::NaiveSampler, clock,
            time, updated, rng)
    key=clock
    if updated==:Fired || updated==:Disabled
        heap_handle=propagator.transition_entry[key]
        # We store distributions in order to calculate remaining hazard
        # which will happen AFTER the state has changed.
        update!(propagator.firing_queue, heap_handle,
            NRTransition(key, -1.))
        todelete=pop!(propagator.firing_queue)
        delete!(propagator.transition_entry, key)
        push!(propagator.disabled, clock)

    elseif updated==:Enabled
        # if haskey(propagator.disabled, clock)
        #     error("Cannot re-enable a transition with this sampler.")
        # end
        when_fire=Sample(clock.intensity, time, rng)
        heap_handle=push!(propagator.firing_queue,
                NRTransition(key, when_fire))
        propagator.transition_entry[key]=heap_handle

    elseif updated==:Modified
        # if haskey(propagator.disabled, clock)
        #     error("Cannot modify a transition with this sampler.")
        # end
        when_fire=Sample(clock.intensity, time, rng)
        heap_handle=propagator.transition_entry[key]
        update!(propagator.firing_queue, heap_handle,
                NRTransition(key, when_fire))
    else
        assert(false)
    end
end


text:
import JSON

abstract type AbstractDataFrame end

abstract type StatVariate end

abstract type AbstractFactor{T<:Unsigned} <: StatVariate end

abstract type AbstractCovariate{T<:AbstractFloat} <: StatVariate end

abstract type AbstractBoolVariate <: StatVariate end

struct DataFrame <: AbstractDataFrame
    length::Int64
    factors::AbstractVector{<:AbstractFactor}
    covariates::AbstractVector{<:AbstractCovariate}
    boolvariates::AbstractVector{<:AbstractBoolVariate}
end

mutable struct CovariateStats
    obscount::Int64
    nancount::Int64
    nanpcnt::Float64
    sum::Float64
    sum2::Float64
    mean::Float64
    std::Float64
    min::Float64
    max::Float64
end

mutable struct LevelStats
    level::String
    freq::Int64
    freqpcnt::Float64
end

mutable struct FactorStats
    obscount::Int64
    missingfreq::Int64
    missingpcnt::Float64
    levelstats::Vector{LevelStats}
end

function Base.length(dataframe::AbstractDataFrame)
    dataframe.length
end

getlevels(factor::AbstractFactor{T}) where {T<:Unsigned} = factor.levels

getname(factor::AbstractFactor{T}) where {T<:Unsigned} = factor.name

getname(covariate::AbstractCovariate{T}) where {T<:AbstractFloat} = covariate.name

getname(boolvar::AbstractBoolVariate) = boolvar.name

function widenfactors(factors::Vector{<:AbstractFactor})
    if all(map((factor -> issubtype(typeof(factor), AbstractFactor{UInt8})), factors))
        factors
    elseif all(map((factor -> issubtype(typeof(factor), AbstractFactor{UInt8}) || issubtype(typeof(factor), AbstractFactor{UInt16})), factors))
        [issubtype(typeof(factor), AbstractFactor{UInt16}) ? factor : WiderFactor{UInt8, UInt16}(factor) for factor in factors]
    else
        [issubtype(typeof(factor), AbstractFactor{UInt32}) ? factor : (issubtype(typeof(factor), AbstractFactor{UInt16}) ? WiderFactor{UInt16, UInt32}(factor) : WiderFactor{UInt8, UInt32}(factor)) for factor in factors]
    end
end

function DataFrame(path::String; preload::Bool = true)
    path = abspath(path)
    headerpath = isfile(path) ? path : joinpath(path, "header.txt")
    headerjson = open(headerpath) do f
        readstring(f)
    end
    header = JSON.parse(headerjson)
    factors = Vector{AbstractFactor}()
    covariates = Vector{AbstractCovariate}()
    datacols = header["datacolumns"]
    len = 0
    for datacol in datacols
        datatype = datacol["datatype"]
        len = datacol["length"]
        name = datacol["name"]
        datpath = joinpath(dirname(headerpath), datacol["filename"])

        if datatype == "Float32"
            if preload
                push!(covariates, Covariate{Float32}(name, len, datpath))
            else
                push!(covariates, FileCovariate{Float32}(name, len, datpath))
            end
        end

        if datatype == "UInt8"
            levels = [string(level) for level in datacol["levels"]]
            if length(levels) == 0
                levels = Vector{String}()
            end
            if preload
                push!(factors, Factor{UInt8}(name, len, levels, datpath))
            else
                push!(factors, FileFactor{UInt8}(name, len, levels, datpath))
            end
        end

        if datatype == "UInt16"
            levels = [string(level) for level in datacol["levels"]]
            if length(levels) == 0
                levels = Vector{String}()
            end
            if preload
                push!(factors, Factor{UInt16}(name, len, levels, datpath))
            else
                push!(factors, FileFactor{UInt16}(name, len, levels, datpath))
            end
        end

        if datatype == "UInt32"
            levels = [string(level) for level in datacol["levels"]]
            if length(levels) == 0
                levels = Vector{String}()
            end
            if preload
                push!(factors, Factor{UInt32}(name, len, levels, datpath))
            else
                push!(factors, FileFactor{UInt32}(name, len, levels, datpath))
            end
        end
    end
    DataFrame(len, factors, covariates, AbstractBoolVariate[])
end

function Base.getindex(df::AbstractDataFrame, name::String)
    for factor in df.factors
        if getname(factor) == name
            return factor
        end
    end
    for cov in df.covariates
        if getname(cov) == name
            return cov
        end
    end
end

function Base.summary(factor::AbstractFactor{T}) where {T<:Unsigned}
    io = IOBuffer()
    factorstats = getstats(factor)
    println(io, @sprintf("%-16s%15d", "Obs Count", factorstats.obscount))
    println(io, @sprintf("%-16s%15s%15s", "Level", "Frequency", "Frequency(%)"))
    println(io, @sprintf("%-16s%15d%15G", MISSINGLEVEL, factorstats.missingfreq, factorstats.missingpcnt))
    for levelstats in factorstats.levelstats
        println(io, @sprintf("%-16s%15d%15G", levelstats.level, levelstats.freq, levelstats.freqpcnt))
    end
    print(String(take!(io)))
end

function Base.summary(boolvar::AbstractBoolVariate)
    io = IOBuffer()
    len = length(boolvar)
    slices = slice(factor, 1, len, SLICELENGTH)
    truefreq = fold(0, slices) do acc, slice
        res = acc
        for v in slice
            if v
                res += 1
            end
        end
        res
    end
    println(io, @sprintf("%-15s%12d", "Obs Count", len))
    println(io, @sprintf("%-15s%12d", "True Freq", truefreq))
    println(io, @sprintf("%-15s%12G", "True Freq (%)", 100.0 * truefreq / len))
    print(String(take!(io)))
end

function getstats(factor::AbstractFactor{T}) where {T<:Unsigned}
    len = length(factor)
    levels = getlevels(factor)
    levelcount = length(levels)
    init = zeros(Int64, levelcount + 1)
    slices = slice(factor, 1, len, SLICELENGTH)
    freq = fold(init, slices) do frq, slice
        for levelindex in slice
            frq[levelindex + 1] += 1
        end
        frq
    end
    missingfreq = freq[1]
    missingpcnt = 100.0 * missingfreq / len
    levelstats = [LevelStats(levels[i], freq[i + 1], 100.0 * freq[i + 1] / len) for i in 1:levelcount]
    FactorStats(len, missingfreq, missingpcnt, levelstats)
end

function getstats(covariate::AbstractCovariate{T}) where {T<:AbstractFloat}
    len = length(covariate)
    init = CovariateStats(0, 0, NaN64, NaN64, NaN64, NaN64, NaN64, NaN64, NaN64)
    slices = slice(covariate, 1, len, SLICELENGTH)
    stats = fold(init, slices) do s, slice
        for v in slice
            if isnan(v)
                s.nancount += 1
            else
                if isnan(s.sum)
                    s.sum = v
                    s.sum2 = v * v
                    s.min = v
                    s.max = v
                else
                    s.sum += v
                    s.sum2 += v * v
                    if v < s.min
                        s.min = v
                    end
                    if v > s.max
                        s.max = v
                    end
                end
            end
        end
        s
    end
    stats.obscount = len
    stats.nanpcnt = 100.0 * stats.nancount / len
    stats.mean = stats.sum / (len - stats.nancount)
    stats.std = sqrt(((stats.sum2 - stats.sum * stats.sum / (len - stats.nancount)) / (len - stats.nancount - 1)))
    stats
end

function Base.summary(covariate::AbstractCovariate{T}) where {T<:AbstractFloat}
    io = IOBuffer()
    stats = getstats(covariate)
    println(io, @sprintf("%-15s%12d", "Obs Count", stats.obscount))
    println(io, @sprintf("%-15s%12d", "NaN Freq", stats.nancount))
    println(io, @sprintf("%-15s%12G", "NaN %", stats.nanpcnt))
    println(io, @sprintf("%-15s%12G", "Min", stats.min))
    println(io, @sprintf("%-15s%12G", "Max", stats.max))
    println(io, @sprintf("%-15s%12G", "Mean", stats.mean))
    println(io, @sprintf("%-15s%12G", "Std", stats.std))
    print(String(take!(io)))
end

function Base.show(io::IO, covariate::AbstractCovariate{T}) where {T<:AbstractFloat}
    slices = slice(covariate, 1, HEADLENGTH, HEADLENGTH)
    slice1, _ = tryread(slices)
    len = length(covariate)
    if !isnull(slice1)
        datahead = join([isnan(v) ? "." : string(v) for v in get(slice1)], " ")
        dataend = len > HEADLENGTH ? "  ..." : ""
        println(io, "Covariate $(getname(covariate)) with $(len) obs: $(datahead)$dataend")
    else
        println(io, "Covariate $(getname(covariate)) with $(len) obs")
    end
end

function Base.show(io::IO, boolvar::AbstractBoolVariate)
    slices = slice(boolvar, 1, HEADLENGTH, HEADLENGTH)
    slice1, _ = tryread(slices)
    len = length(boolvar)
    if !isnull(slice1)
        datahead = join([string(v) for v in get(slice1)], " ")
        dataend = len > HEADLENGTH ? "  ..." : ""
        println(io, "BoolVar $(getname(boolvar)) with $(len) obs: $(datahead)$dataend")
    else
        println(io, "BoolVar $(getname(boolvar)) with $(len) obs")
    end
end

function Base.show(io::IO, factor::AbstractFactor{T}) where {T<:Unsigned}
    slices = slice(factor, 1, HEADLENGTH, HEADLENGTH)
    slice1, _ = tryread(slices)
    len = length(factor)
    levels = getlevels(factor)
    levelcount = length(levels)
    if !isnull(slice1)
        datahead = join([index == 0 ? MISSINGLEVEL : levels[index] for index in get(slice1)], " ")
        dataend = len > HEADLENGTH ? "  ..." : ""
        println(io, "Factor $(getname(factor)) with $(len) obs and $(levelcount) levels: $(datahead)$dataend")
    else
        println(io, "Factor $(getname(factor)) with $(len) obs and $(levelcount) levels")
    end
end

function Base.convert(::Type{Vector{T}}, covariate::AbstractCovariate{T}) where {T<:AbstractFloat}
    v, _ = tryread(slice(covariate, 1, length(covariate), length(covariate)))
    get(v)
end

function isordinal(factor::AbstractFactor{T}) where {T<:Unsigned}
    false
 end

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
using Test
using Quon: QuonConst, Zero, Pi, HalfPi, NegHalfPi, InfZero, InfPi
using Quon
using Yao

to_Rz(p::Complex) = Rz(p/im)
function to_Rz(p::QuonConst)
    p === Pi && return Rz(π)
    p === Zero && return Rz(0.0)
    p === HalfPi && return Rz(π/2)
    p === NegHalfPi && return Rz(-π/2)
    p === InfZero && return matblock(projector(0))
    p === InfPi && return matblock(projector(1))
end
to_Rx(p::Complex) = Rx(p/im)
function to_Rx(p::QuonConst)
    p === Pi && return Rx(π)
    p === Zero && return Rx(0.0)
    p === HalfPi && return Rx(π/2)
    p === NegHalfPi && return Rx(-π/2)
    p === InfZero && return matblock((1/sqrt(2) + 0.0im)*[1 1; 1 1])
    p === InfPi && return matblock((1/sqrt(2) + 0.0im)*[1 -1; -1 1])
end

function test_quon_const(a, b, c)
    circ1 = chain(1, to_Rz(a), to_Rx(b), to_Rz(c))
    p, q, r = yang_baxter_param(a, b, c)
    circ2 = chain(1, to_Rx(p), to_Rz(q), to_Rx(r))
    return operator_fidelity(matblock(mat(circ1)), matblock(mat(circ2)))/
        sqrt(operator_fidelity(matblock(mat(circ1)), matblock(mat(circ1))) * operator_fidelity(matblock(mat(circ2)), matblock(mat(circ2))))
end

AllConsts = (Zero, Pi, HalfPi, NegHalfPi, InfZero, InfPi)
FiniteConsts = (Zero, Pi, HalfPi, NegHalfPi)

test_quon_const(Pi, Pi, Pi)

test_123 = [test_quon_const(p1, p2, p3) for p1 in AllConsts, p2 in AllConsts, p3 in AllConsts]
test_12 = [test_quon_const(p1, p2, Complex(p3)) for p1 in AllConsts, p2 in AllConsts, p3 in FiniteConsts]
test_23 = [test_quon_const(Complex(p1), p2, p3) for p1 in FiniteConsts, p2 in AllConsts, p3 in AllConsts]
test_13 = [test_quon_const(p1, Complex(p2), p3) for p1 in AllConsts, p2 in FiniteConsts, p3 in AllConsts]
test_1 = [test_quon_const(p1, Complex(p2), Complex(p3)) for p1 in AllConsts, p2 in FiniteConsts, p3 in FiniteConsts]
test_2 = [test_quon_const(Complex(p1), p2, Complex(p3)) for p1 in FiniteConsts, p2 in AllConsts, p3 in FiniteConsts]
test_3 = [test_quon_const(Complex(p1), Complex(p2), p3) for p1 in FiniteConsts, p2 in FiniteConsts, p3 in AllConsts]

@test sum(test_123 .≈ 1) == length(test_123) - 4
@test sum(test_12 .≈ 1) == length(test_12)
@test sum(test_13 .≈ 1) == length(test_13) - 4
@test sum(test_23 .≈ 1) == length(test_23)
@test sum(test_1 .≈ 1) == length(test_1)
@test sum(test_2 .≈ 1) == length(test_2)
@test sum(test_3 .≈ 1) == length(test_3)


text:
function centroid(points)
    com = sum(points)
    return com/size(points, 1)
end

function covmat(points, com)
    @assert size(points,1) > 1
    covmat = zeros(3,3)
    for i in eachindex(points)
        pt = points[i]-com

        covmat[2,2] += pt[2]*pt[2]
        covmat[2,3] += pt[2]*pt[3]
        covmat[3,3] += pt[3]*pt[3]

        covmat[1,1] += pt[1]*pt[1]
        covmat[1,2] += pt[1]*pt[2]
        covmat[1,3] += pt[1]*pt[3]
    end
    covmat[2,1] = covmat[1,2]
    covmat[3,1] = covmat[1,3]
    covmat[3,2] = covmat[2,3]
    return covmat
end

function covmat(points)
    c = centroid(points)
    return covmat(points, c)
end

"""
    normedcovmat(points, centr)

Compute the normed covariance matrix with the given centroid.
"""
function normedcovmat(points, c)
    m = covmat(points, c)
    return m/size(points,1)
end

"""
    normedcovmat(points)

Compute the normed covariance matrix of the points.
"""
function normedcovmat(points)
    c = centroid(points)
    m = covmat(points, c)
    return m/size(points,1)
end

"""
    vec2homvec(ps)

Convert 3D vectors to homogeneous vectors: append 1 to each vector.
Return a vector of SVectors.
"""
function vec2homvec(ps)
    return [SVector{4,Float64}(p..., 1) for p in ps]
end

"""
    homvec2vec(ps)

Convert 3D homogeneous vectors to "normal" vectors: drop the 1 from the end.
Return a vector of SVectors.
"""
function homvec2vec(ps)
    return [SVector{3,Float64}(p[1:3]) for p in ps]
end

"""
    fit3pointimplplane(points)

Fit an `ImplicitPlane` to 3 points. First is used as one point of the plane.
"""
function fit3pointimplplane(points)
    v1 = normalize(points[2]-points[1])
    v2 = normalize(points[3]-points[1])
    n = normalize(cross(v1, v2))
    return ImplicitPlane(points[1], n)
end

"""
    findOBB(points)

Compute the oriented bounding box of points (in 3D).
Return the cornerpoints and the 6 planes (as `ImplicitPlane`s).
"""
function findOBB(points)
    c = centroid(points)
    m = normedcovmat(points, c)
    evs = eigvecs(m)
    zv = cross(evs[:,1], evs[:,2])
    if abs(dot(zv, evs[:,3])) < cosd(5)
        @warn "evs[3]: $(evs[3]) and zv: $zv are not parallel!"
    end
    evs[:,3] = zv
    # transformation matrix
    rmat = zeros(4,4)
    rmat[4,4] = 1
    rmat[1:3,1:3] = evs'
    rmat[1:3,4] = -1 .* evs' * c
    trm = SMatrix{4,4,Float64}(rmat)

    p_transformed = [trm*p for p in vec2homvec(points)]
    aabb_tr = homvec2vec(findAABB(p_transformed))
    meanDiagonal = sum(aabb_tr)/2

    # size of oriented bounding box
    bbs = aabb_tr[2]-aabb_tr[1]
    nv = -bbs/2
    # box around origo
    bi = 0:1
    corners_ = [nv+bbs .*[i,j,k] for i in bi for j in bi for k in bi]
    corners = vec2homvec(corners_)

    # final transformation
    finalTR = zeros(4,4)
    finalTR[4,4] = 1
    finalTR[1:3,1:3] = evs
    finalTR[1:3,4] = evs*meanDiagonal+c
    ftr = SMatrix{4,4,Float64}(finalTR)

    obbcorners_ = [ftr*p for p in corners]
    obbcorners = homvec2vec(obbcorners_)

    # points are selected so that the normals point outside from the cube
    planetris = [[1,5,2], [1,2,3], [1,3,5],
                [8,4,6], [8,6,7], [8,7,4]]
    planes = [fit3pointimplplane(obbcorners[i]) for i in planetris]
    return obbcorners, planes
end


text:
module GraphingUtilities
    using Statistics, Plots, ..Representation
    function plot_all(fitness_values::Array{Array{Float64, 1}, 1},
                        num_generations::Int, population::Representation.Population,
                        true_values::Array{Float64, 1}, x_values::Array{Float64, 1},
                        hit_ratios::Array{Array{Float64, 1}, 1},
                        population_size::Int64,
                        save_name::Union{Nothing, AbstractString})

        error_plot = plot(xlabel="Generation", ylabel="Error")
        plot_average_fitness!(fitness_values, num_generations)
        plot_max_fitness!(fitness_values, num_generations)
        func_plot = plot_both_functions(population, true_values, x_values)
        hit_ratio_plot = plot_metric(hit_ratios, num_generations, population_size,"Hit Ratio")
        if save_name != nothing
            savefig(error_plot, save_name * "_error.png")
            savefig(func_plot, save_name * "_funcs.png")
            savefig(hit_ratio_plot, save_name * "_hr.png")
        else
            display(error_plot)
            display(func_plot)
            display(hit_ratio_plot)
        end
    end

    function plot_metric(metric_values::Array{Array{Float64, 1}, 1},
                            num_generations::Int, population_size::Int64,
                            metric_name::AbstractString)
        averages = [mean(generation_values) for generation_values in metric_values]
        averages = averages / (population_size + 1)
        plot(1:num_generations, averages, lw=3,
                xlabel="Generation", ylabel=metric_name, label=metric_name)
    end

    function plot_average_fitness!(fitness_values::Array{Array{Float64, 1}, 1},
                                    num_generations::Int)
        averages = [mean(generation_fitvalues) for generation_fitvalues in fitness_values]
        plot!(1:num_generations, averages, lw=3,
                xlabel="Generation", ylabel="Error", label="Mean Fitness")
    end

    function plot_max_fitness!(fitness_values::Array{Array{Float64, 1}, 1},
                                num_generations::Int)
        minimums = [minimum(generation_fitvalues) for generation_fitvalues in fitness_values]
        plot!(1:num_generations, minimums, lw=3, label="Max fitness")
    end

    function plot_both_functions(population::Representation.Population,
                                 true_values::Array{Float64, 1},
                                 x_values::Array{Float64, 1})
        func_plot = plot(xlabel="X", ylabel="Y")
        fittest = sort(population.members, by=i->i.fitness)[end]
        pred_values = [kernels_to_values(x, fittest.gauss_kernels) for x in true_values]
        plot!(x_values, pred_values, lw=2, label="Pred function")
        plot!(x_values, true_values, lw=2, label="True function")
        return func_plot
    end

    function kernels_to_values(x::Float64,
                                kernels::Array{Tuple{Float64, Float64, Float64}, 1})
        sum = Float32(0.0)
        for kernel in kernels
            (weight, c, γ) = kernel
            sum += weight * exp(-γ*(c - x)^2)
        end
        return sum
    end

end


text:
using Utils
using Test

@testset "@something_nothing" begin
    @testset "doesn't throw" begin
        @test isnothing(@something_nothing nothing)
        @test isnothing(@something_nothing nothing nothing)
        @test isnothing(@something_nothing nothing nothing nothing)
    end

    @testset "unwrapping" begin
        @test (@something_nothing Some(0)) == 0
        @test (@something_nothing nothing Some(0)) == 0
        @test (@something_nothing nothing nothing Some(0)) == 0
        @test (@something_nothing nothing nothing Some(0) nothing) == 0
        @test isnothing(@something_nothing Some(nothing) 1)
        @test isnothing(@something_nothing nothing Some(nothing) 1)
        @test isnothing(@something_nothing nothing Some(nothing) 1 nothing)
    end

    @testset "wrapped nothing at the end" begin
        @test isnothing(@something_nothing Some(nothing))
        @test isnothing(@something_nothing nothing Some(nothing))
        @test isnothing(@something_nothing nothing nothing Some(nothing))
    end

    @testset "escaping" begin
        val = 2

        @test_throws DomainError (@something_nothing sqrt(1 - val))
        @test (@something_nothing nothing (1 - val)) == -1
    end
end

max_stars_repo_name: CmdQ/AoC


text:
using BinaryBuilder, Pkg

name = "Couenne"
version = v"0.5.8"

sources = [
    GitSource("https://github.com/coin-or/Couenne.git",
              "7154f7a9b3cd84be378d02b483d090b76fc79ce8"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/Couenne

atomic_patch -p1 ${WORKSPACE}/srcdir/patches/fixes.patch

if [[ ${target} == *mingw* ]]; then
    # Disable user-interrupts on Windows.
    sed -i s/'#define SIGNAL'// ${WORKSPACE}/srcdir/Couenne/Couenne/src/main/CouenneBab.cpp
    # Export symbols from DLL
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/export.patch
fi

# asl.h defines a macro called `filename`. Very unhelpful.
sed -i s/'#define filename'/'#define __asl_filename'/ /workspace/destdir/include/asl.h
sed -i s/'filename;'/'__asl_filename;'/ /workspace/srcdir/Couenne/Couenne/src/readnl/readnl.cpp

# Remove misleading libtool files
rm -f ${libdir}/*.la
update_configure_scripts

# old and custom autoconf
sed -i s/elf64ppc/elf64lppc/ configure

export CPPFLAGS="${CPPFLAGS} -I${includedir} -I${includedir}/coin -I${includedir}/coin-or"
export CXXFLAGS="${CXXFLAGS} -std=c++11"

if [[ ${target} == *mingw* ]]; then
    export LDFLAGS="-L${libdir}"
    export LT_LDFLAGS="-no-undefined"
elif [[ ${target} == *apple* ]]; then
    # Turn off the annoying, but harmless, -Wdeprecated-register.
    export CPPFLAGS="${CPPFLAGS} -Wno-deprecated-register"
fi

./configure \
    --prefix=$prefix \
    --build=${MACHTYPE} \
    --host=${target} \
    --enable-static \
    --enable-shared \
    lt_cv_deplibs_check_method=pass_all \
    --with-asl-lib="-lasl -lipoptamplinterface" \
    --with-bonmin-lib="-lbonminampl -lbonmin -lipoptamplinterface -lipopt -lCbc -lCgl -lOsiClp -lClp -lOsi -lCoinUtils -lasl -lopenblas"

make -j${nproc}
make install
"""

platforms = supported_platforms()
platforms = filter!(!Sys.isfreebsd, platforms)
platforms = expand_cxxstring_abis(platforms)
platforms = filter(x -> cxxstring_abi(x) != "cxx03", platforms)

# The products that we will ensure are always built
products = [
    LibraryProduct("libCouenne", :libcouenne),
    ExecutableProduct("couenne", :amplexe),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("ASL_jll", v"0.1.2"),
    Dependency("Bonmin_jll", v"1.8.8"),
    Dependency("Cbc_jll", v"2.10.5"),
    Dependency("Ipopt_jll", v"3.13.4"),
]

# Note: for obscure reasons I miss, `Clp_jll` built for
# `x86_64-linux-gnu-libgfortran3-cxx11` with GCC 6.1 provides the symbol
#     virtual thunk to OsiClpSolverInterface::getRowName(int, unsigned int) const
# instead of
#     virtual thunk to OsiClpSolverInterface::getRowName[abi:cxx11](int, unsigned int) const
# Thus, if we build Couenne with GCC < 7 (=> libgfortran3) we'd have the symbol
# without the `abi:cxx11` tag.  Auditor doesn't complain because at the end of
# the build everything has a consistent ABI for the libgfortran3 runtime.
# However, when we load the package at runtime we'd have libgfortran5 and the
# corresponding Clp provides the `abi:cxx11`-tagged symbol, causing the error
#    /tmp/jl_KKxJ1h/artifacts/ecc837f417130ef88b6288424541354733c387b1/lib/libCouenne.so: undefined symbol: _ZTv0_n720_NK21OsiClpSolverInterface10getRowNameEij
# The solution is to build Couenne with GCC 7 (=> libgfortran4) so that
# `Clp_jll` used during the build has the correctly tagged symbol.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"7")

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
@symbol_func function B_V(cur_reactor::AbstractReactor)
  cur_B = 0.1

  cur_B *= cur_reactor.I_P

  cur_B /= cur_reactor.R_0

  cur_B *= ( K_VI(cur_reactor) + beta_P(cur_reactor) )

  cur_B
end


text:
"""
    getobjectiveconstraints(model::AbstractModel)

Returns an instance of `AbstractFunction` that when called with input `x` returns the objective and constraint values of `model` at `x` in a single vector.
"""
getobjectiveconstraints

"""
    getobjective(model::AbstractModel)

Returns a that when called with input `x` returns the objective value of `model` at `x`.
"""
getobjective

"""
    getineqconstraint(model::AbstractModel, i::Integer)

Returns a that when called with input `x` returns the constraint violation value of the `i`^th inequality constraint of `model` at `x`.
"""
getineqconstraint

"""
    geteqconstraint(model::AbstractModel, i::Integer)

Returns a that when called with input `x` returns the constraint violation value of the `i`^th equality constraint of `model` at `x`.
"""
geteqconstraint

"""
    getineqconstraints(model::AbstractModel)

Returns a that when called with input `x` returns the constraint violation values of all the inequality constraints of `model` at `x`.
"""
getineqconstraints

"""
    geteqconstraints(model::AbstractModel)

Returns a that when called with input `x` returns the constraint violation values of all the equality constraints of `model` at `x`.
"""
geteqconstraints

"""
    get_objective_multiple(model::AbstractModel)

Returns the factor by which the objective of `model` was scaled. This is 1 by default.
"""
get_objective_multiple

"""
    set_objective_multiple!(model::AbstractModel, multiple::Number)

Set the objective's scaling factor of `model` to `multiple`.
"""
set_objective_multiple!

"""
    getnconstraints(model::AbstractModel)

Returns the total number of constraints in `model`.
"""
getnconstraints

"""
    getnineqconstraints(model::AbstractModel)

Returns the number of inequality constraints in `model`.
"""
getnineqconstraints

"""
    getneqconstraints(model::AbstractModel)

Returns the number of equality constraints in `model`.
"""
getneqconstraints

"""
    getnvars(model::AbstractModel)

Returns the number of variables in `model`.
"""
getnvars

"""
    getmin(model::AbstractModel)

Returns the variables' lower bounds of `model`.
"""
getmin

"""
    getmax(model::AbstractModel)

Returns the variables' upper bounds of `model`.
"""
getmax

"""
    setmin!(model::AbstractModel, min::Union{Number, AbstractVector})
    setmin!(model::AbstractModel, i::Integer, min::Number)

Sets the variables' lower bounds in `model` to `min`. Or if `i` is given, it sets the lower bound of the `i`^th variable in `model` to `min`.
"""
setmin!

"""
    setmax!(model::AbstractModel, max::Union{Number, AbstractVector})
    setmax!(model::AbstractModel, i::Integer, max::Number)

Sets the variables' upper bound in `model` to `max`. Or if `i` is given, it sets the upper bound of the `i`^th variable in `model` to `max`.
"""
setmax!

"""
    setbox!(model::AbstractModel, min::Union{Number, AbstractVector}, max::Union{Number, AbstractVector})
    setbox!(model::AbstractModel, i::Integer, min::Number, max::Number)

Sets the variables' lower and upper bounds in `model` to `min` and `max` respectively. Or if `i` is given, it sets the lower and upper bounds of the `i`^th variable in `model` to `min` and `max` respectively.
"""
setbox!

"""
    addvar!(model::AbstractModel, lb::Number, ub::Number; init = lb, integer = false)

Adds a new variable to `model` with lower and upper bounds `lb` and `ub`, initial value `init` and `integer = true` if the variable is integer.
"""
addvar!

"""
    add_ineq_constraint!(model::AbstractModel, f::IneqConstraint)
    add_ineq_constraint!(model::AbstractModel, fs::Vector{<:IneqConstraint})

Adds a new constraint `f` or multiple new constraints `fs` to `model`.
"""
add_ineq_constraint!

"""
    init!(model::AbstractModel)

Initializes the model.
"""
init!


text:
struct LazyArrayStyle{N} <: AbstractArrayStyle{N} end
LazyArrayStyle(::Val{N}) where N = LazyArrayStyle{N}()
LazyArrayStyle{M}(::Val{N}) where {N,M} = LazyArrayStyle{N}()

struct BroadcastArray{T, N, BRD<:Broadcasted} <: AbstractArray{T, N}
    broadcasted::BRD
end

BroadcastArray{T,N}(bc::BRD) where {T,N,BRD<:Broadcasted} = BroadcastArray{T,N,BRD}(bc)
BroadcastArray{T}(bc::Broadcasted{<:Union{Nothing,BroadcastStyle},<:Tuple{Vararg{Any,N}},<:Any,<:Tuple}) where {T,N} =
    BroadcastArray{T,N}(bc)

_broadcast2broadcastarray(a, b...) = tuple(a, b...)
_broadcast2broadcastarray(a::Broadcasted, b...) = tuple(BroadcastArray(a), b...)

_BroadcastArray(bc::Broadcasted) = BroadcastArray{combine_eltypes(bc.f, bc.args)}(bc)
BroadcastArray(bc::Broadcasted{S}) where S =
    _BroadcastArray(instantiate(Broadcasted{S}(bc.f, _broadcast2broadcastarray(bc.args...))))
BroadcastArray(b::BroadcastArray) = b
BroadcastArray(f, A, As...) = BroadcastArray(broadcasted(f, A, As...))

axes(A::BroadcastArray) = axes(A.broadcasted)
size(A::BroadcastArray) = map(length, axes(A))

IndexStyle(::BroadcastArray{<:Any,1}) = IndexLinear()

@propagate_inbounds getindex(A::BroadcastArray, kj::Int...) = A.broadcasted[kj...]

@propagate_inbounds _broadcast_getindex_range(A::Union{Ref,AbstractArray{<:Any,0},Number}, I) = A[] # Scalar-likes can just ignore all indices
# Everything else falls back to dynamically dropping broadcasted indices based upon its axes
@propagate_inbounds _broadcast_getindex_range(A, I) = A[I]

getindex(B::BroadcastArray{<:Any,1}, kr::AbstractVector{<:Integer}) =
    BroadcastArray(B.broadcasted.f, map(a -> _broadcast_getindex_range(a,kr), B.broadcasted.args)...)

copy(bc::Broadcasted{<:LazyArrayStyle}) = BroadcastArray(bc)

BroadcastStyle(::Type{<:BroadcastArray{<:Any,N}}) where N = LazyArrayStyle{N}()
BroadcastStyle(L::LazyArrayStyle{N}, ::StaticArrayStyle{N}) where N = L
BroadcastStyle(::StaticArrayStyle{N}, L::LazyArrayStyle{N})  where N = L

## scalar-range broadcast operations ##
# Ranges already support smart broadcasting
for op in (+, -, big)
    @eval begin
        broadcasted(::LazyArrayStyle{1}, ::typeof($op), r::AbstractRange) =
            broadcast(DefaultArrayStyle{1}(), $op, r)
    end
end

for op in (-, +, *, /)
    @eval broadcasted(::LazyArrayStyle{1}, ::typeof($op), r::AbstractRange, x::Real) =
        broadcast(DefaultArrayStyle{1}(), $op, r, x)
end

for op in (-, +, *, \)
    @eval broadcasted(::LazyArrayStyle{1}, ::typeof($op), x::Real, r::AbstractRange) =
        broadcast(DefaultArrayStyle{1}(), $op, x, r)
end

broadcasted(::LazyArrayStyle{N}, op, r::AbstractFill{T,N}) where {T,N} =
    broadcast(DefaultArrayStyle{N}(), op, r)
broadcasted(::LazyArrayStyle{N}, op, r::AbstractFill{T,N}, x::Number) where {T,N} =
    broadcast(DefaultArrayStyle{N}(), op, r, x)
broadcasted(::LazyArrayStyle{N}, op, x::Number, r::AbstractFill{T,N}) where {T,N} =
    broadcast(DefaultArrayStyle{N}(), op, x, r)
broadcasted(::LazyArrayStyle{N}, op, r1::AbstractFill{T,N}, r2::AbstractFill{V,N}) where {T,V,N} =
    broadcast(DefaultArrayStyle{N}(), op, r1, r2)

broadcasted(::LazyArrayStyle{N}, ::typeof(*), a::Zeros{T,N}, b::Zeros{V,N}) where {T,V,N} =
    broadcast(DefaultArrayStyle{N}(), *, a, b)
broadcasted(::LazyArrayStyle{N}, ::typeof(*), a::AbstractArray{T,N}, b::Zeros{V,N}) where {T,V,N} =
    broadcast(DefaultArrayStyle{N}(), *, a, b)
broadcasted(::LazyArrayStyle{N}, ::typeof(*), a::Zeros{T,N}, b::AbstractArray{V,N}) where {T,V,N} =
    broadcast(DefaultArrayStyle{N}(), *, a, b)


text:
@everywhere begin

using Dagger

struct FakeProc <: Dagger.Processor
    owner::Int
end
FakeProc() = FakeProc(myid())
Dagger.get_parent(proc::FakeProc) = OSProc(proc.owner)
Base.:(==)(proc1::FakeProc, proc2::FakeProc) = proc1.owner == proc2.owner

struct FakeVal
    x
end
fakesum(xs...) = FakeVal(sum(map(y->y.x, xs)))

Dagger.iscompatible_func(proc::FakeProc, opts, f) = true
Dagger.iscompatible_arg(proc::FakeProc, opts, x::Integer) = true
Dagger.iscompatible_arg(proc::FakeProc, opts, x::FakeVal) = true
Dagger.move(ctx, from_proc::OSProc, to_proc::FakeProc, x::Integer) = FakeVal(x)
Dagger.move(ctx, from_proc::FakeProc, to_proc::OSProc, x::Vector) =
    map(y->y.x, x)
Dagger.move(ctx, from_proc::FakeProc, to_proc::OSProc, x::FakeVal) =
    x.x
Dagger.execute!(proc::FakeProc, func, args...) = FakeVal(42+func(args...).x)

end


text:
struct ADTestFunction
    name::String
    f::Function
    x::Vector
end
struct DistSpec{Tθ<:Tuple, Tx}
    name::Union{Symbol, Expr}
    θ::Tθ
    x::Tx
end

vectorize(v::Number) = [v]
vectorize(v) = vec(v)
pack(vals...) = reduce(vcat, vectorize.(vals))
function unpack(x, vals...)
    unpacked = []
    i = 1
    for v in vals
        if v isa Number
            push!(unpacked, x[i])
            i += 1
        elseif v isa Vector
            push!(unpacked, x[i:i+length(v)-1])
            i += length(v)
        elseif v isa Matrix
            push!(unpacked, reshape(x[i:(i+length(v)-1)], size(v)))
            i += length(v)
        else
            throw("Unsupported argument")
        end
    end
    @assert i == length(x) + 1
    return (unpacked...,)
end
function get_function(dist::DistSpec, inds, val)
    syms = []
    args = []
    for (i, a) in enumerate(dist.θ)
        if i in inds
            sym = gensym()
            push!(syms, sym)
            push!(args, sym)
        else
            push!(args, a)
        end
    end
    if val
        sym = gensym()
        push!(syms, sym)
        expr = :(($(syms...),) -> logpdf($(dist.name)($(args...)), $(sym)))
        if length(inds) == 0
            f = x -> Base.invokelatest(eval(expr), unpack(x, dist.x)...)
            return ADTestFunction(string(expr), f, pack(dist.x))
        else
            f = x -> Base.invokelatest(eval(expr), unpack(x, dist.θ[inds]..., dist.x)...)
            return ADTestFunction(string(expr), f, pack(dist.θ[inds]..., dist.x))
        end
    else
        @assert length(inds) > 0
        expr = :(($(syms...),) -> logpdf($(dist.name)($(args...)), $(dist.x)))
        f = x -> Base.invokelatest(eval(expr), unpack(x, dist.θ[inds]...)...)
        return ADTestFunction(string(expr), f, pack(dist.θ[inds]...))
    end
end
function get_all_functions(dist::DistSpec, continuous=false)
    fs = []
    if length(dist.θ) == 0
        push!(fs, get_function(dist, (), true))
    else
        for inds in combinations(1:length(dist.θ))
            push!(fs, get_function(dist, inds, false))
            if continuous
                push!(fs, get_function(dist, inds, true))
            end
        end
    end
    return fs
end

function test_ad(f, at = 0.5; rtol = 1e-8, atol = 1e-8)
    isarr = isa(at, AbstractArray)
    reverse = Tracker.data(Tracker.gradient(f, at)[1])
    if isarr
        forward = ForwardDiff.gradient(f, at)
        @test isapprox(reverse, forward, rtol=rtol, atol=atol)
    else
        forward = ForwardDiff.derivative(f, at)
        finite_diff = central_fdm(5,1)(f, at)
        @test isapprox(reverse, forward, rtol=rtol, atol=atol)
        @test isapprox(reverse, finite_diff, rtol=rtol, atol=atol)
    end
end

"""
    test_reverse_mode_ad(forward, f, ȳ, x...; rtol=1e-8, atol=1e-8)

Check that the reverse-mode sensitivities produced by an AD library are correct for `f`
at `x...`, given sensitivity `ȳ` w.r.t. `y = f(x...)` up to `rtol` and `atol`.
`forward` should be either `Tracker.forward` or `Zygote.forward`.
"""
function test_reverse_mode_ad(forward, f, ȳ, x...; rtol=1e-8, atol=1e-8)

    # Perform a regular forwards-pass.
    y = f(x...)

    # Use tracker to compute reverse-mode sensitivities.
    y_tracker, back = forward(f, x...)
    x̄s_tracker = back(ȳ)

    # Use finite differencing to compute reverse-mode sensitivities.
    x̄s_fdm = FDM.j′vp(central_fdm(5, 1), f, ȳ, x...)
    if length(x) == 1
        x̄s_fdm = (x̄s_fdm,)
    end

    # Check that forwards-pass produces the correct answer.
    @test y ≈ y_tracker

    # Check that reverse-mode sensitivities are correct.
    @test all([x̄_tracker ≈ x̄_fdm for (x̄_tracker, x̄_fdm) in zip(x̄s_tracker, x̄s_fdm)])
end

# See `test_reverse_mode_ad` for details.
function test_tracker_ad(f, ȳ, x...; rtol=1e-8, atol=1e-8)
    return test_reverse_mode_ad(Tracker.forward, f, ȳ, x...; rtol=rtol, atol=atol)
end

_to_cov(B) = B * B' + Matrix(I, size(B)...)


text:
struct Reporter{T,F<:Function}
    name::String
    call::F
end

"""
    Reporter(name, datatype, call)

Return a `Reporter` that specifies a piece of data that should be collected.

# Arguments
* `name::String`: specifies the name of the piece of data.
* `datatype::DataType`: specifies the datatype of the piece of data.
* `call::Function`: specifies a Function that returns the corresponding value when called as
    `f(args...)`.
"""
function Reporter(name::String, datatype::DataType, call::F) where F<:Function
    return Reporter{datatype,F}(name, call)
end

struct Collector{T,C,F1<:Function,F2<:Function,F3<:Function,F4<:Function}
    names::Vector{String} # Holds the reporter names.
    calls::C # Holds the reporter calls.
    iter::F1 # Returns an iterable when called as iter(model, agents, exp).
    condition::F2 # Specifies the conditions under which calls to update the
                  # collector should proceed. Continues if and only if `condition(model,
                  # agents, exp)` returns `true`.
    prepare::F3 # Specifies algorithms to run right before the collector collects data.
                # Called as `prepare(model, agents, exp)`.
    finish::F4 # Specifies algorithms that run right after the collector collects data.
               # Called as `finish(model, agents, exp)`.
end

"""
    Collector(reporters...; iter = passthrough, condition = always_true, prepare = always_nothing,
              finish = always_nothing)

Return a `Collector` that collects data for the reporters.

# Arguments
* `reporters::Reporter...`: one or more `Reporter`s that specify data that should be collected.
* `iter::Function`: specifies a Function that provides an iterable that provides the arguments
    with which to call the reporters. The function is called as `iter(args...)`. Defaults to a
    function that passes the arguments directly.
* `condition::Function`: specifies the conditions under which the collector should collect data.
    Calls to collect data with the collector will only proceed when `condition(args...)` returns
    `true`. Defaults to a function that always returns `true`.
* `prepare::Function`: specifies a function to run right before the collector collects data, but
    after the check that data should be collected. The function is called as `prepare(args...)`.
    Defaults to an empty function.
* `finish::Function`: specifies a function to run right after the collector collected data, but
    only if data was actually collected. The function is called as `finish(args...)`. Defaults to
    an empty function.
"""
function Collector(reporters::Reporter...;
                   iter::F1 = passthrough, condition::F2 = always_true,
                   prepare::F3 = always_nothing, finish::F4 = always_nothing) where {F1,F2,F3,F4}
    T = Tuple{(typeof(reporter).parameters[1] for reporter in reporters)...}
    C = Tuple{(typeof(reporter).parameters[2] for reporter in reporters)...}
    names = [reporter.name for reporter in reporters]
    calls = ((reporter.call for reporter in reporters)..., )
    return Collector{T,C,F1,F2,F3,F4}(names, calls, iter, condition, prepare, finish)
end


text:
using PolynomialZeros
const AGCD = PolynomialZeros.AGCD
const MultRoot = PolynomialZeros.MultRoot
using Polynomials
using Test

poly_coeffs = PolynomialZeros.poly_coeffs
multroot=PolynomialZeros.MultRoot.multroot
pejroot=PolynomialZeros.MultRoot.pejroot
agcd = PolynomialZeros.AGCD.agcd
identify_z0s_ls = PolynomialZeros.MultRoot.identify_z0s_ls

x = variable(Float64)
bx = variable(BigFloat)
_poly(zs,ls,x=x) = prod((x-z)^l for (z,l) in zip(zs, ls))

@testset "agcd-Float64" begin

    x = variable()

    p = (x-1)^3
    u,v,w,err = AGCD.agcd(p)
    @test Polynomials.degree(v) == 1

    p = prod(x-i for i in 1:6)
    u,v,w,err = AGCD.agcd(p)
    @test Polynomials.degree(v) == Polynomials.degree(p)

    n = 4
    p = prod((x-i)^i for i in 1:n)
    u,v,w,err = AGCD.agcd(p)
    @test Polynomials.degree(v) == n

    # can fails
    n = 6
    p = prod((x-i)^i for i in 1:n)
    u,v,w,err = AGCD.agcd(p)
    @test Polynomials.degree(v) >= n

    # use big
    n = 6
    p = prod((bx-i)^i for i in 1:n)
    u,v,w,err = AGCD.agcd(p)
    @test Polynomials.degree(v) == n

    T = Float64
    x = variable(Complex{T})
    p = (x-im)^2 * (x-1)^2 * (x+2im)^2
    u,v,w,err = AGCD.agcd(p)
    @test length(v)-1 == 3

end

@testset "agcd-othertypes" begin

    T = Float32
    zs, ls = [1,2,-3],[2,2,3]
    p = _poly(zs, ls, variable(T))
    ps = poly_coeffs(p)
    u,v,w,err = AGCD.agcd(p, θ=1e-4)

    @test degree(v) == 3

    # big float
    p = _poly(zs, ls, bx)
    u,v,w,err = AGCD.agcd(p^15)

    @test degree(v) == 3 # fails for Float64

end

@testset "pejroot" begin

    zs, ls = [1.0,2,3,4], [4,3,2,1]
    p = _poly(zs, ls,x)

    delta = 0.1 # works
    z0 = zs + delta*[1,-1,1,-1]
    z1 = pejroot(p, z0, ls)
    @test !all(sort(z0) .== sort(z1))

    delta = 0.2 # fails
    z0 = zs + delta*[1,-1,1,-1]
    z1 = pejroot(p, z0, ls)
    @test all(sort(z0) .== sort(z1))

    ls = [30,20,10,5]
    p = _poly(zs, ls,x)

    delta = 0.01 # works
    z0 = zs + delta*[1,-1,1,-1]
    z1 = pejroot(p, z0, ls)
    @test !all(sort(z0) .== sort(z1))

end

@testset "identify_ls" begin

    zs, ls = [1.0,2,3,4], [4,3,2,1]
    p = _poly(zs, ls,x)
    _zs, _ls = identify_z0s_ls(poly_coeffs(p))
    @test all(sort(ls) .== sort(_ls))

    p = _poly(zs, 3*ls,x)
    _zs, _ls = identify_z0s_ls(poly_coeffs(p))
    @test length(ls) != length(_ls) # fails w/o preconditioning

    n = 4
    zs, ls = cumsum(ones(n)), cumsum(ones(Int, n))
    p = _poly(zs, ls,x)
    _zs, _ls = identify_z0s_ls(poly_coeffs(p))
    @test all(sort(ls) .== sort(_ls))

    n = 5
    zs, ls = cumsum(ones(n)), cumsum(ones(Int, n))
    p = _poly(zs, ls,x)
    _zs, _ls = identify_z0s_ls(poly_coeffs(p))
    @test !(length(ls) == length(_ls))

end

@testset "multroot" begin

    zs, ls = [1.0,2,3,4], [4,3,2,1]
    _zs, _ls = multroot(_poly(zs, ls,x))
    @test all(sort(ls) .== sort(_ls))

    _zs, _ls = multroot(_poly(zs, 3ls,x))
    @test !all(sort(2ls) .== sort(_ls)) # XXX fails!

    _zs, _ls = multroot(_poly(big.(zs), 3ls,bx))
    @test all(sort(3ls) .== sort(_ls)) # passes

    delta = 0.01
    zs, ls = [1-delta, 1, 1+delta], [5,4,3]
    _zs, _ls = multroot(_poly(zs, ls,x))
    @test all(sort(ls) .== sort(_ls))

    # should work, but giving issue with CI
    # n = 20
    # zs,ls = collect(1.0:n), ones(Int, n)
    # _zs, _ls = multroot(_poly(zs, ls,x))
    # @test all(sort(ls) .== sort(_ls))

    _zs, _ls = multroot(_poly(zs, 2ls,x))
    @test !(length(ls) == length(_ls)) ##XXX fails!

    n = 10
    zs,ls = collect(1.0:n), ones(Int, n)
    _zs, _ls = multroot(_poly(zs, 2ls,x)) #fails *badly*
    @test !(length(_ls) == length(ls))

    n = 10
    zs,ls = collect(1.0:n), ones(Int, n)
    _zs, _ls = multroot(_poly(big.(zs), 2ls,x)) # works now
    @test all(sort(2ls) .== sort(_ls))

    n = 5
    zs,ls = collect(1.0:n), ones(Int, n)
    _zs, _ls = multroot(_poly(zs, 4ls,x))
    @test all(sort(4ls) .== sort(_ls))

end

@testset "examples" begin
    # examples from paper
    # our agcd implementation (and identify_z0s_ls) do not perform
    # as advertised therein.

    # 3.6.1
    zs, ls = [1.0,2,3,4], [4,3,2,1]
    p1 = _poly(zs, ls, x)
    ps = poly_coeffs(p1)
    z0 = zs + 0.1 * [1,-1,1,-1]
    @test maximum(abs.(sort(pejroot(ps, z0, ls)) - [1,2,3,4])) <= 1e-10

    ls = [40,30,20,10]
    p2 = _poly(zs, ls, x)
    ps = poly_coeffs(p2)
    z0 = zs + 0.01 * [1,-1,1,-1]
    @test maximum(abs.(sort(pejroot(ps, z0, ls)) - [1,2,3,4])) <= 1e-13

    zs, ls = [sqrt(2), sqrt(3)], [20,10]
    p3 = _poly(zs, ls, x)
    ps = poly_coeffs(p3)
    zs_, ls_ = multroot(ps)
    @test maximum(abs.(sort(zs_) - sort(zs))) <= 1e-14

    # 3.62
    delta = 0.1
    zs, ls = [1,1,1]-delta*[-1,0,1], [9,5,2]
    p4 = _poly(zs, ls, x)
    ps = poly_coeffs(p4)
    zs_,ls_ = multroot(ps)
    @test  maximum(abs.(sort(zs_) - sort(zs))) <= 1e-13

    zs, ls = [1,1,1]+delta*[-1,0,1], 2 * [9,5,2] # from paper
    p4 = _poly(zs, ls, x)
    ps = poly_coeffs(p4)
    zs_,ls_ = multroot(ps)
    # @test  maximum(abs.(sort(zs_) - zs)) <= 1e-14 # fails, not even correct size for zs

    # pejroots
    z0 = [0.89999999993, 0.9999999993, 1.0999999998]
    zs_ = pejroot(ps, z0, ls)
    @test maximum(abs.(sort(zs_) - sort(zs))) <= 1e-10

    # will work with BigFloat
    zs, ls = [1,1,1]+delta*[-1,0,1], 2 * [9,5,2] # from paper
    p5 = _poly(zs, ls, bx)
    ps = poly_coeffs(p5)
    zs_,ls_ = multroot(ps)
    @test  maximum(abs.(sort(zs_) - zs)) <= 1e-50

    zs, ls = [0.3 + 0.6im, 0.1 + 0.7im, 0.7+0.5im, 0.3 + 0.4im], [100, 200, 300, 400]
    p7 = _poly(zs, ls, variable(Complex{Float64}))
    ps = poly_coeffs(p7)
    ps = ps + rand(-1:1, length(ps)) * 1e-6
    z0s = [.289+.601im, .100 + .702im, .702 + .498im, .301 + .399im]
    zs_ = pejroot(ps, z0s, ls)
    @test maximum(abs.(sort(abs.(zs_)) - sort(abs.(zs)))) <= 1e-5

    # 4.6
    #
    # k = 5 is in paper, we can't even get k=4. In fact k=3 not
    # working without preconditioning
    k=3
    zs, ls = [1,2,3,4], k*[4,3,2,1]
    p8 = _poly(zs, ls, x)
    ps = poly_coeffs(p8)
    zs_, ls_ = multroot(ps)
    @test all(sort(ls_) .== sort(ls))

    zs, ls = [10/11, 20/11, 30/11], [5,5,5]
    p9 = _poly(zs, ls, x)
    ps = poly_coeffs(p9)
    zs_, ls_ = multroot(ps)
    @test all(sort(ls_) .== sort(ls))

    # 5.2 -- these fail
    for ε in [1e-1, 1e-2, 1e-3]
        zs, ls = [1-ε, 1, 1/2], [20,20,5]
        p10 = _poly(zs, ls, x)
        zs_, ls_ = multroot(p10)
#        @test all(sort(ls) .== sort(ls_))
    end

    ε = .0001
     zs, ls = [1-ε, 1, 1/2], [20,20,5]
    p10 = _poly(zs, ls, x)
    zs_, ls_ = multroot(p10,  ϕ = 5)
#    @test all(sort(ls) .== sort(ls_))

end


text:
# ------------------------------------------------------------------------------------------
# ## 上手 Jupyter notebook
#
# ### 目录
# - [运行代码块（cell）](#运行代码块（cell）)
# - [查看 Julia 的帮助文档](#使用-shell-命令)
# - [使用 shell 命令](#使用-shell-命令)
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### 运行代码块（cell）
# 选中一个代码块之后（选中的代码块会被绿色矩形框起来），有以下几种方法可以运行它：
# 1. 按下组合键 `Shift` + `Enter`
# 2. 点击上方工具栏中的 **运行（Run）** 按钮
# 3. 按下组合键 `Ctrl` + `Enter`
# ------------------------------------------------------------------------------------------

1 + 1
2 + 2

# ------------------------------------------------------------------------------------------
# 新手注意！执行一个代码块时默认打印最后一行。想要什么也不输出，在最后一行的行尾加上分号 `;` 就行了。如下所示
# ------------------------------------------------------------------------------------------

1 + 1
2 + 2;

# ------------------------------------------------------------------------------------------
# ### 查看 Julia 的帮助文档
#
# 遇到不熟悉的 Julia 函数或语法，在函数名或表达式前加个问号 `?` 就可以查询对应的帮助文档。（在Julia的REPL中也好使哦！）
# ------------------------------------------------------------------------------------------

?println

# ------------------------------------------------------------------------------------------
# ### 使用 shell 命令
#
# 在 shell 命令前加上分号 `;` 即可，像这样：
#
# > 【译注】以下命令仅在 macOS 和 *nix 系统下可用
# ------------------------------------------------------------------------------------------

;ls

;pwd

# ------------------------------------------------------------------------------------------
# shell 命令在 Julia 的 REPL 中也好使哦！
#
# > 【译注】对于 Windows 系统，你需要使用 Windows cmd 命令，像这样：
# ------------------------------------------------------------------------------------------

;cd

# ------------------------------------------------------------------------------------------
# 但并不是所有的 cmd 命令都可以使用。如 `dir` 命令不可用，会报错
# ------------------------------------------------------------------------------------------

;dir

# ------------------------------------------------------------------------------------------
# 中文社区的讨论贴
# - [Windows 中 Julia 的 shell 模式与 `;dir` 报错 - 综合讨论区 / 心得体会 -
# Julia中文社区](https://discourse.juliacn.com/t/topic/2890)
# ------------------------------------------------------------------------------------------

max_stars_repo_name: NahsiN/JuliaTutorials

max_issues_repo_name: NahsiN/JuliaTutorials

max_forks_repo_name: NahsiN/JuliaTutorials


text:
# code for Video 5.5:
julia> using RDatasets
julia> iris = dataset("datasets", "iris")
julia> labels = iris[:Species]
150-element DataArrays.PooledDataArray{ASCIIString,UInt8,1}:
 "setosa"
 "setosa"
 "setosa"
 ...
 julia> features = convert(Array, iris[:, 1:4])
150x4 Array{Float64,2}:
 5.1  3.5  1.4  0.2
 4.9  3.0  1.4  0.2
 4.7  3.2  1.3  0.2
 4.6  3.1  1.5  0.2
 5.0  3.6  1.4  0.2

julia> Pkg.add("LIBSVM")
julia> using LIBSVM

features_train, features_test = features[1:2:end, :], features[2:2:end, :]
(
75x4 Array{Float64,2}:
 5.1  3.5  1.4  0.2
 4.7  3.2  1.3  0.2
 5.0  3.6  1.4  0.2
 4.6  3.4  1.4  0.3
 4.4  2.9  1.4  0.2
 5.4  3.7  1.5  0.2
 4.8  3.0  1.4  0.1
 5.8  4.0  1.2  0.2
 5.4  3.9  1.3  0.4
 ...
 )
75x4 Array{Float64,2}:
 4.9  3.0  1.4  0.2
 4.6  3.1  1.5  0.2
 5.4  3.9  1.7  0.4
 5.0  3.4  1.5  0.2
 4.9  3.1  1.5  0.1
 4.8  3.4  1.6  0.2
 4.3  3.0  1.1  0.1
 ...
 )

julia> labels_train, labels_test = labels[1:2:end], labels[2:2:end]
(ASCIIString["setosa","setosa","setosa","setosa","setosa","setosa","setosa","seto
sa","setosa","setosa"  …  "virginica","virginica","virginica","virginica","virgin
ica","virginica","virginica","virginica","virginica","virginica"],ASCIIString["se
tosa","setosa","setosa","setosa","setosa","setosa","setosa","setosa","setosa","se
tosa"  …  "virginica","virginica","virginica","virginica","virginica","virginica"
,"virginica","virginica","virginica","virginica"])

julia> features_train' # we need a transposed matrix for our SVM model
4x75 Array{Float64,2}:
 5.1  4.7  5.0  4.6  4.4  5.4  4.8  5.8  5.4  5.7  …  7.4  6.4  6.1  6.3  6.0  6.7  5.8  6.7  6.3  6.2
 3.5  3.2  3.6  3.4  2.9  3.7  3.0  4.0  3.9  3.8     2.8  2.8  2.6  3.4  3.0  3.1  2.7  3.3  2.5  3.4
 1.4  1.3  1.4  1.4  1.4  1.5  1.4  1.2  1.3  1.7     6.1  5.6  5.6  5.6  4.8  5.6  5.1  5.7  5.0  5.4
 0.2  0.2  0.2  0.3  0.2  0.2  0.1  0.2  0.4  0.3     1.9  2.2  1.4  2.4  1.8  2.4  1.9  2.5  1.9  2.3
# Rows: 1) SepalLength, 2) SepalWidth, 3) PetalLength, 4) PetalWidth

julia> model = svmtrain(labels_train, features_train')
LIBSVM.SVMModel{ASCIIString}(Ptr{Void} @0x00000000044e9620,[LIBSVM.SVMParameter(0,2,3,0.25,0.0,100.0,0.001,1.0,0,Ptr{Int32} @0x00007ef9f467d800,Ptr{Float64} @0x00007ef9f467d840,0.5,0.1,1,0)],[LIBSVM.SVMProblem(75,Ptr{Float64} @0x00007ef9f3f37a00,Ptr{Ptr{LIBSVM.SVMNode}} @0x00007ef9f4bc0040)],5x75 Array{LIBSVM.SVMNode,2}:
 LIBSVM.SVMNode(1,5.1)   LIBSVM.SVMNode(1,4.7)   …  LIBSVM.SVMNode(1,6.2)
 LIBSVM.SVMNode(2,3.5)   LIBSVM.SVMNode(2,3.2)      LIBSVM.SVMNode(2,3.4)
 LIBSVM.SVMNode(3,1.4)   LIBSVM.SVMNode(3,1.3)      LIBSVM.SVMNode(3,5.4)
 LIBSVM.SVMNode(4,0.2)   LIBSVM.SVMNode(4,0.2)      LIBSVM.SVMNode(4,2.3)
 LIBSVM.SVMNode(-1,NaN)  LIBSVM.SVMNode(-1,NaN)     LIBSVM.SVMNode(-1,NaN),[Ptr{LIBSVM.SVMNode} @0x00000000044e7710,Ptr{LIBSVM.SVMNode} @0x00000000044e7760,Ptr{LIBSVM.SVMNode} @0x00000000044e77b0,Ptr{LIBSVM.SVMNode} @0x00000000044e7800,Ptr{LIBSVM.SVMNode} @0x00000000044e7850,Ptr{LIBSVM.SVMNode} @0x00000000044e78a0,Ptr{LIBSVM.SVMNode} @0x00000000044e78f0,Ptr{LIBSVM.SVMNode} @0x00000000044e7940,Ptr{LIBSVM.SVMNode} @0x00000000044e7990,Ptr{LIBSVM.SVMNode} @0x00000000044e79e0  …  Ptr{LIBSVM.SVMNode} @0x00000000044e8b60,Ptr{LIBSVM.SVMNode} @0x00000000044e8bb0,Ptr{LIBSVM.SVMNode} @0x00000000044e8c00,Ptr{LIBSVM.SVMNode} @0x00000000044e8c50,Ptr{LIBSVM.SVMNode} @0x00000000044e8ca0,Ptr{LIBSVM.SVMNode} @0x00000000044e8cf0,Ptr{LIBSVM.SVMNode} @0x00000000044e8d40,Ptr{LIBSVM.SVMNode} @0x00000000044e8d90,Ptr{LIBSVM.SVMNode} @0x00000000044e8de0,Ptr{LIBSVM.SVMNode} @0x00000000044e8e30],ASCIIString["setosa","versicolor","virginica"],Int32[],Float64[],4,false)

julia> (predicted_labels, decision_values) = svmpredict(model, features_test')

(ASCIIString["setosa","setosa","setosa","setosa","setosa","setosa","setosa","setosa","setosa","setosa"  …  "virginica","versicolor","virginica","virginica","virginica","virginica","virginica","virginica","virginica","virginica"],
3x75 Array{Float64,2}:
 1.07954   1.07629   1.05125   1.16964   …  -0.809238  -0.897749  -0.936477
 1.10042   1.07671   1.03685   1.1533       -1.00333   -1.0383    -1.07345
 0.263692  0.254304  0.358478  0.282262     -0.810783  -0.571849  -0.315681)

julia> @printf "Accuracy: %.2f%%\n" mean((predicted_labels .== labels_test))*100
Accuracy: 93.33%

julia> using DecisionTree

julia> confusion_matrix(convert(Array, labels_test), predicted_labels)
Classes:  ASCIIString["setosa","versicolor","virginica"]
Matrix:   3x3 Array{Int64,2}:
 25   0   0
  0  24   1
  0   4  21

Accuracy: 0.9333333333333333
Kappa:    0.9

### Take 80% of records (120) for training and 20% of records (30) for testing
julia> julia> n_test = Int(length(labels) * 0.2) # 30

julia> collect(1:length(labels)) .> n_test
150-element BitArray{1}:
 false   # first 30 boolean values are false
 false
 false
 false
 false
 false
 false
 false
 false
 false
 false
 false
     ⋮
  true
  true

julia> train_rows = shuffle(collect(1:length(labels)) .> n_test) # shuffle randomizes the true and false values
150-element BitArray{1}:
  true  # shuffle randomizes the true and false values
  true
  true
  true
  true
 false
  true
  ...

julia> features_train, features_test = features[train_rows, :], features[!train_rows, :]
julia> labels_train, labels_test = labels[train_rows], labels[!train_rows]

julia> model = svmtrain(labels_train, features_train')
julia> (predictions, decision_values) = svmpredict(model, features_test')

julia> confusion_matrix(labels_test, predictions)
Classes:  ASCIIString["setosa","versicolor","virginica"]
Matrix:   3x3 Array{Int64,2}:
 15  0  0
  0  9  1
  0  0  5
Accuracy: 0.9666666666666667
Kappa:    0.9459459459459458

max_stars_repo_path: Section 5/5_5.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 5/5_5.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 5/5_5.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
# Copyright 2017,2018 Eric S. Tellez
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

import SimilaritySearch: optimize!
export NaiveBayesClassifier, NaiveBayesKernel, kernel_prob, predict, predict_proba, optimize!

abstract type NaiveBayesKernel end

mutable struct NaiveBayesClassifier{ItemType,LabelType}
    kernel::NaiveBayesKernel
    probs::Vector{Float64}
    le::LabelEncoder{LabelType}
end

include("nbgaussian.jl")
include("nbmultinomial.jl")

function NaiveBayesClassifier(X::AbstractVector{ItemType}, y::AbstractVector{LabelType}; kernel=GaussianKernel) where {ItemType,LabelType}
    le = fit(LabelEncoder, y)
    # y_ = [transform(le, l) for l in y]
    y_ = transform.(le, y)
    probs = Float64[freq/length(y) for freq in le.freqs]
    kernel_ = kernel(X, y_, length(le.labels))
    NaiveBayesClassifier{ItemType,LabelType}(kernel_, probs, le)
end

function predict(nbc::NaiveBayesClassifier{ItemType,LabelType}, vector)::Vector{LabelType} where {ItemType,LabelType}
    y = Vector{LabelType}(undef, length(vector))
    for i in 1:length(vector)
        y[i] = predict_one(nbc, vector[i])
    end

    y
end

function predict_one_proba(nbc::NaiveBayesClassifier{ItemType,LabelType}, x) where {ItemType,LabelType}
    w = kernel_prob(nbc, nbc.kernel, x)
    ws = sum(w)
    @inbounds for i in 1:length(w)
        w[i] /= ws
    end
    w
end

function predict_one(nbc::NaiveBayesClassifier{ItemType,LabelType}, x) where {ItemType,LabelType}
    p, i = findmax(kernel_prob(nbc, nbc.kernel, x))
    inverse_transform(nbc.le, i)
end

function optimize!(nbc::NaiveBayesClassifier{ItemType,LabelType}, X::AbstractVector{ItemType}, y::AbstractVector{LabelType}, scorefun::Function; runs=3, trainratio=0.5, testratio=0.5, folds=0, shufflefolds=true) where {ItemType,LabelType}
    @info "optimizing nbc $(typeof(nbc))"
    # y::Vector{Int} = transform.(nbc.le, y)
    mem = Dict{Any,Float64}()
    function f(train_X, train_y, test_X, test_y)
        tmp = NaiveBayesClassifier(train_X, train_y, kernel=MultinomialKernel)
        for smoothing in [0.0, 0.1, 0.3, 1.0]
            tmp.kernel.smoothing = smoothing
            pred_y = predict(tmp, test_X)
            score = scorefun(test_y, pred_y)
            mem[(MultinomialKernel, smoothing)] = get(mem, (MultinomialKernel, smoothing), 0.0) + score
        end

        tmp.kernel = GaussianKernel(train_X, transform.(tmp.le, train_y), length(tmp.le.labels))
        pred_y = predict(tmp, test_X)
        score = scorefun(test_y, pred_y)
        mem[(GaussianKernel, -1)] = get(mem, (GaussianKernel, -1), 0.0) + score
        0
    end

    if folds > 1
        kfolds(f, X, y, folds=folds, shuffle=shufflefolds)
        bestlist = [(score/folds, conf) for (conf, score) in mem]
    else
        montecarlo(f, X, y, runs=runs, trainratio=trainratio, testratio=testratio)
        bestlist = [(score/runs, conf) for (conf, score) in mem]
    end

    sort!(bestlist, by=x -> (-x[1], x[2][2]))
    best = bestlist[1][2]
    if best[1] == GaussianKernel
        nbc.kernel = GaussianKernel(X, transform.(nbc.le, y), length(nbc.le.labels))
    else
        nbc.kernel = MultinomialKernel(X, transform.(nbc.le, y), length(nbc.le.labels), smoothing=best[2])
    end

    bestlist
end


text:
apply(A::AbstractMatrix, x::AbstractVector) = A * x
apply(f, x) = f(x)

function apply(operator, x, α₀, α₁)
    y = apply(operator, x)
    if α₀ != zero(α₀) || α₁ != one(α₁)
        axpby!(α₀, x, α₁, y)
    end

    return y
end

# GKL, SVD, LSMR
apply_normal(A::AbstractMatrix, x::AbstractVector) = A * x
apply_adjoint(A::AbstractMatrix, x::AbstractVector) = A' * x
apply_normal((f,fadjoint)::Tuple{Any, Any}, x) = f(x)
apply_adjoint((f,fadjoint)::Tuple{Any, Any}, x) = fadjoint(x)
apply_normal(f, x) = f(x, Val(false))
apply_adjoint(f, x) = f(x, Val(true))

# generalized eigenvalue problem
genapply((A,B)::Tuple{Any, Any}, x) = (apply(A, x), apply(B, x))
genapply(f, x) = f(x)


text:
#Read/write DataMD datasets to/from HDF5
#-------------------------------------------------------------------------------

#==Useful constants
===============================================================================#

#Possible T element types of DataHR{T} & DataRS{T}:
const MAP_STR2TELEM = Dict{String, Type}(
	"DataRS" => DataRS,
	"DataF1" => DataF1,
	"INT64" => Int64,
	"INT32" => Int32,
	"FLOAT64" => Float64,
	"FLOAT32" => Float32,
	"CPLX128" => Complex{Float64},
	"BOOL" => Bool,
)
const MAP_TELEM2STR = Dict{Type, String}(v=>k for (k,v) in MAP_STR2TELEM)

#==Helper functions
===============================================================================#

#Read/write dataset subtype:
#-------------------------------------------------------------------------------
_write_elemtype_attr(grp::HDF5.Group, ::DataHR{T}) where T = HDF5.write_attribute(grp, "ELEMTYPE", MAP_TELEM2STR[T])
_write_elemtype_attr(grp::HDF5.Group, ::DataRS{T}) where T = HDF5.write_attribute(grp, "ELEMTYPE", MAP_TELEM2STR[T])

_write_datatype_attr(grp::HDF5.Group, d::DataF1) = _write_datatype_attr(grp, DataF1)
function _write_datatype_attr(grp::HDF5.Group, d::DataHR)
	_write_datatype_attr(grp, DataHR)
	_write_elemtype_attr(grp, d)
end
function _write_datatype_attr(grp::HDF5.Group, d::DataRS)
	_write_datatype_attr(grp, DataRS)
	_write_elemtype_attr(grp, d)
end

#Read DataMD datatype, as required
function _read_mddatatype_attr(grp::HDF5.Group)
	basetype = _read_datatype_attr(grp)
	if !(basetype <: DataMD)
		path = HDF5.name(grp)
		throw(Meta.ParseError("Expected T<:DataMD, got $basetype:\n$path"))
	elseif (basetype <: DataHR) || (basetype <: DataRS)
		elemtype = MAP_STR2TELEM[HDF5.read_attribute(grp, "ELEMTYPE")]
		return basetype{elemtype}
	end
	return basetype
end

#Read/write PSweep (DataHR):
#-------------------------------------------------------------------------------
function writepsweep(grp::HDF5.Group, sweeps::Vector{PSweep})
	HDF5.write_attribute(grp, "sweepnames", names(sweeps))
	for i in 1:length(sweeps)
		HDF5.write_attribute(grp, "sweep$i", sweeps[i].v)
	end
end

function readpsweep(grp::HDF5.Group)
	sweepnames = HDF5.read_attribute(grp, "sweepnames")
	result = PSweep[]
	for i in 1:length(sweepnames)
		push!(result, PSweep(sweepnames[i], HDF5.read_attribute(grp, "sweep$i")))
	end
	return result
end

#==Main EasyDataHDF5 dataset read/write functions
===============================================================================#

#Read/write DataF1:
#-------------------------------------------------------------------------------
function _write_typed(grp::HDF5.Group, d::DataF1)
	_write_datatype_attr(grp, d)
	grp["x"] = d.x
	grp["y"] = d.y
end

function __read(::Type{DataF1}, grp::HDF5.Group)
	return DataF1(HDF5.read_dataset(grp, "x"), HDF5.read_dataset(grp, "y"))
end

#Read/write DataRS{T<:Number}:
#-------------------------------------------------------------------------------
function _write_typed(grp::HDF5.Group, d::DataRS{T}) where T<:Number
	_write_datatype_attr(grp, d)
	grp["data"] = d.elem
end

function __read(::Type{DataRS{T}}, grp::HDF5.Group) where T<:Number
	sweeps = readpsweep(grp)
	if length(sweeps) != 1
		path = HDF5.name(grp)
		throw(Meta.ParseError("Each DataRS node must have exactly one sweep:\n$path"))
	end
	data = grp["data"]
	return DataRS{T}(sweeps[1], data)
end

#Read/write DataHR{T<:Number}:
#-------------------------------------------------------------------------------
function _write_typed(grp::HDF5.Group, d::DataHR{T}) where T<:Number
	_write_datatype_attr(grp, d)
	writepsweep(grp, d.sweeps)
	grp["data"] = d.elem
end

function __read(::Type{DataHR{T}}, grp::HDF5.Group) where T<:Number
	sweeps = readpsweep(grp)
	data = HDF5.read_dataset(grp, "data")
	return DataHR{T}(sweeps, data)
end

#Read/write DataRS{DataF1}:
#-------------------------------------------------------------------------------
function _write_typed(grp::HDF5.Group, d::DataRS{DataF1})
	_write_datatype_attr(grp, d)
	writepsweep(grp, PSweep[d.sweep])

	for i in 1:length(d.elem)
		_write_typed(grp, "$i", d.elem[i])
	end
end

function __read(::Type{DataRS{DataF1}}, grp::HDF5.Group)
	sweeps = readpsweep(grp)
	if length(sweeps) != 1
		path = HDF5.name(grp)
		throw(Meta.ParseError("Each DataRS node must have exactly one sweep:\n$path"))
	end
	data = DataRS{DataF1}(sweeps[1])

	for i in 1:length(data.elem)
		data.elem[i] = _read_datamd(grp, "$i")
	end

	return data
end

#Read/write DataHR{DataF1}:
#-------------------------------------------------------------------------------
function _write_typed(grp::HDF5.Group, d::DataHR{DataF1})
	_write_datatype_attr(grp, d)
	writepsweep(grp, d.sweeps)

	for inds in subscripts(d)
		subpath = join(inds, "/")
		_write_typed(grp, subpath, d.elem[inds...])
	end
end

function __read(::Type{DataHR{DataF1}}, grp::HDF5.Group)
	sweeps = readpsweep(grp)
	data = DataHR{DataF1}(sweeps)

	for inds in subscripts(data)
		subpath = join(inds, "/")
		data.elem[inds...] = _read_datamd(grp, "$subpath")
	end

	return data
end

#Read/write DataRS{DataRS}:
#-------------------------------------------------------------------------------
function _write_typed(grp::HDF5.Group, d::DataRS{DataRS})
	_write_datatype_attr(grp, d)
	writepsweep(grp, PSweep[d.sweep])

	for i in 1:length(d.elem)
		_write_typed(grp, "$i", d.elem[i])
	end
end

function __read(::Type{DataRS{DataRS}}, grp::HDF5.Group)
	sweeps = readpsweep(grp)
	if length(sweeps) != 1
		throw(Meta.ParseError("Each DataRS node must have exactly one sweep:\n$path"))
	end
	data = DataRS{DataRS}(sweeps[1])

	for i in 1:length(data.elem)
		data.elem[i] = _read_datamd(grp, "$i")
	end

	return data
end

#Catch-all for _write_typed() interface:
#-------------------------------------------------------------------------------
function _write_typed(grp::HDF5.Group, name::String, d::DataMD)
	objgrp = creategrp(grp, name)
	_write_typed(objgrp, d)
end

#Traps for typed _read_typed() function (because we need element type for DataRS & DataHR):
#-------------------------------------------------------------------------------
function __read(::Union{Type{DataHR},Type{DataRS}}, grp::HDF5.Group)
	dtype = _read_mddatatype_attr(grp)
		path = HDF5.name(grp)
	return __read(dtype, grp) #Call more specific version
end

#Read DataMD types (start point):
#(Slightly more direct than calling _read_typed(grp, name) because it knows
#to look for element type as well)
#-------------------------------------------------------------------------------
function _read_datamd(grp::HDF5.Group, name::String)
	dgrp = opengrp(grp, name)
	dtype = _read_mddatatype_attr(dgrp)
	return __read(dtype, dgrp) #Read appropriate data structure
end

#==User-level interface:
===============================================================================#

function Base.write(w::EasyDataWriter, d::DataMD, name::String)
	grp = opengrp(w, hdf5dataroot)
	return _write_typed(grp, name, d)
end
function readdata(r::EasyDataReader, name::String) #Avoids specifying ::Type{DataMD}
	grp = opengrp(r, hdf5dataroot)
	return _read_datamd(grp, name)
end
Base.read(::Type{DataMD}, r::EasyDataReader, name::String) = readdata(r, name)

#Last Line