text:
export getdoughertylist4copyvars3steps

const doughertylist4copyvars3steps = [
    2 4 2 1 3 1 0 2 0;
    2 3 3 1 5 2 0 2 0;
    3 6 3 1 6 3 0 3 0;
    2 4 2 1 2 0 0 2 3;
    2 3 3 2 2 0 0 2 0;
    4 6 4 3 4 2 1 4 0;
    2 5 2 1 2 0 0 2 0;
    2 4 3 1 2 0 0 2 0;
    2 4 1 2 2 3 0 2 0;
    3 7 4 1 4 1 0 3 0;
    4 6 11 3 6 2 0 4 0;
    3 6 3 1 4 1 0 3 5;
    7 8 12 12 7 5 5 7 0;
    5 14 9 1 7 2 0 5 0;
    6 7 11 11 6 3 3 6 0;
    3 4 6 3 6 2 0 3 0;
    11 23 28 3 11 7 5 11 0;
    5 6 8 7 5 3 2 5 0;
    6 12 6 3 6 4 3 6 0;
    4 5 16 4 10 6 0 4 0;
    3 6 5 1 5 3 0 3 0;
    4 13 7 1 4 2 2 4 0;
    4 5 7 6 4 1 1 4 0;
    4 8 4 1 10 6 0 4 0;
    5 16 13 1 5 1 1 5 0;
    5 6 11 11 5 1 1 5 0;
    2 3 4 1 4 5 0 2 0;
    4 5 6 4 4 2 4 4 0;
    4 7 4 2 4 1 1 4 7;
    2 3 3 1 2 3 2 2 0;
    4 10 2 3 9 9 0 4 0;
    4 7 3 4 5 5 0 4 0;
    3 5 4 2 3 0 0 3 4;
    5 14 11 1 5 2 2 5 0;
    6 15 10 2 6 0 0 6 11;
    11 31 18 3 13 4 0 11 0;
    18 38 31 6 18 6 6 18 0;
    4 9 3 2 8 4 0 4 0;
    5 12 3 3 10 9 0 5 0;
    8 19 6 4 9 14 0 8 0;
    3 5 4 2 4 1 0 3 0;
    7 19 8 1 9 8 2 7 0;
    6 16 2 9 6 11 0 6 0;
    7 8 11 11 7 7 7 7 0;
    5 8 10 3 5 1 1 5 0;
    4 10 10 1 4 1 4 4 0;
    8 9 14 14 8 8 8 8 0;
    3 5 4 1 8 5 0 3 0;
    6 11 10 2 6 6 9 6 0;
    7 19 11 2 7 5 10 7 0;
    6 13 10 2 6 2 2 6 0;
    9 12 16 7 9 13 11 9 0;
    7 8 16 16 7 3 3 7 0;
    5 9 4 4 5 3 1 5 0;
    8 17 7 3 16 8 0 8 0;
    3 9 2 2 3 0 0 3 0;
    9 17 15 3 9 15 5 9 0;
    3 4 5 4 4 1 0 3 0;
    7 16 12 2 7 4 3 7 0;
    3 5 6 2 3 0 0 3 0;
    6 11 10 2 6 9 4 6 0;
    6 7 13 13 6 2 2 6 0;
    10 23 16 3 10 5 5 10 0;
    4 5 9 6 6 3 0 4 0;
    9 19 8 4 15 7 0 9 0;
    10 11 22 22 10 9 9 10 0;
    5 20 15 1 5 0 0 5 0;
    4 8 11 1 4 4 4 4 0;
    7 16 12 2 7 3 4 7 0;
    4 9 8 1 4 3 2 4 0;
    6 12 5 5 6 3 0 6 0;
    7 12 10 3 14 6 0 7 0;
    13 24 16 5 26 10 0 13 0;
    6 14 4 3 12 11 0 6 0;
    10 18 11 7 13 3 0 10 0;
    12 22 14 5 23 9 0 12 0;
    4 14 9 1 4 0 0 4 0;
    6 11 10 2 6 8 5 6 0;
    5 11 4 2 10 8 0 5 0;
    8 14 13 4 12 8 0 8 0;
    10 20 9 6 12 7 0 10 0;
    8 17 9 5 8 0 0 8 0;
    3 8 1 4 6 7 0 3 0;
    8 9 15 15 8 6 6 8 0;
    7 8 21 21 7 2 2 7 0;
    5 6 15 15 5 0 0 5 0;
    5 6 5 5 5 10 10 5 0;
    9 10 19 19 9 7 7 9 0;
    6 7 17 17 6 1 1 6 0;
    10 17 11 11 10 0 0 10 0;
    4 9 3 2 5 7 0 4 0;
    4 6 5 4 4 4 0 4 0;
    4 5 10 9 4 0 0 4 0;
    3 4 5 5 3 0 0 3 0;
    6 13 5 3 8 0 2 6 10;
    3 4 4 5 4 1 0 3 0;
    7 25 18 1 13 6 0 7 0;
    10 17 8 9 10 4 2 10 0;
    9 16 14 4 15 9 0 9 0;
    9 10 18 18 9 8 8 9 0;
    7 15 6 3 12 7 0 7 0;
    3 6 5 1 3 1 3 3 0;
    5 9 6 2 10 5 0 5 0;
    7 35 28 1 11 4 0 7 0;
    4 7 5 1 13 8 0 4 0;
    7 19 2 11 7 12 1 7 0;
    4 13 9 1 4 0 0 4 7;
    4 13 10 1 4 0 0 4 0;
    3 8 5 1 3 0 0 3 0;
    3 6 2 2 4 0 1 3 5;
    4 9 8 1 4 2 3 4 0;
    5 7 11 6 5 2 0 5 0;
    5 10 5 1 15 10 0 5 0;
    4 11 7 1 4 4 5 4 0;
    8 17 7 4 10 0 2 8 12;
    7 8 9 9 7 14 14 7 0;
    4 5 10 6 7 2 0 4 0;
    10 26 3 16 10 18 1 10 0;
    9 12 22 10 13 10 0 9 0;
    3 7 2 2 3 0 0 3 2;
    7 9 16 10 10 2 0 7 0;
    6 9 6 4 6 3 4 6 0;
    4 8 2 5 4 0 0 4 4;
    13 37 22 3 17 6 0 13 0;
    8 10 23 12 8 1 4 8 0;
    6 7 7 7 6 11 11 6 0;
    11 29 3 19 11 21 2 11 0;
    9 12 8 13 9 9 2 9 0;
    4 7 5 2 8 2 0 4 0;
    4 11 7 1 4 1 7 4 0;
    8 14 6 7 8 4 2 8 0;
    6 8 6 5 6 4 5 6 0;
    7 13 5 6 7 5 1 7 0;
    12 16 25 18 12 3 0 12 0;
    8 12 9 5 8 4 8 8 0;
    5 6 14 9 5 1 2 5 0;
    15 27 45 5 15 27 11 15 0;
    10 18 12 4 21 7 0 10 0;
    9 17 9 6 11 4 0 9 0;
    11 29 16 3 15 6 0 11 0;
    18 22 51 30 18 3 6 18 0;
    3 6 6 1 3 1 1 3 0;
    18 32 13 16 18 13 2 18 0;
    6 14 9 2 6 2 3 6 0;
    14 24 13 14 14 11 0 14 0;
    5 9 5 4 6 2 0 5 0;
    6 17 2 9 9 3 0 6 0;
    4 11 7 1 4 2 6 4 0;
    18 30 20 12 37 15 0 18 0;
    8 10 22 13 8 1 2 8 0;
    5 8 3 5 5 5 2 5 0;
    9 19 7 5 12 8 0 9 0;
    4 5 4 4 4 6 5 4 0;
    14 28 10 13 14 0 0 14 11;
    7 15 6 4 8 6 0 7 0;
    7 18 26 2 7 2 1 7 0;
    11 19 34 4 11 20 8 11 0;
    4 5 5 4 4 3 3 4 0;
    5 13 8 1 8 3 0 5 0;
    6 7 15 12 6 2 3 6 0;
    14 21 17 9 14 6 12 14 0;
    3 5 2 4 3 2 0 3 0;
    4 7 4 3 6 2 0 4 0;
    19 32 21 12 40 17 0 19 0;
    10 20 9 7 11 5 0 10 0;
    8 21 34 2 8 3 2 8 0;
    6 8 6 7 6 2 1 6 0;
    4 10 2 3 5 1 0 4 4;
    4 5 9 7 5 1 0 4 0;
    5 8 5 3 5 4 4 5 0;
    4 13 1 9 4 7 0 4 0;
    9 25 12 3 11 4 0 9 0;
    5 6 7 6 5 4 4 5 0;
    6 9 9 5 8 2 0 6 0;
    7 11 6 6 7 3 2 7 0;
    4 9 5 1 7 3 0 4 0;
    9 18 7 6 12 9 0 9 0;
    3 7 2 2 3 4 0 3 0;
    7 14 6 4 10 7 0 7 0;
    6 8 8 5 6 2 1 6 0;
    9 24 31 2 9 12 3 9 0;
    9 24 43 2 9 4 3 9 0;
    5 9 7 2 5 3 2 5 0;
    13 20 8 14 13 12 6 13 0;
    5 7 7 4 5 1 2 5 0;
    21 32 24 13 21 10 18 21 0;
    5 10 3 6 5 0 0 5 2;
    6 11 6 3 12 6 0 6 0;
    20 34 17 18 20 5 4 20 0;
    24 40 17 24 24 13 6 24 0;
    5 11 14 1 5 6 6 5 0;
    7 18 2 21 7 1 7 7 3;
    8 10 14 13 8 2 1 8 0;
    5 8 6 5 5 2 0 5 0;
    19 29 22 12 19 8 16 19 0;
    4 6 5 4 6 2 0 4 0;
    26 43 25 21 26 6 7 26 0;
    11 19 28 4 11 20 12 11 0;
    5 10 19 1 5 8 8 5 0;
    7 13 6 7 7 3 0 7 0;
    36 61 35 29 36 6 11 36 0;
    6 9 8 4 6 2 5 6 0;
    10 12 18 17 10 4 3 10 0;
    14 36 45 4 14 17 2 14 0;
    14 22 10 5 16 20 7 14 0;
    3 6 5 1 6 2 0 3 0;
    11 19 30 4 11 20 10 11 0;
    20 34 20 16 20 3 6 20 0;
    16 42 51 4 16 19 4 16 0;
    15 25 15 12 15 3 4 15 0;
    8 21 2 25 8 2 9 8 1;
    47 79 37 41 50 15 9 47 0;
    25 66 81 6 25 30 7 25 0;
    20 30 14 15 24 8 7 20 0]

const doughertycoefs = [
    (1,2,0),
    (1,2,3),
    (1,3,2),
    (2,3,1),
    (1,2,4),
    (1,4,2),
    (2,4,1),
    (3,4,0),
    (3,4,1)]

function mysubmod(jki::Tuple{Int,Int,Int}, σ)
    J = mymap(σ, set(jki[1]), 4)
    K = mymap(σ, set(jki[2]), 4)
    I = mymap(σ, jki[3] == 0 ? emptyset() : set(jki[3]), 4)
    submodular(4, J, K, I)
end

function doughertyrow2dualentropy(x, σ=collect(1:4))
    h = -mysubmod(doughertycoefs[1], σ) * x[1]
    for i in 2:9
        h += mysubmod(doughertycoefs[i], σ) * x[i]
    end
    h
end

