text:
import SymmetricDominance
using Base.Test

const sd = SymmetricDominance

const gdb = sd.GeneDB()

for _ = 1:2
    sd.insert!(sd.WithNewAllele, gdb, 1, 0)
end
for i = 1:2, _ = 1:2
    sd.insert!(sd.WithNewAllele, gdb, 2, i)
end
for i = 3:6, _ = 1:2
    sd.insert!(sd.WithNewAllele, gdb, 3, i)
end
@test gdb.currentid == 14
@test gdb.currentstate == 14
sd.insert!(sd.WithoutNewAllele, gdb, 4, 1, sd.select(gdb, 1, :state)[1])
@test gdb.currentid == 15
@test gdb.currentstate == 14

@test sd.isidbystate(gdb, 1, 1) == true
@test sd.isidbystate(gdb, 1, 15) == true
@test sd.isidbystate(gdb, 15, 1) == true
@test sd.isidbystate(gdb, 1, 2) != true
@test sd.isidbystate(gdb, 2, 1) != true
@test sd.isidbystate(gdb, 2, 6) != true
@test sd.isidbystate(gdb, 6, 2) != true

@test sd.getancestor(gdb, 1, 1) == 1
@test sd.getancestor(gdb, 2, 1) == 2
@test sd.getancestor(gdb, 3, 1) == 1
@test sd.getancestor(gdb, 4, 1) == 1
@test sd.getancestor(gdb, 5, 1) == 2
@test sd.getancestor(gdb, 6, 1) == 2
@test sd.getancestor(gdb, 7, 1) == 1
@test sd.getancestor(gdb, 8, 1) == 1
@test sd.getancestor(gdb, 9, 1) == 1
@test sd.getancestor(gdb, 10, 1) == 1
@test sd.getancestor(gdb, 11, 1) == 2
@test sd.getancestor(gdb, 12, 1) == 2
@test sd.getancestor(gdb, 13, 1) == 2
@test sd.getancestor(gdb, 14, 1) == 2
@test sd.getancestor(gdb, 15, 1) == 1
@test sd.getancestor(gdb, 3, 2) == 3
@test sd.getancestor(gdb, 4, 2) == 4
@test sd.getancestor(gdb, 5, 2) == 5
@test sd.getancestor(gdb, 6, 2) == 6
@test sd.getancestor(gdb, 7, 2) == 3
@test sd.getancestor(gdb, 8, 2) == 3
@test sd.getancestor(gdb, 9, 2) == 4
@test sd.getancestor(gdb, 10, 2) == 4
@test sd.getancestor(gdb, 11, 2) == 5
@test sd.getancestor(gdb, 12, 2) == 5
@test sd.getancestor(gdb, 13, 2) == 6
@test sd.getancestor(gdb, 14, 2) == 6
@test sd.getancestor(gdb, 15, 2) == 1
@test sd.getancestor(gdb, 7, 3) == 7
@test sd.getancestor(gdb, 8, 3) == 8
@test sd.getancestor(gdb, 9, 3) == 9
@test sd.getancestor(gdb, 10, 3) == 10
@test sd.getancestor(gdb, 11, 3) == 11
@test sd.getancestor(gdb, 12, 3) == 12
@test sd.getancestor(gdb, 13, 3) == 13
@test sd.getancestor(gdb, 14, 3) == 14
@test sd.getancestor(gdb, 15, 3) == 1

const o = sd.Organism(3)
@test sd.nloci(o) == 3

const pop = sd.Population(2, 1)
pop[1, 1, 1] = 7
pop[2, 1, 1] = 8
pop[1, 1, 2] = 9
pop[2, 1, 2] = 10
@test sd.hascoalesced(gdb, pop, 1, 1)
@test !sd.hascoalesced(gdb, pop, 1, 2)

for i = 1:15
    gdb.state[i] += 10
end

sdata3 = Array(Int, 2, 1, 2)
sdata3[1, 1, 1] = 7
sdata3[2, 1, 1] = 8
sdata3[1, 1, 2] = 9
sdata3[2, 1, 2] = 10
@test sd.toarray(gdb, pop, :id) == sdata3
sdata3[1, 1, 1] = 17
sdata3[2, 1, 1] = 18
sdata3[1, 1, 2] = 19
sdata3[2, 1, 2] = 20
@test sd.toarray(gdb, pop, :state) == sdata3
sdata3[1, 1, 1] = 3
sdata3[2, 1, 1] = 3
sdata3[1, 1, 2] = 4
sdata3[2, 1, 2] = 4
@test sd.toarray(gdb, pop, :parent) == sdata3
sdata3[1, 1, 1] = 3
sdata3[2, 1, 1] = 3
sdata3[1, 1, 2] = 3
sdata3[2, 1, 2] = 3
@test sd.toarray(gdb, pop, :epoch) == sdata3

const acs, gcs, hcs = sd.counts(gdb, pop)
@test acs == [Dict{Int,Int}(17=>1, 18=>1, 19=>1, 20=>1)]
@test gcs == [Dict{Tuple{Int,Int}, Int}((17, 19)=>1, (18, 20)=>1)]
@test hcs == Dict{Tuple{Int}, Int}((17,)=>1, (18,)=>1, (19,)=>1, (20,)=>1)

const afs, gfs, hfs = sd.spectra(gdb, pop)
@test afs == [Dict{Int,Int}(1=>4)]
@test gfs == [Dict{Int,Int}(1=>2)]
@test hfs == Dict{Int,Int}(1=>4)

@test sd.history(gdb, 1, 1) == [1]
@test sd.history(gdb, 2, 1) == [2]
@test sd.history(gdb, 3, 1) == [1, 3]
@test sd.history(gdb, 4, 1) == [1, 4]
@test sd.history(gdb, 5, 1) == [2, 5]
@test sd.history(gdb, 6, 1) == [2, 6]
@test sd.history(gdb, 7, 1) == [1, 3, 7]
@test sd.history(gdb, 8, 1) == [1, 3, 8]
@test sd.history(gdb, 9, 1) == [1, 4 ,9]
@test sd.history(gdb, 10, 1) == [1, 4, 10]
@test sd.history(gdb, 11, 1) == [2, 5, 11]
@test sd.history(gdb, 12, 1) == [2, 5, 12]
@test sd.history(gdb, 13, 1) == [2, 6, 13]
@test sd.history(gdb, 14, 1) == [2, 6, 14]
@test sd.history(gdb, 15, 1) == [1, 15]
@test sd.history(gdb, 3, 2) == [3]
@test sd.history(gdb, 4, 2) == [4]
@test sd.history(gdb, 5, 2) == [5]
@test sd.history(gdb, 6, 2) == [6]
@test sd.history(gdb, 7, 2) == [3, 7]
@test sd.history(gdb, 8, 2) == [3, 8]
@test sd.history(gdb, 9, 2) == [4 ,9]
@test sd.history(gdb, 10, 2) == [4, 10]
@test sd.history(gdb, 11, 2) == [5, 11]
@test sd.history(gdb, 12, 2) == [5, 12]
@test sd.history(gdb, 13, 2) == [6, 13]
@test sd.history(gdb, 14, 2) == [6, 14]
@test sd.history(gdb, 15, 2) == [1, 15]
@test sd.history(gdb, 7, 3) == [7]
@test sd.history(gdb, 8, 3) == [8]
@test sd.history(gdb, 9, 3) == [9]
@test sd.history(gdb, 10, 3) == [10]
@test sd.history(gdb, 11, 3) == [11]
@test sd.history(gdb, 12, 3) == [12]
@test sd.history(gdb, 13, 3) == [13]
@test sd.history(gdb, 14, 3) == [14]
@test sd.history(gdb, 15, 3) == [1, 15]

ds = sd.distances(gdb, pop, 1)
@test ds == [1 17 18 2;
             1 17 19 4;
             1 17 20 4;
             1 18 19 4;
             1 18 20 4;
             1 19 20 2]
ds = sd.distances(gdb, pop, 2)
@test ds == [1 17 18 2;
             1 17 19 -1;
             1 17 20 -1;
             1 18 19 -1;
             1 18 20 -1;
             1 19 20 2]
ds = sd.distances(gdb, pop, 3)
@test ds == [1 17 18 -1;
             1 17 19 -1;
             1 17 20 -1;
             1 18 19 -1;
             1 18 20 -1;
             1 19 20 -1]

@test sd.nsegsites(gdb, pop, 1) == [6]
@test sd.nsegsites(gdb, pop, 2) == [-1]
@test sd.nsegsites(gdb, pop, 3) == [-1]

@test sd.tmrca(gdb, pop) == [1]


text:
module TestReduce

using FLoops
using MicroCollections
using Test
using Transducers: SplitBy

function test_sum()
    @floop for (x, y) in zip(1:3, 1:2:6)
        a = x + y
        b = x - y
        @reduce(s += a, t += b)
    end
    @test (s, t) === (15, -3)
end

function test_sum_with_init()
    @floop for (x, y) in zip(1:3, 1:2:6)
        a = x + y
        b = x - y
        @reduce(s = 0im + a, t = 0im + b)
    end
    @test (s, t) === (15 + 0im, -3 + 0im)
end

function test_sum_nested_loop()
    @floop for x in 1:3, y in 1:2:6
        a = x + y
        b = x - y
        @reduce(s += a, t += b)
    end
    @test (s, t) === (45, -9)
end

function test_findmax()
    @floop for (i, v) in pairs([0, 1, 3, 2])
        y = 2v
        @reduce() do (ymax = -Inf; y), (imax = 0; i)
            if isless(ymax, y)
                ymax = y
                imax = i
            end
        end
    end
    @test (ymax, imax) == (6, 3)
end

function test_findminmax()
    @floop for (i, v) in pairs([0, 1, 3, 2])
        y = 2v
        @reduce() do (ymax = -Inf; y), (imax = 0; i)
            if isless(ymax, y)
                ymax = y
                imax = i
            end
        end
        @reduce() do (ymin = Inf; y), (imin = 0; i)
            if isless(y, ymin)
                ymin = y
                imin = i
            end
        end
    end
    @test (ymax, imax) == (6, 3)
    @test (ymin, imin) == (0, 1)
end

function sum_arrays_broadcast(arrays, ex = nothing)
    @floop ex for x in arrays
        @reduce s .+= x
    end
    try
        return s
    catch err
        return err
    end
end

function test_simple_broadcast()
    @test sum_arrays_broadcast([[1], [2], [3]]) == [sum(1:3)]
    @test sum_arrays_broadcast([[1], [2], [3]], SequentialEx()) == [sum(1:3)]
    @test sum_arrays_broadcast([]) == UndefVarError(:s)
    @test sum_arrays_broadcast([], SequentialEx()) == UndefVarError(:s)
end

function fused_broadcast(xs)
    ys = nothing
    vs = 1:2:11
    @floop for x in xs
        @reduce ys .+= vs .== x
    end
    return ys
end

function test_fused_broadcast()
    function desired(n)
        m = cld(n, 2)
        return [ones(Int, m); zeros(Int, 6 - m)]
    end
    @test fused_broadcast(1:0) === nothing
    @testset for n in 1:11
        @test fused_broadcast(1:n) == desired(n)
    end
end

⊕(a, b) = a .+ b

function mixed_broadcasts(xs)
    @floop for x in xs
        @reduce(
            a .+= x,              # symbol
            b .+= isodd.(x),      # dot call
            c .+= .√(x),          # unary dot op
            d .+= x .- 1,         # binary dot op
            e .+= prod(x),        # normal call
            f .+= x ⊕ 1,          # normal binary op
            g .+= -x,             # normal unary op
        )
    end
    try
        return (a = a, b = b, c = c, d = d, e = e, f = f, g = g)
    catch err
        return err
    end
end

function mixed_broadcasts2(xs)
    @floop for x in xs
        @reduce(
            a .= 0 .+ x,          # symbol
            b .= 0 .+ isodd.(x),  # dot call
            c .= 0 .+ .√(x),      # unary dot op
            d .= 0 .+ (x .- 1),   # binary dot op
            e .= 0 .+ prod(x),    # normal call
            f .= 0 .+ (x ⊕ 1),    # normal binary op
            g .= 0 .+ -x,         # normal unary op
        )
    end
    try
        return (a = a, b = b, c = c, d = d, e = e, f = f, g = g)
    catch err
        return err
    end
end

function test_mixed_broadcasts()
    @test mixed_broadcasts([[4], [9]]) ==
          mixed_broadcasts2([[4], [9]]) ==
          (a = [13], b = [1], c = [5], d = [11], e = 13, f = [15], g = [-13])
    @test mixed_broadcasts([[4]]) ==
          mixed_broadcasts2([[4]]) ==
          (a = [4], b = [0], c = [2], d = [3], e = 4, f = [5], g = [-4])
    @test mixed_broadcasts(1:0) == UndefVarError(:a)
    @test mixed_broadcasts2(1:0) == (a = 0, b = 0, c = 0, d = 0, e = 0, f = 0, g = 0)
end

function sum_onehot(indices, ex = nothing)
    l, h = extrema(indices)
    n = h - l + 1
    @floop ex for i in indices
        @inbounds @reduce h .+= OneHotVector(i - l + 1 => 1, n)
    end
    return h
end

function sum_prealloc_onehot(indices, ex = nothing)
    l, h = extrema(indices)
    n = h - l + 1
    @floop ex for i in indices
        @reduce h .= zeros(Int, n) .+ OneHotVector(i - l + 1 => 1, n)
    end
    return h
end

function test_onehot()
    @testset "$(repr(ex))" for ex in [SequentialEx(), nothing]
        @test sum_onehot(1:3, ex) == [1, 1, 1]
        @test sum_onehot([1, 2, 4, 1], ex) == [2, 1, 0, 1]
        @test sum_prealloc_onehot(1:3, ex) == [1, 1, 1]
        @test sum_prealloc_onehot([1, 2, 4, 1], ex) == [2, 1, 0, 1]
    end
end

function test_break()
    @floop for x in 1:10
        @reduce(s += x)
        x == 3 && break
    end
    @test s == 6
end

