text:
using Base.Test
using DataFrames

macro test_da_pda(da, code)
    esc(quote
        let $da = copy($da)
            $code
        end
        let $da = PooledDataArray($da)
            $code
        end
    end)
end

let
    unary_operators = [(+), (-), (!)]

    numeric_unary_operators = [(+), (-)]

    logical_unary_operators = [(!)]

    elementary_functions = [abs, sign, acos, acosh, asin,
                            asinh, atan, atanh, sin, sinh,
                            cos, cosh, tan, tanh, ceil, floor,
                            round, trunc, exp, exp2, expm1, log, log10, log1p,
                            log2, exponent, sqrt, gamma, lgamma, digamma,
                            erf, erfc]

    special_comparison_operators = [isless]

    comparison_operators = [(==), (.==), (!=), (.!=),
                            (>), (.>), (>=), (.>=), (<), (.<),
                            (<=), (.<=)]

    scalar_comparison_operators = [(==), (!=), (>), (>=),
                                   (<), (<=)]

    array_comparison_operators = [(.==), (.!=), (.>), (.>=), (.<), (.<=)]

    binary_operators = [(+), (.+), (-), (.-), (*), (.*), (/), (./),
                        (.^), (div), (mod), (fld), (rem),
                        (&), (|), ($)]

    induced_binary_operators = [(^)]

    arithmetic_operators = [(+), (.+), (-), (.-), (*), (.*), (/), (./),
                            (.^), (div), (mod), (fld), (rem)]

    induced_arithmetic_operators = [(^)]

    biscalar_operators = [(max), (min)]

    scalar_arithmetic_operators = [(+), (-), (*), (/),
                                   (div), (mod), (fld), (rem)]

    induced_scalar_arithmetic_operators = [(^)]

    array_arithmetic_operators = [(+), (.+), (-), (.-), (.*), (./), (.^)]

    bit_operators = [(&), (|), ($)]

    unary_vector_operators = [minimum, maximum, prod, sum, mean, median, std,
                              var, mad, norm, skewness, kurtosis]

    pairwise_vector_operators = [diff, percent_change]

    cumulative_vector_operators = [cumprod, cumsum, cumsum_kbn, cummin, cummax]

    ffts = [fft]

    binary_vector_operators = [dot, cor, cov, cor_spearman]

    columnar_operators = [colmins, colmaxs, colprods, colsums,
                          colmeans, colmedians, colstds, colvars,
                          colnorms]

    boolean_operators = [any, all]

    # All unary operators return NA when evaluating NA
    for f in unary_operators
        @assert isna(f(NA))
    end

    # All elementary functions return NA when evaluating NA
    for f in elementary_functions
        @assert isna(f(NA))
    end

    # All comparison operators return NA when comparing NA with NA
    # All comparison operators return NA when comparing scalars with NA
    # All comparison operators return NA when comparing NA with scalars
    for f in comparison_operators
        @assert isna(f(NA, NA))
        @assert isna(f(NA, 1))
        @assert isna(f(1, NA))
    end

    # All arithmetic operators return NA when operating on two NA's
    # All arithmetic operators return NA when operating on a scalar and an NA
    # All arithmetic operators return NA when operating on an NA and a scalar
    for f in arithmetic_operators
        @assert isna(f(NA, NA))
        @assert isna(f(1, NA))
        @assert isna(f(NA, 1))
    end

    # All bit operators return NA when operating on two NA's
    # All bit operators return NA when operating on a scalar and an NA
    # All bit operators return NA when operating on an NA and a scalar
    for f in bit_operators
        @assert isna(f(NA, NA))
        @assert isna(f(1, NA))
        @assert isna(f(NA, 1))
    end

    # Unary operators on DataVector's should be equivalent to elementwise
    # application of those same operators
    N = 5
    dv = dataones(N)
    @test_da_pda dv begin
        for f in numeric_unary_operators
            for i in 1:length(dv)
                @assert f(dv)[i] == f(dv[i])
            end
        end
    end
    dv = datatrues(N)
    @test_da_pda dv begin
        for f in logical_unary_operators
            for i in 1:length(dv)
                @assert f(dv)[i] == f(dv[i])
            end
        end
    end

    # Unary operators on DataFrame's should be equivalent to elementwise
    # application of those same operators
    df = DataFrame(quote
                       A = [1, 2, 3, 4]
                       B = [1.0, pi, pi, e]
                   end)
    for f in numeric_unary_operators
        for i in 1:nrow(df)
            for j in 1:ncol(df)
                @assert f(df)[i, j] == f(df[i, j])
            end
        end
    end
    df = DataFrame(quote
                       A = [true, false, true, false]
                   end)
    for f in logical_unary_operators
        for i in 1:nrow(df)
            for j in 1:ncol(df)
                @assert f(df)[i, j] == f(df[i, j])
            end
        end
    end

    # Elementary functions on DataVector's
    N = 5
    dv = dataones(N)
    @test_da_pda dv begin
        for f in elementary_functions
            for i in 1:length(dv)
                @assert f(dv)[i] == f(dv[i])
            end
        end
    end

    # Elementary functions on DataFrames's
    N = 5
    df = DataFrame(quote
                       A = dataones($(N))
                       B = dataones($(N))
                   end)
    for f in elementary_functions
        for i in 1:nrow(df)
            for j in 1:ncol(df)
                  @assert f(df)[i, j] == f(df[i, j])
            end
        end
    end

    # Broadcasting operations between NA's and DataVector's
    N = 5
    dv = dataones(N)
    @test_da_pda dv begin
        for f in arithmetic_operators
            for i in 1:length(dv)
                @assert isna(f(dv, NA)[i])
                @assert isna(f(NA, dv)[i])
            end
        end
    end

    # Broadcasting operations between NA's and DataVector's
    N = 5
    dv = dataones(N)
    @test_da_pda dv begin
        for f in arithmetic_operators
            for i in 1:length(dv)
                @assert isna(f(dv, NA)[i])
                @assert isna(f(NA, dv)[i])
            end
        end
    end

    # Broadcasting operations between scalars and DataVector's
    N = 5
    dv = dataones(N)
    @test_da_pda dv begin
        for f in arithmetic_operators
            for i in 1:length(dv)
                @assert f(dv, 1)[i] == f(dv[i], 1)
                @assert f(1, dv)[i] == f(1, dv[i])
            end
        end
    end
    dv = DataVector[false, true, false, true, false]
    for f in bit_operators
        for i in 1:length(dv)
            @assert f(dv, true)[i] == f(dv[i], true)
            @assert f(true, dv)[i] == f(true, dv[i])
        end
    end

    # Broadcasting operations between NA's and DataFrames's
    N = 5
    df = DataFrame(quote
                       A = dataones($(N))
                       B = dataones($(N))
                   end)
    for f in arithmetic_operators
        for i in 1:nrow(df)
            for j in 1:ncol(df)
                @assert isna(f(df, NA)[i, j])
                @assert isna(f(NA, df)[i, j])
            end
        end
    end

    # Broadcasting operations between scalars and DataFrames's
    N = 5
    df = DataFrame(quote
                       A = dataones($(N))
                       B = dataones($(N))
                   end)
    for f in arithmetic_operators
        for i in 1:nrow(df)
            for j in 1:ncol(df)
                @assert f(df, 1)[i, j] == f(df[i, j], 1)
                @assert f(1, df)[i, j] == f(1, df[i, j])
            end
        end
    end

    # Binary operations on (DataVector, Vector) or (Vector, DataVector)
    N = 5
    v = ones(N)
    dv = dataones(N)
    dv[1] = NA
    @test_da_pda dv begin
        for f in array_arithmetic_operators
            for i in 1:length(dv)
                @assert isna(f(v, dv)[i]) && isna(dv[i]) ||
                        f(v, dv)[i] == f(v[i], dv[i])
                @assert isna(f(dv, v)[i]) && isna(dv[i]) ||
                        f(dv, v)[i] == f(dv[i], v[i])
            end
        end
    end

    # Binary operations on pairs of DataVector's
    N = 5
    dv = dataones(N)
    dv[1] = NA
    @test_da_pda dv begin
        for f in array_arithmetic_operators
            for i in 1:length(dv)
                @assert isna(f(dv, dv)[i]) && isna(dv[i]) ||
                        f(dv, dv)[i] == f(dv[i], dv[i])
            end
        end
    end

    # Binary operations on pairs of DataFrame's
    # TODO: Test in the presence of in-operable types like Strings
    N = 5
    df = DataFrame(quote
                       A = dataones($(N))
                       B = dataones($(N))
                   end)
    for f in array_arithmetic_operators
        for i in 1:nrow(df)
            for j in 1:ncol(df)
                @assert isna(f(df, df)[i, j]) && isna(df[i, j]) ||
                        f(df, df)[i, j] == f(df[i, j], df[i, j])
            end
        end
    end

    # Unary vector operators on DataVector's
    N = 5
    dv = dataones(5)
    for f in unary_vector_operators
        @assert isequal(f(dv), f(dv.data))
    end
    dv[1] = NA
    for f in unary_vector_operators
        @assert isna(f(dv))
    end

    # Pairwise vector operators on DataVector's
    N = 5
    dv = DataVector[911, 269, 835.0, 448, 772]
    for f in pairwise_vector_operators
        @assert isequal(f(dv), f(dv.data))
    end
    dv = DataVector[NA, 269, 835.0, 448, 772]
    for f in pairwise_vector_operators
        v = f(dv)
        @assert isna(v[1])
        @assert isequal(v[2:4], f(dv.data)[2:4])
    end
    dv = DataVector[911, NA, 835.0, 448, 772]
    for f in pairwise_vector_operators
        v = f(dv)
        @assert isna(v[1])
        @assert isna(v[2])
        @assert isequal(v[3:4], f(dv.data)[3:4])
    end
    dv = DataVector[911, 269, 835.0, 448, NA]
    for f in pairwise_vector_operators
        v = f(dv)
        @assert isna(v[4])
        @assert isequal(v[1:3], f(dv.data)[1:3])
    end

    # Cumulative vector operators on DataVector's
    N = 5
    dv = dataones(N)
    for f in cumulative_vector_operators
        for i in 1:length(dv)
            @assert f(dv)[i] == f(dv.data)[i]
        end
    end
    dv[4] = NA
    for f in cumulative_vector_operators
        for i in 1:3
            @assert f(dv)[i] == f(dv.data)[i]
        end
        for i in 4:N
            @assert isna(f(dv)[i])
        end
    end

    # FFT's on DataVector's
    N = 5
    dv = dataones(5)
    for f in ffts
        @assert f(dv) == f(dv.data)
    end
    dv[1] = NA
    for f in ffts
        @assert isna(f(dv))
    end

    # Binary vector operators on DataVector's
    N = 5
    dv = dataones(5)
    for f in binary_vector_operators
        @assert f(dv, dv) == f(dv.data, dv.data) ||
                (isnan(f(dv, dv)) && isnan(f(dv.data, dv.data)))
    end
    dv[1] = NA
    for f in binary_vector_operators
        @assert isna(f(dv, dv))
    end

    # TODO: Columnar operators on DataFrame's

    # Boolean operators on DataVector's
    N = 5
    @assert any(datafalses(N)) == false
    @assert any(datatrues(N)) == true
    @assert all(datafalses(N)) == false
    @assert all(datatrues(N)) == true
    @assert any(PooledDataArray(datafalses(N))) == false
    @assert any(PooledDataArray(datatrues(N))) == true
    @assert all(PooledDataArray(datafalses(N))) == false
    @assert all(PooledDataArray(datatrues(N))) == true

    dv = datafalses(N)
    dv[3] = true
    @test_da_pda dv begin
        @assert any(dv) == true
        @assert all(dv) == false
    end

    dv = datafalses(N)
    dv[2] = NA
    dv[3] = true
    @test_da_pda dv begin
        @assert any(dv) == true
        @assert all(dv) == false
    end

    dv = datafalses(N)
    dv[2] = NA
    @test_da_pda dv begin
        @assert isna(any(dv))
        @assert all(dv) == false
    end

    dv = datafalses(1)
    dv[1] = NA
    @test_da_pda dv begin
        @assert isna(any(dv))
        @assert isna(all(dv))
    end

    # Boolean operators on DataFrames's
    N = 5
    df = DataFrame(quote
                       A = datafalses($(N))
                   end)
    @assert any(df) == false
    @assert any(!df) == true
    @assert all(df) == false
    @assert all(!df) == true

    df = DataFrame(quote
                       A = datafalses($(N))
                   end)
    df[3, 1] = true
    @assert any(df) == true
    @assert all(df) == false

    df = DataFrame(quote
                       A = datafalses($(N))
                   end)
    df[2, 1] = NA
    df[3, 1] = true
    @assert any(df) == true
    @assert all(df) == false

    df = DataFrame(quote
                       A = datafalses($(N))
                   end)
    df[2, 1] = NA
    @assert isna(any(df))
    @assert all(df) == false

    df = DataFrame(quote
                     A = datafalses($(N))
                   end)
    df[1, 1] = NA
    @assert isna(any(dv))
    @assert isna(all(dv))

    # Is this a genuine special case?
    @assert isna(NA ^ 2.0)

    #
    # Equality tests
    #

    dv = DataVector[1, NA]
    alt_dv = DataVector[2, NA]
    pdv = PooledDataArray(DataVector[1, NA])
    alt_pdv = PooledDataArray(DataVector[2, NA])
    df = DataFrame({dv})
    alt_df = DataFrame({alt_dv})

    @assert isna(NA == NA)
    @assert isna(NA != NA)

    # @assert isna(dv == dv) # SHOULD RAISE ERROR
    # @assert isna(dv != dv) # SHOULD RAISE ERROR

    # @assert isna(df == df) # SHOULD RAISE ERROR
    # @assert isna(df != df) # SHOULD RAISE ERROR

    @assert isequal(dv, dv)
    @assert isequal(pdv, pdv)
    @assert isequal(df, df)

    @assert !isequal(dv, alt_dv)
    @assert !isequal(pdv, alt_pdv)
    @assert !isequal(df, alt_df)

    @assert isequal(DataVector[1, NA] .== DataVector[1, NA], DataVector[true, NA])
    @assert isequal(PooledDataVector[1, NA] .== PooledDataVector[1, NA], DataVector[true, NA])
    @assert isequal(DataFrame({dv}) .== DataFrame({dv}), DataFrame({DataVector[true, NA]}))

    @assert all(isna(NA .== dataones(5)))
    @assert all(isna(dataones(5) .== NA))
    @assert all(isna(NA .== PooledDataArray(dataones(5))))
    @assert all(isna(PooledDataArray(dataones(5)) .== NA))

    @assert all(isna(NA .== df))
    @assert all(isna(df .== NA))

    # Run length encoding
    dv = dataones(5)
    dv[3] = NA

    v, l = rle(dv)
    @assert isequal(v, DataVector[1.0, NA, 1.0])
    @assert (l == [2, 1, 2])

    rdv = inverse_rle(v, l)
    @assert isequal(dv, rdv)
