text:
include("../src/read/read_datasets.jl")
include("../src/utils.jl")
include("../src/pq/PQ.jl");
include("../src/linscan/Linscan.jl");

function demo_pq(
  dataset_name="SIFT1M",
  nread::Integer=Int(1e4)) # Increase this to 1e5 to use the full dataset

  # === Hyperparams ===
  m       = 8
  h       = 256
  verbose = true
  nquery  = Int(1e4)
  knn     = Int(1e3) # Compute recall up to
  b       = Int( log2(h) * m )

  # === Train ===
  x_train           = read_dataset(dataset_name, nread )
  C, B, train_error = train_pq(x_train, m, h, verbose)
  @printf("Error in training is %e\n", train_error)

  # === Encode the base set ===
  nread_base = Int(1e6)
  x_base     = read_dataset(dataset_name * "_base", nread_base )
  B_base     = quantize_pq( x_base, C, verbose )
  base_error = qerror_pq( x_base, B_base, C )
  @printf("Error in base is %e\n", base_error)

  # === Compute recall ===
  x_query = read_dataset( dataset_name * "_query", nquery, verbose )
  gt      = read_dataset( dataset_name * "_groundtruth", nquery, verbose )
  if dataset_name == "SIFT1M" || dataset_name == "GIST1M"
    gt = gt + 1;
  end
  gt = convert( Vector{UInt32}, gt[1,1:nquery] )
  B_base  = convert( Matrix{UInt8}, B_base-1 )

  print("Querying m=$m ... ");
  @time dists, idx = linscan_pq( B_base, x_query[:,1:nquery], C, b, knn )
  println("done");

  rec = eval_recall( gt, idx, knn );

end # function train_process

# train
demo_pq()

max_issues_repo_name: xiaosu-zhu/local-search-quantization

max_forks_repo_name: xiaosu-zhu/local-search-quantization


text:
using HTTP

function fetch_cap(url)
    servicename = "git-upload-pack"
    r1 = HTTP.get("$url/info/refs?service=$servicename",["Git-Protocol" => "version=1"]; verbose=2)
    @assert r1.status in (200, 304)

    h1 = Dict(r1.headers)
    @assert h1["Content-Type"] == "application/x-$servicename-advertisement"

    b1 = IOBuffer(r1.body)
    pkt1 = readpktline(b1)
    @assert pkt1 == "# service=$servicename"

    pkt2 = readpktline(b1)
    @assert pkt2 === nothing

    fetch_refs(b1)
end

ext: jl

max_stars_repo_name: simonbyrne/GitX

max_issues_repo_name: simonbyrne/GitX

max_forks_repo_name: simonbyrne/GitX


text:
module TestOneDimensionalRangeIterators

using Test
using MLJBase
using Random
Random.seed!(123)

mutable struct DummyModel <: Deterministic
    K::Int
    metric::Float64
    kernel::Char
end

dummy_model = DummyModel(4, 9.5, 'k')

mutable struct SuperModel <: Deterministic
    lambda::Float64
    model1::DummyModel
    model2::DummyModel
end

dummy1 = DummyModel(1, 9.5, 'k')
dummy2 = DummyModel(2, 9.5, 'k')
super_model = SuperModel(0.5, dummy1, dummy2)

z1 = range(dummy_model, :K, lower=1, upper=10)
z2 = range(dummy_model, :K, lower=10, origin=10^6, unit=10^5)
z3 = range(dummy_model, :K, upper=-10, origin=-10^6, unit=10^5)
z4 = range(super_model, :lambda, lower=1, upper=10)
z5 = range(dummy_model, :K, origin=10, unit=20)
p1 = range(dummy_model, :K, lower=1, upper=10, scale=:log10)
p2 = range(dummy_model, :kernel, values=['c', 'd'])
p3 = range(super_model, :lambda, lower=0.1, upper=1, scale=:log2)
p4 = range(dummy_model, :K, lower=1, upper=3, scale=x->2x)

@testset "scale transformations" begin
    @test transform(MLJBase.Scale, scale(:log), ℯ) == 1
    @test inverse_transform(MLJBase.Scale, scale(:log), 1) == float(ℯ)
end

@testset "inspecting scales of ranges" begin
    @test scale(p1) == :log10
    @test scale(p2) == :none
    @test scale(p3) == :log2
    @test scale(p4) == :custom
    @test scale(sin) === sin
end

@testset "iterators" begin

    @test iterator(p1, 5)  == [1, 2, 3, 6, 10]
    @test iterator(p2) == collect(p2.values)
    u = 2^(log2(0.1)/2)
    @test iterator(p3, 3) ≈ [0.1, u, 1]
    @test iterator(p4, 3) == [2, 4, 6]

    # semi-unbounded ranges:
    v = Int.(round.(exp.([(1-t)*log(10) + t*log(10+2e5)
                     for t in 0:(1/3):1]))) |> unique
    @test iterator(z2, 4) == v
    @test iterator(z3, 4) == reverse(-v)

    # doubly-unbounded ranges:
    @test iterator(z5, 4) ==
        iterator(range(Int, :foo, lower=-10, upper=30), 4)

    # truncated nominal range iteration:
    rr = range(Char, :foo, values=collect("abcdefg"))
    @test iterator(rr, 3) == ['a', 'b', 'c']

    # random:
    rng = MersenneTwister(123)
    @test iterator(p1, 5, rng) == [10, 6, 1, 2, 3]
    @test iterator(rr, rng) == ['b', 'a', 'g', 'c', 'e', 'f', 'd']
    @test iterator(rr, 3, rng) == ['b', 'c', 'a']

end

end
true


text:
module SymbolicRegression

# Types
export Population,
    PopMember,
    HallOfFame,
    Options,
    Node,

    #Functions:
    EquationSearch,
    s_r_cycle,
    calculate_pareto_frontier,
    count_nodes,
    compute_complexity,
    copy_node,
    print_tree,
    string_tree,
    eval_tree_array,
    eval_diff_tree_array,
    eval_grad_tree_array,
    differentiable_eval_tree_array,
    node_to_symbolic,
    symbolic_to_node,
    combine_operators,
    gen_random_tree,
    gen_random_tree_fixed_size,

    #Operators
    plus,
    sub,
    mult,
    square,
    cube,
    pow,
    div,
    log_abs,
    log2_abs,
    log10_abs,
    log1p_abs,
    acosh_abs,
    sqrt_abs,
    neg,
    greater,
    relu,
    logical_or,
    logical_and,

    # special operators
    gamma,
    erf,
    erfc,
    atanh_clip

using Distributed
using JSON3: JSON3
import Printf: @printf, @sprintf
using Pkg: Pkg
import Random: seed!, shuffle!
using Reexport

# for symbolic constraints
using PyCall

@reexport import LossFunctions:
    MarginLoss,
    DistanceLoss,
    SupervisedLoss,
    ZeroOneLoss,
    LogitMarginLoss,
    PerceptronLoss,
    HingeLoss,
    L1HingeLoss,
    L2HingeLoss,
    SmoothedL1HingeLoss,
    ModifiedHuberLoss,
    L2MarginLoss,
    ExpLoss,
    SigmoidLoss,
    DWDMarginLoss,
    LPDistLoss,
    L1DistLoss,
    L2DistLoss,
    PeriodicLoss,
    HuberLoss,
    EpsilonInsLoss,
    L1EpsilonInsLoss,
    L2EpsilonInsLoss,
    LogitDistLoss,
    QuantileLoss,
    LogCoshLoss

include("Core.jl")
include("Recorder.jl")
include("Utils.jl")
include("EquationUtils.jl")
include("EvaluateEquation.jl")
include("EvaluateEquationDerivative.jl")
include("CheckConstraints.jl")
include("MutationFunctions.jl")
include("LossFunctions.jl")
include("PopMember.jl")
include("ConstantOptimization.jl")
include("Population.jl")
include("HallOfFame.jl")
include("InterfaceSymbolicUtils.jl")
include("SimplifyEquation.jl")
include("Mutate.jl")
include("RegularizedEvolution.jl")
include("SingleIteration.jl")
include("ProgressBars.jl")

import .CoreModule:
    CONST_TYPE,
    MAX_DEGREE,
    BATCH_DIM,
    FEATURE_DIM,
    RecordType,
    Dataset,
    Node,
    copy_node,
    Options,
    plus,
    sub,
    mult,
    square,
    cube,
    pow,
    div,
    log_abs,
    log2_abs,
    log10_abs,
    log1p_abs,
    sqrt_abs,
    acosh_abs,
    neg,
    greater,
    greater,
    relu,
    logical_or,
    logical_and,
    gamma,
    erf,
    erfc,
    atanh_clip,
    SRConcurrency,
    SRSerial,
    SRThreaded,
    SRDistributed,
    string_tree,
    print_tree
import .UtilsModule:
    debug, debug_inline, is_anonymous_function, recursive_merge, next_worker, @sr_spawner
import .EquationUtilsModule:
    count_nodes,
    compute_complexity,
    get_constants,
    set_constants,
    index_constants,
    NodeIndex
import .EvaluateEquationModule: eval_tree_array, differentiable_eval_tree_array
import .EvaluateEquationDerivativeModule: eval_diff_tree_array, eval_grad_tree_array
import .CheckConstraintsModule: check_constraints
import .MutationFunctionsModule:
    gen_random_tree,
    gen_random_tree_fixed_size,
    random_node,
    random_node_and_parent,
    crossover_trees
import .LossFunctionsModule: eval_loss, loss, score_func
import .PopMemberModule: PopMember, copy_pop_member
import .PopulationModule: Population, best_sub_pop, record_population, best_of_sample
import .HallOfFameModule:
    HallOfFame, calculate_pareto_frontier, string_dominating_pareto_curve
import .SingleIterationModule: s_r_cycle, optimize_and_simplify_population
import .InterfaceSymbolicUtilsModule: node_to_symbolic, symbolic_to_node
import .SimplifyEquationModule: combine_operators, simplify_tree
import .ProgressBarsModule: ProgressBar, set_multiline_postfix
import .RecorderModule: @recorder, find_iteration_from_record

include("Configure.jl")
include("Deprecates.jl")

StateType{T} = Tuple{
    Union{Vector{Vector{Population{T}}},Matrix{Population{T}}},
    Union{HallOfFame,Vector{HallOfFame}},
}

