text:
getdim(A::AbstractArray, d::Integer, i) = getindex(A, Base._setindex(i, d, axes(A)...)...)

function compute_scores(model, x; chunksize = 10000)
    x_obs = ndims(x)
    n = size(x, x_obs)
    scores = zeros(eltype(x), n)

    for inds in partition(1:n, chunksize)
        scores[inds] .= model(getdim(x, x_obs, inds))[:]
    end
    return scores
end

# batch provider
function batch_provider(d::Data, t::Train, x, y)
    Random.seed!(t.seed)
    if length(y) == d.batchsize || d.batchsize == 0
        return ((x, y) for iter in 1:t.iters)
    else

        neg = findall(vec(y) .== 0)
        pos = findall(vec(y) .== 1)

        n_neg = d.batchsize ÷ 2
        n_pos = d.batchsize - n_neg

        x_obs = ndims(x)
        y_obs = ndims(y)

        function make_batch()
            inds = vcat(
                sample(neg, n_neg; replace = length(neg) < n_neg),
                sample(pos, n_pos; replace = length(pos) < n_pos),
            )
            shuffle!(inds)
            return (getdim(x, x_obs, inds), getdim(y, y_obs, inds))
        end
        return (make_batch() for iter in 1:t.iters)
    end
end

# train!
function custom_train!(loss, ps, data, opt; cb = (args...) -> ())
    ps = Flux.Zygote.Params(ps)
    cb = Flux.Optimise.runall(cb)
    local loss_val

    for d in data
        try
            gs = Flux.Zygote.gradient(ps) do
                loss_val = loss(Flux.Optimise.batchmemaybe(d)...)
                return loss_val
            end
            Flux.Optimise.update!(opt, ps, gs)
            cb(loss_val, Flux.Optimise.batchmemaybe(d)...)
        catch ex
            if ex isa Flux.Optimise.StopException
                break
            else
                rethrow(ex)
            end
        end
    end
end

# ------------------------------------------------------------------------------------------
# callback and savemodel
# ------------------------------------------------------------------------------------------
Base.@kwdef mutable struct CallBack
    iters::Int
    title::String = "Training:"
    bar::Progress = Progress(iters, 5, title)
    saveat::Int = 1000
    savefunc::Function = (args...) -> nothing
    counter::Int = 0
    loss = Float32[]
end

CallBack(iters; kwargs...) = CallBack(; iters, kwargs...)

function (c::CallBack)(loss_val, x, y)
    c.counter += 1
    push!(c.loss, loss_val)
    mod(c.counter, c.saveat) == 0 && c.savefunc(c, x, y)

    tm = round((c.bar.tlast - c.bar.tfirst)/c.bar.counter; sigdigits = 2)

    next!(c.bar; showvalues = vcat(
        ("Iteration", string(c.counter, "/", c.iters)),
        ("Time/iteration", string(tm, "s")),
        ("Loss0", c.loss[1]),
        ("Loss", c.loss[end]),
    ))
    return
end

# savemodel
function savemodel(
    c::CallBack,
    d::Data,
    t::Train,
    m::Model,
    model,
    x,
    y,
)
    file = modeldir(string(d), string(t), string(m, ".bson"))
    dict = isfile(file) ? BSON.load(file) : Dict{Symbol, Any}()

    s, y = compute_scores(model, x), vec(y) |> cpu

    get!(dict, :dataset_settings, todict(d))
    get!(dict, :train_settings, todict(t))
    get!(dict, :model_settings, todict(m))
    get!(dict, :iterations, Dict{Int, Any}())

    # upadte loss
    get!(dict, :loss, c.loss)
    dict[:loss] = c.loss

    # update time
    tm = (c.bar.tlast - c.bar.tfirst)
    get!(dict, :time_per_iter, 1.)
    get!(dict, :time_total, 1.)
    dict[:time_per_iter] = tm/c.bar.counter
    dict[:time_total] = tm

    # add new iter
    dict[:iterations][c.counter] = Dict(
        :minibatch => Dict(:targets => y, :scores => s),
        :model => deepcopy(cpu(model)),
    )

    mkpath(dirname(file))
    BSON.bson(file, dict)
    return
end

function loadmodel(d::Data, t::Train, m::Model)
    file = modeldir(string(d), string(t), string(m, ".bson"))
    if isfile(file)
        return BSON.load(file)
    else
        error("model not trained")
        return
    end
end

# ------------------------------------------------------------------------------------------
# sumilations
# ------------------------------------------------------------------------------------------
function trainmodel(
    d::Data,
    t::Train,
    m::Model{T},
    x,
    y;
    force = false,
    runon = cpu,
) where {T}

    file = modeldir(string(d), string(t), string(m, ".bson"))
    isfile(file) && !force && return
    isfile(file) && rm(file)

    # create model
    model, pars = build_network(d, t) |> runon
    objective = build_loss(m)
    loss(x, y) = objective(x, y, model, pars)

    # create callback
    savefunc(c, x, y) = savemodel(c, d, t, m, model, x, y)

    cb = CallBack(;
        title = string(nameof(T), ": "),
        iters = t.iters,
        saveat = t.saveat,
        savefunc = savefunc,
    )

    # training
    Random.seed!(t.seed)
    batches = batch_provider(d, t, x, y) |> runon
    opt = t.optimiser(t.step)

    custom_train!(loss, pars, batches, opt; cb = cb)
    return
end

# run simulations
function run_simulations(
    Datasets,
    Trains,
    Models;
    runon = cpu,
    force = false,
)
    Datasets = isa(Datasets, Data) ? [Datasets] : Datasets
    Trains = isa(Trains, Train) ? [Trains] : Trains
    Models = isa(Models, Model) ? [Models] : Models

    for dataset in Datasets
        @info dataset
        train, valid, test = loaddata(dataset) |> runon

        for t in Trains, m in Models
            trainmodel(dataset, t, m, train...; force, runon)

            dict = loadmodel(dataset, t, m)
            update_results!(dataset, t, m, dict, train..., :train; runon)
            update_results!(dataset, t, m, dict, valid..., :valid; runon)
            update_results!(dataset, t, m, dict, test..., :test; runon)

            file = modeldir(string(dataset), string(t), string(m, ".bson"))
            BSON.bson(file, dict)
        end
    end
    return
end

function update_results!(
    d::Data,
    t::Train,
    m::Model,
    dict::Dict,
    x,
    y,
    key::Symbol;
    runon = cpu
)

    model0, pars0 = build_network(d, t; zero = true) |> runon
    for iter in keys(dict[:iterations])
        model = deepcopy(dict[:iterations][iter][:model]) |> runon
        pars = params(model)
        delete!(pars, model.b)
        loss = build_loss(m)

        get!(dict[:iterations][iter], key, Dict{Symbol, Any}())
        get!(dict[:iterations][iter][key], :targets, cpu(vec(y)))
        get!(dict[:iterations][iter][key], :scores, cpu(compute_scores(model, x)))
        get!(dict[:iterations][iter][key], :loss, cpu(loss(x, y, model, pars)))
        get!(dict[:iterations][iter][key], :loss_zero, cpu(loss(x, y, model0, pars0)))
    end
end

# ------------------------------------------------------------------------------------------
# benchmarks
# ------------------------------------------------------------------------------------------
function benchmarkmodel(
    d::Data,
    t::Train,
    m::Model{T},
    x,
    y;
    runon = cpu,
) where {T}

    # create model
    model, pars = build_network(d, t) |> runon
    objective = build_loss(m)
    loss(x, y) = objective(x, y, model, pars)

    # training
    if d.batchsize == 0
        batches = [(x, y)]
    else
        batches = Flux.Data.DataLoader((x, y); batchsize = 512, shuffle = true, partial = false)
    end
    batches = batches |> runon
    opt = t.optimiser(t.step)

    bench = @benchmark Flux.train!($loss, $pars, $batches, $opt)

    dict = merge(todict(d), todict(m))
    tms = bench.times .* 1e-9
    dict[:time] = measurement(mean(tms), std(tms))

    return select(DataFrame(dict), [:dataset, :model, :time])
end

function run_benchmarks(Datasets; runon = cpu, force = false)

    file = datadir("results", "benchmark.csv")
    isfile(file) && !force && return CSV.read(file, DataFrame; header = true)

    Datasets = isa(Datasets, Data) ? [Datasets] : Datasets
    Models = vcat(
        Model(PatMat; τ = 0.05, β = 1, surrogate = hinge),
        Model(PatMatNP; τ = 0.05, β = 1, surrogate = hinge),
        Model(TopPush; λ = 1),
        Model(TopPushK; K = 5),
        Model(τFPL; τ = 0.05, λ = 1, surrogate = hinge),
        Model(TopMean; τ = 0.05, λ = 1, surrogate = hinge),
        Model(Grill; τ = 0.05, λ = 1, surrogate = hinge),
        Model(GrillNP; τ = 0.05, λ = 1, surrogate = hinge),
    )
    t = Train(; seed = 1, iters = 1, optimiser = ADAM, step = 0.01)

    rows = []
    runs = length(Datasets) * length(Models)
    bar = Progress(runs, 5, "Benchmarks")
    for dataset in Datasets
        train, valid, test = loaddata(dataset) |> runon
        for m in Models
            row = benchmarkmodel(dataset, t, m, train...; runon)
            next!(bar)
            push!(rows, row)
        end
    end
    table = unstack(reduce(vcat, rows), :model, :time)

    mkpath(dirname(file))
    CSV.write(file, table)
    return table
end


text:
module Algorithms
export create_ncu_circuit, create_grover_circuit

include("algs/NCU.jl")
include("algs/Oracle.jl")
include("algs/Diffusion.jl")
include("algs/Grover.jl")

using QuantExQASM.GateOps
using QuantExQASM.Circuit
using PicoQuant

function create_ncu_circuit(num_qubits::Int, use_aux::Bool=false, gate_label::GateOps.GateLabel=GateOps.GateLabel(:x))
    # Create empty circuit with qiven qubit count
    cct = Circuit.Circ(num_qubits)

    # Initialise default intermediate gates for use in NCU
    NCU.init_intermed_gates(cct, num_qubits-1)

    try # check if accessing Gate exists
        U = Circuit.gate_cache[gate_label]
    catch e
        println("GateLabel $(gate_label) does not exist.")
        throw(e)
    end

    if use_aux == false || num_qubits < 9 # If insufficient qubits, use unoptimised NCU
        ctrl = collect(0:num_qubits-2)
        aux = Int[]
        tgt = num_qubits-1
    else
        num_ctrl = Int(ceil(num_qubits/2))
        num_aux = num_qubits - num_ctrl -1
        ctrl = collect(0:num_ctrl-1)
        aux = collect(num_ctrl:num_ctrl+num_aux-1)
        tgt = num_qubits-1
    end

    for i in ctrl
        Circuit.add_gatecall!(cct, GateOps.pauli_x(i) )
    end

    NCU.apply_ncu!(cct, ctrl, aux, tgt, gate_label);

    # Convert circuit to QASM populated buffer
    Circuit.to_qasm(cct, true, "out.qasm")
    cct_s = String( Circuit.to_qasm(cct, true) )
    qq_cct = PicoQuant.load_qasm_as_circuit( cct_s )
    return qq_cct
end

function create_grover_circuit(num_qubits::Int, use_aux::Bool=false, bit_pattern::Int=11)
    # Create empty circuit with qiven qubit count
    cct = Circuit.Circ(num_qubits)

    # Initialise default intermediate gates for use in NCU
    NCU.init_intermed_gates(cct, num_qubits-1)

    if use_aux == false || num_qubits < 9 # If insufficient qubits, use unoptimised NCU
        Grover.run_grover!(cct, collect(0:num_qubits-1), bit_pattern)
    else
        num_ctrl = Int(ceil(num_qubits/2))
        num_aux = num_qubits - num_ctrl -1
        ctrl = collect(0:num_ctrl-1)
        aux = collect(num_ctrl:num_ctrl+num_aux-1)
        tgt = num_qubits-1

        Grover.run_grover!(cct, vcat(ctrl, tgt), aux, bit_pattern)
    end

    #bit_pattern_s = reverse(join(digits(bit_pattern, base=2, pad=num_qubits)))

    cct_s = String(Circuit.to_qasm(cct, true) )
    qq_cct = PicoQuant.load_qasm_as_circuit( cct_s )
    return qq_cct
end

end

ext: jl


text:
#
# Factors Iteration
#
# Iterating over indices, potential, etc.

Base.start(ϕ::Factor) = 1
Base.next(ϕ::Factor, i) = (sub2at(ϕ, ind2sub(ϕ, i)...), i+1)
Base.done(ϕ::Factor, i) = (i == length(ϕ))

"""
    pattern(ϕ::Factor, [dim])

Return the index of the states of `dim`, given its current position in ϕ.
"""
pattern(ϕ::Factor, dim::Symbol) = pattern(ϕ, [dim])

