text:
for f in (:sum, :prod, :maximum, :minimum, :extrema, :all, :any, :findmax)
    @eval begin
        function Base.$(f)(a::IndicesArray; dims=:, kwargs...)
            d = to_dims(a, dims)
            return rebuild(
                a,
                Base.$(f)(parent(a); dims=d, kwargs...),
                reduce_axes(a, d)
               )
        end
    end
end

for f in (:cumsum, :cumprod, :sort, :sort!)
    @eval begin
        function Base.$(f)(a::IndicesArray; dims, kwargs...)
            return rebuild(a, $(f)(parent(a), dims=to_dims(a, dims), kwargs...), axes(a))
        end

        function Base.$(f)(a::IndicesVector; kwargs...)
            return rebuild(a, $(f)(parent(a); kwargs...), axes(a))
        end
    end
end

for (A,B,fa,fb) in ((:AbstractMatrix, :IndicesMatrix,  identity, parent),
                    (:AbstractMatrix, :IndicesVector,  identity, parent),
                    (:IndicesMatrix,  :AbstractVector, parent,   parent),
                    (:IndicesMatrix,  :AbstractMatrix, parent,   identity),
                    (:IndicesVector,  :AbstractMatrix, parent,   identity),
                    (:AbstractVector, :IndicesMatrix,  identity, parent),
                    (:IndicesMatrix,  :IndicesVector,  parent,   parent),
                    (:IndicesMatrix,  :IndicesMatrix,  parent,   parent),
                    (:IndicesVector,  :IndicesMatrix,  parent,   parent))
    @eval begin
        function Base.:*(a::$A, b::$B)
            rebuild(rebuild_rule(a, b), *($(fa)(a), $(fb)(b)), matmul_axes(a, b))
        end
    end
end

for A in (Adjoint{<:Any, <:AbstractVector}, Transpose{<:Real, <:AbstractVector{<:Real}})
    @eval begin
        Base.:*(a::$A, b::IndicesVector) where {T,I} = *(a, parent(b))
        function Base.:*(b::IndicesVector, a::$A) where {T,I}
            return rebuild(b, *(parent(b), a), matmul_axes(b, a))
        end
    end
end

# Between arrays
function Base.:+(a::IndicesArray, b::IndicesArray)
    return rebuild(
        promote_rule(typeof(a), typeof(b)),
        +(parent(a), parent(b)),
        combine_axes(a, b)
       )
end
function Base.:+(a::AbstractArray, b::IndicesArray)
    return rebuild(
        promote_rule(typeof(a), typeof(b)),
        +(a, parent(b)),
        combine_axes(a, b)
       )
end
function Base.:+(a::IndicesArray, b::AbstractArray)
    return rebuild(
        promote_rule(typeof(a), typeof(b)),
        +(parent(a), b),
        combine_axes(a, b)
       )
end

function Base.:-(a::IndicesArray, b::IndicesArray)
    return rebuild(
        promote_rule(typeof(a), typeof(b)),
        -(parent(a), parent(b)),
        combine_axes(a, b)
       )
end
function Base.:-(a::AbstractArray, b::IndicesArray)
    return rebuild(
        promote_rule(typeof(a), typeof(b)),
        -(a, parent(b)),
        combine_axes(a, b)
       )
end
function Base.:-(a::IndicesArray, b::AbstractArray)
    return rebuild(
        promote_rule(typeof(a), typeof(b)),
        -(parent(a), b),
        combine_axes(a, b)
       )
end
function Base.:-(a::IndicesArray, b::IndicesArray, c...)
    return -(rebuild(a, -(parent(a), parent(b)), combine_indices(a, b)), c...)
end
function Base.:-(a::AbstractArray, b::IndicesArray, c...)
    return -(rebuild(a, -(a, parent(b)), axes(b)), c...)
end
function Base.:-(a::IndicesArray, b::AbstractArray, c...)
    return -(rebuild(a, -(parent(a), b), axes(a)), c...)
end

Base.:*(a::Number, b::IndicesArray) = rebuild(b, *(a, parent(b)), axes(b))
Base.:*(a::IndicesArray, b::Number) = rebuild(a, *(parent(a), b), axes(a))

Base.:\(a::Number, b::IndicesArray) = rebuild(b, \(a, parent(b)), axes(b))
Base.:/(a::IndicesArray, b::Number) = rebuild(a, /(parent(a), b), axes(a))


text:
using DrWatson, Test

# Test commit function
com = current_commit(@__DIR__)
@test com === nothing

com = current_commit(dirname(@__DIR__))
@test com !== nothing
@test typeof(com) == String

# tag!
d1 = Dict(:x => 3, :y => 4)
d2 = Dict("x" => 3, "y" => 4)
for d in (d1, d2)
    d = tag!(d, dirname(@__DIR__))

    @test haskey(d, keytype(d)(:commit))
    @test d[keytype(d)(:commit)] |> typeof == String
end

# @tag!
for d in (d1, d2)
    d = @tag!(d, @__DIR__)
    @test !haskey(d, keytype(d)(:commit))

    d = @tag!(d, dirname(@__DIR__))
    @test d[keytype(d)(:commit)] |> typeof == String
    @test d[keytype(d)(:script)][1:4] == "test"
end

# Test dictionary expansion
c = Dict(:a => [1, 2], :b => 4);
c1 = [ Dict(:a=>1,:b=>4)
 Dict(:a=>2,:b=>4)]

v1 = dict_list(c)
for el in c1
    @test el ∈ v1
end
@test keytype(eltype(v1)) == Symbol

c[:c] = "test"; c[:d] = ["lala", "lulu"];
c2 = [ Dict(:a=>1,:b=>4,:d=>"lala",:c=>"test")
 Dict(:a=>2,:b=>4,:d=>"lala",:c=>"test")
 Dict(:a=>1,:b=>4,:d=>"lulu",:c=>"test")
 Dict(:a=>2,:b=>4,:d=>"lulu",:c=>"test")]

v2 = dict_list(c)
for el in c2
    @test el ∈ v2
end
@test keytype(eltype(v2)) == Symbol

c[:e] = [[1, 2], [3, 5]];
c3 = [
Dict(:a=>1,:b=>4,:d=>"lala",:e=>[1, 2],:c=>"test")
Dict(:a=>2,:b=>4,:d=>"lala",:e=>[1, 2],:c=>"test")
Dict(:a=>1,:b=>4,:d=>"lulu",:e=>[1, 2],:c=>"test")
Dict(:a=>2,:b=>4,:d=>"lulu",:e=>[1, 2],:c=>"test")
Dict(:a=>1,:b=>4,:d=>"lala",:e=>[3, 5],:c=>"test")
Dict(:a=>2,:b=>4,:d=>"lala",:e=>[3, 5],:c=>"test")
Dict(:a=>1,:b=>4,:d=>"lulu",:e=>[3, 5],:c=>"test")
Dict(:a=>2,:b=>4,:d=>"lulu",:e=>[3, 5],:c=>"test")
]

v3 = dict_list(c)
@test dict_list_count(c) == length(c3) == length(v3)
for el in c3
    @test el ∈ v3
end
@test keytype(eltype(v3)) == Symbol

v4 = dict_list(Dict(:a => 1, :b => 2.0)) # both non-iterable
@test keytype(eltype(v4)) == Symbol

v5 = dict_list(Dict(:a => [1], :b => 2.0)) # one non-iterable
@test keytype(eltype(v5)) == Symbol

v6 = dict_list(Dict(:a => [1], :b => [2.0])) # both iterable
@test keytype(eltype(v6)) == Symbol

### tmpsave ###
tmpdir = joinpath(@__DIR__, "tmp")
ret = tmpsave(v3, tmpdir)
for r in ret
    @test isfile(joinpath(tmpdir, r))
    a = load(joinpath(tmpdir, r))
    @test a ∈ v3
end
rm(tmpdir, force = true, recursive = true)
@test !isdir(tmpdir)

ext: jl


text:
# Just a work on the

Ae = [3.0 1.0; 1.0 3.0]
Ae2 = [3.0 1.0; -1.0 3.0]
function test_sym_non_per(Ae, nel, nb, nbslv)

    m = zeros(Int, 2, nel)
    for e in 1:nel
        m[1, e] = e
        m[2, e] = e+1
    end

    Ag = HPFEM.BBMatrix{Float64}(nb, nbslv)

    # Assemble the matrix:
    for e = 1:nel
        HPFEM.assemble!(Ag, Ae, m[:,e])
    end
    A = copy(Ag.A)

    HPFEM.trf!(Ag)

    f = rand(nbslv)

    x0 = A\copy(f)

    x1 = HPFEM.trs!(Ag, copy(f))

    Ag_st = HPFEM.BBSymTri{Float64}(nb, nbslv)
    for e = 1:nel
        HPFEM.assemble!(Ag_st, Ae, m[:,e])
    end

    HPFEM.trf!(Ag_st)
    x2 = HPFEM.trs!(Ag_st, copy(f))

    Ag_tri = HPFEM.BBTri{Float64}(nb, nbslv)

    for e = 1:nel
        HPFEM.assemble!(Ag_tri, Ae, m[:,e])
    end

    HPFEM.trf!(Ag_tri)
    x3 = HPFEM.trs!(Ag_tri, copy(f))

    Ag_tp = HPFEM.BBTriP{Float64}(nb, nbslv, false)

    for e = 1:nel
        HPFEM.assemble!(Ag_tp, Ae, m[:,e])
    end
    HPFEM.trf!(Ag_tp)
    x4 = HPFEM.trs!(Ag_tp, copy(f))

    hcat(x0, x1-x0, x2-x0, x3-x0, x4-x0)
end

function test_non_per(Ae, nel, nb, nbslv)

    m = zeros(Int, 2, nel)
    for e in 1:nel
        m[1, e] = e
        m[2, e] = e+1
    end

    Ag = HPFEM.BBMatrix{Float64}(nb, nbslv)

    # Assemble the matrix:
    for e = 1:nel
        HPFEM.assemble!(Ag, Ae, m[:,e])
    end
    A = copy(Ag.A)

    HPFEM.trf!(Ag)

    f = rand(nbslv)

    x0 = A\copy(f)

    x1 = HPFEM.trs!(Ag, copy(f))

    Ag_tri = HPFEM.BBTri{Float64}(nb, nbslv)

    for e = 1:nel
        HPFEM.assemble!(Ag_tri, Ae, m[:,e])
    end

    HPFEM.trf!(Ag_tri)
    x2 = HPFEM.trs!(Ag_tri, copy(f))

    Ag_tp = HPFEM.BBTriP{Float64}(nb, nbslv, false)

    for e = 1:nel
        HPFEM.assemble!(Ag_tp, Ae, m[:,e])
    end
    HPFEM.trf!(Ag_tp)
    x3 = HPFEM.trs!(Ag_tp, copy(f))

    hcat(x0, x1-x0, x2-x0, x3-x0)
end

function test_per(Ae, nel)

    m = zeros(Int, 2, nel)
    for e in 1:nel
        m[1, e] = e
        m[2, e] = e+1
    end
    m[2,nel] = 1
    nb = nel
    nbslv = nb
    Ag = HPFEM.BBMatrix{Float64}(nb, nbslv)

    # Assemble the matrix:
    for e = 1:nel
        HPFEM.assemble!(Ag, Ae, m[:,e])
    end
    A = copy(Ag.A)

    HPFEM.trf!(Ag)

    f = rand(nbslv)

    x0 = A\copy(f)

    x1 = HPFEM.trs!(Ag, copy(f))

    Ag_tp = HPFEM.BBTriP{Float64}(nb, nbslv, true)

    for e = 1:nel
        HPFEM.assemble!(Ag_tp, Ae, m[:,e])
    end
    HPFEM.trf!(Ag_tp)
    x2 = HPFEM.trs!(Ag_tp, copy(f))

    hcat(x0, x1-x0, x2-x0)
end