"""
    EquationSearch(X, y[; kws...])

Perform a distributed equation search for functions `f_i` which
describe the mapping `f_i(X[:, j]) ≈ y[i, j]`. Options are
configured using SymbolicRegression.Options(...),
which should be passed as a keyword argument to options.
One can turn off parallelism with `numprocs=0`,
which is useful for debugging and profiling.

# Arguments
- `X::AbstractMatrix{T}`:  The input dataset to predict `y` from.
    The first dimension is features, the second dimension is rows.
- `y::Union{AbstractMatrix{T}, AbstractVector{T}}`: The values to predict. The first dimension
    is the output feature to predict with each equation, and the
    second dimension is rows.
- `niterations::Int=10`: The number of iterations to perform the search.
    More iterations will improve the results.
- `weights::Union{AbstractMatrix{T}, AbstractVector{T}, Nothing}=nothing`: Optionally
    weight the loss for each `y` by this value (same shape as `y`).
- `varMap::Union{Array{String, 1}, Nothing}=nothing`: The names
    of each feature in `X`, which will be used during printing of equations.
- `options::Options=Options()`: The options for the search, such as
    which operators to use, evolution hyperparameters, etc.
- `numprocs::Union{Int, Nothing}=nothing`:  The number of processes to use,
    if you want `EquationSearch` to set this up automatically. By default
    this will be `4`, but can be any number (you should pick a number <=
    the number of cores available).
- `procs::Union{Array{Int, 1}, Nothing}=nothing`: If you have set up
    a distributed run manually with `procs = addprocs()` and `@everywhere`,
    pass the `procs` to this keyword argument.
- `multithreading::Bool=false`: Whether to use multithreading. Otherwise,
    will use multiprocessing. Multithreading uses less memory, but multiprocessing
    can handle multi-node compute.
- `runtests::Bool=true`: Whether to run (quick) tests before starting the
    search, to see if there will be any problems during the equation search
    related to the host environment.
- `saved_state::Union{StateType, Nothing}=nothing`: If you have already
    run `EquationSearch` and want to resume it, pass the state here.
    To get this to work, you need to have stateReturn=true in the options,
    which will cause `EquationSearch` to return the state. Note that
    you cannot change the operators or dataset, but most other options
    should be changeable.
- `addprocs_function::Union{Function, Nothing}=nothing`: If using distributed
    mode (`multithreading=false`), you may pass a custom function to use
    instead of `addprocs`. This function should take a single positional argument,
    which is the number of processes to use, as well as the `lazy` keyword argument.
    For example, if set up on a slurm cluster, you could pass
    `addprocs_function = addprocs_slurm`, which will set up slurm processes.

# Returns
- `hallOfFame::HallOfFame`: The best equations seen during the search.
    hallOfFame.members gives an array of `PopMember` objects, which
    have their tree (equation) stored in `.tree`. Their score (loss)
    is given in `.score`. The array of `PopMember` objects
    is enumerated by size from `1` to `options.maxsize`.
"""
function EquationSearch(
    X::AbstractMatrix{T},
    y::AbstractMatrix{T};
    niterations::Int=10,
    weights::Union{AbstractMatrix{T},AbstractVector{T},Nothing}=nothing,
    varMap::Union{Array{String,1},Nothing}=nothing,
    options::Options=Options(),
    numprocs::Union{Int,Nothing}=nothing,
    procs::Union{Array{Int,1},Nothing}=nothing,
    multithreading::Bool=false,
    runtests::Bool=true,
    saved_state::Union{StateType{T},Nothing}=nothing,
    addprocs_function::Union{Function,Nothing}=nothing,
) where {T<:Real}
    nout = size(y, FEATURE_DIM)
    if weights !== nothing
        weights = reshape(weights, size(y))
    end
    datasets = [
        Dataset(
            X,
            y[j, :];
            weights=(weights === nothing ? weights : weights[j, :]),
            varMap=varMap,
        ) for j in 1:nout
    ]

    return EquationSearch(
        datasets;
        niterations=niterations,
        options=options,
        numprocs=numprocs,
        procs=procs,
        multithreading=multithreading,
        runtests=runtests,
        saved_state=saved_state,
        addprocs_function=addprocs_function,
    )
end

function EquationSearch(
    X::AbstractMatrix{T1}, y::AbstractMatrix{T2}; kw...
) where {T1<:Real,T2<:Real}
    U = promote_type(T1, T2)
    return EquationSearch(
        convert(AbstractMatrix{U}, X), convert(AbstractMatrix{U}, y); kw...
    )
end

function EquationSearch(
    X::AbstractMatrix{T1}, y::AbstractVector{T2}; kw...
) where {T1<:Real,T2<:Real}
    return EquationSearch(X, reshape(y, (1, size(y, 1))); kw...)
end

function EquationSearch(
    datasets::Array{Dataset{T},1};
    niterations::Int=10,
    options::Options=Options(),
    numprocs::Union{Int,Nothing}=nothing,
    procs::Union{Array{Int,1},Nothing}=nothing,
    multithreading::Bool=false,
    runtests::Bool=true,
    saved_state::Union{StateType{T},Nothing}=nothing,
    addprocs_function::Union{Function,Nothing}=nothing,
) where {T<:Real}

    # Population(datasets[j], baselineMSEs[j], npop=options.npop,
    #            nlength=3, options=options, nfeatures=datasets[j].nfeatures),
    # HallOfFame(options),

    noprocs = (procs === nothing && numprocs == 0)
    someprocs = !noprocs

    concurrency = if multithreading
        @assert procs === nothing && numprocs in [0, nothing]
        SRThreaded()
    elseif someprocs
        SRDistributed()
    else #noprocs, multithreading=false
        SRSerial()
    end

    return _EquationSearch(
        concurrency,
        datasets;
        niterations=niterations,
        options=options,
        numprocs=numprocs,
        procs=procs,
        runtests=runtests,
        saved_state=saved_state,
        addprocs_function=addprocs_function,
    )
end

