text:
using PedigreeBase
using LinearAlgebra
using Test

include("functions.jl")

@testset "Reading pedigree" begin
   # character file
   pedlist,idtable = read_ped("ped.henderson.char")
   @test pedlist == [0 0 1 1 3 1 5; 0 0 0 2 4 4 6]
   @test idtable == Dict("A"=>1,"B"=>2,"C"=>3,"D"=>4,"E"=>5,"F"=>6,"G"=>7,"0"=>0)
   try
      pedlist = read_ped("ped.henderson.char",integer=true)
      @test false
   catch
      # must be failed
      @test true
   end
   # integer file
   pedlist = read_ped("ped.henderson.int",integer=true)
   @test pedlist == [0 0 1 1 3 1 5; 0 0 0 2 4 4 6]
   # shuffled file
   pedlist,idtable = read_ped("ped.henderson.char.2",order=vec([2 4 5]))
   @test pedlist==[2 0 0 5 2 4 2; 3 0 0 1 0 7 1]
   @test idtable == Dict("A"=>2,"B"=>3,"C"=>5,"D"=>1,"E"=>4,"F"=>7,"G"=>6,"0"=>0)
   pedlist = read_ped("ped.henderson.int.2",order=vec([2 3 1]),integer=true)
   @test pedlist==[0 0 1 1 3 1 5; 0 0 0 2 4 4 6]
end

@testset "Wrong file format" begin
   # wrong file name
   @test_throws ErrorException read_ped("ThereIsNotSuchAFile.")
   # invalid order
   @test_throws BoundsError read_ped("ped.henderson.char.2",order=vec([2 8 5]))
   @test_throws ErrorException read_ped("ped.henderson.char.2",order=vec([2 4 0]))
   # missing columns
   @test_throws BoundsError read_ped("ped.henderson.char.error1")
   @test_throws BoundsError read_ped("ped.henderson.int.error1",integer=true)
   # non-integer data
   @test_throws ErrorException read_ped("ped.henderson.int.error2",integer=true)
   # duplication
   println("The following warning is expected.")
   @test_logs (:warn,"ID: 'A' appears in this file again. The last definision will be used.") read_ped("ped.henderson.char.error3")
   @test_logs (:warn,"ID: '1' appears in this file again. The last definision will be used.") read_ped("ped.henderson.int.error3",integer=true)
end

@testset "Logical error in pedigree" begin
   # okay
   pedlist,idtable = read_ped("ped.henderson.char")
   @test check_ped(pedlist)
   pedlist = read_ped("ped.henderson.int",integer=true)
   @test check_ped(pedlist)
   # mixed sire/dam
   pedlist,idtable = read_ped("ped.henderson.char.error4")
   #@test !check_ped(pedlist)
   @test_logs (:warn,"sire/dam found in the alternative sex: 1") check_ped(pedlist)
   pedlist = read_ped("ped.henderson.int.error4",integer=true)
   @test_logs (:warn,"sire/dam found in the alternative sex: 1") check_ped(pedlist)
   # loops
   pedlist,idtable = read_ped("ped.henderson.char.error5")
   @test_logs (:warn,"pedigree loops found: 3") check_ped(pedlist)
   pedlist = read_ped("ped.henderson.int.error5",integer=true)
   @test_logs (:warn,"pedigree loops found: 3") check_ped(pedlist)
   # progeny precede ancestor
   pedlist,idtable = read_ped("ped.henderson.char.3")
   @test_logs (:warn,"pedigree not chronologically sorted: 1") check_ped(pedlist, parentsfirst=true)
   pedlist = read_ped("ped.henderson.int.3",integer=true)
   @test_logs (:warn,"pedigree not chronologically sorted: 1") check_ped(pedlist, parentsfirst=true)
end

@testset "Permutation of pedigree" begin
   # character pedigree
   pedlist,idtable = read_ped("ped.henderson.char.3")
   perm,invp = find_ped_order(pedlist)
   permute_ped!(invp,pedlist,idtable)
   @test check_ped(pedlist, parentsfirst=true)
   # integer pedigree
   pedlist = read_ped("ped.henderson.int.3", integer=true)
   perm,invp = find_ped_order(pedlist)
   permute_ped!(invp,pedlist)
   @test check_ped(pedlist, parentsfirst=true)
   # make a subset
   pedlist,idtable = read_ped("ped.henderson.char.3")
   perm,invp = find_ped_order(pedlist,vec([7 3]))
   newpedlist,newidtable = extract_ped(invp,pedlist,idtable)
   @test check_ped(newpedlist, parentsfirst=true)
   # the same but different starting animals
   pedlist,idtable = read_ped("ped.henderson.char.3")
   perm,invp = find_ped_order(pedlist,vec([5 2 7]))
   newpedlist,newidtable = extract_ped(invp,pedlist,idtable)
   @test check_ped(newpedlist, parentsfirst=true)
end

@testset "Pedigree with UPG" begin
   pedlist,idtable = read_ped("ped.henderson.char.4")
   upgtable = Dict("x"=>1,"y"=>2)
   pedlist,idtable = set_upg(pedlist,idtable,upgtable)
   perm,invp = find_ped_order(pedlist)
   permute_ped!(invp,pedlist,idtable)
   @test check_ped(pedlist, parentsfirst=true)

   pedlist,idtable = read_ped("ped.henderson.char.4")
   upgtable = Dict("x"=>1,"y"=>2)
   pedlist,idtable = set_upg(pedlist,idtable,upgtable)
   perm,invp = find_ped_order(pedlist,vec([6]))
   newpedlist,newidtable = extract_ped(invp,pedlist,idtable)
   @test check_ped(newpedlist, parentsfirst=true)
end

@testset "Inbreeding coefficients" begin
   # reference
   pedlist1,idtable1 = read_ped("ped.henderson.char")
   f1 = get_inb(pedlist1, check=true)
   # shuffled data
   pedlist2,idtable2 = read_ped("ped.henderson.char.2",order=vec([2 4 5]))
   perm,invp = find_ped_order(pedlist2)
   permute_ped!(invp,pedlist2,idtable2)
   f2 = get_inb(pedlist2, check=true)
   @test f1 ≈ f2
   # reordered data
   pedlist3,idtable3 = read_ped("ped.henderson.char.3")
   perm,invp = find_ped_order(pedlist3)
   permute_ped!(invp,pedlist3,idtable3)
   f3 = get_inb(pedlist3, check=true)
   @test f1 ≈ f3

   # integer case
   pedlist=read_ped("ped.henderson.int",integer=true)
   s = pedlist[1,:]
   d = pedlist[2,:]
   f1a = get_inb(s,d)
   f1b = get_inb(pedlist)
   @test f1a ≈ f1b
end

@testset "Inb/upg code" begin
   # reference
   pedlist=read_ped("ped.henderson.int",integer=true)
   s = pedlist[1,:]
   d = pedlist[2,:]
   f = get_inb(s,d)
   refcode = [1000,1000,1333,2000,2000,2000,2462]
   @test all( get_inbupg_code(s,d,f) .== refcode )
   @test all( get_inbupg_code(pedlist,f) .== refcode )
   n = size(s,1)
   for i=1:n
      @test get_inbupg_code(s[i],d[i],f) == refcode[i]
   end
end

@testset "Inbreeding coefficients with sorting" begin
   # reference
   pedlist1,idtable1 = read_ped("ped.henderson.char")
   f1 = get_inb(pedlist1, check=true)
   # try again Case 3 keeping the original order
   pedlist3,idtable3 = read_ped("ped.henderson.char.3")
   perm,invp = find_ped_order(pedlist3)
   permped3 = copy(pedlist3)
   permute_ped!(invp,permped3)
   f3 = get_inb(permped3, check=true)
   f3 .= f3[invp]
   nerr=0
   for k in sort(collect(keys(idtable1)))
      v=idtable1[k]
      if v>0
         #println((k,v,f1[v]))
         if f1[v]!=f3[idtable3[k]]
            nerr = nerr + 1
         end
      end
   end
   @test nerr==0

   # built-in sorting
   f3s = get_inb(pedlist3, check=true, sorthere=true)
   @test f3s ≈ f3
end

@testset "Numerator relationship matrix" begin
   # reference
   pedlist1,idtable1 = read_ped("ped.henderson.char")
   A1 = get_nrm(pedlist1)
   f1 = get_inb(pedlist1, check=true)
   n = size(A1,1)
   p1,ip1 = sort_by_dict(idtable1,n)
   A1 = A1[p1,p1]
   @test f1 ≈ diag(A1) .- 1.0
   # shuffled data
   pedlist2,idtable2 = read_ped("ped.henderson.char.2",order=vec([2 4 5]))
   perm,invp = find_ped_order(pedlist2)
   permute_ped!(invp,pedlist2,idtable2)
   A2 = get_nrm(pedlist2)
   p2,ip2 = sort_by_dict(idtable2,n)
   A2 = A2[p2,p2]
   @test A1 ≈ A2
   # shuffled with UPGs
   pedlist3,idtable3 = read_ped("ped.henderson.char.4")
   upgtable3 = Dict("x"=>1,"y"=>2)
   pedlist3,idtable3 = set_upg(pedlist3,idtable3,upgtable3)
   perm,invp = find_ped_order(pedlist3)
   permute_ped!(invp,pedlist3,idtable3)
   A3 = get_nrm(pedlist3)
   p3,ip3 = sort_by_dict(idtable3,n)
   A3 = A3[p3,p3]
   @test A1 ≈ A3
end

@testset "A-inverse" begin
   # reference
   pedlist1,idtable1 = read_ped("ped.henderson.char")
   A1 = get_nrm(pedlist1)
   f1 = get_inb(pedlist1, check=true)
   n = size(A1,1)
   p1,ip1 = sort_by_dict(idtable1,n)
   A1 = A1[p1,p1]
   Ainv1 = Matrix( get_nrminv(pedlist1, f1) )
   Ainv1 = Ainv1[p1,p1]
   @test Ainv1 ≈ inv(A1)
   # shuffled data
   pedlist2,idtable2 = read_ped("ped.henderson.char.2",order=vec([2 4 5]))
   perm,invp = find_ped_order(pedlist2)
   permute_ped!(invp,pedlist2,idtable2)
   A2 = get_nrm(pedlist2)
   f2 = get_inb(pedlist2, check=true)
   p2,ip2 = sort_by_dict(idtable2,n)
   A2 = A2[p2,p2]
   Ainv2 = Matrix( get_nrminv(pedlist2, f=f2) )
   Ainv2 = Ainv2[p2,p2]
   @test Ainv2 ≈ inv(A2)
   # shuffled with UPGs
   pedlist3,idtable3 = read_ped("ped.henderson.char.4")
   upgtable3 = Dict("x"=>1,"y"=>2)
   pedlist3,idtable3 = set_upg(pedlist3,idtable3,upgtable3)
   perm,invp = find_ped_order(pedlist3)
   permute_ped!(invp,pedlist3,idtable3)
   A3 = get_nrm(pedlist3)
   f3 = get_inb(pedlist3, check=true)
   p3,ip3 = sort_by_dict(idtable3,n)
   A3 = A3[p3,p3]
   Ainv3 = Matrix( get_nrminv(pedlist3, f=f3) )
   Ainv3 = Ainv3[p3,p3]
   @test Ainv3 ≈ inv(A3)
end

@testset "A-inverse with UPG" begin
   pedlist=read_ped("ped.quaas.int",integer=true,hasupg=true)
   f = get_inb(pedlist)
   Ai = Matrix(get_nrminv(pedlist,f=f))
   M = zeros(6,6)
   M[1,1]=4/3; M[1,2]=-2/3; M[1,5]=-1/2; M[1,6]=-1/6;
   M[2,2]=11/6; M[2,3]=1/2; M[2,4]=-1; M[2,6]=-2/3;
   M[3,3]=3/2; M[3,4]=-1; M[3,5]=-1/2; M[3,6]=-1/2;
   M[4,4]=2; M[5,5]=1/2; M[5,6]=1/2; M[6,6]=5/6;
   #M = M + M' - diagm(diag(M))
   M = M + M' - Matrix(Diagonal(diag(M)))
   @test M ≈ Ai
   upgs = get_upg_range(pedlist)
   @test upgs == 5:6
end

@testset "Q matrix" begin
   # regular process
   pedlist,idtable = read_ped("ped.quaas.char")
   upgtable = Dict("a"=>1,"b"=>2,"c"=>2,"d"=>1,"e"=>2)
   pedlist,idtable = set_upg(pedlist,idtable,upgtable)
   perm,invp = find_ped_order(pedlist)
   newpedlist,newidtable = extract_ped(invp,pedlist,idtable)
   P1 = get_pm(newpedlist)
   Q1 = get_qm(newpedlist)
   cp = [newidtable["a"]-4, newidtable["b"]-4]
   rp,rip = sort_by_dict(newidtable,size(pedlist,2))

   # sorted by the original ID
   pedlist,idtable = read_ped("ped.quaas.char")
   @test_logs (:warn,"pedigree with 0 index; not all animals assined to UPG") get_pm(pedlist)
   P0 = get_pm(pedlist, verbose=false)
   p,ip = sort_by_dict(idtable,size(pedlist,2))
   Pb = P0[p,p][1:4,5:9]
   P = P0[p,p][1:4,1:4]
   Qb = get_qb(upgtable)
   Q0 = inv(Matrix(1.0I,4,4)-P)*Pb*Qb
   @test inv(Matrix(1.0I,4,4)-P) ≈ Matrix(1.0I,4,4) + P + P*P
   @test Q0 ≈ Q1[rp,cp]
end

@testset "Subset pedigree" begin
   #           1  2  3  4  5  6  7  8  9 10 11  12 13 14 15
   pedlist = [ 0  0  0  0  0  0  2  1  2  7  7  11 11  9 11;
               0  0  0  0  0  0  5  4  3  6  4   8 10 13 10]
   idlist = [14]
   perm,invp,subpedlist = subset_ped(pedlist,idlist)
   f = get_inb(pedlist)
   subf = get_inb(subpedlist)
   @test all( map(x->f[perm[x]],findall(perm.>0)) .== subf )
   @test all( map(x->subf[invp[x]],findall(invp.>0)) .== subf )
end