function getdoughertylist4copyvars3steps(perm=true)
    m = size(doughertylist4copyvars3steps, 1)
    if perm
        list = Vector{DualEntropy{15,Int}}(24*m)
        for i in 1:m
            coefs = doughertylist4copyvars3steps[i,:]
            σ = collect(1:4)
            for j in 1:24
                list[(i-1)*24+j] = doughertyrow2dualentropy(coefs, σ)
                nextp = next_perm(σ)
                @assert (j == 24) $ nextp
            end
        end
        sort!(list)
        keep = BitSet()
        for i in 1:length(list)
            if i == 1 || list[i] != list[i-1]
                push!(keep, i)
            end
        end
        list[collect(keep)]
    else
        DualEntropy{15,Int}[doughertyrow2dualentropy(doughertylist4copyvars3steps[i,:]) for i in 1:m]
    end
end


text:
# Unit test functions
function unit_test_stack()
    s = Stack(Vector{Float64}(undef, 0))
    @test length(s.data) == 0

    return nothing
end

function unit_test_stack_isempty()
    s = Stack(Vector{Float64}(undef, 0))
    @test isempty(s)

    return nothing
end

function unit_test_stack_push()
    s = Stack(Vector{Float64}(undef, 0))
    val = 1.0
    push!(s, val)
    @test s.data[1] == val

    return nothing
end

function unit_test_stack_pop()
    s = Stack(Vector{Float64}(undef, 0))
    val = 1.0
    push!(s, val)
    @test pop!(s) == val
    @test isempty(s)

    return nothing
end

function unit_test_stack_peek()
    s = Stack(Vector{Float64}(undef, 0))
    val = 1.0
    push!(s, val)
    @test peek(s) == val

    return nothing
end

function unit_test_stack_size()
    s = Stack(Vector{Float64}(undef, 0))
    val = 1.0
    push!(s, val)
    @test size(s) == 1
    push!(s, val)
    @test length(s) == 2

    return nothing
end

function unit_test_stack_length()
    s = Stack(Vector{Float64}(undef, 0))
    val = 1.0
    push!(s, val)
    @test length(s) == 1
    push!(s, val)
    @test length(s) == 2

    return nothing
end

# Run the unit tests
@testset "Stack" begin
    unit_test_stack()
    unit_test_stack_isempty()
    unit_test_stack_push()
    unit_test_stack_pop()
    unit_test_stack_peek()
    unit_test_stack_size()
    unit_test_stack_length()
end

ext: jl


text:
using CPUTime
using Plots
using DifferentialEquations
using LinearAlgebra
using MatrixEquations

# function ARE(A, B, Q, R)
#     """有本ポッターで解く（未完成）"""
#     ℋ = [
#         A' -B*inv(R)*B'
#         -Q -A
#     ]  #  ハミルトン行列
#     println("ℋ = ", ℋ)

#     λ_ = eigvals(ℋ)  # ハミルトン行列の固有値
#     ω_ = eigvecs(ℋ)  # ハミルトン行列の固有ベクトル
#     println("ℋ の固有値 = ", λ_)
#     println("ℋ の固有ベクトル", ω_)

#     # ハミルトン行列の固有値が負のものを探す
#     index = []
#     for i in 1:size(λ_)[1]
#         print(i)
#         if real(λ_[i]) < 0
#             println("実部が負")
#             push!(index, i)
#         else
#             println("実部が正")
#         end
#     end
#     println("ℋ の固有値の実部が負なのは, ", index)
#     n = size(index)[1]

#     # Y,Zを計算
#     ω = Matrix{ComplexF64}(undef, size(ω_)[1], n)  # 未初期化のMatrix

#     for i in 1:n
#         global ω[:, i] .= ω_[:, index[i]]
#     end

#     Y = ω[1:n, :]
#     Z = ω[n+1:end, :]

#     P = Z * inv(Y)
#     println("P = ", P)
#     return P
# end

function AREbyMatrixEquations(A, B, Q, R)
    """ライブラリ使用
    https://juliapackages.com/p/matrixequations
    """
    S = zero(B)
    P, E, F = arec(A, B, R, Q, S)
    println("P = ", P)
    return P
end

function simu(A, B, C, D, Q, R, x0)
    """解を描写"""
    #P = ARE(A, B, Q, R)  # リカッチ方程式を解く
    P = AREbyMatrixEquations(A, B, Q, R)
    K = inv(R) * B' * P  # 最適フィードバックゲイン
    Ā = A - B * K
    println("Ā = ", Ā)
    println(eigvals(Ā))

    # 解
    x(t) = exp(Ā .* t) * x0  # 状態方程式の解
    u(t) = -K * x(t)  # 最適入力
    y(t) = C * x(t) + D * u(t)  # 出力方程式の解

    # Plot．システムの次元に合わせてplotの数を変更．
    plot(xlims=(0, 10), xaxis="time", )
    plot!(t -> x.(t)[1], label="x1")
    plot!(t -> x.(t)[2], label="x2")
    plot!(t -> x.(t)[3], label="x3")
end

### システムの例1 ###
A = [
    1.1 2.0 3.0
    0 0.95 1.20
    1.2 0.01 10.5
]
B = [
    1.0
    0.0
    0.847
]
C = Matrix(I, 2, 2)
D = [1 1 1]
# 重み行列
Q = diagm(0 => [1,1,1]) .* 1000
R = ones((1,1))
x0 = [1, 2, 3]

# ### システムの例2 ###
# A = [0.0 1.0; -10.0 -1.0]
# B = [0.0, 1.0]
# Q = [300 0; 0 60]
# R = ones((1, 1))
# C = ones((2, 2))
# D = B
# x0 = [1, 2]

@time simu(A, B, C, D, Q, R, x0)  # 実行


text:
struct SafetyMask{M <: MDP, P <: Policy}
    mdp::M
    policy::P
    threshold::Float64
end

function safe_actions(mask::M, s) where M <: SafetyMask
    A = actiontype(mask.mdp)
    safe_acts = A[]
    vals = actionvalues(mask.policy, s)
    safe = maximum(vals) > mask.threshold ? true : false
    if !safe # follow safe controller
        push!(safe_acts, actions(mask.mdp)[argmax(vals)])
    else
        for (j, a) in enumerate(actions(mask.mdp))
            if vals[j] > mask.threshold
                push!(safe_acts, a)
            end
        end
    end
    return safe_acts
end

# """
# extract a vector P of size |S| where P(s) is the probability of satisfying a property
# when starting in state s
# `get_proba(mdp::MDP, result::ModelCheckingResult)`
# """
# function get_proba(mdp::MDP, result::ModelCheckingResult)
#     P = zeros(n_states(mdp))
#     for (i, val) in enumerate(result.result[:get_values]())
#         P[i] = val
#     end
#     return P
# end

# """
# Returns a matrix of dimension |S|x|A| where each element is the probability of satisfying an LTL formula for a given state action pair.
# This algorithm is mathematically sound only for basic LTL property like "!a U b" !!! (should technically be over the product MDP)
# Arguments:
# - `mdp::MDP` the MDP model
# - `result::ModelCheckingResult` result from model checking
#  `get_state_action_proba(mdp::MDP, P::Vector{Float64})`
# """
# function get_state_action(mdp::MDP, result::ModelCheckingResult)
#     P = get_proba(mdp, result)
#     P_sa = get_state_action_proba(mdp, P)
#     return P_sa
# end
# function get_state_action_proba(mdp::MDP, P::Vector{Float64})
#     P_sa = zeros(n_states(mdp), n_actions(mdp))
#     states = ordered_states(mdp)
#     actions = ordered_actions(mdp)
#     for (si, s) in enumerate(states)
#         P[si] == 0. ? continue : nothing
#         for (ai, a) in enumerate(actions)
#             dist = transition(mdp, s, a)
#             for (sp, p) in  weighted_iterator(dist)
#                 p == 0.0 ? continue : nothing # skip if zero prob
#                 spi = stateindex(mdp, sp)
#                 P_sa[si, ai] += p * P[spi]
#             end
#         end
#     end
#     return P_sa
# end

# """
# Extract a scheduler as a `VectorPolicy` from the model checking result
# """
# function POMDPToolbox.VectorPolicy{S, A}(mdp::MDP{S, A}, result::ModelCheckingResult{S})
#     P = get_proba(mdp, result)
#     P_sa = get_state_action_proba(mdp, P)
#     Pmax, amax = fargmax(P_sa, 2)
#     amax = ind2sub.((size(P_sa),), amax)
#     actions = ordered_actions(mdp)
#     act_vec = Vector{A}(n_states(mdp))
#     map!(x -> actions[x[2]], act_vec, amax)
#     return VectorPolicy(mdp, act_vec)
# end

# """
# Build a safety mask: a `SafetyMask` object contains infromation on the "safety" of each state action pair of the MDP. Safety is measure as the probability of satisfying the desired LTL formula. It takes as input the mdp model, the result from model checking and a threshold on the probability of success. State action pairs for which the probability of success is below the threshold are considered as unsafe.
# `SafetyMask{S, A}(mdp::MDP{S, A}, result::ModelCheckingResult, threshold::Float64)`
# """
# struct SafetyMask{M <: MDP, A}
#     mdp::M
#     threshold::Float64
#     risk_vec::Vector{Float64}
#     risk_mat::Array{Float64, 2}
#     actions::Vector{A}
# end

# function SafetyMask{S, A}(mdp::MDP{S, A}, result::ModelCheckingResult, threshold::Float64)
#     P = get_proba(mdp, result)
#     P_sa = get_state_action_proba(mdp, P)
#     return SafetyMask(mdp, threshold, P, P_sa, actions(mdp))
# end

# """
# Returns a vector of safe actions to execute in state s
# An action is safe if the probability of success is above the threshold of the safety mask.
# `safe_actions{M, A, S}(mask::SafetyMask{M,A}, s::S)`
# """
# function safe_actions{M, A, S}(mask::SafetyMask{M,A}, s::S)
#     safe_acts = A[]
#     sizehint!(safe_acts, n_actions(mask.mdp))
#     si = stateindex(mask.mdp, s)
#     safe = mask.risk_vec[si] > mask.threshold ? true : false
#     if !safe # follow safe controller
#         push!(safe_acts, mask.actions[argmax(mask.risk_mat[si, :])])
#     else
#         for (j, a) in enumerate(mask.actions)
#             if mask.risk_mat[si, j] > mask.threshold
#                 push!(safe_acts, a)
#             end
#         end
#     end
#     return safe_acts
# end

# """
# returns the initial probability of satisfying \$\phi\$:
# \$\Sigma_s b(s)P(s \models \phi)\$

#     initial_probability(mdp::MDP, result::ModelCheckingResult)
# """
# function initial_probability(mdp::MDP, result::ModelCheckingResult)
#     P = get_proba(mdp, result)
#     p_init = 0.
#     d0 = initialstate_distribution(mdp)
#     for (s, p) in weighted_iterator(d0)
#         si = stateindex(mdp, s)
#         p_init += p*P[si]
#     end
#     return p_init
# end


text:
using Potato
using Test
using JuMP
using GLPK

@testset "Potato" begin
    @test 1 == 1
end


text:
module CartesianGrids

import Base: @propagate_inbounds, show, summary, fill!

#using Compat
using FFTW
using SpecialFunctions
using Statistics

using LinearAlgebra
using SparseArrays
using Interpolations

import LinearAlgebra: mul!, ldiv!, cross, ×, dot, ⋅

import Base: parentindices
const GAMMA = MathConstants.γ

export Primal, Dual, ScalarGridData, VectorGridData, TensorGridData, GridData,
       CollectedGridData,
       PointData, ScalarData, VectorData, TensorData,
       celltype, griddatatype, indexshift,
       diff!,grid_interpolate!,
       curl, curl!, Curl, divergence, divergence!, Divergence,
       grad, grad!, Grad,
       laplacian, laplacian!, laplacian_symm!, plan_laplacian, plan_laplacian!,
       helmholtz, helmholtz!, plan_helmholtz, plan_helmholtz!,
       plan_intfact,plan_intfact!,Identity,
       product, product!, ∘,
       magsq!,magsq,mag!,mag,
       directional_derivative!, directional_derivative_conserve!, curl_cross!,
       convective_derivative!, convective_derivative_rot!,
       DDF, GradDDF,
       Regularize, RegularizationMatrix, InterpolationMatrix,
       CircularConvolution,
       AbstractSpatialField, Gaussian, DGaussian, radius, center, strength,
       EmptySpatialField,
       SpatialGaussian,GeneratedField, datatype, grid, PulseField, ModulatedField