function _EquationSearch(
    ::ConcurrencyType,
    datasets::Array{Dataset{T},1};
    niterations::Int=10,
    options::Options=Options(),
    numprocs::Union{Int,Nothing}=nothing,
    procs::Union{Array{Int,1},Nothing}=nothing,
    runtests::Bool=true,
    saved_state::Union{StateType{T},Nothing}=nothing,
    addprocs_function::Union{Function,Nothing}=nothing,
) where {T<:Real,ConcurrencyType<:SRConcurrency}
    if options.deterministic
        if ConcurrencyType != SRSerial
            error("Determinism is only guaranteed for serial mode.")
        end
    end

    can_read_input = true
    try
        Base.start_reading(stdin)
        bytes = bytesavailable(stdin)
        if bytes > 0
            # Clear out initial data
            read(stdin, bytes)
        end
    catch err
        if isa(err, MethodError)
            can_read_input = false
        else
            throw(err)
        end
    end

    # Redefine print, show:
    @eval begin
        function Base.print(io::IO, tree::Node)
            return print(io, string_tree(tree, $options; varMap=$(datasets[1].varMap)))
        end
        function Base.show(io::IO, tree::Node)
            return print(io, string_tree(tree, $options; varMap=$(datasets[1].varMap)))
        end
    end

    example_dataset = datasets[1]
    nout = size(datasets, 1)

    if runtests
        test_option_configuration(T, options)
        # Testing the first output variable is the same:
        test_dataset_configuration(example_dataset, options)
    end

    if example_dataset.weighted
        avgys = [
            sum(dataset.y .* dataset.weights) / sum(dataset.weights) for dataset in datasets
        ]
        baselineMSEs = [
            loss(dataset.y, ones(T, dataset.n) .* avgy, dataset.weights, options) for
            dataset in datasets, avgy in avgys
        ]
    else
        avgys = [sum(dataset.y) / dataset.n for dataset in datasets]
        baselineMSEs = [
            loss(dataset.y, ones(T, dataset.n) .* avgy, options) for
            (dataset, avgy) in zip(datasets, avgys)
        ]
    end

    if options.seed !== nothing
        seed!(options.seed)
    end
    # Start a population on every process
    #    Store the population, hall of fame
    allPopsType = if ConcurrencyType == SRSerial
        Tuple{Population,HallOfFame,RecordType,Float64}
    elseif ConcurrencyType == SRDistributed
        Future
    else
        Task
    end

    allPops = [allPopsType[] for j in 1:nout]
    init_pops = [allPopsType[] for j in 1:nout]
    # Set up a channel to send finished populations back to head node
    if ConcurrencyType in [SRDistributed, SRThreaded]
        if ConcurrencyType == SRDistributed
            channels = [
                [RemoteChannel(1) for i in 1:(options.npopulations)] for j in 1:nout
            ]
        else
            channels = [[Channel(1) for i in 1:(options.npopulations)] for j in 1:nout]
        end
        tasks = [Task[] for j in 1:nout]
    end

    # This is a recorder for populations, but is not actually used for processing, just
    # for the final return.
    returnPops = [
        [
            Population(
                datasets[j],
                baselineMSEs[j];
                npop=1,
                options=options,
                nfeatures=datasets[j].nfeatures,
            ) for i in 1:(options.npopulations)
        ] for j in 1:nout
    ]
    # These initial populations are discarded:
    bestSubPops = [
        [
            Population(
                datasets[j],
                baselineMSEs[j];
                npop=1,
                options=options,
                nfeatures=datasets[j].nfeatures,
            ) for i in 1:(options.npopulations)
        ] for j in 1:nout
    ]
    if saved_state === nothing
        hallOfFame = [HallOfFame(options) for j in 1:nout]
    else
        hallOfFame = saved_state[2]::Union{HallOfFame,Vector{HallOfFame}}
        if !isa(hallOfFame, Vector{HallOfFame})
            hallOfFame = [hallOfFame]
        end
        hallOfFame::Vector{HallOfFame}
    end
    actualMaxsize = options.maxsize + MAX_DEGREE

    # Make frequencyComplexities a moving average, rather than over all time:
    window_size = 100000
    smallest_frequency_allowed = 1
    # 3 here means this will last 3 cycles before being "refreshed"
    # We start out with even numbers at all frequencies.
    frequencyComplexities = [
        ones(T, actualMaxsize) * window_size / actualMaxsize for i in 1:nout
    ]

    curmaxsizes = [3 for j in 1:nout]
    record = RecordType("options" => "$(options)")

    if options.warmupMaxsizeBy == 0.0f0
        curmaxsizes = [options.maxsize for j in 1:nout]
    end

    # Records the number of evaluations:
    # Real numbers indicate use of batching.
    num_evals = [[0.0 for i in 1:(options.npopulations)] for j in 1:nout]

    we_created_procs = false
    ##########################################################################
    ### Distributed code:
    ##########################################################################
    if ConcurrencyType == SRDistributed
        if addprocs_function === nothing
            addprocs_function = addprocs
        end
        if numprocs === nothing && procs === nothing
            numprocs = 4
            procs = addprocs_function(4; lazy=false)
            we_created_procs = true
        elseif numprocs === nothing
            numprocs = length(procs)
        elseif procs === nothing
            procs = addprocs_function(numprocs; lazy=false)
            we_created_procs = true
        end

        if we_created_procs
            project_path = splitdir(Pkg.project().path)[1]
            activate_env_on_workers(procs, project_path, options)
            import_module_on_workers(procs, @__FILE__, options)
        end
        move_functions_to_workers(procs, options, example_dataset)
        if runtests
            test_module_on_workers(procs, options)
        end

        if runtests
            test_entire_pipeline(procs, example_dataset, options)
        end
    end
    # Get the next worker process to give a job:
    worker_assignment = Dict{Tuple{Int,Int},Int}()

    for j in 1:nout
        for i in 1:(options.npopulations)
            worker_idx = next_worker(worker_assignment, procs)
            if ConcurrencyType == SRDistributed
                worker_assignment[(j, i)] = worker_idx
            end
            if saved_state === nothing
                new_pop = @sr_spawner ConcurrencyType worker_idx (
                    Population(
                        datasets[j],
                        baselineMSEs[j];
                        npop=options.npop,
                        nlength=3,
                        options=options,
                        nfeatures=datasets[j].nfeatures,
                    ),
                    HallOfFame(options),
                    RecordType(),
                    Float64(options.npop),
                )
                # This involves npop evaluations, on the full dataset:
            else
                is_vector = typeof(saved_state[1]) <: Vector{Vector{Population{T}}}
                cur_saved_state = is_vector ? saved_state[1][j][i] : saved_state[1][j, i]

                if length(cur_saved_state.members) >= options.npop
                    new_pop = @sr_spawner ConcurrencyType worker_idx (
                        cur_saved_state, HallOfFame(options), RecordType(), 0.0
                    )
                else
                    # If population has not yet been created (e.g., exited too early)
                    println(
                        "Warning: recreating population (output=$(j), population=$(i)), as the saved one only has $(length(cur_saved_state.members)) members.",
                    )
                    new_pop = @sr_spawner ConcurrencyType worker_idx (
                        Population(
                            datasets[j],
                            baselineMSEs[j];
                            npop=options.npop,
                            nlength=3,
                            options=options,
                            nfeatures=datasets[j].nfeatures,
                        ),
                        HallOfFame(options),
                        RecordType(),
                        Float64(options.npop),
                    )
                end
            end
            push!(init_pops[j], new_pop)
        end
    end
    # 2. Start the cycle on every process:
    for j in 1:nout
        dataset = datasets[j]
        baselineMSE = baselineMSEs[j]
        frequencyComplexity = frequencyComplexities[j]
        curmaxsize = curmaxsizes[j]
        for i in 1:(options.npopulations)
            @recorder record["out$(j)_pop$(i)"] = RecordType()
            worker_idx = next_worker(worker_assignment, procs)
            if ConcurrencyType == SRDistributed
                worker_assignment[(j, i)] = worker_idx
            end

            # TODO - why is this needed??
            # Multi-threaded doesn't like to fetch within a new task:
            updated_pop = @sr_spawner ConcurrencyType worker_idx let
                in_pop = if ConcurrencyType in [SRDistributed, SRThreaded]
                    fetch(init_pops[j][i])[1]
                else
                    init_pops[j][i][1]
                end

                cur_record = RecordType()
                @recorder cur_record["out$(j)_pop$(i)"] = RecordType(
                    "iteration0" => record_population(in_pop, options)
                )
                tmp_num_evals = 0.0
                tmp_pop, tmp_best_seen, evals_from_cycle = s_r_cycle(
                    dataset,
                    baselineMSE,
                    in_pop,
                    options.ncyclesperiteration,
                    curmaxsize,
                    copy(frequencyComplexity) / sum(frequencyComplexity);
                    verbosity=options.verbosity,
                    options=options,
                    record=cur_record,
                )
                tmp_num_evals += evals_from_cycle
                tmp_pop, evals_from_optimize = optimize_and_simplify_population(
                    dataset, baselineMSE, tmp_pop, options, curmaxsize, cur_record
                )
                tmp_num_evals += evals_from_optimize
                if options.batching
                    for i_member in 1:(options.maxsize + MAX_DEGREE)
                        score, result_loss = score_func(
                            dataset,
                            baselineMSE,
                            tmp_best_seen.members[i_member].tree,
                            options,
                        )
                        tmp_best_seen.members[i_member].score = score
                        tmp_best_seen.members[i_member].loss = result_loss
                        tmp_num_evals += 1
                    end
                end
                (tmp_pop, tmp_best_seen, cur_record, tmp_num_evals)
            end
            push!(allPops[j], updated_pop)
        end
    end

    debug(options.verbosity > 0 || options.progress, "Started!")
    start_time = time()
    total_cycles = options.npopulations * niterations
    cycles_remaining = [total_cycles for j in 1:nout]
    sum_cycle_remaining = sum(cycles_remaining)
    if options.progress && nout == 1
        #TODO: need to iterate this on the max cycles remaining!
        progress_bar = ProgressBar(1:sum_cycle_remaining; width=options.terminal_width)
        cur_cycle = nothing
        cur_state = nothing
    end

    last_print_time = time()
    num_equations = 0.0
    print_every_n_seconds = 0
    equation_speed = Float32[]

    if ConcurrencyType in [SRDistributed, SRThreaded]
        for j in 1:nout
            for i in 1:(options.npopulations)
                # Start listening for each population to finish:
                t = @async put!(channels[j][i], fetch(allPops[j][i]))
                push!(tasks[j], t)
            end
        end
    end

    # Randomly order which order to check populations:
    # This is done so that we do work on all nout equally.
    all_idx = [(j, i) for j in 1:nout for i in 1:(options.npopulations)]
    shuffle!(all_idx)
    kappa = 0
    head_node_time = Dict("occupied" => 0, "start" => time())
    while sum(cycles_remaining) > 0
        kappa += 1
        if kappa > options.npopulations * nout
            kappa = 1
        end
        # nout, npopulations:
        j, i = all_idx[kappa]

        # Check if error on population:
        if ConcurrencyType in [SRDistributed, SRThreaded]
            if istaskfailed(tasks[j][i])
                fetch(tasks[j][i])
                error("Task failed for population")
            end
        end
        # Non-blocking check if a population is ready:
        population_ready =
            ConcurrencyType in [SRDistributed, SRThreaded] ? isready(channels[j][i]) : true
        # Don't start more if this output has finished its cycles:
        # TODO - this might skip extra cycles?
        population_ready &= (cycles_remaining[j] > 0)
        if population_ready
            head_node_start_work = time()
            # Take the fetch operation from the channel since its ready
            (cur_pop, best_seen, cur_record, cur_num_evals) =
                if ConcurrencyType in [SRDistributed, SRThreaded]
                    take!(channels[j][i])
                else
                    allPops[j][i]
                end
            returnPops[j][i] = cur_pop
            cur_pop::Population
            best_seen::HallOfFame
            cur_record::RecordType
            cur_num_evals::Float64
            bestSubPops[j][i] = best_sub_pop(cur_pop; topn=options.topn)
            @recorder record = recursive_merge(record, cur_record)
            num_evals[j][i] += cur_num_evals

            dataset = datasets[j]
            baselineMSE = baselineMSEs[j]
            curmaxsize = curmaxsizes[j]

            #Try normal copy...
            bestPops = Population([
                member for pop in bestSubPops[j] for member in pop.members
            ])

            ###################################################################
            # Hall Of Fame updating ###########################################
            for (i_member, member) in enumerate(
                Iterators.flatten((cur_pop.members, best_seen.members[best_seen.exists]))
            )
                part_of_cur_pop = i_member <= length(cur_pop.members)
                size = compute_complexity(member.tree, options)

                if part_of_cur_pop
                    if size <= options.maxsize
                        frequencyComplexities[j][size] += 1
                    end
                end
                actualMaxsize = options.maxsize + MAX_DEGREE

                valid_size = size < actualMaxsize
                if valid_size
                    already_filled = hallOfFame[j].exists[size]
                    better_than_current = member.score < hallOfFame[j].members[size].score
                    if !already_filled || better_than_current
                        hallOfFame[j].members[size] = copy_pop_member(member)
                        hallOfFame[j].exists[size] = true
                    end
                end
            end
            ###################################################################

            # Dominating pareto curve - must be better than all simpler equations
            dominating = calculate_pareto_frontier(dataset, hallOfFame[j], options)
            hofFile = options.hofFile
            if nout > 1
                hofFile = hofFile * ".out$j"
            end
            open(hofFile, "w") do io
                println(io, "Complexity|MSE|Equation")
                for member in dominating
                    println(
                        io,
                        "$(compute_complexity(member.tree, options))|$(member.loss)|$(string_tree(member.tree, options, varMap=dataset.varMap))",
                    )
                end
            end
            cp(hofFile, hofFile * ".bkup"; force=true)

            ###################################################################
            # Migration #######################################################
            # Try normal copy otherwise.
            if options.migration
                for k in rand(
                    1:(options.npop), round(Int, options.npop * options.fractionReplaced)
                )
                    to_copy = rand(1:size(bestPops.members, 1))

                    # Explicit copy here resets the birth.
                    cur_pop.members[k] = PopMember(
                        copy_node(bestPops.members[to_copy].tree),
                        copy(bestPops.members[to_copy].score),
                        copy(bestPops.members[to_copy].loss);
                        ref=copy(bestPops.members[to_copy].ref),
                        parent=copy(bestPops.members[to_copy].parent),
                        deterministic=options.deterministic,
                    )
                    # TODO: Clean this up using copy_pop_member.
                end
            end

            if options.hofMigration && size(dominating, 1) > 0
                for k in rand(
                    1:(options.npop), round(Int, options.npop * options.fractionReplacedHof)
                )
                    # Copy in case one gets used twice
                    to_copy = rand(1:size(dominating, 1))
                    cur_pop.members[k] = PopMember(
                        copy_node(dominating[to_copy].tree),
                        copy(dominating[to_copy].score),
                        copy(dominating[to_copy].loss);
                        ref=copy(dominating[to_copy].ref),
                        parent=copy(dominating[to_copy].parent),
                        deterministic=options.deterministic,
                    )
                    # TODO: Clean this up with copy_pop_member.
                end
            end
            ###################################################################

            cycles_remaining[j] -= 1
            if cycles_remaining[j] == 0
                break
            end
            worker_idx = next_worker(worker_assignment, procs)
            if ConcurrencyType == SRDistributed
                worker_assignment[(j, i)] = worker_idx
            end
            @recorder begin
                key = "out$(j)_pop$(i)"
                iteration = find_iteration_from_record(key, record) + 1
            end

            allPops[j][i] = @sr_spawner ConcurrencyType worker_idx let
                cur_record = RecordType()
                @recorder cur_record[key] = RecordType(
                    "iteration$(iteration)" => record_population(cur_pop, options)
                )
                tmp_num_evals = 0.0
                tmp_pop, tmp_best_seen, evals_from_cycle = s_r_cycle(
                    dataset,
                    baselineMSE,
                    cur_pop,
                    options.ncyclesperiteration,
                    curmaxsize,
                    copy(frequencyComplexities[j]) / sum(frequencyComplexities[j]);
                    verbosity=options.verbosity,
                    options=options,
                    record=cur_record,
                )
                tmp_num_evals += evals_from_cycle
                tmp_pop, evals_from_optimize = optimize_and_simplify_population(
                    dataset, baselineMSE, tmp_pop, options, curmaxsize, cur_record
                )
                tmp_num_evals += evals_from_optimize

                # Update scores if using batching:
                if options.batching
                    for i_member in 1:(options.maxsize + MAX_DEGREE)
                        if tmp_best_seen.exists[i_member]
                            score, result_loss = score_func(
                                dataset,
                                baselineMSE,
                                tmp_best_seen.members[i_member].tree,
                                options,
                            )
                            tmp_best_seen.members[i_member].score = score
                            tmp_best_seen.members[i_member].loss = result_loss
                            tmp_num_evals += 1
                        end
                    end
                end

                (tmp_pop, tmp_best_seen, cur_record, tmp_num_evals)
            end
            if ConcurrencyType in [SRDistributed, SRThreaded]
                tasks[j][i] = @async put!(channels[j][i], fetch(allPops[j][i]))
            end

            cycles_elapsed = total_cycles - cycles_remaining[j]
            if options.warmupMaxsizeBy > 0
                fraction_elapsed = 1.0f0 * cycles_elapsed / total_cycles
                if fraction_elapsed > options.warmupMaxsizeBy
                    curmaxsizes[j] = options.maxsize
                else
                    curmaxsizes[j] =
                        3 + floor(
                            Int,
                            (options.maxsize - 3) * fraction_elapsed /
                            options.warmupMaxsizeBy,
                        )
                end
            end
            num_equations += options.ncyclesperiteration * options.npop / 10.0

            if options.progress && nout == 1
                # set_postfix(iter, Equations=)
                equation_strings = string_dominating_pareto_curve(
                    hallOfFame[j], baselineMSE, datasets[j], options, avgys[j]
                )
                load_string =
                    @sprintf(
                        "Head worker occupation: %.1f",
                        100 * head_node_time["occupied"] /
                            (time() - head_node_time["start"])
                    ) * "%\n"
                # TODO - include command about "q" here.
                load_string *= @sprintf(
                    "Press 'q' and then <enter> to stop execution early.\n"
                )
                equation_strings = load_string * equation_strings
                set_multiline_postfix(progress_bar, equation_strings)
                if cur_cycle === nothing
                    (cur_cycle, cur_state) = iterate(progress_bar)
                else
                    (cur_cycle, cur_state) = iterate(progress_bar, cur_state)
                end
                sum_cycle_remaining = sum(cycles_remaining)
            end
            head_node_end_work = time()
            head_node_time["occupied"] += (head_node_end_work - head_node_start_work)

            # Move moving window along frequencyComplexities:
            cur_size_frequency_complexities = sum(frequencyComplexities[j])
            if cur_size_frequency_complexities > window_size
                difference_in_size = cur_size_frequency_complexities - window_size
                # We need frequencyComplexities to be positive, but also sum to a number.
                num_loops = 0
                # TODO: Clean this code up. Should not have to have
                # loop catching.
                while difference_in_size > 0
                    indices_to_subtract = findall(
                        frequencyComplexities[j] .> smallest_frequency_allowed
                    )
                    num_remaining = size(indices_to_subtract, 1)
                    amount_to_subtract = min(
                        difference_in_size / num_remaining,
                        min(frequencyComplexities[j][indices_to_subtract]...) -
                        smallest_frequency_allowed,
                    )
                    frequencyComplexities[j][indices_to_subtract] .-= amount_to_subtract
                    total_amount_to_subtract = amount_to_subtract * num_remaining
                    difference_in_size -= total_amount_to_subtract
                    num_loops += 1
                    if num_loops > 1000 || total_amount_to_subtract < 1e-6
                        # Sometimes, total_amount_to_subtract can be a very very small number.
                        break
                    end
                end
            end
        end
        sleep(1e-6)

        ################################################################
        ## Printing code
        elapsed = time() - last_print_time
        #Update if time has passed, and some new equations generated.
        if elapsed > print_every_n_seconds && num_equations > 0.0
            # Dominating pareto curve - must be better than all simpler equations
            current_speed = num_equations / elapsed
            average_over_m_measurements = 10 #for print_every...=5, this gives 50 second running average
            push!(equation_speed, current_speed)
            if length(equation_speed) > average_over_m_measurements
                deleteat!(equation_speed, 1)
            end
            if (options.verbosity > 0) || (options.progress && nout > 1)
                @printf("\n")
                average_speed = sum(equation_speed) / length(equation_speed)
                @printf("Cycles per second: %.3e\n", round(average_speed, sigdigits=3))
                @printf("Time since last print: %.3f seconds\n", elapsed)
                @printf(
                    "Head worker occupation: %.1f%%\n",
                    100 * head_node_time["occupied"] / (time() - head_node_time["start"])
                )
                cycles_elapsed = total_cycles * nout - sum(cycles_remaining)
                @printf(
                    "Progress: %d / %d total iterations (%.3f%%)\n",
                    cycles_elapsed,
                    total_cycles * nout,
                    100.0 * cycles_elapsed / total_cycles / nout
                )

                @printf("==============================\n")
                for j in 1:nout
                    if nout > 1
                        @printf("Best equations for output %d\n", j)
                    end
                    equation_strings = string_dominating_pareto_curve(
                        hallOfFame[j], baselineMSEs[j], datasets[j], options, avgys[j]
                    )
                    print(equation_strings)
                    @printf("==============================\n")

                    # Debugging code for frequencyComplexities:
                    # for size_i=1:actualMaxsize
                    #     @printf("frequencyComplexities size %d = %.2f\n", size_i, frequencyComplexities[j][size_i])
                    # end
                end
                @printf("Press 'q' and then <enter> to stop execution early.\n")
            end
            last_print_time = time()
            num_equations = 0.0
        end
        ################################################################

        ################################################################
        ## Early stopping code
        if options.earlyStopCondition !== nothing
            # Check if all nout are below stopping condition.
            all_below = true
            for j in 1:nout
                dominating = calculate_pareto_frontier(datasets[j], hallOfFame[j], options)
                # Check if zero size:
                if length(dominating) == 0
                    all_below = false
                elseif !any([
                    options.earlyStopCondition(
                        member.loss, compute_complexity(member.tree, options)
                    ) for member in dominating
                ])
                    # None of the equations meet the stop condition.
                    all_below = false
                end

                if !all_below
                    break
                end
            end
            if all_below # Early stop!
                break
            end
        end
        ################################################################

        ################################################################
        ## Signal stopping code
        if can_read_input
            bytes = bytesavailable(stdin)
            if bytes > 0
                # Read:
                data = read(stdin, bytes)
                control_c = 0x03
                quit = 0x71
                if length(data) > 1 && (data[end] == control_c || data[end - 1] == quit)
                    break
                end
            end
        end
        ################################################################

        ################################################################
        ## Timeout stopping code
        if options.timeout_in_seconds !== nothing
            if time() - start_time > options.timeout_in_seconds
                break
            end
        end
        ################################################################

        ################################################################
        ## num_evals stopping code
        if options.max_evals !== nothing
            if options.max_evals <= sum(sum, num_evals)
                break
            end
        end
        ################################################################
    end
    if can_read_input
        Base.stop_reading(stdin)
    end
    if we_created_procs
        rmprocs(procs)
    end
    # TODO - also stop threads here.

    ##########################################################################
    ### Distributed code^
    ##########################################################################

    @recorder begin
        open(options.recorder_file, "w") do io
            JSON3.write(io, record; allow_inf=true)
        end
    end

    if options.stateReturn
        state = (returnPops, (nout == 1 ? hallOfFame[1] : hallOfFame))
        state::StateType{T}
        return state
    else
        if nout == 1
            return hallOfFame[1]
        else
            return hallOfFame
        end
    end
end

end #module SR


text:
using DemoPackagePM
using Documenter

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

makedocs(;
    modules=[DemoPackagePM],
    authors="Pingchuan Ma <mapingchuan@wustl.edu> and contributors",
    repo="https://github.com/pcmaaaaa/DemoPackagePM.jl/blob/{commit}{path}#{line}",
    sitename="DemoPackagePM.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://pcmaaaaa.github.io/DemoPackagePM.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/pcmaaaaa/DemoPackagePM.jl",
    devbranch="main",
)


text:
import Base.+, Base.-, Base./, Base.*

function +(spl::Schumaker, num::Real)
    new_coefficients_ = hcat(spl.coefficient_matrix_[:,1:2],  spl.coefficient_matrix_[:,3] .+ num)
    return Schumaker(spl.IntStarts_, new_coefficients_)
end
function -(spl::Schumaker, num::Real)
    new_coefficients_ = hcat(spl.coefficient_matrix_[:,1:2],  spl.coefficient_matrix_[:,3] .- num)
    return Schumaker(spl.IntStarts_, new_coefficients_)
end
function *(spl::Schumaker, num::Real)
    new_coefficients_ = spl.coefficient_matrix_  .* num
    return Schumaker(spl.IntStarts_, new_coefficients_)
end
function /(spl::Schumaker, num::Real)
    new_coefficients_ = spl.coefficient_matrix_  ./ num
    return Schumaker(spl.IntStarts_, new_coefficients_)
end
function +(num::Real, spl::Schumaker)
    return +(spl,num)
end
function -(num::Real, spl::Schumaker)
    return (-1)*spl + num
end
function *(num::Real, spl::Schumaker)
    return *(spl,num)
end


text:
using Pkg

Pkg.activate(@__DIR__)
Pkg.instantiate()

using PackageCompiler

cd(@__DIR__)

create_sysimage([:RoME,:RoMEPlotting], sysimage_path="RoMEPlottingSysimage.so", precompile_execution_file="precompile_triggers_plotting.jl")

## to use RoME and RoMEPlotting with the newly created sysimage, start julia with:
# julia -J RoMEPlottingSysimage.so
# February 2021: starting julia with the sysimage crashes


text:
# Copyright (c) 2014: SCS.jl contributors
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module TestSCS

using Test
using MathOptInterface
import SCS

const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

test_DirectSolver() = _test_runtests(SCS.DirectSolver)

test_IndirectSolver() = _test_runtests(SCS.IndirectSolver)

function _test_runtests(linear_solver)
    optimizer = SCS.Optimizer()
    MOI.set(
        optimizer,
        MOI.RawOptimizerAttribute("linear_solver"),
        linear_solver,
    )
    MOI.set(optimizer, MOI.RawOptimizerAttribute("eps_abs"), 1e-6)
    MOI.set(optimizer, MOI.Silent(), true)
    model = MOI.Bridges.full_bridge_optimizer(
        MOI.Utilities.CachingOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
            optimizer,
        ),
        Float64,
    )
    MOI.Test.runtests(
        model,
        MOI.Test.Config(
            atol = 1e-2,
            exclude = Any[
                MOI.ConstraintBasisStatus,
                MOI.VariableBasisStatus,
                MOI.ConstraintName,
                MOI.VariableName,
                MOI.ObjectiveBound,
            ],
        ),
        exclude = String[
            # Unexpected failures:
            #   TODO(odow): looks like a tolerance issue?
            "test_linear_add_constraints",
            # Expected test failures:
            #   TODO(odow): get not supported for primal/dual starts
            "test_model_ModelFilter_AbstractConstraintAttribute",
            #   Unimplemented feature
            "test_attribute_SolverVersion",
            #   ArgumentError: The number of constraints must be greater than 0
            "test_attribute_RawStatusString",
            "test_attribute_SolveTimeSec",
            "test_objective_ObjectiveFunction_blank",
            "test_solve_TerminationStatus_DUAL_INFEASIBLE",
            #   Problem is a nonconvex QP
            "test_basic_ScalarQuadraticFunction_EqualTo",
            "test_basic_ScalarQuadraticFunction_GreaterThan",
            "test_basic_ScalarQuadraticFunction_Interval",
            "test_basic_VectorQuadraticFunction_",
            "test_quadratic_SecondOrderCone_basic",
            "test_quadratic_nonconvex_",
            #   MathOptInterface.jl issue #1431
            "test_model_LowerBoundAlreadySet",
            "test_model_UpperBoundAlreadySet",
        ],
    )
    return
end

function test_RawOptimizerAttribute()
    model = SCS.Optimizer()
    MOI.set(model, MOI.RawOptimizerAttribute("eps_abs"), 1.0)
    @test MOI.get(model, MOI.RawOptimizerAttribute("eps_abs")) == 1.0
    @test MOI.get(model, MOI.RawOptimizerAttribute("eps_abs")) == 1.0
    MOI.set(model, MOI.RawOptimizerAttribute("eps_abs"), 2.0)
    @test MOI.get(model, MOI.RawOptimizerAttribute("eps_abs")) == 2.0
    @test MOI.get(model, MOI.RawOptimizerAttribute("eps_abs")) == 2.0
end

function test_constrained_variables()
    model = MOI.Bridges.full_bridge_optimizer(
        MOI.Utilities.CachingOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
            SCS.Optimizer(),
        ),
        Float64,
    )
    @test MOI.supports_constraint(
        model,
        MOI.VectorOfVariables,
        MOI.PositiveSemidefiniteConeTriangle,
    )
    x = MOI.add_variables(model, 6)
    f = MOI.VectorOfVariables(x)
    s = MOI.PositiveSemidefiniteConeTriangle(3)
    @test isa(
        MOI.add_constraint(model, f, s),
        MOI.ConstraintIndex{typeof(f),typeof(s)},
    )
    return
end

function test_unsupported()
    model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
    optimizer = SCS.Optimizer()
    x = MOI.add_variable(model)
    MOI.add_constraint(model, 1.0x, MOI.EqualTo(1.0))
    err = MOI.UnsupportedConstraint{
        MOI.ScalarAffineFunction{Float64},
        MOI.EqualTo{Float64},
    }()
    @test_throws err MOI.optimize!(optimizer, model)
    MOI.empty!(model)
    x = MOI.add_variable(model)
    MOI.set(model, MOI.Test.UnknownVariableAttribute(), x, 1.0)
    err = MOI.UnsupportedAttribute{MOI.Test.UnknownVariableAttribute}
    @test_throws err MOI.optimize!(optimizer, model)
    MOI.empty!(model)
    x = MOI.add_variable(model)
    c = MOI.add_constraint(model, MOI.Utilities.vectorize([1.0x]), MOI.Zeros(1))
    MOI.set(model, MOI.Test.UnknownConstraintAttribute(), c, 1.0)
    err = MOI.UnsupportedAttribute{MOI.Test.UnknownConstraintAttribute}
    @test_throws err MOI.optimize!(optimizer, model)
end

function test_empty_problem()
    model = MOI.Utilities.Model{Float64}()
    scs = SCS.Optimizer()
    MOI.optimize!(scs, model)
    @test MOI.get(scs, MOI.TerminationStatus()) == MOI.INVALID_MODEL
    @test MOI.get(scs, MOI.PrimalStatus()) == MOI.NO_SOLUTION
    @test MOI.get(scs, MOI.DualStatus()) == MOI.NO_SOLUTION
    return
end

function test_conic_no_variables()
    model = MOI.Utilities.Model{Float64}()
    scs = SCS.Optimizer()
    f = MOI.VectorAffineFunction(
        MOI.VectorAffineTerm{Float64}[],
        [1.0, 0.5, 0.5],
    )
    MOI.add_constraint(model, f, MOI.SecondOrderCone(3))
    MOI.optimize!(scs, model)
    @test MOI.get(scs, MOI.TerminationStatus()) == MOI.INVALID_MODEL
    @test MOI.get(scs, MOI.PrimalStatus()) == MOI.NO_SOLUTION
    @test MOI.get(scs, MOI.DualStatus()) == MOI.NO_SOLUTION
    return
end

function test_Name_skip()
    model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}())
    MOI.set(model, MOI.Name(), "My problem")
    x = MOI.add_variables(model, 3)
    f = MOI.Utilities.operate(vcat, Float64, 1.0 .* x...)
    MOI.add_constraint(model, f, MOI.Nonnegatives(3))
    scs = SCS.Optimizer()
    MOI.optimize!(scs, model)
    @test MOI.get(scs, MOI.TerminationStatus()) == MOI.OPTIMAL
    return
end

end  # module

TestSCS.runtests()


