text:
printstyled(stdout,"  merge! behavior and intent\n", color=:light_green)
nx = 100
fs = 100.0
Δ = round(Int64, 1000000/fs)

id = "VV.STA1.00.EHZ"
loc = GeoLoc(lat=46.8523, lon=121.7603, el=4392.0)
loc1 = GeoLoc(lat=rand(), lon=rand(), el=rand())
resp1 = fctoresp(0.2, 1.0)
resp2 = fctoresp(2.0, 1.0)
units = "m/s"
src = "test"
t0 = 0

mkC1() = SeisChannel( id = id, name = "Channel 1",
                  loc = loc, fs = fs, resp = deepcopy(resp1), units = units,
                  src = "test channel 1",
                  notes = [tnote("New channel 1.")],
                  misc = Dict{String,Any}( "P" => 6.1 ),
                  t = [1 t0; nx 0],
                  x = randn(nx) )
mkC2() = SeisChannel( id = id, name = "Channel 2",
                  loc = loc, fs = fs, resp = deepcopy(resp1), units = units,
                  src = "test channel 2",
                  notes = [tnote("New channel 2.")],
                  misc = Dict{String,Any}( "S" => 11.0 ),
                  t = [1 t0+nx*Δ; nx 0],
                  x = randn(nx) )
mkC4() = (C = mkC2(); C.loc = loc1;
                  C.name = "Channel 4";
                  C.src = "test channel 4";
                  C.notes = [tnote("New channel 4.")]; C)
mkC5() = (C = mkC2(); C.loc = loc1;
                  C.name = "Channel 5";
                  C.src = "test channel 5";
                  C.notes = [tnote("New channel 5.")];
                  C.t = [1 t0+nx*2Δ; nx 0];
                  C)
C2_ov() = (nov = 3; C = mkC2(); C.t = [1 t0+(nx-nov)*Δ; nx 0]; C)
C3_ov() = (nov = 3; C = mkC2(); C.t = [1 t0+2*(nx-nov)*Δ; nx 0]; C)

function prandSC(c::Bool)
  if c == true
    C = randSeisChannel(c=true, nx=1000)
  else
    C = randSeisChannel(c=false, nx=10000)
  end
  C.name = randstring(20)
  C.misc = Dict{String,Any}()
  return C
end

function mk_tcat(T::Array{Array{Int64,2},1}, fs::Float64)
  L = length(T)
  τ = Array{Array{Int64,1},1}(undef,L)
  for i = 1:L
    τ[i] = t_expand(T[i], fs)
  end
  tt = sort(unique(vcat(τ...)))
  return t_collapse(tt, fs)
end

# ===========================================================================
printstyled(stdout,"    xtmerge!\n", color=:light_green)
δ = 20000
tmax = div(typemax(Int64),2)-2δ
x = randn(12)
t = sort(rand(0:2δ:tmax, 12))
while (length(unique(t)) < 12)
  push!(t, rand(0:2δ:tmax))
  t = sort(unique(t))
end
x = vcat(x, x[1:6])
t = vcat(t, t[1:6])
xtmerge!(t, x, div(δ,2))
@test length(t) == 12
@test length(x) == 12

# ===========================================================================
printstyled(stdout,"    removal of traces with no data or time info\n", color=:light_green)
S = SeisData(prandSC(false), prandSC(false), prandSC(false), prandSC(false))
S.x[2] = Float64[]
S.x[3] = Float64[]
S.t[4] = Array{Int64,2}(undef,0,0)
merge!(S, v=1)
basic_checks(S)
sizetest(S, 1)

printstyled(stdout,"    ability to handle irregularly-sampled data\n", color=:light_green)
C = prandSC(true)
namestrip!(C)
S = SeisData(C, prandSC(true), prandSC(true))
namestrip!(S)
for i = 2:3
  S.id[i] = identity(S.id[1])
  S.resp[i] = copy(S.resp[1])
  S.loc[i] = deepcopy(S.loc[1])
  S.units[i] = identity(S.units[1])
end
T = merge(S, v=1)
basic_checks(T)
sizetest(T, 1)

# ===========================================================================
printstyled(stdout,"    simple merges\n", color=:light_green)
printstyled(stdout,"      three channels, two w/same params, no overlapping data\n", color=:light_green)
S = SeisData(mkC1(), mkC2(), prandSC(false))

# Is the group merged correctly?
U = deepcopy(S)
merge!(S)
basic_checks(S)
@test S.n == 2
i = findid(id, S)
@test vcat(U.x[1], U.x[2])==S.x[i]
@test mk_tcat(U.t[1:2], fs) == S.t[i]

# Do the notes log the extra source?
@test findfirst([occursin("New channel 1", i) for i in S.notes[i]]) != nothing
@test findfirst([occursin("test channel 1", i) for i in S.notes[i]]) != nothing

# Are dictionaries merging correctly?
@test haskey(S.misc[i], "P")
@test haskey(S.misc[i], "S")

printstyled(stdout, "      \"zipper\" merge I: two channels, staggered time windows, no overlap\n", color=:light_green)
S = SeisData(mkC1(), mkC2())
W = Array{Int64,2}(undef, 8, 2);
for i = 1:8
  W[i,:] = t0 .+ [(i-1)*nx*Δ ((i-1)*nx + nx-1)*Δ]
end
w1 = W[[1,3,5,7], :]
w2 = W[[2,4,6,8], :]
S.t[1] = w_time(w1, Δ)
S.t[2] = w_time(w2, Δ)
S.x[1] = randn(4*nx)
S.x[2] = randn(4*nx)
U = deepcopy(S)
merge!(S)
basic_checks(S)
@test S.n == 1
@test S.t[1] == [1 t0; 8*nx 0]

# ===========================================================================
# (II) as (I) with another channel at a new location
printstyled(stdout,"    channels must have identical :fs, :loc, :resp, and :units to merge \n", color=:light_green)
S = SeisData(mkC1(), mkC2(), prandSC(false), mkC4())

U = deepcopy(S)
merge!(S)
@test S.n == 3
basic_checks(S)
# Are there two channels with the same ID?
@test length(findall(S.id.==id)) == 2

# Is the subgroup merged correctly?
i = findfirst(S.src.=="test channel 2")
@test S.loc[i] == loc
@test mk_tcat(U.t[1:2], fs) == S.t[i]
@test vcat(U.x[1], U.x[2]) == S.x[i]

# Is the odd channel left alone?
j = findfirst(S.src.=="test channel 4")
@test S.loc[j] == loc1
for i in datafields
  if i != :notes
    @test getfield(U, i)[4] == getfield(S, i)[j]
  end
end

# with a second channel at the new location
S = SeisData(deepcopy(U), mkC5())
U = deepcopy(S)
merge!(S)
@test S.n == 3
basic_checks(S)

# Are there two channels with the same ID?
@test length(findall(S.id.==id)) == 2

# Is the first subgroup merged correctly?
i = findfirst(S.src.=="test channel 2")
@test S.loc[i] == loc
@test mk_tcat(U.t[1:2], fs) == S.t[i]
@test vcat(U.x[1], U.x[2]) == S.x[i]

# Is the second subgroup merged correctly?
j = findfirst(S.src.=="test channel 5")
@test S.loc[j] == loc1
@test mk_tcat(U.t[4:5], fs) == S.t[j]
@test vcat(U.x[4], U.x[5]) == S.x[j]

# with resp, not loc
S = deepcopy(U)
S.loc[4] = deepcopy(loc)
S.loc[5] = deepcopy(loc)
S.resp[4] = deepcopy(resp2)
S.resp[5] = deepcopy(resp2)

U = deepcopy(S)
merge!(S)
@test S.n == 3
basic_checks(S)

# Are there two channels with the same ID?
@test length(findall(S.id.==id)) == 2

# Is the first subgroup merged correctly?
i = findfirst(S.src.=="test channel 2")
@test S.loc[i] == loc
@test mk_tcat(U.t[1:2], fs) == S.t[i]
@test vcat(U.x[1], U.x[2]) == S.x[i]

# Is the second subgroup merged correctly?
j = findfirst(S.src.=="test channel 5")
@test S.resp[j] == resp2
@test mk_tcat(U.t[4:5], fs) == S.t[j]
@test vcat(U.x[4], U.x[5]) == S.x[j]

# ===========================================================================
printstyled(stdout,"    one merging channel with a time gap\n", color=:light_green)
S = SeisData(mkC1(), mkC2(), prandSC(false))
S.x[2] = rand(2*nx)
S.t[2] = vcat(S.t[2][1:1,:], [nx 2*Δ], [2*nx 0])

# Is the group merged correctly?
U = deepcopy(S)
merge!(S)
basic_checks(S)
@test S.n == 2
i = findid(id, S)
@test U.x[1] == S.x[i][1:nx]
@test U.x[2][1:nx] == S.x[i][nx+1:2nx]
@test U.x[2][nx+1:2nx] == S.x[i][2nx+1:3nx]
@test mk_tcat(U.t[1:2], fs) == S.t[i]

printstyled(stdout,"      merge window is NOT the first\n", color=:light_green)
os = 2
S = SeisData(mkC1(), mkC2(), prandSC(false))
S.x[1] = rand(2*nx)
S.t[1] = vcat(S.t[1][1:1,:], [nx os*Δ], [2*nx 0])
S.t[2][1,2] += (os+nx)*Δ
U = deepcopy(S)
merge!(S)
basic_checks(S)
@test S.n == 2
i = findid(id, S)
@test U.x[1] == S.x[i][1:2nx]
@test U.x[2] == S.x[i][2nx+1:3nx]
@test mk_tcat(U.t[1:2], fs) == S.t[i]

# ===========================================================================
printstyled(stdout,"    one merge group has non-duplication time overlap\n", color=:light_green)
printstyled(stdout,"      check for averaging\n", color=:light_green)
nov = 3
S = SeisData(mkC1(), C2_ov(), prandSC(false))

# Is the group merged correctly?
U = deepcopy(S)
merge!(S)
basic_checks(S)

@test S.n == 2
i = findid(id, S)
@test S.x[i][1:nx-nov] == U.x[1][1:nx-nov]
@test S.x[i][nx-nov+1:nx] == 0.5*(U.x[1][nx-nov+1:nx] + U.x[2][1:nov])
@test S.x[i][nx+1:2nx-nov] == U.x[2][nov+1:nx]
@test S.t[i] == [1 U.t[1][1,2]; 2nx-nov 0]

# Do the notes log the extra source?
@test findfirst([occursin("New channel 1", i) for i in S.notes[i]]) != nothing
@test findfirst([occursin("test channel 1", i) for i in S.notes[i]]) != nothing

# Are dictionaries merging correctly?
@test haskey(S.misc[i], "P")
@test haskey(S.misc[i], "S")

printstyled(stdout,"      src overlap window is NOT first\n", color=:light_green)
os = 2
S = SeisData(mkC1(), C2_ov(), prandSC(false))
S.x[1] = rand(2*nx)
S.t[1] = vcat(S.t[1][1:1,:], [nx os*Δ], [2*nx 0])
S.t[2][1,2] += (os+nx)*Δ
U = deepcopy(S)
merge!(S)
basic_checks(S)
@test S.n == 2
i = findid(id, S)
@test S.x[i][1:2nx-nov] == U.x[1][1:2nx-nov]
@test S.x[i][2nx-nov+1:2nx] == 0.5*(U.x[1][2nx-nov+1:2nx] + U.x[2][1:nov])
@test S.x[i][2nx+1:3nx-nov] == U.x[2][nov+1:nx]
@test S.t[i] == vcat(U.t[1][1:2,:], [length(S.x[i]) 0])

printstyled(stdout,"      dest overlap window is NOT first\n", color=:light_green)
nov = 3
S = SeisData(mkC1(), C2_ov(), prandSC(false))
S.x[2] = rand(2*nx)
S.t[2] = [1 t0-nx*Δ; nx+1 Δ*(nx-nov); 2*nx 0]

U = deepcopy(S)
merge!(S)
basic_checks(S)
@test S.n == 2
i = findid(id, S)
@test S.x[i][1:nx] == U.x[2][1:nx]
@test S.x[i][nx+1:2nx-nov] == U.x[1][1:nx-nov]
@test S.x[i][2nx-nov+1:2nx] == 0.5*(U.x[1][nx-nov+1:nx] + U.x[2][nx+1:nx+nov])
@test S.x[i][2nx+1:3nx-nov] == U.x[2][nx+nov+1:2nx]
@test mk_tcat(U.t[1:2], fs) == S.t[i]

# ===========================================================================
printstyled(stdout,"    overlap with time mismatch\n", color=:light_green)

#= mk_tcat stops working here as merge shifts one window back in
time one sample to account for the intentional one-Δ time mismatch =#

printstyled(stdout,"      data overlap one sample off of S.t\n", color=:light_green)

# (a) 3_sample overlap with wrong time (C2[1:2] == C1[99:100])
nov = 2
S = SeisData(mkC1(), C2_ov())
S.x[2] = vcat(copy(S.x[1][nx-nov+1:nx]), rand(nx-nov))
U = deepcopy(S)
merge!(S)
basic_checks(S)
@test S.n == 1
i = findid(id, S)
@test length(S.x[i]) == 2nx-nov
@test S.x[i][1:nx-nov] == U.x[1][1:nx-nov]
@test S.x[i][nx-nov+1:nx] == U.x[1][nx-nov+1:nx] == U.x[2][1:nov]
@test S.x[i][nx+1:2nx-nov] == U.x[2][nov+1:nx]
@test S.t[i] == [1 U.t[1][1,2]-Δ; 2nx-nov 0]

printstyled(stdout,"      src overlap window is NOT first\n", color=:light_green)
S = deepcopy(U)
S.x[1] = rand(2*nx)
S.t[1] = vcat(S.t[1][1:1,:], [nx os*Δ], [2*nx 0])
S.x[2] = vcat(copy(S.x[1][2nx-nov+1:2nx]), rand(nx-nov))
S.t[2][1,2] += (os+nx)*Δ
U = deepcopy(S)
merge!(S)
basic_checks(S)
@test S.n == 1
i = findid(id, S)
@test S.x[i][1:2nx-nov] == U.x[1][1:2nx-nov]
@test S.x[i][2nx-nov+1:2nx] == U.x[1][2nx-nov+1:2nx] == U.x[2][1:nov]
@test S.x[i][2nx+1:3nx-nov] == U.x[2][nov+1:nx]
@test S.t[i] == [1 0; U.t[1][2,1] U.t[1][2,2]-Δ; 3nx-nov 0]

#= mk_tcat starts working again here as the time shift is now
applied to the second window, rather than the first. =#

printstyled(stdout,"      dest overlap window is NOT first\n", color=:light_green)
S = SeisData(mkC1(), C2_ov())
S.t[2] = [1 t0-nx*Δ; nx+1 Δ*(nx-nov); 2*nx 0]
S.x[2] = vcat(randn(nx), copy(S.x[1][nx-nov+1:nx]), randn(nx-nov))
U = deepcopy(S)
merge!(S)
basic_checks(S)
@test S.n == 1
i = findid(id, S)
@test S.x[i][1:nx] == U.x[2][1:nx]
@test S.x[i][nx+1:2nx-nov] == U.x[1][1:nx-nov]
@test S.x[i][2nx-nov+1:2nx] == U.x[1][nx-nov+1:nx] == U.x[2][nx+1:nx+nov]
@test S.x[i][2nx+1:3nx-nov] == U.x[2][nx+nov+1:2nx]
@test mk_tcat(U.t[1:2], fs) == S.t[i]

# ===========================================================================
printstyled(stdout,"    multichannel merge with overlap\n", color=:light_green)
breakpt_1 = nx-nov        # 97
breakpt_2 = 2*(nx-nov)-1  # 195
breakpt_3 = 2*(nx-nov)+2  # 198

