text:
module Test_klf

using Random
using LinearAlgebra
using MatrixPencils
using Test

Random.seed!(2351);

@testset "Matrix Pencils Utilities" begin

@testset "klf_rlsplit" begin

fast = true
for fast in (false,true)

M2 = zeros(0,0); N2 = zeros(0,0);
M = copy(M2); N = copy(N2);

@time M1, N1, Q, Z, ν, μ, n, m, p  = klf_rlsplit(M, N, fast = fast)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      ν == [] && μ == [] && n == 0 && m == 0  && p == 0

M2 = zeros(3,0); N2 = zeros(3,0);
M = copy(M2); N = copy(N2);

@time M1, N1, Q, Z, ν, μ, n, m, p  = klf_rlsplit(M, N, fast = fast, finite_infinite = false)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      ν == [] && μ == [] && n == 0 && m == 0  && p == 3

M2 = zeros(3,0); N2 = zeros(3,0);
M = copy(M2); N = copy(N2);

@time M1, N1, Q, Z, ν, μ, n, m, p  = klf_rlsplit(M, N, fast = fast, finite_infinite = true)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      ν == [3] && μ == [0] && n == 0 && m == 0  && p == 0

M2 = zeros(0,3); N2 = zeros(0,3);
M = copy(M2); N = copy(N2);

@time M1, N1, Q, Z, ν, μ, n, m, p  = klf_rlsplit(M, N, fast = fast, finite_infinite = false)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      ν == [0] && μ == [3] && n == 0 && m == 0  && p == 0

M2 = zeros(0,3); N2 = zeros(0,3);
M = copy(M2); N = copy(N2);

@time M1, N1, Q, Z, ν, μ, n, m, p  = klf_rlsplit(M, N, fast = fast, finite_infinite = true)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      ν == [] && μ == [] && n == 0 && m == 3  && p == 0

M2 = ones(1,1); N2 = zeros(1,1);
M = copy(M2); N = copy(N2);

@time M1, N1, Q, Z, ν, μ, n, m, p  = klf_rlsplit(M, N, fast = fast, finite_infinite = false)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      ν == [1] && μ == [1] && n == 0 && m == 0  && p == 0

M2 = ones(1,1); N2 = zeros(1,1);
M = copy(M2); N = copy(N2);

@time M1, N1, Q, Z, ν, μ, n, m, p  = klf_rlsplit(M, N, fast = fast, finite_infinite = true)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      ν == [1] && μ == [1] && n == 0 && m == 0  && p == 0

M2 = zeros(1,1); N2 = ones(1,1);
M = copy(M2); N = copy(N2);

@time M1, N1, Q, Z, ν, μ, n, m, p  = klf_rlsplit(M, N, fast = fast, finite_infinite = false)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      ν == [] && μ == [] && n == 1 && m == 0  && p == 0

M2 = zeros(1,1); N2 = ones(1,1);
M = copy(M2); N = copy(N2);

@time M1, N1, Q, Z, ν, μ, n, m, p  = klf_rlsplit(M, N, fast = fast, finite_infinite = true)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      ν == [] && μ == [] && n == 1 && m == 0  && p == 0

# Test Suite 1 (Kagstrom)
M2 = [  22  34  31   31  17
        45  45  42   19  29
        39  47  49   26  34
        27  31  26   21  15
        38  44  44   24  30 ];

N2 = [   13  26  25  17  24
        31  46  40  26  37
        26  40  19  25  25
        16  25  27  14  23
        24  35  18  21  22  ];

M = copy(M2); N = copy(N2);

@time M1, N1, Q, Z, ν, μ, n, m, p  = klf_rlsplit(M, N, fast = fast, finite_infinite = false)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      ν == [1] && μ == [2] && n == 3 && m == 0  && p == 1

M = copy(M2); N = copy(N2);

@time M1, N1, Q, Z, ν, μ, n, m, p  = klf_rlsplit(M, N, fast = fast, finite_infinite = true)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      ν == [2] && μ == [1] && n == 3 && m == 1  && p == 0

# Test data for MB04UD
M2 = [
   2.0  0.0  2.0 -2.0
   0.0 -2.0  0.0  2.0
   2.0  0.0 -2.0  0.0
   2.0 -2.0  0.0  2.0];
N2 = [
   1.0  0.0  1.0 -1.0
   0.0 -1.0  0.0  1.0
   1.0  0.0 -1.0  0.0
   1.0 -1.0  0.0  1.0];

M = copy(M2); N = copy(N2);

@time M1, N1, Q, Z, ν, μ, n, m, p  = klf_rlsplit(M, N, fast = fast, finite_infinite = false)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      ν == [] && μ == [] && n == 4 && m == 0  && p == 0

M = copy(M2); N = copy(N2);

@time M1, N1, Q, Z, ν, μ, n, m, p  = klf_rlsplit(M, N, fast = fast, finite_infinite = true)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      ν == [] && μ == [] && n == 4 && m == 0  && p == 0

# Test Suite 2 (Example 2.2.1, Beelen)
M = zeros(14,16);
i1 = [1 2 3 6 7 8 9 10 11 12 13 13 14];
j1 = [4 6 7 8 9 10 11 12 13 14 15 16 16];
aval = [1 1 1 1 1 1 1 1 1 2 3 1 3];
for k = 1:length(i1)
    M[i1[k],j1[k]] = aval[k]
end
N = zeros(14,16);
i1 = [1 2 3 5 6 7 10 12 13 14]; j1 = [3 5 6 8 9 10 13 14 15 16];
for k = 1:length(i1)
    N[i1[k],j1[k]] = 1
end

Q = qr(rand(14,14)).Q; Z = qr(rand(16,16)).Q;
M2 = Q*M*Z; N2 = Q*N*Z;

M = copy(M2); N = copy(N2);

@time M1, N1, Q, Z, ν, μ, n, m, p  = klf_rlsplit(M, N, atol1 = 1.e-7, atol2 = 1.e-7, fast = fast, finite_infinite = false)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      ν == [4, 2, 0 ] && μ == [6, 3, 1]  && n == 6 && m == 0  && p == 2

M = copy(M2); N = copy(N2);