function pattern(ϕ::Factor, dims::Vector{Symbol})
    isempty(dims) && return Int[]

    _check_dims_valid(dims, ϕ)

    inds = indexin(dims, ϕ)
    len = length(ϕ)
    p = Array{Int}(len, length(dims))

    for i in 1:len
        t = ind2sub(ϕ, i)
        p[i, :] = [t[inds]...]
    end

    return p
end

function pattern(ϕ::Factor)
    sz = size(ϕ)

    len = length(ϕ)
    p = Array{Int}(len, ndims(ϕ))

    for i in 1:len
        p[i, :] = [ind2sub(sz, i)...]
    end

    return p
end

"""
    pattern_states(ϕ::Factor, [dim])

Returns the states of `dim`, given its current position in ϕ, for all possible states.
"""
pattern_states(ϕ::Factor, d::Symbol) = getdim(ϕ, d)[pattern(ϕ, d)]

function pattern_states(ϕ::Factor, dims::Vector{Symbol})
    isempty(dims) && return []

    pat = pattern(ϕ, dims)
    dims = getdim(ϕ, dims)

    return hcat([d[pat[:, i]] for (i, d) in enumerate(dims)]...)
end

function pattern_states(ϕ::Factor)
    pat = pattern(ϕ)

    return hcat([d[pat[:, i]] for (i, d) in enumerate(ϕ.dimensions)]...)
end


text:
num = 600851475143
sqr = sqrt(num)
c = 3
max = num
while c <= sqr
    while num % c == 0
        num /= c
        max = c
    end
    c += 2
end
println(max)

max_stars_repo_path: 003 - Largest Prime Factor/003 - Largest_Prime_Factor.jl

max_stars_repo_name: LanguageAdventures/Project-Euler-Julia

max_issues_repo_name: LanguageAdventures/Project-Euler-Julia

max_forks_repo_name: LanguageAdventures/Project-Euler-Julia


text:
##################################################################################
#
#  ORTHOFILTER TRANSFORMS
#  Periodic boundaries, Orthogonal
#
##################################################################################

# DWT
# 1-D
# writes to y
function _dwt!(y::AbstractVector{Ty}, x::AbstractVector{Tx},
                filter::OrthoFilter, L::Integer, fw::Bool) where {Tx<:Number, Ty<:Number}
    T = promote_type(Tx, Ty)
    si = Vector{T}(undef, length(filter)-1) # tmp filter vector
    scfilter, dcfilter = WT.makereverseqmfpair(filter, fw, T)
    return _dwt!(y, x, filter, L, fw, dcfilter, scfilter, si)
end
function _dwt!(y::AbstractVector{<:Number}, x::AbstractVector{<:Number},
                filter::OrthoFilter, L::Integer, fw::Bool,
                dcfilter::Vector{T}, scfilter::Vector{T},
                si::Vector{T}, snew::Vector{T} = Vector{T}(undef, ifelse(L>1, length(x)>>1, 0))) where T<:Number
    n = length(x)
    size(x) == size(y) ||
        throw(DimensionMismatch("in and out array size must match"))
    0 <= L ||
        throw(ArgumentError("L must be positive"))
    sufficientpoweroftwo(y, L) ||
        throw(ArgumentError("size must have a sufficient power of 2 factor"))
    y === x &&
        throw(ArgumentError("in array is out array"))
    L > 1 && length(snew) != n>>1 &&
        throw(ArgumentError("length of snew incorrect"))

    if L == 0
        return copyto!(y,x)
    end
    s = x                           # s is current scaling coefs location
    filtlen = length(filter)

    lrange = 1:L
    !fw && (lrange = reverse(lrange))

    for l in lrange
        if fw
            # detail coefficients
            filtdown!(dcfilter, si, y, detailindex(n,l,1), detailn(n,l), s, 1,-filtlen+1, true)
            # scaling coefficients
            filtdown!(scfilter, si, y, 1, detailn(n,l), s, 1, 0, false)
        else
            # scaling coefficients
            filtup!(false, scfilter, si, y, 1, detailn(n,l-1), s, 1, -filtlen+1, false)
            # detail coefficients
            filtup!(true,  dcfilter, si, y, 1, detailn(n,l-1), x, detailindex(n,l,1), 0, true)
        end
        # if not final iteration: copy to tmp location
        l != lrange[end] && copyto!(snew,1,y,1,detailn(n, fw ? l : l-1))
        L > 1 && (s = snew)
    end
    return y
end
function unsafe_dwt1level!(y::AbstractVector{<:Number}, x::AbstractVector{<:Number},
                            filter::OrthoFilter, fw::Bool,
                            dcfilter::FVector{T}, scfilter::FVector{T},
                            si::FVector{T}) where T<:Number
    n = length(x)
    l = 1
    filtlen = length(filter)

    if fw
        # detail coefficients
        filtdown!(dcfilter, si, y, detailindex(n,l,1), detailn(n,l), x, 1,-filtlen+1, true)
        # scaling coefficients
        filtdown!(scfilter, si, y, 1, detailn(n,l), x, 1, 0, false)
    else
        # scaling coefficients
        filtup!(false, scfilter, si, y, 1, detailn(n,l-1), x, 1, -filtlen+1, false)
        # detail coefficients
        filtup!(true,  dcfilter, si, y, 1, detailn(n,l-1), x, detailindex(n,l,1), 0, true)
    end
    return y
end

function dwt_transform_strided!(y::AbstractArray{<:Number}, x::AbstractArray{<:Number},
                            msub::Int, nsub::Int, stride::Int, idx_func::Function,
                            tmpvec::FVector{T}, tmpvec2::FVector{T},
                            filter::OrthoFilter, fw::Bool,
                            dcfilter::FVector{T}, scfilter::FVector{T}, si::FVector{T}) where T<:Number
    for i=1:msub
        xi = idx_func(i)
        stridedcopy!(tmpvec, x, xi, stride, nsub)
        unsafe_dwt1level!(tmpvec2, tmpvec, filter, fw, dcfilter, scfilter, si)
        stridedcopy!(y, xi, stride, tmpvec2, nsub)
    end
end

function dwt_transform_cols!(y::AbstractArray{<:Number}, x::AbstractArray{<:Number},
                            msub::Int, nsub::Int, idx_func::Function,
                            tmpvec::FVector{T},
                            filter::OrthoFilter, fw::Bool,
                            dcfilter::FVector{T}, scfilter::FVector{T}, si::FVector{T}) where T<:Number
    for i=1:nsub
        xi = idx_func(i)
        copyto!(tmpvec, 1, x, xi, msub)
        ya = unsafe_vectorslice(y, xi, msub)
        unsafe_dwt1level!(ya, tmpvec, filter, fw, dcfilter, scfilter, si)
    end
end

# 2-D
# writes to y
function _dwt!(y::AbstractMatrix{Ty}, x::AbstractMatrix{Tx},
                filter::OrthoFilter, L::Integer, fw::Bool) where {Tx<:Number, Ty<:Number}
    m, n = size(x)
    T = promote_type(Tx, Ty)
    si = Vector{T}(undef, length(filter)-1)       # tmp filter vector
    tmpbuffer = Vector{T}(undef, max(n<<1, m))    # tmp storage vector
    scfilter, dcfilter = WT.makereverseqmfpair(filter, fw, T)

    return _dwt!(y, x, filter, L, fw, dcfilter, scfilter, si, tmpbuffer)
end
function _dwt!(y::AbstractMatrix{<:Number}, x::AbstractMatrix{<:Number},
                filter::OrthoFilter, L::Integer, fw::Bool,
                dcfilter::Vector{T}, scfilter::Vector{T},
                si::Vector{T}, tmpbuffer::Vector{T}) where T<:Number

    m, n = size(x)
    size(x) == size(y) ||
        throw(DimensionMismatch("in and out array size must match"))
    0 <= L ||
        throw(ArgumentError("L must be positive"))
    sufficientpoweroftwo(y, L) ||
        throw(ArgumentError("size must have a sufficient power of 2 factor"))
    y === x &&
        throw(ArgumentError("in array is out array"))
    length(tmpbuffer) >= max(n<<1,m) ||
        throw(ArgumentError("length of tmpbuffer incorrect"))

    if L == 0
        return copyto!(y,x)
    end
    row_stride = m
    #s = x

    if fw
        lrange = 1:L
        nsub = n
        msub = m
    else
        lrange = L:-1:1
        nsub = div(n,2^(L-1))
        msub = div(m,2^(L-1))
        copyto!(y,x)
    end

    inputArray = x

    row_idx_func = i -> row_idx(i, m)
    col_idx_func = i -> col_idx(i, m)
    for l in lrange
        tmpvec  = unsafe_vectorslice(tmpbuffer, 1, msub)
        tmpvec2 = unsafe_vectorslice(tmpbuffer, 1, nsub)
        tmpvec3 = unsafe_vectorslice(tmpbuffer, nsub+1, nsub)
        if fw
            # rows
            dwt_transform_strided!(y, inputArray, msub, nsub, row_stride, row_idx_func,
                tmpvec2, tmpvec3, filter, fw, dcfilter, scfilter, si)
            l == lrange[1] && (inputArray = y)

            # columns
            dwt_transform_cols!(y, y, msub, nsub, col_idx_func,
                tmpvec, filter, fw, dcfilter, scfilter, si)
        else
            # columns
            dwt_transform_cols!(y, inputArray, msub, nsub, col_idx_func,
                tmpvec, filter, fw, dcfilter, scfilter, si)
            l == lrange[1] && (inputArray = y)

            # rows
            dwt_transform_strided!(y, y, msub, nsub, row_stride, row_idx_func,
                tmpvec2, tmpvec3, filter, fw, dcfilter, scfilter, si)
        end
        msub = (fw ? msub>>1 : msub<<1)
        nsub = (fw ? nsub>>1 : nsub<<1)
    end
    return y
end

# 3-D
# writes to y
function _dwt!(y::AbstractArray{Ty, 3}, x::AbstractArray{Tx, 3},
               filter::OrthoFilter, L::Integer, fw::Bool) where {Tx<:Number, Ty<:Number}
    m, n, d = size(x)
    T = promote_type(Tx, Ty)
    si = Vector{T}(undef, length(filter)-1)            # tmp filter vector
    tmpbuffer = Vector{T}(undef, max(m, n<<1, d<<1))   # tmp storage vector
    scfilter, dcfilter = WT.makereverseqmfpair(filter, fw, T)

    return _dwt!(y, x, filter, L, fw, dcfilter, scfilter, si, tmpbuffer)
end
function _dwt!(y::AbstractArray{<:Number, 3}, x::AbstractArray{<:Number, 3},
                filter::OrthoFilter, L::Integer, fw::Bool,
                dcfilter::Vector{T}, scfilter::Vector{T},
                si::Vector{T}, tmpbuffer::Vector{T}) where T<:Number

    m, n, d = size(x)
    size(x) == size(y) ||
        throw(DimensionMismatch("in and out array size must match"))
    0 <= L ||
        throw(ArgumentError("L must be positive"))
    sufficientpoweroftwo(y, L) ||
        throw(ArgumentError("size must have a sufficient power of 2 factor"))
    y === x &&
        throw(ArgumentError("in array is out array"))
    length(tmpbuffer) >= n<<1 ||
        throw(ArgumentError("length of tmpbuffer incorrect"))

    if L == 0
        return copyto!(y,x)
    end
    row_stride = m
    plane_stride = m*n

    if fw
        lrange = 1:L
        msub = m
        nsub = n
        dsub = d
    else
        lrange = L:-1:1
        msub = div(m,2^(L-1))
        nsub = div(n,2^(L-1))
        dsub = div(d,2^(L-1))
        copyto!(y,x)
    end

    inputArray = x

    for l in lrange
        tmpcol  = unsafe_vectorslice(tmpbuffer, 1, msub)
        tmprow  = unsafe_vectorslice(tmpbuffer, 1, nsub)
        tmprow2 = unsafe_vectorslice(tmpbuffer, nsub+1, nsub)
        tmphei  = unsafe_vectorslice(tmpbuffer, 1, dsub)
        tmphei2 = unsafe_vectorslice(tmpbuffer, dsub+1, dsub)
        if fw
            # planes
            for j in 1:nsub
                plane_idx_func = i -> plane_idx(i, j, m)
                dwt_transform_strided!(y, inputArray, msub, dsub, plane_stride, plane_idx_func,
                    tmphei, tmphei2, filter, fw, dcfilter, scfilter, si)
            end
            l == lrange[1] && (inputArray = y)

            # rows
            for j in 1:dsub
                row_idx_func = i -> row_idx(i, j, m, n)
                dwt_transform_strided!(y, y, msub, nsub, row_stride, row_idx_func,
                    tmprow, tmprow2, filter, fw, dcfilter, scfilter, si)
            end
            # columns
            for j in 1:dsub
                col_idx_func = i -> col_idx(i, j, m, n)
                dwt_transform_cols!(y, y, msub, nsub, col_idx_func,
                    tmpcol, filter, fw, dcfilter, scfilter, si)
            end
        else
            # columns
            for j in 1:dsub
                col_idx_func = i -> col_idx(i, j, m, n)
                dwt_transform_cols!(y, inputArray, msub, nsub, col_idx_func,
                    tmpcol, filter, fw, dcfilter, scfilter, si)
            end
            l == lrange[1] && (inputArray = y)

            # rows
            for j in 1:dsub
                row_idx_func = i -> row_idx(i, j, m, n)
                dwt_transform_strided!(y, y, msub, nsub, row_stride, row_idx_func,
                    tmprow, tmprow2, filter, fw, dcfilter, scfilter, si)
            end
            # planes
            for j in 1:nsub
                plane_idx_func = i -> plane_idx(i, j, m)
                dwt_transform_strided!(y, y, msub, dsub, plane_stride, plane_idx_func,
                    tmphei, tmphei2, filter, fw, dcfilter, scfilter, si)
            end
        end
        msub = (fw ? msub>>1 : msub<<1)
        nsub = (fw ? nsub>>1 : nsub<<1)
        dsub = (fw ? dsub>>1 : dsub<<1)
    end
    return y