function test_findmax_with_filtering()
    @floop for (i, v) in pairs([0, 1, 4, 3, 1])
        if isodd(v)
            @reduce() do (vmax = -Inf; v), (imax = 0; i)
                if isless(vmax, v)
                    vmax = v
                    imax = i
                end
            end
        end
    end
    @test (vmax, imax) == (3, 4)
end

function test_nonunique_arguments()
    @testset "`@reduce(a += x, b *= x)`" begin
        xs = 1:10
        @floop for x in xs
            @reduce(a += x, b *= x)
        end
        @test (sum(xs), prod(xs)) == (a, b)
    end

    @testset "`@reduce(a = 0 + x, b = 1 * x)`" begin
        xs = 1:10
        @floop for x in xs
            @reduce(a = 0 + x, b = 1 * x)
        end
        @test (sum(xs), prod(xs)) == (a, b)
    end

    @testset "`@reduce a = 0 + x; @reduce b = 1 * x`" begin
        xs = 1:10
        @floop for x in xs
            @reduce a = 0 + x
            @reduce b = 1 * x
        end
        @test (sum(xs), prod(xs)) == (a, b)
    end
end

function floop_with_init_binops(xs, ex = nothing)
    @floop ex for x in xs
        @reduce(s = 0 + x, c = 0 + 1)
    end
    return f() = (s, c)
end

function two_floops(xs, ex = nothing)
    @floop ex for x in xs
        @reduce(s = 0 + x, c = 0 + 1)
    end
    @floop ex for x in xs
        @reduce(v = 0 + (x - (s / c))^2)
    end
    return v / c
end

function test_empty()
    @test !FLoops.has_boxed_variables(floop_with_init_binops(1:10))
    @test two_floops(1:8) == 5.25
end

function sum_arrays(arrays, ex = nothing)
    @floop ex for x in arrays
        @reduce() do (s = zero(x); x)
            s .+= x
        end
    end
    return s
end

function test_at_reduce_init_scope()
    @test sum_arrays([[1], [2], [3]]) == [sum(1:3)]
    @test sum_arrays([[1], [2], [3]], SequentialEx()) == [sum(1:3)]
    @testset "let" begin
        local y = 0
        @floop for x in 1:10
            let y = 123
            end
            @reduce(s = y + x)
        end
        @test s == sum(1:10)
    end
end

function sum_arrays2(arrays, ex = nothing)
    @floop ex for x in arrays
        @reduce(s = zero(x) + x)
    end
    try
        return s
    catch err
        return err
    end
end

function test_init_in_loop()
    @test sum_arrays2([[1], [2], [3]]) == [sum(1:3)]
    @test sum_arrays2([[1], [2], [3]], SequentialEx()) == [sum(1:3)]
    @test sum_arrays2([]) == UndefVarError(:s)
    @test sum_arrays2([], SequentialEx()) == UndefVarError(:s)
end

function maximum_partition_length(f, xs, ex = nothing)
    @floop ex for chunk in SplitBy(f)(xs)
        @reduce m = max(typemin(Int), length(chunk))
    end
    return m
end

function test_SplitBy()
    @test maximum_partition_length(isodd, 1:10) == 1
    @test maximum_partition_length(isodd, 1:10, SequentialEx()) == 1
    @test maximum_partition_length(==(7), 1:10) == 6
    @test maximum_partition_length(==(7), 1:10, SequentialEx()) == 6
    xs = rand(0:9, 1000)
    @test maximum_partition_length(iszero, xs) ==
          maximum_partition_length(iszero, xs, SequentialEx())
end

function sum_halved_arrays(arrays, ex = nothing)
    @floop ex for x in arrays
        @init y = zero(x)
        y .= x .÷ 2
        @reduce(s = 0 + sum(y))
    end
    return s
end

function two_inits(arrays, ex = nothing)
    @floop ex for x in arrays
        @init begin
            y = zero(x)
            z = similar(y)
        end
        y .= x .÷ 2
        z .= 2 .* y
        r = sum(a * b for (a, b) in zip(y, z))
        @reduce(s = 0 + r)
    end
    return s
end

function two_inits2(arrays, ex = nothing)
    @floop ex for x in arrays
        @init y = zero(x)
        @init z = similar(y)
        y .= x .÷ 2
        z .= 2 .* y
        r = sum(a * b for (a, b) in zip(y, z))
        @reduce(s = 0 + r)
    end
    return s
end

function test_at_init()
    arrays = [[1, 2], [3, 4], [5, 6], [7, 8]]
    desired = sum(sum(x .÷ 2) for x in arrays)
    @test sum_halved_arrays(arrays) == desired
    @test sum_halved_arrays(arrays, SequentialEx()) == desired

    desired = sum((c ÷ 2) * (2 * (c ÷ 2)) for x in arrays for c in x)
    @test two_inits(arrays) == desired
    @test two_inits(arrays, SequentialEx()) == desired
    @test two_inits2(arrays) == desired
    @test two_inits2(arrays, SequentialEx()) == desired
end

function test_just_at_init()
    n = 10
    dest = zeros(Int, n)
    @floop for (i, x) in pairs(1:n)
        @init y = zeros(Int, 3)
        y .= (x, 2x, 3x)
        dest[i] = y[1]
    end
    @test dest == 1:n
end

function probe_init(xs, ex = nothing)
    ninit = Threads.Atomic{Int}(0)
    @floop ex for x in xs
        @init y = begin
            Threads.atomic_add!(ninit, 1)
            nothing
        end
    end
    return ninit[]
end

function test_at_init_called_once()
    @testset "default" begin
        @test probe_init(1:(Threads.nthreads()*100)) == Threads.nthreads()
    end
    @testset "SequentialEx" begin
        @test probe_init(1:(Threads.nthreads()*100), SequentialEx()) == 1
    end
end

function init_is_private(xs, ex = nothing)
    @floop ex for x in xs
        @init p = []
        push!(p, x)
        y = pop!(p)
        @reduce s += y
    end
    return (sum = s, isdefined_p = @isdefined(p))
end

function test_init_is_private()
    @test init_is_private(1:10) == (sum = sum(1:10), isdefined_p = false)
    @test init_is_private(1:10, SequentialEx()) == (sum = sum(1:10), isdefined_p = false)
end

function test_unprocessed_at_reduce()
    err = try
        @reduce(s += y, p *= y)
        nothing
    catch err
        err
    end
    @test err isa FLoops.ReduceOpSpec
    @test occursin("used outside `@floop`", sprint(showerror, err))
end

function test_unprocessed_at_init()
    err = try
        @init x = 0
        nothing
    catch err
        err
    end
    @test err isa FLoops.InitSpec
    @test occursin("used outside `@floop`", sprint(showerror, err))
end

function test_invalid_at_init()
    @testset "toplevel" begin
        err = try
            @eval @init(a)
            nothing
        catch err
            err
        end
        @test err isa Exception
        @test occursin("requires an assignment", sprint(showerror, err))
    end
    @testset "non assignment (+=)" begin
        err = try
            @eval @init a += 1
            nothing
        catch err
            err
        end
        @test err isa Exception
        @test occursin("requires an assignment", sprint(showerror, err))
    end
    @testset "non assignment (constant)" begin
        err = try
            @eval @init begin
                a = 1
                1
            end
            nothing
        catch err
            err
        end
        @test err isa Exception
        @test occursin("requires an assignment", sprint(showerror, err))
    end
end

function test_duplicated_accumulators()
    err = try
        @eval @macroexpand @floop for x in xs
            @reduce(y += x, y += x)
        end
        nothing
    catch err
        err
    end
    @test err isa Exception
    @test occursin("`y` used 2 times", sprint(showerror, err))
end

function test_heterogeneous_dot_update_syntax()
    err = try
        @eval @macroexpand @floop for x in xs
            @reduce(a .+= x, b .= 0 .+ x)
        end
        nothing
    catch err
        err
    end
    @test err isa Exception
    @test occursin("e.g., don't mix `a .+= x` and `a .= 0 .+ x`", sprint(showerror, err))
end

function test_non_dot_call()
    err = try
        @eval @macroexpand @floop for x in xs
            @reduce a .= 0 + x
        end
        nothing
    catch err
        err
    end
    @test err isa Exception
    @test occursin("requires a binary dot call", sprint(showerror, err))
end

end  # module

ext: jl


text:
test_volume_results = [0, 1, 1, 8, 1, 1, 1, 1, 2, 1, 1, 1, 3, 2, 1, 9, 3, 2, 2, 6]

black_surface = zeros(Int, 5, 5)
white_surface = ones(Int, 5, 5)

separated_voxels = [0 0 0 0 0; 0 1 0 1 0; 0 0 0 0 0; 0 1 0 1 0; 0 0 0 0 0]
checkerboard_1 = [1 0 1 0 1; 0 1 0 1 0; 1 0 1 0 1; 0 1 0 1 0; 1 0 1 0 1]
checkerboard_2 = [0 1 0 1 0; 1 0 1 0 1; 0 1 0 1 0; 1 0 1 0 1; 0 1 0 1 0]
hollow_structure = [0 0 0 0 0; 0 1 1 1 0; 0 1 0 1 0; 0 1 1 1 0; 0 0 0 0 0]
cutoff_hollow = [0 0 0 0 0; 0 1 1 1 0; 0 1 0 1 0; 0 1 0 1 0; 0 1 0 1 0]
two_blocks = [0 0 0 0 0; 1 1 1 1 0; 0 0 0 0 0; 0 1 1 1 1; 0 0 0 0 0]
single_hollow_voxel = [1 0 1 0 1; 1 1 1 1 1; 1 1 0 1 1; 1 1 1 1 1; 0 1 1 1 0]
enclosed_hollow = [1 0 1 0 1; 1 1 1 1 1; 1 0 0 0 1; 1 1 1 1 1; 0 1 1 1 0]
multiple_single_hollow_voxel = [1 0 1 0 1; 1 1 1 1 1; 1 0 1 0 1; 1 1 1 1 1; 0 1 0 1 0]
border_voxels = [1 1 1 1 1; 1 0 0 0 1; 1 0 0 0 1; 1 0 0 0 1; 1 1 1 1 1]

# Black voxels
test_volume_1 = zeros(Int, 5, 5, 5)

# White voxels
test_volume_2 = ones(Int, 5, 5, 5)

# Single voxel
_test_volume_3 = zeros(Int, 5, 5, 5)
_test_volume_3[3, 3, 3] = 1
test_volume_3 = _test_volume_3

# Multiple single voxels
test_volume_4 = cat(black_surface, separated_voxels, black_surface, separated_voxels, black_surface; dims = 3)

# Checkerboard 1
test_volume_5 = cat(checkerboard_1, checkerboard_2, checkerboard_1, checkerboard_2, checkerboard_1; dims = 3)

# Checkerboard 2
test_volume_6 = cat(checkerboard_2, checkerboard_1, checkerboard_2, checkerboard_1, checkerboard_2; dims = 3)

# Hollow structure
test_volume_7 = cat(black_surface, hollow_structure, hollow_structure, hollow_structure, black_surface; dims = 3)

# Cut-off hollow structure
test_volume_8 = cat(black_surface, cutoff_hollow, cutoff_hollow, cutoff_hollow, cutoff_hollow; dims = 3)

# Two blocks
test_volume_9 = cat(black_surface, two_blocks, two_blocks, two_blocks, black_surface; dims = 3)

# Single hollow voxel
test_volume_10 = cat(single_hollow_voxel, single_hollow_voxel, single_hollow_voxel, single_hollow_voxel, single_hollow_voxel; dims = 3)

# Enclosed hollow structure
test_volume_11 = cat(enclosed_hollow, enclosed_hollow, enclosed_hollow, black_surface, black_surface; dims = 3)

# Multiple single hollow voxels
test_volume_12 = cat(black_surface, black_surface, multiple_single_hollow_voxel, black_surface, black_surface; dims = 3)

# Border voxels
test_volume_13 = cat(border_voxels, black_surface, border_voxels, black_surface, border_voxels; dims = 3)

# Cube in cube
_test_volume_14 = cat(white_surface, border_voxels, border_voxels, border_voxels, white_surface; dims = 3)
_test_volume_14[3, 3, 3] = 1
test_volume_14 = _test_volume_14

# Unfixed hollow voxel sizes
test_volume_15 = cat(border_voxels, hollow_structure, hollow_structure, border_voxels; dims = 3)

# Cubes in cube
_white_surface_7 = ones(Int, 7, 7)
_border_voxels_7 = [1 1 1 1 1 1 1; 1 0 0 0 0 0 1; 1 0 0 0 0 0 1; 1 0 0 0 0 0 1; 1 0 0 0 0 0 1; 1 0 0 0 0 0 1; 1 1 1 1 1 1 1]
_cube_slices = [1 1 1 1 1 1 1; 1 0 0 0 0 0 1; 1 0 1 0 1 0 1; 1 0 0 0 0 0 1; 1 0 1 0 1 0 1; 1 0 0 0 0 0 1; 1 1 1 1 1 1 1]
test_volume_16 = cat(_white_surface_7, _border_voxels_7, _cube_slices, _border_voxels_7, _cube_slices, _border_voxels_7, _white_surface_7; dims = 3)

# Cube in cube in cube
_white_surface_9 = ones(Int, 9, 9)
_border_voxels_9 = [1 1 1 1 1 1 1 1 1; 1 0 0 0 0 0 0 0 1; 1 0 0 0 0 0 0 0 1; 1 0 0 0 0 0 0 0 1; 1 0 0 0 0 0 0 0 1; 1 0 0 0 0 0 0 0 1; 1 0 0 0 0 0 0 0 1; 1 0 0 0 0 0 0 0 1; 1 1 1 1 1 1 1 1 1]
_test_volume_17 = cat(_white_surface_9, _border_voxels_9, _border_voxels_9, _border_voxels_9, _border_voxels_9, _border_voxels_9, _border_voxels_9, _border_voxels_9, _white_surface_9; dims = 3)
_test_volume_17[3:7, 3:7, 3:7] = test_volume_14
test_volume_17 = _test_volume_17