@time M1, N1, Q, Z, ν, μ, n, m, p  = klf_rlsplit(M, N, fast = fast, finite_infinite = true, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      ν == [1, 1, 2, 4] && μ == [0, 1, 2, 3] && n == 6 && m == 4  && p == 0

Ty = Float64
for Ty in (Float64, Complex{Float64})

abstol = sqrt(eps(one(real(Ty))))

mr = 2; nr = mr+1; ni = 4; nf = 10; nl = 4; ml = nl+1;
mM = mr+ni+nf+ml;
nM = nr+ni+nf+nl;
M2 = [rand(Ty,mr,nM);
     zeros(Ty,ni,nr) Matrix{Ty}(I,ni,ni) rand(Ty,ni,nf+nl);
     zeros(Ty,nf,nr+ni) rand(Ty,nf,nf+nl);
     zeros(Ty,ml,nr+ni+nf) rand(Ty,ml,nl)]
N2 = [rand(Ty,mr,nM);
     zeros(Ty,ni,nr) diagm(1 => ones(Ty,ni-1)) rand(Ty,ni,nf+nl);
     zeros(Ty,nf,nr+ni) rand(Ty,nf,nf+nl);
     zeros(Ty,ml,nr+ni+nf) rand(Ty,ml,nl)]
Q = qr(rand(Ty,mM,mM)).Q;
Z = qr(rand(Ty,nM,nM)).Q;
M2 = Q*M2*Z;
N2 = Q*N2*Z;

M = copy(M2); N = copy(N2);
atol1 = 1.e-7; atol2 = 1.e-7; rtol = 1.e-7;
@time M1, N1, Q1, Z1, νr, μr, νi, nfe, νl, μl = klf(M, N, finite_infinite=true,fast = true, atol1 = 0*atol1, atol2 = 0*atol2, rtol = rtol)
@test norm(Q1'*M2*Z1-M1) < atol1 &&
      norm(Q1'*N2*Z1-N1) < atol2 &&
      νr == [ones(Int,mr); 0] && μr == ones(Int,nr) &&
      νi == ones(Int,ni) && νl == ones(Int,ml) && μl == [0; ones(Int,nl)] && nfe == nf

M = copy(M2); N = copy(N2);

@time M1, N1, Q, Z, ν, μ, n, m, p  = klf_rlsplit(M, N, atol1 = 1.e-7, atol2 = 1.e-7, fast = fast, finite_infinite = false)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      ν == [2, 2, 1, 1] && μ == [2, 2, 2, 1]  && n == 14 && m == 0  && p == 1

M = copy(M2); N = copy(N2);

@time M1, N1, Q, Z, ν, μ, n, m, p  = klf_rlsplit(M, N, atol1 = 1.e-7, atol2 = 1.e-7, fast = fast, finite_infinite = true)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      ν == [1, 2, 2, 2, 2] && μ == [0, 2, 2, 2, 2]  && n == 12 && m == 1  && p == 0
end

end
end

@testset "klf_left and klf_right" begin

fast = true
for fast in (false,true)

M2 = zeros(0,0); N2 = zeros(0,0);
M = copy(M2); N = copy(N2);

@time M1, N1, Q, Z, ν, μ, nf, νl, μl = klf_left(M, N, fast = fast)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      ν == [] && μ == [] && νl == [] && μl == []  && nf == 0

M2 = zeros(0,0); N2 = zeros(0,0);
M = copy(M2); N = copy(N2);

@time M1, N1, Q, Z, νr, μr, nf, ν, μ  = klf_right(M, N, fast = fast)
@test norm(Q'*M2*Z-M) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N) < sqrt(eps(1.)) &&
      νr == [] && μr == [] && ν == [] && μ == []  && nf == 0

M2 = zeros(3,0); N2 = zeros(3,0);
M = copy(M2); N = copy(N2);

@time M1, N1, Q, Z, ν, μ, nf, νl, μl = klf_left(M, N, fast = fast)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      ν == [] && μ == [] && νl == [3] && μl == [0]  && nf == 0

M2 = zeros(3,0); N2 = zeros(3,0);
M = copy(M2); N = copy(N2);

@time M1, N1, Q, Z, νr, μr, nf, ν, μ  = klf_right(M, N, fast = fast)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      νr == [] && μr == [] && ν == [3] && μ == [0] && nf == 0

M2 = zeros(0,3); N2 = zeros(0,3);
M = copy(M2); N = copy(N2);

@time M1, N1, Q, Z, ν, μ, nf, νl, μl = klf_left(M, N, fast = fast)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      ν == [0] && μ == [3] && νl == [] && μl == [] && nf == 0

M2 = zeros(0,3); N2 = zeros(0,3);
M = copy(M2); N = copy(N2);
@time M1, N1, Q, Z, νr, μr, nf, ν, μ  = klf_right(M, N, fast = fast)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      νr == [0] && μr == [3] && ν == [] && μ == [] && nf == 0

N2 = zeros(2,1); M2 = [zeros(1,1); ones(1,1)];
M = copy(M2); N = copy(N2);
@time M1, N1, Q, Z, ν, μ, nf, νl, μl = klf_left(M, N, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      ν == [1] && μ == [1] && νl == [1] && μl == [0] && nf == 0

N2 = zeros(2,1); M2 = [zeros(1,1); ones(1,1)];
M = copy(M2); N = copy(N2);
@time M1, N1, Q, Z, νr, μr, nf, ν, μ  = klf_right(M, N, fast = fast)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      νr == [] && μr == [] && ν == [2] && μ == [1] && nf == 0

M2 = zeros(1,1); N2 = zeros(1,1);
M = copy(M2); N = copy(N2);
@time M1, N1, Q, Z, ν, μ, nf, νl, μl = klf_left(M, N, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      ν == [0] && μ == [1]  && νl == [1] && μl == [0] && nf == 0

M2 = zeros(1,1); N2 = zeros(1,1);
M = copy(M2); N = copy(N2);
@time M1, N1, Q, Z, νr, μr, nf, ν, μ  = klf_right(M, N, fast = fast)
@test norm(Q'*M2*Z-M1) < sqrt(eps(1.)) &&
      norm(Q'*N2*Z-N1) < sqrt(eps(1.)) &&
      νr == [0] && μr == [1] && ν == [1] && μ == [0] && nf == 0

# Test Suite 1 (Kagstrom)
M2 = [  22  34  31   31  17
        45  45  42   19  29
        39  47  49   26  34
        27  31  26   21  15
        38  44  44   24  30 ];

N2 = [   13  26  25  17  24
        31  46  40  26  37
        26  40  19  25  25
        16  25  27  14  23
        24  35  18  21  22  ];

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, ν, μ, nf, νl, μl = klf_left(M, N, fast = fast)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      ν == [1] && μ == [2] && νl == [1] && μl == [0] && nf == 3

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, νr, μr, nf, ν, μ = klf_right(M, N, fast = fast)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [0] && μr == [1] && ν == [2] && μ == [1] && nf == 3

# Test data for MB04UD
M2 = [
   2.0  0.0  2.0 -2.0
   0.0 -2.0  0.0  2.0
   2.0  0.0 -2.0  0.0
   2.0 -2.0  0.0  2.0];
N2 = [
   1.0  0.0  1.0 -1.0
   0.0 -1.0  0.0  1.0
   1.0  0.0 -1.0  0.0
   1.0 -1.0  0.0  1.0];

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, ν, μ, nf, νl, μl = klf_left(M, N, fast = fast)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      ν == [ ] && μ == [ ]  && νl == [] && μl == [] && nf == 4

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, νr, μr, nf, ν, μ  = klf_right(M, N, fast = fast)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      ν == [ ] && μ == [ ]  && νl == [] && μl == [] && nf == 4

# Test Suite 2 (Example 2.2.1, Beelen)
M = zeros(14,16);
i1 = [1 2 3 6 7 8 9 10 11 12 13 13 14];
j1 = [4 6 7 8 9 10 11 12 13 14 15 16 16];
aval = [1 1 1 1 1 1 1 1 1 2 3 1 3];
for k = 1:length(i1)
    M[i1[k],j1[k]] = aval[k]
end
N = zeros(14,16);
i1 = [1 2 3 5 6 7 10 12 13 14]; j1 = [3 5 6 8 9 10 13 14 15 16];
for k = 1:length(i1)
    N[i1[k],j1[k]] = 1
end

Q = qr(rand(14,14)).Q; Z = qr(rand(16,16)).Q;
M2 = Q*M*Z; N2 = Q*N*Z;

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, ν, μ, nf, νl, μl = klf_left(M, N, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      ν == [4, 2, 0 ] && μ == [6, 3, 1] && νl == [1, 1, 1, 2] && μl == [0, 1, 1, 1] && nf == 3

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, ν, μ, nf, νl, μl = klf_left(M, N, fast = fast, ut = true, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      ν == [4, 2, 0 ] && μ == [6, 3, 1] && νl == [1, 1, 1, 2] && μl == [0, 1, 1, 1] && nf == 3

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, n, m, νi, νl, μl = klf_leftinf(M, N, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      n == 6 && m == 4 && νl == [1, 1, 1, 2] && μl == [0, 1, 1, 1] && νi == [2, 1]

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, n, m, νi, νl, μl = klf_leftinf(M, N, fast = fast, ut = true, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      n == 6 && m == 4 && νl == [1, 1, 1, 2] && μl == [0, 1, 1, 1] && νi == [2, 1]

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, νr, μr, nf, ν, μ  = klf_right(M, N, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [2, 1, 0] && μr == [4, 2, 1]  && ν == [1, 1, 2, 4] && μ == [0, 1, 2, 3] && nf == 3

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, νr, μr, nf, ν, μ  = klf_right(M, N, fast = fast, ut = true, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [2, 1, 0] && μr == [4, 2, 1]  && ν == [1, 1, 2, 4] && μ == [0, 1, 2, 3] && nf == 3

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, νr, μr, νi, n, p  = klf_rightinf(M, N, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [2, 1, 0] && μr == [4, 2, 1]  && νi == [1, 2] && n == 6 && p == 2

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, νr, μr, νi, n, p  = klf_rightinf(M, N, fast = fast, ut = true, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [2, 1, 0] && μr == [4, 2, 1]  && νi == [1, 2] && n == 6 && p == 2

M = copy(M2'); N = copy(N2');
@time M1, N1, Q1, Z1, ν, μ, nf, νl, μl = klf_left(M, N, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2'*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2'*Z1-N1) < sqrt(eps(1.)) &&
      ν == [3, 2, 1, 0] && μ == [4, 2, 1, 1] && νl == [1, 2, 4] && μl == [0, 1, 2] && nf == 3

M = copy(M2'); N = copy(N2');
@time M1, N1, Q1, Z1, νr, μr, nf, ν, μ  = klf_right(M, N, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2'*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2'*Z1-N1) < sqrt(eps(1.)) &&
      νr == [1, 1, 1, 0] && μr == [2, 1, 1, 1] && ν == [1, 3, 6] && μ == [0, 2, 4] && nf == 3

# intermediary computation
M2 = [
    1.0  -1.14463e-18  -1.47479e-15  4.87132e-15  6.96416e-31  -2.32652e-31  -1.50216e-16   1.69785e-16   1.42103e-16   7.77705e-17
    0.0   1.0           0.0          0.0          0.0           0.0          -2.34988e-16  -3.55922e-16  -6.84622e-17  -6.11787e-17
    0.0   0.0           1.0          0.0          0.0           0.0           4.90878e-16   3.80428e-16   3.62292e-16   3.74223e-16
    0.0   0.0           0.0          1.0          0.0           0.0           8.12392e-16   5.94559e-16   1.17495e-16   6.02872e-16
    0.0   0.0           0.0          0.0          0.0           0.0           1.0           5.10024e-16  -2.95823e-31  -2.69854e-16
    0.0   0.0           0.0          0.0          0.0           0.0           0.0           1.0           0.0           5.41115e-16
]
N2 = [
    0.0  0.0  0.0  0.0  0.0  0.0  -1.29895e-15   4.68167e-15   1.51642e-15  -2.88014e-17
    0.0  0.0  0.0  0.0  0.0  0.0   0.689011     -0.0434045     0.72345      -3.72429e-16
    0.0  0.0  0.0  0.0  0.0  0.0   0.723189      0.10666      -0.682364     -7.51122e-18
    0.0  0.0  0.0  0.0  0.0  0.0  -0.0475455     0.993348      0.10488       1.43459e-17
    0.0  0.0  0.0  0.0  0.0  0.0   0.0           0.0           0.0          -5.10024e-16
    0.0  0.0  0.0  0.0  0.0  0.0   0.0           0.0           0.0          -1.0
]

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, νr, μr, nf, νl, μl  = klf_right(M, N, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [2, 1, 0] && μr == [4, 2, 1] &&  νl == [1, 2] && μl == [1, 2] && nf == 0

# intermediary result
M2 = [
1.0  -1.51959e-15   4.9582e-15   -7.05588e-17   6.96416e-31  -2.32652e-31   2.38551e-16  -2.02364e-16   5.79165e-17   7.77705e-17
0.0   1.0           2.77556e-17  -1.73472e-18  -2.76763e-47  -1.35446e-47  -1.43825e-16  -3.50328e-16  -5.85653e-16   3.74223e-16
0.0   0.0          -1.0           2.77556e-17   4.26419e-48   9.29779e-48   5.70351e-16   5.60542e-16   6.19106e-16  -6.02872e-16
0.0   0.0           0.0          -1.0           7.98872e-48   1.74576e-47  -1.64948e-16  -3.84582e-16  -1.94459e-16   6.11787e-17
0.0   0.0           0.0           0.0           0.0           0.0           0.731013      0.0979        0.675304      2.52908e-16
0.0   0.0           0.0           0.0           0.0           0.0           0.0          -0.989654      0.143472      5.77433e-16
]
N2 = [
    0.0  0.0  0.0  0.0  0.0  0.0   1.40857e-15  -4.84718e-15  -4.04657e-17  -4.16334e-17
    0.0  0.0  0.0  0.0  0.0  0.0  -1.0           2.21239e-7    5.83691e-7    2.19061e-17
    0.0  0.0  0.0  0.0  0.0  0.0   0.0           1.0           2.90233e-7   -1.24119e-17
    0.0  0.0  0.0  0.0  0.0  0.0   0.0           0.0           1.0           3.44402e-16
    0.0  0.0  0.0  0.0  0.0  0.0   0.0           0.0           0.0           0.145907
    0.0  0.0  0.0  0.0  0.0  0.0   0.0           0.0           0.0          -0.989298
]
M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, νr, μr, nf, νl, μl  = klf_right(M, N, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [2, 1, 0] && μr == [4, 2, 1] &&  νl == [1, 2] && μl == [1, 2] && nf == 0

# Example  Karathanasi and Karampetakis IMA J. Math. Contr Inform. (2019)
N2 = [
1. 1 0 0 1 0 -1 0
1 1 1 1 0 0 0 0
-1 1 0 0 0 0 0 0
1 0 0 0 0 -1 0 0
0 -2 0 0 0 1 1 0
1 0 0 0 0 -1 0 0
]
M2 = [
2. 1 0 0 0 -2 0 0
0 2 1 0 0 0 -2 0
-2 1 0 0 0 1 -1 1
3 -1 1 1 0 -2 2 0
0 -2 0 0 0 1 2 -1
2 -1 0 0 0 -2 1 0
]

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, ν, μ, nf, νl, μl = klf_left(M, N, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      ν == [3, 1 ] && μ == [3, 3] && νl == [] && μl == [] && nf == 2

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, νr, μr, nf, ν, μ  = klf_right(M, N, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [2, 0 ] && μr == [2, 2] && ν == [1, 1] && μ == [1, 1] && nf == 2

for Ty in (Float64, Complex{Float64})

abstol = sqrt(eps(one(real(Ty))))

# generic cases
m = 1; n = 2;
M2 = rand(Ty,m,n)
N2 = rand(Ty,m,n)
M = copy(M2); N = copy(N2);
@time M1, N1, Q, Z, ν, μ, nf, νl, μl = klf_left(M, N, fast = fast)
@test norm(Q'*M2*Z-M1) < abstol &&
      norm(Q'*N2*Z-N1) < abstol &&
      ν == [1, 0 ] && μ == [1, 1] && νl == [] && μl == [] && nf == 0

M = copy(M2); N = copy(N2);
@time M1, N1, Q, Z, νr, μr, nf, ν, μ = klf_right(M, N, fast = fast)
@test norm(Q'*M2*Z-M1) < abstol &&
      norm(Q'*N2*Z-N1) < abstol &&
      νr == [1, 0 ] && μr == [1, 1] && ν == [] && μ == [] && nf == 0

m = 3; n = 5;
M2 = rand(Ty,m,n)
N2 = rand(Ty,m,n)
M = copy(M2); N = copy(N2);
@time M1, N1, Q, Z, ν, μ, nf, νl, μl = klf_left(M, N, fast = fast)
@test norm(Q'*M2*Z-M1) < abstol &&
      norm(Q'*N2*Z-N1) < abstol &&
      ν == [2, 1, 0] && μ == [2, 2, 1] && νl == [] && μl == [] && nf == 0

m = 3; n = 5;
M2 = rand(Ty,m,n)
N2 = rand(Ty,m,n)
M = copy(M2); N = copy(N2);
@time M1, N1, Q, Z, νr, μr, nf, ν, μ = klf_right(M, N, fast = fast)
@test norm(Q'*M2*Z-M1) < abstol &&
      norm(Q'*N2*Z-N1) < abstol &&
      νr == [2, 1, 0] && μr == [2, 2, 1] && ν == [] && μ == [] && nf == 0

m = 5; n = 3;
M2 = rand(Ty,m,n)
N2 = rand(Ty,m,n)
M = copy(M2); N = copy(N2);
@time M1, N1, Q, Z, ν, μ, nf, νl, μl = klf_left(M, N, fast = fast);
@test norm(Q'*M2*Z-M1) < abstol &&
      norm(Q'*N2*Z-N1) < abstol &&
      ν == [] && μ == [] && νl == [1, 2, 2] && μl == [0, 1, 2] && nf == 0

m = 5; n = 3;
M2 = rand(Ty,m,n)
N2 = rand(Ty,m,n)
M = copy(M2); N = copy(N2);
@time M1, N1, Q, Z, νr, μr, nf, ν, μ = klf_right(M, N, fast = fast);
@test norm(Q'*M2*Z-M1) < abstol &&
      norm(Q'*N2*Z-N1) < abstol &&
      νr == [] && μr == [] && ν == [1, 2, 2] && μ == [0, 1, 2] && nf == 0

m = 3; n = 5; r = 2;
M2 = rand(Ty,m,n)
N2 = rand(Ty,m,m)*[rand(Ty,r,r) zeros(Ty,r,n-r); zeros(Ty,m-r,n)]*rand(Ty,n,n)
M = copy(M2); N = copy(N2);
@time M1, N1, Q, Z, ν, μ, nf, νl, μl = klf_left(M, N; fast = fast);
@test norm(Q'*M2*Z-M1) < abstol &&
      norm(Q'*N2*Z-N1) < abstol &&
      ν == [3, 0] && μ == [3, 2] && νl == [] && μl == [] && nf == 0

m = 3; n = 5; r = 2;
M2 = rand(Ty,m,n)
N2 = rand(Ty,m,m)*[rand(Ty,r,r) zeros(Ty,r,n-r); zeros(Ty,m-r,n)]*rand(Ty,n,n)
M = copy(M2); N = copy(N2);
@time M1, N1, Q, Z, νr, μr, nf, ν, μ = klf_right(M, N; fast = fast)
@test norm(Q'*M2*Z-M1) < abstol &&
      norm(Q'*N2*Z-N1) < abstol &&
      νr == [2, 0] && μr == [2, 2] && ν == [1] && μ == [1] && nf == 0

m = 3; n = 5; r = 2;
M2 = rand(Ty,m,n)
N2 = rand(Ty,m,m)*[rand(Ty,r,r) zeros(Ty,r,n-r); zeros(Ty,m-r,n)]*rand(Ty,n,n)
M = copy(M2); N = copy(N2);
@time M1, N1, Q, Z, νr, μr, nf, ν, μ = klf_right(M, N; fast = fast)
@test norm(Q'*M2*Z-M1) < abstol &&
      norm(Q'*N2*Z-N1) < abstol &&
      νr == [2, 0] && μr == [2, 2] && ν == [1] && μ == [1] && nf == 0

m = 3; n = 5; r = 2;
M2 = rand(Ty,m,n)
N2 = rand(Ty,m,m)*[rand(Ty,r,r) zeros(Ty,r,n-r); zeros(Ty,m-r,n)]*rand(Ty,n,n)
M = copy(M2); N = copy(N2);
@time M1, N1, Q, Z, ν, μ, nf, νl, μl = klf_left(M, N; fast = fast)
@test norm(Q'*M2*Z-M1) < abstol &&
      norm(Q'*N2*Z-N1) < abstol &&
      ν == [3, 0] && μ == [3, 2] && νl == [] && μl == [] && nf == 0

M2 = rand(Ty,7,7)
N2 = [zeros(Ty,4,3) rand(Ty,4,4); zeros(Ty,3,7)]
M = copy(M2); N = copy(N2);
@time M1, N1, Q, Z, ν, μ, nf, νl, μl = klf_left(M, N, fast = fast);
@test norm(Q'*M2*Z-M1) < abstol &&
      norm(Q'*N2*Z-N1) < abstol &&
      ν == [3] && μ == [3] && νl == [] && μl == [] && nf == 4

M2 = rand(Ty,7,7)
N2 = [zeros(Ty,4,3) rand(Ty,4,4); zeros(Ty,3,7)]
M = copy(M2); N = copy(N2);
@time M1, N1, Q, Z, νr, μr, nf, ν, μ = klf_right(M, N, fast = fast);
@test norm(Q'*M2*Z-M1) < abstol &&
      norm(Q'*N2*Z-N1) < abstol &&
      νr == [] && μr == [] && ν == [3] && μ == [3] && nf == 4

end
end

@testset "klf" begin

fast = true
for fast in (false,true)

M2 = zeros(0,0); N2 = zeros(0,0);
M = copy(M2); N = copy(N2);

@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl  = klf(M, N, fast = fast)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [] && μr == [] && νi == [] && νl == [] && μl == []  && nf == 0

M = copy(M2); N = copy(N2);

@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl  = klf(M, N, fast = fast, finite_infinite = true)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [] && μr == [] && νi == [] && νl == [] && μl == []  && nf == 0

M2 = zeros(3,0); N2 = zeros(3,0);
M = copy(M2); N = copy(N2);

@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl  = klf(M, N, fast = fast)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [] && μr == [] && νi == [] && νl == [3] && μl == [0]  && nf == 0

M2 = zeros(0,3); N2 = zeros(0,3);
M = copy(M2); N = copy(N2);

@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl  = klf(M, N, fast = fast)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [0] && μr == [3] && νi == [] && νl == [] && μl == []  && nf == 0

M2 = ones(1,1); N2 = zeros(1,1);
M = copy(M2); N = copy(N2);

@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl  = klf(M, N, fast = fast)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [] && μr == [] && νi == [1] && νl == [] && μl == []  && nf == 0

M2 = zeros(1,1); N2 = ones(1,1);
M = copy(M2); N = copy(N2);

@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl  = klf(M, N, fast = fast)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [] && μr == [] && νi == [] && νl == [] && μl == []  && nf == 1

N2 = zeros(2,1); M2 = [zeros(1,1); ones(1,1)];
M = copy(M2); N = copy(N2);

@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl  = klf(M, N, fast = fast)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [] && μr == [] && νi == [1] && νl == [1] && μl == [0]  && nf == 0

M2 = zeros(1,1); N2 = zeros(1,1);
M = copy(M2); N = copy(N2);

@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl  = klf(M, N, fast = fast)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [0] && μr == [1] && νi == [] && νl == [1] && μl == [0]  && nf == 0

# Test Suite 1 (Kagstrom)
M2 = [  22  34  31   31  17
        45  45  42   19  29
        39  47  49   26  34
        27  31  26   21  15
        38  44  44   24  30 ];

N2 = [   13  26  25  17  24
        31  46  40  26  37
        26  40  19  25  25
        16  25  27  14  23
        24  35  18  21  22  ];

M = copy(M2); N = copy(N2);

@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl  = klf(M, N, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [0] && μr == [1] && νi == [1] && νl == [1] && μl == [0]  && nf == 3

M = copy(M2); N = copy(N2);

@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl  = klf(M, N, fast = fast, finite_infinite = true)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [0] && μr == [1] && νi == [1] && νl == [1] && μl == [0]  && nf == 3

M = copy(M2); N = copy(N2);

@time N1, M1, Q1, Z1, νr, μr, νi, nf, νl, μl  = klf(N, M, fast = fast)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [0] && μr == [1] && νi == [1, 1] && νl == [1] && μl == [0]  && nf == 2

M = copy(M2); N = copy(N2);

@time N1, M1, Q1, Z1, νr, μr, νi, nf, νl, μl  = klf(N, M, fast = fast, finite_infinite = true)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [0] && μr == [1] && νi == [1, 1] && νl == [1] && μl == [0]  && nf == 2

# Test data for MB04UD
M2 = [
   2.0  0.0  2.0 -2.0
   0.0 -2.0  0.0  2.0
   2.0  0.0 -2.0  0.0
   2.0 -2.0  0.0  2.0];
N2 = [
   1.0  0.0  1.0 -1.0
   0.0 -1.0  0.0  1.0
   1.0  0.0 -1.0  0.0
   1.0 -1.0  0.0  1.0];

M = copy(M2); N = copy(N2);

@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl  = klf(M, N, fast = fast)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [] && μr == [] && νi == [] && νl == [] && μl == []  && nf == 4

# Test data for MB04VD
M2 = [
   1.0  0.0 -1.0  0.0
   1.0  1.0  0.0 -1.0];
N2 = [
   0.0 -1.0  0.0  0.0
   0.0 -1.0  0.0  0.0];

M = copy(M2); N = copy(N2);

@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl  = klf(M, N, fast = fast)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [1, 0] && μr == [2, 1] && νi == [1] && νl == [] && μl == []  && nf == 0

# Test Suite 2 (Example 2.2.1, Beelen)
M = zeros(14,16);
i1 = [1 2 3 6 7 8 9 10 11 12 13 13 14];
j1 = [4 6 7 8 9 10 11 12 13 14 15 16 16];
aval = [1 1 1 1 1 1 1 1 1 2 3 1 3];
for k = 1:length(i1)
    M[i1[k],j1[k]] = aval[k]
end
N = zeros(14,16);
i1 = [1 2 3 5 6 7 10 12 13 14]; j1 = [3 5 6 8 9 10 13 14 15 16];
for k = 1:length(i1)
    N[i1[k],j1[k]] = 1
end

Q = qr(rand(14,14)).Q; Z = qr(rand(16,16)).Q;
M2 = Q*M*Z; N2 = Q*N*Z;

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = klf(M, N, fast = fast, ut = false, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [2, 1, 0] && μr == [4, 2, 1] && νi == [1, 2] && νl == [1, 1, 1, 2] && μl == [0, 1, 1, 1] && nf == 3

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = klf(M, N, fast = fast, ut = true, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [2, 1, 0] && μr == [4, 2, 1] && νi == [1, 2] && νl == [1, 1, 1, 2] && μl == [0, 1, 1, 1] && nf == 3

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl  = klf(M, N, fast = fast, finite_infinite = true, ut = false, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [2, 1, 0] && μr == [4, 2, 1] && νi == [2, 1] && νl == [1, 1, 1, 2] && μl == [0, 1, 1, 1] && nf == 3

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl  = klf(M, N, fast = fast, finite_infinite = true, ut = true, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [2, 1, 0] && μr == [4, 2, 1] && νi == [2, 1] && νl == [1, 1, 1, 2] && μl == [0, 1, 1, 1] && nf == 3

M = copy(M2'); N = copy(N2');
@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = klf(M, N, fast = fast, ut = false, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2'*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2'*Z1-N1) < sqrt(eps(1.)) &&
      νr == [1, 1, 1, 0] && μr == [2, 1, 1, 1] && νi == [1, 2] && νl == [1, 2, 4] && μl == [0, 1, 2] && nf == 3

M = copy(M2'); N = copy(N2');
@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = klf(M, N, fast = fast, ut = true, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2'*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2'*Z1-N1) < sqrt(eps(1.)) &&
      νr == [1, 1, 1, 0] && μr == [2, 1, 1, 1] && νi == [1, 2] && νl == [1, 2, 4] && μl == [0, 1, 2] && nf == 3

M = copy(M2'); N = copy(N2');
@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = klf(M, N, fast = fast, finite_infinite = true, ut = false, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2'*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2'*Z1-N1) < sqrt(eps(1.)) &&
      νr == [1, 1, 1, 0] && μr == [2, 1, 1, 1] && νi == [2, 1] && νl == [1, 2, 4] && μl == [0, 1, 2] && nf == 3

M = copy(M2'); N = copy(N2');
@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = klf(M, N, fast = fast, finite_infinite = true, ut = true, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2'*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2'*Z1-N1) < sqrt(eps(1.)) &&
      νr == [1, 1, 1, 0] && μr == [2, 1, 1, 1] && νi == [2, 1] && νl == [1, 2, 4] && μl == [0, 1, 2] && nf == 3

# intermediary computation
M2 = [
    1.0  -1.14463e-18  -1.47479e-15  4.87132e-15  6.96416e-31  -2.32652e-31  -1.50216e-16   1.69785e-16   1.42103e-16   7.77705e-17
    0.0   1.0           0.0          0.0          0.0           0.0          -2.34988e-16  -3.55922e-16  -6.84622e-17  -6.11787e-17
    0.0   0.0           1.0          0.0          0.0           0.0           4.90878e-16   3.80428e-16   3.62292e-16   3.74223e-16
    0.0   0.0           0.0          1.0          0.0           0.0           8.12392e-16   5.94559e-16   1.17495e-16   6.02872e-16
    0.0   0.0           0.0          0.0          0.0           0.0           1.0           5.10024e-16  -2.95823e-31  -2.69854e-16
    0.0   0.0           0.0          0.0          0.0           0.0           0.0           1.0           0.0           5.41115e-16
]
N2 = [
    0.0  0.0  0.0  0.0  0.0  0.0  -1.29895e-15   4.68167e-15   1.51642e-15  -2.88014e-17
    0.0  0.0  0.0  0.0  0.0  0.0   0.689011     -0.0434045     0.72345      -3.72429e-16
    0.0  0.0  0.0  0.0  0.0  0.0   0.723189      0.10666      -0.682364     -7.51122e-18
    0.0  0.0  0.0  0.0  0.0  0.0  -0.0475455     0.993348      0.10488       1.43459e-17
    0.0  0.0  0.0  0.0  0.0  0.0   0.0           0.0           0.0          -5.10024e-16
    0.0  0.0  0.0  0.0  0.0  0.0   0.0           0.0           0.0          -1.0
]

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = klf(M, N, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [2, 1, 0] && μr == [4, 2, 1] && νi == [1, 2] && νl == [] && μl == [] && nf == 0

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = klf(M', N', fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2'*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2'*Z1-N1) < sqrt(eps(1.)) &&
      μl == [0, 1, 2] && νl == [1, 2, 4] && νi == [1, 2] && νr == [] && μr == [] && nf == 0

# Example  Karathanasi and Karampetakis IMA J. Math. Contr Inform. (2019)
N2 = [
1. 1 0 0 1 0 -1 0
1 1 1 1 0 0 0 0
-1 1 0 0 0 0 0 0
1 0 0 0 0 -1 0 0
0 -2 0 0 0 1 1 0
1 0 0 0 0 -1 0 0
]
M2 = [
2. 1 0 0 0 -2 0 0
0 2 1 0 0 0 -2 0
-2 1 0 0 0 1 -1 1
3 -1 1 1 0 -2 2 0
0 -2 0 0 0 1 2 -1
2 -1 0 0 0 -2 1 0
]

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = klf(M, N, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [2, 0] && μr == [2, 2] && νi == [1, 1] && νl == [] && μl == [] && nf == 2

M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = klf(M', N', fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2'*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2'*Z1-N1) < sqrt(eps(1.)) &&
      νl == [2, 2] && μl == [0, 2] && νi == [1, 1] && νr == [] && μr == [] && nf == 2

Ty = Float64;
#for Ty in (Float64, Float32, Complex{Float64},  Complex{Float32})
for Ty in (Float64, Complex{Float64})

abstol = sqrt(eps(one(real(Ty))))

# generic cases
m = 1; n = 2;
M2 = rand(Ty,m,n)
N2 = rand(Ty,m,n)
M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = klf(M, N, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [1, 0] && μr == [1, 1] && νi == [] && νl == [] && μl == [] && nf == 0

m = 3; n = 5;
M2 = rand(Ty,m,n)
N2 = rand(Ty,m,n)
M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = klf(M, N, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [2, 1, 0] && μr == [2, 2, 1] && νi == [] && νl == [] && μl == [] && nf == 0

m = 5; n = 3;
M2 = rand(Ty,m,n)
N2 = rand(Ty,m,n)
M = copy(M2); N = copy(N2);
@time M1, N1, Q1, Z1, νr, μr, νi, nf, νl, μl = klf(M, N, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q1'*M2*Z1-M1) < sqrt(eps(1.)) &&
      norm(Q1'*N2*Z1-N1) < sqrt(eps(1.)) &&
      νr == [] && μr == [] && νi == [] && νl == [1, 2, 2] && μl == [0, 1, 2] && nf == 0

m = 3; n = 5; r = 2;
M2 = rand(Ty,m,n)
N2 = rand(Ty,m,m)*[rand(Ty,r,r) zeros(Ty,r,n-r); zeros(Ty,m-r,n)]*rand(Ty,n,n)
M = copy(M2); N = copy(N2);
@time M1, N1, Q, Z, νr, μr, νi, nf, νl, μl = klf(M, N, fast = fast, atol1 = 1.e-7, atol2 = 1.e-7)
@test norm(Q'*M2*Z-M1) < abstol &&
      norm(Q'*N2*Z-N1) < abstol &&
      νr == [2, 0] && μr == [2, 2] && νi == [1] && νl == [] && μl == [] && nf == 0

mr = 2; nr = mr+1; ni = 4; nf = 10; nl = 4; ml = nl+1;
mM = mr+ni+nf+ml;
nM = nr+ni+nf+nl;
M2 = [rand(Ty,mr,nM);
     zeros(Ty,ni,nr) Matrix{Ty}(I,ni,ni) rand(Ty,ni,nf+nl);
     zeros(Ty,nf,nr+ni) rand(Ty,nf,nf+nl);
     zeros(Ty,ml,nr+ni+nf) rand(Ty,ml,nl)]
N2 = [rand(Ty,mr,nM);
     zeros(Ty,ni,nr) diagm(1 => ones(Ty,ni-1)) rand(Ty,ni,nf+nl);
     zeros(Ty,nf,nr+ni) rand(Ty,nf,nf+nl);
     zeros(Ty,ml,nr+ni+nf) rand(Ty,ml,nl)]
# Q = qr(rand(Ty,mM,mM)).Q;
# Z = qr(rand(Ty,nM,nM)).Q;
# M2 = Q*M2*Z;
# N2 = Q*N2*Z;

M = copy(M2); N = copy(N2);
atol1 = 1.e-6; atol2 = 1.e-7; rtol = 1.e-7;
@time M1, N1, Q1, Z1, νr, μr, νi, nfe, νl, μl = klf(M, N, finite_infinite=true,fast = fast, atol1 = atol1, atol2 = 0*atol2, rtol = rtol)
@test norm(Q1'*M2*Z1-M1) < atol1 &&
      norm(Q1'*N2*Z1-N1) < atol2 &&
      νr == [ones(Int,mr); 0] && μr == ones(Int,nr) &&
      νi == ones(Int,ni) && νl == ones(Int,ml) && μl == [0; ones(Int,nl)] && nfe == nf

M = copy(M2); N = copy(N2);
atol1 = 1.e-6; atol2 = 1.e-7; rtol = 1.e-7;
@time M1, N1, Q1, Z1, νr, μr, νi, nfe, νl, μl = klf(M, N, finite_infinite=false,fast = fast, atol1 = atol1, atol2 = 0*atol2, rtol = rtol)
@test norm(Q1'*M2*Z1-M1) < atol1 &&
      norm(Q1'*N2*Z1-N1) < atol2 &&
      νr == [ones(Int,mr); 0] && μr == ones(Int,nr) &&
      νi == ones(Int,ni) && νl == ones(Int,ml) && μl == [0; ones(Int,nl)] && nfe == nf

end
end

end

end
end
end


text:
using Seismic,PyCall
@pyimport matplotlib.pyplot as plt

nt = 256;
nx = convert(Int64,128);
h = Array(Header,nx);
for ix = 1:nx
  h[ix] = Seismic.InitSeisHeader();
  h[ix].tracenum = ix;
  h[ix].d1 = 0.004;
  h[ix].n1 = nt;
  h[ix].imx = ix;
  h[ix].imy = 0;
  h[ix].ihx = 0;
  h[ix].ihy = 0;
end

dt=0.004;
dx=25;
f_niquest=int(floor(1/dt))
param=["dt"=>dt, "dx"=>dx]

data=rand(nt,nx)

f1 = plt.figure(1)
plt.imshow(data[1:nt,:],cmap="PuOr",aspect="auto",vmin=-1,vmax=1)
plt.title("data before fk filter")
plt.xlabel("X (m)")
plt.ylabel("T (s)")
plt.show()

m=fft(data,1)/sqrt(size(data,1));
m=fft(m,2)/sqrt(size(data,2))
m=fftshift(m)

m=m[int(floor(f_niquest/2)):f_niquest,:]

f11 = plt.figure(11)
plt.imshow(abs(m),cmap="PuOr",aspect="auto",vmin=-1,vmax=1)
plt.title("fk spectrum before filter")
plt.xlabel("K (1/m)")
plt.ylabel("W (1/s)")
plt.show()

# apply FK filter
data_fk,fk_spec,h = FKFilter(data,h,param)

f2 = plt.figure(2)
plt.imshow(data_fk[1:nt,:],cmap="PuOr",aspect="auto",vmin=-1,vmax=1)
plt.title("data after fk filter")
plt.xlabel("X (m)")
plt.ylabel("T (s)")
plt.show()

f21 = plt.figure(21)
plt.imshow(fk_spec,cmap="PuOr",aspect="auto",vmin=-1,vmax=1)
plt.title("fk spectrum after filter")
plt.xlabel("K (1/m)")
plt.ylabel("W (1/s)")
plt.show()

ext: jl


text:
@testset "Career progression report" begin

@testset "function generateCareerProgression" begin
    report = generateCareerProgression( mpSim, "Test", "Sim1", "Sim100",
        "Sim500", "Sim1000" )
    @test all( haskey.( Ref( report ), ["Sim1", "Sim100", "Sim500"] ) ) &&
        !any( haskey.( Ref( report ), ["Test", "Sim1000"] ) )
    @test all( [report["Sim1"][1] == 240, report["Sim100"][1] == 240,
        report["Sim500"][1] == 240] )
    @test all( [
        all( report["Sim1"][2][:, "timeIndex"] .== [0, 24] ),
        all( report["Sim100"][2][:, "timeIndex"] .== [36, 60, 156] ),
        all( report["Sim500"][2][:, "timeIndex"] .==
            [192, 216, 252] )] )
    @test all( [
        all( report["Sim1"][2][:, "transition"] .== ["EW", "B-"] ),
        all( report["Sim100"][2][:, "transition"] .==
            ["EW", "Promotion", "PE"] ),
        all( report["Sim500"][2][:, "transition"] .==
            ["EW", "Promotion", "Promotion"] )] )
end  # @testset "function generateCareerProgression"

end  # @testset "Career progression report"


text:
#import MAGMA
using MAGMA

#using MAGMA:gesvd!
using MAGMA: MagmaAllVec, gesvd!, libmagma

### some JuliaGPU packages, maybe useful (who knows)
using CUDAdrv
using CUDAapi
using CUDAnative
using CuArrays

#
using Test, LinearAlgebra

matrixToTest = rand(Float64, 2, 2)

right_answer = svd(matrixToTest).S
S = right_answer

jobu = MagmaAllVec
jobvt = MagmaAllVec

ldu=2
ldvt=2
lwork=134
success=magmaInit()

U, s, VT, work, info = gesvd!(jobu,jobvt,matrixToTest,ldu,ldvt,lwork)

diff = S .- s
error_value = norm(diff)

@test error_value < 1e-7


text:
# sample from the Muller potential using EM
using Plots
using Printf
using Random
using LinearAlgebra
using ForwardDiff
using StaticArrays

push!(LOAD_PATH,"../src/")

using JuBasicMD

include("potentials.jl")

β = 0.1;
x₀ = [-0.75, 1.5];
Δt = 1e-3;
n_iters = 10^6;

V = x->Muller(x);
cfg = ForwardDiff.GradientConfig(V, x₀);
gradV! = (gradV, x)-> ForwardDiff.gradient!(gradV, V, x, cfg);

sampler = EM(gradV!, β, Δt);

Random.seed!(100);
X₀ = copy(x₀);
sample_trajectory!(X₀, sampler, options=Options(n_iters=n_iters));
@printf("In Place X after %d iterations: (%g, %g)\n",n_iters, X₀[1], X₀[2])

Random.seed!(100);
Xvals = sample_trajectory(x₀, sampler, options=Options(n_iters=n_iters,n_save_iters=n_iters));
X = Xvals[end];
@printf("X after %d iterations: (%g, %g)\n",n_iters, X[1], X[2])

Random.seed!(100);
X_vals = sample_trajectory(x₀, sampler, options=Options(n_iters=n_iters));
histogram2d([X[1] for X in X_vals], [X[2] for X in X_vals],normalize=true,color=:viridis)
xlims!(-1.5,1.5)
ylims!(-0.5, 2.0)
xlabel!("x")
ylabel!("y")

max_stars_repo_name: cortner/JuBasicMD

max_issues_repo_name: cortner/JuBasicMD

max_forks_repo_name: cortner/JuBasicMD


text:
DLIB_NAME="dlib-19.4"

# download dlib

run(`rm -fr $(DLIB_NAME).tar.bz2`)
run(`rm -fr $(DLIB_NAME)`)
run(`wget http://dlib.net/files/$(DLIB_NAME).tar.bz2`)
run(`tar -vxjf $(DLIB_NAME).tar.bz2`)
run(`mv $(DLIB_NAME) dlib`)
run(`rm -fr $(DLIB_NAME).tar.bz2`)

# build dlib

cd("dlib/dlib")
mkdir("build")
cd("build")
run(`cmake ..`)
run(`make`)


text:
struct ScaledPerm{T, I} <: GroupsCore.AbstractAlgebra.AbstractPerm
    indices::Vector{Pair{I, T}}
end
function trace_matrix_representative(p::ScaledPerm)
    return sum(pair.second for (i, pair) in enumerate(p.indices) if pair.first == i)
end
function Base.inv(p::ScaledPerm{T,I}) where {T,I}
    indices = Vector{Pair{I,T}}(undef, length(p.indices))
    for i in eachindex(indices)
        pair = p.indices[i]
        indices[pair.first] = i => inv(pair.second)
    end
    return ScaledPerm(indices)
end

function PermutationGroups.order(p::ScaledPerm)
    cur = p
    _one = one(p)
    o = 1
    while cur != _one
        o += 1
        cur *= p
    end
    return o
end
Base.one(p::ScaledPerm{T}) where {T} = ScaledPerm([i => one(T) for i in eachindex(p.indices)])
Base.:(==)(p::ScaledPerm, q::ScaledPerm) = p.indices == q.indices
Base.hash(p::ScaledPerm, u::UInt64) = hash(p.indices, u)
SymbolicWedderburn.degree(p::ScaledPerm) = length(p.indices)
Base.:^(i::Integer, p::ScaledPerm) = p.indices[i]
function SymbolicWedderburn.add_inverse_permutation!(result, val, i::Int, j::Pair)
    result[i, j.first] += val / j.second
end
function SymbolicWedderburn.action_character(conjugacy_cls::AbstractVector{<:AbstractOrbit{<:ScaledPerm}})
    vals = Int[trace_matrix_representative(first(cc)) for cc in conjugacy_cls]
    return SymbolicWedderburn.Character(vals, conjugacy_cls)
end
function Base.:*(p::ScaledPerm, q::ScaledPerm)
    return ScaledPerm(map(eachindex(q.indices)) do i
        pair_q = q.indices[i]
        pair_p = p.indices[pair_q.first]
        pair_p.first => pair_p.second * pair_q.second
    end)
end
Base.:^(p::ScaledPerm, k::Integer) = Base.power_by_squaring(p, k)

function SymbolicWedderburn.permutation(ehom::SymbolicWedderburn.ExtensionHomomorphism{<:AbstractMonomial}, terms::Vector{<:AbstractTerm})
    return ScaledPerm([ehom[monomial(term)] => MP.coefficient(term) for term in terms])
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct ProvidersApi <: SwaggerApi
    client::Swagger.Client
end

"""

Gets the specified resource provider.
Param: resourceProviderNamespace::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: expand::String
Return: Provider
"""
function providersGet(_api::ProvidersApi, resourceProviderNamespace::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Provider, "/subscriptions/{subscriptionId}/providers/{resourceProviderNamespace}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceProviderNamespace", resourceProviderNamespace)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "\$expand", expand)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Gets all resource providers for a subscription.
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: top::Int32
Param: expand::String
Return: ProviderListResult
"""
function providersList(_api::ProvidersApi, api_version::String, subscriptionId::String; top=nothing, expand=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", ProviderListResult, "/subscriptions/{subscriptionId}/providers", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "\$top", top)  # type Int32
    Swagger.set_param(_ctx.query, "\$expand", expand)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Registers a subscription with a resource provider.
Param: resourceProviderNamespace::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Provider
"""
function providersRegister(_api::ProvidersApi, resourceProviderNamespace::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Provider, "/subscriptions/{subscriptionId}/providers/{resourceProviderNamespace}/register", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceProviderNamespace", resourceProviderNamespace)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Unregisters a subscription from a resource provider.
Param: resourceProviderNamespace::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Provider
"""
function providersUnregister(_api::ProvidersApi, resourceProviderNamespace::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Provider, "/subscriptions/{subscriptionId}/providers/{resourceProviderNamespace}/unregister", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceProviderNamespace", resourceProviderNamespace)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export providersGet, providersList, providersRegister, providersUnregister


text:
module Deploy

module Docker

import Revise
import Genie, Genie.FileTemplates

"""
    dockerfile(path::String = "."; user::String = "genie", env::String = "dev",
              filename::String = "Dockerfile", port::Int = 8000, dockerport::Int = 80, force::Bool = false)

Generates a `Dockerfile` optimised for containerizing Genie apps.

# Arguments
- `path::String`: where to generate the file
- `filename::String`: the name of the file (default `Dockerfile`)
- `user::String`: the name of the system user under which the Genie app is run
- `env::String`: the environment in which the Genie app will run
- `host::String`: the local IP of the Genie app inside the container
- `port::Int`: the port of the Genie app inside the container
- `dockerport::Int`: the port to use on the host (used by the `EXPOSE` directive)
- `force::Bool`: if the file already exists, when `force` is `true`, it will be overwritten
"""
function dockerfile(path::String = "."; filename::String = "Dockerfile", user::String = "genie", env::String = "dev",
                    host = "0.0.0.0", port::Int = 8000, dockerport::Int = 80, force::Bool = false,
                    websockets_port::Int = 8001, websockets_dockerport::Int = 8001)
  filename = normpath(joinpath(path, filename))
  isfile(filename) && force && rm(filename)
  isfile(filename) && error("File $(filename) already exists. Use the `force = true` option to overwrite the existing file.")

  open(filename, "w") do io
    write(io, FileTemplates.dockerfile(user = user, env = env, filename = filename, host = host,
                                        port = port, dockerport = dockerport,
                                        websockets_port = websockets_port, websockets_dockerport = websockets_dockerport))
  end

  "Docker file successfully written at $(abspath(filename))" |> println
end

"""
    build(path::String = "."; appname = "genie")

Builds the Docker image based on the `Dockerfile`
"""
function build(path::String = "."; appname::String = "genie")
  `docker build -t "$appname" $path` |> Base.run

  "Docker container successfully built" |> println
end

"""
    run(; containername::String = "genieapp", hostport::Int = 80, containerport::Int = 8000, appdir::String = "/home/genie/app",
        mountapp::Bool = false, image::String = "genie", command::String = "bin/server", rm::Bool = true, it::Bool = true)

Runs the Docker container named `containername`, binding `hostport` and `containerport`.

# Arguments
- `containername::String`: the name of the container of the Genie app
- `hostport::Int`: port to be used on the host for accessing the app
- `containerport::Int`: the port on which the app is running inside the container
- `appdir::String`: the folder where the app is stored within the container
- `mountapp::String`: if true the app from the host will be mounted so that changes on the host will be reflected when accessing the app in the container (to be used for dev)
- `image::String`: the name of the Docker image
- `command::String`: what command to run when starting the app
- `rm::Bool`: removes the container upon exit
- `it::Bool`: runs interactively
"""
function run(; containername::String = "genieapp", hostport::Int = 80, containerport::Int = 8000, appdir::String = "/home/genie/app",
                mountapp::Bool = false, image::String = "genie", command::String = "bin/server", rm::Bool = true, it::Bool = true,
                websockets_hostport::Int = 8001, websockets_containerport::Int = 8001)
  options = []
  it && push!(options, "-it")
  rm && push!(options, "--rm")
  push!(options, "-p")
  push!(options, "$hostport:$containerport")
  push!(options, "-p")
  push!(options, "$websockets_hostport:$websockets_containerport")
  push!(options, "--name")
  push!(options, "$containername")
  if mountapp
    push!(options, "-v")
    push!(options,  "$(pwd()):$appdir")
  end
  push!(options, image)
  isempty(command) || push!(options, command)

  "Starting docker container with `docker run $(join(options, " "))`" |> println

  `docker run $options` |> Base.run
end

end # end module Docker

########

module Heroku

function createapp(appname::String; region::String = "us")
  `heroku create $(lowercase(appname)) --region $region` |> Base.run
end

function push(appname::String; apptype::String = "web")
  `heroku container:push $apptype -a $appname` |> Base.run
end

function release(appname::String; apptype::String = "web")
  `heroku container:release $apptype -a $appname` |> Base.run
end

function open(appname::String)
  `heroku open -a $appname` |> Base.run
end

function login()
  `heroku container:login` |> Base.run
end

function logs(appname::String; lines::Int = 1_000)
  `heroku logs --tail -a $appname -n $lines` |> Base.run
end

end # end module Heroku

end # end module Deploy


text:
using CachedFunctions
using Test

#strain test, different types, a dict.
@testset "basics" begin
    function f!(y,x)
        y[:p] = x[1]
        y[:t] = x[2]+ y[:p]
        for i = 1:10
            y[:p] = cos(y[:p]+1.5)
        end
        y
    end
    y0 = Dict(:t=>0.0,:p=>0.0)
    x = rand(2)
    ff = CachedFunction(f!,x,y0)
    ff(x)
    @test length(cached_methods(ff)) == 1
    @test calls(ff) == 1
    allocate!(ff,BigFloat)
    @test length(cached_methods(ff)) == 2
    ff(big.(rand(2)))
    @test length(cached_methods(ff)) == 2 #check if no new caches were created.
    @test calls(ff) == 2
end


text:
"""
    allele_freq(genotype::Union{Missing, NTuple{N,<:Integer}}) where N
Calculate allele frequency for a single locus of a single sample. Returns a
`Dict` of alleles and their frequencies.
"""
function allele_freq(genotype::Union{Missing, NTuple{N,<:Integer}}) where N
    x === missing && return missing
    d = Dict()
    for allele in genotype
        # sum up alleles
        d[allele] = get!(d, allele, 0) + 1
    end
    total = values(d) |> sum    # sum of all non-missing alleles
    [d[i] = d[i] / total for i in keys(d)]  # allele count / sum
    return d
end

"""
    allele_freq(locus::Vector{<:Union{Missing, NTuple{N,<:Integer}}}) where N
Calculate allele counts for a single locus of a `PopObj`. Returns a `Dict` of
allele's and their frequencies.
"""
function allele_freq(locus::Vector{<:Union{Missing, NTuple{N,<:Integer}}}) where N
    d = Dict()
    for genotype in skipmissing(locus)
        # sum up alleles
        for allele in genotype
            d[allele] = get!(d, allele, 0) + 1
        end
    end
    total = values(d) |> sum    # sum of all non-missing alleles
    [d[i] = d[i] / total for i in keys(d)]  # allele count / sum
    return d
end

"""
    allele_freq(locus::Vector{<:Union{Missing, Tuple{Vararg}}}
Calculate allele counts for a single locus of a `PopObj` with unequal ploidy across
samples. Returns a `Dict` of allele's and their frequencies.
"""
function allele_freq(locus::Vector{<:Union{Missing, Tuple{Vararg}}})
    d = Dict()
    for genotype in skipmissing(locus)
        # sum up alleles
        for allele in genotype
            d[allele] = get!(d, allele, 0) + 1
        end
    end
    total = values(d) |> sum    # sum of all non-missing alleles
    [d[i] = d[i] / total for i in keys(d)]  # allele count / sum
    return d
end

"""
    allele_freq(locus::SubArray{<:Union{Missing,NTuple{N,<:Integer}}}) where N
Calculate allele counts for a single locus of a `PopObj` split by population
using `group()`. Returns a `Dict` of allele's and their frequencies.
"""
function allele_freq(locus::SubArray{<:Union{Missing,NTuple{N,<:Integer}}}) where N
    d = Dict()
    for genotype in locus
        genotype === missing && continue
        # sum up alleles
        for allele in genotype
            d[allele] = get!(d, allele, 0) + 1
        end
    end
    total = values(d) |> sum    # sum of all non-missing alleles
    [d[i] = d[i] / total for i in keys(d)]  # allele count / sum
    return d
end

"""
    allele_freq(locus::SubArray{<:Union{Missing,Tuple{Vararg}}})
Calculate allele counts for a single locus of a `PopObj` split by population
using `group()` with unequal ploidy across samples. Returns a `Dict` of
allele's and their frequencies.
"""
function allele_freq(locus::SubArray{<:Union{Missing,Tuple{Vararg}}})
    d = Dict()
    for genotype in locus
        genotype === missing && continue
        # sum up alleles
        for allele in genotype
            d[allele] = get!(d, allele, 0) + 1
        end
    end
    total = values(d) |> sum    # sum of all non-missing alleles
    [d[i] = d[i] / total for i in keys(d)]  # allele count / sum
    return d
end

"""
    allele_matrix(data::PopObj, allele::Int)
Convert a DataFrame of genotypes from a PopObj into a matrix of a single
allele. Use `allele` to specify which allele you want isolated. Intended
to be used to create separate allele matrices for comparing alleles 1 and
2, such as assessing heterozygosity.
"""
function allele_matrix(data::T, allele::Int) where T <: DataFrame
    geno_matrix = Matrix(data)
    map(geno_matrix) do i
        i === missing && return missing
        i[allele]
    end
end

"""
    geno_freq(locus::Vector{<:Union{Missing, NTuple{N,<:Integer}}}) where N
Calculate genotype frequencies of all loci in a `PopObj`. Returns a `Dict` of
genotypes and their frequencies.
"""
function geno_freq(locus::Vector{<:Union{Missing, NTuple{N,<:Integer}}}) where N
    d = Dict()
    # conditional testing if all genos are missing
    all(ismissing.(locus)) == true && return missing
    for genotype in skipmissing(locus)
        # sum up non-missing genotypes
        d[genotype] = get!(d, genotype, 0) + 1
    end
    total = values(d) |> sum    # sum of all non-missing genotypes
    [d[i] = d[i] / total for i in keys(d)] # genotype count/total
    return d
end

"""
    geno_freq(locus::Vector{<:Union{Missing, Tuple{Vararg}}})
Calculate genotype frequencies of all loci in a `PopObj` with unequal ploidy across
samples. Returns a `Dict` of genotypes and their frequencies.
"""
function geno_freq(locus::Vector{<:Union{Missing, Tuple{Vararg}}})
    d = Dict()
    # conditional testing if all genos are missing
    all(ismissing.(locus)) == true && return missing
    for genotype in skipmissing(locus)
        # sum up non-missing genotypes
        d[genotype] = get!(d, genotype, 0) + 1
    end
    total = values(d) |> sum    # sum of all non-missing genotypes
    [d[i] = d[i] / total for i in keys(d)] # genotype count/total
    return d
end

"""
    geno_freq(locus::SubArray{<:Union{Missing, NTuple{N,<:Integer}}}) where N
Calculate genotype frequencies of all loci in `PopObj` split
by population using `group()`. Returns a `Dict` of genotypes and their frequencies.
"""
function geno_freq(locus::SubArray{<:Union{Missing, NTuple{N,<:Integer}}}) where N
    d = Dict()
    # conditional testing if all genos are missing
    all(ismissing.(locus)) == true && return missing
    for genotype in skipmissing(locus)
        # sum up non-missing genotypes
        d[genotype] = get!(d, genotype, 0) + 1
    end
    total = values(d) |> sum    # sum of all non-missing genotypes
    [d[i] = d[i] / total for i in keys(d)] # genotype count/total
    return d
end

"""
    geno_freq(locus::SubArray{<:Union{Missing, Tuple{Vararg}}})
Calculate genotype frequencies of all loci in `PopObj` split by population
using `group()` with unequal ploidy across samples. Returns a `Dict` of
genotypes and their frequencies.
"""
function geno_freq(locus::SubArray{<:Union{Missing, Tuple{Vararg}}})
    d = Dict()
    # conditional testing if all genos are missing
    all(ismissing.(locus)) == true && return missing
    for genotype in skipmissing(locus)
        # sum up non-missing genotypes
        d[genotype] = get!(d, genotype, 0) + 1
    end
    total = values(d) |> sum    # sum of all non-missing genotypes
    [d[i] = d[i] / total for i in keys(d)] # genotype count/total
    return d
end

"""
    get_genotype(data::PopObj; sample::String, locus::String)
Return the genotype of a specific sample for specific locus in a `PopObj`.
This is a barebones variant to `isolate_genotypes` that is ~1000x faster and
suitable for development.

`get_genotype(nancycats, sample = "N115" , locus = "fca8")`

"""
function get_genotype(data::PopObj; sample::String, locus::String)
    idx = findfirst(i -> i == sample, data.samples.name)
    return getindex(data.loci[!, Symbol(locus)], idx)
end

"""
    get_sample_genotypes(data::PopObj, sample::String)
Return all the genotypes of a specific sample in a `PopObj`.
This is an extension for the internal function `get_genotype`.

`get_sample_genotypes(nancycats(), "N115")`
"""
function get_sample_genotypes(data::PopObj, sample::String)
    idx = findfirst(i -> i == sample, data.samples.name)
    map(i -> i[idx] , eachcol(data.loci))
end


text:
using CITIfile: rdciti
using Test

s,f = rdciti("hfss_terret_inner_perp.cit")
@testset "RI-style file" begin
    @test size(s) == (2,2,13)
    @test size(f) == (13,)
    @test s[1,1,1] ≈ -0.07517356 - 0.2106189im
    @test s[1,2,4] ≈ 0.9124461 - 0.3167198im
    @test s[end,end,end] ≈ -0.1570225 - 0.3573307im
end


text:
reload("src/master.jl")

#Setup Task Parameters
task         = H5Task("work.h5", "input/", "results/", H5Job[], Dict{Any,Any}())
task["cols"] = 18
script       = "example/example-job.jl"
func         = "job_function"

#Define Individual Jobs
job1         = H5Job(task, script, func, ["cols", "rows"])
job1["rows"] = 3
job2         = H5Job(task, script, func, ["cols", "rows"])
job2["rows"] = 4
job3         = H5Job(task, script, func, ["cols", "rows"])
job3["rows"] = 2
job4         = H5Job(task, script, func, ["cols", "rows"])
job4["rows"] = 1
job5         = H5Job(task, script, func, ["cols", "rows"])
job5["rows"] = 30
job6         = H5Job(task, script, func, ["cols", "rows"])
job6["rows"] = 3

#Run The Jobs
run(task, 8)

#View Some Results
read(job6)


text:
include("header.jl")
include("dataset.jl")
include("network.jl")
include("callback.jl")

# pyplot()
plot_loss(l_loss_train, l_loss_val; yscale = :log10)

function loss_neuralode_res(p)
    l_loss_exp = zeros(n_exp)
    for i_exp in 1:n_exp
        exp_data = l_exp_data[i_exp]
        pred = Array(pred_n_ode(p, i_exp, exp_data))
        masslist = sum(clamp.(@view(pred[1:end-1, :]), 0, Inf), dims = 1)'
        l_loss_exp[i_exp] = mae(masslist, @view(exp_data[1:length(masslist), 3]))
    end
    return l_loss_exp
end

println("results after pruning")
maxiters = 1e5
p_cutoff = 0.185
loss_epoch = zeros(Float64, n_exp);
for i_exp = 1:n_exp
    loss_epoch[i_exp] = loss_neuralode(p, i_exp)
end
loss_train = mean(loss_epoch[l_train])
loss_val = mean(loss_epoch[l_val])
@printf(
    "Loss train: %.2f val: %.2f p_cutoff: %.2e",
    log10(loss_train),
    log10(loss_val),
    p_cutoff
)

display_p(p)
for i_exp in randperm(n_exp)
    cbi(p, i_exp)
end

l_loss_exp = loss_neuralode_res(p)

l_exp_data = []
l_exp_info = zeros(Float64, length(l_exp), 3)
for (i_exp, value) in enumerate(l_exp)
    filename = string("exp_data/expdata_no", string(value), ".txt")
    exp_data = Float64.(load_exp(filename))
    push!(l_exp_data, exp_data)
    l_exp_info[i_exp, 1] = exp_data[1, 2] # initial temperature, K
end
l_exp_info[:, 2] = readdlm("exp_data/beta.txt")[l_exp]
l_exp_info[:, 3] = log.(readdlm("exp_data/ocen.txt")[l_exp] .+ llb)

l_loss_exp = loss_neuralode_res(p)

conc_end = zeros(n_exp, ns)

# Plot TGA data
list_plt = []
for i_exp = 1:14
    T0, beta, ocen = l_exp_info[i_exp, :]
    exp_data = l_exp_data[i_exp]
    sol = pred_n_ode(p, i_exp, exp_data)
    Tlist = similar(sol.t)
    T0, beta, ocen = l_exp_info[i_exp, :]
    for (i, t) in enumerate(sol.t)
        Tlist[i] = getsampletemp(t, T0, beta)
    end

    conc_end[i_exp, :] .= sol[:, end]

    if beta < 100
        plt = plot(
            Tlist,
            exp_data[:, 3],
            seriestype = :scatter,
            label = "Exp-$i_exp",
            framestyle = :box,
        )
        plot!(
            plt,
            Tlist,
            sum(sol[1:end-1, :], dims = 1)',
            lw = 2,
            legend = false,
        )
        xlabel!(plt, "Temperature [K]")

    else
        plt = plot(
            exp_data[:, 1] / 60.0,
            exp_data[:, 3],
            seriestype = :scatter,
            label = "Exp-$i_exp",
            framestyle = :box,
        )
        plot!(
            plt,
            sol.t / 60,
            sum(sol[1:end-1, :], dims = 1)',
            lw = 2,
            legend = false,
        )
        xlabel!(plt, "Time [min]")
    end
    exp_cond = @sprintf("T0:%.0f K", T0)
    if beta < 100.0
        exp_cond *= @sprintf("\nbeta:%.0f K/min", beta)
    end
    if exp(ocen) * 100.0 > 0.001
        exp_cond *= @sprintf("\n[O2]:%.2f", exp(ocen))
    else
        exp_cond *= @sprintf("\ninert")
    end

    ann_loc = [0.2, 0.3]
    if i_exp in [6, 7, 14]
        ann_loc = [0.7, 0.4]
    end
    annotate!(plt, xlims(plt)[1] + (xlims(plt)[end] - xlims(plt)[1]) * ann_loc[1],
                   ylims(plt)[1] + (ylims(plt)[end] - ylims(plt)[1]) * ann_loc[2],
                   text(exp_cond, 11))

    ylabel!(plt, "Mass [-] (No. $i_exp)")
    ylims!(plt, (0.0, 1.0))
    plot!(
        plt,
        xtickfontsize = 11,
        ytickfontsize = 11,
        xguidefontsize = 12,
        yguidefontsize = 12,
    )

    push!(list_plt, plt)
end
plt_all = plot(list_plt..., layout = (7, 2))
plot!(plt_all, size = (800, 1200))
png(plt_all, string(fig_path, "/TGA_mass_summary"))

# https://github.com/DENG-MIT/Biomass.jl/blob/main/backup/crnn_cellulose_ocen_test.jl
varnames = ["Cellu", "S2", "S3", "Vola"]
for i_exp in 1:n_exp
    T0, beta, ocen = l_exp_info[i_exp, :]
    exp_data = l_exp_data[i_exp]
    sol = pred_n_ode(p, i_exp, exp_data)
    Tlist = similar(sol.t)
    T0, beta, ocen = l_exp_info[i_exp, :]
    for (i, t) in enumerate(sol.t)
        Tlist[i] = getsampletemp(t, T0, beta)
    end
    value = l_exp[i_exp]
    list_plt = []
    scale_factor = 1 ./ maximum(sol[:, :], dims=2)
    scale_factor .= 1.0
    plt = plot(Tlist, clamp.(sol[1, :], 0, Inf), label = varnames[1])
    for i in 2:ns
        if scale_factor[i] > 1.1
            _label =  @sprintf("%s x %.2e", varnames[i], scale_factor[i])
        else
            _label = varnames[i]
        end
        plot!(plt, Tlist, clamp.(sol[i, :], 0, Inf) * scale_factor[i], label = _label)
    end
    xlabel!(plt, "Temperature [K]");
    ylabel!(plt, "Mass (-)");
    plot!(plt, size=(350, 350), legend=:topleft, framestyle=:box)
    plot!(
        plt,
        xtickfontsize = 11,
        ytickfontsize = 11,
        xguidefontsize = 12,
        yguidefontsize = 12,
    )
    png(plt, string(fig_path, "/pred_S_exp_$value"))
end

include("dataset.jl")
gr()

max_stars_repo_name: DENG-MIT/CRNN-Pyrolysis


text:
module LimnoSES

using Accessors
using Reexport
@reexport using Agents, DiffEqNoiseProcess
using StatsBase
using BlackBoxOptim
using DataInterpolations
using Parameters
import OrdinaryDiffEq
import Random

include("config.jl")

include("models/scheffer.jl")
include("models/martin.jl")
include("lake_state.jl")

include("initialise.jl")
include("evolve.jl")
include("reporting.jl")

include("decisions.jl")

end # module


text:
module CATKEVerticalDiffusivities

using Adapt
using KernelAbstractions: @kernel, @index

using Oceananigans.Architectures
using Oceananigans.Grids
using Oceananigans.Utils
using Oceananigans.Fields
using Oceananigans.BoundaryConditions: default_prognostic_field_boundary_condition
using Oceananigans.BoundaryConditions: BoundaryCondition, FieldBoundaryConditions, DiscreteBoundaryFunction, FluxBoundaryCondition
using Oceananigans.BuoyancyModels: ∂z_b, top_buoyancy_flux
using Oceananigans.Operators: ℑzᵃᵃᶜ

using Oceananigans.TurbulenceClosures:
    get_closure_ij,
    AbstractTurbulenceClosure,
    ExplicitTimeDiscretization,
    VerticallyImplicitTimeDiscretization

import Oceananigans.BoundaryConditions: getbc
import Oceananigans.Utils: with_tracers
import Oceananigans.TurbulenceClosures:
    add_closure_specific_boundary_conditions,
    calculate_diffusivities!,
    DiffusivityFields,
    z_viscosity,
    z_diffusivity,
    viscous_flux_ux,
    viscous_flux_uy,
    viscous_flux_uz,
    viscous_flux_vx,
    viscous_flux_vy,
    viscous_flux_vz,
    viscous_flux_wx,
    viscous_flux_wy,
    viscous_flux_wz,
    diffusive_flux_x,
    diffusive_flux_y,
    diffusive_flux_z

function hydrostatic_turbulent_kinetic_energy_tendency end

struct CATKEVerticalDiffusivity{TD, CD, CL, CQ} <: AbstractTurbulenceClosure{TD}
    Cᴰ :: CD
    mixing_length :: CL
    surface_TKE_flux :: CQ
end

function CATKEVerticalDiffusivity{TD}(Cᴰ:: CD,
                                      mixing_length :: CL,
                                      surface_TKE_flux :: CQ) where {TD, CD, CL, CQ}

    return CATKEVerticalDiffusivity{TD, CD, CL, CQ}(Cᴰ, mixing_length, surface_TKE_flux)
end

const CATKEVD = CATKEVerticalDiffusivity

# Support for "ManyIndependentColumnMode"
const CATKEVDArray = AbstractArray{<:CATKEVD}

function with_tracers(tracer_names, closure::Union{CATKEVD, CATKEVDArray})
    :e ∈ tracer_names || error("Tracers must contain :e to represent turbulent kinetic energy for `CATKEVerticalDiffusivity`.")
    return closure
end

include("mixing_length.jl")
include("surface_TKE_flux.jl")
include("turbulent_kinetic_energy_equation.jl")

for S in (:MixingLength, :SurfaceTKEFlux)
    @eval @inline convert_eltype(::Type{FT}, s::$S) where FT = $S{FT}(; Dict(p => getproperty(s, p) for p in propertynames(s))...)
    @eval @inline convert_eltype(::Type{FT}, s::$S{FT}) where FT = s
end

"""
    CATKEVerticalDiffusivity(FT=Float64;
                             Cᴰ = 2.91,
                             mixing_length = MixingLength{FT}(),
                             surface_TKE_flux = SurfaceTKEFlux{FT}(),
                             time_discretization::TD = VerticallyImplicitTimeDiscretization())

Returns the `CATKEVerticalDiffusivity` turbulence closure for vertical mixing by
small-scale ocean turbulence based on the prognostic evolution of subgrid
Turbulent Kinetic Energy (TKE).
"""
function CATKEVerticalDiffusivity(FT=Float64;
                                  Cᴰ = 2.91,
                                  mixing_length = MixingLength{FT}(),
                                  surface_TKE_flux = SurfaceTKEFlux{FT}(),
                                  warning = true,
                                  time_discretization::TD = VerticallyImplicitTimeDiscretization()) where TD

    if warning
        @warn "CATKEVerticalDiffusivity is an experimental turbulence closure that \n" *
              "is unvalidated and whose default parameters are not calibrated for \n" *
              "realistic ocean conditions or for use in a three-dimensional \n" *
              "simulation. Use with caution and report bugs and problems with physics \n" *
              "to https://github.com/CliMA/Oceananigans.jl/issues."
    end

    Cᴰ = convert(FT, Cᴰ)
    mixing_length = convert_eltype(FT, mixing_length)
    surface_TKE_flux = convert_eltype(FT, surface_TKE_flux)

    return CATKEVerticalDiffusivity{TD}(Cᴰ, mixing_length, surface_TKE_flux)
end

#####
##### Diffusivities and diffusivity fields utilities
#####

function DiffusivityFields(grid, tracer_names, bcs, closure::Union{CATKEVD, CATKEVDArray})

    default_diffusivity_bcs = (Kᵘ = FieldBoundaryConditions(grid, (Center, Center, Center)),
                               Kᶜ = FieldBoundaryConditions(grid, (Center, Center, Center)),
                               Kᵉ = FieldBoundaryConditions(grid, (Center, Center, Center)))

    bcs = merge(default_diffusivity_bcs, bcs)

    Kᵘ = CenterField(grid, boundary_conditions=bcs.Kᵘ)
    Kᶜ = CenterField(grid, boundary_conditions=bcs.Kᶜ)
    Kᵉ = CenterField(grid, boundary_conditions=bcs.Kᵉ)

    return (; Kᵘ, Kᶜ, Kᵉ)
end

function calculate_diffusivities!(diffusivities, closure::Union{CATKEVD, CATKEVDArray}, model)

    arch = model.architecture
    grid = model.grid
    velocities = model.velocities
    tracers = model.tracers
    buoyancy = model.buoyancy
    clock = model.clock
    top_tracer_bcs = NamedTuple(c => tracers[c].boundary_conditions.top for c in propertynames(tracers))

    event = launch!(arch, grid, :xyz,
                    calculate_CATKE_diffusivities!,
                    diffusivities, grid, closure, velocities, tracers, buoyancy, clock, top_tracer_bcs,
                    dependencies = device_event(arch))

    wait(device(arch), event)

    return nothing
end

@kernel function calculate_CATKE_diffusivities!(diffusivities, grid, closure::Union{CATKEVD, CATKEVDArray}, args...)
    i, j, k, = @index(Global, NTuple)

    # Ensure this works with "ensembles" of closures, in addition to ordinary single closures
    closure_ij = get_closure_ij(i, j, closure)

    @inbounds begin
        diffusivities.Kᵘ[i, j, k] = Kuᶜᶜᶜ(i, j, k, grid, closure_ij, args...)
        diffusivities.Kᶜ[i, j, k] = Kcᶜᶜᶜ(i, j, k, grid, closure_ij, args...)
        diffusivities.Kᵉ[i, j, k] = Keᶜᶜᶜ(i, j, k, grid, closure_ij, args...)
    end
end

@inline turbulent_velocity(i, j, k, grid, e) = @inbounds sqrt(max(zero(eltype(grid)), e[i, j, k]))

@inline function Kuᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, top_tracer_bcs)
    u★ = turbulent_velocity(i, j, k, grid, tracers.e)
    ℓu = momentum_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, top_tracer_bcs)
    return ℓu * u★
end

@inline function Kcᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, top_tracer_bcs)
    u★ = turbulent_velocity(i, j, k, grid, tracers.e)
    ℓc = tracer_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, top_tracer_bcs)
    return ℓc * u★
end

@inline function Keᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, top_tracer_bcs)
    u★ = turbulent_velocity(i, j, k, grid, tracers.e)
    ℓe = TKE_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, clock, top_tracer_bcs)
    return ℓe * u★
end

#####
##### Viscous flux, diffusive fluxes, plus shenanigans for diffusive fluxes of TKE (eg TKE "transport")
#####

# Special "index type" alternative to Val for dispatch
struct TKETracerIndex{N} end

@inline TKETracerIndex(N) = TKETracerIndex{N}()

@inline function viscous_flux_uz(i, j, k, grid, closure::CATKEVD, clock, velocities, diffusivities, args...)
    Ku = ℑxzᶠᵃᶠ(i, j, k, grid, diffusivities.Kᵘ)
    return - Ku * ∂zᵃᵃᶠ(i, j, k, grid, velocities.u)
end

@inline function viscous_flux_vz(i, j, k, grid, closure::CATKEVD, clock, velocities, diffusivities, args...)
    Kv = ℑyzᵃᶠᶠ(i, j, k, grid, diffusivities.Kᵘ)
    return - Kv * ∂zᵃᵃᶠ(i, j, k, grid, velocities.v)
end

@inline function viscous_flux_wz(i, j, k, grid, closure::CATKEVD, clock, velocities, diffusivities, args...)
    @inbounds Kw = diffusivities.Kᵘ[i, j, k]
    return - Kw * ∂zᵃᵃᶜ(i, j, k, grid, velocities.w)
end

@inline function diffusive_flux_z(i, j, k, grid, closure::CATKEVD, c, tracer_index, clock, diffusivities, args...)
    Kcᶜᶜᶠ = ℑzᵃᵃᶠ(i, j, k, grid, diffusivities.Kᶜ)
    return - Kcᶜᶜᶠ * ∂zᵃᵃᶠ(i, j, k, grid, c)
end

# Diffusive flux of TKE!
@inline function diffusive_flux_z(i, j, k, grid, closure::CATKEVD, e, ::TKETracerIndex, clock, diffusivities, args...)
    Keᶜᶜᶠ = ℑzᵃᵃᶠ(i, j, k, grid, diffusivities.Kᵉ)
    return - Keᶜᶜᶠ * ∂zᵃᵃᶠ(i, j, k, grid, e)
end

# "Translations" for diffusive transport by non-CATKEVD closures
@inline diffusive_flux_x(i, j, k, grid, closure, e, ::TKETracerIndex{N}, args...) where N = diffusive_flux_x(i, j, k, grid, closure, e, Val(N), args...)
@inline diffusive_flux_y(i, j, k, grid, closure, e, ::TKETracerIndex{N}, args...) where N = diffusive_flux_y(i, j, k, grid, closure, e, Val(N), args...)
@inline diffusive_flux_z(i, j, k, grid, closure, e, ::TKETracerIndex{N}, args...) where N = diffusive_flux_z(i, j, k, grid, closure, e, Val(N), args...)

# Shortcuts --- CATKEVD incurs no horizontal transport
@inline viscous_flux_ux(i, j, k, grid, ::CATKEVD, args...) = zero(eltype(grid))
@inline viscous_flux_uy(i, j, k, grid, ::CATKEVD, args...) = zero(eltype(grid))
@inline viscous_flux_vx(i, j, k, grid, ::CATKEVD, args...) = zero(eltype(grid))
@inline viscous_flux_vy(i, j, k, grid, ::CATKEVD, args...) = zero(eltype(grid))
@inline viscous_flux_wx(i, j, k, grid, ::CATKEVD, args...) = zero(eltype(grid))
@inline viscous_flux_wy(i, j, k, grid, ::CATKEVD, args...) = zero(eltype(grid))
@inline diffusive_flux_x(i, j, k, grid, ::CATKEVD, args...) = zero(eltype(grid))
@inline diffusive_flux_y(i, j, k, grid, ::CATKEVD, args...) = zero(eltype(grid))

# Disambiguate
@inline diffusive_flux_x(i, j, k, grid, ::CATKEVD, e, ::TKETracerIndex, args...) = zero(eltype(grid))
@inline diffusive_flux_y(i, j, k, grid, ::CATKEVD, e, ::TKETracerIndex, args...) = zero(eltype(grid))

#####
##### Support for VerticallyImplicitTimeDiscretization
#####

const VITD = VerticallyImplicitTimeDiscretization

@inline z_viscosity(closure::Union{CATKEVD, CATKEVDArray}, diffusivities, args...) = diffusivities.Kᵘ

@inline function z_diffusivity(closure::Union{CATKEVD, CATKEVDArray}, ::Val{tracer_index},
                               diffusivities, tracers, args...) where tracer_index

    tke_index = findfirst(name -> name === :e, keys(tracers))

    if tracer_index === tke_index
        return diffusivities.Kᵉ
    else
        return diffusivities.Kᶜ
    end
end

const VerticallyBoundedGrid{FT} = AbstractGrid{FT, <:Any, <:Any, <:Bounded}

@inline diffusive_flux_z(i, j, k, grid, ::VITD, closure::CATKEVD, args...) = zero(eltype(grid))
@inline viscous_flux_uz(i, j, k, grid, ::VITD, closure::CATKEVD, args...) = zero(eltype(grid))
@inline viscous_flux_vz(i, j, k, grid, ::VITD, closure::CATKEVD, args...) = zero(eltype(grid))

@inline function diffusive_flux_z(i, j, k, grid::VerticallyBoundedGrid{FT}, ::VITD, closure::CATKEVD, args...) where FT
    return ifelse(k == 1 || k == grid.Nz+1, # on boundaries, calculate fluxes explicitly
                  diffusive_flux_z(i, j, k, grid, ExplicitTimeDiscretization(), closure, args...),
                  zero(FT))
end

@inline function viscous_flux_uz(i, j, k, grid::VerticallyBoundedGrid{FT}, ::VITD, closure::CATKEVD, args...) where FT
    return ifelse(k == 1 || k == grid.Nz+1, # on boundaries, calculate fluxes explicitly
                  viscous_flux_uz(i, j, k, grid, ExplicitTimeDiscretization(), closure, args...),
                  zero(FT))
end

@inline function viscous_flux_vz(i, j, k, grid::VerticallyBoundedGrid{FT}, ::VITD, closure::CATKEVD, args...) where FT
    return ifelse(k == 1 || k == grid.Nz+1, # on boundaries, calculate fluxes explicitly
                  viscous_flux_vz(i, j, k, grid, ExplicitTimeDiscretization(), closure, args...),
                  zero(FT))
end

@inline function viscous_flux_wz(i, j, k, grid::VerticallyBoundedGrid{FT}, ::VITD, closure::CATKEVD, args...) where FT
    return ifelse(k == 1 || k == grid.Nz+1, # on boundaries, calculate fluxes explicitly
                  viscous_flux_wz(i, j, k, grid, ExplicitTimeDiscretization(), closure, args...),
                  zero(FT))
end

end


text:
using ReSHOP
using Pkg
using Compat.Test

pkgs = Pkg.installed()

@testset "Common basic tests" begin
	include("nl_convert.jl")
	include("nl_linearity.jl")
end

if pkgs["JuMP"] >= v"0.19"
	include("runtests_moi.jl")
	include("export_gms.jl")
	export_gms(mktempdir())
else
	include("runtests_mbp.jl")
end


text:
using DiffusionDefinition
using Test
using StaticArrays, Random
DD = DiffusionDefinition

@testset "DiffusionDefinition.jl" begin
    @diffusion_process TestDiffusion begin
        :dimensions
        process --> 3
        wiener --> 5

        :parameters
        _ --> (3, Float64)
        param --> Float64
        stem --> (2, Int64)
        (theta, alpha) --> (2, Int32)
        (beta, gamma) --> Float32
        (yota, zeta) --> (Float64, Int64)
        v --> SArray{Tuple{4},Float64,1,4}

        :conjugate
        phi(t, u) --> (
            (p1*u[2]*(u[2]-1),u[1],u[1]),
            (p3*u[5]*t,u[2],u[4]),
        )
        nonhypo(x) --> x
        num_non_hypo --> 3

        :additional
        statespace --> LowerBoundedStateSpace((1,2,3), (0.0, 1.0, 0.0))
        constdiff --> false
        linear --> false
        eltype --> Float64
    end a b c d

    P = TestDiffusion(1.0, 2.0, 3.0, 4.0, 5, 6, Int32(7), Int32(8),
                      Float32(9.0), Float32(10.0), 11.0, 12,
                      @SVector [1.0, 2.0, 3.0, 4.0])

    @test DD.nonhypo(P, :x) == :x
    @testset "satisfying bounds" begin
        @test !DD.bound_satisfied(P, [0.1, 2.0, 0.0])
        @test DD.bound_satisfied(P, [0.1, 2.0, 0.1])
        @test !DD.bound_satisfied(P, [0.1, 1.0, 0.1])
    end
    @test !DD.constdiff(P)
    @test typeof(P) <: DD.DiffusionProcess
    @test !(typeof(P) <: DD.LinearDiffusion)
    @test eltype(P) == Float64
end

@testset "State space restrictions" begin
    lb = LowerBoundedStateSpace((1,3), (-2.0, 3.0))
    ub = UpperBoundedStateSpace((1,3), (4.0, 10.0))
    b = BoundedStateSpace(lb, ub)
    nb = UnboundedStateSpace()

    a = [-1,0,1.0,4.0]
    v = [-1,0,1.0,20.0]
    @test DD.bound_satisfied(nb, a)
    @test !DD.bound_satisfied(lb, a)
    @test DD.bound_satisfied(ub, a)
    @test !DD.bound_satisfied(b, a)
    @test DD.bound_satisfied(nb, v)
    @test !DD.bound_satisfied(lb, v)
    @test DD.bound_satisfied(ub, v)
    @test !DD.bound_satisfied(b, v)
end

@testset "utility functions" begin
    @test DD.remove_curly(Vector{Float64}) == Array
    @test DD.get_curly(Vector{Float64}) == (:Float64, 1)
end

@testset "Buffers" begin
    d, m = 4,5
    el = Float64
    sb = DD.StandardEulerBuffer{el}(d, m)
    @test size(sb) == (d + d*m + m,)
    @test length(sb.b) == d
    @test size(sb.σ) == (d, m)
    @test length(sb.dW) == m
    @test eltype(sb) == el
    sb2 = similar(sb)
    @test sb2 == sb
    @test sb2 !== sb
    sb3 = similar(sb, SVector{3,Float64})
    @test size(sb3) == (d + d*m + m,)
    @test sb3.data == zeros(SVector{3,Float64}, d+d*m+m)

    @testset "indexing StandardEulerBuffer" begin
        sb[1] = 1.0
        sb[:] == sb.data
        sb[2:5] .= [2.0, 3.0, 4.0, 5.0]
        sb[[6,8]] .= [6.0, 8.0]
        sb[15] = 15.0
        sb[25] = 25.0

        @test sb.b == [1.0, 2.0, 3.0, 4.0]
        @test sb.σ == [
            5.0 0.0 0.0 0.0 0.0;
            6.0 0.0 0.0 0.0 0.0;
            0.0 0.0 15.0 0.0 0.0;
            8.0 0.0 0.0 0.0 0.0
        ]
        @test sb.dW == [25.0, 0.0, 0.0, 0.0, 0.0]
    end
    lb = DD.LinearDiffBuffer{el}(d, m)
    @test size(lb) == (d + d*m + d*d + m,)
    @test length(lb.b) == d
    @test size(lb.B) == (d,d)
    @test size(lb.σ) == (d, m)
    @test length(lb.dW) == m
end

@testset "Sampling trajectories" begin
    Random.seed!(3)
    N, dt, T = 30, 0.0001, 10.0
    tt = collect(0.0:dt:T)
    quad_var(x) = sum(abs2.(diff(x)))

    tr = trajectory(tt, rand(SVector{N,Float64},length(tt)))
    wr = Wiener()
    rand!(tr, wr)

    tr_mut = trajectory(tt, [rand(N) for _ in 1:length(tt)])
    wr_mut = Wiener(N)
    rand!(tr_mut, wr_mut)

    @test maximum([abs(quad_var(map(x->x[i], tr.x))/T-1.0) for i in 1:N])<0.05
    @test maximum([abs(quad_var(map(x->x[i], tr_mut.x))/T-1.0) for i in 1:N])<0.05

    # an example of a pre-defined diffusion
    @load_diffusion LotkaVolterraAux

    α, β, γ, δ, σ1, σ2 = 2.0/3.0, 4.0/3.0, 1.0, 1.0, 0.2, 0.3
    lv_aux = LotkaVolterraAux(
        α, β, γ, δ, σ1, σ2,
        0.0, 1.0, zero(DD.ℝ{2}), zero(DD.ℝ{2})
    )
    N, dt, T = 2, 0.0001, 1000.0
    tt = collect(0.0:dt:T)

    WW = trajectory(tt, rand(SVector{N,Float64},length(tt)))
    wr = Wiener()
    rand!(WW, wr)

    v0 = @SVector [1.0, 2.0]
    XX = trajectory(tt, rand(SVector{N,Float64},length(tt)))
    DD.solve!(XX, WW, lv_aux, v0)

    mean(x) = sum(x)/length(x)
    x1 = map(x->x[1], XX.x)
    x2 = map(x->x[2], XX.x)

    @test abs(mean(x1)-γ/δ) < 0.05
    @test abs(mean(x2)-α/β) < 0.05
    @test abs(quad_var(x1)/T-σ1^2) < 0.01
    @test abs(quad_var(x2)/T-σ2^2) < 0.01

    # and in-place simulation
    function DD.B!(buffer, t, P::LotkaVolterraAux)
        buffer.B[1,1] = 0.0
        buffer.B[1,2] = -P.β*P.γ/P.δ
        buffer.B[2,1] = P.α*P.δ/P.β
        buffer.B[2,2] = 0.0
    end

    function DD.β!(buffer, t, P::LotkaVolterraAux)
        buffer.b[1] = P.γ/P.δ*P.α
        buffer.b[2] = -P.α/P.β*P.γ
    end

    function DD.σ!(buffer, t, P::LotkaVolterraAux)
        buffer.σ[1,1] = P.σ1
        buffer.σ[1,2] = 0.0
        buffer.σ[2,1] = 0.0
        buffer.σ[2,2] = P.σ2
    end

    WW = trajectory(tt, [rand(N) for _ in 1:length(tt)])
    wr = Wiener(N)
    rand!(WW, wr)

    XX = trajectory(tt, [rand(N) for _ in 1:length(tt)])
    v0 = [1.0, 2.0]
    buffer = DD.LinearDiffBuffer{Float64}(DD.dimension(lv_aux)...)
    DD.solve!(DD.EulerMaruyama(), XX, WW, lv_aux, v0, buffer)

    x1 = map(x->x[1], XX.x)
    x2 = map(x->x[2], XX.x)
    @test abs(mean(x1)-γ/δ) < 0.05
    @test abs(mean(x2)-α/β) < 0.05
    @test abs(quad_var(x1)/T-σ1^2) < 0.01
    @test abs(quad_var(x2)/T-σ2^2) < 0.01
end


text:
using TestPkg2
using Test

@testset "TestPkg2.jl" begin
    # Write your tests here.
end


text:
using BinaryBuilder

name = "licensecheck"
version = v"0.3.1"

sources = [
    GitSource("https://github.com/google/licensecheck",
              "16aaea36649f556bae5a5ee972c247f58a0de1c4"),
    DirectorySource("./bundled")

]

# Bash recipe for building across all platforms
script = raw"""
install_license ${WORKSPACE}/srcdir/licensecheck/LICENSE
cd $WORKSPACE/srcdir/licensecheck/
mkdir clib
cp $WORKSPACE/srcdir/main.go clib/main.go
mkdir -p ${libdir}
CGO_ENABLED=1 go build -buildmode=c-shared -o ${libdir}/licensecheck.${dlext} clib/main.go
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("licensecheck", :licensecheck),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; compilers=[:c, :go])

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
# Forward/Backward

# {forward,backward}(a, A, L)
for f in (:forward, :backward)
    f!  = Symbol("$(f)!")    # forward!
    fl! = Symbol("$(f)log!") # forwardlog!

    @eval begin
        """
            $($f)(a, A, L; logl) -> (Vector, Float)

        Compute $($f) probabilities using samples likelihoods.
        See [Forward-backward algorithm](https://en.wikipedia.org/wiki/Forward–backward_algorithm).

        **Output**
        - `Vector{Float64}`: $($f) probabilities.
        - `Float64`: log-likelihood of the observed sequence.
        """
        function $(f)(a::AbstractVector, A::AbstractMatrix, L::AbstractMatrix; logl = false)
            m = Matrix{Float64}(undef, size(L))
            c = Vector{Float64}(undef, size(L)[1])
            if logl
                $(fl!)(m, c, a, A, L)
            else
                warn_logl(L)
                $(f!)(m, c, a, A, L)
            end
            m, sum(log.(c))
        end
    end
end

# {forward,backward}(hmm, observations)
for f in (:forward, :backward)
    f!  = Symbol("$(f)!")    # forward!
    fl! = Symbol("$(f)log!") # forwardlog!

    @eval begin
        """
            $($f)(hmm, observations; logl, robust) -> (Vector, Float)

        Compute $($f) probabilities of the `observations` given the `hmm` model.

        **Output**
        - `Vector{Float64}`: $($f) probabilities.
        - `Float64`: log-likelihood of the observed sequence.

        **Example**
        ```julia
        using Distributions, HMMBase
        hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0,1), Normal(10,1)])
        y = rand(hmm, 1000)
        probs, tot = $($f)(hmm, y)
        ```
        """
        function $(f)(hmm::AbstractHMM, observations; robust = false, kwargs...)
            L = likelihoods(hmm, observations; robust = robust, kwargs...)
            $(f)(hmm.a, hmm.A, L; kwargs...)
        end
    end
end

# Posteriors

"""
    posteriors(α, β) -> Vector

Compute posterior probabilities from `α` and `β`.

**Arguments**
- `α::AbstractVector`: forward probabilities.
- `β::AbstractVector`: backward probabilities.
"""
function posteriors(α::AbstractMatrix, β::AbstractMatrix)
    γ = Matrix{Float64}(undef, size(α))
    posteriors!(γ, α, β)
    γ
end

"""
    posteriors(a, A, L; logl) -> Vector

Compute posterior probabilities using samples likelihoods.
"""
function posteriors(a::AbstractVector, A::AbstractMatrix, L::AbstractMatrix; kwargs...)
    α, _ = forward(a, A, L; kwargs...)
    β, _ = backward(a, A, L; kwargs...)
    posteriors(α, β)
end

"""
    posteriors(hmm, observations; logl, robust) -> Vector

Compute posterior probabilities using samples likelihoods.

**Example**
```julia
using Distributions, HMMBase
hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0,1), Normal(10,1)])
y = rand(hmm, 1000)
γ = posteriors(hmm, y)
```
"""
function posteriors(hmm::AbstractHMM, observations; robust = false, kwargs...)
    L = likelihoods(hmm, observations; robust = robust, kwargs...)
    posteriors(hmm.a, hmm.A, L; kwargs...)
end

# Likelihood

"""
    loglikelihood(hmm, observations; logl, robust) -> Float64

Compute the log-likelihood of the observations under the model.
This is defined as the sum of the log of the normalization coefficients in the forward filter.

**Output**
- `Float64`: log-likelihood of the observations sequence under the model.

**Example**
```jldoctest
using Distributions, HMMBase
hmm = HMM([0.9 0.1; 0.1 0.9], [Normal(0,1), Normal(10,1)])
loglikelihood(hmm, [0.15, 0.10, 1.35])
# output
-4.588183811489616
```
"""
function loglikelihood(hmm::AbstractHMM, observations; logl = false, robust = false)
    forward(hmm, observations, logl = logl, robust = robust)[2]
end

max_stars_repo_name: Shu-Lea-Lai/module-MultiHMMBase