@testset "MGS pedigree" begin
   # Henderson (1975)
   # for all pedigree including females
   #           1  2 X_3 Y_4 Z_5 3_6 4_7 5_8 6_9
   pedlist = [ 0  0  1   2   2   2   1   0   7
               0  0  0   0   0   3   4   5   0]
   mgslist = mgs_ped(pedlist)
   ref_mgs = [ 0  0  1   2   2   2   1   0   7
               0  0  0   0   0   1   2   2   0]
   @test all(mgslist .== ref_mgs)

   # for all pedigree including males only
   #           1  2 X_3 Y_4 Z_5 3_6 4_7 5_8 6_9
   pedlist = [ 0  0  1   2   2   2   1   0   7
               0  0  0   0   0   3   4   5   0]
   males =   [true, true, false, false, false, true, true, true, true]
   mgslist = mgs_ped(pedlist, males)
   ref_mgs = [ 0  0  0   0   0   2   1   0   7
               0  0  0   0   0   1   2   2   0]
   @test all(mgslist .== ref_mgs)

   Ainv = get_mgsnrminv(mgslist, males)
   idx = [1,2,6,7,8,9]
   ref_Ainv = [
      21.8182 5.4545  -5.4545  -10.9091    0 0
      5.4545  22.8182 -10.9091 -5.4545  -4.0 0
     -5.4545 -10.9091  21.8182  0        0   0
     -10.9091 -5.4545  0        26.8182  0  -10.0
     0        -4.0     0        0       16.0 0
     0         0       0        -10.0   0    20.0
   ]
   @test isapprox(ref_Ainv,Ainv[idx,idx]*15,atol=1e-3)
end


text:
j1 = Jump(0.5,"a","b")
j2 = Jump(0.5,"a","b")
j3 = Jump(0.6,"a","b")
@test j1.t==0.5
@test j2.from == "a"
@test j2.to == "b"
@test j1==j2
@test j1!=j3

j1 = DurationJump(0.5,"a","b",1.3)
j2 = DurationJump(0.5,"a","b",1.3)
j3 = DurationJump(0.5,"a","b",2.2)
@test j1.t==0.5
@test j2.from == "a"
@test j2.to == "b"
@test j1.δ == 1.3
@test duration(j2) == 1.3
@test j1==j2
@test j1!=j3

ext: jl


text:
@testset "743.network-delay-time.jl" begin
    @test network_delay_time([[2,1,1],[2,3,1],[3,4,1]], 4, 2) == 2
    @test network_delay_time([[1,2,1]], 2, 1) == 1
    @test network_delay_time([[1,2,1]], 2, 2) == -1
end


text:
include("../console.jl")

arr = [21, -22, [23, 24], "owl", Dict("eagle" => [25])];
map = Dict("cat" => [31, 32], "dog" => 33, 34 => "fox")
console.log("Julia:", 1337, arr, map)

max_forks_repo_name: lgarron/console


text:
"""
Numeric representation of SAT problem in DIMACS-CNF form

Input:
- `vars` is a list of symbolic Boolean variables
- `clauses` are integer vectors; negative numbers indicate negated variables
"""

"""SAT problem represented numerically in DIMACS-CNF form:
clauses are integer vectors with negative numbers indicating negated literals
"""
struct SATProblem
    num_variables::Int
    clauses::Vector{Vector{Int}}
end

max_var(clause) = maximum(abs.(clause))

"Constructor which calculates the number of variables automatically"
function SATProblem(clauses::Vector{Vector{Int}})
    num_variables = maximum(max_var.(clauses))

    return SATProblem(num_variables, clauses)
end

push!(p::SATProblem, clause) = push!(p.clauses, clause)
ban!(p::SATProblem, clause) = push!(p, .-(clause))

function all_solutions(p::SATProblem)

    solns = []
    status = :sat

    while true
        status, results = solve(p)
        status == :unsat && break

        push!(solns, results)
        ban!(p, results)

    end

    return solns
end


text:
#const CIS_LUT = SVector{32}(complex.(floor.(Int8, cos.(2π .* (0:31) ./ 32) .* (1 << 4)), floor.(Int8, sin.(2π .* (0:31) ./ 32) .* (1 << 4))))
#cis_fast(x) = @inbounds CIS_LUT[x & 31 + 1]
"""
$(SIGNATURES)

Get sine approximation very fast. It is valid between -π and π.
```julia-repl
julia> sin_vfast(π / 4)
```
"""
@inline function sin_vfast(x::T) where T <: Real
    B = T(4 / π)
    C = T(-4 / (π * π))
    B * x + C * x * abs(x)
end

"""
$(SIGNATURES)

Get cosine approximation very fast. It is valid between -π and π.
```julia-repl
julia> cos_vfast(π / 4)
```
"""
@inline function cos_vfast(x::T) where T <: Real
    x += T(π / 2)
    x -= (x > T(π)) * T(2π)
    sin_vfast(x)
end

"""
$(SIGNATURES)

Get sine approximation fast. It is valid between -π and π.
```julia-repl
julia> sin_fast(π / 4)
```
"""
@inline function sin_fast(x::T) where T <: Real
    P = T(0.225)
    y = sin_vfast(x)
    P * (y * abs(y) - y) + y
end

"""
$(SIGNATURES)

Get cosine approximation fast. It is valid between -π and π.
```julia-repl
julia> cos_fast(π / 4)
```
"""
@inline function cos_fast(x::T) where T <: Real
    x += T(π / 2)
    x -= (x > T(π)) * T(2π)
    sin_fast(x)
end

"""
$(SIGNATURES)

Get \exp(iz) approximation fast. It is valid between -π and π.
```julia-repl
julia> cos_fast(π / 4)
```
"""
function cis_fast(x)
    complex(cos_fast(x), sin_fast(x))
end

"""
$(SIGNATURES)

Get \exp(iz) approximation very fast. It is valid between -π and π.
```julia-repl
julia> cis_vfast(π / 4)
```
"""
function cis_vfast(x)
    complex(cos_vfast(x), sin_vfast(x))
end

"""
$(SIGNATURES)

Get \exp(iz) approximation fast. It is valid between -π and π.
```julia-repl
julia> get_carrier_fast_unsafe(π / 4)
```
"""
@inline function get_carrier_fast_unsafe(x)
    cis_fast(x)
end

"""
$(SIGNATURES)

Get \exp(iz) approximation very fast. It is valid between -π and π.
```julia-repl
julia> get_carrier_vfast_unsafe(π / 4)
```
"""
@inline function get_carrier_vfast_unsafe(x)
    cis_vfast(x)
end

"""
$(SIGNATURES)

Fixed point carrier phases
"""
function fpcarrier_phases!(
    phases::Vector{T},
    carrier_frequency,
    sample_frequency,
    start_phase::AbstractFloat;
    start_sample::Integer = 1,
    num_samples::Integer = length(phases),
    bits::Val{N} = Val(5)
) where {T <: Integer, N}
    n = N + 2
    fixed_point = 32 - n - 2
    delta = floor(Int32, carrier_frequency * 1 << (fixed_point + n) / sample_frequency)
    fixed_point_start_phase = floor(Int32, start_phase * 1 << (fixed_point + n))
    fixed_point_phase = fixed_point_start_phase - delta
    @inbounds for i = start_sample:num_samples + start_sample - 1
        fixed_point_phase = fixed_point_phase + delta
        phases[i] = T(fixed_point_phase >> fixed_point)
    end
    phases
end

"""
$(SIGNATURES)

Fixed point carrier
"""
function fpcarrier!(
    carrier_sin::VT,
    carrier_cos::VT,
    phases::VT;
    start_sample::Integer = 1,
    num_samples::Integer = length(phases),
    bits::Val{N} = Val(5)
) where {VT <: Vector{Int16}, N}
    @avx unroll = 6 for i = start_sample:num_samples + start_sample - 1
        carrier_sin[i] = fpsin(phases[i], bits)
        carrier_cos[i] = fpcos(phases[i], bits)
    end
end

"""
$(SIGNATURES)

Fixed point carrier
"""
function fpcarrier!(
    carrier::StructArray{Complex{Int16}},
    carrier_frequency,
    sample_frequency,
    start_phase::AbstractFloat;
    start_sample::Integer = 1,
    num_samples::Integer = length(carrier),
    bits::Val{N} = Val(5)
) where N
    fpcarrier_phases!(
        carrier.re,
        carrier_frequency,
        sample_frequency,
        start_phase,
        start_sample = start_sample,
        num_samples = num_samples,
        bits = bits
    )
    fpcarrier!(
        carrier.im,
        carrier.re,
        carrier.re,
        start_sample = start_sample,
        num_samples = num_samples,
        bits = bits
    )
    carrier
end


text:
# There's not much to do as you slowly descend to the bottom of the ocean. The
# submarine computer challenges you to a nice game of Dirac Dice.
#
# This game consists of a single die, two pawns, and a game board with a
# circular track containing ten spaces marked 1 through 10 clockwise. Each
# player's starting space is chosen randomly (your puzzle input). Player 1 goes
# first.
#
# Players take turns moving. On each player's turn, the player rolls the die
# three times and adds up the results. Then, the player moves their pawn that
# many times forward around the track (that is, moving clockwise on spaces in
# order of increasing value, wrapping back around to 1 after 10). So, if a
# player is on space 7 and they roll 2, 2, and 1, they would move forward 5
# times, to spaces 8, 9, 10, 1, and finally stopping on 2.
#
# After each player moves, they increase their score by the value of the space
# their pawn stopped on. Players' scores start at 0. So, if the first player
# starts on space 7 and rolls a total of 5, they would stop on space 2 and add
# 2 to their score (for a total score of 2). The game immediately ends as a win
# for any player whose score reaches at least 1000.
#
# Since the first game is a practice game, the submarine opens a compartment
# labeled deterministic dice and a 100-sided die falls out. This die always
# rolls 1 first, then 2, then 3, and so on up to 100, after which it starts
# over at 1 again. Play using this die.
#
# For example, given these starting positions:
#
# Player 1 starting position: 4
# Player 2 starting position: 8
# This is how the game would go:
#
# Player 1 rolls 1+2+3 and moves to space 10 for a total score of 10.
# Player 2 rolls 4+5+6 and moves to space 3 for a total score of 3.
# Player 1 rolls 7+8+9 and moves to space 4 for a total score of 14.
# Player 2 rolls 10+11+12 and moves to space 6 for a total score of 9.
# Player 1 rolls 13+14+15 and moves to space 6 for a total score of 20.
# Player 2 rolls 16+17+18 and moves to space 7 for a total score of 16.
# Player 1 rolls 19+20+21 and moves to space 6 for a total score of 26.
# Player 2 rolls 22+23+24 and moves to space 6 for a total score of 22.
# ...after many turns...
#
# Player 2 rolls 82+83+84 and moves to space 6 for a total score of 742.
# Player 1 rolls 85+86+87 and moves to space 4 for a total score of 990.
# Player 2 rolls 88+89+90 and moves to space 3 for a total score of 745.
# Player 1 rolls 91+92+93 and moves to space 10 for a final score, 1000.
#
# Since player 1 has at least 1000 points, player 1 wins and the game ends. At
# this point, the losing player had 745 points and the die had been rolled a
# total of 993 times; 745 * 993 = 739785.
#
# Play a practice game using the deterministic 100-sided die. The moment either
# player wins, what do you get if you multiply the score of the losing player
# by the number of times the die was rolled during the game?

# read input
positions = map(l -> l[end] - '0', eachline("input.txt"))

# do eet
rolls = Iterators.Stateful(Iterators.cycle(1:100))
scores = zeros(Int, length(positions))
turn = 0
while all(s -> s < 1000, scores)
  player = (turn & 1) + 1
  position = positions[player] = (positions[player] + sum(Iterators.take(rolls, 3)) - 1) % 10 + 1
  scores[player] += position
  global turn += 1
end

println("Result: ", scores[scores .< 1000][1] * turn * 3)


text:
# ---
# title: 324. Wiggle Sort II
# id: problem324
# author: Indigo
# date: 2021-06-10
# difficulty: Medium
# categories: Sort
# link: <https://leetcode.com/problems/wiggle-sort-ii/description/>
# hidden: true
# ---
#
# Given an unsorted array `nums`, reorder it such that `nums[0] < nums[1] >
# nums[2] < nums[3]...`.
#
# **Example 1:**
#
#
#
#     Input:nums = [1, 5, 1, 1, 6, 4]
#     Output: One possible answer is [1, 4, 1, 5, 1, 6].
#
# **Example 2:**
#
#
#
#     Input:nums = [1, 3, 2, 2, 3, 1]
#     Output: One possible answer is [2, 3, 1, 3, 1, 2].
#
# **Note:**
# You may assume all input has valid answer.
#
# **Follow Up:**
# Can you do it in O(n) time and/or in-place with O(1) extra space?
#
#
## @lc code=start
using LeetCode

function wiggle_sort!(nums::Vector{Int})
    len = length(nums)
    md = len ÷ 2
    partialsort!(nums, md)
    nums[1], nums[md] = nums[md], nums[1]
    i, j = 2, len
    iseven(len) && (j -= 1)
    while i <= j
        nums[i], nums[j] = nums[j], nums[i]
        i += 2
        j -= 2
    end
    return nums
end
## @lc code=end


text:
using AbstractSphericalHarmonics
using LinearAlgebra: norm
using Random
using StaticArrays
using Test

const bitsign = AbstractSphericalHarmonics.bitsign

chop(x) = abs2(x) < 100eps(x) ? zero(x) : x
chop(x::Complex) = Complex(chop(real(x)), chop(imag(x)))
chop(x::SArray) = chop.(x)

function integrate(f::AbstractMatrix, g::AbstractMatrix, lmax::Integer)
    sz = ash_grid_size(lmax)
    @assert size(f) == size(g) == sz

    s = zero(eltype(f)) * zero(eltype(g)) * zero(Float64)
    for pt in CartesianIndices(sz)
        theta, phi = ash_point_coord(pt, lmax)
        dtheta, dphi = ash_point_delta(pt, lmax)
        s += conj(f[pt]) * g[pt] * sin(theta) * dtheta * dphi
    end

    return s
end

################################################################################

# Half angle formulae:
#     sin(θ/2)^2 = (1-cos(θ))/2
#     cos(θ/2)^2 = (1+cos(θ))/2

