text:
using GNSSTools
using ProgressMeter
using Distributions
using JLD
using PyPlot
pygui(true)
include("kalman_filter_gain_error_functions.jl")

tcxo_h_parms = [1.9e-6, 2e-6, 9.7e-8, 9.1e-11, 9.9e-15]./16.3676e6^2  # TCXO
ocxo_h_parms = [2.2e-9, 1.6e-10, 1.6e-10, 7.1e-13, 7.2e-16] ./ 10e6^2  # OCXO
h_parms = tcxo_h_parms
h₋₂ = h_parms[1]
h₋₁ = h_parms[2]
h₀ = h_parms[3]
sig_freq = L1_freq
t_length = 1  # s
phi_init = 0  # rad
f_d_init = 1000  # Hz
fd_rate_init = -0.7  # Hz/s

#----------------------------------SIMULATION----------------------------------#

σω = 1
q_a = 0.1
T = 1e-3
M = floor(Int, t_length/T)
t = calctvector(M, 1/T)

A = calcA(T, 3)
C = calcC(T, 3)
Q = calcQ(T, h₀, h₋₂, q_a, sig_freq, 3)
R = [phase_noise_variance(CN0, T)]

x⁺_init = [phi_init; 2π*f_d_init; 2π*fd_rate_init]
x⁺ = deepcopy(x⁺_init)
x_truth_noiseless = zeros(3, M)
ϕ_truth_noiseless = zeros(M)
x_truth = zeros(3, M)
ϕ_truth = zeros(M)
Δϕ_truth = zeros(M)

# Generate truth
x_truth_noiseless[:,1] = x⁺_init
x⁺ = deepcopy(x⁺_init)
for i in 2:M
    x⁺ = A*x⁺
    x_truth_noiseless[:,i] = x⁺
    ϕ_truth_noiseless[i] = (C*x⁺)[1]
end

# Simulate states
x_truth[:,1] = x⁺_init
x⁺ = deepcopy(x⁺_init)
for i in 2:M
    x⁺ = A*x⁺ + rand(MvNormal(Q))
    x_truth[:,i] = x⁺
    ϕ_truth[i] = (C*x⁺ + rand(MvNormal(R)))[1]
    if i > 1
        Δϕ_truth[i] = (C*(x_truth_noiseless[:,i] .- x_truth[:,i]))[1]
    end
end

fig = figure(figsize=(8, 5))

ax1 = fig.add_subplot(3,2,1)
ax1.plot(t, x_truth[1,:], "k-")
xlabel("Time (s)")
ylabel("ϕ State (rad)")

ax2 = fig.add_subplot(3,2,2)
ax2.plot(t, ϕ_truth, "k-")
xlabel("Time (s)")
ylabel("ϕ Meas. (rad)")

ax3 = fig.add_subplot(3,2,3)
ax3.plot(t, x_truth[2,:]./2π, "k-")
xlabel("Time (s)")
ylabel("f_d (Hz)")

ax4 = fig.add_subplot(3,2,4)
ax4.plot(t, Δϕ_truth.*(180/π), "k-")
xlabel("Time (s)")
ylabel("Δϕ Meas. (degrees)")

ax5 = fig.add_subplot(3,2,5)
ax5.plot(t, x_truth[3,:]./2π, "k-")
xlabel("Time (s)")
ylabel("f_d rate (Hz/s)")

ax6 = fig.add_subplot(3,2,6)
ax6.plot(t, (ϕ_truth_noiseless .- ϕ_truth).*(180/π), "k-")
# ax6.plot(t, ϕ_truth_noiseless.*(180/π), "k:")
xlabel("Time (s)")
ylabel("ϕ Truth (rad)")

subplots_adjust(wspace=0.4, hspace=0.6, bottom=0.15, left=0.12, right=0.88,
                top=0.93)

#------------------------------------FILTER------------------------------------#

CN0 = 45
σω = 1
q_a = 0.1
T = 1e-3
state_num = 3

A = calcA(T, state_num)
C = calcC(T, state_num)
Q = calcQ(T, h₀, h₋₂, q_a, sig_freq, state_num)
R = [phase_noise_variance(CN0, T)]

if state_num == 3
    P⁺ = [phase_noise_variance(B_init, CN0, T, h₀, h₋₁, h₋₂) 0 0;
          0 (2π/T)^2 0;
          0 0 σω^2]
    x_init = [phi_init; 2π*f_d_init; 2π*fd_rate_init]# + rand(MvNormal(P⁺))
else
    P⁺ = [phase_noise_variance(B_init, CN0, T, h₀, h₋₁, h₋₂) 0;
          0 (2π/T)^2]
    x_init = [phi_init; 2π*f_d_init]# + rand(MvNormal(P⁺))
end