# Cuboid 1 (4 x 5 x 5)
_test_volume_18_1 = [0 1 1 1 0; 1 1 0 1 1; 1 0 0 1 1; 1 1 1 1 0]
_test_volume_18_2 = [0 1 1 1 0; 1 1 1 1 1; 1 1 1 1 1; 1 1 1 0 0]
_test_volume_18_3 = [0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0]
_test_volume_18_4 = [0 1 1 1 0; 1 1 0 1 1; 1 1 0 1 1; 1 1 0 0 0]
_test_volume_18_5 = [0 1 1 1 0; 1 1 0 1 1; 1 1 0 1 1; 1 1 0 0 0]
test_volume_18 = cat(_test_volume_18_1, _test_volume_18_2, _test_volume_18_3, _test_volume_18_4, _test_volume_18_5; dims = 3)

# Cuboid 2 (7 x 5 x 3)
_test_volume_19_1 = [1 1 1 1 1; 1 0 1 0 1; 0 0 1 0 0; 1 0 1 0 1; 1 1 1 1 1; 0 0 0 0 0; 0 0 0 0 0]
_test_volume_19_2 = [0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0]
_test_volume_19_3 = [1 0 0 0 0; 0 1 0 0 0; 0 0 1 0 0; 0 0 0 1 0; 0 0 0 0 1; 0 0 0 1 0; 0 0 1 0 0]
test_volume_19 = cat(_test_volume_19_1, _test_volume_19_2, _test_volume_19_3; dims = 3)

# Cuboid 3 (4 x 6 x 2)
_test_volume_20_1 = [1 1 0 0 0 0 0; 0 0 0 0 0 1 0; 0 0 0 0 0 0 0; 0 0 1 0 0 0 0]
_test_volume_20_2 = [0 0 0 1 0 0 0; 0 0 0 0 0 0 0; 1 0 0 0 0 0 0; 0 0 0 0 1 1 0]
test_volume_20 = cat(_test_volume_20_1, _test_volume_20_2; dims = 3)


text:
#=
Created on Saturday 26 December 2020
Last update: -

@author: Bram De Jaegher
bram.de.jaegher@gmail.com
=#

# --- Types --- #

abstract type AbstractQuestion end
abstract type AbstractDifficulty end
abstract type AbstractQuestionBlock end
struct NoDiff <: AbstractDifficulty end
struct Easy <: AbstractDifficulty end
struct Intermediate <: AbstractDifficulty end
struct Hard <: AbstractDifficulty end

"""
Defines a mandatory question.

# Arguments

- description (md_str): a markdown string to be displayed above the validators
- validators (Array{Bool}): an array of booleans with the tests the answer to the question should solve.
- status (md_str): a markdown string used to change/display the question state (correct, missing, incorrect).
	the default value is missing and should probably not be changed.
"""
mutable struct Question <: AbstractQuestion
	description::Markdown.MD
	validators::Any
	status::Markdown.MD

	Question(;description=md"",
						validators=[],
						status=md"") = new(description, validators, status)
end

"""
Defines an optional question.

# Arguments

- `description` (md_str): a markdown string to be displayed above the validators
- `validators` (Array{Bool}): an array of booleans with the tests the answer to the question should solve.
- `status` (md_str): a markdown string used to change/display the question state (correct, missing, incorrect).
	the default value is missing and should probably not be changed.

Difficulty types
------
- `T ∈ [NoDiff, Easy, Intermediate, Hard]`
- At this point purely easthetic since the difficulty is only used as display
"""
mutable struct QuestionOptional{T<:AbstractDifficulty}  <: AbstractQuestion
	description::Markdown.MD
	validators::Any
	status::Markdown.MD
	difficulty::T

	QuestionOptional{T}(;description=md"",
						validators=[],
						status=md"") where {T<:AbstractDifficulty} = new{T}(description, validators, status, T())
end

"""
Defines a Question block.

# Arguments
- `title` (md_str): a markdown string to be displayed as title
- `description` (md_str): the general descriptions, this is diplayed directly below the title.
- `hints` (Array{md_str}): an array of markdown strings for hint admonitions but this can me any markdown and will be displayed in the \"hints\" section.
- `questions` (Array{AbstractQuestion}): an array of questions. Currently exactly one mandatory question is expected but 0-∞ optional questions can be defined.

# Examples
```julia
q₁ = Question(;
	description=md\"\"\"
	Complete the function `myclamp(x)` that clamps a number `x` between 0 and 1.

	Open assignments always return `missing`.
	\"\"\",
	validators= @safe[myclamp(-1)==0, myclamp(0.3)==0.3, myclamp(1.1)==1.0]
)

q₂ = QuestionOptional{Easy}(;
	description=md\"\"\"
	Try to make the clamping also work for arrays.
	\"\"\",
	validators= @safe[myclamp([2.0, 0.3])==[1.0, 0.3]]
)

q₃ = QuestionOptional{Intermediate}(;
	description=md\"\"\"
	This is an intermediate question. Surely you can complete this
	\"\"\",
	validators= @safe[true]
)

q₄ = QuestionOptional{Hard}(;
	description=md\"\"\"
	I admit, this one is definitely harder
	\"\"\",
	validators= @safe[false]
)

 qb = QuestionBlock(;
title=md"### Question 1.0: What a crazy exercise",
description=md\"\"\"
	Some additional general kind off description and all.
	Anything markdowny. Just make sure to use the triple accolades.

	\"\"\",
questions = [q₁, q₂, q₃, q₄],
hints=[	hint(md"Have you tried this?"),
		hint(md"Have you tried switching it on and off again?")]
);
```
"""
mutable struct QuestionBlock <: AbstractQuestionBlock
	title::Markdown.MD
	description::Markdown.MD
	hints::Array{Markdown.MD}
	questions::Array{T} where {T<:AbstractQuestion}

	QuestionBlock(;title=md"",
									description=md"",
									hints = Markdown.MD[],
									questions = [Question()]) = new(title, description, hints, questions)
end

# --- Rendering --- #
Base.show(io::IO, ::MIME"text/html", q::AbstractQuestionBlock) = print(io::IO, tohtml(q))

function tohtml(q::QuestionBlock)

	hint_string = ""
	if length(q.hints) > 0
		hint_string = "

Hints:
"
		for hint in q.hints
			hint_string *= "
" * html(hint) * "
"
		end
	end

	N_mandatory = sum(isa.(q.questions,Question))
	state_string = ""
	for index in 1:N_mandatory
		if q.questions[index].description !== ""
			state_string *= "
 $(html(q.questions[index].description))
"
		end
		state_string *= "
 $(html(q.questions[index].status))
"
	end

	opt_state_string = ""
	if length(q.questions) > 1
		for opt_question in q.questions[N_mandatory+1:end]
			if opt_question.difficulty !== ""
			opt_state_string *= "
  Optional ($(split(string(typeof(opt_question.difficulty)), ".")[2])):
"
			else
			opt_state_string *= "
  Optional:
"
			end
			if opt_question.description !== ""
				opt_state_string *= "
 $(html(opt_question.description))
"
			end
			opt_state_string *= "
 $(html(opt_question.status))
"
		end
	end

	out = """

			$(html(q.title))

 $(html(q.description))

			$state_string
			$opt_state_string
			$hint_string

		$question_css
	"""
	return out
end

# --- Macro(s) --- #
"""
The @safe macro is a hidden try-catch statement to avoid the Markdown admonitions to crash when
the user introduces an error in the exercise functions.

Usage:
-----
@safe test_func(input) == correct_output
@safe[test_func1(input) == correct_output1, test_func2(input) == correct_output2]

"""
macro safe(ex)
	safe_ex = quote
		try $(esc(ex))
		catch e
			false
		end
	end
	return safe_ex
end

# --- Validation --- #
function check_answer(q::AbstractQuestion)
	validators = q.validators
	all_valid = all(validators)
	some_valid = any(validators)

	if length(validators) < 1 && return md""
	elseif ismissing(all_valid)
		status = still_missing()
	elseif some_valid && !all_valid
		status = keep_working()
	elseif !all_valid
		status = keep_working()
	elseif all_valid
		status = correct()
	end
	return status
end

function validate(q::QuestionBlock)
	for (index, question) in enumerate(q.questions)
		q.questions[index].status = check_answer(question)
	end
	return q
end


text:
mutable struct ConicForm{T, AT, VT, C} <: MOI.ModelLike
    num_rows::Vector{Int}
    dimension::Dict{Int, Int}
    sense::MOI.OptimizationSense
    objective_constant::T # The objective
    A::Union{Nothing, AT} # The constraints
    b::VT          # `b - Ax in cones`
    c::VT          # `sense c'x + objective_constant`
    cone_types::C
    cone_types_dict::Dict{DataType, Int}

    function ConicForm{T, AT, VT}(cone_types) where {T, AT, VT}
        model = new{T, AT, VT, typeof(cone_types)}()
        model.cone_types = cone_types
        model.cone_types_dict = Dict{DataType, Int}(
            s => i for (i, s) in enumerate(cone_types)
        )
        model.num_rows = zeros(Int, length(cone_types))
        model.dimension = Dict{Int, Int}()
        model.A = nothing
        return model
    end
end

_set_type(::MOI.ConstraintIndex{F,S}) where {F,S} = S
cons_offset(conic::MOI.Zeros) = conic.dimension
cons_offset(conic::MOI.Nonnegatives) = conic.dimension
cons_offset(conic::MOI.SecondOrderCone) = conic.dimension
cons_offset(conic::MOI.PositiveSemidefiniteConeTriangle) = Int64((conic.side_dimension*(conic.side_dimension+1))/2)

function get_conic_form(::Type{T}, model::M, con_idx) where {T, M <: MOI.AbstractOptimizer}
    # extract cones
    cones = _set_type.(con_idx)
    cones = unique(cones)

    conic = ConicForm{T, SparseMatrixCSRtoCSC{Int64}, Array{T, 1}}(Tuple(cones))

    idxmap = MOI.copy_to(conic, model)

    # fix optimization sense
    if MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
        conic.c = -conic.c
    end

    return conic
end

MOI.is_empty(model::ConicForm) = model.A === nothing
function MOI.empty!(model::ConicForm{T}) where T
    empty!(model.dimension)
    fill!(model.num_rows, 0)
    model.A = nothing
    model.sense = MOI.FEASIBILITY_SENSE
    model.objective_constant = zero(T)
end

function MOI.supports_constraint(
    model::ConicForm,
    ::Type{MOI.VectorAffineFunction{Float64}},
    ::Type{S}) where S <: MOI.AbstractVectorSet
    return haskey(model.cone_types_dict, S)
end

function _allocate_variables(model::ConicForm{T, AT, VT}, vis_src, idxmap) where {T, AT, VT}
    model.A = AT(length(vis_src))
    for (i, vi) in enumerate(vis_src)
        idxmap[vi] = MOI.VariableIndex(i)
    end
    return
end

function rows(model::ConicForm, ci::CI{MOI.VectorAffineFunction{Float64}})
    return ci.value .+ (1:model.dimension[ci.value])
end

function MOI.set(::ConicForm, ::MOI.VariablePrimalStart,
                 ::MOI.VariableIndex, ::Nothing)
end
function MOI.set(model::ConicForm, ::MOI.VariablePrimalStart,
                 vi::MOI.VariableIndex, value::Float64)
    model.primal[vi.value] = value
end
function MOI.set(::ConicForm, ::MOI.ConstraintPrimalStart,
                 ::MOI.ConstraintIndex, ::Nothing)
end
function MOI.set(model::ConicForm, ::MOI.ConstraintPrimalStart,
                 ci::MOI.ConstraintIndex, value)
    offset = constroffset(model, ci)
    model.slack[rows(model, ci)] .= value
end
function MOI.set(::ConicForm, ::MOI.ConstraintDualStart,
                  ::MOI.ConstraintIndex, ::Nothing)
end
function MOI.set(model::ConicForm, ::MOI.ConstraintDualStart,
                  ci::MOI.ConstraintIndex, value)
    offset = constroffset(model, ci)
    model.dual[rows(model, ci)] .= value
end
function MOI.set(model::ConicForm, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    model.sense = sense
end
variable_index_value(t::MOI.ScalarAffineTerm) = t.variable_index.value
variable_index_value(t::MOI.VectorAffineTerm) = variable_index_value(t.scalar_term)
function MOI.set(model::ConicForm, ::MOI.ObjectiveFunction,
                 f::MOI.ScalarAffineFunction{Float64})
    c = Vector(sparsevec(variable_index_value.(f.terms), MOI.coefficient.(f.terms),
                         model.A.n))
    model.objective_constant = f.constant
    model.c = c
    return nothing
end

function _allocate_constraint(model::ConicForm, src, indexmap, cone_id, ci)
    # TODO use `CanonicalConstraintFunction`
    func = MOI.get(src, MOI.ConstraintFunction(), ci)
    func = MOIU.is_canonical(func) ? func : MOI.Utilities.canonical(func)
    allocate_terms(model.A, indexmap, func)
    offset = model.num_rows[cone_id]
    model.num_rows[cone_id] = offset + MOI.output_dimension(func)
    return ci, offset, func
end

function _allocate_constraints(model::ConicForm, src, indexmap, cone_id, ::Type{S}) where S
    cis = MOI.get(src, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, S}())
    return map(cis) do ci
        _allocate_constraint(model, src, indexmap, cone_id, ci)
    end
end

function _load_variables(model::ConicForm, nvars::Integer)
    m = sum(model.num_rows)
    model.A.m = m
    model.b = zeros(m)
    model.c = zeros(model.A.n)
    allocate_nonzeros(model.A)
end

function _load_constraints(model::ConicForm, src, indexmap, cone_offset, i, cache)
    for (ci_src, offset_in_cone, func) in cache
        offset = cone_offset + offset_in_cone
        set = MOI.get(src, MOI.ConstraintSet(), ci_src)
        load_terms(model.A, indexmap, func, offset)
        copyto!(model.b, offset + 1, func.constants)
        model.dimension[offset] = MOI.output_dimension(func)
        indexmap[ci_src] = typeof(ci_src)(offset)
    end