end


text:
using ConfigEnv
using AlphaVantage
using DataFrames, CSV
using Dates
using Indicators
using Setfield
using Plots
using BenchmarkTools

dotenv();

########################################
# Downloading data
########################################
df = digital_currency_daily("BTC", datatype = "csv", parser = x -> CSV.File(x.body)) |> DataFrame

df = filter(:timestamp => x -> Date(2019, 1, 1) <= x <= Date(2020, 1, 1), df)
price = df."close (USD)"

sma1 = sma(price; n = 10)
sma2 = sma(price; n = 20)

########################################
# Converting raw data to signal generation
########################################

struct MAH{T}
    ma::T
    mah::T
end

function build_mah(v)
    T = eltype(v)
    res = Vector{MAH{T}}(undef, length(v))
    @inbounds res[1] = MAH(v[1], T(NaN))
    @inbounds for i in 2:length(v)
        res[i] = MAH(v[i], v[i - 1])
    end
    return res
end

sma1 = build_mah(sma1)
sma2 = build_mah(sma2)

function cross_signal(mah1, mah2)
    (mah1.mah < mah2.mah) & (mah1.ma > mah2.ma) && return 1
    (mah1.mah > mah2.mah) & (mah1.ma < mah2.ma) && return -1
    return 0
end

########################################
# Broker and single simulation run
########################################

struct Broker
    cash::Float64
    shares::Int
end

function run(broker, price, sma1, sma2)
    @inbounds for i in eachindex(sma1)
        signal = cross_signal(sma1[i], sma2[i])
        if signal > 0
            p = price[i]
            shares = floor(Int, broker.cash/p)
            @set! broker.shares = shares
            @set! broker.cash -= shares*p
        elseif signal < 0
            broker.shares > 0 || continue # no shorting
            p = price[i]
            @set! broker.cash += broker.shares*p
            @set! broker.shares = 0
        end
    end

    # Cash out
    p = price[end]
    @set! broker.cash += broker.shares*p
    @set! broker.shares = 0

    return broker
end

broker = Broker(1_000_000, 0)

@btime run($broker, $price, $sma1, $sma2)
  # 352.389 ns (0 allocations: 0 bytes)

########################################
# Grid run
########################################

smas = map(n -> sma(price; n = n), 10:100) .|> build_mah

struct SimInfo
    i::Int
    j::Int
    x::Float64
end

function run(broker, price, smas)
    res = SimInfo[]
    for i in 1:length(smas)
        for j in i+1:length(smas)
            sma1 = smas[i]
            sma2 = smas[j]
            out = run(broker, price, sma1, sma2)
            push!(res, SimInfo(i, j, out.cash))
        end
    end
    return res
end

@btime run($broker, $price, $smas);
  # 1.635 ms (12 allocations: 192.55 KiB)

########################################
# Final heatmap
########################################

res = run(broker, price, smas)

res_dual = map(x -> SimInfo(x.j, x.i, x.x), res)

res_diag = map(i -> SimInfo(i, i, broker.cash), 1:length(smas))

res = vcat(res, res_dual, res_diag)
sort!(res, by = x -> (x.i, x.j))
vals = map(x -> x.x/broker.cash - 1, reshape(res, :, length(smas)))
heatmap(vals, c=cgrad([:blue,:gray,:red]))
savefig("BTC_DMAC.png")

ext: jl


text:
dop = DOP(2.3468360693153523, 0, 0, 0, 0)
projected_position_ECI2ECEF = PVTSolution(
    ECEF([4.0186869577910933e6, 427050.24116486, 4.91825870249732e6]), -2.133209701701266e7, dop, 1:length(satellite_states), [], [])
prev_PVT = PVTSolution()
min_elevation =  0

@testset "User position ECI2ECEF" begin
    positions = map(x -> sat_position_ECI_2_ECEF(x), satellite_states)
    pseudo_ranges = PositionVelocityTime.pseudo_ranges(satellite_states)
    out = user_position(prev_PVT.pos, positions, pseudo_ranges)
    @test out.pos ≈ projected_position_ECI2ECEF.pos
    @test out.receiver_time_correction ≈ projected_position_ECI2ECEF.receiver_time_correction
    @test get_gdop(out) ≈ get_gdop(projected_position_ECI2ECEF)

    out = false
    try
        calc_PVT([sv1_struct], prev_PVT, min_elevation)
    catch e
        println(e)
        if typeof(e) == PositionVelocityTime.SmallData
            out = true
        end
    end
    @test out == true
end

projected_position_ECEF = PVTSolution(
    ECEF([4.0186869577910956e6, 427050.2411648569, 4.918258702497325e6]),
    -2.133209701701266e7, DOP(2.346836069315351, 0, 0, 0, 0), 1:length(satellite_states), [], [])
@testset "User position ECEF" begin
    positions = map( x -> sat_position_ECEF(x), satellite_states)
    pseudo_ranges = PositionVelocityTime.pseudo_ranges(satellite_states)
    out = user_position(prev_PVT.pos, positions, pseudo_ranges)
    @test out.pos ≈ projected_position_ECEF.pos
    @test out.receiver_time_correction ≈ projected_position_ECEF.receiver_time_correction
    @test get_gdop(out) ≈ get_gdop(projected_position_ECI2ECEF)

    out = false
    try
        calc_PVT([sv1_struct], prev_PVT, min_elevation)
    catch e
        if typeof(e) == PositionVelocityTime.SmallData
            out = true
        end
    end
    @test out == true
end


text:
using TableShowUtils
using Test
using DataValues
using Dates

@testset "TableShowUtils" begin

source = [(a=1,b="A"),(a=2,b="B")]

@test sprint(TableShowUtils.printtable, source, "foo file") == """
2x2 foo file
a │ b
──┼──
1 │ A
2 │ B"""

@test sprint(TableShowUtils.printHTMLtable, source) == """

a	b

1	"A"

2	"B"

"""

@test sprint((stream) -> TableShowUtils.printtable(stream, source, "foo file", force_unknown_rows = true)) == "?x2 foo file\na │ b\n──┼──\n1 │ A\n2 │ B\n... with more rows"

source_with_many_columns = [(a0=1,b0=1,c0=1,a1=1,b1=1,c1=1,a2=1,b2=1,c2=1,a3=1,b3=1,c3=1,a4=1,b4=1,c4=1,a5=1,b5=1,c5=1,a6=1,b6=1,c6=1,a7=1,b7=1,c7=1,a8=1,b8=1,c8=1,a9=1,b9=1,c9=1,a10=1,b10=1,c10=1)]
@test sprint(TableShowUtils.printtable, source_with_many_columns, "foo file") == "1x33 foo file\na0 │ b0 │ c0 │ a1 │ b1 │ c1 │ a2 │ b2 │ c2 │ a3 │ b3 │ c3 │ a4 │ b4 │ c4 │ a5\n───┼────┼────┼────┼────┼────┼────┼────┼────┼────┼────┼────┼────┼────┼────┼───\n1  │ 1  │ 1  │ 1  │ 1  │ 1  │ 1  │ 1  │ 1  │ 1  │ 1  │ 1  │ 1  │ 1  │ 1  │ 1 \n... with 17 more columns: b5, c5, a6, b6, c6, a7, b7, c7, a8, b8, c8, a9, b9, c9, a10, b10, c10"

source_with_NA = [(a=1,b="A"),(a=2,b=NA)]
@test sprint(TableShowUtils.printtable, source_with_NA, "foo file") == "2x2 foo file\na │ b  \n──┼────\n1 │ A  \n2 │ #NA"

@test sprint((stream) -> TableShowUtils.printHTMLtable(stream, source, force_unknown_rows = true)) == "

a	b

1	"A"

2	"B"

⋮	⋮

... with more rows.
"

@test sprint(TableShowUtils.printdataresource, source) == "{\"schema\":{\"fields\":[{\"name\":\"a\",\"type\":\"integer\"},{\"name\":\"b\",\"type\":\"string\"}]},\"data\":[{\"a\":1,\"b\":\"A\"},{\"a\":2,\"b\":\"B\"}]}"

@test sprint(TableShowUtils.printdataresource, source_with_NA) == "{\"schema\":{\"fields\":[{\"name\":\"a\",\"type\":\"string\"},{\"name\":\"b\",\"type\":\"string\"}]},\"data\":[{\"a\":1,\"b\":\"A\"},{\"a\":2,\"b\":null}]}"

@test TableShowUtils.julia_type_to_schema_type(AbstractFloat) == "number"
@test TableShowUtils.julia_type_to_schema_type(Bool) == "boolean"
@test TableShowUtils.julia_type_to_schema_type(Dates.Time) == "time"
@test TableShowUtils.julia_type_to_schema_type(Dates.Date) == "date"
@test TableShowUtils.julia_type_to_schema_type(Dates.DateTime) == "datetime"
@test TableShowUtils.julia_type_to_schema_type(DataValues.DataValue{Integer}) == "integer"

end


text:
# tensornetwork.jl
#
# Definition and implementation of a TensorNetwork type that
# contains a list of tensors together with a specification of how they
# should be contracted, along with methods to perform the actual contraction
# and to optimize the contraction order.
#
type TensorNetwork{S,P,T,L} <: AbstractTensorNetwork{S,P,T}
    tensors::Vector{AbstractTensor{S,P,T}}
    labels::Vector{Vector{L}}
    conjs::Vector{Char}
    contracttree::Tuple
    outputorder::Vector{L}
    outputspaces::Dict{L,S}
end
function TensorNetwork{S,P,T,L}(tensors::Vector{AbstractTensor{S,P,T}},labels::Vector{Vector{L}},conjs::Vector{Char},contracttree::Tuple,outputorder::Vector{L}=Array(L,0))
    # basic concistency check on arguments
    numtensors=length(tensors)
    numtensors>1 || throw(Argument("number of tensors should be larger than one for a proper tensor network"))
    length(labels)==numtensors || throw(ArgumentError("number of label lists does not equal number of tensors"))
    length(conjs)==numtensors || throw(ArgumentError("number of conjugation identifiers does not equal number of tensors"))
    for k=1:numtensors
        numind(tensors[k])==length(labels[k]) || throw(ArgumentError("number of labels should equal number of indices for tensor $k"))
        conjs[k]=='N' || conjs[k]=='C' || throw(ArgumentError("conjugation identifiers should equal 'N' or 'C'"))
    end

    # analyse contraction network
    alllabels=vcat(labels...)
    numindices=length(alllabels)
    alllabels=unique(alllabels)
    numlabels=length(alllabels)

    table=zeros(Int,(numlabels,2))
    for n=1:numtensors
        length(labels[n])==length(unique(labels[n])) || throw(ArgumentError("duplicate labels for tensor $n: handle inner contraction first with tensortrace"))
        indn=findin(alllabels,labels[n])
        for i in indn
            if table[i,1]==0
                table[i,1]=n
            elseif table[i,2]==0
                table[i,2]=n
            else
                throw(ArgumentError("no label should appear more than two times"))
            end
        end
    end

    # build lookup table for output spaces
    outputspaces=Dict{L,S}()
    sizehint!(outputspaces,2*numlabels-numindices)
    for i=1:numindices
        if table[i,2]==0
            n=table[i,1]
            l=alllabels[i]
            s=space(tensors[n],findfirst(labels[n],l))
            outputspaces[l]=conjs[n]=='C' ? conj(s) : s
        end
    end

    # check / create output order
    outputlabels=collect(keys(outputspaces))
    if isempty(outputorder)
        outputorder=sort!(outputlabels)
    else
        isperm(indexin(outputorder,outputlabels)) || throw(ArgumentError("output order does not match with uncontracted indices"))
    end

    # check contraction tree to visit every tensor once
    visited=zeros(Int,numtensors)
    visit(node::Tuple)=(visit(node[1]);visit(node[2]))
    visit(leaf::Int)=(visited[leaf]+=1)
    visit(contracttree)
    all(visited.==1) || throw(ArgumentError("invalid contractrion tree"))

    TensorNetwork{S,P,T,L}(tensors,labels,conjs,contracttree,outputorder,outputspaces)