# These sYlm are taken from Wikipedia and black-holes.org:
# <https://en.wikipedia.org/wiki/Spherical_harmonics> and
# <https://en.wikipedia.org/wiki/Spin-weighted_spherical_harmonics>
# <https://www.black-holes.org/SpinWeightedSphericalHarmonics.nb>
function sYlm(s::Integer, l::Integer, m::Integer, θ::Real, ϕ::Real)
    @assert abs(s) ≤ l
    @assert -l ≤ m ≤ l

    # Parity:
    s < 0 && return bitsign(s + m) * conj(sYlm(-s, l, -m, θ, ϕ))

    (s, l, m) == (0, 0, 0) && return sqrt(1 / 4π)
    (s, l, m) == (0, 1, -1) && return sqrt(3 / 8π) * sin(θ) * cis(-ϕ)
    (s, l, m) == (0, 1, 0) && return sqrt(3 / 4π) * cos(θ)
    (s, l, m) == (0, 1, +1) && return -sqrt(3 / 8π) * sin(θ) * cis(ϕ)
    (s, l, m) == (0, 2, -2) && return sqrt(15 / 2π) * cos(θ / 2)^2 * sin(θ / 2)^2 * cis(-2ϕ)
    (s, l, m) == (0, 2, -1) && return -sqrt(15 / 2π) * cos(θ / 2) * sin(θ / 2) * (-cos(θ / 2)^2 + sin(θ / 2)^2) * cis(-ϕ)
    (s, l, m) == (0, 2, 0) && return sqrt(5 / 4π) * (cos(θ / 2)^4 - 4 * cos(θ / 2)^2 * sin(θ / 2)^2 + sin(θ / 2)^4)
    (s, l, m) == (0, 2, +1) && return -sqrt(15 / 2π) * cos(θ / 2) * sin(θ / 2) * (cos(θ / 2)^2 - sin(θ / 2)^2) * cis(ϕ)
    (s, l, m) == (0, 2, +2) && return sqrt(15 / 2π) * cos(θ / 2)^2 * sin(θ / 2)^2 * cis(2ϕ)

    (s, l, m) == (+1, 1, -1) && return -sqrt(3 / 16π) * (1 + cos(θ)) * cis(-ϕ)
    (s, l, m) == (+1, 1, 0) && return sqrt(3 / 8π) * sin(θ)
    (s, l, m) == (+1, 1, +1) && return -sqrt(3 / 16π) * (1 - cos(θ)) * cis(ϕ)
    (s, l, m) == (+1, 2, -2) && return -sqrt(5 / π) * cos(θ / 2)^3 * sin(θ / 2) * cis(-2ϕ)
    (s, l, m) == (+1, 2, -1) && return -sqrt(5 / 4π) * cos(θ / 2)^2 * (cos(θ / 2)^2 - 3 * sin(θ / 2)^2) * cis(-ϕ)
    (s, l, m) == (+1, 2, 0) && return sqrt(15 / 2π) * cos(θ / 2) * sin(θ / 2) * (cos(θ / 2)^2 - sin(θ / 2)^2)
    (s, l, m) == (+1, 2, +1) && return sqrt(5 / 4π) * sin(θ / 2)^2 * (-3 * cos(θ / 2)^2 + sin(θ / 2)^2) * cis(ϕ)
    (s, l, m) == (+1, 2, +2) && return sqrt(5 / π) * cos(θ / 2) * sin(θ / 2)^3 * cis(2ϕ)

    (s, l, m) == (+2, 2, -2) && return sqrt(5 / 4π) * cos(θ / 2)^4 * cis(-2ϕ)
    (s, l, m) == (+2, 2, -1) && return -sqrt(5 / π) * cos(θ / 2)^3 * sin(θ / 2) * cis(-ϕ)
    (s, l, m) == (+2, 2, 0) && return sqrt(15 / 2π) * cos(θ / 2)^2 * sin(θ / 2)^2
    (s, l, m) == (+2, 2, +1) && return -sqrt(5 / π) * cos(θ / 2) * sin(θ / 2)^3 * cis(ϕ)
    (s, l, m) == (+2, 2, +2) && return sqrt(5 / 4π) * sin(θ / 2)^4 * cis(2ϕ)

    @assert false
end

# ð F = - (sin θ)^s (∂_θ + i / sin(θ) ∂_ϕ) (sin θ)^-s F
# (Calculated manually from sYlm above)
function ðsYlm(s::Integer, l::Integer, m::Integer, θ::Real, ϕ::Real)
    @assert abs(s) ≤ l
    @assert -l ≤ m ≤ l

    (s, l, m) == (0, 0, 0) && return 0
    (s, l, m) == (0, 1, -1) && return -sqrt(3 / 8π) * (1 + cos(θ)) * cis(-ϕ)
    (s, l, m) == (0, 1, 0) && return sqrt(3 / 4π) * sin(θ)
    (s, l, m) == (0, 1, +1) && return -sqrt(3 / 8π) * (1 - cos(θ)) * cis(ϕ)
    (s, l, m) == (0, 2, -2) && return -sqrt(15 / 8π) * (1 + cos(θ)) * sin(θ) * cis(-2ϕ)
    (s, l, m) == (0, 2, -1) && return -sqrt(15 / 2π) * cos(θ / 2)^2 * (-1 + 2 * cos(θ)) * cis(-ϕ)
    (s, l, m) == (0, 2, 0) && return sqrt(45 / 16π) * sin(2θ)
    (s, l, m) == (0, 2, +1) && return -sqrt(15 / 2π) * (1 + 2 * cos(θ)) * sin(θ / 2)^2 * cis(ϕ)
    (s, l, m) == (0, 2, +2) && return -sqrt(15 / 8π) * (-1 + cos(θ)) * sin(θ) * cis(2ϕ)

    (s, l, m) == (+1, 1, -1) && return 0
    (s, l, m) == (+1, 1, 0) && return 0
    (s, l, m) == (+1, 1, +1) && return 0
    (s, l, m) == (+1, 2, -2) && return sqrt(5 / π) * cos(θ / 2)^4 * cis(-2ϕ)
    (s, l, m) == (+1, 2, -1) && return -sqrt(5 / π) * cos(θ / 2)^2 * sin(θ) * cis(-ϕ)
    (s, l, m) == (+1, 2, 0) && return sqrt(15 / 8π) * sin(θ)^2
    (s, l, m) == (+1, 2, +1) && return -sqrt(5 / π) * sin(θ / 2)^2 * sin(θ) * cis(ϕ)
    (s, l, m) == (+1, 2, +2) && return sqrt(5 / π) * sin(θ / 2)^4 * cis(2ϕ)

    (s, l, m) == (+2, 2, -2) && return 0
    (s, l, m) == (+2, 2, -1) && return 0
    (s, l, m) == (+2, 2, 0) && return 0
    (s, l, m) == (+2, 2, +1) && return 0
    (s, l, m) == (+2, 2, +2) && return 0

    @assert false
end

# ð̄ F = - (sin θ)^-s (∂_θ - i / sin(θ) ∂_ϕ) (sin θ)^s F
# (Calculated manually from sYlm above)
function ð̄sYlm(s::Integer, l::Integer, m::Integer, θ::Real, ϕ::Real)
    @assert abs(s) ≤ l
    @assert -l ≤ m ≤ l

    (s, l, m) == (0, 0, 0) && return 0
    (s, l, m) == (0, 1, -1) && return sqrt(3 / 8π) * (1 - cos(θ)) * cis(-ϕ)
    (s, l, m) == (0, 1, 0) && return sqrt(3 / 4π) * sin(θ)
    (s, l, m) == (0, 1, +1) && return sqrt(3 / 8π) * (1 + cos(θ)) * cis(ϕ)
    (s, l, m) == (0, 1, +1) && return -sqrt(3 / 8π) * (1 - cos(θ)) * cis(ϕ)
    (s, l, m) == (0, 2, -2) && return -sqrt(15 / 8π) * (-1 + cos(θ)) * sin(θ) * cis(-2ϕ)
    (s, l, m) == (0, 2, -1) && return sqrt(15 / 2π) * (1 + 2 * cos(θ)) * sin(θ / 2)^2 * cis(-ϕ)
    (s, l, m) == (0, 2, 0) && return sqrt(45 / 16π) * sin(2θ)
    (s, l, m) == (0, 2, +1) && return sqrt(15 / 2π) * cos(θ / 2)^2 * (-1 + 2 * cos(θ)) * cis(ϕ)
    (s, l, m) == (0, 2, +2) && return -sqrt(15 / 8π) * (1 + cos(θ)) * sin(θ) * cis(2ϕ)

    (s, l, m) == (+1, 1, -1) && return -sqrt(3 / 4π) * sin(θ) * cis(-ϕ)
    (s, l, m) == (+1, 1, 0) && return -sqrt(3 / 2π) * cos(θ)
    (s, l, m) == (+1, 1, +1) && return sqrt(3 / 4π) * sin(θ) * cis(ϕ)
    (s, l, m) == (+1, 2, -2) && return -sqrt(45 / 16π) * sin(θ)^2 * cis(-2ϕ)
    (s, l, m) == (+1, 2, -1) && return -sqrt(45 / 16π) * sin(2θ) * cis(-ϕ)
    (s, l, m) == (+1, 2, 0) && return -sqrt(15 / 32π) * (1 + 3 * cos(2θ))
    (s, l, m) == (+1, 2, +1) && return sqrt(45 / 16π) * sin(2θ) * cis(ϕ)
    (s, l, m) == (+1, 2, +2) && return -sqrt(45 / 16π) * sin(θ)^2 * cis(2ϕ)

    (s, l, m) == (+2, 2, -2) && return sqrt(5 / 4π) * (1 + cos(θ)) * sin(θ) * cis(-2ϕ)
    (s, l, m) == (+2, 2, -1) && return sqrt(5 / 4π) * (cos(θ) + cos(2θ)) * cis(-ϕ)
    (s, l, m) == (+2, 2, 0) && return -sqrt(15 / 2π) * cos(θ) * sin(θ)
    (s, l, m) == (+2, 2, +1) && return sqrt(5 / 4π) * (cos(θ) - cos(2θ)) * cis(ϕ)
    (s, l, m) == (+2, 2, +2) && return -sqrt(5 / π) * sin(θ / 2)^2 * sin(θ) * cis(2ϕ)

    @assert false
end

function rand_tensor(::Val{D}, ::Type{T}, lmax::Int) where {D,T}
    sz = ash_grid_size(lmax)
    Dims = Tuple{[2 for d in 1:D]...}
    f = randn(SArray{Dims,T}, sz)
    return f
end

function const_tensor(::Val{D}, ::Type{T}, lmax::Int) where {D,T}
    sz = ash_grid_size(lmax)
    Dims = Tuple{[2 for d in 1:D]...}
    α = randn(SArray{Dims,T})
    f = fill(α, sz)
    return f
end

################################################################################

@testset "Mode indices" begin
    for lmax in 0:20
        nmodes = ash_nmodes(lmax)
        for s in -4:4, l in abs(s):lmax, m in (-l):l
            ind = ash_mode_index(s, l, m, lmax)
            @test length(nmodes) ≡ length(ind)
            @test all(1 ≤ ind[d] ≤ nmodes[d] for d in 1:length(nmodes))
            l′, m′ = ash_mode_numbers(s, ind, lmax)
            @test l′ == l && m′ == m
        end
    end
end

Random.seed!(100)
modes = [
    (name="(s=$s,l=$l,m=$m)", spin=s, el=l, m=m, fun=(θ, ϕ) -> sYlm(s, l, m, θ, ϕ), modes=(l′, m′) -> l′ == l && m′ == m) for
    s in -2:+2 for l in abs(s):2 for m in (-l):l
]
@testset "Simple transforms: $(mode.name)" for mode in modes
    for lmax in (mode.el):20
        sz = ash_grid_size(lmax)
        f = Array{Complex{Float64}}(undef, sz)
        for ij in CartesianIndices(sz)
            θ, ϕ = ash_point_coord(ij, lmax)
            f[ij] = mode.fun(θ, ϕ)
        end
        # function setvalue(ij::CartesianIndex{2})
        #     θ, ϕ = ash_point_coord(ij, lmax)
        #     return Complex{Float64}(mode.fun(θ, ϕ))
        # end
        # f = map(setvalue, CartesianIndices(sz))

        flm = ash_transform(f, mode.spin, lmax)
        @test all(
            isapprox(flm[ash_mode_index(mode.spin, l, m, lmax)], mode.modes(l, m); atol=100eps()) for l in abs(mode.spin):lmax for
            m in (-l):l
        )

        f′ = ash_evaluate(flm, mode.spin, lmax)
        @test isapprox(f′, f; atol=1000eps())
    end
end

Random.seed!(100)
@testset "Parity" begin
    for iter in 1:20
        lmax = rand(0:100)
        sz = ash_grid_size(lmax)
        nmodes = ash_nmodes(lmax)

        s = rand(-4:4)
        abs(s) ≤ lmax || continue
        l = rand(abs(s):lmax)
        m = rand((-l):l)

        flm = zeros(Complex{Float64}, nmodes)
        flm[ash_mode_index(s, l, m, lmax)] = 1

        flm′ = zeros(Complex{Float64}, nmodes)
        flm′[ash_mode_index(-s, l, -m, lmax)] = 1

        f = ash_evaluate(flm, s, lmax)
        f′ = ash_evaluate(flm′, -s, lmax)

        # Phase: conj(sYlm) = (-1)^(s+m) (-s)Yl(-m)
        @test conj(f) ≈ bitsign(s + m) * f′
    end
end

Random.seed!(100)
@testset "Linearity" begin
    for iter in 1:20
        lmax = rand(0:100)
        sz = ash_grid_size(lmax)
        nmodes = ash_nmodes(lmax)

        spin = rand(-4:4)

        f = randn(Complex{Float64}, sz)
        g = randn(Complex{Float64}, sz)
        α = randn(Complex{Float64})
        h = f + α * g

        flm = ash_transform(f, spin, lmax)
        glm = ash_transform(g, spin, lmax)
        hlm = ash_transform(h, spin, lmax)

        if !(flm + α * glm ≈ hlm)
            @show iter lmax sz nmodes spin
            @show α
            @show any(isnan, f) any(isnan, g) any(isnan, h)
            @show any(isnan, flm) any(isnan, glm) any(isnan, hlm)
            @show f[1] g[1] h[1]
            @show flm[1] glm[1] hlm[1]
        end
        @test flm + α * glm ≈ hlm

        flm = randn(Complex{Float64}, nmodes)
        glm = randn(Complex{Float64}, nmodes)
        hlm = flm + α * glm

        f = ash_evaluate(flm, spin, lmax)
        g = ash_evaluate(glm, spin, lmax)
        h = ash_evaluate(hlm, spin, lmax)

        @test f + α * g ≈ h
    end