end

function MOI.copy_to(dest::ConicForm, src::MOI.ModelLike; copy_names::Bool=true)
    MOI.empty!(dest)

    vis_src = MOI.get(src, MOI.ListOfVariableIndices())
    idxmap = MOIU.IndexMap()

    has_constraints = BitSet()
    for (F, S) in MOI.get(src, MOI.ListOfConstraints())
        i = get(dest.cone_types_dict, S, nothing)
        if i === nothing || F != MOI.VectorAffineFunction{Float64}
            throw(MOI.UnsupportedConstraint{F, S}())
        end
        push!(has_constraints, i)
    end

    _allocate_variables(dest, vis_src, idxmap)

    # Allocate constraints
    caches = map(collect(has_constraints)) do i
        _allocate_constraints(dest, src, idxmap, i, dest.cone_types[i])
    end

    # Load variables
    _load_variables(dest, length(vis_src))

    # Set variable attributes
    MOIU.pass_attributes(dest, src, copy_names, idxmap, vis_src)

    # Set model attributes
    MOIU.pass_attributes(dest, src, copy_names, idxmap)

    # Load constraints
    offset = 0
    for (i, cache) in zip(has_constraints, caches)
        _load_constraints(dest, src, idxmap, offset, i, cache)
        offset += dest.num_rows[i]
    end

    return idxmap
end


text:
using POMDPs, POMDPPolicies, POMDPGym, Crux
using Flux, Random

mdp = AtariPOMDP(:Pong, version = :v0)
S = ContinuousSpace((80,80,4), Float32) # Needs to be Float32 or the conv layers through an error on gpu
as = actions(mdp)

Q() = DiscreteNetwork(Chain(x->x ./ 255f0, Conv((8,8), 4=>16, relu, stride = 4), Conv((4,4), 16=>32, relu, stride = 2), flatten, Dense(2048, 256, relu), Dense(256, length(as))), as) |> gpu

𝒮 = DQN(π=Q(), S=S, N=5000000, buffer_size=10000, max_steps=1000, buffer_init=5000, log=(fns=[log_undiscounted_return(1)],), c_opt=(optimizer=Flux.Optimiser(ClipValue(1f0), ADAM(1f-3)),))
solve(𝒮, mdp)

# Plot the learning curve
p = plot_learning(𝒮, title = "Pong Training Curve")

# Produce a gif with the final policy
gif(mdp, 𝒮.π, "pong.gif")


text:
using FractrixZeta
using Test

@testset "Placeholder" begin
    @test true
end

max_stars_repo_name: zoemcc/fractrixzeta

max_issues_repo_name: zoemcc/fractrixzeta

max_forks_repo_name: zoemcc/fractrixzeta


text:
using Documenter, ConstLab
ENV["JULIA_DEBUG"] = "Documenter,LiveServer"
# Load remaining dependencies to avoid leaking
# precompiling ... into the docs.
import PGFPlotsX, LaTeXStrings, Literate, Tensors

if haskey(ENV, "GITHUB_ACTIONS")
    PGFPlotsX.latexengine!(PGFPlotsX.PDFLATEX)
    ENV["GITHUB_EVENT_PATH"] = joinpath(@__DIR__, "actions-event.json")
end

# Generate examples
include("generate.jl")

GENERATEDEXAMPLES = [joinpath("examples", f) for f in (
    "von_mises_plasticity_mixed_hardening.md",
    )]

# Build documentation.
makedocs(
    format = Documenter.HTML(
        prettyurls = haskey(ENV, "GITHUB_ACTIONS"), # disable for local builds
        assets = ["assets/custom.css"],
    ),
    sitename = "ConstLab.jl",
    doctest = false,
    strict = false,
    clean = !isdefined(Main, :LiveServer),
    pages = Any[
        "Home" => "index.md",
        "Examples" => GENERATEDEXAMPLES,
        ],
)

# Deploy built documentation from Travis.
deploydocs(
    repo = "github.com/fredrikekre/ConstLab.jl.git",
    push_preview = true,
)


text: @test get_discretization_counts(LinearDiscretizer([0.0,1.0,2.0]), [0.5,0.5,0.5,1.5,1.5]) == [3,2]


text:
module UnitfulMoles

using Unitful

export @mol, @xmol, @compound, @u_str, @unit, uconvert

"""
    @mol(x::Symbol)

The @mol macro lets you declare specific molar units.

Usage example:
```
julia> @mol N
julia> 2μmolN
2 μmolN
```
"""
macro mol(x::Symbol)
    abbr = "mol$x"
    symb = Symbol(abbr)
    name = Symbol("Moles $x")
    equals = :(1u"mol")
    tf = true

    expr = Expr(:block)
    push!(expr.args, :(@unit $symb $abbr $name $equals $tf))
    return esc(expr)
end

"""
    @mol(x::Symbol, grams::Real)

The @mol macro lets you declare specific molar units, with weight in grams.

Usage example:
```
julia> @mol N 14.007
julia> 2μmolN/u"L"
2 L^-1 μmolN
```
"""
macro mol(x::Symbol, grams::Union{Real,Symbol,Expr})
    grams = grams isa Expr ? esc(grams) : grams
    abbr = "mol$x"
    symb = Symbol(abbr)
    name = Symbol("Moles $x")
    equals = :($grams * u"g")
    tf = true

    expr = Expr(:block)
    push!(expr.args, :(@unit $symb $abbr $name $equals $tf))
    return esc(expr)
end

"""
    @compound(compound::Symbol)

The @compound macro lets you declare a compound moles using their chemical
formula. Gram conversions are generated for free from components.

Usage example:
```
julia> @mol H 1.008
molH
julia> @mol O 15.999
molO
julia> @compound H2O
molH2O
julia>
2molH2O |> u"g"
36.03 g
```
"""
macro compound(name::Symbol)
    elements = parse_compound(string(name))

    sum_expr = Expr(:block, :(weight = 0.0u"g"))
    for (el, n) in elements
        x = "mol$el"
        w = :(Unitful.uconvert(u"g", 1 * (@u_str $x)))
        push!(sum_expr.args, :(weight += $n * $w))
    end
    name = Symbol(subscriptify(string(name)))

    weight_expr = :(weight = ustrip($sum_expr))
    mol_expr = :(@mol $name weight)
    return Expr(:block, esc(weight_expr), esc(mol_expr))
end

"""
    @xmol(base::Symbol, compound::Symbol)

The @mol macro lets you declare a fractional compound mole with a particular
base, such as a C-mol.

Usage example:
```
julia> @xmol C C8H10N4O2
julia> 1CmolC8H10N4O2 |> molC8H10N4O2
0.125 molC8H10N4O2
```
"""
macro xmol(base::Symbol, compound::Symbol)
    elements = parse_compound(string(compound))
    n = 1/sum_base(base, elements)

    # Subscriptify compound
    compound = Symbol(subscriptify(string(compound)))

    symb = Symbol(base, "mol", compound)
    abbr = "$base-mol$compound"
    name = Symbol("$base-moles $compound")
    equals = :($n * $(Symbol("mol$compound")))
    tf = true

    expr = Expr(:block)
    push!(expr.args, :(@compound $compound))
    push!(expr.args, :(@unit $symb $abbr $name $equals $tf))
    esc(expr)
end

function sum_base(base, elements::Array{Pair{T,Int},1}) where {T<:AbstractString}
    number = 0
    for (element, n) in elements
        if string(base) == element;
            number += n
        end
    end
    return number
end

# turn integers in compound name into subscripts
# inspired from Unitful.superscript
subscriptify(s::AbstractString) = map(s) do c
    c == '0' ? '\u2080' :
    c == '1' ? '\u2081' :
    c == '2' ? '\u2082' :
    c == '3' ? '\u2083' :
    c == '4' ? '\u2084' :
    c == '5' ? '\u2085' :
    c == '6' ? '\u2086' :
    c == '7' ? '\u2087' :
    c == '8' ? '\u2088' :
    c == '9' ? '\u2089' :
    c
end
# turn subscript integers in compound name into integer characters
desubscriptify(s::AbstractString) = map(s) do c
    c == '\u2082' ? '2' :
    c == '\u2083' ? '3' :
    c == '\u2084' ? '4' :
    c == '\u2085' ? '5' :
    c == '\u2086' ? '6' :
    c == '\u2087' ? '7' :
    c == '\u2088' ? '8' :
    c == '\u2089' ? '9' :
    c == '\u2080' ? '0' :
    c == '\u2081' ? '1' :
    c
end

parse_compound(s::String) = [element_num_pair(Xn) for Xn in split(desubscriptify(s), r"(?=\p{Lu}(?:\p{Ll})*(?:\d)*)")]
element_num_pair(s::AbstractString) = strip(isdigit, s) => num_element(lstrip(isletter, s))
num_element(s::AbstractString) = isempty(s) ? 1 : Base.parse(Int, s)

include("conventionalmoles.jl")

# Allow precompile, and register mol units with u_str macro.
const localunits = Unitful.basefactors
function __init__()
    merge!(Unitful.basefactors, localunits)
    Unitful.register(UnitfulMoles)
end

end # module


text:
using RNGTest
using Test

include("circularRNG.jl")

@time out = RNGTest.smallcrushJulia(generator)
println("\nSmall Crush p-values:\n")
for i in eachindex(out)
  println(out[i])
end
minpv = 1.0
for i in eachindex(out)
  global minpv = min(minpv, minimum(out[i]))
  @test minimum(out[i]) > 1e-4
end
println("Smallest p-value = ", minpv)


text:
function isdivisible(n)
   (0 == (n % 3)) || (0 == (n % 5))
end

function solve()
   sum(filter(isdivisible, 1:999))
end

result = solve()
println(result)

max_stars_repo_name: Ernulphus/polyglot-euler

max_issues_repo_name: Ernulphus/polyglot-euler

max_forks_repo_name: Ernulphus/polyglot-euler


text:
#TODO remove this module from the SOM app to separate concerns
#
# Module for retrieving and storing SOM data in Mongodb

using Mongoc

client = Mongoc.Client()
database = client["nomisnator-test"]
collection = database["names"]
mynames = Mongoc.find(collection) |> collect
names_short = [ x["name"] for x= mynames if x["stats"]["count"] > 5000 && x["sex"] == "F"]

names_short |> unique!

# println(names_short)

function cache_som(som,train_data,row_labels=[],record_name="default")
   database = client["som-cache"]
   collection = database[record_name]
   # document["codes"] = serialize_2d(som.codes)
   codes = serialize_2d(som.codes)
   for i in 1:length(codes)
      document = Mongoc.BSON()
      document["code"] = [i,codes[i]]
      Mongoc.insert_one(collection,document)
   end
   # document["train_data"] = serialize_2d(train_data)
   # document["row_labels"] = row_labels
end

function serialize_2d(array)
   array = convert(Matrix,array)
   serialized = []
   for i in 1:nrow(array)
      push!(serialized,array[1,:])
   end
   serialized
end

function deserialize2d(array::Array)
   ncols = length(array[1][1])
   out = reshape([],(0,ncols))
   for row in array
      r = reshape(row,(1,ncols))
      out = vcat(out,r)
   end
   out
end

function restore_som(som,set="default")
   database = client["som-cache"]
   collection = database[set]
   # som_object = collection |> Mongoc.find |> collect |> last
   som_object = collection |> Mongoc.find |> collect
   # som_object["codes"]
end

# function cache_test(name)
#    database = client["test"]
#    collection = database["stuff"]
#    document = Mongoc.BSON()
#    document["name"] = name
#    Mongoc.insert_one(collection,document)
# end
#
# function find_test()
#    database = client["test"]
#    collection = database["stuff"]
#    Mongoc.find(collection) |> collect
# end


text:
module AutoDiff
# (c) David Barber, University College London 2015

push!(LOAD_PATH,pwd())
#push!(LOAD_PATH, joinpath(pwd(), "src"))

# push all subdirectories from src
#map(d -> push!(LOAD_PATH, joinpath(pwd(), "src", d)),
#    filter(d -> isdir(joinpath("src", d)), readdir("src")))#
#
#push!(LOAD_PATH, joinpath(pwd(), "Demos"))

include("gpumacros.jl")
#export @gpu, @cpu

include("initfile.jl")

@cpu println("Using CPU")
#@gpu println("Compiling kernels...")
#@gpu include("compile_kernels.jl")
@gpu println("Using GPU")

@gpu (using CUDArt; println("using CUDArt"))
@gpu (using CUBLAS; println("using CUBLAS"))

@gpu (import CUDArt.to_host; to_host(A::Array)=A; export to_host;) # CHECKTHIS
@cpu (to_host(A::Array)=A; export to_host;) # CHECKTHIS

using Reexport
@gpu (@reexport using CUDArt;println("reexport using CUDArt"))

global nodecounter
export nodecounter

global GPU

@gpu ArrayOrCudaArray = Union{Array,CudaArray}
@cpu ArrayOrCudaArray = Array

function StartCode()
    global nodecounter = 0
    global node=[]
end
export StartCode

function EndCode()
    return network()
end
export EndCode

function NodeCounter()
    global nodecounter
    nodecounter
end
export NodeCounter

function Node()
    global node
    node
end
export Node

type ADnode
    index #node index
    parents # node parent indices
    f::Function   # function that the node computes
    f_inplace::Function   # in place version of function
    df::Function  # derivative function
    children::Array{Int,1} # node child indices
    takederivative # whether to take the derivative
    returnderivative::Bool # whether to return the derivative
#    input::Bool # whether this is an input variable
    isconst::Bool
    constval
   #ADnode(f=nx,parents=[];returnderivative=false,isconst=false,constval=nothing)=
       ADnode(f=nx,parents=[];returnderivative=false,isconst=false,constval=nothing)=
        begin
            global nodecounter+=1
            global node
            #if f==nx
            #    input=true
            #else
            #    input=false
            #end
            returnderivative=returnderivative==true
            takederivative=returnderivative