function test_sym_per(Ae, nel)

    m = zeros(Int, 2, nel)
    for e in 1:nel
        m[1, e] = e
        m[2, e] = e+1
    end
    m[2,nel] = 1
    nb = nel
    nbslv = nb
    Ag = HPFEM.BBMatrix{Float64}(nb, nbslv)

    # Assemble the matrix:
    for e = 1:nel
        HPFEM.assemble!(Ag, Ae, m[:,e])
    end
    A = copy(Ag.A)

    HPFEM.trf!(Ag)

    f = rand(nbslv)

    x0 = A\copy(f)

    x1 = HPFEM.trs!(Ag, copy(f))

    Ag_tp = HPFEM.BBTriP{Float64}(nb, nbslv, true)

    for e = 1:nel
        HPFEM.assemble!(Ag_tp, Ae, m[:,e])
    end
    HPFEM.trf!(Ag_tp)
    x2 = HPFEM.trs!(Ag_tp, copy(f))

    Ag_sp = HPFEM.BBSymTriP{Float64}(nb, nbslv, true)

    for e = 1:nel
        HPFEM.assemble!(Ag_sp, Ae, m[:,e])
    end
    HPFEM.trf!(Ag_sp)
    x3 = HPFEM.trs!(Ag_sp, copy(f))

    hcat(x0, x1-x0, x2-x0, x3-x0)
end


text:
module TestProgressMacro

using ProgressLogging: @progress, ProgressLevel
using Test
using Test: collect_test_logs
using OffsetArrays

@testset "@progress" begin
  let i = 0, x
      x = @progress for _ = 1:100
          i += 1
      end
      @test i == 100
      @test x == nothing
  end

  let i = 0, r = -50:10:50, x
      x = @progress for _ in r
          i += 1
      end
      @test i == 11
      @test x == nothing
  end

  let i = 0, x
      x = @progress "named" for _ = 1:100
          i += 1
      end
      @test i == 100
      @test x == nothing
  end

  let i = 0, j = 0, x
      x = @progress for _ = 1:10, __ = 1:20
          i += 1
      end
      @test i == 200
      @test x == nothing
  end

  let i = 0, j = 0, x
      bar = "bar"
      x = @progress "foo $bar" for _ = 1:10
          i += 1
      end
      @test i == 10
      @test x == nothing
  end

  let x, y
      x = @progress y = [i + 3j for i = 1:3, j = 1:4]
      @test y == reshape(4:15, 3, 4)
      @test x == y
  end

  let a = [], x
      x = @progress for i = 1:3, j in [-5, -2, -1, 8]
          j > 0 && continue
          push!(a, (i, j))
          i > 1 && break
      end
      @test a == [(1, -5), (1, -2), (1, -1), (2, -5)]
      @test x == nothing
  end

  # Multi-dimensional arrays in comprehension and offset axes
  let off1 = -2, off2 = 21
    v1 = OffsetArray(2:3, off1)
    v2 = OffsetArray(-1:2, off2)
    logs, _ = collect_test_logs(min_level = ProgressLevel) do
      x = @progress y = [i*j for i in v1, j in v2]
      @test x == y == OffsetArray([-2 0 2 4; -3 0 3 6], off1, off2)
    end
    @test isequal(
        [l.kwargs[:progress] for l in logs],
        [nothing; (1:8)./8; "done"],
    )
    m = OffsetArray(reshape(2:7,2,3), off1, off2)
    x = @progress y = [i*j for i in v1, j in m]
    @test x == y == [i*j for i in v1, j in m]
  end

  # non-indexable iterables with axes
  @testset "non-indexable" for off in (0,10)
    let r = OffsetVector(1:5, off)
      x1 = @progress y1 = [i for i in (x^2 for x in r)]
      x2 = @progress y2 = [i for i in zip(r,r)]
      @test x1 == y1 == r.^2
      @test x2 == y2 == collect(zip(r,r))

      y1, y2 = [], []
      x1 = @progress for i in (x^2 for x in r)
        push!(y1, i)
      end
      x2 = @progress for i in zip(r,r)
        push!(y2, i)
      end
      @test x1 == x2 == nothing
      @test OffsetVector(y1,off) == r.^2
      @test OffsetVector(y2,off) == collect(zip(r,r))
    end
  end
end

end  # module


text:
function colNullities(P::AbstractArray, tol=0)
    A, B = toPencil(P)
    return colNullities(A, B, tol)
end

"""
    e, d, Ar, Br = colNullities(A, B, tol=0)
computes the column nullities of the pencil Bx-A.
e = column nullities
d =  nullities of infinite eigenvalue
x*Br - Ar = reduced pencil containing only finite eigenvalues and row indices
"""
function colNullities(A::AbstractArray{T, N}, B::AbstractArray{T, N}, tol=0) where {T<:Number, N}

    r = rank(B, atol=tol)
    s = size(B, 2) - r
    blocks  = [s]
    ranks = Int[]

    while s != 0
        F = svd(B, full=true);
        B = B*F.V
        A = A*F.V

        F = svd(A[:, r+1:end], full=true)
        Ua = F.U
        ra = rank(A[:, r+1:end], atol=tol)
        push!(ranks, ra)
        P = Matrix(I, size(Ua))
        P = P[[ra+1:end; 1:ra], :]
        A = P*Ua'*A
        B = P*Ua'*B

        A = A[1:end-ra, 1:r]
        B = B[1:end-ra, 1:r]
        r = rank(B, atol=tol)
        s = size(B, 2) - r
        push!(blocks, s)
    end
    e = blocks[1:end-1] - ranks
    d = ranks - blocks[2:end]

    removeTailZeroes!(e)
    removeTailZeroes!(d)
    iszero(d) && (d = Int[])
    iszero(e) && (e = Int[])
    return e, d, A, B
end

function rowNullities(P::AbstractArray, tol=0)
    A, B = toPencil(P)
    return rowNullities(A, B, tol)
end

"""
    e, d, Ar, Br = colNullities(A, B, tol=0)
computes the row nullities of the pencil Bx-A.
e = row nullities
d =  nullities of infinite eigenvalue
x*Br - Ar = reduced pencil containing only finite eigenvalues and column indices
"""
function rowNullities(A::AbstractArray{T, N}, B::AbstractArray{T, N}, tol=0) where {T<:Number, N}

    r = rank(B, atol=tol)
    s = size(B, 1) - r
    blocks = [s]
    ranks = Int[]
    while s != 0
        F = svd(B, full=true)
        Ub = F.U
        P = Matrix(I, size(Ub))
        P = P[[r+1:end; 1:r], :]
        B = P*Ub'*B
        A = P*Ub'*A
        F = svd(A[1:s, :], full=true)
        ra = rank(A[1:s, :], atol=tol)
        push!(ranks, ra)
        A = A*F.V
        B = B*F.V
        A = A[s+1:end, ra+1:end]
        B = B[s+1:end, ra+1:end]
        r = rank(B, atol=tol)
        s = size(B, 1) - r
        push!(blocks, s)
    end

    e = blocks[1:end-1] - ranks
    d = ranks - blocks[2:end]

    removeTailZeroes!(e)
    removeTailZeroes!(d)
    iszero(d) && (d = Int[])
    iszero(e) && (e = Int[])
    return e, d, A, B
end

function eigNullities(λ::Number, A::AbstractArray{T, N}, B::AbstractArray{T, N}, tol=0) where {T<:Number, N}
    A = A - λ*B
    r = rank(A, atol=tol)
    s = size(A, 2) - r

    blocks  = [s]

    while s != 0
        F = svd(A, full=true);
        B = B*F.V
        A = A*F.V

        F = svd(B[:, r+1:end], full=true)
        Ua = F.U

        P = Matrix(I, size(Ua))
        P = P[[s+1:end; 1:s], :]
        A = P*Ua'*A
        B = P*Ua'*B

        A = A[1:r, 1:r]
        B = B[1:r, 1:r]
        r = rank(A, atol=tol)
        s = size(A, 2) - r
        push!(blocks, s)
    end
    e = blocks[1:end-1] - blocks[2:end]

    return e
end

function eigNullities(λ::Number, P::AbstractArray, tol=0)
    _, _, Ar, Br = colNullities(P, tol)
    _, _, Ar, Br = rowNullities(Ar, Br, tol)
    return eigNullities(λ, Ar, Br, tol)
end


text:
@inline function cos_kernel(::Union{NS,MV}, x::MCNoGrad, z::Interval{Float64},
                            cv_tp1::Float64, cv_tp2::Float64, cc_tp1::Float64, cc_tp2::Float64)
    xL = x.Intv.lo;    xU = x.Intv.hi;    zL = z.lo;    zU = z.hi
    eps_min, eps_max = cos_arg(xL, xU)
    mcc, cci = mid3(x.cc, x.cv, eps_max)
    mcv, cvi = mid3(x.cc, x.cv, eps_min)
    cc, dcc, cc_tp1, cc_tp2 = cc_cos(mcc, xL, xU, cc_tp1, cc_tp2)
    cv, dcv, cv_tp1, cv_tp2 = cv_cos(mcv, xL, xU, cv_tp1, cv_tp2)
    if zL > zcv
        cvi = 3
        dcv = 0.0
    end
    if zU < zcc
        cci = 3
        dcc = 0.0
    end
    return MCNoGrad(cv, cc, y, x.cnst), cvi, cci, dcv, dcc, cv_tp1, cc_tp1, cv_tp1, cc_tp2
end
@inline function cos(::ANYRELAX, x::MCNoGrad)
    y, _, _, _, _, _, _, _, _, = cos_kernel(x, cos(x.Intv), Inf, Inf, Inf, Inf)
    return y
end
function sin_kernel(::ANYRELAX, x::MCNoGrad, y::Interval{Float64},
                    cv_tp1::Float64, cv_tp2::Float64,
                    cc_tp1::Float64, cc_tp2::Float64)
    cos_kernel(x-pi/2.0, y, cv_tp1, cv_tp2, cc_tp1, cc_tp2)
end
function sin(::ANYRELAX, x::MCNoGrad)
    y, _, _, _, _ = sin_kernel(x, sin(x.Intv), Inf, Inf, Inf, Inf)
    return y
end


text:
using Diagonalizations, LinearAlgebra, PosDefManifold, Test