end

Random.seed!(100)
@testset "Orthonormality" begin
    for iter in 1:20
        lmax = rand(0:100)
        nmodes = ash_nmodes(lmax)

        smax = min(4, lmax ÷ 2)
        spin = rand((-smax):smax)

        flm = zeros(Complex{Float64}, nmodes)
        glm = zeros(Complex{Float64}, nmodes)

        # The produce will have lh = lf + lg
        lf = rand(abs(spin):(lmax ÷ 2))
        mf = rand((-lf):lf)
        lg = rand(abs(spin):(lmax - lf))
        mg = rand((-lg):lg)

        flm[ash_mode_index(spin, lf, mf, lmax)] = 1
        glm[ash_mode_index(spin, lg, mg, lmax)] = 1

        f = ash_evaluate(flm, spin, lmax)
        g = ash_evaluate(glm, spin, lmax)

        @test isapprox(integrate(f, f, lmax), 1; atol=20 / lmax^2)
        @test isapprox(integrate(f, g, lmax), (lf == lg) * (mf == mg); atol=1 / lmax^2)

        h = conj(f) .* f
        hlm = ash_transform(h, 0, lmax)
        @test isapprox(hlm[ash_mode_index(0, 0, 0, lmax)], sqrt(1 / 4π); atol=sqrt(eps()))

        h = conj(f) .* g
        hlm = ash_transform(h, 0, lmax)
        @test isapprox(hlm[ash_mode_index(0, 0, 0, lmax)], (lf == lg) * (mf == mg) * sqrt(1 / 4π); atol=sqrt(eps()))
    end
end

Random.seed!(100)
modes = [
    (
        name="(s=$s,l=$l,m=$m)",
        spin=s,
        el=l,
        fun=(θ, ϕ) -> sYlm(s, l, m, θ, ϕ),
        ðfun=(θ, ϕ) -> ðsYlm(s, l, m, θ, ϕ),
        ð̄fun=(θ, ϕ) -> ð̄sYlm(s, l, m, θ, ϕ),
    ) for s in 0:+2 for l in abs(s):2 for m in (-l):l
]
@testset "Simple derivatives (eth, eth-bar): $(mode.name)" for mode in modes
    for lmax in (mode.el):20
        sz = ash_grid_size(lmax)
        f = Array{Complex{Float64}}(undef, sz)
        ðf₀ = Array{Complex{Float64}}(undef, sz)
        ð̄f₀ = Array{Complex{Float64}}(undef, sz)
        for ij in CartesianIndices(sz)
            θ, ϕ = ash_point_coord(ij, lmax)
            f[ij] = mode.fun(θ, ϕ)
            ðf₀[ij] = mode.ðfun(θ, ϕ)
            ð̄f₀[ij] = mode.ð̄fun(θ, ϕ)
        end

        flm = ash_transform(f, mode.spin, lmax)

        ðflm = ash_eth(flm, mode.spin, lmax)
        ðf = ash_evaluate(ðflm, mode.spin + 1, lmax)
        @test isapprox(ðf, ðf₀; atol=10000eps())

        ð̄flm = ash_ethbar(flm, mode.spin, lmax)
        ð̄f = ash_evaluate(ð̄flm, mode.spin - 1, lmax)
        @test isapprox(ð̄f, ð̄f₀; atol=10000eps())
    end
end

Random.seed!(100)
@testset "Eigenvectors of Laplacian" begin
    for iter in 1:20
        lmax = rand(0:100)
        nmodes = ash_nmodes(lmax)

        flm = zeros(Complex{Float64}, nmodes)

        l = rand(0:lmax)
        m = rand((-l):l)
        flm[ash_mode_index(0, l, m, lmax)] = 1

        ðflm = ash_eth(flm, 0, lmax)
        ð̄ðflm = ash_ethbar(ðflm, +1, lmax)

        ð̄flm = ash_ethbar(flm, 0, lmax)
        ðð̄flm = ash_eth(ð̄flm, -1, lmax)

        f = ash_evaluate(flm, 0, lmax)
        ð̄ðf = ash_evaluate(ð̄ðflm, 0, lmax)
        ðð̄f = ash_evaluate(ðð̄flm, 0, lmax)

        @test isapprox(ð̄ðf, -l * (l + 1) * f; atol=(lmax + 1)^2 * 100eps())
        @test isapprox(ðð̄f, -l * (l + 1) * f; atol=(lmax + 1)^2 * 100eps())
    end
end

################################################################################

Random.seed!(100)
@testset "Arbitrary modes" begin
    for iter in 1:20
        lmax = rand(0:100)
        nmodes = ash_nmodes(lmax)

        flm = zeros(Complex{Float64}, nmodes)

        s = rand(-4:4)
        abs(s) ≤ lmax || continue
        l = rand(abs(s):lmax)
        m = rand((-l):l)
        flm[ash_mode_index(s, l, m, lmax)] = 1

        f = ash_evaluate(flm, s, lmax)

        f′ = [
            begin
                θ, ϕ = ash_point_coord(ij, lmax)
                # We need the increased precision of `BigFloat` for `l >≈ 50`
                # AbstractSphericalHarmonics.sYlm(Val(s), Val(l), Val(m), θ, ϕ)
                Complex{Float64}(AbstractSphericalHarmonics.sYlm(Val(s), Val(l), Val(m), big(θ), big(ϕ)))
            end for ij in CartesianIndices(f)
        ]

        @test f ≈ f′
    end
end

################################################################################

Random.seed!(100)
@testset "Tensors on the sphere (rank $D)" for D in 0:4
    for iter in 1:20
        lmax = rand(0:100)

        f = rand_tensor(Val(D), Complex{Float64}, lmax)

        t = Tensor{D}(f, lmax)
        st = SpinTensor{D}(t)
        t′ = Tensor{D}(st)
        st′ = SpinTensor{D}(t′)
        t″ = Tensor{D}(st′)

        @test st′ ≈ st
        @test t″ ≈ t′
    end
end

Random.seed!(100)
@testset "Linearity of tensors on the sphere (rank $D)" for D in 0:4
    for iter in 1:20
        lmax = rand(0:100)

        f = rand_tensor(Val(D), Complex{Float64}, lmax)
        g = rand_tensor(Val(D), Complex{Float64}, lmax)
        α = randn(Complex{Float64})

        h = f + α * g

        t = Tensor{D}(f, lmax)
        u = Tensor{D}(g, lmax)
        v = Tensor{D}(h, lmax)

        st = SpinTensor{D}(t)
        su = SpinTensor{D}(u)
        sv = SpinTensor{D}(v)

        @test sv ≈ st + α * su
    end
end

@testset "Simple derivatives of tensors on the sphere" begin
    for lmax in 1:20
        sz = ash_grid_size(lmax)

        # 1

        s = Tensor{0}(fill(Scalar{Complex{Float64}}(1), sz), lmax)
        ds₀ = Tensor{1}(fill(SVector{2,Complex{Float64}}(0, 0), sz), lmax)

        s̃ = SpinTensor(s)
        ds̃ = tensor_gradient(s̃)
        ds = Tensor(ds̃)

        @test isapprox(ds, ds₀; atol=(lmax + 1)^2 * 1000eps())

        # x

        x = Tensor{0}([
            begin
                θ, ϕ = ash_point_coord(ij, lmax)
                Scalar{Complex{Float64}}(sin(θ) * cos(ϕ))
            end for ij in CartesianIndices(sz)
        ], lmax)
        dx₀ = Tensor{1}([
            begin
                θ, ϕ = ash_point_coord(ij, lmax)
                SVector{2,Complex{Float64}}(cos(θ) * cos(ϕ), -sin(ϕ))
            end for ij in CartesianIndices(sz)
        ], lmax)

        x̃ = SpinTensor(x)
        dx̃ = tensor_gradient(x̃)
        dx = Tensor(dx̃)

        @test isapprox(dx, dx₀; atol=(lmax + 1)^2 * 100eps())

        # y

        y = Tensor{0}([
            begin
                θ, ϕ = ash_point_coord(ij, lmax)
                Scalar{Complex{Float64}}(sin(θ) * sin(ϕ))
            end for ij in CartesianIndices(sz)
        ], lmax)
        dy₀ = Tensor{1}([
            begin
                θ, ϕ = ash_point_coord(ij, lmax)
                SVector{2,Complex{Float64}}(cos(θ) * sin(ϕ), cos(ϕ))
            end for ij in CartesianIndices(sz)
        ], lmax)

        ỹ = SpinTensor(y)
        dỹ = tensor_gradient(ỹ)
        dy = Tensor(dỹ)

        @test isapprox(dy, dy₀; atol=(lmax + 1)^2 * 100eps())

        # z

        z = Tensor{0}([
            begin
                θ, ϕ = ash_point_coord(ij, lmax)
                Scalar{Complex{Float64}}(cos(θ))
            end for ij in CartesianIndices(sz)
        ], lmax)
        dz₀ = Tensor{1}([
            begin
                θ, ϕ = ash_point_coord(ij, lmax)
                SVector{2,Complex{Float64}}(-sin(θ), 0)
            end for ij in CartesianIndices(sz)
        ], lmax)

        z̃ = SpinTensor(z)
        dz̃ = tensor_gradient(z̃)
        dz = Tensor(dz̃)

        @test isapprox(dz, dz₀; atol=(lmax + 1)^2 * 100eps())

        if lmax ≥ 2

            # grad x
            #
            # q_θθ = 1
            # q_ϕϕ = sin(ϕ)^2
            #
            # Γ^θ_ϕϕ = - cos(θ) sin(θ)
            # Γ^ϕ_θϕ = cos(θ) / sin(θ)
            # Γ^ϕ_ϕθ = cos(θ) / sin(θ)
            #
            # x = sin(θ) cos(ϕ)
            # y = sin(θ) sin(ϕ)
            # z = cos(θ)
            #
            # ∇x_θ = cos(θ) cos(ϕ)
            # ∇x_ϕ = - sin(θ) sin(ϕ)
            #
            # ∇∇x_θθ = - sin(θ) cos(ϕ)
            # ∇∇x_ϕϕ = - sin(θ)^3 cos(ϕ)

            gradx = Tensor{1}([
                begin
                    θ, ϕ = ash_point_coord(ij, lmax)
                    SVector{2,Complex{Float64}}(cos(θ) * cos(ϕ), -sin(ϕ))
                end for ij in CartesianIndices(sz)
            ], lmax)
            dgradx₀ = Tensor{2}(
                [
                    begin
                        θ, ϕ = ash_point_coord(ij, lmax)
                        SMatrix{2,2,Complex{Float64}}(-sin(θ) * cos(ϕ), 0, 0, -sin(θ) * cos(ϕ))
                    end for ij in CartesianIndices(sz)
                ], lmax
            )

            gradx̃ = SpinTensor(gradx)
            dgradx̃ = tensor_gradient(gradx̃)
            dgradx = Tensor(dgradx̃)

            @test isapprox(dgradx, dgradx₀; atol=(lmax + 1)^2 * 100eps())

            ddgradx₀ = Tensor{3}(
                [
                    begin
                        θ, ϕ = ash_point_coord(ij, lmax)
                        SArray{Tuple{2,2,2},Complex{Float64}}(-cos(θ) * cos(ϕ), 0, 0, -cos(θ) * cos(ϕ), sin(ϕ), 0, 0, sin(ϕ))
                    end for ij in CartesianIndices(sz)
                ],
                lmax,
            )

            ddgradx̃ = tensor_gradient(dgradx̃)
            ddgradx = Tensor(ddgradx̃)

            @test isapprox(ddgradx, ddgradx₀; atol=(lmax + 1)^2 * 1000eps())

            # curl x
            #
            # curl x_θ = sin(ϕ)
            # curl x_ϕ = cos(θ) sin(θ) cos(ϕ)
            #
            # ∇ curl x_θθ = - sin(θ) cos(ϕ)
            # ∇ curl x_ϕϕ = - sin(θ)^3 cos(ϕ)

            curlx = Tensor{1}([
                begin
                    θ, ϕ = ash_point_coord(ij, lmax)
                    SVector{2,Complex{Float64}}(sin(ϕ), cos(θ) * cos(ϕ))
                end for ij in CartesianIndices(sz)
            ], lmax)
            dcurlx₀ = Tensor{2}(
                [
                    begin
                        θ, ϕ = ash_point_coord(ij, lmax)
                        SMatrix{2,2,Complex{Float64}}(0, -sin(θ) * cos(ϕ), sin(θ) * cos(ϕ), 0)
                    end for ij in CartesianIndices(sz)
                ], lmax
            )

            curlx̃ = SpinTensor(curlx)
            dcurlx̃ = tensor_gradient(curlx̃)
            dcurlx = Tensor(dcurlx̃)

            @test isapprox(dcurlx, dcurlx₀; atol=(lmax + 1)^2 * 100eps())
        end
    end
end