#            if returnderivative & !input
#                error("cannot return derivative for a node that has parents")
#            end
            if isa(parents,Array{ADnode})
                parents=map(n->n.index,parents)
            end
            if isa(parents,ADnode)
                parents=parents.index
            end
            thisnode=new(nodecounter,collect(parents),f,Inplace[f],Derivative[f],[0],takederivative,returnderivative,isconst,constval)
            if isempty(node)
                node=Array(Any,0)
            end
            push!(node,thisnode)
            return thisnode
        end
end

abstract ADdummy

type ADtrans<: ADdummy # transpose node. Dummy node that can be used for code optimisation
    index # node index
    parent # node parent index
#    input::Bool # we set this to true since this prevents dummy nodes being differentiated
    ADtrans(parentnode)=
        begin
            global nodecounter+=1
            global node
            if isa(parentnode,ADnode)
                parent=parentnode.index
            end
            #thisnode=new(nodecounter,parent,true)
            thisnode=new(nodecounter,parent)
            if isempty(node)
                node=Array(Any,0)
            end
            push!(node,thisnode)
            return thisnode
        end
end
export ADtrans

type ADdiag<:ADdummy# Dummy diag node that can be used for code optimisation
    index # node index
    parent # node parent index
#    input::Bool # we set this to true since this prevents dummy nodes being differentiated
    ADdiag(parentnode)=
        begin
            global nodecounter+=1
            global node
            if isa(parentnode,ADnode)
                parent=parentnode.index
            end
            #thisnode=new(nodecounter,parent,true)
            thisnode=new(nodecounter,parent)
            if isempty(node)
                node=Array(Any,0)
            end
            push!(node,thisnode)
            return thisnode
        end
end
export ADdiag

ArrayADnode=Array{ADnode}
ADnodeOrArrayADnode=Union{ADnode,Array{ADnode}}
export ADnodeOrArrayADnode, ArrayADnode

import Base.getindex
getindex(x::Array,A::ADnode)=getindex(x,A.index)
export getindex

import Base.setindex!
setindex!(x::Array,value,A::ADnode)=setindex!(x,value,A.index)
export setindex!

function setindex!(x::Array,value::Float64,A::ADnode)
    tmp=cArray((1,1))
    fill!(tmp,value)
    setindex!(x,tmp,A.index)
end
export setindex!

type network
    #node::Array{ADnode,1}
    node::Array{Any,1}
    FunctionNode::Int # Node that forms the scalar function (by default the last node in the graph)
    value
    auxvalue
    gradient
    validnodes
    ancestors
    relevantchildren
    ForwardPassList
    parentIDX
    gpu::Bool
    eltype # Float32 or Float64

   function network()
    vn=find(map((x)->( x!=nothing && !isa(x,ADdummy))  ,Node()))
        return new(Node(),NodeCounter(),Array(Any,NodeCounter()),Array(Any,NodeCounter()),Array(Any,NodeCounter()),vn,nothing,nothing,nothing,nothing,PROC=="GPU",Float64)
    end

    function network(node)
        vn=find(map((x)->( x!=nothing && !isa(x,ADdummy))  ,Node()))
        return new(node,NodeCounter(),Array(Any,NodeCounter()),Array(Any,NodeCounter()),Array(Any,NodeCounter()),vn,nothing,nothing,nothing,nothing,PROC=="GPU",Float64)
    end

    function network(node,FunctionNode,value,auxvalue,gradient,anc,relevantchildren,forwardlist=nothing)
    vn=find(map((x)->( x!=nothing && !isa(x,ADdummy))  ,Node()))
        return new(node,NodeCounter(),value,auxvalue,gradient,vn,anc,relevantchildren,forwardlist,nothing,PROC=="GPU",Float64)
    end

end

include("utils.jl")
include("CUDAutils.jl")

include("defs.jl")
include("compile.jl")

include("ADforward!.jl")
include("ADbackward!.jl")

include("gradcheckGPU.jl"); export gradcheckGPU
include("gradcheckCPU.jl"); export gradcheckCPU

include("netutils.jl")

function gradcheck(net;showgrad=false)
    if net.gpu
        gradcheckGPU(net;showgrad=showgrad)
    else
        gradcheckCPU(net;showgrad=showgrad)
    end
end
export gradcheck

export matread, jldopen, matopen
export ADnode, network, compile
export ADforward!, ADbackward!
export gradcheck
export ArrayOrCudaArray

#ADvariable(;returnderivative=true)=ADnode(;returnderivative=returnderivative)
ADvariable(;returnderivative=true)=ADnode(;returnderivative=returnderivative)
export ADvariable

# make the following form an array if constval is a scalar:
ADconst(constval)=ADnode(;returnderivative=false,isconst=true,constval=Float64(constval))
ADint(constval)=ADnode(;returnderivative=false,isconst=true,constval=Int(constval))
export ADconst
export ADint

#export @gpu, @cpu

end

#function MapReduce(f,op,node::Array{ADnode,1})
#    # f is the function mapped onto each node
#    # op is the binary reduction
#    f(node[1])
#    for i=2:length(node)
#        oldcounter=NodeCounter()
#        f(node[i])
#        op(oldcounter, NodeCounter())
#    end
#    return Node()[end]
#end
#export MapReduce

#function endnode(node)
#    d=falses(length(node))
#    for i=1:length(d)
#        d[i]=isdefined(node,i);
#    end
#    return last(find(d))
#end
#export endnode

# Source code generation:
#include("genFcode.jl")
#include("genRcode.jl")
#export genFcode, genRcode

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
function verif(n::Int64,f)
	if n < 0
        return true
    elseif f(n)
        return verif(n-1,f)
    else
        return false
    end
end

function f(x) true end
function g(x) false end
function h(x) x <= 43 end

println(verif(314151692, f))
println(verif(1, f))
println(verif(0, f))
println(verif(-42, f))
println(verif(314151692, g))
println(verif(1, g))
println(verif(0, g))
println(verif(-42, g))
println(verif(314151692, h))
println(verif(1, h))
println(verif(0, h))
println(verif(-42, h))
println(verif(42, h))
println(verif(43, h))

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
function read_fasta(filename::AbstractString, max_gap_fraction::Real, theta::Any, remove_dups::Bool)
    Z = read_fasta_alignment(filename, max_gap_fraction)
    if remove_dups
        Z, _ = remove_duplicate_sequences(Z)
    end
    N, M = size(Z)
    q = round(Int, maximum(Z))
    W, Meff = compute_weights(Z,  theta)
    return W, Z, N, M, q
end

function computep0(var)
    @extract var:W Z q
    p0 = zeros(q)
    for i in 1:length(W)
        p0[Z[1,i]] += W[i]
    end
    p0
end

function compute_empirical_freqs(Z::AbstractArray{Ti,2}, W::AbstractVector{Float64}, q::Ti) where Ti <: Integer
    N, M = size(Z)
    f = zeros(q, N)
    @inbounds for i in 1:N
        for s in 1:M
            f[Z[i,s],i] += W[s]
        end
    end
    f
end

function entropy(Z::AbstractArray{Ti,2}, W::AbstractVector{Float64}) where Ti <: Integer
    N, M = size(Z)
    q = maximum(Z)
    f = compute_empirical_freqs(Z, W, q)
    S = zeros(N)
    @inbounds for i in 1:N
        _s = 0.0
        for a in 1:q
            _s -= f[a,i] > 0 ? f[a,i] * log(f[a,i]) : 0.0
        end
        S[i] = _s
    end
    S
end

function unpack_params(θ, arvar::ArVar)
    @extract arvar:N q

    arrJ = Array{Float64,3}[]
    arrH = Vector{Float64}[]
    ctr = 0
    for site in 1:N - 1
        _arrJ = zeros(q, q, length(1:site))
        for i in 1:site
            for a in 1:q
                for b in 1:q
                    ctr += 1
                    _arrJ[b,a,i] = θ[ctr]
                end
            end
        end
        push!(arrJ, _arrJ)
        _arrH = zeros(q)
        for a in 1:q
            ctr += 1
            _arrH[a] = θ[ctr]
        end
        push!(arrH, _arrH)
    end
    @assert ctr == length(θ)
    computep0(arvar), arrJ, arrH
end

function softmax!(r::Vector{Float64}, x::Vector{Float64})
    n = length(x)
    u = maximum(x)
    s = 0.
    @inbounds @avx for i = 1:n
        r[i] = exp(x[i] - u)
        s += r[i]
    end
    invs =  inv(s)
    @inbounds @avx for i = 1:n
        r[i] *= invs
    end
    r
end

"""
    softmax(x::AbstractArray{<:Real})
Return the [`softmax transformation`](https://en.wikipedia.org/wiki/Softmax_function) applied to `x`
"""
softmax!(x::Vector{Float64}) = softmax!(x, x)
softmax(x::Vector{Float64}) = softmax!(similar(x, Float64), x)

"""
    sample(arnet::ArNet, msamples::Int)
Return a generated alignment in the form of a `N × msamples`  matrix of type `::Matrix{Int}`
# Examples
```
julia> arnet,arvar=ardca("file.fasta",verbose=true,permorder=:ENTROPIC, lambdaJ=0.001,lambdaH=0.001);
julia> Zgen=Zgen=sample(arnet,1000);
```
"""
function sample(arnet::ArNet, msamples::Int)
    @extract arnet:H J p0 idxperm
    q = length(p0)
    N = length(H) # here N is N-1 !!
    backorder = sortperm(idxperm)
    totH = Vector{Float64}(undef, q)
    res = @inbounds @distributed hcat for i in 1:msamples
        sample_z = Vector{Int}(undef, N + 1)
        sample_z[1] = wsample(1:q, p0)
        for site in 1:N
            Js = J[site]
            h = H[site]
            copy!(totH,h)
            @avx for i in 1:site
                for a in 1:q
                    totH[a] += Js[a,sample_z[i],i]
                end
            end
            p = softmax(totH)
            sample_z[site + 1] = wsample(1:q, p)
        end
        sample_z
    end
    permuterow!(res, backorder)
end

"""
    sample_with_weights(arnet::ArNet, msamples::Int)
Return a generated alignment in the form of a `N × msamples`  matrix of type `::Matrix{Int}` and the relative probabilities under the module

# Examples
```
julia> arnet,arvar=ardca("file.fasta",verbose=true,permorder=:ENTROPIC, lambdaJ=0.001,lambdaH=0.001);
julia> Wgen,Zgen=sample_with_weights(arnet,1000);
```
"""
function sample_with_weights(arnet::ArNet, msamples)
    @extract arnet:H J p0 idxperm
    q = length(p0)
    N = length(H) # here N is N-1 !!
    backorder = sortperm(idxperm)
    totH = Vector{Float64}(undef, q)
    W = SharedArray{Float64}(msamples)
    res = @inbounds @distributed hcat for i in 1:msamples
        sample_z = Vector{Int}(undef, N + 1)
        sample_z[1] = wsample(1:q, p0)
        logw = log(p0[sample_z[1]])
        for site in 1:N
            Js = J[site]
            h = H[site]
            copy!(totH,h)
            @avx for i in 1:site
                for a in 1:q
                    totH[a] += Js[a,sample_z[i],i]
                end
            end
            p = softmax(totH)
            sample_z[site + 1] = wsample(1:q, p)
            logw += log(p[sample_z[site + 1]])
        end
        W[i] = exp(logw)
        sample_z
    end
    permuterow!(res, backorder)
    return W , res
end

function permuterow!(x::AbstractMatrix, p::Vector)
    isperm(p) || error("not a permutation")
    for j in 1:size(x, 2)
        vx = @view x[:,j]
        Base.permute!(vx, p)
    end
    x
end

function permuterow!(x::AbstractVector, p::Vector)
    isperm(p) || error("not a permutation")
    Base.permute!(x, p)
end

# warning the gauge of H[:,1] is to be determined !!
function tensorize(arnet::ArNet; tiny::Float64=1e-16)
    @extract arnet:J H idxperm p0
    N = length(idxperm)
    q = length(H[1])
    p0pc = (1.0-tiny)*p0 .+ tiny/q
    outJ = zeros(q, q, N, N)
    outH = zeros(q, N)
    shiftH0 = sum(log.(p0pc)) / q
    outH[:,idxperm[1]] .= log.(p0pc) .- shiftH0
    for i in 1:N - 1
        si = idxperm[i + 1]
        Js = J[i]
        outH[:,si] .= H[i]
        for j in 1:i
            sj = idxperm[j]
            outJ[:,:,si,sj] .= Js[:,:,j]
            outJ[:,:,sj,si] .= Js[:,:,j]'
        end
    end
    outJ, outH
end

"""
    dms_single_site(arnet::ArNet, arvar::ArVar, seqid::Int; pc::Float64=0.1)

Return a `q×L` matrix of containing `-log(P(mut))/log(P(seq))` for all single
site mutants of the reference sequence `seqid`, and a vector of the indices of
the residues of the reference sequence that contain gaps (i.e. the 21
amino-acid) for which the score has no sense and is set by convention to `+Inf`.
A negative value indicate a beneficial mutation, a value 0 indicate
the wild-type amino-acid.
"""
function dms_single_site(arnet::ArNet, arvar::ArVar, seqid::Int; pc::Float64=0.1)
    @extract arnet : H J p0 idxperm
    @extract arvar : Z M N q

    1 ≤ seqid ≤ M || error("seqid=$seqid should be in the interval [1,...,$M]")

    ppc = (1-pc) * p0 + pc * ones(q)/q
    Da = fill(Inf64,q,N)
    xori = Z[:,seqid]
    xmut = copy(xori)
    idxnogap=findall(x->x!=21,xori)
    arlike = zeros(N)
    arlike0 = zeros(N)
    _outputarnet!(arlike0,xori,J,H,ppc,N,q)
    ll0 = -sum(log.(arlike0))

    @inbounds for i in idxnogap
        if xori[i] == 21
            continue
        end
        for a in 1:q
            if a != xori[i]
                xmut[i] = a
                _outputarnet!(arlike,xmut, J, H, ppc, N, q)
                Da[a,i] = -sum(log.(arlike)) - ll0
            else
                Da[a,i] = 0.0
            end
        end
        xmut[i] = xori[i] #reset xmut to the original velue
    end
    invperm = sortperm(idxperm)
    return Da[:,invperm],sort!(idxperm[setdiff(1:N,idxnogap)])