end
function TensorNetwork{S,P,T,L}(tensors::Vector{AbstractTensor{S,P,T}},labels::Vector{Vector{L}},conjs::Vector{Char},contractorder::Vector{L}=Array(L,0),outputorder::Vector{L}=Array(L,0))
    # basic concistency check on arguments
    numtensors=length(tensors)
    numtensors>1 || throw(Argument("number of tensors should be larger than one for a proper tensor network"))
    length(labels)==numtensors || throw(ArgumentError("number of label lists does not equal number of tensors"))
    length(conjs)==numtensors || throw(ArgumentError("number of conjugation identifiers does not equal number of tensors"))
    for k=1:numtensors
        numind(tensors[k])==length(labels[k]) || throw(ArgumentError("number of labels should equal number of indices for tensor $k"))
        conjs[k]=='N' || conjs[k]=='C' || throw(ArgumentError("conjugation identifiers should equal 'N' or 'C'"))
    end

    # analyse contraction network
    alllabels=vcat(labels...)
    numindices=length(alllabels)
    alllabels=unique(alllabels)
    numlabels=length(alllabels)

    table=zeros(Int,(numlabels,2))
    for n=1:numtensors
        length(labels[n])==length(unique(labels[n])) || throw(ArgumentError("duplicate labels for tensor $n: handle inner contraction first with tensortrace"))
        indn=findin(alllabels,labels[n])
        for i in indn
            if table[i,1]==0
                table[i,1]=n
            elseif table[i,2]==0
                table[i,2]=n
            else
                throw(ArgumentError("no label should appear more than two times"))
            end
        end
    end

    # build lookup table for output spaces
    outputspaces=Dict{L,S}()
    sizehint!(outputspaces,2*numlabels-numindices)
    for i=1:numindices
        if table[i,2]==0
            n=table[i,1]
            l=alllabels[i]
            s=space(tensors[n],findfirst(labels[n],l))
            outputspaces[l]=conjs[n]=='C' ? conj(s) : s
        end
    end

    # check / create output order
    outputlabels=collect(keys(outputspaces))
    if isempty(outputorder)
        outputorder=sort!(outputlabels)
    else
        isperm(indexin(outputorder,outputlabels)) || throw(ArgumentError("output order does not match with uncontracted indices"))
    end

    # check / create contraction order
    contractlabels=setdiff(alllabels,outputlabels)
    if isempty(contractorder)
        contractorder=sort!(contractlabels,rev=true)
    else
        isperm(indexin(contractorder,contractlabels)) || throw(ArgumentError("contraction order does not match with uncontracted indices"))
    end

    # create corresponding contraction tree
    indexorder=indexin(contractorder,alllabels)
    lookup=(Int=>Any)[n=>n for n=1:numtensors]
    for i in indexorder
        n1,n2=table[i,1],table[i,2]
        node1=lookup[n1]
        node2=lookup[n2]
        if node1!=node2
            node=tuple(node1,node2)
            lookup[n1]=node
            lookup[n2]=node
        end
    end
    nodes=unique(values(lookup))
    contracttree=nodes[1]
    for k=2:length(nodes)
        contracttree=tuple(tree,nodes[k])
    end

    TensorNetwork{S,P,T,L}(tensors,labels,conjs,contracttree,outputorder,outputspaces)
end

# Convenience constructor
#-------------------------
function network(tensors::Vector,labels::Vector,conjs::Vector{Char}=fill('N',length(tensors));contractorder::Vector=[],outputorder::Vector=[],contracttree::Tuple=())
    length(tensors)>1 || throw(ArgumentError("number of tensors should be larger than one for a proper tensor network"))
    length(tensors)==length(labels) || throw(ArgumentError("number of label lists does not equal number of tensors"))
    length(tensors)==length(conjs) || throw(ArgumentError("number of conjugation identifiers does not equal number of tensors"))
    for n=1:length(tensors)
        isa(tensors[n],AbstractTensor) || throw(ArgumentError("not a list of tensors"))
        isa(labels[n],Vector) || throw(ArgumentError("not a list of labels"))
    end
    Slist=map(spacetype,tensors)
    S=Slist[1]
    all(Slist .== S) || throw(ArgumentError("not a homogeneous list of tensors"))
    Plist=map(tensortype,tensors)
    P=Plist[1]
    all(Plist .== P) || throw(ArgumentError("not a homogeneous list of tensors"))
    T = eltype(tensors[1])
    L = eltype(labels[1])
    for n=2:length(tensors)
        T=promote_type(T,eltype(tensors[n]))
        L=promote_type(L,eltype(labels[n]))
    end

    if isempty(contracttree)
        return TensorNetwork(convert(Vector{AbstractTensor{S,P,T}},tensors),convert(Vector{Vector{L}},labels),convert(Vector{L},contractorder),convert(Vector{L},outputorder))
    elseif isempty(contractorder)
        return TensorNetwork(convert(Vector{AbstractTensor{S,P,T}},tensors),convert(Vector{Vector{L}},labels),contracttree,convert(Vector{L},outputorder))
    else
        throw(ArgumentError("provide contractorder or contracttree, not both"))
    end
end

# Basic properties
#------------------
Base.length(network::TensorNetwork)=length(network.tensors)

# Basic methods for characterising the output tensor:
#-----------------------------------------------------
numind(network::TensorNetwork)=length(network.outputorder)
space(network::TensorNetwork,ind::Integer)=network.outputspaces[network.outputorder[ind]]
space(network::TensorNetwork)=ntuple(ind->space(network,ind),numind(network))

# Compute the output tensor:
#----------------------------
function tensor(network::TensorNetwork)
    tree=network.contracttree
    t1,c1,labels1=_contract(network,tree[1])
    t2,c2,labels2=_contract(network,tree[2])
    return tensorcontract(t1,labels1,c1,t2,labels2,c2,network.outputorder)
end
tensorcontract(network::TensorNetwork)=tensor(network)

# Internal methods to perform the contraction
#---------------------------------------------
_contract(network::TensorNetwork,n::Int)=(network.tensors[n],network.conjs[n],network.labels[n])
function _contract(network::TensorNetwork,tree::Tuple)
    t1,conj1,labels1=_contract(network,tree[1])
    t2,conj2,labels2=_contract(network,tree[2])

    newlabels=symdiff(labels1,labels2)
    if conj1==conj2
        newconj=conj1
        conj1='N'
        conj2='N'
    else
        newconj='N'
    end

    return tensorcontract(t1,labels1,conj1,t2,labels2,conj2,newlabels),newconj,newlabels
end


text:
using VortexDistributions, Plots

Nx=1000; Ny=1000;
θx = LinRange(0,4*pi,Nx); θy = LinRange(0,4*pi,Ny);
psix = exp.(im*θx) + 0.1*(randn(Nx)+im*randn(Nx)); psiy = exp.(im*θy) + 0.1*(randn(Ny)+im*randn(Ny));
ϕx = angle.(psix); ϕy = angle.(psiy);
ϕ = [ϕx ϕy]
ϕu = unwrap(ϕ)
ϕui = similar(ϕu)
unwrap!(ϕui,ϕ)
@test ϕui == ϕu

ϕ = zeros(2, Nx)
ϕ[1, :] .= ϕx
ϕ[2, :] .=ϕy
ϕu = unwrap(ϕ,2)
ϕui = similar(ϕu)
unwrap!(ϕui,ϕ,2)
@test ϕui == ϕu


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdspin

"""
    'feffery'_antdspin(;kwargs...)
    'feffery'_antdspin(children::Any;kwargs...)
    'feffery'_antdspin(children_maker::Function;kwargs...)

An AntdSpin component.

Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): The content of the tab - will only be displayed if this tab is selected
- `id` (String; optional)
- `className` (String; optional)
- `debug` (Bool; optional)
- `delay` (Real; optional)
- `excludeProps` (Array of Strings; optional)
- `includeProps` (Array of Strings; optional)
- `listenPropsMode` (a value equal to: 'default', 'exclude', 'include'; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `size` (a value equal to: 'small', 'middle', 'large'; optional)
- `spinning` (Bool; optional)
- `style` (Dict; optional)
- `text` (String; optional)
- `wrapperClassName` (String; optional)
"""
function 'feffery'_antdspin(; kwargs...)
        available_props = Symbol[:children, :id, :className, :debug, :delay, :excludeProps, :includeProps, :listenPropsMode, :loading_state, :size, :spinning, :style, :text, :wrapperClassName]
        wild_props = Symbol[]
        return Component("'feffery'_antdspin", "AntdSpin", "feffery_antd_components", available_props, wild_props; kwargs...)
end

'feffery'_antdspin(children::Any; kwargs...) = 'feffery'_antdspin(;kwargs..., children = children)
'feffery'_antdspin(children_maker::Function; kwargs...) = 'feffery'_antdspin(children_maker(); kwargs...)

max_stars_repo_name: limmyforget/feffery-antd-components

max_issues_repo_name: limmyforget/feffery-antd-components

max_forks_repo_name: limmyforget/feffery-antd-components


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module Foo
import Bar, Baz, Qux
this = "Foo1"
which = "depot"
end

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
function solve(
    prob::DiffEqBase.AbstractODEProblem{uType,tupType,isinplace},
    alg::AlgType,
    timeseries=[], ts=[], ks=[];
    verbose=true,
    save_timeseries=nothing,
    saveat=eltype(tupType)[], reltol = 1e-5, abstol = 1e-8,
    save_everystep=isempty(saveat),
    dense = save_everystep && isempty(saveat),
    save_start = save_everystep || isempty(saveat) || typeof(saveat) <: Number ?
                 true : prob.tspan[1] in saveat,
    callback=nothing,
    dtmin = abs(prob.tspan[2]-prob.tspan[1])/1e-9,
    dtmax = abs(prob.tspan[2]-prob.tspan[1])/2.5,
    timeseries_errors=true, dense_errors=false,
    dt = 0.0, norm = LinearAlgebra.norm,
    kwargs...) where {uType,tupType,isinplace,AlgType<:ODEjlAlgorithm}

    tType = eltype(tupType)

    if verbose
        warned = !isempty(kwargs) && check_keywords(alg, kwargs, warnlist)
        if !(typeof(prob.f) <: DiffEqBase.AbstractParameterizedFunction) &&
            typeof(alg) <: ode23s
            if DiffEqBase.has_tgrad(prob.f)
                @warn("Explicit t-gradient given to this stiff solver is ignored.")
                warned = true
            end
            if DiffEqBase.has_jac(prob.f)
                @warn("Explicit Jacobian given to this stiff solver is ignored.")
                warned = true
            end
        end
        warned && warn_compat()
    end

    if save_timeseries !== nothing
        verbose && @warn("save_timeseries is deprecated. Use save_everystep instead")
        save_everystep = save_timeseries
    end

    if prob.f.mass_matrix != I
        error("This solver is not able to use mass matrices.")
    end

    if callback !== nothing || :callback ∈ keys(prob.kwargs)
        error("ODE is not compatible with callbacks.")
    end

    tspan = prob.tspan

    if typeof(saveat) <: Number
      if (tspan[1]:saveat:tspan[end])[end] == tspan[end]
        saveat_vec = convert(Vector{tType},collect(tType,tspan[1]+saveat:saveat:tspan[end]))
      else
        saveat_vec = convert(Vector{tType},collect(tType,tspan[1]+saveat:saveat:(tspan[end]-saveat)))
      end
    else
        saveat_vec = convert(Vector{tType}, collect(saveat))
    end

    if !isempty(saveat_vec) && saveat_vec[end] == tspan[2]
        pop!(saveat_vec)
    end

    if !isempty(saveat_vec) && saveat_vec[1] == tspan[1]
        Ts = unique([saveat_vec;tspan[2]])
    else
        Ts = unique([tspan[1];saveat_vec;tspan[2]])
    end

    if save_everystep
        points = :all
    else
        points = :specified
    end

    sizeu = size(prob.u0)
    p = prob.p

    if isinplace
        f = (t,u) -> (du = zero(u); prob.f(du,u,p,t); vec(du))
    elseif uType <: AbstractArray
        f = (t,u) -> vec(prob.f(reshape(u,sizeu),p,t))
    else
        f = (t,u) -> prob.f(u,p,t)
    end

    u0 = uType <: AbstractArray ? vec(prob.u0) : prob.u0

    # Calling the solver, i.e. if the algorithm is ode45,
    # then AlgType(...) is ode45(...)
    ts, timeseries_tmp = AlgType(f,u0,Ts;
                                 norm = norm,
                                 abstol=abstol,
                                 reltol=reltol,
                                 maxstep=dtmax,
                                 minstep=dtmin,
                                 initstep=dt,
                                 points=points)

    if save_start
        start_idx = 1 # The index to start making the timeseries from
    else
        start_idx = 2
        ts = ts[2:end]
    end

    # Reshape the result if needed
    if uType <: AbstractArray
        timeseries = Vector{uType}()
        for i=start_idx:length(timeseries_tmp)
            push!(timeseries,reshape(timeseries_tmp[i],sizeu))
        end
    else
        timeseries = timeseries_tmp
    end

    DiffEqBase.build_solution(prob,alg,ts,timeseries,
                   timeseries_errors = timeseries_errors,
                   dense_errors = dense_errors,
                   retcode = :Succss)
end # function solve


text:
module SpikingNeuralNetworks

export SNN
const SNN = SpikingNeuralNetworks
using Requires
using UnPack

using LinearAlgebra
using SparseArrays
try
    using Reexport
    using Parameters
    using Requires

catch
    import Pkg; Pkg.add("Reexport")
    import Pkg; Pkg.add("Parameters")
    import Pkg; Pkg.add("Requires")
    using Reexport
    using Parameters
    using Requires
end
#using Unitful
#using Unitful.DefaultSymbols
#@reexport using Utils

