text:
module Test_psutils

using PeriodicSystems
using DescriptorSystems
using Symbolics
using Test
using LinearAlgebra
using LinearAlgebra: BlasInt
using ApproxFun
#using BenchmarkTools

println("Test_psutils")

@testset "test_psutils" begin

# MB03BD example
A1 = Matrix{Float64}(I,3,3); A2 = [   1.0   2.0   0.0; 4.0  -1.0   3.0; 0.0   3.0   1.0]; A3 = Matrix{Float64}(I,3,3);
E1 =  [2.0   0.0   1.0; 0.0  -2.0  -1.0; 0.0   0.0   3.0]; E2 = Matrix{Float64}(I,3,3);
E3 = [ 1.0   0.0   1.0; 0.0   4.0  -1.0; 0.0   0.0  -2.0];
ev = eigvals(inv(E1)*A2*inv(E3))

# using the SLICOT wrapper
A = reshape([E1 A2 E3],3,3,3);
KMAX = 3
NMAX = 3
LDA1 = NMAX
LDA2 = NMAX
LDQ1 = NMAX
LDQ2 = NMAX
LDWORK = KMAX + max( 2*NMAX, 8*KMAX )
LIWORK = 2*KMAX + NMAX
QIND = Array{BlasInt,1}(undef, KMAX)
S = [-1,1,-1];
Q = Array{Float64,3}(undef, LDQ1, LDQ2, KMAX)
ALPHAR = Array{Float64,1}(undef, NMAX)
ALPHAI = Array{Float64,1}(undef, NMAX)
BETA = Array{Float64,1}(undef, NMAX)
SCAL = Array{BlasInt,1}(undef, NMAX)
IWORK = Array{BlasInt,1}(undef, LIWORK)
DWORK = Array{Float64,1}(undef, LDWORK)

mb03bd!('T','C','I',QIND,3,3,2,1,3,S,A,Q,ALPHAR, ALPHAI, BETA, SCAL, LIWORK, LDWORK)

poles = (ALPHAR+im*ALPHAI) ./ BETA .* (2. .^SCAL)

@test sort(real(poles)) ≈ sort(real(ev)) &&
      sort(imag(poles)) ≈ sort(imag(ev))

A = reshape([A2 inv(E3) inv(E1)],3,3,3);
mb03wd!('S','I',3,3,1,3,1,3,A,Q,ALPHAR, ALPHAI, LDWORK)

poles = (ALPHAR+im*ALPHAI)

@test sort(real(poles)) ≈ sort(real(ev)) &&
      sort(imag(poles)) ≈ sort(imag(ev))

# MB03VD example
A1 = [1.5 -.7 3.5 -.7; 1.  0.  2.  3.; 1.5 -.7 2.5 -.3; 1.  0.  2.  1.];
A2 = [1.5 -.7 3.5 -.7; 1.  0.  2.  3.; 1.5 -.7 2.5 -.3; 1.  0.  2.  1.];
n = 4; K = 2;
A = reshape([A1 A2], n, n, K);
@time H, Z, ihess = phess(A; rev = false);
@test check_psim(A,Z,H; rev = false) && istriu(H[:,:,ihess],-1) && ihess == 1

@time H1, Z1, ihess1 = phess(A; rev = false, withZ = false);
@test H[1] == H1[1] && H[2] == H1[2]

@time H, Z, ihess = phess(A, hind = 2, rev = false);
@test check_psim(A,Z,H; rev = false) && istriu(H[:,:,ihess],-1) && ihess == 2

@time H1, Z1, ihess1 = phess(A, hind = 2, rev = false, withZ = false);
@test H[1] == H1[1] && H[2] == H1[2]

@time H, Z, ihess = phess(A; rev = true);
@test check_psim(A,Z,H; rev = true) && istriu(H[:,:,ihess],-1) && ihess == 1

@time H1, Z1, ihess1 = phess(A; rev = true, withZ = false);
@test H[1] == H1[1] && H[2] == H1[2]

@time H, Z, ihess = phess(A; hind = 2, rev = true);
@test check_psim(A,Z,H; rev = true) && istriu(H[:,:,ihess],-1) && ihess == 2

@time H1, Z1, ihess1 = phess(A; hind = 2, rev = true, withZ = false);
@test H[1] == H1[1] && H[2] == H1[2]

ev = eigvals(A[:,:,1]*A[:,:,2])
@time S, Z, eigs, ischur, α, γ = pschur(A; rev = false);
@test check_psim(A,Z,S; rev = false) && istriu(S[:,:,ischur],-1) && ischur == 1
@test sort(real(eigs)) ≈ sort(real(ev)) &&
      sort(imag(eigs)) ≈ sort(imag(ev))

@time S1, Z1, eigs1, ischur1, α1, γ1 = pschur(A; rev = false, withZ = false);
@test S[1] == S1[1] && S[2] == S1[2] && eigs == eigs1

@time S, Z, eigs, ischur, α, γ = pschur(A, sind = 2, rev = false);
@test check_psim(A,Z,S; rev = false) && istriu(S[:,:,ischur],-1) && ischur == 2
@test sort(real(eigs)) ≈ sort(real(ev)) &&
      sort(imag(eigs)) ≈ sort(imag(ev))

ev = eigvals(A[:,:,2]*A[:,:,1])
@time S, Z, eigs, ischur, α, γ = pschur(A; rev = true);
@test check_psim(A,Z,S; rev = true) && istriu(S[:,:,ischur],-1) && ischur == 1
@test sort(real(eigs)) ≈ sort(real(ev)) &&
      sort(imag(eigs)) ≈ sort(imag(ev))

@time S, Z, eigs, ischur, α, γ = pschur(A; sind = 2, rev = true);
@test check_psim(A,Z,S; rev = true) && istriu(S[:,:,ischur],-1) && ischur == 2
@test sort(real(eigs)) ≈ sort(real(ev)) &&
      sort(imag(eigs)) ≈ sort(imag(ev))

A1 = [1.5 -.7 3.5 -.7; 1.  0.  2.  3.; 1.5 -.7 2.5 -.3; 1.  0.  2.  1.];
n = 4; K = 19; hind = 10; rev = false;
A = reshape(hcat([rand()*A1 for i in 1:K]...), n, n, K);
@time H, Z, ihess = phess(A; hind, rev);
@test check_psim(A,Z,H; rev) && istriu(H[:,:,ihess],-1) && ihess == hind

@time H, Z, ihess = phess1(A; hind, rev);
@test check_psim(A,Z,H; rev) && istriu(H[:,:,ihess],-1) && ihess == hind

@time S, Z, eigs, ischur, = pschur(A; sind = hind, rev);
@test check_psim(A,Z,S; rev)

@time S, Z, eigs, ischur, = pschur1(A; sind = hind, rev);
@test check_psim(A,Z,S; rev)

# this would fail
# n =4; K = 19; hind = 10; rev = false;
# A = reshape(hcat([rand()*A1 for i in 1:K]...), n, n, K);
# @time S, Z, eigs, ischur, = pschur2(A; sind = hind, rev);
# @test check_psim(A,Z,S; rev)

n =4; K = 15; hind = 10; rev = false;
A = reshape(hcat([rand()*A1 for i in 1:K]...), n, n, K);
@time S, Z, eigs, ischur, = pschur2(A; sind = hind, rev);
@test check_psim(A,Z,S; rev)

hind = 14; rev = true;
@time H, Z, ihess = phess(A; hind, rev);
@test check_psim(A,Z,H; rev) && istriu(H[:,:,ihess],-1) && ihess == hind

@time H, Z, ihess = phess1(A; hind, rev);
@test check_psim(A,Z,H; rev) && istriu(H[:,:,ihess],-1) && ihess == hind

@time S, Z, eigs, ischur, = pschur(A; sind = hind, rev);
@test check_psim(A,Z,S; rev)

@time S, Z, eigs, ischur, = pschur1(A; sind = hind, rev);
@test check_psim(A,Z,S; rev)

# this would fail
# n =4; K = 19; hind = 10; rev = false;
# A = reshape(hcat([rand()*A1 for i in 1:K]...), n, n, K);
# @time S, Z, eigs, ischur, = pschur2(A; sind = hind, rev);
# @test check_psim(A,Z,S; rev)

n =4; K = 15; hind = 10; rev = false;
A = reshape(hcat([rand()*A1 for i in 1:K]...), n, n, K);
@time S, Z, eigs, ischur, = pschur2(A; sind = hind, rev);
@test check_psim(A,Z,S; rev)

A1 = [1.5 -.7 3.5 -.7; 1.  0.  2.  3.; 1.5 -.7 2.5 -.3];
A2 = [1.5 -.7 3.5; 1.  0.  2.; 1.5 -.7 2.5 ; 1.  0.  2.];
Ai=[A1,A2];
A = [Ai[j] for j = 1:2,  i in 1:9][:]

rev = false;
@time S, Z, eigs, ischur, = pschur(A; rev);
@test check_psim(A,Z,S; rev)

rev = true;
@time S, Z, eigs, ischur, = pschur(A; rev);
@test check_psim(A,Z,S; rev)

# @time S1, Z1, eigs1, ischur1 = pschurw(A; rev = false);
# @test check_psim(A,Z1,S1; rev = false)

A1 = [1.5 -.7 3.5 -.7; 1.  0.  2.  3.; 1.5 -.7 2.5 -.3; 1.  0.  2.  1.];
ev1 = eigvals(A1)
n = 4; K = 18;

k = 1; N = div(K,k)
Ap = PeriodicArray(reshape(hcat([A1 for i in 1:N]...), n, n, N), K; nperiod = k)
@time ev = pseig(Ap)
@test sort(real(ev)) ≈ sort(real(ev1.^K)) && sort(imag(ev)) ≈ sort(imag(ev1.^K))

k = 9; N = div(K,k)
Ap = PeriodicArray(reshape(hcat([A1 for i in 1:N]...), n, n, N),K; nperiod = k)
@time ev = pseig(Ap)
@test sort(real(ev)) ≈ sort(real(ev1.^K)) && sort(imag(ev)) ≈ sort(imag(ev1.^K))

k = 18; N = div(K,k)
Ap = PeriodicArray(reshape(hcat([A1 for i in 1:N]...), n, n, N),K; nperiod = k)
@time ev = pseig(Ap)
@test sort(real(ev)) ≈ sort(real(ev1.^K)) && sort(imag(ev)) ≈ sort(imag(ev1.^K))

# modified MB03VD example
A1 = [1.5 -.7 3.5 -.7; 1.  0.  2.  3.; 1.5 -.7 2.5 -.3; 1.  0.  2.  1.];
A2 = [1.5 -.7 3.5 -.7; 1.  0.  2.  3.; 1.5 -.7 2.5 -.3; 1.  0.  2.  1.] .+ 1;
A3 = [1.5 -.7 3.5 -.7; 1.  0.  2.  3.; 1.5 -.7 2.5 -.3; 1.  0.  2.  1.] .- 1;

n = 4; K = 3; sind = 2;
A = reshape([A1 A2 A3], n, n, K);
#AH, Z = phess(A; h);
@time S, Z, eigs, ischur = pschur(A; sind);
@test check_psim(A, Z, S) && istriu(S[:,:,ischur],-1)

ev = eigvals(A[:,:,1]*A[:,:,2]*A[:,:,3])
@test sort(real(eigs)) ≈ sort(real(ev)) &&
      sort(imag(eigs)) ≈ sort(imag(ev))

i1 = reverse(1:K)
@time a, e = psreduc_reg(view(A,:,:,i1)); eigs = eigvals(a,e);
@test sort(real(eigs)) ≈ sort(real(ev)) &&
      sort(imag(eigs)) ≈ sort(imag(ev))

# modified MB03VD example
A1 = [1.5 -.7 3.5 -.7; 1.  0.  2.  3.; 1.5 -.7 2.5 -.3; 1.  0.  2.  1.];
A2 = [1.5 -.7 3.5 -.7; 1.  0.  2.  3.; 1.5 -.7 2.5 -.3; 1.  0.  2.  1.] .+ 1;
A3 = [1.5 -.7 3.5 -.7; 1.  0.  2.  3.; 1.5 -.7 2.5 -.3; 1.  0.  2.  1.] .- 1;
n = 4; K = 3;
A = reshape([A1 A2 A3], n, n, K);
M = PeriodicArray(A,K);
MA = PeriodicMatrix([A1,A2,A3],K);

@time eigs = pseig(A); ev = pseig(A,fast = true);
@test sort(real(eigs)) ≈ sort(real(ev)) &&
      isapprox(sort(imag(eigs)),sort(imag(ev)), atol = 1.e-7)
@test sort(eigs, by = real) ≈ sort(ev, by = real)
@test sort(eigs, by = abs) ≈ sort(ev, by = abs)
@test sort(psceig(M),by = real) ≈ sort(psceig(MA),by = real)

@time eigs = pseig(A); ev = pseig(M,fast = true);
@test sort(real(eigs)) ≈ sort(real(ev)) &&
      isapprox(sort(imag(eigs)),sort(imag(ev)), atol = 1.e-7)
@test sort(eigs, by = real) ≈ sort(ev, by = real)
@test sort(eigs, by = abs) ≈ sort(ev, by = abs)

@time AH, Z, ihess = phess(A; rev = false);
@test check_psim(A,Z,AH; rev = false) && istriu(AH[:,:,ihess],-1) && ihess == 1

@time AH, Z, ihess = phess(A, hind = 2, rev = false);
@test check_psim(A,Z,AH; rev = false) && istriu(AH[:,:,ihess],-1) && ihess == 2

@time AH, Z, ihess = phess(A; rev = true);
@test check_psim(A,Z,AH; rev = true) && istriu(AH[:,:,ihess],-1) && ihess == 1

@time AH, Z, ihess = phess(A; hind = 2, rev = true);
@test check_psim(A,Z,AH; rev = true) && istriu(AH[:,:,ihess],-1) && ihess == 2