end

max_issues_repo_name: pagnani/ArDCA

max_forks_repo_name: pagnani/ArDCA


text:
function VectorizationBase.vectorizable(ScA::ScatteredArray{E,M,T,N,Np1}) where {E,M,T,N,Np1}
    ptr_data = pointer(ScA.data)
    VectorizedScatteredArray{E,M,T,N,Np1}(ptr_data, size(ScA.data))
end

const ScatteredVector{E,M,T} = ScatteredArray{E,M,T,1,2}
const ScatteredMatrix{E,M,T} = ScatteredArray{E,M,T,2,3}

function ScA_setindex_quote(N, L)
    q = quote $(Expr(:meta, :inline)) end
    inds = [:(i[$j]) for j ∈ 1:N]
    for j ∈ 1:L
        push!(q.args, Expr(:call, :setindex!, :(ScA.data), :(v[LinearStorage(), $j]), inds..., j))
    end
    q
end

@generated function Base.setindex!(ScA::ScatteredArray{E,M,T,N,Np1}, v::T,
                                   # i::Vararg{<:Integer,N}) where {T,E,M,N,Np1}
                                   i::Vararg{<:Integer,N}) where {E,M,T,N,Np1}
    ScA_setindex_quote(N, type_length(T))
end

@generated function Base.setindex!(ScA::ScatteredArray{E,M,T,N,Np1}, v::T,
                                    # i::CartesianIndex{N}) where {T,E,M,N,Np1}
                                    i::MultiDimIndex{N}) where {E,M,T,N,Np1}
    ScA_setindex_quote(N, type_length(T))
end

@generated function Base.getindex(ScA::ScatteredArray{E,M,T,N,Np1}, i::Vararg{<:Integer,N}) where {E,M,T,N,Np1}
    inds = [:(i[$n]) for n in 1:N]
    ind_expr = [Expr(:call, :getindex, :(ScA.data), inds..., j) for j in 1:type_length(T)]
    quote
        $(Expr(:meta, :inline))
        $(construct_expr(T, ind_expr))
    end
end

@generated function Base.getindex(ScA::ScatteredArray{E,M,T,N,Np1},
                        i::MultiDimIndex{N}) where {E,M,T,N,Np1}
    inds = [:(i[$n]) for n in 1:N]
    ind_expr = [Expr(:call, :getindex, :(ScA.data), inds..., j) for j in 1:type_length(T)]
    quote
        $(Expr(:meta, :inline))
        $(construct_expr(T, ind_expr))
    end
end

function scattered_array_view_inds(V, N)
    inds = Expr[]
    i = 0
    for n ∈ 1:N
        if V.parameters[n] == Colon
            i += 1
            push!(inds, :(i[$i]))
        else
            push!(inds, :(ScA.view[$n]))
        end
    end
    inds
end

@generated function Base.getindex(ScA::ScatteredArrayView{E,M,T,new_N,N,Np1,V},
                        i::Vararg{<:Any,new_N}) where {E,M,T,new_N,N,Np1,V}
    inds = scattered_array_view_inds(V, N)
    # The strategy is to calculate the first index, as well as the stride between indices.
    # After the first index is calculated, the remaining indices can be calcualted much more
    # efficiently via the stide.
    q = quote
        $(Expr(:meta, :inline))
        stride = $(Expr(:call, :*, [:(ScA.full_size[$n]) for n ∈ 1:N]...))
        ind_1 = sub2ind(ScA.full_size, $(Expr(:tuple,inds...,1)))
    end
    for j ∈ 2:type_length(T)
        push!(q.args, :($(Symbol(:ind_,j)) = ind_1 + stride * $(j-1)))
    end
    push!(q.args, construct_expr(T, [:(unsafe_load(ScA.ptr, $(Symbol(:ind_,j)))) for j ∈ 1:type_length(T)]))
    q
end
@generated function Base.setindex!(ScA::ScatteredArrayView{E,M,T,new_N,N,Np1,V}, v,
                        i::Vararg{<:Any,new_N}) where {E,M,T,new_N,N,Np1,V}
    inds = scattered_array_view_inds(V, N)
    q = quote
        $(Expr(:meta, :inline))
        stride = $(Expr(:call, :*, [:(ScA.full_size[$n]) for n ∈ 1:N]...))
        ind_1 = sub2ind(ScA.full_size, $(Expr(:tuple,inds...,1)))
        unsafe_store!(ScA.ptr, v[LinearStorage(), 1], ind_1)
    end
    for j ∈ 2:type_length(T)
        push!(q.args, quote
            unsafe_store!(ScA.ptr, v[LinearStorage(), $j], ind_1 + stride * $(j-1))
        end)
    end
    push!(q.args, :(v))
    q
end
@inline function Base.view(ScA::ScatteredArray{E,M,T,N,Np1}, i::Vararg{<:Any,N}) where {E,M,T,N,Np1}
    ScatteredArrayView(ScA, i)
end
@generated function ScatteredArrayView(ScA::ScatteredArray{E,M,T,N,Np1}, i::V) where {E,M,T,N,Np1,V}
    new_N = 0
    new_size = Expr(:tuple)
    for n ∈ 1:N
        if V.parameters[n] == Colon
            new_N += 1
            push!(new_size.args, :(full_size[$n]))
        end
    end
    quote
        full_size = size(ScA.data)
        ScatteredArrayView{$E,$M,$T,$new_N,$N,$Np1,$V}(pointer(ScA.data), $new_size, full_size, i)
    end
end

# @inline function SArray{Tuple{S},T,N,L}(vs::Vararg{SIMDPirates.SVec{W,T},L}) where {S,T,N,L,W}
#     SArray{S}(vs...)
# end

function scattered_vload_quote(N, W, E, T, V)
    inds = [:(i[$n]) for n in 1:N]
    quote
        $(Expr(:meta,:inline))
        $(construct_expr(T, [Expr(:call, :vload, V, :(vScA.ptr),
            Expr(:call, :sub2ind, :(vScA.size), Expr(:tuple, inds..., j))) for j in 1:type_length(T)]))
    end
end

@generated function SIMDPirates.vload(::Type{SIMDPirates.SVec{W,E}}, vScA::VectorizedScatteredArray{E,M,T,N,Np1},
                                                # i::CartesianIndex{N}) where {E,M,T,N,Np1,W}
                                                i::Vararg{<:Integer,N}) where {E,M,T,N,Np1,W}
    scattered_vload_quote(N, W, E, T, SIMDPirates.SVec{W,E})
end
@generated function SIMDPirates.vload(::Type{SIMDPirates.SVec{W,E}}, vScA::VectorizedScatteredArray{E,M,T,N,Np1},
                                                # i::Vararg{<:Integer,N}) where {E,M,T,N,Np1,W}
                                                i::MultiDimIndex{N}) where {E,M,T,N,Np1,W}
    scattered_vload_quote(N, W, E, T, SIMDPirates.SVec{W,E})
end

@generated function SIMDPirates.vload(::Type{SIMDPirates.SVec{W,E}}, vScA::VectorizedScatteredArray{E,M,T,N,Np1},
                                                # i::Vararg{<:Integer,N}) where {E,M,T,N,Np1,W}
                                                i::Integer) where {E,M,T,N,Np1,W}
    quote
        $(Expr(:meta,:inline))
        $(construct_expr(T, [Expr(:call, :vload, SIMDPirates.SVec{W,E}, :(vScA.ptr),
            :(i + length(vScA)*$(j-1)) ) for j in 1:type_length(T)]))
    end
end
@generated function SIMDPirates.vload(
                ::Type{SIMDPirates.SVec{W,E}}, vScA::VectorizedScatteredArray{E,M,T,N,Np1}
            ) where {E,M,T,N,Np1,W}
    quote
        $(Expr(:meta,:inline))
        $(construct_expr(T, [Expr(:call, :vload, SIMDPirates.SVec{W,E}, :(vScA.ptr + length(vScA)*$(sizeof(E)*(j-1))) ) for j in 1:type_length(T)]))
    end
end
@generated function SIMDPirates.vload(
                ::Type{SIMDPirates.SVec{W,Ef1}}, vScA::VectorizedScatteredArray{Ef2,M,T,N,Np1}
            ) where {Ef1,Ef2,M,T,N,Np1,W}
    quote
        $(Expr(:meta,:inline))
        $(construct_expr(T,
            [Expr(:call, :convert, SVec{W,Ef1},
                Expr(:call, :vload, SIMDPirates.SVec{W,Ef2}, :(vScA.ptr + length(vScA)*$(sizeof(Ef2)*(j-1))) )) for j in 1:type_length(T)
        ])
        )
    end
end

@inline function Base.:+(i::Integer, v::VectorizedScatteredArray{E,M,T,N,Np1}) where {E,M,T,N,Np1}
    VectorizedScatteredArray{E,M,T,N,Np1}(v.ptr + sizeof(E) * i, v.size)
end
@inline function Base.:+(v::VectorizedScatteredArray{E,M,T,N,Np1}, i::Integer) where {E,M,T,N,Np1}
    VectorizedScatteredArray{E,M,T,N,Np1}(v.ptr + sizeof(E) * i, v.size)
end

function Base.copyto!(ScA::ScatteredArray{E,M,T,N,Np1}, A::AbstractArray{T}) where {E,M,T,N,Np1}
    @boundscheck size(ScA) == size(A) || throwdimensionmismatcherror(size(ScA), size(A))
    @inbounds for i in eachindex(A)
        ScA[i] = A[i]
    end
    ScA
end

@generated function ScatteredArray(A::AbstractArray{T,N}) where {E,M,T <: AbstractArray{E,M},N}
    isbitstype(T) || thrownotisbitserror(T)
    TL = type_length(T)
    quote
        ScA = ScatteredArray{$E,$M,$T,$N,$(N+1)}(
            Array{$E}(undef, $(Expr(:tuple, [:(size(A,$n)) for n ∈ 1:N]..., TL)) )
        )
        @inbounds copyto!(ScA, A)
    end
end
@generated function ScatteredArray{E,M,T,N,Np1}(::UndefInitializer, I::Vararg{<:Integer,N}) where {E,M,T <: AbstractArray{E,M},N,Np1}
    isbitstype(T) || thrownotisbitserror(T)
    N + 1 == Np1 || throw("N + 1 = $(N + 1) != Np1 = $Np1")
    # @show T
    TL = type_length(T)
    quote
        ScatteredArray{$E,$M,$T,$N,$Np1}(
            Array{$E}(undef, $(Expr(:tuple, [:(I[$n]) for n ∈ 1:N]..., TL)) )
        )
    end
end

@generated function Base.size(ScA::ScatteredArray{E,M,T,N,Np1}) where {E,M,T,N,Np1}
    quote
        $(Expr(:meta, :inline))
        $(Expr(:tuple, [:(size(ScA,$n)) for n ∈ 1:N]...))
    end
end
@inline Base.size(ScA::ScatteredArray, n::Integer) = size(ScA.data, n)


text:
@testset "DiffOps" begin
  # uniform weights for simple mesh
  points = P2[(0,0), (1,0), (0,1), (1,1), (0.5,0.5)]
  connec = connect.([(1,2,5),(2,4,5),(4,3,5),(3,1,5)], Triangle)
  mesh = SimpleMesh(points, connec)
  L = laplacematrix(mesh, weights=:uniform)
  @test L == [
     -1 1/3 1/3   0 1/3
    1/3  -1   0 1/3 1/3
    1/3   0  -1 1/3 1/3
      0 1/3 1/3  -1 1/3
    1/4 1/4 1/4 1/4  -1
  ]

  # cotangent weights only defined for triangle meshes
  points = P2[(0,0), (1,0), (1,1), (0,1)]
  connec = connect.([(1,2,3,4)], Quadrangle)
  mesh = SimpleMesh(points, connec)
  @test_throws AssertionError laplacematrix(mesh, weights=:cotangent)

  # full Laplace-Beltrami operator
  sphere = Sphere(P3(0,0,0), T(1))
  mesh = triangulate(sphere)
  L = laplacematrix(mesh)
  M = measurematrix(mesh)
  @test issymmetric(L)
  @test issparse(L)
  @test isdiag(M)
end


text:
using ModernGL, GLFW, FixedPointNumbers

include("GLAbstraction/GLAbstraction.jl")
using .GLAbstraction

const atlas_texture_cache = Dict{GLFW.Window, Tuple{Texture{Float16, 2}, Function}}()

function get_texture!(atlas)
    # clean up dead context!
    filter!(atlas_texture_cache) do (ctx, tex_func)
        if GLAbstraction.context_alive(ctx)
            true
        else
            AbstractPlotting.remove_font_render_callback!(tex_func[2])
            false
        end
    end
    tex, func = get!(atlas_texture_cache, GLAbstraction.current_context()) do
        tex = Texture(
                atlas.data,
                minfilter = :linear,
                magfilter = :linear,
                anisotropic = 16f0,
        )
        # update the texture, whenever a new font is added to the atlas
        function callback(distance_field, rectangle)
            ctx = tex.context
            if GLAbstraction.context_alive(ctx)
                prev_ctx = GLAbstraction.current_context()
                GLAbstraction.switch_context!(ctx)
                tex[rectangle] = distance_field
                GLAbstraction.switch_context!(prev_ctx)
            end
        end
        AbstractPlotting.font_render_callback!(callback)
        return (tex, callback)
    end
    tex
end

function Base.display(x::OpenGLBackend, scene::Scene)
    display(global_gl_screen(), scene)
end

include("GLVisualize/GLVisualize.jl")
using .GLVisualize

include("glwindow.jl")
include("screen.jl")
include("rendering.jl")
include("events.jl")
include("drawing_primitives.jl")