const SNNInt = Int32
const SNNFloat = Float32
# srand(1000)
include("unit.jl")
#include("units.jl")
include("main.jl")
include("util.jl")

include("neuron/if.jl")
include("neuron/if2.jl")
include("neuron/noisy_if.jl")
include("neuron/poisson.jl")
include("neuron/iz.jl")
include("neuron/iz_more_param.jl")
include("neuron/adexp.jl")

include("neuron/hh.jl")
include("neuron/rate.jl")

include("synapse/spiking_synapse.jl")
include("synapse/rate_synapse.jl")
include("synapse/fl_synapse.jl")
include("synapse/fl_sparse_synapse.jl")
include("synapse/pinning_synapse.jl")
include("synapse/pinning_sparse_synapse.jl")

function __init__()
    @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("plot.jl")
end

end


text:
using DiffEqBase, Test
using Distributions

@test DiffEqBase.promote_tspan((0.0,1.0)) == (0.0,1.0)
@test DiffEqBase.promote_tspan((0,1.0)) == (0.0,1.0)
@test DiffEqBase.promote_tspan(1.0) == (0.0,1.0)
@test DiffEqBase.promote_tspan(nothing) == (nothing,nothing)
@test DiffEqBase.promote_tspan(Real[0,1.0]) == (0.0,1.0)

prob = ODEProblem((u,p,t)->u,(p,t0)->p[1],(p)->(0.0,p[2]),(2.0,1.0))
prob2 = DiffEqBase.get_concrete_problem(prob,true)

@test prob2.u0 == 2.0
@test prob2.tspan == (0.0,1.0)

prob = ODEProblem((u,p,t)->u,(p,t)->Normal(p,1),(0.0,1.0),1.0)
prob2 = DiffEqBase.get_concrete_problem(prob,true)
@test typeof(prob2.u0) == Float64

kwargs(; kw...) = kw
prob = ODEProblem((u,p,t)->u,1.0,nothing)
prob2 = DiffEqBase.get_concrete_problem(prob,true,tspan=(1.2,3.4))
@test prob2.tspan === (1.2,3.4)

prob = ODEProblem((u,p,t)->u,nothing,nothing)
prob2 = DiffEqBase.get_concrete_problem(prob,true,u0=1.01,tspan=(1.2,3.4))
@test prob2.u0 === 1.01

prob = ODEProblem((u,p,t)->u,1.0,(0,1))
prob2 = DiffEqBase.get_concrete_problem(prob,true)
@test prob2.tspan == (0.0,1.0)

prob = DDEProblem((u, h, p, t) -> - h(p, t - p[1]), (p, t0) -> p[2], (p, t) -> 0,
                  (p) -> (0.0, p[3]), (1.0, 2.0, 3.0); constant_lags = (p) -> [p[1]])
prob2 = DiffEqBase.get_concrete_problem(prob, true)

@test prob2.u0 == 2.0
@test prob2.tspan == (0.0, 3.0)
@test prob2.constant_lags == [1.0]


text:
__precompile__(true)
"""
Regex functions for Str strings

Copyright 2018-2019 Gandalf Software, Inc., Scott P. Jones, and other contributors to the Julia language
Licensed under MIT License, see LICENSE.md
Based in part on julia/base/regex.jl and julia/base/pcre.jl
"""
module StrRegex

import Base.Threads

using ModuleInterfaceTools

const BASE_REGEX_MT = isdefined(Base.PCRE, :PCRE_COMPILE_LOCK)

@api extend! StrBase

@api base Regex, match, compile, eachmatch, findall, count

@api public RegexStr, RegexStrMatch, "@r_str", "@R_str"

const _not_found = StrBase._not_found

using Base: RefValue, replace_err, SubstitutionString

using PCRE2
const PCRE = PCRE2

const deb = RefValue(false)

const MATCH_CONTEXT = [(C_NULL,C_NULL,C_NULL)]

codeunit_index(::Type{UInt8})  = 1
codeunit_index(::Type{UInt16}) = 2
codeunit_index(::Type{UInt32}) = 3

match_context(::Type{T}, tid) where {T<:CodeUnitTypes} =
    @inbounds MATCH_CONTEXT[tid][codeunit_index(T)]

const JIT_STACK_START_SIZE = 32768
const JIT_STACK_MAX_SIZE = 1048576

function __init__()
    if (n = Threads.nthreads()) != 1
        resize!(MATCH_CONTEXT, n)
        fill!(MATCH_CONTEXT, (C_NULL, C_NULL, C_NULL))
    end
end

# UCP, UTF and NO_UTF_CHECK are based on the string type
const DEFAULT_COMPILER_OPTS = PCRE.ALT_BSUX
const DEFAULT_MATCH_OPTS    = 0%UInt32

const Binary_Regex_CSEs = Union{ASCIICSE,BinaryCSE,Text1CSE,Text2CSE,Text4CSE}
const Regex_CSEs = Union{Binary_Regex_CSEs,Latin_CSEs,UTF8_CSEs,UCS2_CSEs,UTF16CSE, UTF32_CSEs}

const _VALID = PCRE.NO_UTF_CHECK
const _UTF   = PCRE.UTF
const _UCP   = PCRE.UCP

_clear_opts(opt) = UInt32(opt) & ~(_VALID | _UTF | _UCP)

const comp_add   = (0%UInt32, _UCP, _UTF, _UCP|_UTF, _UTF|_VALID, _UCP|_UTF|_VALID)
const match_add  = (0%UInt32, 0%UInt32, 0%UInt32, 0%UInt32, _VALID, _VALID)

@noinline _check_compile(options) =
    (options & ~PCRE.COMPILE_MASK) == 0 ? UInt32(options) :
    throw(ArgumentError("invalid regex compile options: $options"))
@noinline _check_match(options) =
    (options & ~PCRE.MATCH_MASK) == 0 ? UInt32(options) :
    throw(ArgumentError("invalid regex match options: $options"))

function finalize! end

fin(exp) = finalizer(finalize!, exp) ; exp

# There are 18 valid combinations of size (8,16,32), UTF/no-UTF, possibly UCP, and NO_UTF_CHECK
#
#  1 - Text1/Binary/ASCII   8
#  2 - Latin/_Latin         8,UCP
#  3                        8,UTF
#  4 - RawUTF8CSE           8,UCP,UTF
#  5                        8,UTF,NO_UTF_CHECK
#  6 - UTF8                 8,UCP,UTF,NO_UTF_CHECK

#  1 - Text2               16
#  2 - UCS2/_UCS2          16,UCP
#  3                       16,UTF
#  4 - RawUTF16CSE         16,UCP,UTF
#  5                       16,UTF,NO_UTF_CHECK
#  6 - UTF16               16,UCP,UTF,NO_UTF_CHECK

#  1 - Text4               32
#  2                       32,UCP
#  3                       32,UTF
#  4 - RawUTF32CSE         32,UCP,UTF
#  5                       32,UTF,NO_UTF_CHECK
#  6 - UTF32/_UTF32        32,UCP,UTF,NO_UTF_CHECK

opt_index(::Type{C}) where {C<:CSE} = 1
# opt_index(::Type{C}) where {C<:Union{ASCIICSE, BinaryCSE, Text1CSE, Text2CSE, Text4CSE}} = 1
opt_index(::Type{C}) where {C<:Union{LatinCSE, _LatinCSE, UCS2CSE, _UCS2CSE}} = 2
opt_index(::Type{C}) where {C<:Union{RawUTF8CSE, RawUTF16CSE}} = 4
opt_index(::Type{C}) where {C<:Union{UTF8CSE, UTF16CSE, UTF32CSE, _UTF32CSE}} = 6

opt_index(::Type{S}) where {S<:AbstractString} = opt_index(cse(S))

# Match tables only need to be allocated for each size (8,16,32) and for each thread
mutable struct MatchTab
    match_data::NTuple{3, Ptr{Cvoid}}
    ovec::NTuple{3, Vector{Csize_t}}
    MatchTab() = new((C_NULL, C_NULL, C_NULL), (Csize_t[], Csize_t[], Csize_t[]))
end

md_free(::Type{T}, md) where {T<:CodeUnitTypes} =
    md == C_NULL || PCRE.match_data_free(T, md)

function finalize!(mt::MatchTab)
    md_free(UInt8,  mt.match_data[1])
    md_free(UInt16, mt.match_data[2])
    md_free(UInt32, mt.match_data[3])
    mt.match_data = (C_NULL, C_NULL, C_NULL)
    mt.ovec = (Csize_t[], Csize_t[], Csize_t[])
end

const empty_table = (C_NULL, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL)

mutable struct RegexStr
    pattern::String
    compile_options::UInt32
    match_options::UInt32
    negated_options::UInt32
    table::Vector{NTuple{6, Ptr{Cvoid}}}
    match::Vector{MatchTab}

    function RegexStr(pattern::AbstractString,
                      compile_options::Integer,
                      match_options::Integer,
                      negated_options::Integer=0)
        re = new(String(pattern),
                 _check_compile(compile_options),
                 _check_match(match_options),
                 _check_compile(negated_options),
                 [empty_table for i=1:3],
                 [MatchTab() for i=1:Threads.nthreads()])
        compile(cse(pattern), pattern, re)
        fin(re)
    end
end

const tabtype = (UInt8,UInt16,UInt32)

function _finalize!(re::RegexStr, i)
    tab = re.table[i]
    typ = tabtype[i]
    for j = 1:6
        (r = tab[j]) == C_NULL || PCRE.code_free(typ, r)
    end
    re.table[i] = empty_table
end

function finalize!(re::RegexStr)
    _finalize!(re, 1)
    _finalize!(re, 2)
    _finalize!(re, 3)
    for m in re.match
        finalize!(m)
    end
end

_update_match(t, v, n) =
    (n == 1 ? (v, t[2], t[3]) : n == 2 ? (t[1], v, t[3]) : (t[1], t[2], v))

_update_table(t, v, n) =
    (n == 1 ? v : t[1], n == 2 ? v : t[2], n == 3 ? v : t[3],
     n == 4 ? v : t[4], n == 5 ? v : t[5], n == 6 ? v : t[6])

const RegexTypes = Union{Regex, RegexStr}

@noinline _mode_err()     = throw(ArgumentError("a and u mode not allowed together"))
@noinline _unknown_err(f) = throw(ArgumentError("unknown regex flag: $f"))

function _add_compile_options(flags)
    negated = 0%UInt32
    options = DEFAULT_COMPILER_OPTS
    for f in flags
        if f == 'a'
            (options & PCRE.UCP) && _mode_err()
            negated |= PCRE.UCP
        elseif f == 'u'
            (negated & PCRE.UCP) && _mode_err()
            options |= PCRE.UCP
        elseif f == 'i' ; options |= PCRE.CASELESS
        elseif f == 'm' ; options |= PCRE.MULTILINE
        elseif f == 's' ; options |= PCRE.DOTALL
        elseif f == 'x' ; options |= PCRE.EXTENDED
        else
            _unknown_err(f)
        end
    end
    options, negated
end

RegexStr(pattern::AbstractString) =
    RegexStr(pattern, DEFAULT_COMPILER_OPTS, 0%UInt32, DEFAULT_MATCH_OPTS)
RegexStr(pattern::AbstractString, flags::AbstractString) =
    RegexStr(pattern, _add_compile_options(flags)..., DEFAULT_MATCH_OPTS)

@static if isdefined(Main, :LETS_BE_PIRATES) && Core.eval(Main, :LETS_BE_PIRATES)
    import Base.@r_str
    macro r_str(pattern::ANY, flags...) ; cmp_all(RegexStr(pattern, flags...)) ; end
    # Yes, this is type piracy, but it is needed to make all string types work together easily
    Base.Regex(pattern::AbstractString, co::Integer, mo::Integer, no::Integer=0) =
        RegexStr(pattern, co, mo, no)
    Base.Regex(pattern::AbstractString, flags::AbstractString) = RegexStr(pattern, flags)
    Base.Regex(pattern::AbstractString) = RegexStr(pattern)
else
    Regex(pattern::MaybeSub{<:Str}, co, mo, no=0) = RegexStr(pattern, co, mo, no)
    Regex(pattern::MaybeSub{<:Str}, flags::AbstractString) = RegexStr(pattern, flags)
    Regex(pattern::MaybeSub{<:Str}) = RegexStr(pattern)
end

function check_compile(::Type{C}, re::RegexStr) where {C<:CSE}
    try
        compile(C, re)
    catch ex
        (ex isa PCRE.PCRE2_Error && ex.errno in (134, 177)) || rethrow()
    end
end

"""Precompile all of the common cases needed for String, UTF8Str, UTF16Str, and UniStr"""
function cmp_all(re::RegexStr)
    pat = re.pattern
    is_bmp(pat)   && check_compile(UCS2CSE,  re)
    is_latin(pat) && check_compile(LatinCSE, re)
    if is_unicode(pat)
        check_compile(RawUTF8CSE, re)
        check_compile(UTF32CSE, re)
        check_compile(UTF16CSE, re)
        check_compile(UTF8CSE,  re)
    end
    is_ascii(pat) && check_compile(ASCIICSE, re)
    re
end

macro R_str(pattern, flags...) cmp_all(RegexStr(pattern, flags...)) end