end

# WPT
# 1-D
# writes to y
function _wpt!(y::AbstractVector{T}, x::AbstractVector{T}, filter::OrthoFilter, tree::BitVector, fw::Bool) where T<:Number
    si = Vector{T}(undef, length(filter)-1)
    ns = ifelse(fw, length(x)>>1, length(x))
    snew = Vector{T}(undef, ns)
    scfilter, dcfilter = WT.makereverseqmfpair(filter, fw, T)

    return _wpt!(y, x, filter, tree, fw, dcfilter, scfilter, si, snew)
end
function _wpt!(y::AbstractVector{T}, x::AbstractVector{T}, filter::OrthoFilter, tree::BitVector, fw::Bool, dcfilter::Vector{T}, scfilter::Vector{T}, si::Vector{T}, snew::Vector{T}) where T<:Number

    size(x) == size(y) ||
        throw(DimensionMismatch("in and out array size must match"))
    y === x &&
        throw(ArgumentError("in array is out array"))
    isvalidtree(y, tree) ||
        throw(ArgumentError("invalid tree"))
    if tree[1] && length(snew) < ifelse(fw, length(x)>>1, length(x))
        throw(ArgumentError("length of snew incorrect"))
    end

    if !tree[1]
        return copyto!(y,x)
    end

    first = true
    n = length(x)
    Lmax = maxtransformlevels(n)
    L = Lmax
    while L > 0
        ix = 1
        k = 1
        Lfw = (fw ? Lmax-L : L-1)
        nj = detailn(n, Lfw)
        treeind = 2^(Lfw)-1
        dx = first ? x : unsafe_vectorslice(snew, 1, nj) # dx will be overwritten if first

        while ix <= n
            if tree[treeind+k]
                dy = unsafe_vectorslice(y, ix, nj)
                if first
                    dx = unsafe_vectorslice(x, ix, nj)
                else
                    copyto!(dx, dy)
                end
                unsafe_dwt1level!(dy, dx, filter, fw, dcfilter, scfilter, si)
            elseif first
                dy = unsafe_vectorslice(y, ix, nj)
                dx = unsafe_vectorslice(x, ix, nj)
                copyto!(dy, dx)
            end
            ix += nj
            k += 1
        end
        L -= 1
        first = false
    end

    return y
end

macro filtermainloop(si, silen, b, val)
    quote
        @inbounds for j=2:$(esc(silen))
            $(esc(si))[j-1] = $(esc(si))[j] + $(esc(b))[j]*$(esc(val))
        end
        @inbounds $(esc(si))[$(esc(silen))] = $(esc(b))[$(esc(silen))+1]*$(esc(val))
    end
end
macro filtermainloopzero(si, silen)
    quote
        @inbounds for j=2:$(esc(silen))
            $(esc(si))[j-1] = $(esc(si))[j]
        end
        @inbounds $(esc(si))[$(esc(silen))] = 0.0
    end
end

# periodic filter and downsampling (by 2)
# apply f to x[ix:ix+nx-1] and write to out[iout:iout+nout-1]
# f : filter
# si: tmp array of length 1 less than f
# out : result gets written to out[iout:iout+nout-1]
# x : filter convolved with x[ix:ix+nx-1], where nx=nout*2 (shifted by shift)
# ss : shift downsampling
# based on Base.filt
function filtdown!(f::AbstractVector{T}, si::AbstractVector{T},
                  out::AbstractVector{<:Number}, iout::Integer, nout::Integer,
                  x::AbstractVector{<:Number}, ix::Integer,
                  shift::Integer=0, ss::Bool=false) where T<:Number
    nx = nout<<1
    silen = length(si)
    flen = length(f)
    @assert length(x) >= ix+nx-1
    @assert length(out) >= iout+nout-1
    @assert (nx-1)>>1 + iout <= length(out)
    @assert silen == flen-1
    @assert shift <= 0

    fill!(si,0.0)
    istart = flen + Int(ss)
    dsshift = (flen%2 + Int(ss))%2  # is flen odd, and shift downsampling

    rout1, rin, rout2 = splitdownrangeper(istart, ix, nx, shift)
    @inbounds begin
        for i in rout1  # rout1 assumed to be in 1:istart-1
            # periodic in the range [ix:ix+nx-1]
            xatind = x[mod(i-1+shift, nx) + ix]
            @filtermainloop(si, silen, f, xatind)
        end
        # rin is inbounds in [ix:ix+nx-1]
        ixsh = -1 + shift + ix
        for i in rin
            xatind = x[i + ixsh]
            # take every other value after istart
            if (i+dsshift)%2 == 0 && i >= istart
                out[(i-istart)>>1 + iout] = si[1] + f[1]*xatind
            end
            @filtermainloop(si, silen, f, xatind)
        end
        for i in rout2
            # periodic in the range [ix:ix+nx-1]
            xatind = x[mod(i-1+shift, nx) + ix]
            # take every other value after istart
            if (i+dsshift)%2 == 0 && i >= istart
                out[(i-istart)>>1 + iout] = si[1] + f[1]*xatind
            end
            @filtermainloop(si, silen, f, xatind)
        end
    end

    return nothing
end
# find part of range which is inbounds for [ix:ix+nx-1] (ixsh = -1 + shift + ix)
# where mod(i-1+shift, nx) + ix == i + ixsh
function splitdownrangeper(istart, ix, nx, shift)
    inxi = 0
    ixsh = -1 + shift + ix
    if mod(shift, nx) + ix == 1 + ixsh   # shift likely 0
        inxi = 1
        iend = nx-1
        while mod(iend - 1 + shift, nx) + ix != iend + ixsh
            iend -= 1
        end
        return (0:-1, inxi:iend, (iend+1):(nx-1+istart))
    elseif mod(istart - 1 + shift, nx) + ix == istart + ixsh
        inxi = istart
        iend = nx-1+istart
        while mod(iend - 1 + shift, nx) + ix != iend + ixsh
            iend -= 1
        end
        return (1:inxi-1, inxi:iend, (iend+1):(nx-1+istart))
    else
        return (0:-1, 0:-1, 1:(nx-1+istart))
    end
end

# periodic filter and upsampling (by 2)
# apply f to x[ix:ix+nx-1] upsampled and write to out[iout:iout+nout-1]
# f : filter
# si: tmp array of length 1 less than f
# out : result gets written to out[iout:iout+nout-1]
# x : filter convolved with x[ix:ix+nx-1] upsampled, where nout==nx*2 (then shifted by shift)
# ss : shift upsampling
# based on Base.filt
function filtup!(add2out::Bool, f::Vector{T}, si::Vector{T},
              out::AbstractVector{<:Number}, iout::Integer, nout::Integer,
              x::AbstractVector{<:Number}, ix::Integer,
              shift::Integer=0, ss::Bool=false) where T<:Number
    nx = nout>>1
    silen = length(si)
    flen = length(f)
    @assert length(x) >= ix+nx-1                # check array size
    @assert length(out) >= iout+nout-1          # check array size
    @assert (nx<<1-1)>>1 + iout <= length(out)  # max for out index
    @assert silen == flen-1
    @assert shift <= 0

    fill!(si,0.0)
    istart = flen - shift%2
    dsshift = Int(ss)%2  # shift upsampling

    rout1, rin, rout2 = splituprangeper(istart, ix, nx, nout, shift)
    @inbounds begin
        xatind = 0.0
        for i in rout1  # rout1 assumed to be in 1:istart-1
            if (i+dsshift)%2 == 0
                @filtermainloopzero(si, silen)
            else
                # periodic in the range [ix:ix+nx-1]
                xindex = mod((i-1)>>1+shift>>1,nx) + ix    #(i-1)>>1 increm. every other
                xatind = x[xindex]
                @filtermainloop(si, silen, f, xatind)
            end
        end
        ixsh = shift>>1 + ix
        for i in rin
            if (i+dsshift)%2 == 0
                xatind = 0.0
            else
                xindex = (i-1)>>1 + ixsh
                xatind = x[xindex]
            end
            if i >= istart
                if add2out
                    out[(i-istart) + iout] += si[1] + f[1]*xatind
                else
                    out[(i-istart) + iout] = si[1] + f[1]*xatind
                end
            end
            if (i+dsshift)%2 == 0
                @filtermainloopzero(si, silen)
            else
                @filtermainloop(si, silen, f, xatind)
            end
        end
        for i in rout2
            if (i+dsshift)%2 == 0
                xatind = 0.0
            else
                xindex = mod((i-1)>>1+shift>>1,nx) + ix
                xatind = x[xindex]
            end
            if i >= istart
                if add2out
                    out[(i-istart) + iout] += si[1] + f[1]*xatind
                else
                    out[(i-istart) + iout] = si[1] + f[1]*xatind
                end
            end
            if (i+dsshift)%2 == 0
                @filtermainloopzero(si, silen)
            else
                @filtermainloop(si, silen, f, xatind)
            end
        end
    end

    return nothing
end
# find part of range which is inbounds for [ix:ix+nx-1] (ixsh = shift>>1 + ix)
# where mod((i-1)>>1+shift>>1, nx) + ix == (i-1)>>1 + ixsh
function splituprangeper(istart, ix, nx, nout, shift)
    inxi = 0
    ixsh = shift>>1 + ix
    if mod(shift>>1, nx) + ix == ixsh   # shift likely 0
        inxi = 1
        iend = nout-1
        while mod((iend-1)>>1+shift>>1, nx) + ix != (iend-1)>>1 + ixsh
            iend -= 1
        end
        return (0:-1, inxi:iend, (iend+1):(nout-1+istart))
    elseif mod((istart-1)>>1+shift>>1, nx) + ix == (istart-1)>>1 + ixsh
        inxi = istart
        iend = nout-1+istart
        while mod((iend-1)>>1+shift>>1, nx) + ix != (iend-1)>>1 + ixsh
            iend -= 1
        end
        return (1:inxi-1, inxi:iend, (iend+1):(nout-1+istart))
    else
        return (0:-1, 0:-1, 1:(nout-1+istart))
    end
end


text:
"Parse Julog terms using Prolog-like syntax."
function parse_term(expr)
    if isa(expr, Union{String,Number,Enum,Bool})
        # Strings, numbers, enums and bools are automatically constants
        return :(Const($expr))
    elseif isa(expr, Symbol)
        if expr == :_
            # Underscores are wildcards, parsed to variables with new names
            return Expr(:call, Var, Meta.quot(gensym(expr)))
        elseif isuppercase(string(expr)[1]) || string(expr)[1] == "_"
            # Initial uppercase / underscored symbols are parsed to variables
            return Expr(:call, Var, Meta.quot(expr))
        else
            # Everything else is parsed to a constant
            return Expr(:call, Const, Meta.quot(expr))
        end
    elseif isa(expr, Expr) && expr.head == :tuple
        # Fully evaluated tuples are valid constants
        return :(Const($(esc(expr))))
    elseif isa(expr, Expr) && expr.head == :vect
        # Parse vector as Prolog style list
        return parse_list(expr.args)
    elseif isa(expr, Expr) && expr.head == :$
        # Evaluate interpolated expression as Const within scope of caller
        val = expr.args[1]
        return :(Const($(esc(val))))
    elseif isa(expr, Expr) && expr.head == :call
        # A compound term comprises its name (functor) and arguments
        name = Meta.quot(expr.args[1])
        args = [parse_term(e) for e in expr.args[2:end]]
        return :(Compound($name, [$(args...)]))
    elseif isa(expr, QuoteNode)
        # Evaluate quoted expression within scope of caller
        val = expr.value
        return :($(esc(val)))
    else
        dump(expr)
        error("syntax error in Julog term at $expr")
    end