text:
####
#### MulticlassLDA
####
"""
    LDA(; kwargs...)

$LDA_DESCR

# Keyword Parameters

- `method::Symbol=:gevd`:  choice of solver, one of `:gevd` or `:whiten` methods
- `cov_w::CovarianceEstimator`=SimpleCovariance: an estimator for the within-class
    covariance (used in computing within-class scatter matrix, Sw), by default set
    to the standard `MultivariateStats.SimpleCovariance()` but
    could be set to any robust estimator from `CovarianceEstimation.jl`.
- `cov_b::CovarianceEstimator`=SimpleCovariance: same as `cov_w` but for the between-class
    covariance (used in computing between-class scatter matrix, Sb)
- `out_dim::Int=0`: the output dimension, i.e dimension of the transformed space,
    automatically set if 0 is given (default).
- `regcoef::Float64=1e-6`: regularization coefficient (default value 1e-6). A positive
    value `regcoef * eigmax(Sw)` where `Sw` is the within-class scatter matrix, is added
    to the diagonal of Sw to improve numerical stability. This can be useful if using
    the standard covariance estimator.
- `dist::SemiMetric=SqEuclidean`: the distance metric to use when performing classification
    (to compare the distance between a new point and centroids in the transformed space),
    an alternative choice can be the `CosineDist`.Defaults to `SqEuclidean`

See also the
[package documentation](https://multivariatestatsjl.readthedocs.io/en/latest/lda.html).
For more information about the algorithm, see the paper by Li, Zhu and Ogihara,
[Using Discriminant Analysis for Multi-class Classification:
An Experimental Investigation](http://citeseerx.ist.psu.edu/viewdoc/
download?doi=10.1.1.89.7068&rep=rep1&type=pdf).
"""
@mlj_model mutable struct LDA <: MMI.Probabilistic
    method::Symbol = :gevd::(_ in (:gevd, :whiten))
    cov_w::CovarianceEstimator = MS.SimpleCovariance()
    cov_b::CovarianceEstimator = MS.SimpleCovariance()
    out_dim::Int = 0::(_ ≥ 0)
    regcoef::Float64 = 1e-6::(_ ≥ 0)
    dist::SemiMetric = SqEuclidean()
end

function MMI.fit(model::LDA, ::Int, X, y)
    Xm_t, yplain, classes_seen, p, n, nc, nclasses, integers_seen, out_dim =
        _check_lda_data(model, X, y)
    core_res = MS.fit(
        MS.MulticlassLDA, nc, Xm_t, Int.(yplain);
        method=model.method,
        outdim=out_dim,
        regcoef=model.regcoef,
        covestimator_within=model.cov_w,
        covestimator_between=model.cov_b
    )
    cache = nothing
    report = (
        classes=classes_seen,
        out_dim=MS.outdim(core_res),
        class_means=MS.classmeans(core_res),
        mean=MS.mean(core_res),
        class_weights=MS.classweights(core_res),
        Sw=MS.withclass_scatter(core_res),
        Sb=MS.betweenclass_scatter(core_res),
        nc=nc
    )
    fitresult = (core_res, classes_seen)
    return fitresult, cache, report
end

function _check_lda_data(model, X, y)
    class_list = MMI.classes(y[1]) # Class list containing entries in pool of y.
    nclasses = length(class_list)
    # Class list containing unique entries in y.
    classes_seen = filter(in(unique(y)), class_list)
    nc = length(classes_seen) # Number of classes in pool of y.
    integers_seen = MMI.int(classes_seen)
    # NOTE: copy/transpose.
    Xm_t = _matrix_transpose(model, X) # Now p x n matrix
    yplain = MMI.int(y) # Vector of n ints in {1,..., nclasses}.
    p, n = size(Xm_t)
    # Recode yplain to be in {1,..., nc}
    nc == nclasses || replace!(yplain, (integers_seen .=> 1:nc)...)
    # Check to make sure we have more than one class in training sample.
    # This is to prevent Sb from being a zero matrix.
    if nc <= 1
        throw(
            ArgumentError(
                "The number of unique classes in "*
                "traning sample has to be greater than one"
            )
        )
    end

    # Check to make sure we have more samples than classes.
    # This is to prevent Sw from being the zero matrix.
    if n <= nc
        throw(
            ArgumentError(
                "The number of training samples `n` has"*
                " to be greater than number of unique classes `nc`"
            )
        )
    end
    # Check output dimension default is min(p, nc-1)
    def_outdim = min(p, nc - 1)
    # If unset (0) use the default; otherwise try to use the provided one
    out_dim = ifelse(model.out_dim == 0, def_outdim, model.out_dim)
    # Check if the given one is sensible
    if out_dim > p
        throw(
            ArgumentError(
                "`out_dim` must not be larger than `p`"*
                "where `p` is the number of features in `X`"
            )
        )
    end
    return Xm_t, yplain, classes_seen, p, n, nc, nclasses, integers_seen, out_dim
end

function MMI.fitted_params(::LDA, (core_res, classes_seen))
    return (
        projected_class_means=MS.classmeans(core_res),
        projection_matrix=MS.projection(core_res)
    )
end

function MMI.predict(m::LDA, (core_res, classes_seen), Xnew)
    # projection of Xnew, XWt is nt x o  where o = number of out dims
    # nt = number ot test samples
    XWt = MMI.matrix(Xnew) * core_res.proj
    # centroids in the transformed space, nc x o
    centroids = transpose(core_res.pmeans)

    # compute the distances in the transformed space between pairs of rows
    # the probability matrix Pr is `n x nc` and normalised accross rows
    Pr = pairwise(m.dist, XWt, centroids, dims=1)
    Pr .= Pr .* -1
    # apply a softmax transformation
    softmax!(Pr)
    return MMI.UnivariateFinite(classes_seen, Pr)
end

metadata_model(LDA,
    input=Table(Continuous),
    target=AbstractVector{<:Finite},
    weights=false,
    output=Table(Continuous),
    descr=LDA_DESCR,
    path="$(PKG).LDA"
)

####
#### BayesianLDA
####

"""
    BayesianLDA(; kwargs...)

$BayesianLDA_DESCR

# Keyword Parameters

- `method::Symbol=:gevd`: choice of solver, one of `:gevd` or `:whiten` methods
- `cov_w::CovarianceEstimator=SimpleCovariance()`: an estimator for the within-class
    covariance (used in computing within-class scatter matrix, Sw), by default set to the
    standard `MultivariateStats.CovarianceEstimator` but could be set to any robust
    estimator from `CovarianceEstimation.jl`.
- `cov_b::CovarianceEstimator=SimpleCovariance()`: same as `cov_w` but for the
    between-class covariance(used in computing between-class scatter matrix, Sb).
- `out_dim::Int=0`: the output dimension, i.e dimension of the transformed space,
    automatically set if 0 is given (default).
- `regcoef::Float64=1e-6`: regularization coefficient (default value 1e-6). A positive
value `regcoef * eigmax(Sw)` where `Sw` is the within-class covariance estimator, is added
    to the diagonal of Sw to improve numerical stability. This can be useful if using the
    standard covariance estimator.
- `priors::Union{Nothing, Vector{Float64}}=nothing`: For use in prediction with Baye's rule. If `priors = nothing` then
    `priors` are estimated from the class proportions in the training data. Otherwise it
    requires a `Vector` containing class probabilities with probabilities specified using
    the order given by `levels(y)` where y is the target vector.

See also the [package documentation](
https://multivariatestatsjl.readthedocs.io/en/latest/lda.html).
For more information about the algorithm, see the paper by Li, Zhu and Ogihara,
[Using Discriminant Analysis for Multi-class Classification: An Experimental Investigation](
http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.89.7068&rep=rep1&type=pdf).
"""
@mlj_model mutable struct BayesianLDA <: MMI.Probabilistic
    method::Symbol = :gevd::(_ in (:gevd, :whiten))
    cov_w::CovarianceEstimator=MS.SimpleCovariance()
    cov_b::CovarianceEstimator=MS.SimpleCovariance()
    out_dim::Int=0::(_ ≥ 0)
    regcoef::Float64=1e-6::(_ ≥ 0)
    priors::Union{Nothing, Vector{Float64}}=nothing
end

function MMI.fit(model::BayesianLDA, ::Int, X, y)
    Xm_t, yplain, classes_seen, p, n, nc, nclasses, integers_seen, out_dim =
        _check_lda_data(model, X, y)
    ## If piors are specified check if they makes sense.
    ## This was put here to through errors much earlier
    if isa(model.priors, Vector)
        priors = _check_lda_priors(model.priors, nc, nclasses, integers_seen)
    end

    core_res = MS.fit(
        MS.MulticlassLDA, nc, Xm_t, Int.(yplain);
        method=model.method,
        outdim=out_dim,
        regcoef=model.regcoef,
        covestimator_within=model.cov_w,
        covestimator_between=model.cov_b
    )

    ## Estimates prior probabilities if specified by user.
    ## Put it here to avoid recomputing as the fitting process does this already
    if isa(model.priors, Nothing)
        weights = MS.classweights(core_res)
        total = core_res.stats.tweight
        priors = weights ./ total
    end
    cache     = nothing
    report    = (
        classes=classes_seen,
        out_dim=MS.outdim(core_res),
        class_means=MS.classmeans(core_res),
        mean=MS.mean(core_res),
        class_weights=MS.classweights(core_res),
        Sw=MS.withclass_scatter(core_res),
        Sb=MS.betweenclass_scatter(core_res),
        nc=nc
    )

    fitresult = (core_res, classes_seen, priors, n)
    return fitresult, cache, report
end

function _matrix_transpose(model::Union{LDA, BayesianLDA}, X)
    return MMI.matrix(X; transpose=true)
end

function _check_lda_priors(priors, nc, nclasses, integers_seen)
    if length(priors) != nclasses
        throw(ArgumentError("Invalid size of `priors`."))
    end
    if !isapprox(sum(priors), 1)
        throw(ArgumentError("probabilities specified in `priors` must sum to 1"))
    end
    if any(model.priors .< 0)
        throw(ArgumentError("probabilities specified in `priors` must non-negative"))
    end
    # Select priors for unique classes in `y` (For resampling purporses).
    priors_ = nc == nclasses ? model.priors : @view model.priors[integers_seen]
    return priors_
end

function MMI.fitted_params(::BayesianLDA, (core_res, classes_seen, priors, n))
   return (
       projected_class_means=MS.classmeans(core_res),
       projection_matrix=MS.projection(core_res),
       priors=priors
    )
end

function MMI.predict(m::BayesianLDA, (core_res, classes_seen, priors, n), Xnew)
     # projection of Xnew, XWt is nt x o  where o = number of out dims
    # nt = number ot test samples
    XWt = MMI.matrix(Xnew) * core_res.proj
    # centroids in the transformed space, nc x o
    centroids = transpose(core_res.pmeans)

    # The discriminant matrix `Pr` is of dimension `nt x nc`
    # Pr[i,k] = -0.5*(xᵢ −  µₖ)ᵀ(Σw⁻¹)(xᵢ −  µₖ) + log(priorsₖ) where (Σw = Sw/n)
    # In the transformed space this becomes
    # Pr[i,k] = -0.5*(Pᵀxᵢ −  Pᵀµₖ)ᵀ(PᵀΣw⁻¹P)(Pᵀxᵢ −  Pᵀµₖ) + log(priorsₖ)
    # But PᵀSw⁻¹P = I and PᵀΣw⁻¹P = n*I due to the nature of the projection_matrix, P
    # Giving Pr[i,k] = -0.5*n*(Pᵀxᵢ −  Pᵀµₖ)ᵀ(Pᵀxᵢ −  Pᵀµₖ) + log(priorsₖ)
    # with (Pᵀxᵢ −  Pᵀµₖ)ᵀ(Pᵀxᵢ −  Pᵀµₖ) being the SquaredEquclidean distance between
    # pairs of rows in the transformed space
    Pr = pairwise(SqEuclidean(), XWt, centroids, dims=1)
    Pr .*= (-0.5*n)
    Pr .+= log.(priors)'

    # apply a softmax transformation to convert Pr to a probability matrix
    softmax!(Pr)
    return MMI.UnivariateFinite(classes_seen, Pr)
end

function MMI.transform(m::T, (core_res,), X) where T<:Union{LDA, BayesianLDA}
    # projection of X, XWt is nt x o  where o = out dims
    proj = core_res.projw * core_res.projLDA #proj is the projection_matrix
    XWt = MMI.matrix(X) * proj
    return MMI.table(XWt, prototype = X)
end

metadata_model(
    BayesianLDA,
    input=Table(Continuous),
    target= AbstractVector{<:Finite},
    weights=false,
    output=Table(Continuous),
    descr=BayesianLDA_DESCR,
    path="$(PKG).BayesianLDA"
)

####
#### SubspaceLDA
####

"""
    SubspaceLDA(; kwargs...)

$SubspaceLDA_DESCR

# Keyword Parameters

- `normalize=true`: Option to normalize the between class variance for the number of
    observations in each class, one of `true` or `false`.
- `out_dim`: the dimension of the transformed space to be used by `predict` and
    `transform` methods, automatically set if 0 is given (default).
- `dist=SqEuclidean`: the distance metric to use when performing classification
    (to compare the distance between a new point and centroids in the transformed space),
    an alternative choice can be the `CosineDist`.

See also the [package documentation](
https://multivariatestatsjl.readthedocs.io/en/latest/lda.html).
For more information about the algorithm, see the paper by Howland & Park (2006),
"Generalizing discriminant analysis using the generalized singular value decomposition",
IEEE Trans. Patt. Anal. & Mach. Int., 26: 995-1006.
"""
@mlj_model mutable struct SubspaceLDA <: MMI.Probabilistic
    normalize::Bool = true
    out_dim::Int = 0::(_ ≥ 0)
    dist::SemiMetric = SqEuclidean()
end