####  Create data for testing the case k=1, m>1
# `t` is the number of samples,
# `m` is the number of datasets,
# `n` is the number of variables,
# `noise` must be smaller than 1.0. The smaller the noise,
#  the more data are correlated.
function getData(t, m, n, noise)
    # create m identical data matrices and rotate them by different
    # random orthogonal matrices V_1,...,V_m
    𝐕=[randU(n) for i=1:m] # random orthogonal matrices
    X=randn(n, t)  # data common to all subjects
    # each subject has this common part plus a random part
    𝐗=[𝐕[i]'*((1-noise)*X + noise*randn(n, t)) for i=1:m]
    return 𝐗
end

function getData(::Type{Complex{T}}, t, m, n, noise) where {T<:AbstractFloat}
    # create m identical data matrices and rotate them by different
    # random orthogonal matrices V_1,...,V_m
    𝐕=[randU(ComplexF64, n) for i=1:m] # random orthogonal matrices
    X=randn(ComplexF64, n, t)  # data common to all subjects
    # each subject has this common part plus a random part
    𝐗=[𝐕[i]'*((1-noise)*X + noise*randn(ComplexF64, n, t)) for i=1:m]
    return 𝐗
end

# REAL data: check that for the case m=2 gMCA gives the same result as MCA
t, m, n, noise = 20, 2, 6, 0.1
Xset=getData(t, m, n, noise)
Cx=(Xset[1]*Xset[1]')/t
Cy=(Xset[2]*Xset[2]')/t
Cxy=(Xset[1]*Xset[2]')/t

gm=gmca(Xset; simple=true)
m=mca(Cxy; simple=true)

@test (m.F[1]'*Cxy*m.F[2]) ≈ (gm.F[1]'*Cxy*gm.F[2])
# the following must be the identity matrix out of a possible sign ambiguity
@test abs.(m.F[1]'*gm.F[1]) ≈ I
@test abs.(m.F[2]'*gm.F[2]) ≈ I

# COMPLEX data: check that for the case m=2 gMCA gives the same result as MCA
t, m, n, noise = 20, 2, 6, 0.1
Xcset=getData(ComplexF64, t, m, n, noise)
Ccx=(Xcset[1]*Xcset[1]')/t
Ccy=(Xcset[2]*Xcset[2]')/t
Ccxy=(Xcset[1]*Xcset[2]')/t

gmc=gmca(Xcset; simple=true)
mc=mca(Ccxy; simple=true)

# for complex data just do a sanity check as the order of vectors
# is arbitrary
@test spForm(mc.F[1]'gmc.F[1])<0.01
@test spForm(mc.F[2]'gmc.F[2])<0.01

# REAL data: m>2 case
t, m, n, noise = 20, 4, 6, 0.1
Xset=getData(t, m, n, noise)

# ... selecting subspace dimension allowing an explained variance = 0.9
gm=gmca(Xset, eVar=0.9)

# name of the filter
gm.name

𝒞=Array{Matrix}(undef, 1, m, m)
for i=1:m, j=1:m 𝒞[1, i, j]=(Xset[i]*Xset[j]')/t end

using Plots
# plot regularized accumulated eigenvalues
plot(gm.arev)

# plot the original cross-covariance matrices and the rotated
# cross-covariance matrices

# Get all products 𝐔[i]' * 𝒞[l, i, j] * 𝐔[j]
function _rotate_crossCov(𝐔, 𝒞, m, k)
    𝒮=Array{Matrix}(undef, k, m, m)
    @inbounds for l=1:k, i=1:m, j=1:m 𝒮[l, i, j]=𝐔[i]'*𝒞[l, i, j]*𝐔[j] end
    return 𝒮
end

# Put all cross-covariances in a single matrix of dimension m*n x m*n for visualization
function 𝒞2Mat(𝒞::AbstractArray, m, k)
    n=size(𝒞[1, 1, 1], 1)
    C=Matrix{Float64}(undef, m*n, m*n)
    for i=1:m, j=1:m, x=1:n, y=1:n C[i*n-n+x, j*n-n+y]=𝒞[k, i, j][x, y] end
    return C
end

 C=𝒞2Mat(𝒞, m, 1)
 Cmax=maximum(abs.(C));
 h1 = heatmap(C, clim=(-Cmax, Cmax), yflip=true, c=:bluesreds, title="all cross-covariances")
 𝒮=_rotate_crossCov(gm.F, 𝒞, m, 1)
 S=𝒞2Mat(𝒮, m, 1)
 Smax=maximum(abs.(S));
 h2 = heatmap(S, clim=(0, Smax), yflip=true, c=:amp, title="all rotated cross-covariances")
 📈=plot(h1, h2, size=(700,300))
# savefig(📈, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FiggMCA.png")

# COMPLEX data: m>2 case
t, m, n, noise = 20, 4, 6, 0.1
Xcset=getData(ComplexF64, t, m, n, noise)

# ... selecting subspace dimension allowing an explained variance = 0.9
gmc=gmca(Xcset, eVar=0.9)


text:
using ColorLab,Test,YAML,Makie,StatsMakie

"Get digital RGB color spectral measured from a specific display"
function RGBSpectral(measurement)
    C = map(i->parse.(Float64,split(i)),measurement["Color"])
    λ = measurement["WL"]
    I = measurement["Spectral"]
    return C,λ,I
end
cd(@__DIR__)

## plot cone fundamentals
Plots.plot(sscone2le[:,1],sscone2le[:,2:end],linewidth=2,color=["Red" "Green" "Blue"],xlabel="Wavelength (nm)",ylabel="Sensitivity",label=["L" "M" "S"],title="Cone Fundamentals(2deg)")
foreach(i->savefig("Cone Fundamentals(2deg)$i"),[".png",".svg"])

Plots.plot(sscone10le[:,1],sscone10le[:,2:end],linewidth=2,color=["Red" "Green" "Blue"],xlabel="Wavelength (nm)",ylabel="Sensitivity",label=["L" "M" "S"],title="Cone Fundamentals(10deg)")
foreach(i->savefig("Cone Fundamentals(10deg)$i"),[".png",".svg"])

# plot rgb color matching funcionts
Plots.plot(sbrgb2[:,1],sbrgb2[:,2:end],linewidth=2,color=["Red" "Green" "Blue"],xlabel="Wavelength (nm)",ylabel="Tristimulus Value",label=["r" "g" "b"],title="Color Matching Functions_rgb(2deg)")
Plots.scatter!(sbrgb_primary',zeros(3)',markersize=8,markerstrokewidth=0,markercolor=["Red" "Green" "Blue"],legend=false)
foreach(i->savefig("Color Matching Functions_rgb(2deg)$i"),[".png",".svg"])

Plots.plot(sbrgb10[:,1],sbrgb10[:,2:end],linewidth=2,color=["Red" "Green" "Blue"],xlabel="Wavelength (nm)",ylabel="Tristimulus Value",label=["r" "g" "b"],title="Color Matching Functions_rgb(10deg)")
Plots.scatter!(sbrgb_primary',zeros(3)',markersize=8,markerstrokewidth=0,markercolor=["Red" "Green" "Blue"],legend=false)
foreach(i->savefig("Color Matching Functions_rgb(10deg)$i"),[".png",".svg"])

# plot xyz color matching funcionts
xyz2 = newcmf(sscone2le,LMSToXYZ2)
Plots.plot(xyz2[:,1],xyz2[:,2:end],linewidth=2,color=["Red" "Green" "Blue"],xlabel="Wavelength (nm)",ylabel="Tristimulus Value",label=["x" "y" "z"],title="Color Matching Functions_xyz(2deg)")
foreach(i->savefig("Color Matching Functions_xyz(2deg)$i"),[".png",".svg"])

xyz10 = newcmf(sscone2le,LMSToXYZ10)
Plots.plot(xyz10[:,1],xyz10[:,2:end],linewidth=2,color=["Red" "Green" "Blue"],xlabel="Wavelength (nm)",ylabel="Tristimulus Value",label=["x" "y" "z"],title="Color Matching Functions_xyz(10deg)")
foreach(i->savefig("Color Matching Functions_xyz(10deg)$i"),[".png",".svg"])

# ASUS ROG Swift PG279Q IPS LCD
displayname = "ROGPG279Q"
# ViewSonic VX3276mhd IPS LCD
displayname = "VX3276mhd"
# Sony Trinitron CRT
displayname = "Trinitron"

resultdir = "./$displayname"
mkpath(resultdir)

display = YAML.load_file("DisplayMeasurement.yaml")
RGBToLMS,LMSToRGB = RGBLMSMatrix(RGBSpectral(display[displayname]["SpectralMeasurement"])...)
RGBToXYZ,XYZToRGB = RGBXYZMatrix(RGBSpectral(display[displayname]["SpectralMeasurement"])...)

# unit cube digital colors
ur=0:0.01:1
uc=[[i,j,k] for i=ur,j=ur,k=ur][:]
ucm=hcat(uc...)
uqcm=quavectors(ucm)

s=Makie.scatter(ucm[1,:],ucm[2,:],ucm[3,:],color=[RGB(i...) for i in uc],markersize=0.01,transparency=true)
record(s,"Unit Color Space.mp4",1:360/5,framerate=12) do i
    rotate_cam!(s,5pi/180,0.0,0.0)
end

# Transformation between specific display RGB and LMS Spaces
lms_rgb = RGBToLMS*uqcm
rgb_lms = LMSToRGB*uqcm

s=Makie.scatter(ucm[1,:],ucm[2,:],ucm[3,:],color=[RGBA(i...,1) for i in uc],markersize=0.01,transparency=true)
Makie.scatter!(lms_rgb[1,:],lms_rgb[2,:],lms_rgb[3,:],color=[RGBA(i...,1) for i in uc],markersize=0.01,transparency=true)
s.center=false
s[Axis][:names,:axisnames]=("R / L","G / M","B / S")
record(s,joinpath(resultdir,"Unit RGB To LMS Space.mp4"),1:360/5,framerate=12) do i
    rotate_cam!(s,5pi/180,0.0,0.0)
end

s=Makie.scatter(ucm[1,:],ucm[2,:],ucm[3,:],color=[RGBA(i...,1) for i in uc],markersize=0.01,transparency=true)
Makie.scatter!(rgb_lms[1,:],rgb_lms[2,:],rgb_lms[3,:],color=[RGBA(i...,1) for i in uc],markersize=0.01,transparency=true)
s.center=false
s[Axis][:names,:axisnames]=("R / L","G / M","B / S")
record(s,joinpath(resultdir,"Unit LMS To RGB Space.mp4"),1:360/5,framerate=12) do i
    rotate_cam!(s,5pi/180,0.0,0.0)
end

# Transformation between specific display RGB and CIE XYZ Spaces
XYZ_rgb = RGBToXYZ*uqcm
rgb_XYZ = XYZToRGB*uqcm

s=Makie.scatter(ucm[1,:],ucm[2,:],ucm[3,:],color=[RGBA(i...,1) for i in uc],markersize=0.01,transparency=true)
Makie.scatter!(XYZ_rgb[1,:],XYZ_rgb[2,:],XYZ_rgb[3,:],color=[RGBA(i...,1) for i in uc],markersize=0.01,transparency=true)
s.center=false
s[Axis][:names,:axisnames]=("R / X","G / Y","B / Z")
record(s,joinpath(resultdir,"Unit RGB To XYZ Space.mp4"),1:360/5,framerate=12) do i
    rotate_cam!(s,5pi/180,0.0,0.0)
end

s=Makie.scatter(ucm[1,:],ucm[2,:],ucm[3,:],color=[RGBA(i...,1) for i in uc],markersize=0.01,transparency=true)
Makie.scatter!(rgb_XYZ[1,:],rgb_XYZ[2,:],rgb_XYZ[3,:],color=[RGBA(i...,1) for i in uc],markersize=0.01,transparency=true)
s.center=false
s[Axis][:names,:axisnames]=("R / X","G / Y","B / Z")
record(s,joinpath(resultdir,"Unit XYZ To RGB Space.mp4"),1:360/5,framerate=12) do i
    rotate_cam!(s,5pi/180,0.0,0.0)
end

# Transformation between LMS and Cone Contrast
bg = [0.5,0.5,0.5]
LMSToContrast,ContrastToLMS = LMSContrastMatrix(bg)
cc_lms = LMSToContrast*uqcm
lms_cc = ContrastToLMS*uqcm

s=Makie.scatter(ucm[1,:],ucm[2,:],ucm[3,:],color=[RGBA(i...,1) for i in uc],markersize=0.01,transparency=true)
Makie.scatter!(cc_lms[1,:],cc_lms[2,:],cc_lms[3,:],color=[RGBA(i...,0.1) for i in uc],markersize=0.01,transparency=true)

s=Makie.scatter(ucm[1,:],ucm[2,:],ucm[3,:],color=[RGBA(i...,0.1) for i in uc],markersize=0.01,transparency=true)
Makie.scatter!(lms_cc[1,:],lms_cc[2,:],lms_cc[3,:],color=[RGBA(i...,1) for i in uc],markersize=0.01,transparency=true)

bg = RGBToLMS*[0.5,0.5,0.5,1]
LMSToContrast,ContrastToLMS = LMSContrastMatrix(bg)
cc_lms = LMSToContrast*lms_rgb

s=Makie.scatter(lms_rgb[1,:],lms_rgb[2,:],lms_rgb[3,:],color=[RGBA(i...,1) for i in uc],markersize=0.01,transparency=true)
Makie.scatter!(cc_lms[1,:],cc_lms[2,:],cc_lms[3,:],color=[RGBA(i...,1) for i in uc],markersize=0.01,transparency=true)

# Transformation between LMS and DKL Spaces
bg = [0.5,0.5,0.5]
LMSToDKL,DKLToLMS = LMSDKLMatrix(bg,isnorm=true)
dkl_lms = LMSToDKL*uqcm
lms_dkl = DKLToLMS*uqcm

s=Makie.scatter(ucm[1,:],ucm[2,:],ucm[3,:],color=[RGBA(i...,1) for i in uc],markersize=0.01,transparency=true)
Makie.scatter!(dkl_lms[1,:],dkl_lms[2,:],dkl_lms[3,:],color=[RGBA(i...,0.1) for i in uc],markersize=0.01,transparency=true)
s.center=false
s[Axis][:names,:axisnames]=("L / L+M","M / L-M","S / S-(L+M)")
record(s,"Unit LMS To DKL$bg Space.mp4",1:360/5,framerate=12) do i
    rotate_cam!(s,5pi/180,0.0,0.0)
end

s=Makie.scatter(ucm[1,:],ucm[2,:],ucm[3,:],color=[RGBA(i...,1) for i in uc],markersize=0.01,transparency=true)
Makie.scatter!(lms_dkl[1,:],lms_dkl[2,:],lms_dkl[3,:],color=[RGBA(i...,1) for i in uc],markersize=0.01,transparency=true)
s.center=false
s[Axis][:names,:axisnames]=("L / L+M","M / L-M","S / S-(L+M)")

bg = RGBToLMS*[0.5,0.5,0.5,1]
LMSToDKL,DKLToLMS = LMSDKLMatrix(bg,isnorm=true)
dkl_lms = LMSToDKL*lms_rgb

s=Makie.scatter(lms_rgb[1,:],lms_rgb[2,:],lms_rgb[3,:],color=[RGBA(i...,1) for i in uc],markersize=0.01,transparency=true)
Makie.scatter!(dkl_lms[1,:],dkl_lms[2,:],dkl_lms[3,:],color=[RGBA(i...,1) for i in uc],markersize=0.01,transparency=true)
s.center=false
s[Axis][:names,:axisnames]=("L / L+M","M / L-M","S / S-(L+M)")
record(s,joinpath(resultdir,"LMS To DKL$bg Space.mp4"),1:360/5,framerate=12) do i
    rotate_cam!(s,5pi/180,0.0,0.0)
end

# Cone Isolating Modulation Range
lms_gray= RGBToLMS*[0.5,0.5,0.5,1]

s=Makie.surface(StatsMakie.density,lms_rgb[2,:],lms_rgb[3,:],colormap=:reds)
lines!([0,0].+lms_gray[2],[0,0].+lms_gray[3],[0,1],linewidth=3,color=:gray)
s.center=false
s[Axis][:names,:axisnames]=("M","S","Liso Strength")

s=Makie.surface(StatsMakie.density,lms_rgb[1,:],lms_rgb[3,:],colormap=:greens)
lines!([0,0].+lms_gray[1],[0,0].+lms_gray[3],[0,1],linewidth=3,color=:gray)
s.center=false
s[Axis][:names,:axisnames]=("L","S","Miso Strength")

s=Makie.surface(StatsMakie.density,lms_rgb[1,:],lms_rgb[2,:],colormap=:blues)
lines!([0,0].+lms_gray[1],[0,0].+lms_gray[2],[0,1],linewidth=3,color=:gray)
s.center=false
s[Axis][:names,:axisnames]=("L","M","Siso Strength")

# Cone Isolating RGB through a color
th = [0.5,0.5,0.5]
# each column of LMSToRGB is the cone isolating RGB direction
ConeIsoRGBVec = trimatrix(LMSToRGB)
# Scale RGB Vector into Unit Cube
ConeIsoRGBVec./=maximum(abs.(ConeIsoRGBVec),dims=1)
# since through color is the center of RGB cube, the line intersects at two symmatric points on the faces of unit cube
minc = th.-0.5*ConeIsoRGBVec;maxc = th.+0.5*ConeIsoRGBVec
ConeIsoRGB =map((i,j)->collect(i.+j), minc,[0:0.001:1].*(maxc.-minc))

Liso = hcat(ConeIsoRGB[:,1]...)'
Lisoc = [RGB(Liso[:,i]...) for i in 1:size(Liso,2)]
Lisolms=RGBToLMS*quavectors(Liso)

Miso = hcat(ConeIsoRGB[:,2]...)'
Misoc = [RGB(Miso[:,i]...) for i in 1:size(Miso,2)]
Misolms=RGBToLMS*quavectors(Miso)

Siso = hcat(ConeIsoRGB[:,3]...)'
Sisoc = [RGB(Siso[:,i]...) for i in 1:size(Siso,2)]
Sisolms=RGBToLMS*quavectors(Siso)

s=Makie.scatter(lms_rgb[1,:],lms_rgb[2,:],lms_rgb[3,:],color=[RGBA(i...,0.01) for i in uc],markersize=0.01,transparency=true)
Makie.scatter!(Lisolms[1,:],Lisolms[2,:],Lisolms[3,:],color=Lisoc,markersize=0.01,transparency=true)
Makie.scatter!(Misolms[1,:],Misolms[2,:],Misolms[3,:],color=Misoc,markersize=0.01,transparency=true)
Makie.scatter!(Sisolms[1,:],Sisolms[2,:],Sisolms[3,:],color=Sisoc,markersize=0.01,transparency=true)
s.center=false
s[Axis][:names,:axisnames]=("L","M","S")
record(s,joinpath(resultdir,"ConeIsolating_ThroughGray_LMSSpace.mp4"),1:360/5,framerate=12) do i
    rotate_cam!(s,5pi/180,0.0,0.0)
end

s=Makie.scatter(ucm[1,:],ucm[2,:],ucm[3,:],color=[RGBA(i...,0.01) for i in uc],markersize=0.01,transparency=true)
Makie.scatter!(Liso[1,:],Liso[2,:],Liso[3,:],color=Lisoc,markersize=0.01,transparency=true)
Makie.scatter!(Miso[1,:],Miso[2,:],Miso[3,:],color=Misoc,markersize=0.01,transparency=true)
Makie.scatter!(Siso[1,:],Siso[2,:],Siso[3,:],color=Sisoc,markersize=0.01,transparency=true)
s.center=false
s[Axis][:names,:axisnames]=("R","G","B")
record(s,joinpath(resultdir,"ConeIsolating_ThroughGray_RGBSpace.mp4"),1:360/5,framerate=12) do i
    rotate_cam!(s,5pi/180,0.0,0.0)
end

# DKL Isolating RGB through a background color
bg = RGBToLMS*[th;1]
LMSToDKL,DKLToLMS = LMSDKLMatrix(bg,isnorm=true)
# convert each column of DKLToLMS which are the DKL isolating LMS to RGB direction
DKLIsoRGBVec = trimatrix(LMSToRGB*DKLToLMS)
# Scale RGB Vector into Unit Cube
DKLIsoRGBVec./=maximum(abs.(DKLIsoRGBVec),dims=1)
# since through color is the center of RGB cube, the line intersects at two symmatric points on the faces of unit cube
minc = th.-0.5*DKLIsoRGBVec;maxc = th.+0.5*DKLIsoRGBVec
DKLIsoRGB =map((i,j)->collect(i.+j), minc,[0:0.001:1].*(maxc.-minc))

Lumiso = hcat(DKLIsoRGB[:,1]...)'
Lumisoc = [RGB(Lumiso[:,i]...) for i in 1:size(Lumiso,2)]
Lumisolms=RGBToLMS*quavectors(Lumiso)
Lumisodkl = LMSToDKL*Lumisolms

LMiso = hcat(DKLIsoRGB[:,2]...)'
LMisoc = [RGB(LMiso[:,i]...) for i in 1:size(LMiso,2)]
LMisolms=RGBToLMS*quavectors(LMiso)
LMisodkl = LMSToDKL*LMisolms

SLMiso = hcat(DKLIsoRGB[:,3]...)'
SLMisoc = [RGB(SLMiso[:,i]...) for i in 1:size(SLMiso,2)]
SLMisolms=RGBToLMS*quavectors(SLMiso)
SLMisodkl = LMSToDKL*SLMisolms

s=Makie.scatter(dkl_lms[1,:],dkl_lms[2,:],dkl_lms[3,:],color=[RGBA(i...,0.01) for i in uc],markersize=0.01,transparency=true)
Makie.scatter!(Lumisodkl[1,:],Lumisodkl[2,:],Lumisodkl[3,:],color=Lumisoc,markersize=0.01,transparency=true)
Makie.scatter!(LMisodkl[1,:],LMisodkl[2,:],LMisodkl[3,:],color=LMisoc,markersize=0.01,transparency=true)
Makie.scatter!(SLMisodkl[1,:],SLMisodkl[2,:],SLMisodkl[3,:],color=SLMisoc,markersize=0.01,transparency=true)
s.center=false
s[Axis][:names,:axisnames]=("L+M","L-M","S-(L+M)")
record(s,joinpath(resultdir,"DKLIsolating_ThroughGray_DKLSpace.mp4"),1:360/5,framerate=12) do i
    rotate_cam!(s,5pi/180,0.0,0.0)
end

s=Makie.scatter(lms_rgb[1,:],lms_rgb[2,:],lms_rgb[3,:],color=[RGBA(i...,0.01) for i in uc],markersize=0.01,transparency=true)
Makie.scatter!(Lumisolms[1,:],Lumisolms[2,:],Lumisolms[3,:],color=Lumisoc,markersize=0.01,transparency=true)
Makie.scatter!(LMisolms[1,:],LMisolms[2,:],LMisolms[3,:],color=LMisoc,markersize=0.01,transparency=true)
Makie.scatter!(SLMisolms[1,:],SLMisolms[2,:],SLMisolms[3,:],color=SLMisoc,markersize=0.01,transparency=true)
s.center=false
s[Axis][:names,:axisnames]=("L","M","S")
record(s,joinpath(resultdir,"DKLIsolating_ThroughGray_LMSSpace.mp4"),1:360/5,framerate=12) do i
    rotate_cam!(s,5pi/180,0.0,0.0)
end

s=Makie.scatter(ucm[1,:],ucm[2,:],ucm[3,:],color=[RGBA(i...,0.01) for i in uc],markersize=0.01,transparency=true)
Makie.scatter!(Lumiso[1,:],Lumiso[2,:],Lumiso[3,:],color=Lumisoc,markersize=0.01,transparency=true)
Makie.scatter!(LMiso[1,:],LMiso[2,:],LMiso[3,:],color=LMisoc,markersize=0.01,transparency=true)
Makie.scatter!(SLMiso[1,:],SLMiso[2,:],SLMiso[3,:],color=SLMisoc,markersize=0.01,transparency=true)
s.center=false
s[Axis][:names,:axisnames]=("R","G","B")
record(s,joinpath(resultdir,"DKLIsolating_ThroughGray_RGBSpace.mp4"),1:360/5,framerate=12) do i
    rotate_cam!(s,5pi/180,0.0,0.0)
end

# DKL Isoluminance plane
DKLToRGB = LMSToRGB*DKLToLMS
RGBToDKL = LMSToDKL*RGBToLMS
lum = 0
# rotate l-m direction around l+m axis within the Isoluminance plane
DKLIsoLumRGBVec = trivectors(hcat(map(i->DKLToRGB*RotateXYZMatrix(deg2rad(i),dims=1)*[0,1,0,0], 0:0.05:179.95)...))
# find Intersections of Isoluminance directions with faces of unit RGB cube
minmaxc = hcat(map(i->intersectlineunitorigincube((DKLToRGB*[lum,0,0,1])[1:3],DKLIsoLumRGBVec[:,i]),1:size(DKLIsoLumRGBVec,2))...)
minc=minmaxc[:,1:2:end];maxc=minmaxc[:,2:2:end]
DKLIsoLumRGB =map((i,j)->collect(i.+j), minc,[0:0.001:1].*(maxc.-minc))

IsoLum=hcat(map(i->hcat(DKLIsoLumRGB[:,i]...)',1:size(DKLIsoLumRGB,2))...)
IsoLumc = [RGB(IsoLum[:,i]...) for i in 1:size(IsoLum,2)]
IsoLumlms=RGBToLMS*quavectors(IsoLum)
IsoLumdkl = LMSToDKL*IsoLumlms

s=Makie.scatter(IsoLumdkl[2,:],IsoLumdkl[3,:],color=IsoLumc,markersize=0.01,scale_plot=false)
s[Axis][:names,:axisnames]=("L-M","S-(L+M)")
s=title(s,"DKL Isoluminance Plane through $lum")
s.center=false
Makie.save(joinpath(resultdir,"DKL Isoluminance Plane through $lum.png"),s)

# combine several Isoluminance planes
IsoLum = hcat(IsoLum,hcat(map(i->hcat(DKLIsoLumRGB[:,i]...)',1:size(DKLIsoLumRGB,2))...))

IsoLumc = [RGB(IsoLum[:,i]...) for i in 1:size(IsoLum,2)]
IsoLumlms=RGBToLMS*quavectors(IsoLum)
IsoLumdkl = LMSToDKL*IsoLumlms

s=Makie.scatter(dkl_lms[1,:],dkl_lms[2,:],dkl_lms[3,:],color=[RGBA(i...,0.01) for i in uc],markersize=0.01,transparency=true)
Makie.scatter!(IsoLumdkl[1,:],IsoLumdkl[2,:],IsoLumdkl[3,:],color=IsoLumc,markersize=0.01,transparency=true)
s.center=false
s[Axis][:names,:axisnames]=("L+M","L-M","S-(L+M)")
record(s,joinpath(resultdir,"DKLIsoluminancePlane_DKLSpace.mp4"),1:360/5,framerate=12) do i
    rotate_cam!(s,5pi/180,0.0,0.0)
end

s=Makie.scatter(lms_rgb[1,:],lms_rgb[2,:],lms_rgb[3,:],color=[RGBA(i...,0.01) for i in uc],markersize=0.01,transparency=true)
Makie.scatter!(IsoLumlms[1,:],IsoLumlms[2,:],IsoLumlms[3,:],color=IsoLumc,markersize=0.01,transparency=true)
s.center=false
s[Axis][:names,:axisnames]=("L","M","S")
record(s,joinpath(resultdir,"DKLIsoluminancePlane_LMSSpace.mp4"),1:360/5,framerate=12) do i
    rotate_cam!(s,5pi/180,0.0,0.0)
end

s=Makie.scatter(ucm[1,:],ucm[2,:],ucm[3,:],color=[RGBA(i...,0.01) for i in uc],markersize=0.01,transparency=true)
Makie.scatter!(IsoLum[1,:],IsoLum[2,:],IsoLum[3,:],color=IsoLumc,markersize=0.01,transparency=true)
s.center=false
s[Axis][:names,:axisnames]=("R","G","B")
record(s,joinpath(resultdir,"DKLIsoluminancePlane_RGBSpace.mp4"),1:360/5,framerate=12) do i
    rotate_cam!(s,5pi/180,0.0,0.0)
end

# rgb color matching chromaticity
ploc = divsum(hcat([matchlambda(w,sbrgb10) for w in sbrgb_primary]...))
wloc = divsum(hcat([matchlambda(w,sbrgb10) for w in 390:830]...))

s=Makie.lines(wloc[1,:],wloc[2,:],wloc[3,:],color=:black,linewidth=3,transparency=true)
Makie.scatter!(ploc[1,:],ploc[2,:],ploc[3,:],color=[:red,:green,:blue],markersize=0.05,transparency=true)
Makie.lines!(wloc[1,:],wloc[2,:],color=:gray40,linewidth=3,transparency=true)
Makie.scatter!(ploc[1,:],ploc[2,:],color=[:red,:green,:blue],markersize=0.05,transparency=true)

aw =  [450,500,550,600]
awloc = divsum(hcat([matchlambda(w,sbrgb10) for w in aw]...))
Makie.scatter!(awloc[1,:],awloc[2,:],color=:gray20,markersize=0.03,transparency=true)
annotations!(string.(" ",aw),Point2.(awloc[1,:],awloc[2,:]),textsize=0.08)
s[Axis][:names,:axisnames]=("r","g","b")
record(s,"Color Matching Chromaticity_rgb.mp4",1:360/5,framerate=12) do i
    rotate_cam!(s,5pi/180,0.0,0.0)
end

# xyz color matching chromaticity
d = 250
wloc = divsum(hcat([matchlambda(w,xyz10) for w in 390:0.15:830]...))
aw =  [450,500,550,600,650]
awloc = divsum(hcat([matchlambda(w,xyz10) for w in aw]...))
realxy = hcat([linepoints(wloc[1:2,1],wloc[1:2,i],d=d) for i in 2:size(wloc,2)]...)
realxyz = vcat(realxy,1 .- sum(realxy,dims=1))

xyz_rgb_primary = divsum(trimatrix(RGBToXYZ))
xyz_white = divsum(trimatrix(RGBToXYZ)*[1,1,1])
xyzToRGB = inv(xyz_rgb_primary)
xyzToRGB ./= xyzToRGB*xyz_white
rgb_realxyz = desaturate2gamut!(xyzToRGB*realxyz)
gamuttriangle = [xyz_rgb_primary xyz_rgb_primary[:,1:1]]

s = Makie.scatter(realxyz[1,:],realxyz[2,:],realxyz[3,:],color=[RGBA(rgb_realxyz[:,i]...,1) for i in 1:size(rgb_realxyz,2)],markersize=0.005,transparency=true)

s = Makie.scatter(realxy[1,:],realxy[2,:],color=[RGBA(rgb_realxyz[:,i]...,1) for i in 1:size(rgb_realxyz,2)],markersize=0.005,transparency=true,limits = FRect(0,0,0.8,0.8))
Makie.lines!(wloc[1,:],wloc[2,:],color=:black,linewidth=4)
Makie.scatter!(awloc[1,:],awloc[2,:],color=:black,markersize=0.01)
annotations!(string.(aw),[Point2(awloc[1,i]>0.3 ? awloc[1,i]+0.01 : awloc[1,i]-0.05,awloc[2,i]) for i in 1:length(aw)],textsize=0.02)
Makie.lines!(gamuttriangle[1,:],gamuttriangle[2,:],gamuttriangle[3,:],linewidth=0.5,scale_plot=false)
s=title(s,"RGB gamut in CIE xy chromaticity")
s.center=false
Makie.save(joinpath(resultdir,"RGB gamut in CIE xy chromaticity.png"),s)

# CIECAM16 Uniform Color Space
vc = cam16view(Surround=:Dark)
cam = XYZ2CAM16(100*trivectors(XYZ_rgb);vc...)
camucs = CAM16UCS(cam.J,cam.M,cam.h,form=:cartesian)

s=Makie.scatter(camucs[1,:],camucs[2,:],camucs[3,:],color=[RGBA(i...,1) for i in uc],markersize=0.01,transparency=true)
s.center=false
s[Axis][:names,:axisnames]=("J′","a′","b′")
record(s,"CAM16 Uniform Color Space.mp4",1:360/5,framerate=12) do i
    rotate_cam!(s,5pi/180,0.0,0.0)
end

s = Makie.scatter(camucs[1,:],camucs[2,:],camucs[3,:],color=[RGBA(rgb_realxyz[:,i]...,1) for i in 1:size(rgb_realxyz,2)],markersize=0.005,transparency=true)

cami = CAM16UCSinv(camucs,form=:cartesian)
t = CAM162XYZ(cami.J,cami.M,cami.h;vc...)

s=Makie.lines(StatsMakie.histogram(nbins=200),cam02.J,colormap=:reds)

lines(StatsMakie.histogram,camucs[3,:])

lines(StatsMakie.histogram,cami.h)

extrema(cam.h' .-cami.h)

t.-trivectors(100*XYZ_rgb)
sum(t.-trivectors(100*XYZ_rgb))

xyz_rgb = divsum(trivectors(XYZ_rgb))

s = Makie.scatter(realxyz[1,:],realxyz[2,:],realxyz[3,:],color=[RGBA(rgb_realxyz[:,i]...,1) for i in 1:size(rgb_realxyz,2)],markersize=0.005,transparency=true)

s.center=false
s[Axis][:names,:axisnames]=("R / L","G / M","B / S")
record(s,joinpath(resultdir,"Unit RGB To LMS Space.mp4"),1:360/5,framerate=12) do i
    rotate_cam!(s,5pi/180,0.0,0.0)
end

tIsoLumdkl = TranslateXYZMatrix(x=-0.5)*IsoLumdkl
Makie.scatter!(tIsoLumdkl[1,:],tIsoLumdkl[2,:],tIsoLumdkl[3,:],color=IsoLumc,markersize=0.01,transparency=true)
tIsoLumrgb = trivectors(DKLToRGB*tIsoLumdkl)
tIsoLumrgb = clamp.(tIsoLumrgb,0,1)
tIsoLumrgb./=maximum(abs.(tIsoLumrgb),dims=1)
Makie.scatter!(tIsoLumrgb[1,:],tIsoLumrgb[2,:],tIsoLumrgb[3,:],color=IsoLumc,markersize=0.01,transparency=true)

DKLIsoLumRGBVec = trivectors(hcat([DKLToRGB*RotateXYZMatrix(deg2rad(i>90 ? -j : j),dims=3)*RotateXYZMatrix(deg2rad(i),dims=1)*[0,1,0,0] for i=0:1:180,j=[75]][:]...))

DKLIsoLumRGBVec./=maximum(abs.(DKLIsoLumRGBVec),dims=1)
minc = th.-0.5*DKLIsoLumRGBVec;maxc = th.+0.5*DKLIsoLumRGBVec
DKLIsoLumRGB =map((i,j)->collect(i.+j), minc,[0:0.01:1].*(maxc.-minc))

IsoLum=hcat(map(i->hcat(DKLIsoLumRGB[:,i]...)',1:size(DKLIsoLumRGB,2))...)
IsoLumc = [RGB(IsoLum[:,i]...) for i in 1:size(IsoLum,2)]
IsoLumlms=RGBToLMS*quavectors(IsoLum)
IsoLumdkl = LMSToDKL*IsoLumlms

s=Makie.scatter(dkl_lms[1,:],dkl_lms[2,:],dkl_lms[3,:],color=[RGBA(i...,0.01) for i in uc],markersize=0.01,transparency=true)
Makie.scatter!(IsoLumdkl[1,:],IsoLumdkl[2,:],IsoLumdkl[3,:],color=IsoLumc,markersize=0.01,transparency=true)
s.center=false
s[Axis][:names,:axisnames]=("L+M","L-M","S-(L+M)")

s=Makie.scatter(lms_rgb[1,:],lms_rgb[2,:],lms_rgb[3,:],color=[RGBA(i...,0.01) for i in uc],markersize=0.01,transparency=true)
Makie.scatter!(IsoLumlms[1,:],IsoLumlms[2,:],IsoLumlms[3,:],color=IsoLumc,markersize=0.01,transparency=true)
s.center=false
s[Axis][:names,:axisnames]=("L","M","S")

s=Makie.scatter(ucm[1,:],ucm[2,:],ucm[3,:],color=[RGBA(i...,0.01) for i in uc],markersize=0.01,transparency=true)
Makie.scatter!(IsoLum[1,:],IsoLum[2,:],IsoLum[3,:],color=IsoLumc,markersize=0.01,transparency=true)
s.center=false
s[Axis][:names,:axisnames]=("R","G","B")

t=Makie.scatter(lms_rgb[1,:],lms_rgb[2,:],lms_rgb[3,:],color=[RGBA(i...,0.01) for i in uc],markersize=0.01,transparency=true)
Makie.scatter!(IsoLumlms[1,:],IsoLumlms[2,:],IsoLumlms[3,:],color=IsoLumc,markersize=0.01,transparency=true)

s=Makie.scatter(minc[1,:],minc[2,:],minc[3,:],color=[RGBA(minc[:,i]...,1) for i in size(minc,2)],markersize=0.01,transparency=true)
Makie.scatter!(maxc[1,:],maxc[2,:],maxc[3,:],color=[RGBA(maxc[:,i]...,1) for i in size(maxc,2)],markersize=0.01,transparency=false)

s=Makie.scatter(minc[1,:],minc[2,:],minc[3,:],color=:black,markersize=0.01,transparency=true)

Makie.scatter!(s,maxc[1,:],maxc[2,:],maxc[3,:],color=:red,markersize=0.01,transparency=true)

t=mapreduce(i->hcat(minc[:,i],maxc[:,i]),hcat,1:size(minc,2))

linesegments!(s,t[1,:],t[2,:],t[3,:])

s=Makie.scatter(maxc[1,:],maxc[2,:],maxc[3,:],color=[RGB(maxc[:,i]...) for i in size(maxc,2)],markersize=10,transparency=false)

[RGB(minc[:,i]...) for i in 1:size(minc,2)]

[RGB(maxc[:,i]...) for i in 1:size(maxc,2)]

max_stars_repo_name: Experica/ColorLab

max_issues_repo_name: Experica/ColorLab

max_forks_repo_name: Experica/ColorLab


text:
module test_errors_error

using Test

@testset "error testset" begin

@test true

g() = must_be_an_error
f() = g()
f()

end

end # module test_errors_error


text:
import Base: length, iterate, eltype
export Monic

struct Monic{T<:QuotientRing,X}
    n::Int
    Monic(::Type{P}, n) where {X,T,P<:UnivariatePolynomial{T,X}} = new{T,X}(n)
end
eltype(::Type{Z}) where Z<:Ring = Z
length(::Type{Z}) where Z<:Ring = order(Z)

function iterate(::Type{Z}) where Z<:QuotientRing
    z = zero(Z)
    (z, z)
end
function iterate(::Type{Z}, s) where Z<:ZZmod
    v = Z(s.val + 1)
    iszero(v) ? nothing : (v, v)
end
function iterate(::Type{Q}, s) where {Z<:ZZmod,P<:UnivariatePolynomial{Z},Q<:Quotient{P}}
    c = copy(s.val.coeff)
    m = length(c)
    n = deg(modulus(Q))
    for i = 1:m
        ci = next(c[i])
        if ci !== nothing
            c[i] = ci
            z = Q(c)
            return (z, z)
        else
            c[i] = zero(Z)
        end
    end
    if m < n
        resize!(c, m+1)
        c[m+1] = one(Z)
        m >= 1 && (c[m] = zero(Z))
        z = Q(c)
        return (z, z)
    end
    nothing
end

function next(c)
    cs = iterate(typeof(c), c)
    cs === nothing ? nothing : cs[1]
end

function iterate(::Type{G}) where G<:GaloisField
    G[0], 0
end

function iterate(::Type{G}, s::Integer) where G<:GaloisField
    s += 1
    s >= order(G) && return nothing
    G[s], s
end

function next(g::G) where G<:GaloisField
    s = gettypevar(G).exptable[g.val+1] + 1
    s >= order(G) ? nothing : G[s]
end

Base.length(mo::Monic{Z,X}) where {X,Z<:Ring} = length(Z)^mo.n
Base.eltype(::Monic{Z,X}) where {X,Z<:Ring} = Z[X]
function Base.iterate(mo::Monic{Z,X}) where {X,Z<:Ring}
    p0 = monom(Z[X], mo.n)
    p0, p0
end
function Base.iterate(mo::Monic{Z,X}, s) where {X,Z<:Ring}
    c = copy(s.coeff)
    n = deg(s)
    for i = 1:n
        ci = next(c[i])
        if ci !== nothing
            c[i] = ci
            z = Z[X](c)
            return (z, z)
        else
            c[i] = zero(Z)
        end
    end
    nothing
end

isqrt2(i::T) where T<:Integer = T(floor(sqrt(8*i+1) - 1)) ÷ 2
function ipair(i::Integer)
    m = isqrt2(i)
    c = m * (m + 1) ÷ 2
    m - i + c, i - c
end

# assume n1 == 0, n2 == 0 means no limitation
function index(i::T, n1::T, n2::T) where T<:Integer
    (n1 >= 0 && n2 >= 0) || throw(ArgumentError("(n1, n2) = $((n1, n2)) >= 0 required"))
    0 <= i || throw(ArgumentError("no negative index: $i"))
    if n1 > 0 && n2 > 0
        i = n1 * n2 <= 0 ? i : mod(i, n1 * n2)
    end
    m1, m2 = (n2 == 0) || (0 < n1 <= n2) ? (n1, n2) : (n2, n1)
    m = (m1 + 1) * m1 ÷ 2
    if m1 == 0 || i < m
        ipair(i)
    else
        j = i - m
        if m2 == 0 || j < (m2 - m1) * m1
            j1, j2 = fldmod(j, m1)
            s = n1 == m1 ? n1 - j2 - 1 : j1 + m1 - j2
            t = j1 + m1 - s
            s, t
        else
            s, t = index(n1 * n2 - 1 - i, n1, n2)
            n1 - s - 1, n2 - t - 1
        end
    end
end

function indexv(i::T, nn::Vector{T}) where T<:Integer
    m = length(nn)
    m == 0 && return T[]
    n1 = nn[1]
    m == 1 && return [n1 == 0 ? i : mod(i, n1)]
    m == 2 && return [index(i, n1, nn[2])...]
    m2 = m >> 1
    n1 = prod(nn[1:m2])
    n2 = prod(nn[m2+1:m])
    s, t = index(i, n1, n2)
    vcat(indexv(s, nn[1:m2]), indexv(t, nn[m2+1:m]))
end

len(::Type, d...) = 0
len(T::Type{<:ZZmod}, d...) = modulus(T)
function len(T::Type{<:FractionRing{S}}, d...) where S
    n = len(S, d...)
    n == 0 ? 0 : (n-1)^2 + 1
end
len(T::Type{UnivariatePolynomial{S}}, d::Integer) where S = intpower(len(S), d)
len(T::Type{<:QuotientRing{S}}) where S<:UnivariatePolynomial = len(S, deg(modulus(T)-1))
len(T::Type{<:GaloisField}) = order(T)

ofindex(a::Integer, T::Type{<:Unsigned}) = T(a)
ofindex(a::Integer, T::Type{<:Signed}) = iseven(a) ? -(T(a) >> 1) : T(a+1) >> 1
ofindex(a::Integer, T::Type{ZZ{S}}) where S = T(ofindex(a, S))
ofindex(a::Integer, T::Type{<:ZZmod{m,S}}) where {m,S} = T(ofindex(a, unsigned(S)))
function ofindex(a::Integer, T::Type{<:QuotientRing{S}}) where {S<:UnivariatePolynomial}
    d = deg(modulus(T))
    T(ofindex(a, S, d) - monom(S,d))
end

function ofindex(a::Integer, T::Type{<:FractionRing{S}}) where S
    a == 0 && return zero(T)
    s, t = index(a - 1, len(T), len(T))
    T(ofindex(t+1, S), S(ofindex(s + 1, unsigned(S))))
end
function ofindex(a::Integer, P::Type{<:UnivariatePolynomial{S}}, d::Integer) where S
    P([ofindex.(indexv(a, fill(oftype(a, len(S)), d)), S); 1])
end
function ofindex(a::Integer, G::Type{<:GaloisField})
    G[mod(a, order(G))]
end

struct Factors{T<:Integer,P}
   f::P
   Factors(x::V) where {T,V<:AbstractVector{<:Pair{T}}} = new{T,V}(x)
end
Base.length(fi::Factors) = isempty(fi.f) ? 1 : prod(x+1 for x in _values(fi.f))

_values(d::AbstractDict) = values(d)
_values(d::Vector{<:Pair}) = (last(x) for x in d)

Base.iterate(fi::Factors) = (1, zeros(Int,length(fi.f)))
function Base.iterate(fi::Factors, s::Array{Int})
    f = fi.f
    n = length(f)
    n == 0 && return nothing
    for i = 1:n
        si = s[i]
        if si < last(f[i])
            s[i] = si + 1
            break
        else
            i == n && return nothing
            s[i] = 0
        end
    end
    p = prod(first(f[i])^s[i] for i = 1:n)
    p, s
end
"""
    factors(n)
    factors(Primes.factor(n))

Return an iterator generating all factors of integer `n`.
"""
factors(f::Primes.Factorization) = Factors(f.pe)
factors(f::AbstractVector) = Factors(f)
factors(n::Integer) = factors(factor(n).pe)


text:
###
# Courtesy of @shashi

struct SortMergeItr{S, T, F<:Function}
    a::S
    b::T
    lt::F
end

SortMergeItr(a, b) = SortMergeItr(a, b, isless)

import Base: IsInfinite, SizeUnknown, HasLength, HasShape,
             IteratorSize, IteratorEltype, HasEltype, EltypeUnknown

function IteratorSize(::Type{SortMergeItr{S,T, F}}) where {S,T, F}
    iS = IteratorSize(S)
    iT = IteratorSize(T)
    if all(isa.((iS, iT), Union{HasLength, HasShape}))
        return HasLength()
    elseif IsInfinite() in (iS, iT)
        return IsInfinite()
    elseif SizeUnknown() in (iS, iT)
        return SizeUnknown()
    end
end
Base.length(itr::SortMergeItr) = length(itr.a) + length(itr.b)
Base.eltype(::Type{SortMergeItr{S,T, F}}) where {S,T, F} = Union{eltype(S), eltype(T)}
function IteratorEltype(::Type{SortMergeItr{S,T, F}}) where {S,T,F}
    if IteratorEltype(S) == HasEltype() && IteratorEltype(T) == HasEltype()
        HasEltype()
    else
        EltypeUnknown()
    end
end

function _pick(a, b, lt)
    if a === nothing && b === nothing
        return nothing
    elseif a === nothing
        return b[1], (2, a, b)
    elseif b === nothing
        return a[1], (1, a, b)
    else
        lt(a[1], b[1]) ? (a[1], (1, a, b)) : (b[1], (2, a, b))
    end
end

function Base.iterate(s::SortMergeItr)
    a = iterate(s.a)
    b = iterate(s.b)
    _pick(a, b, s.lt)
end

function Base.iterate(s::SortMergeItr, (select, a, b,))
    if select == 1
        a = iterate(s.a, a[2])
    else
        b = iterate(s.b, b[2])
    end

    _pick(a, b, s.lt)
end

###

function hassmallerkey(a, b)
    a_key, _ = a
    b_key, _ = b
    isless(a_key, b_key)
end

function simplexiterator(cplx)
    vals = map(t -> t.values, cplx)
    itrs = map(finduniquesorted, vals)
    aug_itrs = map(enumerate(itrs)) do (n, itr)
        (i[1] => (n, i[2]) for i in itr)
    end

    foldl((a, b) -> SortMergeItr(b, a, hassmallerkey), aug_itrs)
end


text:
abstract type AbstractMetric end

struct UnitEuclideanMetric <: AbstractMetric
    dim::Int
end

# Create a `UnitEuclideanMetric`; required for an unified interface
(ue::UnitEuclideanMetric)(::Nothing) = UnitEuclideanMetric(ue.dim)

function _string_diag(d, n_chars::Int=32) :: String
    s_diag = string(d)
    l = length(s_diag)
    s_dots = " ..."
    n_diag_chars = n_chars - length(s_dots)
    return s_diag[1:min(n_diag_chars,end)] * (l > n_diag_chars ? s_dots : "")
end

Base.show(io::IO, uem::UnitEuclideanMetric) = print(io, _string_diag(ones(uem.dim)))

struct DiagEuclideanMetric{A<:AbstractVector} <: AbstractMetric
    # Diagnal of the inverse of the mass matrix
    M⁻¹     ::  A
    # Sqare root of the inverse of the mass matrix
    sqrtM⁻¹ ::  A
    # Pre-allocation for intermediate variables
    _temp   ::  A
end

function DiagEuclideanMetric(M⁻¹::AbstractVector{T}) where {T<:Real}
    return DiagEuclideanMetric(M⁻¹, sqrt.(M⁻¹), Vector{T}(undef, size(M⁻¹, 1)))
end
DiagEuclideanMetric(D::Int) = DiagEuclideanMetric(ones(Float64, D))

# Create a `DiagEuclideanMetric` with a new `M⁻¹`
(dem::DiagEuclideanMetric)(M⁻¹::AbstractVector{<:Real}) = DiagEuclideanMetric(M⁻¹)

Base.show(io::IO, dem::DiagEuclideanMetric) = print(io, _string_diag(dem.M⁻¹))

function Base.getproperty(dem::DiagEuclideanMetric, d::Symbol)
    return d === :dim ? size(getfield(dem, :M⁻¹), 1) : getfield(dem, d)
end

struct DenseEuclideanMetric{
    AV<:AbstractVector,
    AM<:AbstractMatrix,
    TcholM⁻¹<:UpperTriangular,
} <: AbstractMetric
    # Inverse of the mass matrix
    M⁻¹::AM
    # U of the Cholesky decomposition of the mass matrix
    cholM⁻¹::TcholM⁻¹
    # Pre-allocation for intermediate variables
    _temp::AV
end

function DenseEuclideanMetric(M⁻¹::AbstractMatrix{T}) where {T<:Real}
    _temp = Vector{T}(undef, size(M⁻¹, 1))
    return DenseEuclideanMetric(M⁻¹, cholesky(Symmetric(M⁻¹)).U, _temp)
end
DenseEuclideanMetric(D::Int) = DenseEuclideanMetric(Matrix{Float64}(I, D, D))

# Create a `DenseEuclideanMetric` with a new `M⁻¹`
(dem::DenseEuclideanMetric)(M⁻¹::AbstractMatrix{<:Real}) = DenseEuclideanMetric(M⁻¹)

Base.show(io::IO, dem::DenseEuclideanMetric) = print(io, _string_diag(diag(dem.M⁻¹)))

function Base.getproperty(dem::DenseEuclideanMetric, d::Symbol)
    return d === :dim ? size(getfield(dem, :M⁻¹), 1) : getfield(dem, d)
end


text:
@testset "Mutability" begin
    @testset "setindex!" begin
        s = dna"ACGT"
        s[1] = DNA_A
        @test s == dna"ACGT"
        s[1] = DNA_C
        @test s == dna"CCGT"
        s[2] = DNA_T
        @test s == dna"CTGT"
        s[4] = DNA_A
        @test s == dna"CTGA"
        @test_throws BoundsError s[0]
        @test_throws BoundsError s[5]

        s = dna"ACGTACGT"
        s[3:5] = DNA_A
        @test s == dna"ACAAACGT"
        s[3:5] = DNA_T
        @test s == dna"ACTTTCGT"
        s[1:8] = DNA_C
        @test s == dna"CCCCCCCC"
        s[:] = DNA_G
        @test s == dna"GGGGGGGG"
        @test_throws BoundsError s[0:3]
        @test_throws BoundsError s[5:10]

        s = dna"ACGTACGT"
        s[[1,2,3]] = DNA_T
        @test s == dna"TTTTACGT"
        s[[1,5,8]] = DNA_C
        @test s == dna"CTTTCCGC"
        @test_throws BoundsError s[[3,9]] = DNA_A

        s = dna"ACGT"
        s[[true, false, false, true]] = DNA_G
        @test s == dna"GCGG"
        s[trues(4)] = DNA_A
        @test s == dna"AAAA"
        @test_throws BoundsError s[[true, false, false]] = DNA_G
        @test_throws BoundsError s[[true, false, false, false, true]] = DNA_G

        s = dna"ACGTACGT"
        s[2:3] = dna"AA"
        @test s == dna"AAATACGT"
        s[7:8] = dna"CC"
        @test s == dna"AAATACCC"
        s[:] = dna"AACCGGTT"
        @test s == dna"AACCGGTT"
        @test_throws BoundsError       s[0:1] = dna"AA"
        @test_throws DimensionMismatch s[3:4] = dna"A"

        s = dna"ACGTACGT"
        s[[1,4]] = dna"TA"
        @test s == dna"TCGAACGT"
        s[[2,3,5]] = dna"CAT"
        @test s == dna"TCAATCGT"
        @test_throws BoundsError       s[[1,2,9]] = dna"AAA"
        @test_throws DimensionMismatch s[[1,2,8]] = dna"AA"

        s = dna"ACGT"
        s[[true,false,true,false]] = dna"TT"
        @test s == dna"TCTT"
        s[trues(4)] = dna"AAAA"
        @test s == dna"AAAA"
        @test_throws BoundsError       s[[true,false,true]] = dna"TT"
        @test_throws DimensionMismatch s[[true,false,true,true]] = dna"TT"
    end

    @testset "resize!" begin
        seq = dna""
        resize!(seq, 100)
        @test length(seq) == 100
        resize!(seq, 200)
        @test length(seq) == 200
        seq1 = seq[3:198]
        resize!(seq1, 55)
        @test length(seq1) == 55
        resize!(seq,  10)
        @test length(seq) == 10
        @test_throws ArgumentError resize!(seq, -1)
    end

    @testset "empty!" begin
        seq = dna"ACG"
        @test empty!(seq) == dna""
        @test length(seq) == 0
    end

    @testset "push!" begin
        seq = dna""
        @test push!(seq, DNA_A) == dna"A"
        @test push!(seq, DNA_C) == dna"AC"
        @test seq == dna"AC"
    end

    @testset "pushfirst!" begin
        seq = dna""
        @test pushfirst!(seq, DNA_A) == dna"A"
        @test pushfirst!(seq, DNA_C) == dna"CA"
        @test seq == dna"CA"
    end

    @testset "pop!" begin
        seq = dna"ACGT"
        @test pop!(seq) === DNA_T
        @test seq == dna"ACG"
        @test pop!(seq) === DNA_G
        @test seq == dna"AC"
        @test_throws ArgumentError pop!(dna"")
    end

    @testset "popfirst!" begin
        seq = dna"ACGT"
        @test popfirst!(seq) === DNA_A
        @test seq == dna"CGT"
        @test popfirst!(seq) === DNA_C
        @test seq == dna"GT"
        @test_throws ArgumentError popfirst!(dna"")
    end

    @testset "insert!" begin
        seq = dna"ACGT"
        @test insert!(seq, 2, DNA_G) == dna"AGCGT"
        @test insert!(seq, 5, DNA_A) == dna"AGCGAT"
        @test_throws BoundsError insert!(seq, 10, DNA_T)
    end

    @testset "deleteat!" begin
        seq = dna"ACGT"
        @test deleteat!(seq, 1) == dna"CGT"
        @test deleteat!(seq, 2) == dna"CT"
        @test_throws BoundsError deleteat!(seq, 10)

        seq = dna"ACGTACGT"
        @test deleteat!(seq, 3:5) == dna"ACCGT"
        @test_throws BoundsError deleteat!(seq, 10:12)
    end

    @testset "append!" begin
        seq = dna""
        @test append!(seq, dna"A") == dna"A"
        @test append!(seq, dna"ACG") == dna"AACG"
    end

    @testset "copyto!" begin
        seq = dna"GGG"
        @test copyto!(seq, dna"ACG") == dna"ACG"
        @test copyto!(seq, dna"TTA") == dna"TTA"

        seq = dna"TCCC"
        @test copyto!(seq, 2, dna"TT", 1, 2) == dna"TTTC"
        seq = dna"TCCC"
        @test copyto!(seq, 2, dna"TT", 1, 1) == dna"TTCC"

        seq = dna"ACGT"
        @test copyto!(seq, seq) == dna"ACGT"
        @test copyto!(seq, 1, seq, 3, 2) == dna"GTGT"
        seq = dna"ACGT"
        @test copyto!(seq, 3, seq, 1, 2) == dna"ACAC"
    end
end


text:
include("../src/Interfaces.jl")
using Interfaces
using Base.Test

type Foo
    x::Int
end

type Bar
    x::ASCIIString
end

@interface MyInterface begin
    func(self::MyInterface, x::Int)
end
common_func(obj::MyInterface) = func(obj, 2)

@interface OtherInterface begin
    run(self::OtherInterface, x::Int)
end

@test_throws(UndefVarError, methods_exist(MyInterface, Foo, current_module()))
func(self::Foo) = self.x
@test_throws(ErrorException, methods_exist(MyInterface, Foo, current_module()))
func(self::Foo, x::Int) = self.x + x

implements(Foo, MyInterface)
@test(common_func(Foo(4)) == 6)

func(self::Bar, x::Int) = string(self.x, ",", x)
run(self::Bar, x::Int) = println("Running...")
@implements Bar <: (MyInterface, OtherInterface)


text:
module RocketZipObservableTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "ZipObservable" begin

    println("Testing: zipped")

    @testset begin
        @test_throws ErrorException zipped()
    end

    @testset begin
        source = zipped(of(1), of(0.0))
        io = IOBuffer()

        show(io, source)

        printed = String(take!(io))

        @test occursin("ZipObservable", printed)
        @test occursin(string(eltype(source)), printed)
        @test eltype(source) === Tuple{Int, Float64}
    end

    run_testset([
        (
            source      = zipped(of(1), of(0.0)),
            values      = @ts([ (1, 0.0), c ]),
            source_type = Tuple{Int, Float64}
        ),
        (
            source      = zipped(completed(Int), of(0.0)),
            values      = @ts(c),
            source_type = Tuple{Int, Float64}
        ),
        (
            source      = zipped(of(1), completed(Int)),
            values      = @ts(c),
            source_type = Tuple{Int, Int}
        ),
        (
            source      = zipped(completed(Int), completed(Int)),
            values      = @ts(c),
            source_type = Tuple{Int, Int}
        ),
        (
            source      = zipped(faulted("e1"), of(0.0)),
            values      = @ts(e("e1"))
        ),
        (
            source      = zipped(of(1), faulted("e2")),
            values      = @ts(e("e2"))
        ),
        (
            source      = zipped(faulted("e1"), faulted("e2")),
            values      = @ts(e("e1"))
        ),
        (
            source      = zipped(never(), of(0.0)),
            values      = @ts()
        ),
        (
            source      = zipped(of(1), never()),
            values      = @ts()
        ),
        (
            source      = zipped(never(), never()),
            values      = @ts()
        ),
        (
            source      = zipped(from(1:2), from(1:5)),
            values      = @ts([ (1, 1), (2, 2), c ]),
            source_type = Tuple{Int, Int}
        ),
        (
            source      = zipped(from(1:2) |> async(0), from(1:5)),
            values      = @ts([ (1, 1) ] ~ [ (2, 2) ] ~ c),
            source_type = Tuple{Int, Int}
        ),
        (
            source      = zipped(from(1:2), from(1:5) |> async(0)),
            values      = @ts([ (1, 1) ] ~ [ (2, 2), c ]),
            source_type = Tuple{Int, Int}
        ),
        (
            source      = zipped(from(1:2) |> async(0), from(1:5) |> async(0)),
            values      = @ts([ (1, 1) ] ~ [ (2, 2), c ]),
            source_type = Tuple{Int, Int}
        ),
        (
            source      = zipped(from(1:2), from(1:5), of('c')),
            values      = @ts([ (1, 1, 'c'), c ]),
            source_type = Tuple{Int, Int, Char}
        ),
        (
            source      = zipped(from(1:2), from(1:2)),
            values      = @ts([ (1, 1), (2, 2), c ]),
            source_type = Tuple{Int, Int}
        ),
        (
            source      = zipped(from(1:3), from(1:3) |> ignore(1), from(1:3) |> ignore(2)),
            values      = @ts([ (1, 2, 3), c ]),
            source_type = Tuple{Int, Int, Int}
        ),
    ])

    @testset begin
        s1 = Subject(Int)
        s2 = Subject(Float64)
        s3 = Subject(String)

        latest = zipped(s1, s2, s3)
        values = Vector{Any}()
        subscription = subscribe!(latest, lambda(
            on_next     = (d) -> push!(values, d),
            on_error    = (e) -> push!(values, e),
            on_complete = ()  -> push!(values, "completed")
        ))

        @test values == [ ]

        next!(s1, 1)

        @test values == [ ]

        next!(s2, 2.0)

        @test values == [ ]

        next!(s3, "Hello")

        @test values == [ (1, 2.0, "Hello") ]

        next!(s2, 3.0)

        @test values == [ (1, 2.0, "Hello") ]

        next!(s3, "Hello, world!")

        @test values == [ (1, 2.0, "Hello") ]

        next!(s1, 2)

        @test values == [ (1, 2.0, "Hello"), (2, 3.0, "Hello, world!") ]

        next!(s1, 3)

        @test values == [ (1, 2.0, "Hello"), (2, 3.0, "Hello, world!") ]

        next!(s1, 4)

        @test values == [ (1, 2.0, "Hello"), (2, 3.0, "Hello, world!") ]

        next!(s2, 4.0)

        @test values == [ (1, 2.0, "Hello"), (2, 3.0, "Hello, world!") ]

        next!(s3, "")

        @test values == [ (1, 2.0, "Hello"), (2, 3.0, "Hello, world!"), (3, 4.0, "") ]

        complete!(s1);

        @test values == [ (1, 2.0, "Hello"), (2, 3.0, "Hello, world!"), (3, 4.0, "")]

        next!(s2, 5.0)

        @test values == [ (1, 2.0, "Hello"), (2, 3.0, "Hello, world!"), (3, 4.0, "") ]

        complete!(s2)

        @test values == [ (1, 2.0, "Hello"), (2, 3.0, "Hello, world!"), (3, 4.0, "") ]

        next!(s3, "last")

        @test values == [ (1, 2.0, "Hello"), (2, 3.0, "Hello, world!"), (3, 4.0, ""), (4, 5.0, "last"), "completed" ]

        unsubscribe!(subscription)
    end

end

end


text:
module BinaryTraits

using MacroTools: rmlines, postwalk
using ExprTools: splitdef, combinedef

export BinaryTrait
export @trait, @assign, @implement, @check, @traitfn
export check, is_trait, traits, required_contracts
export init_traits

include("types.jl")
include("misc.jl")
include("utils.jl")
include("trait.jl")
include("assignment.jl")
include("interface.jl")
include("traitfn.jl")

end # module


text:
# Replacing `isfile(path...)` could be an issue due to its use of Varargs

@testset "isfile" begin
    null_file = string(@__FILE__, ".null")  # Note: tempfile() on Windows creates a file
    @test isfile(null_file) == false

    # Can't use `null_file` since closures don't currently work.
    # Note: @__FILE__ is resolved in the context of this file.
    patch = @patch isfile(p...) = first(p) == string(@__FILE__, ".null")
    apply(patch) do
        @test (@mock isfile(null_file)) == true
    end
end


text:
export isdecomposable, issmooth, isdeterministic

function check_decomposable(ns::AbstractVector{<:AbstractNode})
    d = true
    for i in eachindex(ns)
        l = ns[i]
        scope_l = scope(l)
        for j in (i+1):length(ns)
            r = ns[j]
            d &= isdisjoint(scope_l, scope(r))
        end
    end
    return d
end

function check_smooth(ns::AbstractVector{<:AbstractNode})
    s = true
    n1 = first(ns)
    for ni in ns[2:end]
        s &= scope(ni) == scope(n1)
    end
    return s
end

function check_determinism(ns::AbstractVector{<:AbstractNode})
    d = true
    for i in eachindex(ns)
        l = ns[i]
        sl = support(l)
        for j in (i+1):length(ns)
            r = ns[j]
            if scope(l) == scope(r)
                d &= isdisjointsupport(sl, support(r))
            else
                d &= false
            end
        end
    end
    return d
end

"""
    isdecomposable(n)

Return true if the node is decomposable and false otherwise.
"""
isdecomposable(n::Node{V,S,T,P,ProductNode}) where {V,S,T,P<:Tuple{Decomposability}} = true
isdecomposable(n::AbstractNode) = false

"""
    issmooth(n)

Return true if the node is smooth and false otherwise.
"""
issmooth(n::Node{V,S,T,P,SumNode}) where {V,S,T,P<:Tuple{Smoothness}} = true
issmooth(n::Node{V,S,T,P,SumNode}) where {V,S,T,P<:Tuple{Smoothness,Determinism}} = true
issmooth(n::AbstractNode) = false

"""
    isdeterministic(n)

Return true if the node is deterministic and false otherwise.
"""
isdeterministic(n::Node{V,S,T,P,SumNode}) where {V,S,T,P<:Tuple{Determinism}} = true
isdeterministic(n::Node{V,S,T,P,SumNode}) where {V,S,T,P<:Tuple{Smoothness,Determinism}} = true
isdeterministic(n::AbstractNode) = false


text:
using BenchmarkTools
include("forest_fire.jl")

function counter(model::ABM)
  on_fire = 0
  green = 0
  burned = 0
  for tree in values(model.agents)
    if tree.status == 1
      green += 1
    elseif tree.status == 2
      on_fire += 1
    else
      burned += 1
    end
  end
  return green, on_fire, burned
end

height=100
d=0.6
nsteps = 100
when = 1:nsteps

aggregators = Dict(:model => [counter])
size_range = 100:100:1000
agent_properties = [:model]
results = Float64[]
for width in size_range
  b = @benchmarkable data=step!(forest, tree_step!, nsteps,
  aggregators, when=when) setup=(forest=model_initiation(d=d,
  griddims=($width, height), seed=2))

  j = run(b)
  push!(results, minimum(j.times)/1e9)  # convert them to seconds
end

results

#=
0.0623855
0.1409112
0.2434189
0.3797918
0.8206902
0.7178797
0.9326405
1.364939
1.5469985
2.051536
=#


text:
# # Tutorials for Trixi.jl

# The tutorial section for [Trixi.jl](https://github.com/trixi-framework/Trixi.jl) also contains
# interactive step-by-step explanations via [Binder](https://mybinder.org).

# Right now, you are using the classic documentation. The corresponding interactive notebooks can
# be opened in [Binder](https://mybinder.org/) and viewed in [nbviewer](https://nbviewer.jupyter.org/)
# via the icons ![](https://mybinder.org/badge_logo.svg) and ![](https://raw.githubusercontent.com/jupyter/design/master/logos/Badges/nbviewer_badge.svg)
# in the respective tutorial.
# You can download the raw notebooks from GitHub via ![](https://camo.githubusercontent.com/aea75103f6d9f690a19cb0e17c06f984ab0f472d9e6fe4eadaa0cc438ba88ada/68747470733a2f2f696d672e736869656c64732e696f2f62616467652f646f776e6c6f61642d6e6f7465626f6f6b2d627269676874677265656e).
# **Note:** To improve responsiveness via caching, the notebooks are updated only once a week. They are only
# available for the latest stable release of Trixi at the time of caching.
# There are tutorials for the following topics:

# ### 1 Adding a new equation
# #### [1.1 Scalar conservation law](@ref cubic_conservation_law)
#-
# This tutorial explains how to add a new physics model using the example of the cubic conservation
# law. First, we define the equation using a `struct` `CubicEquation` and the physical flux. Then,
# the corresponding standard setup in Trixi.jl (`mesh`, `solver`, `semi` and `ode`) is implemented
# and the ODE problem is solved by OrdinaryDiffEq's `solve` method.

# #### [1.2 Nonconservative advection](@ref nonconservative_advection)
#-
# In this part, another physics model is implemented, the nonconservative linear advection equation.
# We run two different simulations with different levels of refinement and compare the resulting errors.

# ### [2 Differentiable programming](@ref differentiable_programming)
#-
# This part deals with some basic differentiable programming topics. For example, a Jacobian, its
# eigenvalues and a curve of total energy (through the simulation) are calculated and plotted for
# a few semidiscretizations. Moreover, we calculate an example for propagating errors with Measurement.jl
# at the end.

# ### [3 Unstructured meshes with HOHQMesh.jl](@ref hohqmesh_tutorial)
#-
# The purpose of this tutorial is to demonstrate how to use the [`UnstructuredMesh2D`](@ref)
# functionality of Trixi.jl. This begins by running and visualizing an available unstructured
# quadrilateral mesh example. Then, the tutorial will demonstrate how to conceptualize a problem
# with curved boundaries, generate a curvilinear mesh using the available [HOHQMesh](https://github.com/trixi-framework/HOHQMesh)
# software in the Trixi.jl ecosystem, and then run a simulation using Trixi.jl on said mesh.


text:
path=pwd()

##
import Pkg; Pkg.add("Documenter")
using Documenter
Pkg.activate(".")
using WavesAndEigenvalues
##
cd("./docs")
makedocs(format = Documenter.HTML(
        prettyurls = false,
        ),
        sitename="WavesAndEigenvalues.jl",
        authors = "Georg A. Mensah, Alessandro Orchini, and contributors.",
        modules =[WavesAndEigenvalues],
        pages = [
        "Home" => "index.md",
        "Install" => "installation.md",
        "Examples" => Any["tutorial_00_NLEVP.md",
                        "tutorial_01_rijke_tube.md",
                        "tutorial_02_global_eigenvalue_solver.md",
                        "tutorial_03_local_eigenvalue_solver.md",
                        "tutorial_04_perturbation_theory.md",
                        "tutorial_05_mesh_refinement.md",
                        "tutorial_06_second_order_elements.md",
                        "tutorial_07_Bloch_periodicity.md",
                        "tutorial_08_custom_FTF.md",
                        "tutorial_09_forcing.md",
                        ],
        "API" => Any["NLEVP.md","Meshutils.md", "Helmholtz.md"]
        ]
        )
cd(path)

deploydocs(
    repo = "github.com/JulHoltzDevelopers/WavesAndEigenvalues.jl.git",
    push_preview=true
    )


text:
# EP simulations with i.i.d. patterns and noiseless labels

using EPDilutedClassifier

N=128; ρ=0.25; α=6.0;
epsconv=1e-4; maxiter=50000; η=1.0; damp=0.9995;
w=generate_weights(N,ρ);
report=Reporter(w);

run_instance_gaussian_iid(w,ρ,α,η,callback=report,ρ0=ρ,
   λ0=1.0,δρ=0.0,η0=η,δη=0.0,damp=damp,epsconv=epsconv,maxiter=maxiter)

max_stars_repo_name: mpieropan/EPDilutedClassifier

max_issues_repo_name: mpieropan/EPDilutedClassifier

max_forks_repo_name: mpieropan/EPDilutedClassifier


text:
"""
    lagFn(x,n=1)

Create a matrix or vector of lagged values.

# Input
- `x::Array`: T Vector or Txk matrix
- `n::Int`:   scalar, order of lag. For instance, 2 puts x[t-2,:] on row t

# Output
- `z::Array`:  Txk matrix of lags

"""
function lagFn(x,n=1)

  (T,k) = (size(x,1),size(x,2))

  if n < 0                                    #leads
    z = [ x[1-n:T,:]; fill(NaN,(abs(n),k)) ]
  elseif n == 0                               #no lag or lead
    z = copy(x)
  elseif n > 0                                #lags
    z = [ fill(NaN,(n,k)); x[1:T-n,:] ]
  end

  isa(x,AbstractVector) && (z = z[:,1])               #z should be vector if x is

  return z

end

max_stars_repo_name: duzaichuan/FinancialEconometrics

max_issues_repo_name: duzaichuan/FinancialEconometrics

max_forks_repo_name: duzaichuan/FinancialEconometrics


text:
immutable QueryableSelect{T,Provider} <: Queryable{T,Provider}
    source
    f::Expr
end

function select{TS,Provider}(source::Queryable{TS,Provider}, f::Function, f_expr::Expr)
    T = Base.return_types(f, (TS,))[1]
    return QueryableSelect{T,Provider}(source, f_expr)
end


text:
module AutocorrelationTime

import StatsBase:autocor
export integrated_autocorrelation_time

"""
Find the integrated autocorrelation time of a vector of measurements

Based on the tips given in http://www.stat.unc.edu/faculty/cji/Sokal.pdf
"""
function integrated_autocorrelation_time(samples::Vector)
  c = 6.0 #Some constant determining the tradeoff between large bias (small c) and large variance (large c). Okay value should be around c ∈ [4,10]
  Δtstart = 1
  Δtstop = 30
  ac = autocor(samples, Δtstart:Δtstop)

  τint = 1/2
  M = Δtstart #Window width (to be increased until minimum width is reached)

  for i=1:30
    if c*τint < M
      break
    end
    τint += ac[M - Δtstart + 1]
    M+=1

    if M>Δtstop #More autocorrelation calculation needed
      Δtstart = Δtstop + 1
      Δtstop = 3Δtstop - Δtstart + 1
      ac = autocor(samples, Δtstart:Δtstop)
    end
  end

  return τint
end

end # module


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.4.2
#   kernelspec:
#     display_name: Julia 1.4.0
#     language: julia
#     name: julia-1.4
# ---

# # 情報幾何学のお勉強のぉと

# ## 導入
#
# $p$,$q$ を正の実数とする. $y=px^2$, $x=qy^2$ という $xy$ - 平面上の曲線を考える．この二つの曲線は原点ともう一つの点 $P=(x_0,y_0)$ と交わる．簡単な計算によってそれは $p$, $q$ を用いると $P = (x_0,y_0) = (p^{-2/3}q^{-1/3}, p^{-1/3}q^{-2/3})$ で与えられることがわかる．これをプログラミング言語 Julia を用いて確認する．

using Plots

using Interact

function vis_crosspoint(p::Float64,q::Float64)
    x = range(0,3,length=30) |> collect
    y = range(0,3,length=30) |> collect
    plot!(x,p*x.^2,color=:blue)
    plot!(q*y.^2,y,color=:orange)
    scatter!([p^(-2/3)*q^(-1/3)],[p^(-1/3)*q^(-2/3)])
end

# ## 試しに描画

plot(xlim=(0.,3.), ylim=(0.,3.), legend=:false)
vis_crosspoint(1.,2.4)

# ## Interact.jl で $p$ と $q$ を変動させた様子を描画

# +
p_slider = slider(0.001:0.1:3,label="p")
q_slider = slider(0.001:0.1:3,label="q")

plot(xlim=(0.,3.), ylim=(0.,3.), legend=:false)
cross_pt = map(
    (p,q)->vis_crosspoint(p,q),
    map(observe,[p_slider,q_slider])...,
)

[p_slider,q_slider,cross_pt] .|> display
# -
# ## $q$ を固定して $p$ だけ動かす

# +
plot(xlim=(0.,3.), ylim=(0.,3.), legend=:false)
anim = @animate for p in range(0.1,1,length=10)
    p=vis_crosspoint(p,1.4)
    plot!(p)
end

gif(anim, "/tmp/anim_fps15.gif", fps = 3)