end

"Parse arguments of vector using Prolog-style list syntax."
function parse_list(args)
    if (length(args) > 0 && isa(args[end], Expr) &&
        args[end].head == :call && args[end].args[1] == :|)
        # Handle [... | Tail] syntax for last element
        tail = parse_term(args[end].args[3])
        pretail = parse_term(args[end].args[2])
        tail = :(Compound(:cons, [$pretail, $tail]))
        args = args[1:end-1]
    else
        # Initialize tail to empty list
        tail = :(Compound(:cend, []))
    end
    # Recursively build list using :cons
    elts = [parse_term(a) for a in args]
    for e in reverse(elts)
        tail = :(Compound(:cons, [$e, $tail]))
    end
    return tail
end

"Parse body of Julog clause using Prolog-like syntax, with '&' replacing ','."
function parse_body(expr)
    if expr == true
        return []
    elseif isa(expr, Symbol) || isa(expr, QuoteNode) ||
           isa(expr, Expr) && expr.head in [:tuple, :vect, :$]
        return [parse_term(expr)]
    elseif isa(expr, Expr) && expr.head == :call
        if expr.args[1] == :&
            # '&' is left-associative, so we descend the parse tree accordingly
            return [parse_body(expr.args[2]); parse_term(expr.args[3])]
        else
            return [parse_term(expr)]
        end
    else
        dump(expr)
        error("syntax error in body of Julog clause at $expr")
    end
end

"Parse Julog expression using Prolog-like syntax. '<<=' replaces ':-' in clauses."
function parse_julog(expr)
    if !isa(expr, Expr)
        return parse_term(expr)
    elseif expr.head == :<<=
        head = parse_term(expr.args[1])
        body = parse_body(expr.args[2])
        return :(Clause($head, [$(body...)]))
    elseif expr.head == Symbol("'")
        head = parse_term(expr.args[1])
        return :(Clause($head, []))
    elseif expr.head == :vect
        exprs = [parse_julog(a) for a in expr.args]
        return :([$(exprs...)])
    elseif expr.head == :ref && expr.args[1] == :Clause
        ty = expr.args[1]
        exprs = [parse_julog(a) for a in expr.args[2:end]]
        return :(Clause[$(exprs...)])
    elseif expr.head == :ref && expr.args[1] in [:Term, :Const, :Var, :Compound]
        ty = expr.args[1]
        exprs = [parse_term(a) for a in expr.args[2:end]]
        return :($ty[$(exprs...)])
    elseif expr.head == :ref && expr.args[1] == :list
        return parse_list(expr.args[2:end])
    else
        return parse_term(expr)
    end
end

"""
    @julog expr

Parse and return Julog expressions.

- `@julog <term>` parses a Prolog-style term.
- `@julog <term> <<= true` or `@julog <term>'` parses a fact (body-less clause).
- `@julog <head> <<= <body>` parses a definite clause.
- `@julog [<term|clause>, ...]` parses a vector of terms or clauses
- `@julog <T>[<term>, ...]` parses to a vector of type `T` (e.g. `Const`)
- `@julog list[<term>, ...]` parses a Prolog-style list directly to a term.

Additionally, the `\$` operator can be used to interpolate regular Julia
expressions as Julog constants, while the `:` operator can be used to
interpolate variables referring to pre-constructed Julog terms into
another Julog term or clause.
"""
macro julog(expr)
    return parse_julog(expr)
end

"Parse Julog substitutions, e.g. {X => hello, Y => world}."
macro varsub(expr)
    if !(isa(expr, Expr) && expr.head == :braces)
        error("Invalid format for Julog substitutions.")
    end
    vars = [Var(a.args[2]) for a in expr.args]
    terms = [eval(parse_term(a.args[3])) for a in expr.args]
    return Subst(v => t for (v,t) in zip(vars, terms))
end

"Convert Prolog string to list of Julog strings."
function convert_prolog_to_julog(str::String)
    clauses = String[]
    # Match each clause (being careful to handle periods in lists + floats)
    for m in eachmatch(r"((?:\d\.\d+|[^\.]|\.\()*)\.\s*", str)
        clause = m.captures[1]
        # Replace cuts, negations and implications
        clause = replace(clause, "!" => "cut")
        clause = replace(clause, "\\+" => "!")
        clause = replace(clause, "->" => "=>")
        clause = replace(clause, ".(" => "cons(")
        # Try to match to definite clause
        m = match(r"(.*):-(.*)", clause)
        if isnothing(m)
            # Push fact on to list of clauses
            clause = clause * " <<= true"
            push!(clauses, clause)
            continue
        end
        # Handle definite clauses
        head, body = m.captures[1:2]
        # Handle disjunctions within the body
        for bd in split(body, ";")
            # Handle conjuctions within body (find commas outside of brackets)
            br_count = 0
            commas = [0]
            for (idx, chr) in enumerate(bd)
                if (chr == '(') br_count += 1
                elseif (chr == ')') br_count -= 1
                elseif (chr == ',') && br_count == 0 push!(commas, idx) end
            end
            push!(commas, length(bd)+1)
            terms = [strip(bd[a+1:b-1]) for (a, b) in
                     zip(commas[1:end-1], commas[2:end])]
            subclause = strip(head) * " <<= " * join(terms, " & ")
            push!(clauses, subclause)
        end
    end
    return clauses
end

"Parse Julog expression from string using standard Prolog syntax."
function parse_prolog(str::String)
    strs = convert_prolog_to_julog(str)
    exprs = [@julog($(Meta.parse(s))) for s in strs]
    return exprs
end

"Write list of Julog clauses to Prolog string."
function write_prolog(clauses::Vector{Clause})
    str = ""
    for clause in clauses
        clause = repr(clause)
        clause = replace(clause, "!" => "\\+")
        clause = replace(clause, "=>" => "->")
        clause = replace(clause, " &" => ",")
        clause = replace(clause, "<<=" => ":-")
        str = str * clause * ".\n"
    end
    return str
end

"Parse Prolog program as a string and return a list of Julog clauses."
macro prolog(str::String)
    strs = convert_prolog_to_julog(str)
    exprs = [parse_julog(Meta.parse(s)) for s in strs]
    return :([$(exprs...)])
end


text:
##
# Provides access to map results
type MapResultReader <: MapInputReader
    url::AbstractString
    jid::Int64
    results::Any

    function MapResultReader(url::AbstractString="")
        r = new("", 0, nothing)
        isempty(url) || reset_pos(r, url)
        r
    end
end

function reset_pos(r::MapResultReader, url::AbstractString)
    (url == r.url) && return

    # new source
    r.jid = parse(Int, url)
    j = ((myid() == 1) ? _def_wrkr_job_store : _job_store)[r.jid]
    r.results = j.info.results
    r.url = url
    nothing
end

##
# Input for map
type MRMapInput <: MRInput
    job_list::Vector{JobId}
    reader_fn::Function

    function MRMapInput(job_list, reader_fn::Function)
        jl = JobId[]
        for job_id in job_list
            push!(jl, convert(JobId, job_id))
        end
        new(jl, reader_fn)
    end
end

input_reader_type(inp::MRMapInput) = (MRMapInput, nothing)
get_input_reader(::Type{MRMapInput}, ::Any) = MapResultReader()

max_stars_repo_name: tanmaykm/HadoopBlocks

max_issues_repo_name: tanmaykm/HadoopBlocks

max_forks_repo_name: tanmaykm/HadoopBlocks


text:
struct DiscretizeQuantile <: DiscretizationAlgorithm
    nbins::Int
    trim::Bool # Whether to automatically remove extremely narrow bins caused by duplicate entries
end

# Set trim to false by default
DiscretizeQuantile(nbins::Int) = DiscretizeQuantile(nbins, false)

function binedges(alg::DiscretizeQuantile, data::AbstractArray{N}) where {N<:Union{Int, AbstractFloat }}
    nbins = alg.nbins

    n = length(data)
    n ≥ nbins || error("too many bins requested")

    qs = range(0, 1; length=(nbins+1))
    retval = quantile(data, qs)

    # trim bins that are too small, caused by duplicate entiries
    keep_bin = diff(retval) .> 1e-8

    if alg.trim # trim bins that are too small
        # Ref: https://github.com/scikit-learn/scikit-learn/blob/main/sklearn/preprocessing/_discretization.py#L288
        # Always keep the first bin edge
        insert!(keep_bin, 1, true)
        retval = retval[keep_bin] # NOTE: this reduces our bin count
        length(retval) > 1 || error("no bins remaining")
    else
        # As with other methods, throw an error if there are non-unique bin edges
        !all(keep_bin) && error("binedges non-unique")
    end

    retval
end


text:
using PDFToc
using Test

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


text:
# ---
# title: 404. Sum of Left Leaves
# id: problem404
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Tree
# link: <https://leetcode.com/problems/sum-of-left-leaves/description/>
# hidden: true
# ---
#
# Find the sum of all left leaves in a given binary tree.
#
# **Example:**
#
#
#
#         3
#        / \
#       9  20
#         /  \
#        15   7
#
#     There are two left leaves in the binary tree, with values **9** and **15** respectively. Return **24**.
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#=
Created on 09/01/2021 15:00:09
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Implement the information entropy:
$$
H(p) = -\sum_i p_i \log p_i\,,
$$

Take caution that, by definition,$0 \log 0$ is evaluated as 0 (i.e., its limit value).

Have an optional flag to specify the bag of log, default is base 2, making the output in bits.

Make a plot of the entropy of a Bernouilli-distributed variable for varying probability of success
$p$, i.e., the entropy of probability vector $[p, 1-p]$ for $p\in [0, 1]$.
=#

using Plots

entropy(p; base=2) = p .|> (pᵢ -> pᵢ > 0.0 ? -pᵢ * log(pᵢ) : 0.0) |> sum |> (H -> H / log(base))

entropy([0, 1, 0])  # 0.0
entropy([0.5, 0.5])  # 1.0
entropy([0.5, 0.5], base=10)  # 0.30102999566398114

plot(p->entropy((p, 1-p)), 0, 1, xlabel="p", label="entropy")


text:
mutable struct NumFieldEmbNfAbsNS <: NumFieldEmb{NfAbsNS}
  field::NfAbsNS
  index::Vector{Int}
  absolute_index::Int
  isreal::Bool
  roots::Vector{acb}
  conjugate::Int
end

_absolute_index(P::NumFieldEmbNfAbsNS) = P.absolute_index

number_field(P::NumFieldEmbNfAbsNS) = P.field

embedding_type(::Type{NfAbsNS}) = NumFieldEmbNfAbsNS

embedding_type(::NfAbsNS) = NumFieldEmbNfAbsNS

isreal(P::NumFieldEmbNfAbsNS) = P.isreal

isimaginary(P::NumFieldEmbNfAbsNS) = !P.isreal

conj(P::NumFieldEmbNfAbsNS) = complex_embeddings(number_field(P))[P.conjugate]

function complex_embeddings(K::NfAbsNS; conjugates::Bool = true)
  emb = get_attribute!(K, :complex_embeddings) do
    _complex_embeddings(K)
  end::Vector{NumFieldEmbNfAbsNS}
  if conjugates
    return emb
  else
    r, s = signature(K)
    return emb[1:r + s]
  end
end

function _complex_embeddings(K::NfAbsNS)
  c = conjugate_data_arb_roots(K, 32, copy = false)

  r, s = signature(K)
  res = Vector{NumFieldEmbNfAbsNS}(undef, degree(K))

  l = ngens(K)

  j = 1

  for v in c[2]
    res[j] = NumFieldEmbNfAbsNS(K, v, j, true, acb[c[1][i].roots[v[i]] for i in 1:l], j)
    j += 1
  end

  for v in c[3]
    res[j] = NumFieldEmbNfAbsNS(K, v, j, false, acb[c[1][i].roots[v[i]] for i in 1:l], j + s)
    res[j + s] = NumFieldEmbNfAbsNS(K, v, j + s, false, acb[conj(c[1][i].roots[v[i]]) for i in 1:l], j)
    j += 1
  end

  return res
end

function Base.show(io::IO, ::MIME"text/plain", f::NumFieldEmbNfAbsNS)
  print(io, "Embedding of\n")
  print(IOContext(io, :compact => true), number_field(f))
  print(io, "\n")
  print(io, "with roots ≈ ")
  print(io, "[ ")
  for i in 1:length(f.roots)
    _print_acb_neatly(io, f.roots[i])
    if i < length(f.roots)
      print(io, ", ")
    end
  end
  print(io, "]")
end

function Base.show(io::IO, f::NumFieldEmbNfAbsNS)
  print(io, "Embedding corresponding to ≈ ")
  print(io, "[ ")
  for i in 1:length(f.roots)
    _print_acb_neatly(io, f.roots[i])
    if i < length(f.roots)
      print(io, ", ")
    end
  end
  print(io, "]")