function show(io::IO, re::RegexStr)
    imsx = PCRE.CASELESS|PCRE.MULTILINE|PCRE.DOTALL|PCRE.EXTENDED|PCRE.UCP
    opts = re.compile_options
    neg  = re.negated_options
    if (opts & ~imsx) == DEFAULT_COMPILER_OPTS || neg != 0
        print(io, 'R')
        Base.print_quoted_literal(io, re.pattern)
        (opts & PCRE.CASELESS ) == 0 || print(io, 'i')
        (opts & PCRE.MULTILINE) == 0 || print(io, 'm')
        (opts & PCRE.DOTALL   ) == 0 || print(io, 's')
        (opts & PCRE.EXTENDED ) == 0 || print(io, 'x')
        (opts & PCRE.UCP      ) == 0 || print(io, 'u')
        (neg  & PCRE.UCP      ) == 0 || print(io, 'a')
    else
        print(io, "RegexStr(")
        show(io, re.pattern)
        print(io, ',')
        show(io, opts)
        neg == 0 || print(io, "," , neg)
        print(io, ')')
    end
end

struct RegexStrMatch{T<:AbstractString}
    match::SubString{T}
    captures::Vector{Union{Nothing,SubString{T}}}
    offset::Int
    offsets::Vector{Int}
    regex::RegexStr
end

get_regex(re::RegexStrMatch{T}) where {T<:AbstractString} =
    (C = cse(T) ; re.regex.table[codeunit_index(codeunit(C))][opt_index(C)])
get_regex(re::RegexMatch) = re.regex.regex

function show(io::IO, m::RegexStrMatch{T}) where {T}
    print(io, "RegexStrMatch{$T}(")
    show(io, m.match)
    idx_to_capture_name = PCRE.capture_names(codeunit(T), get_regex(m))
    if !is_empty(m.captures)
        print(io, ", ")
        for i = 1:length(m.captures)
            # If the capture group is named, show the name.
            # Otherwise show its index.
            print(io, get(idx_to_capture_name, i, i), "=")
            show(io, m.captures[i])
            i < length(m.captures) && print(io, ", ")
        end
    end
    print(io, ")")
end

# Capture group extraction
getindex(m::RegexStrMatch, idx::Integer) = m.captures[idx]
function getindex(m::RegexStrMatch{T}, name::Symbol) where {T}
    idx = PCRE.substring_number_from_name(codeunit(T), get_regex(m), name)
    idx <= 0 && error("no capture group named $name found in regex")
    m[idx]
end

function outtab(tab)
    for rt in tab
        print(rt == C_NULL ? "0 " : "0x" * outhex(reinterpret(UInt, rt)) * " ")
    end
end

getindex(m::RegexStrMatch, name::AbstractString) = m[Symbol(name)]

get_comp_ind(opt) =
    ifelse((opt & _UTF) == 0, 1, ifelse((opt & _VALID) == 0, 3, 5)) + ((opt & _UCP) != 0)

function _make_match(CU, re, mt, cu_index)
    md = PCRE.match_data_create_from_pattern(CU, re, C_NULL)
    ov = PCRE.get_ovec(CU, md)
    mt.match_data = _update_match(mt.match_data, md, cu_index)
    mt.ovec       = _update_match(mt.ovec, ov, cu_index)
end

function compile(::Type{C}, pattern, regex::RegexStr) where {C<:Regex_CSEs}
    deb[] && println("compile(::Type{$C}, \"$pattern\", $regex)")
    CU = codeunit(C)
    cu_index = codeunit_index(CU)
    retab = regex.table[cu_index]
    cvtcomp = (regex.compile_options | comp_add[opt_index(C)]) & ~regex.negated_options
    ind = get_comp_ind(cvtcomp)
    # Get index based on actual options set
    re = retab[ind]
    if re == C_NULL
        pat = convert(Str{C,Nothing,Nothing,Nothing}, pattern)
        if PCRE.PCRE_LOCK === nothing
            re = PCRE.compile(pat, cvtcomp)
            deb[] && println("_update_table: ", cu_index, ", ", ind, ", ", retab, ", ", re)
            regex.table[cu_index] = _update_table(retab, re, ind)
            deb[] && println(" => ", regex.table[cu_index])
            PCRE.jit_compile(CU, re)
        else
            l = PCRE.PCRE_LOCK::Threads.SpinLock
            lock(l)
            try
                # Check again while locked if some other thread has compiled this
                if retab[ind] == C_NULL
                    re = PCRE.compile(pat, cvtcomp)
                    deb[] && println("_update_table: ", cu_index, ", ", ind, ", ", retab, ", ", re)
                    regex.table[cu_index] = _update_table(retab, re, ind)
                    deb[] && println(" => ", regex.table[cu_index])
                    PCRE.jit_compile(CU, re)
                end
            finally
                unlock(l)
            end
        end
    end
    tid = Threads.threadid()
    mt = regex.match[tid]
    mt.match_data[cu_index] == C_NULL || return regex
    if PCRE.PCRE_LOCK === nothing
        _make_match(CU, re, mt, cu_index)
    else
        l = PCRE.PCRE_LOCK::Threads.SpinLock
        lock(l)
        try
            # check again under lock
            mt.match_data[cu_index] == C_NULL && _make_match(CU, re, mt, cu_index)
        finally
            unlock(l)
        end
    end
    regex
end

compile(::Type{C}, regex::RegexStr) where {C<:Regex_CSEs} = compile(C, regex.pattern, regex)

function compile(::Type{C}, regex::Regex) where {C<:Regex_CSEs}
    ind = opt_index(C)
@static if BASE_REGEX_MT
    # Keep old match type in regex.extra (which doesn't seem to ever be used)
    oldind = reinterpret(UInt, regex.extra)%Int
    ind == oldind && return regex

    re = regex.regex
    regex.compile_options = cvtcomp = _clear_opts(regex.compile_options) | comp_add[ind]
    if oldind != 0
        S = tabtype[oldind]
        re == C_NULL || (PCRE.code_free(S, re); regex.regex = C_NULL)
    else
        S = Nothing
    end
    regex.regex = re = PCRE.compile(convert(Str{C,Nothing,Nothing,Nothing}, regex.pattern),
                                    cvtcomp)
    T = codeunit(C)
    PCRE.jit_compile(T, re)
    if S !== T
        S === Nothing || md_free(S, regex.match_data)
        regex.match_data = md = PCRE.match_data_create_from_pattern(T, re, C_NULL)
        regex.ovec = PCRE.get_ovec(T, md)
    end
    regex.extra = reinterpret(Ptr{Cvoid}, ind)
else
    cvtcomp = _clear_opts(regex.compile_options) | comp_add[ind]
    deb[] && println("ind = ", ind, ", old co = ", regex.compile_options, ", co = ", cvtcomp)
    regex
end
end

"""Get a thread-specific match context"""
function get_match_context(::Type{T}, tid) where {T<:CodeUnitTypes}
    cu_index = codeunit_index(T)
    if (mc = MATCH_CONTEXT[tid][cu_index]) == C_NULL
        mc = PCRE.match_context_create(T, C_NULL)
        js = PCRE.jit_stack_create(T, JIT_STACK_START_SIZE, JIT_STACK_MAX_SIZE, C_NULL)
        PCRE.jit_stack_assign(T, mc, C_NULL, js)
        MATCH_CONTEXT[tid] = _update_match(MATCH_CONTEXT[tid], mc, cu_index)
    end
    mc
end

_exec_err(rc) = error("StrRegex.exec error: $(PCRE.err_message(rc))")

function exec(C, re::RegexStr, subject, offset, options)
    deb[] && print("exec($re, \"$subject\", $offset, $options, match_data)")
    @preserve subject begin
        pnt = pointer(subject)
        siz = ncodeunits(subject)
        CU = eltype(pnt)
        0 <= offset <= siz || boundserr(subject, offset)
        tid = Threads.threadid()
        cu_index = codeunit_index(CU)
        tab = re.table[cu_index]
        cvtcomp = (re.compile_options | comp_add[opt_index(C)]) & ~re.negated_options
        ind = get_comp_ind(cvtcomp)
        deb[] && println(" => $tab, $ind")
        rc = PCRE.match(CU, tab[ind], pnt, siz, offset,
                        re.match_options | match_add[ind] | _check_match(options),
                        re.match[tid].match_data[cu_index], get_match_context(CU, tid))
        # rc == -1 means no match, -2 means partial match.
        #dump(re)
        rc < -2 && _exec_err(rc)
        rc >= 0
    end
end

function exec(C, re::Regex, subject, offset, options)
    @preserve subject begin
        pnt = pointer(subject)
        siz = ncodeunits(subject)
        T = eltype(pnt)
        #loc = bytoff(T, offset)
        0 <= offset <= siz || boundserr(subject, offset)
        opts = re.match_options | match_add[opt_index(C)] | _check_match(options)
        rc = PCRE.match(T, re.regex, pnt, siz, offset, opts, re.match_data,
                        get_match_context(T, Threads.threadid()))
        # rc == -1 means no match, -2 means partial match.
        rc < -2 && _exec_err(rc)
        rc >= 0
    end
end

comp_exec(C, re, subject, offset, options) = exec(C, compile(C, re), subject, offset, options)

get_range(ov, str, i = 0) = Int(ov[2*i+1]+1) : prevind(str, Int(ov[2*i+2]+1))

@static BASE_REGEX_MT || (get_ovec(::Type{<:Any}, re::Regex) = re.ovec)
get_ovec(::Type{T}, re::RegexStr) where {T<:CodeUnitTypes} =
     re.match[Threads.threadid()].ovec[codeunit_index(T)]
get_ovec(::Type{C}, re::RegexStr) where {C<:CSE} = get_ovec(codeunit(C), re)

function _match(::Type{C}, re, str, idx, opts) where {C<:CSE}
    comp_exec(C, re, str, idx - 1, opts) || return nothing
    ov = get_ovec(C, re)
    n = div(length(ov),2) - 1
    rng = get_range(ov, str)
    mat = SubString(str, rng)
    cap = Union{Nothing,SubString{typeof(str)}}[ov[2*i+1] == PCRE.UNSET ? nothing :
                                                SubString(str, get_range(ov, str, i)) for i=1:n]
    RegexStrMatch(mat, cap, rng.start, Int[ ov[2*i+1]+1 for i=1:n ], re)
end

regex_type_error(T, S) =
    throw(ArgumentError("$T matching is not supported for $S; use UniStr(s) to convert"))

match(re::Regex, str::MaybeSub{<:Str}, idx::Integer, add_opts=0) =
    regex_type_error(Regex, typeof(str))

match(re::Regex, str::MaybeSub{<:Str{<:Regex_CSEs}}, idx::Integer, add_opts::UInt32=UInt32(0)) =
    _match(basecse(str), re, str, Int(idx), add_opts)

match(r::RegexStr, str::AbstractString, idx::Integer, add_opts=0) =
    regex_type_error(RegexStr, typeof(str))

match(re::RegexStr, str::MaybeSub{<:Str}, idx::Integer, add_opts::UInt32=UInt32(0)) =
    _match(basecse(str), re, str, Int(idx), add_opts)

match(re::RegexStr, str::MaybeSub{String}, idx::Integer, add_opts::UInt32=UInt32(0)) =
    _match(RawUTF8CSE, re, str, Int(idx), add_opts)

match(re::Regex, str::MaybeSub{<:Str})   = match(re, str, 1)
match(re::RegexStr, str::AbstractString) = match(re, str, 1)

@inline __find(::Type{C}, re, str, idx) where {C} =
    comp_exec(C, re, str, idx, 0) ? get_range(get_ovec(C, re), str) : _not_found

@inline _find(::Type{C}, re, str) where {C} = __find(C, re, str, 0)

@inline _find(::Type{C}, re, str, idx) where {C} =
    (idx-1 <= ncodeunits(str)
     ? __find(C, re, str, idx-1)
     : (@boundscheck boundserr(str, idx) ; return _not_found))

find(::Type{Fwd}, re::RegexTypes, str::AbstractString, idx::Integer) =
    regex_type_error(typeof(re), typeof(str))

find(::Type{Fwd}, re::RegexTypes, str::MaybeSub{<:Str{C}}, idx::Integer) where {C<:Regex_CSEs} =
    _find(C, re, str, idx)
find(::Type{Fwd}, re::RegexTypes, str::MaybeSub{<:Str{_LatinCSE}}, idx::Integer) =
    _find(LatinCSE, re, str, idx)
find(::Type{Fwd}, re::RegexTypes, str::MaybeSub{String}, idx::Integer) =
    _find(RawUTF8CSE, re, str, idx)

find(::Type{First}, re::RegexTypes, str::AbstractString) = find(Fwd, re, str, 1)

find(::Type{First}, re::RegexTypes, str::MaybeSub{<:Str{C}}) where {C<:Regex_CSEs} =
    __find(C, re, str, 0)
find(::Type{First}, re::RegexTypes, str::MaybeSub{<:Str{_LatinCSE}}) =
    __find(LatinCSE, re, str, 0)
find(::Type{First}, re::RegexTypes, str::MaybeSub{String}) =
    __find(RawUTF8CSE, re, str, 0)

find(::Type{All}, re::RegexTypes, str::AbstractString, idx::Integer) =
    findall(re, str, idx)
find(::Type{All}, re::RegexTypes, str::AbstractString) = find(All, re, str, 1)

# Borrow idea from @dalum on GitHub (sakse on Julia Discourse), PR #29790
starts_with(s::AbstractString, r::RegexStr) =
    comp_exec(UTF8CSE, r, UTF8Str(s), 0, PCRE.ANCHORED)
starts_with(s::MaybeSub{<:Str{C}}, r::RegexStr) where {C<:Regex_CSEs} =
    comp_exec(C, r, s, 0, PCRE.ANCHORED)
ends_with(s::AbstractString, r::RegexStr) =
    comp_exec(UTF8CSE, r, UTF8Str(s), 0, PCRE.ENDANCHORED)