Random.seed!(100)
@testset "Derivatives of tensors on the sphere (rank $D)" for D in 0:3
    for iter in 1:20
        lmax = rand(0:100)

        f = rand_tensor(Val(D), Complex{Float64}, lmax)
        g = rand_tensor(Val(D), Complex{Float64}, lmax)
        α = randn(Complex{Float64})

        # Linear combination
        s = f + α * g

        tf = Tensor{D}(f, lmax)
        tg = Tensor{D}(g, lmax)
        ts = Tensor{D}(s, lmax)

        stf = SpinTensor{D}(tf)
        stg = SpinTensor{D}(tg)
        sts = SpinTensor{D}(ts)

        dstf = tensor_gradient(stf)
        dstg = tensor_gradient(stg)
        dsts = tensor_gradient(sts)

        dtf = Tensor{D + 1}(dstf)
        dtg = Tensor{D + 1}(dstg)
        dts = Tensor{D + 1}(dsts)

        @test dts ≈ dtf + α * dtg

        dstf′ = zero(dstf)
        dstg′ = zero(dstg)
        dsts′ = zero(dsts)
        tensor_gradient!(dstf′, stf)
        tensor_gradient!(dstg′, stg)
        tensor_gradient!(dsts′, sts)
        @test dstf′ ≈ dstf
        @test dstg′ ≈ dstg
        @test dsts′ ≈ dsts

        if D == 0
            # Constant function (derivative is zero)
            c = const_tensor(Val(D), Complex{Float64}, lmax)
            # Product of two functions
            p = map(.*, f, g)

            tc = Tensor{D}(c, lmax)
            tp = Tensor{D}(p, lmax)

            stc = SpinTensor{D}(tc)
            stp = SpinTensor{D}(tp)

            dstc = tensor_gradient(stc)
            dstp = tensor_gradient(stp)

            dtc = Tensor{D + 1}(dstc)
            dtp = Tensor{D + 1}(dstp)

            dtc′ = Tensor{D + 1}(zero(dtc.values), lmax)
            @test isapprox(dtc, dtc′; atol=sqrt(eps()))
            dtp′ = Tensor{D + 1}(
                map((dx, x) -> dx * x[], dtf.values, tg.values) + map((x, dx) -> x[] * dx, tf.values, dtg.values), lmax
            )
            # TODO: This needs smooth input, not noise
            # if !isapprox(dtp, dtp′; atol=1 / lmax^2)
            #     @show D iter lmax
            #     @show maximum(map(x -> maximum(abs.(x)), dtp.values - dtp′.values))
            # end
            @test_skip isapprox(dtp, dtp′; atol=1 / lmax^2)
        end
    end
end

Random.seed!(100)
@testset "Filtering tensors on the sphere D=$D" for D in 0:4
    for lmax in 0:20
        sz = ash_grid_size(lmax)
        t = Tensor{D}(rand_tensor(Val(D), Complex{Float64}, lmax), lmax)
        st = SpinTensor(t)
        st′ = filter_modes(st)
        # This is only a weak test
        lfilter = lmax * 2 ÷ 3
        for (ij, cs) in zip(CartesianIndices(st′.coeffs), st′.coeffs)
            s = count(Tuple(ij) .== 1) - count(Tuple(ij) .== 2)
            for l in abs(s):lmax, m in (-l):l
                if l ≤ lfilter
                    @test cs[ash_mode_index(s, l, m, lmax)] ≠ 0
                else
                    @test abs(cs[ash_mode_index(s, l, m, lmax)]) ≤ 1000eps()
                end
            end
        end
        # Ideas:
        # - test effect of filtering on particular modes
        # - test that filtering and derivatives commute
        # - test that filtering is linear
    end
end

Random.seed!(100)
@testset "Calculate Ricci scalar" begin
    for iter in 1:20
        lmax = rand(2:100)
        sz = ash_grid_size(lmax)

        # q_ab = (1, 0, 0, 1)
        q = Tensor{2}(fill(SMatrix{2,2,Complex{Float64}}(1, 0, 0, 1), sz), lmax)
        q̃ = SpinTensor(q)
        # qu^ab
        qu = q
        # q_ab,c
        dq̃ = tensor_gradient(q̃)
        dq = Tensor(dq̃)
        # Γ^a_bc = q^ad (q_dc,b + q_bd,c - q_bc,d) / 2
        Γ = Tensor{3}(
            [
                SArray{Tuple{2,2,2}}(
                    sum(qu[a, d] * (dq[d, c, b] + dq[b, d, c] - dq[b, c, d]) / 2 for d in 1:2) for a in 1:2, b in 1:2, c in 1:2
                ) for (qu, dq) in zip(qu.values, dq.values)
            ],
            lmax,
        )
        # dΓ^a_bc,d
        Γ̃ = SpinTensor(Γ)
        dΓ̃ = tensor_gradient(Γ̃)
        dΓ = Tensor(dΓ̃)
        # R_ab = dΓ^c_ab,c - dΓ^c_ac,b + Γ^c_ab Γ^d_cd - Γ^c_ad Γ^d_bc
        R = Tensor{2}(
            [
                SMatrix{2,2}(
                    sum(dΓ[c, a, b, c] - dΓ[c, a, c, b] for c in 1:2) +
                    sum(Γ[c, a, b] * Γ[d, c, d] - Γ[c, a, d] * Γ[d, b, c] for c in 1:2, d in 1:2) for a in 1:2, b in 1:2
                ) for (Γ, dΓ) in zip(Γ.values, dΓ.values)
            ],
            lmax,
        )
        Rsc = Tensor{0}([Scalar(sum(qu[a, b] * R[a, b] for a in 1:2, b in 1:2)) for (qu, R) in zip(qu.values, R.values)], lmax)

        @test isapprox(map(x -> x[], Rsc.values), zeros(sz); atol=(lmax + 1)^4 * 10eps())

        # println("q:")
        # display(map(x -> chop.(x), q.values))
        # println()
        # println("Γ:")
        # display(map(x -> chop.(x), Γ.values))
        # println()
        # println("R:")
        # display(map(x -> chop.(x), R.values))
        # println()
        # println("Rsc:")
        # display(map(x -> chop.(x)[], Rsc.values))
        # println()
    end
end


text:
# Here's the problem:
#
# You have linearly spaced buckets for allocation.
# You try to allocate a bucket of size `N`, only to find that the freelist for size `N`
# is empty.
#
# How do you efficiently find the first non-empty bucket of size greater than `N`?
#
# One option is to linearly scan the freelists ... but that's linear.
#
# The approach here is use a `FindNextTree`.
#
# The idea is that we use a bitmask UInt64 to encode groups of 64 buckets.
# A 1 in the mask represents that a bucket is non-empty.
# A quick lookup in the mask can be done using an `and` and a call to `trailing_zeros`.
#
# Then, we encode a next level in the tree, also using UInt64 masks which will then
# cover 4096 buckets - using a 1 to indicate if any entry in the sub-tree is non empty.
# We can continue building trees as long as we want.
struct Mask
    val::UInt64
end

Base.:|(a::Mask, b::UInt64) = Mask(a.val | b)
Base.:&(a::Mask, b::UInt64) = Mask(a.val & b)

"""
    firstentry(m::Mask, [i::Integer])

Return the index of the first non-zero entry in `m`.
If `i` is provided, return the index of the first non-zero entry in `m` greater than `i`.

Behavior is undefined if `i > 64` - so don't so that.
"""
firstentry(m::Mask, i::Integer = 0) = trailing_zeros(m.val >> i) + i + 1

"""
    isvalidentry(x::Integer) -> Bool

Return `true` if the value returned by `firstentry` is a valid entry.
"""
isvalidentry(x::Integer) = x <= 64

"""
    hasentry(m::Mask, [i::Integer]) -> Bool

Return `true` is `m` has at least one non-zero entry.
If `i` is provided, return `true` is `m` has at least one non-zero entry greater than `i`.
"""
hasentry(m::Mask, i::Integer = 0) = (m.val >> i) >= one(UInt64)
hasentryat(m::Mask, i::Integer) = isodd(m.val >> (i-1))

"""
    numentries(m::Mask) -> Integer

Return the number of set bits in `m`.
"""
numentries(m::Mask) = count_ones(m.val)

cdiv(x::T, y::T) where {T} = one(T) + div(x - one(T), y)

struct FindNextTree
    # runs[1] is the highest-level
    # higher indices descend through the tree.
    # runs[N] indicate the leaves.
    runs::Vector{Vector{Mask}}
    length::Int
end

function FindNextTree(len::Integer)
    original_length = len
    runs = Vector{Vector{Mask}}()

    # Keep adding entries until we've covered the whole tree.
    while true
        thislen = cdiv(len, 64)
        pushfirst!(runs, [Mask(0) for _ in 1:thislen])
        thislen == 1 && break
        len = thislen
    end
    return FindNextTree(runs, original_length)
end

# Faster mod for powers of 64
@inline mod64(i::T) where {T} = i & T(63)
@inline div64(i::T) where {T} = i >> 6

@inline function divrem64(x::T) where {T}
    mod = mod64(x)
    div = div64(x)
    return iszero(mod) ? (div, T(64)) : (div+1, mod)
end

Base.@propagate_inbounds function isset(M::FindNextTree, index::Integer)
    @boundscheck begin
        index > M.length && throw(BoundsError(M, index))
    end

    run = @inbounds M.runs[length(M.runs)]
    newindex, modindex = divrem64(index)
    return hasentryat(@inbounds(run[newindex]), modindex)
end

# Setting and clearing entries.
Base.@propagate_inbounds function setentry!(M::FindNextTree, index::Integer)
    @boundscheck begin
        index > M.length && throw(BoundsError(M, index))
    end

    level = length(M.runs)
    while true
        run = @inbounds M.runs[level]

        newindex, modindex = divrem64(index)
        @inbounds update = !hasentry(run[newindex])
        @inbounds run[newindex] |= (one(UInt64) << (modindex-1))

        # Break out if either we don't need an update, or this is the top level.
        (update && level > 1) || break
        index = newindex
        level -= 1
    end
    return nothing
end

Base.@propagate_inbounds function clearentry!(M::FindNextTree, index::Integer)
    @boundscheck begin
        index > M.length && throw(BoundsError(M, index))
    end

    level = length(M.runs)
    while true
        run = M.runs[level]

        newindex, modindex = divrem64(index)

        @inbounds run[newindex] &= ~(one(UInt64) << (modindex-1))
        @inbounds update = !hasentry(run[newindex])
        (update && level > 1) || break
        index = newindex
        level -= 1
    end
    return nothing
end

Base.@propagate_inbounds function Base.findnext(M::FindNextTree, index::Integer)
    index > 64 * length(M.runs[end]) && return nothing
    level = length(M.runs)

    # Go up the tree until we find a non-empty entry.
    @inbounds while true
        run = M.runs[level]
        newindex, modindex = divrem64(index)
        mask = run[newindex]

        # Find the first entry for this index.
        # If it is valid, then we're good to go.
        entry = firstentry(mask, modindex)
        if isvalidentry(entry)
            # Go into this level.
            index = ((newindex - 1) << 6) + entry
            level += 1
            break
        end

        index = newindex
        # If we're at the top level and haven't found anything, we won't find anything.
        level == 1 && return nothing
        level -= 1
    end

    # Now, we have to go back down the tree to get the final index
    @inbounds while level <= length(M.runs)
        run = M.runs[level]
        index = ((index - 1) << 6) + firstentry(run[index])
        level += 1
    end
    return index
end


text:
acknowledge help from Tommy Hofmann

****************************
linearcode.jl
****************************
include("linearcode.jl");

# To create a linear code, need to make generator matrix. We start by creating a field for this.
F, _ = FiniteField(2, 1, "α");
# one may either create a matrix space and then embedd the matrix into it as in
M = MatrixSpace(F, 4, 7);
G = M([1 0 0 0 0 1 1;
       0 1 0 0 1 0 1;
       0 0 1 0 1 1 0;
       0 0 0 1 1 1 1]);
# or use either of the two shorthand notations
G2 = matrix(F, [1 0 0 0 0 1 1;
       0 1 0 0 1 0 1;
       0 0 1 0 1 1 0;
       0 0 0 1 1 1 1]);
G3 = matrix(F, 4, 7, [1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1]);
G == G2
G == G3