end

function (f::NumFieldEmbNfAbsNS)(a::NfAbsNSElem, prec::Int = 32)
  K = parent(a)
  wprec = prec
  pol_a = data(a)
  r, s = signature(K)
  i = _absolute_index(f)
  if i > r + s
    return conj(conj(f)(a, prec))
  end
  while true
    conjs, ind_real, ind_complex = conjugate_data_arb_roots(K, wprec)
    if i <= r
      res = _evaluate(pol_a, acb[conjs[j].roots[ind_real[i][j]] for j = 1:ngens(K)])
    else
      i = i - r
      ev = acb[conjs[j].roots[ind_complex[i][j]] for j = 1:ngens(K)]
      res = _evaluate(pol_a, ev)
    end
    if !isfinite(res) || !radiuslttwopower(res, -prec)
      wprec = 2 * wprec
    else
      return res
    end
  end
end

function evaluate(a::NfAbsNSElem, P::NumFieldEmbNfAbsNS, prec::Int)
  return P(a, prec)
end

function restrict(e::NumFieldEmb, f::NumFieldMor{<: NfAbsNS, <: Any, <: Any})
  @req number_field(e) === codomain(f) "Number fields do not match"
  L = domain(f)
  emb = complex_embeddings(L)
  cn = Bool[ all(overlaps.(ee.roots, e.(f.(gens(L))))) for ee in emb]
  @assert count(cn) == 1
  i = findfirst(cn)
  return emb[i]
end

################################################################################
#
#  Easier creation of embeddings
#
################################################################################

function _find_nearest_embedding(K::NfAbsNS, x::Vector{<:Union{BigFloat, Float64}})
  r = real_embeddings(K)
  l = length(K.pol)
  @assert length(x) == l
  gK = gens(K)
  diffs = [[e(gK[i]) - x[i] for i in 1:l] for e in r]
  t = [maximum([abs(w) for w in z]) for z in diffs]
  _, i = findmin(t)
  return r[i]
end

function real_embedding(K::NfAbsNS, x::Vector)
  return _find_nearest_embedding(K, x)
end

function _find_nearest_embedding(K::NfAbsNS, x::Vector{<:Tuple})
  r = real_embeddings(K)
  p = 32
  gK = gens(K)
  fls = [all(_is_contained_in_interval(real(i(gK[j], p)), x[j]) for j in 1:length(gK)) for i in r]
  while count(fls) != 1
    p = 2 * p
    fls = [all(_is_contained_in_interval(real(i(gK[j], p)), x[j]) for j in 1:length(gK)) for i in r]
    if count(fls) > 1
      possible = [[ Float64(real(x)) for x in e.roots] for e in r]
      s = IOBuffer()
      for k in 1:length(possible)
        for i in 1:length(possible[k])
          @printf s "%.2f" possible[k][i]
          if i < length(possible[k])
            print(s, ", ")
          end
        end
        if k < length(possible)
          print(s, "\n")
        end
      end
      ss = String(take!(s))
      error("""Interval contains more than one root. Possible roots are:
            $ss
            """)
    end
    p > 2^8 && error("Something wrong!")
  end
  i = findfirst(fls)
  @assert i !== nothing
  return r[i]
end

function real_embedding(K::NfAbsNS, x::Vector{<:Tuple})
  return _find_nearest_embedding(K, x)
end


text:
using DrWatson
@quickactivate
using ArgParse
using GenerativeAD
using GenerativeAD.Models
import StatsBase: fit!, predict
using StatsBase
using BSON
# because of vae and 2stage
using DataFrames
using CSV
using ValueHistories
using Flux
using ConditionalDists
using GenerativeModels
import GenerativeModels: VAE
using Distributions
using DistributionsAD
using Distances
using Statistics

s = ArgParseSettings()
@add_arg_table! s begin
   "max_seed"
		required = true
		arg_type = Int
		help = "seed"
	"dataset"
		required = true
		arg_type = String
		help = "dataset"
	"contamination"
    	arg_type = Float64
    	help = "contamination rate of training data"
    	default = 0.0
end
parsed_args = parse_args(ARGS, s)
@unpack dataset, max_seed, contamination = parsed_args

# we select the encoding of the model with lowest loss
tab_name = "vae_LOSS_tabular"
mi = 2

#######################################################################################
################ THIS PART IS TO BE PROVIDED FOR EACH MODEL SEPARATELY ################
# get criterion
sp = split(tab_name, "_")
enc = sp[1]
criterion = lowercase(sp[2])

modelname = "$(enc)_ocsvm"

# default hyperparams
function median_l2_dist(X)
	dists = pairwise(Euclidean(), X)
	# take only the upper diagonal
	ds = []
	for i in 1:size(dists,1)
		for j in (i+1):size(dists,1)
			push!(ds, dists[i,j])
		end
	end
	ml=median(ds)
end
function set_params(data)
	D, N = size(data[1][1])
	gamma = 1/median_l2_dist(data[1][1])
	return (gamma=gamma, kernel="rbf", nu=0.5)
end

function fit(data, parameters, aux_info)
	# construct model - constructor should only accept kwargs
	model = GenerativeAD.Models.OCSVM(;parameters...)

	# sumbsample and fit train data
	tr_data = data[1][1]
	M,N = size(tr_data)
	tr_data = tr_data[:, sample(1:N, min(10000, N), replace=false)]
	try
		global info, fit_t, _, _, _ = @timed fit!(model, data[1][1])
	catch e
		# return an empty array if fit fails so nothing is computed
		return (fit_t = NaN,), []
	end

	# construct return information - put e.g. the model structure here for generative models
	training_info = (
		fit_t = fit_t,
		model = nothing
		)

	# now return the different scoring functions
	training_info, [(x->predict(model, x), merge(parameters, aux_info))]
end

####################################################################
################ THIS PART IS COMMON FOR ALL MODELS ################
# set a maximum for parameter sampling retries
cont_string = (contamination == 0.0) ? "" : "_contamination-$contamination"
for seed in 1:max_seed
	savepath = datadir("experiments/tabular_clean_val_default/$(modelname)/$(dataset)/seed=$(seed)")
	mkpath(savepath)

	aux_info = (model_index=mi, criterion=criterion)
	# get data
	data = GenerativeAD.load_data(dataset, seed=seed, contamination=contamination)
	output, encode_t, _, _, _ = @timed GenerativeAD.Models.load_encoding(tab_name, data, dataset=dataset, seed=seed, model_index=mi)
	data, encoding_name, encoder_params, encoder_fit_t = output

	# set parameters
	parameters = set_params(data)

	@info "Trying to fit $modelname on $dataset with parameters $(parameters)..."
	# check if a combination of parameters and seed alread exists
	# fit
	training_info, results = fit(data, parameters, aux_info)
	# here define what additional info should be saved together with parameters, scores, labels and predict times
	save_entries = merge(training_info, (modelname = modelname,
										 seed = seed,
										 dataset = dataset,
										 encoder=encoding_name,
										 encoder_fit_t = encoder_fit_t,
										 encoder_params=encoder_params,
										 encode_t = encode_t,
										 model_index=mi,
										 criterion=criterion,
										 contamination = contamination))

	# now loop over all anomaly score funs
	for result in results
		GenerativeAD.experiment(result..., data, savepath; save_entries...)
	end
end


text:
module TracerAdvDiff

export
   Problem,
   set_c!,
   updatevars!

using
  FFTW,
  Reexport

@reexport using FourierFlows

import LinearAlgebra: mul!, ldiv!

# --
# Problems
# --

"""
    Problem(; parameters...)

Construct a constant diffusivity problem with steady or time-varying flow.
"""
noflow(args...) = 0.0 # used as defaults for u, v functions in Problem()

function Problem(;
    nx = 128,
    Lx = 2π,
    ny = nx,
    Ly = Lx,
  grid = TwoDGrid(nx, Lx, ny, Ly),
   kap = 0.1,
   eta = kap,
     u = noflow,
     v = noflow,
    dt = 0.01,
  stepper = "RK4",
  steadyflow = false
  )

  if steadyflow; pr = ConstDiffSteadyFlowParams(eta, kap, u, v, grid)
  else;          pr = ConstDiffParams(eta, kap, u, v)
  end

  vs = Vars(grid)
  eq = Equation(pr, grid)

  FourierFlows.Problem(eq, stepper, dt, grid, vs, pr)
end

# --
# Params
# --

abstract type AbstractTracerParams <: AbstractParams end
abstract type AbstractConstDiffParams <: AbstractParams end
abstract type AbstractSteadyFlowParams <: AbstractParams end

"""
    ConstDiffParams(eta, kap, kaph, nkaph, u, v)
    ConstDiffParams(eta, kap, u, v)

Returns the params for constant diffusivity problem with time-varying flow.
"""
struct ConstDiffParams{T} <: AbstractConstDiffParams
  eta::T                 # Constant isotropic horizontal diffusivity
  kap::T                 # Constant isotropic vertical diffusivity
  kaph::T                # Constant isotropic hyperdiffusivity
  nkaph::Int             # Constant isotropic hyperdiffusivity order
  u::Function            # Advecting x-velocity
  v::Function            # Advecting y-velocity
end
ConstDiffParams(eta, kap, u, v) = ConstDiffParams(eta, kap, 0eta, 0, u, v)

"""
    ConstDiffSteadyFlowParams(eta, kap, kaph, nkaph, u, v, g)
    ConstDiffSteadyFlowParams(eta, kap, u, v, g)

Returns the params for constant diffusivity problem with time-steady flow.
"""
struct ConstDiffSteadyFlowParams{T,A} <: AbstractSteadyFlowParams
  eta::T       # Constant horizontal diffusivity
  kap::T       # Constant vertical diffusivity
  kaph::T      # Constant isotropic hyperdiffusivity
  nkaph::Int   # Constant isotropic hyperdiffusivity order
  u::A         # Advecting x-velocity
  v::A         # Advecting y-velocity
end

function ConstDiffSteadyFlowParams(eta, kap, kaph, nkaph, u::Function, v::Function, g)
  x, y = gridpoints(g)
  ugrid = u.(x, y)
  vgrid = v.(x, y)
  ConstDiffSteadyFlowParams(eta, kap, kaph, nkaph, ugrid, vgrid)
end

ConstDiffSteadyFlowParams(eta, kap, u, v, g) = ConstDiffSteadyFlowParams(eta, kap, 0eta, 0, u, v, g)

# --
# Equations
# --

"""
    Equation(p, g)

Returns the equation for constant diffusivity problem with params p and grid g.
"""
function Equation(p::ConstDiffParams, g::AbstractGrid{T}) where T
  LC = zero(g.Krsq)
  @. LC = -p.eta*g.kr^2 - p.kap*g.l^2 - p.kaph*g.Krsq^p.nkaph
  FourierFlows.Equation(LC, calcN!, g)
end

function Equation(p::ConstDiffSteadyFlowParams, g::AbstractGrid{T}) where T
  LC = zero(g.Krsq)
  @. LC = -p.eta*g.kr^2 - p.kap*g.l^2 - p.kaph*g.Krsq^p.nkaph
  FourierFlows.Equation(LC, calcN_steadyflow!, g)
end

# --
# Vars
# --

# Construct Vars types
 physicalvars = [:c, :cx, :cy]
transformvars = [:ch, :cxh, :cyh]

eval(FourierFlows.structvarsexpr(:Vars, physicalvars, transformvars))

"""
    Vars(g)

Returns the vars for constant diffusivity problem on grid g.
"""
function Vars(g; T=typeof(g.Lx))
  @createarrays T (g.nx, g.ny) c cx cy
  @createarrays Complex{T} (g.nkr, g.nl) ch cxh cyh
  Vars(c, cx, cy, ch, cxh, cyh)
end

# --
# Solvers
# --

"""
    calcN!(N, sol, t, cl, v, p, g)

Calculate the advective terms for a tracer equation with constant diffusivity and time-varying flow.
"""
function calcN!(N, sol, t, cl, v, p::AbstractConstDiffParams, g)
  @. v.cxh = im*g.kr*sol
  @. v.cyh = im*g.l*sol

  ldiv!(v.cx, g.rfftplan, v.cxh) # destroys v.cxh when using fftw
  ldiv!(v.cy, g.rfftplan, v.cyh) # destroys v.cyh when using fftw

  x, y = gridpoints(g)
  @. v.cx = -p.u(x, y, cl.t)*v.cx - p.v(x, y, cl.t)*v.cy # copies over v.cx so v.cx = N in physical space
  mul!(N, g.rfftplan, v.cx)
  nothing
end