K = zeros(state_num, M)
P = zeros(state_num, M)
x = zeros(state_num, M)
x[:,1] = x_init
P[:,1] = diag(P⁺)
pᵩ = zeros(M)
x⁺ = deepcopy(x_init)
for i in 2:M
    x⁻ = A*x⁺
    P⁻ = A*P⁺*A' + Q
    Δϕ = Δϕ_truth[i]
    Kᵢ = (P⁻*C')*pinv(C*P⁻*C' + R)
    P⁺ = (I - Kᵢ*C)*P⁻
    x⁺ = x⁻ + Kᵢ*[Δϕ]
    K[:,i] = Kᵢ
    P[:,i] = diag(P⁺)
    x[:,i] = x⁺
    pᵩ[i] = (C*P⁺*C')[1]
end

fig = figure(figsize=(8, 5))

ax1 = fig.add_subplot(3,2,1)
ax1.plot(t, x_truth_noiseless[1,:] .- x[1,:], "k-")
# ax1.plot(t, x_truth_noiseless[1,:], "k:")
xlabel("Time (s)")
ylabel("ϕ State (rad)")

ax2 = fig.add_subplot(3,2,2)
ax2.plot(t, sqrt.(pᵩ).*(180/π), "k-")
xlabel("Time (s)")
ylabel("pᵩ (degrees)")

ax3 = fig.add_subplot(3,2,3)
ax3.plot(t, x_truth_noiseless[2,:]./2π .- x[2,:]./2π, "k-")
# ax3.plot(t, x_truth_noiseless[2,:]./2π, "k:")
xlabel("Time (s)")
ylabel("f_d (Hz)")

ax4 = fig.add_subplot(3,2,4)
ax4.plot(t, Δϕ_truth.*(180/π), "k-")
xlabel("Time (s)")
ylabel("Δϕ Meas. (degrees)")

if state_num == 3
    ax5 = fig.add_subplot(3,2,5)
    ax5.plot(t, x_truth_noiseless[3,:]./2π .- x[3,:]./2π, "k-")
    # ax1.plot(t, x_truth_noiseless[3,:]./2π, "k:")
    xlabel("Time (s)")
    ylabel("f_d rate (Hz/s)")
end

ax6 = fig.add_subplot(3,2,6)
ax6.plot(t, K[1,:], label="k₁")
ax6.plot(t, K[2,:], label="k₂")
if state_num == 3
    ax6.plot(t, K[3,:], label="k₃")
end
xlabel("Time (s)")
ylabel("KF Gain")

subplots_adjust(wspace=0.4, hspace=0.6, bottom=0.15, left=0.12, right=0.88,
                top=0.93)

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

function test_KF(CN0, T)
    f_s = 2.5e6
    f_if = 0
    sig_freq = L1_freq
    B = 2*l1ca_chipping_rate
    Tsys = 535
    nADC = 8
    t_length = 1
    tcxo_h_parms = [1.9e-6, 2e-6, 9.7e-8, 9.1e-11, 9.9e-15]./16.3676e6^2  # TCXO
    h₋₂ = tcxo_h_parms[1]
    h₋₁ = tcxo_h_parms[2]
    h₀ = tcxo_h_parms[3]
    f_d0 = 1000
    f_dd0 = 0
    phi0 = 1.2
    state_num = 3
    qₐ = 1

    carrier_amp = sqrt(2*GNSSTools.k*Tsys)*10^(CN0/20)  # for sinusoid
    noise_amp = sqrt(GNSSTools.k*B*Tsys)
    adc_scale = 2^(nADC-1)
    N = floor(Int, t_length*f_s)
    t = calctvector(N, f_s)
    phase_noise = real.(generate_phase_noise(t_length, f_s, sig_freq, tcxo_h_parms))
    thermal_noise = noise_amp .* randn(ComplexF64, N)
    signal = carrier_amp.*cis.((2π).*(f_if .+ f_d0 .+ 0.5 .* f_dd0.*t).*t .+ phi0)
    phi_true = (2π).*(f_if .+ f_d0 .+ 0.5 .* f_dd0.*t.^2).*t .+ phi0
    data = signal.*cis.(phase_noise) .+ thermal_noise
    # data = signal
    sigmax = sqrt(maximum(abs2, data))
    # data = round.(data.*(adc_scale/sigmax))
    for i in 1:N
        data[i] = round(data[i]*adc_scale/sigmax)
        # Check that the max positive value for I/Q channels is
        # (2^(nADC - 1) - 1)
        if real(data[i]) == adc_scale
            data[i] = data[i] - 1
        end
        if imag(data[i]) == adc_scale
            data[i] = data[i] - 1im
        end
    end

    M = floor(Int, t_length/T)
    M_over_2 = floor(Int, M/2)
    Mₙ = floor(Int, T*f_s)
    σ²ᵩ_unfiltered = phase_noise_variance(CN0, T)
    A = calcA(T, state_num)
    C = calcC(T, state_num)
    R = [σ²ᵩ_unfiltered]
    Q = calcQ(T, h₀, h₋₂, qₐ, sig_freq, state_num)
    x = zeros(state_num, M)
    Δx = zeros(state_num, M)
    P = zeros(state_num, M)
    K = zeros(state_num, M)
    Δϕs = zeros(M)
    ZP = zeros(ComplexF64, M)
    integrated_phase_noise = zeros(M)
    if state_num == 3
        P⁻ = diagm([sqrt(σ²ᵩ_unfiltered), 2π/(2T), 2π*5]).^2
        x⁻ = [phi0, 2π*f_d0, 2π*f_dd0] .+ rand(MvNormal(P⁻))
    elseif state_num == 2
        P⁻ = diagm([sqrt(σ²ᵩ_unfiltered), 2π/(2T)]).^2
        x⁻ = [phi0, 2π*f_d0] .+ rand(MvNormal(P⁻))
    else
        @error("Invalid state_num value of $(state_num). Only 2 and 3 are valid.")
    end
    x⁻ = x⁻ .+ rand(MvNormal(P⁻))
    track_t = calctvector(M, 1/T)

    for i in 1:M
        if state_num == 3
            ϕ, ω, ω_dot = x⁻
        else
            ϕ, ω = x⁻
            ω_dot = 0
        end
        subdata = data[(i-1)*Mₙ+1:i*Mₙ]
        replica = cis.(-(2π.*f_if.*t[1:Mₙ] .+ ω.*t[1:Mₙ] .+ ω_dot.*t[1:Mₙ].^2 ./ 2 .+ ϕ))
        zp = sum(subdata.*replica)/Mₙ
        # Δϕ = atan(imag(zp)/real(zp))
        Δϕ = atan(imag(zp), real(zp))
        Kᵢ = (P⁻*C')*pinv(C*P⁻*C' + R)
        P⁺ = (I - Kᵢ*C)*P⁻*(I - Kᵢ*C)' + Kᵢ*R*Kᵢ'
        correction = Kᵢ.*Δϕ
        x⁺ = x⁻ + correction
        x[:,i] = x⁺
        Δx[:,i] = correction
        P[:,i] = diag(P⁺)
        K[:,i] = Kᵢ
        Δϕs[i] = Δϕ
        ZP[i] = zp
        sub_phase_noise = phase_noise[(i-1)*Mₙ+1:i*Mₙ]
        integrated_phase_noise[i] = mean(sub_phase_noise)
        x⁻ = A*x⁺
    end
    return (track_t, x, P, K, Δϕs, Δx, ZP, integrated_phase_noise,
            N, M, Mₙ, M_over_2, σ²ᵩ_unfiltered, phi_true)
end

CN0 = 45
T = 5e-3

track_t, x, P, K, Δϕs, Δx, ZP, integrated_phase_noise,
N, M, Mₙ, M_over_2, σ²ᵩ_unfiltered, phi_true = test_KF(CN0, T)

fig = figure()
ax1 = fig.add_subplot(3, 1, 1)
ax1.plot(track_t, Δϕs[:].*(180/π), "k.", label=L"\Delta\phi")
ax1.plot(track_t, Δx[1,:].*(180/π), label=L"\Delta\hat{\phi}")
ax1.plot(track_t, sqrt.(P[1,:]).*(180/π),
         linestyle=":", color="r")
ax1.plot(track_t, -sqrt.(P[1,:]).*(180/π), linestyle=":", color="r")
unfiltered_sigma_phi = std(Δϕs[M_over_2:end])
filtered_sigma_phi = std(Δx[1,M_over_2:end])
ax1.plot([track_t[1], track_t[end]], [unfiltered_sigma_phi, unfiltered_sigma_phi].*(180/π),
          linestyle=":", color="g")
ax1.plot([track_t[1], track_t[end]],
         [-unfiltered_sigma_phi, -unfiltered_sigma_phi].*(180/π),
         linestyle=":", color="g")
ax1.plot([track_t[1], track_t[end]],
         [filtered_sigma_phi, filtered_sigma_phi].*(180/π),
          linestyle=":", color="gray")
ax1.plot([track_t[1], track_t[end]],
         [-filtered_sigma_phi, -filtered_sigma_phi].*(180/π),
         linestyle=":", color="gray")
ax1.plot([track_t[1], track_t[end]],
         [sqrt(σ²ᵩ_unfiltered), sqrt(σ²ᵩ_unfiltered)].*(180/π),
          linestyle=":", color="orange")
ax1.plot([track_t[1], track_t[end]],
         [-sqrt(σ²ᵩ_unfiltered), -sqrt(σ²ᵩ_unfiltered)].*(180/π),
         linestyle=":", color="orange")
legend(ncol=1)
xlabel("Time (seconds)")
ylabel("Phase (deg)")

ax2 = fig.add_subplot(3, 1, 2)
ax2.plot(track_t, (phi_true[1:Mₙ:N] .+ integrated_phase_noise .- x[1,:]).*(180/π),
         label=L"\phi + \phi_{oscillator} - \hat{\phi}")
ax2.plot(track_t, sqrt.(P[1,:]).*(180/π), label=L"\sqrt{P_{1,1}}",
         linestyle=":", color="r")
ax2.plot(track_t, -sqrt.(P[1,:]).*(180/π), linestyle=":", color="r")
unfiltered_sigma_phi = std(Δϕs[M_over_2:end])
filtered_sigma_phi = std(Δx[1,M_over_2:end])
ax2.plot([track_t[1], track_t[end]], [unfiltered_sigma_phi, unfiltered_sigma_phi].*(180/π),
          label=L"\sigma_{\Delta\phi}", linestyle=":", color="g")
ax2.plot([track_t[1], track_t[end]],
         [-unfiltered_sigma_phi, -unfiltered_sigma_phi].*(180/π),
         linestyle=":", color="g")
ax2.plot([track_t[1], track_t[end]],
         [filtered_sigma_phi, filtered_sigma_phi].*(180/π),
          label=L"\sigma_{\Delta\hat{\phi}}", linestyle=":", color="gray")
ax2.plot([track_t[1], track_t[end]],
         [-filtered_sigma_phi, -filtered_sigma_phi].*(180/π),
         linestyle=":", color="gray")
ax2.plot([track_t[1], track_t[end]],
         [sqrt(σ²ᵩ_unfiltered), sqrt(σ²ᵩ_unfiltered)].*(180/π),
          label=L"\sqrt{\sigma^2_\phi}}", linestyle=":", color="orange")
ax2.plot([track_t[1], track_t[end]],
         [-sqrt(σ²ᵩ_unfiltered), -sqrt(σ²ᵩ_unfiltered)].*(180/π),
         linestyle=":", color="orange")
legend()
xlabel("Time (seconds)")
ylabel("Phase (deg)")

ax3 = fig.add_subplot(3, 1, 3)
ax3.plot(track_t, real.(ZP), label="real(ZP)")
ax3.plot(track_t, imag.(ZP), label="imag(ZP)")
legend()
xlabel("Time (seconds)")
ylabel("ZP")

function test_KF_errs(CN0, T, iterations)
    P₁₁ = 0.
    σ²_Δϕ_unfilt = 0.
    σ²_Δϕ_filt = 0.
    σ²ᵩ_unfiltered = phase_noise_variance(CN0, T)
    for iteration in 1:iterations
        track_t, x, P, K, Δϕs, Δx, ZP, integrated_phase_noise,
        N, M, Mₙ, M_over_2, σ²ᵩ_unfiltered, phi_true = test_KF(CN0, T)
        σ²_Δϕ_unfilt += var(Δϕs[M_over_2:end])
        σ²_Δϕ_filt += var(Δx[1,M_over_2:end])
        P₁₁ += mean(P[1,M_over_2:end])
    end
    P₁₁_sqrt = sqrt(P₁₁/iterations)
    σ_Δϕ_unfilt = sqrt(σ²_Δϕ_unfilt/iterations)
    σ_Δϕ_filt = sqrt(σ²_Δϕ_filt/iterations)
    σᵩ_unfiltered = sqrt(σ²ᵩ_unfiltered)
    return (P₁₁_sqrt, σ_Δϕ_unfilt, σ_Δϕ_filt, σᵩ_unfiltered)
end

CN0s = Array(range(20, 50, step=1))
Ts = [1e-3, 2e-3]
# Ts = [1e-3]
iterations = 5

P₁₁_sqrts = []
σ_Δϕ_unfilts = []
σ_Δϕ_filts = []
σᵩ_unfiltereds = []

for j in 1:length(Ts)
    T = Ts[j]
    N = length(CN0s)
    P₁₁_sqrts_ = zeros(N)
    σ_Δϕ_unfilts_ = zeros(N)
    σ_Δϕ_filts_ = zeros(N)
    σᵩ_unfiltereds_ = zeros(N)
    p = Progress(N, 1, "Processing...")
    for i in 1:N
        CN0 = CN0s[i]
        P₁₁_sqrt, σ_Δϕ_unfilt, σ_Δϕ_filt, σᵩ_unfiltered = test_KF_errs(CN0, T, iterations)
        # println([P₁₁_sqrt, σ_Δϕ_unfilt, σ_Δϕ_filt, σᵩ_unfiltered])
        P₁₁_sqrts_[i] = P₁₁_sqrt
        σ_Δϕ_unfilts_[i] = σ_Δϕ_unfilt
        σ_Δϕ_filts_[i] = σ_Δϕ_filt
        σᵩ_unfiltereds_[i] = σᵩ_unfiltered
        next!(p)
    end
    push!(P₁₁_sqrts, P₁₁_sqrts_)
    push!(σ_Δϕ_unfilts, σ_Δϕ_unfilts_)
    push!(σ_Δϕ_filts, σ_Δϕ_filts_)
    push!(σᵩ_unfiltereds, σᵩ_unfiltereds_)
end

fig = figure(figsize=(7.5, 5))
ax1 = fig.add_subplot(2, 2, 1)
ax1.plot(CN0s, P₁₁_sqrts[1].*(180/π), label=L"\sqrt{P_{11}}")
ax1.plot(CN0s, σ_Δϕ_filts[1].*(180/π), label=L"\sigma_{\Delta\hat{\phi}}")
hlines(y=10log10(15), xmin=CN0s[1], xmax=CN0s[end], color="grey",
       linestyle=":", label=string(L"15^\circ", " Tracking Limit"))
ax1.set_yscale("log")
legend()
xlabel("C/N₀ (dB⋅Hz)")
ylabel("Carrier Phase Error (deg)")

ax2 = fig.add_subplot(2, 2, 2)
ax2.plot(CN0s, σᵩ_unfiltereds[1].*(180/π), label=L"\sigma_{\phi}")
ax2.plot(CN0s, σ_Δϕ_unfilts[1].*(180/π), label=L"\sigma_{\Delta\phi}")
ax2.set_yscale("log")
legend()
xlabel("C/N₀ (dB⋅Hz)")
ylabel("Carrier Phase Error (deg)")

ax3 = fig.add_subplot(2, 2, 3)
ax3.plot(CN0s, P₁₁_sqrts[2].*(180/π), label=L"\sqrt{P_{11}}")
ax3.plot(CN0s, σ_Δϕ_filts[2].*(180/π), label=L"\sigma_{\Delta\hat{\phi}}")
hlines(y=10log10(15), xmin=CN0s[1], xmax=CN0s[end], color="grey",
       linestyle=":", label=string(L"15^\circ", " Tracking Limit"))
ax3.set_yscale("log")
legend()
xlabel("C/N₀ (dB⋅Hz)")
ylabel("Carrier Phase Error (deg)")

ax4 = fig.add_subplot(2, 2, 4)
ax4.plot(CN0s, σᵩ_unfiltereds[2].*(180/π), label=L"\sigma_{\phi}")
ax4.plot(CN0s, σ_Δϕ_unfilts[2].*(180/π), label=L"\sigma_{\Delta\phi}")
ax4.set_yscale("log")
legend()
xlabel("C/N₀ (dB⋅Hz)")
ylabel("Carrier Phase Error (deg)")

subplots_adjust(wspace=0.3, bottom=0.15, left=0.08, right=0.92, top=0.92)


text:
"Simulator that records the state trajectory, including the start state."
@kwdef struct StateRecorder <: Simulator
    max_steps::Union{Int,Nothing} = nothing
end

function simulate(sim::StateRecorder, domain::Domain, state::State, actions)
    trajectory = [state]
	for (t, act) in enumerate(actions)
        state = transition(domain, state, act)
		push!(trajectory, state)
		sim.max_steps !== nothing && t >= sim.max_steps && break
    end
    return trajectory
end

ext: jl


text:
struct Chain<:AbstractLattice{1, 2}
    length::Int
    Chain(length::Int) = new(length)
end

nsite(c::Chain) = c.length
neighbors(c::Chain, i::Int) = [mod(i-2, c.length)+1, mod(i, c.length)+1]
neighbors(ci::Pair{Chain, Int}) = [mod(ci.second-2, nsite(ci.first))+1, mod(ci.second, nsite(ci.first))+1]
Base.getindex(c::Chain, i::Int) = c=>i
nunitcells(c::Chain) = 1


text:
export translate, transpose, louden, randomnotes, subdivision
export velocities, positions, pitches, durations, relpos
export timesort, timesort!, combine

velocities(notes::Notes) = [Int(x.velocity) for x in notes]
positions(notes::Notes) = [Int(x.position) for x in notes]
pitches(notes::Notes) = [Int(x.pitch) for x in notes]
durations(notes::Notes) = [Int(x.duration) for x in notes]

"""
    randomnotes(n::Int, tpq = 960)
Generate some random notes that start sequentially.
"""
function randomnotes(n::Int, tpq = 960, randchannel = false)
    notes = Note[]
    prevpos = 0
    durran = 1:tpq
    posran = 0:4*tpq
    for i in 1:n
        if randchannel
            note = Note(rand(UInt8), rand(0:0x7F), prevpos + rand(posran), rand(durran), rand(0:127))
        else
            note = Note(rand(UInt8), rand(0:127), prevpos + rand(posran), rand(durran))
        end
        push!(notes, note)
        prevpos = note.position
    end
    return Notes(notes, tpq)
end

"""
    translate(notes, ticks)
Translate the `notes` for the given amount of `ticks`.
Also works for a single note.
"""
translate(notes::Notes, ticks) = Notes(translate(notes.notes, ticks), notes.tpq)
translate(notes::Vector{N}, ticks) where {N<:AbstractNote} =
[N(n.pitch, n.velocity, n.position + ticks, n.duration, n.channel) for n in notes]
translate(n::N, ticks) where {N<:AbstractNote} =
N(n.pitch, n.velocity, n.position+ticks, n.duration, n.channel)

"""
    translate!(notes, ticks)
In-place version of [`translate`](@ref).
"""
function translate!(notes::Notes, ticks)
    for note in notes
        note.position += ticks
    end
end

"""
    transpose(notes, semitones)
Transpose the `notes` by the given amount of `semitones`.
Also works for a single note.
"""
Base.transpose(notes::Notes, semitones) =
Notes(transpose(notes.notes, semitones), notes.tpq)
Base.transpose(notes::Vector{N}, semitones) where {N<:AbstractNote} =
[Note(n.pitch + semitones, n.velocity, n.position, n.duration, n.channel) for n in notes]
Base.transpose(n::N, ticks) where {N<:AbstractNote} =
N(n.pitch+ticks, n.velocity, n.position, n.duration, n.channel)

"""
    louden(notes, v::Int)
Change the velocity of the notes by `v` (which could also be negative).
Also works for a single note.
"""
louden(notes::Notes, v) =
Notes(louden(notes.notes, v), notes.tpq)
louden(notes::Vector{N}, v) where {N<:AbstractNote} =
[Note(n.pitch, n.velocity + v, n.position, n.duration, n.channel) for n in notes]
louden(n::N, ticks) where {N<:AbstractNote} =
N(n.pitch, n.velocity+ticks, n.position, n.duration, n.channel)

"""
    timesort!(notes::Notes)
In-place sort the `notes` by their temporal position.
Use `timesort` for a non-mutating version.
"""
function timesort!(notes::Notes)
    sort!(notes.notes, by = x -> x.position)
    return notes
end
timesort(notes) = timesort!(copy(notes))

"""
    combine(note_container) -> notes
Combine the given container (either `Array{Notes}` or `Dict{Any, Notes}`) into
a single `Notes` instance. In the process, sort the notes by position in the
final container.
"""
function combine(notearray::AbstractArray{<:Notes}, tsort = true)
    notes = copy(notearray[1])
    for i in 2:length(notearray)
        append!(notes, notearray[i])
    end
    tsort && timesort!(notes)
    return notes
end

function combine(notearray::AbstractArray{<:AbstractArray{<:AbstractNote}}, tsort=true)
    notes = copy(notearray[1])
    for i in 2:length(notearray)
        append!(notes, notearray[i])
    end
    tsort && timesort!(notes)
    return notes
end

function combine(notedict::Dict{<:Any, Notes{N}}, tsort = true) where {N}
    n = Notes(N[], first(values(notedict)).tpq)
    for (k, v) in notedict
        append!(n, v)
    end
    tsort && timesort!(n)
    return n
end

"""
    relpos(notes::Notes, grid)
Return the *relative* positions of the notes with respect to the current
`grid`, i.e. all notes are brought within one quarter note.
"""
function relpos(notes::Notes, grid)
    tpq = notes.tpq
    c = (grid[end-1] + (1 - grid[end-1])/2) * tpq
    rpos = zeros(Int, length(notes))
    for (i, n) in enumerate(notes)
        m = mod1(Int(n.position), tpq)
        if m ≥ c
            rpos[i] = m - tpq
        else
            rpos[i] = m
        end
    end
    return rpos
end

"""
    repeat(notes, i = 1)
Repeat the `notes` `i` times, by successively adding duplicates of `notes`
shifted by the total duration of `notes`. Return a single `Notes` container
for convenience.

The function assumes that notes are `timesort`ed.
"""
Base.repeat(n::Notes, i::Int = 1) = Notes(repeat(n.notes, i), n.tpq)
function Base.repeat(n::Vector{<:AbstractNote}, i::Int = 1)
    maxdur = maximum(a.position + a.duration for a in n)
    r = [copy(n)]
    for j in 1:i
        push!(r, translate(r[end], maxdur))
    end
    return combine(r, false)
end


text:
using GaussQuadrature, LinearAlgebra

function cuby_sphere(ξ¹, ξ², ξ³; radii = nothing, rectangle = nothing, face = 1)
    if isnothing(radii)
        r², r¹ = 1.0, 0.8
    else
        r², r¹ = radii.r², radii.r¹
    end
    if isnothing(rectangle)
        bˣ, aˣ = 1.0, -1.0
        bʸ, aʸ = 1.0, -1.0
    else
        aˣ, bˣ, aʸ, bʸ = rectangle.aˣ, rectangle.bˣ, rectangle.aʸ, rectangle.bʸ
    end

    R = (ξ³ + 1) / 2 * (r² - r¹) + r¹
    η² = (ξ¹ + 1) / 2 * (bˣ - aˣ) + aˣ
    η³ = (ξ² + 1) / 2 * (bʸ - aʸ) + aʸ
    facewarp(x) = tan(π / 4 * x) # x # asin(x) * 2/ π # x + 1.0*sin(π *x)
    ρ = sqrt(1 + facewarp(η²)^2 + facewarp(η³)^2)
    x¹ = R / ρ
    x² = R * facewarp(η²) / ρ
    x³ = R * facewarp(η³) / ρ
    if face == 1
        return x¹, x², x³
    elseif face == 2
        return -x¹, -x², -x³
    elseif face == 3
        return x², x¹, x³
    elseif face == 4
        return -x², -x¹, -x³
    elseif face == 5
        return x³, x², x¹
    elseif face == 6
        return -x³, -x², -x¹
    else
        return nothing
    end
    return nothing
end

N = 6
ξ¹, ω¹ = GaussQuadrature.legendre(N, GaussQuadrature.both)
ξ² = copy(ξ¹)
ξ³ = copy(ξ¹)

ξ¹ = reshape(ξ¹, (N, 1, 1))
ξ² = reshape(ξ¹, (1, N, 1))
ξ³ = reshape(ξ¹, (1, 1, N))

cuby = cuby_sphere.(ξ¹, ξ², ξ³)[:]

cuby1 = cuby_sphere.(ξ¹, ξ², ξ³, face = 1)[:]
cuby2 = cuby_sphere.(ξ¹, ξ², ξ³, face = 2)[:]

cuby3 = cuby_sphere.(ξ¹, ξ², ξ³, face = 3)[:]
cuby4 = cuby_sphere.(ξ¹, ξ², ξ³, face = 4)[:]

cuby5 = cuby_sphere.(ξ¹, ξ², ξ³, face = 5)[:]
cuby6 = cuby_sphere.(ξ¹, ξ², ξ³, face = 6)[:]

cuby = [cuby1..., cuby2...]

cuby = [cuby..., cuby3..., cuby4...]
cuby = [cuby..., cuby5..., cuby6...]

vec = ξ¹
ξp = [vec[i]^(j - 1) for i = 1:length(vec), j = 1:length(vec)]
dξp = [(j - 1) * vec[i]^(j - 2) for i = 1:length(vec), j = 1:length(vec)]
∂ξ¹ = dξp / ξp
∂ξ² = copy(∂ξ¹)
∂ξ³ = copy(∂ξ¹)

ω¹ = reshape(ω¹, (N, 1, 1))
ω² = reshape(ω¹, (1, N, 1))
ω³ = reshape(ω¹, (1, 1, N))

##
#using GLMakie
#fig, ax, sc = scatter(cuby)

##
# Check Metric terms
scale = 1
rectangle = (; aˣ = -1.0 / scale, bˣ = 1.0 / scale, aʸ = -1.0 / scale, bʸ = 1.0 / scale)
face1 = cuby_sphere.(ξ¹, ξ², ξ³, face = 1, rectangle = rectangle)
L, M, N = size(face1)
x¹ = [face1[i, j, k][1] for i = 1:L, j = 1:M, k = 1:N]
x² = [face1[i, j, k][2] for i = 1:L, j = 1:M, k = 1:N]
x³ = [face1[i, j, k][3] for i = 1:L, j = 1:M, k = 1:N]

∂x¹∂ξ¹ = copy(x¹) .* 0.0
∂x¹∂ξ² = copy(x¹) .* 0.0
∂x¹∂ξ³ = copy(x¹) .* 0.0

∂x²∂ξ¹ = copy(x¹) .* 0.0
∂x²∂ξ² = copy(x¹) .* 0.0
∂x²∂ξ³ = copy(x¹) .* 0.0

∂x³∂ξ¹ = copy(x¹) .* 0.0
∂x³∂ξ² = copy(x¹) .* 0.0
∂x³∂ξ³ = copy(x¹) .* 0.0

jacobian = zeros((3, 3, size(x¹)...)) # aᵢ
ijacobian = copy(jacobian) # aⁱ
detjacobian = copy(x¹) .* 0.0

kopriva = zeros(3, 3, size(x¹)...) # Jaⁱ

yyz = zeros(size(x¹))

for i = 1:L, j = 1:M, k = 1:N
    for ii = 1:L
        ∂x¹∂ξ¹[i, j, k] += ∂ξ¹[i, ii] * x¹[ii, j, k]
        ∂x²∂ξ¹[i, j, k] += ∂ξ¹[i, ii] * x²[ii, j, k]
        ∂x³∂ξ¹[i, j, k] += ∂ξ¹[i, ii] * x³[ii, j, k]
    end
    for jj = 1:M
        ∂x¹∂ξ²[i, j, k] += ∂ξ²[j, jj] * x¹[i, jj, k]
        ∂x²∂ξ²[i, j, k] += ∂ξ²[j, jj] * x²[i, jj, k]
        ∂x³∂ξ²[i, j, k] += ∂ξ²[j, jj] * x³[i, jj, k]
    end
    for kk = 1:N
        ∂x¹∂ξ³[i, j, k] += ∂ξ³[k, kk] * x¹[i, j, kk]
        ∂x²∂ξ³[i, j, k] += ∂ξ³[k, kk] * x²[i, j, kk]
        ∂x³∂ξ³[i, j, k] += ∂ξ³[k, kk] * x³[i, j, kk]
    end
    jacobian[1, 1, i, j, k] = ∂x¹∂ξ¹[i, j, k]
    jacobian[2, 1, i, j, k] = ∂x²∂ξ¹[i, j, k]
    jacobian[3, 1, i, j, k] = ∂x³∂ξ¹[i, j, k]

    jacobian[1, 2, i, j, k] = ∂x¹∂ξ²[i, j, k]
    jacobian[2, 2, i, j, k] = ∂x²∂ξ²[i, j, k]
    jacobian[3, 2, i, j, k] = ∂x³∂ξ²[i, j, k]

    jacobian[1, 3, i, j, k] = ∂x¹∂ξ³[i, j, k]
    jacobian[2, 3, i, j, k] = ∂x²∂ξ³[i, j, k]
    jacobian[3, 3, i, j, k] = ∂x³∂ξ³[i, j, k]

    detjacobian[i, j, k] = det(jacobian[:, :, i, j, k])
    ijacobian[:, :, i, j, k] = inv(jacobian[:, :, i, j, k])

end

# Kopriva modification to compute Jaⁱ
# ξ derivative in kopriva is ξ¹
# η derivative in kopriva is ξ²
# ζ derivative is ξ³
# need all threads to finish since we loop over different ∂xⁱ∂ξʲ combos
for i = 1:L, j = 1:M, k = 1:N
    for ii in 1:L
        # Ja², ξ derivatives
        kopriva[1, 2, i, j, k] += ∂ξ¹[i, ii] * ∂x²∂ξ³[ii, j, k] * x³[ii, j, k]
        kopriva[2, 2, i, j, k] += ∂ξ¹[i, ii] * ∂x³∂ξ³[ii, j, k] * x¹[ii, j, k]
        kopriva[3, 2, i, j, k] += ∂ξ¹[i, ii] * ∂x¹∂ξ³[ii, j, k] * x²[ii, j, k]
        # Ja³, ξ derivatives
        kopriva[1, 3, i, j, k] -= ∂ξ¹[i, ii] * ∂x²∂ξ²[ii, j, k] * x³[ii, j, k]
        kopriva[2, 3, i, j, k] -= ∂ξ¹[i, ii] * ∂x³∂ξ²[ii, j, k] * x¹[ii, j, k]
        kopriva[3, 3, i, j, k] -= ∂ξ¹[i, ii] * ∂x¹∂ξ²[ii, j, k] * x²[ii, j, k]
    end

    for jj in 1:M
        # Ja², η derivatives
        kopriva[1, 3, i, j, k] += ∂ξ²[j, jj] * ∂x²∂ξ¹[i, jj, k] * x³[i, jj, k]
        kopriva[2, 3, i, j, k] += ∂ξ²[j, jj] * ∂x³∂ξ¹[i, jj, k] * x¹[i, jj, k]
        kopriva[3, 3, i, j, k] += ∂ξ²[j, jj] * ∂x¹∂ξ¹[i, jj, k] * x²[i, jj, k]
        # Ja³, η derivatives
        kopriva[1, 1, i, j, k] -= ∂ξ²[j, jj] * ∂x²∂ξ³[i, jj, k] * x³[i, jj, k]
        kopriva[2, 1, i, j, k] -= ∂ξ²[j, jj] * ∂x³∂ξ³[i, jj, k] * x¹[i, jj, k]
        kopriva[3, 1, i, j, k] -= ∂ξ²[j, jj] * ∂x¹∂ξ³[i, jj, k] * x²[i, jj, k]
    end

    for kk in 1:N
        # Ja², ζ derivatives
        kopriva[1, 1, i, j, k] += ∂ξ³[k, kk] * ∂x²∂ξ²[i, j, kk] * x³[i, j, kk]
        kopriva[2, 1, i, j, k] += ∂ξ³[k, kk] * ∂x³∂ξ²[i, j, kk] * x¹[i, j, kk]
        kopriva[3, 1, i, j, k] += ∂ξ³[k, kk] * ∂x¹∂ξ²[i, j, kk] * x²[i, j, kk]
        # Ja³, ζ derivatives
        kopriva[1, 2, i, j, k] -= ∂ξ³[k, kk] * ∂x²∂ξ¹[i, j, kk] * x³[i, j, kk]
        kopriva[2, 2, i, j, k] -= ∂ξ³[k, kk] * ∂x³∂ξ¹[i, j, kk] * x¹[i, j, kk]
        kopriva[3, 2, i, j, k] -= ∂ξ³[k, kk] * ∂x¹∂ξ¹[i, j, kk] * x²[i, j, kk]
    end
end

if scale == 1
    volume_numerical = sum(detjacobian .* ω¹ .* ω² .* ω³)
    volume_exact = 4 / 3 * π * (1^3 - 0.8^3) / 6 # 1/6 since one of the faces
    volume_error = abs(abs(volume_numerical) - volume_exact) / abs(volume_exact)
    println("The relative error in computing the volume is ", volume_error)
end

##
# Check Free-Stream
divu = zeros(size(x¹))
c = [1, 1, 1] # vector field in cartesian space
# rows of ijacobian[:,:, i,j,k] are the contravariant basis vectors
for i = 1:L, j = 1:M, k = 1:N
    for ii = 1:L
        local u¹ = dot(ijacobian[1, :, ii, j, k], c)
        divu[i, j, k] += ∂ξ¹[i, ii] * detjacobian[ii, j, k] * u¹
    end
    for jj = 1:M
        local u² = dot(ijacobian[2, :, i, jj, k], c)
        divu[i, j, k] += ∂ξ²[j, jj] * detjacobian[i, jj, k] * u²
    end
    for kk = 1:N
        local u³ = dot(ijacobian[3, :, i, j, kk], c)
        divu[i, j, k] += ∂ξ³[k, kk] * detjacobian[i, j, kk] * u³
    end
end
divu = divu ./ detjacobian
divu_error = maximum(abs.(divu))
println("The error in computing the free-stream property is ", divu_error)

##
# Now use Kopriva metrics
divu = zeros(size(x¹))
c = [1, 1, 1] # vector field in cartesian space
# rows of ijacobian[:,:, i,j,k] are the contravariant basis vectors
for i = 1:L, j = 1:M, k = 1:N
    for ii = 1:L
        local Ju¹ = dot(kopriva[:, 1, ii, j, k], c)
        divu[i, j, k] += ∂ξ¹[i, ii] * Ju¹
    end
    for jj = 1:M
        local Ju² = dot(kopriva[:, 2, i, jj, k], c)
        divu[i, j, k] += ∂ξ²[j, jj] * Ju²
    end
    for kk = 1:N
        local Ju³ = dot(kopriva[:, 3, i, j, kk], c)
        divu[i, j, k] += ∂ξ³[k, kk] * Ju³
    end
end
divu = divu ./ detjacobian
divu_error = maximum(abs.(divu))
println("The error in computing the free-stream property with kopriva is ", divu_error)

##
# Compare metrics
index = (1, 1, 1)
detjacobian[index...] .* ijacobian[:, :, index...]'
kopriva[:, :, index...]

# Jacobian inflation
inflation_factor = det(kopriva[:, :, index...] ./ detjacobian[index...]) * detjacobian[index...]
println("The jacobian has inflated by a factor, ", inflation_factor)

##
# Could also construct metrics like this
II = 0 * ∂ξ¹ + I
tmp3 = kron(∂ξ¹, II, II)
tmp2 = kron(II, ∂ξ¹, II)
tmp1 = kron(II, II, ∂ξ¹)
div = [tmp1 tmp2 tmp3]
incompressible_subspace = nullspace(div)
inc_sub = incompressible_subspace
maximum(abs.(div * inc_sub))

Fˣ = (ξ¹.+0*ξ².+0*ξ³)[:]
Fʸ = (0*ξ¹.+1*ξ².+0*ξ³)[:]
Fᶻ = (0*ξ¹.+0*ξ².+1*ξ³)[:]
F⃗ = [Fˣ; Fʸ; Fᶻ]
divF = div * F⃗
all(divF .≈ 3.0)
maximum(abs.(tmp3 * Fˣ))
# need to project all throw
tmp = zeros(L, M, N, 3, 3)
ktmp = zeros(L, M, N, 3, 3)
aftertmp = zeros(L, M, N, 3, 3)

for c in 1:3, s in 1:3, i in 1:L, j in 1:M, k in 1:N
    tmp[i, j, k, c, s] = detjacobian[i, j, k] * ijacobian[s, c, i, j, k]
    ktmp[i, j, k, c, s] = kopriva[c, s, i, j, k]
end

chol_sub = cholesky((inc_sub' * inc_sub)) # this is just the identity matrix
# nullspace yields thing that are already orthonormalized (apparantely)
for c in 1:3
    flattened_tmp = tmp[:, :, :, c, :][:]
    projected_tmp = inc_sub * (chol_sub \ (inc_sub' * flattened_tmp))
    vtmp = view(tmp, :, :, :, c, :)
    vtmp[:] .= projected_tmp
    for s in 1:3, i in 1:L, j in 1:M, k in 1:N
        aftertmp[i, j, k, c, s] = vtmp[i, j, k, s]
    end
end

aftertmp[index..., :, :,]
kopriva[:, :, index...]
detjacobian[index...] .* ijacobian[:, :, index...]'
# orthogonal project should be closer to where we start
norm(aftertmp[index..., :, :]' - detjacobian[index...] .* ijacobian[:, :, index...])
norm(kopriva[:, :, index...]' - detjacobian[index...] .* ijacobian[:, :, index...])
kerr1 = div * (ktmp[:, :, :, 1, :][:])
kerr2 = div * (ktmp[:, :, :, 2, :][:])
kerr3 = div * (ktmp[:, :, :, 3, :][:])

println("Kopriva free-stream yields ", norm(kerr1, Inf))
println("Kopriva free-stream yields ", norm(kerr2, Inf))
println("Kopriva free-stream yields ", norm(kerr3, Inf))

err1 = div * (aftertmp[:, :, :, 1, :][:])
err2 = div * (aftertmp[:, :, :, 2, :][:])
err3 = div * (aftertmp[:, :, :, 3, :][:])

println("Projection free-stream yields ", norm(err1, Inf))
println("Projection free-stream yields ", norm(err2, Inf))
println("Projection free-stream yields ", norm(err3, Inf))

# P = inc_sub * inc_sub'# this is indeed a projection operator, eigvals = {1, 0}
# λP = eigvals(P)

max_stars_repo_name: sandreza/HeldSuarezVisualizationScripts

max_issues_repo_name: sandreza/HeldSuarezVisualizationScripts

max_forks_repo_name: sandreza/HeldSuarezVisualizationScripts


text:
# Generates the Codata module automatically from NIST data available at
# http://physics.nist.gov/cuu/Constants/Table/allascii.txt
module CODATAgen

export GenCodataDefault,GenCodataSIUnits,autogen_codata_module
const units_match=r"(mol|kg|m|Pa|W|MeV|GeV|eV|J|u|ohm|C|F|S|N|T)(\^\-?[0-9])?"
const ignore_unit=r"E_h|C_90"

function parse_to_codata(name,value,units,unc)
    if name == "bohr_magneton_in_ev_per_t"
        "# FIXME $name\n"
    elseif name == "boltzmann_constant_in_ev_per_k"
        "# FIXME $name\n"
    elseif name == "natural_unit_of_action_in_ev_s"
        "# FIXME $name\n"
    elseif name == "nuclear_magneton_in_ev_per_t"
        "# FIXME $name\n"
    elseif name == "planck_constant_in_ev_s"
        "# FIXME $name\n"
    elseif name == "planck_constant_over_2_pi_in_ev_s"
        "# FIXME $name\n"
    else
        if contains(units, "E_h")
            "# FIXME $name ($units)\n"
        elseif contains(units, "sr")
            "# FIXME $name ($units)\n"
        elseif contains(units, "Wb")
            "# FIXME $name ($units)\n"
        elseif contains(units, "/c")
            "# FIXME $name ($units)\n"
        elseif units == "ohm"
            units = "Ohm"
            "const $name = PhysicalConstant(\"$name\",$value$units,$unc$units,\"$units\")\n"
        elseif contains(units, "C_90")
            "# FIXME $name ($units)\n"
        else
            "const $name = PhysicalConstant(\"$name\",$value$units,$unc$units,\"$units\")\n"
        end
    end
end
function getTmpName()
    let makeTmpName=x->ccall((:tmpnam,"libc"),Ptr{Cchar},(Ptr{Uint8},),x)
        tmpname=makeTmpName(0)
        bytestring(tmpname)
    end
end
function get_codata_and_gen_module(write_exports=true,
                                   gen_name=parse_to_codata,
                                   write_header=write_header_default,
                                   name="CODATA.jl")
    tmpName=getTmpName()
    download("http://physics.nist.gov/cuu/Constants/Table/allascii.txt",tmpName)
    autogen_codata_module(tmpName,write_exports,gen_name,write_header,name)
end
function autogen_codata_module(filename :: String,write_exports=true,
                               gen_name=parse_to_codata,
                               write_header=write_header_default,
                               name="Codata.jl")
    # open the input file file and get all the lines.
    ifile = open(filename)
    lines = readlines(ifile)
    close(ifile)

    # now open the target for writing.
    exports = IOBuffer()
    consts = IOBuffer()
    ofile = open(name,false, true, true, true, false)
    # write the header and footer on the file
    write_header(ofile)

    # first 10 lines are the header, who cares.
    data = lines[11:end]

    # now we just iterate over the dataset.  a few times, unfortunately.
    for l in data
        name = strip(l[1:58])
        value = strip(l[59:83])
        unc = strip(l[84:107])
        units = strip(l[108:end])

        # replace spaces in name with underscores, and a bunch of other stuff.
        name = replace(name, ' ', '_')
        name = replace(name, ',', '_')
        name = replace(name, '-', '_')
        name = replace(name, '/', "_per_")
        name = replace(name, '.', "")
        name = replace(name, "{220}_lattice_spacing_of_silicon",
                       "si_220_lattice_spacing")
        name = replace(name, '(', '_')
        name = replace(name, ')', "")
        name = lowercase(name)

        # replace space in values and uncertainties with nothing
        value, unc = replace(value, ' ', ""), replace(unc, ' ', "")
        unc = replace(unc, "(exact)", "0.0")
        value = replace(value,"...","")
        units = replace(replace(units, "  ", ' '), ' ', '*')
        if(gen_name(name,value,units,unc) != Nothing())
            write(consts, gen_name(name,value,units,unc))
        end
    end
    if(write_exports)
    write(ofile,takebuf_string(exports))
    end
    close(exports)
    write(ofile,takebuf_string(consts))
    close(consts)
    # write the footer
    write_footer(ofile)

    # all done
    close(ofile)
end

function write_header_default(hndl :: IOStream)
    write(hndl, "# Codata.jl\n")
    write(hndl, "# This file has been automatically generated.\n\n")
    write(hndl, "module CODATA\n")
    write(hndl, "using SIUnits\n\n")
    write(hndl, "using SIUnits.ShortUnits\n\n")
    write(hndl, "# The PhysicalConstant type encapsulates all of the\n")
    write(hndl, "# information about a particular constant - the value,\n")
    write(hndl, "# the uncertainty in that value, the units, and so on.\n")
    write(hndl, "type PhysicalConstant\n")
    write(hndl, "\tname :: String\n")
    write(hndl, "\tvalue\n")
    write(hndl, "\tuncertainty\n")
    write(hndl, "\tunits_string :: String\n")
    write(hndl, "end\n\n")
    write(hndl, "# autogenerated non-SI units\n")
    write(hndl, "const MeV = SIUnits.Mega*eV # Mega eV\n")
    write(hndl, "const GeV = SIUnits.Giga*eV # Giga eV\n")
    write(hndl, "const fm = SIUnits.Femto*m # femtometers\n")
    write(hndl, "export AtomicMassUnit\n")
    write(hndl, "const u = SIUnits.NonSIUnit{typeof(kg), :u}() # AMU\n")
    write(hndl, "convert(::Type{SIUnits.SIQuantity},::typeof(u))=9.1093826kg\n ")
    write(hndl, "const T = V*s*m^2\n # Tesla")
    write(hndl, "# begin constants\n\n")
end

function write_footer(hndl :: IOStream)
    write(hndl, "end # end of module definition\n")
end
function GenCodataDefault()
    autogen_codata_module("allascii.txt")
end
function GenCodataSIUnits()
    autogen_codata_module("allascii.txt",true,parse_to_SIUnits,write_header_SIUnits)
end
end


text:
# Dense example
using DynamicPolynomials
using SumOfSquares
using MosekTools

@polyvar x[1:2]

m = SOSModel(Mosek.Optimizer)
@variable m t
@objective m Max t
@variable m sos SOSPoly([1, x[1], x[2]])
@constraint m sum(x.^2) - t == sos
optimize!(m)

# Exploit Symmetry
using SymbolicWedderburn
using PermutationGroups
using Cyclotomics

msym = SOSModel(Mosek.Optimizer)
@variable   msym t
@objective  msym Max t

C2 = PermGroup([perm"(1,2)"])

# @constraint msym sum(x.^2) - t in SOSCone() symmetry_group = C2

using MultivariateBases
using MultivariatePolynomials
const MP = MultivariatePolynomials

mvec = monomials(x, 0:1)
dict = Dict(m => i for (i, m) in enumerate(mvec))

function Base.:^(m::MP.AbstractMonomial, p::PermutationGroups.Perm)
    return prod(MP.variables(m).^(MP.exponents(m)^p))
end

G = PermGroup([Perm([dict[m] for m in mvec.^Ref(g)]) for g in gens(C2)])
ccG = conjugacy_classes(G)
chars = SymbolicWedderburn.characters_dixon(ccG)

c1, U1 = SymbolicWedderburn.central_projection(chars[1])
c2, U2 = SymbolicWedderburn.central_projection(chars[2])

#Int64 will not work for all Groups
R1, ids1 = SymbolicWedderburn.row_echelon_form(Int64.(U1))
R2, ids2 = SymbolicWedderburn.row_echelon_form(Int64.(U2))

@variable msym sos1 SOSPoly(FixedPolynomialBasis(R1[1:length(ids1),:]*mvec))
@variable msym sos2 SOSPoly(FixedPolynomialBasis(R2[1:length(ids2),:]*mvec))
@constraint msym sum(x.^2) - t == sos1 + sos2
optimize!(msym)

# do we scale or not?

# if we have a complex character, we need to find its conjugate
#


text:
"""
    pressure(mom::JuSwalbe.Macroquant; γ = 0.01, θ = 1/9)

Film pressure of the thin film equation.

This is just the summation of the laplacian of the height `Δh` and the disjoining potential `Π`.
Since this function uses Π one has to be careful of θ as it should be an array.
However it can be either a single element array or as large as the whole domain.
See also: [`Δh`](@ref), [`Π`](@ref)

# Math
The film pressure is simply

`` p_{film} = -\\gamma\\Delta h + \\Pi(h) ``

The sign is not 100% fixed and can change from paper to paper.
Here I choose to be in agreement with Thiele et al., which is known to be a good theoretician.

# Example
```jldoctest
julia> using JuSwalbe

julia> height = reshape([i for i in 1.0:1.0:16.0],4,4)
4×4 Array{Float64,2}:
 1.0  5.0   9.0  13.0
 2.0  6.0  10.0  14.0
 3.0  7.0  11.0  15.0
 4.0  8.0  12.0  16.0

julia> p = pressure(height)
4×4 Array{Float64,2}:
 -0.200016  -0.0400001   -0.04        0.12
 -0.160002  -7.44536e-8  -1.6082e-8   0.16
 -0.160001  -4.68862e-8  -1.20826e-8  0.16
 -0.12       0.04         0.04        0.2

julia> mom = JuSwalbe.Macroquant(height=height, velocity=JuSwalbe.Twovector(x=zeros(4,4), y=zeros(4,4)), pressure=zeros(4,4), energy=zeros(4,4))
JuSwalbe.Macroquant{Array{Float64,2},JuSwalbe.Twovector{Array{Float64,2}}}
  height: Array{Float64}((4, 4)) [1.0 5.0 9.0 13.0; 2.0 6.0 10.0 14.0; 3.0 7.0 11.0 15.0; 4.0 8.0 12.0 16.0]
  velocity: JuSwalbe.Twovector{Array{Float64,2}}
  pressure: Array{Float64}((4, 4)) [0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0]
  energy: Array{Float64}((4, 4)) [0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0]

julia> mom.pressure
4×4 Array{Float64,2}:
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

julia> pressure(mom) # Pressure is stored in the moments struct!
4×4 Array{Float64,2}:
 -0.200016  -0.0400001   -0.04        0.12
 -0.160002  -7.44536e-8  -1.6082e-8   0.16
 -0.160001  -4.68862e-8  -1.20826e-8  0.16
 -0.12       0.04         0.04        0.2

julia> mom.pressure
4×4 Array{Float64,2}:
 -0.200016  -0.0400001   -0.04        0.12
 -0.160002  -7.44536e-8  -1.6082e-8   0.16
 -0.160001  -4.68862e-8  -1.20826e-8  0.16
 -0.12       0.04         0.04        0.2

```
So what does this mean?
A positive pressure is force that drives the film down in height.
While a negative pressure generates a flux towards that location.
In the example a linear increasing height field was used, an equilibrium though is reached when the film is flat.

# References
## Recent (short)
- [Signatures of slip in dewetting polymer films](https://www.pnas.org/content/116/19/9275)

## Review
- [Dynamics and stability of thin liquid films](https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.81.1131)
"""
function pressure(mom::JuSwalbe.Macroquant{Matrix{T}, JuSwalbe.Twovector{Matrix{T}}}; γ::T=T(0.01), θ::Matrix{T}=ones(T,(1,1))*T(1/9)) where {T<:Number}
    # Dummy array to store the result
    p = zeros(T, size(mom.height))
    # All calculation needed here
    p = -γ * Δh(mom) .+ Π(mom, γ=γ, θ=θ)
    # Write the solution into the struct as well
    mom.pressure = p

    return p
end

function pressure(height::Array{T,2}; γ::T=0.01, θ::Matrix{T}=ones(T,(1,1))*T(1/9)) where {T<:Number}
    # Dummy array to store the result
    p = zeros(size(height))
    # All calculation needed here
    p = -γ * Δh(height) .+ Π(height, γ=γ, θ=θ)

    return p
end
# One dimensional implementation
function pressure(mom::JuSwalbe.Macroquant{Vector{T}, Vector{T}}; γ::T=T(0.01), θ::Vector{T}=ones(T,1)*T(1/9)) where {T<:Number}
    # Dummy array to store the result
    p = zeros(T, length(mom.height))
    # All calculation needed here
    p = -γ * Δh(mom) .+ Π(mom, γ=γ, θ=θ)
    # Write the solution into the struct as well
    mom.pressure = p

    return p
end

function pressure(height::Array{T,1}; γ::T=0.01, θ::Vector{T}=ones(T,1)*T(1/9)) where {T<:Number}
    # Dummy array to store the result
    p = zeros(size(height))
    # All calculation needed here
    p = -γ * Δh(height) .+ Π(height, γ=γ, θ=θ)

    return p
end

"""
    ∇p(mom::JuSwalbe.Macroquant)

Computation of the pressure gradient in one and two spatial dimensions.

Given a struct `Macroquant` with an `.pressure` entry performs a finite difference approximation of the spatial derivative.
A five-point stencil is used for one spatial dimension.
For two spatial dimensions a weighted nine-point stencil is applied.
An perodicity around the edges is assumed (ring, torus).

# Math
Since this is essentially discrete operators I will not go into details.
## Two spatial dimensions

`` \\nabla p = 3 \\sum_{i=0}^8 w_i \\mathbf{c}_i p(\\mathbf{r}+\\mathbf{c}_i) +O(\\nabla^3). ``

Where `` \\mathbf{c}_i, w_i `` are the standard D2Q9 lattice velocities and weights.
Of course this vector `` \\nabla p `` needs to be multiple with the height field to yield the correct force.

## One spatial dimension

`` \\partial_x p = \\frac{p(x-2\\Delta x) - 8p(x-1\\Delta x) + 0p(x) + 8p(x+1\\Delta x) - p(x+2\\Delta x)}{12\\Delta x} + O(\\partial_x^5). ``

Where the step size is limited by the resolution of the lattice Boltzmann grid `` \\Delta x ``.

# Example
## Two spatial dimensions
```jldoctest
julia> using JuSwalbe

julia> zerovelocity = JuSwalbe.Twovector(x=zeros(5,5), y=zeros(5,5)) # Thanks to @with_kw macro
JuSwalbe.Twovector{Array{Float64,2}}
  x: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  y: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]

julia> mom = JuSwalbe.Macroquant(height=fill(1.0, (5,5)), velocity=zerovelocity, pressure=zerovelocity.x, energy=zerovelocity.y)
JuSwalbe.Macroquant{Array{Float64,2},JuSwalbe.Twovector{Array{Float64,2}}}
  height: Array{Float64}((5, 5)) [1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0; … ; 1.0 1.0 … 1.0 1.0; 1.0 1.0 … 1.0 1.0]
  velocity: JuSwalbe.Twovector{Array{Float64,2}}
  pressure: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]
  energy: Array{Float64}((5, 5)) [0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0]

julia> mom.height[2,3] = 2.0; mom.height[3,2] = 2.0; mom.height[4,3] = 2.0; mom.height[3,4] = 2.0
2.0

julia> p = pressure(mom)
5×5 Array{Float64,2}:
 -1.6082e-5   -0.00168275  -0.00668275  -0.00168275  -1.6082e-5
 -0.00168275  -0.0133494    0.029998    -0.0133494   -0.00168275
 -0.00668275   0.029998    -0.0266827    0.029998    -0.00668275
 -0.00168275  -0.0133494    0.029998    -0.0133494   -0.00168275
 -1.6082e-5   -0.00168275  -0.00668275  -0.00168275  -1.6082e-5

julia> mom.pressure
5×5 Array{Float64,2}:
 -1.6082e-5   -0.00168275  -0.00668275  -0.00168275  -1.6082e-5
 -0.00168275  -0.0133494    0.029998    -0.0133494   -0.00168275
 -0.00668275   0.029998    -0.0266827    0.029998    -0.00668275
 -0.00168275  -0.0133494    0.029998    -0.0133494   -0.00168275
 -1.6082e-5   -0.00168275  -0.00668275  -0.00168275  -1.6082e-5

julia> f = JuSwalbe.Forces(slip=zerovelocity, h∇p=zerovelocity, thermal=zerovelocity, bathymetry=zerovelocity)
JuSwalbe.Forces{JuSwalbe.Twovector{Array{Float64,2}}}
  slip: JuSwalbe.Twovector{Array{Float64,2}}
  h∇p: JuSwalbe.Twovector{Array{Float64,2}}
  bathymetry: JuSwalbe.Twovector{Array{Float64,2}}
  thermal: JuSwalbe.Twovector{Array{Float64,2}}

julia> gradp = - ∇p(mom, f)
5×5×2 Array{Float64,3}:
[:, :, 1] =
  0.00166667    0.00097105  -0.0102825    0.00097105   0.00166667
  0.000137716  -0.00833802   0.00277309  -0.00833802   0.000137716
 -0.0          -0.0         -0.0         -0.0         -0.0
 -0.000137716   0.00833802  -0.00277309   0.00833802  -0.000137716
 -0.00166667   -0.00097105   0.0102825   -0.00097105  -0.00166667

[:, :, 2] =
  0.00166667   0.000137716  -0.0  -0.000137716  -0.00166667
  0.00097105  -0.00833802   -0.0   0.00833802   -0.00097105
 -0.0102825    0.00277309   -0.0  -0.00277309    0.0102825
  0.00097105  -0.00833802   -0.0   0.00833802   -0.00097105
  0.00166667   0.000137716  -0.0  -0.000137716  -0.00166667

julia> f.h∇p.x # The force comes with sign
5×5 Array{Float64,2}:
  0.00166667    0.00097105  -0.0102825    0.00097105   0.00166667
  0.000137716  -0.00833802   0.00277309  -0.00833802   0.000137716
 -0.0          -0.0         -0.0         -0.0         -0.0
 -0.000137716   0.00833802  -0.00277309   0.00833802  -0.000137716
 -0.00166667   -0.00097105   0.0102825   -0.00097105  -0.00166667

julia> f.h∇p.y # Same as for the x values
5×5 Array{Float64,2}:
  0.00166667   0.000137716  -0.0  -0.000137716  -0.00166667
  0.00097105  -0.00833802   -0.0   0.00833802   -0.00097105
 -0.0102825    0.00277309   -0.0  -0.00277309    0.0102825
  0.00097105  -0.00833802   -0.0   0.00833802   -0.00097105
  0.00166667   0.000137716  -0.0  -0.000137716  -0.00166667

```

## One spatial dimension
```jldoctest
julia> using JuSwalbe

julia> mom = JuSwalbe.Macroquant(height=ones(10), velocity=zeros(10), pressure=zeros(10), energy=zeros(10))
JuSwalbe.Macroquant{Array{Float64,1},Array{Float64,1}}
  height: Array{Float64}((10,)) [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
  velocity: Array{Float64}((10,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
  pressure: Array{Float64}((10,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
  energy: Array{Float64}((10,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

julia> f = JuSwalbe.Forces(slip=zeros(10), h∇p=zeros(10), thermal=zeros(10), bathymetry=zeros(10))
JuSwalbe.Forces{Array{Float64,1}}
  slip: Array{Float64}((10,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
  h∇p: Array{Float64}((10,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
  bathymetry: Array{Float64}((10,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
  thermal: Array{Float64}((10,)) [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

julia> mom.height[4] = 0.5; mom.height[5] = 2.0; mom.height[6] = 0.5;

julia> p = pressure(mom, θ=zeros(1), γ=1.0)
10-element Array{Float64,1}:
 -0.0
 -0.0
  0.5
 -2.0
  3.0
 -2.0
  0.5
 -0.0
 -0.0
 -0.0

julia> gradp = ∇p(mom, f)
10-element Array{Float64,1}:
 -0.041666666666666664
  0.5
 -1.5833333333333333
  0.9166666666666666
  0.0
 -0.9166666666666666
  1.5833333333333333
 -0.5
  0.041666666666666664
  0.0

```

# References
## Two spatial dimensions
- [Isotropic discrete Laplacian operators from lattice hydrodynamics](https://www.sciencedirect.com/science/article/pii/S0021999112004226)
- [Discretizations for the Incompressible Navier--Stokes Equations Based on the Lattice Boltzmann Method](https://epubs.siam.org/doi/abs/10.1137/S1064827599357188)

The second reference is the older and more mathematically one, while the first reference is short and great for a fast look up.

## One spatial dimension
- [Numerical Differentiation](http://www2.math.umd.edu/~dlevy/classes/amsc466/lecture-notes/differentiation-chap.pdf)

Any lecture notes or book on discrete mathematics should do for the one dimensional case.
"""
function ∇p(mom::JuSwalbe.Macroquant{Vector{T},Vector{T}}, force::JuSwalbe.Forces{Vector{T}}) where {T<:Number}
    # A (deep) copy of the pressure field
    p = deepcopy(mom.pressure)
    len = length(p)
    gradp = zeros(T, len)
    # Periodic extension of pressure, for higher order differential
    pushfirst!(p, p[len]) # Add last element on first spot
    pushfirst!(p, p[len]) # Add second last element on first spot
    push!(p, p[3])
    push!(p, p[4])
    ∇pout = zeros(T, length(p))
    # Five point first derivative
    for i in 3:len+2
        ∇pout[i] = T(1.0/12.0) * (p[i-2] + 8*p[i+1] - 8*p[i-1] - p[i+2])
    end
    gradp = ∇pout[3:len+2]
    # Write force to force struct
    force.h∇p = - mom.height .* gradp
    # Return result
    return  mom.height .* gradp
end

function ∇p(mom::JuSwalbe.Macroquant{Matrix{T},JuSwalbe.Twovector{Matrix{T}}}, force::JuSwalbe.Forces{JuSwalbe.Twovector{Matrix{T}}}) where {T<:Number}
    # A (deep) copy of the pressure field
    p = deepcopy(mom.pressure)
    width, thick = size(p)
    ∇poutx = zeros(T, (width, thick))
    ∇pouty = zeros(T, (width, thick))
    output = zeros(T, (width, thick, 2))
    # Periodic extension of pressure, for higher order differential
    wrapper = 1
    pperiodic = padarray(p, Pad(:circular, wrapper, wrapper))
    # Five point first derivative
    @inbounds for j in 1:thick, i in 1:width
        ∇poutx[i,j] = T(3.0) * (T(1.0/9.0) * (pperiodic[i+1,j] - pperiodic[i-1,j])
                              + T(1.0/36.0) * (pperiodic[i+1,j+1] - pperiodic[i-1,j+1] - pperiodic[i-1,j-1] + pperiodic[i+1,j-1]))
        ∇pouty[i,j] = T(3.0) * (T(1.0/9.0) * (pperiodic[i,j+1] - pperiodic[i,j-1] )
                              + T(1.0/36.0) * (pperiodic[i+1,j+1] + pperiodic[i-1,j+1] - pperiodic[i-1,j-1] - pperiodic[i+1,j-1]))
    end

    # Write force to force struct
    force.h∇p = JuSwalbe.Twovector(x = .- mom.height .* ∇poutx, y = .- mom.height .* ∇pouty)
    output[:, :, 1] = mom.height .* ∇poutx
    output[:, :, 2] = mom.height .* ∇pouty
    # Return result
    return  output
end

"""
    Π(mom::JuSwalbe.Macroquant; h_star = 0.1, exponents = [9,3], γ = 0.01, θ = 1.0/9.0)

Calculates the disjoining pressure for a given surface tension and contact angle at every lattice point.

The disjoing pressure potential enables the simulation to study phenomena like dewetting of thin liquid films.
Only with this potential it is possible to "dewett" without using highly sophisticated boundary conditions at `` h\\rightarrow 0 ``.
The functional shape of this potential can vary, ranging from the powerlaw used here to exponentials to a mixture of both.
Using the powerlaw shape the exponents (9,3) mimic the Lenard-Jones potential for the liquid substrate interaction.

# Arguments
- `mom::JuSwalbe.Macroquant` : Macroscopic moments of the simulation, important here `.height` which contains the height field. Allows for array input as well.
- `h_star::T` : Minimum of the pontential. T needs to be a subtype of `Number`!
- `exponents::Array{Int64,1}` : Exponents for the powerlaw potential
- `γ::T` : Surface tension
- `θ::T` : Equilibrium contact angle in multiple of π (for accuarcy reasons). Needs to be an array when doing patterning. For simple substrate use a one entry array with correct dimension.

# Math
The potential can be derived from the assumption that a given surface energy demands an equilibrium contact angle for the fluid.
For the exact derivation take a look at the references, in principle it is the derivative of the interfacial potential with respect to `h`.

`` \\Phi'(h) = \\Pi(h) ``

`` \\Pi(h) = \\kappa f(h) = (1-\\cos\\theta)\\frac{(n-1)(m-1)}{(n-m)h_{\\ast}}\\Big[\\Big(\\frac{h_{\\ast}}{h}\\Big)^n - \\Big(\\frac{h_{\\ast}}{h}\\Big)^m\\Big] ``

`` \\Pi(h_{\\ast}) = 0 ``

# Example
```jldoctest
julia> using JuSwalbe

julia> n,m = (4,4)
(4, 4)

julia> moment = JuSwalbe.Macroquant(height=ones(n,m), velocity=JuSwalbe.Twovector(x=zeros(n,m), y=zeros(n,m)), pressure=zeros(n,m), energy=zeros(n,m))
JuSwalbe.Macroquant{Array{Float64,2},JuSwalbe.Twovector{Array{Float64,2}}}
  height: Array{Float64}((4, 4)) [1.0 1.0 1.0 1.0; 1.0 1.0 1.0 1.0; 1.0 1.0 1.0 1.0; 1.0 1.0 1.0 1.0]
  velocity: JuSwalbe.Twovector{Array{Float64,2}}
  pressure: Array{Float64}((4, 4)) [0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0]
  energy: Array{Float64}((4, 4)) [0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0]

julia> p = Π(moment)
4×4 Array{Float64,2}:
 -1.6082e-5  -1.6082e-5  -1.6082e-5  -1.6082e-5
 -1.6082e-5  -1.6082e-5  -1.6082e-5  -1.6082e-5
 -1.6082e-5  -1.6082e-5  -1.6082e-5  -1.6082e-5
 -1.6082e-5  -1.6082e-5  -1.6082e-5  -1.6082e-5

julia> p = Π(moment, θ=zeros(1,1)) # Fully wetting substrate
4×4 Array{Float64,2}:
 -0.0  -0.0  -0.0  -0.0
 -0.0  -0.0  -0.0  -0.0
 -0.0  -0.0  -0.0  -0.0
 -0.0  -0.0  -0.0  -0.0

julia> h = reshape([i for i in 1.0:(n*m)],n,m)
4×4 Array{Float64,2}:
 1.0  5.0   9.0  13.0
 2.0  6.0  10.0  14.0
 3.0  7.0  11.0  15.0
 4.0  8.0  12.0  16.0

julia> p = Π(h)
4×4 Array{Float64,2}:
 -1.6082e-5   -1.28656e-7  -2.20603e-8  -7.31997e-9
 -2.01025e-6  -7.44536e-8  -1.6082e-8   -5.86078e-9
 -5.95628e-7  -4.68862e-8  -1.20826e-8  -4.76503e-9
 -2.51281e-7  -3.14101e-8  -9.30669e-9  -3.92626e-9

```

# References
To get a good understanding:
- [Long-scale evolution of thin liquid films](https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.69.931)
- [Wetting and spreading](https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.81.739)
- [Dynamics and stability of thin liquid films](https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.81.1131)

A rather recent new setup for the shape of `Π` can be found in
- [Signatures of slip in dewetting polymer films](https://www.pnas.org/content/116/19/9275)
"""
function Π(mom::JuSwalbe.Macroquant{Matrix{T}, JuSwalbe.Twovector{Matrix{T}}}; h_star::T=T(0.1), exponents=[9,3], γ::T=T(0.01), θ::Matrix{T}=ones(T,(1,1))*T(1/9)) where {T<:Number}
    # Theoretical minimum of the wetting pontential
    hbyhstar_n = ones(T, size(mom.height))
    hbyhstar_m = ones(T, size(mom.height))
    Π_h = zeros(T, size(mom.height))

    hbyhstar_n = power(h_star ./ mom.height, exponents[1])
    hbyhstar_m = power(h_star ./ mom.height, exponents[2])

    # Actual formular of the disjoining potential, long range attracion short range repulsion.
    Π_h = (γ .* (1 .- cospi.(θ))
          .* (exponents[1] - 1)*(exponents[2] - 1)/((exponents[1] - exponents[2])*h_star)
          .* (hbyhstar_n .- hbyhstar_m))

    return Π_h
end

function Π(height::Array{T,2}; h_star::T=T(0.1), exponents=[9,3], γ::T=T(0.01), θ::Matrix{T}=ones(T,(1,1))*T(1/9)) where {T<:Number}
    # Theoretical minimum of the wetting potential, two dimensional on array
    Π_h = zeros(T, size(height))

    hbyhstar_n = power(h_star ./ height, exponents[1])
    hbyhstar_m = power(h_star ./ height, exponents[2])
    # Actual formular of the disjoining potential, long range attracion short range repulsion.
    Π_h = (γ .* (1 .- cospi.(θ))
          .* (exponents[1] - 1)*(exponents[2] - 1) / ((exponents[1] - exponents[2])*h_star)
          .* (hbyhstar_n .- hbyhstar_m))

    return Π_h
end
# One dimensional implementation
function Π(mom::JuSwalbe.Macroquant{Vector{T}, Vector{T}}; h_star::T=T(0.1), exponents=[9,3], γ::T=T(0.01), θ::Vector{T}=ones(T,1)*T(1/9)) where {T<:Number}
    # Theoretical minimum of the wetting pontential
    Π_h = zeros(T, length(mom.height))

    hbyhstar_n = power(h_star ./ mom.height, exponents[1])
    hbyhstar_m = power(h_star ./ mom.height, exponents[2])
    # Actual formular of the disjoining potential, long range attracion short range repulsion.
    Π_h = (γ .* (1 .- cospi.(θ))
          .* (exponents[1] - 1)*(exponents[2] - 1) / ((exponents[1] - exponents[2])*h_star)
          .* (hbyhstar_n .- hbyhstar_m))

    return Π_h
end

function Π(height::Array{T,1}; h_star::T=T(0.1), exponents=[9,3], γ::T=T(0.01), θ::Vector{T}=ones(T,1)*T(1/9)) where {T<:Number}
    # Theoretical minimum of the wetting potential, one dimensional on array
    Π_h = zeros(T, length(height))

    hbyhstar_n = power(h_star ./ height, exponents[1])
    hbyhstar_m = power(h_star ./ height, exponents[2])
    # Actual formular of the disjoining potential, long range attracion short range repulsion.
    Π_h = (γ .* (1 .- cospi.(θ))
          .* (exponents[1] - 1)*(exponents[2] - 1) / ((exponents[1] - exponents[2])*h_star)
          .* (hbyhstar_n .- hbyhstar_m))

    return Π_h
end

"""
  power(array, n)

Computes the `n`th power elementwise of an `array`.

Turns out this is about **two** times faster than using `.^`, the build in power function.

# Example
```jldoctest
julia> using JuSwalbe

julia> h = fill(0.1, (5,5))
5×5 Array{Float64,2}:
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1
 0.1  0.1  0.1  0.1  0.1

julia> JuSwalbe.power(h, 3)
5×5 Array{Float64,2}:
 0.001  0.001  0.001  0.001  0.001
 0.001  0.001  0.001  0.001  0.001
 0.001  0.001  0.001  0.001  0.001
 0.001  0.001  0.001  0.001  0.001
 0.001  0.001  0.001  0.001  0.001

```
"""
function power(array::Array{T, 2}, n::Int) where {T<:Number}
  temp = ones(T, size(array))
  for i = 1:n
    temp .*= array
  end
  return temp
end
function power(array::Array{T, 1}, n::Int) where {T<:Number}
  temp = ones(T, length(array))
  for i = 1:n
    temp .*= array
  end
  return temp
end

"""
    Δh(mom::JuSwalbe.Macroquant)

Calculates the laplacian of the height field with periodic boundaries.

The calculation of the laplacian is central for the thin film evolution.
Only the pressure gradient will induce a flow, at least for the non-fluctuating version.
Therefore it is fairly important to have an accurate computation of the laplacian.

# Math
The Laplace equation is given by

`` \\Delta \\rho = 0 ``

The laplace operator is simply Δ = ∇ . ∇.
Since the film pressure has a contribution Δh we need a discreticed laplace operator.

## Two spatial dimensions

In two dimensions I follow the paper from Santosh and Succi.
A nine-point stencil is used and the equation goes as follow

`` \\Delta h_{i,j} = \\frac{1}{6}[4\\sum_{nn}h_{i,j} + \\sum_{diag}h_{i,j} - 20 h_{i,j}] ``

Where `i` and `j` are the x and y coordinates.
The index `nn` relates to the nearest neighbors, all four elements which are exactly Δx away from (i,j).
On the other hand the diagonal elements are those four which are have a distance √2Δx to (i,j).
Periodicity is taken care of by a circular padding of the height array.

## One spatial dimension

For a single spatial dimension it is mostly sufficient to use the central difference approach.
This approach is given by

`` \\partial_x^2 h_{i} = h_{i-1} - 2h_{i} + h_{i+1}

with `i` being the index along the spatial dimension.

# Example
```jldoctest
julia> using JuSwalbe

julia> height = reshape(collect(1.0:16.0), (4,4))
4×4 Array{Float64,2}:
 1.0  5.0   9.0  13.0
 2.0  6.0  10.0  14.0
 3.0  7.0  11.0  15.0
 4.0  8.0  12.0  16.0

julia> velx = zeros(Float64, (4,4))
4×4 Array{Float64,2}:
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

julia> vely = zeros(Float64, (4,4))
4×4 Array{Float64,2}:
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0

julia> vel = JuSwalbe.Twovector(x=velx, y=vely)
JuSwalbe.Twovector{Array{Float64,2}}
  x: Array{Float64}((4, 4)) [0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0]
  y: Array{Float64}((4, 4)) [0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0]

julia> moment = JuSwalbe.Macroquant(height=height, velocity=vel, pressure=zeros(4,4), energy=zeros(4,4))
JuSwalbe.Macroquant{Array{Float64,2},JuSwalbe.Twovector{Array{Float64,2}}}
  height: Array{Float64}((4, 4)) [1.0 5.0 9.0 13.0; 2.0 6.0 10.0 14.0; 3.0 7.0 11.0 15.0; 4.0 8.0 12.0 16.0]
  velocity: JuSwalbe.Twovector{Array{Float64,2}}
  pressure: Array{Float64}((4, 4)) [0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0]
  energy: Array{Float64}((4, 4)) [0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0]

julia> moment.height
4×4 Array{Float64,2}:
 1.0  5.0   9.0  13.0
 2.0  6.0  10.0  14.0
 3.0  7.0  11.0  15.0
 4.0  8.0  12.0  16.0

julia> laplace = Δh(moment)
4×4 Array{Float64,2}:
 20.0   4.0   4.0  -12.0
 16.0   0.0   0.0  -16.0
 16.0   0.0   0.0  -16.0
 12.0  -4.0  -4.0  -20.0

julia> height = collect(1.0:16.0)
16-element Array{Float64,1}:
  1.0
  2.0
  3.0
  4.0
  5.0
  6.0
  7.0
  8.0
  9.0
 10.0
 11.0
 12.0
 13.0
 14.0
 15.0
 16.0

julia> Δh(height)
16-element Array{Float64,1}:
  16.0
   0.0
   0.0
   0.0
   0.0
   0.0
   0.0
   0.0
   0.0
   0.0
   0.0
   0.0
   0.0
   0.0
   0.0
 -16.0

```

# References
## Two spatial dimensions
There are plenty of papers concerning discret differentail operators.
Many of them are good, although definitly not the best I go here with:
- [Isotropic discrete Laplacian operators from lattice hydrodynamics](https://www.sciencedirect.com/science/article/pii/S0021999112004226)

## One spatial dimension
Almost any reference on discrete differentiation is good enough.
- [Numerical Differentiation](http://www2.math.umd.edu/~dlevy/classes/amsc466/lecture-notes/differentiation-chap.pdf)

"""
function Δh(mom::JuSwalbe.Macroquant{Matrix{T}, Twovector{Matrix{T}}}) where {T<:Number}
    # Get the size of the problem
    width, thick = size(mom.height)
    # Build a wrap around to mimic periodic boundaries
    wrapper = 1
    toruswrapped = padarray(mom.height, Pad(:circular, wrapper, wrapper))
    # Generate a dummy array for the solution
    Δ = zeros(T, size(toruswrapped))
    # Exclude the wrapped boundaries
    for j in 1:width, i in 1:thick
        @inbounds Δ[i,j] = 1.0/6.0 * (4.0 * (toruswrapped[i+1, j] + toruswrapped[i-1, j] + toruswrapped[i, j+1] + toruswrapped[i, j-1])
                                     + 1.0 * (toruswrapped[i+1, j+1] + toruswrapped[i+1, j-1] + toruswrapped[i-1, j+1] + toruswrapped[i-1, j-1])
                                     - 20.0 * toruswrapped[i,j])
    end
    # Return the result without the wrapped boundaries
    return Δ[1:(end-2), 1:(end-2)]
end

function Δh(height::Array{T,2}) where {T<:Number}
    # Get the size of the problem
    width, thick = size(height)
    # Build a wrap around to mimic periodic boundaries
    wrapper = 1
    toruswrapped = padarray(height, Pad(:circular, wrapper, wrapper))
    # Generate a dummy array for the solution
    Δ = zeros(T, size(toruswrapped))
    # Exclude the wrapped boundaries
    for j in 1:width, i in 1:thick
        @inbounds Δ[i,j] = 1.0/6.0 * (4.0 * (toruswrapped[i+1, j] + toruswrapped[i-1, j] + toruswrapped[i, j+1] + toruswrapped[i, j-1])
                                     + 1.0 * (toruswrapped[i+1, j+1] + toruswrapped[i+1, j-1] + toruswrapped[i-1, j+1] + toruswrapped[i-1, j-1])
                                     - 20.0 * toruswrapped[i,j])
    end
    # Return the result without the wrapped boundaries
    return Δ[1:(end-2), 1:(end-2)]
end

# The one dimensonal cases
function Δh(mom::JuSwalbe.Macroquant{Vector{T}, Vector{T}}) where {T<:Number}
    # Get the size of the problem
    len = length(mom.height)
    # Build a wrap around to mimic periodic boundaries
    wrapper = 1
    # Ensures stability of height
    copyheight = deepcopy(mom.height)
    # Extend the one dimensional array periodic
    pushfirst!(copyheight, last(copyheight))
    push!(copyheight, copyheight[2])
    # Generate a dummy array for the solution
    Δ = zeros(T, size(copyheight))
    # Exclude the wrapped boundaries
    for i in 2:len+1
        @inbounds Δ[i] = copyheight[i-1] - 2*copyheight[i] + copyheight[i+1]
    end
    # Return the result without the wrapped boundaries
    return Δ[2:(end-1)]
end

function Δh(height::Array{T,1}) where {T<:Number}
    # Get the size of the problem
    len = length(height)
    # Build a wrap around to mimic periodic boundaries
    wrapper = 1
    # Ensures stability of height
    copyheight = deepcopy(height)
    # Extend the one dimensional array periodic
    pushfirst!(copyheight, last(copyheight))
    push!(copyheight, copyheight[2])
    # Generate a dummy array for the solution
    Δ = zeros(T, size(copyheight))
    # Exclude the wrapped boundaries
    for i in 2:len+1
        @inbounds Δ[i] = copyheight[i-1] - 2*copyheight[i] + copyheight[i+1]
    end
    # Return the result without the wrapped boundaries
    return Δ[2:(end-1)]
end

max_stars_repo_name: Zitzeronion/JuSwalbe

max_issues_repo_name: Zitzeronion/JuSwalbe


text:
include("$(@__DIR__)/../src/FiniteGroups.jl")
using .FiniteGroups, LinearAlgebra, Test
import .FiniteGroups: check_rep

for gi = 1:32
    g = pointgroup(gi)
    χ = charactertable(g)
    R = irreps(g, χ)
    @testset verbose = true "$(name(g)) Test" begin
        @testset "Character" begin
            # Test dimensions
            @test size(χ) == (length(class(g)), length(class(g)))
            # Test Orthogonality
            for i = 1:length(class(g)), j = 1:length(class(g))
                dt = sum(conj(χ[i,k]) * χ[j,k] * mult(g, k) for k = 1:length(class(g)))/order(g)
                dt -= (i==j)
                @test abs(dt) < 1e-5
            end
        end
        @testset "Irreps" begin
            for j = 1:size(χ, 1)
                rep = R[j]
                chi = χ[j, :]
                for k = 1:order(g)
                    @test chi[inclass(g, k)] ≈ tr(rep[k]) atol=1e-7     # check trace
                    @test norm(rep[k] * rep[k]' - I) < 1e-7             # check unitary
                end
                @test check_rep(g, rep)                                 #test multiplication
            end
        end
    end
end

for n = 2:5
    g = permutationgroup(n)
    χ = charactertable(g)
    R = irreps(g, χ)
    @testset verbose = true "$(name(g)) Test" begin
        @testset "Character" begin
            # Test dimensions
            @test size(χ) == (length(class(g)), length(class(g)))
            # Test Orthogonality
            for i = 1:length(class(g)), j = 1:length(class(g))
                dt = sum(conj(χ[i,k]) * χ[j,k] * mult(g, k) for k = 1:length(class(g)))/order(g)
                dt -= (i==j)
                @test abs(dt) < 1e-5
            end
        end
        @testset "Irreps" begin
            for j = 1:size(χ, 1)
                rep = R[j]
                chi = χ[j, :]
                for k = 1:order(g)
                    @test chi[inclass(g, k)] ≈ tr(rep[k]) atol=1e-7     # check trace
                    @test norm(rep[k] * rep[k]' - I) < 1e-7             # check unitary
                end
                @test check_rep(g, rep)                                 #test multiplication
            end
        end
    end
end


text:
function _circular_neighbor_offsets(npoints::Integer, radius::Real)
    radius >= 1 || throw(ArgumentError("radius >= 1.0 is expected."))

    ntuple(npoints) do p
        θ=2π*(p-1)/npoints
        pos = (-1, 1) .* radius .* sincos(θ)
        # rounding the positions for more stable result when applying interpolation
        round.(pos; digits=8)
    end
end

# A helper function to let `interpolation` keyword correctly understands `Degree` inputs
@inline wrap_BSpline(itp::Interpolations.InterpolationType) = itp
@inline wrap_BSpline(degree::Interpolations.Degree) = BSpline(degree)

@inline _inbounds_getindex(A, I) = @inbounds A[I...]
@inline _inbounds_getindex(A::Interpolations.AbstractInterpolation, I) = @inbounds A(I...)


text:
using HDF5

export save_hdf, load_hdf

function save_hdf(hdf_file_path, model; group_name::String="")
    h5open(hdf_file_path, "w") do file
        to_hdf(file, group_name, model)
    end
end

function load_hdf(hdf_file_path; group_name::String="")
    model = h5open(hdf_file_path, "r") do file
        matfac_from_hdf(file, group_name)
    end
    return model
end

####################################
# MatFacModel
####################################

function to_hdf(hdf_file, path::String, model::MatFacModel)

    # Factors
    write(hdf_file, string(path,"/X"), model.X)
    write(hdf_file, string(path,"/Y"), model.Y)

    # covariate coefficients
    if model.instance_covariate_coeff != nothing
        write(hdf_file, string(path,"/instance_covariate_coeff"),
                        model.instance_covariate_coeff)
    end
    to_hdf(hdf_file, string(path,"/instance_covariate_coeff_reg"),
                     model.instance_covariate_coeff_reg)

    # Offsets
    write(hdf_file, string(path,"/instance_offset"), model.instance_offset)
    write(hdf_file, string(path,"/feature_offset"), model.feature_offset)

    # Offset Regularizers
    to_hdf(hdf_file, string(path,"/instance_offset_reg"), model.instance_offset_reg)
    to_hdf(hdf_file, string(path,"/feature_offset_reg"), model.feature_offset_reg)

    # Precisions
    write(hdf_file, string(path,"/feature_precision"), model.feature_precision)

    # Losses
    loss_strs = String[string(typeof(loss)) for loss in model.losses]
    write(hdf_file, string(path, "/loss_types"), loss_strs)
    loss_scales = Float64[loss.scale for loss in model.losses]
    write(hdf_file, string(path, "/loss_scales"), loss_scales)

    # Instance regularizer matrices
    for (i, regmat) in enumerate(model.instance_reg_mats)
        to_hdf(hdf_file, string(path, "/inst_reg/mat_", i), regmat)
    end

    # Feature regularizer matrices
    for (i, regmat) in enumerate(model.feature_reg_mats)
        to_hdf(hdf_file, string(path, "/feat_reg/mat_", i), regmat)
    end

    write(hdf_file, string(path,"/factor_importances"), model.factor_importances)

end

function matfac_from_hdf(hdf_file, path::String)
    # Factors
    X = hdf_file[string(path,"/X")][:,:]
    Y = hdf_file[string(path,"/Y")][:,:]

    # Covariate coefficients
    if "instance_covariate_coeff" in keys(hdf_file[string(path,"/")])
        covariate_coeff = hdf_file[string(path,"/instance_covariate_coeff")][:,:]
    end
    covariate_coeff_reg = spmat_from_hdf(hdf_file, string(path,"/instance_covariate_coeff_reg"))

    # Offsets
    instance_offset = hdf_file[string(path,"/instance_offset")][:]
    feature_offset = hdf_file[string(path,"/feature_offset")][:]

    # Offset Regularizers
    instance_offset_reg = spmat_from_hdf(hdf_file, string(path, "/instance_offset_reg"))
    feature_offset_reg = spmat_from_hdf(hdf_file, string(path, "/feature_offset_reg"))

    # Precisions
    feature_precision = hdf_file[string(path,"/feature_precision")][:]

    # Losses
    loss_names = hdf_file[string(path,"/loss_types")][:]
    loss_scales = hdf_file[string(path,"/loss_scales")][:]
    losses = Loss[eval(Meta.parse(lname))(lscale) for (lname, lscale) in zip(loss_names, loss_scales)]

    # Instance regularizers
    inst_reg_gp = hdf_file[string(path, "/inst_reg")]
    inst_reg_mats = SparseMatrixCSC[spmat_from_hdf(inst_reg_gp, k) for k in sort(keys(inst_reg_gp))]

    # Feature regularizers
    feat_reg_gp = hdf_file[string(path, "/feat_reg")]
    feat_reg_mats = SparseMatrixCSC[spmat_from_hdf(inst_reg_gp, k) for k in sort(keys(feat_reg_gp))]

    # Factor importances
    factor_importances = hdf_file[string(path,"/factor_importances")][:]

    return MatFacModel(X, Y, inst_reg_mats, feat_reg_mats,
                             instance_offset, instance_offset_reg,
                             feature_offset, feature_offset_reg,
                             feature_precision,
                             covariate_coeff, covariate_coeff_reg,
                             losses, factor_importances)
end

##################################
# SparseMatrixCSC
##################################

function to_hdf(hdf_file, path::String, mat::SparseMatrixCSC)

    write(hdf_file, string(path, "/colptr"), mat.colptr)
    write(hdf_file, string(path, "/m"), mat.m)
    write(hdf_file, string(path, "/n"), mat.n)
    write(hdf_file, string(path, "/nzval"), mat.nzval)
    write(hdf_file, string(path, "/rowval"), mat.rowval)

end

function spmat_from_hdf(hdf_file, path::String)
    colptr = read(hdf_file, string(path, "/colptr"))
    m = read(hdf_file, string(path, "/m"))
    n = read(hdf_file, string(path, "/n"))
    nzval = read(hdf_file, string(path, "/nzval"))
    rowval = read(hdf_file, string(path, "/rowval"))

    return SparseMatrixCSC(m, n, colptr, rowval, nzval)
end


text:
## definition of the only built-in function in Aetherlang - USE
function lib_use(line::Ref{UInt}, ns::Namespace, arg1::AetherlangObject, arg2=AetherlangObject(""))::AetherlangObject
    lib::String = arg1.dataref
    include("Libs/"*lib*"/"*lib*".jl")
    nsmod::Namespace = eval(Meta.parse("$(lib)_namespace_modify"))
    prefix::String = arg2.dataref
    for (name, obj) in nsmod
        _fullname::String = prefix*name
        if !haskey(ns, _fullname)
            ns[_fullname] = obj
        else
            throw(AetherlangError("Library `$lib` tries to redefine an existing name `$_fullname`. Use prefix"))
        end
    end
    AetherlangObject()
end

max_stars_repo_name: Aetherlang/aetherlang

max_issues_repo_name: gojakuch/aetherlang

max_forks_repo_name: gojakuch/aetherlang


text:
function plot_results(lw::AbstractDataFrame)
    lw[!,:xname] = ["# $(Int64(round(r.rank_zscore_total))) $(r.name)" for r in eachrow(lw)]
    lw[!,:notey] = [(r.zscore_total >= 0.0 ? r.zscore_total + 0.5 : r.zscore_total - 0.5) for r in eachrow(lw)]
    lw[!,:note] =[(r.rank_zscore_total, r.notey, Plots.text("$(round(r.zscore_total,digits=1))",:center,10)) for r in eachrow(lw)]
    sort!(lw,:rank_zscore_total)
    add_colors!(lw)
    fig = plot(ylabel = "Zscore",
        xlabel = "Ranking",
        guidefont = (18, :darkgreen),
        bottom_margin = 100px,
        left_margin = 40px,
        size=(1000,1000))
    bar!(fig,lw.rank_zscore_total,lw.zscore_total,
        xticks = (1:nrow(lw),lw.xname),
        tickfont = 12,
        yticks = -20:1:20,
        xrotation = 45,
        ylims = (minimum(lw.zscore_total)-1.2,maximum(lw.zscore_total)+1.2),
        fillcolor = lw.heat,
        legend = false)
    annotate!(lw.note)
end

function plot_results(clan::AbstractDict; remove = ["player_names"])
    res = update_results(clan)
    lw = lastweek(res)
    filt = filter(row -> !in(row.name,remove),lw)
    p = AOM.plot_results(filt)
    filename = clan[:name]*"_"*string(today())*".png"
    dir = "/Users/dariosarra/Documents/House/AOM/Stats/"
    savefig(p,joinpath(dir,filename))
    return p
end

function add_colors!(lw::AbstractDataFrame)
    lw[!,:color] .= RGB(61/255,164/255,77/255)
    for r in eachrow(lw)
        if r.rank_zscore_total > 24
            r.color = RGB(172/255,141/255,24/255)
        elseif r.rank_zscore_total < 6
            r.color = RGB(0/255,153/255,249/255)
        end
    end
    pal = range(HSL(colorant"green"), stop=HSL(colorant"red"), length=nrow(lw))
    lw[!,:heat] = collect(pal)
end


text:
# http://en.wikipedia.org/wiki/Function_representation
#----------------------------------

"""
 Overloading CoordinateTransformations.Translation for vectorized implementation.
"""
function (trans::Translation{V})(x::AbstractMatrix) where {V}
    x .+ trans.translation
end

#function FRep(p::MapContainer{N,T,P}, v) where {N,T,P}
#    FRep(p.primitive, p.inv(v))
#end

function FRep(p::MapContainer{N,T,P}, v...) where {N,T,P}
    FRep(p.primitive, p.inv(vcat(v...)))
end
#----------------------------------

function FRep(u::CSGUnion, v...)
    min.(FRep(u.left, v...),FRep(u.right, v...))
end

function FRep(u::CSGDiff, v...)
    max.(FRep(u.left, v...), -FRep(u.right, v...))
end

function FRep(u::CSGIntersect, v...)
    max.(FRep(u.left, v...), FRep(u.right, v...))
end
#----------------------------------

function FRep(p::Sphere,v::AbstractVector)
    norm(v) - p.radius
end

function FRep(p::Sphere,x::AbstractMatrix,y::AbstractMatrix,z::AbstractMatrix)
    @assert size(x) == size(y) == size(z)
    @. sqrt(x*x + y*y + z*z) - p.radius
end

function FRep(p::Sphere,xyz::AbstractMatrix)
    @assert size(xyz)[1] == 3
    r2 = sum(xyz .* xyz,dims=1)
    @. sqrt(r2) - p.radius
end
#----------------------------------

function FRep(p::Cylinder,v::AbstractVector)
    x = v[1]
    y = v[2]
    z = v[3]
    max(-z+p.bottom, z-p.height-p.bottom, sqrt(x*x + y*y) - p.radius)
end

function FRep(p::Cylinder,x::AbstractMatrix,y::AbstractMatrix,z::AbstractMatrix)
    @assert size(x) == size(y) == size(z)
    r = @. sqrt(x*x + y*y) - p.radius
    a = @. -z+p.bottom
    b = @.  z-(p.height+p.bottom)
    max.(a, b, r)
end

function FRep(p::Cylinder,xyz::AbstractMatrix)
    @assert size(xyz)[1] == 3

    x = xyz[1,:]'
    y = xyz[2,:]'
    z = xyz[3,:]'
    fr = FRep(p,x,y,z)
end
#----------------------------------

function FRep(p::Circle, v::AbstractVector)
    norm(v) - p.radius
end

function FRep(p::Circle,x::AbstractMatrix,y::AbstractMatrix)
    @assert size(x) == size(y)
    @. sqrt(x*x + y*y) - p.radius
end

function FRep(p::Circle,xy::AbstractMatrix)
    @assert size(xy)[1] == 2
    r2 = sum(xy .* xy,dims=1)
    @. sqrt(r2) - p.radius
end
#----------------------------------

function FRep(p::Cuboid, v::AbstractVector)
    x = v[1]
    y = v[2]
    z = v[3]
    dx, dy, dz = p.dimensions
    lbx, lby,lbz = p.lowercorner
    max(-x+lbx, x-dx-lbx,
        -y+lby, y-dy-lby,
        -z+lbz, z-dz-lbz)
end

function FRep(p::Cuboid,x::AbstractMatrix,y::AbstractMatrix,z::AbstractMatrix)
    @assert size(x) == size(y) == size(z)
    dx, dy, dz = p.dimensions
    lbx, lby,lbz = p.lowercorner
    a = @. -x + lbx
    b = @.  x - (dx+lbx)
    c = @. -y + lby
    d = @.  y - (dy+lby)
    e = @. -z + lbz
    f = @.  z - (dz+lbz)
    max.(a,b,c,d,e,f)
end

function FRep(p::Cuboid,xyz::AbstractMatrix)
    @assert size(xyz)[1] == 3

    x = xyz[1,:]'
    y = xyz[2,:]'
    z = xyz[3,:]'
    fr = FRep(p,x,y,z)
end
#----------------------------------

function FRep(p::Square, v::AbstractVector)
    x = v[1]
    y = v[2]
    dx, dy = p.dimensions
    lbx, lby = p.lowercorner
    max(-x+lbx, x-dx-lbx,
        -y+lby, y-dy-lby)
end

function FRep(p::Square,x::AbstractMatrix,y::AbstractMatrix)
    @assert size(x) == size(y)
    dx, dy = p.dimensions
    lbx, lby = p.lowercorner
    a = @. -x + lbx
    b = @.  x - (dx+lbx)
    c = @. -y + lby
    d = @.  y - (dy+lby)

    max.(a,b,c,d)
end

function FRep(p::Square,xy::AbstractMatrix)
    @assert size(xy)[1] == 2

    x = xy[1,:]'
    y = xy[2,:]'
    fr = FRep(p,x,y)
end
#----------------------------------

function FRep(s::Shell, v)
    r = FRep(s.primitive, v)
    max.(r, -r .- s.distance)
end

#----------------------------------
function _radius(a::Real,b::Real,r::Real)
    if abs(a-b) >= r
        return min(a,b)
    else
        return b+r*sin(pi/4+asin((a-b)/(r*sqrt(2))))-r
    end
end

function _radius(a::AbstractVector,b::AbstractVector,c::AbstractVector)
    cond = @. (abs(a-b) >= r) # bit-vector
    ans1 = @. min(a,b)
    ans2 = @. b+r*sin(pi/4+asin((a-b)/(r*sqrt(2))))-r

    @. cond*ans1 + (1-cond)*ans2
end

function FRep(u::RadiusedCSGUnion, v)
    a = FRep(u.left, v)
    b = FRep(u.right, v)
    r = u.radius
    _radius(a,b,r)
end
#----------------------------------

function FRep(p::Piping{T}, v::AbstractVector) where {T}
    num_pts = length(p.points)

    val = typemax(T)

    for i = 1:num_pts-1
        e1 = p.points[i]
        e2 = p.points[i+1]
        v = e2 - e1
        w = v - e1
        if dot(w,v) <= 0
            nv = norm(v - e1)
        elseif dot(v,v) <= dot(w,v)
            nv = norm(v - e2)
        else
            nv = norm(cross(v-e1,v-e2))/norm(e2-e1)
        end
        val = min(nv, val)
    end
    val - p.radius
end

function FRep(p::Piping{T}, x::AbstractMatrix, y::AbstractMatrix, z::AbstractMatrix) where {T}

    xyz = vcat(x,y,z)
    collect(FRep(p,xyz[:,i]) for i in 1:size(xyz)[2])'
end

function FRep(p::Piping{T}, xyz::AbstractMatrix) where {T}
    @warn "F-Reps for Piping is not vectorized. switching to non-vectorized implementation"

    collect(FRep(p,xyz[:,i]) for i in 1:size(xyz)[2])'
end
#----------------------------------

function FRep(p::LinearExtrude, v::AbstractVector)
    x = v[1]
    y = v[2]
    z = v[3]
    r = FRep(p.primitive, [x,y])
    max.(max.(-z,z-p.distance), r)
end

function FRep(p::LinearExtrude,x::AbstractMatrix,y::AbstractMatrix,z::AbstractMatrix)
    @assert size(x) == size(y) == size(z)

    xy = vcat(x,y)
    r = FRep(p.primitive,xy)
    @. max(max(-z,z-p.distance), r)
end

function FRep(p::LinearExtrude,xyz::AbstractMatrix)
    @assert size(xyz)[1] == 3

    x = xyz[1,:]'
    y = xyz[2,:]'
    z = xyz[3,:]'
    fr = FRep(p,x,y,z)

end
#----------------------------------


text:
using Tk

const frameinterval = 0.12 # partial seconds between change on screen display

function windowanim(stepinterval::Float64)
    wind = Window("Animation", 300, 100)
    frm = Frame(wind)
    hello = "Hello World!                                           "
    but = Button(frm, width=30, text=hello)
    rightward = true
    callback(s) = (rightward = !rightward)
    bind(but, "command", callback)
    pack(frm, expand=true, fill = "both")
    pack(but, expand=true, fill = "both")
    permut = [hello[i:end] * hello[1:i-1] for i in length(hello)+1:-1:2]
    ppos = 1
    pmod = length(permut)
    while true
        but[:text] = permut[ppos]
        sleep(stepinterval)
        if rightward
            ppos += 1
            if ppos > pmod
                ppos = 1
            end
        else
            ppos -= 1
            if ppos < 1
                ppos = pmod
            end
        end
    end
end

windowanim(frameinterval)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
import StatsBase.sample
function softmax_greedy(X, Y, stepsize, W, n_epoch, topk, exact=1)
    # SGD
    N, D = size(X)
    K = size(Y, 1)

    #n_epoch = 50
    #W = zeros(D, K)
    #tic()
    #XW = X*W
    #@show(toc())
    XT = X'
    for epoch = 1:n_epoch
        #grad = zeros(D, K)
        #XW = X*W

        perm = randperm(N)
        tic()
        prec = 0.0
        prec_rand = 0.0
        for i = 1:1000
            grad = zeros(D, K)
            sample_idx = perm[i]
            class = find(Y[:,sample_idx])[1]
            Wx = W'*XT[:,sample_idx]
            #get topk
            if exact == 1
                #exact topk
                top_idx = sortperm(Wx, rev=true)[1:topk]
            elseif exact == 2
                # greedy mips topk
                maxWx = zeros(K)
                for j in 1:K
                    maxWx[j] = maximum( W[:,j] .* XT[:,sample_idx] )
                end
                top_idx = sortperm(maxWx, rev=true)[1:topk]
            else
                # random topk
                top_idx = sample([1:K;], topk)
            end

            # check greedy mips
            maxWx = zeros(K)
            for j in 1:K
                maxWx[j] = maximum( W[:,j] .* XT[:,sample_idx] )
            end
            top_idx_gmips = sortperm(maxWx, rev=true)[1:topk]
            top_idx_true = sortperm(Wx, rev=true)[1:topk]
            top_idx_rand = sample([1:K;], topk)
            #println( sortperm(Wx, rev=true) )
            #if i == 500
                #@show(top_idx_gmips)
                #@show(top_idx_true)
                #@show( intersect(top_idx_gmips, top_idx_true) )
            #end
            prec += length(intersect(top_idx_gmips, top_idx_true[1:10]) )/10
            prec_rand +=length(intersect(top_idx_rand, top_idx_true[1:10]) )/10

            for kk in 1:topk
                if class == top_idx[kk]
                    deleteat!(top_idx, kk)
                    break
                end
            end

            noise_sample = []
            while length(noise_sample) < 10
                noise = rand(1:K)
                if noise in top_idx || noise == class
                    continue
                else
                    append!(noise_sample, noise)
                end
            end

            tmp = exp.(Wx[top_idx])
            deno_estimate = sum(tmp) + exp(Wx[class]) + (K - length(top_idx) - 1)/length(noise_sample) * sum( exp.(Wx[noise_sample]) )

            # for j = 1:K
            #     grad[:,j] = tmp[j]/deno_estimate * XT[:,sample_idx]
            # end
            # grad[:, class] -= XT[:,sample_idx]

            #update
            for ii in top_idx
                W[:,ii] -= stepsize* exp(Wx[ii])/deno_estimate * XT[:,sample_idx]
            end

            W[:,class] -= stepsize*(exp(Wx[class])/deno_estimate - 1) * XT[:,sample_idx]

        end
        @show(toq())
        @printf "%2.1e epoch complete \n" epoch

        # calculate loss
        prec = prec / 1000
        prec_rand = prec_rand / 1000
        @printf "prec@10 for gmips: %2.6e, prec@10 for random: %2.6e \n" prec prec_rand
        loss = 0
        risk = 0.0
        XW = X*W
        for i = 1:N
            Wxi = XW[i,:]
            class = find(Y[:,i])[1]
            deno = sum( exp.(Wxi) )
            loss += -log( exp(Wxi[class]) / deno )

            ii = indmax(Wxi)
            if ii != class
                risk += 1
            end
        end
        loss = loss / N
        loss = loss / log(2)
        risk = risk / N
        @printf "loss: %3.3e, risk %3.3e \n" loss risk

    end
    return W
end

# softmax
function softmax(X, Y, stepsize, W, n_epoch, topk)
    # SGD
    N, D = size(X)
    K = size(Y, 1)

    #n_epoch = 50
    #W = zeros(D, K)
    #tic()
    #XW = X*W
    #@show(toc())
    XT = X'
    for epoch = 1:n_epoch
        #grad = zeros(D, K)
        #XW = X*W

        perm = randperm(N)
        tic()
        for i = 1:1000
            #println(i)
            grad = zeros(D, K)
            sample_idx = perm[i]
            class = find(Y[:,sample_idx])[1]
            Wx = W'*XT[:,sample_idx]
            tmp = exp.(Wx)
            deno = sum(tmp)
            for j = 1:K
                grad[:,j] = tmp[j]/deno * XT[:,sample_idx]
            end
            grad[:, class] -= XT[:,sample_idx]

            #update
            W -= stepsize*grad
        end
        @show(toq())
        @printf "%2.1e epoch complete \n" epoch

        # show topk statistics
        # randomly choose 3 samples
        for i = 1:5
            sample_idx = rand(1:N)
            class = find(Y[:,sample_idx])[1]
            Wx = W'*XT[:,sample_idx]
            tmp = exp.(Wx)
            grad = zeros(D, K)
            deno = sum(tmp)
            for j = 1:K
                grad[:,j] = tmp[j]/deno * XT[:,sample_idx]
            end
            grad[:, class] -= XT[:,sample_idx]

            norm_grad = zeros(K)
            for j = 1:K
                norm_grad[j] = norm( grad[:,j], 1 )
            end

            sorted_tmp = sort(norm_grad, rev=true)
            total = sum(sorted_tmp)
            total_topk = sum(sorted_tmp[1:topk])
            true_tmp = norm_grad[class]
            # @show(norm_grad[class])
            # @show(sum(norm_grad))
            # @printf "    top %i propotion: %1.4e\n" topk total_topk/total
            # @printf "    true class propotion: %1.4e\n" true_tmp/total
        end

        # calculate loss
        loss = 0
        risk = 0.0
        XW = X*W
        for i = 1:N
            Wxi = XW[i,:]
            class = find(Y[:,i])[1]
            deno = sum( exp.(Wxi) )
            loss += -log( exp(Wxi[class]) / deno )

            ii = indmax(Wxi)
            if ii != class
                risk += 1
            end
        end
        loss = loss / N
        loss = loss / log(2)
        risk = risk / N
        @printf "loss: %3.3e, risk %3.3e \n" loss risk

    end
    return W
end

# pairwise rank loss
function train_pairwise(X, Y, stepsize, n_epoch, W)
    N, D = size(X)
    K = size(Y, 1)

    #W = zeros(D, K)
    XT = X'
    for epoch in 1:n_epoch
        perm = randperm(N)
        tic()
        for i = 1:N
            rand_class = rand(1:K)
            sample_idx = perm[i]
            class = find(Y[:,sample_idx])[1]
            if rand_class == class
                continue
            end
            tmp_class = dot(W[:,class], XT[:, sample_idx])
            tmp_other = dot(W[:,rand_class], XT[:, sample_idx])
            grad = exp(tmp_other - tmp_class) / (1 + exp(tmp_other - tmp_class)) * XT[:, sample_idx]

            # update
            W[:, class] += stepsize * grad
            W[:, rand_class] -= stepsize * grad

        end

        @show(toq())
        @printf "%2.1e epoch complete \n" epoch

        # calculate loss
        loss = 0
        risk = 0.0
        XW = X*W
        for i = 1:N
            Wxi = XW[i,:]
            class = find(Y[:,i])[1]
            #deno = sum( exp.(Wxi) )
            #loss += -log( exp(Wxi[class]) / deno )
            tmp_class = dot(W[:,class], XT[:, i])
            #for j = 1:K
            #    if j == class
            #        continue
            #    end
            #    tmp_other = XW[i,j]
            #    loss += log( 1 + exp(tmp_other - tmp_class ) )
            #end

            ii = indmax(Wxi)
            if ii != class
                risk += 1
            end
        end
        loss = loss / N
        loss = loss / log(2)
        risk = risk / N
        @printf "loss: %3.3e, risk %3.3e \n" loss risk

    end

    return W
end

import StatsBase.sample
# negative sampling
function train_neg_sampling(X, Y, stepsize, n_epoch, W, k)
    N, D = size(X)
    K = size(Y, 1)

    #W = zeros(D, K)
    XT = X'
    # set # neg samples to be 5
    #k = 5
    for epoch in 1:n_epoch
        perm = randperm(N)
        tic()
        for i = 1:N
            sample_idx = perm[i]
            # choose negative samples
            class = find(Y[:,sample_idx])[1]
            neg_samples = sample([1:K;], k)
            grad = 1 / (1 + exp( dot(W[:,class], XT[:,sample_idx]) )) * XT[:, sample_idx]
            # update
            W[:,class] += stepsize*grad
            for ss in neg_samples
                if ss == class
                    continue
                end
                grad = 1 / (1 + exp(-dot(W[:,ss], XT[:,sample_idx]))) * XT[:, sample_idx]
                # update
                W[:,ss] -= stepsize*grad
            end

        end

        @show(toq())
        @printf "%2.1e epoch complete \n" epoch

        # calculate loss
        loss = 0
        risk = 0.0
        XW = X*W
        for i = 1:N
            Wxi = XW[i,:]
            class = find(Y[:,i])[1]
            #deno = sum( exp.(Wxi) )
            #loss += -log( exp(Wxi[class]) / deno )
            tmp_class = dot(W[:,class], XT[:, i])
            #for j = 1:K
            #    if j == class
            #        continue
            #    end
            #    tmp_other = XW[i,j]
            #    loss += log( 1 + exp(tmp_other - tmp_class ) )
            #end

            ii = indmax(Wxi)
            if ii != class
                risk += 1
            end
        end
        loss = loss / N
        loss = loss / log(2)
        risk = risk / N
        @printf "loss: %3.3e, risk %3.3e \n" loss risk

    end

    return W
end


text:
module TestComputeFP

using QuantEcon
using Base.Test
using FactCheck

# set up
mu_1 = 0.2  # 0 is unique fixed point forall x_0 \in [0, 1]

# (4mu - 1)/(4mu) is a fixed point forall x_0 \in [0, 1]
mu_2 = 0.3

# starting points on (0, 1)
unit_inverval = [0.1, 0.3, 0.6, 0.9]

# arguments for compute_fixed_point
kwargs = Dict{Symbol,Any}(:err_tol => 1e-5, :max_iter => 200,
                                  :verbose => true, :print_skip => 30)

rough_kwargs = Dict{Symbol,Any}(:atol => 1e-4)

T(x, mu) = 4.0 * mu * x * (1.0 - x)

# shorthand
abs_fp(f, i, other) = abs(compute_fixed_point(f, i; kwargs...) - other)

facts("Testing compute_fp.jl") do
    # convergence inside interval of convergence
    let f(x) = T(x, mu_1)
        for i in unit_inverval
            @fact abs_fp(f, i, 0.0) --> roughly(0.0; rough_kwargs...)
        end
    end

    # non convergence outside interval of convergence
    let f(x) = T(x, mu_2)
        for i in unit_inverval
            # none of these should converge to 0
            @fact abs_fp(f, i, 0.0) < 1e-4 --> false
        end
    end

    # convergence inside interval of convergence
    let f(x) = T(x, mu_2), fp = (4 * mu_2 - 1) / (4 * mu_2)
        for i in unit_inverval
            @fact abs_fp(f, i, fp) --> roughly(0.0; rough_kwargs...)
        end
    end

    # non convergence outside interval of convergence
    let f(x) = T(x, mu_1), fp = (4 * mu_1 - 1) / (4 * mu_1)
        for i in unit_inverval
            # none of these should converge to 0
            @fact abs_fp(f, i, fp) < 1e-4 --> false
        end
    end

end  # facts
end  # module


text:
movie_data = dataset("ggplot2", "movies")
function get_movie_minibatch(x, y, batch_size)
    x_out = [CuArray{Float32}([NaN for _ in 1:length(x[1, :])]) for _ in 1:batch_size]
    y_out = [NaN for _ in 1:batch_size]
    for i in 1:batch_size
        rn = rand(1:length(y))
        while (any(x->ismissing(x), x[rn, :]) || any(x->ismissing(x), y[rn]))
            rn = rand(1:length(y))
        end

        x_out[i] = CuArray{Float32}(collect(Iterators.flatten(x[rn, :])))
        y_out[i] = y[rn]
    end
    return x_out, y_out
end
convert_names(titles) = begin
    out = zeros(length(titles), maximum(titles .|> length))
    for i in eachindex(titles)
        for j in eachindex(titles[i])
            out[i, j] = titles[i][j] |> Int
        end
    end
    return out
end

TRAIN_PERCENTAGE = 70 # how much data will be used for training

titles = convert_names(movie_data.Title)
normalized_budget = movie_data.Budget ./ (sum(skipmissing(movie_data.Budget)) / length(movie_data.Budget))
x_data = cat(normalized_budget, movie_data.Length, titles; dims=2)
y_data = movie_data.Rating

# print number of rows with missing values in x_data
no_missing_rows = map((x->any(y->ismissing(y), x)), [x_data[i,:] for i in 1:length(x_data[:, 1])]) |> sum
println("number of rows with missing values: ", no_missing_rows)
println("usable number of rows: ", length(y_data) - no_missing_rows)

DATA_DIVIDION_INDEX = (length(x_data[:, 1])/100) * TRAIN_PERCENTAGE |> round |> Int
train_x = x_data[1:DATA_DIVIDION_INDEX, :]
train_y = y_data[1:DATA_DIVIDION_INDEX]
test_x = x_data[DATA_DIVIDION_INDEX+1:end, :]
test_y = y_data[DATA_DIVIDION_INDEX+1:end]


text:
# Load actual module (in case not installed)
pkgpath = joinpath(dirname(@__FILE__), "../src/")
push!(LOAD_PATH,pkgpath)
using GeneralQP

# Load test module
using Test

# actual tests
include("runtests.jl")


text:
using AbstractFFTs

for f in [
    :(AbstractFFTs.fft),
    :(AbstractFFTs.bfft),
    :(AbstractFFTs.ifft),
    :(AbstractFFTs.fft!),
    :(AbstractFFTs.bfft!),
    :(AbstractFFTs.ifft!),
    :(AbstractFFTs.rfft),
]
    # TODO: should we try to alter the image headers to change the units?
    @eval ($f)(img::AstroImage, args...; kwargs...) = copyheader(img, $f(parent(img), args...; kwargs...))
end

for f in [
    :(AbstractFFTs.fftshift),
]
    # TODO: should we try to alter the image headers to change the units?
    @eval ($f)(img::AstroImage, args...; kwargs...) = shareheader(img, $f(parent(img), args...; kwargs...))
end

# AbstractFFTs.complexfloat(img::AstroImage{T,N,D,R,StridedArray{T}}) where {T<:Complex{<:BlasReal}} = img
# AbstractFFTs.realfloat(img::AstroImage{T,N,D,R,StridedArray{T}}) where {T<:BlasReal} = img