S = SeisData(mkC1(), C2_ov(), C3_ov(), prandSC(false))
S.x[2] = vcat(copy(S.x[1][nx-nov+1:nx]), rand(nx-nov))
U = deepcopy(S)
merge!(S)
basic_checks(S)
i = findid(id, S)
@test S.n == 2
@test S.x[i][1:breakpt_1] == U.x[1][1:nx-nov]
@test S.x[i][breakpt_1+1:nx] == U.x[1][nx-nov+1:nx] == U.x[2][1:nov]
@test S.x[i][nx+1:breakpt_2] == U.x[2][nov+1:nx-nov-1]
@test S.x[i][breakpt_2+1:breakpt_3] == 0.5*(U.x[2][nx-nov:nx] + U.x[3][1:nov+1])
@test S.x[i][breakpt_3+1:end] == U.x[3][nov+2:end]
@test S.t[i] == [1 U.t[1][1,2]-Δ; 3nx-2nov-1 0]

# What happens when there's a gap in one trace?
breakpt_4 = length(S.x[i])
S = deepcopy(U)
S.x[3] = rand(2*nx)
S.t[3] = vcat(S.t[3][1:1,:], [nx 2*Δ], [2*nx 0])
@test w_time(t_win(S.t[3], Δ), Δ)  == S.t[3]
U = deepcopy(S)
merge!(S)
basic_checks(S)
i = findid(id, S)
@test S.t[i] == [1 t0-Δ; 295 2*Δ; 395 0]
@test S.x[i][1:breakpt_1] == U.x[1][1:nx-nov]
@test S.x[i][breakpt_1+1:nx] == U.x[1][nx-nov+1:nx] == U.x[2][1:nov]
@test S.x[i][nx+1:breakpt_2] == U.x[2][nov+1:nx-nov-1]
@test S.x[i][breakpt_2+1:breakpt_3] == 0.5*(U.x[2][nx-nov:nx] + U.x[3][1:nov+1])
@test S.x[i][breakpt_3+1:breakpt_4] == U.x[3][nov+2:nx]
@test S.x[i][breakpt_4+1:end] == U.x[3][nx+1:end]

# ===========================================================================
printstyled(stdout,"    \"zipper\" merges II\n", color=:light_green)
printstyled(stdout,"      two traces with staggered time windows, some with overlap\n", color=:light_green)

# (a) One overlap in a late window should not shift back previous windows
nov = 2
S = SeisData(mkC1(), mkC2())
W = Array{Int64,2}(undef, 8, 2);
for i = 1:8
  W[i,:] = t0 .+ [(i-1)*nx*Δ ((i-1)*nx + nx-1)*Δ]
end
W[8,1] -= nov*Δ
W[8,2] -= nov*Δ
w1 = W[[1,3,5,7], :]
w2 = W[[2,4,6,8], :]
S.t[1] = w_time(w1, Δ)
S.t[2] = w_time(w2, Δ)
S.x[1] = randn(4*nx)
S.x[2] = randn(4*nx)
U = deepcopy(S)
merge!(S)
basic_checks(S)
i = findid(id, S)
@test S.n == 1
@test S.t[i] == [1 t0; 8*nx-nov 0]

# These should be untouched
for j = 1:6
  si = 1 + nx*(j-1)
  ei = nx*j
  k = div(j,2)
  if isodd(j)
    usi = k*nx + 1
    uei = usi + nx - 1
    @test S.x[i][si:ei] == U.x[1][usi:uei]
  else
    usi = (k-1)*nx + 1
    uei = usi + nx - 1
    @test S.x[i][si:ei] == U.x[2][usi:uei]
  end
end

# The only overlap should be here:
@test S.x[i][1+6nx:7nx-nov] == U.x[1][1+3nx:4nx-nov]
@test S.x[i][7nx-nov+1:7nx] == 0.5*(U.x[1][4nx-nov+1:4nx] .+ U.x[2][3nx+1:3nx+nov])
@test S.x[i][7nx+1:8nx-nov] == U.x[2][3nx+nov+1:4nx]

printstyled(stdout,"      one overlap, late window, time mismatch\n", color=:light_green)
nov = 3
true_nov = 2
S = SeisData(mkC1(), mkC2())
W = Array{Int64,2}(undef, 8, 2);
for i = 1:8
  W[i,:] = t0 .+ [(i-1)*nx*Δ ((i-1)*nx + nx-1)*Δ]
end
W[8,1] -= nov*Δ
W[8,2] -= nov*Δ
w1 = W[[1,3,5,7], :]
w2 = W[[2,4,6,8], :]
S.t[1] = w_time(w1, Δ)
S.t[2] = w_time(w2, Δ)
S.x[1] = randn(4*nx)
S.x[2] = randn(4*nx)
S.x[2][3nx+1:3nx+2] = copy(S.x[1][4nx-1:4nx])
U = deepcopy(S)
merge!(S)
basic_checks(S)
@test S.n == 1
for j = 1:5
  si = 1 + nx*(j-1)
  ei = nx*j
  k = div(j,2)
  if isodd(j)
    usi = k*nx + 1
    uei = usi + nx - 1
    @test S.x[i][si:ei] == U.x[1][usi:uei]
  else
    usi = (k-1)*nx + 1
    uei = usi + nx - 1
    @test S.x[i][si:ei] == U.x[2][usi:uei]
  end
end
δi = nov - true_nov

# Because we moved the start time of 2:3 back by 1, we expect:
@test S.x[i][1+5nx:6nx-δi] == U.x[2][1+2nx:3nx-δi]

# Thus there's a one-point overlap (that gets resolved by averaging) at:
@test S.x[i][6nx-δi+1:6nx] == 0.5*(U.x[2][3nx-δi+1:3nx] + U.x[1][3nx+1:3nx+δi])

# Proceeding through the merged time series, we expect:
@test S.x[i][6nx+1:7nx-nov] == U.x[1][3nx+δi+1:4nx-true_nov]
@test S.x[i][7nx-nov+1:7nx-δi] == U.x[1][4nx-true_nov+1:4nx] == U.x[2][3nx+1:3nx+true_nov]
@test S.x[i][7nx-δi+1:end] == U.x[2][3nx+true_nov+1:end]

# ============================================================================
printstyled(stdout,"    distributivity: S1*S3 + S2*S3 == (S1+S2)*S3\n", color=:light_green)
imax = 10
printstyled("      trial ", color=:light_green)
for i = 1:imax
  if i > 1
    print("\b\b\b\b\b")
  end
  printstyled(string(lpad(i, 2), "/", imax), color=:light_green)
  S1 = randSeisData()
  S2 = randSeisData()
  S3 = randSeisData()
  # M1 = (S1+S2)*S3
  # M2 = S1*S3 + S2*S3
  @test ((S1+S2)*S3) == (S1*S3 + S2*S3)
  if i == imax
    println("")
  end
end

# ============================================================================
printstyled(stdout,"    checking (formerly-breaking) end-member cases\n", color=:light_green)
printstyled(stdout,"      time windows not in chronological order\n", color=:light_green)
C1 = mkC1()
C2 = mkC2()
C3 = deepcopy(C1)
C3.t = [1 0; 101 -2000000; 200 0]
append!(C3.x, randn(100))
S = SeisData(C3, C2, prandSC(false))

# Is the group merged correctly?
U = deepcopy(S)
merge!(S)
basic_checks(S)
@test S.n == 2
i = findid(id, S)
@test S.t[i] == [1 -1000000; 300 0]
@test vcat(U.x[1][101:200], U.x[1][1:100], U.x[2]) == S.x[i]
@test mk_tcat(U.t[1:2], fs) == S.t[i]

# Do the notes log the extra source?
@test findfirst([occursin("New channel 1", j) for j in S.notes[i]]) != nothing
@test findfirst([occursin("test channel 1", j) for j in S.notes[i]]) != nothing

# Are dictionaries merging correctly?
@test haskey(S.misc[i], "P")
@test haskey(S.misc[i], "S")

printstyled(stdout,"      sequential one-sample windows\n", color=:light_green)
C1 = mkC1()
C2 = mkC2()
C2.t = [1 1200000; 99 90059; 100 90210]
C2.t = [1 1000000; 99 90059; 100 90210]
S = SeisData(C1, C2, prandSC(false))

# Is the group merged correctly?
U = deepcopy(S)
merge!(S)
basic_checks(S)
@test S.n == 2
i = findid(id, S)
@test S.t[i] == [1  0; 199  90059; 200  90210]
@test vcat(U.x[1], U.x[2])==S.x[i]
@test mk_tcat(U.t[1:2], fs) == S.t[i]

# Do the notes log the extra source?
@test findfirst([occursin("New channel 1", j) for j in S.notes[i]]) != nothing
@test findfirst([occursin("test channel 1", j) for j in S.notes[i]]) != nothing

# Are dictionaries merging correctly?
@test haskey(S.misc[i], "P")
@test haskey(S.misc[i], "S")

# ============================================================================
printstyled(stdout,"  merge! and new/extended methods\n", color=:light_green)
fs1 = 50.0
Δ = round(Int64, sμ/fs1)

# t1 = round(Int64,time()/μs)
t1 = 0

# s1 and s2 represent data from a fictitious channel
# s2 begins 1 second later, but has a gap of 1s after sample 25
s1 = SeisChannel(fs = fs1, gain = 10.0, name = "DEAD.STA.EHZ", id = "DEAD.STA..EHZ",
  t = [1 t1; 100 0], x=randn(100))
s2 = SeisChannel(fs = fs1, gain = 10.0, name = "DEAD.STA.EHZ", id = "DEAD.STA..EHZ",
  t = [1 t1+1000000; 26 1000000; 126 200000; 150 0], x=randn(150))

s3 = SeisChannel(fs = 100.0, gain = 5.0, name = "DEAD.STA.EHE", id = "DEAD.STA..EHE",
  t = [1 t1; 100 0], x=rand(100).-0.5)
s4 = SeisChannel(fs = 100.0, gain = 50.0, name = "DEAD.STA.EHE", id = "DEAD.STA..EHE",
  t = [1 t1+1000000; 100 1000000; 150 0], x=randn(150))

# We expect:
# (1) LAST 25-26 points in s1 will be averaged with first 25 points in s3
# (2) s3 will be ungapped with exactly 50 samples of junk
# (3) after sync of a seisobj formed from [s1+s3, s2], s2 will be padded with 0.5s
# at start and 3s at end

S = SeisData()
for i in fieldnames(typeof(S))
  if i != :n
    @test ≈(isempty([]), isempty(getfield(S,i)))
  end
end

printstyled("    fixing data gaps\n", color=:light_green)
s2u = ungap(s2)
@test length(s2.x) / s2.fs + μs * sum(s2.t[2:end - 1, 2]) ≈ length(s2u.x) / s2u.fs

printstyled("    channel add and simple merges\n", color=:light_green)
S = S + (s1 * s2)

# Do in-place operations only change leftmost variable?
@test length(s1.x) ≈ 100
@test length(s2.x) ≈ 150

# expected behavior for S = s1+s2
# (0) length[s.x[1]] = 250
# (1) S.x[1][1:50] = s1.x[1:50]
# (2) S.x[1][51:75] = 0.5.*(s1.x[51:75] + s2.x[1:25])
# (3) S.x[1][76:100] = s1.x[76:100]
# (4) S.x[1][101:125] = mean(S.x[1])
# (5) S.x[1][126:250] = s2.x[126:250]

# Basic merge ops
@test S.n == 1
@test length(S.fs) == 1
@test ==(S.fs[1], s1.fs)
@test ==(length(S.gain), 1)
@test ==(S.gain[1], s1.gain)
@test ==(length(S.name), 1)
@test ==(length(S.t),1)
@test ==(length(S.x),1)
@test S.id[1]=="DEAD.STA..EHZ"
U = deepcopy(S)
ungap!(S, m=false, tap=false) # why do I have to force type here
@test ==(length(S.x[1]), 260)
@test ==(S.x[1][1:50], s1.x[1:50])
@test ==(S.x[1][51:75], 0.5.*(s1.x[51:75] .+ s2.x[1:25]))
@test ==(S.x[1][76:100], s1.x[76:100])
@test minimum(isnan.(S.x[1][101:125]))
@test ==(S.x[1][126:225], s2.x[26:125])
@test minimum(isnan.(S.x[1][226:235]))
@test ==(S.x[1][236:260], s2.x[126:150])

# Auto-tapering after a merge
T = deepcopy(S)
ii = findall(isnan.(S.x[1]))
nanfill!(S)
taper!(S, N_min = 0)
@test length(findall(isnan.(S.x[1])))==0          # No more NaNs?
@test sum(diff(S.x[1][ii]))==0                    # All NaNs filled w/same val?
@test ≈(T.x[1][15:90], S.x[1][15:90])             # Un-windowed vals untouched?

printstyled("    channel delete\n", color=:light_green)
S -= 1
for i in fieldnames(typeof(S))
  if i != :n
    @test ≈(isempty([]), isempty(getfield(S,i)))
  end
end
@test isempty(S)

printstyled("    a more difficult merge\n", color=:light_green)
S = SeisData()
S *= (s1 * s3)
S *= (s2 * s4)
@test ≈(S.n, 2)
@test ≈(length(S.fs), 2)
@test ≈(length(S.gain), 2)
@test ≈(length(S.name), 2)
@test ≈(length(S.t),2)
@test ≈(length(S.x),2)
i = findid("DEAD.STA..EHE", S)
j = findid("DEAD.STA..EHZ", S)
@test ≈(S.fs[i], 100.0)
@test ≈(S.gain[j], 10.0)
ungap!(S, m=false, tap=false)
@test ≈(length(S.x[j]), 260)
@test ≈(length(S.x[i]), 350)
@test ≈(S.x[i][101]/S.gain[i], s4.x[1]/s4.gain)

printstyled("      repeating data\n", color=:light_green)
loc1 = GeoLoc(lat=45.28967, lon=-121.79152, el=1541.0)
loc2 = GeoLoc(lat=48.78384, lon=-121.90093, el=1676.0)
s1 = SeisChannel(fs = fs1, gain = 10.0, name = "DEAD.STA.EHZ", id = "DEAD.STA..EHZ",
t = [1 t1; 100 0], x=randn(100))
s2 = SeisChannel(fs = fs1, gain = 10.0, name = "DEAD.STA.EHZ", id = "DEAD.STA..EHZ",
  t = [1 t1+1000000; 150 0], x=vcat(s1.x[51:100], randn(100)))
C = (s1 * s2)[1]
@test length(C.x) == 200
@test C.x[1:100] == s1.x
@test C.x[101:200] == s2.x[51:150]

printstyled("    overlapping times\n", color=:light_green)
os = 50
nx = length(s1.x)
lx = length(s1.x)/s1.fs
τ = round(Int, sμ*(2.0-os/fs1))
s2 = SeisChannel(fs = fs1, gain = 10.0, name = "DEAD.STA.EHZ", id = "DEAD.STA..EHZ",
  t = [1 t1+τ; 150 0], x=randn(150))
C = (s1 * s2)[1]
@test length(C.x) == 250-os
@test C.x[1:nx-os] == s1.x[1:os]
@test C.x[nx-os+1:nx] == 0.5.*(s1.x[nx-os+1:nx]+s2.x[1:os])
@test C.x[nx+1:200] == s2.x[os+1:150]

printstyled("    two-sample offset\n", color=:light_green)
os = 2
nx = length(s1.x)
lx = length(s1.x)/s1.fs
τ = round(Int, sμ*(2.0-(os-1)/fs1))
s2 = SeisChannel(fs = fs1, gain = 10.0, name = "DEAD.STA.EHZ", id = "DEAD.STA..EHZ",
  t = [1 t1+τ; 150 0], x=vcat(copy(s1.x[nx-os+1:nx]), randn(150-os)))
@test s1.x[nx] == s2.x[os]
C = (s1 * s2)[1]
@test length(C.x) == 250-os+1
@test C.x[1:99] == s1.x[1:99]
@test C.x[100] == 0.5*(s2.x[1] + s1.x[100])
@test C.x[101:249] == s2.x[2:150]

# Ensure merge works correctly with traces separated in time
printstyled("    operator \"*\"\n", color=:light_green)
s5 = SeisChannel(fs = 100.0, loc=loc1, gain = 32.0, name = "DEAD.STA.EHE", id = "DEAD.STA..EHE",
  t = [1 t1; 100 0], x=randn(100))
s6 = SeisChannel(fs = 100.0, loc=loc1, gain = 16.0, name = "UNNAMED", id = "DEAD.STA..EHE",
  t = [1 t1+30000000; 200 0], x=randn(200))