ev = eigvals(A[:,:,1]*A[:,:,2]*A[:,:,3])
@time S, Z, eigs, ischur, α, γ = pschur(A; rev = false);
@test check_psim(A,Z,S; rev = false) && istriu(S[:,:,ischur],-1) && ischur == 1
@test sort(real(eigs)) ≈ sort(real(ev)) &&
      sort(imag(eigs)) ≈ sort(imag(ev))

@time S, Z, eigs, ischur, α, γ = pschur(A, sind = 2, rev = false);
@test check_psim(A,Z,S; rev = false) && istriu(S[:,:,ischur],-1) && ischur == 2
@test sort(real(eigs)) ≈ sort(real(ev)) &&
      sort(imag(eigs)) ≈ sort(imag(ev))

ev = eigvals(A[:,:,3]*A[:,:,2]*A[:,:,1])
@time S, Z, eigs, ischur, α, γ = pschur(A; rev = true);
@test check_psim(A,Z,S; rev = true) && istriu(S[:,:,ischur],-1) && ischur == 1
@test sort(real(eigs)) ≈ sort(real(ev)) &&
      sort(imag(eigs)) ≈ sort(imag(ev))

@time S, Z, eigs, ischur, α, γ = pschur(A; sind = 2, rev = true);
@test check_psim(A,Z,S; rev = true) && istriu(S[:,:,ischur],-1) && ischur == 2
@test sort(real(eigs)) ≈ sort(real(ev)) &&
      sort(imag(eigs)) ≈ sort(imag(ev))

ev = eigvals(A[:,:,3]*A[:,:,2]*A[:,:,1])
a, e = psreduc_reg(A); eigs = eigvals(a,e);
@test sort(real(eigs)) ≈ sort(real(ev)) &&
      sort(imag(eigs)) ≈ sort(imag(ev))

A = [A1, A2, A3];
@time eigs = pseig(A); ev = pseig(A,fast = true);
@test sort(real(eigs)) ≈ sort(real(ev)) &&
      isapprox(sort(imag(eigs)),sort(imag(ev)), atol = 1.e-7)

@time eigs = pseig(A; rev = false); ev = pseig(A; fast = true, rev = false);
@test sort(real(eigs)) ≈ sort(real(ev)) &&
      isapprox(sort(imag(eigs)),sort(imag(ev)), atol = 1.e-7)

@time eigs = pseig(A, 2; rev = false);
@time ev = pseig(A, 2; fast = true, rev = false);
@test sort(real(eigs)) ≈ sort(real(ev)) &&
      isapprox(sort(imag(eigs)),sort(imag(ev)), atol = 1.e-7)

A3 = rand(5,4); A2 = rand(4,3); A1 = rand(3,5);
#A = [A1, A2, A3];
K = 3;
A = PeriodicMatrix([A1, A2, A3],K);
@time eigs = pseig(A); ev = pseig(A,fast = true);
@test sort(real(eigs)) ≈ sort(real(ev)) &&
      isapprox(sort(imag(eigs)),sort(imag(ev)), atol = 1.e-7)
@test sort(eigs, by = real) ≈ sort(ev, by = real)
@test sort(eigs, by = abs) ≈ sort(ev, by = abs)

ev = eigvals(prod(A.M[K:-1:1])); nmin = minimum(size.(A.M,1))
@time S, Z, eigs, ischur, α, γ = pschur(A.M; rev = true);
@test check_psim(A.M,Z,S; rev = true) && istriu(S[ischur][1:nmin,1:nmin],-1) && eigs == α.*γ
@test isapprox(sort(real(eigs)),sort(real(ev)), atol = 1.e-7) &&
      isapprox(sort(imag(eigs)),sort(imag(ev)), atol = 1.e-7)

A = [A3, A2, A1];
ev = eigvals(prod(A)); nmin = minimum(size.(A,1))

@time S, Z, eigs, ischur, α, γ = pschur(A; rev = false);
@test check_psim(A,Z,S; rev = false) && istriu(S[ischur][1:nmin,1:nmin],-1) && eigs == α.*γ
@test sort(real(eigs)) ≈ sort(real(ev)) &&
      isapprox(sort(imag(eigs)),sort(imag(ev)), atol = 1.e-7)

# periodic time-varying matrices from IFAC2005 paper
@variables t
A = PeriodicSymbolicMatrix([cos(t) 1; 1 1-sin(t)],2*pi);
println("A(t) = $(A.F)")
B = HarmonicArray([0;1],[[1;0]],[[1;-1]],2*pi);
println("B(t) = $(convert(PeriodicSymbolicMatrix,B).F)")
C = PeriodicFunctionMatrix(t-> [sin(2*t)+cos(2*t) 1],pi);
println("C(t) = $(C.f(t))")
D = PeriodicFunctionMatrix(zeros(1,1),pi)
K = 100
Δ = 2*pi/K
ts = (0:K-1)*Δ

psys = ps(A,B,C,D);
psys = ps(A,B,C);
ps(A);
ps(B);
ps(C);
ps(D);
psys = ps(A,B,C,D,4*pi);
psys = ps(A,B,C,4*pi);
psys = ps(HarmonicArray,A,B,C,D,4*pi);
psys = ps(PeriodicSymbolicMatrix, A,B,C,4*pi);

psys = ps(rss(4,3,2),10)
@test islti(psys)
psys=ps(rand(2,2),rand(2),rand(3,2),rand(3),10)
psys=ps(rand(2,2),rand(2),rand(3,2),rand(3),10;Ts = 1)
psys=ps(rand(2,2),rand(2),rand(3,2),10;Ts = 1)

psys = PeriodicStateSpace(convert(PeriodicFunctionMatrix,A),
                          convert(PeriodicFunctionMatrix,B),
                          convert(PeriodicFunctionMatrix,C),
                          convert(PeriodicFunctionMatrix,D));
print(psys);
convert(PeriodicStateSpace{PeriodicSymbolicMatrix},psys);
convert(PeriodicStateSpace{HarmonicArray},psys);
convert(PeriodicStateSpace{PeriodicTimeSeriesMatrix},psys);
convert(PeriodicStateSpace{FourierFunctionMatrix},psys);

psys = PeriodicStateSpace(convert(PeriodicFunctionMatrix{:c,BigFloat},A),
                          convert(PeriodicFunctionMatrix{:c,BigFloat},B),
                          convert(PeriodicFunctionMatrix{:c,BigFloat},C),
                          convert(PeriodicFunctionMatrix{:c,BigFloat},D));

psys = ps(convert(PeriodicFunctionMatrix,A),
          convert(PeriodicFunctionMatrix,B),
          convert(PeriodicFunctionMatrix,C),
          convert(PeriodicFunctionMatrix,D));

psys = ps(convert(PeriodicFunctionMatrix,A),
          convert(PeriodicFunctionMatrix,B),
          convert(PeriodicFunctionMatrix,C));

psys = PeriodicStateSpace(convert(PeriodicSymbolicMatrix,A),
                          convert(PeriodicSymbolicMatrix,B),
                          convert(PeriodicSymbolicMatrix,C),
                          convert(PeriodicSymbolicMatrix,D));
print(psys);
convert(PeriodicStateSpace{PeriodicFunctionMatrix},psys);
convert(PeriodicStateSpace{HarmonicArray},psys);
convert(PeriodicStateSpace{PeriodicTimeSeriesMatrix},psys);
convert(PeriodicStateSpace{FourierFunctionMatrix},psys);

psys = ps(convert(PeriodicSymbolicMatrix,A),
          convert(PeriodicSymbolicMatrix,B),
          convert(PeriodicSymbolicMatrix,C),
          convert(PeriodicSymbolicMatrix,D));

psys = ps(convert(PeriodicSymbolicMatrix,A),
          convert(PeriodicSymbolicMatrix,B),
          convert(PeriodicSymbolicMatrix,C));

psys = PeriodicStateSpace(convert(HarmonicArray,A),
                          convert(HarmonicArray,B),
                          convert(HarmonicArray,C),
                          convert(HarmonicArray,D));
print(psys);
convert(PeriodicStateSpace{PeriodicSymbolicMatrix},psys);
convert(PeriodicStateSpace{PeriodicTimeSeriesMatrix},psys);
convert(PeriodicStateSpace{PeriodicFunctionMatrix},psys);
convert(PeriodicStateSpace{FourierFunctionMatrix},psys);

psys = ps(convert(HarmonicArray,A),
          convert(HarmonicArray,B),
          convert(HarmonicArray,C),
          convert(HarmonicArray,D));

psys = ps(convert(HarmonicArray,A),
          convert(HarmonicArray,B),
          convert(HarmonicArray,C));

At = PeriodicTimeSeriesMatrix(convert(PeriodicFunctionMatrix,A).f.(ts),A.period);
Bt = PeriodicTimeSeriesMatrix(convert(PeriodicFunctionMatrix,B).f.(ts),B.period);
Ct = PeriodicTimeSeriesMatrix(convert(PeriodicFunctionMatrix,C).f.(ts),C.period);
Dt = PeriodicTimeSeriesMatrix(convert(PeriodicFunctionMatrix,D).f.(ts),D.period);
psys = PeriodicStateSpace(At, Bt, Ct, Dt);
print(psys);
convert(PeriodicStateSpace{PeriodicFunctionMatrix},psys);
convert(PeriodicStateSpace{PeriodicSymbolicMatrix},psys);
convert(PeriodicStateSpace{HarmonicArray},psys);
convert(PeriodicStateSpace{FourierFunctionMatrix},psys);

psys = ps(At, Bt, Ct, Dt);
psys = ps(At, Bt, Ct);

psys = PeriodicStateSpace(convert(FourierFunctionMatrix,A),
                          convert(FourierFunctionMatrix,B),
                          convert(FourierFunctionMatrix,C),
                          convert(FourierFunctionMatrix,D));
print(psys);
convert(PeriodicStateSpace{PeriodicFunctionMatrix},psys);
convert(PeriodicStateSpace{HarmonicArray},psys);
convert(PeriodicStateSpace{PeriodicTimeSeriesMatrix},psys);
convert(PeriodicStateSpace{PeriodicSymbolicMatrix},psys);

psys = ps(convert(FourierFunctionMatrix,A),
          convert(FourierFunctionMatrix,B),
          convert(FourierFunctionMatrix,C),
          convert(FourierFunctionMatrix,D));
psys = ps(convert(FourierFunctionMatrix,A),
          convert(FourierFunctionMatrix,B),
          convert(FourierFunctionMatrix,C));

psys = ps(FourierFunctionMatrix,A,B,C,D);

# #new tests
# AF = ffm2hr(psys.A)
# convert(PeriodicSymbolicMatrix,AF)

# constant dimensions
Ad = PeriodicMatrix([[1. 0; 0 0], [1 1;1 1], [0 1; 1 0]], 6, nperiod = 2);
Bd = PeriodicMatrix( [[ 1; 0 ], [ 1; 1]] ,2);
Cd = PeriodicMatrix( [[ 1 1 ], [ 1 0]] ,2);
Dd = PeriodicMatrix( [[ 1 ]], 1);
psys = PeriodicStateSpace(Ad,Bd,Cd,Dd);
convert(PeriodicStateSpace{PeriodicArray},psys);
convert(PeriodicStateSpace{PeriodicTimeSeriesMatrix},psys);

psys = ps(Ad,Bd,Cd,Dd);
psys = ps(Ad,Bd,Cd);

Ad = PeriodicMatrix([[1. 0], [1;1]],2);
Bd = PeriodicMatrix( [[ 1 ], [ 1; 1]] ,2);
Cd = PeriodicMatrix( [[ 1 1 ], [ 1 ]] ,2);
Dd = PeriodicMatrix( [[ 1 ]], 1);
psys = PeriodicStateSpace(Ad,Bd,Cd,Dd);
print(psys);
convert(PeriodicStateSpace{PeriodicArray},psys);
convert(PeriodicStateSpace{PeriodicTimeSeriesMatrix},psys);

psys = ps(Ad,Bd,Cd,Dd);
psys = ps(Ad,Bd,Cd);

Ad = PeriodicArray(rand(Float32,2,2,10),10);
Bd = PeriodicArray(rand(2,1,2),2);
Cd = PeriodicArray(rand(1,2,3),3);
Dd = PeriodicArray(rand(1,1,1), 1);
psys = PeriodicStateSpace(Ad,Bd,Cd,Dd);
print(psys);
convert(PeriodicStateSpace{PeriodicMatrix},psys)
convert(PeriodicStateSpace{PeriodicTimeSeriesMatrix},psys)

psys = ps(Ad,Bd,Cd,Dd);
psys = ps(Ad,Bd,Cd);
psys = ps(convert(PeriodicMatrix,Ad), Bd, convert(PeriodicMatrix,Cd), Dd);
psys = ps(convert(PeriodicMatrix,Ad), Bd, convert(PeriodicMatrix,Cd));

# symbolic periodic
@variables t
A = [cos(t) 1; 1 1-sin(t)];
B = [cos(t)+sin(t); 1-sin(t)];
C = [sin(t)+cos(2*t) 1];
Ap = PeriodicSymbolicMatrix(A,2*pi);
Bp = PeriodicSymbolicMatrix(B,2*pi);
Cp = PeriodicSymbolicMatrix(C,2*pi);

# functional expressions
tA(t::Real) = [cos(t) 1; 1 1-sin(t)];
tB(t::Real) = [cos(t)+sin(t); 1-sin(t)];
tC(t::Real) = [sin(t)+cos(2*t) 1];
# store snapshots as 3d arrays
N = 200;
tg = collect((0:N-1)*2*pi/N);
time = (0:N-1)*2*pi/N;
# At = reshape(hcat(A.(t)...),2,2,N);
# Bt = reshape(hcat(B.(t)...),2,1,N);
# Ct = reshape(hcat(C.(t)...),1,2,N);

# time series expressions
At = PeriodicTimeSeriesMatrix(tA.(time),2*pi);
Bt = PeriodicTimeSeriesMatrix(tB.(time),2*pi);
Ct = PeriodicTimeSeriesMatrix(tC.(time),2*pi);

# harmonic expressions
@time Ahr = ts2hr(At);
@test Ahr.values[:,:,1] ≈ [0. 1; 1 1] && Ahr.values[:,:,2] ≈ [1. 0; 0 -im]
@time Bhr = ts2hr(Bt);
@test Bhr.values[:,:,1] ≈ [0.; 1] && Bhr.values[:,:,2] ≈ [1.0+im; -im]
@time Chr = ts2hr(Ct);
@test Chr.values[:,:,1] ≈ [0. 1] && Chr.values[:,:,2] ≈ [im 0] && Chr.values[:,:,3] ≈ [1 0]

#@time Affm = ts2ffm(At);

nperiod = 24
time1 = (0:N-1)*2*pi*nperiod/N;
At1 = PeriodicTimeSeriesMatrix(tA.(time1),2*pi*nperiod);
Ahr1 = ts2hr(At1);
@test convert(PeriodicFunctionMatrix,Ahr1).f(1) ≈ convert(PeriodicFunctionMatrix,Ahr).f(1)

@test iszero(hr2psm(Ahr,1:1) + hr2psm(Ahr,0:0) - hr2psm(Ahr))
@test iszero(hr2psm(Ahr1,1:1) + hr2psm(Ahr1,0:0) - hr2psm(Ahr1))

# harmonic vs. symbolic
@test norm(substitute.(convert(PeriodicSymbolicMatrix,Ahr).F - A, (Dict(t => rand()),))) < 1e-15
@test norm(substitute.(convert(PeriodicSymbolicMatrix,Bhr).F - B, (Dict(t => rand()),))) < 1e-15
@test norm(substitute.(convert(PeriodicSymbolicMatrix,Chr).F - C, (Dict(t => rand()),))) < 1e-15

# harmonic vs. time series
@test all(norm.(tvmeval(Ahr,tg).-At.values) .< 1.e-7)
@test all(norm.(tvmeval(Bhr,tg).-Bt.values) .< 1.e-7)
@test all(norm.(tvmeval(Chr,tg).-Ct.values) .< 1.e-7)

# check time values on the grid
for method in ("constant", "linear", "quadratic", "cubic")
      @test all(norm.(tvmeval(At,tg; method).-At.values) .< 1.e-7)
      @test all(norm.(tvmeval(Bt,tg; method).-Bt.values) .< 1.e-7)
      @test all(norm.(tvmeval(Ct,tg; method).-Ct.values) .< 1.e-7)
end
# check interpolated values: time series vs. harmonic
tt = rand(10)*2pi;
for method in ("linear", "quadratic", "cubic")
    @test all(norm.(tvmeval(At,tt; method).-tvmeval(Ahr,tt; exact = true)) .< 1.e-3)
    @test all(norm.(tvmeval(Bt,tt; method).-tvmeval(Bhr,tt; exact = false)) .< 1.e-3)
    @test all(norm.(tvmeval(Ct,tt; method).-tvmeval(Chr,tt; exact = true)) .< 1.e-3)
end

# check conversion to function form
Amat = convert(PeriodicFunctionMatrix,Ahr);
@test all(norm.(tvmeval(At,tt; method = "linear").-tvmeval(Ahr,tt)) .< 1.e-3)
@test iszero(convert(PeriodicSymbolicMatrix,Amat).F-Ap.F)
@test isperiodic(Amat)

Amat =  convert(PeriodicFunctionMatrix,At);
@test all(norm.(tvmeval(At,tt; method = "linear").-Amat.f.(tt)) .< 1.e-3)
#@test iszero(convert(PeriodicSymbolicMatrix,Amat).F-Ap.F)
@test isperiodic(Amat)

Amat = PeriodicFunctionMatrix(tA,2pi);
@test all(norm.(tvmeval(At,tt; method = "linear").-tvmeval(Amat,tt)) .< 1.e-3)
@test iszero(convert(PeriodicSymbolicMatrix,Amat).F-Ap.F)
@test isperiodic(Amat)

Amat =  convert(PeriodicFunctionMatrix,Ap);
@test all(norm.(tvmeval(At,tt; method = "linear").-tvmeval(Ap,tt)) .< 1.e-3)
@test iszero(convert(PeriodicSymbolicMatrix,Amat).F-Ap.F)
@test isperiodic(Amat) && isperiodic(Ap) && size(Amat) == size(Ap)

for method in ("constant", "linear", "quadratic", "cubic")
      Amat = ts2pfm(At; method);
      @test all(norm.(At.values.-Amat.f.(tg)) .< 1.e-10)
      Bmat = ts2pfm(Bt; method);
      @test all(norm.(Bt.values.-Bmat.f.(tg)) .< 1.e-10)
      Cmat = ts2pfm(Ct; method);
      @test all(norm.(Ct.values.-Cmat.f.(tg)) .< 1.e-10)
end

# example of Colaneri
at(t) = [0 1; -10*cos(t) -24-10*sin(t)];
Afun=PeriodicFunctionMatrix(at,2pi);
ev = pseig(Afun; solver = "non-stiff", reltol = 1.e-10, abstol = 1.e-10)
# @time Φ = tvstm(Afun.f, Afun.period; reltol = 1.e-10)
cvals = log.(complex(ev))/2/pi
println("cvals = $cvals  -- No one digit accuracy!!!")
@test maximum(abs.(ev)) ≈ 1

# using ApproxFun
s = Fourier(0..2π)
At = FourierFunctionMatrix(Fun(t -> [0 1; -10*cos(t) -24-10*sin(t)],s), 2pi)
Atfun = convert(PeriodicFunctionMatrix,At)
Ahrfun = convert(PeriodicFunctionMatrix,pfm2hr(Afun))

@time cvals = psceig(Afun, 500; solver = "non-stiff", reltol = 1.e-10, abstol = 1.e-10)
@time cvals1 = psceig(Atfun, 500; solver = "non-stiff", reltol = 1.e-10, abstol = 1.e-10)
@time cvals2 = psceig(Ahrfun, 500; solver = "non-stiff", reltol = 1.e-10, abstol = 1.e-10)
@test cvals ≈ cvals1 ≈ cvals2
# Tt = Fun(t -> [12+5*sin(t) 1/2; 1 0],s)
# Tinvt=inv(Tt)
# Atilde=Tt*At.M*Tinvt+Tt'*Tinvt
# Aref = Fun(t -> [0 0; 2 -24-10*sin(t)],s)
# @test norm(Aref-Atilde) < 1.e-10

# example Floquet analysis from ApproxFun.jl
a=0.15
at1(t) = -[0 -1 0 0; (2+a*cos(2t)) 0 -1 0; 0 0 0 -1; -1 0 (2+a*cos(2t)) 0]
Afun1=PeriodicFunctionMatrix(at1,pi);
ev1 = pseig(Afun1; solver = "non-stiff", reltol = 1.e-10, abstol = 1.e-10)
cvals1 = log.(complex(ev1))/pi

a=0.15
at2(t) = -[0 -1 0 0; (2+a*cos(2t)) 0 -1 0; 0 0 0 -1; -1 0 (2+a*cos(2t)) 0]
Afun2=PeriodicFunctionMatrix(at2,2*pi;nperiod=2);
ev2 = pseig(Afun2; solver = "non-stiff", reltol = 1.e-10, abstol = 1.e-10)
cvals2 = log.(complex(ev2))/(2pi)
@test ev1.^2 ≈ ev2 && real(cvals1) ≈ real(cvals2)

# full accuracy characteristic exponents
# solver = "symplectic"
# for solver in ("non-stiff", "stiff", "linear", "symplectic", "noidea")
#     @time M = monodromy(Afun, 500; solver, reltol = 1.e-10, abstol = 1.e-10);
#     cvals = log.(complex(pseig(M)))/2/pi
#     #println("solver = $solver cvals = $cvals")
#     @test isapprox(cvals, [0; -24], atol = 1.e-7)
# end

solver = "non-stiff"
for solver in ("non-stiff", "stiff", "linear", "symplectic", "noidea")
    @time cvals = psceig(Afun, 500; solver, reltol = 1.e-10, abstol = 1.e-10)
    @test isapprox(cvals, [0; -24], atol = 1.e-7)
end

# Vinograd example: unstable periodic system with all A(t) stable
#at(t) = [3.5 6;-6 -5.5]+[-4.5 6; 6 4.5]*cos(12*t)+[6 4.5;4.5 -6]*sin(12*t); T = pi/6; # does not work
function at(t::Real)
    [-1-9*(cos(6*t))^2+12*sin(6*t)*cos(6*t) 12*(cos(6*t))^2+9*sin(6*t)*cos(6*t);
    -12*(sin(6*t))^2+9*sin(6*t)*cos(6*t) -1-9*(sin(6*t))^2-12*sin(6*t)*cos(6*t)]
end
@test eigvals(at(rand())) ≈ [-10,-1]
T = pi/3;
Afun=PeriodicFunctionMatrix(at,T);
for solver in ("non-stiff", "stiff", "linear", "symplectic", "noidea")
      @time cvals = psceig(Afun, 500; solver, reltol = 1.e-10)
      @test cvals ≈ [2; -13]
end

# #using BenchmarkTools
# solver = "non-stiff"
# at1(t) = [0 1; -10*cos(t) -24-10*sin(t)];
# Afun1=PeriodicFunctionMatrix(at1,2pi);
# solvers = ("non-stiff", "stiff", "linear", "symplectic", "noidea")
# global k = 1
# for i = 1:5
#       println("solver = $(solvers[k])")
#       @btime  monodromy(Afun1, 500; solver = solvers[k], reltol = 1.e-10, abstol = 1.e-10);
#       k +=1
#       # @btime psceig(Afun, 500; solver, reltol = 1.e-10)
# end

# at1(t) = [0 1; -10*cos(t) -24-10*sin(t)];
# Afun1=PeriodicFunctionMatrix(at1,2pi);
# println("solver = non-stiff")
# @btime  monodromy(Afun1, 500; solver = "non-stiff", reltol = 1.e-10, abstol = 1.e-10);
# println("solver = stiff")
# @btime  monodromy(Afun1, 500; solver = "stiff", reltol = 1.e-10, abstol = 1.e-10);
# println("solver = linear")
# @btime  monodromy(Afun1, 500; solver = "linear", reltol = 1.e-10, abstol = 1.e-10);
# println("solver = symplectic")
# @btime  monodromy(Afun1, 500; solver = "symplectic", reltol = 1.e-10, abstol = 1.e-10);
# println("solver = noidea")
# @btime  monodromy(Afun1, 500; solver = "noidea", reltol = 1.e-10, abstol = 1.e-10);

# at(t) = [BigFloat(0) 1; -10*cos(t) -24-10*sin(t)];
# Afun=PeriodicFunctionMatrix(at,2pi);

# f(u,p,t) = at(t)*u;
# u0 = Matrix{BigFloat}(I,2,2);
# tspan = (BigFloat(0.),2BigFloat(pi))
# prob = ODEProblem(f, u0, tspan)
# sol = solve(prob, Vern9(), rtol = 1.e-50, atol = 1.e-50, save_evrystep = false)
# cvals = log.(complex(eigvals(sol[end])))/2/pi

# at(t) = [BigFloat(0) 1; -10*cos(t) -24-10*sin(t)];
# f!(du,u,p,t) = mul!(du,at(t),u)
# u0 = Matrix{BigFloat}(I,2,2);
# tspan = (BigFloat(0.),2BigFloat(pi))
# prob1 = ODEProblem(f!, u0, tspan)
# sol1 = solve(prob1, Vern9(), rtol = 1.e-50, atol = 1.e-50, save_evrystep = false)
# cvals1 = log.(complex(eigvals(sol1[end])))/2/pi

# at(t) = [BigFloat(0) 1; -10*cos(t) -24-10*sin(t)];
# Afun=PeriodicFunctionMatrix(at,2BigFloat(pi));
# @time Φ = tvstm(Afun, Afun.period; reltol = 1.e-50, abstol = 1.e-50)  # running forever???
# cvals = log.(complex(eigvals(Φ)))/2/pi
# println("cvals = $cvals  -- No one digit accuracy")
# @test maximum(abs.(eigvals(Φ))) > 1

# @variables x
# Cx(x::Real) = [ sin(x)+cos(2*x) 1]

# Anum = PeriodicSymbolicMatrix(Cx(x),2pi)

# # using the lifted periodic system
# ap = [A3,A2,A1]; ep = [E3,E2,E1]; bp = [rand(3,0), rand(3,0), rand(3,0)];
# cp = [rand(0,3), rand(0,3), rand(0,3)]; dp = [rand(0,0), rand(0,0), rand(0,0)];
# psys = PeriodicDiscreteDescriptorStateSpace{Float64}(ap,ep,bp,cp,dp,-1)
# poles = gpole(dss(psys))

# M, N = psreduc_fast([A3,A2,A1],[E3,E2,E1],atol = 1.e-7);
# poles = pzeros(M,N)[1]

# A1 = rand(3,4); A2 = rand(4,3); A = [A1,A2];
# E1 = rand(3,3); E2 = rand(4,4); E = [E1,E2];

# M, N = psreduc_fast(A,E,atol = 1.e-7)
# #E1 = eye(Float64,3); E2 = rand(4,4); E = [E1,E2];
# B = [rand(3,2), rand(4,1)];
# C = [rand(2,4), rand(3,3)];
# D = [rand(2,2), rand(3,1)];
# psys = PeriodicDiscreteDescriptorStateSpace{Float64}(A,E,B,C,D,-1)
# sys = dss(psys)

# ev = eigvals(inv(E2)*A2*inv(E1)*A1);
# poles = gpole(sys);

# @test sort(real(poles)) ≈ sort(real(ev)) &&
#       sort(imag(poles)) ≈ sort(imag(ev))