abstract type CellType end
abstract type Primal <: CellType end
abstract type Dual <: CellType end

abstract type GridData{NX,NY,T} <: AbstractMatrix{T} end
abstract type ScalarGridData{NX,NY,T} <: GridData{NX,NY,T} end
abstract type VectorGridData{NX,NY,T} <: GridData{NX,NY,T} end
abstract type TensorGridData{NX,NY,T} <: GridData{NX,NY,T} end
CollectedGridData = Union{VectorGridData,TensorGridData}

abstract type PointData{N,T} <: AbstractVector{T} end

# List of scalar grid types. Each pair of numbers specifies
# the number of grid points in each direction for this data type, relative
# to the reference grid. The two pairs of numbers correspond to Primal
# and Dual versions of this grid data type.
const SCALARLIST = [ :Nodes, (-1,-1), (0,0)],
                   [ :XEdges, (0,-1), (-1,0)],
                   [ :YEdges, (-1,0), (0,-1)]

# List of collection grid types. This specifies the list of different types
# to create regularization and coordinates routines for.
const VECTORLIST = [:Edges, [:Primal]],
                   [:Edges, [:Dual]],
                   [:NodePair, [:Primal,:Dual]],
                   [:NodePair, [:Dual,:Primal]]

const TENSORLIST = [:EdgeGradient, [:Dual,:Primal]],
                   [:EdgeGradient, [:Primal,:Dual]]

function othertype end

macro othertype(celltype, k)
    esc(quote
        othertype(::$celltype) = $k
        othertype(::Type{$celltype}) = $k
    end)
end

@othertype Primal Dual
@othertype Dual Primal
@othertype CellType CellType

unpack(bc::Base.Broadcast.Broadcasted) = unpack(bc.args)
unpack(args::Tuple) = unpack(unpack(args[1]), Base.tail(args))
unpack(x) = x
unpack(::Tuple{}) = nothing
unpack(a::GridData, rest) = a
unpack(a::PointData, rest) = a
unpack(::Any, rest) = unpack(rest)

@inline unpack_data(bc::Broadcast.Broadcasted, i) = Broadcast.Broadcasted(bc.f, unpack_data_args(i, bc.args))
unpack_data(x,::Any) = x
unpack_data(x::GridData, ::Nothing) = x.data
unpack_data(x::PointData, ::Nothing) = x.data

@inline unpack_data_args(i, args::Tuple) = (unpack_data(args[1], i), unpack_data_args(i, Base.tail(args))...)
unpack_data_args(i, args::Tuple{Any}) = (unpack_data(args[1], i),)
unpack_data_args(::Any, args::Tuple{}) = ()

#@wraparray ScalarGridData data 2

include("fields/fieldmacros.jl")
include("fields/scalargrid.jl")

# Generate the scalar grid field types and associated functions
for (wrapper,primaldn,dualdn) in SCALARLIST
    @eval @scalarfield $wrapper $primaldn $dualdn
end

include("fields/collections.jl")

include("gridoperations/basicoperations.jl")
include("gridoperations/innerproducts.jl")
include("gridoperations/convolution.jl")
include("gridoperations/lgf.jl")
include("gridoperations/lgf-helmholtz.jl")
include("gridoperations/laplacian.jl")
include("gridoperations/helmholtz.jl")
include("gridoperations/intfact.jl")
include("gridoperations/diffcalculus.jl")
include("gridoperations/differencing1d.jl")
include("gridoperations/interpolation1d.jl")
include("gridoperations/shift.jl")
include("gridoperations/nlcalculus.jl")

include("points/points.jl")
include("points/ddf.jl")
include("points/regularization.jl")

include("physical/physicalgrid.jl")
include("physical/interpolation.jl")
include("physical/generatefields.jl")

#== Plot Recipes ==#

include("plot_recipes.jl")

end

max_stars_repo_name: JuliaIBPM/CartesianGrids

max_issues_repo_name: JuliaIBPM/CartesianGrids

max_forks_repo_name: JuliaIBPM/CartesianGrids


text:
"""
    pack(x)

Serialize `x` to MessagePack format and return the resulting `Vector{UInt8}`.

This function uses [`msgpack_type`](@ref) and [`to_msgpack`](@ref) to determine
the appropriate translation of the `value` into MessagePack format.

See also: [`unpack`](@ref)
"""
function pack(x)
    io = IOBuffer(UInt8[]; append = true)
    pack(io, x)
    return take!(io)
end

"""
    pack(io::IO, x)

Like `pack(x)`, but write the resulting bytes to `io`.

See also: [`unpack`](@ref)
"""
function pack(io::IO, x)
    pack_type(io, msgpack_type(typeof(x)), x)
    return io
end

#####
##### `AnyType`
#####

# This function might've been reached via another `pack_type` method that
# relied on reflection and ended up calling `msgpack_type` on a type that
# didn't have a non-`AnyType` mapping, even if the underlying value it's trying
# serialize does (e.g. an object field with a `Union` type). Thus, before giving
# up, we first attempt to resolve the issue by calling `msgpack_type(typeof(x))`
# directly.
function pack_type(io, t::AnyType, x)
    tx = msgpack_type(typeof(x))
    if tx isa AnyType
        error("no non-`AnyType` MsgPack mapping found for ", typeof(x), "; please ",
              "overload `msgpack_type` for this type.")
    end
    return pack_type(io, tx, x)
end

#####
##### `ImmutableStructType` + `MutableStructType`
#####

function pack_type(io,
                   t::Union{ImmutableStructType,MutableStructType},
                   x::T) where {T}
    N = fieldcount(T)
    if N <= 15
        write(io, magic_byte_min(MapFixFormat) | UInt8(N))
    elseif N <= typemax(UInt16)
        write(io, magic_byte(Map16Format))
        write(io, hton(UInt16(N)))
    elseif N <= typemax(UInt32)
        write(io, magic_byte(Map32Format))
        write(io, hton(UInt32(N)))
    else
        invalid_pack(io, t, x)
    end
    Base.@nexprs 32 i -> begin
        F_i = fieldtype(T, i)
        pack_type(io, StringType(), fieldname(T, i))
        pack_type(io, msgpack_type(fieldtype(T, i)), getfield(x, i))
        N == i && return nothing
    end
    for i in 33:N
        F_i = fieldtype(T, i)
        pack_type(io, StringType(), fieldname(T, i))
        pack_type(io, msgpack_type(fieldtype(T, i)), getfield(x, i))
    end
    return nothing
end

#####
##### `IntegerType`
#####

function pack_type(io, t::IntegerType, x)
    x = to_msgpack(t, x)
    if x < 0
        x >= -32 && return pack_format(io, IntFixNegativeFormat(Int8(x)))
        x >= typemin(Int8) && return pack_format(io, Int8Format(), x)
        x >= typemin(Int16) && return pack_format(io, Int16Format(), x)
        x >= typemin(Int32) && return pack_format(io, Int32Format(), x)
        x >= typemin(Int64) && return pack_format(io, Int64Format(), x)
    else
        x <= 127 && return pack_format(io, IntFixPositiveFormat(UInt8(x)))
        x <= typemax(UInt8) && return pack_format(io, UInt8Format(), x)
        x <= typemax(UInt16) && return pack_format(io, UInt16Format(), x)
        x <= typemax(UInt32) && return pack_format(io, UInt32Format(), x)
        x <= typemax(UInt64) && return pack_format(io, UInt64Format(), x)
    end
    invalid_pack(io, t, x)
end

pack_format(io, f::Union{IntFixNegativeFormat,IntFixPositiveFormat}) = write(io, f.byte)
pack_format(io, f::Int8Format, x) = _pack_integer(io, f, Int8, x)
pack_format(io, f::Int16Format, x) = _pack_integer(io, f, Int16, x)
pack_format(io, f::Int32Format, x) = _pack_integer(io, f, Int32, x)
pack_format(io, f::Int64Format, x) = _pack_integer(io, f, Int64, x)
pack_format(io, f::UInt8Format, x) = _pack_integer(io, f, UInt8, x)
pack_format(io, f::UInt16Format, x) = _pack_integer(io, f, UInt16, x)
pack_format(io, f::UInt32Format, x) = _pack_integer(io, f, UInt32, x)
pack_format(io, f::UInt64Format, x) = _pack_integer(io, f, UInt64, x)

function _pack_integer(io, ::F, ::Type{T}, x) where {F,T}
    y = hton(T(x))
    write(io, magic_byte(F))
    write(io, y)
end

#####
##### `NilType`
#####

pack_type(io, ::NilType, x) = pack_format(io, NilFormat(), x)

pack_format(io, ::NilFormat, ::Any) = write(io, magic_byte(NilFormat))

#####
##### `BooleanType`
#####

function pack_type(io, t::BooleanType, x)
    x = to_msgpack(t, x)
    x == true && return pack_format(io, TrueFormat(), x)
    x == false && return pack_format(io, FalseFormat(), x)
    invalid_pack(io, t, x)
end

pack_format(io, ::TrueFormat, ::Any) = write(io, magic_byte(TrueFormat))
pack_format(io, ::FalseFormat, ::Any) = write(io, magic_byte(FalseFormat))

#####
##### `FloatType`
#####

function pack_type(io, t::FloatType, x)
    x = to_msgpack(t, x)
    x isa Float32 && return pack_format(io, Float32Format(), x)
    x isa Float64 && return pack_format(io, Float64Format(), x)
    invalid_pack(io, t, x)
end

function pack_format(io, ::Float32Format, x)
    y = Float32(x)
    write(io, magic_byte(Float32Format))
    write(io, hton(y))
end

function pack_format(io, ::Float64Format, x)
    y = Float64(x)
    write(io, magic_byte(Float64Format))
    write(io, hton(y))
end

#####
##### `StringType`
#####

function pack_type(io, t::StringType, x)
    x = to_msgpack(t, x)
    n = x isa PointerString ? x.len : sizeof(x)
    n <= 31 && return pack_format(io, StrFixFormat(magic_byte_min(StrFixFormat) | UInt8(n)), x)
    n <= typemax(UInt8) && return pack_format(io, Str8Format(), x)
    n <= typemax(UInt16) && return pack_format(io, Str16Format(), x)
    n <= typemax(UInt32) && return pack_format(io, Str32Format(), x)
    invalid_pack(io, t, x)
end

function pack_format(io, f::StrFixFormat, x)
    write(io, f.byte)
    write(io, x)
end

function pack_format(io, ::Str8Format, x)
    write(io, magic_byte(Str8Format))
    write(io, UInt8(sizeof(x)))
    write(io, x)
end

function pack_format(io, ::Str16Format, x)
    write(io, magic_byte(Str16Format))
    write(io, hton(UInt16(sizeof(x))))
    write(io, x)
end

function pack_format(io, ::Str32Format, x)
    write(io, magic_byte(Str32Format))
    write(io, hton(UInt32(sizeof(x))))
    write(io, x)
end

#####
##### `BinaryType`
#####

function pack_type(io, t::BinaryType, x)
    x = to_msgpack(t, x)
    n = length(x)
    n <= typemax(UInt8) && return pack_format(io, Bin8Format(), x)
    n <= typemax(UInt16) && return pack_format(io, Bin16Format(), x)
    n <= typemax(UInt32) && return pack_format(io, Bin32Format(), x)
    invalid_pack(io, t, x)
end

function pack_format(io, ::Bin8Format, x)
    write(io, magic_byte(Bin8Format))
    write(io, UInt8(length(x)))
    write(io, x)
end

function pack_format(io, ::Bin16Format, x)
    write(io, magic_byte(Bin16Format))
    write(io, hton(UInt16(length(x))))
    write(io, x)
end

function pack_format(io, ::Bin32Format, x)
    write(io, magic_byte(Bin32Format))
    write(io, hton(UInt32(length(x))))
    write(io, x)
end

#####
##### `ArrayType`
#####

function pack_type(io, t::ArrayType, x)
    x = to_msgpack(t, x)
    n = length(x)
    n <= 15 && return pack_format(io, ArrayFixFormat(magic_byte_min(ArrayFixFormat) | UInt8(n)), x)
    n <= typemax(UInt16) && return pack_format(io, Array16Format(), x)
    n <= typemax(UInt32) && return pack_format(io, Array32Format(), x)
    invalid_pack(io, t, x)
end

function pack_format(io, f::ArrayFixFormat, x)
    write(io, f.byte)
    for i in x
        pack(io, i)
    end
end

function pack_format(io, ::Array16Format, x)
    write(io, magic_byte(Array16Format))
    write(io, hton(UInt16(length(x))))
    for i in x
        pack(io, i)
    end
end

function pack_format(io, ::Array32Format, x)
    write(io, magic_byte(Array32Format))
    write(io, hton(UInt32(length(x))))
    for i in x
        pack(io, i)
    end
end

#####
##### `MapType`
#####

function pack_type(io, t::MapType, x)
    x = to_msgpack(t, x)
    n = length(x)
    n <= 15 && return pack_format(io, MapFixFormat(magic_byte_min(MapFixFormat) | UInt8(n)), x)
    n <= typemax(UInt16) && return pack_format(io, Map16Format(), x)
    n <= typemax(UInt32) && return pack_format(io, Map32Format(), x)
    invalid_pack(io, t, x)
end

function pack_format(io, f::MapFixFormat, x)
    write(io, f.byte)
    for (k, v) in x
        pack(io, k)
        pack(io, v)
    end
end

function pack_format(io, ::Map16Format, x)
    write(io, magic_byte(Map16Format))
    write(io, hton(UInt16(length(x))))
    for (k, v) in x
        pack(io, k)
        pack(io, v)
    end
end

function pack_format(io, ::Map32Format, x)
    write(io, magic_byte(Map32Format))
    write(io, hton(UInt32(length(x))))
    for (k, v) in x
        pack(io, k)
        pack(io, v)
    end
end

#####
##### `ExtensionType`
#####

function pack_type(io, t::ExtensionType, x)
    x = to_msgpack(t, x)::Extension
    n = length(x.data)
    n == 1 && return pack_format(io, ExtFix1Format(), x)
    n == 2 && return pack_format(io, ExtFix2Format(), x)
    n == 4 && return pack_format(io, ExtFix4Format(), x)
    n == 8 && return pack_format(io, ExtFix8Format(), x)
    n == 16 && return pack_format(io, ExtFix16Format(), x)
    n <= typemax(UInt8) && return pack_format(io, Ext8Format(), x)
    n <= typemax(UInt16) && return pack_format(io, Ext16Format(), x)
    n <= typemax(UInt32) && return pack_format(io, Ext32Format(), x)
    invalid_pack(io, t, x)
end

const ExtFixFormat = Union{ExtFix1Format,ExtFix2Format,ExtFix4Format,ExtFix8Format,ExtFix16Format}

function pack_format(io, ::F, x::Extension) where {F<:ExtFixFormat}
    write(io, magic_byte(F))
    write(io, x.type)
    write(io, x.data)
end

function pack_format(io, ::Ext8Format, x::Extension)
    write(io, magic_byte(Ext8Format))
    write(io, UInt8(length(x.data)))
    write(io, x.type)
    write(io, x.data)
end

function pack_format(io, ::Ext16Format, x::Extension)
    write(io, magic_byte(Ext16Format))
    write(io, hton(UInt16(length(x.data))))
    write(io, x.type)
    write(io, x.data)
end

function pack_format(io, ::Ext32Format, x::Extension)
    write(io, magic_byte(Ext32Format))
    write(io, hton(UInt32(length(x.data))))
    write(io, x.type)
    write(io, x.data)
end

#####
##### utilities
#####

@noinline function invalid_pack(io, t, x)
    error("cannot serialize Julia value $(x) as MsgPack type $(t) to $(io)")
end


text:
import Trebuchet as Trebuchets

println(names(Trebuchets))

@doc Trebuchets.shoot

Trebuchets.shoot(5, 0.25pi, 500)

max_stars_repo_name: lucaferranti/julia-novice

max_issues_repo_name: lucaferranti/julia-novice

max_forks_repo_name: lucaferranti/julia-novice


text:
using TimerOutputs, Printf

using Oceananigans
using Oceananigans.TurbulenceClosures

include("benchmark_utils.jl")

#####
##### Benchmark setup and parameters
#####

const timer = TimerOutput()

Nt = 10  # Number of iterations to use for benchmarking time stepping.

# Run benchmark across these parameters.
            Ns = [(32, 32, 32), (128, 128, 128)]
   float_types = [Float64]       # Float types to benchmark.
         archs = [CPU()]         # Architectures to benchmark on.
@hascuda archs = [CPU(), GPU()]  # Benchmark GPU on systems with CUDA-enabled GPUs.
      closures = [ConstantIsotropicDiffusivity, ConstantAnisotropicDiffusivity, SmagorinskyLilly,
	              VerstappenAnisotropicMinimumDissipation]

#####
##### Run benchmarks
#####

for arch in archs, FT in float_types, N in Ns, Closure in closures
	grid = RegularCartesianGrid(FT, size=N, extent=(1, 1, 1))
    model = IncompressibleModel(architecture=arch, float_type=FT, grid=grid, closure=Closure(FT))

    time_step!(model, 1)  # precompile

    bn =  benchmark_name(N, string(Closure), arch, FT)
    @printf("Running benchmark: %s...\n", bn)
    for i in 1:Nt
        @timeit timer bn time_step!(model, 1)
    end
end

#####
##### Print benchmark results
#####

println()
print_benchmark_info()
print_timer(timer, title="Turbulence closure benchmarks", sortby=:name)
println()


text:
using HartreeFock
using Test

@testset "HartreeFock.jl" begin
    include("scf.jl")
    include("utils.jl")
end


text:
import SDtoolbox: ct, cell_size, CANTERA_MECH_FILETYPE

P₁ = 1e5
T₁ = 300
ϕ = 1

X₁ = "H2:$(42*ϕ), O2:21,N2:79"
mech = "gri30."*CANTERA_MECH_FILETYPE

λ = cell_size(T₁,P₁,X₁,mech)

println("P₁ = $(P₁*1e-5) bar, T₁ = $T₁ K, ϕ = $ϕ : λ = $(λ*1e3) mm")


text:
# ------ promote_op with dimensions and units ------

for op in (.+, .-, +, -)
    @eval function promote_op{S<:Dimensions,T<:Dimensions}(::typeof($op),
        ::Type{S}, ::Type{T})
        if S==T   # add or subtract same dimension, get same dimension
            x
        else
            throw(DimensionError())
        end
    end
    @eval function promote_op{S<:DimensionedUnits,T<:DimensionedUnits}(
        ::typeof($op), ::Type{S}, ::Type{T})
        if dimension(S())==dimension(T())
            promote_type(S,T)
        else
            throw(DimensionError())
        end
    end
end

function promote_op{S<:Unitlike,T<:Unitlike}(op, ::Type{S}, ::Type{T})
    typeof(op(S(), T()))
end

# ------ promote_op with quantities ------

# quantity, quantity
function promote_op{T1,D1,U1,T2,D2,U2}(op, x::Type{Quantity{T1,D1,U1}},
    y::Type{Quantity{T2,D2,U2}})
    # figuring out numeric type can be subtle if D1 == D2 but U1 != U2.
    # in particular, consider adding 1m + 1cm... the numtype is not Int.
    unittype = promote_op(op, U1, U2)
    numtype = if D1 == D2
        promote_type(T1, T2, typeof(convfact(U1(),U2())))
    else
        promote_type(T1, T2)
    end
    if unittype == typeof(NoUnits)
        numtype
    else
        dimtype = typeof(dimension(unittype()))
        Quantity{numtype, dimtype, unittype}
    end
end

# dim'd, quantity
promote_op{T2,D1,D2,U}(op, ::Type{DimensionedQuantity{D1}},
    ::Type{Quantity{T2,D2,U}}) = DimensionedQuantity{promote_op(op,D1,D2)}
promote_op{T2,D1,D2,U}(op, x::Type{Quantity{T2,D2,U}},
    y::Type{DimensionedQuantity{D1}}) = DimensionedQuantity{promote_op(op,D2,D1)}

# number, quantity
function promote_op{R<:Number,S,D,U}(op, ::Type{R}, ::Type{Quantity{S,D,U}})
    unittype = promote_op(op, typeof(NoUnits), U)
    numtype = if D == Dimensions{()}
        promote_type(R, S, typeof(convfact(NoUnits,U())))
    else
        promote_type(R, S)
    end
    if unittype == typeof(NoUnits)
        numtype
    else
        dimtype = typeof(dimension(unittype()))
        Quantity{numtype, dimtype, unittype}
    end
end

function promote_op{R<:Number,S,D,U}(op, x::Type{Quantity{S,D,U}}, y::Type{R})
    unittype = promote_op(op, U, typeof(NoUnits))
    numtype = if D == Dimensions{()}
        promote_type(R, S, typeof(convfact(U(),NoUnits)))
    else
        promote_type(R, S)
    end
    if unittype == typeof(NoUnits)
        numtype
    else
        dimtype = typeof(dimension(unittype()))
        Quantity{numtype, dimtype, unittype}
    end
end

# dim'd, dim'd
promote_op{D1,D2}(op, ::Type{DimensionedQuantity{D1}},
    ::Type{DimensionedQuantity{D2}}) = DimensionedQuantity{promote_op(op,D1,D2)}

# dim'd, number
promote_op{D}(op, ::Type{DimensionedQuantity{D}}, ::Type{Number}) = Number
promote_op{D}(op, ::Type{Number}, ::Type{DimensionedQuantity{D}}) = Number

# ------ promote_op with units ------

# units, quantity
function promote_op{R<:Units,S,D,U}(op, ::Type{Quantity{S,D,U}}, ::Type{R})
    numtype = S
    unittype = typeof(op(U(), R()))
    if unittype == typeof(NoUnits)
        numtype
    else
        dimtype = typeof(dimension(unittype()))
        Quantity{numtype, dimtype, unittype}
    end
end
function promote_op{R<:Units,S,D,U}(op, x::Type{R}, y::Type{Quantity{S,D,U}})
    numtype = S
    unittype = typeof(op(R(), U()))
    if unittype == typeof(NoUnits)
        numtype
    else
        dimtype = typeof(dimension(unittype()))
        Quantity{numtype, dimtype, unittype}
    end
end

# units, number
function promote_op{R<:Number,S<:Units}(op, ::Type{R}, ::Type{S})
    unittype = typeof(op(NoUnits, S()))
    if unittype == typeof(NoUnits)
        R
    else
        dimtype = typeof(dimension(unittype()))
        Quantity{R, dimtype, unittype}
    end
end
function promote_op{R<:Number,S<:Units}(op, ::Type{S}, ::Type{R})
    unittype = typeof(op(S(), NoUnits))
    if unittype == typeof(NoUnits)
        R
    else
        dimtype = typeof(dimension(unittype()))
        Quantity{R, dimtype, unittype}
    end
end
# ------ promote_rule ------

# quantity, quantity (different dims)
promote_rule{S1,S2,D1,D2,U1,U2}(::Type{Quantity{S1,D1,U1}},
    ::Type{Quantity{S2,D2,U2}}) = Number

# quantity, quantity (same dims)
function promote_rule{S1,S2,D,U1,U2}(::Type{Quantity{S1,D,U1}},
    ::Type{Quantity{S2,D,U2}})

    numtype = promote_type(S1,S2,typeof(convfact(U1(),U2())))
    if promote_type(U1,U2) != typeof(NoUnits)
        Quantity{numtype, D, promote_type(U1,U2)}
    else
        numtype
    end
end

# quantity, quantity (same dims, same units)
promote_rule{S1,S2,D,U}(::Type{Quantity{S1,D,U}}, ::Type{Quantity{S2,D,U}}) =
    Quantity{promote_type(S1,S2),D,U}

# dim'd, quantity (different dims)
promote_rule{S2,D1,D2,U}(::Type{DimensionedQuantity{D1}},
    ::Type{Quantity{S2,D2,U}}) = Number

# dim'd, quantity (same dims)
promote_rule{S2,D,U}(::Type{DimensionedQuantity{D}},
    ::Type{Quantity{S2,D,U}}) = DimensionedQuantity{D}

# number, quantity
function promote_rule{S,T<:Number,D,U}(::Type{Quantity{S,D,U}}, ::Type{T})
    if D == Dimensions{()}
        promote_type(S,T,typeof(convfact(U(),NoUnits)))
    else
        Number
    end
end

# dim'd, dim'd (different dims)
promote_rule{D1,D2}(::Type{DimensionedQuantity{D1}},
    ::Type{DimensionedQuantity{D2}}) = Number

# dim'd, dim'd (same dims)
promote_rule{D}(::Type{DimensionedQuantity{D}},
    ::Type{DimensionedQuantity{D}}) = DimensionedQuantity{D}

# dim'd, number
promote_rule{D,T<:Number}(::Type{DimensionedQuantity{D}}, ::Type{T}) = Number


text: # Write your solution here


text:
module ExploreExpression
export compute_symbolreferences, compute_usings, SymbolsState

import Base: union, ==