ends_with(s::MaybeSub{<:Str{C}}, r::RegexStr) where {C<:Regex_CSEs} =
    comp_exec(C, r, s, 0, PCRE.ENDANCHORED)

function _write_capture(io, ::Type{C}, group::Integer, md) where {C<:CSE}
    T = codeunit(C)
    len = PCRE.sub_length_bynumber(T, md, group)
    buf, out = _allocate(T, len)
    PCRE.sub_copy_bynumber(T, md, group, out, len+1)
    print(io, Str(C, buf))
end

const SUB_CHAR = '\\'
const GROUP_CHAR = 'g'
const KEEP_ESC = [SUB_CHAR, GROUP_CHAR, '0':'9'...]
const LBRACKET = '<'
const RBRACKET = '>'

function Base._replace(io, repl_s::SubstitutionString,
                       str::T, r, re::RegexStr) where {T<:AbstractString}
    tid = Threads.threadid()
    C = cse(T)
    CU = codeunit(C)
    cu_index = codeunit_index(CU)
    md = re.match[tid].match_data[cu_index]
    regex = re.table[cu_index][opt_index(C)]
    @static if VERSION < v"1.3"
        repl = repl_s.string
    else
        repl = Base.unescape_string(repl_s.string, KEEP_ESC)
    end
    pos = 1
    lst = lastindex(repl)
    # This needs to be careful with writes!
    while pos <= lst
        ch = repl[pos]
        if ch == SUB_CHAR
            nxt = nextind(repl, pos)
            nxt > lst && replace_err(repl)
            ch = repl[nxt]
            if ch == SUB_CHAR
                print(io, SUB_CHAR)
                pos = nextind(repl, nxt)
            elseif is_digit(ch)
                group = parse(Int, ch)
                pos = nextind(repl, nxt)
                while pos <= lst && (ch = repl[pos]; is_digit(ch))
                    group = 10 * group + parse(Int, ch)
                    pos = nextind(repl, pos)
                end
                _write_capture(io, C, group, md)
            elseif ch == GROUP_CHAR
                pos = nextind(repl, nxt)
                (pos > lst || repl[pos] != LBRACKET) && replace_err(repl)
                pos = nextind(repl, pos)
                pos > lst && replace_err(repl)
                groupstart = pos
                while repl[pos] != RBRACKET
                    pos = nextind(repl, pos)
                    pos > pos && replace_err(repl)
                end
                #  TODO: avoid this allocation
                groupname = SubString(repl, groupstart, prevind(repl, pos))
                if all(isdigit, groupname)
                    _write_capture(io, C, parse(Int, groupname), md)
                else
                    gn = convert(T, groupname)
                    group = PCRE.substring_number_from_name(CU, regex, gn)
                    group < 0 && replace_err("Group $groupname not found in regex $re")
                    _write_capture(io, C, group, md)
                end
                pos = nextind(repl, pos)
            else
                replace_err(repl)
            end
        else
            print(io, ch)
            pos = nextind(repl, pos)
        end
    end
end

struct RegexStrMatchIterator{T<:AbstractString}
    regex::RegexStr
    string::T
    overlap::Bool
end
RegexStrMatchIterator(r::RegexStr, s::AbstractString) = RegexStrMatchIterator(r, s, false)

compile(itr::RegexStrMatchIterator{T}) where {T} = (compile(cse(T), itr.regex); itr)
eltype(::Type{RegexStrMatchIterator{T}}) where {T} = RegexStrMatch{T}
firstindex(itr::RegexStrMatchIterator) = match(itr.regex, itr.string, 1, UInt32(0))
IteratorSize(::Type{RegexStrMatchIterator{T}}) where {T<:AbstractString} = Base.SizeUnknown()

function iterate(itr::RegexStrMatchIterator, (offset,prevempty)=(1,false))
    opts_nonempty = UInt32(PCRE.ANCHORED | PCRE.NOTEMPTY_ATSTART)
    while true
        mat = match(itr.regex, itr.string, offset, prevempty ? opts_nonempty : UInt32(0))

        if mat === nothing
            if prevempty && offset <= sizeof(itr.string)
                offset = nextind(itr.string, offset)
                prevempty = false
                continue
            else
                break
            end
        else
            if itr.overlap
                offset = isempty(mat.match) ? mat.offset : nextind(itr.string, mat.offset)
            else
                offset = mat.offset + ncodeunits(mat.match)
            end
            return (mat, (offset, isempty(mat.match)))
        end
    end
    nothing
end

eachmatch(re::RegexStr, str::AbstractString; overlap = false) =
    RegexStrMatchIterator(re, str, overlap)

const MS_Str    = MaybeSub{<:Str}
const MS_String = MaybeSub{String}

split(str::MS_Str, splitter::Regex;
      limit::Integer=0, keepempty::Bool=true, keep::Union{Nothing,Bool}=nothing) =
    __split(str, splitter, limit, checkkeep(keepempty, keep, :split), splitarr(str))
split(str::MS_Str, splitter::RegexStr;
      limit::Integer=0, keepempty::Bool=true, keep::Union{Nothing,Bool}=nothing) =
    __split(str, splitter, limit, checkkeep(keepempty, keep, :split), splitarr(str))
split(str::MS_String, splitter::RegexStr;
      limit::Integer=0, keepempty::Bool=true, keep::Union{Nothing,Bool}=nothing) =
    __split(str, splitter, limit, checkkeep(keepempty, keep, :split), splitarr(str))

rsplit(str::MS_Str, splitter::Regex;
      limit::Integer=0, keepempty::Bool=true, keep::Union{Nothing,Bool}=nothing) =
    __rsplit(str, splitter, limit, checkkeep(keepempty, keep, :rsplit), splitarr(str))
rsplit(str::MS_Str, splitter::RegexStr;
      limit::Integer=0, keepempty::Bool=true, keep::Union{Nothing,Bool}=nothing) =
    __rsplit(str, splitter, limit, checkkeep(keepempty, keep, :rsplit), splitarr(str))
rsplit(str::MS_String, splitter::RegexStr;
      limit::Integer=0, keepempty::Bool=true, keep::Union{Nothing,Bool}=nothing) =
    __rsplit(str, splitter, limit, checkkeep(keepempty, keep, :rsplit), splitarr(str))

replace(str::MS_Str, pat_repl::Pair{Regex}; count::Integer=typemax(Int)) =
    __replace(str, pat_repl; count=count)
replace(str::MS_Str, pat_repl::Pair{RegexStr}; count::Integer=typemax(Int)) =
    __replace(str, pat_repl; count=count)
replace(str::String, pat_repl::Pair{RegexStr}; count::Integer=typemax(Int)) =
    __replace(str, pat_repl; count=count)
replace(str::SubString{String}, pat_repl::Pair{RegexStr}; count::Integer=typemax(Int)) =
    __replace(str, pat_repl; count=count)

## comparison ##

==(a::RegexTypes, b::RegexTypes) =
    a.pattern == b.pattern &&
    a.compile_options == b.compile_options &&
    a.match_options == b.match_options

## hash ##
hash(r::RegexStr, h::UInt) =
    hash(r.match_options,
         hash(r.compile_options,
              hash(r.pattern, h + UInt === UInt64 ? 0x67e195eb8555e72d : 0xe32373e4)))

_occurs_in(r::RegexTypes, s::AbstractString, off::Integer) =
    comp_exec(UTF8CSE, r, UTF8Str(s), off, 0)
_occurs_in(r::RegexTypes, s::MaybeSub{<:Str{C}}, off::Integer) where {C<:Regex_CSEs} =
    comp_exec(C, r, s, off, 0)

occurs_in(needle::RegexStr, hay::AbstractString; off::Integer=0) = _occurs_in(needle, hay, off)
occurs_in(needle::Regex, hay::MaybeSub{<:Str}; off::Integer=0)   = _occurs_in(needle, hay, off)

findnext(a::RegexStr, b::AbstractString, i) = StrBase.nothing_sentinel(find(Fwd, a, b, i))

@static isdefined(Base, :isnothing) || (isnothing(r) = (r === nothing))

# Copied from julia/base/regex.jl
# Regex really should be moved to stdlib, and an AbstractRegex or AbstractPattern type added
#
"""
    findall(
        pattern::Union{AbstractString,Regex,RegexStr},
        string::AbstractString;
        overlap::Bool = false,
    )

Return a `Vector{UnitRange{Int}}` of all the matches for `pattern` in `string`.
Each element of the returned vector is a range of indices where the
matching sequence is found, like the return value of [`findnext`](@ref).

If `overlap=true`, the matching sequences are allowed to overlap indices in the
original string, otherwise they must be from disjoint character ranges.
"""
function findall(t::RegexStr, s::AbstractString; overlap::Bool=false)
    found = UnitRange{Int}[]
    i, e = firstindex(s), lastindex(s)
    while true
        r = findnext(t, s, i)
        isnothing(r) && break
        push!(found, r)
        j = overlap || isempty(r) ? first(r) : last(r)
        j > e && break
        @inbounds i = nextind(s, j)
    end
    return found
end

"""
    count(
        pattern::Union{AbstractString,Regex,RegexStr},
        string::AbstractString;
        overlap::Bool = false,
    )

Return the number of matches for `pattern` in `string`. This is equivalent to
calling `length(findall(pattern, string))` but more efficient.

If `overlap=true`, the matching sequences are allowed to overlap indices in the
original string, otherwise they must be from disjoint character ranges.
"""
function count(t::RegexStr, s::AbstractString; overlap::Bool=false)
    n = 0
    i, e = firstindex(s), lastindex(s)
    while true
        r = findnext(t, s, i)
        isnothing(r)  && break
        n += 1
        j = overlap || isempty(r) ? first(r) : last(r)
        j > e && break
        @inbounds i = nextind(s, j)
    end
    return n
end

@api freeze

end # module StrRegex


text:
module ManualMemory

mutable struct MemoryBuffer{N,T}
  data::NTuple{N,T}
  @inline function MemoryBuffer{N,T}(::UndefInitializer) where {N,T}
    @assert Base.allocatedinline(T)
    new{N,T}()
  end
end
@inline Base.unsafe_convert(::Type{Ptr{T}}, m::MemoryBuffer) where {T} = Ptr{T}(pointer_from_objref(m))
@inline Base.pointer(m::MemoryBuffer{N,T}) where {N,T} = Ptr{T}(pointer_from_objref(m))

"""
    PseudoPtr(data, position=firstindex(data))

Provides a convenient wrapper that functions like `pointer(data)` for instances where `data`
cannot produce a viable pointer.
"""
struct PseudoPtr{T,D} <: Ref{T}
  data::D
  position::Int

  PseudoPtr(data::D, position) where {D} = new{eltype(D),D}(data, position)
  PseudoPtr(data) = PseudoPtr(data, firstindex(data))
end
Base.:(+)(x::PseudoPtr, y::Int) = PseudoPtr(getfield(x, :data), getfield(x, :position) + y)
Base.:(+)(x::Int, y::PseudoPtr) = y + x

@inline load(x::PseudoPtr) = @inbounds(getindex(getfield(x, :data), getfield(x, :position)))
@generated function load(p::Ptr{T}) where {T}
  if Base.allocatedinline(T)
    Expr(:block, Expr(:meta,:inline), :(unsafe_load(p)))
  else
    Expr(:block, Expr(:meta,:inline), :(ccall(:jl_value_ptr, Ref{$T}, (Ptr{Cvoid},), unsafe_load(Base.unsafe_convert(Ptr{Ptr{Cvoid}}, p)))))
  end
end
@inline load(p::Ptr{UInt}, ::Type{T}) where {T} = load(p, T, 0)[2]

@inline function store!(x::PseudoPtr, val)
    @inbounds(setindex!(getfield(x, :data), val, getfield(x, :position)))
end
@generated function store!(p::Ptr{T}, v::T) where {T}
  if Base.allocatedinline(T)
    Expr(:block, Expr(:meta,:inline), :(unsafe_store!(p, v); return nothing))
  else
    Expr(:block, Expr(:meta,:inline), :(unsafe_store!(Base.unsafe_convert(Ptr{Ptr{Cvoid}}, p), Base.pointer_from_objref(v)); return nothing))
  end
end
@generated offsetsize(::Type{T}) where {T} = Base.allocatedinline(T) ? sizeof(T) : sizeof(Int)

@inline store!(p::Ptr{T}, v) where {T} = store!(p, convert(T, v))

mutable struct Reference{T}
  data::T
  Reference{T}() where {T} = new{T}()
  Reference(x) = new{typeof(x)}(x)
end
@inline Base.pointer(r::Reference{T}) where {T} = Ptr{T}(pointer_from_objref(r))
@inline load(p::Ptr{Reference{T}}) where {T} = getfield(ccall(:jl_value_ptr, Ref{Reference{T}}, (Ptr{Cvoid},), unsafe_load(Base.unsafe_convert(Ptr{Ptr{Cvoid}}, p))), :data)
@inline dereference(r::Reference) = getfield(r, :data)
@inline dereference(x) = x

"""
    LazyPreserve(x, ptrcall=nothing)

Used by [`preserve`](@ref) to identify arguments that will be unwrapped with
[`preserve_buffer`](@ref), which is in turn converted to a pointer. If `ptrcall` is
specified conversion to a pointer occurs through a call equivalent to
`ptrcall(preserve_buffer(x), x)`. Otherwise, a call equivalent to
`pointer(preserve_buffer(x))` occurs.
"""
struct LazyPreserve{A,P}
  arg::A
  ptrcall::P
end
LazyPreserve(arg) = LazyPreserve(arg, nothing)
(p::LazyPreserve)(x) = p.ptrcall(x, p.arg)
(p::LazyPreserve{A,Nothing})(x) where {A} = pointer(x)

"""
    preserve_buffer(x)

For structs wrapping arrays, using `GC.@preserve` can trigger heap allocations.
`preserve_buffer` attempts to extract the heap-allocated part. Isolating it by itself
will often allow the heap allocations to be elided. For example:

```julia
julia> using StaticArrays, BenchmarkTools
julia> # Needed until a release is made featuring https://github.com/JuliaArrays/StaticArrays.jl/commit/a0179213b741c0feebd2fc6a1101a7358a90caed
       Base.elsize(::Type{<:MArray{S,T}}) where {S,T} = sizeof(T)
julia> @noinline foo(A) = unsafe_load(A,1)
foo (generic function with 1 method)
julia> function alloc_test_1()
           A = view(MMatrix{8,8,Float64}(undef), 2:5, 3:7)
           A[begin] = 4
           GC.@preserve A foo(pointer(A))
       end
alloc_test_1 (generic function with 1 method)
julia> function alloc_test_2()
           A = view(MMatrix{8,8,Float64}(undef), 2:5, 3:7)
           A[begin] = 4
           pb = parent(A) # or `LoopVectorization.preserve_buffer(A)`; `perserve_buffer(::SubArray)` calls `parent`
           GC.@preserve pb foo(pointer(A))
       end
alloc_test_2 (generic function with 1 method)
julia> @benchmark alloc_test_1()
BenchmarkTools.Trial:
  memory estimate:  544 bytes
  allocs estimate:  1
  --------------
  minimum time:     17.227 ns (0.00% GC)
  median time:      21.352 ns (0.00% GC)
  mean time:        26.151 ns (13.33% GC)
  maximum time:     571.130 ns (78.53% GC)
  --------------
  samples:          10000
  evals/sample:     998
julia> @benchmark alloc_test_2()
BenchmarkTools.Trial:
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     3.275 ns (0.00% GC)
  median time:      3.493 ns (0.00% GC)
  mean time:        3.491 ns (0.00% GC)
  maximum time:     4.998 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1000
```
"""
@inline preserve_buffer(x::LazyPreserve) = preserve_buffer(x.arg)
@inline preserve_buffer(x) = x
@inline preserve_buffer(A::AbstractArray) = _preserve_buffer(A, parent(A))
@inline _preserve_buffer(a::A, p::P) where {A,P<:AbstractArray} = _preserve_buffer(p, parent(p))
@inline _preserve_buffer(a::A, p::A) where {A<:AbstractArray} = a
@inline _preserve_buffer(a::A, p::P) where {A,P} = p

function load_aggregate(::Type{T}, offset::Int) where {T}
  numfields = fieldcount(T)
  call = (T <: Tuple) ? Expr(:tuple) : Expr(:new, T)
  for f ∈ 1:numfields
    TF = fieldtype(T, f)
    if Base.issingletontype(TF)
      push!(call.args, TF.instance)
    elseif (fieldcount(TF) ≡ 0) || (TF <: Reference)
      ptr = :(p + (offset + $offset))
      ptr = TF === UInt ? ptr : :(reinterpret(Ptr{$TF}, $ptr))
      push!(call.args, :(load($ptr)))
      offset += offsetsize(TF)
    else
      arg, offset = load_aggregate(TF, offset)
      push!(call.args, arg)
    end
  end
  return call, offset
end
@generated function load(p::Ptr{UInt}, ::Type{T}, offset::Int) where {T}
  if Base.issingletontype(T)
    call = Expr(:tuple, :offset, T.instance)
  elseif (fieldcount(T) ≡ 0) || (T <: Reference)
    ptr = :(p + offset)
    ptr = T === UInt ? ptr : :(reinterpret(Ptr{$T}, $ptr))
    call = :(((offset + $(offsetsize(T)), load($ptr))))
  else
    call, off = load_aggregate(T, 0)
    call = Expr(:tuple, :(offset + $off), call)
  end
  Expr(:block, Expr(:meta,:inline), call)
end

function store_aggregate!(q::Expr, sym, ::Type{T}, offset::Int) where {T}
  gf = GlobalRef(Core,:getfield)
  for f ∈ 1:fieldcount(T)
    TF = fieldtype(T, f)
    Base.issingletontype(TF) && continue
    gfcall = Expr(:call, gf, sym, f)
    if (fieldcount(TF) ≡ 0) || (TF <: Reference)
      ptr = :(p + (offset + $offset))
      ptr = TF === UInt ? ptr : :(reinterpret(Ptr{$TF}, $ptr))
      push!(q.args, :(store!($ptr, $gfcall)))
      offset += offsetsize(TF)
    else
      newsym = gensym(sym)
      push!(q.args, Expr(:(=), newsym, gfcall))
      offset = store_aggregate!(q, newsym, TF, offset)
    end
  end
  return offset
end
@generated function store!(p::Ptr{UInt}, x::T, offset::Int) where {T}
  Base.issingletontype(T) && return :offset
  body = Expr(:block, Expr(:meta,:inline))
  if (fieldcount(T) ≡ 0) || (T <: Reference)
    ptr = :(p + offset)
    ptr = T === UInt ? ptr : :(reinterpret(Ptr{$T}, $ptr))
    push!(body.args, :(store!($ptr, x)))
    off = offsetsize(T)
  else
    off = store_aggregate!(body, :x, T, 0)
  end
  push!(body.args, Expr(:call, +, :offset, off))
  return body
end

"""
    preserve(op, args...; kwargs...)

Searches through `args` and `kwargs` for instances of [`LazyPreserve`](@ref), which are
unwrapped using [`preserve_buffer`](@ref) and preserved from garbage collection
(`GC.@preserve`). the resulting buffers are converted to pointers and passed in order to `op`.

# Examples

```julia
julia> using ManualMemory: store!, preserve, LazyPreserve

julia> x = [0 0; 0 0];

julia> preserve(store!, LazyPreserve(x), 1)

julia> x[1]
1

```
"""
preserve(op, args...; kwargs...) = _preserve(op, args, kwargs.data)
@generated function _preserve(op, args::A, kwargs::NamedTuple{syms,K}) where {A,syms,K}
  _preserve_expr(A, syms, K)
end
function _preserve_expr(::Type{A}, syms::Tuple{Vararg{Symbol}}, ::Type{K}) where {A,K}
  body = Expr(:block, Expr(:meta,:inline))
  call = Expr(:call, :op)
  pres = :(GC.@preserve)
  @inbounds for i in 1:length(A.parameters)
    arg_i = _unwrap_preserve(body, pres, :(getfield(args, $i)), A.parameters[i])
    push!(call.args, arg_i)
  end
  if length(syms) > 0
    kwargs = Expr(:parameters)
    @inbounds for i in 1:length(syms)
      arg_i = _unwrap_preserve(body, pres, :(getfield(kwargs, $i)), K.parameters[i])
      push!(call.args, Expr(:kw, syms[i], arg_i))
    end
    push!(call.args, kwargs)
  end
  push!(pres.args, call)
  push!(body.args, pres)
  return body
end
function _unwrap_preserve(body::Expr, pres::Expr, argexpr::Expr, argtype::Type)
  if argtype <: LazyPreserve
    bufsym = gensym()
    push!(body.args, Expr(:(=), bufsym, Expr(:call, :preserve_buffer, argexpr)))
    push!(pres.args, bufsym)
    return :($argexpr($bufsym))
  else
    return argexpr
  end
end

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "ICU"
version = v"69.1"

# Collection of sources required to build ICU
sources = [
    ArchiveSource("https://github.com/unicode-org/icu/releases/download/release-$(version.major)-$(version.minor)/icu4c-$(version.major)_$(version.minor)-src.tgz",
                  "4cba7b7acd1d3c42c44bb0c14be6637098c7faf2b330ce876bc5f3b915d09745"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/icu/

# Apply patch to link `libicudata` against the default standard libraries
# to avoid toolchain weirdness when you have a dynamic library that has
# _no_ dependencies (not even `libc`).  See this bug report for more:
# https://bugs.debian.org/cgi-bin/bugreport.cgi?bug=653457
atomic_patch -p1 $WORKSPACE/srcdir/patches/yes_stdlibs.patch

# Do the native build
(
    cp -r source/ native_build/
    cd native_build
    CC="${CC_BUILD}"
    CXX="${CXX_BUILD}"
    AR="${AR_BUILD}"
    LD="${LD_BUILD}"
    RANLIB="${RANLIB_BUILD}"

    # See https://git.alpinelinux.org/aports/tree/main/icu/APKBUILD?id=334ebffde9dec34becdd628ad56007699e98ea81
    update_configure_scripts
    sed -i -e 's,DU_HAVE_STRTOD_L=1,DU_HAVE_STRTOD_L=0,' configure.ac
    sed -i -e 's,DU_HAVE_STRTOD_L=1,DU_HAVE_STRTOD_L=0,' configure
    for x in ARFLAGS CFLAGS CPPFLAGS CXXFLAGS FFLAGS LDFLAGS; do
        sed -i -e "/^${x} =.*/s:@${x}@::" "config/Makefile.inc.in"
    done

    ./configure --prefix=$prefix --build=${MACHTYPE} \
        ac_cv_prog_ac_ct_AR=${AR} \
        ac_cv_prog_ac_ct_RANLIB=${RANLIB}
    make -j${nproc}
)

# Do the cross build
cd source/

if [[ "${target}" == *-apple-* ]]; then
    # Do not append `-c` flag to ar, which isn't supported by LLVM's ar
    atomic_patch -p1 $WORKSPACE/srcdir/patches/argflags-no--c.patch
    export LDFLAGS="-headerpad_max_install_names"
fi

update_configure_scripts
./configure --prefix=$prefix --host=$target \
    --with-cross-build="/workspace/srcdir/icu/native_build"
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(; experimental=true))

# The products that we will ensure are always built
products = [
    LibraryProduct(["libicudata", "icudt$(version.major)"], :libicudata; dont_dlopen=true),
    LibraryProduct(["libicui18n", "icuin$(version.major)"], :libicui18n),
    LibraryProduct(["libicuio", "icuio$(version.major)"], :libicuio),
    LibraryProduct(["libicutest", "icutest$(version.major)"], :libicutest),
    LibraryProduct(["libicutu", "icutu$(version.major)"], :libicutu),
    LibraryProduct(["libicuuc", "icuuc$(version.major)"], :libicuuc),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"7", julia_compat="1.6")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# # Eady turbulence example
#
# In this example, we initialize a random velocity field and observe its viscous,
# turbulent decay in a two-dimensional domain. This example demonstrates:
#
#   * How to use a tuple of turbulence closures
#   * How to use hyperdiffusivity
#   * How to implement background velocity and tracer distributions
#   * How to use `ComputedField`s for output
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add Oceananigans, JLD2, Plots"
# ```

# ## The Eady problem
#
# The "Eady problem" simulates the baroclinic instability problem proposed by Eric Eady in
# the classic paper
# ["Long waves and cyclone waves," Tellus (1949)](https://doi.org/10.3402/tellusa.v1i3.8507).
# The Eady problem is a simple, canonical model for the generation of mid-latitude
# atmospheric storms and the ocean eddies that permeate the world sea.
#
# In the Eady problem, baroclinic motion and turublence is generated by the interaction
# between turbulent motions and a stationary, geostrophically-balanced basic state that
# is unstable to baroclinic instability. In this example, the baroclinic generation of
# turbulence due to extraction of energy from the geostrophic basic state
# is balanced by a bottom boundary condition that extracts momentum from turbulent motions
# and serves as a crude model for the drag associated with an unresolved and small-scale
# turbulent bottom boundary layer.
#
# ### The geostrophic basic state
#
# The geostrophic basic state in the Eady problem is represented by the streamfunction,
#
# ```math
# ψ(y, z) = - α y (z + L_z) \, ,
# ```
#
# where ``α`` is the geostrophic shear and ``L_z`` is the depth of the domain.
# The background buoyancy includes both the geostrophic flow component,
# ``f ∂_z ψ``, where ``f`` is the Coriolis parameter, and a background stable stratification
# component, ``N^2 z``, where ``N`` is the buoyancy frequency:
#
# ```math
# B(y, z) = f ∂_z ψ + N^2 z = - α f y + N^2 z \, .
# ```
#
# The background velocity field is related to the geostrophic streamfunction via
# ``U = - ∂_y ψ`` such that
#
# ```math
# U(z) = α (z + L_z) \, .
# ```
#
# ### Boundary conditions
#
# All fields are periodic in the horizontal directions.
# We use "insulating", or zero-flux boundary conditions on the buoyancy perturbation
# at the top and bottom. We thus implicitly assume that the background vertical density
# gradient, ``N^2 z``, is maintained by a process external to our simulation.
# We use free-slip, or zero-flux boundary conditions on ``u`` and ``v`` at the surface
# where ``z=0``. At the bottom, we impose a momentum flux that extracts momentum and
# energy from the flow.
#
# #### Bottom boundary condition: quadratic bottom drag
#
# We model the effects of a turbulent bottom boundary layer on the eddy momentum budget
# with quadratic bottom drag. A quadratic cottom drag is introduced by imposing a vertical flux
# of horizontal momentum that removes momentum from the layer immediately above: in other words,
# the flux is negative (downwards) when the velocity at the bottom boundary is positive, and
# positive (upwards) with the velocity at the bottom boundary is negative.
# This drag term is "quadratic" because the rate at which momentum is removed is proportional
# to ``\boldsymbol{u} |\boldsymbol{u}|``, where ``\boldsymbol{u} = u \boldsymbol{\hat{x}} +
# v \boldsymbol{\hat{y}}`` is the horizontal velocity.
#
# The ``x``-component of the quadratic bottom drag is thus
#
# ```math
# \tau_{xz}(z=L_z) = - c^D u \sqrt{u^2 + v^2} \, ,
# ```
#
# while the ``y``-component is
#
# ```math
# \tau_{yz}(z=L_z) = - c^D v \sqrt{u^2 + v^2} \, ,
# ```
#
# where ``c^D`` is a dimensionless drag coefficient and ``\tau_{xz}(z=L_z)`` and ``\tau_{yz}(z=L_z)``
# denote the flux of ``u`` and ``v`` momentum at ``z = L_z``, the bottom of the domain.
#
# ### Vertical and horizontal viscosity and diffusivity
#
# Vertical and horizontal viscosities and diffusivities are required
# to stabilize the Eady problem and can be idealized as modeling the effect of
# turbulent mixing below the grid scale. For both tracers and velocities we use
# a Laplacian vertical diffusivity ``κ_z ∂_z^2 c`` and a horizontal
# hyperdiffusivity ``ϰ_h (∂_x^4 + ∂_y^4) c``.
#
# ### Eady problem summary and parameters
#
# To summarize, the Eady problem parameters along with the values we use in this example are
#
# | Parameter name | Description | Value | Units |
# |:--------------:|:-----------:|:-----:|:-----:|
# | ``f``          | Coriolis parameter | ``10^{-4}`` | ``\mathrm{s^{-1}}`` |
# | ``N``          | Buoyancy frequency (square root of ``\partial_z B``) | ``10^{-3}`` | ``\mathrm{s^{-1}}`` |
# | ``\alpha``     | Background vertical shear ``\partial_z U`` | ``10^{-3}`` | ``\mathrm{s^{-1}}`` |
# | ``c^D``        | Bottom quadratic drag coefficient | ``10^{-4}`` | none |
# | ``κ_z``        | Laplacian vertical diffusivity | ``10^{-2}`` | ``\mathrm{m^2 s^{-1}}`` |
# | ``ϰ_h``        | Biharmonic horizontal diffusivity | ``10^{-2} \times \Delta x^4 / \mathrm{day}`` | ``\mathrm{m^4 s^{-1}}`` |
#
# We start off by importing `Oceananigans`, `Printf`, and some convenient constants
# for specifying dimensional units:

using Printf
using Oceananigans
using Oceananigans.Units: hours, day, days

# ## The grid
#
# We use a three-dimensional grid with a depth of 4000 m and a
# horizontal extent of 1000 km, appropriate for mesoscale ocean dynamics
# with characteristic scales of 50-200 km.

grid = RectilinearGrid(size=(48, 48, 16), x=(0, 1e6), y=(0, 1e6), z=(-4e3, 0))

# ## Rotation
#
# The classical Eady problem is posed on an ``f``-plane. We use a Coriolis parameter
# typical to mid-latitudes on Earth,

coriolis = FPlane(f=1e-4) # [s⁻¹]

# ## The background flow
#
# We build a `NamedTuple` of parameters that describe the background flow,

basic_state_parameters = ( α = 10 * coriolis.f, # s⁻¹, geostrophic shear
                           f = coriolis.f,      # s⁻¹, Coriolis parameter
                           N = 1e-3,            # s⁻¹, buoyancy frequency
                          Lz = grid.Lz)         # m, ocean depth

# and then construct the background fields ``U`` and ``B``

## Background fields are defined via functions of x, y, z, t, and optional parameters
U(x, y, z, t, p) = + p.α * (z + p.Lz)
B(x, y, z, t, p) = - p.α * p.f * y + p.N^2 * z

U_field = BackgroundField(U, parameters=basic_state_parameters)
B_field = BackgroundField(B, parameters=basic_state_parameters)

# ## Boundary conditions
#
# The boundary conditions prescribe a quadratic drag at the bottom as a flux
# condition.

cᴰ = 1e-4 # quadratic drag coefficient

@inline drag_u(x, y, t, u, v, cᴰ) = - cᴰ * u * sqrt(u^2 + v^2)
@inline drag_v(x, y, t, u, v, cᴰ) = - cᴰ * v * sqrt(u^2 + v^2)

drag_bc_u = FluxBoundaryCondition(drag_u, field_dependencies=(:u, :v), parameters=cᴰ)
drag_bc_v = FluxBoundaryCondition(drag_v, field_dependencies=(:u, :v), parameters=cᴰ)

u_bcs = FieldBoundaryConditions(bottom = drag_bc_u)
v_bcs = FieldBoundaryConditions(bottom = drag_bc_v)

# ## Turbulence closures
#
# We use a horizontal hyperdiffusivity and a Laplacian vertical diffusivity
# to dissipate energy in the Eady problem.
# To use both of these closures at the same time, we set the keyword argument
# `closure` to a tuple of two closures.

κ₂z = 1e-2 # [m² s⁻¹] Laplacian vertical viscosity and diffusivity
κ₄h = 1e-1 / day * grid.Δxᶜᵃᵃ^4 # [m⁴ s⁻¹] horizontal hyperviscosity and hyperdiffusivity

Laplacian_vertical_diffusivity = AnisotropicDiffusivity(νh=0, κh=0, νz=κ₂z, κz=κ₂z)
biharmonic_horizontal_diffusivity = AnisotropicBiharmonicDiffusivity(νh=κ₄h, κh=κ₄h)

# ## Model instantiation
#
# We instantiate the model with the fifth-order WENO advection scheme, a 3rd order
# Runge-Kutta time-stepping scheme, and a `BuoyancyTracer`.

model = NonhydrostaticModel(
                   grid = grid,
              advection = WENO5(),
            timestepper = :RungeKutta3,
               coriolis = coriolis,
                tracers = :b,
               buoyancy = BuoyancyTracer(),
      background_fields = (b=B_field, u=U_field),
                closure = (Laplacian_vertical_diffusivity, biharmonic_horizontal_diffusivity),
    boundary_conditions = (u=u_bcs, v=v_bcs)
)

# ## Initial conditions
#
# We seed our initial conditions with random noise stimulate the growth of
# baroclinic instability.

## A noise function, damped at the top and bottom
Ξ(z) = randn() * z/grid.Lz * (z/grid.Lz + 1)

## Scales for the initial velocity and buoyancy
Ũ = 1e-1 * basic_state_parameters.α * grid.Lz
B̃ = 1e-2 * basic_state_parameters.α * coriolis.f

uᵢ(x, y, z) = Ũ * Ξ(z)
vᵢ(x, y, z) = Ũ * Ξ(z)
bᵢ(x, y, z) = B̃ * Ξ(z)

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

# We subtract off any residual mean velocity to avoid exciting domain-scale
# inertial oscillations. We use a `sum` over the entire `parent` arrays or data
# to ensure this operation is efficient on the GPU (set `architecture = GPU()`
# in `NonhydrostaticModel` constructor to run this problem on the GPU if one
# is available).

ū = sum(model.velocities.u.data.parent) / (grid.Nx * grid.Ny * grid.Nz)
v̄ = sum(model.velocities.v.data.parent) / (grid.Nx * grid.Ny * grid.Nz)

model.velocities.u.data.parent .-= ū
model.velocities.v.data.parent .-= v̄
nothing # hide

# ## Simulation set-up
#
# We set up a simulation that runs for 10 days with a `JLD2OutputWriter` that saves the
# vertical vorticity and divergence every 2 hours. We limit the time-step to
# the maximum allowable due either to diffusion, internal waves, or advection by the background flow.

## Calculate absolute limit on time-step using diffusivities and
## background velocity.
Ū = basic_state_parameters.α * grid.Lz

max_Δt = min(grid.Δxᶜᵃᵃ / Ū, grid.Δxᶜᵃᵃ^4 / κ₄h, grid.Δxᶜᵃᵃ^2 / κ₂z, 1 / basic_state_parameters.N)

simulation = Simulation(model, Δt = max_Δt, stop_time = 8days)

# ### The `TimeStepWizard`
#
# The TimeStepWizard manages the time-step adaptively, keeping the
# Courant-Freidrichs-Lewy (CFL) number close to `1.0` while ensuring
# the time-step does not increase beyond the maximum allowable value
# for numerical stability given the specified background flow, internal wave
# time scales, and diffusion time scales.

wizard = TimeStepWizard(cfl=0.85, max_change=1.1, max_Δt=max_Δt)

simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))

# ### A progress messenger
#
# We add a callback that prints out a helpful progress message while the simulation runs.

start_time = time_ns()

progress(sim) = @printf("i: % 6d, sim time: % 10s, wall time: % 10s, Δt: % 10s, CFL: %.2e\n",
                        sim.model.clock.iteration,
                        prettytime(sim.model.clock.time),
                        prettytime(1e-9 * (time_ns() - start_time)),
                        prettytime(sim.Δt),
                        AdvectiveCFL(sim.Δt)(sim.model))

simulation.callbacks[:progress] = Callback(progress, IterationInterval(10))

# ### Output
#
# To visualize the baroclinic turbulence ensuing in the Eady problem,
# we use `ComputedField`s to diagnose and output vertical vorticity and divergence.
# Note that `ComputedField`s take "AbstractOperations" on `Field`s as input:

u, v, w = model.velocities # unpack velocity `Field`s

## Vertical vorticity [s⁻¹]
ζ = ComputedField(∂x(v) - ∂y(u))

## Horizontal divergence, or ∂x(u) + ∂y(v) [s⁻¹]
δ = ComputedField(-∂z(w))

# With the vertical vorticity, `ζ`, and the horizontal divergence, `δ` in hand,
# we create a `JLD2OutputWriter` that saves `ζ` and `δ` and add them to
# `simulation`.

simulation.output_writers[:fields] = JLD2OutputWriter(model, (ζ=ζ, δ=δ),
                                                      schedule = TimeInterval(4hours),
                                                        prefix = "eady_turbulence",
                                                         force = true)
nothing # hide

# All that's left is to press the big red button:

run!(simulation)

# ## Visualizing Eady turbulence
#
# We animate the results by opening the JLD2 file, extracting data for
# the iterations we ended up saving at, and ploting slices of the saved
# fields. We prepare for animating the flow by creating coordinate arrays,
# opening the file, building a vector of the iterations that we saved
# data at, and defining a function for computing colorbar limits:

using JLD2, Plots

## Coordinate arrays
xζ, yζ, zζ = nodes(ζ)
xδ, yδ, zδ = nodes(δ)

## Open the file with our data
file = jldopen(simulation.output_writers[:fields].filepath)

## Extract a vector of iterations
iterations = parse.(Int, keys(file["timeseries/t"]))

# This utility is handy for calculating nice contour intervals:

function nice_divergent_levels(c, clim, nlevels=31)
    levels = range(-clim, stop=clim, length=nlevels)
    cmax = maximum(abs, c)
    clim < cmax && (levels = vcat([-cmax], levels, [cmax]))
    return levels
end
nothing # hide

# Now we're ready to animate.

@info "Making an animation from saved data..."

anim = @animate for (i, iter) in enumerate(iterations)

    ## Load 3D fields from file
    t = file["timeseries/t/$iter"]
    R_snapshot = file["timeseries/ζ/$iter"] ./ coriolis.f
    δ_snapshot = file["timeseries/δ/$iter"]

    surface_R = R_snapshot[:, :, grid.Nz]
    surface_δ = δ_snapshot[:, :, grid.Nz]

    slice_R = R_snapshot[:, 1, :]
    slice_δ = δ_snapshot[:, 1, :]

    Rlim = 0.5 * maximum(abs, R_snapshot) + 1e-9
    δlim = 0.5 * maximum(abs, δ_snapshot) + 1e-9

    Rlevels = nice_divergent_levels(R_snapshot, Rlim)
    δlevels = nice_divergent_levels(δ_snapshot, δlim)

    @info @sprintf("Drawing frame %d from iteration %d: max(ζ̃ / f) = %.3f \n",
                   i, iter, maximum(abs, surface_R))

    xy_kwargs = (xlims = (0, grid.Lx), ylims = (0, grid.Lx),
                 xlabel = "x (m)", ylabel = "y (m)",
                 aspectratio = 1,
                 linewidth = 0, color = :balance, legend = false)

    xz_kwargs = (xlims = (0, grid.Lx), ylims = (-grid.Lz, 0),
                 xlabel = "x (m)", ylabel = "z (m)",
                 aspectratio = grid.Lx / grid.Lz * 0.5,
                 linewidth = 0, color = :balance, legend = false)

    R_xy = contourf(xζ, yζ, surface_R'; clims=(-Rlim, Rlim), levels=Rlevels, xy_kwargs...)
    δ_xy = contourf(xδ, yδ, surface_δ'; clims=(-δlim, δlim), levels=δlevels, xy_kwargs...)
    R_xz = contourf(xζ, zζ, slice_R'; clims=(-Rlim, Rlim), levels=Rlevels, xz_kwargs...)
    δ_xz = contourf(xδ, zδ, slice_δ'; clims=(-δlim, δlim), levels=δlevels, xz_kwargs...)

    plot(R_xy, δ_xy, R_xz, δ_xz,
           size = (1000, 800),
           link = :x,
         layout = Plots.grid(2, 2, heights=[0.5, 0.5, 0.2, 0.2]),
          title = [@sprintf("ζ(t=%s) / f", prettytime(t)) @sprintf("δ(t=%s) (s⁻¹)", prettytime(t)) "" ""])

    iter == iterations[end] && close(file)
end

mp4(anim, "eady_turbulence.mp4", fps = 8) # hide