"""
    calcN_steadyflow!(N, sol, t, cl, v, p, g)

Calculate the advective terms for a tracer equation with constant diffusivity and time-constant flow.
"""
function calcN_steadyflow!(N, sol, t, cl, v, p::AbstractSteadyFlowParams, g)
  @. v.cxh = im*g.kr*sol
  @. v.cyh = im*g.l*sol

  ldiv!(v.cx, g.rfftplan, v.cxh) # destroys v.cxh when using fftw
  ldiv!(v.cy, g.rfftplan, v.cyh) # destroys v.cyh when using fftw

  @. v.cx = -p.u*v.cx - p.v*v.cy # copies over v.cx so v.cx = N in physical space
  mul!(N, g.rfftplan, v.cx)
  nothing
end

# --
# Helper functions
# --

"""
    updatevars!(prob)

Update the vars in v on the grid g with the solution in sol.
"""
function updatevars!(prob)
  v, g, sol = prob.vars, prob.grid, prob.sol
  v.ch .= sol
  ch1 = deepcopy(v.ch)
  ldiv!(v.c, g.rfftplan, ch1)
  nothing
end

"""
    set_c!(prob, c)
    set_c!(prob, c::Function)

Set the solution sol as the transform of c and update variables v
on the grid g.
"""
function set_c!(prob, c)
  sol, v, g = prob.sol, prob.vars, prob.grid

  mul!(sol, g.rfftplan, c)
  updatevars!(prob)
  nothing
end

function set_c!(prob, c::Function)
  sol, v, g = prob.sol, prob.vars, prob.grid

  x, y = gridpoints(g)
  cgrid = c.(x, y)
  mul!(sol, g.rfftplan, cgrid)
  updatevars!(prob)
  nothing
end

end # module


text:
function read_obsheader(f::IOStream,
                obs_hstrings_set::Set{String}=obs_hstrings_set,
        )::Dict{String,String}
    header_hs = Dict{String,String}()
    if get_header!(header_hs, f, obs_hstrings_set)
        haskey(header_hs, ohsVersion) ||
            throw(error("header error, not $ohsVersion"))
        return header_hs
    end
    throw(error("read header error!"))
end

function read_obsheader(filename::String)
    open(filename, "r") do f
        return read_obsheader(f)
    end
end

function parse_obsheader(header_hs::Dict{String,String})
    Version           = parse_ohsVersion(get(header_hs, ohsVersion, ""))
    RunBy             = parse_ohsRunBy(get(header_hs, ohsRunBy, ""))
    Comment           = parse_ohsComment(get(header_hs, ohsComment, ""))
    MarkerName        = parse_ohsMarkerName(get(header_hs, ohsMarkerName, ""))
    MarkerNumber      = parse_ohsMarkerNumber(get(header_hs, ohsMarkerNumber, ""))
    MarkerType        = parse_ohsMarkerType(get(header_hs, ohsMarkerType, ""))
    Observer          = parse_ohsObserver(get(header_hs, ohsObserver, ""))
    Receiver          = parse_ohsReceiver(get(header_hs, ohsReceiver, ""))
    AntennaType       = parse_ohsAntennaType(get(header_hs, ohsAntennaType, ""))
    AntennaPosition   = parse_ohsAntennaPosition(get(header_hs, ohsAntennaPosition, ""))
    AntennaDeltaHEN   = parse_ohsAntennaDeltaHEN(get(header_hs, ohsAntennaDeltaHEN, ""))
    AntennaDeltaXYZ   = parse_ohsAntennaDeltaXYZ(get(header_hs, ohsAntennaDeltaXYZ, ""))
    AntennaPhaseCtr   = parse_ohsAntennaPhaseCtr(get(header_hs, ohsAntennaPhaseCtr, ""))
    AntennaBsightXYZ  = parse_ohsAntennaBsightXYZ(get(header_hs, ohsAntennaBsightXYZ, ""))
    AntennaZeroDirAzi = parse_ohsAntennaZeroDirAzi(get(header_hs, ohsAntennaZeroDirAzi, ""))
    AntennaZeroDirXYZ = parse_ohsAntennaZeroDirXYZ(get(header_hs, ohsAntennaZeroDirXYZ, ""))
    CenterOfMass      = parse_ohsCenterOfMass(get(header_hs, ohsCenterOfMass, ""))
    NumObs            = parse_ohsNumObs(get(header_hs, ohsNumObs, ""))
    if Version[1] > 2.99
        SystemNumObs  = parse_ohsSystemNumObs3(get(header_hs, ohsSystemNumObs, ""))
    else
        SystemNumObs  = parse_ohsSystemNumObs2(get(header_hs, ohsNumObs, ""), Version[3])
    end
    WaveFact          = parse_ohsWaveFact(get(header_hs, ohsWaveFact, ""))
    SigStrengthUnit   = parse_ohsSigStrengthUnit(get(header_hs, ohsSigStrengthUnit, ""))
    Interval          = parse_ohsInterval(get(header_hs, ohsInterval, ""))
    FirstTime         = parse_ohsFirstTime(get(header_hs, ohsFirstTime, ""))
    LastTime          = parse_ohsLastTime(get(header_hs, ohsLastTime, ""))
    ReceiverOffset    = parse_ohsReceiverOffset(get(header_hs, ohsReceiverOffset, ""))
    SystemDCBSapplied = parse_ohsSystemDCBSapplied(get(header_hs, ohsSystemDCBSapplied, ""))
    SystemPCVSapplied = parse_ohsSystemPCVSapplied(get(header_hs, ohsSystemPCVSapplied, ""))
    SystemScaleFac    = parse_ohsSystemScaleFac(get(header_hs, ohsSystemScaleFac, ""))
    SystemPhaseShift  = parse_ohsSystemPhaseShift(get(header_hs, ohsSystemPhaseShift, ""))
    GlonassSlotFreqNo = parse_ohsGlonassSlotFreqNo(get(header_hs, ohsGlonassSlotFreqNo, ""))
    GlonassCodPhsBias = parse_ohsGlonassCodPhsBias(get(header_hs, ohsGlonassCodPhsBias, ""))
    LeapSeconds       = parse_ohsLeapSeconds(get(header_hs, ohsLeapSeconds, ""))
    NumSats           = parse_ohsNumSats(get(header_hs, ohsNumSats, ""))
    PrnObs            = parse_ohsPrnObs(get(header_hs, ohsPrnObs, ""))
    EoH               = parse_ohsEoH(get(header_hs, ohsEoH, ""))
    RinexObsHeader(
        Version, RunBy, Comment,
        MarkerName, MarkerNumber, MarkerType,
        Observer, Receiver,
        AntennaType, AntennaPosition, AntennaDeltaHEN,
        AntennaDeltaXYZ, AntennaPhaseCtr, AntennaBsightXYZ,
        AntennaZeroDirAzi, AntennaZeroDirXYZ,
        CenterOfMass,
        NumObs, SystemNumObs, WaveFact,
        SigStrengthUnit, Interval,
        FirstTime, LastTime, ReceiverOffset,
        SystemDCBSapplied, SystemPCVSapplied,
        SystemScaleFac, SystemPhaseShift,
        GlonassSlotFreqNo, GlonassCodPhsBias,
        LeapSeconds, NumSats, PrnObs,
        EoH, haskey(header_hs, rchsVersion),
    )
end

function read_obsdata(oh::RinexObsHeader, f::IOStream)
    v = oh.version[1]
    sysNumObs = oh.systemNumObs
    interval = isnan(oh.interval) ? 10.0 : oh.interval
    toUTC =  oh.firstTime[2]

    recommend_length = round(Int64, interval * 160)
    if v > 2.99
        return read_3obsdata(f, sysNumObs, toUTC, recommend_length)
    else
        return read_2obsdata(f, sysNumObs, toUTC, recommend_length)
    end
end

function load_obs(filename::String)
    open(filename, "r") do f
        ohd = read_obsheader(f)
        oh  = parse_obsheader(ohd)
        if haskey(ohd, rchsVersion)
            @warn "crx, convert to rnx first!"
            return nothing
        end
        od  = read_obsdata(oh, f)
        return RinexObs(basename(filename), oh, od)
    end
end

function merge_obs(obss::Array{RinexObs,1})
    len_obss = length(obss)
    tp = Array{Tuple{DateTime,ECEF{Float64}}}(undef, len_obss)
    mheaders = Array{RinexObsHeader,1}()

    for (i, iobs) in enumerate(obss)
        firstT = iobs.header.firstTime[1]
        pos = iobs.header.antennaPosition
        tp[i] = (firstT, pos)
    end
    sort!(tp, by=x -> x[1])
    center_pos = tp[ceil(Int64, len_obss / 2)][2]

    data1 = Dict{SatID,Tuple{Array{ObsID,1},Array{SatelliteObsData,1}}}()
    lens = Dict{SatID,Int64}()
    for iobs in obss
        this_pos = iobs.header.antennaPosition
        if 10 < euclidean_distance(this_pos,
                            center_pos)
            @warn "same sation different position?:\n$center_pos; $this_pos"
            continue
        end
        push!(mheaders, iobs.header)
        for (isat, sod) in iobs.data.obs
            if haskey(data1, isat)
                intersect!(data1[isat][1], sod.obsids)
                if length(data1[isat][1]) > 4
                    lens[isat] += length(sod.dts)
                else
                    lens[isat] = 0
                end
                push!(data1[isat][2], sod)
            else
                data1[isat] = (sod.obsids, [sod])
                sizehint!(data1[isat][2], len_obss)
                lens[isat] = length(sod.dts)
                continue
            end
        end
    end
    data = Dict{SatID,SatelliteObsData}()
    for (isat, (ioids, sods)) in data1
        len = lens[isat]
        len_oids = length(ioids)
        dts = Array{DateTime,1}()
        sizehint!(dts, len)
        obs = Array{Array{RinexDatum,1},1}(undef, len_oids)
        for i in 1:len_oids
            obs[i] = Array{RinexDatum,1}()
            sizehint!(obs[i], len)
        end
        for isod in sods
            append!(dts, isod.dts)
            this_oids = isod.obsids
            this_obss = isod.obs
            for (i, ioid) in enumerate(ioids)
                j = findfirst(x -> x == ioid, this_oids)
                append!(obs[i], this_obss[j])
            end
        end
        p = sortperm(dts)
        dts = dts[p]
        for i in 1:len_oids
            obs[i] = obs[i][p]
        end
        data[isat] = SatelliteObsData(isat, ioids, dts, obs)
    end
    mheaders, data, center_pos
end

function load_obss(obs_files::Array{String,1})
    obss = Array{RinexObs,1}(undef, length(obs_files))
    for (i, this_nav_file) in enumerate(obs_files)
        obss[i] = load_obs(this_nav_file)
    end
    merge_obs(obss)
end

function get_rnx3_filename(file::String, roh::RinexObsHeader)

    bf = basename(file)
    if bf[10] == '-'
        return  bf
    end

    sn = uppercase(bf[1:4])
    mr = "00"
    cc = HTrg.reverse_geocode(roh.antennaPosition)[4][1]
    s = "R"
    ft = roh.firstTime[1] - roh.firstTime[2]
    fts = string(@sprintf("%04d",year(ft)),
                    @sprintf("%03d",dayofyear(ft)),
                    Dates.format(ft, "HHMM") )
    fp = "01D" # TODO
    fre = roh.interval > 0.99 ?
                @sprintf("%02.0fS",roh.interval) :
                "00U"
    if length(roh.systemNumObs) == 1
        type = HTgnss.mapSS2C[collect(keys(roh.systemNumObs))[1]]
    else
        type = 'M'
    end

    fmt = roh.iscrx ? "crx" : "rnx"
    filename = string(sn, mr,cc,"_",s,"_",fts,"_",fp,"_",
                        fre,"_",type,"O",".",fmt )
    # NAIN00CAN_R_20112520000_01D_30S_GO.crx.gz
    # MILP00LHT_R_20151150000_01D_01S_GO.crx
    return filename
end


text:
export Multirate

"""
    Multirate(fast, slow)

A multirate Runge--Kutta scheme, combining `fast` and `slow` algorithms

`slow` can be any algorithm providing methods for the following functions
  - `init_inner(prob, outercache, dt)`
  - `update_inner!(innerinteg, outercache, f_slow, u, p, t, dt, stage)`
Algorithms which currently support this are:
 - [`LowStorageRungeKutta2N`](@ref)
 - [`MultirateInfinitesimalStep`](@ref)
 - [`WickerSkamarockRungeKutta`](@ref)
"""
struct Multirate{F,S} <: DistributedODEAlgorithm
    fast::F
    slow::S
end

struct MultirateCache{OC,II}
    outercache::OC
    innerinteg::II
end

function cache(
    prob::DiffEqBase.AbstractODEProblem,
    alg::Multirate;
    dt, fast_dt, kwargs...)

    @assert prob.f isa DiffEqBase.SplitFunction

    # subproblems
    outerprob = DiffEqBase.remake(prob; f=prob.f.f2)
    outercache = cache(outerprob, alg.slow)

    innerfun = init_inner(prob, outercache, dt)
    innerprob = DiffEqBase.remake(prob; f=innerfun)
    innerinteg = DiffEqBase.init(innerprob, alg.fast; dt=fast_dt, kwargs...)
    return MultirateCache(outercache, innerinteg)