T = (s5 * s6)
ungap!(T)
@test ≈(length(T.x[1]),3200)

printstyled("    one repeating segment in two channels\n", color=:light_green)
fs1 = 50.0
Δ = round(Int64, sμ/fs1)
g1 = 10.0
ns = 2
n1 = "DEAD.STA..EHZ"
id = "DEAD.STA..EHZ"
t_in = [1 0; 300 0]
nx = round(Int64, ns*fs1)
seg_1 = randn(nx)
seg_2 = randn(nx)
seg_3 = randn(nx)
t = t_expand(t_in, fs1)
i1 = vcat(collect(1:100), collect(201:300))
i2 = collect(101:300)
t1 = t_collapse(t[i1], fs1)
t2 = t_collapse(t[i2], fs1)
s1 = SeisChannel(fs = fs1, gain = g1, name = n1, id = id, t = t1, x = vcat(seg_1, seg_3))
s2 = SeisChannel(fs = fs1, gain = g1, name = n1, id = id, t = t2, x = vcat(seg_2, seg_3))
C = (s1 * s2)[1]
@test length(C.x) == 300
@test C.x == vcat(seg_1, seg_2, seg_3)
@test C.t == t_in

printstyled(stdout,"    SeisData * SeisChannel ==> SeisData\n", color=:light_green)
(S,T) = mktestseis()
merge!(S,T)
C = prandSC(false)
D = prandSC(false)
U = merge(S,C)
U = sort(U)
V = merge(C,S)
V = sort(V)
@test U == V
U = merge(C,D)

(S,T) = mktestseis()
A = deepcopy(S[5])
B = deepcopy(T[4])
T*=S[1]
sizetest(T, 5)

printstyled(stdout,"    SeisData merge tests\n", color=:light_green)
printstyled(stdout,"      one common channel (fast, no gaps)\n", color=:light_green)
(S,T) = mktestseis()
t = t_expand(S.t[4], S.fs[4])
T.t[3][1,2] = t[end] + t[end]-t[end-1]
merge!(T, S[4])

printstyled(stdout,"      one common channel (slow, has gaps)\n", color=:light_green)
(S,T) = mktestseis()
A = deepcopy(S[5])
B = deepcopy(T[4])
T*=S[1]
sizetest(T, 5)

T = sort(T)
S*=T[2]
sizetest(S, 5)
i = findid(A, S)
@test ≈(S.t[i][2,1], 1+length(A.x))
@test ≈(S.t[i][2,2], (5-1/S.fs[i])*sμ)

printstyled(stdout,"      common channels and \"splat\" notation (mseis!)\n", color=:light_green)
(S,T) = mktestseis()
U = merge(S,T)
sizetest(U, 7)

V = SeisData(S,T)
merge!(V)
@test U == V

mseis!(S,T)
@test S == V

printstyled(stdout,"      mseis! with Types from SeisIO.Quake\n", color=:light_green)
S = randSeisData()
Ev = randSeisEvent(2)
C = ungap(deepcopy(Ev.data[1]))
Ev.data[1] = deepcopy(C)
C.t[1,2] += 1000000
C.x = randn(eltype(C.x), length(C.x))
C.az = 0.0
C.baz = 0.0
C.dist = 0.0
Ev.data[2] = deepcopy(C)
mseis!(S,   randSeisChannel(),
            convert(EventChannel, randSeisChannel()),
            rand(Float64, 23),  # should warn
            convert(EventTraceData, randSeisData()),
            Ev,
            randSeisEvent())

printstyled(stdout,"      two independent channels ==> same as \"+\"\n", color=:light_green)
(S,T) = mktestseis()
U = S[1] * T[2]
sizetest(U, 2)
@test U == S[1]+T[2]

printstyled(stdout,"      two identical channel IDs\n", color=:light_green)
U = S[4] * T[3]
@test typeof(U)==SeisData
@test U.id[1]==S.id[4]
@test U.id[1]==T.id[3]

printstyled(stdout,"    pull\n", color=:light_green)
ux = Float64.(randperm(800))
C = pull(S,4)
C.x = ux[1:600]
@test C.name=="Longmire"
@test S.n==4
@test length(findall(S.name.=="Longmire"))==0

printstyled(stdout,"      notes still faithfully track names in modified objects\n", color=:light_green)
str1 = "ADGJALMGFLSFMGSLMFLSeptember Lobe sucks"
str2 = "HIGH SNR ON THIS CHANNEL"
note!(S,str2)
note!(S,str1)
@test (findall([maximum([occursin(str1, S.notes[i][j]) for j=1:length(S.notes[i])]) for i = 1:S.n]) == [4])
@test (length(findall([maximum([occursin(str2, S.notes[i][j]) for j = 1:length(S.notes[i])]) for i = 1:S.n]))==S.n)

# merge test: when we merge, does each field have exactly 7 entries?
printstyled(stdout,"      merge! after pull doesn't affect pulled channel\n", color=:light_green)
T.x[3] = ux[601:800]
merge!(S,T)
n_targ = 7
@test ≈(S.n, n_targ)
@test ≈(maximum([length(getfield(S,i)) for i in datafields]), n_targ)
@test ≈(minimum([length(getfield(S,i)) for i in datafields]), n_targ)
@test any(maximum([C.x.==i for i in S.x[1]])) == false

# Old: (Test Passed, 1.42978256,  154864097, 0.038663895, Base.GC_Diff(154864097, 109, 0, 1845386, 4165, 0,  38663895, 7, 0))
# New: (Test Passed, 1.263168574, 128490661, 0.108295874, Base.GC_Diff(128490661,  81, 0, 1324714, 3857, 0, 108295874, 6, 1))

printstyled(stdout,"    purge!\n", color=:light_green)
(S,T) = mktestseis()
S.t[5] = Array{Int64,2}(undef,0,2)
S.x[5] = Array{Float32,1}(undef,0)
U = purge(S)
purge!(S)
@test S == U
@test S.n == 4

printstyled(stdout,"    issue 29\n", color=:light_green)
S = read_data(path * "/SampleFiles/SEED/CIRIO__BHE___2017101.mseed")
merge!(S)
C2 = S[1]
ungap!(C2)

#=  In Issue #29, ObsPy handled the merge without repeated data, so this is the
    basis for comparison
=#
io = open(path * "/SampleFiles/SEED/obspy.dat", "r")
X2 = Array{Int32,1}(undef, length(C2.x))
read!(io, X2)
close(io)
X = map(Float32, X2)
@test isapprox(C2.x, X)

printstyled(stdout,"      merge_ext! on SeisData == nothing\n", color=:light_green)
C = randSeisChannel(s=true)
S = randSeisData(3, s=1.0)
S.id = [C.id, C.id, C.id]
S.fs = [C.fs, C.fs, C.fs]
U = deepcopy(S)
Ω = 3
rest = [1,2]
@test SeisIO.merge_ext!(S, Ω, rest) == nothing
@test S==U


text:
module UnitSystems

#   This file is part of UnitSystems.jl. It is licensed under the MIT license
#   UnitSystems Copyright (C) 2020 Michael Reed

import Base: @pure, length, time

const Systems = (:Metric,:SI2019,:CODATA,:Conventional,:MTS,:English,:EnglishUS,:IAU,:SI1976,:Mixed,:ESU2019,:EMU2019,:EMU,:ESU,:Gauss,:LorentzHeaviside,:Thomson,:Kennelly,:Planck,:PlanckGauss,:Stoney,:Hartree,:Rydberg,:Schrodinger,:Electronic,:Natural,:NaturalGauss,:QCD,:QCDGauss,:QCDoriginal)
const Constants = (:hyperfine,:lightspeed,:planck,:planckreduced,:electronmass,:molarmass,:boltzmann,:permeability,:rationalization,:lorentz,:luminousefficacy)
const Physics = (:atomicmass,:protonmass,:planckmass,:newton,:einstein,:hartree,:rydberg,:bohr,:bohrreduced,:electronradius,:avogadro,:universal,:stefan,:radiationdensity,:permittivity,:coulomb,:ampere,:biotsavart,:charge,:faraday,:impedance,:conductance,:klitzing,:josephson,:magneticflux,:magneton)
const Kinematic = (:time,:length,:area,:volume,:wavenumber,:fuelefficiency,:frequency,:frequencydrift,:speed,:acceleration,:jerk,:snap,:volumeflow)
const Mechanical = (:mass,:massflow,:lineardensity,:areadensity,:density,:specificvolume,:force,:stiffness,:pressure,:compressibility,:viscosity,:diffusivity,:rotationalinertia,:momentum,:angularmomentum,:yank,:energy,:specificenergy,:action,:fluence,:power,:powerdensity,:intensity,:spectralflux,:soundexposure,:impedance,:specificimpedance,:admittance,:compliance,:inertance)
const Electromagnetic = (:charge,:chargedensity,:linearchargedensity,:exposure,:mobility,:current,:currentdensity,:resistance,:conductance,:resistivity,:conductivity,:capacitance,:inductance,:reluctance,:permeance,:permittivity,:permeability,:susceptibility,:specificsusceptibility,:demagnetizingfactor,:vectorpotential,:electricpotential,:magneticpotential,:electricfield,:magneticfield,:electricflux,:magneticflux,:electricfluxdensity,:magneticfluxdensity,:electricdipolemoment,:magneticdipolemoment,:electricpolarizability,:magneticpolarizability,:magneticmoment,:magnetizability,:magnetization,:specificmagnetization,:rigidity,:polestrength)
const Thermodynamic = (:temperature,:entropy,:specificentropy,:volumeheatcapacity,:thermalconductivity,:thermalconductance,:thermalresistance,:thermalexpansion,:lapserate)
const Molar = (:molarmass,:molality,:mole,:molarity,:molarvolume,:molarentropy,:molarenergy,:molarconductivity,:molarsusceptibility,:catalysis,:specificity)
const Photometric = (:luminousflux,:luminance,:luminousenergy,:luminousexposure,:luminousefficacy)
const Mechanics = [Kinematic...,Mechanical...]
const Convert = [Mechanics...,Electromagnetic...,Thermodynamic...,Molar...,Photometric...]

listext(x) = join(x,"`, `")

@pure measure(x) = x

# universe

"""
    Coupling{αG,α,μₑᵤ,μₚᵤ}

Specification of `Universe` with the dimensionless `Coupling` constants `coupling`, `finestructure`, `electronunit`, `protonunit`, and `protonelectron`.
Alterations to these values can be facilitated and quantified using parametric polymorphism.
Due to the `Coupling` interoperability, the `MeasureSystems` package is made possible to support calculations with `Measurements` having error standard deviations.
"""
struct Coupling{αG,α,μₑᵤ,μₚᵤ} end
@pure coupling(U::Coupling{αG}) where αG = measure(αG)
@pure finestructure(U::Coupling{αG,α}) where {αG,α} = measure(α)
@pure electronunit(U::Coupling{αG,α,μₑᵤ}) where {αG,α,μₑᵤ} = measure(μₑᵤ)
@pure protonunit(U::Coupling{αG,α,μₑᵤ,μₚᵤ}) where {αG,α,μₑᵤ,μₚᵤ} = measure(μₚᵤ)
@pure protonelectron(U::Coupling) = protonunit(U)/electronunit(U)
Base.display(U::Coupling) = println("Coupling{αG=$(coupling(U)),α=$(finestructure(U)),μₑᵤ=$(electronunit(U)),μₚᵤ=$(protonunit(U))}")

# unit systems

"""
    UnitSystem{kB,ħ,𝘤,μ₀,mₑ,λ,αL}

Fundamental constants of physics are: `kB` Boltzmann's constant, `ħ` reduced Planck's constant, `𝘤` speed of light, `μ₀` vacuum permeability, `mₑ` electron rest mass, `λ` Gauss rationalization, and `αL` Lorentz's constant.
Primarily the `Metric` SI unit system is used in addition to the historic `English` engineering unit system.
These constants induce derived values for `avogadro`, `boltzmann`, `universal`, `planck`, `planckreduced`, `lightspeed`, `planckmass`, `atomicmass`, `protonmass`, `electronmass`, `newton`, `einstein`, `permeability`, `permittivity`, `coulomb`, and
additional constants `molarmass`, `hyperfine`, `luminousefficacy`, `stefan`, `radiationdensity`, `ampere`, `lorentz`, `biotsavart`, `rationalization`, `impedance`, `charge`, `magneton`, `conductance`, `faraday`, `magneticflux`, `josephson`, `klitzing`, `hartree`, `rydberg`, `bohr`, and `bohrreduced`.

Additional reference `UnitSystem` variants: `EMU`, `ESU`, `Gauss`, `LorentzHeaviside`, `MTS`, `SI2019`, `CODATA`, `Conventional`, `IAU`, `EnglishUS`; and natural atomic units based on gravitational coupling `αG` and the fine structure `1/αinv` constant (`Planck`, `PlanckGauss`, `Stoney`, `Hartree`, `Rydberg`, `Schrodinger`, `Electronic`, `Natural`, `NaturalGauss`, `QCD`, `QCDGauss`, and `QCDoriginal`).

**Derived unit conversions:**

Mechanics: `$(listext(Kinematic))`, `$(listext(Mechanical))`;
Electromagnetics: `$(listext(Electromagnetic))`;
Thermodynamics: `$(listext(Thermodynamic))`,
`$(listext(Molar))`, `$(listext(Photometric))`.
""" #`Rᵤ,mᵤ,σ,ħ,μ₀,ε₀,kₑ,𝘦,𝔉,RK,Z₀,G₀`
struct UnitSystem{kB,ħ,𝘤,μ₀,mₑ,λ,αL} end
@pure UnitSystem{k,ħ,𝘤,μ,m,λ}() where {k,ħ,𝘤,μ,m,λ} = UnitSystem{k,ħ,𝘤,μ,m,λ,1}()
@pure UnitSystem{k,ħ,𝘤,μ,m}() where {k,ħ,𝘤,μ,m} = UnitSystem{k,ħ,𝘤,μ,m,1}()
@pure boltzmann(::UnitSystem{k}) where k = measure(k)
@pure planckreduced(::UnitSystem{k,ħ}) where {k,ħ} = measure(ħ)
@pure lightspeed(::UnitSystem{k,ħ,𝘤}) where {k,ħ,𝘤} = measure(𝘤)
@pure permeability(::UnitSystem{k,ħ,𝘤,μ}) where {k,ħ,𝘤,μ} = measure(μ)
@pure electronmass(::UnitSystem{k,ħ,𝘤,μ,m}) where {k,ħ,𝘤,μ,m} = measure(m)
@pure rationalization(::UnitSystem{k,ħ,𝘤,μ,m,λ}) where {k,ħ,𝘤,μ,m,λ} = measure(λ)
@pure lorentz(::UnitSystem{k,ħ,𝘤,μ,m,λ,α}) where {k,ħ,𝘤,μ,m,λ,α} = measure(α)
# ΔνCs:s⁻¹, c:m⋅s⁻¹, h:kg⋅m²⋅s⁻¹, kB:kg⋅m²⋅s⁻²⋅K⁻¹, NA:mol⁻¹, Kcd: cd⋅sr⋅s³⋅kg⁻¹⋅m⁻²

isrationalized(U::UnitSystem) = rationalization(U) ≠ 4π

Base.display(U::UnitSystem) = println("UnitSystem{kB=$(boltzmann(U)),ħ=$(planckreduced(U)),𝘤=$(lightspeed(U)),μ₀=$(permeability(U)),mᵤ=$(electronmass(U)),λ=$(isrationalized(U) ? rationalization(U) : "4π"),αL=$(lorentz(U))}")

@pure universe(::UnitSystem) = Universe
@pure unit(x,y=1) = isapprox(y,x,rtol=eps()^0.9) ? y : x
@pure mass(U::UnitSystem,S::UnitSystem) = electronmass(U,S)
@pure electronmass(𝘩::Float64) = αinv^2*R∞*2𝘩/𝘤
@pure electronmass(𝘩::Float64,C::Coupling) = inv(finestructure(C))^2*R∞*2𝘩/𝘤
@pure planckmass(U::UnitSystem,C::Coupling=universe(U)) = electronmass(U,C)/√coupling(C)
@pure planck(U::UnitSystem,C::Coupling=universe(U)) = 2π*planckreduced(U,C)
@pure newton(U::UnitSystem,C::Coupling=universe(U)) = lightspeed(U,C)*planckreduced(U,C)/planckmass(U,C)^2
@pure charge(U::UnitSystem,C::Coupling=universe(U)) = sqrt(2planck(U)/(permeability(U)/finestructure(U))/(lightspeed(U)*rationalization(U)*lorentz(U)^2))

for unit ∈ (:coupling,:finestructure,:electronunit,:protonunit,:protonelectron)
    @eval @pure $unit(U::UnitSystem) = $unit(universe(U))
end
for unit ∈ (:boltzmann,:planckreduced,:lightspeed,:permeability,:electronmass,:molarmass)
    @eval @pure $unit(U::UnitSystem,C::Coupling) = $unit(U)
end
for unit ∈ Constants
    @eval @pure $unit(U::UnitSystem,S::UnitSystem) = unit($unit(S)/$unit(U))
end
for unit ∈ Convert
    @eval begin
        @pure @inline $unit(v::Real,U::UnitSystem) = $unit(v,U,Metric)
        @pure @inline $unit(v::Real,U::UnitSystem,S::UnitSystem) = (u=$unit(U,S);isone(u) ? v : v/u)
        @pure @inline $unit(v::Real,U::UnitSystem{kB,ħ,𝘤,μ₀,mₑ},S::UnitSystem{kB,ħ,𝘤,μ₀,mₑ}) where {kB,ħ,𝘤,μ₀,mₑ} = v
        @pure @inline $unit(U::UnitSystem{kB,ħ,𝘤,μ₀,mₑ},S::UnitSystem{kB,ħ,𝘤,μ₀,mₑ}) where {kB,ħ,𝘤,μ₀,mₑ} = 1
    end
    if unit ∉ (Constants...,:permittivity,:charge,:magneticflux,:impedance,:conductance)
        @eval @pure @inline $unit(U::UnitSystem) = $unit(U,Metric)
    end
end
for unit ∈ (Systems...,Constants...,Physics...,Convert...)
    @eval export $unit
end

function (U::UnitSystem)(JK=1,Js=1,ms=1,Hm=1,kg=1)
    kB = boltzmann(U)*JK
    ħ = planckreduced(U)*Js
    𝘤 = lightspeed(U)*ms
    μ₀ = permeability(U)*Hm
    mₑ = electronmass(U)*kg
    λ = rationalization(U)
    αL = lorentz(U)
    UnitSystem{kB,ħ,𝘤,μ₀,mₑ,λ,isone(αL) ? αL : αL/ms}()
end

# common conversion factors

const g₀,ft,ftUS = 9.80665,0.3048,1200/3937; const lbm,lbmUS = g₀/ft,g₀/ftUS
const slug,slugUS,rankine,kelvin,atm = 0.45359237lbm,0.45359237lbmUS,5/9,9/5,101325.0
const kcalₜₕ,kcal₄,kcal₁₀,kcal₂₀,kcalₘ,kcalᵢₜ = 4184,4204,4185.5,4182,4190,4186.8
const calₜₕ,cal₄,cal₁₀,cal₂₀,calₘ,calᵢₜ = (kcalₜₕ,kcal₄,kcal₁₀,kcal₂₀,kcalₘ,kcalᵢₜ)./1e3
const kcal = kcalₜₕ; const cal = kcal/1000 # calₜₕ thermal calorie

# fundamental constants, αinv = (34259-1/4366.8123)/250 # 137.036 exactly?

const ΔνCs,Kcd,mP = 9192631770.0,683*555.016/555,2.176434e-8 # planck mass (kg)
const NA,kB,𝘩,𝘤,𝘦 = 6.02214076e23,1.380649e-23,6.62607015e-34,299792458.,1.602176634e-19
const μₑᵤ,μₚᵤ,αinv,R∞ = 1/1822.888486209,1.007276466621,137.035999084,10973731.5681601
const mₑ,μ₀ = electronmass(𝘩),2𝘩/𝘤/αinv/𝘦^2 # ≈ 4π*(1e-7+5.5e-17), exact charge
const ħ,δμ₀,μₚₑ,Rᵤ,αL,αG = 𝘩/2π,μ₀-4π*1e-7,μₚᵤ/μₑᵤ,NA*kB,0.01/𝘤,(mₑ/mP)^2
const RK1990,RK2014,KJ1990,KJ2014 = 25812.807,25812.8074555,4.835979e14,4.835978525e14
const ħ1990,ħ2014 = 2/RK1990/KJ1990^2/π,2/RK2014/KJ2014^2/π
const mₑ1990,mₑ2014 = electronmass(2π*ħ1990),electronmass(2π*ħ2014)

# engineering units # Thomson: αL = 1/2

const Universe = Coupling{αG,1/αinv,μₑᵤ,μₚᵤ}()
const Gauss = UnitSystem{1e10*Rᵤ*mₑ/μₑᵤ,1e7*ħ,100𝘤,1,1000mₑ,4π,0.01/𝘤}()
const LorentzHeaviside = UnitSystem{1e10*Rᵤ*mₑ/μₑᵤ,1e7*ħ,100𝘤,1,1000mₑ,1,0.01/𝘤}()
const Thomson = UnitSystem{1e10*Rᵤ*mₑ/μₑᵤ,1e7*ħ,100𝘤,1,1000mₑ,4π,1/2}()
const Kennelly = UnitSystem{Rᵤ*mₑ/μₑᵤ/0.001,ħ,𝘤,1e-7,mₑ,4π}() # ?
const ESU = UnitSystem{1e10*Rᵤ*mₑ/μₑᵤ,1e7*ħ,100𝘤,(100𝘤)^-2,1000mₑ,4π}()
const ESU2019 = UnitSystem{1e7*kB,1e7*ħ,100𝘤,1e3*μ₀/𝘤^2,1000mₑ}()
const EMU = UnitSystem{1e10*Rᵤ*mₑ/μₑᵤ,1e7*ħ,100𝘤,1,1000mₑ,4π}()
const EMU2019 = UnitSystem{1e7*kB,1e7*ħ,100𝘤,1e7*μ₀,1000mₑ}()
const MTS = UnitSystem{1e6*Rᵤ*mₑ/μₑᵤ,1000ħ,𝘤,4π/1e4,mₑ/1000}()
const Mixed = UnitSystem{Rᵤ*mₑ/μₑᵤ/0.001,ħ,𝘤,μ₀,mₑ}()
const Metric = UnitSystem{Rᵤ*mₑ/μₑᵤ/0.001,ħ,𝘤,4π*1e-7,mₑ}()
const SI1976 = UnitSystem{8.31432mₑ/μₑᵤ/0.001,ħ,𝘤,4π*1e-7,mₑ}()
const SI2019 = UnitSystem{kB,ħ,𝘤,μ₀,mₑ}()
const CODATA = UnitSystem{Rᵤ*mₑ2014/μₑᵤ/0.001,ħ2014,𝘤,2RK2014/𝘤/αinv,mₑ2014}()
const Conventional = UnitSystem{Rᵤ*mₑ1990/μₑᵤ/0.001,ħ1990,𝘤,2RK1990/𝘤/αinv,mₑ1990}()
const English = UnitSystem{kB*rankine/slug/ft^2,ħ/slug/ft^2,𝘤/ft,4π,mₑ/slug}()
const EnglishUS = UnitSystem{1000Rᵤ*mₑ/μₑᵤ*rankine/slug/ftUS^2,ħ/slug/ftUS^2,𝘤/ftUS,4π,mₑ/slug}()

# astronomical units

const GMsun,GMearth,GMjupiter = 1.32712442099e20,398600441.8e6,1.26686534e17
const au,LD,day = 149597870.7e3,384402e3,60^2*24
const pc,ly,GG = au*648000/π,365.25𝘤*day,𝘤*ħ/mP^2
const mₛ = GMsun/GG; const Jₛ = mₛ*au^2/day^2; export mₛ,Jₛ,au,day
const IAU = UnitSystem{Rᵤ*mₑ/μₑᵤ/0.001/Jₛ,ħ/day/Jₛ,day*𝘤/au,4π*1e-7*day^2/Jₛ,mₑ/mₛ}()

# aliased & humorous units

const mf = mass(90/lbm,Metric,English); const Jf = mf*(201.168/14day)^2
const FFF = UnitSystem{1000Rᵤ*mₑ/μₑᵤ*rankine/Jf,ħ/14day/Jf,14day*𝘤/201.168,0,mₑ/mf}()
const units, US, SI, MKS, temp = UnitSystem, UnitSystem, SI2019, Metric, temperature
const CGS, CGS2019, CGSm, CGSe, HLU = Gauss, EMU2019, EMU, ESU, LorentzHeaviside

# natural units

const Planck = UnitSystem{1,1,1,1,√(4π*αG)}()
const PlanckGauss = UnitSystem{1,1,1,4π,√αG}()
const Stoney = UnitSystem{1,αinv,1,4π,√(αG*αinv)}()
const Hartree = UnitSystem{1,1,αinv,4π/αinv^2,1}()
const Rydberg = UnitSystem{1,1,2αinv,π/αinv^2,1/2}()
const Schrodinger = UnitSystem{1,1,αinv,4π/αinv^2,√(αG*αinv)}()
const Electronic = UnitSystem{1,αinv,1,4π,1}()
const Natural = UnitSystem{1,1,1,1,1}()
const NaturalGauss = UnitSystem{1,1,1,4π,1}()
const QCD = UnitSystem{1,1,1,1,1/μₚₑ}()
const QCDGauss = UnitSystem{1,1,1,4π,1/μₚₑ}()
const QCDoriginal = UnitSystem{1,1,1,4π/αinv,1/μₚₑ}()

# physical constants

@pure electronmass(U::typeof(Planck),C::Coupling) = sqrt(4π*coupling(C))
@pure electronmass(U::typeof(PlanckGauss),C::Coupling) = sqrt(coupling(C))
@pure electronmass(U::UnitSystem{kB,ħ,𝘤,μ₀,√(αG*αinv)},C::Coupling) where {kB,ħ,𝘤,μ₀} = sqrt(coupling(C)/finestructure(C))
@pure electronmass(U::UnitSystem{kB,ħ,𝘤,μ₀,1/μₚₑ},C::Coupling) where {kB,ħ,𝘤,μ₀} = 1/protonelectron(C)
@pure permeability(U::UnitSystem{kB,ħ,𝘤,4π/αinv^2},C::Coupling) where {kB,ħ,𝘤} = 4π*finestructure(C)^2
@pure permeability(U::UnitSystem{kB,ħ,𝘤,π/αinv^2},C::Coupling) where {kB,ħ,𝘤} = π*finestructure(C)^2
@pure lightspeed(U::UnitSystem{kB,ħ,αinv},C::Coupling) where {kB,ħ} = 1/finestructure(C)
@pure lightspeed(U::UnitSystem{kB,ħ,2αinv},C::Coupling) where {kB,ħ} = 2/finestructure(C)
@pure planckreduced(U::UnitSystem{kB,αinv},C::Coupling) where kB = 1/finestructure(C)

@pure electronmass(U::UnitSystem{kB,ħ,𝘤,μ₀,mₑ},C::Coupling) where {kB,ħ,μ₀} = electronmass(planck(U),C)
@pure electronmass(U::UnitSystem{kB,ħ,100𝘤,μ₀,1000mₑ},C::Coupling) where {kB,ħ,μ₀} = 1000electronmass(SI,C)
@pure electronmass(U::UnitSystem{kB,ħ,𝘤,μ₀,mₑ/1000},C::Coupling) where {kB,ħ,μ₀} = electronmass(SI,C)/1000
@pure electronmass(U::UnitSystem{kB,ħ,𝘤,μ₀,mₑ2014},C::Coupling) where {kB,ħ,μ₀} = electronmass(planck(U),C)
@pure electronmass(U::UnitSystem{kB,ħ,𝘤,μ₀,mₑ1990},C::Coupling) where {kB,ħ,μ₀} = electronmass(planck(U),C)
@pure electronmass(U::UnitSystem{kB,ħ,𝘤/ftUS,μ₀,mₑ/slug},C::Coupling) where {kB,ħ,μ₀} = electronmass(SI,C)/slug
@pure permeability(U::UnitSystem{kB,ħ,𝘤,μ₀},C::Coupling) where {kB,ħ,𝘤} = finestructure(C)*2𝘩/𝘤/𝘦^2
@pure permeability(U::typeof(ESU2019),C::Coupling) = 1e3*permeability(SI,C)/𝘤^2
@pure permeability(U::typeof(EMU2019),C::Coupling) = 1e7*permeability(SI,C)
@pure permeability(U::typeof(CODATA),C::Coupling) = 2RK2014*finestructure(C)/𝘤
@pure permeability(U::typeof(Conventional),C::Coupling) = 2RK1990*finestructure(C)/𝘤

@pure molarmass(U::UnitSystem{1}) = 1
@pure molarmass(U::UnitSystem{kB},C::Coupling=universe(U)) = NA*electronmass(U,C)/electronunit(C)
@pure molarmass(U::UnitSystem{1e7*kB},C::Coupling=universe(U)) = 1000molarmass(SI2019,C)
@pure molarmass(U::UnitSystem{1e3*kB},C::Coupling=universe(U)) = molarmass(SI2019,C)/1000
@pure molarmass(U::UnitSystem{kB}) where kB = molarmass(CGS)/1000
@pure molarmass(U::UnitSystem{boltzmann(MTS)}) = molarmass(CGS)/1e6
@pure molarmass(U::UnitSystem{boltzmann(CGS)}) = molarmass(Natural)
@pure molarmass(U::UnitSystem{boltzmann(FFF)}) = molarmass(Natural)
@pure molarmass(U::UnitSystem{boltzmann(English)},C::Coupling=universe(U)) = 1000molarmass(SI2019,C)
@pure molarmass(U::UnitSystem{boltzmann(EnglishUS)}) = molarmass(Natural)
@pure molarmass(U::UnitSystem{boltzmann(IAU)}) = 1/1000mₛ

@pure luminousefficacy(U::UnitSystem{1}) = 1
@pure luminousefficacy(U::UnitSystem) = power(Kcd,SI2019,U)

include("kinematic.jl")
include("electromagnetic.jl")
include("thermodynamic.jl")
include("physics.jl")
include("systems.jl")

end # module

ext: jl


text:
using LightGraphs
using GraphIO
using ParserCombinator
using JuMP
using BioSequences
using CSV

export read_from_gfa, find_graph_component, get_summary, get_terminus, get_gfa, find_longest_path, get_fasta

include("struct.jl")
include("utils.jl")
include("read_gfa.jl")
include("graph_component.jl")
include("longest_path.jl")


text:
using ForwardAutodiff
using Test: @test, @testset