# This is the generator matrix for the [7, 4, 3] binary Hamming code in standard form
C = LinearCode(G)
# We can get the basic information
length(C)
dimension(C)
cardinality(C)
rate(C)
# since we gave it a full rank matrix, the rank should equal the dimension
rank(G) == dimension(C)
# no minimum distance calculations were performed so the object so far does not know
# this value and the current state of the variable is set to 'missing'
# functions which rely on this value will report an error in this case
minimumdistance(C)
relativedistance(C)
genus(C)
# since we know the minimum distance of this code is 3, we can go ahead and set that
# no checking is done to verify the user's input
# since setting the minimum distance is a desirable thing to be able to do, all code
# objects are mutable (and hence live on the heap)
setminimumdistance!(C, 3)
minimumdistance(C)
relativedistance(C)
genus(C)
isMDS(C)
numbercorrectableerrors(C)
# can check the generator and parity check matrices and check their defining properties
GC = generatormatrix(C)
G == GC
H = paritycheckmatrix(C)
iszero(G * H')
iszero(H * G')
# if G is full rank, the matrix is stored directly as the generator matrix of C as written
# and H is computed as the right kernel of G
# the standard form of G and H may be retrieved by passing true into the optional parameter
Gstand = generatormatrix(C, true)
Hstand = paritycheckmatrix(C, true)
# G in the example was already in standard form, so we can check the computations
Gstand == GC
Hstand == H
# finally, we may compute the dual
# in order to save computational time, this command merely flipps the roles of G and H
# in the LinearCode constructor, and no new computation is done
dual(C)

# a vector v is in the code C if it has zero syndrome
r = generatormatrix(C)[1, :]
iszero(syndrome(r, C))
# a code C1 is a subset of C2 if every row of the generator matrix of C1 is in C2
C ⊆ C
C ⊆ dual(C)
issubcode(C, dual(C))
# two codes C1 and C2 are equivalent if C1 ⊆ C2 and C2 ⊆ C1
isequivalent(C, dual(C))
# a code is self dual if it is equivalent to its dual
isselfdual(C)
# code is self orthogonal if it is a subcode of its dual
isselforthogonal(C)

# finally, one may encode a length k vector by C using
encode(C.G[:, 1], C)
# here the vector to be encoded is an element of MatrixSpace(Nemo.GF(2), 4, 1)
# allowing for proper matrix multiplication between it and the generator matrix
# one may also pass in a Vector{Int64} of the proper length and the appropriate
# matrix space will automatically be created. this may be slow if this is done
# repeatedly as this object will be built from scratch each time. the same is
# true for syndrome(). either a (1, C.k) or a (C.k, 1) dimensional vector may
# be used as the code will automatically apply the transpose if necessary
# it will not however accept objects of type adjoint such as v'
generatormatrix(C)[:, 1]
v = [1, 0, 0, 0];
encode(v, C)
v2 = [1; 0; 0; 0];
encode(v2, C)
generatormatrix(C)[1, :]
v = [1, 0, 0, 0, 0, 1, 1];
syndrome(v, C)
v = [1; 0; 0; 0; 0; 1; 1];
syndrome(v, C)
# might want to remove these transpose examples and remove from code since these have same size

# need examples, not full rank
# passing in H
# passing in G and H directly

# we can preform several standard methods of building new codes from old ones
# direct sum
D = C ⊕ C;
generatormatrix(D)
D = directsum(C, C);
generatormatrix(D)
# direct product
D = C ⊗ C;
generatormatrix(D)
D = directproduct(C, C);
generatormatrix(D)
D = tensorproduct(C, C);
generatormatrix(D)
# the most common form of extending a code is to add parity check column to right of generator matrix
D = extend(C);
length(D)
generatormatrix(D)
F3, _ = FiniteField(3, 1, "α");
# M2 = MatrixSpace(Nemo.GF(3), 2, 4); #- ex page 15
G2 = matrix(F3, [1 0 1 1; 0 1 1 -1]);
C2 = LinearCode(G2);
generatormatrix(extend(C2))
paritycheckmatrix(extend(C2))
originalgeneratormatrix(C2) == G2
# in cases such as this it is often useful to keep track of the original generator matrix which
# was used to generate the code. can also do this for the parity check matrix. in cases where this
# does not make sense, the normal generator or parity check matrix is returned
originalgeneratormatrix(D)
originalgeneratormatrix(D) == generatormatrix(D)[:, 1:end-1]
paritycheckmatrix(D)
originalparitycheckmatrix(D)
# puncturing deletes columns from a generator matrix
C2 = puncture(D, [length(D)]);
isequivalent(C, C2)
generatormatrix(C) == generatormatrix(C2)
originalgeneratormatrix(C2) == generatormatrix(C)
# to expurgate a code is to delete rows from generator matrix and then remove any potentially resulting zero columns
D = expurgate(C, [dimension(C)]);
length(D)
dimension(D)
generatormatrix(D)
originalgeneratormatrix(D) == generatormatrix(C)
# to augment is to add rows to the generator matrix
# in this example we are going to first add the all 1's vector to the code, which is common
# we note that the all 1's vector is actually already in the code so the code
# shouldn't change
# M = MatrixSpace(Nemo.GF(2), 1, 7);
v = matrix(F, [1 for i in 1:7]')
v ∈ C
D = augment(C, v);
generatormatrix(D) == generatormatrix(C)
v = matrix(F, [1 0 1 0 1 1 1])
v ∈ C
D = augment(C, v);
generatormatrix(D)
originalgeneratormatrix(D) == generatormatrix(C)
# to shorten is to expurgate then puncture
D = shorten(C, [length(C)]);
generatormatrix(D)
originalgeneratormatrix(D) == generatormatrix(C)
# M2 = MatrixSpace(Nemo.GF(2), 3, 6);
G2 = matrix(F, [1 0 0 1 1 1; 0 1 0 1 1 1; 0 0 1 1 1 1]);
C2 = LinearCode(G2);
generatormatrix(shorten(C2, [5, 6]))
# the most common form of lengthening a code it to augment the all 1's vector and then extend
D = lengthen(C);
generatormatrix(D)
originalgeneratormatrix(D) == generatormatrix(C)
# the (u | u + v)- or Plotkin construction
# M2 = MatrixSpace(Nemo.GF(2), 3, 4);
# M3 = MatrixSpace(Nemo.GF(2), 1, 4);
G2 = matrix(F, [1 0 1 0; 0 1 0 1; 0 0 1 1]);
G3 = matrix(F, [1 1 1 1]);
C2 = LinearCode(G2);
C3 = LinearCode(G3);
generatormatrix(uuplusv(C2, C3)) #- ex p. 19

SingletonBound

****************************
cyclotomic.jl
****************************
include("cyclotomic.jl")

# ord_n(q)
q = 2;
n = 15;
b = 3;
δ = 4;
ord(n, q)
# at the moment this comes back sorted so one can easily record a coset representative
coset = cyclotomiccoset(3, q, n, false)
cyclotomiccoset(3, q, n, true)
allcyclotomiccosets(q, n, true)
dualqcosets(q, n, [coset])
complementqcosets(q, n, [coset])
# pairings
qcosetpairings(q, n)
# find all in a range
qcosettable(10, 13, q)

****************************
cycliccode.jl
****************************
include("cycliccode.jl")
q = 2;
n = 15;
b = 3;
δ = 4;
cosets = definingset([i for i = b:(b + δ - 2)], q, n, false)
C = CyclicCode(q, n, cosets)
basefield(C)
splittingfield(C)
polynomialring(C)
primitiveroot(C)
offset(C)
designdistance(C)
qcosets(C)
qcosetsreps(C)
definingset(C)
generatorpolynomial(C)
paritycheckpolynomial(C)
idempotent(C)
# CyclicCode is a subtype of LinearCode so the functions of the previous section also apply here
originalgeneratormatrix(C)
G = matrix(F, [1 0 0 0 0 1 1;
       0 1 0 0 1 0 1;
       0 0 1 0 1 1 0;
       0 0 0 1 1 1 1]);
D = LinearCode(G);
CD = directsum(C, D);
generatormatrix(CD)
# the generator matrix of a cyclic code is kept in cyclic form but the standard form is of course always available
generatormatrix(C)
generatormatrix(C, true)
# note how the previous code auto detected that it was a BCH code and returned an object of type BCHCode
complement(C)
# the above is equivalent to
B = BCHCode(q, n, δ, b)
C == B
# the optional offset parameter is set to 0 by default
BCHCode(q, n, δ)
# for cyclic codes, we can
C ∩ B == C
B = BCHCode(q, n, δ - 1, b + 4)
D = C ∩ B # check later that this is == repetition code
C + B

# Reed Solomon codes are BCH codes with n = q^m - 1
# this forces the cyclotomic cosets to each have size one
q = 16;
n = 15;
b = 3;
δ = 4;
allcyclotomiccosets(q, n, true)
# as with BCH codes, it will auto detect whether or not the code is Reed Solomon and call the appropriate constructor
cosets = definingset([i for i = b:(b + δ - 2)], q, n, false);
CyclicCode(q, n, cosets)
BCHCode(q, n, δ, b)

############## this doesn't work
ReedSolomonCode(q, δ, b)

come up with better examples from the book but basis functionality is working so far and interaction with LinearCode is good
fix C ∩ complement(C) example where defining set is an extreme (Vector{Any} so probably empty)
2nd form of CyclicCode constructor

****************************
ReedMuller.jl
****************************
include("ReedMuller.jl")

# so far only binary implemented
# fields are created via F, _ = FiniteField(2, 1, "α") by default
# codes are created using the recursive form of the generator matrix
q = 2;
ReedMullergeneratormatrix(q, 0, 3)
ReedMullergeneratormatrix(q, 2, 2)
ReedMullerCode(q, 1, 2)
C = ReedMullerCode(q, 1, 3)
ReedMullerCode(q, 2, 3)
isselfdual(C)

max_stars_repo_name: esabo/CodingTheory

max_issues_repo_name: esabo/CodingTheory

max_forks_repo_name: esabo/CodingTheory


text:
# if it's one symbol, set the guide and return the column
function handle_dfs(df::AbstractDataFrame, d::KW, letter, sym::Symbol)
    get!(d, Symbol(letter * "guide"), string(sym))
    collect(df[sym])
end

# if it's an array of symbols, set the labels and return a Vector{Any} of columns
function handle_dfs(df::AbstractDataFrame, d::KW, letter, syms::AbstractArray{Symbol})
    get!(d, :label, reshape(syms, 1, length(syms)))
    Any[collect(df[s]) for s in syms]
end

# for anything else, no-op
function handle_dfs(df::AbstractDataFrame, d::KW, letter, anything)
    anything
end

# handle grouping by DataFrame column
function Plots.extractGroupArgs(group::Symbol, df::AbstractDataFrame, args...)
    Plots.extractGroupArgs(collect(df[group]))
end

# if a DataFrame is the first arg, lets swap symbols out for columns
@recipe function f(df::AbstractDataFrame, args...)
    # if any of these attributes are symbols, swap out for the df column
    for k in (:fillrange, :line_z, :marker_z, :markersize, :ribbon, :weights, :xerror, :yerror)
        if haskey(d, k) && isa(d[k], Symbol)
            d[k] = collect(df[d[k]])
        end
    end

    # return a list of new arguments
    tuple(Any[handle_dfs(df, d, (i==1 ? "x" : i==2 ? "y" : "z"), arg) for (i,arg) in enumerate(args)]...)
end


text:
@testset "AlexNet" begin
  model = alexnet()
  @test size(model(rand(Float32, 256, 256, 3, 50))) == (1000, 50)
end


text:
using Terming
using Crayons

function main()
    Terming.cshow(false)
    Terming.displaysize(20, 75); Terming.clear()

    julia_str = [
    "     _       _ _       ",
    "    | |_   _| (_) __ _ ",
    " _  | | | | | | |/ _` |",
    "| |_| | |_| | | | (_| |",
    " \\___/ \\__,_|_|_|\\__,_|"
    ]

    Terming.buffered() do buffer
        Terming.cmove(2, 2)
        for str in julia_str
            Terming.csave(buffer)
            Terming.print(buffer, Crayon(foreground=(175, 122, 197), bold=false), str)
            Terming.crestore(buffer)
            Terming.cmove_down(buffer)
        end
        Terming.print(buffer, "-"^73)

        Terming.cmove_line_down(buffer, 2); Terming.cmove_col(buffer, 2)
        Terming.print(buffer, Crayon(foreground=(93, 173, 226), bold=true), "[noun]")

        Terming.cmove_line_down(buffer, 1); Terming.cmove_col(buffer, 6)
        def_str = "The Julia Language: A fresh approach to technical computing."
        Terming.print(buffer, Crayon(foreground=(88, 214, 141), bold=false), def_str)

        Terming.cmove_line_down(buffer, 2); Terming.cmove_col(buffer, 2)
        Terming.print(buffer, Crayon(foreground=(93, 173, 226), bold=true), "[verb]")

        Terming.cmove_line_down(buffer, 1); Terming.cmove_col(buffer, 6)
        def_str = "To write code in Julia."
        Terming.print(buffer, Crayon(foreground=(88, 214, 141), bold=false), def_str)

        Terming.cmove_line_last(buffer)
    end

    readline()

    Terming.cshow(true)

    return
end

main()


text:
# !! auto-generated code, do not edit !!
module nav_msgs
using RobotOSData.Messages
using RobotOSData.StdMsgs
using RobotOSData.CommonMsgs
struct MapMetaData <: Readable
    map_load_time::ROSTime
    resolution::Float32
    width::UInt32
    height::UInt32
    origin::geometry_msgs.Pose
end
struct GridCells <: Readable
    header::Header
    cell_width::Float32
    cell_height::Float32
    cells::Vector{geometry_msgs.Point}
end
struct Odometry <: Readable
    header::Header
    child_frame_id::String
    pose::geometry_msgs.PoseWithCovariance
    twist::geometry_msgs.TwistWithCovariance
end
struct OccupancyGrid <: Readable
    header::Header
    info::MapMetaData
    data::Vector{Int8}
end
struct Path <: Readable
    header::Header
    poses::Vector{geometry_msgs.PoseStamped}
end
end


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

using Roots
import DataFrames

# luminance (cd/m2)	Multiple	Value	Item
# 10−6	µcd/m2	1 µcd/m2	Absolute threshold of vision[1]
# 10−5
# 10−4		400 µcd/m2	Darkest sky[2]
# 10−3	mcd/m2	1 mcd/m2	Night sky[3]
# 1.4   mcd/m2	Typical photographic scene lit by full moon[4]
# 5     mcd/m2	Approximate scotopic/mesopic threshold[5]
# 10−2		40 mcd/m2	Phosphorescent markings on a watch dial after 1 h in the dark[6][7]
# 10−1
# 100	cd/m2	2 cd/m2	Floodlit buildings, monuments, and fountains[8]
# 5     cd/m2	Approximate mesopic/photopic threshold[5]
# 101		25 cd/m2	Typical photographic scene at sunrise or sunset[4]
# 30    cd/m2	Green electroluminescent source[2]
# 102		250 cd/m2	Peak luminance of a typical LCD monitor[10][11]
# 700   cd/m2	Typical photographic scene on overcast day[4][8][11]
# 103	kcd/m2	2 kcd/m2	Average cloudy sky[2]
# 5     kcd/m2	Typical photographic scene in full sunlight[4][8]

"""This function returns the radius of the longest basis vector of the lattice cluster. Most lattices defined in this project have symmetric basis vectors so the radii of all basis vectors will be identical. This function is used in determing which clusters "fit" within the eye pupil."""
function latticediameter(basismatrix::SMatrix)
    maximum(norm.([basismatrix[:,i] for i in 1:size(basismatrix)[2]]))
end
export latticediameter

latticediameter(a::Repeat.AbstractLatticeCluster) =   latticediameter(Repeat.basismatrix(Repeat.clusterbasis(a)))
latticediameter(a::Repeat.AbstractBasis) =  latticediameter(Repeat.basismatrix(a))

const hex3latticeclusterbasis = [2 // 1 -1 // 1;-1 // 1 2 // 1]
export hex3latticeclusterbasis

"""Only will work properly if lattice basis matrix contains only integer or rational terms. Returns the integer lattice coords of point in the given basis if the point is in the span of latticebasis. Otherwise returns nothing"""
function latticepoint(latticebasis::AbstractMatrix, origin, point)
    Ainv = inv(Rational.(latticebasis))
    b = [(point .- origin)...]
    x = Ainv * b
    if reduce(&, (1, 1) .== denominator.(x))
        return Integer.(x)
    else
        return nothing
    end
end
export latticepoint

colorbasis(::Repeat.HexBasis1) = SMatrix{2,2}(2, -1, 1, 1)
colorbasis(::Repeat.HexBasis3) = SMatrix{2,2}(2, -1, 1, 1)
colororigins(::Repeat.HexBasis1) = ((0, 0), (-1, 0), (-1, 1))
colororigins(::Repeat.HexBasis3) = ((0, 0), (0, -1), (1, -1))

"""To reduce color channel cross talk it may be useful to arrange the color lenslets in a hexagonal lattice so that each color is surrounded by lenslets of a different color. With appropriate color filtering crosstalk between immediately neighboring lenslets can be reduced. This function computes the color which should be assigned to any lattice point to ensure this properly holds."""
function pointcolor(point, cluster::Repeat.AbstractLatticeCluster)
    latticematrix = colorbasis(Repeat.elementbasis(cluster))
    origins = colororigins(Repeat.elementbasis(cluster))
    colors = zip(origins, ("red", "green", "blue"))
    for (origin, color) in colors
        if nothing !== latticepoint(latticematrix, origin, point)
            return color
        end
    end
end
export pointcolor

defaultclusterproperties() = (mtf = .2, minfnumber = 2.0, cyclesperdegree = 11, λ = 530nm, pixelpitch = .9μm)
export defaultclusterproperties

"""maximum allowable display for visibility reasons"""
maxlensletdisplaysize() = (250μm, 250μm)

const ρ_quartervalue = 2.21509 # value of ρ at which the airy disk function has magnitude .25
export ρ_quartervalue
const ρ_zerovalue = 3.832 # value of ρ at which the airy disk function has magnitude 0

"""given pixelpitch and angular subtense (in degrees) of pixel returns focal length"""
focallength(pixelpitch::Unitful.Length,θ) = uconvert(mm, .5 * pixelpitch / tand(θ / 2))
export focallength

"""returns the diffraction limit frequency in cycles/deg

focal length = 𝒇𝒍
diffraction cutoff frequency,fc, in cycles/mm = 1/λF# = diameter/λ*𝒇𝒍
cutoff wavelength, Wc, = 1/cutoff frequency = λ*𝒇𝒍/diameter

angular wavelength, θc, radians/cycle, corresponding to cutoff wavelength:

θ ≈ tanθ for small θ
θc corresponding to Wc:  θc ≈ tanθ = Wc/𝒇𝒍
Wc = θc*𝒇𝒍

from equation for Wc:

λ*𝒇𝒍/diameter = Wc = θc*𝒇𝒍
θc = λ/diameter
cycles/rad = 1/θc = diameter/λ

"""
cyclesperdegree(diameter::Unitful.Length,λ::Unitful.Length) = uconvert(Unitful.NoUnits, diameter / (rad2deg(1) * λ))
export cyclesperdegree

"""diffraction limited response for a circular aperture, normalized by maximum cutoff frequency"""
function mtfcircular(freq, freqcutoff)
    s =  freq / freqcutoff
    return 2 / π * (acos(s) - s * ((1 - s^2)^.5))
end
export mtfcircular

"""returns the diffraction limit frequency in cycles/degree. At this frequeny the response of the system is zero"""
diffractionlimit(λ::Unitful.Length,diameter::Unitful.Length) = uconvert(Unitful.NoUnits, diameter / λ) / rad2deg(1)
export diffractionlimit

"""computes the diameter of a diffraction limited lens that has response mtf at frequency cyclesperdeg"""
function diameter_for_cycles_deg(mtf, cyclesperdeg, λ::Unitful.Length)
    cyclesperrad = rad2deg(1) * cyclesperdeg
    f(x) = mtfcircular(x, 1.0) - mtf
    normalizedfrequency = find_zero(f, (0.0, 1.0))
    fc = cyclesperrad / normalizedfrequency
    return uconvert(mm, λ * fc)
end
export diameter_for_cycles_deg

airydisk(ρ) = (2 * SpecialFunctions.besselj1(ρ) / ρ)^2

"""The f# required for the first zero of the airy diffraction disk to be at the next sample point"""
diffractionfnumber(λ::Unitful.Length,pixelpitch::Unitful.Length,indexofrefraction) = uconvert(Unitful.NoUnits, pixelpitch / (2.44λ / indexofrefraction))
export diffractionfnumber

"""returns ρ, the normalized distance, at which the airy disk will have the value airyvalue"""
function ρatairyvalue(airyvalue)
    ρ = 1e-5 # start at small angle
    while airydisk(ρ) > airyvalue
        ρ += 1e-5 # numerical precision not an issue here
    end
    return ρ
end
export ρatairyvalue

""" Spacing between lenslets which guarantess that for any pixel visible in all lenslets every point on the eyebox plane is covered. This is closest packing of circles."""
closestpackingdistance(pupildiameter::Unitful.Length) = pupildiameter * cosd(30)
export closestpackingdistance

"""Tries clusters of various sizes to choose the largest one which fits within the eye pupil. Larger clusters allow for greater reduction of the fov each lenslet must cover so it returns the largest feasible cluster"""
function choosecluster(pupildiameter::Unitful.Length, lensletdiameter::Unitful.Length)
    clusters = (hex3RGB(), hex4RGB(), hex7RGB(), hex9RGB(), hex12RGB(), hex19RGB())
    # cdist = closestpackingdistance(pupildiameter)
    cdist = pupildiameter
    maxcluster = clusters[1]
    ratio = 0.0

    for cluster in clusters
        temp = cdist / (lensletdiameter * latticediameter(cluster))
            if temp >= 1.0
            ratio = temp
            maxcluster = cluster
        else
            break
        end
    end

    @assert ratio ≥ 1.0 "ratio $ratio cdist $cdist lensletdiameter $lensletdiameter latticediameter $(latticediameter(maxcluster)) scaled=$(lensletdiameter * latticediameter(maxcluster))"

    return (cluster = maxcluster, lensletdiameter = lensletdiameter * ratio, diameteroflattice = latticediameter(maxcluster) / ratio, packingdistance = cdist * ustrip(mm, lensletdiameter))
end
export choosecluster

"""Computes the largest feasible cluster size meeting constraints"""
function choosecluster(pupildiameter::Unitful.Length, λ::Unitful.Length, mtf, cyclesperdeg::T) where {T <: Real}
    diam = diameter_for_cycles_deg(mtf, cyclesperdeg, λ)
    return choosecluster(pupildiameter, diam) # use minimum diameter for now.
end

"""Computes display size assuming lenslet normal to eyebox plane passes through lenslet. This is an approximation but for the narrow fov we are considering it is accurate enough to estimate pixel redundandcy, etc."""
function sizeoflensletdisplay(eyerelief::T, eyebox::AbstractVector{T}, ppd, pixelpitch::S) where {T <: Unitful.Length,S <: Unitful.Length}
    θ = eyeboxangles(eyebox, eyerelief)
    displaysize = @.  θ * ppd * pixelpitch
    return uconvert.(mm, displaysize)
end
export sizeoflensletdisplay

"""This computes the approximate size of the entire display, not the individual lenslet displays."""
sizeofdisplay(fov,eyerelief) = @. 2 * tand(fov / 2) * eyerelief
export sizeofdisplay

"""Computes the number of lenslets required to create a specified field of view at the given eyerelief"""
function numberoflenslets(fov, eyerelief::Unitful.Length, lensletdiameter::Unitful.Length)
    lensletarea = π * (lensletdiameter / 2)^2
    dispsize = sizeofdisplay(fov, eyerelief)
     return dispsize[1] * dispsize[2] / lensletarea
end
export numberoflenslets

"""given the angles each lenslet has to cover compute the corresponding display size"""
sizeoflensletdisplay(angles,ppd,pixelpitch::Unitful.Length) = @. angles * ppd * pixelpitch

"""angular size of the eyebox when viewed from distance eyerelief"""
eyeboxangles(eyebox,eyerelief) = @. atand(uconvert(Unitful.NoUnits, eyebox / eyerelief))
export eyeboxangles

"""computes how the fov can be subdivided among lenslets based on cluster size. Assumes the horizontal size of the eyebox is larger than the vertical so the larger number of subdivisions will always be the first number in the returns Tuple."""
function anglesubdivisions(pupildiameter::Unitful.Length, λ::Unitful.Length, mtf, cyclesperdegree;RGB=true)
    cluster, _ = choosecluster(pupildiameter, λ, mtf, cyclesperdegree)
    numelements = Repeat.clustersize(cluster)
    if numelements == 19
        return RGB ? (3, 2) : (5, 3)
    elseif numelements == 12
        return RGB ? (2, 2) : (4, 3)
   elseif numelements == 9
        return RGB ? (3, 1) : (3, 3)
    elseif numelements == 7
        return RGB ? (3, 1) : (3, 2)
    elseif numelements == 4 || numelements == 3
        return RGB ? (1, 1) : (3, 1)
    else throw(ErrorException("this should never happen"))
    end
end
export anglesubdivisions

"""Computes the approximate fov required of each lenslet for the given constraints. This is strictly correct only for a lenslet centered in front of the eyebox, but the approximation is good enough for high level analysis"""
function lensletangles(eyerelief::Unitful.Length, eyebox::NTuple{2,Unitful.Length}, pupildiameter::Unitful.Length, ppd; clusterproperties=defaultclusterproperties(), RGB=true)
    cyclesperdegree = ppd / 2.0
    return eyeboxangles(eyebox, eyerelief) ./ anglesubdivisions(pupildiameter, clusterproperties.λ, clusterproperties.mtf, clusterproperties.cyclesperdegree, RGB=RGB)
end
export lensletangles

testangles() = lensletangles(18mm, (10mm, 6mm), 4mm, 45)
export testangles

lensletresolution(angles,ppd) = angles .* ppd

"""Multilens displays tradeoff pixel redundancy for a reduction in total track of the display, by using many short focal length lenses to cover the eyebox. This function computes the ratio of pixels in the multilens display vs. a conventional display of the same nominal resolution"""
function pixelredundancy(fov, eyerelief::Unitful.Length, eyebox::NTuple{2,Unitful.Length}, pupildiameter::Unitful.Length, ppd; RGB=true)
    lensprops = defaultclusterproperties()
    clusterdata = choosecluster(pupildiameter, lensprops.λ, lensprops.mtf, lensprops.cyclesperdegree)
    nominalresolution = fov .* ppd
    angles = lensletangles(eyerelief, eyebox, pupildiameter, ppd, RGB=RGB)
    pixelsperlenslet = angles .* ppd
    numlenses = numberoflenslets(fov, eyerelief, clusterdata.lensletdiameter)
    return (numlenses * pixelsperlenslet[1] * pixelsperlenslet[2]) / ( nominalresolution[1] * nominalresolution[2])
end
export pixelredundancy

testpixelredundancy() = pixelredundancy((55, 35), 18mm, (10mm, 6mm), 4mm, 45, RGB=false)
export testpixelredundancy

label(color) = color ? "RGB" : "Monochrome"

"""generates a contour plot showing pixel redundancy as a function of ppd and pupil diameter"""
function redundancy_ppdvspupildiameter()
    x = 20:2:45
    y = 3.0:.05:4

    RGB = true
    Plots.plot(Plots.contour(x, y, (x, y) -> pixelredundancy((50, 35), 18mm, (10mm, 6mm), y * mm, x, RGB=RGB), fill=true, xlabel="pixels per degree", ylabel="pupil diameter", legendtitle="pixel redundancy", title="$(label(RGB)) lenslets"))
end
export redundancy_ppdvspupildiameter

"""computes lenslet display size to match the design constraints"""
function lensletdisplaysize(fov, eyerelief::Unitful.Length, eyebox::NTuple{2,Unitful.Length}, pupildiameter::Unitful.Length, ppd; RGB=true)
    lensprops = defaultclusterproperties()
    angles = lensletangles(eyerelief, eyebox, pupildiameter, ppd, RGB=RGB)
    return @. angles * ppd * lensprops.pixelpitch
end
export lensletdisplaysize

testlensletdisplaysize() = lensletdisplaysize((55, 35), 18mm, (10mm, 6mm), 4mm, 30, RGB=true)
export testlensletdisplaysize

"""generates a contour plot of lenslet display size as a function of ppd and pupil diameter"""
function displaysize_ppdvspupildiameter()
    x = 20:2:45
    y = 3.0:.05:4
    RGB = true

    Plots.plot(Plots.contour(x, y, (x, y) -> maximum(ustrip.(μm, lensletdisplaysize((50, 35), 18mm, (10mm, 6mm), y * mm, x, RGB=RGB))), fill=true, xlabel="pixels per degree", ylabel="pupil diameter", legendtitle="display size μm", title="$(label(RGB)) lenslets"))
end
export displaysize_ppdvspupildiameter

""" Compute high level properties of a lenslet display system. Uses basic geometric and optical constraints to compute approximate values which should be within 10% or so of a real physical system.

`eyerelief` is in mm: `18mm` not `18`.

`eyebox` is a 2 tuple of lengths, also in mm, that specifies the size of the rectangular eyebox, assumed to lie on vertex of the cornea when the eye is looking directly forward.

`fov` is the field of view of the display as seen by the eye, in degrees (no unit type necessary here).

`pupildiameter` is the diameter of the eye pupil in mm.

`mtf` is the MTF response of the system at the angular frequency `cyclesperdegree`.

`pixelsperdegree` is the number of pixels per degree on the display, as seen by the eye through the lenslets. This is not a specification of the optical resolution of the system (use `mtf` and `cyclesperdegree` for that). It is a physical characteristic of the display.

Example:
```
julia> systemproperties(18mm,(10mm,9mm),(55°,45°),4.0mm,.2,11,30)
(lenslet_diameter = 0.7999999999999999 mm, diffraction_limit = 26.344592482933276, fnumber = 1.9337325040249893, focal_length = 1.5469860032199914 mm, display_size = (18.740413819862866 mm, 14.911688245431423 mm), lenslet_display_size = (261.4914368916943 μm, 358.6281908905529 μm), total_silicon_area = 52.1360390897504 mm^2, number_lenslets = 555.9505147668694, pixel_redundancy = 28.895838544429424 °^-2, eyebox_angles = (29.054604099077146, 26.56505117707799), lenslet_fov = (9.684868033025715, 13.282525588538995), subdivisions = (3, 2))
```
"""
function systemproperties(eyerelief::Unitful.Length, eyebox::NTuple{2,Unitful.Length}, fov, pupildiameter::Unitful.Length, mtf, cyclesperdegree,pixelsperdegree; minfnumber=2.0,RGB=true,λ=530nm,pixelpitch=.9μm)
    diameter = diameter_for_cycles_deg(mtf, cyclesperdegree, λ)
    clusterdata = choosecluster(pupildiameter, diameter)
    difflimit = diffractionlimit(λ, clusterdata.lensletdiameter)
    numlenses = numberoflenslets(fov, eyerelief, clusterdata.lensletdiameter)
    redundancy = pixelredundancy(fov, eyerelief, eyebox, pupildiameter, difflimit, RGB=RGB)
    subdivisions = anglesubdivisions(pupildiameter, λ, mtf, cyclesperdegree, RGB=RGB)
    eyebox_angles = eyeboxangles(eyebox,eyerelief)
    angles = lensletangles(eyerelief, eyebox, pupildiameter, difflimit, clusterproperties=(mtf = mtf, minfnumber = minfnumber, cyclesperdegree = cyclesperdegree, λ = λ, pixelpitch = pixelpitch))
    dispsize = lensletdisplaysize(angles, eyerelief, eyebox, pupildiameter, pixelsperdegree, RGB=RGB)
    focal_length = focallength(pixelpitch,1/pixelsperdegree)
    fnumber = focal_length/clusterdata.lensletdiameter
    siliconarea = uconvert(mm^2,numlenses  * dispsize[1]*dispsize[2])
    fulldisplaysize = sizeofdisplay(fov,eyerelief)

    return (lenslet_diameter = clusterdata.lensletdiameter, diffraction_limit = difflimit, fnumber = fnumber, focal_length = focal_length, display_size = fulldisplaysize, lenslet_display_size = dispsize, total_silicon_area = siliconarea, number_lenslets = numlenses, pixel_redundancy = redundancy, eyebox_angles = eyebox_angles, lenslet_fov = angles, subdivisions = subdivisions)
end
export systemproperties

"""prints system properties nicely"""
function printsystemproperties(eyerelief::Unitful.Length, eyebox::NTuple{2,Unitful.Length}, fov, pupildiameter::Unitful.Length, mtf, cyclesperdegree,pixelsperdegree; minfnumber=2.0,RGB=true,λ=530nm,pixelpitch=.9μm)
    println("eye relief = $eyerelief")
    println("eye box = $eyebox")
    println("fov = $(fov)°")
    println("pupil diameter = $pupildiameter")
    println("mtf = $mtf @ $cyclesperdegree cycles/degree")
    for (key,value) in pairs(systemproperties(eyerelief, eyebox, fov, pupildiameter, mtf, cyclesperdegree,pixelsperdegree, minfnumber = minfnumber,RGB=RGB,λ=λ,pixelpitch=pixelpitch))
        if key == :diffraction_limit
            println("$key = $value cycles/°")
        else
            println("$key = $value")
        end
    end
end
export printsystemproperties


text:
using ChangePointMean
using Documenter

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

makedocs(;
    modules=[ChangePointMean],
    authors="Daymond Ling",
    repo="https://github.com/DaymondLing/ChangePointMean.jl/blob/{commit}{path}#{line}",
    sitename="ChangePointMean.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://DaymondLing.github.io/ChangePointMean.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "User's Guide" => [
            "Single change" => "man/singlechg.md",
            "Multiple changes" => "man/multichg.md",
            "Plot routine" => "man/plot.md",
            ],
        "Function Reference" => "Reference.md",
    ],
)

deploydocs(;
    repo = "github.com/DaymondLing/ChangePointMean.jl.git",
    devbranch = "main"
)


text:
module OrderedCollections

    import Base: <, <=, ==, convert, length, isempty, iterate, delete!,
                 show, dump, empty!, getindex, setindex!, get, get!,
                 in, haskey, keys, merge, copy, cat,
                 push!, pop!, insert!,
                 union!, delete!, empty, sizehint!,
                 isequal, hash,
                 map, map!, reverse,
                 first, last, eltype, getkey, values, sum,
                 merge, merge!, lt, Ordering, ForwardOrdering, Forward,
                 ReverseOrdering, Reverse, Lt,
                 isless,
                 union, intersect, symdiff, setdiff, setdiff!, issubset,
                 searchsortedfirst, searchsortedlast, in,
                 filter, filter!, ValueIterator, eachindex, keytype,
                 valtype, lastindex, nextind,
                 copymutable, emptymutable, dict_with_eltype

    export OrderedDict, OrderedSet, LittleDict
    export freeze

    include("dict_support.jl")
    include("ordered_dict.jl")
    include("little_dict.jl")
    include("ordered_set.jl")
    include("dict_sorting.jl")

    import Base: similar
    @deprecate similar(d::OrderedDict) empty(d)
    @deprecate similar(s::OrderedSet) empty(s)
end


text:
module DesignModels

## Import Automa.jl for parsing stl strings
import Automa
import Automa.RegExp: @re_str
# Define the notation for RegExp in Automa.jl
const re = Automa.RegExp
# Concatenation with spaces: \times (tab completion)
× = (x::re.RE, y::re.RE) -> x * re"[\t\n\r ]+" * y

# Import JSON3
import JSON3
import JSON3.StructTypes

# Exports
export import_design, export_design, convert_design

# Includes
include("design.jl")
include("stl.jl")
include("io.jl")

"""
    import_design(path::String; format = :stl)
Import a design to internal format. Default to `.stl` designs.
"""
function import_design(path::String; format=:stl)
    if format == :stl
        return parse_stl(read(path, String))[1]
    end
end

function export_design(design::Design; format=:json)
    if format == :json
        return JSON3.write(design)
    end
end

function convert_design(path::String; from=:stl, to=:json, file="")
    str = export_design(import_design(path; format=from); format=to)
    if file != ""
        write(file, str)
    end
    return str
end

end # module


text:
if VERSION >= v"0.7.0-DEV.2437"
    using Base.Meta.parse
end

const VERSIONS = Dict("Julia" => string(VERSION),
                      "BenchmarkTools" => string(BENCHMARKTOOLS_VERSION))

# TODO: Add any new types as they're added
const SUPPORTED_TYPES = [Benchmark, BenchmarkGroup, Parameters, TagFilter, Trial,
                         TrialEstimate, TrialJudgement, TrialRatio]

for T in SUPPORTED_TYPES
    @eval function JSON.lower(x::$T)
        d = Dict{String,Any}()
        for i = 1:nfields(x)
            name = String(fieldname($T, i))
            field = getfield(x, i)
            value = typeof(field) in SUPPORTED_TYPES ? JSON.lower(field) : field
            push!(d, name => value)
        end
        [string(typeof(x)), d]
    end
end

function recover(x::Vector)
    length(x) == 2 || throw(ArgumentError("Expecting a vector of length 2"))
    typename = x[1]::String
    fields = x[2]::Dict
    T = eval(parse(typename))::Type
    fc = fieldcount(T)
    xs = Vector{Any}(uninitialized, fc)
    for i = 1:fc
        ft = fieldtype(T, i)
        fn = String(fieldname(T, i))
        xs[i] = if ft in SUPPORTED_TYPES
            recover(fields[fn])
        else
            convert(ft, fields[fn])
        end
        if T == BenchmarkGroup && xs[i] isa Dict
            for (k, v) in xs[i]
                if v isa Vector && length(v) == 2 && v[1] isa String
                    xs[i][k] = recover(v)
                end
            end
        end
    end
    T(xs...)
end

function badext(filename)
    noext, ext = splitext(filename)
    msg = if ext == ".jld"
        "JLD serialization is no longer supported. Benchmarks should now be saved in\n" *
        "JSON format using `save(\"$noext\".json, args...)` and loaded from JSON using\n" *
        "using `load(\"$noext\".json, args...)`. You will need to convert existing\n" *
        "saved benchmarks to JSON in order to use them with this version of BenchmarkTools."
    else
        "Only JSON serialization is supported."
    end
    throw(ArgumentError(msg))
end

function save(filename::AbstractString, args...)
    endswith(filename, ".json") || badext(filename)
    isempty(args) && throw(ArgumentError("Nothing to save"))
    goodargs = Any[]
    for arg in args
        if arg isa String
            warn("Naming variables in serialization is no longer supported.\nThe name " *
                 "will be ignored and the object will be serialized in the order it appears " *
                 "in the input.")
            continue
        elseif !any(T->arg isa T, SUPPORTED_TYPES)
            throw(ArgumentError("Only BenchmarkTools types can be serialized."))
        end
        push!(goodargs, arg)
    end
    isempty(goodargs) && error("Nothing to save")
    open(filename, "w") do io
        JSON.print(io, [VERSIONS, goodargs])
    end
end

function load(filename::AbstractString, args...)
    endswith(filename, ".json") || badext(filename)
    if !isempty(args)
        throw(ArgumentError("Looking up deserialized values by name is no longer supported, " *
                            "as names are no longer saved."))
    end
    parsed = open(JSON.parse, filename, "r")
    if !isa(parsed, Vector) || length(parsed) != 2 || !isa(parsed[1], Dict) || !isa(parsed[2], Vector)
        error("Unexpected JSON format. Was this file originally written by BenchmarkTools?")
    end
    versions = parsed[1]::Dict
    values = parsed[2]::Vector
    map!(recover, values, values)
end


text:
"""
    node_clustering(s,v)

Return the clustering coefficient of the given node v in stream s:
```math
cc \\left(v \\right) = \\frac{ \\sum_{uw \\in V \\otimes V} \\left| T_{vu} \\cap T_{vw} \\cap T_{uw} \\right| }{\\sum_{uw \\in V \\otimes V} \\left| T_{vu} \\cap T_{vw} \\right|}
```
Note: If node v never has two neighbors present at the same
time in the stream, this function returns 0.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function node_clustering(s::AbstractUndirectedStream, v::AbstractString)
    nomin=0
    denom=0
    N=Set(AbstractString[])
    for n in keys(neighborhood(s,v))
        push!(N,n)
    end
    for (u,w) in N ⊗ N
        nomin+=length((times(s,v,u) ∩ times(s,v,w)) ∩ times(s,u,w))
        denom+=length(times(s,v,u) ∩ times(s,v,w))
    end
    denom != 0 ? nomin/denom : 0.0
end

"""
    node_clustering(s,v,t)

Return the instantaneous clustering coefficient of the given node v
at the given time t in the given stream:
```math
cc_t \\left(v \\right)= \\frac{\\sum_{uw} vu_t vw_t uw_t}{\\sum_{uw} vu_t vw_t}
```
Note: If node v does not have two neighbors at the given
time t, this function returns 0.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function node_clustering(s::AbstractUndirectedStream,v::AbstractString,t::Float64)
    Nt=Set(AbstractString[])
    for n in neighborhood(s,v,t)
        push!(Nt,n)
    end
    nom=sum([1.0 for (u,w) in Nt ⊗ Nt if ((t ∈ times(s,v,u)) & (t ∈ times(s,v,w)) & (t ∈ times(s,u,w)))])
    denom=sum([1.0 for (u,w) in Nt ⊗ Nt if ((t ∈ times(s,v,u)) & (t ∈ times(s,v,w)))])
    denom != 0 ? nom/denom : 0.0
end

"""
    node_clustering(s)

Return the node clustering coefficient of the given stream:
```math
cc \\left( V \\right)= \\frac{1}{n} \\sum_{v \\in V} n_v cc \\left( v \\right) = \\sum_{v \\in V} \\frac{\\left| T_v \\right|}{\\left|W \\right|} cc \\left(v \\right)
```
Where n is the number of nodes in the stream.
Note: If there is no node in the stream, this function returns 0.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function node_clustering(s::AbstractUndirectedStream)
    if length(s.V)>0
        1.0/length(s.V)*sum([contribution(s,v)*node_clustering(s,v) for v in s.V])
    else
        0.0
    end
end

"""
    time_clustering(s,t)

Return time clustering coefficient of a given time instant in the stream:
```math
cc \\left( t \\right)= \\frac{\\sum_{v}cc_t\\left(v \\right) \\sum_{uw} vu_t vw_t}{\\sum_{v} \\sum_{uw} vu_t vw_t}
```

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function time_clustering(s::AbstractUndirectedStream,t::Float64)
    nom::Float64=0.0
    denom::Float64=0.0
    for v in s.V
        acc::Float64=0.0
        for (u,w) in s.V ⊗ s.V
            if ((t ∈ times(s,v,u)) & (t ∈ times(s,v,w)))
                acc+=1.0
            end
        end
        nom+=node_clustering(s,v,t) * acc
        denom+=acc
    end
    denom != 0 ? nom/denom : 0.0
end

"""
    time_clustering(s)

Return the time clustering coefficient of the given stream:
```math
cc \\left( T \\right) = \\int_{t} \\frac{\\left| V_t \\right|}{\\left| W \\right|} cc \\left( t \\right) dt
```
Note: If nodes are never present in the stream, this function returns 0.

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function time_clustering(s::AbstractUndirectedStream)
    τ = times(s)
    dW = duration(s.W)
    dW != 0 ? 1.0/dW*sum([time_clustering(s,0.5*(t[2]+t[1])) * length(nodes(s,0.5*(t[2]+t[1]))) for t in zip(τ[1:end-1],τ[2:end])]) : 0.0
end

"""
    time_clustering(s)

Return the time clustering coefficient of the given stream using a TimeCursor.
"""
function time_clustering(s::AbstractUndirectedStream,tc::TimeCursor)
    throw("Not Implemented")
end

"""
    clustering(s)

Return the clustering coefficient of the given stream:
```math
cc \\left( S \\right) = \\int_t \\frac{1}{\\left|T \\right|} \\sum_{v} \\frac{cc_t\\left(v \\right)}{\\left|V \\right|} dt
```

### Reference
- Matthieu Latapy, Tiphaine Viard and Clémence Magnien Social
  Networks Analysis and Mining, 8: 61, 2018. "Stream Graphs and
  Link Streams for the Modeling of Interactions over Time".
  [(arXiv)](https://arxiv.org/pdf/1710.04073.pdf)
"""
function clustering(s::AbstractUndirectedStream)
    card_TV = duration(s)*length(s.V)
    if card_TV==0
        return 0.0
    end
    τ = times(s)
    acc::Float64 = 0.0
    for v in s.V
        for t in zip(τ[1:end-1],τ[2:end])
            acc+=node_clustering(s,v,0.5*(t[1]+t[2]))*(t[2]-t[1])
        end
    end
    1.0/card_TV*acc
end


text:
# ComponentSol{T,N,C} = SciMLBase.AbstractODESolution{T,N,C} where C<:AbstractVector{<:ComponentArray}

# Plotting stuff
function SciMLBase.getsyms(sol::SciMLBase.AbstractODESolution{T,N,C}) where {T,N,C<:AbstractVector{<:ComponentArray}}
    if SciMLBase.has_syms(sol.prob.f)
        return sol.prob.f.syms
    else
        return labels(sol.u[1])
    end
end

# A little bit of type piracy. Should probably make this a PR to DiffEqBase
SciMLBase.cleansyms(syms::AbstractArray{<:String}) = SciMLBase.cleansyms.(syms)
SciMLBase.cleansyms(syms::String) = syms