end

function step_u!(int, cache::MultirateCache)
    outercache = cache.outercache
    tab = outercache.tableau

    u = int.u
    p = int.p
    dt = int.dt
    t = int.t

    innerinteg = cache.innerinteg
    fast_dt = innerinteg.dt

    N = nstages(outercache)
    for stage in 1:N

        update_inner!(innerinteg, outercache, int.prob.f.f2, u, p, t, dt, stage)

        # solve inner problem
        #  dv/dτ .= B[s]/(C[s+1] - C[s]) .* du .+ f_fast(v,τ) τ ∈ [τ0,τ1]

        # TODO: make this more generic
        # there are 2 strategies we can use here:
        #  a. use same fast_dt for all slow stages, use `adjustfinal=true`
        #     - problems for ARK (e.g. requires expensive LU factorization)
        #  b. use different fast_dt, cache expensive ops

        innerinteg.adjustfinal = true
        DiffEqBase.solve!(innerinteg)
        innerinteg.dt = fast_dt # reset
    end
end


text:
using Pkg
Pkg.activate("c:/git/CuCountmap")

using CuArrays
CuArrays.allowscalar(false)
using CUDAnative

using CUDAnative: atomic_add!

using StatsBase
import StatsBase: countmap

function cucountmap!(buffer, v)
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = gridDim().x * blockDim().x

    for j = i:stride:length(v)
        b = Int(v[j]) + Int(32769)
        @atomic buffer[b] = buffer[b] + UInt32(1)
    end
    return
end

function countmap(v::CuArray{T, N, NN}; threads = 256, blocks = 1024) where {T, N, NN}
    buffer = CuArrays.zeros(UInt32, 2^16)
    CuArrays.@sync @cuda threads = threads blocks = blocks cucountmap!(buffer, v)
    buffer
end

function cucountmap2!(buffer, v)
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = gridDim().x * blockDim().x

    for j = i:stride:length(v)
         b = Int(v[j]) + Int(32769)
         @atomic buffer[b, blockIdx().x] = buffer[b, blockIdx().x] + UInt32(1)
    end

    return
end

function countmap2(v::CuArray{T, N, NN}; threads = 256, blocks = 1024) where {T, N, NN}
    buffer = CuArrays.zeros(UInt32, 2^16, blocks)
    CuArrays.@sync @cuda threads = threads blocks = blocks cucountmap2!(buffer, v)

    for i in 2:blocks
        buffer[:, 1] .+= buffer[:, i]
    end
    buffer[:, 1]
end

v = CuArray(rand(Int16, 100_000_000))

################################################################
# basic check v2
################################################################

buffer = CuArrays.zeros(UInt32, 2^16, 1024)

CuArrays.@sync @cuda threads = 512 blocks = 1024 cucountmap2!(buffer, v)

@time countmap2(v)
@time countmap(v)
vc = collect(v)

################################################################
# basic check
################################################################

buffer = CuArrays.zeros(UInt32, 2^16)

CuArrays.@sync @cuda threads = 512 blocks = 1024 cucountmap!(buffer, v)

countmap(v)
vc = collect(v)

################################################################
# detail testing
################################################################
res = @time countmap(collect(v))

buffer_check = zeros(Int, 2^16)
for (k, i) in collect(res)
    buffer_check[k+32769] = i
end

collect(buffer) == buffer_check

#@device_code_warntype cucountmap!(buffer, v)

################################################################
# benchmark
################################################################
using BenchmarkTools
vc = collect(v)

@benchmark countmap($vc)
@benchmark countmap($v)


text:
# This subroutine discretize an AR(1) process with Tauchen (86)
# λ' = (1 - ρ)μ_λ + ρλ + ϵ, where ϵ ~ N(0, (σ_ϵ)^2)
# Inputs:
# Outputs:

using Distributions

function Tauchen(μ_λ, ρ, σ_ϵ, n; m = 3)
    σ_λ = σ_ϵ/(1-ρ^2)^(1/2);

    λ_min = μ_λ - m*σ_λ;
    λ_max = μ_λ + m*σ_λ;
    λ_vals = range(λ_min, λ_max, length = n); # can also use LinRange(λ_min, λ_max, n)

    w = λ_vals[2] - λ_vals[1];

    d = Normal();
    P = zeros(n,n);

    for i in 1:n
        for j = 1:n
            if j == 1
                P[i,j] = cdf(d, (λ_vals[1] + w/2 - (1-ρ)*μ_λ - ρ*λ_vals[i])/σ_ϵ);
            elseif j == n
                P[i,j] = 1 - cdf(d, (λ_vals[n] - w/2 - (1-ρ)*μ_λ - ρ*λ_vals[i])/σ_ϵ);
            else
                P[i,j] = cdf(d, (λ_vals[j] + w/2 - (1-ρ)*μ_λ - ρ*λ_vals[i])/σ_ϵ) - cdf(d, (λ_vals[j] - w/2 - (1-ρ)*μ_λ - ρ*λ_vals[i])/σ_ϵ);
            end

        end
    end
    return λ_vals, P
end


text:
function derivatives(y)

# evaluates the matrices for the first and second derivatives
# with respect to y
#
#
# input
#  y -  meridional coordinate

#
# output
# d1, d2 - matrices for first and second derivative

dy = y[2]-y[1];
ngrid = length(y);
a = Vector{Float64}(undef,ngrid);
b = Vector{Float64}(undef,ngrid-1);
c = Vector{Float64}(undef,ngrid-1);

#
# second derivative
a[1:ngrid] .= -2.0/dy^2;
b[1:ngrid-1] .= 1.0/dy^2;
d2 = SymTridiagonal(a,b);

# first derivative
c[1:ngrid-1] .= 1.0/(2*dy);
d1 = Tridiagonal(-c,zeros(ngrid),c);

return d1,d2

end

max_stars_repo_name: bbuffett/MagRossby

max_issues_repo_name: bbuffett/MagRossby

max_forks_repo_name: bbuffett/MagRossby


text:
printStackTrace(ex) = Base.showerror(stdout, ex, catch_backtrace())

println("Start Server...")
while true
  try
    run(`java -jar cywwtaip.jar`)
  catch ex
    printStackTrace(ex)
  end
  println("Start next session in 3 sec...")
  sleep(3)
end

max_stars_repo_name: Gilga/AI-Cywwtaip

max_issues_repo_name: Gilga/AI-Cywwtaip

max_forks_repo_name: Gilga/AI-Cywwtaip


text:
export SimpleDeterministicCompositeModel

using MLJBase

"""
    SimpleDeterministicCompositeModel(;regressor=ConstantRegressor(),
                              transformer=FeatureSelector())

Construct a composite model consisting of a transformer
(`Unsupervised` model) followed by a `Deterministic` model. Mainly
intended for internal testing .

"""
mutable struct SimpleDeterministicCompositeModel{L<:Deterministic,
                             T<:Unsupervised} <: DeterministicComposite
    model::L
    transformer::T

end

function SimpleDeterministicCompositeModel(;
                      model=DeterministicConstantRegressor(),
                      transformer=FeatureSelector())

    composite =  SimpleDeterministicCompositeModel(model, transformer)

    message = MLJBase.clean!(composite)
    isempty(message) || @warn message

    return composite

end

MLJBase.is_wrapper(::Type{<:SimpleDeterministicCompositeModel}) = true

function MLJBase.fit(composite::SimpleDeterministicCompositeModel,
                     verbosity::Integer, Xtrain, ytrain)
    X = source(Xtrain) # instantiates a source node
    y = source(ytrain)

    t = machine(composite.transformer, X)
    Xt = transform(t, X)

    l = machine(composite.model, Xt, y)
    yhat = predict(l, Xt)

    mach = machine(Deterministic(), X, y; predict=yhat)
    fit!(mach, verbosity=verbosity)

    return mach()
end

MLJBase.load_path(::Type{<:SimpleDeterministicCompositeModel}) =
    "MLJBase.SimpleDeterministicCompositeModel"
MLJBase.package_uuid(::Type{<:SimpleDeterministicCompositeModel}) = ""
MLJBase.package_url(::Type{<:SimpleDeterministicCompositeModel}) =
    "https://github.com/alan-turing-institute/MLJBase.jl"
MLJBase.is_pure_julia(::Type{<:SimpleDeterministicCompositeModel}) = true
MLJBase.input_scitype(::Type{<:SimpleDeterministicCompositeModel{L,T}}) where {L,T} =
    MLJBase.input_scitype(T)
MLJBase.target_scitype(::Type{<:SimpleDeterministicCompositeModel{L,T}}) where {L,T} =
    MLJBase.target_scitype(L)


text:
@testset "NoiseGrid" begin

using DiffEqNoiseProcess, DiffEqBase, Test

t = 0:0.001:1
grid = exp.(t)
W = NoiseGrid(t,grid)

dt = 0.1
calculate_step!(W,dt,nothing,nothing)

for i in 1:10
  accept_step!(W,dt,nothing,nothing)
end

W = NoiseGrid(t,grid)
prob = NoiseProblem(W,(0.0,1.0))
sol = solve(prob;dt=0.1)

@test !sol.step_setup
@test_throws ErrorException accept_step!(sol,dt,nothing,nothing)

t = 0:0.001:1
grid = [[exp.(ti) for i in 1:8] for ti in t]
W = NoiseGrid(t,grid)
prob = NoiseProblem(W,(0.0,1.0))
sol = solve(prob;dt=0.1)

dt = 0.001
t = 0:dt:1
brownian_values = cumsum([0;[sqrt(dt)*randn() for i in 1:length(t)-1]])
W = NoiseGrid(t,grid)

end


text:
abstract type AbstractGrid{N,T,Names} end
struct LatticeGrid{N} end
struct ParticleGrid end

include("latticegrids.jl")
include("particlegrids.jl")
include("abstractgrid.jl")
include("utils.jl")
include("units.jl")


text:
#=
test:
- Julia version:
- Author: bramb
- Date: 2019-02-05
=#
import Pkg; Pkg.add("PyCall")

using PyCall

# start julia in HRGParser/src, where TAGMLLexer.py can be found
pushfirst!(PyVector(pyimport("sys")["path"]), @__DIR__) # add current directory to python path

Pkg.add("Conda")
using Conda

Conda.add_channel("conda-forge")
Conda.add("antlr-python-runtime")

function __init__()
    py"""
    from antlr4 import *
    from TAGMLLexer import TAGMLLexer

    def tokenize(filename):
        file = open(filename, "r")
        text = file.read()
        file.close()
        lexer = TAGMLLexer(InputStream(text))
        return CommonTokenStream(lexer)
    """
end

cs = py"tokenize"("test.tag")
keys(cs) # shows available functions on the commontokenstream PyObject
token0 = cs[:get](0)
token1 = cs[:get](1)
cs[:LT](-2)[:text]

# walk the stream until the end
cs[:fill]()
t=cs[:getTokens](0,25)
map(x->x[:text],t)

max_stars_repo_name: brambg/HyperGraphTools

max_issues_repo_name: brambg/HyperGraphTools

max_forks_repo_name: brambg/HyperGraphTools


text:
function get_indentation(nspaces::Int)
    res = ""
    for i = 1:nspaces
        res *= " "
    end
    return res
end

function is_leafnode(obj)
    if isa(obj, Array{Any, 1})
        n = length(obj)
        non_dict = true
        i = 1
        while non_dict
            if !isa(obj[i], String) && !(typeof(obj[i]) <: Number)
                non_dict = false
                res = false
            elseif i == n
                res = true
                break
            end
            i += 1
        end
    else
        res = false
    end
    return res
end

is_multidict(obj) = isa(obj, MultiDict)

function any_multidicts(v::Array{Any, 1})
    res = any(x -> isa(x, MultiDict), v)
    return res
end

function show_key_structure(xmldict_obj, nspaces = 4)
    if is_multidict(xmldict_obj)
        keys_array = collect(keys(xmldict_obj))
        indent_str = get_indentation(nspaces)
        for k in keys_array
            println(indent_str, "-", k)
            if is_multidict(xmldict_obj[k]) || typeof(xmldict_obj[k]) == Array{Any, 1}
                num_elem = length(xmldict_obj[k])
                for i = 1:num_elem
                    show_key_structure(xmldict_obj[k][i], nspaces + 4)
                    if i < num_elem
                        println(indent_str, "-", k)
                    end
                end
            end
        end
    end
end


text:
log_range(lb,ub,n) = exp.(range(log(lb), log(ub); length = n))

quadratic_range(ub,n; reverse=false) = let f(x) = sign(x) * x^2 * ub
    l, r = reverse ? (1,-1) : (-1,1)
    f.(range(l, r; length=n))
end