@testset "FowardAutodiff.jl tests" begin
    @testset "Transcendental Functions" begin
        x = randn()
        @test isapprox(D(exp)(x), exp(x))
        @test isapprox(D(sin)(x), cos(x))
        @test isapprox(D(cos)(x), -sin(x))
        @test isapprox(D(tan)(x), sec(x)*sec(x))
    end

    @testset "sqrt" begin
        x = randn()

        # We need to be in the complex domain explicitly, or sqrt(x < 0) will error out
        @test isapprox(D(sqrt)(x + 0im), 1/(2*sqrt(x + 0im)))
    end

    @testset "Type promotion" begin
        @test isapprox(D(sin)(2//3), cos(2//3)) # Rationals
        @test isapprox(D(sin)(1+2im), cos(1+2im)) # Complex integers!
        @test isapprox(D(cos)(1.0+3.5im), -sin(1.0+3.5im)) # Boring complex floating point numbers
        @test isapprox(D(sqrt)(2), 1/(2*sqrt(2.0))) # Automatically promotes to the correct type.
    end

    @testset "compound function" begin
        freq = sqrt(2)

        function f(x)
            return exp(-x/pi)*sin(2.0*pi*freq*x)
        end
        fprime = D(f)

        function laborious_fprime(x)
            return exp(-x/pi)*(2*pi*freq*cos(2*pi*freq*x) - sin(2*pi*freq*x)/pi)
        end

        @test isapprox(fprime(3), laborious_fprime(3))
    end

    @testset "partials of simple function" begin
        function linear(x, m, b)
            return m*x+b
        end
        dldx = D(1, linear)
        dldslope = D(2, linear)
        dldintercept = D(3, linear)

        x = randn()
        m = randn()
        b = randn()
        @test isapprox(dldx(x,m,b), m)
        @test isapprox(dldslope(x,m,b), x)
        @test isapprox(dldintercept(x,m,b), 1)
    end

    @testset "Cycloid" begin
        function cycloid(t)
            return [cos(2*pi*t), sin(2*pi*t), t]
        end

        @test isapprox(D(cycloid)(0), [0.0, 2*pi, 1])

        t = randn()
        @test isapprox(D(cycloid)(t), [-2*pi*sin(2*pi*t), 2*pi*cos(2*pi*t), 1])
    end

    @testset "rotations in 2D" begin
        function R2(theta)
            return [[cos(theta), sin(theta)], [-sin(theta), cos(theta)]]
        end

        @test isapprox(D(R2)(0), [[0.0, 1.0], [-1.0, 0.0]])
    end
end


text:
module NCCL

import CUDAdrv
import CuArrays
using Printf

import CUDAdrv: CuPtr, CuStream_t
import CuArrays: CuArray

const ext = joinpath(dirname(@__DIR__), "deps", "ext.jl")
isfile(ext) || error("NCCL.jl has not been built, please run Pkg.build(\"NCCL\").")
include(ext)
if !configured
    # default (non-functional) values for critical variables,
    # making it possible to _load_ the package at all times.
    const libnccl = nothing
end

include("base.jl")

function version()
    version = Ref{Cint}()
    @apicall(:ncclGetVersion, (Ref{Cint},), version)
    version[]
end

include("communicator.jl")
include("group.jl")
include("collectives.jl")

end


text:
using OrdinaryDiffEq, CuArrays, LinearAlgebra
u0 = cu(rand(1000))
A  = cu(randn(1000,1000))
f(du,u,p,t)  = mul!(du,A,u)
prob = ODEProblem(f,u0,(0.0f0,1.0f0)) # Float32 is better on GPUs!
sol = solve(prob,Tsit5())


text:
using Optim, LineSearches
using LinearAlgebra: I, norm
using TimerOutputs

"""
    energy(h, ipeps; χ, tol, maxiter)

return the energy of the `ipeps` 2-site hamiltonian `h` and calculated via a
ctmrg with parameters `χ`, `tol` and `maxiter`.
"""
function energy(h, model::HamiltonianModel, ipeps::IPEPS; χ::Int, tol::Real, maxiter::Int, verbose = false)
    ipeps = indexperm_symmetrize(ipeps)  # NOTE: this is not good
    D = getd(ipeps)^2
    s = gets(ipeps)
    ap = ein"abcdx,ijkly -> aibjckdlxy"(ipeps.bulk, conj(ipeps.bulk))
    ap = reshape(ap, D, D, D, D, s, s)
    a = ein"ijklaa -> ijkl"(ap)

    folder = "./data/$(model)/"
    mkpath(folder)
    chkp_file = folder*"vumps_env_D$(D)_chi$(χ).jld2"
    if isfile(chkp_file)
        rt = SquareVUMPSRuntime(a, chkp_file, χ; verbose = verbose)
    else
        rt = SquareVUMPSRuntime(a, Val(:random), χ; verbose = verbose)
    end

    env = vumps(rt; tol=tol, maxiter=maxiter, verbose = verbose)
    save(chkp_file, "env", env)
    e = expectationvalue(h, ap, env)
    return e
end

"""
    expectationvalue(h, ap, env)

return the expectationvalue of a two-site operator `h` with the sites
described by rank-6 tensor `ap` each and an environment described by
a `SquareCTMRGRuntime` `env`.
"""
function expectationvalue(h, ap, env::SquareVUMPSRuntime)
    AL,C,AR,FL,FR = env.AL,env.C,env.AR,env.FL,env.FR
    ap /= norm(ap)

    l = ein"abc,cde,anm,ef,ml,bnodpq -> lofpq"(FL,AL,conj(AL),C,conj(C), ap)
    r = ein"fgh,lkj,hij,okigrs -> folrs"(AR,conj(AR),FR,ap)
    e = ein"lofpq,folrs,pqrs -> "(l,r,h)[]
    n = ein"lofpp,folrr -> "(l,r)[]

    # AC = ein"asc,cb -> asb"(AL,C)
    # _, FL4 = bigleftenv(AL, M)
    # _, FR4 = bigrightenv(AR, M)
    # e2 = ein"dcba,def,aji,fghi,ckgepq,bjhkrs,pqrs -> "(FL4,AC,conj(AC),FR4,ap,ap,h)[]
    # n2 = ein"dcba,def,aji,fghi,ckgepq,bjhkrs -> pqrs"(FL4,AC,conj(AC),FR4,ap,ap)
    # n2 = ein"pprr -> "(n2)[]
    # @show e/n e2/n2 (e/n+e2/n2)/2

    return e/n
end

"""
    init_ipeps(model::HamiltonianModel; D::Int, χ::Int, tol::Real, maxiter::Int)

Initial `ipeps` and give `key` for use of later optimization. The key include `model`, `D`, `χ`, `tol` and `maxiter`.
The iPEPS is random initial if there isn't any calculation before, otherwise will be load from file `/data/model_D_chi_tol_maxiter.jld2`
"""
function init_ipeps(model::HamiltonianModel; D::Int, χ::Int, tol::Real, maxiter::Int, verbose = true)
    folder = "./data/$(model)/"
    mkpath(folder)
    key = (model, D, χ, tol, maxiter)
    chkp_file = folder*"$(model)_D$(D)_chi$(χ)_tol$(tol)_maxiter$(maxiter).jld2"
    if isfile(chkp_file)
        bulk = load(chkp_file)["ipeps"]
        verbose && println("load iPEPS from $chkp_file")
    else
        bulk = rand(D,D,D,D,2)
        verbose && println("random initial iPEPS")
    end
    ipeps = SquareIPEPS(bulk)
    ipeps = indexperm_symmetrize(ipeps)
    return ipeps, key
end

"""
    optimiseipeps(ipeps, h; χ, tol, maxiter, optimargs = (), optimmethod = LBFGS(m = 20))

return the tensor `bulk'` that describes an ipeps that minimises the energy of the
two-site hamiltonian `h`. The minimization is done using `Optim` with default-method
`LBFGS`. Alternative methods can be specified by loading `LineSearches` and
providing `optimmethod`. Other options to optim can be passed with `optimargs`.
The energy is calculated using vumps with key include parameters `χ`, `tol` and `maxiter`.
"""
function optimiseipeps(ipeps::IPEPS{LT}, key; f_tol = 1e-6, verbose= false, optimmethod = LBFGS(m = 20)) where LT
    model, D, χ, tol, maxiter = key
    h = hamiltonian(model)
    to = TimerOutput()
    f(x) = @timeit to "forward" real(energy(h, model, IPEPS{LT}(x); χ=χ, tol=tol, maxiter=maxiter, verbose=verbose))
    ff(x) = real(energy(h, model, IPEPS{LT}(x); χ=χ, tol=tol, maxiter=maxiter, verbose=verbose))
    g(x) = @timeit to "backward" Zygote.gradient(ff,x)[1]
    res = optimize(f, g,
        ipeps.bulk, optimmethod,inplace = false,
        Optim.Options(f_tol=f_tol,
        extended_trace=true,
        callback=os->writelog(os, key)),
        )
    println(to)
    return res
end

"""
    writelog(os::OptimizationState, key=nothing)

return the optimise infomation of each step, including `time` `iteration` `energy` and `g_norm`, saved in `/data/model_D_chi_tol_maxiter.log`. Save the final `ipeps` in file `/data/model_D_chi_tol_maxiter.jid2`
"""
function writelog(os::OptimizationState, key=nothing)
    message = "$(round(os.metadata["time"],digits=2))   $(os.iteration)   $(os.value)   $(os.g_norm)\n"

    printstyled(message; bold=true, color=:red)
    flush(stdout)

    model, D, χ, tol, maxiter = key
    if !(key === nothing)
        logfile = open("./data/$(model)/$(model)_D$(D)_chi$(χ)_tol$(tol)_maxiter$(maxiter).log", "a")
        write(logfile, message)
        close(logfile)
        save("./data/$(model)/$(model)_D$(D)_chi$(χ)_tol$(tol)_maxiter$(maxiter).jld2", "ipeps", os.metadata["x"])
    end
    return false
end


text:
__precompile__(true)

module Estuaries

using Reexport

using DataFrames
@reexport using DataStreams

import Base.convert
import Base: length, size
import Base: getindex, setindex!
import Base.copy
import Base.eltype
import Base.empty!, Base.delete!
import Base: insert!, merge!
import Base.isnull
import Base.show
import Base.haskey

import DataFrames.index
import DataFrames.eltypes
import DataFrames: ncol, nrow
import DataFrames: head, tail

include("abstracts.jl")
include("utils.jl")
include("estuary.jl")
include("datastreams.jl")

end # module


text:
"""
    SimpleRDArray{T,N} <: AbstractRDArray{T,N}

A simple implementation of resizable dense array.
"""
struct SimpleRDArray{T,N} <: AbstractRDArray{T,N}
    parent::Vector{T}
    sz::Size{N}
end
function SimpleRDArray(A::AbstractArray)
    vec = Vector{eltype(A)}(undef, length(A))
    copyto!(vec, A)
    sz = Size(A)
    return SimpleRDArray(vec, sz)
end
Base.IndexStyle(::Type{<:SimpleRDArray}) = IndexLinear()
Base.@propagate_inbounds Base.getindex(A::SimpleRDArray, i::Int) = parent(A)[i]
Base.@propagate_inbounds Base.getindex(A::SimpleRDArray, I::AbstractArray) = parent(A)[I]
Base.@propagate_inbounds Base.setindex!(A::SimpleRDArray, v, i::Int) = parent(A)[i] = v
Base.@propagate_inbounds Base.setindex!(A::SimpleRDArray, v, i::AbstractArray) =
    parent(A)[i] = v
Base.parent(A::SimpleRDArray) = A.parent
ArrayInterface.parent_type(::Type{<:SimpleRDArray{T,N}}) where {T,N} = Vector{T}
size_type(::Type{<:SimpleRDArray{T,N}}) where {T,N} = Size{N}
@inline getsize(A::SimpleRDArray) = A.sz


text:
@testset "UniformScaling rules" begin

    @testset "constructor" begin
        test_rrule(UniformScaling, rand())
    end

    @testset "+" begin
        # Forward
        test_frule(+, rand(3,3), I * rand(ComplexF64))
        test_frule(+, I, rand(3,3))

        # Reverse
        test_rrule(+, rand(3,3), I)
        test_rrule(+, rand(3,3), I * rand(ComplexF64))
        test_rrule(+, I, rand(3,3))
        test_rrule(+, I * rand(), rand(ComplexF64, 3,3))
    end

    @testset "-" begin
        # Forward
        test_frule(-, rand(3,3), I * rand(ComplexF64))
        test_frule(-, I, rand(3,3))

        # Reverse
        test_rrule(-, rand(3,3), I)
        test_rrule(-, rand(3,3), I * rand(ComplexF64))
        test_rrule(-, I, rand(3,3))
        test_rrule(-, I * rand(), rand(ComplexF64, 3,3))
    end

    @testset "Matrix" begin
        test_rrule(Matrix, I, (2, 2))
        test_rrule(Matrix{ComplexF64}, rand()*I, (3, 3))
        test_rrule(Matrix, rand(ComplexF64)*I, (2, 4))
    end

end


text:
module test_81_SaveFigure

import ModiaResult
ModiaResult.@usingModiaPlot

println("\n... test test_81_SaveFigure.jl:\n")

t = range(0.0, stop=10.0, length=100)
result = Dict{String,Any}("time" => t, "phi" => sin.(t))
info   = ModiaResult.resultInfo(result)
ModiaResult.printResultInfo(result)

plot(result, "phi", figure=2)
saveFigure(2, "test_saveFigure2.png")

end

max_stars_repo_name: ModiaSim/ModiaResult

max_issues_repo_name: ModiaSim/ModiaResult

max_forks_repo_name: ModiaSim/ModiaResult


text:
# Practical Machine learning
# Reinforcement learning - Q learning example
# Chapter 12

using DeepQLearning, NNGraph
# reload("DeepQLearning")
xs = linspace(0,360,100)
ys = round(sin(deg2rad(xs)),3)
deg2rad(xs)
# plot(x=xs,y=ys)

m = DQN(2,100,2)
alpha=0.0001; t_alpha =0.15
epsilon = 0.2; t_epsilon =0.45

init = [0. 0.]
s0 = NNMatrix(init'); a0 = 1; r0 = 0.
t = 0
for epoch = 1:1000 #0000
    t += 1
    avgReward = 0
    m.epsilon = epsilon * 1/t^t_epsilon
    m.alpha = alpha * 1/t^t_alpha
    for i = 2:length(xs)
        x, x2, y = xs[i],xs[i-1], ys[i]
        s = [x x2]
        s1 = NNMatrix(s')
        a1 = act(m,s1)
        r1 = (a1==1?-1:1) * sign(y)
        avgReward += r1
        if i > 2 learn(m,s0,a0,r1,s1) end
        s0 = s1; a0 = a1; r0 = r1
    end
    avgReward = avgReward / (length(xs)-2)
    if epoch % 100 == 0 println("$t $epoch avgReward = $(round(avgReward,3))   m.alpha=$(round(m.alpha,6))  m.epsilon=$(round(m.epsilon,6))") end
end


text: BaseTopicModel = Union{AbstractLDA, AbstractfLDA, AbstractCTM, AbstractfCTM}


text:
push!(LOAD_PATH, pwd())
using Irradiance, Sockets
const subnet = ip"192.168.0.255"
run_app(true, 8080, 37322, subnet)


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

""" Elasticity equations.

Field equation is:

    m∂²u/∂t² = ∇⋅σ - b

Weak form is: find u∈U such that ∀v in V

    δW := ∫ρ₀∂²u/∂t²⋅δu dV₀ + ∫S:δE dV₀ - ∫b₀⋅δu dV₀ - ∫t₀⋅δu dA₀ = 0

where

    ρ₀ = density
    b₀ = displacement load
    t₀ = displacement traction

Formulations
------------
plane stress, plane strain, 3D

References
----------

https://en.wikipedia.org/wiki/Linear_elasticity
https://en.wikipedia.org/wiki/Finite_strain_theory
https://en.wikipedia.org/wiki/Stress_measures
https://en.wikipedia.org/wiki/Mooney%E2%80%93Rivlin_solid
https://en.wikipedia.org/wiki/Strain_energy_density_function
https://en.wikipedia.org/wiki/Plane_stress
https://en.wikipedia.org/wiki/Hooke's_law

"""
mutable struct Elasticity <: FieldProblem
    # these are found from problem.properties for type Problem{Elasticity}
    formulation :: Symbol
    finite_strain :: Bool
    geometric_stiffness :: Bool
    store_fields :: Vector{Symbol}
end
function Elasticity()
    # formulations: plane_stress, plane_strain, continuum
    return Elasticity(:continuum, false, false, [])
end

function get_unknown_field_name(problem::Problem{Elasticity})
    return "displacement"
end

function get_formulation_type(problem::Problem{Elasticity})
    return :incremental
end

"""
    assemble!(assembly:Assembly, problem::Problem{Elasticity}, elements, time)

Start finite element assembly procedure for Elasticity problem.

Function groups elements to arrays by their type and assembles one element type
at time. This makes it possible to pre-allocate matrices common to same type
of elements.
"""
function assemble!(assembly::Assembly, problem::Problem{Elasticity},
                   elements::Vector{Element}, time)
    formulation = Val{problem.properties.formulation}
    for (element_type, elements_subset) in group_by_element_type(elements)
        assemble!(assembly, problem, elements_subset, time, formulation)
    end
end

include("problems_elasticity_2d.jl")

const Elasticity3DSurfaceElements = Union{Poi1,Tri3,Tri6,Quad4,Quad8,Quad9}
const Elasticity3DVolumeElements = Union{Tet4, Pyr5, Wedge6, Wedge15, Hex8, Tet10, Hex20, Hex27}

function initialize_internal_params!(params, ip, type_) #::Type{Val{:type_2d}})
    param_keys = keys(params)
    all_keys = ip.fields.keys
    ip_fields = filter(x->isassigned(all_keys, x), collect(1:length(all_keys)))

    if !("params_initialized" in ip_fields)
        for key in param_keys
            update!(ip, key, 0.0 => params[key])
        end
        if type_ == Val{:type_2d}
            update!(ip, "stress", 0.0 => [0.0,0.0,0.0])
            update!(ip, "strain", 0.0 => [0.0,0.0,0.0])
        elseif type_ == Val{:type_3d}
            update!(ip, "stress", 0.0 => [0.0,0.0,0.0,0.0,0.0,0.0])
            update!(ip, "strain", 0.0 => [0.0,0.0,0.0,0.0,0.0,0.0])
        else
            error("daa")
        end
        update!(ip, "prev_time", 0.0 => 0.0)
        update!(ip, "params_initialized", 0.0 => true)
    end
end

""" Assemble 3d continuum elements in general solid mechanics problem. """
function assemble!(assembly::Assembly,
                   problem::Problem{Elasticity},
                   elements::Vector{Element{El}},
                   time, ::Type{Val{:continuum}}) where El<:Elasticity3DVolumeElements
    props = problem.properties
    dim = get_unknown_field_dimension(problem)

    nnodes = length(El)
    ndofs = dim*nnodes
    BL = zeros(6, ndofs)
    BNL = zeros(9, ndofs)
    Km = zeros(ndofs, ndofs)
    Kg = zeros(ndofs, ndofs)
    f_int = zeros(ndofs)
    f_ext = zeros(ndofs)
    bi = BasisInfo(El)
    gradu = zeros(dim, dim)
    strain = zeros(dim, dim)
    strain_vec = zeros(6)
    stress_vec = zeros(6)
    F = zeros(dim, dim)
    D = zeros(6, 6)
    Dtan = zeros(6, 6)

    Bt_mul_D = zeros(ndofs, 6)
    Bt_mul_D_mul_B = zeros(ndofs, ndofs)
    Bt_mul_S = zeros(ndofs)

    for element in elements

        u = element("displacement", time)
        X = element("geometry", time)

        fill!(Km, 0.0)
        fill!(Kg, 0.0)
        fill!(f_int, 0.0)
        fill!(f_ext, 0.0)

        for ip in get_integration_points(element)
            eval_basis!(bi, X, ip)
            w = ip.weight*bi.detJ
            N = bi.N
            dN = bi.grad  # deriatives of basis functions w.r.t. X, i.e. ∂N/∂X
            grad!(bi, gradu, u)  # displacement gradient ∇u

            # calculate strain tensor and deformation gradient
            fill!(strain, 0.0)
            fill!(F, 0.0)
            F[:,:] += I
            if props.finite_strain
                strain[:,:] = 1/2 * (gradu + gradu' + gradu'*gradu)
                F[:,:] += gradu
            else
                strain[:,:] = 1/2 * (gradu + gradu')
            end

            strain_vec[1] = strain[1,1]
            strain_vec[2] = strain[2,2]
            strain_vec[3] = strain[3,3]
            strain_vec[4] = 2.0*strain[1,2]
            strain_vec[5] = 2.0*strain[2,3]
            strain_vec[6] = 2.0*strain[1,3]

            # material stiffness start

            fill!(BL, 0.0)
            if props.finite_strain
                for i=1:nnodes
                    BL[1, 3*(i-1)+1] = F[1,1]*dN[1,i]
                    BL[1, 3*(i-1)+2] = F[2,1]*dN[1,i]
                    BL[1, 3*(i-1)+3] = F[3,1]*dN[1,i]
                    BL[2, 3*(i-1)+1] = F[1,2]*dN[2,i]
                    BL[2, 3*(i-1)+2] = F[2,2]*dN[2,i]
                    BL[2, 3*(i-1)+3] = F[3,2]*dN[2,i]
                    BL[3, 3*(i-1)+1] = F[1,3]*dN[3,i]
                    BL[3, 3*(i-1)+2] = F[2,3]*dN[3,i]
                    BL[3, 3*(i-1)+3] = F[3,3]*dN[3,i]
                    BL[4, 3*(i-1)+1] = F[1,1]*dN[2,i] + F[1,2]*dN[1,i]
                    BL[4, 3*(i-1)+2] = F[2,1]*dN[2,i] + F[2,2]*dN[1,i]
                    BL[4, 3*(i-1)+3] = F[3,1]*dN[2,i] + F[3,2]*dN[1,i]
                    BL[5, 3*(i-1)+1] = F[1,2]*dN[3,i] + F[1,3]*dN[2,i]
                    BL[5, 3*(i-1)+2] = F[2,2]*dN[3,i] + F[2,3]*dN[2,i]
                    BL[5, 3*(i-1)+3] = F[3,2]*dN[3,i] + F[3,3]*dN[2,i]
                    BL[6, 3*(i-1)+1] = F[1,3]*dN[1,i] + F[1,1]*dN[3,i]
                    BL[6, 3*(i-1)+2] = F[2,3]*dN[1,i] + F[2,1]*dN[3,i]
                    BL[6, 3*(i-1)+3] = F[3,3]*dN[1,i] + F[3,1]*dN[3,i]
                end
            else
                for i=1:nnodes
                    BL[1, 3*(i-1)+1] = dN[1,i]
                    BL[2, 3*(i-1)+2] = dN[2,i]
                    BL[3, 3*(i-1)+3] = dN[3,i]
                    BL[4, 3*(i-1)+1] = dN[2,i]
                    BL[4, 3*(i-1)+2] = dN[1,i]
                    BL[5, 3*(i-1)+2] = dN[3,i]
                    BL[5, 3*(i-1)+3] = dN[2,i]
                    BL[6, 3*(i-1)+1] = dN[3,i]
                    BL[6, 3*(i-1)+3] = dN[1,i]
                end
            end

            # calculate stress

            fill!(D, 0.0)
            E = element("youngs modulus", ip, time)::Float64
            nu = element("poissons ratio", ip, time)::Float64
            la = E*nu/((1.0+nu)*(1.0-2.0*nu))
            mu = E/(2.0*(1.0+nu))
            D[1,1] = D[2,2] = D[3,3] = 2*mu + la
            D[4,4] = D[5,5] = D[6,6] = mu
            D[1,2] = D[2,1] = D[2,3] = D[3,2] = D[1,3] = D[3,1] = la

            # determine material model

            material_model = :linear_elasticity
            if haskey(element, "plasticity")
                material_model = :ideal_plasticity
            end

            # calculate stress vector based on material model

            if material_model == :linear_elasticity
                Dtan[:,:] = D[:,:]
                stress_vec[:] = Dtan * strain_vec
            end

            if material_model == :ideal_plasticity
                plastic_def = element("plasticity")[ip.id]

                calculate_stress! = plastic_def["type"]
                yield_surface_ = plastic_def["yield_surface"]
                params = plastic_def["params"]

                initialize_internal_params!(params, ip, Val{:type_3d})

                if time == 0.0
                    error("Given step time = $(time). Please select time > 0.0")
                end

                t_last = ip("prev_time", time)
                update!(ip, "prev_time", time => t_last)
                dt = time - t_last
                stress_last = ip("stress", t_last)
                strain_last = ip("strain", t_last)
                dstrain_vec = strain_vec - strain_last
                fill!(stress_vec, 0.0)
                fill!(Dtan, 0.0)
                plastic_strain = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
                calculate_stress!(stress_vec, stress_last, dstrain_vec, plastic_strain, D, params, Dtan, yield_surface_, time, dt, Val{:type_3d})

            end

            :strain in props.store_fields && update!(ip, "strain", time => strain_vec)
            :stress in props.store_fields && update!(ip, "stress", time => stress_vec)
            :stress11 in props.store_fields && update!(ip, "stress11", time => stress_vec[1])
            :stress22 in props.store_fields && update!(ip, "stress22", time => stress_vec[2])
            :stress33 in props.store_fields && update!(ip, "stress33", time => stress_vec[3])
            :stress12 in props.store_fields && update!(ip, "stress12", time => stress_vec[4])
            :stress23 in props.store_fields && update!(ip, "stress23", time => stress_vec[5])
            :stress13 in props.store_fields && update!(ip, "stress13", time => stress_vec[6])
            :plastic_strain in props.store_fields && update!(ip, "plastic_strain", time => plastic_strain)

            #Km += w*BL'*Dtan*BL
            mul!(Bt_mul_D, transpose(BL), Dtan)
            mul!(Bt_mul_D_mul_B, Bt_mul_D, BL)
            rmul!(Bt_mul_D_mul_B, w)
            for i=1:ndofs^2
                @inbounds Km[i] += Bt_mul_D_mul_B[i]
            end

            # material stiffness end

            if props.geometric_stiffness
                # take geometric stiffness into account

                fill!(BNL, 0.0)

                for i=1:size(dN, 2)
                    BNL[1, 3*(i-1)+1] = dN[1,i]
                    BNL[2, 3*(i-1)+1] = dN[2,i]
                    BNL[3, 3*(i-1)+1] = dN[3,i]
                    BNL[4, 3*(i-1)+2] = dN[1,i]
                    BNL[5, 3*(i-1)+2] = dN[2,i]
                    BNL[6, 3*(i-1)+2] = dN[3,i]
                    BNL[7, 3*(i-1)+3] = dN[1,i]
                    BNL[8, 3*(i-1)+3] = dN[2,i]
                    BNL[9, 3*(i-1)+3] = dN[3,i]
                end

                S3 = zeros(3*dim, 3*dim)
                S3[1,1] = stress_vec[1]
                S3[2,2] = stress_vec[2]
                S3[3,3] = stress_vec[3]
                S3[1,2] = S3[2,1] = stress_vec[4]
                S3[2,3] = S3[3,2] = stress_vec[5]
                S3[1,3] = S3[3,1] = stress_vec[6]
                S3[4:6,4:6] = S3[7:9,7:9] = S3[1:3,1:3]

                Kg += w*BNL'*S3*BNL

            end

            # internal load
            mul!(Bt_mul_S, transpose(BL), stress_vec)
            rmul!(Bt_mul_S, w)
            for i=1:ndofs
                @inbounds f_int[i] += Bt_mul_S[i]
            end

            # external load start

            if haskey(element, "displacement load")
                T = element("displacement load", ip, time)
                f_ext += w*vec(T*N)
            end

            for i=1:dim
                if haskey(element, "displacement load $i")
                    b = element("displacement load $i", ip, time)
                    f_ext[i:dim:end] += w*vec(b*N)
                end
            end

            # external load end

        end

        gdofs = get_gdofs(problem, element)

        # add contributions to K, Kg, f
        add!(assembly.K, gdofs, gdofs, Km)

        if props.geometric_stiffness
            add!(assembly.Kg, gdofs, gdofs, Kg)
        end

        add!(assembly.f, gdofs, f_ext - f_int)

    end

    return nothing
end

""" Elasticity equations, surface traction for continuum formulation. """
function assemble!(assembly::Assembly,
                   problem::Problem{Elasticity},
                   elements::Vector{Element{El}},
                   time, ::Type{Val{:continuum}}) where El<:Elasticity3DSurfaceElements

    props = problem.properties
    dim = get_unknown_field_dimension(problem)

    for element in elements
        nnodes = size(element, 2)
        f = zeros(dim*nnodes)

        has_concentrated_forces = false
        for ip in get_integration_points(element)
            detJ = element(ip, time, Val{:detJ})
            w = ip.weight*detJ
            N = element(ip, time)
            if haskey(element, "displacement traction force")
                T = element("displacement traction force", ip, time)
                f += w*vec(T*N)
            end
            for i in 1:dim
                if haskey(element, "displacement traction force $i")
                    T = element("displacement traction force $i", ip, time)
                    f[i:dim:end] += w*vec(T*N)
                end
                if haskey(element, "concentrated force $i")
                    has_concentrated_forces = true
                    T = element("concentrated force $i", ip, time)
                    f[i:dim:end] += w*vec(T*N)
                end
            end
            if haskey(element, "surface pressure")
                J = element(ip, time, Val{:Jacobian})'
                n = cross(J[:,1], J[:,2])
                n /= norm(n)
                # sign convention, positive pressure is towards surface
                p = -element("surface pressure", ip, time)
                f += w*p*vec(n*N)
            end
        end
        if has_concentrated_forces
            update!(element, "concentrated force", time => Any[f])
        end

        gdofs = get_gdofs(problem, element)
        add!(assembly.f, gdofs, f)

    end
end

""" Return strain tensor. """
function get_strain_tensor(problem, element, ip, time)
    gradu = element("displacement", ip, time, Val{:Grad})
    eps = 0.5*(gradu' + gradu)
    return eps
end

""" Return stress tensor. """
function get_stress_tensor(problem, element, ip, time)
    eps = get_strain_tensor(problem, element, ip, time)
    E = element("youngs modulus", ip, time)
    nu = element("poissons ratio", ip, time)
    mu = E/(2.0*(1.0+nu))
    la = E*nu/((1.0+nu)*(1.0-2.0*nu))
    S = la*tr(eps)*I + 2.0*mu*eps
    return S
end

""" Return stain vector in "ABAQUS" order 11, 22, 33, 12, 23, 13. """
function get_strain_vector(problem, element, ip, time)
    eps = get_strain_tensor(problem, element, ip, time)
    return [eps[1,1], eps[2,2], eps[3,3], eps[1,2], eps[2,3], eps[1,3]]
end

""" Return stress vector in "ABAQUS" order 11, 22, 33, 12, 23, 13. """
function get_stress_vector(problem, element, ip, time)
    S = get_stress_tensor(problem, element, ip, time)
    return [S[1,1], S[2,2], S[3,3], S[1,2], S[2,3], S[1,3]]
end

""" Make least squares fit for some field to nodes. """
function lsq_fit(problem, elements, field, time)
    A = SparseMatrixCOO()
    b = SparseMatrixCOO()
    volume = 0.0
    for element in elements
        gdofs = get_connectivity(element)
        for ip in get_integration_points(element)
            detJ = element(ip, time, Val{:detJ})
            w = ip.weight*detJ
            N = element(ip, time)
            f = field(problem, element, ip, time)
            add!(A, gdofs, gdofs, w*kron(N', N))
            for i=1:length(f)
                add!(b, gdofs, w*f[i]*N, i)
            end
            volume += w
        end
    end
    A = sparse(A)
    b = sparse(b)
    A = 1/2*(A + A')

    nz = get_nonzero_rows(A)
    F = ldlt(A[nz,nz])

    x = F \ b[nz, :]

    nodal_values = Dict(node_id => Vector(x[idx, :]) for (idx, node_id) in enumerate(nz))
    return nodal_values
end

""" Postprocessing, extrapolate strain to nodes using least-squares fit. """
function postprocess!(problem::Problem{Elasticity}, time::Float64, ::Type{Val{:strain}})
    elements = get_elements(problem)
    strain = lsq_fit(problem, elements, get_strain_vector, time)
    update!(elements, "strain", time => strain)
end

function postprocess!(problem::Problem{Elasticity}, time::Float64, ::Type{Val{:stress}})
    elements = get_elements(problem)
    stress = lsq_fit(problem, elements, get_stress_vector, time)
    update!(elements, "stress", time => stress)
end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*--------------------------------------------------------------------
#* Author: Gerson Balbueno Bicca
#* $Id: HairpinIncr.mso								$
#*------------------------------------------------------------------
type Properties_Average
	Properties_Average()=begin
		N=outers.N
		new(
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of zones",
				:Default=>2
			)),
			molweight(Dict{Symbol,Any}(
				:Brief=>"Average Mol Weight",
				:Default=>75,
				:Lower=>1,
				:Upper=>1e8
			)),
			fill(temperature(Dict{Symbol,Any}(
				:Brief=>"Average  Incremental Temperature",
				:Lower=>50
			)),(N)),
			fill(pressure(Dict{Symbol,Any}(
				:Brief=>"Average  Incremental Pressure",
				:Default=>1,
				:Lower=>1e-10,
				:Upper=>2e4,
				:DisplayUnit=>"kPa"
			)),(N)),
			fill(dens_mass(Dict{Symbol,Any}(
				:Brief=>"Stream Incremental Density" ,
				:Default=>1000,
				:Lower=>1e-3,
				:Upper=>5e5,
				:Symbol=>"\\rho"
			)),(N)),
			fill(viscosity(Dict{Symbol,Any}(
				:Brief=>"Stream Incremental Viscosity",
				:Lower=>0.0001,
				:Symbol=>"\\mu"
			)),(N)),
			fill(cp_mol(Dict{Symbol,Any}(
				:Brief=>"Stream Incremental Molar Heat Capacity",
				:Upper=>1e10
			)),(N)),
			fill(conductivity(Dict{Symbol,Any}(
				:Brief=>"Stream Incremental Thermal Conductivity",
				:Default=>1.0,
				:Lower=>1e-5,
				:Upper=>500
			)),(N)),
			[:N,],
			[:Mw,:T,:P,:rho,:Mu,:Cp,:K,]
		)
	end
	N::DanaInteger
	Mw::molweight
	T::Array{temperature}
	P::Array{pressure}
	rho::Array{dens_mass}
	Mu::Array{viscosity}
	Cp::Array{cp_mol}
	K::Array{conductivity}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Properties_Average
function atributes(in::Properties_Average,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Brief]="Average incremental physical properties of the streams."
	fields[:Info]="to be documented."
	drive!(fields,_)
	return fields
end
Properties_Average(_::Dict{Symbol,Any})=begin
	newModel=Properties_Average()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
const errors = CEnum.namemap(MDAL_Status)

function check_MDAL(p::Ptr{Nothing})
    if p != C_NULL
        return p
    else
        state = errors[Int(MDAL.MDAL_LastStatus())]
        if state == :None
            return p
        else
            error(state)
        end
    end
end

function check_MDAL(p::Any)
    state = errors[Int(MDAL.MDAL_LastStatus())]
    if state == :None
        p
    else
        error(state)
    end
end


text:
using Test, Flux, MLToolkit.Neural

@testset "Neural" begin
    @warn "`Neural.jl` is not tested."

    @testset "Architecture" begin
        σ, σlast = x -> sin.(x), x -> x
        Dout, B = 10, 5

        @testset "DenseNet" begin
            Din = 50
            mlp = DenseNet(Din, Dout, σ, σlast)
            Dhs = (40, 30, 20)
            mlp = DenseNet(Din, Dhs, Dout, σ, σlast)
            @test_throws DimensionMismatch mlp(randn(Din - 1, B))
            x = randn(Din, B)
            @test size(mlp(x)) == (Dout, B)
            x = randn(5, 5, 2, B)
            @test size(mlp(x)) == (Dout, B)
        end

        @testset "ConvNet" begin
            W, H, C = 28, 28, 1
            @test_throws ErrorException ConvNet((W - 1, H - 1, C), Dout, σ, σlast)
            convnet = ConvNet((W, H, C), Dout, σ, σlast)

            @test_throws DimensionMismatch convnet(randn(W - 1, H, C, B))
            @test_throws DimensionMismatch convnet(randn(W, H - 1, C, B))
            @test_throws DimensionMismatch convnet(randn(W, H, C + 1, B))
            for x in [
                randn(W, H, C, B),
                randn(W * H * C, B)
            ]
                @test size(convnet(x)) == (Dout, B)
            end
        end
    end
end


text:
"""
    struct ComplexSpace <: EuclideanSpace{ℂ}

A standard complex vector space ℂ^d with Euclidean inner product and no additional
structure. It is completely characterised by its dimension and whether its the normal space
or its dual (which is canonically isomorphic to the conjugate space).
"""
struct ComplexSpace <: EuclideanSpace{ℂ}
  d::Int
  dual::Bool
end
ComplexSpace(d::Integer; dual = false) = ComplexSpace(Int(d), dual)

# convenience constructor
Base.:^(::ComplexNumbers, d::Int) = ComplexSpace(d)
Base.getindex(::ComplexNumbers) = ComplexSpace
Base.getindex(::ComplexNumbers, d::Int) = ComplexSpace(d)

# Corresponding methods:
#------------------------
dim(V::ComplexSpace) = V.d
isdual(V::ComplexSpace) = V.dual
Base.axes(V::ComplexSpace) = Base.OneTo(dim(V))

Base.conj(V::ComplexSpace) = ComplexSpace(dim(V), !isdual(V))

Base.oneunit(::Type{ComplexSpace}) = ComplexSpace(1)
⊕(V1::ComplexSpace, V2::ComplexSpace) = isdual(V1) == isdual(V2) ?
    ComplexSpace(dim(V1)+dim(V2), isdual(V1)) :
    throw(SpaceMismatch("Direct sum of a vector space and its dual does not exist"))
fuse(V1::ComplexSpace, V2::ComplexSpace) = ComplexSpace(V1.d*V2.d)
flip(V::ComplexSpace) = dual(V)

Base.min(V1::ComplexSpace, V2::ComplexSpace) = isdual(V1) == isdual(V2) ?
    ComplexSpace(min(dim(V1), dim(V2)), isdual(V1)) :
    throw(SpaceMismatch("Mininum of space and dual space does not exist"))
Base.max(V1::ComplexSpace, V2::ComplexSpace) = isdual(V1) == isdual(V2) ?
    ComplexSpace(max(dim(V1), dim(V2)), isdual(V1)) :
    throw(SpaceMismatch("Maximum of space and dual space does not exist"))

Base.show(io::IO, V::ComplexSpace) = print(io, isdual(V) ? "(ℂ^$(V.d))'" : "ℂ^$(V.d)")
Base.show(io::IO, ::Type{ComplexSpace}) = print(io, "ComplexSpace")


text:
abstract type MHV2RuleModel <: MixingRule end

struct MHV2Rule{γ} <: MHV2RuleModel
    components::Array{String,1}
    activity::γ
end

@registermodel MHV2Rule
export MHV2Rule
function MHV2Rule(components::Vector{String}; activity = Wilson, userlocations::Vector{String}=String[],activity_userlocations::Vector{String}=String[], verbose::Bool=false)
    init_activity = activity(components;userlocations = activity_userlocations,verbose)

    model = MHV2Rule(components, init_activity)
    return model
end

function mixing_rule(model::RKModel,V,T,z,mixing_model::MHV2RuleModel,α,a,b,c)
    n = sum(z)
    x = z./n
    invn2 = (one(n)/n)^2
    g_E = excess_gibbs_free_energy(mixing_model.activity,1e5,T,z) / n
    b̄ = dot(z,Symmetric(b),z) * invn2
    c̄ = dot(z,c)/n

    ᾱ = a.*sqrt.(α.*α')./(b*R̄*T)

    q1 = -0.4783
    q2 = -0.0047
    c  = -q1*sum(x[i]*ᾱ[i,i] for i ∈ @comps)-q2*sum(x[i]*ᾱ[i,i]^2 for i ∈ @comps)-g_E/(R̄*T)-sum(x[i]*log(b̄/b[i,i]) for i ∈ @comps)

    ā = b̄*R̄*T*(-q1-sqrt(q1^2-4*q2*c))/(2*q2)
    return ā,b̄,c̄
end

function mixing_rule(model::PRModel,V,T,z,mixing_model::MHV2RuleModel,α,a,b,c)
    n = sum(z)
    x = z./n
    invn2 = (one(n)/n)^2
    g_E = excess_gibbs_free_energy(mixing_model.activity,1e5,T,z) / n
    b̄ = dot(z,Symmetric(b),z) * invn2
    c̄ = dot(z,c)/n

    ᾱ = a.*sqrt.(α.*α')./(b*R̄*T)

    q1 = -0.4347
    q2 = -0.003654
    c  = -q1*sum(x[i]*ᾱ[i,i] for i ∈ @comps)-q2*sum(x[i]*ᾱ[i,i]^2 for i ∈ @comps)-g_E/(R̄*T)-sum(x[i]*log(b̄/b[i,i]) for i ∈ @comps)

    ā = b̄*R̄*T*(-q1-sqrt(q1^2-4*q2*c))/(2*q2)
    return ā,b̄,c̄
end


text:
### A Pluto.jl notebook ###
# v0.12.4

using Markdown
using InteractiveUtils

# ╔═╡ bb6cae46-4108-11eb-2309-3d11f13a19fa
begin
	using Pkg
	Pkg.activate(mktempdir())
end

# ╔═╡ d8b334ae-4108-11eb-2546-99e75efc18d6
begin
	Pkg.add(["Plots", "PlutoUI","StatsBase", "DataStructures", "Unitful", "PhysicalConstants", "UnitfulEquivalences"])

	using Plots
	plotly()
	using PlutoUI
	using StatsBase
	using DataStructures
	using Unitful
	using PhysicalConstants.CODATA2018: c_0, h, ħ
	using UnitfulEquivalences
end

# ╔═╡ 56127994-4129-11eb-1b74-e92215af5b82
energy1 = 6.2u"keV"

# ╔═╡ 54d9f070-4129-11eb-13f6-2b138f584c42
wavelength1 = uconvert(u"pm", energy1, Spectral())

# ╔═╡ 83bc9488-4129-11eb-3c64-d71221b14aeb
freq1 = uconvert(u"Hz", energy1, Spectral())

# ╔═╡ 5990a14a-410b-11eb-086a-f129b7fdf82c
refracted_speed2 = uconvert(u"Mm/s", c_0 / 1.33)

# ╔═╡ 6209ea56-410c-11eb-0bb3-ff80947cdcb0
wavelength2 = 1.5u"μm"

# ╔═╡ a78b9c3a-4109-11eb-26ca-1556cfeb2b8c
quantum_energy2 = uconvert(u"eV", wavelength2, Spectral())

# ╔═╡ 84c21a6c-4109-11eb-166e-151c460d28c8
md"""
## E2

Light, which possesses in vacuum the frequency of 500 THz and quantum energy of $quantum_energy2, passes through
sheet of glass with the refractive index of 1.5. The speed of such light in glass is $refracted_speed2, the frequency is
500 THz, the wavelength is $wavelength2 and quantum energy is 2.07 eV.
"""

# ╔═╡ 1fd66092-412a-11eb-3468-bf3c5c496686
freq2 = uconvert(u"THz", wavelength2, Spectral())

# ╔═╡ fdae38fe-410c-11eb-18dd-9bae51fb4e00
freq3 = 350u"THz"

# ╔═╡ 220ec448-410d-11eb-17c3-d1682a1588c9
refractive3 = 1.5

# ╔═╡ fdbe0894-410c-11eb-0ace-3d269b2942aa
n = refractive3 ^2

# ╔═╡ 74fa7748-410c-11eb-0688-a37608618c96
md"""
## E3.1

The electric field strength that oscillates with the frequency of **$freq3** is reduced **$n** times in some substance
with respect to the vacuum. The refractive index of this substance is **$refractive3**
"""

# ╔═╡ 31bce1b8-410d-11eb-008b-bf04b66bf555
slit_length4 = 0.03u"mm"

# ╔═╡ 2e0551ea-410d-11eb-176c-650b9c498d25
distance_from_slit = 2u"m"

# ╔═╡ 0aaf6748-410e-11eb-035a-8be0950091d6
minima_distance = 8u"cm"

# ╔═╡ 0aa8bb82-410e-11eb-163c-6de6eda9e616
minima_angle = uconvert(u"rad", atan(minima_distance/2, distance_from_slit))

# ╔═╡ 70a82a76-410e-11eb-2f0e-295b06d5d8fb
wavelength4 = uconvert(u"nm", sin(minima_angle) * slit_length4)

# ╔═╡ 3bb15c6c-410d-11eb-1b1b-4b13708ceb70
md"""
## E4.1

The monochromatic light is penetrating a single **$slit_length4** wide vertical slit. On the screen positioned perpendicular to the light propagating direction at the distance of **$distance_from_slit** behind the slit, a diffraction pattern appears in which the two first minima are located at the distance of **$minima_distance** from each other. The first minimum is located at the angle of **$minima_angle** with respect to the diffraction pattern centre. The wavelength of the light used was **$wavelength4**
"""

# ╔═╡ 0d5b2ca4-412f-11eb-2dec-d1d78f201ff5
wavelength42 = 600u"nm"

# ╔═╡ 0d6c3466-412f-11eb-2ed7-6311848972ed
distance_from_grate = 4u"m"

# ╔═╡ 564061be-412f-11eb-13e1-894885ffb2be
maxima_distance = 24u"cm"

# ╔═╡ 0c352744-412f-11eb-37ae-ff4b3812f6c2
md"""
## E4.2
Red light, which possesses in vacuum the wavelength of **$wavelength42** is penetrating a diffraction grating. On the screen positioned perpendicular to the light propagating direction at the distance of **$distance_from_grate** from the grating, the two first order maxima in the diffraction pattern are located at the distance of **$maxima_distance** from each other.

"""

# ╔═╡ 97321f46-412f-11eb-2cda-b7e915c51a5b
maxima_angle = uconvert(u"rad", atan(maxima_distance/2, distance_from_grate))

# ╔═╡ cb110218-4130-11eb-1284-0fce3810c11f
grating_period = uconvert(u"mm", wavelength42/sin(maxima_angle))

# ╔═╡ c08c3f10-4117-11eb-059e-1fa7c75e4240
absorption_coef_blue = 0.02u"cm^-1"

# ╔═╡ ffd6d630-4117-11eb-2200-a147d74af92f
absorbtion_coef_green = 4u"cm^-1"

# ╔═╡ 38b97e80-4118-11eb-34c5-65de18666aed
absorption_coef_red = 7u"cm^-1"

# ╔═╡ 45bb6c38-4118-11eb-0ba5-f55cde355ccd
green_layer_thickness = uconvert(u"cm", 1/absorbtion_coef_green)

# ╔═╡ 40d8a62c-4118-11eb-15c3-e9e0c5d87609
n_intensity_reduced = 7

# ╔═╡ 793b338a-4119-11eb-1916-6513f6eed696
ℯ

# ╔═╡ afcc4710-4119-11eb-28db-fdf5aa042d8e
substance_color1 = "blue"

# ╔═╡ d0340d2e-4119-11eb-209d-7dcf7fa4c2e7
substance_color2 = "black"

# ╔═╡ c094cb3a-4117-11eb-168a-df762ad5739e
md"""
## E5
For some substance, the absorption coefficient for the blue light is **$absorption_coef_blue**. For the green light, it is **$absorbtion_coef_green** and for the red light **$absorption_coef_red**. This means that the intensity of the green light will be reduced Euler number (e = 2.7183..) times by passing through the **$green_layer_thickness** thick layer of the substance. The layer with the thickness of 1 cm reduces the intensity of green light e **$n_intensity_reduced** times. If we illuminate this substance with the white light, the colour of the substance will be **$substance_color1**. If we illuminate the substance with the blue light, the colour of the substance will be **$substance_color2**.

"""

# ╔═╡ 80572be0-411b-11eb-3f15-fbf4e6ed7cf3
electron_speed = 364.7u"km/s"

# ╔═╡ 5788010a-411e-11eb-2c69-bd73c5cd6d3f
ground_speed = 2188u"km/s"

# ╔═╡ 7fe77a98-411b-11eb-1835-75b0687dedfd
n6 = round(ground_speed / electron_speed)

# ╔═╡ 9f78c474-412b-11eb-1cf0-d7f1a152ffcf
reduced_planck = 1.05E-34u"kg*m^2*s^-1"

# ╔═╡ 7fc27112-411b-11eb-3b4b-97599735cd31
angular_momentum = uconvert(u"kg*m^2*s^-1", reduced_planck * n6)

# ╔═╡ 224bba9e-411f-11eb-2901-87f766e95a56
ground_frequency = 6580u"THz"

# ╔═╡ 3ed8da6c-411e-11eb-3a35-4164dfd2aa54
imagined_rotation_frequency = uconvert(u"THz", ground_frequency / n6 ^3)

# ╔═╡ 91b1a86c-411f-11eb-0700-a30281bd1f40
ground_energy = -13.60u"eV"

# ╔═╡ 8a2e2e76-411f-11eb-0da9-ed7c9382599f
net_energy = ground_energy / n6^2

# ╔═╡ 0396e5f2-4120-11eb-1502-f358ec34acd9
ground_radius = 0.529u"Å"

# ╔═╡ 03d6d638-4120-11eb-20d0-0fa5cffc6d3d
mean_distance = ground_radius * n6^2

# ╔═╡ 296b84ce-4121-11eb-2c07-1b1ff4cac46c
emitted_photon_wavelength = 410.6u"nm"

# ╔═╡ f74d6ff0-4120-11eb-2de1-9d7c40c9354c
quantum_energy6 = uconvert(u"eV", emitted_photon_wavelength, Spectral())

# ╔═╡ 66942310-4121-11eb-063c-37022902216f
energy_after_photon = net_energy - quantum_energy6

# ╔═╡ 81bc1422-4121-11eb-0cb5-215e6fb5c308
qn = round(sqrt(ground_energy/energy_after_photon))

# ╔═╡ 8238f47e-4121-11eb-06a7-a3f1bf8887cb
imag_rot_freq2 = uconvert(u"THz", ground_frequency / qn ^3)

# ╔═╡ 0e7c4e68-4122-11eb-3b1a-8f9808df6829
mean_distance2 = ground_radius * qn^2

# ╔═╡ 7fec90b4-411b-11eb-04db-a913940a9f22
md"""
## E6
In the atom of hydrogen, the electron moves according to the Bohr model. Its speed is **$electron_speed**. The angular momentum of this electron is **$angular_momentum**, the frequency of imagined rotation is **$imagined_rotation_frequency**, the net energy is **$net_energy** and the mean distance from the nucleus is **$mean_distance**. The electron comes closer to the nucleus and a photon with the wavelength of **$emitted_photon_wavelength** is emitted. The quantum energy of this radiation is **$quantum_energy6**. At the end of the transition, the frequency of the imagined rotation of the electron is **$imag_rot_freq2** and its mean distance from the nucleus is **$mean_distance2**.
"""

# ╔═╡ 46652b7e-4122-11eb-1576-29e34caf1b03
mass = 20u"mg"

# ╔═╡ 78e0e354-4122-11eb-1a11-0b1d36a38b70
diam = 5u"mm"

# ╔═╡ 46798970-4122-11eb-31a8-ab967162b1ae
h_bar = 1E-34u"J*s"

# ╔═╡ 4691ff34-4122-11eb-1614-b5df5b328525
speed_increase = uconvert(u"m/s", h_bar / (mass * diam))

# ╔═╡ c9767010-4125-11eb-3529-43e03fc0ed7b
time_for_diam = uconvert(u"s", diam / speed_increase)

# ╔═╡ 63ef8d68-412e-11eb-15fd-79cc263180f7
seconds_in_year = (365 * 24 * 60 * 60)u"s"

# ╔═╡ c9680bce-4125-11eb-1a30-93e346ab5ef8
time2 = uconvert(u"s", time_for_diam)

# ╔═╡ 451edba2-4122-11eb-0ea9-b3e5a7021bd1
md"""
# E7

The uranium atom (a microscopic object) with the mass of **$mass** and diameter of **$diam** moves in a certain direction. The location of the atom is fixed by its passing through the slit possessing the width almost equal to the atom diameter. The screen containing this slit is positioned perpendicular to the velocity of the atom. As a result of penetrating this slit, the speed of the atom receives some unpredictable increase perpendicular to the direction of previous motion. Using the reduced Planck constant 10–34 J*s as a typical value of impact uncertainty, we can conclude that the increase of the speed is **$speed_increase**. After **$time_for_diam** seconds, the atom deviates from its classical straight trajectory by its own diameter. This result means that we can predict the motion of uranium atom in a classical way for maximally **$time2** nanoseconds.

"""

# ╔═╡ 571de63e-412e-11eb-3b48-d9192483a1de
time2_years = time2 / seconds_in_year

# ╔═╡ aa066c66-4126-11eb-1909-713aeda74edc
mass8 = 9E-31u"kg"

# ╔═╡ c4375a50-4126-11eb-2667-7b693682e9a8
well_width = 2u"Å"

# ╔═╡ aa323b5c-4126-11eb-18fa-a5aa100d0678
qn8 = 5

# ╔═╡ aa487444-4126-11eb-306c-1d98f391cfc6
planck8 = 6.6E-34u"J*s"

# ╔═╡ aa708326-4126-11eb-2129-df2f5666aa93
wavelength = 2*well_width / qn8

# ╔═╡ aa9ab420-4126-11eb-206b-2bd429f38b13
lin_momentum = uconvert(u"kg*m/s", planck8 / wavelength)

# ╔═╡ aa805792-4126-11eb-20eb-b127d6ae2898
energy8 = uconvert(u"eV", lin_momentum^2/(2*mass8))

# ╔═╡ c98691d4-4125-11eb-3e71-d1ffa4fc0f3a
md"""
## E8

An electron with the mass of **$mass8** is located in a single-dimensional potential well with the width of **$well_width** (approximate diameter of the atom). Let us describe this electron according to the laws of quantum mechanics as a standing wave in the string, whereas the length of the string is equal to the width of the potential well. In the state described by the quantum number **$qn8**, the wavelength of the electron is **$wavelength**, its
linear momentum is **$lin_momentum** and its energy is **$energy8**.

"""

# ╔═╡ aa5f27d4-4126-11eb-19db-9b02f89cf0f0
md"""

## E9

An incandescent lamp has the nominal power of 60 W and the nominal voltage of 220 V. The cylindrical
filament of this lamp operating at the absolute temperature of 2650 K is made from tungsten (W). The
resistivity of tungsten at 2650 K is 0.64 μΩ*m. The electric field strength inside the filament is 500 V/m.
The resistance of the filament is 807 ohms, its length is 44 cm, the cross-sectional area of the filament is
3.49 *10–10 m2 and its radius is 10.5 micrometers. The cylinder surface area of the filament is 29.0 mm2.
The surface radiance of the filament is 2.07 MW/m2. If we describe the filament as a blackbody, its radiance
will be 2.80 MW/m2. The maximum of the spectral radiance of the lamp will be located in the spectrum at
the wavelength of 1094 nm.

"""

# ╔═╡ aa1d3ab8-4126-11eb-11d9-ab8a72c7e68f
power = 21u"W"

# ╔═╡ 244d7438-4131-11eb-1b84-c37c9c688b08
voltage = 12u"V"

# ╔═╡ 2d25a12a-4131-11eb-0cf6-7f77bfcf2071
temp = 2750u"K"

# ╔═╡ 3860b1f6-4131-11eb-198e-f3f101300a20
resistivity = 0.67u"μΩ*m"

# ╔═╡ 53aad108-4131-11eb-393a-3b7c4707017a
efs = 350u"V/m"

# ╔═╡ 624d2ab2-4131-11eb-0d4d-1f55c87e0abb
resistance = uconvert(u"Ω", voltage^2/power)

# ╔═╡ 7be2c68c-4131-11eb-3e4f-41b7e535bd69
length = uconvert(u"cm", voltage / efs)

# ╔═╡ b418c7fc-4131-11eb-36ad-0bcaf71c41b3
cross_section = uconvert(u"m^2", (resistivity * length) /  resistance)

# ╔═╡ 038ff9c2-4132-11eb-3ccb-eb745f9d9827
radius = uconvert(u"μm", sqrt(cross_section/pi))

# ╔═╡ 3c07c514-4132-11eb-0a3e-8b13abce99bb
silinder_surface = uconvert(u"mm^2", 2*pi*radius*length)

# ╔═╡ 9cbd7b60-4132-11eb-0692-67048ff2a347
wien_displacement = 2896u"μm*K"

# ╔═╡ 103478d8-4132-11eb-116e-6137e6dde875
radiant_exitance = uconvert(u"MW/m^2", power / silinder_surface)

# ╔═╡ Cell order:
# ╠═bb6cae46-4108-11eb-2309-3d11f13a19fa
# ╠═d8b334ae-4108-11eb-2546-99e75efc18d6
# ╠═56127994-4129-11eb-1b74-e92215af5b82
# ╠═54d9f070-4129-11eb-13f6-2b138f584c42
# ╠═83bc9488-4129-11eb-3c64-d71221b14aeb
# ╟─84c21a6c-4109-11eb-166e-151c460d28c8
# ╠═a78b9c3a-4109-11eb-26ca-1556cfeb2b8c
# ╠═5990a14a-410b-11eb-086a-f129b7fdf82c
# ╠═6209ea56-410c-11eb-0bb3-ff80947cdcb0
# ╠═1fd66092-412a-11eb-3468-bf3c5c496686
# ╟─74fa7748-410c-11eb-0688-a37608618c96
# ╠═fdae38fe-410c-11eb-18dd-9bae51fb4e00
# ╠═fdbe0894-410c-11eb-0ace-3d269b2942aa
# ╠═220ec448-410d-11eb-17c3-d1682a1588c9
# ╠═3bb15c6c-410d-11eb-1b1b-4b13708ceb70
# ╠═31bce1b8-410d-11eb-008b-bf04b66bf555
# ╠═2e0551ea-410d-11eb-176c-650b9c498d25
# ╠═0aaf6748-410e-11eb-035a-8be0950091d6
# ╠═0aa8bb82-410e-11eb-163c-6de6eda9e616
# ╠═70a82a76-410e-11eb-2f0e-295b06d5d8fb
# ╠═0c352744-412f-11eb-37ae-ff4b3812f6c2
# ╠═0d5b2ca4-412f-11eb-2dec-d1d78f201ff5
# ╠═0d6c3466-412f-11eb-2ed7-6311848972ed
# ╠═564061be-412f-11eb-13e1-894885ffb2be
# ╠═97321f46-412f-11eb-2cda-b7e915c51a5b
# ╠═cb110218-4130-11eb-1284-0fce3810c11f
# ╠═c094cb3a-4117-11eb-168a-df762ad5739e
# ╠═c08c3f10-4117-11eb-059e-1fa7c75e4240
# ╠═ffd6d630-4117-11eb-2200-a147d74af92f
# ╠═38b97e80-4118-11eb-34c5-65de18666aed
# ╠═45bb6c38-4118-11eb-0ba5-f55cde355ccd
# ╠═40d8a62c-4118-11eb-15c3-e9e0c5d87609
# ╠═793b338a-4119-11eb-1916-6513f6eed696
# ╠═afcc4710-4119-11eb-28db-fdf5aa042d8e
# ╠═d0340d2e-4119-11eb-209d-7dcf7fa4c2e7
# ╟─7fec90b4-411b-11eb-04db-a913940a9f22
# ╠═80572be0-411b-11eb-3f15-fbf4e6ed7cf3
# ╠═5788010a-411e-11eb-2c69-bd73c5cd6d3f
# ╠═7fe77a98-411b-11eb-1835-75b0687dedfd
# ╠═9f78c474-412b-11eb-1cf0-d7f1a152ffcf
# ╠═7fc27112-411b-11eb-3b4b-97599735cd31
# ╠═224bba9e-411f-11eb-2901-87f766e95a56
# ╠═3ed8da6c-411e-11eb-3a35-4164dfd2aa54
# ╠═91b1a86c-411f-11eb-0700-a30281bd1f40
# ╠═8a2e2e76-411f-11eb-0da9-ed7c9382599f
# ╠═0396e5f2-4120-11eb-1502-f358ec34acd9
# ╠═03d6d638-4120-11eb-20d0-0fa5cffc6d3d
# ╠═296b84ce-4121-11eb-2c07-1b1ff4cac46c
# ╠═f74d6ff0-4120-11eb-2de1-9d7c40c9354c
# ╠═66942310-4121-11eb-063c-37022902216f
# ╠═81bc1422-4121-11eb-0cb5-215e6fb5c308
# ╠═8238f47e-4121-11eb-06a7-a3f1bf8887cb
# ╠═0e7c4e68-4122-11eb-3b1a-8f9808df6829
# ╟─451edba2-4122-11eb-0ea9-b3e5a7021bd1
# ╠═46652b7e-4122-11eb-1576-29e34caf1b03
# ╠═78e0e354-4122-11eb-1a11-0b1d36a38b70
# ╠═46798970-4122-11eb-31a8-ab967162b1ae
# ╠═4691ff34-4122-11eb-1614-b5df5b328525
# ╠═c9767010-4125-11eb-3529-43e03fc0ed7b
# ╠═571de63e-412e-11eb-3b48-d9192483a1de
# ╠═63ef8d68-412e-11eb-15fd-79cc263180f7
# ╠═c9680bce-4125-11eb-1a30-93e346ab5ef8
# ╠═c98691d4-4125-11eb-3e71-d1ffa4fc0f3a
# ╠═aa066c66-4126-11eb-1909-713aeda74edc
# ╠═c4375a50-4126-11eb-2667-7b693682e9a8
# ╠═aa323b5c-4126-11eb-18fa-a5aa100d0678
# ╠═aa487444-4126-11eb-306c-1d98f391cfc6
# ╠═aa708326-4126-11eb-2129-df2f5666aa93
# ╠═aa9ab420-4126-11eb-206b-2bd429f38b13
# ╠═aa805792-4126-11eb-20eb-b127d6ae2898
# ╠═aa5f27d4-4126-11eb-19db-9b02f89cf0f0
# ╠═aa1d3ab8-4126-11eb-11d9-ab8a72c7e68f
# ╠═244d7438-4131-11eb-1b84-c37c9c688b08
# ╠═2d25a12a-4131-11eb-0cf6-7f77bfcf2071
# ╠═3860b1f6-4131-11eb-198e-f3f101300a20
# ╠═53aad108-4131-11eb-393a-3b7c4707017a
# ╠═624d2ab2-4131-11eb-0d4d-1f55c87e0abb
# ╠═7be2c68c-4131-11eb-3e4f-41b7e535bd69
# ╠═b418c7fc-4131-11eb-36ad-0bcaf71c41b3
# ╠═038ff9c2-4132-11eb-3ccb-eb745f9d9827
# ╠═3c07c514-4132-11eb-0a3e-8b13abce99bb
# ╠═9cbd7b60-4132-11eb-0692-67048ff2a347
# ╠═103478d8-4132-11eb-116e-6137e6dde875


text:
using BEAST
using CompScienceMeshes
using Test

dir = point(0,0,1)
width = 1.0
delay = 1.5
scaling = 1.0
sig = creategaussian(width, delay, scaling)

pw = BEAST.planewave(dir, 1.0, sig)

CompScienceMeshes.cartesian(p::typeof(dir)) = p
CompScienceMeshes.cartesian(p::Number) = p
x = point(1,0,0)
t = 1.0
@test pw(x,t) ≈ sig(t-dot(dir,x))

pw2 = BEAST.gradient(pw)
@test pw2.direction ≈ dir
@test pw2.polarisation ≈ -dir
@test pw2.speedoflight ≈ pw.speed_of_light

dsig = derive(sig)
@test pw2(x,t) ≈ -dir*dsig(t-dot(dir,x))

trc = dot(BEAST.n,pw2)
ch = simplex(
    point(0,0,0),
    point(1,0,0),
    point(0,1,0))
ctr = center(ch)
val = trc(ctr,t)
@test val ≈ -dsig(t-dot(dir,x))


text:
using BinaryBuilder, Pkg

include("../../../fancy_toys.jl")

name = "CUDA_loader"
version = v"0.1"

cuda_versions = [v"9.0", v"9.2", v"10.0", v"10.1", v"10.2", v"11.0", v"11.1", v"11.2", v"11.3"]
for cuda_version in cuda_versions
    cuda_tag = "$(cuda_version.major).$(cuda_version.minor)"
    include("build_$(cuda_tag).jl")

    any(should_build_platform.(triplet.(platforms))) || continue
    build_tarballs(ARGS, name, version, [], script, platforms, products, dependencies;
                   lazy_artifacts=true)
end

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil