text:
using Statistics
using Plots
using Flux
using OceanParameterizations
using Oceananigans.Grids
using BSON
using OrdinaryDiffEq, DiffEqSensitivity
using WindMixing
using JLD2
using FileIO
using WindMixing: prepare_parameters_NDE_training
using WindMixing: prepare_BCs
using BenchmarkTools
using LinearAlgebra

PATH = joinpath(pwd(), "extracted_training_output")
DATA_NAME = "NDE_training_mpp_10sim_windcooling_windheating_diffusivity_1e-1_Ri_1e-1_weights_divide1f5_gradient_smallNN_scale_5e-3_rate_1e-4"
DATA_PATH = joinpath(PATH, "$(DATA_NAME)_extracted.jld2")

# file = jldopen(DATA_PATH, "r")
file = jldopen(DATA_PATH, "r")

train_files = file["training_info/train_files"]
train_parameters = file["training_info/parameters"]

𝒟train = WindMixing.data(train_files, scale_type=ZeroMeanUnitVarianceScaling, enforce_surface_fluxes=false)

uw_NN = file["neural_network/uw"]
vw_NN = file["neural_network/vw"]
wT_NN = file["neural_network/wT"]

close(file)
[uw_NN(rand(96)) uw_NN(rand(96))]

constants, scalings, derivatives, NN_constructions, weights, NN_sizes, NN_ranges, filters = prepare_parameters_NDE_training(
    𝒟train, uw_NN, vw_NN, wT_NN, 1f-4, 32, 9.81, 1.67f-4, 1f-4, 1f-1, 0.25f0, 0.1f0, 1f0, 10, (modified_pacanowski_philander=true, something=false))

test_files = ["-1e-3"]
𝒟test = WindMixing.data(test_files, scale_type=ZeroMeanUnitVarianceScaling, enforce_surface_fluxes=false)
BCs = prepare_BCs(𝒟test, scalings)
uvT₀ = [scalings.u.(𝒟test.u.coarse[:,1]); scalings.v.(𝒟test.v.coarse[:,1]); scalings.T.(𝒟test.T.coarse[:,1])]
trange = 1:1:1153
ts = 𝒟test.t[trange] ./ constants.τ
timestepper = ImplicitEuler(autodiff=false)

Threads.nthreads() = 6
BLAS.set_num_threads(4)

@info "CPU"
sol = solve_NDE_mutating(uw_NN, vw_NN, wT_NN, scalings, constants, BCs, derivatives, uvT₀, ts, timestepper)
# @btime solve_NDE_mutating(uw_NN, vw_NN, wT_NN, scalings, constants, BCs, derivatives, uvT₀, ts, timestepper)
# @btime solve_NDE_mutating(uw_NN, vw_NN, wT_NN, scalings, constants, BCs, derivatives, uvT₀, ts, timestepper)

@info "CPU 40 threads"
BLAS.set_num_threads(40)
# sol = solve_NDE_mutating(uw_NN, vw_NN, wT_NN, scalings, constants, BCs, derivatives, uvT₀, ts, timestepper=timestepper)
@btime solve_NDE_mutating(uw_NN, vw_NN, wT_NN, scalings, constants, BCs, derivatives, uvT₀, ts, timestepper=timestepper)
@btime solve_NDE_mutating(uw_NN, vw_NN, wT_NN, scalings, constants, BCs, derivatives, uvT₀, ts, timestepper=timestepper)

# plot(sol[65:end, end], -31:1:0)

using CUDA
using WindMixing: solve_NDE_mutating_GPU
CUDA.allowscalar(false)

test_files = ["-1e-3"]
𝒟test = WindMixing.data(test_files, scale_type=ZeroMeanUnitVarianceScaling, enforce_surface_fluxes=false)
BCs = prepare_BCs(𝒟test, scalings)
uvT₀_gpu = [scalings.u.(𝒟test.u.coarse[:,1]); scalings.v.(𝒟test.v.coarse[:,1]); scalings.T.(𝒟test.T.coarse[:,1])] |> gpu
trange = 1:1:1153
ts = 𝒟test.t[trange] ./ constants.τ
tspan_gpu = (ts[1], ts[end])
ts_gpu = ts
timestepper=ROCK4()

file = jldopen(DATA_PATH, "r")
uw_NN = file["neural_network/uw"] |> gpu
vw_NN = file["neural_network/vw"] |> gpu
wT_NN = file["neural_network/wT"] |> gpu
close(file)

uw_NN_gpu = uw_NN |> gpu
vw_NN_gpu = vw_NN |> gpu
wT_NN_gpu = wT_NN |> gpu

derivatives_gpu = (face=derivatives.face |> gpu, cell=derivatives.cell |> gpu)

@info "GPU"
# sol_gpu = solve_NDE_mutating_GPU(uw_NN_gpu, vw_NN_gpu, wT_NN_gpu, scalings, constants, BCs, derivatives_gpu, uvT₀_gpu, ts_gpu, tspan_gpu, timestepper=timestepper)
@btime solve_NDE_mutating_GPU(uw_NN_gpu, vw_NN_gpu, wT_NN_gpu, scalings, constants, BCs, derivatives_gpu, uvT₀_gpu, ts_gpu, tspan_gpu, timestepper=timestepper)
@btime solve_NDE_mutating_GPU(uw_NN_gpu, vw_NN_gpu, wT_NN_gpu, scalings, constants, BCs, derivatives_gpu, uvT₀_gpu, ts_gpu, tspan_gpu, timestepper=timestepper)


text:
using Random
import SharedArrays
import Base.Threads.@threads

function execute!(model::GebModel)
    # log births of the initial organisms in the first execution step
    if model.time == 0
        for organism in model.organisms
            log_birth(model.logger, model, organism)
        end
    end

    shuffle(model.organisms)
    organisms_in_batch = [(i, o) for (i, o) in enumerate(model.organisms)]

    activations_in_batch = SharedArrays.SharedArray{Float64}(length(organisms_in_batch), 4*MAX_NEURONS)
    activation_size_in_batch = SharedArrays.SharedArray{UInt8}(length(organisms_in_batch))

    for i in 1:length(organisms_in_batch)
        for j in 1:MAX_NEURONS
            activations_in_batch[i, j] = 0.0
        end
        activation_size_in_batch[i] = 0
    end

    for (i, organism) in organisms_in_batch
        develop_nodes!(organism.network, organism.rules)
    end

    for (i, organism) in organisms_in_batch
        activations = determine_input_activations(model, organism)
        num_outputs = length(activations)

        activations_in_batch[i, 1:num_outputs] = activations
        activation_size_in_batch[i] = num_outputs
    end

    for (i, organism) in organisms_in_batch
        fire!(organism.network, activations_in_batch[i, 1:activation_size_in_batch[i]])
    end

    for (_, organism) in organisms_in_batch
        if !(organism in model.organisms) continue end

        perform!(model, organism)
        organism.age += 1
    end

    model.time += 1
    log_step(model.logger, model)
end

ext: jl

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
"""
Module: SymbolicDiff (Symbolic Operation for Arithmetic)
"""

"""
seval(f, (dvar1, dver2), env, cache)
Return the second derivative of expr f with respect to dvar1 and dver2
"""

function seval(f, dvar::Tuple{Symbol,Symbol})
    seval(f, dvar, globalenv, SymbolicCache())
end

function seval(f, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv)
    seval(f, dvar, env, SymbolicCache())
end

function seval(f, dvar::Tuple{Symbol,Symbol}, cache::SymbolicCache)
    seval(f, dvar, globalenv, cache)
end

###

function seval(f, dvar::Tuple{SymbolicVariable{Tv},SymbolicVariable{Tv}}) where Tv
    seval(f, (dvar[1].var, dvar[2].var))
end

function seval(f, dvar::Tuple{SymbolicVariable{Tv},SymbolicVariable{Tv}}, env::SymbolicEnv) where Tv
    seval(f, (dvar[1].var, dvar[2].var), env)
end

function seval(f, dvar::Tuple{SymbolicVariable{Tv},SymbolicVariable{Tv}}, cache::SymbolicCache) where Tv
    seval(f, (dvar[1].var, dvar[2].var), cache)
end

function seval(f, dvar::Tuple{SymbolicVariable{Tv},SymbolicVariable{Tv}}, env::SymbolicEnv, cache::SymbolicCache) where Tv
        seval(f, (dvar[1].var, dvar[2].var), env, cache)
end

###

function seval(f::SymbolicValue{Tv}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    0
end

function seval(f::SymbolicVariable{Tv}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    0
end

function seval(f::AbstractNumberSymbolic{Tv}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    (dvar[1] in f.params) || (dvar[2] in f.params) || return Tv(0)
    get(cache, (f,dvar)) do
        retval = _eval(Val(f.op), f, dvar, env, cache)
        cache[(f,dvar)] = retval
    end
end

"""
_eval(::Val{xx}, dvar, f, env, cache)

Dispached function to evaluate the second derivative of f
"""

function _eval(::Val{:+}, f::SymbolicExpression{Tv}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    args = [seval(x, dvar, env, cache) for x = f.args]
    +(args...)
end

function _eval(::Val{:-}, f::SymbolicExpression{Tv}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    args = [seval(x, dvar, env, cache) for x = f.args]
    -(args...)
end

function _eval(::Val{:*}, f::SymbolicExpression{Tv}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    args = [seval(x, env, cache) for x = f.args]
    dargs_a = [seval(x, dvar[1], env, cache) for x = f.args]
    dargs_b = [seval(x, dvar[2], env, cache) for x = f.args]
    dargs_ab = [seval(x, dvar, env, cache) for x = f.args]

    ret = dargs_ab[1]
    s = args[1]
    s_a = dargs_a[1]
    s_b = dargs_b[1]
    for i = 2:length(args)
        ret *= args[i]
        ret += s * dargs_ab[i]
        ret += s_a * dargs_b[i]
        ret += s_b * dargs_a[i]
        (i == length(args)) && break
        s_a *= args[i]
        s_b *= args[i]
        s_a += s * dargs_a[i]
        s_b += s * dargs_b[i]
        s *= args[i]
    end
    ret
end

function _eval(::Val{:/}, f::SymbolicExpression{Tv}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    x,y = [seval(x, env, cache) for x = f.args]
    dx_a,dy_a = [seval(x, dvar[1], env, cache) for x = f.args]
    dx_b,dy_b = [seval(x, dvar[2], env, cache) for x = f.args]
    dx_ab,dy_ab = [seval(x, dvar, env, cache) for x = f.args]
    ((dx_ab * y - dx_a * dy_b  - dx_b * dy_a - x * dy_ab) * y + 2 * x * dy_a * dy_b) / y^3
end

function _eval(::Val{:^}, f::SymbolicExpression{Tv}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    x,y = [seval(x, env, cache) for x = f.args]
    dx_a,dy_a = [seval(x, dvar[1], env, cache) for x = f.args]
    dx_b,dy_b = [seval(x, dvar[2], env, cache) for x = f.args]
    dx_ab,dy_ab = [seval(x, dvar, env, cache) for x = f.args]
    f = x^y
    f_a = f * (dy_a * x * log(x) + y * dx_a) / x
    f_b = f * (dy_b * x * log(x) + y * dx_b) / x
    f_ab = (f_b * (dy_a * x * log(x) + y * dx_a) - f_a * dx_b + f * (dy_ab * x * log(x) + dy_a * dx_b * (1 + log(x)) + dx_a * dy_b + y * dx_ab)) / x
    f_ab
end

function _eval(::Val{:exp}, f::SymbolicExpression{Tv}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    x, = [seval(x, env, cache) for x = f.args]
    dx_a, = [seval(x, dvar[1], env, cache) for x = f.args]
    dx_b, = [seval(x, dvar[2], env, cache) for x = f.args]
    dx_ab, = [seval(x, dvar, env, cache) for x = f.args]
    exp(x) * (dx_b * dx_a + dx_ab)
end

function _eval(::Val{:log}, f::SymbolicExpression{Tv}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    x, = [seval(x, env, cache) for x = f.args]
    dx_a, = [seval(x, dvar[1], env, cache) for x = f.args]
    dx_b, = [seval(x, dvar[2], env, cache) for x = f.args]
    dx_ab, = [seval(x, dvar, env, cache) for x = f.args]
    (dx_ab * x - dx_a * dx_b) / x^2
end

function _eval(::Val{:sqrt}, f::SymbolicExpression{Tv}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    x, = [seval(x, env, cache) for x = f.args]
    dx_a, = [seval(x, dvar[1], env, cache) for x = f.args]
    dx_b, = [seval(x, dvar[2], env, cache) for x = f.args]
    dx_ab, = [seval(x, dvar, env, cache) for x = f.args]
    sqrt(x) * (dx_ab * 2 - dx_a * dx_b) / (4*x)
end

function _eval(::Val{:sum}, f::SymbolicExpression{Tv}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    dx_ab, = [seval(x, dvar, env, cache) for x = f.args]
    sum(dx_ab)
end

function _eval(::Val{:dot}, f::SymbolicExpression{Tv}, dvar::Tuple{Symbol,Symbol}, env::SymbolicEnv, cache::SymbolicCache)::Tv where Tv
    x,y = [seval(x, env, cache) for x = f.args]
    dx_a,dy_a = [seval(x, dvar[1], env, cache) for x = f.args]
    dx_b,dy_b = [seval(x, dvar[2], env, cache) for x = f.args]
    dx_ab,dy_ab = [seval(x, dvar, env, cache) for x = f.args]
    dot(x,dy_ab) + dot(dx_b,dy_a) + dot(dx_a,dy_b) + dot(dx_ab,y)
end


text:
# Global strings
function _unsafe_mj_string(A::AbstractArray{Ptr{Cchar}})
    map(A) do x
        replace(Base.unsafe_string(x), "&" => "")
    end
end

for (name, dims) in (
    (:mjDISABLESTRING, (mjNDISABLE, )),
    (:mjENABLESTRING, (mjNENABLE, )),
    (:mjTIMERSTRING, (mjNTIMER, )),
    (:mjLABELSTRING, (mjNLABEL, )),
    (:mjFRAMESTRING, (mjNFRAME, )),
    (:mjRNDSTRING, (3, mjNRNDFLAG)),
    (:mjVISSTRING, (3, mjNVISFLAG)),
)
    # TODO is this safe to do at module initialization time?
    dims = map(Int, dims)
    T = SArray{Tuple{dims...}, Ptr{Cchar}, length(dims), prod(dims)}
    cref = CRef{name, T}()
    A = _unsafe_mj_string(cref[])
    @eval const $name = $A
end

# Callbacks
for name in (
    :mjcb_passive,
    :mjcb_control,
    :mjcb_contactfilter,
    :mjcb_sensor,
    :mju_user_error,
    :mju_user_warning,
    :mjcb_time,
    :mjcb_act_dyn,
    :mjcb_act_gain,
    :mjcb_act_bias,
)
    @eval const $name = CRef{$(QuoteNode(name)), Ptr{Cvoid}}()
end

mju_user_warning_cb(msg::Cstring) = (@warn unsafe_string(msg); nothing)
mju_user_error_cb(msg::Cstring) = (@mjerrormsg msg; exit(1))

const mjCOLLISIONFUNC = let N = Int(mjNGEOMTYPES), L = N*N
    MMatrix{N,N,Ptr{Cvoid},L}(ntuple(_->C_NULL, Val(L)))
end

function _init_globals()
    warncb = @cfunction(mju_user_warning_cb, Cvoid, (Cstring,))
    errcb = @cfunction(mju_user_error_cb, Cvoid, (Cstring,))
    MJCore.mju_user_warning[] = warncb
    MJCore.mju_user_error[] = errcb

    ptr = cglobal((:mjCOLLISIONFUNC, libmujoco), typeof(mjCOLLISIONFUNC))
    A = Base.unsafe_load(ptr)
    for i in eachindex(A, mjCOLLISIONFUNC)
        mjCOLLISIONFUNC[i] = A[i]
    end
    nothing
end


text:
##
# layer 3 access
# read data as records, in the following forms:
# - plain Julia types (use isdefined to check if member is set, but can't do this check for members that are primitive type)
# - Thrift types (use Thrift.isfilled to check if member is set)
# - ProtoBuf types (use ProtoBuf.isfilled to check if member is set)
#
# Builder implementations do the appropriate translation from column chunks to types.
# RecordCursor uses Builders for the result types.

abstract type AbstractBuilder{T} end
#

##
# Row cursor iterates through row numbers of a column
mutable struct RowCursor
    par::ParFile

    rows::UnitRange{Int}            # rows to scan over
    row::Int                        # current row

    rowgroups::Vector{RowGroup}     # row groups in range
    rg::Union{Int,Nothing}          # current row group
    rgrange::Union{UnitRange{Int},Nothing} # current rowrange

    function RowCursor(par::ParFile, rows::UnitRange{Int64}, col::AbstractString, row::Signed=first(rows))
        rgs = rowgroups(par, col, rows)
        cursor = new(par, rows, row, rgs, nothing, nothing)
        setrow(cursor, row)
        cursor
    end
end

function setrow(cursor::RowCursor, row::Signed)
    cursor.row = row
    cursor.rgrange!==nothing && (row in cursor.rgrange) && return
    startrow = 1
    rgs = cursor.rowgroups
    for rg in 1:length(rgs)
        rgrange = startrow:(startrow + rgs[rg].num_rows)
        if row in rgrange
            cursor.row = row
            cursor.rg = rg
            cursor.rgrange = rgrange
            return
        end
        startrow = last(rgrange) + 1
    end
    throw(BoundsError(par.path, row))
end

rowgroup_offset(cursor::RowCursor) = cursor.row - first(cursor.rgrange)

function _start(cursor::RowCursor)
    row = first(cursor.rows)
    setrow(cursor, row)
    row
end
_done(cursor::RowCursor, row::Signed) = (row > last(cursor.rows))
function _next(cursor::RowCursor, row::Signed)
    setrow(cursor, row)
    row, (row+1)
end

function Base.iterate(cursor::RowCursor, state)
    _done(cursor, state) && return nothing
    return _next(cursor, state)
end

function Base.iterate(cursor::RowCursor)
    r = iterate(x, _start(x))
    return r
end

##
# Column cursor iterates through all values of the column, including null values.
# Each iteration returns the value (as a Union{T,Nothing}), definition level, and repetition level for each value.
# Row can be deduced from repetition level.
mutable struct ColCursor{T}
    row::RowCursor
    colname::AbstractString
    maxdefn::Int

    colchunks::Union{Vector{ColumnChunk},Nothing}
    cc::Union{Int,Nothing}
    ccrange::Union{UnitRange{Int},Nothing}

    vals::Vector{T}
    valpos::Int
    #valrange::UnitRange{Int}

    defn_levels::Vector{Int}
    repn_levels::Vector{Int}
    levelpos::Int
    levelrange::UnitRange{Int}

    function ColCursor{T}(row::RowCursor, colname::AbstractString) where T
        maxdefn = max_definition_level(schema(row.par), colname)
        new{T}(row, colname, maxdefn,nothing,nothing,nothing)
    end
end

function ColCursor(par::ParFile, rows::UnitRange{Int64}, colname::AbstractString, row::Signed=first(rows))
    rowcursor = RowCursor(par, rows, colname, row)

    rg = rowcursor.rowgroups[rowcursor.rg]
    colchunks = columns(par, rg, colname)
    ctype = coltype(colchunks[1])
    T = PLAIN_JTYPES[ctype+1]
    if (ctype == _Type.BYTE_ARRAY) || (ctype == _Type.FIXED_LEN_BYTE_ARRAY)
        T = Vector{T}
    end

    cursor = ColCursor{T}(rowcursor, colname)
    setrow(cursor, row)
    cursor
end

function setrow(cursor::ColCursor{T}, row::Signed) where {T}
    par = cursor.row.par
    rg = cursor.row.rowgroups[cursor.row.rg]
    ccincr = (row - cursor.row.row) == 1 # whether this is just an increment within the column chunk
    setrow(cursor.row, row) # set the row cursor
    cursor.colchunks!==nothing || (cursor.colchunks = columns(par, rg, cursor.colname))

    # check if cursor is done
    if _done(cursor.row, row)
        cursor.cc = length(cursor.colchunks) + 1
        cursor.ccrange = row:(row-1)
        cursor.vals = Array{T}(undef, 0)
        cursor.repn_levels = cursor.defn_levels = Int[]
        cursor.valpos = cursor.levelpos = 0
        cursor.levelrange = 0:-1 #cursor.valrange = 0:-1
        return
    end

    # find the column chunk with the row
    if cursor.ccrange===nothing || !(row in cursor.ccrange)
        offset = rowgroup_offset(cursor.row) # the offset of row from beginning of current rowgroup
        colchunks = cursor.colchunks

        startrow = row - offset
        for cc in 1:length(colchunks)
            vals, defn_levels, repn_levels = values(par, colchunks[cc])
            if isempty(repn_levels)
                nrowscc = length(vals) # number of values is number of rows
            else
                nrowscc = length(repn_levels) - length(find(repn_levels))   # number of values where repeation level is 0
            end
            ccrange = startrow:(startrow + nrowscc)

            if row in ccrange
                cursor.cc = cc
                cursor.ccrange = ccrange
                cursor.vals = vals
                cursor.defn_levels = defn_levels
                cursor.repn_levels = repn_levels
                cursor.valpos = cursor.levelpos = 0
                ccincr = false
                break
            else
                startrow = last(ccrange) + 1
            end
        end
    end

    # find the starting positions for values and levels
    ccrange = cursor.ccrange
    defn_levels = cursor.defn_levels
    repn_levels = cursor.repn_levels

    # compute the level and value pos for row
    if isempty(repn_levels)
        # no repetitions, so each entry corresponds to one full row
        levelpos = row - first(ccrange) + 1
        levelrange = levelpos:levelpos
    else
        # multiple entries may constitute one row
        idx = first(ccrange)
        levelpos = findfirst(repn_levels, 0) # NOTE: can start from cursor.levelpos to optimize, but that will prevent using setrow to go backwards
        while idx < row
            levelpos = findnext(repn_levels, 0, levelpos+1)
            idx += 1
        end
        levelend = max(findnext(repn_levels, 0, levelpos+1)-1, length(repn_levels))
        levelrange = levelpos:levelend
    end

    # compute the val pos for row
    if isempty(defn_levels)
        # all entries are required, so there must be a corresponding value
        valpos = levelpos
        #valrange = levelrange
    else
        maxdefn = cursor.maxdefn
        if ccincr
            valpos = cursor.valpos
        else
            valpos = sum(view(defn_levels, 1:(levelpos-1)) .== maxdefn) + 1
        end
        #nvals = sum(sub(defn_levels, levelrange) .== maxdefn)
        #valrange = valpos:(valpos+nvals-1)
    end

    cursor.levelpos = levelpos
    cursor.levelrange = levelrange
    cursor.valpos = valpos
    #cursor.valrange = valrange
    nothing
end

function _start(cursor::ColCursor)
    row = _start(cursor.row)
    setrow(cursor, row)
    row, cursor.levelpos
end
function _done(cursor::ColCursor, rowandlevel::Tuple{Int,Int})
    row, levelpos = rowandlevel
    (levelpos > last(cursor.levelrange)) && _done(cursor.row, row)
end
function _next(cursor::ColCursor{T}, rowandlevel::Tuple{Int,Int}) where {T}
    # find values for current row and level in row
    row, levelpos = rowandlevel
    (levelpos == cursor.levelpos) || throw(InvalidStateException("Invalid column cursor state", :levelpos))

    maxdefn = cursor.maxdefn
    defn_level = isempty(cursor.defn_levels) ? maxdefn : cursor.defn_levels[levelpos]
    repn_level = isempty(cursor.repn_levels) ? 0 : cursor.repn_levels[levelpos]
    cursor.levelpos += 1
    if defn_level == maxdefn
        val = (cursor.vals[cursor.valpos])::Union{Nothing,T}
        cursor.valpos += 1
    else
        val = (nothing)::Union{Nothing,T}
    end

    # advance row
    if cursor.levelpos > last(cursor.levelrange)
        row += 1
        setrow(cursor, row)
    end

    (val, defn_level, repn_level), (row, cursor.levelpos)
end

function Base.iterate(cursor::ColCursor, state)
    _done(cursor, state) && return nothing
    return _next(cursor, state)
end

function Base.iterate(cursor::ColCursor)
    r = iterate(x, _start(x))
    return r
end

##
# Record cursor iterates over multiple columns and returns rows as records
mutable struct RecCursor{T}
    colnames::Vector{AbstractString}
    colcursors::Vector{ColCursor}
    builder::T

    colstates::Vector{Tuple{Int,Int}}
    #colfilters::Vector{Function}
    #recfilter::Function
end

function RecCursor(par::ParFile, rows::UnitRange{Int64}, colnames::Vector{AbstractString}, builder::T, row::Signed=first(rows)) where {T <: AbstractBuilder}
    colcursors = [ColCursor(par, rows, colname, row) for colname in colnames]
    RecCursor{T}(colnames, colcursors, builder, Array{Tuple{Int,Int}}(undef, length(colcursors)))
end

function state(cursor::RecCursor)
    col1state = cursor.colstates[1]
    col1state[1] # return row as state
end

function _start(cursor::RecCursor)
    cursor.colstates = [_start(colcursor) for colcursor in cursor.colcursors]
    state(cursor)
end
_done(cursor::RecCursor, row::Signed) = _done(cursor.colcursors[1].row, row)

function _next(cursor::RecCursor{T}, row::Signed) where {T}
    states = cursor.colstates
    cursors = cursor.colcursors
    builder = cursor.builder

    row = init(cursor.builder)
    for colid in 1:length(states)                               # for each column
        colcursor = cursors[colid]
        colval, colstate = _next(colcursor, states[colid])       # for each value, defn level, repn level in column
        val, def, rep = colval
        update(builder, row, colcursor.colname, val, def, rep)  # update record
        states[colid] = colstate                                # set last state to states
    end
    row, state(cursor)
end

function Base.iterate(cursor::RecCursor, state)
    _done(cursor, state) && return nothing
    return _next(cursor, state)
end

function Base.iterate(cursor::RecCursor)
    r = iterate(cursor, _start(cursor))
    return r
end

##
# JuliaBuilder creates a plain Julia object
function default_init(::Type{T}) where {T}
    if T <: Array
        Array{eltype(T)}(0)
    else
        ccall(:jl_new_struct_uninit, Any, (Any,), T)::T
    end
end

mutable struct JuliaBuilder{T} <: AbstractBuilder{T}
    par::ParFile
    rowtype::Type{T}
    initfn::Function
    col_repeat_state::Dict{AbstractString,Int}
end
JuliaBuilder(par::ParFile, T::DataType, initfn::Function=default_init) = JuliaBuilder{T}(par, T, initfn, Dict{AbstractString,Int}())

function init(builder::JuliaBuilder{T}) where {T}
    empty!(builder.col_repeat_state)
    builder.initfn(T)::T
end

function update(builder::JuliaBuilder{T}, row::T, fqcolname::AbstractString, val, defn_level::Signed, repn_level::Signed) where {T}
    #@debug("updating $fqcolname")
    nameparts = split(fqcolname, '.')
    sch = builder.par.schema
    F = row  # the current field corresponding to the level in fqcolname
    Fdefn = 0
    Frepn = 0

    # for each name part of colname (a field)
    for idx in 1:length(nameparts)
        colname = join(nameparts[1:idx], '.')
        #@debug("updating part $colname of $fqcolname isnull:$(val === nothing), def:$(defn_level), rep:$(repn_level)")
        leaf = nameparts[idx]
        symleaf = Symbol(leaf)

        required = isrequired(sch, colname)         # determine whether field is optional and repeated
        repeated = isrepeated(sch, colname)
        required || (Fdefn += 1)                    # if field is optional, increment defn level
        repeated && (Frepn += 1)                    # if field can repeat, increment repn level

        defined = (val === nothing) ? isdefined(F, symleaf) : false
        mustdefine = defn_level >= Fdefn
        mustrepeat = repeated && (repn_level == Frepn)
        repkey = fqcolname * ":" * colname
        repidx = get(builder.col_repeat_state, repkey, 0)
        if mustrepeat
            repidx += 1
            builder.col_repeat_state[repkey] = repidx
        end
        nreps = defined ? length(getfield(F, symleaf)) : 0

        #@debug("repeat:$mustrepeat, nreps:$nreps, repidx:$repidx, defined:$defined, mustdefine:$mustdefine")
        if mustrepeat && (nreps < repidx)
            if !defined && mustdefine
                Vrep = builder.initfn(fieldtype(typeof(F), symleaf))
                setfield!(F, symleaf, Vrep)
            else
                Vrep = getfield(F, symleaf)
            end
            if length(Vrep) < repidx
                resize!(Vrep, repidx)
                if !isbits(eltype(Vrep))
                    Vrep[repidx] = builder.initfn(eltype(Vrep))
                end
            end
            F = Vrep[repidx]
        elseif !defined && mustdefine
            if idx == length(nameparts)
                V = val
            else
                V = builder.initfn(fieldtype(typeof(F), symleaf))
            end
            setfield!(F, symleaf, V)
            F = V
        end
    end
    nothing
end


text:
function show(io::IO, x::Mag)
    str = string(x)
    print(io, str)
end

function show(io::IO, x::ArbFloat{P}; midpoint::Bool=false) where {P}
    str = string(x, midpoint=midpoint)
    print(io, str)
end
showall(io::IO, x::ArbFloat{P}) where {P} = show(io, x, midpoint=true)

function show(io::IO, x::ArbReal{P}; midpoint::Bool=false, radius::Bool=false) where {P}
    str = string(x, midpoint=midpoint, radius=radius)
    print(io, str)
end
showall(io::IO, x::ArbReal{P}; radius::Bool=true) where {P} = show(io, x, midpoint=true, radius=radius)

function show(io::IO, x::ArbComplex{P}; midpoint::Bool=false, radius::Bool=false) where {P}
    str = string(x, midpoint=midpoint, radius=radius)
    print(io, str)
end
showall(io::IO, x::ArbComplex{P}; radius::Bool=true) where {P} = show(io, x, midpoint=true, radius=radius)

function show(x::ArbFloat{P}; midpoint::Bool=false) where {P}
    str = string(x, midpoint=midpoint)
    print(Base.stdout, str)
end
showall(x::ArbFloat{P}) where {P} = show(x, midpoint=true)

function show(x::ArbReal{P}; midpoint::Bool=false, radius::Bool=false) where {P}
    str = string(x, midpoint=midpoint, radius=radius)
    print(Base.stdout, str)
end
showall(x::ArbReal{P}; radius::Bool=true) where {P} = show(x, midpoint=true, radius=radius)

function show(x::ArbComplex{P}; midpoint::Bool=false, radius::Bool=false) where {P}
    str = string(x, midpoint=midpoint, radius=radius)
    print(Base.stdout, str)
end
showall(x::ArbComplex{P}; radius::Bool=true) where {P} = show(x, midpoint=true, radius=radius)


text:
### Interface to iSAM ###
### Ted Steiner ###
### January 2014   ###

# Note: To use covariances, the isam executable must be modified to
# compute these and output them to the file. One of the isam examples provides
# and example of how this can be done.
# The iSAM library is available at people.csail.mit.edu/kaess/isam/
function isam( isam_path::String,
               fname_graph::String;
               fname_result::String="NA",
               fname_cov::String="NA",
               fname_cov_rt::String="NA",
               batch::Bool=false,
               quiet::Bool=false,
               powell::Bool=true )

    args = String[]

    if batch
        push!(args,"-B")
    end
    if quiet
        push!(args,"-q")
    end
    if powell
        push!(args,"-P")
    end
    if fname_cov != "NA"
        push!(args,"-C -U $fname_cov")
    end
    if fname_cov_rt != "NA"
        push!(args,"-K $fname_cov_rt")
    end
    if fname_result != "NA"
        push!(args,"-W $fname_result")
    end

    c = `$isam_path $args $fname_graph`
    run(`$isam_path $args $fname_graph`)
    return c
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

struct Relation <: AbstractRelation
    iid::String
    type::RelationType
end

function Relation(t::Proto.Thing)
    iid = bytes2hex(t.iid)
    isempty(iid) && throw(TypeDBClientException(CONCEPT_MISSING_IID))
    return Relation(iid, instantiate(t._type))
end

# TODO seems unnecessary?
as_relation(r::Relation) = r

# Remote functions

# TODO depends on concepts: RoleType, Thing


text:
# ---
# title: 560. Subarray Sum Equals K
# id: problem560
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Hash Table
# link: <https://leetcode.com/problems/subarray-sum-equals-k/description/>
# hidden: true
# ---
#
# Given an array of integers `nums` and an integer `k`, return _the total number
# of continuous subarrays whose sum equals to`k`_.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,1,1], k = 2
#     Output: 2
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1,2,3], k = 3
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 2 * 104`
#   * `-1000 <= nums[i] <= 1000`
#   * `-107 <= k <= 107`
#
#
## @lc code=start
using LeetCode

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


text:
"""
A single point or an array of points on the cells of a Triangulation
CellField objects can be evaluated efficiently at CellPoint instances.
"""
struct CellPoint{DS,A,B,C} <: CellDatum
  cell_ref_point::A
  cell_phys_point::B
  trian::C
  domain_style::DS
end

function CellPoint(
  cell_ref_point::AbstractArray,
  trian::Triangulation,
  domain_style::ReferenceDomain)

  cell_map = get_cell_map(trian)
  cell_phys_point = lazy_map(evaluate,cell_map,cell_ref_point)
  CellPoint(cell_ref_point,cell_phys_point,trian,domain_style)
end

function CellPoint(
  cell_phys_point::AbstractArray,
  trian::Triangulation,
  domain_style::PhysicalDomain)
  cell_map = get_cell_map(trian)
  cell_invmap = lazy_map(inverse_map,cell_map)
  cell_ref_point = lazy_map(evaluate,cell_invmap,cell_phys_point)
  CellPoint(cell_ref_point,cell_phys_point,trian,domain_style)
end

function get_data(f::CellPoint)
  if DomainStyle(f) == ReferenceDomain()
    f.cell_ref_point
  else
    f.cell_phys_point
  end
end

get_triangulation(f::CellPoint) = f.trian
DomainStyle(::Type{CellPoint{DS,A,B,C}}) where {DS,A,B,C} = DS()

function change_domain(a::CellPoint,::ReferenceDomain,::PhysicalDomain)
  CellPoint(a.cell_ref_point,a.cell_phys_point,a.trian,PhysicalDomain())
end

function change_domain(a::CellPoint,::PhysicalDomain,::ReferenceDomain)
  CellPoint(a.cell_ref_point,a.cell_phys_point,a.trian,ReferenceDomain())
end

# Possibly with a different name
"""
"""
function get_cell_points(trian::Triangulation)
  cell_ref_coords = get_cell_ref_coordinates(trian)
  cell_phys_coords = get_cell_coordinates(trian)
  CellPoint(cell_ref_coords,cell_phys_coords,trian,ReferenceDomain())
end

"""
"""
abstract type CellField <: CellDatum end

function similar_cell_field(f::CellField,cell_data,trian,ds)
  GenericCellField(cell_data,trian,ds)
end

function Base.show(io::IO,::MIME"text/plain",f::CellField)
  show(io,f)
  print(io,":")
  print(io,"\n num_cells: $(num_cells(f))")
  print(io,"\n DomainStyle: $(DomainStyle(f))")
  print(io,"\n Triangulation: $(get_triangulation(f))")
  print(io,"\n Triangulation id: $(objectid(get_triangulation(f)))")
end

function CellField(f::Function,trian::Triangulation,domain_style::DomainStyle)
  s = size(get_cell_map(trian))
  cell_field = Fill(GenericField(f),s)
  GenericCellField(cell_field,trian,PhysicalDomain())
end

function CellField(f::Number,trian::Triangulation,domain_style::DomainStyle)
  s = size(get_cell_map(trian))
  cell_field = Fill(constant_field(f),s)
  GenericCellField(cell_field,trian,domain_style)
end

function CellField(f::AbstractArray{<:Number},trian::Triangulation,domain_style::DomainStyle)
  @check length(f)==num_cells(trian)  """\n
  You are trying to build a CellField from an array of length $(length(f))
  on a Triangulation with $(num_cells(trian)) cells. The length of the given array
  and the number of cells should match.
  """
  cell_field = lazy_map(constant_field,f)
  GenericCellField(cell_field,trian,domain_style)
end

function CellField(f::CellField,trian::Triangulation,domain_style::DomainStyle)
  change_domain(f,trian,domain_style)
end

function CellField(f,trian::Triangulation)
  CellField(f,trian,ReferenceDomain())
end

function get_normal_vector(trian::Triangulation)
  cell_normal = get_facet_normal(trian)
  get_normal_vector(trian,cell_normal)
end

function get_normal_vector(trian::Triangulation,cell_normal::AbstractArray)
  GenericCellField(cell_normal,trian,ReferenceDomain())
end

evaluate!(cache,f::Function,x::CellPoint) = CellField(f,get_triangulation(x))(x)

function change_domain(a::CellField,::ReferenceDomain,::PhysicalDomain)
  trian = get_triangulation(a)
  cell_map = get_cell_map(trian)
  cell_invmap = lazy_map(inverse_map,cell_map)
  cell_field_ref = get_data(a)
  cell_field_phys = lazy_map(Broadcasting(∘),cell_field_ref,cell_invmap)
  similar_cell_field(a,cell_field_phys,trian,PhysicalDomain())
end

function change_domain(a::CellField,::PhysicalDomain,::ReferenceDomain)
  trian = get_triangulation(a)
  cell_map = get_cell_map(trian)
  cell_field_phys = get_data(a)
  cell_field_ref = lazy_map(Broadcasting(∘),cell_field_phys,cell_map)
  similar_cell_field(a,cell_field_ref,trian,ReferenceDomain())
end

function change_domain(a::CellField,target_trian::Triangulation,target_domain::DomainStyle)
  change_domain(a,DomainStyle(a),target_trian,target_domain)
end

function change_domain(a::CellField,::ReferenceDomain,ttrian::Triangulation,::ReferenceDomain)
  msg = """\n
  We cannot move the given CellField to the reference domain of the requested triangulation.
  Make sure that the given triangulation is either the same as the triangulation on which the
  CellField is defined, or that the latter triangulation is the background of the former.
  """
  strian = get_triangulation(a)
  if strian === ttrian
    return a
  end
  @assert is_change_possible(strian,ttrian) msg
  D = num_cell_dims(strian)
  sglue = get_glue(strian,Val(D))
  tglue = get_glue(ttrian,Val(D))
  change_domain_ref_ref(a,ttrian,sglue,tglue)
end

function change_domain(a::CellField,::PhysicalDomain,ttrian::Triangulation,::PhysicalDomain)
  msg = """\n
  We cannot move the given CellField to the physical domain of the requested triangulation.
  Make sure that the given triangulation is either the same as the triangulation on which the
  CellField is defined, or that the latter triangulation is the background of the former.
  """
  strian = get_triangulation(a)
  if strian === ttrian
    return a
  end
  @assert is_change_possible(strian,ttrian) msg
  D = num_cell_dims(strian)
  sglue = get_glue(strian,Val(D))
  tglue = get_glue(ttrian,Val(D))
  change_domain_phys_phys(a,ttrian,sglue,tglue)
end

function change_domain(a::CellField,::PhysicalDomain,trian::Triangulation,::ReferenceDomain)
  a_trian = change_domain(a,trian,PhysicalDomain())
  change_domain(a_trian,ReferenceDomain())
end

function change_domain(a::CellField,::ReferenceDomain,trian::Triangulation,::PhysicalDomain)
  a_phys = change_domain(a,PhysicalDomain())
  change_domain(a_phys,trian,PhysicalDomain())
end

function change_domain_ref_ref(
  a::CellField,ttrian::Triangulation,sglue::FaceToFaceGlue,tglue::FaceToFaceGlue)
  sface_to_field = get_data(a)
  mface_to_sface = sglue.mface_to_tface
  tface_to_mface = tglue.tface_to_mface
  tface_to_mface_map = tglue.tface_to_mface_map
  mface_to_field = extend(sface_to_field,mface_to_sface)
  tface_to_field_s = lazy_map(Reindex(mface_to_field),tface_to_mface)
  tface_to_field_t = lazy_map(Broadcasting(∘),tface_to_field_s,tface_to_mface_map)
  similar_cell_field(a,tface_to_field_t,ttrian,ReferenceDomain())
end

function change_domain_phys_phys(
  a::CellField,ttrian::Triangulation,sglue::FaceToFaceGlue,tglue::FaceToFaceGlue)
  sface_to_field = get_data(a)
  mface_to_sface = sglue.mface_to_tface
  tface_to_mface = tglue.tface_to_mface
  mface_to_field = extend(sface_to_field,mface_to_sface)
  tface_to_field = lazy_map(Reindex(mface_to_field),tface_to_mface)
  similar_cell_field(a,tface_to_field,ttrian,PhysicalDomain())
end

"""
"""
struct GenericCellField{DS} <: CellField
  cell_field::AbstractArray
  trian::Triangulation
  domain_style::DS
  function GenericCellField(
    cell_field::AbstractArray,
    trian::Triangulation,
    domain_style::DomainStyle)

    DS = typeof(domain_style)
    new{DS}(Fields.MemoArray(cell_field),trian,domain_style)
  end
end

get_data(f::GenericCellField) = f.cell_field
get_triangulation(f::GenericCellField) = f.trian
DomainStyle(::Type{GenericCellField{DS}}) where DS = DS()
function similar_cell_field(f::GenericCellField,cell_data,trian,ds)
  GenericCellField(cell_data,trian,ds)
end

"""
   dist = distance(polytope::ExtrusionPolytope,
                   inv_cmap::Field,
                   x::Point)

Calculate distance from point `x` to the polytope. The polytope is
given by its type and by the inverse cell map, i.e. by the map from
the physical to the reference space.

Positive distances are outside the polytope, negative distances are
inside the polytope.

The distance is measured in an unspecified norm, currently the L∞
norm.
"""
function distance(polytope::ExtrusionPolytope, inv_cmap::Field, x::Point)
  extrusion = polytope.extrusion
  isempty(extrusion) && return zero(eltype(x))
  p = inv_cmap(x)
  if all(e == HEX_AXIS for e in extrusion)
    # Boundaries are at `a=0` and `a=1` in each direction
    return maximum(max(0 - a, a - 1) for a in p)
  elseif all(e == TET_AXIS for e in extrusion)
    # Calculate barycentric coordinates
    λ = Point(p..., 1 - sum(p))
    return maximum(-λ)
  else
    @notimplemented "Only hypercubes and simplices are implemented so far"
  end
end

function _point_to_cell!(cache, x::Point)
  searchmethod, kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map, table_cache = cache

  # Loop over the first m.num_nearest_vertex
  for (id,dist) in zip(knn(kdtree, SVector(Tuple(x)), searchmethod.num_nearest_vertices, true)...)

    # Find all neighbouring cells
    cells = getindex!(table_cache,vertex_to_cells,id)
    @assert !isempty(cells)

    # Calculate the distance from the point to all the cells. Without
    # round-off, and with non-overlapping cells, the distance would be
    # negative for exactly one cell and positive for all other ones. Due
    # to round-off, the distance can be slightly negative or slightly
    # positive for points near cell boundaries, in particular near
    # vertices. In this case, choose the cell with the smallest
    # distance, and check that the distance (if positive) is at most at
    # round-off level.
    T = eltype(dist)
    function cell_distance(cell::Integer)
      ctype = cell_to_ctype[cell]
      polytope = ctype_to_polytope[ctype]
      cmap = cell_map[cell]
      inv_cmap = inverse_map(cmap)
      return distance(polytope, inv_cmap, x)
    end
    # findmin, without allocating an array
    cell = zero(eltype(cells))
    dist = T(Inf)
    for jcell in cells
      jdist = cell_distance(jcell)
      if jdist < dist
        cell = jcell
        dist = jdist
      end
    end

    dist ≤ 1000eps(T) && return cell

  end

  # Output error message if cell not found
  @check false "Point $x is not inside any active cell"
end

function evaluate!(cache,f::CellField,x::Point)
  cache1,cache2 = cache
  cell_f_cache, f_cache, cell_f, f₀ = cache2
  @check f === f₀ "Wrong cache"

  cell = _point_to_cell!(cache1, x)
  cf = getindex!(cell_f_cache, cell_f, cell)
  fx = evaluate!(f_cache, cf, x)
  return fx
end

return_cache(f::CellField,xs::AbstractVector{<:Point}) = return_cache(f,testitem(xs))

# # Simple version:
# function evaluate!(cache,f::CellField,xs::AbstractVector{<:Point})
#   return map(x->evaluate!(cache,f,x), xs)
# end

function make_inverse_table(i2j::AbstractVector{<:Integer},nj::Int)
  ni = length(i2j)
  @assert nj≥0

  p = sortperm(i2j)
  # i2j[p] is a sorted array of js

  data = p
  ptrs = Array{Int}(undef, nj+1)
  i2lis = Array{Int}(undef, ni)
  jprev = zero(eltype(i2j))
  liprev = 0
  for (n,i) in enumerate(p)
    j = i2j[i]
    @assert jprev≤j≤nj
    li = (j==jprev ? liprev : 0) + 1
    ptrs[jprev+1:j] .= n
    i2lis[i] = li
    jprev = j
    liprev = li
  end
  ptrs[jprev+1:nj+1] .= ni+1
  j2is = Table(data,ptrs)

  return j2is,i2lis
end

# Efficient version:
function evaluate!(cache,f::CellField,point_to_x::AbstractVector{<:Point})
  cache1,cache2 = cache
  searchmethod, kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map = cache1
  cell_f_cache, f_cache, cell_f, f₀ = cache2
  @check f === f₀ "Wrong cache"

  ncells = length(cell_map)
  x_to_cell(x) = _point_to_cell!(cache1,x)
  point_to_cell = map(x_to_cell,point_to_x)
  cell_to_points,point_to_lpoint = make_inverse_table(point_to_cell,ncells)
  cell_to_xs = lazy_map(Broadcasting(Reindex(point_to_x)),cell_to_points)
  cell_to_f = get_array(f)
  cell_to_fxs = lazy_map(evaluate,cell_to_f,cell_to_xs)
  point_to_fxs = lazy_map(Reindex(cell_to_fxs),point_to_cell)
  point_to_fx = lazy_map(getindex,point_to_fxs,point_to_lpoint)
  collect(point_to_fx)          # Collect into a plain array
end

function compute_cell_points_from_vector_of_points(xs::AbstractVector{<:Point}, trian::Triangulation, domain_style::PhysicalDomain)
    searchmethod = KDTreeSearch()
    cache1 = _point_to_cell_cache(searchmethod,trian)
    x_to_cell(x) = _point_to_cell!(cache1, x)
    point_to_cell = map(x_to_cell, xs)
    ncells = num_cells(trian)
    cell_to_points, point_to_lpoint = make_inverse_table(point_to_cell, ncells)
    cell_to_xs = lazy_map(Broadcasting(Reindex(xs)), cell_to_points)
    cell_point_xs = CellPoint(cell_to_xs, trian, PhysicalDomain())
end

(a::CellField)(x) = evaluate(a,x)

function evaluate!(cache,f::CellField,x::CellPoint)
  _f, _x = _to_common_domain(f,x)
  cell_field = get_data(_f)
  cell_point = get_data(_x)
  lazy_map(evaluate,cell_field,cell_point)
end

function _to_common_domain(f::CellField,x::CellPoint)
  trian_f = get_triangulation(f)
  trian_x = get_triangulation(x)
  f_on_trian_x = change_domain(f,trian_x,DomainStyle(x))
  f_on_trian_x, x
end

# Gradient

function gradient(a::CellField)
  cell_∇a = lazy_map(Broadcasting(∇),get_data(a))
  if DomainStyle(a) == PhysicalDomain()
    g = cell_∇a
  else
    cell_map = get_cell_map(get_triangulation(a))
    g = lazy_map(Broadcasting(push_∇),cell_∇a,cell_map)
  end
  similar_cell_field(a,g,get_triangulation(a),DomainStyle(a))
end

function DIV(a::CellField)
  DIVa = DIV(get_data(a))
  if DomainStyle(a) == PhysicalDomain()
    @notimplemented
  end
  similar_cell_field(a,DIVa,get_triangulation(a),DomainStyle(a))
end

function ∇∇(a::CellField)
  cell_∇∇a = lazy_map(Broadcasting(∇∇),get_data(a))
  if DomainStyle(a) == PhysicalDomain()
    h = cell_∇∇a
  else
    cell_map = get_cell_map(get_triangulation(a))
    h = lazy_map(Broadcasting(push_∇∇),cell_∇∇a,cell_map)
  end
  similar_cell_field(a,h,get_triangulation(a),DomainStyle(a))
end

# This function has to be removed when ∇⋅∇(a) is implemented
laplacian(a::CellField) = tr(∇∇(a))

# Operations between CellField

function evaluate!(cache,k::Operation,a::CellField...)
  _operate_cellfields(k,a...)
end

function evaluate!(cache,k::Operation,a::Union{Function,CellField}...)
  _operate_cellfields(k,_convert_to_cellfields(a...)...)
end

function evaluate!(cache,k::Operation,a::Union{Number,CellField}...)
  _operate_cellfields(k,_convert_to_cellfields(a...)...)
end

function evaluate!(cache,k::Operation,a::Union{AbstractArray{<:Number},CellField}...)
  _operate_cellfields(k,_convert_to_cellfields(a...)...)
end

# Why julia hangs with this method????
#
#function evaluate!(cache,k::Operation,a::Union{Function,Number,AbstractArray,CellField}...)
#  b = _convert_to_cellfields(a...)
#  _operate_cellfields(k,b...)
#end

struct OperationCellField{DS} <: CellField
  op::Operation
  args::Tuple
  trian::Triangulation
  domain_style::DS
  memo::Dict{Any,Any}
  function OperationCellField(op::Operation,args::CellField...)

    @assert length(args) > 0
    trian = get_triangulation(first(args))
    domain_style = DomainStyle(first(args))
    @check all( map(i->DomainStyle(i)==domain_style,args) )
    #@check all( map(i->get_triangulation(i)===trian,args) )

    if num_cells(trian)>0
      x = _get_cell_points(args...)
      try
         ax = map(i->i(x),args)
         axi = map(first,ax)
         r = Fields.BroadcastingFieldOpMap(op.op)(axi...)
      catch
        @unreachable """\n
        It is not possible to perform the operation "$(op.op)" on the given cell fields.

        See the caught error for more information. (If you are using the Visual
          Studio Code REPL you might not see the caught error, please use the
          command-line REPL instead).
        """
      end
    end

    new{typeof(domain_style)}(op,args,trian,domain_style,Dict())
  end
end

function _get_cell_points(args::CellField...)
  k = findfirst(i->isa(i,CellState),args)
  if k === nothing
    j = findall(i->isa(i,OperationCellField),args)
    if length(j) == 0
      _get_cell_points(first(args))
    else
      _get_cell_points(args[j]...)
    end
  else
    args[k].points
  end
end

function _get_cell_points(a::CellField)
  trian = get_triangulation(a)
  get_cell_points(trian)
end

function _get_cell_points(a::OperationCellField...)
  b = []
  for ai in a
    for i in ai.args
      push!(b,i)
    end
  end
  _get_cell_points(b...)
end

function _get_cell_points(a::OperationCellField)
  _get_cell_points(a.args...)
end

function get_data(f::OperationCellField)
  a = map(get_data,f.args)
  lazy_map(Broadcasting(f.op),a...)
end
get_triangulation(f::OperationCellField) = f.trian
DomainStyle(::Type{OperationCellField{DS}}) where DS = DS()

function evaluate!(cache,f::OperationCellField,x::CellPoint)
  #key = (:evaluate,objectid(x))
  #if ! haskey(f.memo,key)
  #  ax = map(i->i(x),f.args)
  #  f.memo[key] = lazy_map(Fields.BroadcastingFieldOpMap(f.op.op),ax...)
  #end
  #f.memo[key]
  ax = map(i->i(x),f.args)
  lazy_map(Fields.BroadcastingFieldOpMap(f.op.op),ax...)
end

function change_domain(f::OperationCellField,target_trian::Triangulation,target_domain::DomainStyle)
  args = map(i->change_domain(i,target_trian,target_domain),f.args)
  OperationCellField(f.op,args...)
end

function _operate_cellfields(k::Operation,a...)
  b = _to_common_domain(a...)
  OperationCellField(k,b...)
end

function _convert_to_cellfields(a...)
  a1 = filter(i->isa(i,CellField),a)
  a2 = _to_common_domain(a1...)
  target_domain = DomainStyle(first(a2))
  target_trian = get_triangulation(first(a2))
  map(i->CellField(i,target_trian,target_domain),a)
end

function _to_common_domain(a::CellField...)

  # Find a suitable domain style
  if any( map(i->DomainStyle(i)==ReferenceDomain(),a) )
    target_domain = ReferenceDomain()
  else
    target_domain = PhysicalDomain()
  end

  # Find a suitable triangulation
  msg = """\n
  You are trying to operate CellField objects defined on incompatible triangulations.

  Make sure that all CellField objects are defined on the background triangulation
  or that the number of different sub-triangulations is equal to one.

  For instace:

  - 3 cell fields 2, two them on the same Neumann boundary and the other on the background mesh is OK.

  - 2 cell fields defined on 2 different Neumann boundaries is NOT OK.
  """
  trian_candidates = unique(objectid,map(get_triangulation,a))
  if length(trian_candidates) == 1
    target_trian = first(trian_candidates)
  elseif length(trian_candidates) == 2
    trian_a, trian_b = trian_candidates
    sa_tb = is_change_possible(trian_a,trian_b)
    sb_ta = is_change_possible(trian_b,trian_a)
    if sa_tb && sb_ta
      target_trian = best_target(trian_a,trian_b)
    elseif !sa_tb && sb_ta
      target_trian = trian_a
    elseif sa_tb && !sb_ta
      target_trian = trian_b
    else
      @unreachable msg
    end
  else
    m = """\n
    Cannote operate cellfields defined over more than 2 different
    triangulations at this moment.
    """
    @notimplemented
  end
  map(i->change_domain(i,target_trian,target_domain),a)
end

# Composition (this replaces the @law macro)

Base.:(∘)(f::Function,g::CellField) = Operation(f)(g)
Base.:(∘)(f::Function,g::Tuple{Vararg{CellField}}) = Operation(f)(g...)
Base.:(∘)(f::Function,g::Tuple{Vararg{Union{AbstractArray{<:Number},CellField}}}) = Operation(f)(g...)
Base.:(∘)(f::Function,g::Tuple{Vararg{Union{Function,CellField}}}) = Operation(f)(g...)

# Define some of the well known arithmetic ops

# Unary ops

for op in (:symmetric_part,:inv,:det,:abs,:abs2,:+,:-,:tr,:transpose,:adjoint,:grad2curl,:real,:imag,:conj)
  @eval begin
    ($op)(a::CellField) = Operation($op)(a)
  end
end

# Binary ops

for op in (:inner,:outer,:double_contraction,:+,:-,:*,:cross,:dot,:/)
  @eval begin
    ($op)(a::CellField,b::CellField) = Operation($op)(a,b)
    ($op)(a::CellField,b::Number) = Operation($op)(a,b)
    ($op)(a::Number,b::CellField) = Operation($op)(a,b)
    ($op)(a::CellField,b::Function) = Operation($op)(a,b)
    ($op)(a::Function,b::CellField) = Operation($op)(a,b)
    ($op)(a::CellField,b::AbstractArray{<:Number}) = Operation($op)(a,b)
    ($op)(a::AbstractArray{<:Number},b::CellField) = Operation($op)(a,b)
  end
end

Base.broadcasted(f,a::CellField,b::CellField) = Operation((i,j)->f.(i,j))(a,b)
Base.broadcasted(f,a::Number,b::CellField) = Operation((i,j)->f.(i,j))(a,b)
Base.broadcasted(f,a::CellField,b::Number) = Operation((i,j)->f.(i,j))(a,b)
Base.broadcasted(f,a::Function,b::CellField) = Operation((i,j)->f.(i,j))(a,b)
Base.broadcasted(f,a::CellField,b::Function) = Operation((i,j)->f.(i,j))(a,b)
Base.broadcasted(::typeof(*),::typeof(∇),f::CellField) = Operation(Fields._extract_grad_diag)(∇(f))
Base.broadcasted(::typeof(*),s::Fields.ShiftedNabla,f::CellField) = Operation(Fields._extract_grad_diag)(s(f))

dot(::typeof(∇),f::CellField) = divergence(f)
function (*)(::typeof(∇),f::CellField)
  msg = "Syntax ∇*f has been removed, use ∇⋅f (\\nabla \\cdot f) instead"
  error(msg)
end
outer(::typeof(∇),f::CellField) = gradient(f)
outer(f::CellField,::typeof(∇)) = transpose(gradient(f))
cross(::typeof(∇),f::CellField) = curl(f)

"""
    get_physical_coordinate(trian::Triangulation)

In contrast to get_cell_map, the returned object:
- is a [`CellField`](@ref)
- its gradient is the identity tensor
"""
function get_physical_coordinate(trian::Triangulation)
  CellField(_phys_coord,trian)
end

_phys_coord(x) = x

_phys_coord_grad(x) = one(typeof(outer(x,x)))

gradient(::typeof(_phys_coord)) = _phys_coord_grad

# Skeleton related Operations

function Base.getproperty(x::CellField, sym::Symbol)
  if sym in (:⁺,:plus)
    CellFieldAt{:plus}(x)
  elseif sym in (:⁻, :minus)
    CellFieldAt{:minus}(x)
  else
    getfield(x, sym)
  end
end

function Base.propertynames(x::CellField, private::Bool=false)
  (fieldnames(typeof(x))...,:⁺,:plus,:⁻,:minus)
end

struct CellFieldAt{T,F} <: CellField
  parent::F
  CellFieldAt{T}(parent::CellField) where T = new{T,typeof(parent)}(parent)
end

get_data(f::CellFieldAt) = get_data(f.parent)
get_triangulation(f::CellFieldAt) = get_triangulation(f.parent)
DomainStyle(::Type{CellFieldAt{T,F}}) where {T,F} = DomainStyle(F)
gradient(a::CellFieldAt{P}) where P = CellFieldAt{P}(gradient(a.parent))
∇∇(a::CellFieldAt{P}) where P = CellFieldAt{P}(∇∇(a.parent))
function similar_cell_field(f::CellFieldAt{T},cell_data,trian,ds) where T
  parent = similar_cell_field(f.parent,cell_data,trian,ds)
  CellFieldAt{T}(parent)
end

function CellFieldAt{T}(parent::OperationCellField) where T
  args = map(i->CellFieldAt{T}(i),parent.args)
  OperationCellField(parent.op,args...)
end

function get_normal_vector(trian::Triangulation,cell_normal::SkeletonPair)
  plus = get_normal_vector(trian,cell_normal.plus)
  minus = get_normal_vector(trian,cell_normal.minus)
  SkeletonPair(plus,minus)
end

for op in (:outer,:*,:dot)
  @eval begin
    ($op)(a::CellField,b::SkeletonPair{<:CellField}) = Operation($op)(a,b)
    ($op)(a::SkeletonPair{<:CellField},b::CellField) = Operation($op)(a,b)
  end
end

function evaluate!(cache,k::Operation,a::CellField,b::SkeletonPair{<:CellField})
  plus = k(a.plus,b.plus)
  minus = k(a.minus,b.minus)
  SkeletonPair(plus,minus)
end

function evaluate!(cache,k::Operation,a::SkeletonPair{<:CellField},b::CellField)
  plus = k(a.plus,b.plus)
  minus = k(a.minus,b.minus)
  SkeletonPair(plus,minus)
end

jump(a::CellField) = a.⁺ - a.⁻
jump(a::SkeletonPair{<:CellField}) = a.⁺ + a.⁻ # a.⁻ results from multiplying by n.⁻. Thus we need to sum.

mean(a::CellField) = Operation(_mean)(a.⁺,a.⁻)
_mean(x,y) = 0.5*x + 0.5*y

# This is the fundamental part to make operations on the skeleton work.

for fun in (:change_domain_ref_ref,:change_domain_phys_phys)
  @eval begin

    function $fun(
      a::CellField,ttrian::Triangulation,sglue::FaceToFaceGlue,tglue::SkeletonPair)
      msg = """\n
        It is not possible to use the given CellField on a SkeletonTriangulation.
        Make sure that you are specifying which of the two possible traces,
        either plus (aka ⁺) or minus (aka ⁻) you want to use.
        """
      # If the underlying array is Fill it does not matter if we restrict if from
      # the plus or minus side.
      @check isa(get_array(a.cell_field),Fill)
      plus = $fun(a,ttrian,sglue,tglue.plus)
      plus
    end

    function $fun(
      a::CellFieldAt,ttrian::Triangulation,sglue::FaceToFaceGlue,tglue::SkeletonPair)
      if isa(a,CellFieldAt{:plus})
        $fun(a,ttrian,sglue,tglue.plus)
      elseif isa(a,CellFieldAt{:minus})
        $fun(a,ttrian,sglue,tglue.minus)
      else
        @unreachable
      end
    end

    function $fun(
      f::OperationCellField,ttrian::Triangulation,sglue::FaceToFaceGlue,tglue::SkeletonPair)
      args = map(i->$fun(i,ttrian,sglue,tglue),f.args)
      OperationCellField(f.op,args...)
    end

  end
end

# Just to provide more meaningful error messages
function (a::SkeletonPair{<:CellField})(x)
  @unreachable """\n
  You are trying to evaluate a CellField on a mesh skeleton but you have not specified which of the
  two sides i.e. plus (aka ⁺) or minus (aka ⁻) you want to select.

  For instance, if you have extracted the normal vector and the cell points from a SkeletonTriangulation

      x = get_cell_points(strian)
      n = get_normal_vector(strian)

  Evaluating `n(x)` is not allowed. You need to call either `n.⁺(x)` or `n.⁻(x)`.
  """
end

# Interpolable struct
struct KDTreeSearch
  num_nearest_vertices::Int
  function KDTreeSearch(;num_nearest_vertices=1)
    new(num_nearest_vertices)
  end
end

struct Interpolable{M,A} <: Function
  uh::A
  tol::Float64
  searchmethod::M
  function Interpolable(uh; tol=1e-6, searchmethod=KDTreeSearch())
    new{typeof(searchmethod),typeof(uh)}(uh, tol,searchmethod)
  end
end

(a::Interpolable)(x) = evaluate(a,x)
evaluate!(cache,a::Interpolable,x::Point) = evaluate!(cache,a.uh,x)
return_cache(f::CellField,x::Point) = return_cache(Interpolable(f),x)

function return_cache(a::Interpolable,x::Point)
  f = a.uh
  trian = get_triangulation(f)
  cache1 = _point_to_cell_cache(a.searchmethod,trian)

  cell_f = get_array(f)
  cell_f_cache = array_cache(cell_f)
  cf = testitem(cell_f)
  f_cache = return_cache(cf,x)
  cache2 = cell_f_cache, f_cache, cell_f, f

  return cache1,cache2
end

function _point_to_cell_cache(searchmethod::KDTreeSearch,trian::Triangulation)
  model = get_active_model(trian)
  topo = get_grid_topology(model)
  vertex_coordinates = Geometry.get_vertex_coordinates(topo)
  kdtree = KDTree(map(nc -> SVector(Tuple(nc)), vertex_coordinates))
  D = num_cell_dims(trian)
  vertex_to_cells = get_faces(topo, 0, D)
  cell_to_ctype = get_cell_type(trian)
  ctype_to_reffe = get_reffes(trian)
  ctype_to_polytope = map(get_polytope, ctype_to_reffe)
  cell_map = get_cell_map(trian)
  table_cache = array_cache(vertex_to_cells)
  cache1 = searchmethod, kdtree, vertex_to_cells, cell_to_ctype, ctype_to_polytope, cell_map, table_cache
end


text:
# ---------------------------------- Type of datasets ---------------------------------- */

abstract type DatasetName end
abstract type Tabular <: DatasetName end

abstract type ImageType end
abstract type MLImage <: ImageType end

abstract type Image{T <: ImageType} <: DatasetName end
abstract type GrayImage{T <: ImageType} <: Image{T} end
abstract type ColorImage{T <: ImageType} <: Image{T} end

# -------------------------------------------------------------------------------------- */
abstract type File end
abstract type Labels <: File end
abstract type Header <: File end

abstract type Problem end
abstract type Classification <: Problem end
abstract type Regression <: Problem end

abstract type Normalization end
abstract type L2 <: Normalization end
abstract type Std <: Normalization end
abstract type MinMax <: Normalization end

abstract type Split end
abstract type Train <: Split end
abstract type Valid <: Split end
abstract type Test <: Split end

max_stars_repo_name: GruAna/PreprocessData

max_issues_repo_name: GruAna/PreprocessData

max_forks_repo_name: GruAna/PreprocessData


text:
import TPT

using Base.Test

using Calculus
using Plots
pyplot()

println("--- Unit AHS ---")

#
# Prepare a directory to store the results
#
!isdir("results") && mkdir("results")

resdir = joinpath("results", "unit_ahs")
!isdir(resdir) && mkdir(resdir)

# default settings for plots
default(xlabel = "r", ylabel = "g(r)", label = "")

#
# Tests of cavity functions
# Ref: E. W. Grundke and D. Henderson: Mol. Phys. 24 (1972) 269-281.
#
ahs = TPT.AHS(η = 0.5068, σ = [1.0, 3.0], c = [0.5, 0.5])
ρ₀ = TPT.totalnumberdensity(ahs)
y = TPT.cavityfunction(ahs)
y1 = y[1,1](0)
y2 = y[2,2](0)
logy = [ r -> log10(y[i,j](r)) for i in 1:2, j in 1:2 ]
plot([logy[1,1], logy[1,2], logy[2,2]], 0, 5, ylims=(0, 5))
png(joinpath(resdir, "y"))

#
# A binary hard-sphere mixture
# Ref: S. B. Yuste et al. J. Chem. Phys. 108 (1998) 3683-3693
#
ahs1 = TPT.AHS(η = 0.49, σ = [1.0, 0.3], c = [1/16, 15/16], approx=:RFA)

α1 = ahs1.α
g1 = TPT.paircorrelation(ahs1)

# Fig.1(a)
plot(g1[1,1], 1, 1.25, ylims = (0, 10), yticks = 0:1:10)
png(joinpath(resdir, "g11a"))

# Fig.1(b)
plot(g1[1,1], 1, 2.25, xticks = 1:0.2:2.2, ylims = (0.55, 1.25),
     yticks = 0.0:0.1:1.3)
png(joinpath(resdir, "g11b"))

# Fig.2(a)
plot(g1[1,2], 0.65, 0.90, ylims = (0.8, 4.8), yticks = 0:0.5:10)
png(joinpath(resdir, "g12a"))

# Fig.2(b)
plot(g1[1,2], 0.7, 1.9, xticks = (0:0.2:2), ylims = (0.88, 1.18),
     yticks = 0:0.05:2)
png(joinpath(resdir, "g12b"))

# Fig.3(a)
plot(g1[2,2], 0.30, 0.55, xticks = (0.3:0.05:0.6), ylims = (0.8, 3.8), yticks = 0:0.5:4)
png(joinpath(resdir, "g22a"))

# Fig.3(b)
plot(g1[2,2], 0.4, 1.5, xticks = (0:0.2:2), ylims = (0.93, 1.15),
     yticks = 0:0.05:2)
png(joinpath(resdir, "g22b"))

#
# Tests for α values
#
ahs2 = TPT.AHS(η = 0.49, σ = [1.0, 0.3], c = [1/2, 1/2], approx=:RFA)
α2 = ahs2.α

ahs3 = TPT.AHS(η = 0.49, σ = [1.0, 0.3, 0.1], c = [1/102, 1/102, 100/102],
               approx=:RFA)
α3 = ahs3.α

@testset "Unit AHS" begin
  @test isapprox(y1, 24.0, rtol=1e-1)
  @test isapprox(y2, 4.45e+9, rtol=1e-3)
  @test isapprox(α1, 0.01886, atol=1e-5)
  @test isapprox(α2, 0.02784, atol=1e-5)
  @test isapprox(α3, 0.01837, atol=1e-5)
end


text:
export czt, iczt

"""
    czt_1d(xin , scaled , d)

Chirp z transform along a single direction d of an ND array `xin` into the ND array 'xout'.
Note that xin and xout can be the same array for inplace operations.
Note that the result type is defined by `eltype(xin)` and not by `scales`.
This code is based on a 2D Matlab version of the CZT, written by H. Gross et al.

#References: Rabiner, Schafer, Rader, The Cirp z-Transform Algorithm, IEEE Trans AU 17(1969) p. 86
"""
function czt_1d(xin, scaled, d)
    sz=size(xin)
    # returns the real datatype
    rtype = real(eltype(xin))
    scaled = rtype(scaled)
    dsize = sz[d]
    nn = (0:dsize-1)
    # 2N +1 new k-space positions
    kk = ((-dsize+1):(dsize-1))
    kk2 = rtype.((kk .^ 2) ./ 2)

    # scalar factor
    w = cispi(-2/(dsize*scaled))
    # this is needed to center the array according to the center pixel convention
    half_pix_shift = (2*(dsize÷2))/dsize
    # scalar factor. The correction factor to the right was introduced to be centered correctly on the pixel as ft is.
    a = cispi(-1/scaled * half_pix_shift)
    ww = w .^ kk2
    aa = a .^ (-nn)
    # is a 1d list of factors. This defines the shift in Fourier space (centering of frequencies)
    aa = aa .* ww[dsize .+ nn]
    to_fft = NDTools.select_region(1 ./ ww[1:(2*dsize-1)], new_size=(2*dsize,), center=(dsize+1,))
    # return tofft
    # is always 1d (small array)
    fv = fft(to_fft)

    y = xin .* reorient(aa,d)
    # twice the size along direction d
    nsz = sz .* NDTools.single_dim_size(d,2,length(sz))
    to_fft = NDTools.select_region(y, new_size=nsz, center=nsz.÷2 .+1)
    # convolve on a larger grid along one dimension
    g = ifft(fft(to_fft, d) .* reorient(fv,d), d)
    # return g
    oldctr = sz[d]÷2 + 1
    newctr = size(g) .÷ 2 .+1
    ctr = ntuple(md -> (md == d) ? newctr[md] + oldctr - 2 : newctr[md], length(newctr))

    # This is to deal with a strange phase shift appearing for odd-sized arrays
    if isodd(dsize)
        extra_phase = (2*dsize-2)/(2*dsize) # 5: 12 / 15, 7: 12/14, 9: 16/18, 11: 20/22
    else
        extra_phase = 1
    end
    # is a 1d list of factors
    fak =  ww[dsize:(2*dsize-1)] .* cispi.(ramp(rtype,1,dsize, scale=1/scaled * extra_phase))
    # return select_region(g, new_size=sz,center=ctr)
    xout = select_region(g, new_size=sz,center=ctr) .* reorient(fak,d)

    # this is a fix to deal with the problem that imaginary numbers are appearing for even-sized arrays, caused by the first entry
    if iseven(dsize) && (scaled>1.0)
        # a `scaled` of one meas this is an ordinary Fourier-transformation without zoom, which needs to keep the highes frequency value
        midp = dsize÷2+1
        for o in 1+mod(midp,2):2:dsize
            NDTools.slice(xout,d, o) .-= NDTools.slice(xin,d, 1) .* (1im).^mod(o-midp,4)
        end
    end
    return xout
    # xout .= g[dsize:(2*dsize-1)] .* reorient(fak,d)
end

"""
    czt(xin , scale, dims=1:length(size(xin)))
Chirp z transform of the ND array `xin`
This code is based on a 2D Matlab version of the CZT, written by H. Gross.
The tuple `scale` defines the zoom factors in the Fourier domain. Each has to be bigger than one.

#See also: iczt, czt_1d

#References: Rabiner, Schafer, Rader, The Cirp z-Transform Algorithm, IEEE Trans AU 17(1969) p. 86

#Example:
```jdoctest
julia> using IndexFunArrays

julia> sz = (10,10);

julia> xin = disc(sz,4)
10×10 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  1.0  1.0  1.0  1.0  0.0  0.0
 0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0
 0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0
 0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0
 0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0
 0.0  0.0  0.0  1.0  1.0  1.0  1.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0

julia> xft = czt(xin,(1.2,1.3));

julia> zoomed = real.(ift(xft))
10×10 Matrix{Float64}:
 -0.0197423    0.0233008  -0.0449251   0.00295724   0.205593  -0.166546   0.205593   0.00295724  -0.0449251   0.0233008
  0.0239759   -0.028264    0.0541186  -0.0116475   -0.261294   0.312719  -0.261294  -0.0116475    0.0541186  -0.028264
 -0.0569       0.0666104  -0.122277    0.140354     0.78259    1.34381    0.78259    0.140354    -0.122277    0.0666104
  0.00540611  -0.0117886   0.0837357   1.30651      1.8283     1.50127    1.8283     1.30651      0.0837357  -0.0117886
  0.11892     -0.147731    0.368046    1.76537      1.33218    1.66119    1.33218    1.76537      0.368046   -0.147731
 -0.00389861   0.0145979   1.21842     1.52989      1.67375    1.543      1.67375    1.52989      1.21842     0.0145979
  0.11892     -0.147731    0.368046    1.76537      1.33218    1.66119    1.33218    1.76537      0.368046   -0.147731
  0.00540611  -0.0117886   0.0837357   1.30651      1.8283     1.50127    1.8283     1.30651      0.0837357  -0.0117886
 -0.0569       0.0666104  -0.122277    0.140354     0.78259    1.34381    0.78259    0.140354    -0.122277    0.0666104
  0.0239759   -0.028264    0.0541186  -0.0116475   -0.261294   0.312719  -0.261294  -0.0116475    0.0541186  -0.028264
```
"""
function czt(xin::Array{T,N}, scale, dims=1:length(size(xin)))::Array{complex(T),N} where {T,N}
    xout = xin
    for d in dims
        xout = czt_1d(xout, scale[d], d)
    end
    return xout
end

"""
    iczt(xin , scale, dims=1:length(size(xin)))
Inverse chirp z transform of the ND array `xin`
This code is based on a 2D Matlab version of the CZT, written by H. Gross.
The tuple `scale` defines the zoom factors in the Fourier domain. Each has to be bigger than one.

#References: Rabiner, Schafer, Rader, The Cirp z-Transform Algorithm, IEEE Trans AU 17(1969) p. 86

#See also: czt, czt_1d

#Example:
```jdoctest
julia> using IndexFunArrays

julia> sz = (10,10);

julia> xin = disc(sz,4)
10×10 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  1.0  1.0  1.0  1.0  0.0  0.0
 0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0
 0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0
 0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0
 0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0
 0.0  0.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.0
 0.0  0.0  0.0  1.0  1.0  1.0  1.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0

julia> xft = ft(xin);

julia> iczt(xft,(1.2,1.3))
10×10 Matrix{ComplexF64}:
 0.00648614+0.0213779im  0.0165456+0.0357733im  0.0389356+0.0482465im  -0.235491-0.156509im    …  0.178484-0.0730099im  -0.245418-5.88331e-5im  0.0471654-0.0814548im  0.0141525+0.0734892im
  -0.104602-0.160481im   -0.163859-0.125535im    0.495205+0.135059im    0.660335+0.00736904im     0.764045-0.0497007im    0.67753+0.263814im      0.48095-0.0775406im  -0.159713-0.0637132im
   0.139304+0.111348im    0.454973+0.106869im    0.970263-0.0249785im    1.25999-0.166495im        1.07328-0.0481437im    1.24013-0.14664im      0.986722-0.0414382im   0.450186+0.111656im
  -0.035645-0.0311352im    1.03899-0.0589268im     1.1463-0.0940003im   0.790545+0.283668im       0.994255+0.134865im     0.80774-0.0124851im     1.13205+0.151519im     1.04314-0.130321im
   0.292575+0.0853233im   0.929883+0.0687029im    1.06514-0.0649952im   0.989483-0.019913im        1.02311+0.018235im    0.979555-0.136654im      1.07337+0.0317868im    0.92749+0.0405597im
    1.12254-0.0464723im    1.03467-0.0239316im    0.92709+0.0822984im     1.0521-0.0992709im   …  0.983655-0.0663123im     1.0521+0.0992709im     0.92709-0.0822984im    1.03467+0.0239316im
   0.287928-0.0306724im    0.92749-0.0405597im    1.07337-0.0317868im   0.979555+0.136654im        1.01648+0.0597475im   0.989483+0.019913im      1.06514+0.0649952im   0.929883-0.0687029im
 -0.0275957+0.169775im     1.04314+0.130321im     1.13205-0.151519im     0.80774+0.0124851im       1.00574+0.0629632im   0.790545-0.283668im       1.1463+0.0940003im    1.03899+0.0589268im
   0.130009-0.120643im    0.450186-0.111656im    0.986722+0.0414382im    1.24013+0.14664im         1.06002+0.0348813im    1.25999+0.166495im     0.970263+0.0249785im   0.454973-0.106869im
 -0.0965531+0.0404296im  -0.159713+0.0637132im    0.48095+0.0775406im    0.67753-0.263814im        0.77553-0.121603im    0.660335-0.00736904im   0.495205-0.135059im   -0.163859+0.125535im
 ```
"""
function iczt( xin , scale, dims=1:length(size(xin)))
    conj(czt(conj(xin), scale, dims)) / prod(size(xin))
end


text:
using Random, Test

function regression_test()
    rng = MersenneTwister(123);
    @testset "" begin
        nr = 4;
        nRegr = 3;
        coeffM = randn(rng, nr, nRegr) .* 10.0;
        seM = 0.1 .* coeffM;
        coeffNameV = ["beta$j"  for j = 1 : nr];
        regrNameV = ["Regr$j" for j = 1 : nRegr];
        for onePlot = [false, true]
            p = plot_regressions(coeffNameV, regrNameV, coeffM, seM;
                onePlot = onePlot);
            p = plot_regressions(Symbol.(coeffNameV), Symbol.(regrNameV),
                coeffM, seM; onePlot = onePlot);

            fName = "plot_regressions";
            if onePlot
                fName = fName * "_oneplot";
            end
            delete_test_figure(fName);
            save_test_figure(p,  fName);
            @test isfile(test_fig_path(fName))
        end
	end
end

@testset "Regression plots" begin
    regression_test();
end

# -------------


text:
using MultiConvex
using Base.Test

# write your own tests here
include("test_multiconvex.jl")


text: include("backstepping_position_controller_static_allocator_multicopter.jl")


text:
module MAFE

using DataFrames, DocStringExtensions, StatFiles

# Types
export Country, Congo, Ghana, Senegal
export Spell, SpellHistory, n_spells, get_spell, start_year, end_year, spell_years, spell_from_year_range
export MigrationHistory, Migration
export JobHistory, Job

# Stata files
export read_activity_file, read_migration_file, read_person_file, read_stata_file, clean_spells

# Making histories
export make_migration_history, make_job_history, make_job_histories, find_country

# Main
export count_wage_gains

include("constants.jl")
include("types.jl")
include("stata_files.jl")
include("main.jl")

end # module


text:
# # Model Parameters
@with_kw mutable struct CRPar
    r = 2.0
    k = 3.0
    a = 1.1
    h = 0.8
    e = 0.7
    m = 0.4
end

## CR exp resource, type I consumer
function cr_expI_model(du, u, p, t)
    @unpack r, a, e, m = p
    R, C, = u
    du[1] = r * R - a * R * C
    du[2] = e * a * R * C - m * C
    return
end

function cr_expI_model(u, p, t)
    du = similar(u)
    cr_expI_model(du, u, p, t)
    return du
end

## CR log resource, type I consumer
function cr_logI_model(du, u, p, t)
    @unpack r, k, a, e, m = p
    R, C, = u
    du[1] = r * R * (1 - R / k) - a * R * C
    du[2] = e * a * R * C - m * C
    return
end

function cr_logI_model(u, p, t)
    du = similar(u)
    cr_logI_model(du, u, p, t)
    return du
end

## CR log resource, type II consumer
function cr_logII_model(du, u, p, t)
    @unpack r, k, a, e, h, m = p
    R, C, = u
    du[1] = r * R * (1 - R / k) - a * R * C / (1 + a * h * R)
    du[2] = e * a * R * C / (1 + a * h * R) - m * C
    return
end

function cr_logII_model(u, p, t)
    du = similar(u)
    cr_logII_model(du, u, p, t)
    return du
end

## CR exp resource, type II consumer
function cr_expII_model(du, u, p, t)
    @unpack r, a, e, h, m = p
    R, C, = u
    du[1] = r * R  - a * R * C / (1 + a * h * R)
    du[2] = e * a * R * C / (1 + a * h * R) - m * C
    return
end

function cr_expII_model(u, p, t)
    du = similar(u)
    cr_expII_model(du, u, p, t)
    return du
end

#### Code that is used in multiple scripts ####

function expII_resiso(R, p)
    @unpack r, a, h = p
    return R * h * r + r / a
end

function expII_coniso(p)
    @unpack m, a, e, h = p
    return m / (a * (e - h * m))
end

# Function to calculate the jacobian at any point (with any model)
function jac(u, model, p)
    ForwardDiff.jacobian(u -> model(u, p, NaN), u)
end

# Function to calculate maximum eigenvalue (real part)
λ_stability(M) = maximum(real.(eigvals(M)))

# Function to find the interior equilibrium value for any efficiency value (for CR_logII model)
function logII_eq(p)
    @unpack r, k, a, h, m, e = p
    R = m / (a * (e - h * m))
    C = r * (-R^2 * a * h + R * a * h * k - R + k) / (a * k)
    return [R, C]
end

max_stars_repo_name: McCannLab/TheoryClub

max_issues_repo_name: McCannLab/TheoryClub

max_forks_repo_name: McCannLab/TheoryClub


text:
using LinearAlgebra
using RoME
# , IncrementalInference, TransformUtils, Distributions
using Test
# import  IncrementalInference: getSample

mutable struct RotationTest <: IncrementalInference.FunctorPairwise
  z::MvNormal
end

# 3 dimensional line, z = [a b][x y]' + c
function (rt::RotationTest)(res::Vector{Float64}, userdata, idx, meas, var1, var2)
  z = view(meas[1],:,idx)
  dq = convert(Quaternion, Euler(z...))
  s1 = so3(var1[:,idx])
  s2 = so3(var2[:,idx])
  q1 = convert(Quaternion, s1)
  q2 = convert(Quaternion, s2)
  q12 = q1*q_conj(q2)
  qq = dq*q_conj(q12)
  vee!(res, convert(so3, qq))
  nothing
end

global rr = RotationTest(MvNormal(zeros(3), 0.001*Matrix{Float64}(LinearAlgebra.I, 3,3)))

@testset "Increased dimension root finding test" begin

# known rotations

global eul = zeros(3,1)

global R1 = zeros(3,1)
global R2 = zeros(3,1)

global res = randn(3)

rr(res, nothing, 1, (zeros(3,1),), R1, R2)
@test norm(res) < 1e-10

# random rotations
global eul = 0.25*randn(3, 1)

global R1 = rand(3,1)
global R2 = rand(3,1)

global res = zeros(3)

rr(res, nothing, 1, (zeros(3),), R1, R2)

@test norm(res) > 1e-3

end

@testset "test CommonConvWrapper functions" begin

global N = 10

global i=1
for i in 1:5

global eul = 0.25*randn(3, N)

# res = zeros(3)
# @show rotationresidual!(res, eul, (zeros(0),x0))
# @show res
# gg = (res, x) -> rotationresidual!(res, eul, (zeros(0),x))
global x0 = 0.1*randn(3)
global res = zeros(3)
# @show gg(res, x0)
# @show res

global A = rand(3,N)
global B = rand(3,N)
global At = deepcopy(A)
global t = Array{Array{Float64,2},1}()
push!(t,A)
push!(t,B)
global rr = RotationTest(MvNormal(zeros(3), 0.001*Matrix{Float64}(LinearAlgebra.I, 3,3)))

global zDim = 3

global ccw = CommonConvWrapper(rr, t[1], zDim, t, measurement=(eul,)) # old bug where measurement is not patched through fixed in IIF v0.3.9

@test ccw.xDim == 3

# TODO remove
ccw.measurement = (eul,)

ccw(res, x0)

# and return complete fr/gwp
global n = 1
for n in 1:N

ccw.cpt[Threads.threadid()].particleidx = n
numericRootGenericRandomizedFnc!( ccw )

# test the result
global qq = convert(Quaternion, Euler(eul[:,ccw.cpt[Threads.threadid()].particleidx]...))
global q1 = convert(Quaternion, so3(ccw.cpt[Threads.threadid()].Y))
global q2 = convert(Quaternion, so3(B[:,ccw.cpt[Threads.threadid()].particleidx]))
@test TransformUtils.compare(q1*q_conj(q2), qq, tol=1e-8)

end # particle for

end # i for

end # testset


text:
export Car

mutable struct Car <: Robot
  btCollisionObject
end
Car() = Car(BulletCollision.convex_hull([zeros(3)]))


text:
###
### **arb_poly.h** -- polynomials over the real numbers
###

### Types, macros and constants

### Memory management
arbcall"void arb_poly_init(arb_poly_t poly)"
arbcall"void arb_poly_clear(arb_poly_t poly)"
arbcall"void arb_poly_fit_length(arb_poly_t poly, slong len)"
arbcall"void _arb_poly_set_length(arb_poly_t poly, slong len)"
arbcall"void _arb_poly_normalise(arb_poly_t poly)"
arbcall"slong arb_poly_allocated_bytes(const arb_poly_t x)"

### Basic manipulation
#mo arbcall"slong arb_poly_length(const arb_poly_t poly)" # clashes with Base.length
arbcall"slong arb_poly_degree(const arb_poly_t poly)"
arbcall"int arb_poly_is_zero(const arb_poly_t poly)"
arbcall"int arb_poly_is_one(const arb_poly_t poly)"
arbcall"int arb_poly_is_x(const arb_poly_t poly)"
arbcall"void arb_poly_zero(arb_poly_t poly)"
arbcall"void arb_poly_one(arb_poly_t poly)"
arbcall"void arb_poly_set(arb_poly_t dest, const arb_poly_t src)"
arbcall"void arb_poly_set_round(arb_poly_t dest, const arb_poly_t src, slong prec)"
arbcall"void arb_poly_set_trunc(arb_poly_t dest, const arb_poly_t src, slong n)"
arbcall"void arb_poly_set_trunc_round(arb_poly_t dest, const arb_poly_t src, slong n, slong prec)"
arbcall"void arb_poly_set_coeff_si(arb_poly_t poly, slong n, slong c)"
arbcall"void arb_poly_set_coeff_arb(arb_poly_t poly, slong n, const arb_t c)"
arbcall"void arb_poly_get_coeff_arb(arb_t v, const arb_poly_t poly, slong n)"
arbcall"void _arb_poly_shift_right(arb_ptr res, arb_srcptr poly, slong len, slong n)"
arbcall"void arb_poly_shift_right(arb_poly_t res, const arb_poly_t poly, slong n)"
arbcall"void _arb_poly_shift_left(arb_ptr res, arb_srcptr poly, slong len, slong n)"
arbcall"void arb_poly_shift_left(arb_poly_t res, const arb_poly_t poly, slong n)"
arbcall"void arb_poly_truncate(arb_poly_t poly, slong n)"
arbcall"slong arb_poly_valuation(const arb_poly_t poly)"

### Conversions
#ni arbcall"void arb_poly_set_fmpz_poly(arb_poly_t poly, const fmpz_poly_t src, slong prec)"
#ni arbcall"void arb_poly_set_fmpq_poly(arb_poly_t poly, const fmpq_poly_t src, slong prec)"
arbcall"void arb_poly_set_si(arb_poly_t poly, slong src)"

### Input and output
arbcall"void arb_poly_printd(const arb_poly_t poly, slong digits)"
#ns arbcall"void arb_poly_fprintd(FILE * file, const arb_poly_t poly, slong digits)"

### Random generation
#ns arbcall"void arb_poly_randtest(arb_poly_t poly, flint_rand_t state, slong len, slong prec, slong mag_bits)"

### Comparisons
arbcall"int arb_poly_contains(const arb_poly_t poly1, const arb_poly_t poly2)"
#ni arbcall"int arb_poly_contains_fmpz_poly(const arb_poly_t poly1, const fmpz_poly_t poly2)"
#ni arbcall"int arb_poly_contains_fmpq_poly(const arb_poly_t poly1, const fmpq_poly_t poly2)"
arbcall"int arb_poly_equal(const arb_poly_t A, const arb_poly_t B)"
arbcall"int _arb_poly_overlaps(arb_srcptr poly1, slong len1, arb_srcptr poly2, slong len2)"
arbcall"int arb_poly_overlaps(const arb_poly_t poly1, const arb_poly_t poly2)"
#ni arbcall"int arb_poly_get_unique_fmpz_poly(fmpz_poly_t z, const arb_poly_t x)"

### Bounds
arbcall"void _arb_poly_majorant(arb_ptr res, arb_srcptr poly, slong len, slong prec)"
arbcall"void arb_poly_majorant(arb_poly_t res, const arb_poly_t poly, slong prec)"

### Arithmetic
#mo arbcall"void _arb_poly_add(arb_ptr C, arb_srcptr A, slong lenA, arb_srcptr B, slong lenB, slong prec)" # clashes with _arb_vec_add
arbcall"void arb_poly_add(arb_poly_t C, const arb_poly_t A, const arb_poly_t B, slong prec)"
arbcall"void arb_poly_add_si(arb_poly_t C, const arb_poly_t A, slong B, slong prec)"
#mo arbcall"void _arb_poly_sub(arb_ptr C, arb_srcptr A, slong lenA, arb_srcptr B, slong lenB, slong prec)" # clashes with _arb_vec_sub
arbcall"void arb_poly_sub(arb_poly_t C, const arb_poly_t A, const arb_poly_t B, slong prec)"
arbcall"void arb_poly_add_series(arb_poly_t C, const arb_poly_t A, const arb_poly_t B, slong len, slong prec)"
arbcall"void arb_poly_sub_series(arb_poly_t C, const arb_poly_t A, const arb_poly_t B, slong len, slong prec)"
arbcall"void arb_poly_neg(arb_poly_t C, const arb_poly_t A)"
arbcall"void arb_poly_scalar_mul_2exp_si(arb_poly_t C, const arb_poly_t A, slong c)"
arbcall"void arb_poly_scalar_mul(arb_poly_t C, const arb_poly_t A, const arb_t c, slong prec)"
arbcall"void arb_poly_scalar_div(arb_poly_t C, const arb_poly_t A, const arb_t c, slong prec)"
arbcall"void _arb_poly_mullow_classical(arb_ptr C, arb_srcptr A, slong lenA, arb_srcptr B, slong lenB, slong n, slong prec)"
arbcall"void _arb_poly_mullow_block(arb_ptr C, arb_srcptr A, slong lenA, arb_srcptr B, slong lenB, slong n, slong prec)"
arbcall"void _arb_poly_mullow(arb_ptr C, arb_srcptr A, slong lenA, arb_srcptr B, slong lenB, slong n, slong prec)"
arbcall"void arb_poly_mullow_classical(arb_poly_t C, const arb_poly_t A, const arb_poly_t B, slong n, slong prec)"
arbcall"void arb_poly_mullow_ztrunc(arb_poly_t C, const arb_poly_t A, const arb_poly_t B, slong n, slong prec)"
arbcall"void arb_poly_mullow_block(arb_poly_t C, const arb_poly_t A, const arb_poly_t B, slong n, slong prec)"
arbcall"void arb_poly_mullow(arb_poly_t C, const arb_poly_t A, const arb_poly_t B, slong n, slong prec)"
arbcall"void _arb_poly_mul(arb_ptr C, arb_srcptr A, slong lenA, arb_srcptr B, slong lenB, slong prec)"
arbcall"void arb_poly_mul(arb_poly_t C, const arb_poly_t A, const arb_poly_t B, slong prec)"
arbcall"void _arb_poly_inv_series(arb_ptr Q, arb_srcptr A, slong Alen, slong len, slong prec)"
arbcall"void arb_poly_inv_series(arb_poly_t Q, const arb_poly_t A, slong n, slong prec)"
arbcall"void _arb_poly_div_series(arb_ptr Q, arb_srcptr A, slong Alen, arb_srcptr B, slong Blen, slong n, slong prec)"
arbcall"void arb_poly_div_series(arb_poly_t Q, const arb_poly_t A, const arb_poly_t B, slong n, slong prec)"
arbcall"void _arb_poly_div(arb_ptr Q, arb_srcptr A, slong lenA, arb_srcptr B, slong lenB, slong prec)"
arbcall"void _arb_poly_rem(arb_ptr R, arb_srcptr A, slong lenA, arb_srcptr B, slong lenB, slong prec)"
arbcall"void _arb_poly_divrem(arb_ptr Q, arb_ptr R, arb_srcptr A, slong lenA, arb_srcptr B, slong lenB, slong prec)"
arbcall"int arb_poly_divrem(arb_poly_t Q, arb_poly_t R, const arb_poly_t A, const arb_poly_t B, slong prec)"
arbcall"void _arb_poly_div_root(arb_ptr Q, arb_t R, arb_srcptr A, slong len, const arb_t c, slong prec)"

### Composition
arbcall"void _arb_poly_taylor_shift_horner(arb_ptr g, const arb_t c, slong n, slong prec)"
arbcall"void arb_poly_taylor_shift_horner(arb_poly_t g, const arb_poly_t f, const arb_t c, slong prec)"
arbcall"void _arb_poly_taylor_shift_divconquer(arb_ptr g, const arb_t c, slong n, slong prec)"
arbcall"void arb_poly_taylor_shift_divconquer(arb_poly_t g, const arb_poly_t f, const arb_t c, slong prec)"
arbcall"void _arb_poly_taylor_shift_convolution(arb_ptr g, const arb_t c, slong n, slong prec)"
arbcall"void arb_poly_taylor_shift_convolution(arb_poly_t g, const arb_poly_t f, const arb_t c, slong prec)"
arbcall"void _arb_poly_taylor_shift(arb_ptr g, const arb_t c, slong n, slong prec)"
arbcall"void arb_poly_taylor_shift(arb_poly_t g, const arb_poly_t f, const arb_t c, slong prec)"
arbcall"void _arb_poly_compose_horner(arb_ptr res, arb_srcptr poly1, slong len1, arb_srcptr poly2, slong len2, slong prec)"
arbcall"void arb_poly_compose_horner(arb_poly_t res, const arb_poly_t poly1, const arb_poly_t poly2, slong prec)"
arbcall"void _arb_poly_compose_divconquer(arb_ptr res, arb_srcptr poly1, slong len1, arb_srcptr poly2, slong len2, slong prec)"
arbcall"void arb_poly_compose_divconquer(arb_poly_t res, const arb_poly_t poly1, const arb_poly_t poly2, slong prec)"
arbcall"void _arb_poly_compose(arb_ptr res, arb_srcptr poly1, slong len1, arb_srcptr poly2, slong len2, slong prec)"
arbcall"void arb_poly_compose(arb_poly_t res, const arb_poly_t poly1, const arb_poly_t poly2, slong prec)"
arbcall"void _arb_poly_compose_series_horner(arb_ptr res, arb_srcptr poly1, slong len1, arb_srcptr poly2, slong len2, slong n, slong prec)"
arbcall"void arb_poly_compose_series_horner(arb_poly_t res, const arb_poly_t poly1, const arb_poly_t poly2, slong n, slong prec)"
arbcall"void _arb_poly_compose_series_brent_kung(arb_ptr res, arb_srcptr poly1, slong len1, arb_srcptr poly2, slong len2, slong n, slong prec)"
arbcall"void arb_poly_compose_series_brent_kung(arb_poly_t res, const arb_poly_t poly1, const arb_poly_t poly2, slong n, slong prec)"
arbcall"void _arb_poly_compose_series(arb_ptr res, arb_srcptr poly1, slong len1, arb_srcptr poly2, slong len2, slong n, slong prec)"
arbcall"void arb_poly_compose_series(arb_poly_t res, const arb_poly_t poly1, const arb_poly_t poly2, slong n, slong prec)"
arbcall"void _arb_poly_revert_series_lagrange(arb_ptr h, arb_srcptr f, slong flen, slong n, slong prec)"
arbcall"void arb_poly_revert_series_lagrange(arb_poly_t h, const arb_poly_t f, slong n, slong prec)"
arbcall"void _arb_poly_revert_series_newton(arb_ptr h, arb_srcptr f, slong flen, slong n, slong prec)"
arbcall"void arb_poly_revert_series_newton(arb_poly_t h, const arb_poly_t f, slong n, slong prec)"
arbcall"void _arb_poly_revert_series_lagrange_fast(arb_ptr h, arb_srcptr f, slong flen, slong n, slong prec)"
arbcall"void arb_poly_revert_series_lagrange_fast(arb_poly_t h, const arb_poly_t f, slong n, slong prec)"
arbcall"void _arb_poly_revert_series(arb_ptr h, arb_srcptr f, slong flen, slong n, slong prec)"
arbcall"void arb_poly_revert_series(arb_poly_t h, const arb_poly_t f, slong n, slong prec)"

### Evaluation
arbcall"void _arb_poly_evaluate_horner(arb_t y, arb_srcptr f, slong len, const arb_t x, slong prec)"
arbcall"void arb_poly_evaluate_horner(arb_t y, const arb_poly_t f, const arb_t x, slong prec)"
arbcall"void _arb_poly_evaluate_rectangular(arb_t y, arb_srcptr f, slong len, const arb_t x, slong prec)"
arbcall"void arb_poly_evaluate_rectangular(arb_t y, const arb_poly_t f, const arb_t x, slong prec)"
arbcall"void _arb_poly_evaluate(arb_t y, arb_srcptr f, slong len, const arb_t x, slong prec)"
arbcall"void arb_poly_evaluate(arb_t y, const arb_poly_t f, const arb_t x, slong prec)"
arbcall"void _arb_poly_evaluate_acb_horner(acb_t y, arb_srcptr f, slong len, const acb_t x, slong prec)"
arbcall"void arb_poly_evaluate_acb_horner(acb_t y, const arb_poly_t f, const acb_t x, slong prec)"
arbcall"void _arb_poly_evaluate_acb_rectangular(acb_t y, arb_srcptr f, slong len, const acb_t x, slong prec)"
arbcall"void arb_poly_evaluate_acb_rectangular(acb_t y, const arb_poly_t f, const acb_t x, slong prec)"
arbcall"void _arb_poly_evaluate_acb(acb_t y, arb_srcptr f, slong len, const acb_t x, slong prec)"
arbcall"void arb_poly_evaluate_acb(acb_t y, const arb_poly_t f, const acb_t x, slong prec)"
arbcall"void _arb_poly_evaluate2_horner(arb_t y, arb_t z, arb_srcptr f, slong len, const arb_t x, slong prec)"
arbcall"void arb_poly_evaluate2_horner(arb_t y, arb_t z, const arb_poly_t f, const arb_t x, slong prec)"
arbcall"void _arb_poly_evaluate2_rectangular(arb_t y, arb_t z, arb_srcptr f, slong len, const arb_t x, slong prec)"
arbcall"void arb_poly_evaluate2_rectangular(arb_t y, arb_t z, const arb_poly_t f, const arb_t x, slong prec)"
arbcall"void _arb_poly_evaluate2(arb_t y, arb_t z, arb_srcptr f, slong len, const arb_t x, slong prec)"
arbcall"void arb_poly_evaluate2(arb_t y, arb_t z, const arb_poly_t f, const arb_t x, slong prec)"
arbcall"void _arb_poly_evaluate2_acb_horner(acb_t y, acb_t z, arb_srcptr f, slong len, const acb_t x, slong prec)"
arbcall"void arb_poly_evaluate2_acb_horner(acb_t y, acb_t z, const arb_poly_t f, const acb_t x, slong prec)"
arbcall"void _arb_poly_evaluate2_acb_rectangular(acb_t y, acb_t z, arb_srcptr f, slong len, const acb_t x, slong prec)"
arbcall"void arb_poly_evaluate2_acb_rectangular(acb_t y, acb_t z, const arb_poly_t f, const acb_t x, slong prec)"
arbcall"void _arb_poly_evaluate2_acb(acb_t y, acb_t z, arb_srcptr f, slong len, const acb_t x, slong prec)"
arbcall"void arb_poly_evaluate2_acb(acb_t y, acb_t z, const arb_poly_t f, const acb_t x, slong prec)"

### Product trees
arbcall"void _arb_poly_product_roots(arb_ptr poly, arb_srcptr xs, slong n, slong prec)"
arbcall"void arb_poly_product_roots(arb_poly_t poly, arb_srcptr xs, slong n, slong prec)"
arbcall"void _arb_poly_product_roots_complex(arb_ptr poly, arb_srcptr r, slong rn, acb_srcptr c, slong cn, slong prec)"
arbcall"void arb_poly_product_roots_complex(arb_poly_t poly, arb_srcptr r, slong rn, acb_srcptr c, slong cn, slong prec)"
#ni arbcall"arb_ptr * _arb_poly_tree_alloc(slong len)"
#ni arbcall"void _arb_poly_tree_free(arb_ptr * tree, slong len)"
#ni arbcall"void _arb_poly_tree_build(arb_ptr * tree, arb_srcptr roots, slong len, slong prec)"

### Multipoint evaluation
arbcall"void _arb_poly_evaluate_vec_iter(arb_ptr ys, arb_srcptr poly, slong plen, arb_srcptr xs, slong n, slong prec)"
arbcall"void arb_poly_evaluate_vec_iter(arb_ptr ys, const arb_poly_t poly, arb_srcptr xs, slong n, slong prec)"
#ni arbcall"void _arb_poly_evaluate_vec_fast_precomp(arb_ptr vs, arb_srcptr poly, slong plen, arb_ptr * tree, slong len, slong prec)"
arbcall"void _arb_poly_evaluate_vec_fast(arb_ptr ys, arb_srcptr poly, slong plen, arb_srcptr xs, slong n, slong prec)"
arbcall"void arb_poly_evaluate_vec_fast(arb_ptr ys, const arb_poly_t poly, arb_srcptr xs, slong n, slong prec)"

### Interpolation
arbcall"void _arb_poly_interpolate_newton(arb_ptr poly, arb_srcptr xs, arb_srcptr ys, slong n, slong prec)"
arbcall"void arb_poly_interpolate_newton(arb_poly_t poly, arb_srcptr xs, arb_srcptr ys, slong n, slong prec)"
arbcall"void _arb_poly_interpolate_barycentric(arb_ptr poly, arb_srcptr xs, arb_srcptr ys, slong n, slong prec)"
arbcall"void arb_poly_interpolate_barycentric(arb_poly_t poly, arb_srcptr xs, arb_srcptr ys, slong n, slong prec)"
#ni arbcall"void _arb_poly_interpolation_weights(arb_ptr w, arb_ptr * tree, slong len, slong prec)"
#ni arbcall"void _arb_poly_interpolate_fast_precomp(arb_ptr poly, arb_srcptr ys, arb_ptr * tree, arb_srcptr weights, slong len, slong prec)"
arbcall"void _arb_poly_interpolate_fast(arb_ptr poly, arb_srcptr xs, arb_srcptr ys, slong len, slong prec)"
arbcall"void arb_poly_interpolate_fast(arb_poly_t poly, arb_srcptr xs, arb_srcptr ys, slong n, slong prec)"

### Differentiation
arbcall"void _arb_poly_derivative(arb_ptr res, arb_srcptr poly, slong len, slong prec)"
arbcall"void arb_poly_derivative(arb_poly_t res, const arb_poly_t poly, slong prec)"
arbcall"void _arb_poly_integral(arb_ptr res, arb_srcptr poly, slong len, slong prec)"
arbcall"void arb_poly_integral(arb_poly_t res, const arb_poly_t poly, slong prec)"

### Transforms
arbcall"void _arb_poly_borel_transform(arb_ptr res, arb_srcptr poly, slong len, slong prec)"
arbcall"void arb_poly_borel_transform(arb_poly_t res, const arb_poly_t poly, slong prec)"
arbcall"void _arb_poly_inv_borel_transform(arb_ptr res, arb_srcptr poly, slong len, slong prec)"
arbcall"void arb_poly_inv_borel_transform(arb_poly_t res, const arb_poly_t poly, slong prec)"
arbcall"void _arb_poly_binomial_transform_basecase(arb_ptr b, arb_srcptr a, slong alen, slong len, slong prec)"
arbcall"void arb_poly_binomial_transform_basecase(arb_poly_t b, const arb_poly_t a, slong len, slong prec)"
arbcall"void _arb_poly_binomial_transform_convolution(arb_ptr b, arb_srcptr a, slong alen, slong len, slong prec)"
arbcall"void arb_poly_binomial_transform_convolution(arb_poly_t b, const arb_poly_t a, slong len, slong prec)"
arbcall"void _arb_poly_binomial_transform(arb_ptr b, arb_srcptr a, slong alen, slong len, slong prec)"
arbcall"void arb_poly_binomial_transform(arb_poly_t b, const arb_poly_t a, slong len, slong prec)"

### Powers and elementary functions
arbcall"void _arb_poly_pow_ui_trunc_binexp(arb_ptr res, arb_srcptr f, slong flen, ulong exp, slong len, slong prec)"
arbcall"void arb_poly_pow_ui_trunc_binexp(arb_poly_t res, const arb_poly_t poly, ulong exp, slong len, slong prec)"
arbcall"void _arb_poly_pow_ui(arb_ptr res, arb_srcptr f, slong flen, ulong exp, slong prec)"
arbcall"void arb_poly_pow_ui(arb_poly_t res, const arb_poly_t poly, ulong exp, slong prec)"
arbcall"void _arb_poly_pow_series(arb_ptr h, arb_srcptr f, slong flen, arb_srcptr g, slong glen, slong len, slong prec)"
arbcall"void arb_poly_pow_series(arb_poly_t h, const arb_poly_t f, const arb_poly_t g, slong len, slong prec)"
arbcall"void _arb_poly_pow_arb_series(arb_ptr h, arb_srcptr f, slong flen, const arb_t g, slong len, slong prec)"
arbcall"void arb_poly_pow_arb_series(arb_poly_t h, const arb_poly_t f, const arb_t g, slong len, slong prec)"
arbcall"void _arb_poly_sqrt_series(arb_ptr g, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_sqrt_series(arb_poly_t g, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_rsqrt_series(arb_ptr g, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_rsqrt_series(arb_poly_t g, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_log_series(arb_ptr res, arb_srcptr f, slong flen, slong n, slong prec)"
arbcall"void arb_poly_log_series(arb_poly_t res, const arb_poly_t f, slong n, slong prec)"
arbcall"void _arb_poly_log1p_series(arb_ptr res, arb_srcptr f, slong flen, slong n, slong prec)"
arbcall"void arb_poly_log1p_series(arb_poly_t res, const arb_poly_t f, slong n, slong prec)"
arbcall"void _arb_poly_atan_series(arb_ptr res, arb_srcptr f, slong flen, slong n, slong prec)"
arbcall"void arb_poly_atan_series(arb_poly_t res, const arb_poly_t f, slong n, slong prec)"
arbcall"void _arb_poly_asin_series(arb_ptr res, arb_srcptr f, slong flen, slong n, slong prec)"
arbcall"void arb_poly_asin_series(arb_poly_t res, const arb_poly_t f, slong n, slong prec)"
arbcall"void _arb_poly_acos_series(arb_ptr res, arb_srcptr f, slong flen, slong n, slong prec)"
arbcall"void arb_poly_acos_series(arb_poly_t res, const arb_poly_t f, slong n, slong prec)"
arbcall"void _arb_poly_exp_series_basecase(arb_ptr f, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_exp_series_basecase(arb_poly_t f, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_exp_series(arb_ptr f, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_exp_series(arb_poly_t f, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_sin_cos_series_basecase(arb_ptr s, arb_ptr c, arb_srcptr h, slong hlen, slong n, slong prec, int times_pi)"
arbcall"void arb_poly_sin_cos_series_basecase(arb_poly_t s, arb_poly_t c, const arb_poly_t h, slong n, slong prec, int times_pi)"
arbcall"void _arb_poly_sin_cos_series_tangent(arb_ptr s, arb_ptr c, arb_srcptr h, slong hlen, slong n, slong prec, int times_pi)"
arbcall"void arb_poly_sin_cos_series_tangent(arb_poly_t s, arb_poly_t c, const arb_poly_t h, slong n, slong prec, int times_pi)"
arbcall"void _arb_poly_sin_cos_series(arb_ptr s, arb_ptr c, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_sin_cos_series(arb_poly_t s, arb_poly_t c, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_sin_series(arb_ptr s, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_sin_series(arb_poly_t s, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_cos_series(arb_ptr c, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_cos_series(arb_poly_t c, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_tan_series(arb_ptr g, arb_srcptr h, slong hlen, slong len, slong prec)"
arbcall"void arb_poly_tan_series(arb_poly_t g, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_sin_cos_pi_series(arb_ptr s, arb_ptr c, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_sin_cos_pi_series(arb_poly_t s, arb_poly_t c, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_sin_pi_series(arb_ptr s, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_sin_pi_series(arb_poly_t s, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_cos_pi_series(arb_ptr c, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_cos_pi_series(arb_poly_t c, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_cot_pi_series(arb_ptr c, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_cot_pi_series(arb_poly_t c, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_sinh_cosh_series_basecase(arb_ptr s, arb_ptr c, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_sinh_cosh_series_basecase(arb_poly_t s, arb_poly_t c, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_sinh_cosh_series_exponential(arb_ptr s, arb_ptr c, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_sinh_cosh_series_exponential(arb_poly_t s, arb_poly_t c, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_sinh_cosh_series(arb_ptr s, arb_ptr c, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_sinh_cosh_series(arb_poly_t s, arb_poly_t c, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_sinh_series(arb_ptr s, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_sinh_series(arb_poly_t s, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_cosh_series(arb_ptr c, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_cosh_series(arb_poly_t c, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_sinc_series(arb_ptr s, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_sinc_series(arb_poly_t s, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_sinc_pi_series(arb_ptr s, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_sinc_pi_series(arb_poly_t s, const arb_poly_t h, slong n, slong prec)"

### Lambert W function
arbcall"void _arb_poly_lambertw_series(arb_ptr res, arb_srcptr z, slong zlen, int flags, slong len, slong prec)"
arbcall"void arb_poly_lambertw_series(arb_poly_t res, const arb_poly_t z, int flags, slong len, slong prec)"

### Gamma function and factorials
arbcall"void _arb_poly_gamma_series(arb_ptr res, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_gamma_series(arb_poly_t res, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_rgamma_series(arb_ptr res, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_rgamma_series(arb_poly_t res, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_lgamma_series(arb_ptr res, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_lgamma_series(arb_poly_t res, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_digamma_series(arb_ptr res, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_digamma_series(arb_poly_t res, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_rising_ui_series(arb_ptr res, arb_srcptr f, slong flen, ulong r, slong trunc, slong prec)"
arbcall"void arb_poly_rising_ui_series(arb_poly_t res, const arb_poly_t f, ulong r, slong trunc, slong prec)"

### Zeta function
arbcall"void arb_poly_zeta_series(arb_poly_t res, const arb_poly_t s, const arb_t a, int deflate, slong n, slong prec)"
arbcall"void _arb_poly_riemann_siegel_theta_series(arb_ptr res, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_riemann_siegel_theta_series(arb_poly_t res, const arb_poly_t h, slong n, slong prec)"
arbcall"void _arb_poly_riemann_siegel_z_series(arb_ptr res, arb_srcptr h, slong hlen, slong n, slong prec)"
arbcall"void arb_poly_riemann_siegel_z_series(arb_poly_t res, const arb_poly_t h, slong n, slong prec)"

### Root-finding
arbcall"void _arb_poly_root_bound_fujiwara(mag_t bound, arb_srcptr poly, slong len)"
arbcall"void arb_poly_root_bound_fujiwara(mag_t bound, arb_poly_t poly)"
arbcall"void _arb_poly_newton_convergence_factor(arf_t convergence_factor, arb_srcptr poly, slong len, const arb_t convergence_interval, slong prec)"
arbcall"int _arb_poly_newton_step(arb_t xnew, arb_srcptr poly, slong len, const arb_t x, const arb_t convergence_interval, const arf_t convergence_factor, slong prec)"
arbcall"void _arb_poly_newton_refine_root(arb_t r, arb_srcptr poly, slong len, const arb_t start, const arb_t convergence_interval, const arf_t convergence_factor, slong eval_extra_prec, slong prec)"

### Other special polynomials
arbcall"void _arb_poly_swinnerton_dyer_ui(arb_ptr poly, ulong n, slong trunc, slong prec)"
arbcall"void arb_poly_swinnerton_dyer_ui(arb_poly_t poly, ulong n, slong prec)"