# sys = dss(psys, kstart = 2, compacted = true, atol = 1.e-7)

# ev = eigvals(inv(E1)*A1*inv(E2)*A2);
# poles = gpole(sys);

# @test sort(real(poles)) ≈ sort(real(ev)) &&
#       sort(imag(poles)) ≈ sort(imag(ev))

# psys = PeriodicDiscreteDescriptorStateSpace{Float64}(A,I,B,C,D,-1)
# sys = dss(psys)

# ev = eigvals(A2*A1);
# poles = gpole(sys);

# @test sort(real(poles)) ≈ sort(real(ev)) &&
#       sort(imag(poles)) ≈ sort(imag(ev))

# sys = dss(psys, kstart = 2, compacted = true, atol = 1.e-7)

# ev = eigvals(A1*A2);
# poles = gpole(sys);

# @test sort(real(poles)) ≈ sort(real(ev)) &&
#       sort(imag(poles)) ≈ sort(imag(ev))

end
end


text:
#-------------------------------------------------------------------------------
# Combination of signals from multiple coils
#

"""
Functor for weighted combination of signal from multiple coils.  Create using
`pca_channel_combiner()`.

Channels combination weights should have
* The correct phase so that signals constructively interfere
* The correct amplitude to achieve the best signal to noise ratio

These factors depend on the relative geometry of the sample vs coil for each
channel, and may be inferred from one or more acquisitions, depending on the
experiment.
"""
struct ChannelCombiner{R}
    weights::Vector{ComplexF64}
    channels::R
end

ChannelCombiner(w, c) = ChannelCombiner(ComplexF64.(w), c)

(combiner::ChannelCombiner)(data) = combine_channels(combiner, data)

"""
    combine_channels(combiner::ChannelCombiner, data)

Combine channels data residing in a `NxC` matrix `data`, where `N` is number of
temporal samples and `C` the number of channels.  Alternatively, `data` can be
an acquisition.

See `pca_channel_combiner()` to create a ChannelCombiner functor.
"""
combine_channels(combiner::ChannelCombiner, data::AbstractMatrix) =
    data[:,combiner.channels] * combiner.weights
combine_channels(combiner::ChannelCombiner, acq::Acquisition) =
    acq.data[:,combiner.channels] * combiner.weights