text:
module Rasterization
    using Dates: length
    using Base.Threads
    using LinearAlgebra
    using Dates
    using Printf
    using NetCDF
    using Base
    using SeisSolXDMF
    using ..Util
    using ..NC
    using ..Args

    export rasterize, ZRange

    @enum ZRange z_all=-1 z_floor=0 z_surface=1

    const EDGE_TOL = .0      # Tolerance for whether to consider points that are close to a simplex edge/face to be on that edge/face
    const Z_RANGE_TOL = .001 # Tolerance for points close to the z_range-limit. Points over the limit by less than Z_RANGE_TOL are not discarded

    # Aliases for common array indices
    const X = 1
    const Y = 2
    const Z = 3
    const MIN = 1
    const MAX = 2

    struct HesseNormalForm
        n0  :: NTuple{3, Float64}
        d   :: Float64
    end

    """
    Contains all static information necessary for rasterization, as well as some computed values for convenience.
    """
    struct RasterizationContext
        simplices                   :: AbstractArray{Integer, 2}
        points                      :: AbstractArray{AbstractFloat, 2}

        z_range                     :: ZRange
        domain                      :: NTuple{3, NTuple{2, Float64}} # (x, y, z) with (min, max) respectively
        samples                     :: NTuple{3, UInt}
        sampling_rate               :: NTuple{3, Float64}
        dyn_var_mapping             :: VarMapping
        stat_var_mapping            :: VarMapping
        in_vars_dyn                 :: Array{String, 1}
        out_vars_stat               :: Array{String, 1}
        out_vars_dyn                :: Array{String, 1}

        n_threads                   :: Int
        n_simplex_points            :: Int
        n_dims                      :: Int
        n_simplices                 :: Int
        n_in_vars                   :: Int
        n_out_vars_dyn              :: Int
        n_out_vars_stat             :: Int
        n_timesteps                 :: Int
        n_timesteps_per_iteration   :: Int
        n_iterations                :: Int

        simplices_per_thread        :: Float64

        tanioka                     :: Bool
        water_height                :: Float64

        t_begin                     :: UInt
        t_end                       :: UInt
        times                       :: AbstractArray{Float64, 1}

        out_filename                :: AbstractString
    end

    """
    Contains buffers that are rewritten during rasterization, as well as binning info.
    """
    struct IterationBuffers
        bin_ids
        bin_counts
        out_grids_dyn               :: Dict{String, Array{Float64, 3}}
        out_grids_stat              :: Dict{String, Array{Float64, 2}}
        out_sample_counts           :: Array{UInt16, 2}
        prefetched_vars             :: Dict{String, Array{Float64, 2}}
    end

    """
    Returns the name of the simplex type currently being rasterized (triangle or tetrahedron), either in plural (default) or singular form.
    """
    @inline function simplex_name(ctx, plural=true)
        return ctx.n_dims == 3 ? (plural ? "simplices" : "simplex") : (plural ? "triangles" : "triangle")
    end

    function RasterizationContext(simplices, points, sampling_rate, mem_limit, dyn_var_mapping, stat_var_mapping, z_range,
                                  tanioka, water_height, t_begin, t_end, times, out_filename, custom_domain; max_threads=typemax(Int))

        custom_domain                   = (custom_domain..., (-Inf, Inf)) # Add z-range
        domain                          = Tuple( (max(minimum(points[i, :]), custom_domain[i][MIN]),
                                                  min(maximum(points[i, :]), custom_domain[i][MAX]))                for i in X:Z)

        samples                         = Tuple( ceil(UInt, (domain[i][MAX] - domain[i][MIN]) / sampling_rate[i])   for i in X:Z)

        #============================================#
        # Calculate problem size & hardware limitations
        #============================================#

        n_threads                       = min(max_threads, nthreads())

        n_simplex_points                = size(simplices, 1)
        n_dims                          = n_simplex_points - 1
        n_simplices                     = size(simplices, 2)
        simplices_per_thread            = n_simplices / n_threads # float on purpose

        b_mem_points                    = sizeof(Float64) * 3 * size(points, 2)
        b_mem_simps                     = sizeof(Int32) * size(simplices, 1) * size(simplices, 2)
        b_mem_cnt                       = sizeof(UInt16)  * samples[X] * samples[Y]
        b_mem_misc                      = 16 * 1024^2 * n_threads # 16 MiB per thread
        b_mem_misc                     += samples[X] * sizeof(Float64)
        b_mem_misc                     += n_simplices * sizeof(UInt16)
        b_mem_per_out_var_and_timestep  = sizeof(Float64) * samples[X] * samples[Y]
        b_mem_per_in_var_and_timestep   = sizeof(Float64) * n_simplices

        # static vars are currently only bathymetry and thus extracted from grid geometry. Thus, no in_vars_stat.
        in_vars_dyn                     = collect(keys(dyn_var_mapping))

        # Tanioka's method also needs horizontal displacements, add them if not present yet.
        if tanioka
            for var_name ∈ ("U", "V")
                var_name ∈ in_vars_dyn || push!(in_vars_dyn, var_name)
            end
        end

        n_in_vars                       = length(in_vars_dyn)

        out_vars_stat                   = z_range != z_all ? collect(values(stat_var_mapping)) : []
        out_vars_dyn                    = collect(values(dyn_var_mapping))

        n_out_vars_dyn                  = length(dyn_var_mapping)
        n_out_vars_stat                 = length(out_vars_stat)

        n_timesteps                     = t_end - t_begin + 1
        n_timesteps_per_iteration       = ((mem_limit - b_mem_points - b_mem_simps - b_mem_cnt - b_mem_misc - n_out_vars_stat * b_mem_per_out_var_and_timestep)
                                        ÷  (n_in_vars * b_mem_per_in_var_and_timestep + n_out_vars_dyn * b_mem_per_out_var_and_timestep))
        n_timesteps_per_iteration       = max(1, n_timesteps_per_iteration) # Process at least 1 timestep

        n_iterations                    = ceil(Int, n_timesteps / n_timesteps_per_iteration)

        return RasterizationContext(simplices, points,
                                    z_range, domain, samples, sampling_rate,
                                    dyn_var_mapping, stat_var_mapping,
                                    in_vars_dyn, out_vars_stat, out_vars_dyn,
                                    n_threads, n_simplex_points,
                                    n_dims, n_simplices,
                                    n_in_vars, n_out_vars_dyn,
                                    n_out_vars_stat, n_timesteps,
                                    n_timesteps_per_iteration, n_iterations,
                                    simplices_per_thread,
                                    tanioka, water_height, t_begin, t_end, times,
                                    out_filename)
    end

    function IterationBuffers(bin_ids, bin_counts, ctx :: RasterizationContext)
        # These are the grids that will be written to the NetCDF output later on.
        out_grids_dyn = Dict{String, Array{Float64, 3}}()
        for var_name ∈ ctx.out_vars_dyn
            out_grids_dyn[var_name] = Array{Float64, 3}(undef, (ctx.samples[X:Y]..., ctx.n_timesteps_per_iteration))
        end

        out_grids_stat = Dict{String, Array{Float64, 2}}()
        for var_name ∈ ctx.out_vars_stat
            out_grids_stat[var_name] = Array{Float64, 2}(undef, ctx.samples[X:Y])
        end

        # This matrix counts the total number of samples in each output grid cell
        # (read: the total number of sampled points in z-direction for each xy-index of the output grids)
        # This stays constant across all vars and timesteps because the grid geometry is constant.
        out_sample_counts = Array{UInt16, 2}(undef, ctx.samples[X:Y])

        # Variable prefetch buffers, get rewritten every iteration
        prefetched_vars = Dict{String, Array{Float64, 2}}()
        for var_name ∈ ctx.in_vars_dyn
            prefetched_vars[var_name] = Array{Float64, 2}(undef, (ctx.n_simplices, ctx.n_timesteps_per_iteration))
        end

        return IterationBuffers(bin_ids, bin_counts, out_grids_dyn, out_grids_stat, out_sample_counts, prefetched_vars)
    end

    """
        Rasterize variables stored in an unstructured grid into a regular grid. Write the results to a NetCDF file.

    # Arguments
    - `xdmf`:                   The XDMF file handle of the file to be rasterized
    - `dyn_var_mapping`:        Mapping of input to output names for dynamic variables in the current rasterization step.
    - `stat_var_mapping`:       Analogous, but for static variables.
    - `times`:                  An array of timestamps corresponding to the timesteps
    - `sampling_rate`:          A tuple (Δx, Δy, Δz) defining the regular grid cell size.
    - `out_filename`:           The path/name of the output NetCDF file
    - `mem_limit`:              The soft memory (in bytes) the software is allowed to utilize (soft limit!)
    - `z_range`:                A filter for discarding simplices in a certain mesh region. `z_floor` only keeps bathymetry, `z_surface` only keeps the sea surface,
                                `z_all` keeps all simplices.
    - `t_begin`:                The first timestep to be rasterized
    - `t_end`:                  The last timestep to be rasterized
    - `create_nc_dyn_vars`:     If not empty: create output file containing all value sides of the var mapping as variables. Existing file is overwritten.
                                Else: use existing file. If given, should contain all variable mappings of all rasterization passes.
    - `create_nc_stat_vars`:    Analogous to above.
    - `water_height`:           The offset to translate the water level to be at `z=0`.
    - `tanioka`:                Whether to apply Tanioka's method to the bathymetry. Requires U, V, W displacements to be present in the input XDMF file.
    """
    function rasterize(
        xdmf                :: XDMFFile,
        dyn_var_mapping     :: VarMapping,
        stat_var_mapping    :: VarMapping,
        times               :: AbstractArray{Float64, 1},
        sampling_rate       :: NTuple{3, Float64},
        out_filename        :: AbstractString,
        mem_limit           :: Int;
        z_range             :: ZRange  = z_all,
        t_begin             :: Integer = 1,
        t_end               :: Integer = length(times),
        create_nc_dyn_vars  :: Array{VarMapping, 1} = Array{VarMapping, 1}(),
        create_nc_stat_vars :: Array{VarMapping, 1} = Array{VarMapping, 1}(),
        water_height        :: Float64 = 0.,
        tanioka             :: Bool    = false,
        domain              :: DomainSize = ((-Inf, Inf),(-Inf, Inf)))

        simplices, points = grid_of(xdmf)

        ctx = RasterizationContext(simplices, points, sampling_rate, mem_limit, dyn_var_mapping, stat_var_mapping, z_range, tanioka,
                                   water_height, t_begin, t_end, times, out_filename, domain)

        #============================================#
        # Print rasterization settings
        #============================================#

        if (ctx.z_range == z_all) != (ctx.n_dims == 3)
            error("Z-Range 'all' can only be used with tetrahedra; Z-Ranges 'floor' and 'surface' can only be used with triangles.")
        end

        println("Domain is $(ctx.domain[X]) × $(ctx.domain[Y]) × $(ctx.domain[Z]).")

        if ctx.z_range != z_all
            s_region = ctx.z_range == z_floor ? "floor" : "surface"
            println("Only processing triangles on the $s_region.")
        end

        if ctx.n_dims == 3
            println("Sampling into $(ctx.samples[X]) × $(ctx.samples[Y]) × $(ctx.samples[Z]) cells.")
        else
            println("Sampling into $(ctx.samples[X]) × $(ctx.samples[Y]) cells.")
        end

        #============================================#
        # Set up NetCDF output file
        #============================================#

        if !isempty(create_nc_dyn_vars) || !isempty(create_nc_stat_vars)
            create_netcdf(ctx.out_filename, [ctx.domain[Y][MIN] + i * ctx.sampling_rate[Y] for i ∈ 1:ctx.samples[Y]],
                                  [ctx.domain[X][MIN] + i * ctx.sampling_rate[X] for i ∈ 1:ctx.samples[X]],
                                  ctx.times[ctx.t_begin:ctx.t_end],
                                  create_nc_stat_vars, create_nc_dyn_vars)
        else
            get_netcdf(ctx.out_filename)
        end

        #============================================#
        # Put simplices into bins to avoid thread
        # synchronization & lock contention later on.
        #============================================#

        bin_ids, bin_counts = bin(ctx)

        #============================================#
        # Setup output grids in memory
        #============================================#

        itbuf = IterationBuffers(bin_ids, bin_counts, ctx)

        #============================================#
        # Sample simplices and write to output
        #============================================#

        println("Processing timesteps ", ctx.t_begin-1, "-", ctx.t_end-1, " of 0-", length(times)-1, ".")
        println("RAM limit: ", human_readable_size(mem_limit), "; Can work on ",
                ctx.n_timesteps_per_iteration, " timesteps at once (therefore needing ",
                ctx.n_iterations , " iterations).")

        #============================================#
        # Iterate over and rasterize timesteps
        #============================================#

        iterate(ctx, itbuf, xdmf)

    end # function rasterize

    function bin(ctx)
        println("Binning $(ctx.n_simplices) $(simplex_name(ctx)) into $(ctx.n_threads)+$(ctx.n_threads - 1)+1 buckets...")

        n_excluded_domain = 0
        n_excluded_z_range = 0

        l_bin_start_idxs = [1 + round(Int, ctx.samples[X] / ctx.n_threads * bucket_id) for bucket_id ∈ 0:(ctx.n_threads-1)]

        #============================================#
        # Perform binning
        #============================================#

        # For each tetrahedron, assign a bin with an ID ∈ 1:2*n_threads based on its location
        l_bin_ids = Array{UInt16, 1}(undef, ctx.n_simplices)

        Threads.@threads for thread_id ∈ 1:ctx.n_threads
            thread_range_min = floor(Int, (thread_id - 1) * ctx.simplices_per_thread) + 1
            thread_range_max = thread_id == ctx.n_threads ? ctx.n_simplices : floor(Int, thread_id * ctx.simplices_per_thread)

            m3n_simp_points = Array{Float64, 2}(undef, (3, ctx.n_simplex_points))

            for simp_id ∈ thread_range_min:thread_range_max
                tetrahedron = (@view ctx.simplices[:, simp_id])
                for i ∈ 1:ctx.n_simplex_points
                    m3n_simp_points[:, i] = ctx.points[:, tetrahedron[i]]
                end

                coord_x_min = minimum(m3n_simp_points[X, :])
                coord_x_max = maximum(m3n_simp_points[X, :])
                coord_y_min = minimum(m3n_simp_points[Y, :])
                coord_y_max = maximum(m3n_simp_points[Y, :])

                if coord_x_max <= ctx.domain[X][MIN] || coord_x_min >= ctx.domain[X][MAX] || coord_y_max <= ctx.domain[Y][MIN] || coord_y_min >= ctx.domain[Y][MAX]
                    l_bin_ids[simp_id] = ctx.n_threads * 2 + 1
                    n_excluded_domain += 1
                    continue
                end

                # If simplex is not in the z-domain (only possible if z_range ∈ [z_floor, z_surface]), ignore it.
                if ctx.z_range != z_all
                    is_bath = is_bathy(m3n_simp_points, ctx.water_height)

                    if (ctx.z_range == z_floor) != is_bath
                        l_bin_ids[simp_id] = ctx.n_threads * 2 + 1 # Unused bin, will be ignored later
                        n_excluded_z_range += 1
                        continue
                    end
                else # This would be were tetrahedra in the volume mesh that belong to the ground are filtered out. For now: cutoff below -2000m
                    # TODO: correct z filter for tetrahedra
                    coord_z_min = minimum(m3n_simp_points[Z, :])
                    coord_z_max = maximum(m3n_simp_points[Z, :])
                    if coord_z_max <= 0.00001
                        l_bin_ids[simp_id] = ctx.n_threads * 2 + 1 # Unused bin, will be ignored later
                        n_excluded_z_range += 1
                        continue
                    end
                end

                idx_x_min = floor(Int,(coord_x_min - ctx.domain[X][MIN]) / ctx.sampling_rate[X]) + 1
                idx_x_max =  ceil(Int,(coord_x_max - ctx.domain[X][MIN]) / ctx.sampling_rate[X])

                bin_id_l = findlast(bin_start_idx -> bin_start_idx ≤ idx_x_min, l_bin_start_idxs)
                bin_id_r = findlast(bin_start_idx -> idx_x_max ≥ bin_start_idx, l_bin_start_idxs)
                if bin_id_r > ctx.n_threads; bin_id_r = ctx.n_threads; end # The right border of the last bin is inclusive while the others are not.

                # n_bins := 2*n_threads
                # bin i ∈ 1:n_bins contains the tets completely inside the column of bin i
                # bin j ∈ n_bins+1:2*n_bins-1 contains tets that are overlapping across bin i and i+1
                # bin k = 2*n_bins contains tets overlapping over more than 2 adjascent columns
                if bin_id_l == bin_id_r
                    l_bin_ids[simp_id] = bin_id_l
                elseif bin_id_l == bin_id_r - 1
                    l_bin_ids[simp_id] = ctx.n_threads + bin_id_l
                else
                    l_bin_ids[simp_id] = 2 * ctx.n_threads
                end
            end # for tet_id
        end # for thread_id

        #============================================#
        # Calculate number of simps in each bin,
        # print statistics
        #============================================#

        # The number of tetrahedra in each respective bin (only the first 2*n_threads bins have been filled)
        l_bin_counts = zeros(Int, 2 * ctx.n_threads + 1)

        for bin_id ∈ l_bin_ids
            l_bin_counts[bin_id] += 1
        end

        max_chars_per_bin = floor(Int, max(0, log10(maximum(l_bin_counts)))) + 1
        chars_per_line = 120
        bins_per_line = max(1, (chars_per_line ÷ (max_chars_per_bin + 1)) ÷ 8 * 8)

        println("Bin Counts:")
        for bin_id ∈ 1:length(l_bin_counts) - 1
            if bin_id != 1 && (bin_id - 1) % bins_per_line == 0
                println()
            end
            bin_chars = floor(Int, max(0, log10(l_bin_counts[bin_id]))) + 1
            print(' '^(max_chars_per_bin-bin_chars+1), l_bin_counts[bin_id])
        end
        println()

        if n_excluded_domain != 0
            println(n_excluded_domain, " ", simplex_name(ctx), " are ignored due to domain filter.")
        end

        if n_excluded_z_range != 0
            println(n_excluded_z_range, " ", simplex_name(ctx), " are ignored due to z_range filter.")
        end

        println("Done.")
        return (l_bin_ids, l_bin_counts)
    end

    function iterate(ctx, itbuf, xdmf)
        # In each iteration, process ALL simplices for ALL variables for as many timesteps as possible*
        # *possible means that the maximum memory footprint set above cannot be exceeded. This limits the number of
        # timesteps that can be stored in l_iter_output_grids.
        for iteration ∈ 1:ctx.n_iterations
            GC.gc();GC.gc();GC.gc();GC.gc();GC.gc()

            t_start = ctx.t_begin + (iteration - 1) * ctx.n_timesteps_per_iteration
            n_times = min(ctx.n_timesteps - (iteration - 1) * ctx.n_timesteps_per_iteration, ctx.n_timesteps_per_iteration)
            t_stop = t_start + n_times - 1

            #============================================#
            # Prefetch the values of all vars in all
            # timesteps for the current iteration.
            # This eliminates random accesses later on
            # and D R A S T I C A L L Y improves runtime.
            #============================================#

            for var_name ∈ ctx.in_vars_dyn, t ∈ 1:n_times
                buf = data_of(xdmf, t_start + t - 1, var_name)
                copyto!(itbuf.prefetched_vars[var_name], 1 + (t-1) * ctx.n_simplices, buf, 1, ctx.n_simplices)
            end

            # Reset the output values to 0. Only do so for the timesteps actually processed in this iteration.
            for grid ∈ values(itbuf.out_grids_dyn)
                grid[:,:, 1:n_times] .= 0
            end

            # Static vars only need to be read once. Do that in the first iteration.
            if iteration == 1
                for grid ∈ values(itbuf.out_grids_stat)
                    grid .= 0
                end
            end

            itbuf.out_sample_counts .= 0

            #============================================#
            # Binning
            #============================================#
            # First, process each of the n_threads normal
            # bins. (P-partitioning)
            #
            # Then, process the n_threads-1 bins on the
            # borders between normal bins. (B-partitioning)
            #
            # Lastly, process the 1 bin of tetrahedra
            # that overlapped more than 2 adjacent bins.
            # (R-partitioning)
            #============================================#

            # P-partitioning
            Threads.@threads for thread_id ∈ 1:ctx.n_threads
                iterate_simps!(thread_id, ctx, itbuf, t_start, t_stop, print_progress = (thread_id == ctx.n_threads ÷ 2))
            end # for thread_id

            println("Processing $(simplex_name(ctx)) on bucket borders...")

            # B-partitioning
            Threads.@threads for thread_id ∈ 1:ctx.n_threads-1
                iterate_simps!(ctx.n_threads + thread_id, ctx, itbuf, t_start, t_stop, print_progress = (thread_id == ctx.n_threads ÷ 2))
            end # for thread_id

            # R-partitioning
            iterate_simps!(2 * ctx.n_threads, ctx, itbuf, t_start, t_stop, print_progress = true)

            println("Done.")

            #============================================#
            # Write NetCDF outputs for this iteration
            #============================================#

            println("Writing outputs for timesteps $(t_start-1)-$(t_stop-1)...")

            start = [1, 1, (t_start-ctx.t_begin)+1]
            cnt = [-1, -1, n_times]

            for name ∈ ctx.out_vars_dyn
                ncwrite(itbuf.out_grids_dyn[name], ctx.out_filename, name, start=start, count=cnt)
            end

            if iteration == 1 # Only write static outputs once
                for name ∈ ctx.out_vars_stat
                    if name == ctx.stat_var_mapping["b"] # Water height correction
                        itbuf.out_grids_stat[name] .-= ctx.water_height
                    end
                    ncwrite(itbuf.out_grids_stat[name], ctx.out_filename, name, start=[1, 1], count=[-1, -1])
                end
            end

            println("Done.")
        end # for iteration
    end

    function iterate_simps!(bin_id, ctx:: RasterizationContext, itbuf:: IterationBuffers, t_start, t_stop; print_progress=false)
        n_times = t_stop - t_start + 1

        #============================================#
        # Pre-allocate memory used throughout
        # rasterization by this thread
        #============================================#

        m43_simp_points     = Array{Float64, 2}(undef, (3, ctx.n_simplex_points))
        m32_tet_aabb        = Array{Float64, 2}(undef, (3, 2))
        l_simp_points       = Array{SubArray, 1}(undef, ctx.n_simplex_points)
        l_simp_hnfs         = Array{HesseNormalForm, 1}(undef, ctx.n_simplex_points)

        v_p                 = Array{Float64, 1}(undef, 3)
        simp_sample_counts  = Array{UInt16, 2}(undef, (64, 64))

        d_start             = now()
        d_last_printed      = d_start
        print_interval      = Second(2)
        n_bin_rasterized    = 0
        n_bin_total         = itbuf.bin_counts[bin_id]

        for simp_id ∈ 1:ctx.n_simplices
            # Only process tetrahedra in thread's own bin
            if itbuf.bin_ids[simp_id] != bin_id
                continue
            end

            if print_progress
                d_now = now()

                if (d_now - d_last_printed) >= print_interval && n_bin_rasterized > 0
                    d_last_printed = d_now

                    etr  = (d_now-d_start) * (n_bin_total - n_bin_rasterized) ÷ (n_bin_rasterized)
                    hh   = floor(etr, Dates.Hour)
                    etr -= floor(hh,  Dates.Millisecond)
                    mm   = floor(etr, Dates.Minute)
                    etr -= floor(mm,  Dates.Millisecond)
                    ss   = floor(etr, Dates.Second)

                    @printf("Working on simplex %7d of %7d in timesteps %d-%d (%02.2f%% done). ETR: %02d:%02d:%02d\n",
                            n_bin_rasterized, n_bin_total,
                            t_start - 1,           # Tools like ParaView work 0-indexed. Avoid confusion by outputting 0-indexed here.
                            t_stop  - 1,           # Tools like ParaView work 0-indexed. Avoid confusion by outputting 0-indexed here.
                            n_bin_rasterized / n_bin_total * 100,
                            hh.value, mm.value, ss.value)
                end
            end # if print_progress

            # Read simplex points. m43 is only filled with 3x3 values for triangles
            for i ∈ 1:ctx.n_simplex_points
                m43_simp_points[:, i] = ctx.points[:, ctx.simplices[i, simp_id]]
            end

            for dim ∈ X:Z
                m32_tet_aabb[dim, MIN] = minimum(m43_simp_points[dim, :])
                m32_tet_aabb[dim, MAX] = maximum(m43_simp_points[dim, :])
            end

            #============================================#
            # A note on global & local indices
            #============================================#
            # In the following, idx_g_ will refer to an
            # index of the global sampling domain while
            # idx_l_ refers to an index in the simp-local
            # domain.
            #============================================#

            idx_g_min = Tuple(
                        (floor(Int32, (m32_tet_aabb[dim,MIN] - ctx.domain[dim][MIN]) / ctx.sampling_rate[dim]) + 1) for dim ∈ X:Z)

            idx_g_max = ( ceil(Int32, (m32_tet_aabb[  X,MAX] - ctx.domain[  X][MIN]) / ctx.sampling_rate[  X]),
                          ceil(Int32, (m32_tet_aabb[  Y,MAX] - ctx.domain[  Y][MIN]) / ctx.sampling_rate[  Y]),
                         floor(Int32, (m32_tet_aabb[  Z,MAX] - ctx.domain[  Z][MIN]) / ctx.sampling_rate[  Z]) + 1)

            idx_g_min = Tuple(max(1               , idx_g_min[dim]) for dim ∈ X:Z)
            idx_g_max = Tuple(min(ctx.samples[dim], idx_g_max[dim]) for dim ∈ X:Z)

            n_current_cells = idx_g_max .- idx_g_min .+ 1

            # If the current simp has a bounding box larger than anticipated, enlarge the sample counts array accordingly
            # By allocating double the size of the bounding box, future simps that are just a bit larger than the current one
            # will not trigger a resize operation.
            if size(simp_sample_counts, X) < n_current_cells[X] || size(simp_sample_counts, Y) < n_current_cells[Y]
                simp_sample_counts = Array{UInt16, 2}(undef, (n_current_cells[X] * 2, n_current_cells[Y] * 2))
            end

            # For the bounding box of the current tet, reset all sample counts to 0
            simp_sample_counts[1:n_current_cells[X], 1:n_current_cells[Y]] .= 0

            for i ∈ 1:ctx.n_simplex_points
                l_simp_points[i] = @view m43_simp_points[:,i]
            end

            calculate_hnfs!(ctx, l_simp_points, l_simp_hnfs)

            #============================================#
            # Sample each cuboid that is overlapping the
            # current simplex's AABB.
            # If its CENTER POINT lies within the simp,
            # increase the simp's sample count for the
            # cuboid's (x, y) index by 1.
            #============================================#

            for idx_g_x ∈ idx_g_min[X]:idx_g_max[X]
                v_p[X] = (idx_g_x-1) * ctx.sampling_rate[X] + ctx.domain[X][MIN]

                for idx_g_y ∈ idx_g_min[Y]:idx_g_max[Y]
                    v_p[Y] = (idx_g_y-1) * ctx.sampling_rate[Y] + ctx.domain[Y][MIN]

                    if ctx.n_simplex_points == 4 # processing tets, increment counter for intersected cells

                        for idx_g_z ∈ idx_g_min[Z]:idx_g_max[Z]
                            v_p[Z] = (idx_g_z-1) * ctx.sampling_rate[Z] + ctx.domain[Z][MIN]

                            if is_in_simplex(ctx, l_simp_hnfs, v_p) # accept cell and increment its 2D cell counter for the current tetrahedron by 1
                                idx_l_x = idx_g_x - idx_g_min[X] + 1
                                idx_l_y = idx_g_y - idx_g_min[Y] + 1
                                itbuf.out_sample_counts[idx_l_x, idx_l_y] += 1
                            end
                        end # for z
                    else # processing triangles, immediately rasterize vars
                        if is_in_simplex(ctx, l_simp_hnfs, v_p)
                            idx_l_x = idx_g_x - idx_g_min[X] + 1
                            idx_l_y = idx_g_y - idx_g_min[Y] + 1

                            # Only calculate static vars in the first timestep (read: once)
                            if ctx.n_out_vars_stat != 0 && t_start == ctx.t_begin
                                v_λ = bary_coords_2d(v_p, m43_simp_points)

                                # Calculate bathymetry height from triangle's z-coords and interpolate between them
                                b = sum(m43_simp_points[Z,1:3] .* v_λ)
                                itbuf.out_grids_stat[ctx.stat_var_mapping["b"]][idx_g_x, idx_g_y] = b
                            end

                            if ctx.tanioka
                                (∂b∂x, ∂b∂y) = get_tanioka_db(ctx, itbuf, v_p, m43_simp_points, idx_g_x, idx_g_y)
                            else
                                (∂b∂x, ∂b∂y) = (0., 0.)
                            end

                            for in_name ∈ keys(ctx.dyn_var_mapping)
                                out_name = ctx.dyn_var_mapping[in_name]
                                for t ∈ 1:n_times
                                    if in_name == "W" && ctx.tanioka
                                        tanioka_offset = (∂b∂x * itbuf.prefetched_vars["U"][simp_id, t] +  ∂b∂y * itbuf.prefetched_vars["V"][simp_id, t])
                                    else
                                        tanioka_offset = 0.
                                    end

                                    itbuf.out_grids_dyn[out_name][idx_g_x, idx_g_y, t] = itbuf.prefetched_vars[in_name][simp_id, t] - tanioka_offset
                                end
                            end

                            itbuf.out_sample_counts[idx_l_x, idx_l_y] = 1 # 2D surface is only sampled once in z-direction
                        end # if is_in_simplex()
                    end # if ctx.n_simplex_points == 4
                end # for y
            end # for x

            #============================================#
            # All cuboids overlapping the tet's domain
            # have been sampled,
            # now calculate the running average of each
            # (x, y) output grid cell with the newly gained
            # samples.
            # The weight of each tetrahedron in the
            # average is the number of samples of the grid
            # cell lying within the tet divided by the total
            # number of samples in the cell.
            #============================================#

            if ctx.z_range == z_all # Build running average of variables in volume grid
                for idx_l_y ∈ 1:n_current_cells_y, idx_l_x ∈ 1:n_current_cells_x
                    num_samples = itbuf.out_sample_counts[idx_l_x, idx_l_y]

                    if num_samples > 0
                        idx_g_x = idx_l_x + idx_g_min[X] - 1
                        idx_g_y = idx_l_y + idx_g_min[Y] - 1
                        total_samples = (itbuf.out_sample_counts[idx_g_x, idx_g_y] += num_samples)

                        for in_name ∈ keys(ctx.dyn_var_mapping)
                            out_name = ctx.dyn_var_mapping[in_name]

                            for t ∈ 1:n_times
                                avg = itbuf.out_grids_dyn[out_name][idx_g_x, idx_g_y, t]
                                itbuf.out_grids_dyn[out_name][idx_g_x, idx_g_y, t] =
                                    avg + (itbuf.prefetched_vars[in_name][simp_id, t] - avg) * (num_samples / total_samples)
                            end
                        end
                    end
                end
            end # if z_range

            n_bin_rasterized += 1
        end # for tet_id
    end

    """
        Calculates all 3 / 4 Hesse Normal Forms for the given simplex points.

    # Arguments
    - `ctx` The rasterization context
    - `l_simp_points` a list of ctx.n_simplex_points column vectors.
    - `l_simp_hnfs` a HesseNormalForm-array of size ctx.n_simplex_points. Will be populated with the calculated HNFs.
    """
    @inline function calculate_hnfs!(ctx, l_simp_points, l_simp_hnfs)
        v_n = Array{Float64, 1}(undef, 3)

        for i ∈ 1:ctx.n_simplex_points
            # Exclude point i from the tetrahedron and consider the plane defined by the other 3 points
            # If and only if (x,y,z) lies on the "positive" side (in normal direction) of all 4 planes,
            # it is inside the tetrahedron.
            # Analogous for triangles but the fourth point is placed straight above the third one which
            # allows the usage of the same point-in-simplex code afterwards.
            v_p_excl = l_simp_points[i]
            v_p1     = l_simp_points[( i    % ctx.n_simplex_points) + 1]
            v_p2     = l_simp_points[((i+1) % ctx.n_simplex_points) + 1]

            v_p3 =
                if ctx.n_simplex_points == 4; l_simp_points[((i+2) % 4) + 1]
                elseif ctx.n_simplex_points == 3; [v_p1[1], v_p1[2], v_p1[3] + 1.]
                else throw(DimensionMismatch("Only simplices with 3 or 4 points (read: triangles / tetrahedra) are supported."))
                end

            # Calculate Hesse normal form of the plane defined by p1, p2, p3
            # Normal vector
            cross3!(v_p2-v_p1, v_p3-v_p1, v_n)
            ctx.n_simplex_points == 3 && @assert (v_n[3] == 0.) "Normal vector of 2D line has non-zero z component!"
            normalize!(v_n)

            # n1x1 + n2x2 + n3x3 - d = 0; solve for d
            d = dot3(v_p1, v_n)
            dist_p_excl = dot3(v_p_excl, v_n) - d

            # Ensure that p_excl is considered to be on the "positive" side of the plane
            if dist_p_excl < 0
                v_n = -v_n
                d = -d
            end

            l_simp_hnfs[i] = HesseNormalForm((v_n[1], v_n[2], v_n[3]), d)
        end # for i
    end

    @inline function is_in_simplex(ctx, l_tet_faces, point)
        # For each tetrahedron face, filter out cells that are not in the tetrahedron
        for i ∈ 1:ctx.n_simplex_points
            face = l_tet_faces[i]
            # p_excl is ALWAYS on the "positive" side of the plane.
            # So, if p_excl and p lie on the same side, p is inside the tetrahedron
            dist_p = dot3(point, face.n0) - face.d

            # If signs of distance vectors do not match, p and p_excl are on different sides of the plane (reject)
            # As stated above, the normal vector n is chosen so that p_excl is on the POSITIVE side of the plane
            # Therefore, we only need to check if dist_p is negative to reject p.
            if dist_p < -EDGE_TOL
                return false
            end

            # The point lies approx. ON the face, reject in the tet on the lexicographically
            # more negative side of the face. This ensures that for two neighboring tetrahedra, the point is
            # only accepted for one of them.
            # Def.: lexicographical order:
            #       a ≤ b ⟺ (a.x ≤ b.x) ∨ (a.x = b.x ∧ a.y ≤ b.y) ∨ (a.x = b.x ∧ a.y = b.y ∧ a.z ≤ b.z)
            # This is a total ordering on all vectors in R^3 and will therefore always be larger for the "positive"
            # normal vector in contrast to its negated counterpart.
            if (abs(dist_p) <= EDGE_TOL)
                is_positive_side =
                    (face.n0[X]  > 0.) ||
                    (face.n0[X] == 0. && face.n0[Y]  > 0.) ||
                    (face.n0[X] == 0. && face.n0[Y] == 0. && face.n0[Z] > 0.)

                # Reject if in tetrahedron on face's lex. negative side
                if !is_positive_side
                    return false
                end
            end
        end # for i

        return true
    end

    """
    Checks if the given triangle belongs to the bathymetry mesh.
        This is done as follows:
        If the triangle is not parallel to the xy-plane (different z-coords for at least one of the points), it belongs to bathymetry.
        If the triangle is parallel to the xy-plane but not at water height, then it also belongs to bathymetry.
        Otherwise, it does not belong to bathymetry.

    # Arguments
    - `triangle_points`:    The 2D array of triangle point coords (each column consists of x, y, and z coord)
    - `water_height`:       The water height as a number
    """
    @inline function is_bathy(triangle_points, water_height = 0.)
        return any(abs.(triangle_points[Z,:] .- water_height) .> 1e-8)
    end

    @inline function get_tanioka_db(ctx, itbuf, v_p, m43_simp_points, idx_g_x, idx_g_y)
        # ∂b/∂x
        v_λ = bary_coords_2d(v_p .+ (1., 0., 0.), m43_simp_points)
        b = sum(m43_simp_points[Z,1:3] .* v_λ)
        ∂b∂x = b - itbuf.out_grids_stat[ctx.stat_var_mapping["b"]][idx_g_x, idx_g_y]

        # ∂b/∂y
        v_λ = bary_coords_2d(v_p .+ (0., 1., 0.), m43_simp_points)
        b = sum(m43_simp_points[Z,1:3] .* v_λ)
        ∂b∂y = b - itbuf.out_grids_stat[ctx.stat_var_mapping["b"]][idx_g_x, idx_g_y]

        return (∂b∂x, ∂b∂y)
    end

    """
    In-place cross product to avoid heaps of heap allocations
    """
    @inline function cross3!(a, b, ret)
        ret[1] = (a[2]*b[3] - a[3]*b[2])
        ret[2] = (a[3]*b[1] - a[1]*b[3])
        ret[3] = (a[1]*b[2] - a[2]*b[1])
    end

    """
    Explicit dot product for 3D vectors to avoid huge call stacks
    """
    @inline function dot3(a, b)
        return a[1]*b[1] + a[2]*b[2] + a[3]*b[3]
    end

    """
    Calculate barycentric coordinates of `p` in the triangle defined by `ps`.
    """
    @inline function bary_coords_2d(p :: Array{Float64, 1}, ps :: Array{Float64, 2}) :: NTuple{3, Float64}
        @assert length(p) >= 2
        @assert size(ps, 1) >= 2
        @assert size(ps, 2) >= 3

        T11 = ps[X, 1] - ps[X, 3]
        T12 = ps[X, 2] - ps[X, 3]
        T21 = ps[Y, 1] - ps[Y, 3]
        T22 = ps[Y, 2] - ps[Y, 3]

        det = T11 * T22 - T12 * T21

        λ1 = (T22*(p[X]-ps[X,3]) - T12*(p[Y]-ps[Y,3])) / det
        λ2 = (T11*(p[Y]-ps[Y,3]) - T21*(p[X]-ps[X,3])) / det
        λ3 = 1. - λ1 - λ2

        return (λ1, λ2, λ3)
    end
end # module Rasterization

max_stars_repo_name: SeisSol/SAMPLER

max_issues_repo_name: SeisSol/SAMPLER

max_forks_repo_name: SeisSol/SAMPLER


text:
#various utilities to help with differentiation.

d(::Function) = throw(MethodError(d, (Function)))
dxasy(::Function) = throw(MethodError(dbyx, (Function)))

max_stars_repo_name: interplanetary-robot/GenML

max_issues_repo_name: interplanetary-robot/GenML

max_forks_repo_name: interplanetary-robot/GenML


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/09/2014

module AbstractDynamicModelImpl

export AbstractDynamicModel

abstract AbstractDynamicModel

end

module AbstractDynamicModelInterfaces

export initializeDynamicModel,
       simulateDynamicModel

function initializeDynamicModel()
end

function simulateDynamicModel()
end

end