function MMI.fit(model::SubspaceLDA, ::Int, X, y)
    Xm_t, yplain, classes_seen, p, n, nc, nclasses, integers_seen, out_dim =
        _check_lda_data(model, X, y)

    core_res = MS.fit(
        MS.SubspaceLDA, Xm_t, Int.(yplain), nc;
        normalize = model.normalize
    )
    # λ is a (nc -1) x 1 vector containing the eigen values sorted in descending order.
    λ = core_res.λ
    explained_variance_ratio = λ ./ sum(λ) #proportions of variance

    cache = nothing
    report = (
        explained_variance_ratio=explained_variance_ratio,
        classes=classes_seen,
        class_means=MS.classmeans(core_res),
        mean=MS.mean(core_res),
        class_weights=MS.classweights(core_res),
        nc=nc
    )
    fitresult = (core_res, out_dim, classes_seen)
    return fitresult, cache, report
end

function MMI.fitted_params(::SubspaceLDA, (core_res, _))
    return (class_means=MS.classmeans(core_res), projection_matrix=MS.projection(core_res))
end

function MMI.predict(m::SubspaceLDA, (core_res, out_dim, classes_seen), Xnew)
     # projection of Xnew, XWt is nt x o  where o = number of out dims
    # nt = number ot test samples
    proj = core_res.projw * view(core_res.projLDA, :, 1:out_dim) #proj is the projection_matrix
    XWt = MMI.matrix(Xnew) * proj
    # centroids in the transformed space, nc x o
    centroids = transpose(core_res.cmeans) * proj

    # compute the distances in the transformed space between pairs of rows
    # the probability matrix is `nt x nc` and normalised accross rows
    Pr = pairwise(m.dist, XWt, centroids, dims=1)
    Pr .= Pr .* -1
    # apply a softmax transformation
    softmax!(Pr)
    return MMI.UnivariateFinite(classes_seen, Pr)
end

metadata_model(
    SubspaceLDA,
    input=Table(Continuous),
    target=AbstractVector{<:Finite},
    weights=false,
    output=Table(Continuous),
    descr=SubspaceLDA_DESCR,
    path="$(PKG).SubspaceLDA"
)

####
#### BayesianSubspaceLDA
####

"""
    BayesianSubspaceLDA(; kwargs...)

$BayesianSubspaceLDA_DESCR

# Keyword Parameters

- `normalize::Bool=true`: Option to normalize the between class variance for the number of
    observations in each class, one of `true` or `false`.
- `out_dim::Int=0`: the dimension of the transformed space to be used by `predict` and
    `transform` methods, automatically set if 0 is given (default).
- `priors::Union{Nothing, Vector{Float64}}=nothing`: For use in prediction with Baye's
    rule. If `priors = nothing` then `priors` are estimated from the class proportions
    in the training data. Otherwise it requires a `Vector` containing class
    probabilities with probabilities specified using the order given by `levels(y)`
    where y is the target vector.

For more information about the algorithm, see the paper by Howland & Park (2006),
"Generalizing discriminant analysis using the generalized singular value decomposition"
,IEEE Trans. Patt. Anal. & Mach. Int., 26: 995-1006.
"""
@mlj_model mutable struct BayesianSubspaceLDA <: MMI.Probabilistic
    normalize::Bool=false
    out_dim::Int= 0::(_ ≥ 0)
    priors::Union{Nothing, Vector{Float64}}=nothing
end

function MMI.fit(model::BayesianSubspaceLDA, ::Int, X, y)
    Xm_t, yplain, classes_seen, p, n, nc, nclasses, integers_seen, out_dim =
        _check_lda_data(model, X, y)
    ## If piors are specified check if they makes sense.
    ## This was put here to through errors much earlier
    if isa(model.priors, Vector)
        priors = _check_lda_priors(model.priors, nc, nclasses, integers_seen)
    end

    core_res = MS.fit(
        MS.SubspaceLDA, Xm_t, Int.(yplain), nc;
        normalize = model.normalize
    )
    # λ is a (nc -1) x 1 vector containing the eigen values sorted in descending order.
    λ = core_res.λ
    explained_variance_ratio = λ ./ sum(λ) #proportions of variance
    mult = model.normalize ? n : 1 #used in prediction

    ## Estimates prior probabilities if specified by user.
    ## Put it here to avoid recomputing as the fitting process does this already
    if isa(model.priors, Nothing)
        weights = MS.classweights(core_res)
        priors = weights ./ n
    end

    cache = nothing
    report = (
        explained_variance_ratio=explained_variance_ratio,
        classes=classes_seen,
        class_means=MS.classmeans(core_res),
        mean=MS.mean(core_res),
        class_weights=MS.classweights(core_res),
        nc=nc
    )
    fitresult = (core_res, out_dim, classes_seen, priors, n, mult)
    return fitresult, cache, report
end

function _matrix_transpose(model::Union{SubspaceLDA, BayesianSubspaceLDA}, X)
    return transpose(MMI.matrix(X))
end

function MMI.fitted_params(::BayesianSubspaceLDA, (core_res, _, _, priors,_))
    return (
        projected_class_means=MS.classmeans(core_res),
        projection_matrix=MS.projection(core_res),
        priors=priors
    )
end

function MMI.predict(
    m::BayesianSubspaceLDA,
    (core_res, out_dim, classes_seen, priors, n, mult),
    Xnew
)
    # projection of Xnew, XWt is nt x o  where o = number of out dims
    # nt = number ot test samples
    #proj is the projection_matrix
    proj = core_res.projw * view(core_res.projLDA, :, 1:out_dim)
    XWt = MMI.matrix(Xnew) * proj

    # centroids in the transformed space, nc x o
    centroids = transpose(core_res.cmeans) * proj
    nc = length(classes_seen)

    # compute the distances in the transformed space between pairs of rows
    # The discriminant matrix `Pr` is of dimension `nt x nc`
    # Pr[i,k] = -0.5*(xᵢ −  µₖ)ᵀ(Σw⁻¹)(xᵢ −  µₖ) + log(priorsₖ) where Σw = Sw/(n-nc)
    # In the transformed space this becomes
    # Pr[i,k] = -0.5*(Pᵀxᵢ −  Pᵀµₖ)ᵀ(PᵀΣw⁻¹P)(Pᵀxᵢ −  Pᵀµₖ) + log(priorsₖ)
    # But PᵀSw⁻¹P = (1/mult)*I and PᵀΣw⁻¹P = (n-nc)/mult*I
    # Giving Pr[i,k] = -0.5*n*(Pᵀxᵢ −  Pᵀµₖ)ᵀ(Pᵀxᵢ −  Pᵀµₖ) + log(priorsₖ)
    # (Pᵀxᵢ −  Pᵀµₖ)ᵀ(Pᵀxᵢ −  Pᵀµₖ) is the SquaredEquclidean distance in the
    # transformed space
    Pr = pairwise(SqEuclidean(), XWt, centroids, dims=1)
    Pr .*= (-0.5 * (n-nc)/mult)
    Pr .+= log.(priors)'

    # apply a softmax transformation to convert Pr to a probability matrix
    softmax!(Pr)
    return MMI.UnivariateFinite(classes_seen, Pr)
end

function MMI.transform(m::T, (core_res, out_dim, _), X) where T<:Union{SubspaceLDA, BayesianSubspaceLDA}
    # projection of X, XWt is nt x o  where o = out dims
    proj = core_res.projw * view(core_res.projLDA, :, 1:out_dim) #proj is the projection_matrix
    XWt = MMI.matrix(X) * proj
    return MMI.table(XWt, prototype = X)
end

metadata_model(
    BayesianSubspaceLDA,
    input=Table(Continuous),
    target=AbstractVector{<:Finite},
    weights=false,
    output=Table(Continuous),
    descr=BayesianSubspaceLDA_DESCR,
    path="$(PKG).BayesianSubspaceLDA"
)


text:
#############################################################################
# abs.jl
# Absolute value of an expression
# All expressions and atoms are subtpyes of AbstractExpr.
# Please read expressions.jl first.
#############################################################################
import Base.abs, Base.abs2

### Absolute Value

struct AbsAtom <: AbstractExpr
    head::Symbol
    id_hash::UInt64
    children::Tuple{AbstractExpr}
    size::Tuple{Int,Int}

    function AbsAtom(x::AbstractExpr)
        children = (x,)
        return new(:abs, hash(children), children, x.size)
    end
end

function sign(x::AbsAtom)
    return Positive()
end

function monotonicity(x::AbsAtom)
    return (Nondecreasing() * sign(x.children[1]),)
end

function curvature(x::AbsAtom)
    return ConvexVexity()
end

function evaluate(x::AbsAtom)
    return abs.(evaluate(x.children[1]))
end

function conic_form!(x::AbsAtom, unique_conic_forms::UniqueConicForms)
    if !has_conic_form(unique_conic_forms, x)
        c = x.children[1]
        t = Variable(size(c))
        objective = conic_form!(t, unique_conic_forms)
        if sign(x.children[1]) == ComplexSign()
            for i in 1:length(vec(t))
                conic_form!(
                    t[i] >= norm2([real(c[i]); imag(c[i])]),
                    unique_conic_forms,
                )
            end
        else
            conic_form!(c <= t, unique_conic_forms)
            conic_form!(c >= -t, unique_conic_forms)
        end
        cache_conic_form!(unique_conic_forms, x, objective)
    end
    return get_conic_form(unique_conic_forms, x)
end

abs(x::AbstractExpr) = AbsAtom(x)
abs2(x::AbstractExpr) = square(abs(x))


text:
# Memory
# ======

"""
A contiguous memory.

This type works like a `Vector` method.
"""
struct Memory
    ptr::Ptr{UInt8}
    size::UInt
end

function Memory(data::ByteData)
    return Memory(pointer(data), sizeof(data))
end

@inline function Base.length(mem::Memory)
    return mem.size
end

@inline function Base.lastindex(mem::Memory)
    return Int(mem.size)
end

@inline function Base.checkbounds(mem::Memory, i::Integer)
    if !(1 ≤ i ≤ lastindex(mem))
        throw(BoundsError(mem, i))
    end
end

@inline function Base.getindex(mem::Memory, i::Integer)
    @boundscheck checkbounds(mem, i)
    return unsafe_load(mem.ptr, i)
end

@inline function Base.setindex!(mem::Memory, val::UInt8, i::Integer)
    @boundscheck checkbounds(mem, i)
    return unsafe_store!(mem.ptr, val, i)
end


text:
#=
prompt.jl
prompt user to hit a key to continue
=#

export prompt

import Plots # gui, isplotnull

#using REPL #: RadioMenu, request
import REPL
using REPL.TerminalMenus

# global prompt state
prompt_state = :prompt

"""
`prompt()`
prompt user to hit enter to continue, after gui()
"""
function prompt(; gui::Bool=true)
	global prompt_state

	prompt_state == :nodraw && return nothing

	gui && !Plots.isplotnull() && Plots.gui()

	prompt_state == :prompt && wait_for_key()
	return nothing

#=
	list = ["continue*", "run", "draw", "quit", "nodraw"]
	menu = RadioMenu(list) #;, pagesize=3)

	preface = ""
	[name, line] = caller_name
	if isempty(name)
		preface = []
	else
		preface = sprintf('%s %d: ', name, line)
	end

	what = "enter to continue (or select)"

	choice = request(preface * what, menu)
	nothing
=#
end

"""
`function wait_for_key(; prompt=?, io=stdin)`
from:
https://discourse.julialang.org/t/wait-for-a-keypress/20218
"""
function wait_for_key(; prompt::String = "press any key: ", io = stdin)
	setraw!(raw) = ccall(:jl_tty_set_mode, Int32, (Ptr{Cvoid},Int32), io.handle, raw)
	print(io, prompt)
	setraw!(true)
	read(io, 1)
	setraw!(false)
	nothing
end

"""
`prompt(key::symbol)`
set prompt state to one of:
- `:prompt` call gui() if possible then prompt user
- `:draw` call gui() if possible then continue
- `:nodraw` do not call gui(), just continue

`prompt(:test)`
self test
"""
function prompt(key::Symbol)
	global prompt_state

	if key == :state
		return prompt_state
	end

	if key == :test
		prompt(:draw)
		@test prompt(:state) == :draw
		prompt()
		return true
	end

	key ∉ (:prompt, :draw, :nodraw) && throw("prompt $prompt")
	prompt_state = key
end


text:
export ZeContext, status

mutable struct ZeContext
    handle::ze_context_handle_t

    driver::ZeDriver

    function ZeContext(drv::ZeDriver)
        desc_ref = Ref(ze_context_desc_t())
        handle_ref = Ref{ze_context_handle_t}()
        zeContextCreate(drv, desc_ref, handle_ref)
        obj = new(handle_ref[], drv)
        finalizer(obj) do obj
            zeContextDestroy(obj)
        end
        obj
    end
end

Base.unsafe_convert(::Type{ze_context_handle_t}, dev::ZeContext) = dev.handle

status(ctx::ZeContext) = zeContextGetStatus(ctx)


text:
using Printf

function quadroots(x::Real, y::Real, z::Real)
    a, b, c = promote(float(x), y, z)
    if a ≈ 0.0 return [-c / b] end
    Δ = b ^ 2 - 4a * c
    if Δ ≈ 0.0 return [-sqrt(c / a)] end
    if Δ < 0.0 Δ = complex(Δ) end
    d = sqrt(Δ)
    if b < 0.0
        d -= b
        return [d / 2a, 2c / d]
    else
        d = -d - b
        return [2c / d, d / 2a]
    end
end

a = [1, 1, 1.0, 10]
b = [10, 2, -10.0 ^ 9, 1]
c = [1, 1, 1, 1]

for (x, y, z) in zip(a, b, c)
    @printf "The roots of %.2fx² + %.2fx + %.2f\n\tx₀ = (%s)\n" x y z join(round.(quadroots(x, y, z), 2), ", ")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Test

function twoSum(nums::Array{Int64,1}, target::Int64)
	d = Dict()
	for (i, x) in enumerate(nums)
		if get(d, target-x, -1) > 0
			return (get(d, target-x, nothing), i)
		end
		d[x] = i
	end
	(0, 0)
end

@test twoSum([3,2,4], 6) == (2,3)
@test twoSum([3,2,4], 9) == (0,0)

max_stars_repo_name: cdluminate/MyNotes

max_issues_repo_name: cdluminate/MyNotes

max_forks_repo_name: cdluminate/MyNotes