"""
    pca_channel_combiner(signals; signal_range=nothing, channels=:)

Compute channel combination object using PCA, assuming that each acquisition in
`signals` has the same relative geometry of sample and coils.  `signals[i]` is
assumed to be an array of shape `num_samp × num_chan`.  We use the first few
samples of all acquisition data together in the same calculation to get a good
estimate of the correct channel weights.

These assumptions should be valid to assess the relative coil SNR for
spectroscopy of a single voxel; for other experiments you may want to use a
different method, or restrict calculation of weights to a different
`signal_range`.

`channels` can be set to a index-like object in order to combine only a subset
of channels.
"""
function pca_channel_combiner(signals; signal_range=nothing, channels=:)
    # Select set of time series data from multiple acquisitions, containing
    # mostly signal rather than noise.
    if signal_range == nothing
        # Crude heuristic: Try to focus on a part of the time series containing
        # mostly signal (ie, the start).
        #
        # Using the start of the signal will fail if it's an echo, so arguably
        # we could do something better here...
        signal_range = 1:100
    end
    signal = ComplexF64.(vcat([a[signal_range,channels] for a in signals]...))

    # Compute PCA via the correlation matrix
    eig = eigen(Hermitian(signal'*signal))
    # Weights arise from the first principle component
    pc1_index = argmax(eig.values)
    weights = eig.vectors[:,pc1_index]

    # TODO:
    # * What's the most appropriate normalization for the weights?
    # * Is there any sensible phase factor we should add here?
    ChannelCombiner((1 ./ sum(abs, weights)) .* weights, channels)
end

###
# Hacks: Need the following to preserve the AxisArray axes :-((
function combine_channels(combiner::ChannelCombiner, data::AxisArray)
    combined = combine_channels(combiner, data.data)
    combined::Array  # Check that we don't double wrap in the future...
    AxisArray(combined, AxisArrays.axes(data,Axis{:time}))
end
function Statistics.mean(a::AbstractArray{<:AxisArray}, args...)
	AxisArray(mean((d.data for d in a), args...), AxisArrays.axes(a[1]))
end

ext: jl

max_stars_repo_name: TRIImaging/MagneticResonanceSignals

max_issues_repo_name: TRIImaging/MagneticResonanceSignals

max_forks_repo_name: TRIImaging/MagneticResonanceSignals


text:
include("silicon_runners.jl")

@testset "Silicon PBE (small, Float64)" begin
    run_silicon_pbe(Float64, Ecut=7, test_tol=0.03, n_ignored=0, grid_size=17)
end

if !isdefined(Main, :FAST_TESTS) || !FAST_TESTS
    @testset "Silicon PBE (large, Float64)" begin
        run_silicon_pbe(Float64, Ecut=25, test_tol=1e-5, n_ignored=0,
                        grid_size=33, scf_tol=1e-7)
    end
end


text:
##########################
### Points Conversions ###
##########################

###############################################
### Conversion from LLA to ECEF coordinates ###
###############################################

function ECEF(lla::OpenStreetMapX.LLA, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84)
    ϕdeg, λdeg, h = lla.lat, lla.lon, lla.alt
    d = datum

    sinϕ, cosϕ = sind(ϕdeg), cosd(ϕdeg)
    sinλ, cosλ = sind(λdeg), cosd(λdeg)

    N = d.a / sqrt(1 - d.e² * sinϕ^2)  # Radius of curvature (meters)

    x = (N + h) * cosϕ * cosλ
    y = (N + h) * cosϕ * sinλ
    z = (N * (1 - d.e²) + h) * sinϕ

    return OpenStreetMapX.ECEF(x, y, z)
end

###############################################
### Conversion from ECEF to LLA coordinates ###
###############################################

function LLA(ecef::OpenStreetMapX.ECEF, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84)
    x, y, z = ecef.x, ecef.y, ecef.z
    d = datum

    p = hypot(x, y)
    θ = atan(z*d.a, p*d.b)
    λ = atan(y, x)
    ϕ = atan(z + d.e′² * d.b * sin(θ)^3, p - d.e²*d.a*cos(θ)^3)

    N = d.a / sqrt(1 - d.e² * sin(ϕ)^2)  # Radius of curvature (meters)
    h = p / cos(ϕ) - N

    return OpenStreetMapX.LLA(rad2deg(ϕ), rad2deg(λ), h)
end

###############################################
### Conversion from ECEF to ENU coordinates ###
###############################################

# Given a reference point for linarization
function ENU(ecef::OpenStreetMapX.ECEF, lla_ref::OpenStreetMapX.LLA, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84)
    ϕdeg, λdeg = lla_ref.lat, lla_ref.lon

    ecef_ref = OpenStreetMapX.ECEF(lla_ref, datum)
    ∂x = ecef.x - ecef_ref.x
    ∂y = ecef.y - ecef_ref.y
    ∂z = ecef.z - ecef_ref.z

    # Compute rotation matrix
    sinλ, cosλ = sind(λdeg), cosd(λdeg)
    sinϕ, cosϕ = sind(ϕdeg), cosd(ϕdeg)

    # R = [     -sinλ       cosλ  0.0
    #      -cosλ*sinϕ -sinλ*sinϕ cosϕ
    #       cosλ*cosϕ  sinλ*cosϕ sinϕ]
    #
    # east, north, up = R * [∂x, ∂y, ∂z]
    east  = ∂x * -sinλ      + ∂y * cosλ       + ∂z * 0.0
    north = ∂x * -cosλ*sinϕ + ∂y * -sinλ*sinϕ + ∂z * cosϕ
    up    = ∂x * cosλ*cosϕ  + ∂y * sinλ*cosϕ  + ∂z * sinϕ

    return OpenStreetMapX.ENU(east, north, up)
end

# Given Bounds object for linearization
ENU(ecef::OpenStreetMapX.ECEF, bounds::OpenStreetMapX.Bounds{OpenStreetMapX.LLA}, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84) = OpenStreetMapX.ENU(ecef, OpenStreetMapX.center(bounds), datum)

###############################################
### Conversion from ENU to ECEF coordinates ###
###############################################

function ECEF(enu::OpenStreetMapX.ENU, lla_ref::OpenStreetMapX.LLA, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84)
    ϕdeg, λdeg = lla_ref.lat, lla_ref.lon

    ecef_ref = OpenStreetMapX.ECEF(lla_ref, datum)

    # Compute rotation matrix
    sinλ, cosλ = sind(λdeg), cosd(λdeg)
    sinϕ, cosϕ = sind(ϕdeg), cosd(ϕdeg)

    # R = [-sinλ -sinϕ*cosλ  cosϕ*cosλ
    #      	cosλ -sinϕ*sinλ  cosϕ*sinλ
    #        0.0       cosϕ       sinϕ]
    #
    # x,y,z = R * [east, north, up] + [ecef_ref.x, ecef_ref.y, ecef_ref.z]

    x  = -sinλ*enu.east + -sinϕ*cosλ*enu.north + cosϕ*cosλ*enu.up + ecef_ref.x
    y = cosλ*enu.east +  -sinϕ*sinλ*enu.north + cosϕ*sinλ*enu.up + ecef_ref.y
    z    = 0.0*enu.east + cosϕ*enu.north + sinϕ*enu.up + ecef_ref.z

    return OpenStreetMapX.ECEF(x, y, z)
end

# Given Bounds object for linearization
ECEF(enu::OpenStreetMapX.ENU, bounds::OpenStreetMapX.Bounds{OpenStreetMapX.LLA}, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84) = OpenStreetMapX.ECEF(enu, OpenStreetMapX.center(bounds), datum)

##############################################
### Conversion from LLA to ENU coordinates ###
##############################################

# Given a reference point for linarization
ENU(lla::OpenStreetMapX.LLA, lla_ref::OpenStreetMapX.LLA, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84) = OpenStreetMapX.ENU(OpenStreetMapX.ECEF(lla, datum), lla_ref, datum)

# Given Bounds object for linearization
ENU(lla::OpenStreetMapX.LLA, bounds::OpenStreetMapX.Bounds{OpenStreetMapX.LLA}, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84) = OpenStreetMapX.ENU(OpenStreetMapX.ECEF(lla, datum), bounds, datum)

##############################################
### Conversion from ENU to LLA coordinates ###
##############################################

# Given a reference point for linarization
LLA(enu::OpenStreetMapX.ENU, lla_ref::OpenStreetMapX.LLA, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84) = OpenStreetMapX.LLA(OpenStreetMapX.ECEF(enu,lla_ref))

# Given Bounds object for linearization
LLA(enu::OpenStreetMapX.ENU, bounds::OpenStreetMapX.Bounds{OpenStreetMapX.LLA}, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84) = OpenStreetMapX.LLA(OpenStreetMapX.ECEF(enu,bounds))

#########################################
### Dictionaries of Nodes Conversions ###
#########################################

###############################################
### Conversion from LLA to ECEF coordinates ###
###############################################

function ECEF(nodes::Dict{Int,OpenStreetMapX.LLA}, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84)
    r = Dict{Int,OpenStreetMapX.ECEF}()
    sizehint!(r, ceil(Int, 1.5*length(nodes)))

    for (key, node) in nodes
        r[key] = OpenStreetMapX.ECEF(node, datum)
    end

    return r
end

###############################################
### Conversion from ECEF to LLA coordinates ###
###############################################

function LLA(nodes::Dict{Int,OpenStreetMapX.ECEF}, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84)
    r = Dict{Int,OpenStreetMapX.LLA}()
    sizehint!(r, ceil(Int, 1.5*length(nodes)))

    for (key, node) in nodes
        r[key] = OpenStreetMapX.LLA(node, datum)
    end

    return r
end

######################################################
### Conversion from LLA or ECEF to ENU coordinates ###
######################################################

# Given a reference point
function ENU(nodes::Dict{Int,T}, lla_ref::OpenStreetMapX.LLA,
            datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84) where T<:Union{OpenStreetMapX.LLA,OpenStreetMapX.ECEF}
    r = Dict{Int,OpenStreetMapX.ENU}()
    sizehint!(r, ceil(Int, 1.5*length(nodes)))

    for (key, node) in nodes
        r[key] = OpenStreetMapX.ENU(node, lla_ref, datum)
    end

    return r
end

# Given Bounds
ENU(nodes::Dict{Int,T}, bounds::OpenStreetMapX.Bounds{OpenStreetMapX.LLA}, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84) where T<:Union{OpenStreetMapX.LLA,OpenStreetMapX.ECEF} = OpenStreetMapX.ENU(nodes, OpenStreetMapX.center(bounds), datum)

###############################################
### Conversion from ENU to ECEF coordinates ###
###############################################

function ECEF(nodes::Dict{Int,OpenStreetMapX.ENU},lla_ref::OpenStreetMapX.LLA , datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84)
    r = Dict{Int,OpenStreetMapX.ECEF}()
    sizehint!(r, ceil(Int, 1.5*length(nodes)))

    for (key, node) in nodes
        r[key] = OpenStreetMapX.ECEF(node, lla_ref, datum)
    end

    return r
end

# Given Bounds
ECEF(nodes::Dict{Int,OpenStreetMapX.ENU}, bounds::Bounds{LLA}, datum::Ellipsoid = WGS84) = ECEF(nodes, OpenStreetMapX.center(bounds), datum)

###############################################
### Conversion from ENU to LLA coordinates ###
###############################################

function LLA(nodes::Dict{Int,OpenStreetMapX.ENU},lla_ref::OpenStreetMapX.LLA , datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84)
    r = Dict{Int,OpenStreetMapX.LLA}()
    sizehint!(r, ceil(Int, 1.5*length(nodes)))

    for (key, node) in nodes
        r[key] = OpenStreetMapX.LLA(node, lla_ref, datum)
    end

    return r
end

# Given Bounds
LLA(nodes::Dict{Int,OpenStreetMapX.ENU}, bounds::OpenStreetMapX.Bounds{OpenStreetMapX.LLA}, datum::OpenStreetMapX.Ellipsoid = OpenStreetMapX.WGS84) = OpenStreetMapX.LLA(nodes, OpenStreetMapX.center(bounds), datum)


text:
function read_entry(io, sdf, name)
    blocks = sdf.blocks
    data_block = getindex(blocks, name)
    cache = get_cache(data_kind(typeof(data_block)), sdf)
    data = cached_read(io, data_block, cache)
    @debug "Reading $name"
    store_entry(data_block, data, io, blocks)
end

function read_selected(file, blocks, name, skip_grid)
    data_block = getindex(blocks, name)
    # don't read grid for particles
    isnothing(get_mesh_id(data_block)) && return nothing
    data = read(file, data_block)

    if get_mesh_id(data_block) in skip_grid
        (data=data, block=data_block)
    else
        store_entry(data_block, data, file, blocks)
    end
end

function get_mesh_id(file, idx)
    block = getindex(file.blocks, idx)
    get_mesh_id(block)
end

get_cache(::Data, sdf) = sdf.cache[]
get_cache(::Grid, sdf) = sdf.particle_cache[]

get_mesh_id(block::AbstractBlockHeader) = hasproperty(block, :mesh_id) ? Symbol(block.mesh_id) : nothing

store_entry(data_block::T, data, file, blocks) where T = store_entry(data_kind(T), data_block, data, file, blocks)

store_entry(::Grid, ::T, data, file, blocks) where T = store_entry(discretization_type(T), data)

function store_entry(::Data, data_block, data, file, blocks)
    mesh_block = getindex(blocks, get_mesh_id(data_block))
    grid = make_grid(mesh_block, data_block, file)

    @debug "Grid for data: $grid"
    store_entry(data_block, data, grid)
end

store_entry(data_block::T, data, grid) where T = store_entry(discretization_type(T), data_block, data, grid)
store_entry(::Variable, data::NTuple) = VectorVariable(data)

function store_entry(::StaggeredField, data_block, data, grid)
    name = nameof(data_block)
    ScalarField(data, grid, name)
end

function store_entry(::Variable, data_block, data, grid)
    name = nameof(data_block)
    ScalarVariable(data, grid, name)
end

make_grid(mesh_block::T, data_block, file) where T =
    make_grid(discretization_type(T), mesh_block, data_block, file)

function make_grid(::StaggeredField, mesh_block, data_block, file)
    units = get_units(mesh_block.units)
    minval = mesh_block.minval .* units
    maxval = mesh_block.maxval .* units
    dims = mesh_block.dims

    original_grid = map(eachindex(dims)) do i
        range(minval[i], maxval[i], length=dims[i])
    end
    @debug "Creating grid from $original_grid"

    stagger = data_block.stagger
    @debug "Staggering: $stagger"
    grid = apply_stagger(original_grid, Val(stagger))

    # Fix grid in the cases where it doesn't match the data. See
    # https://cfsa-pmw.warwick.ac.uk/SDF/SDF_C/-/blob/master/src/sdf_control.c#L775-780
    for i in axes(grid)[1]
        if length(grid[i]) ≠ data_block.dims[i]
            grid = setindex!!(grid, grid[i][begin:end-1], i)
        end
    end

    names = Symbol.(lowercase.(labels(mesh_block)))
    @debug "Axis label names: $names"

    return SparseAxisGrid(grid; names)
end

function make_grid(::Variable, mesh_block, data_block, file; cache=nothing)
    if isnothing(cache)
        @debug "No cache"
        grid = read(file, mesh_block)
    else
        @debug "Got cache"
        grid = cached_read(file, mesh_block, cache)
    end
    units = get_units(mesh_block.units)

    minvals = (mesh_block.minval...,) .* units
    maxvals = (mesh_block.maxval...,) .* units

    names = Symbol.(lowercase.(labels(mesh_block)))
    @debug "Axis label names: $names"

    ParticlePositions(grid; names, mins=MVector(minvals), maxs=MVector(maxvals))
end

apply_stagger(grid, ::Val{CellCentre}) = (midpoints.(grid)...,)

function apply_stagger(grid, ::Val{FaceX})
    n = length(grid)
    if n == 1
        grid
    # TODO: 2D
    else
        (grid[1], midpoints.(grid[2:end])...)
    end
end

function apply_stagger(grid, ::Val{FaceY})
    n = length(grid)
    if n == 1
        grid
    elseif n == 2
        (midpoints(grid[1]), grid[2])
    else
        (midpoints(grid[1]), grid[2], midpoints(grid[3]))
    end
end

function apply_stagger(grid, ::Val{FaceZ})
    n = length(grid)
    if n == 1
        grid
    else
        (midpoints.(grid[1:2])..., grid[3])
    end
end

function apply_stagger(grid, ::Val{EdgeX})
    n = length(grid)
    if n == 1
        grid
    else
        (midpoints(grid[1]), grid[2:end]...)
    end
end

function apply_stagger(grid, ::Val{EdgeY})
    n = length(grid)
    if n == 1
        grid
    else
        (grid[1], midpoints(grid[2]), grid[3])
    end
end

function apply_stagger(grid, ::Val{EdgeZ})
    n = length(grid)
    if n == 1
        grid
    else
        (grid[1:2]..., midpoints(grid[3]))
    end
end

apply_stagger(grid, ::Val{Vertex}) = (grid...,)


text:
#Julia Gadfly Scatterplot

using Gadfly
srand(111)
plot(x=rand(7), y=rand(7))


text:
@testset "Kernel Indexing" begin

function idx_kern(X)
    X[1] = workitemIdx().x
    X[2] = workitemIdx().y
    X[3] = workitemIdx().z

    X[4] = workgroupIdx().x
    X[5] = workgroupIdx().y
    X[6] = workgroupIdx().z

    nothing
end

A = zeros(Int64, 6)
HA = HSAArray(A)
@roc groupsize=(1,2,3) gridsize=(4,5,6) idx_kern(HA)
A = Array(HA)
@test all(A .> 0)

function dim_kern(X)
    X[1] = workgroupDim().x
    X[2] = workgroupDim().y
    X[3] = workgroupDim().z

    X[4] = gridDim().x
    X[5] = gridDim().y
    X[6] = gridDim().z

    nothing
end

A = zeros(Int64, 6)
HA = HSAArray(A)
wait(@roc groupsize=(1,2,3) gridsize=(4,5,6) dim_kern(HA))
A = Array(HA)
@test A ≈ [1,2,3,4,5,6]

end


text:
# We need some simple way that julia packages can "call home" and indicate use.
# This is needed since our funding agencies wants proof that our software
# and solutions are of use in society.
function async_send_to_server_tcp(message, url = "127.0.0.1", port = 2001)
  @async begin
    socket = connect(url, port)
    println(socket, message)
    close(socket)
  end
end

function start_test_server(port = 2001)
  @async begin
    server = listen(port)
    while true
      socket = accept(server)
      tstr = strftime("%Y%m%d %H:%M.%S", time())
      println("$(tstr), Server received: ", readline(socket), "\n")
    end
  end
end


text:
include("benchmark_utils.jl")

function simple_benchmark(ntimes = 5)
    ret = []
    function _simple_benchmark!(setup_ex, ex, desc)
        stats = @nbenchmark_freshexec ntimes = $(ntimes) $(setup_ex) $(ex)
        @info desc ntimes setup_ex ex stats
        push!(ret, stats)
    end

    let
        setup_ex = :(using JET)
        ex       = :(@profile_call identity(nothing))
        _simple_benchmark!(setup_ex, ex, "benchmark first time performance")
    end

    let
        setup_ex = quote
            using JET
            @profile_call identity(nothing)
        end
        ex       = :(@profile_call sum("julia"))
        _simple_benchmark!(setup_ex, ex, "benchmark easy error reporting")
    end

    let
        setup_ex = quote
            using JET
            @profile_call sum("julia")
        end
        ex       = :(@profile_call sum("julia"))
        _simple_benchmark!(setup_ex, ex, "benchmark analysis for cached frame")
    end

    let
        setup_ex = quote
            using JET
            @profile_call sum("julia")
        end
        ex       = :(@profile_call rand(Bool))
        _simple_benchmark!(setup_ex, ex, "benchmark a bit complex call")
    end

    return ret
end


text:
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "startin"
version = v"0.4.9"

sources = [
    ArchiveSource(
        "https://github.com/evetion/startin/archive/c-interface.zip",  # PR branch
        "d3a987e9ab20ee3504c8c79f0940ce22527e23ed8b029e8d12136831fe996c63"
    )
]

script = raw"""
cd $WORKSPACE/srcdir/startin-c-interface/
if [[ "${target}" == *-darwin* ]] || [[ "${target}" == *-freebsd* ]]; then
    # Fix linker for BSD platforms
    sed -i "s/${rust_target}-gcc/${target}-gcc/" "${CARGO_HOME}/config"
fi
if [[ "${target}" == *-w64-mingw32* ]]; then
    # Fix from https://github.com/rust-lang/rust/issues/32859#issuecomment-573423629, see https://github.com/rust-lang/rust/issues/47048
    cp -f /opt/${target}/${target}/sys-root/lib/{,dll}crt2.o `rustc --print sysroot`/lib/rustlib/${rust_target}/lib
fi
cargo build --features c_api --release -j${nproc}
mkdir ${libdir}
if [[ "${target}" == *-w64-mingw32* ]]; then
    # Windows generates .dlls without the lib prefix
    cp target/${rust_target}/release/startin.dll ${libdir}/libstartin.dll
else
    cp target/${rust_target}/release/libstartin.${dlext} ${libdir}/libstartin.${dlext}
fi
"""

# musl platforms are failing, as is win32
platforms = [
    Platform("x86_64", "freebsd"),
    Platform("aarch64", "linux"; libc="glibc"),
    # Platform("aarch64", "linux"; libc="musl"),
    Platform("armv7l", "linux"; libc="glibc"),
    # Platform("armv7l", "linux"; libc="musl"),
    Platform("i686", "linux"; libc="glibc"),
    # Platform("i686", "linux"; libc="musl"),
    Platform("powerpc64le", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="glibc"),
    # Platform("x86_64", "linux"; libc="musl"),
    Platform("x86_64", "macos"),
    # Platform("i686", "windows"),  # linking error
    Platform("x86_64", "windows"),
]

products = [
    LibraryProduct("libstartin", :libstartin),
]

dependencies = Dependency[]

build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; compilers=[:c, :rust])

max_stars_repo_name: plut/Yggdrasil

max_issues_repo_name: plut/Yggdrasil

max_forks_repo_name: plut/Yggdrasil


text:
@testset "Correct interval: apply" begin
	datain = DataFrame(
			   datetime=collect(DateTime(2010,1,1,3):Dates.Minute(30):DateTime(2010,1,2,12)),
			   grav = zeros(Float64,67) .+ 9.8,
			   pres = zeros(Float64,67) .+ 1000.,
			   temp = zeros(Float64,67) .+ 27.);
    inan = [5,6];
    datain[:grav][inan] .= NaN;
	datain[:pres][inan] .= NaN;
	corrfile = joinpath(dirname(@__DIR__),"test","input","correctTimeInterval_inputFile.txt");
	dataout = correctinterval(datain,corrfile,includetime=false);
	## Interpolated values
	# check if only :grav column was corrected
	for i in inan
		@test isnan(dataout[:pres][i])
		@test isnan(datain[:grav][i]) # should stay NaN as called without !
		@test dataout[:grav][i] ≈ 9.8
	end
	# Inserted NaNs
	r = findall(x-> x .== DateTime(2010,01,01,08,30,00),datain[:datetime])
	for i = r[1]:r[1]+2
		@test isnan(dataout[:pres][i])
		@test !isnan(datain[:pres][i])
	end
	# Remove step
	r = findall(x-> x .== DateTime(2010,01,02,04),datain[:datetime])
	@test dataout[:temp][r] == datain[:temp][r]
	for i in r[1]+1:size(dataout,1)
		@test dataout[:temp][i] ≈ datain[:temp][i]+10.
		@test dataout[:grav][i] == datain[:grav][i];
	end
	# Set parameters directly and overwrite the input dataframe
	corrpar = DataFrame(column=[1,2,3], id = [3,2,1],
						x1 = [DateTime(2010,01,01,04,30,00),
							  DateTime(2010,01,01,08,00,00),
							  DateTime(2010,01,02,04,00,00)],
					  	x2 = [DateTime(2010,01,01,07,00,00),
							  DateTime(2010,01,01,09,30,09),
							  DateTime(2010,01,02,04,00,00)],
						y1 = [NaN,NaN,10.],y2 = [NaN,NaN,0.0]);
	correctinterval!(datain,corrpar,includetime=false);
	for k in names(datain)
		for i in 1:size(datain,1)
			if typeof(datain[k][i]) != DateTime
				if !isnan(datain[k][i])
					@test datain[k][i] == dataout[k][i];
				else
					isnan(dataout[k][i])
				end
			end
		end
	end
	# Replace values within interval including x1/x2 time
	corrpar = DataFrame(column=[2], id = [5],
						x1 = [DateTime(2010,01,01,03)],
					  	x2 = [DateTime(2010,01,01,04,30)],
						y1 = [0.0],y2 = [1.0]);
	correctinterval!(datain,corrpar);
	@test datain[:pres][1:4] == collect(range(0.,stop=1.,length=4));
	#
	datain = DataFrame(
			   datetime=collect(DateTime(2010,1,1,3):Dates.Minute(60):DateTime(2010,1,2,12)),
			   grav = zeros(Float64,34) .+ 9.8,
			   pres = zeros(Float64,34) .+ 1000.,
			   temp = zeros(Float64,34) .+ 27.);
    datain[:temp][5:end] .+= 3;
	datain[:temp][33:end] .+= 1;
	corrfile = joinpath(dirname(@__DIR__),"test","input","correctTimeInterval_inputFile2.txt");
	dataout = correctinterval(datain,corrfile,includetime=true);
	@test sum(dataout[:temp]) == 27.0*length(dataout[:temp])+2.0
	# all others are unchanged
	for i in [:grav,:pres,:datetime]
		@test dataout[i] == datain[i]
	end
end

@testset "Correct interval: prepare" begin
	dfin = DataFrame(Temp = collect(0.:1.:12.),
		 datetime= collect(DateTime(2000,1,1):Dates.Hour(1):DateTime(2000,1,1,12)))
	dfin[:Temp][[6,10,11]] .= NaN;
	corrpar = prepcorrpar(dfin[:Temp],dfin[:datetime],min_gap=2,defcol=:Temp,defid=2)
	@test size(corrpar) == (1,7)
	@test corrpar[:x1][1] == dfin[:datetime][10]
	@test corrpar[:x2][1] == dfin[:datetime][11]
	@test corrpar[:column][1] == :Temp
	@test corrpar[:id][1] == 2
end


text:
module noninteractive

export main

import Dates
import ErrorTypes: Result, Ok, Err, unwrap
import GeoInterface: coordinates
import GRIB
import HTTP
import Proj4
import PyPlot
import Shapefile
import ZipFile

"""
    fetch_natural_earth(layerurl)::Shapefile.Table

Fetch a Shapefile.Table of the zip file located at `layerurl`
"""
function fetch_natural_earth(layerurl)
    # In a real workflow, you'd have Natural Earth saved somewhere
    #   and load it in, rather than getting from the website every time
    name = tempname()
    HTTP.open("GET", layerurl) do req
        open(name, "w") do io
            write(io, req)
        end
    end
    reader = ZipFile.Reader(name)
    tempdir = mktempdir()
    for file in reader.files
        open(joinpath(tempdir, file.name), "w") do io
            write(io, read(file, String))
        end
    end
    basename = split(layerurl, '/')[end]
    Shapefile.Table(joinpath(tempdir, splitext(basename)[1]))
end

"""
    fetch_cfs(time)::ErrorTypes.Result{String, DomainError}

Downloads the CFS analysis for the given time. Returns the path to a temporary
file if successful. Returns a DomainError if the time cannot represent a valid
analysis time (00/06/12/18Z between April 2011 and a few days ago).
"""
function fetch_cfs(time)::Result{String, DomainError}
    base = "https://www.ncei.noaa.gov/data/climate-forecast-system/access/operational-analysis/6-hourly-by-pressure/"
    tail = ".pgrbhanl.grib2"
    if Dates.hour(time) ∉ [0, 6, 12, 18]
        return Err(DomainError(Dates.hour(time), "The hour must be either 0, 6, 12, or 18"))
    end
    if time < Dates.DateTime(2011, 4, 1, 0) || time > Dates.now(Dates.UTC)
        return Err(DomainError(time, "The date is outside of the period of record (2011-04 to several days ago)"))
    end
    center = Dates.format(time, "yyyy/yyyymm/yyyymmdd/")
    full = "$(base)$(center)cdas1.t$(Dates.hour(time))z$(tail)"
    # Be aware that this is a large file (~80 MB)
    name = tempname()
    HTTP.open("GET", full) do req
        open(name, "w") do io
            write(io, req)
        end
    end
    Ok(name)
end

"""
    subset_indices(lats, lons, bbox)

Create a masking array the same size as `lats` and `lons` that will
subset the portion of lats and lons within `bbox`. `bbox` has fields
`ymin`, `ymax`, `xmin`, `xmax`.
"""
function subset_indices(lats, lons, bbox)
    valid_lats = map(l -> bbox.ymin < l < bbox.ymax, lats)
    valid_lons = map(l -> bbox.xmin < l < bbox.xmax, lons)
    valid_lats .& valid_lons
end

"""
    relative_vorticity(u, v, x, y)

Compute the relative vorticity (ζ). It is recommended to have `x` and `y`
be in a CRS that preserves direction like the Mercator projection. Returns
an `Array` that has the same eltype and shape as `u`.
"""
function relative_vorticity(u, v, x, y)
    ζ = zeros(typeof(u), size(u))
    for col in 2:size(ζ, 2)-1, row in 2:size(ζ, 1)-1,
        ∂v = v[row, col+1] + v[row, col-1] - 2*v[row, col]
        ∂x = (x[row, col+1] - x[row, col-1]) / 2
        ∂u = u[row+1, col] + u[row-1, col] - 2*u[row, col]
        ∂y = (y[row+1, col] - y[row-1, col]) / 2
        ζ[row, col] = ∂v / ∂x - ∂u / ∂y
    end
    ζ
end

"""
    mps_to_kt(x)

Convert `x` from meters per second to knots. The type of the output
is the type of applying `*` and `/` with x and Float64s.
"""
mps_to_kt(x) = x * 60. * 60. / 1852.

function main()
    # Download data and wait on CFS
    cfstime = Dates.DateTime(2013, 2, 9, 12)
    cfsfut = Threads.@spawn fetch_cfs(cfstime)
    statesfut = Threads.@spawn fetch_natural_earth("https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/50m/cultural/ne_50m_admin_1_states_provinces_lakes.zip")
    countriesfut = Threads.@spawn fetch_natural_earth("https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/50m/cultural/ne_50m_admin_0_countries_lakes.zip")

    # Do work that doesn't depend on an Internet download while we wait
    merc = Proj4.Projection(Proj4.epsg[3395]) # Use Mercator because to preserves direction
    wgs84 = Proj4.Projection(Proj4.epsg[4326]) # Input isn't really WGS84 but close enough

    # Just throw errors for this non-operational case
    cfsfname = unwrap(fetch(cfsfut))

    # Open the file
    grib = GRIB.GribFile(cfsfname)

    # Filter for the messages we want (500 hPa Z, T, U, and V)
    wanted = collect(Iterators.filter(grib) do msg
        msg["level"] == 500 && msg["shortName"] ∈ ["t", "gh", "u", "v"]
    end)

    # We don't know their order, but we can use filter to disambiguate them
    # Spawn a thread to uncompress the values as we go because that takes a while
    umsg = only(filter(msg -> msg["shortName"] == "u", wanted))
    ufut = Threads.@spawn umsg["values"]
    vmsg = only(filter(msg -> msg["shortName"] == "v", wanted))
    vfut = Threads.@spawn vmsg["values"]
    heightmsg = only(filter(msg -> msg["shortName"] == "gh", wanted))
    heightfut = Threads.@spawn heightmsg["values"]
    tempmsg = only(filter(msg -> msg["shortName"] == "t", wanted))
    tempfut = Threads.@spawn tempmsg["values"]

    # While the messages could have different grids, we know from prior examination that they don't
    lats = reshape(heightmsg["latitudes"], heightmsg["Ni"], heightmsg["Nj"])
    lons = reshape(heightmsg["longitudes"], heightmsg["Ni"], heightmsg["Nj"])

    # Subset to the NW quadrant so we don't have to deal with periodicity
    goodlats = findall(l -> 0. < l < 84., lats[1, 1:end])
    goodlons = findall(l -> 180. < l < 360., lons[1:end, 1])
    lats = permutedims(lats[goodlons, goodlats])
    lons = map(x -> x - 360, permutedims(lons[goodlons, goodlats]))

    # Calculate Mercator projection (it's fine because it preserves direction)
    xs, ys = map([1, 2]) do i
         a = map(x -> Proj4.transform(wgs84, merc, [x[1], x[2]]), zip(lons, lats))
         getindex.(a, i)
    end

    xs = reshape(xs, size(lats))
    ys = reshape(ys, size(lats))

    # Wait on u and v and compute ζ
    us = permutedims(fetch(ufut)[goodlons, goodlats])
    vs = permutedims(fetch(vfut)[goodlons, goodlats])
    ζ = map(value -> value * 1e5, relative_vorticity(us, vs, xs, ys))

    # Plot filled ζ contours
    PyPlot.figure(figsize=(12., 8.))
    fills = PyPlot.contourf(lons, lats, ζ, 5:50)

    # Next come the state lines
    nestates = fetch(statesfut)
    for row in nestates
        if row.iso_a2 ∈ ["US", "CA"]
            linestring = coordinates(coordinates(Shapefile.shape(row))[1])[1]
            shapex, shapey = map([1, 2]) do i
                a = coordinates(linestring)
                getindex.(a, i)
            end
            PyPlot.plot(shapex, shapey, color="brown", linewidth=1, marker=nothing)
        end
    end

    # Then country lines
    necountries = fetch(countriesfut)
    for row in necountries
        if row.CONTINENT == "North America"
            linestring = coordinates(coordinates(Shapefile.shape(row))[1])[1]
            shapex, shapey = map([1, 2]) do i
                a = coordinates(linestring)
                getindex.(a, i)
            end
            PyPlot.plot(shapex, shapey, color="brown", linewidth=1, marker=nothing)
        end
    end

    # Plot contour lines of height
    heights = map(value -> value / 10, permutedims(fetch(heightfut)[goodlons, goodlats]))
    PyPlot.contour(lons, lats, heights, 492:6:600, colors="k")

    # Plot wind barbs
    ukt = map(mps_to_kt, us)
    vkt = map(mps_to_kt, vs)
    barbstep = 6
    PyPlot.barbs(lons[1:barbstep:end, 1:barbstep:end], lats[1:barbstep:end, 1:barbstep:end], ukt[1:barbstep:end, 1:barbstep:end], vkt[1:barbstep:end, 1:barbstep:end])

    # Plot prettiness
    timestr = Dates.format(cfstime, "yyyy-mm-dd HHZ")
    PyPlot.title("CFS 500 hPa Analysis at $timestr\nGeopotential Height (dam, contour), Relative Vorticity (\$10^{-5} s^{-1}\$, shaded), Wind (kt)")
    PyPlot.xlim(-130, -55)
    PyPlot.ylim(20, 60)
    PyPlot.colorbar(fills)
    PyPlot.savefig("output.png")

    # Safely close the file. It's not super important here because it'd be closed at the end of the process
    #   but it is best practice.
    destroy(grib)

end
#main() # Commented this out because precompilation runs it lol

end # module


text:
function foxWolframMoments(particleArray)
    entries = length(particleArray)
    h0 = 0.0
    hd = 0.0
    # init
    h10 = 0.0
    h20 = 0.0
    h30 = 0.0
    h40 = 0.0

    if entries < 2
        println(stderr, "foxWolframMoments: Not enough entries")
        return [-1., -1., -1., -1.]
    end
	# now do the moments
	@inbounds for i = 1:entries-1
	    h0 += particleArray[i][4]
	    hd += sum(particleArray[i][1:3].^2)
	    @inbounds for j = i+1:entries
    		ptemp = particleArray[i][4]*particleArray[j][4]
	    	cosθ = (particleArray[i][1]*particleArray[j][1] +
	    	        particleArray[i][2]*particleArray[j][2] +
                    particleArray[i][3]*particleArray[j][3]
                   ) / ptemp
            h10 += ptemp*cosθ
		    h20 += ptemp*(1.5 * cosθ^2 - 0.5)
		    h30 += ptemp*(2.5 * cosθ^3 - 1.5*cosθ)
		    h40 += ptemp*(4.375*cosθ^4 - 3.75*cosθ^2 + 0.375)
        end
    end
	h0 += particleArray[entries][4]
	hd += sum(particleArray[entries][1:3].^2)
	h0 = h0*h0 # square it
	# normalize the moments
	h10 = (hd + 2*h10) / h0
	h20 = (hd + 2*h20) / h0
	h30 = (hd + 2*h30) / h0
    h40 = (hd + 2*h40) / h0
    return (h10, h20, h30, h40)
end


text:
function LLegend(parent::Scene; bbox = nothing, kwargs...)

    attrs = merge!(Attributes(kwargs), default_attributes(LLegend))

    @extract attrs (
        halign, valign, padding, margin,
        title, titlefont, titlesize, titlealign, titlevisible,
        labelsize, labelfont, labelcolor, labelhalign, labelvalign,
        bgcolor, strokecolor, strokewidth,
        patchsize, # the side length of the entry patch area
        ncols,
        colgap, rowgap, patchlabelgap,
    )

    decorations = Dict{Symbol, Any}()

    sizeattrs = sizenode!(attrs.width, attrs.height)
    alignment = lift(tuple, halign, valign)

    suggestedbbox = create_suggested_bboxnode(bbox)

    autosizenode = Node{NTuple{2, Optional{Float32}}}((nothing, nothing))

    computedsize = computedsizenode!(sizeattrs, autosizenode)

    finalbbox = alignedbboxnode!(suggestedbbox, computedsize, alignment, sizeattrs, autosizenode)

    scenearea = @lift(IRect2D($finalbbox))

    scene = Scene(parent, scenearea, raw = true, camera = campixel!)

    translate!(scene, (0, 0, 10))

    legendrect = @lift(
        BBox($margin[1], width($scenearea) - $margin[2],
             $margin[3], height($scenearea)- $margin[4]))

    frame = poly!(scene,
        @lift(enlarge($legendrect, repeat([-$strokewidth/2], 4)...)),
        color = bgcolor, strokewidth = strokewidth,
        strokecolor = strokecolor, raw = true)[end]

    entries = Node(LegendEntry[])

    maingrid = GridLayout(bbox = legendrect, alignmode = Outside(20))
    manipulating_grid = Ref(false)

    onany(maingrid.needs_update, margin) do _, margin
        if manipulating_grid[]
            return
        end
        w = determinedirsize(maingrid, Col())
        h = determinedirsize(maingrid, Row())
        if !any(isnothing.((w, h)))
            autosizenode[] = (w + sum(margin[1:2]), h + sum(margin[3:4]))
        end
    end

    entrytexts = LText[]
    entryplots = [AbstractPlot[]]
    entryrects = LRect[]

    titletext = maingrid[1, 1] = LText(scene, text = title, textsize = titlesize,
        halign = titlealign, visible = titlevisible)

    labelgrid = maingrid[2, 1] = GridLayout()

    function relayout()
        manipulating_grid[] = true

        n_entries = length(entries[])
        nrows = ceil(Int, n_entries / ncols[])

        # the grid has twice as many columns as ncols, because of labels and patches
        ncols_with_symbolcols = 2 * min(ncols[], n_entries) # if fewer entries than cols, not so many cols

        for (i, lt) in enumerate(entrytexts)
            irow = (i - 1) ÷ ncols[] + 1
            icol = (i - 1) % ncols[] + 1
            labelgrid[irow, icol * 2] = lt
        end

        for (i, rect) in enumerate(entryrects)
            irow = (i - 1) ÷ ncols[] + 1
            icol = (i - 1) % ncols[] + 1
            labelgrid[irow, icol * 2 - 1] = rect
        end

        for i in labelgrid.nrows : -1 : max((nrows + 1), 2) # not the last row
            deleterow!(labelgrid, i)
        end

        for i in labelgrid.ncols : -1 : max((ncols_with_symbolcols + 1), 2) # not the last col
            deletecol!(labelgrid, i)
        end

        for i in 1:(labelgrid.ncols - 1)
            if i % 2 == 1
                colgap!(labelgrid, i, Fixed(patchlabelgap[]))
            else
                colgap!(labelgrid, i, Fixed(colgap[]))
            end
        end

        for i in 1:(labelgrid.nrows - 1)
            rowgap!(labelgrid, i, Fixed(rowgap[]))
        end

        if titlevisible[]
            if maingrid.nrows == 1
                maingrid[0, 1] = titletext
            end
        else
            if maingrid.nrows == 2
                deleterow!(maingrid, 1)
            end
        end

        manipulating_grid[] = false
        maingrid.needs_update[] = true

        translate!(scene, (0, 0, 10))
    end

    onany(ncols, rowgap, colgap, patchlabelgap, titlevisible) do _, _, _, _, _
        relayout()
    end

    on(entries) do entries

        # first delete all existing labels and patches

        delete!.(entrytexts)
        empty!(entrytexts)

        delete!.(entryrects)
        empty!(entryrects)

        # delete patch plots
        for eplots in entryplots
            # each entry can have a vector of patch plots
            delete!.(scene, eplots)
        end
        empty!(entryplots)

        # the attributes for legend entries that the legend itself carries
        # these serve as defaults unless the legendentry gets its own value set
        preset_attrs = extractattributes(attrs, LegendEntry)

        for (i, e) in enumerate(entries)

            # fill missing entry attributes with those carried by the legend
            merge!(e.attributes, preset_attrs)

            # create the label
            push!(entrytexts, LText(scene,
                text = e.label, textsize = e.labelsize, font = e.labelfont,
                color = e.labelcolor, halign = e.labelhalign, valign = e.labelvalign
                ))

            # create the patch rectangle
            rect = LRect(scene, color = e.patchcolor, strokecolor = e.patchstrokecolor,
                strokewidth = e.patchstrokewidth,
                width = lift(x -> x[1], e.patchsize),
                height = lift(x -> x[2], e.patchsize))
            push!(entryrects, rect)

            # plot the symbols belonging to this entry
            symbolplots = AbstractPlot[
                legendsymbol!(scene, element, rect.layoutnodes.computedbbox, e.attributes)
                for element in e.elements]

            push!(entryplots, symbolplots)
        end
        relayout()
    end

    # entries[] = [
    #     LegendEntry("entry 1", AbstractPlot[]),
    # ]

    # no protrusions
    protrusions = Node(RectSides(0f0, 0f0, 0f0, 0f0))

    # trigger suggestedbbox
    suggestedbbox[] = suggestedbbox[]

    layoutnodes = LayoutNodes{LLegend, GridLayout}(suggestedbbox, protrusions, computedsize, autosizenode, finalbbox, nothing)

    LLegend(scene, entries, layoutnodes, attrs, decorations, entrytexts, entryplots)
end

function legendsymbol!(scene, element::MarkerElement, bbox, defaultattrs::Attributes)
    merge!(element.attributes, defaultattrs)
    attrs = element.attributes

    fracpoints = attrs.markerpoints
    points = @lift(fractionpoint.(Ref($bbox), $fracpoints))
    scatter!(scene, points, color = attrs.color, marker = attrs.marker,
        markersize = attrs.markersize,
        strokewidth = attrs.markerstrokewidth,
        strokecolor = attrs.strokecolor, raw = true)[end]
end

function legendsymbol!(scene, element::LineElement, bbox, defaultattrs::Attributes)
    merge!(element.attributes, defaultattrs)
    attrs = element.attributes

    fracpoints = attrs.linepoints
    points = @lift(fractionpoint.(Ref($bbox), $fracpoints))
    lines!(scene, points, linewidth = attrs.linewidth, color = attrs.color,
        linestyle = attrs.linestyle,
        raw = true)[end]
end

function legendsymbol!(scene, element::PolyElement, bbox, defaultattrs::Attributes)
    merge!(element.attributes, defaultattrs)
    attrs = element.attributes

    fracpoints = attrs.polypoints
    points = @lift(fractionpoint.(Ref($bbox), $fracpoints))
    poly!(scene, points, strokewidth = attrs.polystrokewidth, color = attrs.color,
        strokecolor = attrs.strokecolor,
        raw = true)[end]
end

function Base.getproperty(lentry::LegendEntry, s::Symbol)
    if s in fieldnames(LegendEntry)
        getfield(lentry, s)
    else
        lentry.attributes[s]
    end
end

function Base.setproperty!(lentry::LegendEntry, s::Symbol, value)
    if s in fieldnames(LegendEntry)
        setfield!(lentry, s, value)
    else
        lentry.attributes[s][] = value
    end
end

function Base.propertynames(lentry::LegendEntry)
    [fieldnames(T)..., keys(lentry.attributes)...]
end

function LegendEntry(label::String, plots::Vararg{AbstractPlot}; kwargs...)
    attrs = Attributes(label = label)

    kwargattrs = Attributes(kwargs)
    merge!(attrs, kwargattrs)

    elems = vcat(legendelements.(plots)...)
    LegendEntry(elems, attrs)
end

function LegendEntry(label::String, elements::Vararg{LegendElement}; kwargs...)
    attrs = Attributes(label = label)

    kwargattrs = Attributes(kwargs)
    merge!(attrs, kwargattrs)

    elems = LegendElement[elements...]
    LegendEntry(elems, attrs)
end

function LineElement(;kwargs...)
    LineElement(Attributes(kwargs))
end

function MarkerElement(;kwargs...)
    MarkerElement(Attributes(kwargs))
end
function PolyElement(;kwargs...)
    PolyElement(Attributes(kwargs))
end

function legendelements(plot::Union{Lines, LineSegments})
    LegendElement[LineElement(color = plot.color, linestyle = plot.linestyle)]
end

function legendelements(plot::Scatter)
    LegendElement[MarkerElement(
        color = plot.color, marker = plot.marker,
        strokecolor = plot.strokecolor)]
end

function legendelements(plot::Poly)
    LegendElement[PolyElement(color = plot.color, strokecolor = plot.strokecolor)]
end

function legendelements(plot::T) where T
    error("""
        There is no `legendelements` method defined for plot type $T. This means
        that you can't automatically generate a legend entry for this plot type.
        You can overload this method for your plot type that returns a `LegendElement`
        vector, or manually construct a legend entry from those elements.
    """)
end

function Base.getproperty(legendelement::T, s::Symbol) where T <: LegendElement
    if s in fieldnames(T)
        getfield(legendelement, s)
    else
        legendelement.attributes[s]
    end
end

function Base.setproperty!(legendelement::T, s::Symbol, value) where T <: LegendElement
    if s in fieldnames(T)
        setfield!(legendelement, s, value)
    else
        legendelement.attributes[s][] = value
    end
end

function Base.propertynames(legendelement::T) where T <: LegendElement
    [fieldnames(T)..., keys(legendelement.attributes)...]
end

function Base.push!(legend::LLegend, entry::LegendEntry)
    legend.entries[] = [legend.entries[]; entry]
    nothing
end

function Base.pushfirst!(legend::LLegend, entry::LegendEntry)
    legend.entries[] = [entry; legend.entries[]]
    nothing
end

function Base.push!(legend::LLegend, label::String, plots::Vararg{AbstractPlot}; kwargs...)
    entry = LegendEntry(label, plots...; kwargs...)
    push!(legend, entry)
end

function Base.pushfirst!(legend::LLegend, label::String, plots::Vararg{AbstractPlot}; kwargs...)
    entry = LegendEntry(label, plots...; kwargs...)
    pushfirst!(legend, entry)
end

function LLegend(scene, plots::AbstractArray{<:AbstractPlot}, labels::AbstractArray{String}; kwargs...)
    legend = LLegend(scene; kwargs...)
    if length(plots) != length(labels)
        error("Legend received $(length(plots)) plots but $(length(labels)) labels.")
    end
    legend.entries[] = [LegendEntry(label, plot) for (plot, label) in zip(plots, labels)]
    legend
end

function LLegend(scene, plotgroups::AbstractArray{<:AbstractArray{<:AbstractPlot}}, labels::AbstractArray{String}; kwargs...)
    legend = LLegend(scene; kwargs...)
    if length(plotgroups) != length(labels)
        error("Legend received $(length(plots)) plotgroups but $(length(labels)) labels.")
    end
    legend.entries[] = [LegendEntry(label, plotgroup...) for (plotgroup, label) in zip(plotgroups, labels)]
    legend
end


text:
function ipm_step!(solv::Solver{T}) where T <: Real
    iter = solv.iter
    lp   = solv.lp
    res  = solv.res

    nc  = lp.nc # number of constraints
    nv  = lp.nv # number of variables
    nu  = lp.nu # number of upper bounded variables
    A   = lp.A
    b   = lp.b
    c   = lp.c
    ubi = lp.upb_ind
    ubv = lp.upb_val

    # scaling for problem data
    θxs = iter.s ./ iter.x
    θvw = iter.w ./ iter.v
    @views θxs[ubi] .+= θvw

    # regularization
    r_min = sqrt(eps(T)) # approx 1e-8
    solv.regP .= max.(r_min, solv.regP ./ 10)
    solv.regD .= max.(r_min, solv.regD ./ 10)
    solv.regG  = max(r_min, solv.regG / 10)

    # factorization
    # make three attempts, after increasing regularization
    # see Tulip.jl for details on factorization
    attempt = 0
    while attempt < 4
        try
            update_linear_solver!(solv.ls, θxs, solv.regP, solv.regD)
            break
        catch e
            if !(isa(e, PosDefException) || isa(e, ZeroPivotException))
                # unexpected error
                rethrow(e)
            end

            attempt += 1
            solv.regP .*= 100
            solv.regD .*= 100
            solv.regG  *= 100
        end
    end

    # factorization may not have been complete
    if attempt >= 3
        throw(PosDefException(0))
    end

    # predictor search direction
    Δ  = Iterate{T}(nc, nv, nu)
    # Δc = Iterate{T}(nc, nv, nu)

    δx = zeros(T, nv)
    δy = zeros(T, nc)
    δz = zeros(T, nu)
    solve_augsys!(δx, δy, δz, solv.ls, θvw, ubi, b, c, ubv)
    δ0 = solv.regG + iter.κ / iter.τ - dot(c, δx) + dot(b, δy) - dot(ubv, δz)

    # affine-scaling search direction
    solve_newton_system!(
        Δ, solv.ls, θvw, b, c, ubi, ubv, δx, δy, δz, δ0, iter,
        res.rp, res.ru, res.rd, res.rg,
        - iter.x .* iter.s, # ξxs
        - iter.v .* iter.w, # ξvw
        - iter.τ  * iter.κ  # ξτκ
    )

    # step length for affine-scaling direction
    # γ = (1 - α)^2 * min(β, 1 - α)
    # η = 1 - γ
    α = max_alpha(iter, Δ)
    γ = (oneunit(T) - α)^2 * min(T(1e-1), oneunit(T) - α)
    η = 1 - γ

    # apply mehrotra corrector
    solve_newton_system!(
        Δ, solv.ls, θvw, b, c, ubi, ubv, δx, δy, δz, δ0, iter,
        η .* res.rp, η .* res.ru, η .* res.rd, η .* res.rg,
        - iter.x .* iter.s .+ γ * iter.μ .- Δ.x .* Δ.s,
        - iter.v .* iter.w .+ γ * iter.μ .- Δ.v .* Δ.w,
        - iter.τ  * iter.κ  + γ * iter.μ  - Δ.τ  * Δ.κ
    )

    # calculate new step length, and dampen alpha
    α = max_alpha(iter, Δ)
    α *= T(9995 / 10_000)

    # update iterate
    update_iterate!(iter, Δ, α)

    nothing
end

function solve!(solv::Solver{T}, maxiter::Int) where T <: Real
    slp = solv.lp

    while true
        # update residuals
        update_residuals!(
            solv.res, solv.iter,
            slp.A, slp.b, slp.c,
            slp.upb_ind, slp.upb_val
        )

        # update duality measure
        duality!(solv.iter)

        # check solver status -- exit if needed
        check_status!(
            solv, solv.iter, solv.res, solv.tols,
            slp.A, slp.b, slp.c, slp.upb_ind, slp.upb_val
        )
        if (
            solv.status == SolverOptimal          ||
            solv.status == SolverPrimalInfeasible ||
            solv.status == SolverDualInfeasible
        )
            break
        elseif solv.niter >= maxiter
            solv.status = SolverExceededIterations
            break
        end

        # take next ipm_step
        try
            ipm_step!(solv)
        catch e
            if isa(e, OutOfMemoryError)
                solv.status = SolverExceededMemory
            elseif isa(e, PosDefException) || isa(e, SingularException)
                solv.status = SolverNumericalInstability
            else
                solv.status = SolverUndefined
                rethrow(e)
            end

            break
        end

        solv.niter += 1
    end

    nothing
end


text:
get_input() = map(s -> parse(Int, s), readlines(joinpath(INPUT_DIR, "day01.txt")))

fuelForMass(mass::Int) = div(mass, 3) - 2

function solve(::Puzzle(2019, 1, 1), input)
    map(fuelForMass, get_input()) |> sum
end

recursiveFuel(mass::Int) = begin
    fuel = fuelForMass(mass)

    if fuel <= 0
        return 0
    end

    return fuel + recursiveFuel(fuel)
end

function solve(::Puzzle(2019, 1, 2), input)
    map(recursiveFuel, input) |> sum
end


text:
# Tests for the PWDict type.
# Author: Noah C. Benson <n@nben.net>

@testset "PWDict" begin
    intels = [gensym() => (x,rand(Float64)) for x in 1:2000]
    els = [:b => (2,20.0), :c => (3,30.0), :a => (1,10.0), :d => (4,40.0)]
    p = PWDict{Symbol,Int,Float64}(els...)
    for el in els
        p = push(p, el)
    end
    (lst,mst) = (first(p), pop(p))
    @test lst == (:d => 4)
    (lst,mst) = (first(mst), pop(mst))
    @test lst == (:c => 3)
    (lst,mst) = (first(mst), pop(mst))
    @test lst == (:b => 2)
    (lst,mst) = (first(mst), pop(mst))
    @test lst == (:a => 1)
    @test isempty(mst)
    # make many samples and make sure they resemble the distribution
    counts = Dict(:a => 0.0, :b => 0.0, :c => 0.0, :d => 0.0)
    vcounts = [0.0,0.0,0.0,0.0]
    for _ in 1:100000
        (k,v) = rand(p)
        counts[k] += 1/1000
        vcounts[v] += 1/1000
    end
    @test  8.5 <  counts[:a] < 12.5
    @test 18.5 <  counts[:b] < 22.5
    @test 28.5 <  counts[:c] < 32.5
    @test 38.5 <  counts[:d] < 42.5
    @test  8.5 < vcounts[1] < 12.5
    @test 18.5 < vcounts[2] < 22.5
    @test 28.5 < vcounts[3] < 32.5
    @test 38.5 < vcounts[4] < 42.5
end


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using ContactImplicitMPC
using LinearAlgebra
using Quaternions

# ## Visualizer
vis = ContactImplicitMPC.Visualizer()
ContactImplicitMPC.render(vis)

@show Threads.nthreads()

# ## Simulation
s = get_simulation("centroidal_quadruped_wall", "flat_3D_lc", "flat")
model = s.model
env = s.env

q1_sim = zeros(model.nq)
v1_sim = zeros(model.nq)

v1_sim[1] = 10.0
v1_sim[1] = 10.0
v1_sim[7] = 10.0
v1_sim[10] = 10.0
v1_sim[13] = 10.0
v1_sim[16] = 10.0

q1_sim[3] = 1.0
q1_sim[9] = 1.0
q1_sim[12] = 1.0
q1_sim[15] = 1.0
q1_sim[18] = 1.0

# ## Simulator
sim = simulator(s, 100, h=0.1)

# ## Simulate
q1_sim0 = deepcopy(q1_sim)
RoboDojo.simulate!(sim, q1_sim0, v1_sim)

# ## Visualize
set_light!(vis)
set_floor!(vis, grid=true)
set_background!(vis)
anim = visualize!(vis, model, sim.traj.q; Δt=0.1)

sim.ip.z
sim.ip.θ

sim.traj.q[1][16]
ϕ_func(model, env, sim.traj.q[end])


text:
using OperatorFlux
using Test
using LinearAlgebra

@testset "FourierTransform" begin
    # test constructors
    trafo = FourierTransform(modes = (3, 7))
    @test trafo.modes == (3, 7)
    @test trafo.even == true

    # test forward & inverse
    for mode in [(4,), (4, 6), (4, 7, 9), (5,), (5, 6), (5, 7, 9)]
        z = rand(3, mode..., 5)
        even = Bool(mod(mode[1], 2)) ? false : true
        trafo = FourierTransform(modes = mode, even = even)
        @test norm(
            OperatorFlux.inverse(trafo, OperatorFlux.forward(trafo, z)) - z,
        ) ≤ 20 * eps(norm(z))
    end

    # test truncate_modes & pad_modes
    M = 32
    trafo = FourierTransform(modes = (3, 7, 6))
    c = rand(7, M, M, M, 3)
    inds = [
        [1, 2, 3],
        [1, 2, 3, 4, 5, 6, 7, 26, 27, 28, 29, 30, 31, 32],
        [1, 2, 3, 4, 5, 6, 27, 28, 29, 30, 31, 32],
    ]
    @test all(OperatorFlux.truncate_modes(trafo, c) .== c[:, inds..., :])

    trafo = FourierTransform(modes = (3, 6))
    c = rand(7, M, M, 3)
    inds = [[1, 2, 3], [1, 2, 3, 4, 5, 6, 27, 28, 29, 30, 31, 32]]
    @test all(OperatorFlux.truncate_modes(trafo, c) .== c[:, inds..., :])

    trafo = FourierTransform(modes = (3,))
    c = rand(7, M, 3)
    inds = [[1, 2, 3]]
    @test all(OperatorFlux.truncate_modes(trafo, c) .== c[:, inds..., :])

    trafo = FourierTransform(modes = (12, 5))
    c = rand(27, 16, 14, 1)
    inds = [
        [(1:(div(m, 2) + 1))..., ((s - div(m, 2) + 2):s)...]
        for (s, m) in zip((M, 14), size(c)[2:(end - 1)])
    ]
    inds[1] = collect(1:16)
    @test all(OperatorFlux.pad_modes(trafo, c, (32, 14))[:, inds..., :] .== c)

    # rountrip integration tests
    for M in [32, 33]
        a, b = -π, π
        x = (b - a) .* collect(0:(M - 1)) / M .+ a
        y = copy(x)
        z = sin.(2x) * cos.(2y)' + sin.(2x) * cos.(3y)'
        z = z[:, :, :, :] # emulate channels and batches
        z = permutedims(z, (3, 1, 2, 4))
        even = Bool(mod(M, 2)) ? false : true
        trafo = FourierTransform(modes = (12, 17), even = even)

        a = OperatorFlux.forward(trafo, z)
        b = OperatorFlux.truncate_modes(trafo, a)
        c = OperatorFlux.pad_modes(trafo, b, size(a)[2:(end - 1)])
        d = OperatorFlux.inverse(trafo, c)
        @test norm(d - z) ≤ 3 * eps(norm(z))
    end

    # test Base extensions
    trafo = FourierTransform(modes = (5, 11))
    @test ndims(trafo) == 2
    @test eltype(trafo) == ComplexF32
    @test size(trafo) == [5, 22]
end


text:
#
# Load DFTK-compatible structural information from an external file
# Relies on ASE and other external libraries to do the parsing
#
using PyCall

function load_lattice(file::AbstractString)
    ase = pyimport_e("ase")
    ispynull(ase) && error("Install ASE to load data from exteral files")
    load_lattice(pyimport("ase.io").read(file))
end

function load_atoms(file::AbstractString)
    ase = pyimport_e("ase")
    ispynull(ase) && error("Install ASE to load data from exteral files")
    load_atoms(pyimport("ase.io").read(file))
end

function load_magnetic_moments(file::AbstractString)
    ase = pyimport_e("ase")
    ispynull(ase) && error("Install ASE to load data from exteral files")
    load_magnetic_moments(pyimport("ase.io").read(file))
end


text:
export Model, Parameter, Measurement

abstract type Model end

@doc doc"""
Input parameter of simulation
"""
const Parameter = Dict{String,Any}

const Measurement = Dict{String, Any}

include("lattice.jl")
include("update.jl")
include("estimator.jl")
include("parameter.jl")
include("postproc.jl")


text:
module adapters_jdbc_structured_test

using Test
using Octo.Adapters.JDBC

struct User
end
Schema.model(User, table_name="users")

@test to_sql([FROM User]) == "FROM users"

u = from(User)
@test to_sql([FROM u]) == "FROM users"

end # module adapters_jdbc_structured_test


text:
using SparseArrays
using LinearAlgebra
using MatrixDepot
using BenchmarkTools
using UnicodePlots
using PrettyTables
using ChainPartitioners
using Cthulhu
using Profile

for mtx in [
            #"HB/can_292",
            #"HB/662_bus",
            #"HB/fs_680_2",
            #"HB/can_292",
            #"HB/plat362",
            #"LPnetlib/lp_bandm",
            #"LPnetlib/lp_etamacro",
            #"Pajek/Erdos991",
            "Boeing/ct20stif",
            #"DIMACS10/chesapeake",
            "Schmid/thermal1",
            "Rothberg/3dtube",
           ]
    A = permutedims(1.0 * sparse(mdopen(mtx).A))
    (m, n) = size(A)

    y = zeros(m)
    x = rand(n)
    ref_time = @belapsed(mul!($y, $A, $x))

    mdl = AffineConnectivityModel{Int64}(0, 0, 0, 1)

    #println()
    #println()
    #println(spy(A, maxwidth=50, maxheight=50, title="$mtx"))
    println(mtx)
    rows = []
    for (key, method) in [
        ("dynamic16", DynamicTotalChunker(ConstrainedCost(mdl, VertexCount(), fld(n, 4)))),
        ("quadrangle16", ConvexTotalChunker(ConstrainedCost(mdl, VertexCount(), fld(n, 4)))),
    ]
        setup_time = @belapsed(pack_stripe($A, $method))/ref_time
        Φ = pack_stripe(A, method)
        #@profile pack_stripe(A, method)
        #Profile.print()
        comm = total_value(A, Φ, mdl)
        @assert issorted(Φ.spl)
        @assert Φ.spl[1] == 1
        @assert Φ.spl[end] == n + 1

        push!(rows, [key setup_time comm])
    end
    pretty_table(vcat(rows...), ["method", "setuptime", "communication"])
end


text:
module TestAutodetect
using Test
using UrlDownload
using UrlDownload: autodetect_uri_type
using UrlDownload: File, URL
using UrlDownload: @f_str, @u_str
using DataFrames

@testset "files" begin
    url = "/tmp/data.csv"
    @test typeof(autodetect_uri_type(url)) == File

    url = "~/data.csv"
    @test typeof(autodetect_uri_type(url)) == File
end

@testset "urls" begin
    url = "https://example.com/data.csv"
    @test typeof(autodetect_uri_type(url)) == URL

    url = "http://example.com/data.csv"
    @test typeof(autodetect_uri_type(url)) == URL
end

@testset "download files or urls" begin
    df = DataFrame(x = [1, 3], y = [2, 4])

    url = "https://raw.githubusercontent.com/Arkoniak/UrlDownload.jl/master/data/ext.csv"
    res = urldownload(url) |> DataFrame
    @test df == res

    url = joinpath(@__DIR__, "..", "data", "ext.csv")
    res = urldownload(url) |> DataFrame
    @test df == res

    url = URL("https://raw.githubusercontent.com/Arkoniak/UrlDownload.jl/master/data/ext.csv")
    res = urldownload(url) |> DataFrame
    @test df == res

    url = u"https://raw.githubusercontent.com/Arkoniak/UrlDownload.jl/master/data/ext.csv"
    res = urldownload(url) |> DataFrame
    @test df == res

    url = File(joinpath(@__DIR__, "..", "data", "ext.csv"))
    res = urldownload(url) |> DataFrame
    @test df == res

    url = joinpath(@__DIR__, "..", "data", "ext.csv")
    url = @f_str url
    res = urldownload(url) |> DataFrame
    @test df == res
end

end # module


text:
using Evolutionary

println("Use Genetic Algorithm to fit the parameter")
#         objfun: Objective fitness function
#              N: Search space dimensionality
# initPopulation: Search space dimension ranges as a vector, or initial population values as matrix,
#                 or generation function which produce individual population entities.
# populationSize: Size of the population
#  crossoverRate: The fraction of the population at the next generation, not including elite children,
#                 that is created by the crossover function.
#   mutationRate: Probability of chromosome to be mutated
#              ɛ: Positive integer specifies how many individuals in the current generation
#                 are guaranteed to survive to the next generation.
#                 Floating number specifies fraction of population.

cost_function = build_loss_objective(prob1,Tsit5(),L2Loss(t,data),
                                     maxiters=10000,verbose=false)
N = 1
result, fitness, cnt = ga(cost_function, N;
                   initPopulation = Float64[1.2],
                   populationSize = 100,
                   ɛ = 0.1,
               selection = sus,
               crossover = intermediate(0.25),
               mutation = domainrange(fill(0.5,N)))
@test result[1] ≈ 1.5 atol=3e-1

cost_function = build_loss_objective(prob2,Tsit5(),L2Loss(t,data),
                                     maxiters=10000)
N = 2
result, fitness, cnt = ga(cost_function, N;
                   initPopulation = Float64[1.2,2.8],
                   populationSize = 500,
                   ɛ = 0.1,
               selection = sus,
               crossover = intermediate(0.25),
               mutation = domainrange(fill(0.5,N)))
@test result ≈ [1.5;3.0] atol=3e-1

cost_function = build_loss_objective(prob3,Tsit5(),L2Loss(t,data),
                                     maxiters=10000)
N = 4
result, fitness, cnt = ga(cost_function, N;
                   initPopulation = Float64[1.3,0.8,2.8,1.2],
                   populationSize = 1000,
                   ɛ = 0.1,
               selection = sus,
               crossover = intermediate(0.25),
               mutation = domainrange(fill(0.5,N)))
@test result ≈ [1.5;1.0;3.0;1.0] atol=3e-1