function build_ui(max_v, max_ψ, formation_names)
    scene, layout = layoutscene(resolution = (2000, 1600))

    function mk_silder!(name, range, parent_layout, unit = "")
        ls = labelslider!(scene, name, range, format = x -> @sprintf("%.3g%s", x, unit))
        parent_layout[end+1,1] = ls.layout
        ls.slider
    end

    ax_view = layout[1,1] = LAxis(
        scene, title = "Scene", aspect = DataAspect(),
        backgroundcolor = RGBf0(0.98, 0.98, 0.98))
    display_layout = layout[1,2] = GridLayout()
    controls_layout = layout[2,1] = GridLayout()
    params_layout = layout[2,2] = GridLayout()
    colsize!(layout, 1, Relative(2/3))
    rowsize!(layout, 1, Relative(4/5))

    dy_noise_control = mk_silder!(
        "dynamics noise", log_range(1e-5, 1e-1, 100), params_layout)
    sensor_noise_control = mk_silder!(
        "sensor noise", log_range(1e-5, 1e-1, 100), params_layout)

    max_ψ = rad2deg(max_ψ)
    slider_ticks = 500
    v_slider = mk_silder!(
        "v̂", range(-max_v, max_v, length = slider_ticks), controls_layout, "m/s")
    ψ_slider = mk_silder!(
        "ψ̂", quadratic_range(max_ψ, slider_ticks, reverse = true), controls_layout, "°")
    formation_control = controls_layout[0, 1] = LMenu(
        scene, direction = :up,
        options = formation_names, selection = formation_names[1])

    function handle_keyboard()
        function change_turn(delta)
            set_close_to!(ψ_slider, ψ_slider.value[] + 4delta * max_ψ / slider_ticks)
        end
        function change_speed(delta)
            set_close_to!(v_slider, v_slider.value[] + 4delta * max_v / slider_ticks)
        end
        function change_formation(i)
            formation_control.i_selected[] = i
        end

        let button = scene.events.keyboardbuttons[]
            ispressed(button, Keyboard.left) && change_turn(1)
            ispressed(button, Keyboard.right) && change_turn(-1)
            ispressed(button, Keyboard.up) && change_speed(1)
            ispressed(button, Keyboard.down) && change_speed(-1)
            ispressed(button, Keyboard._1) && change_formation(1)
            ispressed(button, Keyboard._2) && change_formation(2)
            ispressed(button, Keyboard._3) && change_formation(3)
        end
    end

    (
        scene = scene,
        ax_view = ax_view,
        display_layout = display_layout,
        dy_noise = dy_noise_control,
        sensor_noise = sensor_noise_control,
        v_slider = v_slider,
        ψ_slider = ψ_slider,
        formation_control = formation_control,
        handle_keyboard = handle_keyboard,
    )
end

function tryset_ylims(scene, (lower, upper), margin = 0.1)
    Δ = upper - lower
    (Δ != 0) && ylims!(scene, (lower - margin * Δ, upper + margin * Δ))
end

function draw_view(ui, freq, init, central, dy_model, form_from_id)
    @unpack ax_view, display_layout, formation_control, scene = ui

    N = length(init)
    xs_node = Node((p -> p[1]).(init))
    traj_len = round_ceil(freq * 5)
    traj_qs = [Node(constant_queue(init[i][1], traj_len)) for i in 1:N]
    colors = [get(ColorSchemes.Spectral_10, i/N) for i in 1:N]

    loss_node = Node(AxisArray(zeros(ℝ, N, 2), id=1:N, comp=[:u, :x]))
    loss_len = traj_len
    loss_q = constant_queue(loss_node[], loss_len)
    u_loss = [Node(zeros(ℝ, loss_len)) for _ in 1:N]
    x_loss = [Node(zeros(ℝ, loss_len)) for _ in 1:N]
    total_loss = Node(zeros(ℝ, loss_len))

    # draw loss history
    on(loss_node) do loss
        pushpop!(loss_q, copy(loss))
        for i in 1:N
            u_loss[i][] = [l[comp=:u, id = i] for l in loss_q]
            x_loss[i][] = [l[comp=:x, id = i] for l in loss_q]
        end
        total_loss[] = map(sum, loss_q)
    end

    function loss_ax(title)
        display_layout[end+1,1] = LAxis(scene, title = title,
            xzoomlock = true, xpanlock = true, yrectzoom = true)
    end
    ax_total, ax_x, ax_u = loss_ax.(["Total Loss", "X Loss", "U Loss"])
    trim!(display_layout)

    for i in 1:N
        lines!(ax_x, 1:loss_len, x_loss[i]; color=colors[i])
        lines!(ax_u, 1:loss_len, u_loss[i]; color=colors[i])
    end
    lines!(ax_total, 1:loss_len, total_loss; color=:red)
    on(total_loss) do loss
        tryset_ylims(ax_total, (minimum(loss), maximum(loss)))
        get_min(losses) = minimum(n -> minimum(n[]), losses)
        get_max(losses) = maximum(n -> maximum(n[]), losses)
        tryset_ylims(ax_x, (get_min(x_loss), get_max(x_loss)))
        tryset_ylims(ax_u, (get_min(u_loss), get_max(u_loss)))
    end

    function draw_refs(ctrl::FormationControl{RefPointTrackControl})
        refpoints = @lift let
            leader = $xs_node[1]
            @_ (form_from_id(formation_control.i_selected[])
                |> map(ref_point(ctrl.K, _),__)
                |> map(to_formation_frame(ctrl, leader), __))
        end
        points = @lift (x -> Point2f0(x[1])).($refpoints)
        directions = @lift (x -> Point2f0(x[2])).($refpoints)
        arrows!(ax_view, points, directions;
            linecolor=:red, arrowcolor=:red, lengthscale=0.2, arrowsize=0.1)
        scatter!(ax_view, points; color=:red)
    end

    function draw_refs(ctrl::FormationControl{ConfigTrackControl})
        ref_cars = @lift let
            leader = $xs_node[1]
            center, ω = center_of_rotation(leader, dy_model.l)
            form = form_from_id(formation_control.i_selected[])
            s′ = form[1]
            rot = rotation2D(leader.θ - s′.θ)
            offset = get_pos(leader) - get_pos(s′)
            transform = s -> rot * s + offset
            @_ (form
                |> map(transform ∘ get_pos, __)
                |> map(rotate_around(center, ω, dy_model.l, _), __))
        end
        for id in 1:N
            car = @lift let
                @unpack x, y, θ = $ref_cars[id]
                car_triangle(x, y, θ)
            end
            poly!(ax_view, car; color = with_alpha(colors[id], 0.3))
        end
    end

    # draw trajectories
    for id in 1:N
        q = traj_qs[id]
        traj_x = @lift [x.x for x in $q]
        traj_y = @lift [x.y for x in $q]
        on(xs_node) do xs
            pushpop!(q[], xs[id])
            q[] = q[]
        end
        c = colors[id]
        cmap = [with_alpha(c, 0.4), with_alpha(c, 1)]
        color = range(0, 1, length = traj_len)
        lines!(ax_view, traj_x, traj_y; color, colormap = cmap, linewidth=3)
    end
    # draw cars
    for id in 1:N
        car = @lift let
            @unpack x, y, θ = $xs_node[id]
            car_triangle(x, y, θ)
        end
        poly!(ax_view, car; color = colors[id], strokecolor=:black, strokewidth=2)
    end
    # draw ref points
    draw_refs(central)

    display(ui.scene)
    xs_node, loss_node
end

function live_demo(; CF::CFName = onevision_cf, seed = 1234)
    X, U = CarX{ℝ}, CarU{ℝ}
    Z = HVec{U, ℕ}
    ΔT = 5
    freq = 20.0 * ΔT
    delta_t = 1 / freq

    max_v = 5.0
    max_ψ = 45°
    formation_names = ["triangle", "horizontal", "vertical"]
    ui = build_ui(max_v, max_ψ, formation_names)

    # Car dynamics parameters
    dy_model = CarDynamics(;delta_t, max_ψ = 45°, integrator_samples = round_ceil(5 / ΔT))
    rng = MersenneTwister(seed)
    function add_noise(x::X, t)::X where X
        x + ui.dy_noise.value[] * CarX(
            x=0.0, y=0.0, θ=0.0, v = randn(rng, ℝ), ψ = randn(rng, ℝ))
    end
    dy_actual = @set dy_model.add_noise = add_noise
    dy_actual_leader = dy_model

    function xs_observer(xs, t)
        (x -> x + randn(rng, X) * ui.sensor_noise.value[]).(xs)
    end

    N = 4
    triangle_formation = let
        l = 0.8
        Δϕ = 360° / (N-1)
        circle = [X(x = l * cos(Δϕ * i), y = l * sin(Δϕ * i), θ = 0) for i in 1:N-1]
        [[zero(X)]; circle]
    end
    horizontal_formation = let
        l = 0.5
        leader_idx = round_ceil(N/2)
        line = [X(x = 0, y = l * (i - leader_idx), θ = 0) for i in 1:N]
        [line[mod1(leader_idx + j - 1, N)] for j in 1:N]
    end
    vertical_formation = let
        l = 0.6
        leader_idx = round_ceil(N/2)
        line = [X(x = l * (i - leader_idx), y = 0, θ = 0) for i in 1:N]
        [line[mod1(leader_idx + j - 1, N)] for j in 1:N]
    end
    formation_map = Dict{String, Formation{ℝ}}(zip(formation_names,
        [triangle_formation, horizontal_formation, vertical_formation]))
    form_from_id(i) = formation_map[formation_names[i]]

    function external_control(x, t)
        CarU{ℝ}(ui.v_slider.value[], deg2rad(ui.ψ_slider.value[]))
    end
    function external_formation(x, t)
        ui.formation_control.i_selected[]
    end
    leader_z_dy = FormationObsDynamics(external_control, external_formation)

    delays_model  = DelayModel(obs = 2, act = 4, com = 6, ΔT = ΔT)
    # delays_model  = DelayModel(obs = 0, act = 0, com = 4ΔT, ΔT = ΔT)
    delays_actual = delays_model
    H = 20

    RefK = RefPointTrackControl(;
        dy = dy_model, ref_pos = dy_model.l, ctrl_interval = delta_t * ΔT,
        kp = 1.0, ki = 0, kd = 0)
    # RefK = ConfigTrackControl(;dy = dy_model, k1 = 2.0, k2 = 0.5, k3 = 1.0)
    avoidance = CollisionAvoidance(scale=1.0, min_r=dy_model.l, max_r=3*dy_model.l)
    central = FormationControl((_, zs, _) -> form_from_id(zs[1].d),
        RefK, dy_model, avoidance)

    init = let
        formation = rotate_formation(form_from_id(1), 0°)
        map(1:N) do i
            x = formation[i]
            x, HVec(zero(U), 1), zero(U)
        end
    end

    world_model = WorldDynamics(fill((dy_model, StaticObsDynamics()), N))
    x_weights = SVector{N}(fill(X(x = 1, y = 1, θ = 1), N))
    u_weights = SVector{N}(fill(U(v̂ = 1, ψ̂ = 1), N))
    loss_model = RegretLossModel(central, world_model, x_weights, u_weights)

    framework = mk_cf(CF, world_model, central, delays_model, loss_model; X, Z, H)

    world = WorldDynamics([
        (dy_actual_leader, leader_z_dy)
        fill((dy_actual, StaticObsDynamics()), N-1)
    ])

    xs_node, loss_node = draw_view(ui, freq, init, central, dy_model, form_from_id)
    start_time = Dates.now()
    last_time = Dates.now()
    function callback(xs,zs,us,loss,t)
        xs_node[] = xs
        loss_node[] = loss
        ui.handle_keyboard()

        current_time = Dates.now()
        time = Dates.canonicalize(current_time - start_time)
        print("Interactive Simulation running ($time)...  \r"); flush(stdout)

        to_sleep = delta_t - (current_time - last_time).value/1000
        sleep(max(0, to_sleep))
        last_time = current_time
    end
    sleep(1) # wait for the UI to appear
    simulate(
        world, delays_actual, framework, init, callback, (1, typemax(𝕋));
        loss_model, xs_observer)
end


text:
using transD_GP.MT1D, PyPlot
##
ntimesperdecade = 10
T = 10 .^(-3:1/ntimesperdecade:6)
f = 1 ./T
## model 1 conductor above
ρ1 = [10, 100]
z = [0, 5000]
h = diff(z)
Z = MT1D.Z_f(f, ρ1, h)
fig = MT1D.plotcurve(T, Z)
## model 2 conductor below
ρ2 = [100, 10]
Z = MT1D.Z_f(f, ρ2, h)
MT1D.plotcurve(T, Z, fig)
## and to do this with the models plotted
fig = MT1D.plotmodelcurve(T, ρ1, z)
MT1D.plotmodelcurve(T, ρ2, z, fig)

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA