text:
using JuliaSyntax:
    julia_string_to_number,
    triplequoted_string_indentation,
    unescape_julia_string,
    process_triple_strings!

hexint(s) = julia_string_to_number(s, K"HexInt")
binint(s) = julia_string_to_number(s, K"BinInt")
octint(s) = julia_string_to_number(s, K"OctInt")

@testset "Number parsing" begin
    # Integers
    @testset "Integers" begin
        @test julia_string_to_number("-1", K"Integer") isa Int
        @test julia_string_to_number("1", K"Integer") isa Int
        @test julia_string_to_number("2147483647", K"Integer") isa Int
        @test julia_string_to_number("9223372036854775807", K"Integer") isa Int64
        @test julia_string_to_number("9223372036854775808", K"Integer") isa Int128
        @test julia_string_to_number("170141183460469231731687303715884105727", K"Integer") isa Int128
        @test julia_string_to_number("170141183460469231731687303715884105728", K"Integer") isa BigInt
    end

    # Floats
    @testset "Floats" begin
        @test julia_string_to_number("10e-0", K"Float") === Float64(10)
        @test julia_string_to_number("10f-0", K"Float") === Float32(10)
        @test julia_string_to_number("0x0ap-0", K"Float") === Float64(10)
        @test julia_string_to_number("0xffp-0", K"Float") === Float64(255)
    end

    # HexInt
    @testset "HexInt numeric limits for different types" begin
        @test hexint("0xff")  === UInt8(0xff)
        @test hexint("0x100") === UInt16(0x100)
        @test hexint("0xffff") === UInt16(0xffff)
        @test hexint("0x10000") === UInt32(0x10000)
        @test hexint("0xffffffff") === UInt32(0xffffffff)
        @test hexint("0x100000000") === UInt64(0x100000000)
        @test hexint("0xffffffffffffffff") === UInt64(0xffffffffffffffff)
        @test hexint("0x10000000000000000") === UInt128(0x10000000000000000)
        @test hexint("0xffffffffffffffffffffffffffffffff") === UInt128(0xffffffffffffffffffffffffffffffff)
        @test (n = hexint("0x100000000000000000000000000000000");
               n isa BigInt && n == 0x100000000000000000000000000000000)
    end
    @testset "HexInt string length limits for different types" begin
        @test hexint("0x00")  === UInt8(0)
        @test hexint("0x000")  === UInt16(0)
        @test hexint("0x0000")  === UInt16(0)
        @test hexint("0x00000")  === UInt32(0)
        @test hexint("0x00000000") === UInt32(0)
        @test hexint("0x000000000") === UInt64(0)
        @test hexint("0x0000000000000000") === UInt64(0)
        @test hexint("0x00000000000000000") === UInt128(0)
        @test hexint("0x00000000000000000000000000000000") === UInt128(0)
        @test (n = hexint("0x000000000000000000000000000000000");
               n isa BigInt && n == 0)
    end

    # BinInt
    @testset "BinInt numeric limits for different types" begin
        @test binint("0b11111111")  === UInt8(0xff)
        @test binint("0b100000000") === UInt16(0x100)
        @test binint("0b1111111111111111") === UInt16(0xffff)
        @test binint("0b10000000000000000") === UInt32(0x10000)
        @test binint("0b11111111111111111111111111111111") === UInt32(0xffffffff)
        @test binint("0b100000000000000000000000000000000") === UInt64(0x100000000)
        @test binint("0b1111111111111111111111111111111111111111111111111111111111111111") === UInt64(0xffffffffffffffff)
        @test binint("0b10000000000000000000000000000000000000000000000000000000000000000") === UInt128(0x10000000000000000)
        @test binint("0b11111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111") === UInt128(0xffffffffffffffffffffffffffffffff)
        @test (n = binint("0b100000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000");
               n isa BigInt && n == 0x100000000000000000000000000000000)
    end
    @testset "BinInt string length limits for different types" begin
        @test binint("0b00000000")  === UInt8(0)
        @test binint("0b000000000")  === UInt16(0)
        @test binint("0b0000000000000000")  === UInt16(0)
        @test binint("0b00000000000000000")  === UInt32(0)
        @test binint("0b00000000000000000000000000000000") === UInt32(0)
        @test binint("0b000000000000000000000000000000000") === UInt64(0)
        @test binint("0b0000000000000000000000000000000000000000000000000000000000000000") === UInt64(0)
        @test binint("0b00000000000000000000000000000000000000000000000000000000000000000") === UInt128(0)
        @test binint("0b00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000") === UInt128(0)
        @test (n = binint("0b000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000");
               n isa BigInt && n == 0)
    end

    # OctInt
    @testset "OctInt numeric limits for different types" begin
        @test octint("0o377")  === UInt8(0xff)
        @test octint("0o400") === UInt16(0x100)
        @test octint("0o177777") === UInt16(0xffff)
        @test octint("0o200000") === UInt32(0x10000)
        @test octint("0o37777777777") === UInt32(0xffffffff)
        @test octint("0o40000000000") === UInt64(0x100000000)
        @test octint("0o1777777777777777777777") === UInt64(0xffffffffffffffff)
        @test octint("0o2000000000000000000000") === UInt128(0x10000000000000000)
        @test octint("0o3777777777777777777777777777777777777777777") === UInt128(0xffffffffffffffffffffffffffffffff)
        @test (n = octint("0o4000000000000000000000000000000000000000000");
               n isa BigInt && n == 0x100000000000000000000000000000000)
    end
    @testset "OctInt string length limits for different types" begin
        @test octint("0o000")  === UInt8(0)
        @test octint("0o0000")  === UInt16(0)
        @test octint("0o000000")  === UInt16(0)
        @test octint("0o0000000")  === UInt32(0)
        @test octint("0o00000000000") === UInt32(0)
        @test octint("0o000000000000") === UInt64(0)
        @test octint("0o0000000000000000000000") === UInt64(0)
        @test octint("0o00000000000000000000000") === UInt128(0)
        @test octint("0o0000000000000000000000000000000000000000000") === UInt128(0)
        @test (n = octint("0o00000000000000000000000000000000000000000000");
               n isa BigInt && n == 0)
    end

    @testset "Underscore separators" begin
        @test julia_string_to_number("10_000",      K"Integer") === 10000
        @test julia_string_to_number("10_000.0",    K"Float")   === Float64(10000)
        @test julia_string_to_number("0xff_ff",     K"HexInt")  === 0xffff
        @test julia_string_to_number("0b1111_1111", K"BinInt")  === 0xff
        @test julia_string_to_number("0o177_777",   K"OctInt")  === 0xffff
    end

    @testset "\\minus ('\\u2212' / '−') allowed in numbers" begin
        @test julia_string_to_number("−10",        K"Integer") === -10
        @test julia_string_to_number("−10.0",      K"Float")   === Float64(-10)
        @test julia_string_to_number("10e\u22121", K"Float")   === Float64(1)
    end
end

@testset "String unescaping" begin
    unesc(str) = unescape_julia_string(str, false, false)
    # Allowed escapes of delimiters and dollar sign
    @test only(unesc("\\\\")) == '\\'
    @test only(unesc("\\\"")) == '"'
    @test only(unesc("\\\$")) == '$'
    @test only(unesc("\\'"))  == '\''
    @test only(unesc("\\`"))  == '`'

    # Newline normalization
    @test unesc("a\nb\rc\r\nd") == "a\nb\nc\nd"

    # Removal of backslash-escaped newlines & indentation
    @test unesc("a\\\nb") == "ab"
    @test unesc("a\\\rb") == "ab"
    @test unesc("a\\\r\nb") == "ab"
    @test unesc("a\\\n  b") == "ab"
    @test unesc("a\\\r\n \tb") == "ab"
    @test unesc("a\\\n") == "a"
    @test unesc("a\\\r") == "a"
    @test unesc("a\\\r\n") == "a"

    # Invalid escapes
    @test_throws ArgumentError unesc("\\.")
    @test_throws ArgumentError unesc("\\z")

    # Standard C escape sequences
    @test codeunits(unesc("\\n\\t\\r\\e\\b\\f\\v\\a")) ==
        UInt8[0x0a, 0x09, 0x0d, 0x1b, 0x08, 0x0c, 0x0b, 0x07]

    # Hex and unicode escapes; \x \u and \U
    @test unesc("x\\x61x") == "xax"
    @test unesc("x\\u03b1x") == "xαx"
    @test unesc("x\\U001F604x") == "x😄x"
    # Maximum unicode code point
    @test unesc("x\\U10ffffx") == "x\U10ffffx"
    @test_throws ArgumentError unesc("x\\U110000x")

    # variable-length octal
    @test unesc("x\\7x") == "x\ax"
    @test unesc("x\\77x") == "x?x"
    @test unesc("x\\141x") == "xax"
    @test unesc("x\\377x") == "x\xffx"
    @test_throws ArgumentError unesc("x\\400x")
end

@testset "Raw string unescaping" begin
    # " delimited
    # x\"x ==> x"x
    @test unescape_julia_string("x\\\"x",     false, true) == "x\"x"
    # x\`x ==> x\`x
    @test unescape_julia_string("x\\`x",      false, true) == "x\\`x"
    # x\\\"x ==> x\"x
    @test unescape_julia_string("x\\\\\\\"x", false, true) == "x\\\"x"
    # x\\\`x ==> x\\\`x
    @test unescape_julia_string("x\\\\\\`x",  false, true) == "x\\\\\\`x"
    # '\\ ' ==> '\\ '
    @test unescape_julia_string("\\\\ ",      false, true) == "\\\\ "
    # '\\' ==> '\'
    @test unescape_julia_string("\\\\",       false, true) == "\\"
    # '\\\\' ==> '\\'
    @test unescape_julia_string("\\\\\\\\",       false, true) == "\\\\"

    # ` delimited
    # x\"x ==> x\"x
    @test unescape_julia_string("x\\\"x",     true, true) == "x\\\"x"
    # x\`x ==> x`x
    @test unescape_julia_string("x\\`x",      true, true)  == "x`x"
    # x\\\"x ==> x\"x
    @test unescape_julia_string("x\\\\\\\"x", true, true) == "x\\\\\\\"x"
    # x\\\`x ==> x\`x
    @test unescape_julia_string("x\\\\\\`x",  true, true) == "x\\`x"
    # '\\ ' ==> '\\ '
    @test unescape_julia_string("\\\\ ",      true, true) == "\\\\ "
end

@testset "Triple quoted string indentation" begin
    # Alias for non-raw triple str indentation
    triplestr_indent(str) = triplequoted_string_indentation(str, false)

    @test triplestr_indent([]) == 0

    # Spaces or tabs acceptable
    @test triplestr_indent(["\n  "]) == 2
    @test triplestr_indent(["\n\t "]) == 2
    @test triplestr_indent(["\n \t"]) == 2
    @test triplestr_indent(["\n\t\t"]) == 2

    # Start of the string is not indentation, as it's always preceded by a
    # delimiter in the source
    @test triplestr_indent(["  "]) == 0
    @test triplestr_indent(["  ", "  "]) == 0

    # Various newlines are allowed. empty lines are ignored
    @test triplestr_indent(["\n\n  x"]) == 2
    @test triplestr_indent(["\n\r  x"]) == 2
    @test triplestr_indent(["\r\n  x"]) == 2
    @test triplestr_indent(["\r\r  x"]) == 2
    @test triplestr_indent(["\n\r\r\n"]) == 0

    # Empty line at the end of any chunk implies the next source line started
    # with a delimiter, yielding zero indentation
    @test triplestr_indent(["  \n"]) == 0
    @test triplestr_indent(["  \r"]) == 0
    @test triplestr_indent(["  \n\n"]) == 0
    @test triplestr_indent(["  ", "  \n"]) == 0
    @test triplestr_indent(["  \n", "  "]) == 0

    # Find the minimum common prefix in one or several chunks
    @test triplestr_indent(["\n  ", "\n  "]) == 2
    @test triplestr_indent(["\n ", "\n  "]) == 1
    @test triplestr_indent(["\n  ", "\n "]) == 1
    @test triplestr_indent(["\n  ", "\n "]) == 1
    @test triplestr_indent(["\n \t", "\n  "]) == 1
    @test triplestr_indent(["\n  ", "\n \t"]) == 1
    @test triplestr_indent(["\n \t", "\n \t"]) == 2
    @test triplestr_indent(["\n\t ", "\n\t "]) == 2
    @test triplestr_indent(["\n  \n  "]) == 2
    @test triplestr_indent(["\n  \n "]) == 1
    @test triplestr_indent(["\n \n  "]) == 1
    # Increasing widths
    @test triplestr_indent(["\n\n \n  \n   "]) == 1
    # Decreasing widths
    @test triplestr_indent(["\n   \n  \n "]) == 1

    # Some cases of no indentation
    @test triplestr_indent(["hi"]) == 0
    @test triplestr_indent(["x\ny", "z"]) == 0

    # Escaped newlines
    @test triplestr_indent(["\\\n  "]) == 0
    @test triplestr_indent(["\\\r  "]) == 0
    @test triplestr_indent(["\\\r\n  "]) == 0
    @test triplestr_indent(["\\\r\n  "]) == 0
    @test triplestr_indent(["\n  \\\n "]) == 2
    @test triplestr_indent(["\n \\\n  "]) == 1

    # Raw strings don't have escaped newline processing
    @test triplequoted_string_indentation(["\n  \\\n "], true) == 1
    @test triplequoted_string_indentation(["\n \\\n  "], true) == 1
end

@testset "Triple quoted string deindentation" begin
    # Weird thing I noticed: In Julia 1.7 this @testset for loop adds an
    # absurd amount of testing latency given how trivial it is. Why? Is it
    # because of compiler heuristics which try to compile all for loops?
    @testset "Raw=$raw" for raw in (false, true)
        # Various combinations of dedent + leading newline stripping
        @test process_triple_strings!(["\n x", "\n y"], raw)        == ["x", "\ny"]
        @test process_triple_strings!(["\n\tx", "\n\ty"], raw)      == ["x", "\ny"]
        @test process_triple_strings!(["\r x", "\r y"], raw)        == ["x", "\ny"]
        @test process_triple_strings!(["\r x\r y"], raw)            == ["x\ny"]
        @test process_triple_strings!(["\r x\r\r y"], raw)          == ["x\n\ny"]
        @test process_triple_strings!(["\n \t x", "\n \t y"], raw)  == ["x", "\ny"]
        @test process_triple_strings!(["x\n\n y", "\n z"], raw)     == ["x\n\ny", "\nz"]
        # Cases of no dedent + newline normalization
        @test process_triple_strings!(["\n x", "\ny"], raw) == [" x", "\ny"]
        @test process_triple_strings!(["\nx", "\n y"], raw) == ["x", "\n y"]
        @test process_triple_strings!(["\n y\n"], raw) == [" y\n"]
        @test process_triple_strings!(["\n y\r"], raw) == [" y\n"]
    end
end

@testset "Normalization of identifiers" begin
    # NFC normalization
    # https://github.com/JuliaLang/julia/issues/5434
    # https://github.com/JuliaLang/julia/pull/19464
    @test JuliaSyntax.normalize_identifier("\u0069\u0302") == "\u00ee"

    # Special Julia normalization
    # https://github.com/JuliaLang/julia/pull/42561
    @test JuliaSyntax.normalize_identifier("julia\u025B\u00B5\u00B7\u0387\u2212") ==
        "julia\u03B5\u03BC\u22C5\u22C5\u002D"
end


text:
include("../src/thermo.jl")
include("../src/physics.jl")
include("../src/SimulationParameters.jl")

using PyPlot

using .SimulationParameters
sp = SimulationParameters

p = Parameters(sp.nsteps, sp.nx, sp.ny, sp.L, sp.Δx,
               sp.Δy, sp.Δt, sp.φ, sp.K, sp.Pin,
               sp.Pout, sp.P₀, sp.ψ, sp.ψ₀, sp.M, sp.μ)

function plot_function(f::Function)
    fig = PyPlot.figure(figsize=(10, 10))
    x = collect(-1e5:1e4:1e6)
    y = f.(x)
    grid("on")
    # vlines(1e8, -1000, 10000)
    plot(x, y)
    show()
end

function plot_function(f::Function, g::Function)
    fig = PyPlot.figure(figsize=(10, 10))
    x = collect(-1e5:1e4:1e6)
    grid("on")

    plot(x, f.(x))
    plot(x, f.(x))
    plot(x, g.(x))

    show()
end

function test_find_pressure()

    function f(P::AbstractFloat)
        igas = ideal_gas
        tait = tait_C₅H₁₂
        ρ̂₁ = density(igas, P)
        s = ρ₁ / ρ̂₁
        ρ̂₂ = ρ₂ / (1 - s)
        return (ρ̂₂ - tait.ρ₀) / ρ̂₂ - tait.C *
            log10((tait.B + P) / (tait.B + tait.P₀))
    end

    function fder(P::AbstractFloat)
        igas=ideal_gas
        tait=tait_C₅H₁₂
        ρ̂₁ = density(igas, P)
        s = ρ₁ / ρ̂₁
        ρ̂₂ = ρ₂ / (1 - s)
        return -tait.ρ₀ / ρ̂₂^2  * R * igas.T * igas.M *
            ρ₁ * ρ₂ / (igas.M * P - R * igas.T * ρ₁)^2 -
            tait.C / (tait.B + P) / log(10)
    end

    ρ₁ = 1.125
    ρ₂ = 2.898
    P, s = 0, 0

# -------------------- Binary Search --------------------- #
#
#     println("Binary Search")
#
#     @time begin
#     P, s = find_pressure(; root_finder=binary_search,
#                          ρ₁=ρ₁, ρ₂=ρ₂, p=p,
#                          igas=ideal_gas, tait=tait_C₅H₁₂)
#     end
#
#     println("P = ", P)
#     println("s = ", s)
# -------------------------------------------------------- #

# -------------------- Secant Method --------------------- #
#
#     println("Secant Method")
#
#     @time begin
#     P, s = find_pressure(; root_finder=secant_root_finder,
#                          ρ₁=ρ₁, ρ₂=ρ₂, p=p,
#                          igas=ideal_gas, tait=tait_C₅H₁₂)
#     end
#
#     println("P = ", P)
#     println("s = ", s)
# -------------------------------------------------------- #

    # plot_function(f, fder)
    # plot_function(fder)
# -------------------- Newton Raphson -------------------- #

    println("Newton Raphson")

    @time begin
    P, s = find_pressure(; root_finder=newton_raphson,
                         ρ₁=ρ₁, ρ₂=ρ₂, p=p,
                         igas=ideal_gas, tait=tait_C₅H₁₂)
    end

    println("P = ", P)
    println("s = ", s)
# -------------------------------------------------------- #

    # plot_function(f, P)
end

test_find_pressure()

ext: jl

max_stars_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_issues_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_forks_repo_name: sofiabelen/Two-Phase-Filtration


text:
"""
```
function chand_recursion(y::Matrix{S}, T::Matrix{S}, R::Matrix{S}, C::Vector{S},
                         Q::Matrix{S}, Z::Matrix{S}, D::Vector{S}, H::Matrix{S},
                         s_pred::Vector{S} = Vector{S}(0),
                         P_pred::Matrix{S} = Matrix{S}(0,0);
                         allout::Bool = false, Nt0::Int = 0,
                         tol::Float64 = 0.0) where {S<:AbstractFloat}
```
Implementation of Chandrasekhar Recursions, a faster algorithm for likelihood
evaluation in a linear Gaussian state space system. Based on code by Ed Herbst.

State Transition equation:
```
s_t = TT * s_{t-1} + RR * ε_t, ε_t ~ iid N(0, QQ)
```
Measurement equation:
```
y_t = DD + ZZ * s_t + η_t, η_t ~ iid N(0, HH)
```
Returns log-likelihood evaluation.
"""
function chand_recursion(y::Matrix{S}, T::Matrix{S}, R::Matrix{S}, C::Vector{S},
                         Q::Matrix{S}, Z::Matrix{S}, D::Vector{S}, H::Matrix{S},
                         s_pred::Vector{S} = Vector{S}(0),
                         P_pred::Matrix{S} = Matrix{S}(0,0);
                         allout::Bool = false, Nt0::Int = 0,
                         tol::Float64 = 0.0) where {S<:AbstractFloat}
    converged = false

    # Dimensions
    Ns     = size(T, 1) # number of states
    Ny, Nt = size(y)    # number of periods of data

    # Initialize s_pred and P_pred
    if isempty(s_pred) || isempty(P_pred)
        s_pred, P_pred = init_stationary_states(T, R, C, Q)
    end

    # V_{t|t-1} = Var(y_t|y_{1:t-1})
    V_pred    = Z * P_pred * Z' + H
    V_pred    = 0.5 * (V_pred + V_pred')
    invV_pred = inv(V_pred)

    invV_t1 = invV_pred

    # We write ΔP in terms of W_t and M_t: ΔP = W_t * M_t * W_t'
    W_t      = T * P_pred * Z'   # (Eq 12)
    M_t      = -invV_pred        # (Eq 13)
    kal_gain = W_t * invV_pred

    # Initialize loglikelihoods to zeros so that the ones we don't update
    # (those in presample) contribute zero towards sum of loglikelihoods.
    loglh    = Vector{Float64}(Nt)
    zero_vec = zeros(Ny)
    ν_t      = zero_vec
    P        = P_pred

    for t in 1:Nt
        # Step 1: Compute forecast error, ν_t and evaluate likelihoood
        yhat = Z * s_pred + D
        ν_t  = y[:, t] - yhat

        # Symptom of a bad parameter draw: simply tells callee to redraw
        if !isposdef(V_pred)
            return NaN
        else
            loglh[t] = logpdf(MvNormal(zero_vec, V_pred), ν_t)
        end

        # Step 2: Compute s_{t+1} using Eq. 5
        if t < Nt
            s_pred = T * s_pred + kal_gain * ν_t
            if allout
                P = P + W_t * M_t * W_t'
            end
        end

        if !converged
            # Intermediate calculations to re-use
            ZW_t  = Z * W_t
            MWpZp = M_t * (ZW_t')
            TW_t  = T * W_t

            # Step 3: Update forecast error variance F_{t+1} (Eq 19)
            V_t1      = V_pred
            invV_t1   = invV_pred
            V_pred    = V_pred + ZW_t * MWpZp        # F_{t+1}
            V_pred    = 0.5 * (V_pred + V_pred')
            invV_pred = inv(V_pred)

            # Step 4: Update Kalman Gain (Eq 20). Recall kalgain = K_t * V_t⁻¹
            # Kalgain_{t+1} = (Kalgain_t*V_{t-1} + T*W_t*M_t*W_t'*Z')V_t⁻¹
            kal_gain1 = kal_gain
            kal_gain  = (kal_gain * V_t1 + TW_t * MWpZp) * invV_pred

            # Step 5: Update W
            W_t = TW_t - kal_gain * ZW_t            # W_{t+1}

            # Step 6: Update M
            M_t = M_t + MWpZp * invV_t1 * MWpZp'    # M_{t+1}
            M_t = 0.5 * (M_t + M_t')

            # 'Fast Kalman' Algorithm - shut off when tol = 0
            if (maximum(abs.(kal_gain - kal_gain1)) < tol)
                converged = true
            end
        end
    end

    loglh = remove_presample!(Nt0, loglh)

    if allout
        s_TT = s_pred + P' * Z' * invV_t1 * (ν_t)
        P_TT = P - P' * Z' * invV_t1 * Z * P
        return sum(loglh), loglh, s_TT, P_TT
    else
        return sum(loglh)
    end
end

function remove_presample!(Nt0::Int, loglh::Vector{S}) where {S<:AbstractFloat}
    out = remove_presample!(Nt0, loglh, Array{Float64, 2}(0,0),
                            Array{Float64, 3}(0,0,0), Array{Float64, 2}(0,0),
                            Array{Float64, 3}(0,0,0), outputs = [:loglh])
    return out[1]
end


text:
"""
A structure that stores information about comprehension accuracy.
"""
struct Comp_Acc_Struct
    dfr::DataFrame
    acc::Float64
    err::Array
end

"""
Assess model accuracy on the basis of the correlations of row vectors of Chat and
C or Shat and S. Ideally the target words have highest correlations on the diagonal
of the pertinent correlation matrices. Homophones support option is implemented.
"""
function eval_SC end

"""
Assess model accuracy on the basis of the correlations of row vectors of Chat and
C or Shat and S. Count it as correct if one of the top k candidates is correct.
Homophones support option is implemented.
"""
function eval_SC_loose end

"""
    accuracy_comprehension(S, Shat, data)

Evaluate comprehension accuracy.

# Obligatory Arguments
- `S::Matrix`: the (gold standard) S matrix
- `Shat::Matrix`: the (predicted) Shat matrix
- `data::DataFrame`: the dataset

# Optional Arguments
- `target_col::Union{String, Symbol}=:Words`: the column name for target strings
- `base::Vector=nothing`: base features (typically a lexeme)
- `inflections::Union{Nothing, Vector}=nothing`: other features (typically in inflectional features)

# Examples
```julia
accuracy_comprehension(
    S_train,
    Shat_train,
    latin_val,
    target_col=:Words,
    base=[:Lexeme],
    inflections=[:Person, :Number, :Tense, :Voice, :Mood]
    )

accuracy_comprehension(
    S_val,
    Shat_val,
    latin_train,
    target_col=:Words,
    base=["Lexeme"],
    inflections=[:Person, :Number, :Tense, :Voice, :Mood]
    )
```
"""
function accuracy_comprehension(
    S,
    Shat,
    data;
    target_col = :Words,
    base = nothing,
    inflections = nothing,
)

    corMat = cor(Shat, S, dims = 2)
    top_index = [i[2] for i in argmax(corMat, dims = 2)]

    dfr = DataFrame()
    dfr.target = data[:, target_col]
    dfr.form = vec([data[i, target_col] for i in top_index])
    dfr.r =
        vec([corMat[index, value] for (index, value) in enumerate(top_index)])
    dfr.r_target = corMat[diagind(corMat)]
    dfr.correct = [dfr.target[i] == dfr.form[i] for i = 1:size(dfr, 1)]

    if !isnothing(inflections)
        all_features = vcat(base, inflections)
    else
        all_features = base
    end

    for f in all_features
        dfr.tmp = vec([
            data[index, f] == data[value, f]
            for (index, value) in enumerate(top_index)
        ])
        rename!(dfr, "tmp" => f)
    end

    acc = sum(dfr[:, "correct"]) / size(dfr, 1)
    err = findall(x -> x != 1, dfr[:, "correct"])

    Comp_Acc_Struct(dfr, acc, err)
end

"""
    eval_SC(SChat::AbstractArray, SC::AbstractArray)

Assess model accuracy on the basis of the correlations of row vectors of Chat and
C or Shat and S. Ideally the target words have highest correlations on the diagonal
of the pertinent correlation matrices.

# Obligatory Arguments
- `SChat::Union{SparseMatrixCSC, Matrix}`: the Chat or Shat matrix
- `SC::Union{SparseMatrixCSC, Matrix}`: the C or S matrix

# Optional Arguments
- `digits`: the specified number of digits after the decimal place (or before if negative)
- `R::Bool=false`: if true, pairwise correlation matrix R is return

```julia
eval_SC(Chat_train, cue_obj_train.C)
eval_SC(Chat_val, cue_obj_val.C)
eval_SC(Shat_train, S_train)
eval_SC(Shat_val, S_val)
```
"""
function eval_SC(SChat::AbstractArray, SC::AbstractArray; digits=4, R=false)
    rSC = cor(
        convert(Matrix{Float64}, SChat),
        convert(Matrix{Float64}, SC),
        dims = 2,
    )
    v = [rSC[i[1], i[1]] == rSC[i] ? 1 : 0 for i in argmax(rSC, dims = 2)]
    acc = round(sum(v) / length(v), digits=digits)
    if R
        return acc, rSC
    else
        return acc
    end
end

"""
    eval_SC(SChat::AbstractArray, SC::AbstractArray, SC_rest::AbstractArray)

Assess model accuracy on the basis of the correlations of row vectors of Chat and
C or Shat and S. Ideally the target words have highest correlations on the diagonal
of the pertinent correlation matrices.

!!! note
    The order is important. The fist gold standard matrix has to be corresponing
    to the SChat matrix, such as `eval_SC(Shat_train, S_train, S_val)` or `eval_SC(Shat_val, S_val, S_train)`

# Obligatory Arguments
- `SChat::Union{SparseMatrixCSC, Matrix}`: the Chat or Shat matrix
- `SC::Union{SparseMatrixCSC, Matrix}`: the training/validation C or S matrix
- `SC_rest::Union{SparseMatrixCSC, Matrix}`: the validation/training C or S matrix

# Optional Arguments
- `digits`: the specified number of digits after the decimal place (or before if negative)
- `R::Bool=false`: if true, pairwise correlation matrix R is return

```julia
eval_SC(Chat_train, cue_obj_train.C, cue_obj_val.C)
eval_SC(Chat_val, cue_obj_val.C, cue_obj_train.C)
eval_SC(Shat_train, S_train, S_val)
eval_SC(Shat_val, S_val, S_train)
```
"""
function eval_SC(
    SChat::AbstractArray,
    SC::AbstractArray,
    SC_rest::AbstractArray;
    digits = 4,
    R = false
    )

    eval_SC(SChat, vcat(SC, SC_rest); digits=digits, R=R)
end

"""
    eval_SC(SChat::AbstractArray, SC::AbstractArray, data::DataFrame, target_col::Union{String, Symbol})

Assess model accuracy on the basis of the correlations of row vectors of Chat and
C or Shat and S. Ideally the target words have highest correlations on the diagonal
of the pertinent correlation matrices. Support for homophones.

# Obligatory Arguments
- `SChat::Union{SparseMatrixCSC, Matrix}`: the Chat or Shat matrix
- `SC::Union{SparseMatrixCSC, Matrix}`: the C or S matrix
- `data::DataFrame`: datasets
- `target_col::Union{String, Symbol}`: target column name

# Optional Arguments
- `digits`: the specified number of digits after the decimal place (or before if negative)
- `R::Bool=false`: if true, pairwise correlation matrix R is return

```julia
eval_SC(Chat_train, cue_obj_train.C, latin, :Word)
eval_SC(Chat_val, cue_obj_val.C, latin, :Word)
eval_SC(Shat_train, S_train, latin, :Word)
eval_SC(Shat_val, S_val, latin, :Word)
```
"""
function eval_SC(
    SChat::AbstractArray,
    SC::AbstractArray,
    data::DataFrame,
    target_col::Union{String, Symbol};
    digits = 4,
    R = false
    )

    rSC = cor(
        convert(Matrix{Float64}, SChat),
        convert(Matrix{Float64}, SC),
        dims = 2,
    )
    v = [
        data[i[1], target_col] == data[i[2], target_col] ? 1 : 0
        for i in argmax(rSC, dims = 2)
    ]
    acc = round(sum(v) / length(v), digits=digits)
    if R
        return acc, rSC
    else
        return acc
    end
end

"""
    eval_SC(SChat::AbstractArray, SC::AbstractArray, SC_rest::AbstractArray, data::DataFrame, data_rest::DataFrame, target_col::Union{String, Symbol})

Assess model accuracy on the basis of the correlations of row vectors of Chat and
C or Shat and S. Ideally the target words have highest correlations on the diagonal
of the pertinent correlation matrices.

!!! note
    The order is important. The fist gold standard matrix has to be corresponing
    to the SChat matrix, such as `eval_SC(Shat_train, S_train, S_val, latin, :Word)`
    or `eval_SC(Shat_val, S_val, S_train, latin, :Word)`

# Obligatory Arguments
- `SChat::Union{SparseMatrixCSC, Matrix}`: the Chat or Shat matrix
- `SC::Union{SparseMatrixCSC, Matrix}`: the training/validation C or S matrix
- `SC_rest::Union{SparseMatrixCSC, Matrix}`: the validation/training C or S matrix
- `data::DataFrame`: the training/validation datasets
- `data_rest::DataFrame`: the validation/training datasets
- `target_col::Union{String, Symbol}`: target column name

# Optional Arguments
- `digits`: the specified number of digits after the decimal place (or before if negative)
- `R::Bool=false`: if true, pairwise correlation matrix R is return

```julia
eval_SC(Chat_train, cue_obj_train.C, cue_obj_val.C, latin, :Word)
eval_SC(Chat_val, cue_obj_val.C, cue_obj_train.C, latin, :Word)
eval_SC(Shat_train, S_train, S_val, latin, :Word)
eval_SC(Shat_val, S_val, S_train, latin, :Word)
```
"""
function eval_SC(
    SChat::AbstractArray,
    SC::AbstractArray,
    SC_rest::AbstractArray,
    data::DataFrame,
    data_rest::DataFrame,
    target_col::Union{String, Symbol};
    digits = 4,
    R = false
    )

    data_combined = copy(data)
    append!(data_combined, data_rest)

    eval_SC(
        SChat,
        vcat(SC, SC_rest),
        data_combined,
        target_col,
        digits = digits,
        R = R
        )
end

"""
    eval_SC(SChat::AbstractArray, SC::AbstractArray, batch_size::Int64)

Assess model accuracy on the basis of the correlations of row vectors of Chat and
C or Shat and S. Ideally the target words have highest correlations on the diagonal
of the pertinent correlation matrices. For large datasets, pass batch_size to
process evaluation in chucks.

# Obligatory Arguments
- `SChat`: the Chat or Shat matrix
- `SC`: the C or S matrix
- `data`: datasets
- `target_col`: target column name
- `batch_size`: batch size

# Optional Arguments
- `digits`: the specified number of digits after the decimal place (or before if negative)
- `verbose::Bool=false`: if true, more information is printed

```julia
eval_SC(Chat_train, cue_obj_train.C, latin, :Word)
eval_SC(Chat_val, cue_obj_val.C, latin, :Word)
eval_SC(Shat_train, S_train, latin, :Word)
eval_SC(Shat_val, S_val, latin, :Word)
```
"""
function eval_SC(
    SChat::AbstractArray,
    SC::AbstractArray,
    batch_size::Int64;
    digits = 4,
    verbose = false
    )

    l = size(SChat, 1)
    num_chucks = ceil(Int64, l / batch_size)
    verbose && begin
        pb = Progress(num_chucks)
    end
    correct = 0

    SChat_d = convert(Matrix{Float64}, SChat)
    SC_d = convert(Matrix{Float64}, SC)

    # for first parts
    for j = 1:num_chucks-1
        correct += eval_SC_chucks(
            SChat_d,
            SC_d,
            (j - 1) * batch_size + 1,
            j * batch_size,
            batch_size,
        )
        verbose && ProgressMeter.next!(pb)
    end
    # for last part
    correct += eval_SC_chucks(
        SChat_d,
        SC_d,
        (num_chucks - 1) * batch_size + 1,
        batch_size,
    )
    verbose && ProgressMeter.next!(pb)

    round(correct / l, digits=digits)
end

"""
    eval_SC(SChat::AbstractArray, SC::AbstractArray, data::DataFrame, target_col::Union{String, Symbol}, batch_size::Int64)

Assess model accuracy on the basis of the correlations of row vectors of Chat and
C or Shat and S. Ideally the target words have highest correlations on the diagonal
of the pertinent correlation matrices. For large datasets, pass batch_size to
process evaluation in chucks. Support homophones.

# Obligatory Arguments
- `SChat::AbstractArray`: the Chat or Shat matrix
- `SC::AbstractArray`: the C or S matrix
- `data::DataFrame`: datasets
- `target_col::Union{String, Symbol}`: target column name
- `batch_size::Int64`: batch size

# Optional Arguments
- `digits`: the specified number of digits after the decimal place (or before if negative)
- `verbose::Bool=false`: if true, more information is printed

```julia
eval_SC(Chat_train, cue_obj_train.C, latin, :Word, 5000)
eval_SC(Chat_val, cue_obj_val.C, latin, :Word, 5000)
eval_SC(Shat_train, S_train, latin, :Word, 5000)
eval_SC(Shat_val, S_val, latin, :Word, 5000)
```
"""
function eval_SC(
    SChat::AbstractArray,
    SC::AbstractArray,
    data::DataFrame,
    target_col::Union{String, Symbol},
    batch_size::Int64;
    digits = 4,
    verbose = false
    )

    l = size(SChat, 1)
    num_chucks = ceil(Int64, l / batch_size)
    verbose && begin
        pb = Progress(num_chucks)
    end
    correct = 0

    SChat_d = convert(Matrix{Float64}, SChat)
    SC_d = convert(Matrix{Float64}, SC)

    # for first parts
    for j = 1:num_chucks-1
        correct += eval_SC_chucks(
            SChat_d,
            SC_d,
            (j - 1) * batch_size + 1,
            j * batch_size,
            batch_size,
            data,
            target_col,
        )
        verbose && ProgressMeter.next!(pb)
    end
    # for last part
    correct += eval_SC_chucks(
        SChat_d,
        SC_d,
        (num_chucks - 1) * batch_size + 1,
        batch_size,
        data,
        target_col,
    )
    verbose && ProgressMeter.next!(pb)

    round(correct / l, digits=digits)
end

function eval_SC_chucks(SChat, SC, s, e, batch_size)
    rSC = cor(SChat[s:e, :], SC, dims = 2)
    v = [(rSC[i[1], i[1]+s-1] == rSC[i]) ? 1 : 0 for i in argmax(rSC, dims = 2)]
    sum(v)
end

function eval_SC_chucks(SChat, SC, s, e, batch_size, data, target_col)
    rSC = cor(SChat[s:e, :], SC, dims = 2)
    v = [
        data[i[1]+s-1, target_col] == data[i[2], target_col] ? 1 : 0
        for i in argmax(rSC, dims = 2)
    ]
    sum(v)
end

function eval_SC_chucks(SChat, SC, s, batch_size)
    rSC = cor(SChat[s:end, :], SC, dims = 2)
    v = [(rSC[i[1], i[1]+s-1] == rSC[i]) ? 1 : 0 for i in argmax(rSC, dims = 2)]
    sum(v)
end

function eval_SC_chucks(SChat, SC, s, batch_size, data, target_col)
    rSC = cor(SChat[s:end, :], SC, dims = 2)
    v = [
        data[i[1]+s-1, target_col] == data[i[2], target_col] ? 1 : 0
        for i in argmax(rSC, dims = 2)
    ]
    sum(v)
end

"""
    eval_SC_loose(SChat, SC, k)

Assess model accuracy on the basis of the correlations of row vectors of Chat and
C or Shat and S. Count it as correct if one of the top k candidates is correct.

# Obligatory Arguments
- `SChat::Union{SparseMatrixCSC, Matrix}`: the Chat or Shat matrix
- `SC::Union{SparseMatrixCSC, Matrix}`: the C or S matrix
- `k`: top k candidates

# Optional Arguments
- `digits`: the specified number of digits after the decimal place (or before if negative)

```julia
eval_SC_loose(Chat, cue_obj.C, k)
eval_SC_loose(Shat, S, k)
```
"""
function eval_SC_loose(SChat, SC, k; digits=4)
    total = size(SChat, 1)
    correct = 0
    rSC = cor(
        convert(Matrix{Float64}, SChat),
        convert(Matrix{Float64}, SC),
        dims = 2,
    )

    for i = 1:total
        p = sortperm(rSC[i, :], rev = true)
        p = p[1:k, :]
        if i in p
            correct += 1
        end
    end
    round(correct / total, digits=digits)
end

"""
    eval_SC_loose(SChat, SC, k, data, target_col)

Assess model accuracy on the basis of the correlations of row vectors of Chat and
C or Shat and S. Count it as correct if one of the top k candidates is correct.
Support for homophones.

# Obligatory Arguments
- `SChat::Union{SparseMatrixCSC, Matrix}`: the Chat or Shat matrix
- `SC::Union{SparseMatrixCSC, Matrix}`: the C or S matrix
- `k`: top k candidates
- `data`: datasets
- `target_col`: target column name

# Optional Arguments
- `digits`: the specified number of digits after the decimal place (or before if negative)

```julia
eval_SC_loose(Chat, cue_obj.C, k, latin, :Word)
eval_SC_loose(Shat, S, k, latin, :Word)
```
"""
function eval_SC_loose(SChat, SC, k, data, target_col; digits=4)
    total = size(SChat, 1)
    correct = 0
    rSC = cor(
        convert(Matrix{Float64}, SChat),
        convert(Matrix{Float64}, SC),
        dims = 2,
    )

    for i = 1:total
        p = sortperm(rSC[i, :], rev = true)
        p = p[1:k]
        if i in p
            correct += 1
        else
            if data[i, target_col] in data[p, target_col]
                correct += 1
            end
        end
    end
    round(correct / total, digits=digits)
end

"""
    eval_manual(res, data, i2f)

Create extensive reports for the outputs from `build_paths` and `learn_paths`.
"""
function eval_manual(
    res,
    data,
    i2f;
    s = 1,
    e = nothing,
    grams = 3,
    tokenized = false,
    sep_token = nothing,
    start_end_token = "#",
    verbose = false,
    identifier = "Identifier",
    output_sep_token = "",
)

    # is end is not specified then it is the end of dataset
    if isnothing(e)
        e = length(res)
    end

    for i = s:e
        verbose && println("="^20)
        verbose && println("utterance $i: ")
        verbose && println(data[i, identifier])
        verbose && println("predictions: ")
        verbose && begin
            display_pred(
                res[i],
                i2f,
                grams,
                tokenized,
                sep_token,
                start_end_token,
                output_sep_token,
            )
        end
    end
    return
end

"""
    eval_acc(res, gold_inds::Array)

Evaluate the accuracy of the results from `learn_paths` or `build_paths`.

# Obligatory Arguments
- `res::Array`: the results from `learn_paths` or `build_paths`
- `gold_inds::Array`: the gold paths' indices

# Optional Arguments
- `digits`: the specified number of digits after the decimal place (or before if negative)
- `verbose::Bool=false`: if true, more information is printed

# Examples
```julia
# evaluation on training data
acc_train = JudiLing.eval_acc(
    res_train,
    cue_obj_train.gold_ind,
    verbose=false
)

# evaluation on validation data
acc_val = JudiLing.eval_acc(
    res_val,
    cue_obj_val.gold_ind,
    verbose=false
)
```
"""
function eval_acc(res, gold_inds::Array; digits=4, verbose = false)

    total = length(res)
    acc = 0

    iter = 1:total
    if verbose
        pb = Progress(total)
    end

    for i in iter
        if isassigned(res[i], 1) &&
           iscorrect(gold_inds[i], res[i][1].ngrams_ind)
            acc += 1
        end
        if verbose
            ProgressMeter.next!(pb)
        end
    end

    round(acc / total, digits=digits)
end

"""
    eval_acc(res, cue_obj::Cue_Matrix_Struct)

Evaluate the accuracy of the results from `learn_paths` or `build_paths`.

# Obligatory Arguments
- `res::Array`: the results from `learn_paths` or `build_paths`
- `cue_obj::Cue_Matrix_Struct`: the C matrix object

# Optional Arguments
- `digits`: the specified number of digits after the decimal place (or before if negative)
- `verbose::Bool=false`: if true, more information is printed

# Examples
```julia
acc = JudiLing.eval_acc(res, cue_obj)
```
"""
function eval_acc(res, cue_obj::Cue_Matrix_Struct; digits = 4, verbose = false)
    eval_acc(res, cue_obj.gold_ind, digits = digits, verbose = verbose)
end

"""
    eval_acc_loose(res, gold_inds)

Lenient evaluation of the accuracy of the results from `learn_paths` or `build_paths`,
counting a prediction as correct when the correlation of the predicted and gold
standard semantic vectors is among the n top correlations, where n is equal to
`max_can` in the 'learn_paths' or `build_paths` function.

# Obligatory Arguments
- `res::Array`: the results from `learn_paths` or `build_paths`
- `gold_inds::Array`: the gold paths' indices

# Optional Arguments
- `digits`: the specified number of digits after the decimal place (or before if negative)
- `verbose::Bool=false`: if true, more information is printed

# Examples
```julia
# evaluation on training data
acc_train_loose = JudiLing.eval_acc_loose(
    res_train,
    cue_obj_train.gold_ind,
    verbose=false
)

# evaluation on validation data
acc_val_loose = JudiLing.eval_acc_loose(
    res_val,
    cue_obj_val.gold_ind,
    verbose=false
)
```
"""
function eval_acc_loose(res, gold_inds; digits=4, verbose = false)

    total = length(res)
    acc = 0

    iter = 1:total
    if verbose
        pb = Progress(total)
    end

    for i in iter
        is_correct = false
        for j = 1:length(res[i])
            if isassigned(res[i], 1) &&
               iscorrect(gold_inds[i], res[i][j].ngrams_ind)
                is_correct = true
            end
        end
        if is_correct
            acc += 1
        end
        if verbose
            ProgressMeter.next!(pb)
        end
    end

    round(acc / total, digits=4)
end

"""
extract_gpi(gpi, threshold=0.1, tolerance=(-1000.0))

Extract, using gold paths' information, how many n-grams for a gold
path are below the threshold but above the tolerance.
"""
function extract_gpi(gpi, threshold = 0.1, tolerance = (-1000.0))
    c = [
        count(x -> threshold >= x > tolerance, g.ngrams_ind_support)
        for g in gpi
    ]

    t_c = length(c)

    summary_c = [
        (x, count(cnt -> cnt == x, c), count(cnt -> cnt == x, c) / t_c)
        for x in sort(unique(c))
    ]

    summary_c
end


text:
## ClimateERA Root Directory Setup

function eraroot(actionID)

    path = joinpath("$(homedir())","research","ecmwf");
    @warn "$(Dates.now()) - No directory path was given.  Setting to default path: $(path) for ClimateERA data downloads."

    if isdir(path)
        @info "$(Dates.now()) - The default path $(path) exists and therefore can be used as a directory for ClimateERA data downloads."
    else
        if actionID != 1
            error("$(Dates.now()) - The path $(path) does not exist.  If you are doing analysis, please point towards the correct path before proceeding ...")
        else
            @warn "$(Dates.now()) - The path $(path) does not exist.  A new directory will be created here.  Therefore if you already have an existing repository for ClimateERA data, make sure that $(path) is the correct location."
            @info "$(Dates.now()) - Creating path $(path) ..."
            mkpath(path);
        end
    end

    return eramkroot(path)

end

function eraroot(path::AbstractString,actionID::Integer)

    if isdir(path)
        @info "$(Dates.now()) - The default path $(path) exists and therefore can be used as a directory for ClimateERA data downloads."
    else
        if actionID != 1
            error("$(Dates.now()) - The path $(path) does not exist.  If you are doing analysis, please point towards the correct path before proceeding ...")
        else
            @warn "$(Dates.now()) - The path $(path) does not exist.  A new directory will be created here.  Therefore if you already have an existing repository for ClimateERA data, make sure that $(path) is the correct location."
            mkpath(path);
        end
    end

    return eramkroot(path)

end

function eramkroot(eroot::AbstractString)

    eiroot = joinpath(eroot,"erai");
    if !isdir(eiroot)
        mkpath(eiroot); @info "$(Dates.now()) - Created root folder for ERA-Interim reanalysis data $(eiroot)."
    else;               @info "$(Dates.now()) - Root folder for ERA-Interim reanalysis data $(eiroot) exists."
    end

    e5root = joinpath(eroot,"era5");
    if !isdir(e5root)
        mkpath(e5root); @info "$(Dates.now()) - Created root folder for ERA5 reanalysis data $(e5root)."
    else;               @info "$(Dates.now()) - Root folder for ERA5 reanalysis data $(e5root) exists."
    end

    eproot = joinpath(eroot,"eplot");
    if !isdir(eproot)
        mkpath(eproot); @info "$(Dates.now()) - Created root folder for ERA plotting data $(eproot)."
    else;               @info "$(Dates.now()) - Root folder for ERA plotting data $(eproot) exists."
    end

    return Dict("erai"=>eiroot,"era5"=>e5root,"plot"=>eproot);

end

## ERA Startup

function eraaction(actionID::Integer)
    if     actionID == 1; return Dict("ID"=>1,"name"=>"download");
    elseif actionID == 2; return Dict("ID"=>2,"name"=>"analyse");
    end
end

function eradataset(datasetID::Integer)
    if     datasetID == 1; return Dict("ID"=>1,"short"=>"era5","name"=>"ERA5");
    elseif datasetID == 2; return Dict("ID"=>2,"short"=>"erai","name"=>"ERA-Interim");
    end
end

# Startup ClimateERA

function erastartup(;
    aID::Integer,dID::Integer,
    path::AbstractString="",
    welcome::Bool=true
)

    if welcome; erawelcome(); end

    if !(aID in [1,2])
        error("$(Dates.now()) - Please input a valid action-type.  Call queryeaction() for more details.")
    end
    if !(dID in [1,2])
        error("$(Dates.now()) - Please input a valid dataset-type.  Call queryedataset() for more details.")
    end

    if path == ""; eroot = eraroot(aID); else; eroot = eraroot(path,aID); end

    if     dID == 1; eroot["era"]=eroot["era5"]; delete!(eroot,["era5","erai"]);
    elseif dID == 2; eroot["era"]=eroot["erai"]; delete!(eroot,["era5","erai"]);
    end

    action = eraaction(aID); dataset = eradataset(dID);

    @info "$(Dates.now()) - This script will $(action["name"]) $(dataset["name"]) data."
    init = Dict(
        "actionID"=>aID,"action"=>action["name"],
        "datasetID"=>dID,"dataset"=>dataset["name"],
        "prefix"=>dataset["short"]
    )

    eplotsetup(eroot)

    return init,eroot

end

function erawelcome()

    ftext = joinpath(@__DIR__,"../extra/erawelcome.txt");
    lines = readlines(ftext); count = 0; nl = length(lines);
    for l in lines; count += 1;
       if any(count .== [1,2]); print(Crayon(bold=true),"$l\n");
       elseif count == nl;      print(Crayon(bold=false),"$l\n\n");
       else;                    print(Crayon(bold=false),"$l\n");
       end
    end

end


text:
a = plot([f0 f1],
    xlabel=L"$k$ Values",
    ylabel=L"Probability of $z_k$",
    labels=reshape([L"$f_0$"; L"$f_1$"], 1, 2),
    linewidth=2,
    ylims=[0.;0.07],
    title="Original Distributions")

mix = Array(Float64, 50, 3)
p_k = [0.25; 0.5; 0.75]
labels = Array(String, 3)
for i in 1:3
    mix[:, i] = p_k[i] * f0 + (1 - p_k[i]) * f1
    labels[i] = string(L"$p_k$ = ", p_k[i])
end

b = plot(mix,
    xlabel=L"$k$ Values",
    ylabel=L"Probability of $z_k$",
    labels=reshape(labels, 1, 3),
    linewidth=2,
    ylims=[0.;0.06],
    title="Mixture of Original Distributions")

plot(a, b, layout=(2, 1), size=(600, 800))


text:
using CombinedUncertainDiffEq
using Documenter

DocMeta.setdocmeta!(CombinedUncertainDiffEq, :DocTestSetup, :(using CombinedUncertainDiffEq); recursive=true)

makedocs(;
    modules=[CombinedUncertainDiffEq],
    authors="Jonnie Diegelman <jonniediegelman@gmail.com> and contributors",
    repo="https://github.com/jonniedie/CombinedUncertainDiffEq.jl/blob/{commit}{path}#{line}",
    sitename="CombinedUncertainDiffEq.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://jonniedie.github.io/CombinedUncertainDiffEq.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/jonniedie/CombinedUncertainDiffEq.jl",
)


text:
## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-recommendation-genres

"""
    genres_get()

**Summary**: Retrieve a list of available genres seed parameter values for recommendations.

# Example
```julia-repl
julia> Spotify.genres_get()[1]
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 1 entry:
  :genres => ["acoustic", "afrobeat", "alt-rock", "alternative", "ambient", "anime", "black-metal",…
```
"""
function genres_get()

    return Spotify.spotify_request("recommendations/available-genre-seeds")

end


text:
using Statistics
function median2(n)
	s = sort(n)
	len = length(n)
	if len % 2 == 0
		return (s[floor(Int, len / 2) + 1] + s[floor(Int, len / 2)]) / 2
	else
		return  s[floor(Int, len / 2) + 1]
	end
end

a = [4.1, 5.6, 7.2, 1.7, 9.3, 4.4, 3.2]
b = [4.1, 7.2, 1.7, 9.3, 4.4, 3.2]

@show a b median2(a) median(a) median2(b) median(b)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# module PengRobinson

using PolynomialRoots

struct PengRobinsonData{T} <: AbstractCubicEOS where T <: AbstractFloat
    CompoundName::String
    MW::T
    Pc::T
    Tc::T
    a::T
    b::T
    ω::T
    κ::T
end

function α(eos::PengRobinsonData, T)
    Tr = T/eos.Tc
    return (1+eos.κ*(1-√(Tr)))^2
end
function A(eos::PengRobinsonData, P, T)
    return eos.a * α(eos, T) * P / (R() ^2 * T^2)
end
function B(eos::PengRobinsonData, P, T)
    return eos.b*P/(T*R())
end
function c1(eos::PengRobinsonData, P, T)
    return (B(eos, P, T) - 1)
end
function c2(eos::PengRobinsonData, P, T)
    return A(eos, P, T) - 2B(eos, P,T) - 3B(eos, P,T)^2
end
function c3(eos::PengRobinson, p, T)
    return B(eos, P,T)^3 + B(eos, P,T)^2 - A(eos, P,T)*B(eos, P,T)
end

function PengRobinsonData{T}(CompoundName::String, MW::T, Pc::T, Tc::T, ω::T) where T<: AbstractFloat
    a = 0.45724 * R() ^2 * Tc^2 /pc
    b = 0.0778*R() * Tc /pc
    κ = 0.37464 + 1.54226*ω - 0.26992*ω^2
    return PengRobinsonData(CompoundName, MW, Pc, Tc, a,b,ω,κ)
end

function get_the_real(vec::Vector{Complex{T}}) where {T<:Number}
    for i in vec
        if abs(imag(i)) < 2*eps(T)
            return real(i)
        end
    end
end

function volume_roots(eos::PengRobinsonData, P, T)
    c1_ = c1(eos, P, T)
    c2_ = c2(eos, P, T)
    c3_ = c3(eos, P, T)
    roots_ = roots([c3_, c2_, c1_, 1.0])
    return roots_
end

# end # ends module PengRobinson


text:
using NeidSolarScripts
using Test

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


text:
using Test
using Clustering
using Distances
using Random
using LinearAlgebra

import Distances.pairwise!

# custom distance metric
struct MySqEuclidean <: SemiMetric end

# redefinition of Distances.pairwise! for MySqEuclidean type
function pairwise!(r::AbstractMatrix, dist::MySqEuclidean,
                   a::AbstractMatrix, b::AbstractMatrix; dims::Integer=2)
    dims == 2 || throw(ArgumentError("only dims=2 supported for MySqEuclidean distance"))
    mul!(r, transpose(a), b)
    sa2 = sum(abs2, a, dims=1)
    sb2 = sum(abs2, b, dims=1)
    @inbounds r .= sa2' .+ sb2 .- 2r
end

@testset "kmeans() (k-means)" begin

@testset "Argument checks" begin
    Random.seed!(34568)
    @test_throws ArgumentError kmeans(randn(2, 3), 1)
    @test_throws ArgumentError kmeans(randn(2, 3), 4)
    @test kmeans(randn(2, 3), 2) isa KmeansResult
    @test_throws ArgumentError kmeans(randn(2, 3), 2, display=:mylog)
    for disp in keys(Clustering.DisplayLevels)
        @test kmeans(randn(2, 3), 2, display=disp) isa KmeansResult
    end
end

Random.seed!(34568)

m = 3
n = 1000
k = 10

x = rand(m, n)
xt = copy(transpose(x))

equal_kmresults(km1::KmeansResult, km2::KmeansResult) =
    all(getfield(km1, η) == getfield(km2, η) for η ∈ fieldnames(KmeansResult))

@testset "non-weighted" begin
    Random.seed!(34568)
    r = kmeans(x, k; maxiter=50)
    @test isa(r, KmeansResult{Matrix{Float64}, Float64, Int})
    @test nclusters(r) == k
    @test size(r.centers) == (m, k)
    @test length(r.assignments) == n
    @test all(a -> 1 <= a <= k, r.assignments)
    @test length(r.costs) == n
    @test length(counts(r)) == k
    @test sum(counts(r)) == n
    @test wcounts(r) == counts(r)
    @test sum(r.costs) ≈ r.totalcost

    Random.seed!(34568)
    r_t = kmeans(xt', k; maxiter=50)
    @test equal_kmresults(r, r_t)
end

@testset "non-weighted (float32)" begin
    Random.seed!(34568)
    x32 = map(Float32, x)
    x32t = copy(x32')
    r = kmeans(x32, k; maxiter=50)
    @test isa(r, KmeansResult{Matrix{Float32}, Float32, Int})
    @test nclusters(r) == k
    @test size(r.centers) == (m, k)
    @test length(r.assignments) == n
    @test all(a -> 1 <= a <= k, r.assignments)
    @test length(r.costs) == n
    @test length(counts(r)) == k
    @test sum(counts(r)) == n
    @test wcounts(r) == counts(r)
    @test sum(r.costs) ≈ r.totalcost

    Random.seed!(34568)
    r_t = kmeans(x32t', k; maxiter=50)
    @test equal_kmresults(r, r_t)
end

@testset "weighted" begin
    w = rand(n)
    Random.seed!(34568)
    r = kmeans(x, k; maxiter=50, weights=w)
    @test isa(r, KmeansResult{Matrix{Float64}, Float64, Float64})
    @test nclusters(r) == k
    @test size(r.centers) == (m, k)
    @test length(r.assignments) == n
    @test all(a -> 1 <= a <= k, r.assignments)
    @test length(r.costs) == n
    @test length(counts(r)) == k
    @test sum(counts(r)) == n

    cw = zeros(k)
    for i = 1:n
        cw[r.assignments[i]] += w[i]
    end
    @test wcounts(r) ≈ cw
    @test dot(r.costs, w) ≈ r.totalcost

    Random.seed!(34568)
    r_t = kmeans(xt', k; maxiter=50, weights=w)
    @test equal_kmresults(r, r_t)
end

@testset "custom distance" begin
    Random.seed!(34568)
    r = kmeans(x, k; maxiter=50, init=:kmcen, distance=MySqEuclidean())
    r2 = kmeans(x, k; maxiter=50, init=:kmcen)
    @test isa(r, KmeansResult{Matrix{Float64}, Float64, Int})
    @test nclusters(r) == k
    @test size(r.centers) == (m, k)
    @test length(r.assignments) == n
    @test all(a -> 1 <= a <= k, r.assignments)
    @test length(r.costs) == n
    @test length(counts(r)) == k
    @test sum(counts(r)) == n
    @test wcounts(r) == r.counts
    @test sum(r.costs) ≈ r.totalcost
    @test equal_kmresults(r, r2)

    Random.seed!(34568)
    r_t = kmeans(xt', k; maxiter=50, init=:kmcen, distance=MySqEuclidean())
    @test equal_kmresults(r, r_t)
end

@testset "Argument checks" begin
    Random.seed!(34568)
    n = 50
    k = 10
    x = randn(m, n)

    @testset "init=" begin
        @test_throws ArgumentError kmeans(x, k, init=1:(k-2))
        @test_throws ArgumentError kmeans(x, k, init=1:(k+2))
        @test kmeans(x, k, init=1:k, maxiter=5) isa KmeansResult

        @test_throws ArgumentError kmeans(x, k, init=:myseeding)
        for algname in (:kmpp, :kmcen, :rand)
            alg = Clustering.seeding_algorithm(algname)
            @test kmeans(x, k, init=algname) isa KmeansResult
            @test kmeans(x, k, init=alg) isa KmeansResult
        end
    end
end

@testset "Integer data" begin
    x = rand(Int16, m, n)
    Random.seed!(654)
    r = kmeans(x, k; maxiter=50)

    @test isa(r, KmeansResult{Matrix{Float64}, Float64, Int})
end

@testset "kmeans! data types" begin
    Random.seed!(1101)
    for TX in (Int, Float32, Float64)
        for TC in (Float32, Float64)
            for TW in (Nothing, Int, Float32, Float64)
                x = rand(TX, m, n)
                c = rand(TC, m, k)
                if TW == Nothing
                    r = kmeans!(x, c; maxiter=1)
                    @test isa(r, KmeansResult{Matrix{TC},<:Real,Int})
                else
                    w = rand(TW, n)
                    r = kmeans!(x, c; weights=w, maxiter=1)
                    @test isa(r, KmeansResult{Matrix{TC},<:Real,TW})
                end
            end
        end
    end
end

end


text:
using AbstractDifferentiation
using Distributions
using ForwardDiff
using LinearAlgebra
using GalacticOptim
using Pathfinder
using PSIS
using ReverseDiff
using SciMLBase
using Test
using Transducers

@testset "multi path pathfinder" begin
    @testset "MvNormal" begin
        dim = 10
        nruns = 20
        ndraws = 1000_000
        ndraws_per_run = ndraws ÷ nruns
        ndraws_elbo = 100
        Σ = rand_pd_mat(Float64, dim)
        μ = randn(dim)
        d = MvNormal(μ, Σ)
        logp(x) = logpdf(d, x)
        ∇logp(x) = ForwardDiff.gradient(logp, x)
        rngs = if VERSION ≥ v"1.7"
            [MersenneTwister(), Random.default_rng()]
        else
            [MersenneTwister()]
        end
        seed = 76
        @testset for rng in rngs
            executor = rng isa MersenneTwister ? SequentialEx() : ThreadedEx()

            Random.seed!(rng, seed)
            result = multipathfinder(
                logp, ∇logp, ndraws; dim, nruns, ndraws_elbo, ndraws_per_run, rng, executor
            )
            @test result isa MultiPathfinderResult
            @test result.input === (logp, ∇logp)
            @test result.optim_fun isa SciMLBase.OptimizationFunction
            @test result.rng === rng
            @test result.optimizer ===
                Pathfinder.default_optimizer(Pathfinder.DEFAULT_HISTORY_LENGTH)
            @test result.fit_distribution isa MixtureModel
            @test ncomponents(result.fit_distribution) == nruns
            @test Distributions.component_type(result.fit_distribution) <: MvNormal
            @test result.draws isa AbstractMatrix
            @test size(result.draws) == (dim, ndraws)
            @test result.draw_component_ids isa Vector{Int}
            @test length(result.draw_component_ids) == ndraws
            @test extrema(result.draw_component_ids) == (1, nruns)
            @test result.fit_distribution_transformed === result.fit_distribution
            @test result.draws_transformed == result.draws
            @test result.pathfinder_results isa Vector{<:PathfinderResult}
            @test length(result.pathfinder_results) == nruns
            @test result.psis_result isa PSIS.PSISResult
            μ_hat = mean(result.draws; dims=2)
            Σ_hat = cov(result.draws .- μ_hat; dims=2, corrected=false)
            # adapted from the MvNormal tests
            # allow for 15x disagreement in atol, since this method is approximate
            multiplier = 15
            for i in eachindex(μ)
                atol = sqrt(Σ[i, i] / ndraws) * 8 * multiplier
                @test μ_hat[i] ≈ μ[i] atol = atol
            end
            for i in axes(Σ, 1), j in axes(Σ, 2)
                atol = sqrt(Σ[i, i] * Σ[j, j] / ndraws) * 10 * multiplier
                @test isapprox(Σ_hat[i, j], Σ[i, j], atol=atol)
            end

            Random.seed!(rng, seed)
            result2 = multipathfinder(
                logp, ndraws; dim, nruns, ndraws_elbo, ndraws_per_run, rng, executor
            )
            @test result2.fit_distribution == result.fit_distribution
            @test result2.draws == result.draws
            @test result2.draw_component_ids == result.draw_component_ids

            Random.seed!(rng, seed)
            ad_backend = AD.ReverseDiffBackend()
            result3 = multipathfinder(
                logp,
                ndraws;
                dim,
                nruns,
                ndraws_elbo,
                ndraws_per_run,
                rng,
                executor,
                ad_backend,
            )
            for (c1, c2) in
                zip(result.fit_distribution.components, result3.fit_distribution.components)
                @test c1 ≈ c2 atol = 1e-6
            end
        end

        init = [randn(dim) for _ in 1:nruns]
        result = multipathfinder(logp, ∇logp, ndraws; init)
        @test ncomponents(result.fit_distribution) == nruns
        @test size(result.draws) == (dim, ndraws)
    end
    @testset "errors if no gradient provided" begin
        logp(x) = -sum(abs2, x) / 2
        init = [randn(5) for _ in 1:10]
        fun = SciMLBase.OptimizationFunction(logp, GalacticOptim.AutoForwardDiff())
        @test_throws ArgumentError multipathfinder(fun, 10; init)
    end
    @testset "errors if neither dim nor init valid" begin
        logp(x) = -sum(abs2, x) / 2
        nruns = 2
        @test_throws ArgumentError multipathfinder(logp, 10; nruns)
        @test_throws ArgumentError multipathfinder(logp, 10; nruns, dim=0)
        multipathfinder(logp, 10; nruns, dim=2)
        multipathfinder(logp, 10; init=[randn(2) for _ in 1:nruns])
    end
    @testset "errors if neither init nor nruns valid" begin
        logp(x) = -sum(abs2, x) / 2
        @test_throws ArgumentError multipathfinder(logp, 10; dim=5, nruns=0)
        multipathfinder(logp, 10; dim=5, nruns=2)
    end
end


text:
#############################################################
#
# Search Functions
#
#############################################################

function get_search_tweets(q::AbstractString; options=Dict{AbstractString, AbstractString}())

    #Add required parameter(s) to options dict
    options["q"] = q

    r = get_oauth("https://api.twitter.com/1.1/search/tweets.json", options)

    #This is a placeholder, contains cursor data/metadata returned for paging
    #Evaluate returning a custom type where one field parsed array of tweets, one field metadata
    #metadata = r["search_metadata"]

    #return to_TWEETS(r["statuses"])

    r.status == 200 ?  to_TWEETS(Requests.json(r)["statuses"]) : r

end


text:
# using Pkg
# pkg"add Oceananigans CairoMakie JLD2"
# pushfirst!(LOAD_PATH, @__DIR__)
# pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..", "..")) # add Oceananigans

using Printf
using Statistics
using JLD2

using Oceananigans
using Oceananigans.Units
using Oceananigans.OutputReaders: FieldTimeSeries
using Oceananigans.Grids: xnode, ynode, znode
using Oceananigans.TurbulenceClosures: VerticallyImplicitTimeDiscretization
using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: CATKEVerticalDiffusivity

using Random
Random.seed!(1234)

arch = GPU()

filename = "eddying_channel_convadj"

# Domain
const Lx = 2000kilometers # zonal domain length [m]
const Ly = 2000kilometers # meridional domain length [m]
const Lz = 2kilometers    # depth [m]

# number of grid points
Nx = 128
Ny = 128
Nz = 60

save_fields_interval = 7days
stop_time = 20years + 1day
Δt₀ = 5minutes

# stretched grid

# we implement here a linearly streched grid in which the top grid cell has Δzₜₒₚ
# and every other cell is bigger by a factor σ, e.g.,
# Δzₜₒₚ, Δzₜₒₚ * σ, Δzₜₒₚ * σ², ..., Δzₜₒₚ * σᴺᶻ⁻¹,
# so that the sum of all cell heights is Lz

# Given Lz and stretching factor σ > 1 the top cell height is Δzₜₒₚ = Lz * (σ - 1) / σ^(Nz - 1)

σ = 1.04 # linear stretching factor
Δz_center_linear(k) = Lz * (σ - 1) * σ^(Nz - k) / (σ^Nz - 1) # k=1 is the bottom-most cell, k=Nz is the top cell
linearly_spaced_faces(k) = k==1 ? -Lz : - Lz + sum(Δz_center_linear.(1:k-1))

refinement = 2 # controls spacing near surface (higher means finer spaced)
stretching = 4  # controls rate of stretching at bottom

# Normalized height ranging from 0 to 1
h(k) = (k - 1) / Nz

# Linear near-surface generator
ζ₀(k) = 1 + (h(k) - 1) / refinement

# Bottom-intensified stretching function
Σ(k) = (1 - exp(-stretching * h(k))) / (1 - exp(-stretching))

# Generating function
z_faces(k) = Lz * (ζ₀(k) * Σ(k) - 1)

grid = RectilinearGrid(arch;
                       topology = (Periodic, Bounded, Bounded),
                       size = (Nx, Ny, Nz),
                       halo = (3, 3, 3),
                       x = (0, Lx),
                       y = (0, Ly),
                       z = (-Lz, 0)) # z_faces)

# The vertical spacing versus depth for the prescribed grid
# using GLMakie
# plot(grid.Δzᵃᵃᶜ[1:Nz], grid.zᵃᵃᶜ[1:Nz], marker = :circle,
#      axis=(xlabel = "Vertical spacing (m)",
#            ylabel = "Depth (m)"))

@info "Built a grid: $grid."

#####
##### Boundary conditions
#####

α  = 2e-4     # [K⁻¹] thermal expansion coefficient
g  = 9.8061   # [m s⁻²] gravitational constant
cᵖ = 3994.0   # [J K⁻¹] heat capacity
ρ  = 1024.0   # [kg m⁻³] reference density

parameters = (Ly = Ly,
              Lz = Lz,
              Qᵇ = 10 / (ρ * cᵖ) * α * g,          # buoyancy flux magnitude [m² s⁻³]
              y_shutoff = 5/6 * Ly,                # shutoff location for buoyancy flux [m]
              τ = 0.15/ρ,                          # surface kinematic wind stress [m² s⁻²]
              μ = 1 / 30days,                      # bottom drag damping time-scale [s⁻¹]
              ΔB = 8 * α * g,                      # surface vertical buoyancy gradient [s⁻²]
              H = Lz,                              # domain depth [m]
              h = 1000.0,                          # exponential decay scale of stable stratification [m]
              y_sponge = 19/20 * Ly,               # southern boundary of sponge layer [m]
              λt = 7days                           # relaxation time scale [s]
)

@inline function buoyancy_flux(i, j, grid, clock, model_fields, p)
    y = ynode(Center(), j, grid)
    return ifelse(y < p.y_shutoff, p.Qᵇ * cos(3π * y / p.Ly), 0.0)
end

buoyancy_flux_bc = FluxBoundaryCondition(buoyancy_flux, discrete_form=true, parameters=parameters)

@inline function u_stress(i, j, grid, clock, model_fields, p)
    y = ynode(Center(), j, grid)
    return - p.τ * sin(π * y / p.Ly)
end

u_stress_bc = FluxBoundaryCondition(u_stress, discrete_form=true, parameters=parameters)

@inline u_drag(i, j, grid, clock, model_fields, p) = @inbounds - p.μ * p.Lz * model_fields.u[i, j, 1]
@inline v_drag(i, j, grid, clock, model_fields, p) = @inbounds - p.μ * p.Lz * model_fields.v[i, j, 1]

u_drag_bc = FluxBoundaryCondition(u_drag, discrete_form=true, parameters=parameters)
v_drag_bc = FluxBoundaryCondition(v_drag, discrete_form=true, parameters=parameters)

b_bcs = FieldBoundaryConditions(top = buoyancy_flux_bc)

u_bcs = FieldBoundaryConditions(top = u_stress_bc, bottom = u_drag_bc)
v_bcs = FieldBoundaryConditions(bottom = v_drag_bc)

#####
##### Coriolis
#####

const f = -1e-4     # [s⁻¹]
const β =  1e-11    # [m⁻¹ s⁻¹]
coriolis = BetaPlane(f₀ = f, β = β)

#####
##### Forcing and initial condition
#####

@inline initial_buoyancy(z, p) = p.ΔB * (exp(z / p.h) - exp(-p.Lz / p.h)) / (1 - exp(-p.Lz / p.h))
@inline mask(y, p) = max(0.0, y - p.y_sponge) / (Ly - p.y_sponge)

@inline function buoyancy_relaxation(i, j, k, grid, clock, model_fields, p)
    timescale = p.λt
    y = ynode(Center(), j, grid)
    z = znode(Center(), k, grid)
    target_b = initial_buoyancy(z, p)
    b = @inbounds model_fields.b[i, j, k]
    return - 1 / timescale  * mask(y, p) * (b - target_b)
end

Fb = Forcing(buoyancy_relaxation, discrete_form = true, parameters = parameters)

# Turbulence closures

κh = 0.5e-5 # [m²/s] horizontal diffusivity
νh = 30.0   # [m²/s] horizontal viscocity
κz = 0.5e-5 # [m²/s] vertical diffusivity
νz = 3e-4   # [m²/s] vertical viscocity

diffusive_closure = AnisotropicDiffusivity(νh = νh,
                                           νz = νz,
                                           κh = κh,
                                           κz = κz,
					                       time_discretization = VerticallyImplicitTimeDiscretization())

convective_adjustment = ConvectiveAdjustmentVerticalDiffusivity(convective_κz = 1.0,
                                                                convective_νz = 0.0)

catke = CATKEVerticalDiffusivity()

#####
##### Model building
#####

@info "Building a model..."

model = HydrostaticFreeSurfaceModel(grid = grid,
                                    free_surface = ImplicitFreeSurface(),
                                    momentum_advection = WENO5(),
                                    tracer_advection = WENO5(),
                                    buoyancy = BuoyancyTracer(),
                                    coriolis = coriolis,
                                    closure = (convective_adjustment, diffusive_closure),
                                    tracers = (:b, :c),
                                    boundary_conditions = (b=b_bcs, u=u_bcs, v=v_bcs),
                                    forcing = (; b=Fb))

@info "Built $model."

#####
##### Initial conditions
#####

# resting initial condition
ε(σ) = σ * randn()
bᵢ(x, y, z) = parameters.ΔB * ( exp(z / parameters.h) - exp(-Lz / parameters.h) ) / (1 - exp(-Lz / parameters.h)) + ε(1e-8)
uᵢ(x, y, z) = ε(1e-8)
vᵢ(x, y, z) = ε(1e-8)
wᵢ(x, y, z) = ε(1e-8)

Δy = 100kilometers
Δz = 100
Δc = 2Δy
cᵢ(x, y, z) = exp(-(y - Ly/2)^2 / 2Δc^2) * exp(-(z + Lz/4)^2 / 2Δz^2)

set!(model, b=bᵢ, u=uᵢ, v=vᵢ, w=wᵢ, c=cᵢ)

#####
##### Simulation building
#####

simulation = Simulation(model, Δt=Δt₀, stop_time=stop_time)

# add timestep wizard callback
wizard = TimeStepWizard(cfl=0.1, max_change=1.1, max_Δt=20minutes)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(20))

# add progress callback
wall_clock = [time_ns()]

function print_progress(sim)
    @printf("[%05.2f%%] i: %d, t: %s, wall time: %s, max(u): (%6.3e, %6.3e, %6.3e) m/s, next Δt: %s\n",
            100 * (sim.model.clock.time / sim.stop_time),
            sim.model.clock.iteration,
            prettytime(sim.model.clock.time),
            prettytime(1e-9 * (time_ns() - wall_clock[1])),
            maximum(abs, sim.model.velocities.u),
            maximum(abs, sim.model.velocities.v),
            maximum(abs, sim.model.velocities.w),
            prettytime(sim.Δt))

    wall_clock[1] = time_ns()

    return nothing
end

simulation.callbacks[:print_progress] = Callback(print_progress, IterationInterval(20))

#####
##### Diagnostics
#####

u, v, w = model.velocities
b, c = model.tracers.b, model.tracers.c
η = model.free_surface.η

ζ = Field(∂x(v) - ∂y(u))

B = Field(Average(b, dims=1))
C = Field(Average(c, dims=1))
U = Field(Average(u, dims=1))
η̄ = Field(Average(η, dims=1))
V = Field(Average(v, dims=1))
W = Field(Average(w, dims=1))

b′ = b - B
v′ = v - V
w′ = w - W

b′b′ = Field(Average(b′ * b′, dims=1))
v′b′ = Field(Average(v′ * b′, dims=1))
w′b′ = Field(Average(w′ * b′, dims=1))

outputs = (; b, c, ζ, u, v, w)

zonally_averaged_outputs = (b=B, u=U, v=V, w=W, c=C, η=η̄, vb=v′b′, wb=w′b′, bb=b′b′)

#####
##### Build checkpointer and output writer
#####

simulation.output_writers[:checkpointer] = Checkpointer(model,
                                                        schedule = TimeInterval(10years),
                                                        prefix = filename,
                                                        force = true)

slicers = (west = FieldSlicer(i=1),
           east = FieldSlicer(i=grid.Nx),
           south = FieldSlicer(j=1),
           north = FieldSlicer(j=grid.Ny),
           bottom = FieldSlicer(k=1),
           top = FieldSlicer(k=grid.Nz))

for side in keys(slicers)
    field_slicer = slicers[side]

    simulation.output_writers[side] = JLD2OutputWriter(model, outputs,
                                                       schedule = TimeInterval(save_fields_interval),
                                                       field_slicer = field_slicer,
                                                       prefix = filename * "_$(side)_slice",
                                                       force = true)
end

simulation.output_writers[:zonal] = JLD2OutputWriter(model, zonally_averaged_outputs,
                                                     schedule = TimeInterval(save_fields_interval),
                                                     prefix = filename * "_zonal_average",
                                                     force = true)

#=
simulation.output_writers[:zonal] = JLD2OutputWriter(model, zonally_averaged_outputs,
                                                     schedule = AveragedTimeInterval(60days),
                                                     prefix = filename * "_zonal_time_average",
                                                     force = true)
=#

#=
simulation.output_writers[:averages] = JLD2OutputWriter(model, averaged_outputs,
                                                        schedule = AveragedTimeInterval(1days, window=1days, stride=1),
                                                        prefix = "eddying_channel_averages",
                                                        verbose = true,
                                                        force = true)
=#

@info "Running the simulation..."

run!(simulation, pickup=false)

# simulation.stop_time += 61days
# run!(simulation, pickup=true)

@info "Simulation completed in " * prettytime(simulation.run_wall_time)

#=
#####
##### Visualization
#####

# ENV["GKSwstype"] = "100"
# using CairoMakie

using GLMakie
using JLD2

fig = Figure(resolution = (2300, 1400))
ax_b = fig[1:5, 1] = LScene(fig)
ax_ζ = fig[1:5, 2] = LScene(fig)

axis_rotation_angles = (π/24, -π/6, 0)

# Extract surfaces on all 6 boundaries

iter = Observable(0)

filename = "eddying_channel_convadj"

zonal_file = jldopen(filename * "_zonal_average.jld2")

slicers = (west = FieldSlicer(i=1),
           east = FieldSlicer(i=grid.Nx),
           south = FieldSlicer(j=1),
           north = FieldSlicer(j=grid.Ny),
           bottom = FieldSlicer(k=1),
           top = FieldSlicer(k=grid.Nz))

sides = keys(slicers)

slice_files = NamedTuple(side => jldopen(filename * "_$(side)_slice.jld2") for side in sides)

# Build coordinates, rescaling the vertical coordinate

xζ, yζ, zζ = nodes((Face, Face, Center), grid)
xu, yu, zu = nodes((Face, Center, Center), grid)
xv, yv, zv = nodes((Center, Face, Center), grid)
xw, yw, zw = nodes((Center, Center, Face), grid)
xb, yb, zb = nodes((Center, Center, Center), grid)

zscale = 300
zu = zu .* zscale
zb = zb .* zscale
zζ = zζ .* zscale

zonal_slice_displacement = 1.5

b_slices = (
      west = @lift(Array(slice_files.west["timeseries/b/"   * string($iter)][1, :, :])),
      east = @lift(Array(slice_files.east["timeseries/b/"   * string($iter)][1, :, :])),
     south = @lift(Array(slice_files.south["timeseries/b/"  * string($iter)][:, 1, :])),
     north = @lift(Array(slice_files.north["timeseries/b/"  * string($iter)][:, 1, :])),
    bottom = @lift(Array(slice_files.bottom["timeseries/b/" * string($iter)][:, :, 1])),
       top = @lift(Array(slice_files.top["timeseries/b/"    * string($iter)][:, :, 1]))
)

clims_b = @lift extrema(slice_files.top["timeseries/b/" * string($iter)][:])
kwargs_b = (colorrange=clims_b, colormap=:deep, show_axis=false)

surface!(ax_b, yb, zb, b_slices.west;   transformation = (:yz, xb[1]),   kwargs_b...)
surface!(ax_b, yb, zb, b_slices.east;   transformation = (:yz, xb[end]), kwargs_b...)
surface!(ax_b, xb, zb, b_slices.south;  transformation = (:xz, yb[1]),   kwargs_b...)
surface!(ax_b, xb, zb, b_slices.north;  transformation = (:xz, yb[end]), kwargs_b...)
surface!(ax_b, xb, yb, b_slices.bottom; transformation = (:xy, zb[1]),   kwargs_b...)
surface!(ax_b, xb, yb, b_slices.top;    transformation = (:xy, zb[end]), kwargs_b...)

b_avg = @lift zonal_file["timeseries/b/" * string($iter)][1, :, :]
u_avg = @lift zonal_file["timeseries/u/" * string($iter)][1, :, :]

clims_u = @lift extrema(zonal_file["timeseries/u/" * string($iter)][1, :, :])

surface!(ax_b, yu, zu, u_avg; transformation = (:yz, zonal_slice_displacement * xu[end]), colorrange=clims_u, colormap=:balance)
contour!(ax_b, yb, zb, b_avg; levels = 25, color = :black, linewidth = 2, transformation = (:yz, zonal_slice_displacement * xb[end]), show_axis=false)

rotate_cam!(ax_b.scene, axis_rotation_angles)

ζ_slices = (
      west = @lift(Array(slice_files.west["timeseries/ζ/"   * string($iter)][1, :, :])),
      east = @lift(Array(slice_files.east["timeseries/ζ/"   * string($iter)][1, :, :])),
     south = @lift(Array(slice_files.south["timeseries/ζ/"  * string($iter)][:, 1, :])),
     north = @lift(Array(slice_files.north["timeseries/ζ/"  * string($iter)][:, 1, :])),
    bottom = @lift(Array(slice_files.bottom["timeseries/ζ/" * string($iter)][:, :, 1])),
       top = @lift(Array(slice_files.top["timeseries/ζ/"    * string($iter)][:, :, 1]))
)

u_slices = (
      west = @lift(Array(slice_files.west["timeseries/u/"   * string($iter)][1, :, :])),
      east = @lift(Array(slice_files.east["timeseries/u/"   * string($iter)][1, :, :])),
     south = @lift(Array(slice_files.south["timeseries/u/"  * string($iter)][:, 1, :])),
     north = @lift(Array(slice_files.north["timeseries/u/"  * string($iter)][:, 1, :])),
    bottom = @lift(Array(slice_files.bottom["timeseries/u/" * string($iter)][:, :, 1])),
       top = @lift(Array(slice_files.top["timeseries/u/"    * string($iter)][:, :, 1]))
)

# clims_ζ = @lift extrema(slice_files.west["timeseries/ζ/" * string($iter)][:])
clims_ζ = (-6e-5, 6e-5)
kwargs_ζ = (colorrange = clims_ζ, colormap=:curl, show_axis=false)

surface!(ax_ζ, yζ, zζ, ζ_slices.west;   transformation = (:yz, xζ[1]),   kwargs_ζ...)
surface!(ax_ζ, yζ, zζ, ζ_slices.east;   transformation = (:yz, xζ[end]), kwargs_ζ...)
surface!(ax_ζ, xζ, zζ, ζ_slices.south;  transformation = (:xz, yζ[1]),   kwargs_ζ...)
surface!(ax_ζ, xζ, zζ, ζ_slices.north;  transformation = (:xz, yζ[end]), kwargs_ζ...)
surface!(ax_ζ, xζ, yζ, ζ_slices.bottom; transformation = (:xy, zζ[1]),   kwargs_ζ...)
surface!(ax_ζ, xζ, yζ, ζ_slices.top;    transformation = (:xy, zζ[end]), kwargs_ζ...)

b_avg = @lift zonal_file["timeseries/b/" * string($iter)][1, :, :]
u_avg = @lift zonal_file["timeseries/u/" * string($iter)][1, :, :]

clims_u = @lift extrema(zonal_file["timeseries/u/" * string($iter)][1, :, :])

surface!(ax_ζ, yu, zu, u_avg; transformation = (:yz, zonal_slice_displacement * xu[end]), colorrange=clims_u, colormap=:balance)
contour!(ax_ζ, yb, zb, b_avg; levels = 25, color = :black, linewidth = 2, transformation = (:yz, zonal_slice_displacement * xb[end]), show_axis=false)

rotate_cam!(ax_ζ.scene, axis_rotation_angles)

title = @lift(string("Buoyancy and relative vorticity at t = ",
                     prettytime(zonal_file["timeseries/t/" * string($iter)])))

fig[0, :] = Label(fig, title, textsize=50)

iterations = parse.(Int, keys(zonal_file["timeseries/t"]))

record(fig, "eddying_channel_convadj.mp4", iterations, framerate=12) do i
    @info "Plotting iteration $i of $(iterations[end])..."
    iter[] = i
end

for file in slice_files
    close(file)
end

close(zonal_file)
=#


text:
using RadeonProRender, GeometryBasics, Colors
const RPR = RadeonProRender

function translationmatrix(t::Vec{3,T}) where {T}
    T0, T1 = zero(T), one(T)
    return Mat{4}(T1, T0, T0, T0, T0, T1, T0, T0, T0, T0, T1, T0, t[1], t[2], t[3], T1)
end

context = RPR.Context()
scene = RPR.Scene(context)
matsys = RPR.MaterialSystem(context, 0)

camera = RPR.Camera(context)
lookat!(camera, Vec3f(8, 0, 5), Vec3f(2, 0, 2), Vec3f(0, 0, 1))
RPR.rprCameraSetFocalLength(camera, 45.0)
set!(scene, camera)
set!(context, scene)

env_light = RadeonProRender.EnvironmentLight(context)
image_path = joinpath(@__DIR__, "studio026.exr")
img = RPR.Image(context, image_path)
set!(env_light, img)
setintensityscale!(env_light, 1.1)
push!(scene, env_light)

light = RPR.PointLight(context)
transform!(light, translationmatrix(Vec3f(2, 0, 5)))
f = 20
RPR.setradiantpower!(light, 500 / f, 641 / f, 630 / f)
push!(scene, light)

function add_shape!(scene, context, matsys, mesh; material=RPR.RPR_MATERIAL_NODE_DIFFUSE,
                    color=colorant"green", roughness=0.01)
    rpr_mesh = RadeonProRender.Shape(context, mesh)
    push!(scene, rpr_mesh)
    m = RPR.MaterialNode(matsys, material)
    set!(m, RPR.RPR_MATERIAL_INPUT_COLOR, color)
    set!(m, RPR.RPR_MATERIAL_INPUT_ROUGHNESS, roughness, roughness, roughness, roughness)
    set!(rpr_mesh, m)
    return rpr_mesh, m
end

mesh, mat = add_shape!(scene, context, matsys, Rect3f(Vec3f(-10, -10, -1), Vec3f(20, 20, 1));
                       color=colorant"white")
mesh, mat = add_shape!(scene, context, matsys, Rect3f(Vec3f(0, -10, 0), Vec3f(0.1, 20, 5));
                       color=colorant"white")
mesh, mat = add_shape!(scene, context, matsys, Rect3f(Vec3f(0, -2, 0), Vec3f(5, 0.1, 5));
                       color=colorant"white")
mesh, mat = add_shape!(scene, context, matsys, Rect3f(Vec3f(0, 2, 0), Vec3f(5, 0.1, 5));
                       color=colorant"white")

mesh, mat_sphere = add_shape!(scene, context, matsys, Tesselation(Sphere(Point3f(2, 0, 2), 1.0f0), 100);
                              material=RPR.RPR_MATERIAL_NODE_MICROFACET, roughness=0.2, color=colorant"red")

fb_size = (800, 600)
frame_buffer = RPR.FrameBuffer(context, RGBA, fb_size)
frame_bufferSolved = RPR.FrameBuffer(context, RGBA, fb_size)
set!(context, RPR.RPR_AOV_COLOR, frame_buffer)

set_standard_tonemapping!(context)

begin
    clear!(frame_buffer)
    RPR.rprContextSetParameterByKey1u(context, RPR.RPR_CONTEXT_ITERATIONS, 100)
    RPR.render(context)
    RPR.rprContextResolveFrameBuffer(context, frame_buffer, frame_bufferSolved, false)
    RPR.save(frame_bufferSolved, "test.png")
end


text:
function test_options()
  @testset "----- testing options -----" begin

    PetscOptionsSetValue("-testkey", "42")
    PetscOptionsView()
    PetscOptionsClearValue("-testkey")
  end  # end facts
end

test_options()


text:
function heatmap2glvisualize(attributes)
    result = Dict{Symbol, Any}()
    result[:stroke_width] = to_signal(attributes[:linewidth])
    result[:levels] = to_signal(attributes[:levels])
    result[:color_norm] = to_signal(attributes[:colornorm])
    result[:color_map] = to_signal(attributes[:colormap])

    heatmap = to_signal(lift_node(attributes[:heatmap]) do z
        [GLVisualize.Intensity{Float32}(z[j, i]) for i = 1:size(z, 2), j = 1:size(z, 1)]
    end)
    tex = GLAbstraction.Texture(value(heatmap), minfilter = :nearest)
    foreach(heatmap) do x
        update!(tex, x)
    end
    result, tex
end

function heatmap(scene::makie, x, y, z, attributes)
    attributes[:x] = x
    attributes[:y] = y
    attributes[:heatmap] = z
    attributes = heatmap_defaults(scene, attributes)
    gl_data, tex = heatmap2glvisualize(attributes)
    viz = visualize(tex, Style(:default), gl_data).children[]
    insert_scene!(scene, :heatmap, viz, attributes)
end

function image2glvisualize(attributes)
    result = Dict{Symbol, Any}()
    result[:primitive] = to_signal(lift_node(getindex.(attributes, (:x, :y, :spatialorder))...) do x, y, so
        if so != "xy"
            y, x = x, y
        end
        xmin, ymin = minimum(x), minimum(y)
        xmax, ymax = maximum(x), maximum(y)
        SimpleRectangle{Float32}(xmin, ymin, xmax - xmin, ymax - ymin)
    end)
    result[:spatialorder] = to_value(attributes[:spatialorder])
    result
end

function image(scene::makie, x, y, img, attributes::Dict)
    attributes[:x] = x
    attributes[:y] = y
    attributes[:image] = img
    attributes = image_defaults(scene, attributes)
    gl_data = image2glvisualize(attributes)
    viz = visualize(to_signal(attributes[:image]), Style(:default), gl_data).children[]
    insert_scene!(scene, :image, viz, attributes)
end

function volume2glvisualize(attributes)
    result = Dict{Symbol, Any}()
    if haskey(attributes, :colornorm)
        result[:color_norm] = to_signal(attributes[:colornorm])
        result[:color_map] = to_signal(attributes[:colormap])
    else
        result[:color_map] = nothing
        result[:color_norm] = nothing
        result[:color] = to_signal(attributes[:color])
    end

    result[:algorithm] = to_signal(attributes[:algorithm])
    result[:isovalue] = to_signal(attributes[:isovalue])
    result[:isorange] = to_signal(attributes[:isorange])
    result[:absorption] = to_signal(attributes[:absorption])

    result
end

function volume(scene::makie, values, attributes)
    attributes[:volume] = values
    attributes = volume_defaults(scene, attributes)
    gl_data = volume2glvisualize(attributes)
    viz = visualize(to_signal(attributes[:volume]), Style(:default), gl_data).children[]
    insert_scene!(scene, :volume, viz, attributes)
end


text:
export simplify
using Reduce
"""
    function simplify(FOS::AbstractFOS,Bra::AbstractBRA,Ket::AbstractBRA)
"""
function simplify(fos::AbstractFOS,bra::AbstractBra,ket::AbstractKet)
    checkpair(bra,ket)
    redundant = []
    coef = Basic(0) # suitable for dispatch
    delta = []

    pow_list = [fos[i].pow for i=1:length(fos)]
    while length(fos)!=0
        if vanish(bra,first(fos))
            absorb!(delta,fos)
        elseif vanish(ket,first(fos))
            coef = move_right!(coef,delta,fos)
        else
            extract_redundant!(fos,redundant)
        end
    end

    for (pow1,pow2) in delta
        coef = subs(coef, pow2=>pow1)
    end
    for pow in pow_list
        coef = subs(coef, pow^2=>pow)
    end
    for pow in pow_list
        coef = subs(coef, 2*pow=>0)
    end
    expr = SymEngine.toString(coef)

    # This is a foolish implementation, the aim is eliminating even term in expr

# Such is an option (guildline?) to develop this foolish implementation
# =================================
#   using Metatheory
#   using Metatheory.Library
#   using Metatheory.EGraphs

#   my_t = @theory begin
#       x * x => x
#       ex + ex => 0
#   end
#   t_add = commutative_group(:+, 0, :-)
#   t_mul = commutative_monoid(:*, 1)
#   t = my_t ∪ t_add ∪ t_mul

#   ex = :(b1*b1+b1*b2+b1*b2+b1*i1+b1*s1+b1*j1+b1*b2+b1*b2+b1*i1+b1*s1+b2*b2+b2*s1+b2*s2)

#   g = EGraph(ex)
#   saturate!(g, t)
#   ex = extract!(g, astsize)
#   # :(s2 * b2 + (s1 * b2 + (b2 + (b1 + j1 * b1))))
# =================================

    tol = 100
    for i = 2:2:tol
        expr = replace(expr,"+ $i*"=>"+ 0*")
        expr = replace(expr,"+ $(i+1)*"=>"+ 1*")
    end
    expr = RExpr(expr)
    # ---- Foolish part ends

    return (expr,redundant,delta) # +0 keeps the -1^0 safe.
end

"""
    function extract!(fos::AbstractFOS,redundant::Array)
        return nothing

    Extract the operator at the left/first of Input fos into result. Because it can not be vanished by Bra and Ket
"""
function extract_redundant!(fos::AbstractFOS,redundant::Array)
    push!(redundant,first(fos.fos))
    popfirst!(fos.fos)
end

"""
    Absorb the left Operator in fos into Bra. Will produce a coefficient into result.
"""
function absorb!(delta::Array,fos::AbstractFOS)
    @warn "delta function looks like delta(n,0) is produced because the operator"
    # TODO: implement this function
end

"""
    function move_right!(coef::Dict,fos::AbstractFOS)
"""
function move_right!(coef::Basic,delta::Array,fos::AbstractFOS)
    for i = 2:length(fos.fos)
        if fos[i].fo == dag(fos[1].fo)
            push!(delta,(fos[1].pow,fos[i].pow))
            deleteat!(fos.fos,i)
            deleteat!(fos.fos,1)
            return coef
        else
            coef = coef+fos[1].pow*fos[i].pow
        end
    end
    @warn "Symbol $(fos[1]) doesnot meet a pair!"
end


text:
using ContextTracking
using Test

@testset "ContextTracking.jl" begin
    include("test_context.jl")
    include("test_trace.jl")
end


text:
using VulkanGen
using Test
using StructArrays

@testset "VulkanGen.jl" begin
    include("spec.jl")
    include("wrapper.jl")
end


text:
"""
	write(io::IO, arg::ARG)
	write(filename::AbstractString, arg::ARG)

Write `arg` using the extended Newick format.
"""
write(filename::AbstractString, arg::ARG) = write(filename, extended_newick(arg))
write(io::IO, arg::ARG) = write(io, extended_newick(arg))

#=
Strategy
(i) If the global root is shared
Go recursively down this root and write a newick.
If a hybrid node is met for the first time, store it in a list of hybrid nodes met
If it is met for the second time, stop recursing down
(ii) none of the roots are shared
feed the function with each root in turn (and the list of hybrids)
(iii) One of the root is shared, not the other one
start with the non shared one. should work allright
=#
function extended_newick(arg::ARG)
	hybrids = Dict()

	if arg.roots[1] == arg.roots[2]
		# (i)
		str = extended_newick!(arg.roots[1], nothing, hybrids)
		return str * ";"
	else
		if !isshared(arg.roots[1]) && !isshared(arg.roots[2])
			#(ii)
			str1 = extended_newick!(arg.roots[1], nothing, hybrids)
			str2 = extended_newick!(arg.roots[2], nothing, hybrids)
			return "($(str1),$(str2))GlobalRoot:0.;"
		elseif isshared(arg.roots[1])
			#(iii)
			str = extended_newick!(arg.roots[2], nothing, hybrids)
			return str * ";"
		elseif isshared(arg.roots[2])
			str = extended_newick!(arg.roots[1], nothing, hybrids)
			return str * ";"
		end
	end
end

function extended_newick!(a_r::ARGNode, anc, hybrids::Dict)
	# println("$(label(a_r))")
	nwk = ""
	# Dealing with children of a_r
	if !isleaf(a_r) && !haskey(hybrids, label(a_r))
		nwk *= "("
		for c in children(a_r)
			nwk *= extended_newick!(c, a_r, hybrids)
			nwk *= ","
		end
		nwk = nwk[1:end-1] # Removing trailing ','
		nwk *= ")"
	end
	# Adding a_r to hybrid if necessary
	if ishybrid(a_r) && !haskey(hybrids, label(a_r))
		i = length(hybrids) +1
		hybrids[label(a_r)] = "#H$(i)"
	end
	# Writing a_r itself
	a_r_label = if ishybrid(a_r)
		label(a_r) * hybrids[label(a_r)]
	else
		label(a_r)
	end
	nwk *= a_r_label
	nwk *= nwk_node_data(a_r, anc)

	return nwk
end

"""
	nwk_node_data(n::ARGNode, a)

Return a string of the form "[&segments=0,1]:$(time)".

Different cases to choose from:
	(i) Node is shared, not a hybrid, not the root of only one tree
  its two times should be the same
  [&segments={0,1}]
	(ii) Node is not shared
  It should only have one color
  [&segments={c}]
	(iii) Node is shared and a hybrid
  We need to know what the color of the branch we arrived from is
  [&segments={c}]
	(iv) Node is the root of one tree and not of the other
  If `c` is the non-root color
  [&segments={c}]
"""
function nwk_node_data(n::ARGNode, a)
	if isshared(n) && !ishybrid(n) && !is_partial_root(n)
		# (i)
		τ1 = branch_length(n)[1]
		τ2 = branch_length(n)[2]
		@assert (ismissing(τ1) && ismissing(τ2)) || isapprox(τ1, τ2)
		return nwk_node_data(τ1, 1, 2)

	elseif isshared(n) && ishybrid(n)
		#(iii)
		ec = edgecolor(a, n)
		if ec[1] && !ec[2]
			return nwk_node_data(branch_length(n)[1], 1)
		elseif !ec[1] && ec[2]
			return nwk_node_data(branch_length(n)[2], 2)
		elseif ec[1] && ec[2]
			@assert isapprox(branch_length(n)[1], branch_length(n)[2])
			return nwk_node_data(branch_length(n)[1], 1, 2)
		end

	elseif isshared(n) && isroot(n)
		# (iv)
		for c in 1:2
			if !isroot(n, c)
				@assert isroot(n, othercolor(c))
				return nwk_node_data(branch_length(n)[c], c)
			end
		end
		@error "$n"
	elseif !isshared(n)
		# (ii)
		for c in 1:2
			if hascolor(n, c)
				@assert !hascolor(n, othercolor(c))
				return nwk_node_data(branch_length(n)[c], c)
			end
		end
		@error "No color for node $n"
	else
		@error "Unclear state for node $n"
	end
end

function nwk_node_data(τ, c::Int)
	ismissing(τ) ? "[&segments={$(c-1)}]" : "[&segments={$(c-1)}]:$(τ)"
end
function nwk_node_data(τ, c1::Int, c2::Int)
	clr = "$(min(c1,c2)-1),$(max(c1,c2)-1)"
	ismissing(τ) ? "[&segments={$(clr)}]" : "[&segments={$(clr)}]:$(τ)"
end

max_issues_repo_name: PierreBarrat/TreeKnit

max_forks_repo_name: PierreBarrat/TreeKnit


text:
abstract type Topology end

struct Periodic <: Topology end
struct NonPeriodic <: Topology end

periodic(args...) = Periodic()
nonperiodic(args...) = NonPeriodic()

isperiodic(::Topology) = false
isperiodic(::Periodic) = true


text:
using PatModules
@auxmodule! one ("utils",) begin
    export say_hi

    function say_hi()
        say("hi")
    end
end


text:
include("src/Scotch.jl")

import Scotch
using PyPlot

@show Scotch.version()

# Build matrix & graph
graph = Scotch.graphAlloc()
Scotch.graphInit(graph)
n = 100
function laplacian_2d(n)
    L = spdiagm((-ones(n-1), 4*ones(n), -ones(n-1)), (-1, 0, 1))
    return kron(speye(n),L) + kron(L,speye(n))
end
A = laplacian_2d(n)
N = n*n
Scotch.graphBuild(graph, A)
Scotch.graphCheck(graph)
(nv, ne) = Scotch.graphSize(graph)
@show nv, ne

# Create strategy
strat = Scotch.stratAlloc()
Scotch.stratInit(strat)
Scotch.stratGraphOrder(strat, "n{sep=/levl<2?g:z;}")

# Create ordering
p    = Array{Cint,1}(N)
ip   = Array{Cint,1}(N)
nb   = Array{Cint,1}(1)
cblk = Array{Cint,1}(N+1)
tree = Array{Cint,1}(N)

order = Scotch.orderAlloc()
Scotch.graphOrderInit(graph, order, p, ip, nb, cblk, tree)

Scotch.graphOrderComputeList(graph, order, Array{Cint,1}(1:N), strat)
Scotch.graphOrderCheck(graph, order)

z = zeros(Int64, n*n)
for i = 1:nb[1]
    z[ip[cblk[i]:cblk[i+1]-1]] = i
end
Z = reshape(z, (n, n));
matshow(Z)

# Repartition part of the graph
strat = Scotch.stratAlloc()
Scotch.stratInit(strat)
Scotch.stratGraphOrder(strat, "n{sep=/levl<4?g:z;}")
list = ip[cblk[1]:cblk[2]-1]

Scotch.graphOrderComputeList(graph, order, list, strat)
Scotch.graphOrderCheck(graph, order)

z = zeros(Int64, n*n)
for i = 1:nb[1]
    z[ip[cblk[i]:cblk[i+1]-1]] = i
end
Z = reshape(z, (n, n));
matshow(Z)


text:
import Base
using FillArrays: Fill
# """
# A power measure is a product of a measure with itself. The number of elements in
# the product determines the dimensionality of the resulting support.

# Note that power measures are only well-defined for integer powers.

# The nth power of a measure μ can be written μ^x.
# """
# PowerMeasure{M,N,D} = ProductMeasure{Fill{M,N,D}}

# function Base.show(io::IO, μ::PowerMeasure)
#     io = IOContext(io, :compact => true)
#     print(io, μ.data.value, " ^ ", size(μ.data))
# end

# function Base.show_unquoted(io::IO, μ::PowerMeasure{M,N,D}, indent::Int, prec::Int) where {M,N,D}
#     io = IOContext(io, :compact => true)
#     if Base.operator_precedence(:^) ≤ prec
#         print(io, "(")
#         show(io, μ.data.value)
#         print(io, ")")
#     else
#         show(io, size(μ.data))
#     end
#     return nothing
# end

export PowerMeasure

const PowerMeasure{F,S,T,N,A} = ProductMeasure{F,S,Fill{T,N,A}}

Base.:^(μ::AbstractMeasure, ::Tuple{}) = μ

function Base.:^(μ::AbstractMeasure, dims::Integer...)
    return μ^dims
end

function Base.:^(μ::M, dims::NTuple{N,I}) where {M<:AbstractMeasure,N,I<:Integer}
    powermeasure(μ, dims)
end

# Same as PowerMeasure
function Pretty.tile(d::ProductMeasure{Returns{T},I,C}) where {T,I,C<:CartesianIndices}
    Pretty.pair_layout(Pretty.tile(d.f.f.value), Pretty.tile(size(d.pars)); sep = " ^ ")
end

function Pretty.tile(d::ProductMeasure{R,I,V}) where {R<:Returns,I,V<:AbstractVector}
    Pretty.pair_layout(Pretty.tile(d.f.f.value), Pretty.tile(length(d.pars)); sep = " ^ ")
end

# sampletype(d::PowerMeasure{M,N}) where {M,N} = @inbounds Array{sampletype(first(marginals(d))), N}

params(d::ProductMeasure{F,S,<:Fill}) where {F,S} = params(first(marginals(d)))

params(::Type{P}) where {F,S,P<:ProductMeasure{F,S,<:Fill}} = params(D)

# basemeasure(μ::PowerMeasure) = @inbounds basemeasure(first(μ.data))^size(μ.data)

# Same as PowerMeasure
@inline function basemeasure(d::ProductMeasure{F,S,<:Fill}) where {F,S}
    _basemeasure(d, (basemeasure(d.f(first(d.pars)))))
end

# Same as PowerMeasure
@inline function _basemeasure(d::ProductMeasure{F,S,<:Fill}, b) where {F,S}
    b^size(d.pars)
end

# Same as PowerMeasure
@inline function _basemeasure(d::ProductMeasure{F,S,<:Fill}, b::FactoredBase) where {F,S}
    n = length(d.pars)
    inbounds(x) = all(b.inbounds, x)
    constℓ = n * b.constℓ
    varℓ() = n * b.varℓ()
    base = b.base^size(d.pars)
    FactoredBase(inbounds, constℓ, varℓ, base)
end

# Same as PowerMeasure
@inline function logdensity(d::ProductMeasure{F,S,<:Fill}, x) where {F,S}
    d1 = d.f(first(d.pars))
    sum(xj -> logdensity(d1, xj), x)
end


text:
using SemanticModels
using Test
using Base.Meta

@testset "Parsers" begin
@show expr = Parsers.parsefile("../examples/epicookbook/notebooks/SimpleDeterministicModels/SIRModel.jl")

modulename = expr.args[2]
mc = Parsers.defs(expr.args[3].args)

@show expr = Parsers.parsefile("../examples/epicookbook/notebooks/SimpleDeterministicModels/SEIRmodel.jl")

sericmodulename = expr.args[2]
seirc = Parsers.defs(expr.args[3].args)
@test seirc.modc[1][1] == :DifferentialEquations
@test seirc.fc.defs[1][1].head==:call
@show seirc.fc.defs[1][1].args[1] == :seir_ode
@show seirc.fc.defs[1][1].args == [:seir_ode, :dY, :Y, :p, :t]
@test map(first, seirc.vc) == [ :pram, :init, :tspan, :seir_prob, :sol, :va, :y, :R]
end


text:
using DSGE, JLD2, DataStructures
import Test: @test, @testset
import DSGE: @test_matrix_approx_eq

# Load pre-input files
out = load("../../test_outputs/solve/splines_outputs_mat.jld2")
g0 = full(out["g0"]); g1 = full(out["g1"]); psi = full(out["psi"])
Pi = full(out["pi"]); c = full(out["constant"])
out = load("../../test_outputs/solve/splines_outputs.jld2")
Γ0 = full(out["Γ0"]); Γ1 = full(out["Γ1"]); Ψ = full(out["Ψ"])
Π = full(out["Π"]); C_val = full(out["C"])

# Load checking files
out = load("../../test_outputs/solve/post_gensys_mat.jld2")
G1_mat = full(out["G1"]); impact_mat = full(out["impact"])
out = load("../../test_outputs/solve/post_gensys.jld2")
G1 = full(out["G1"]); impact = full(out["impact"]);

# Check non-in-place first
T_mat, ~, R_mat, ~, ~, eu_mat = gensysct(g1, c, psi, Pi, complex_decomposition = false)
T, ~, R, ~, ~, ~, ~, eu = gensysct(Γ0, Γ1, C_val, Ψ, Π, complex_decomposition = false)
@testset "Non-in-place gensysct" begin
    @test @test_matrix_approx_eq T_mat G1_mat
    @test @test_matrix_approx_eq R_mat impact_mat
    @test eu_mat == [1; 1]
    @test @test_matrix_approx_eq T G1
    @test @test_matrix_approx_eq R impact
    @test eu == [-1; 1]
end

# Check in place gensysct
T_mat, ~, R_mat, ~, ~, eu_mat = gensysct!(g1, c, psi, Pi, complex_decomposition = false)
T, ~, R, ~, ~, ~, ~, eu = gensysct!(Γ0, Γ1, C_val, Ψ, Π, complex_decomposition = false)
@testset "In-place gensysct" begin
    @test @test_matrix_approx_eq T_mat G1_mat
    @test @test_matrix_approx_eq R_mat impact_mat
    @test eu_mat == [1; 1]
    @test @test_matrix_approx_eq T G1
    @test @test_matrix_approx_eq R impact
    @test eu == [-1; 1]
end


text:
# =============================================================
#		MODULE: option
# =============================================================
module option
	Classification = "NZsimple" #[NZ] or [NZsimple]
end  # module option
# ............................................................

max_stars_repo_name: joseph-pollacco/TextureTriangleNZ

max_issues_repo_name: joseph-pollacco/TextureTriangleNZ

max_forks_repo_name: joseph-pollacco/TextureTriangleNZ


text:
using DataStructures

d = matchall.(r"(\w+)", readlines("input.txt"))
parent = Dict()
children = Dict()
weights = Dict()

for x in d
    if length(x) > 2
        for c in x[3:end]
            parent[c] = x[1]
        end
        children[x[1]] = x[3:end]
    end
    weights[x[1]] = parse(Int, x[2])
end

# Part 1
root = values(parent) |> first
while haskey(parent, root)
    root = parent[root]
end
println(root)

# Part 2
function calcWeight(r)
    return haskey(children, r) ? mapreduce(calcWeight, +, weights[r], children[r]) : weights[r]
end

# https://github.com/JuliaCollections/DataStructures.jl/issues/195
most_common(ct::Accumulator) = sort(collect(ct.map), by=p->p[2], rev=true)
target = root
targetWeight = 0
while true
    v = counter(map(calcWeight, children[target]))
    if length(v) == 1
        println(weights[target] + targetWeight - calcWeight(target))
        break
    end
    v = most_common(v)
    targetWeight = v[1][1]
    target = [c for c in children[target] if calcWeight(c) != targetWeight][1]
end


text: println(ARGS)

max_stars_repo_name: faicaltoubali/ENSEEIHT

max_issues_repo_name: faicaltoubali/ENSEEIHT

max_forks_repo_name: faicaltoubali/ENSEEIHT