# from the source code: https://github.com/JuliaLang/julia/blob/master/src/julia-parser.scm#L9
const modifiers = [:(+=), :(-=), :(*=), :(/=), :(//=), :(^=), :(÷=), :(%=), :(<<=), :(>>=), :(>>>=), :(&=), :(⊻=), :(≔), :(⩴), :(≕)]

"SymbolsState trickles _down_ the ASTree: it carries referenced and defined variables from endpoints down to the root"
mutable struct SymbolsState
    references::Set{Symbol}
    assignments::Set{Symbol}
    funccalls::Set{Symbol}
    funcdefs::Dict{Symbol, SymbolsState}
end

SymbolsState(references, assignments, funccalls) = SymbolsState(references, assignments, funccalls, Dict{Symbol,SymbolsState}())
SymbolsState(references, assignments) = SymbolsState(references, assignments, Set{Symbol}())
SymbolsState() = SymbolsState(Set{Symbol}(), Set{Symbol}())

"ScopeState moves _up_ the ASTree: it carries scope information up towards the endpoints"
mutable struct ScopeState
    inglobalscope::Bool
    exposedglobals::Set{Symbol}
    hiddenglobals::Set{Symbol}
    definedfuncs::Set{Symbol}
end

function union(a::Dict{Symbol,SymbolsState}, b::Dict{Symbol,SymbolsState})
    # TODO: optimise: reuse `a` as `c`
    c = Dict{Symbol,SymbolsState}()
    for (k, v) in a
        c[k] = v
    end
    for (k, v) in b
        if haskey(c, k)
            c[k] = c[k] ∪ v
        else
            c[k] = v
        end
    end
    c
end

function union(a::SymbolsState, b::SymbolsState)
    SymbolsState(a.references ∪ b.references, a.assignments ∪ b.assignments, a.funccalls ∪ b.funccalls, a.funcdefs ∪ b.funcdefs)
end

function union(a::ScopeState, b::ScopeState)
    SymbolsState(a.inglobalscope && b.inglobalscope, a.exposedglobals ∪ b.exposedglobals, a.hiddenglobals ∪ b.hiddenglobals)
end

function ==(a::SymbolsState, b::SymbolsState)
    a.references == b.references && a.assignments == b.assignments&& a.funccalls == b.funccalls && a.funcdefs == b.funcdefs
end

function will_assign_global(assignee::Symbol, scopestate::ScopeState)::Bool
    (scopestate.inglobalscope || assignee ∈ scopestate.exposedglobals) && (assignee ∉ scopestate.hiddenglobals || assignee ∈ scopestate.definedfuncs)
end

function get_global_assignees(assignee_exprs, scopestate::ScopeState)::Set{Symbol}
    global_assignees = Set{Symbol}()
    for ae in assignee_exprs
        if isa(ae, Symbol)
            will_assign_global(ae, scopestate) && push!(global_assignees, ae)
        else
            if ae.head == :(::)
                will_assign_global(ae.args[1], scopestate) && push!(global_assignees, ae.args[1])
            else
                @warn "Unknown assignee expression" ae
            end
        end
    end
    global_assignees
end

# We handle a list of function arguments separately.
function extractfunctionarguments(funcdef::Expr)::Set{Symbol}
    argnames = Set{Symbol}()

    funcdefargs = if funcdef.head == :call
        funcdef.args[2:end]
    else
        funcdef.args
    end

    for a in funcdefargs
        if isa(a, Symbol)
            push!(argnames, a)
        elseif isa(a, Expr)
            if a.head == :(::)
                if length(a.args) == 2
                    push!(argnames, a.args[1])
                else
                    # f(..., ::JustAType, ...)
                end
            elseif a.head == :parameters || a.head == :tuple  # second is for ((a,b),(c,d)) -> a*b*c*d stuff
                push!(argnames, extractfunctionarguments(a)...)
            elseif a.head == :kw || a.head == :(=) # first is for unnamed function arguments, second is for lambdas
                push!(argnames, extractfunctionarguments(a.args[1])...)
            end
        end
    end

    return argnames
end

function extractfunctionarguments(funcdef::Symbol)::Set{Symbol}
    Set([funcdef])
end

# Possible leaf: value
# Like: a = 1
# 1 is a value (Int64)
function explore!(value, scopestate::ScopeState)::SymbolsState
    # includes: LineNumberNode, Int64, String,
    return SymbolsState(Set{Symbol}(), Set{Symbol}(), Set{Symbol}(), Dict{Symbol,SymbolsState}())
end

# Possible leaf: symbol
# Like a = x
# x is a symbol
# We handle the assignment separately, and explore!(:a, ...) will not be called.
# Therefore, this method only handles _references_, which are added to the symbolstate, depending on the scopestate.
function explore!(sym::Symbol, scopestate::ScopeState)::SymbolsState
    return if !(sym in scopestate.hiddenglobals)
        SymbolsState(Set([sym]), Set{Symbol}(), Set{Symbol}(), Dict{Symbol,SymbolsState}())
    else
        SymbolsState(Set{Symbol}(), Set{Symbol}(), Set{Symbol}(), Dict{Symbol,SymbolsState}())
    end
end

# General recursive method. Is never a leaf.
# Modifies the `scopestate`.
function explore!(ex::Expr, scopestate::ScopeState)::SymbolsState
    if ex.head == :(=)
        # Does not create scope

        assignees = if isa(ex.args[1], Symbol)
            # x = 123
            [ex.args[1]]
        elseif isa(ex.args[1], Expr)
            if ex.args[1].head == :tuple
                # (x, y) = (1, 23)
                filter(s->s isa Symbol, ex.args[1].args)
            elseif ex.args[1].head == :(::)
                # TODO: type is referenced
                [ex.args[1].args[1]]
            elseif ex.args[1].head == :ref
                # TODO: what is the desired behaviour here?
                # right now, it registers no reference, and no assignment
                []
            elseif ex.args[1].head == :(.)
                # TODO: what is the desired behaviour here?
                []
            elseif ex.args[1].head == :call
                # f(x, y) = x + y
                # Rewrite to:
                # function f(x, y) x + y end
                return explore!(Expr(:function, ex.args...), scopestate)
            else
                @warn "unknow use of =. Assignee is unrecognised." ex.args[1]
                []
            end
        else
            # When you assign to a datatype like Int, String, or anything bad like that
            # e.g. 1 = 2
            # This is parsable code, so we have to treat it
            []
        end
        val = ex.args[2]

        global_assignees = get_global_assignees(assignees, scopestate)

        # If we are _not_ assigning a global variable
        for assignee in setdiff(assignees, global_assignees)
            # Then this symbol hides any global definition with that name
            scopestate.hiddenglobals = union(scopestate.hiddenglobals, [assignee])
        end

        innersymstate = explore!(val, scopestate)

        symstate = innersymstate

        for assignee in global_assignees
            scopestate.hiddenglobals = union(scopestate.hiddenglobals, [assignee])
            symstate.assignments = union(symstate.assignments, [assignee])
        end

        return symstate
    elseif ex.head in modifiers
        # We change: a[1] += 123
        # to:        a[1] = a[1] + 123
        # We transform the modifier back to its operator
        # for when users redefine the + function

        operator = Symbol(string(ex.head)[1:end - 1])
        expanded_expr = Expr(:(=), ex.args[1], Expr(:call, operator, ex.args[1], ex.args[2]))
        return explore!(expanded_expr, scopestate)
    elseif ex.head == :let || ex.head == :for || ex.head == :while
        # Creates local scope

        # Because we are entering a new scope, we create a copy of the current scope state, and run it through the expressions.
        innerscopestate = deepcopy(scopestate)
        innerscopestate.inglobalscope = false

        return mapfoldl(a -> explore!(a, innerscopestate), ∪, ex.args, init=SymbolsState())
    elseif ex.head == :struct
        # Creates local scope

        structname = assignee = if isa(ex.args[2], Symbol)
            ex.args[2]
        else
            # We have:   struct a <: b
            ex.args[2].args[1]
            # TODO: record reactive reference to type
        end
        structfields = ex.args[3].args

        equiv_func = Expr(:function, Expr(:call, structname, structfields...), Expr(:block, nothing))

        return explore!(equiv_func, scopestate)
    elseif ex.head == :generator
        # Creates local scope

        # In a `generator`, a single expression is followed by the iterator assignments.
        # In a `for`, this expression comes at the end.

        # TODO: this is not the normal form of a `for`.
        return explore!(Expr(:for, ex.args[2:end]..., ex.args[1]), scopestate)
    elseif ex.head == :function
        symstate = SymbolsState()
        # Creates local scope

        funcroot = if ex.args[1].head == :(::)
            # TODO: record reactive reference to type
            ex.args[1].args[1]
        else
            ex.args[1]
        end

        # is either `f` (Symbol) or `f::String` (Expression)
        funcname_expr = assignee = funcroot.args[1]

        # is either [funcname] or []
        global_assignees = get_global_assignees([funcname_expr], scopestate)

        # Because we are entering a new scope, we create a copy of the current scope state, and run it through the expressions.

        innerscopestate = deepcopy(scopestate)
        innerscopestate.hiddenglobals = union(innerscopestate.hiddenglobals, extractfunctionarguments(funcroot))
        innerscopestate.inglobalscope = false

        innersymstate = explore!(Expr(:block, ex.args[2:end]...), innerscopestate)

        if !isempty(global_assignees)
            funcname = first(global_assignees)

            symstate.funcdefs[funcname] = innersymstate
            push!(scopestate.hiddenglobals, funcname)
            push!(scopestate.definedfuncs, funcname)
            push!(symstate.assignments, funcname)
        else
            # The function is not defined globally. However, the function can still modify the global scope or reference globals, e.g.

            # let
            #     function f(x)
            #         global z = x + a
            #     end
            #     f(2)
            # end

            # so we insert the function's inner symbol state here, as if it was a `let` block.
            symstate = innersymstate
        end

        return symstate
    elseif ex.head == :(->)
        # Creates local scope

        tempname = Symbol("anon",rand(UInt64))

        # We will rewrite this to a normal function definition, with a temporary name
        funcroot = ex.args[1]
        args_ex = if funcroot isa Symbol || (funcroot isa Expr && funcroot.head == :(::))
            [funcroot]
        elseif funcroot.head == :tuple
            funcroot.args
        else
            @error "Unknown lambda type"
        end

        equiv_func = Expr(:function, Expr(:call, tempname, args_ex...), ex.args[2])

        return explore!(equiv_func, scopestate)
    elseif ex.head == :global
        symstate = SymbolsState()
        # Does not create scope

        # We have one of:
        # global x;
        # global x = 1;
        # global x += 1;

        # where x can also be a tuple:
        # global a,b = 1,2

        globalisee = ex.args[1]

        if isa(globalisee, Symbol)
            scopestate.exposedglobals = union(scopestate.exposedglobals, [globalisee])
        elseif isa(globalisee, Expr)
            innerscopestate = deepcopy(scopestate)
            innerscopestate.inglobalscope = true
            innersymstate = explore!(globalisee, innerscopestate)
            symstate = innersymstate
        else
            @error "unknow global use"
        end

        return symstate
    elseif ex.head == :local
        symstate = SymbolsState()
        # Does not create scope

        # Logic similar to :global
        localisee = ex.args[1]

        if isa(localisee, Symbol)
            scopestate.hiddenglobals = union(scopestate.hiddenglobals, [localisee])
        elseif isa(localisee, Expr)
            innerscopestate = deepcopy(scopestate)
            innerscopestate.inglobalscope = false
            innersymstate = explore!(localisee, innerscopestate)
            symstate = innersymstate
        else
            @error "unknow local use"
        end

        return symstate
    elseif ex.head == :tuple
        # Does not create scope
        symstate = SymbolsState()

        # Is something like:
        # a,b,c = 1,2,3

        # This parses to:
        # head: Symbol tuple
        # args: a, b, :(c=1), 2, 3

        # 🤔
        # This one is very messy...

        recursers = ex.args

        # We have a tuple assignment if one of the tuple elements is an assignment expression:
        indexoffirstassignment = findfirst(a->isa(a, Expr) && a.head == :(=), ex.args)
        if indexoffirstassignment !== nothing
            recursers = ex.args[indexoffirstassignment:end]

            exposed = get_global_assignees(ex.args[1:indexoffirstassignment - 1], scopestate)

            scopestate.exposedglobals = union(scopestate.exposedglobals, exposed)
            symstate.assignments = exposed
        end

        return mapfoldl(a -> explore!(a, scopestate), ∪, recursers, init=symstate)
    elseif ex.head == :using || ex.head == :import
        if scopestate.inglobalscope
            imports = if ex.args[1].head == :(:)
                ex.args[1].args[2:end]
            else
                ex.args
            end

            packagenames = map(e->e.args[1], imports)

            return SymbolsState(Set{Symbol}(), Set{Symbol}(packagenames))
        else
            return SymbolsState(Set{Symbol}(), Set{Symbol}())
        end
    elseif ex.head == :macrocall && isa(ex.args[1], Symbol) && ex.args[1] == Symbol("@md_str")
        # Does not create scope
        # The Markdown macro treats things differently, so we must too

        parsed_markdown_str = Meta.parse("\"\"\"$(ex.args[3])\"\"\"", raise = false)
        innersymstate = explore!(parsed_markdown_str, scopestate)
        push!(innersymstate.references, Symbol("@md_str"))
        symstate = innersymstate
        return symstate
    elseif ex.head == :module
        # We completely ignore the contents

        return SymbolsState(Set{Symbol}(), Set{Symbol}([ex.args[2]]))
    elseif ex.head == :call && ex.args[1] isa Symbol
        # Does not create scope

        # We change the `call` to a `block` and recurse again (hitting the fallback below).
        # In particular, this adds the called function as a reference, which is what we want.
        symstate = explore!(Expr(:block, ex.args...), scopestate)
        push!(symstate.funccalls, ex.args[1])
        return symstate
    else
        # fallback, includes:
        # begin, block, do, call,
        # (and hopefully much more!)

        # Does not create scope (probably)

        return mapfoldl(a -> explore!(a, scopestate), ∪, ex.args, init=SymbolsState())
    end
end

function compute_symbolreferences(ex)::SymbolsState
    symstate = explore!(ex, ScopeState(true, Set{Symbol}(), Set{Symbol}(), Set{Symbol}()))

    # We do something special to account for recursive functions:
    # If a function `f` calls a function `g`, and both are defined inside this cell, the reference to `g` inside the symstate of `f` will be deleted.
    # the motivitation is that normally, an assignment (or function definition) will add that symbol to a list of 'hidden globals' - any future references to that symbol will be ignored. i.e. the _local definition hides a global_.
    # In the case of functions, you can reference functions and variables that do not yet exist, and so they won't be in the list of hidden symbols when the function definition is analysed.
    # Of course, our method will fail if a referenced function is defined both inside the cell **and** in another cell. However, this will lead to a MultipleDefinitionError before anything bad happens.
    for (func, inner_symstate) in symstate.funcdefs
        inner_symstate.references = setdiff(inner_symstate.references, keys(symstate.funcdefs))
        inner_symstate.funccalls = setdiff(inner_symstate.funccalls, keys(symstate.funcdefs))
    end
    symstate
end

# TODO: this can be done during the `explore` recursion
"Get the set of `using Module` expressions that are contained in this expression."
function compute_usings(ex)::Set{Expr}
    if isa(ex, Expr)
        if ex.head == :using
            Set{Expr}([ex])
        else
            union(Set{Expr}(), compute_usings.(ex.args)...)
        end
    else
        Set{Expr}()
    end
end

function is_pure_html(ex)::Bool
    if !(isa(ex, Expr) && ex.head == :macrocall && (ex.args[1] == Symbol("@md_str") || ex.args[1] == Symbol("@html_str")) && isa(ex.args[3], String))
        return false
    end
    if ex.args[1] == Symbol("@md_str")
        parsed_markdown_str = Meta.parse("\"\"\"$(ex.args[3])\"\"\"", raise = false)
        isa(parsed_markdown_str, String)
    else
        true
    end
end

end


text:
function mark_cross!(r::Robot)
    for side in (Nord,Ost,Sud,West)
        putmarkers!(r,side)
        if isborder(r,side)==true && isborder(r,next(side))==false
            obhod(r,side)
            putmarker!(r)
            putmarkers!(r,side)
        end
        side_inv=inverse(side)
        move_by_marker!(r,side_inv)
    end
putmarker!(r)
end
function obhod(r::Robot,side::HorizonSide)
    k=0
    j=0

    while isborder(r,side)==true && isborder(r,next(side))==false
        move!(r,next(side))
        k+=1
    end
    if isborder(r,side)==true && isborder(r,next(side))==true
        for _ in 1:k
            move!(r,inverse(next(side)))
        end
        k=0
    end
    while isborder(r,side)==true && isborder(r,inverse(next(side)))==false
        move!(r,inverse(next(side)))
        j+=1
    end
    if isborder(r,side)==true && isborder(r,inverse(next(side)))==true
        for _ in 1:j
            move!(r,next(side))
        end
    end
    if k==0
        k=j
        side1 = next(side)
    else
        side1 = inverse(next(side))
    end
    if isborder(r,side)==false
        print(side)
        move!(r,side)
        while (isborder(r,next(inverse(side)))==true &&(isborder(r,next(side))==false || isborder(r,inverse(next(side)))==false))
            move!(r,side)
        end
        for _ in 1:k
            move!(r,side1)
        end
    end
end
function putmarkers!(r::Robot,side::HorizonSide)
    while isborder(r,side)==false
        move!(r,side)
        putmarker!(r)
    end
end
function move_by_marker!(r::Robot,side::HorizonSide)
    while ismarker(r)==true
        if isborder(r,side)==false
            move!(r,side)
        elseif isborder(r,side)==true && isborder(r,next(side))==false
            obhod(r,side)
        end
    end
end
inverse(side::HorizonSide)=HorizonSide(mod(Int(side)+2,4))
next(side::HorizonSide)=HorizonSide(mod(Int(side)+3,4))


text:
using CompScienceMeshes

# G = readmesh("examples/sphere.in")
# edges = skeleton(G,1)
# faces = skeleton(G,2)
# connectMat = connectivity(edges,faces)
# NV, NE, NF = numvertices(G), numcells(edges), numcells(faces)
# nv = NV + NE
# verts = Array{vertextype(G)}(nv)
# G1 = bisecting_refinement(G)
# G2 = Loop_subdivision(G)
# G3 = Loop_subdivision(G2)
# G4 = Loop_subdivision(G3)
# writemesh(G2,"examples/sphere_subd1.in")
# writemesh(G3,"examples/sphere_subd2.in")
# writemesh(G4,"examples/sphere_subd3.in")
fn = joinpath(dirname(@__FILE__),"sphere_subd3.in")
#G = readmesh("/Users/Benjamin/Documents/sphere.in")
G = readmesh(fn)
# G = Loop_subdivision(G)
subdG=GSubdMesh(G)
nelem = length(subdG.elements)
area = 0.0
ng = 10
gpt,wt = gauss_points(ng)
for e = 1:nelem
    C = chart(subdG,e)
    for ig = 1 : ng*ng
        gpta = gpt[ig,1]
        gptb = gpt[ig,2]
        w = wt[ig]
        gpxi = (1+gpta)*(1-gptb)/4.0
        gpeta = (1+gptb)/2.0
        w *= (1-gpeta)/4.0
        g = neighborhood(C,[gpxi,gpeta])
        jac = jacobian(g)
        area += jac * w
    end
end
print(" area = ")
print(area)


text:
using Revise
using ZonalFlow
using JLD2

tspan   = (0.0,1000.0);
tsargs  = (
            dt=0.001,
            adaptive=false,
            progress=true,
            progress_steps=50000,
            save_everystep=false,
            save_start=true,
            # dense=true,
            save_noise=false,
            saveat=5
           );

dn      = "validate/ranktests/repeat/";
eqs     = [GQL(0),CE2()];

## Point Jet
domain  = Domain(extent=(2π,2π),res=(8,8));
coeffs  = Coefficients(Ω=2π,θ=0.0,μ=0.05,ν=0.0,ν₄=0.0);
forcing = PointJet(Ξ=1.0,Δθ=0.1,τ=20.0);
prob    = BetaPlane(domain,coeffs,forcing);

ql_pj      = integrate(prob,eqs[1],tspan;tsargs...)
write(prob,eqs[1],ql_pj,dn=dn,fn="ql_pj_500")

ce2_pj     = integrate(prob,eqs[2],tspan;tsargs...)
write(prob,eqs[2],ce2_pj,dn=dn,fn="ce2_pj_500")

ce2_pj     = integrate(prob,eqs[2],tspan;u0=ql_pj.u[end],tsargs...)
write(prob,eqs[2],ce2_pj,dn=dn,fn="ce2_pj_qlic_500")

## Kolmogorov forcing
domain  = Domain(extent=(2π,2π),res=(8,8));
coeffs  = Coefficients(Ω=0.0,θ=0.0,μ=0.0,ν=0.02,ν₄=0.0); # β=0
forcing = Kolmogorov(A₁=-1.0,A₄=-2.0);
prob    = BetaPlane(domain,coeffs,forcing);

ql_kf      = integrate(prob,eqs[1],tspan;tsargs...);
write(prob,eqs[1],ql_kf,dn=dn,fn="ql_kf_500");

ce2_kf     = integrate(prob,eqs[2],tspan;tsargs...);
write(prob,eqs[2],ce2_kf,dn=dn,fn="ce2_kf_500");

ce2_kf_qlic     = integrate(prob,eqs[2],tspan;u0=ql_kf.u[end],tsargs...);
write(prob,eqs[2],ce2_kf_qlic,dn=dn,fn="ce2_kf_qlic_500");

"""
# Kolmogorov with positivity check
eqs     = [GQL(0),CE2(poscheck=true,poscheckat=10)];

ce2_kf     = integrate(prob,eqs[2],tspan;tsargs...);
write(prob,eqs[2],ce2_kf,dn=dn,fn="ce2_kf_poscheck_500");

ce2_kf_qlic     = integrate(prob,eqs[2],tspan;u0=ql_kf.u[end],tsargs...);
write(prob,eqs[2],ce2_kf_qlic,dn=dn,fn="ce2_kf_poscheck_qlic_500");
"""

"""
# # Kolmogorov forcing lower viscosity
domain  = Domain(extent=(2π,2π),res=(11,11));
coeffs  = Coefficients(Ω=0.0,θ=0.0,μ=0.0,ν=0.01,ν₄=0.0); # β=0
forcing = Kolmogorov(A₁=-1.0,A₄=-2.0);
prob    = BetaPlane(domain,coeffs,forcing);

eqs     = [GQL(0),CE2(poscheck=true,poscheckat=5)];

ql_kf      = integrate(prob,eqs[1],tspan;tsargs...);
write(prob,eqs[1],ql_kf,dn=dn,fn="ql_kf_nu01_500");

ce2_kf     = integrate(prob,eqs[2],tspan;tsargs...);
write(prob,eqs[2],ce2_kf,dn=dn,fn="ce2_kf_nu01_500");

ce2_kf_qlic     = integrate(prob,eqs[2],tspan;u0=ql_kf.u[end],tsargs...);
write(prob,eqs[2],ce2_kf_qlic,dn=dn,fn="ce2_kf_nu01_qlic_500");
"""

# Kolmogorov Forcing: 11x11 grid
domain  = Domain(extent=(2π,2π),res=(9,9));
coeffs  = Coefficients(Ω=0.0,θ=0.0,μ=0.0,ν=0.02,ν₄=0.0); # β=0
forcing = Kolmogorov(A₁=-1.0,A₄=-2.0);
prob    = BetaPlane(domain,coeffs,forcing);

ql_kf_9      = integrate(prob,eqs[1],tspan;tsargs...);
write(prob,eqs[1],ql_kf_9,dn=dn,fn="ql_kf_9x9_500");

ce2_kf_9     = integrate(prob,eqs[2],tspan;tsargs...);
write(prob,eqs[2],ce2_kf_9,dn=dn,fn="ce2_kf_9x9_500");

ce2_kf_9_qlic     = integrate(prob,eqs[2],tspan;u0=ql_kf_9.u[end],tsargs...);
write(prob,eqs[2],ce2_kf_9_qlic,dn=dn,fn="ce2_kf_9x9_qlic_500");

# Kolmogorov Forcing: 11x11 grid
domain  = Domain(extent=(2π,2π),res=(11,11));
coeffs  = Coefficients(Ω=0.0,θ=0.0,μ=0.0,ν=0.02,ν₄=0.0); # β=0
forcing = Kolmogorov(A₁=-1.0,A₄=-2.0);
prob    = BetaPlane(domain,coeffs,forcing);

tspan   = (0.0,1000.0);

ql_kf_11      = integrate(prob,eqs[1],tspan;tsargs...);
write(prob,eqs[1],ql_kf_11,dn=dn,fn="ql_kf_11x11_1000");

tsargs  = (
            dt=0.0005,
            adaptive=false,
            progress=true,
            progress_steps=50000,
            save_everystep=false,
            save_start=true,
            # dense=true,
            save_noise=false,
            saveat=5
           );

ce2_kf_11     = integrate(prob,eqs[2],tspan;tsargs...);
write(prob,eqs[2],ce2_kf_11,dn=dn,fn="ce2_kf_11x11_500");
@save dn*"ce2_kf_11x11_500"*".jld2" ce2_kf_11

ce2_kf_11_qlic     = integrate(prob,eqs[2],tspan;u0=ql_kf_11.u[end],tsargs...);
write(prob,eqs[2],ce2_kf_11_qlic,dn=dn,fn="ce2_kf_11x11_qlic_500");
@save dn*"ce2_kf_qlic_11x11_500"*".jld2" ce2_kf_11_qlic

# Kolmogorov Forcing: 16x16 grid
domain  = Domain(extent=(2π,2π),res=(16,16));
coeffs  = Coefficients(Ω=0.0,θ=0.0,μ=0.0,ν=0.02,ν₄=0.0); # β=0
forcing = Kolmogorov(A₁=-1.0,A₄=-2.0);
prob    = BetaPlane(domain,coeffs,forcing);

eqs     = [GQL(0),CE2(poscheck=true,poscheckat=5)];

tsargs  = (
            dt=0.0005,
            adaptive=false,
            progress=true,
            progress_steps=50000,
            save_everystep=false,
            save_start=true,
            # dense=true,
            save_noise=false,
            saveat=5
           );

ql_kf_16      = integrate(prob,eqs[1],tspan;tsargs...);
write(prob,eqs[1],ql_kf_16,dn=dn,fn="ql_kf_16x16_1000");
@save dn*"ql_kf_16x16_500"*".jld2" ql_kf_16

ce2_kf_16     = integrate(prob,eqs[2],tspan;tsargs...);
write(prob,eqs[2],ce2_kf_16,dn=dn,fn="ce2_kf_16x16_500");
@save dn*"ce2_kf_16x16_500"*".jld2" ce2_kf_16

ce2_kf_16_qlic     = integrate(prob,eqs[2],tspan;u0=ql_kf_16.u[end],tsargs...);
write(prob,eqs[2],ce2_kf_16_qlic,dn=dn,fn="ce2_kf_16x16_qlic_500");
@save dn*"ce2_kf_16x16_qlic_500"*".jld2" ce2_kf_16_qlic

# Kolmogorov Forcing: 24x24 grid
domain  = Domain(extent=(2π,2π),res=(24,24));
coeffs  = Coefficients(Ω=0.0,θ=0.0,μ=0.0,ν=0.02,ν₄=0.0); # β=0
forcing = Kolmogorov(A₁=-1.0,A₄=-2.0);
prob    = BetaPlane(domain,coeffs,forcing);

eqs     = [GQL(0),CE2()];

tsargs  = (
            dt=0.001,
            adaptive=false,
            progress=true,
            progress_steps=5000,
            save_everystep=false,
            save_start=true,
            # dense=true,
            save_noise=false,
            saveat=20
           );

ql_kf_24      = integrate(prob,eqs[1],tspan;tsargs...);
write(prob,eqs[1],ql_kf_24,dn=dn,fn="ql_kf_24x24_500");
@save dn*ql_kf_24x24_500*".jld2" ql_kf_24

ce2_kf_24     = integrate(prob,eqs[2],tspan;tsargs...);
write(prob,eqs[2],ce2_kf_24,dn=dn,fn="ce2_kf_24x24_500");

## Stochastic jet
domain  = Domain(extent=(2π,2π),res=(8,8));
coeffs  = Coefficients(Ω=5.0,θ=0.0,μ=0.01,ν=0.0,ν₄=1.0,linear=true);
forcing = Stochastic(kf=5,dk=2,ε=0.005,τ=0.0,isotropic=false);
prob    = BetaPlane(domain,coeffs,forcing);

ql_sf      = integrate(prob,eqs[1],tspan;tsargs...)
write(prob,eqs[1],ql_sf,dn=dn,fn="ql_sf_test_500")

ce2_sf     = integrate(prob,eqs[2],tspan;tsargs...)
write(prob,eqs[2],ce2_sf,dn=dn,fn="ce2_sf_test_500")

ce2_sf_qlic     = integrate(prob,eqs[2],tspan;u0=ql_sf.u[end],tsargs...)
write(prob,eqs[2],ce2_sf_qlic,dn=dn,fn="ce2_sf_qlic_500")

#
domain  = Domain(extent=(2π,2π),res=(12,12));
coeffs  = Coefficients(Ω=5.0,θ=0.0,μ=0.01,ν=0.0,ν₄=1.0,linear=false);
forcing = Stochastic(kf=5,dk=2,ε=0.005,τ=0.0,isotropic=false);
prob    = BetaPlane(domain,coeffs,forcing);

ql_sf      = integrate(prob,eqs[1],tspan;tsargs...)
write(prob,eqs[1],ql_sf,dn=dn,fn="ql_sf_test_1000")

ce2_sf     = integrate(prob,eqs[2],tspan;tsargs...)
write(prob,eqs[2],ce2_sf,dn=dn,fn="ce2_sf_test_1000")

tspan = (0.0,500.0)
domain  = Domain(extent=(2π,2π),res=(12,12));
coeffs  = Coefficients(Ω=5.0,θ=0.0,μ=0.01,ν=0.0,ν₄=1.0,linear=false);
forcing = Stochastic(kf=8,dk=2,ε=0.01,τ=0.0,isotropic=false);
prob    = BetaPlane(domain,coeffs,forcing);

ql_sf      = integrate(prob,eqs[1],tspan;tsargs...)
write(prob,eqs[1],ql_sf,dn=dn,fn="ql_sf_test_1000")

ce2_sf     = integrate(prob,eqs[2],tspan;tsargs...)
write(prob,eqs[2],ce2_sf,dn=dn,fn="ce2_sf_test_1000")

ce2_sf_qlic     = integrate(prob,eqs[2],tspan;u0=ql_sf.u[end],tsargs...)
write(prob,eqs[2],ce2_sf_qlic,dn=dn,fn="ce2_sf_qlic_test_1000")


text:
using DashBootstrapComponents, DashHtmlComponents

row = dbc_row(
    [
        dbc_col(html_div("One of three columns")),
        dbc_col(html_div("One of three columns")),
        dbc_col(html_div("One of three columns")),
    ],
    no_gutters=true,
);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
function output_writer(u::EvolVars, chart2D::Chart, atlas, tinfo::Jecco.TimeInfo,
                       io::InOut, potential::Potential)
    Nsys = length(atlas)

    # output structures
    out_bulk  = Jecco.Output(io.out_dir, "bulk_", tinfo)
    bulkevols = getbulkevolvedpartition(u)

    bulkevols_fields = ntuple(i -> (
        Jecco.Field("phi c=$i", getphi(bulkevols[i]), atlas[i])
    ), Nsys)

    last_output_bulk_t     = -io.out_bulk_every_t

    function (u::EvolVars)
        bulkevols = getbulkevolvedpartition(u)

        it = tinfo.it
        tt = tinfo.t

        do_output_bulk     = false

        if (io.out_bulk_every > 0 && it % io.out_bulk_every == 0)
            do_output_bulk = true
        end

        if (io.out_bulk_every_t > 0 &&
            tt >= last_output_bulk_t + io.out_bulk_every_t - 1e-12)
            do_output_bulk     = true
            last_output_bulk_t = tt
        end

        if do_output_bulk
            @inbounds for i in 1:Nsys
                bulkevols_fields[i].data = getphi(bulkevols[i])
            end
            # write data
            out_bulk.(bulkevols_fields)
        end

        nothing
    end
end

function checkpoint_writer(u::EvolVars, chart2D::Chart, atlas, tinfo::Jecco.TimeInfo,
                           io::InOut)
    Nsys = length(atlas)

    # output structure
    out  = Jecco.Output(io.checkpoint_dir, "checkpoint_it", tinfo)

    bulkevols = getbulkevolvedpartition(u)

    # output fields
    bulkevols_fields = ntuple(i -> (
        Jecco.Field("phi c=$i", getphi(bulkevols[i]), atlas[i])
    ), Nsys)

    function (u::EvolVars)
        bulkevols = getbulkevolvedpartition(u)

        @inbounds for i in 1:Nsys
            bulkevols_fields[i].data = getphi(bulkevols[i])
        end

        # write data
        out.(bulkevols_fields)

        nothing
    end
end

# restore all bulk evolved fields
function restore!(bulkevols::BulkPartition{Nsys}, ts::OpenPMDTimeSeries,
                  it::Int) where {Nsys}
    for i in 1:Nsys
        phi, chart = get_field(ts, it=it, field="phi c=$i")
        @assert size(phi) == size(bulkevols[i].phi)
        copyto!(bulkevols[i].phi, phi)
    end

    nothing
end

# recover all necessary values
function recover(bulkevols::BulkPartition, boundary::Boundary, systems::SystemPartition,
                 id::InitialData, recovery_dir::String)
    ts = try
        OpenPMDTimeSeries(recovery_dir; prefix="checkpoint_it")
    catch e
        throw(e)
    end
    # grab last iteration of the timeseries, which should be the latest checkpoint
    it = ts.iterations[end]
    println("INFO: Recovering from checkpoint file $(ts.files[end])")

    restore!(bulkevols, ts, it)
    # for the boundary data
    init_data!(boundary, systems[1], id)

    ts.current_iteration, ts.current_t
end


text:
# Soma S Dhavala
# TBD


text:
using Test
using PotentialCalculation.Atoms

@testset "atoms" begin

    a = AtomOnlySymbol("H")

    am = convert(AtomWithMass, a)
    aa = convert(AtomOnlySymbol, "H")

    @test a.id == am.id

    x = AtomWithMass("x", 3.2)

    O = AtomWithMass("O")

    @test am.mass == masses["H"]
    @test O.mass == masses["O"]
end


text:
using Test
using JSON
using Dates
using Configurations
using IBMQClient.Schema

d = JSON.parse("""
{
"name": "cx",
"parameters": [],
"coupling_map": [[0,1],[0,2],[0,3],[1,2],[0,4]],
"qasm_def": "gate cx q1,q2 { CX q1, q2; }",
"conditional": true,
"latency_map": [[1,0,0,0,0],[1,0,0,0,0],
[1,0,0,0,0],[0,1,0,0,0],
[0,0,0,0,1]],
"description": "CNOT gate"
}
""")

info = Configurations.from_dict(GateInfo, d)
@test info.name == "cx"
@test info.parameters == String[]
@test info.coupling_map == [[0,1],[0,2],[0,3],[1,2],[0,4]]
@test info.conditional
@test info.latency_map == [[1,0,0,0,0],[1,0,0,0,0], [1,0,0,0,0],[0,1,0,0,0], [0,0,0,0,1]]
@test info.description == "CNOT gate"

d = JSON.parse("""
{
"name": "swap",
"parameters": [],
"coupling_map": [[2,3]],
"qasm_def": "gate swap a,b { CX a,b; CX b,a; CX a,b; }",
"conditional": false,
"description": "SWAP gate"
}
""")
info = Configurations.from_dict(GateInfo, d)

@test info.name == "swap"
@test info.parameters == String[]
@test info.coupling_map == [[2,3]]
@test info.qasm_def == "gate swap a,b { CX a,b; CX b,a; CX a,b; }"
@test info.conditional == false
@test info.description == "SWAP gate"

Qobj(;
    qobj_id="bell_Qobj_07272018",
    type="QASM",
    schema_version=v"1",
    header=Dict("description"=>"Bell states"),
    config=ExpConfig(shots=1000, memory_slots=2),
    experiments=[
        Experiment(;
            header=Dict("description"=>"|11>+|00> Bell"),
            instructions=[
                Gate(name="u2", qubits=[0], params=[0.0, π]),
                Gate(name="cx", qubits=[0, 1]),
                Measure(qubits=[0, 1], memory=[0, 1]),
            ]
        ),
        Experiment(;
            header=Dict("description"=>"|01>+|10> Bell"),
            instructions=[
                Gate(name="u2", qubits=[0], params=[0.0, π]),
                Gate(name="cx", qubits=[0, 1]),
                Gate(name="u3", qubits=[0], params=[π, 0.0, π]),
                Measure(qubits=[0, 1], memory=[0, 1]),
            ]
        )
    ]
)


text:
module HTSFileFormats

export
    SAM,
    SAMHeader,
    SAMRecord,
    BAM,
    BAMRecord,
    BAI,
    Tabix,
    ismapped,
    refname,
    nextrefname,
    refid,
    nextrefid,
    position,
    nextposition,
    mappingquality,
    flag,
    templatelength,
    seqname,
    cigar,
    sequence,
    qualities,
    header,
    overlapchunks,
    # SAM flags
    SAM_FLAG_PAIRED,
    SAM_FLAG_PROPER_PAIR,
    SAM_FLAG_UNMAP,
    SAM_FLAG_MUNMAP,
    SAM_FLAG_REVERSE,
    SAM_FLAG_MREVERSE,
    SAM_FLAG_READ1,
    SAM_FLAG_READ2,
    SAM_FLAG_SECONDARY,
    SAM_FLAG_QCFAIL,
    SAM_FLAG_DUP,
    SAM_FLAG_SUPPLEMENTARY

import Bio
import Bio: Ragel
import Bio.Seq: DNASequence
import BGZFStreams: BGZFStream, VirtualOffset, virtualoffset
import BufferedStreams   # this is necessary though I don't know why.
import BufferedStreams: BufferedInputStream

include("chunk.jl")
include("sam/sam.jl")
include("bam/bam.jl")
include("tabix/tabix.jl")

end # module


text:
# Utility #
# ========================== #
function evaluate_objective(objective::JuMP.GenericAffExpr, x::AbstractVector)
    val = objective.constant
    for (var, coeff) in objective.terms
        val += coeff*x[index(var).value]
    end
    return val
end

# function calculate_objective_value!(stochasticprogram::StochasticProgram)
#     first_stage = get_stage_one(stochasticprogram)
#     objective_value = eval_objective(first_stage.obj, first_stage.colVal)
#     objective_value += calculate_subobjectives(scenarioproblems(stochasticprogram))
#     first_stage.objVal = objective_value
#     return nothing
# end
# function calculate_subobjectives(scenarioproblems::ScenarioProblems)
#     return sum([(probability(scenario)*eval_objective(subprob.obj,subprob.colVal))::Float64 for (scenario,subprob) in zip(scenarios(scenarioproblems),subproblems(scenarioproblems))])
# end
# function calculate_subobjectives(scenarioproblems::DScenarioProblems)
#     partial_subobjectives = Vector{Float64}(undef, nworkers())
#     @sync begin
#         for (i,w) in enumerate(workers())
#             @async partial_subobjectives[i] = remotecall_fetch((sp) -> calculate_subobjectives(fetch(sp)),
#                                                                w,
#                                                                scenarioproblems[w-1])
#         end
#     end
#     return sum(partial_subobjectives)
# end

function invalidate_cache!(stochasticprogram::StochasticProgram)
    cache = problemcache(stochasticprogram)
    delete!(cache, :evp)
    delete!(cache, :dep)
    cache = solutioncache(stochasticprogram)
    delete!(cache, :solution)
    return nothing
end

function remove_scenarios!(stochasticprogram::StochasticProgram, s::Integer = 2)
    remove_scenarios!(scenarioproblems(stochasticprogram, s))
    return nothing
end

function remove_decision_variables!(stochasticprogram::StochasticProgram, s::Integer)
    remove_decision_variables!(scenarioproblems(stochasticprogram, s))
    return nothing
end

function remove_subproblems!(stochasticprogram::StochasticProgram, s::Integer = 2)
    remove_subproblems!(scenarioproblems(stochasticprogram, s))
    return nothing
end

function transfer_model!(dest::StochasticProgram, src::StochasticProgram)
    empty!(dest.generator)
    merge!(dest.generator, src.generator)
    return dest
end

function supports_zero(types::Vector, provided_def::Bool)
    for vartype in types
        if !hasmethod(zero, (Type{vartype}, ))
            !provided_def && @warn "Zero not defined for $vartype. Cannot generate zero function."
            return false
        end
    end
    return true
end

function supports_expected(types::Vector, provided_def::Bool)
    for vartype in types
        if !hasmethod(+, (vartype, vartype))
            !provided_def && @warn "Addition not defined for $vartype. Cannot generate expectation function."
            return false
        end
        if !hasmethod(*, (Float64, vartype)) || Base.code_typed(*, (Float64, vartype))[1].second != vartype
            !provided_def && @warn "Scalar multiplication with Float64 not defined for $vartype. Cannot generate expectation function."
            return false
        end
    end
    return true
end

problemcache(stochasticprogram::StochasticProgram) = stochasticprogram.problemcache
solutioncache(stochasticprogram::StochasticProgram) = stochasticprogram.solutioncache
function get_problem(stochasticprogram::StochasticProgram, key::Symbol)
    haskey(stochasticprogram.problemcache, key)|| error("No $key in problem cache")
    return stochasticprogram.problemcache[key]
end
function get_stage_one(stochasticprogram::StochasticProgram)
    haskey(stochasticprogram.problemcache, :stage_1) || error("First-stage problem not generated.")
    return stochasticprogram.problemcache[:stage_1]
end
function get_stage(stochasticprogram::StochasticProgram, stage::Integer)
    stage_key = Symbol(:stage_, stage)
    haskey(stochasticprogram.problemcache, stage_key) || error("Stage problem $stage not generated.")
    return stochasticprogram.problemcache[stage_key]
end

typename(dtype::UnionAll) = dtype.body.name.name
typename(dtype::DataType) = dtype.name.name

function add_subscript(src::AbstractString, subscript::Integer)
    return @sprintf("%s%s", src, unicode_subscript(subscript))
end
add_subscript(src::Symbol, subscript::Integer) = add_subscript(String(src), subscript)

function unicode_subscript(subscript::Integer)
    if subscript < 0
        error("$subscript is negative")
    end
    return join('₀'+d for d in reverse(digits(subscript)))
end

function extract_set(expr)
    set = NoSpecifiedConstraint()
    found = false
    quit = false
    new_expr = prewalk(expr) do x
        if @capture(x, var_Symbol)
            if var == :Bin
                set = MOI.ZeroOne()
                found = true
                return :()
            elseif var == :Int
                set = MOI.Integer()
                found = true
                return :()
            end
        elseif @capture(x, var_Symbol = true)
            if var == :binary
                set = MOI.ZeroOne()
                found = true
                return :()
            elseif var == :integer
                set = MOI.Integer()
                found = true
                return :()
            end
        elseif @capture(x, var_Symbol in constrset_) && !found && !quit
            set = constrset
            found = true
            return :($var)
        elseif @capture(x, var_Symbol[ids__] in constrset_) && !found && !quit
            set = constrset
            found = true
            return :($var[$ids])
        elseif @capture(x, set = constrset_) && !found && !quit
            set = constrset
            found = true
            return :()
        elseif @capture(x, var_Symbol[ids__])
            # Break here to prevent indices from being filtered
            quit = true
            return x
        else
            return x
        end
    end
    return set, new_expr, found
end
# ========================== #


text:
module IndependentHypothesisWeighting

using Reexport

using CategoricalArrays
using Distributions
import Distributions: cdf, pdf
using Intervals
import MLDataPattern: kfolds, FoldsView, shuffleobs

@reexport using MultipleTesting
import MultipleTesting: adjust

using Random
using SplitApplyCombine
using Statistics
using StatsBase
import StatsBase: fit, sample, weights
using Tables
using UnPack

include("grenander.jl")
include("weighted_multiple_testing.jl")
include("ihw.jl")
include("categorical_weight_learners.jl")
include("discretization_weight_learners.jl")

export PriorityWeights,
    IHW,
    GBHLearner,
    GrenanderLearner,
    QuantileSlicingDiscretizer,
    DiscretizationWeightLearner
end