text:
# This file is part of the IntervalArithmetic.jl package; MIT licensed

mutable struct IntervalParameters

    precision_type::Type
    precision::Int
    rounding::Symbol
    pi::Interval{BigFloat}

    IntervalParameters() = new(BigFloat, 256, :narrow)  # leave out pi
end

const parameters = IntervalParameters()

## Precision:

"`bigequiv` creates an equivalent `BigFloat` interval to a given `AbstractFloat` interval."
function bigequiv(a::Interval{T}) where {T <: AbstractFloat}
    return guarded_setprecision(Interval, precision(T)) do
        return atomic(Interval{BigFloat}, a)
    end
end

bigequiv(x::AbstractFloat) = guarded_setprecision(precision(x)) do
    return BigFloat(x)
end

# NOTE: prevents multiple threads from calling setprecision() concurrently
const precision_lock = ReentrantLock()

guarded_setprecision(f, ::Type{Interval}, prec::Integer) = lock(precision_lock) do
    return setprecision(f, Interval, prec)
end

guarded_setprecision(f, prec::Integer) = lock(precision_lock) do
    return setprecision(f, prec)
end

setprecision(::Type{Interval}, ::Type{Float64}) = parameters.precision_type = Float64
# does not change the BigFloat precision

function setprecision(::Type{Interval}, ::Type{T}, prec::Integer) where T<:AbstractFloat
    #println("SETTING BIGFLOAT PRECISION TO $precision")
    setprecision(BigFloat, prec)

    parameters.precision_type = T
    parameters.precision = prec
    parameters.pi = atomic(Interval{BigFloat}, pi)

    prec
end

setprecision(::Type{Interval{T}}, prec) where T<:AbstractFloat =
    setprecision(Interval, T, prec)

setprecision(::Type{Interval}, prec::Integer) = setprecision(Interval, BigFloat, prec)

function setprecision(f::Function, ::Type{Interval}, prec::Integer)

    old_precision = precision(Interval)
    setprecision(Interval, prec)

    try
        return f()
    finally
        setprecision(Interval, old_precision)
    end
end

# setprecision(::Type{Interval}, precision) = setprecision(Interval, precision)
setprecision(::Type{Interval}, t::Tuple) = setprecision(Interval, t...)

precision(::Type{Interval}) = (parameters.precision_type, parameters.precision)

function Base.setrounding(f::Function, ::Type{Rational{T}},
    rounding_mode::RoundingMode) where T
    setrounding(f, float(Rational{T}), rounding_mode)
end

float(x::Interval{T}) where T = atomic( Interval{float(T)}, x)  # https://github.com/dpsanders/IntervalArithmetic.jl/issues/174

big(x::Interval) = atomic(Interval{BigFloat}, x)


text:
using RvSpectMLPlots
using Test

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

ext: jl


text:
"""
$(TYPEDSIGNATURES)

```jldoctest
julia> random_design!((Uniform(2, 3), DiscreteUniform(-1, 5), Uniform(5, 10)), 10)
10×3 Array{Float64,2}:
 2.04922  -1.0  6.44508
 2.59117   3.0  5.57183
 2.77148   5.0  8.72759
 2.25659   1.0  9.75865
 2.64968   4.0  5.72029
 2.31523   5.0  8.48192
 2.86526   5.0  5.42517
 2.44753   4.0  6.75815
 2.08284   3.0  8.94993
 2.81201  -1.0  5.89413

```
"""
function random_design!(distributions::Tuple, n::Int)
    return hcat([rand(distributions[d], n) for d in 1:size(distributions, 1)]...)
end


text:
using LinearAlgebra
using Parameters

########## InputAffineQuadraticCostEnv ##########
"""
An example of continuous-time nonlinear dynamical system introduced in
several studies on approximate dynamic programming.
# Reference
[1] K. G. Vamvoudakis and F. L. Lewis, “Online Actor-Critic Algorithm to Solve the Continuous-Time Infinite Horizon Optimal Control Problem,” Automatica, vol. 46, no. 5, pp. 878–888, 2010, doi: 10.1016/j.automatica.2010.02.018.
[2] V. Nevistic and J. A. Primbs, “Constrained Nonlinear Optimal Control: a Converse HJB Approach,” 1996.
"""
@with_kw struct InputAffineQuadraticCostEnv <: Fym
    Q = Matrix(I, 2, 2)
    R = 1
    P = [0.5 0; 0 1]
end

function f(env::InputAffineQuadraticCostEnv, x)
    x1 = x[1]
    x2 = x[2]
    f1 = -x1 + x2
    f2 = -0.5*x1 -0.5*x2*(1-(cos(2*x1)+2)^2)
    return  [f1, f2]
end

function g(env::InputAffineQuadraticCostEnv, x)
    x1 = x[1]
    g1 = 0
    g2 = cos(2*x1) +2
    return [g1, g2]
end

"""
    r(env, x, u)

Calculate running cost, i.e., V = ∫r dt.
"""
function r(env::InputAffineQuadraticCostEnv, x, u)
    x'*env.Q*x + u'*env.R*u
end

function V_optimal(env::InputAffineQuadraticCostEnv, x)
    return x'*env.P*x
end

function u_optimal(env::InputAffineQuadraticCostEnv, x)
    _g = g(env, x)
    return -_g'*env.P*x
end

function ẋ(env::InputAffineQuadraticCostEnv, x, t, u)
    ẋ = f(env, x) + g(env, x)*u
    return ẋ
end

_size = (2,)
_flatten_length = 2
_index = 1:2
size(env::InputAffineQuadraticCostEnv, x) = _size
flatten_length(env::InputAffineQuadraticCostEnv, x) = _flatten_length
index(env::InputAffineQuadraticCostEnv, x) = _index
########## end ##########


text:
function init_from_experiment(experiment, data, labels, res)
    debug(LOGGER, "[INIT] Start initialization of experiment.")
    pools = copy(experiment[:param][:initial_pools])
    split_strategy = experiment[:split_strategy]

    # model
    train_data, train_pools, _ = get_train(split_strategy, data, pools)
    debug(LOGGER, "[INIT] Initializing model '$(experiment[:model])' with $(format_observations(train_data)) observations.")
    model = instantiate(eval(experiment[:model][:type]), train_data, train_pools, experiment[:model][:param])
    initialize!(model, eval(experiment[:model][:init_strategy]))
    set_model_fitted!(res, model)
    set_adjust_K!(model, experiment[:param][:adjust_K])
    solver = experiment[:param][:solver]
    debug(LOGGER, "[INIT] Model solver for this experiment is '$(typeof(solver))'.")

    # query strategy
    debug(LOGGER, "[INIT] Initializing QS '$(experiment[:query_strategy][:type])'.")
    query_data, _, _ = get_query(split_strategy, data, pools)
    qs = QueryStrategies.initialize_qs(eval(experiment[:query_strategy][:type]), model, query_data, experiment[:query_strategy][:param])

    debug(LOGGER, "[INIT] Initializing oracle.")
    oracle = initialize_oracle(eval(experiment[:oracle]), labels)

    info(LOGGER, "[INIT] Initialization done.")
    return (model, pools, solver, qs, split_strategy, oracle)
end

function check_active_learn_args(data, labels)
    MLLabelUtils.islabelenc(labels, OneClassActiveLearning.LABEL_ENCODING) || throw(ArgumentError("Argument labels is in the wrong encoding."))
    return size(data, 2) == length(labels) || throw(ArgumentError("Number of observations ($(size(data,2))) does not equal number of labels $(length(labels))."))
end

function active_learn(experiment::Dict{Symbol, Any})
    data, labels = OneClassActiveLearning.load_data(experiment[:data_file])
    return active_learn(experiment, data, labels)
end

function active_learn(experiment::Dict{Symbol, Any}, data::Array{T, 2}, labels::Vector{Symbol}) where T <: Real
    check_active_learn_args(data, labels)
    log_experiment_info(experiment)

    res = Result(experiment)
    set_worker_info!(res)
    set_data_stats!(res, data, experiment[:split_strategy])

    model, pools, solver, qs, split_strategy, oracle = try
        init_from_experiment(experiment, data, labels, res)
    catch e
        warn(LOGGER, e)
        if isa(e, KDEException)
            res.status[:exit_code] = :kde_error
        elseif isa(e, MissingLabelTypeException)
            res.status[:exit_code] = :missing_label_type
        else
            throw(e)
        end
        return res
    end

    debug(LOGGER, "Start active learning cycle with $(experiment[:param][:num_al_iterations]) queries.")
    for i in 0:experiment[:param][:num_al_iterations]
        info(LOGGER, "Iteration $(i)")
        debug(LOGGER, "Memory consumption $(round(Int, Sys.free_memory() / 2^20)) MB / $(round(Int, Sys.total_memory() / 2^20)) MB")

        train_data, train_pools, _ = get_train(split_strategy, data, pools)
        time_set_data = @elapsed set_data!(model, train_data)
        set_pools!(model, labelmap(train_pools))

        debug(LOGGER, "[FIT] Start fitting model on $(format_observations(train_data)) observations.")
        # Workaround: redirect solver output
        stdout_orig, stderr_orig = stdout, stderr
        redirect_stdout(open("/dev/null", "w")); redirect_stderr(open("/dev/null", "w"))
        status, time_fit, mem_fit = @timed fit!(model, solver)
        redirect_stdout(stdout_orig); redirect_stderr(stderr_orig)
        debug(LOGGER, "[FIT] Fitting done ($(time_fit) s, $(format_bytes(mem_fit))).")

        @trace res.al_history i time_fit mem_fit time_set_data
        if status !== JuMP.MathOptInterface.Success && status != :Optimal
            warn(LOGGER, "Not solved to optimality. Solver status: $status.")
            res.status[:exit_code] = :solver_error
            return res
        end

        test_data, _, test_indices = get_test(split_strategy, data, pools)
        debug(LOGGER, "[TEST] Testing by predicting $(format_observations(test_data)) observations.")
        predictions = SVDD.predict(model, test_data)
        push_evaluation!(res.al_history, i, predictions, labels[test_indices])
        debug(LOGGER, "[TEST] Testing done.")

        if i < experiment[:param][:num_al_iterations]
            query_data, query_pools, query_indices = try
                get_query(split_strategy, data, pools)
            catch e
                warn(LOGGER, e)
                if isa(e, KDEException)
                     res.status[:exit_code] = :kde_error
                elseif isa(e, MissingLabelTypeException)
                    res.status[:exit_code] = :missing_label_type
                else
                    res.status[:exit_code] = :unknown_qs_error
                end
                return res
            end
            if :U ∉ query_pools
                info(LOGGER, "Aborting '$(experiment[:hash])' after $(i) iterations because no more unlabeled observations are left.")
                al_summarize!(res)
                res.status[:exit_code] = :early_abort
                return res
            end
            debug(LOGGER, "[QS] Starting query strategy on $(format_observations(query_data)) observations.")
            if i == 0
                query_id, time_qs, mem_qs = @timed get_query_object(qs, query_data, query_pools, query_indices, Int[])
            else
                query_id, time_qs, mem_qs = @timed get_query_object(qs, query_data, query_pools, query_indices, values(res.al_history, :query_history))
            end
            debug(LOGGER, "[QS] Query strategy scoring done ($(time_qs) s, $(format_bytes(mem_qs))).")
            query_label = ask_oracle(oracle, query_id)
            update_pools!(pools, query_id, query_label)
            push_query!(res.al_history, i, query_id, query_label, time_qs, mem_qs)
            debug(LOGGER, "[QS] Query(id = $(query_id), label = $(labels[query_id]))")
            debug(LOGGER, "[QS] Query strategy done.")

        end
        debug(LOGGER, "Finished iteration $(i).")
    end
    debug(LOGGER, "Finished active learning cycle.")

    al_summarize!(res)
    debug(LOGGER, "Summary done.")
    info(LOGGER, "Finished experiment '$(experiment[:hash])'.")
    res.status[:exit_code] = :success
    return res
end

function update_pools!(pools, query_id, query_label)
    pools[query_id] = query_label == :inlier ? :Lin : :Lout
    return nothing
end

"""
get_query_object(qs::QueryStrategy, data::Array{T, 2}, pools::Vector{Symbol}, global_indices::Vector{Int}, history::Vector{Int})

# Arguments
- `query_data`: Subset of the full data set
- `pools`: Labels for `query_data`
- `global_indices`: Indices of the observations in `query_data` relative to the full data set.
"""
function get_query_object(qs::QueryStrategy, query_data::Array{T, 2}, pools::Vector{Symbol}, global_indices::Vector{Int}, history::Vector{Int})::Int where T <: Real
    pool_map = labelmap(pools)
    haskey(pool_map, :U) || throw(ArgumentError("No more points that are unlabeled."))
    scores = qs_score(qs, query_data, pool_map)
    @assert length(scores) == size(query_data, 2)
    candidates = [i for i in pool_map[:U] if global_indices[i] ∉ history]
    debug(LOGGER, "[QS] Selecting from $(format_number(length(candidates))) candidates.")
    local_query_index = candidates[argmax(scores[candidates])]
    return global_indices[local_query_index]
end

function push_query!(al_history::MVHistory, i, query_id, query_label, time_qs, mem_qs)
    push!(al_history, :query_history, i, query_id)
    @trace al_history i query_label time_qs mem_qs
    return nothing
end

function push_evaluation!(al_history::MVHistory, i, predictions, labels)
    cm = ConfusionMatrix(SVDD.classify.(predictions), labels)
    push!(al_history, :cm, i, cm)
    for e in [:cohens_kappa, :matthews_corr, :f1_score, :tpr, :fpr]
        push!(al_history, e, i, eval(e)(cm))
    end
    push!(al_history, :auc, i, roc_auc(predictions, labels))
    for k in [0.01, 0.02, 0.05, 0.1, 0.2]
        auc_fpr = OneClassActiveLearning.roc_auc(predictions, labels, fpr = k)
        auc_fpr_normalized = OneClassActiveLearning.roc_auc(predictions, labels, fpr = k, normalize = true)
        push!(al_history, Symbol("auc_fpr_$(k)"), i, auc_fpr)
        push!(al_history, Symbol("auc_fpr_normalized_$(k)"), i, auc_fpr_normalized)
    end
    return nothing
end


text:
@testset "Loop utilities" begin

function looper1(N::Int, K::Int)
    x = 0
    for i in 1:N
        if (i < 20) continue end
        for j in 1:K
            x += 1
        end
    end
    return x
end

function looper2(N::Int)
    x = 0
    for i in 1 : N
        if (i > 10)
            x += 1
        else
            x += 2
        end
    end
    return x
end

function looper3(N::Int)
    x = 0
    for i in 1 : N
        for j in 1 : i
            for k in i : j
                for l in 1 : 10
                    x += 1
                end
            end
        end
    end
    return x
end

@testset "Loop detection" begin

ir = IR(typeof(looper1), Int, Int)
l1, l2 = Genify.detectloops!(ir)
@test l2.body ⊆ l1.body
@test looper1(10, 10) == IRTools.evalir(ir, nothing, 10, 10)

ir = IR(typeof(looper2), Int)
loops = Genify.detectloops!(ir)
@test length(loops) == 1
@test looper2(20) == IRTools.evalir(ir, nothing, 20)

ir = IR(typeof(looper3), Int)
l1, l2, l3, l4 = Genify.detectloops!(ir)
@test l4.body ⊆ l3.body ⊆ l2.body ⊆ l1.body
@test looper3(20) == IRTools.evalir(ir, nothing, 20)

end

@testset "Preheader insertion" begin

ir = IR(typeof(looper1), Int, Int)
preheaders = Genify.preheaders!(ir, Genify.detectloops!(ir))
loops = Genify.detectloops(ir)
@test all(IRTools.block(ir, l.header-1) in preheaders for l in loops)
@test looper1(10, 10) == IRTools.evalir(ir, nothing, 10, 10)

ir = IR(typeof(looper2), Int)
preheaders = Genify.preheaders!(ir, Genify.detectloops!(ir))
loops = Genify.detectloops(ir)
@test all(IRTools.block(ir, l.header-1) in preheaders for l in loops)
@test looper2(20) == IRTools.evalir(ir, nothing, 20)

ir = IR(typeof(looper3), Int)
preheaders = Genify.preheaders!(ir, Genify.detectloops!(ir))
loops = Genify.detectloops(ir)
@test all(IRTools.block(ir, l.header-1) in preheaders for l in loops)
@test looper3(20) == IRTools.evalir(ir, nothing, 20)

end

@testset "Loop count insertion" begin

ir = IR(typeof(looper1), Int, Int)
loops, countvars = Genify.loopcounts!(ir)
@test length(loops) == length(countvars) == 2
@test looper1(10, 10) == IRTools.evalir(ir, nothing, 10, 10)

ir = IR(typeof(looper2), Int)
loops, countvars = Genify.loopcounts!(ir)
@test length(loops) == length(countvars) == 1
@test looper2(20) == IRTools.evalir(ir, nothing, 20)

ir = IR(typeof(looper3), Int)
loops, countvars = Genify.loopcounts!(ir)
@test length(loops) == length(countvars) == 4
@test looper3(20) == IRTools.evalir(ir, nothing, 20)

end

end


text:
#!/usr/bin/env julia

using Eigenbroetler

function fibonacci(x::Number, y::Number)
    # cos(30*log(r)+13*(ϕ))**2+cos(20*log(r)-11*(ϕ))**2
    if x == 0 && y == 0
        return 0.0
    end
    r = sqrt(x * x + y * y)
    ϕ = atan(y, x)
    return cos(30.0 * log(r) + 13.0 * ϕ) ^ 2 + cos(20.0 * log(r) - 11.0 * ϕ) ^ 2
end

eb = Eigenbrot(fibonacci, 512, 512)
save("fibonacci.png", eb, ImageSetting(RealPart, LinearScale))


text:
using PerfChecker

using CompositionalNetworks
using ConstraintDomains

title = "Explore, Learn, and Compose"
path = normpath(joinpath(pathof(PerfChecker), "../perf"))
dependencies = [CompositionalNetworks, ConstraintDomains]
targets = [CompositionalNetworks]

domains = [domain([1, 2]) for i in 1:1]
pre_alloc() = foreach(_ -> explore_learn_compose(domains, allunique), 1:1)
alloc() = explore_learn_compose(domains, allunique)

alloc_check(title, dependencies, targets, pre_alloc, alloc; path)

for d in walkdir(@__DIR__), f in d[end]
    splitext(f)[2] == ".mem" && rm(joinpath(d[1], f))
end
rm(joinpath(path, "mallocs.csv"))


text:
# ---
# title: 1072. Flip Columns For Maximum Number of Equal Rows
# id: problem1072
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table
# link: <https://leetcode.com/problems/flip-columns-for-maximum-number-of-equal-rows/description/>
# hidden: true
# ---
#
# Given a `matrix` consisting of 0s and 1s, we may choose any number of columns
# in the matrix and flip **every**  cell in that column.  Flipping a cell
# changes the value of that cell from 0 to 1 or from 1 to 0.
#
# Return the maximum number of rows that have all values equal after some number
# of flips.
#
#
#
# **Example 1:**
#
#
#
#     Input: [[0,1],[1,1]]
#     Output: 1
#     Explanation: After flipping no values, 1 row has all values equal.
#
#
# **Example 2:**
#
#
#
#     Input: [[0,1],[1,0]]
#     Output: 2
#     Explanation: After flipping values in the first column, both rows have equal values.
#
#
# **Example 3:**
#
#
#
#     Input: [[0,0,0],[0,0,1],[1,1,0]]
#     Output: 2
#     Explanation: After flipping values in the first two columns, the last two rows have equal values.
#
#
#
#
# **Note:**
#
#   1. `1 <= matrix.length <= 300`
#   2. `1 <= matrix[i].length <= 300`
#   3. All `matrix[i].length`'s are equal
#   4. `matrix[i][j]` is `0` or `1`
#
#
## @lc code=start
using LeetCode

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


text:
const SIMULATION_STATUS_GAUSSIAN_VARY_RANK_SMALL                 = true
const SIMULATION_STATUS_GAUSSIAN_VARY_RANK_MEDIUM                = false
const SIMULATION_STATUS_GAUSSIAN_VARY_RANK_LARGE                 = nothing

const SIMULATION_STATUS_GAUSSIAN_VARY_MISSING_PERCENTAGE_SMALL   = false
const SIMULATION_STATUS_GAUSSIAN_VARY_MISSING_PERCENTAGE_MEDIUM  = false
const SIMULATION_STATUS_GAUSSIAN_VARY_MISSING_PERCENTAGE_LARGE   = nothing

@testset "$(format("ADMM Algorithm: Small Input Simulation [Gaussian 400 x 400]"))" begin
  if SIMULATION_STATUS_GAUSSIAN_VARY_RANK_SMALL == true
    let
      Random.seed!(65536)
      ROW = 400
      COL = 400
      for input_rank in 1:400
        @printf("small case: rank = %d\n", input_rank)
        dist = Gaussian()
        timer = TimerOutput()
        RESULTS_DIR    = GLOBAL_SIMULATION_RESULTS_DIR * "gaussian/small(400x400)(vary_rank)/" * "rank" * string(input_rank) * "/"
        LOG_FILE_NAME  = "io.log"
        DATA_FILE_NAME = "saved_variables.h5"
        LOG_FILE_PATH  = RESULTS_DIR * LOG_FILE_NAME
        DATA_FILE_PATH = RESULTS_DIR * DATA_FILE_NAME
        Base.Filesystem.mkpath(RESULTS_DIR)
        io = open(LOG_FILE_PATH, "w")
        truth_matrix        = rand([(FixedRankMatrix(Distributions.Gaussian(10, 5), rank = input_rank), ROW, COL)])
        sample_model        = provide(Sampler{BernoulliModel}(), rate = 0.8)
        input_matrix        = sample_model.draw(truth_matrix)
        manual_type_matrix  = Array{Symbol}(undef, ROW, COL)
        manual_type_matrix .= :Gaussian
        @timeit timer  "Gaussian(400x400)" * "| rank=" * string(input_rank) begin
          completed_matrix, type_tracker, tracker = complete(A               = input_matrix,
                                                             maxiter         = 200,
                                                             ρ               = 0.3,
                                                             use_autodiff    = false,
                                                             gd_iter         = 3,
                                                             debug_mode      = false,
                                                             project_rank    = nothing,
                                                             io              = io,
                                                             type_assignment = manual_type_matrix)
        end

        predicted_matrix = predict(MatrixCompletionModel(),
                                   completed_matrix = completed_matrix,
                                   type_tracker     = type_tracker)

        summary_object   = summary(MatrixCompletionModel(),
                                   predicted_matrix = predicted_matrix,
                                   truth_matrix     = truth_matrix,
                                   type_tracker     = type_tracker,
                                   tracker          = tracker)
        pickle(DATA_FILE_PATH,
               "missing_idx"      => type_tracker[:Missing],
               "completed_matrix" => completed_matrix,
               "predicted_matrix" => predicted_matrix,
               "truth_matrix"     => truth_matrix,
               "summary"          => summary_object)
        print(io, JSON.json(summary_object, 4))
        print(io, timer)
        close(io)
      end
    end
  else
    @info("Already Completed Simualtion[Gaussian vary rank][SMALL]")
  end
end

@testset "$(format("ADMM Algorithm: Medium Input Simulation [Gaussian 2000 x 2000]"))" begin
  if SIMULATION_STATUS_GAUSSIAN_VARY_RANK_MEDIUM == true
    let
      Random.seed!(65536)
      ROW = 2000
      COL = 2000
      for input_rank in union(1, collect(10:10:500))
        @printf("medium case: rank = %d\n", input_rank)
        dist = Gaussian()
        timer = TimerOutput()
        RESULTS_DIR    = GLOBAL_SIMULATION_RESULTS_DIR * "gaussian/medium(2000x2000)(vary_rank)/" * "rank" * string(input_rank) * "/"
        LOG_FILE_NAME  = "io.log"
        DATA_FILE_NAME = "saved_variables.h5"
        LOG_FILE_PATH  = RESULTS_DIR * LOG_FILE_NAME
        DATA_FILE_PATH = RESULTS_DIR * DATA_FILE_NAME
        Base.Filesystem.mkpath(RESULTS_DIR)
        io = open(LOG_FILE_PATH, "w")
        truth_matrix        = rand([(FixedRankMatrix(Distributions.Gaussian(10, 5), rank = input_rank), ROW, COL)])
        sample_model        = provide(Sampler{BernoulliModel}(), rate = 0.8)
        input_matrix        = sample_model.draw(truth_matrix)
        manual_type_matrix  = Array{Symbol}(undef, ROW, COL)
        manual_type_matrix .= :Gaussian
        @timeit timer  "Gaussian(2000x2000)" * "| rank=" * string(input_rank) begin
          completed_matrix, type_tracker, tracker = complete(A               = input_matrix,
                                                             maxiter         = 200,
                                                             ρ               = 0.3,
                                                             use_autodiff    = false,
                                                             gd_iter         = 3,
                                                             debug_mode      = false,
                                                             project_rank    = input_rank,
                                                             io              = io,
                                                             type_assignment = manual_type_matrix)
        end

        predicted_matrix = predict(MatrixCompletionModel(),
                                   completed_matrix = completed_matrix,
                                   type_tracker     = type_tracker)

        summary_object   = summary(MatrixCompletionModel(),
                                   predicted_matrix = predicted_matrix,
                                   truth_matrix     = truth_matrix,
                                   type_tracker     = type_tracker,
                                   tracker          = tracker)
        pickle(DATA_FILE_PATH,
               "missing_idx"      => type_tracker[:Missing],
               "completed_matrix" => completed_matrix,
               "predicted_matrix" => predicted_matrix,
               "truth_matrix"     => truth_matrix,
               "summary"          => summary_object)

        print(io, JSON.json(summary_object, 4))
        print(io, timer)
        close(io)
      end
    end
  else
    @info("Already Completed Simualtion[Gaussian vary rank][MEDIUM]")
  end
end


text:
## Nerlove model by OLS (fminunc)
# and fmincon, imposing that factor shares are in (0,1) and sum to 1
using DelimitedFiles, Econometrics
## prepare the data
cd(@__DIR__)
data = readdlm("../Data/nerlove.data")
data = log.(data[:,2:end])
n = size(data,1)
y = data[:,1]
x = [ones(n,1) data[:,2:5]]

## bounds and restriction
lb = [-1e6, -1e6, 0., 0., 0.0]
ub = [1e6, 1e6, 1., 1., 1.]
R = [0. 0. 1. 1. 1.]
r = [1.0]

## define the objective function and start value
obj = theta -> (y-x*theta)'*(y-x*theta)
startval = (ub+lb)/2.0

## OLS
thetahat, objvalue = fminunc(obj, startval)
println("the OLS estimates: obj. value: ", round(objvalue,digits=5))
prettyprint(thetahat)

## restricted LS
thetahat, objvalue_r, flag = fmincon(obj, startval, R, r, lb, ub) # both lower and upper bounds
println("the restricted LS estimates: obj. value: ", round(objvalue_r,digits=5))
prettyprint(thetahat)

println("Exercise: you should construct a qF test using the unrestricted and")
println("restricted objective function values")

max_stars_repo_name: mcreel/EconometricsNotes

max_issues_repo_name: mcreel/EconometricsNotes

max_forks_repo_name: mcreel/EconometricsNotes


text:
function twin01()

   k = 4
   p = poly([-0.2*ones(k); 0.39*ones(k); 0.40*ones(k)])
   z = [-0.2 k; 0.39 k; 0.4 k];

   p, PolyZeros(z)
end


text:
@NLconstraint(graph, x * y <= 1)
@NLconstraint(graph, x^2 / 3 <= 1)

@NLexpression(graph, name, exp(x))
@NLconstraint(graph, name + y == 0)

@NLconstraint(graph, sin(x)*cos(y) == 5)
@NLconstraint(graph, sin(x)*cos(y) == 5)
@NLconstraint(graph, nlconstr[i=1:2], i*x^2 == i)
@NLconstraint(graph, -0.5 <= sin(x) <= 0.5)
ψ(x) = 1, t(x,y) = 2
@NLconstraint(graph, ψ(x) + t(x,y) <= 3)

graph = Model()
@variable(graph, x)
@NLconstraint(graph, x <= sum(0 for i in []) + prod(1 for i in []))

max_stars_repo_name: BoyangSheng/Skill-studing

max_issues_repo_name: BoyangSheng/Skill-studying

max_forks_repo_name: BoyangSheng/Skill-studying


text:
# Practice implementation of code blocks 4.1 through 4.6 in McElreath's Statistical
# Rethinking, 2nd ed.

using Random
using Distributions: Uniform, rand, pdf
using Plots: histogram

# Set up
Random.seed!(17)

# Code block 4.1
distro = Uniform(-1, 1)
steps = rand(distro, 1000, 16)
positions = sum(steps, dims=2)
display(histogram(positions, bins=:40, xlabel="Position", ylabel="Frequency"))

# Code blocks 4.2 and 4.3: Gaussian through multiplication
distro = Uniform(0, 0.1)
effects = rand(distro, 1000, 12) .+ 1
growth = prod(effects, dims=2)
display(histogram(growth, bins=:40, xlabel="Growth", ylabel="Frequency"))

# Code block 4.4
distro = Uniform(0, 0.5)
effects = rand(distro, 1000, 12) .+ 1
big_growth = prod(effects, dims=2)
display(histogram(big_growth, bins=:40, xlabel="Growth", ylabel="Frequency"))

distro = Uniform(0, 0.01)
effects = rand(distro, 1000, 12) .+ 1
small_growth = prod(effects, dims=2)
display(histogram(small_growth, bins=:40, xlabel="Growth", ylabel="Frequency"))

# Code block 4.5
log_big_growth = log.(big_growth)
display(histogram(log_big_growth, bins=:40, xlabel="Log(growth)", ylabel="Frequency"))

# Code block 4.6
w = 6
n = 9
param_grid = 0.0:1/100:1.0
posterior = [pdf(Binomial(n, p), w) * pdf(Uniform(0, 1), p) for p in param_grid]
posterior = posterior / sum(posterior)


text:
push!(LOAD_PATH,"../../src/")
import UNSflow

alphadef = UNSflow.ConstDef(10. *pi/180)

hdef = UNSflow.ConstDef(0.)

udef = UNSflow.ConstDef(1.)

full_kinem = UNSflow.KinemDef(alphadef, hdef, udef)

pvt = 0.25

geometry = "NACA0012"

surf = UNSflow.TwoDSurfThick(geometry, pvt, full_kinem, ndiv=140, naterm=138)

curfield = UNSflow.TwoDFlowField()

dtstar = 0.015#UNSflow.find_tstep(alphadef)

t_tot = 5.

nsteps =Int(round(t_tot/dtstar))+1

startflag = 0

writeflag = 1

writeInterval = t_tot/10.

#delvort = delSpalart(500, 12, 1e-5)
delvort = UNSflow.delNone()

mat, surf, curfield = UNSflow.lautat(surf, curfield, nsteps, dtstar,startflag, writeflag, writeInterval, delvort)

#UNSflow.makeVortPlots2D()

#UNSflow.makeForcePlots2D()

#UNSflow.makeInfoPlots2D()

UNSflow.cleanWrite()


text:
#
# InterpolationKernels.jl --
#
# Kernel functions used for linear filtering, windowing or linear
# interpolation.
#
#------------------------------------------------------------------------------
#
# This file is part of the InterpolationKernels package licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2016-2021, Éric Thiébaut.
#

module InterpolationKernels

export
    BSpline,
    BSplinePrime,
    CardinalCubicSpline,
    CardinalCubicSplinePrime,
    CatmullRomSpline,
    CatmullRomSplinePrime,
    CubicSpline,
    CubicSplinePrime,
    Kernel,
    LanczosKernel,
    LanczosKernelPrime,
    MitchellNetravaliSpline,
    MitchellNetravaliSplinePrime,
    iscardinal,
    isnormalized

import Base: convert, adjoint, values

# List and union of floating-point types which are specially handled.
const FLOATS = (:BigFloat, :Float16, :Float32, :Float64)
const Floats = @eval $(Expr(:curly, :Union, FLOATS...))

#------------------------------------------------------------------------------
# SUPPORTS

abstract type Bound end
struct Open <: Bound end
struct Closed <: Bound end

"""

`InterpolationKernels.Support{T,S,L,R}` is the abstract type for the support of
an interpolation kernel parameterized by the floating-point type `T` used for
computations, the integer size `S` of the support and the types `L` and `R` of
the left and right bounds which can be `InterpolationKernels.Open` or
`InterpolationKernels.Closed`.

"""
abstract type Support{T<:AbstractFloat,S,L<:Bound,R<:Bound} end

Base.length(sup::Support) = length(typeof(sup))
Base.length(::Type{<:Support{T,S}}) where {T,S} = S

Base.eltype(sup::Support) = eltype(typeof(sup))
Base.eltype(::Type{<:Support{T,S}}) where {T,S} = T

"""
    InterpolationKernels.SymmetricSupport{T,S,L,R}() -> sup

yields an instance of a symmetric support parameterized by the floating-point
type `T`, the integer size `S` of the support and the types `L` and `R` of the
left and right bounds which can be `InterpolationKernels.Open` or
`InterpolationKernels.Closed`.

Depending on `L` and `R`, the support is:

    (L,R) = (Open,Open) -------> sup = (-S/2,S/2)
    (L,R) = (Closed,Open) -----> sup = [-S/2,S/2)
    (L,R) = (Open,Closed) -----> sup = (-S/2,S/2]
    (L,R) = (Closed,Closed) ---> sup = [-S/2,S/2]

"""
struct SymmetricSupport{T<:AbstractFloat,S,L,R} <: Support{T,S,L,R} end

"""
    InterpolationKernels.LeftAnchoredSupport{T,S,L,R}(a)

yields an instance of a support with lower bound `a` and parameterized by the
floating-point type `T`, the integer size `S` of the support and the types `L`
and `R` of the left and right bounds which can be `InterpolationKernels.Open`
or `InterpolationKernels.Closed`.

Depending on `L` and `R`, the support is:

    (L,R) = (Open,Open) -------> sup = (a,a+S)
    (L,R) = (Closed,Open) -----> sup = [a,a+S)
    (L,R) = (Open,Closed) -----> sup = (a,a+S]
    (L,R) = (Closed,Closed) ---> sup = [a,a+S]

"""
struct LeftAnchoredSupport{T,S,L,R} <: Support{T,S,L,R}
    a::T
end

"""
    InterpolationKernels.RightAnchoredSupport{T,S,L,R}(b)

yields an instance of a support with upper bound `b` and parameterized by the
floating-point type `T`, the integer size `S` of the support and the types `L`
and `R` of the left and right bounds which can be `InterpolationKernels.Open`
or `InterpolationKernels.Closed`.

Depending on `L` and `R`, the support is:

    (L,R) = (Open,Open) -------> sup = (b-S,b)
    (L,R) = (Closed,Open) -----> sup = [b-S,b)
    (L,R) = (Open,Closed) -----> sup = (b-S,b]
    (L,R) = (Closed,Closed) ---> sup = [b-S,b]

"""
struct RightAnchoredSupport{T,S,L,R} <: Support{T,S,L,R}
    b::T
end

"""
    InterpolationKernels.supremum(sup) -> b

yields the least upper bound `b` of the kernel support `sup`.

"""
supremum(::SymmetricSupport{T,S}) where {T,S} = T(S)/2
supremum(sup::LeftAnchoredSupport{T,S}) where {T,S} = infimum(sup) + T(S)
supremum(sup::RightAnchoredSupport) = sup.b

"""
    InterpolationKernels.infimum(sup) -> a

yields the greatest lower bound `a` of the kernel support `sup`.

"""
infimum(::SymmetricSupport{T,S}) where {T,S} = T(-S)/2
infimum(sup::LeftAnchoredSupport) = sup.a
infimum(sup::RightAnchoredSupport{T,S}) where {T,S} = supremum(sup) - T(S)

#------------------------------------------------------------------------------
# INTERPOLATION KERNELS

"""
# Interpolation Kernels

An interpolation kernel `Kernel{T,S}` is parameterized by the floating-point
type `T` of its coefficients and by the integer size `S` of its support.  For
efficiency reasons, only kernels with (small) finite size supports are
implemented in `InterpolationKernels`.

A kernel `ker` is a callable object which may be used as a function with a real
argument:

    ker(x::Real)

yields kernel value at offset `x`.  Whatever the type of `x`, `ker(x)` is
always of type `T = eltype(ker)` the floating-point type associated with `ker`.
All kernel supports are symmetric; that is `ker(x)` is zero if `abs(x) > S/2`
with `S = length(ker)` the size of the kernel support.

## Kernel floating-point type conversion

Called as a function with a real argument, a given kernel returns a value of
its associated floating-point type.  This has been chosen to have fast
interpolation methods.  Converting a kernel `ker` to use floating-point type
`T` is simply done by one of:

    T(ker)
    Kernel{T}(ker)
    convert(Kernel{T}, ker)

Beware that changing the floating-point type may lead to a loss of precision if
the kernel has numerical parameters.

## Common methods

A few common methods are specialized for any interpolation kernel `ker`:

    eltype(ker) -> T

yields the floating-point type for calculations,

    length(ker) -> S

yield the size the support of `ker` which is also the number of neighbors
involved in an interpolation by this kernel,

    values(ker)

yields a tuple of the parameters of `ker` such that an identical instance can
be built by:

    typeof(ker)(values(ker)...)

finally:

    compute_offset_and_weights(ker, x) -> off, (w1, w2, ..., wS)

yields the offset `off` and an `S`-tuple of interpolation weights to
interpolate an array at coordinate `x` (in fractional index units).

"""
abstract type Kernel{T<:AbstractFloat,S} <: Function end

Base.eltype(ker::Kernel) = eltype(typeof(ker))
Base.eltype(::Type{<:Kernel{T,S}}) where {T,S} = T

Base.length(ker::Kernel) = length(typeof(ker))
Base.length(::Type{<:Kernel{T,S}}) where {T,S} = S

"""
    isnormalized(ker)

yields whether the kernel `ker` has the partition of unity property.  That is,
the sum of the values computed by the kernel `ker` on a unit spaced grid is
equal to one.

""" isnormalized

"""
    iscardinal(ker)

yields whether the kernel `ker` is zero for any non-zero integer argument.
Cardinal kernels are directly suitable for interpolation.

""" iscardinal

"""
    InterpolationKernels.compute_offset_and_weights(ker, x) -> off, wgt

yields the index offset `off` and the weights `wgt` to interpolate with kernel
`ker` at position `x` in fractional index units.  The offset is a scalar and
the weights are an `n`-tuple with `n = length(ker)` the size of the support of
the kernel, all returned values have the same floating point type `eltype(ker)`
as the kernel.

Not taking into account boundary conditions, interpolating a vector `A` at
position `x` would then write:

    off, wgt = InterpolationKernels.compute_offset_and_weights(ker, x)
    k = Int(off)
    result = wgt[1]*A[k+1] + ... + wgt[n]*A[k+n]

Note that 1-based indexing is assumed by `compute_offset_and_weights` to
interpret the position `x` and compute the offset `off`.  If this is not the
case, the code should be:

    j1 = first(axes(A,1)) # first index in A
    off, wgt = InterpolationKernels.compute_offset_and_weights(ker, x - (j1 - 1))
    k = Int(off) + (j1 - 1)
    result = wgt[1]*A[k+1] + ... + wgt[n]*A[k+n]

where expression `x - (j1 - 1)` is assuming that the position `x` is in
fractional index for `A`, that is `x = j1` at the first entry of `A`.

!!! note
    For fast computations, this method should be specialized for specific
    kernel types.  For kernels with symmetric support, the method
    [`InterpolationKernels.compute_weights`](@ref) is called by
    `compute_offset_and_weights` to calculate the interpolation weights; for
    such kernels it is sufficient to specialize `compute_weights` instead of
    `compute_offset_and_weights`.

""" compute_offset_and_weights

@inline compute_offset_and_weights(ker::Kernel{T}, x::T) where {T} =
    compute_offset_and_weights(support(ker), ker, x)

# The following version is specialized for symmetric supports and rely on an
# optimized version of `compute_weights`.
@generated function compute_offset_and_weights(sup::SymmetricSupport{T,S},
                                               ker::Kernel{T,S},
                                               x::T) where {T,S}
    if isodd(S)
        quote
            $(Expr(:meta, :inline))
            round_x = round(x)
            wgts = compute_weights(ker, x - round_x)
            off = round_x - $((S + 1) >> 1)
            return off, wgts
        end
    else
        quote
            $(Expr(:meta, :inline))
            floor_x = floor(x)
            wgts = compute_weights(ker, x - floor_x)
            off = floor_x - $(S >> 1)
            return off, wgts
        end
    end
end

# This version is to call the generic version by default.
@inline function compute_offset_and_weights(sup::Support{T,S},
                                            ker::Kernel{T}, x::T) where {T,S}
    generic_compute_offset_and_weights(sup, ker, x)
end

# The generic version is only called if no optimized version exists.  It can
# also be used to check code implementing an optimized version.
@inline generic_compute_offset_and_weights(ker::Kernel{T}, x::T) where {T} =
    generic_compute_offset_and_weights(support(ker), ker, x)

@generated function generic_compute_offset_and_weights(sup::Support{T,S},
                                                       ker::Kernel{T,S},
                                                       x::T) where {T,S}
    wgt = ntuple(j -> Symbol("w_", j), Val(S))
    exprs = [:($(wgt[j]) = ker(u - $j)) for j in 1:S]
    quote
        $(Expr(:meta, :inline))
        off = offset(sup, x)
        u = x - off
        $(exprs...)
        return off, $(Expr(:tuple, wgt...))
    end
end

"""
    InterpolationKernels.offset(sup, x) -> off

yields the coordinate offset for computing the interpolation weights
at coordinate `x` for a kernel whose support is `sup`:

    off = floor(x - b)       if `sup` is right-open
          ceil(x - b) - 1    if `sup` is left-open

with `b = supremum(sup)` the upper bound of the support.

""" offset

offset(sup::Support{T,S,<:Bound,Open}, x::T) where {T,S} =
    floor(x - supremum(sup))

offset(sup::Support{T,S,Open,Closed}, x::T) where {T,S} =
    ceil(x - (supremum(sup) - 1))

"""
    InterpolationKernels.compute_weights(ker, t) -> wgt

computes the interpolation weights returned by
[`InterpolationKernels.compute_offset_and_weights`](@ref) for kernel `ker` with
symmetric support.  Assuming interpolation is performed at position `x`,
argument `t` is given by:

     t = x - floor(x)     if length(ker) is even, hence t ∈ [0,1)
     t = x - round(x)     if length(ker) is odd,  hence t ∈ [-1/2,+1/2]

The returned weights are then:

     wgt = ntuple(i -> ker(t + k - i), length(ker))

where `k = (length(ker) + 1) >> 1` (i.e., integer division of `length(ker)+1`
by 2).  These conventions have been adopted so that, by specializing the
`compute_weights` method, computing the `length(ker)` weights at the same time
may be done in much fewer operations than calling `ker` as a function for each
weight.

""" compute_weights

# Call generic version by default.  The generic version has a different name so
# that it can be used to check optimized versions.  Must be in-lined.
@inline compute_weights(ker::Kernel{T}, t::T) where {T} =
    generic_compute_weights(ker, t)

@generated function generic_compute_weights(ker::Kernel{T,S}, t::T) where {T,S}
    k = ((S + 1) >> 1)
    exprs = [(j < 0 ? :(ker(t + $(-j))) :
              j > 0 ? :(ker(t - $j)) : :(ker(t))) for j in 1-k:S-k]
    return quote
        $(Expr(:meta, :inline))
        return $(Expr(:tuple, exprs...))
    end
end

"""
    values(ker::InterpolationKernels.Kernel)

yields a tuple of the parameters of the interpolation kernel `ker` such that an
identical instance can be built by:

    typeof(ker)(values(ker)...)

"""
values(::Kernel) = ()

"""
    InterpolationKernels.support(ker) -> sup

yields the support of the interpolation kernel `ker`.

""" support

#------------------------------------------------------------------------------
# B-SPLINES

"""
    BSpline{S,T}()

yields a B-spline (short for *basis spline*) of order `S` that is a piecewise
polynomial function of degree `S - 1` on a support of length `S`.  The
parameter `T` is the floating-point type for computations, `T = Float64` is
assuled if this parameter is not specified.

Fr now, not all B-spline are implemented in `InterpolationKernels`, `S` must
be: `1` (for a **rectangular** B-spline), `2` (for a **linear** B-spline), `3`
(for a **quadratic** B-spline), or `4` (for a **cubic** B-spline).

If `ker` is a B-spline, then `ker'` is its derivative which can also be
directly constructed by calling [`BSplinePrime`](@ref).

!!! warning
    The derivative of B-splines of order `S ≤ 2` is not defined everywhere.  It
    is allowed to take their derivative but it (arbitrarily) yields zero where
    not defined.  Returning `NaN` would have been more correct but it has been
    considered that it would do more harm than good in practice.

""" BSpline

struct BSpline{S,T} <: Kernel{T,S} end

"""
    BSplinePrime{S,T}()

yields the derivative of a B-spline of order `S` for floating-point `T`.

See the caveats in [`BSpline`](@ref) about taking the derivative of B-splines of
order `S ≤ 2`.

""" BSplinePrime

struct BSplinePrime{S,T} <: Kernel{T,S} end

# Outer Constructors.
BSpline{S}() where {S} = BSpline{S,Float64}()
BSplinePrime{S}() where {S} = BSplinePrime{S,Float64}()

# Support is right-open for 1st order splines and open for all others.
support(::BSpline{1,T}) where {T} = SymmetricSupport{T,1,Closed,Open}()
support(::BSplinePrime{1,T}) where {T} = SymmetricSupport{T,1,Closed,Open}()
support(::BSpline{S,T}) where {S,T} = SymmetricSupport{T,S,Open,Open}()
support(::BSplinePrime{S,T}) where {S,T} = SymmetricSupport{T,S,Open,Open}()

# Only the 2 first B-splines are interpolating.
iscardinal(::Union{K,Type{K}}) where {K<:BSpline{1}} = true
iscardinal(::Union{K,Type{K}}) where {K<:BSpline{2}} = true
iscardinal(::Union{K,Type{K}}) where {K<:BSpline} = false
iscardinal(::Union{K,Type{K}}) where {K<:BSplinePrime} = false

# All B-splines are normalized, not their derivative.
isnormalized(::Union{K,Type{K}}) where {K<:BSpline} = true
isnormalized(::Union{K,Type{K}}) where {K<:BSplinePrime} = false

#
# Rectangular B-spline
# --------------------
#
@inline (::BSpline{1,T})(x::T) where {T} =
    ifelse((x < frac(T,-1,2))|(x ≥ frac(T,1,2)), zero(T), one(T))
compute_offset_and_weights(::BSpline{1,T}, x::T) where {T} =
    (round(x) - 1, (one(T),))
compute_weights(::BSpline{1,T}, t::T) where {T} = (one(T),)

(::BSplinePrime{1,T})(x::T) where {T} = zero(T)
compute_offset_and_weights(::BSplinePrime{1,T}, x::T) where {T} =
    (round(x) - 1, (zero(T),))
compute_weights(::BSplinePrime{1,T}, t::T) where {T} = (zero(T),)

#
# Linear B-splines
# ----------------
#
@inline function (::BSpline{2,T})(x::T) where {T}
    abs_x = abs(x)
    return ifelse(abs_x < 1, 1 - abs_x, zero(T))
end

compute_weights(ker::BSpline{2,T}, t::T) where {T} = (1 - t, t)

# The linear B-spline is not C¹ continuous, its left and right derivatives do
# not match at x ∈ (-1,0,1).  Setting h'(x) = NaN for x ∈ (-1,0,1) is correct
# but not useful in practice.  Taking the mean of the left and right
# derivatives yiedls h'(±1) = ∓1/2 and h'(0) = 0 but, then the support would no
# longer be semi-open interval of size 2.  For now, h'(x) = 0 for x ∈ (-1,0,1).
@inline function (::BSplinePrime{2,T})(x::T) where {T}
    if (-1 < x)&(x < 0)
        return one(T)
    elseif (0 < x)&(x < 1)
        return -one(T)
    else
        return zero(T)
    end
end

compute_weights(ker::BSplinePrime{2,T}, t::T) where {T} =
    ifelse(t > 0, (-one(T), one(T)), (zero(T), zero(T)))

#
# Quadratic B-spline
# ------------------
#
@inline function (::BSpline{3,T})(x::T) where {T}
    abs_x = abs(x)
    if abs_x ≥ frac(T,3,2)
        return zero(T)
    elseif abs_x ≤ frac(T,1,2)
        return frac(T,3,4) - abs_x*abs_x
    else
        return square(abs_x - frac(T,3,2))*frac(T,1,2)
    end
end

@inline function compute_weights(ker::BSpline{3,T}, t::T) where {T}
    #
    # Given `t = x - round(x)`, the weights are:
    #
    #     w1 = (1/8)*(1 - 2*t)^2 = (1/2)*(1/sqrt(2) - t)^2
    #     w2 = (3/4) - t^2
    #     w3 = (1/8)*(1 + 2*t)^2 = (1/2)*(1/sqrt(2) + t)^2
    #
    # which can be computed in 9 operations using the following factorized
    # expressions:
    #
    #     w1 = (1/2)*(1/sqrt(2) - t)^2
    #     w2 = (3/4) - t^2
    #     w3 = (1/2)*(1/sqrt(2) + t)^2
    #
    # or in 6 operations:
    #
    ht = t/2
    w = frac(T,1,8) + ht*t
    w1 = w - ht
    w2 = frac(T,3,4) - t*t
    w3 = w + ht
    return (w1, w2, w3)
end

@inline function (::BSplinePrime{3,T})(x::T) where {T}
    if (x ≤ frac(T,-3,2))|(x ≥ frac(T,3,2))
        return zero(T)
    elseif x < frac(T,-1,2)
        return x + frac(T,3,2)
    elseif x ≤ frac(T,1,2)
        return -2x
    else
        return x - frac(T,3,2)
    end
end

@inline function compute_weights(ker::BSplinePrime{3,T}, t::T) where {T}
    # 3 operations
    h = frac(T,1,2)
    return (t - h, -2t, t + h)
end

#
# Cubic B-spline
# --------------
#
@inline function (::BSpline{4,T})(x::T) where {T}
    abs_x = abs(x)
    if abs_x ≥ 2
        return zero(T)
    elseif abs_x ≥ 1
        return cube(2 - abs_x)*frac(T,1,6)
    else
        return (frac(T,1,2)*abs_x - 1)*abs_x*abs_x + frac(T,2,3)
    end
end

@inline function compute_weights(ker::BSpline{4,T}, t::T) where {T}
    #
    # With `t = x - floor(x)`, the weights are given by:
    #
    #     w1 = 1/6 - t/2 + t^2/2 - t^3/6
    #        = 1/6 + (t^2 - t)/2 - t^3/6
    #        = (1 - t)^3/6
    #     w2 = 2/3 - t^2 + t^3/2
    #        = 2/3 + (t/2 - 1)*t^2
    #     w3 = 1/6 + t/2 + t^2/2 - t^3/2
    #        = 1/6 + (t + t^2 - t^3)/2
    #        = 1/6 - ((t - 1)*t - 1)*t/2
    #        = 4/6 - (1 - t)^2*(t + 1)/2
    #     w4 = t^3/6
    #
    # Horner's scheme takes 6 operations per cubic polynomial, 24 operations
    # for the 4 weights.  Precomputing the powers of t, t^2 and t^3, takes 2
    # operations, then 6 operations per cubic polynomial are needed.
    #
    # Using factorizations, I manage to only use 15 operations:
    #
    h = frac(T,1,2)
    p = frac(T,2,3)
    q = frac(T,1,6)
    u = 1 - t
    u2 = u*u
    t2 = t*t
    w1 = q*u2*u
    w2 = p + (h*t - 1)*t2
    w3 = p - (h*u2)*(t + 1)
    w4 = q*t2*t
    return (w1, w2, w3, w4)
end

@inline function (::BSplinePrime{4,T})(x::T) where {T}
    if (x ≤ -2)|(x ≥ 2)
        return zero(T)
    elseif x < -1
        return frac(T,1,2)*square(x + 2)
    elseif x <  0
        return frac(T,-3,2)*x*(x + frac(T,4,3))
    elseif x <  1
        return frac(T,+3,2)*x*(x - frac(T,4,3))
    else
        return frac(T,-1,2)*square(x - 2)
    end
end

@inline function compute_weights(ker::BSplinePrime{4,T}, t::T) where {T}
    w1 = frac(T,-1,2)*(t - 1)^2
    w2 = frac(T, 3,2)*(t - frac(T,4,3))*t
    w3 = frac(T, 1,2) + (1 - frac(T,3,2)*t)*t
    w4 = frac(T, 1,2)*t^2
    return (w1, w2, w3, w4)
end

#------------------------------------------------------------------------------
# GENERIC CUBIC SPLINE

# Abstract types to implement common traits of cubic splines.
abstract type AbstractCubicSpline{T} <: Kernel{T,4} end
abstract type AbstractCubicSplinePrime{T} <: Kernel{T,4} end

support(::AbstractCubicSpline{T}) where {T} =
    SymmetricSupport{T,4,Open,Open}()
support(::AbstractCubicSplinePrime{T}) where {T} =
    SymmetricSupport{T,4,Open,Open}()

"""
    CubicSpline{T}(a, b) -> ker

yields an instance of a generic cubic spline for floating-point type `T` and
parameters `a = ker'(1)` and `b = ker(1)` the slope and the value of the
function `ker(x)` at `x = 1`.

A cubic spline kernel is at least C¹ continuous, the expression `ker'` yields a
kernel instance implementing the 1st derivative of the generic cubic spline
`ker` (see [`CubicSplinePrime`](@ref) to directly build a derivative).

Depending on the values of the parameters `a` and `b`, more specific cubic
spline kernels can be emulated:

* `CubicSpline{T}(-1/2,1/6)` yields a cubic B-spline as built by
  `BSpline{4,T}()`.

* `CubicSpline{T}(a,0)` yields a cardinal cubic spline as built by
  `CardinalCubicSpline{T}(a)`.

* `CubicSpline{T}(-1/2,0)` yields a Catmull-Rom kernel as built by
  `CatmullRomSpline{T}()`.

* `CubicSpline{T}(-b/2-c,b/6)` yields Mitchell & Netravali cubic spline as
  built by `MitchellNetravaliSpline{T}(b,c)`.

Instances of `CubicSpline` are very well optimized and, in practice, they may
be as fast or even faster than these more specialized counterparts.

""" CubicSpline

"""
    CubicSplinePrime{T}(a, b)

yields a kernel instance that is the 1st derivative of the generic cubic spline
of parameters `a` and `b` (see [`CubicSpline`](@ref)) for floating-point
type `T` (`Float64` by default).

""" CubicSplinePrime

struct CubicSpline{T} <: AbstractCubicSpline{T}
    a::T
    b::T
    c0::T
    c1::T
    c2::T
    c3::T
    c4::T
    c5::T
    c6::T
    function CubicSpline{T}(_a::Real, _b::Real) where {T}
        # NOTE: Paramaters must be signed in the expressions of the kernel
        # constants.  Since there is no loss of precision in computing these
        # constants with floating-point arithmetic even though parameters are
        # integers, we convert the parameters to type T.
        a, b = T(_a), T(_b)
        new{T}(a, b,
               1 - 2b,     # c0
               9b - a - 3, # c1
               2 + a - 6b, # c2
               -a - b,     # c3
               a + 2b,     # c4
               3b - 1,     # c5
               a + 3b)     # c6
    end
end

@inline function (ker::CubicSpline{T})(x::T) where {T}
    abs_x = abs(x)
    if abs_x ≤ 1
        return (ker.c2*abs_x + ker.c1)*abs_x^2 + ker.c0
    elseif abs_x < 2
        return (ker.c3 + ker.c4*abs_x)*(abs_x - 2)^2
    else
        return zero(T)
    end
end

@inline function compute_weights(ker::CubicSpline{T}, t::T) where {T}
    #
    # With `t = x - floor(x)`, the weights are given by:
    #
    #     w1 = (1 - t)^2 (b + (a + 2b) t)
    #        = u^2 (b + c4*t)
    #     w2 = (1 - 2b) + t^2 ((-3 - a + 9b) + (2 + a - 6b) t)
    #        = c0 + t^2 (c1 + c2*t)
    #     w3 = (1 - 2b) + (1 - t)^2 ((3b - 1) - (2 + a - 6b) t)
    #        = c0 + u^2 (c5 - c2*t)
    #     w4 = t^2 ((a + 3b) - (a + 2b) t)
    #        = t^2 (c6 - c4 t)
    #
    # where:
    #
    #     u = 1 - t
    #     c5 = 3b - 1
    #     c6 = a + 3b
    #
    # can be done in 15 operations:
    #
    u = 1 - t
    c4t = ker.c4*t
    c2t = ker.c2*t
    u2 = u*u
    t2 = t*t
    w1 = (ker.b  + c4t)*u2
    w2 = (ker.c1 + c2t)*t2 + ker.c0
    w3 = (ker.c5 - c2t)*u2 + ker.c0
    w4 = (ker.c6 - c4t)*t2
    return (w1, w2, w3, w4)
end

struct CubicSplinePrime{T} <: AbstractCubicSplinePrime{T}
    a::T
    b::T
    c1::T
    c2::T
    c3::T
    c4::T
    c5::T
    function CubicSplinePrime{T}(_a::Real, _b::Real) where {T}
        a, b = T(_a), T(_b)
        new{T}(a, b,
               -6 - 2a + 18b, # c1
               6 + 3a - 18b,  # c2
               -4a - 6b,      # c3
               3a + 6b,       # c4
               2a + 6b)       # c5
    end
end

@inline function (ker::CubicSplinePrime{T})(x::T) where {T}
    sign_x, abs_x = signabs(x)
    if abs_x ≤ 1
        return (ker.c1 + ker.c2*abs_x)*x
    elseif abs_x < 2
        return sign_x*(abs_x - 2)*(ker.c3 + ker.c4*abs_x)
    else
        return zero(T)
    end
end

@inline function compute_weights(ker::CubicSplinePrime{T}, t::T) where {T}
    #
    # With `t = x - floor(x)`, the weights are given by:
    #
    #     w1 = (a - (6b + 3a) t) (1 - t)
    #        = (a - c4 t) (1 - t)
    #     w2 = ((-6 - 2a + 18b) + (6 + 3a - 18b) t) t
    #        = (c1 + c2 t) t
    #     w3 = ((6 + 3a - 18b) t - a) (1 - t)
    #        = (c2 t - a) (1 - t)
    #     w4 = ((2a + 6b) - (3a + 6b) t) t
    #        = (c5 - c4 t) t
    #
    # can be computed in 11 operations in total:
    #
    u = 1 - t
    c4t = ker.c4*t
    c2t = ker.c2*t
    w1 = (ker.a - c4t)*u
    w2 = (ker.c1 + c2t)*t
    w3 = (c2t - ker.a)*u
    w4 = (ker.c5 - c4t)*t
    return (w1, w2, w3, w4)
end

# Outer constructors, traits and standard methods.
for K in (:CubicSpline, :CubicSplinePrime)
    @eval begin
        $K(a::Real, b::Real) = $K(promote(a, b)...)
        $K(a::Integer, b::Integer) = $K{Float64}(a, b)
        $K(a::T, b::T) where {T<:AbstractFloat} = $K{T}(a, b)
        values(ker::$K) = (ker.a, ker.b)
        iscardinal(::Union{K,Type{K}}) where {K<:$K} = false
        isnormalized(::Union{K,Type{K}}) where {K<:$K} =
            $(K === :CubicSpline)
    end
end

#------------------------------------------------------------------------------
# KEYS' CARDINAL CUBIC SPLINE

# Abstract types to implement common traits of cardinal cubic splines.
abstract type AbstractCardinalCubicSpline{T} <: AbstractCubicSpline{T} end
abstract type AbstractCardinalCubicSplinePrime{T} <: AbstractCubicSplinePrime{T} end

"""
    CardinalCubicSpline{T}(a)

yields an instance of the Keys family of cardinal cubic splines for
floating-point type `T` and parameter `a = ker'(1)` the slope of the function
`ker(x)` at `x = 1`.

These kernels are C¹ continuous piecewise normalized cardinal cubic spline
which depend on one parameter `a` and defined by:

    ker(x) = 1 + ((2 + a)*|x| - (3 + a))*x^2    if |x| ≤ 1
             a*(|x| - 1)*(|x| - 2)^2            if 1 ≤ |x| ≤ 2
             0                                  if |x| ≥ 2

The expression `ker'` yields a kernel instance which is the 1st derivative of
the Keys kernel `ker` (also see the constructor
[`CardinalCubicSplinePrime`](@ref)).

Reference:

* Keys, Robert, G., "*Cubic Convolution Interpolation for Digital Image
  Processing*", IEEE Trans. Acoustics, Speech, and Signal Processing,
  Vol. ASSP-29, No. 6, December 1981, pp. 1153-1160.

""" CardinalCubicSpline

struct CardinalCubicSpline{T} <: AbstractCardinalCubicSpline{T}
    a::T   # a
    ap2::T # a + 2
    ap3::T # a + 3
    am2::T # a - 2
    function CardinalCubicSpline{T}(_a::Real) where {T}
        a = T(_a)
        new{T}(a, a + 2, a + 3, a - 2)
    end
end

@inline function (ker::CardinalCubicSpline{T})(x::T) where {T}
    abs_x = abs(x)
    if abs_x ≥ 2
        return zero(T)
    elseif abs_x ≤ 1
        return 1 + (ker.ap2*abs_x - ker.ap3)*abs_x^2
    else
        return ker.a*(abs_x - 1)*(abs_x - 2)^2
    end
end

@inline function compute_weights(ker::CardinalCubicSpline{T}, t::T) where {T}
    #
    # Given `t = x - floor(x)`, compute:
    #
    #     w1 = a*(1 - t)^2*t
    #     w2 = (1 + t - (2 + a) t^2)*(1 - t)
    #     w3 = -t (a - (2a + 3)*t + (a + 2)*t^2)
    #     w4 = a*(1 - t)*t^2
    #
    # Noting that, with `u = 1 - t` and `v = a*t*u`(1 - t)`, the following
    # hold:
    #
    #     w1 = v*u
    #     w4 = v*t
    #     w1 + w4 = v
    #     w2 + w3 = 1 - v
    #     w2 = 1 - (1 + (2 + a)*u)*t^2
    #
    # can be done in in 11 operations:
    #
    u = 1 - t
    v = ker.a*t*u
    w = (1 + ker.ap2*u)*(t*t) # ap2 = a + 2
    w1 = v*u
    w2 = 1 - w
    w3 = w - v
    w4 = v*t
    return (w1, w2, w3, w4)
end

"""
    CardinalCubicSplinePrime{T}(a)

yields a kernel instance that is the 1st derivative of the Keys cardinal cubic
spline (see [`CardinalCubicSpline`](@ref)) for floating-point type `T` and
parameter `a`.  This derivative is given by:

    ker′(x) = (3(a + 2)*|x| - 2(a + 3))*x           if |x| ≤ 1
              (3a)*(|x| - 2)*(|x| - 4/3)*sign(x)    if 1 ≤ |x| ≤ 2
              0                                     if |x| ≥ 2

""" CardinalCubicSplinePrime

struct CardinalCubicSplinePrime{T} <: AbstractCardinalCubicSplinePrime{T}
    a ::T
    c1::T # 3(a + 2)
    c2::T # 2(a + 3)
    c3::T # 3a
    CardinalCubicSplinePrime{T}(a::Real) where {T} =
        new{T}(a, 3(a + 2), 2(a + 3), 3a)
end

@inline function (ker::CardinalCubicSplinePrime{T})(x::T) where {T}
    sign_x, abs_x = signabs(x)
    if abs_x ≥ 2
        return zero(T)
    elseif abs_x ≤ 1
        return (ker.c1*abs_x - ker.c2)*x
    else
        return (sign_x*ker.c3)*(abs_x - 2)*(abs_x - frac(T,4,3))
    end
end

# FIXME: optimize compute_weights for CardinalCubicSplinePrime

# Traits and outer constructors.
for K in (:CardinalCubicSpline, :CardinalCubicSplinePrime)
    @eval begin
        $K(a::Real) = $K{Float64}(a)
        $K(a::T) where {T<:AbstractFloat} = $K{T}(a)
        values(ker::$K) = (ker.a,)
        iscardinal(::Union{K,Type{K}}) where {K<:$K} =
            $(K === :CardinalCubicSpline)
        isnormalized(::Union{K,Type{K}}) where {K<:$K} =
            $(K === :CardinalCubicSpline)
    end
end

#------------------------------------------------------------------------------
# CATMULL-ROM INTERPOLATION KERNEL

"""
    CatmullRomSpline{T}()

yields an instance of the Catmull-Rom interpolation kernel for floating-point
type `T` which is assumed to be `Float64` if omitted.

Catmull-Rom interpolation kernel is a cardinal piecewise cubic spline defined
by:

    ker(x) = ((3/2)*|x| - (5/2))*x^2 + 1             if |x| ≤ 1
             (((5/2) - (1/2)*|x|)*|x| - 4)*|x| + 2   if 1 ≤ |x| ≤ 2
             0                                       if |x| ≥ 2

The expression `ker'` yields a kernel instance which is the 1st derivative of
the Catmull-Rom interpolation kernel `ker` (also see the constructor
[`CatmullRomSplinePrime`](@ref)).

""" CatmullRomSpline

struct CatmullRomSpline{T} <: AbstractCardinalCubicSpline{T} end

@inline function (::CatmullRomSpline{T})(x::T) where {T}
    abs_x = abs(x)
    if abs_x ≥ 2
        return zero(T)
    elseif abs_x ≤ 1
        # (((3/2)*x - (5/2))*x^2 + 1)
        return ((T(3)/2)*abs_x - T(5)/2)*abs_x*abs_x + 1
    else
        # (((5/2) - (1/2)*x)*x - 4)*x + 2
        # = (2 - x)^2 (1 - x)/2
        return (2 - abs_x)^2*(1 - abs_x)/2
    end
end

@inline function compute_weights(::CatmullRomSpline{T}, t::T) where {T}
    # 10 operations:
    u = 1 - t
    v = frac(T,-1,2)*t*u
    w1 = v*u
    w4 = v*t
    w = w4 - w1
    w2 = u - w1 + w
    w3 = t - w4 - w
    return (w1, w2, w3, w4)
end

"""
    CatmullRomSplinePrime{T}()

yields a kernel instance that is the 1st derivative of the Catmull-Rom
interpolation kernel (see [`CatmullRomSpline`](@ref)) for floating-point type
`T` which is assumed to be `Float64` if omitted.

The 1st derivative of the Catmull-Rom interpolation kernel is given by:

    ker′(x) = ((9/2)*|x| - 5)*x                      if a = |x| ≤ 1
              (5 - (3/2)*|x|)*x - 4*sign(x)          if 1 ≤ |x| ≤ 2
              0                                      if |x| ≥ 2

""" CatmullRomSplinePrime

struct CatmullRomSplinePrime{T} <: AbstractCardinalCubicSplinePrime{T} end

@inline function (::CatmullRomSplinePrime{T})(x::T) where {T}
    if x < 0
        if x ≤ -2
            return zero(T)
        elseif x ≥ -1
            return (frac(T,-9,2)*x - 5)*x
        else
            return (5 + frac(T,3,2)*x)*x + 4
        end
    else
        if x ≥ 2
            return zero(T)
        elseif x ≤ 1
            return (frac(T,9,2)*x - 5)*x
        else
            return (5 - frac(T,3,2)*x)*x - 4
        end
    end
end

@inline function compute_weights(::CatmullRomSplinePrime{T}, t::T) where {T}
    #
    # Weights are given by (18 operations):
    #
    #   w1 = -(3/2)*t^2 + 2*t - (1/2);
    #   w2 =  (9/2)*t^2 - 5*t;
    #   w3 = -(9/2)*t^2 + 4*t + (1/2);
    #   w4 =  (3/2)*t^2 - t;
    #
    #   w4 = ((3/2)*t - 1)*t;
    #   w1 = t - w4 - (1/2);
    #   w2 = 3*w4 - 2*t;
    #   w3 = t - 3*w4 + (1/2);
    #
    # can be computed in 11 operations:
    #
    w4 = frac(T,3,2)*t*t - t;
    w1 = t - w4 - frac(T,1,2);
    w2 = T(3)*w4 - 2*t;
    w3 = t - T(3)*w4 + frac(T,1,2);
    return (w1, w2, w3, w4)
end

# Traits and outer constructors.
for K in (:CatmullRomSpline, :CatmullRomSplinePrime)
    @eval begin
        $K() = $K{Float64}()
        iscardinal(::Union{K,Type{K}}) where {K<:$K} =
            $(K === :CatmullRomSpline)
        isnormalized(::Union{K,Type{K}}) where {K<:$K} =
            $(K === :CatmullRomSpline)
    end
end

#------------------------------------------------------------------------------
# MITCHELL & NETRAVALI KERNELS

"""
    MitchellNetravaliSpline{T}(b=1/3, c=1/3)

yields an instance of the Mitchell & Netravali family of kernels for
floating-point type `T` and parameters `(b,c)`.

These kernels are cubic splines which depends on 2 parameters, `b` and `c`.
Whatever the values of `(b,c)`, Mitchell & Netravali kernels are normalized,
even and C¹ continuous functions (these kernels and their first derivatives are
continuous).

Taking `b = 0` yields the family of cardinal cubic splines (see
[`CardinalCubicSpline`](@ref)) and is a sufficient and necessary condition to
have Mitchell & Netravali kernels be cardinal functions.

Using the constraint: `b + 2c = 1` yields a cubic filter with, at least,
quadratic order approximation.

Some specific values of `(b,c)` yield other well known kernels:

    (b,c) = (1,0)      --> cubic B-spline
    (b,c) = (0,-a)     --> Keys's cardinal cubic spline CardinalCubicSpline(a)
    (b,c) = (0,1/2)    --> Catmull-Rom kernel CatmullRomSpline()
    (b,c) = (b,0)      --> Duff's tensioned B-spline
    (b,c) = (6β,-α-3β) --> generic cubic spline CubicSpline(α,β)
    (b,c) = (1/3,1/3)  --> recommended by Mitchell-Netravali

The expression `ker'` yields a kernel instance which is the 1st derivative of
the Mitchell & Netravali kernel `ker` (also see the constructor
[`MitchellNetravaliSplinePrime`](@ref)).

Mitchell & Netravali family of kernels are currently instances of
[`CubicSpline`](@ref).

Reference:

* [Mitchell & Netravali, "*Reconstruction Filters in Computer Graphics*", in
  Computer Graphics, Vol. 22, Num. 4
  (1988)](http://www.cs.utexas.edu/users/fussell/courses/cs384g/lectures/mitchell/Mitchell.pdf).

""" MitchellNetravaliSpline

abstract type MitchellNetravaliSpline{T} <: AbstractCubicSpline{T} end

"""
    MitchellNetravaliSplinePrime([T=Float64,] [b=1/3, c=1/3,] B=Flat)

yields a kernel instance that is the 1st derivative of the Mitchell & Netravali
kernel (see [`MitchellNetravaliSpline`](@ref)) for floating-point type `T`,
parameters `b` and `c` and boundary conditions `B`.

""" MitchellNetravaliSplinePrime

abstract type MitchellNetravaliSplinePrime{T} <: AbstractCubicSplinePrime{T} end

for K in (:MitchellNetravaliSpline, :MitchellNetravaliSplinePrime)
    @eval begin
        $K() = $K{Float64}()
        $K{T}() where {T<:AbstractFloat} = $K{T}(one(T)/3, one(T)/3)
        $K(b::Real, c::Real) = $K{floating_point_type(b, c)}(b, c)
        function $K{T}(_b::Real, _c::Real) where {T<:AbstractFloat}
            b, c = T(_b), T(_c)
            return CubicSpline{T}(-(b/2 + c), b/6)
        end
    end
end

floating_point_type(x::Real) = floating_point_type(typeof(x))
floating_point_type(args::Real...) =
    floating_point_type(map(typeof, args)...)

floating_point_type(T::Type{<:AbstractFloat}) = T
floating_point_type(T::Type{<:Real}) = Float64
floating_point_type(types::Type{<:Real}...) =
    floating_point_type(promote_type(types...))

#------------------------------------------------------------------------------
# LANCZOS RESAMPLING KERNEL

"""
    LanczosKernel{S,T}()

yields an instance of a Lanczos re-sampling kernel of support size `S` (which
must be even) and for floating-point type `T`.

The Lanczos re-sampling kernels are even cardinal functions which tend to be
normalized for large support size.  They are defined by:

    ker(x) = S/(2*(π*x)^2)*sin(π*x)*sin(2*π*x/S)     if |x| ≤ S/2
             0                                       if |x| ≥ S/2

The expression `ker'` yields the first derivative of a Lanczos re-sampling
kernel `ker` (also see the constructor [`LanczosKernelPrime`](@ref)).

""" LanczosKernel

"""
    LanczosKernelPrime{S,T}()

yields a kernel instance that is the 1st derivative of the Lanczos re-sampling
kernel (see [`LanczosKernel`](@ref)) of support size `S` and for floating-point
type `T`.

""" LanczosKernelPrime

struct LanczosKernel{S,T} <: Kernel{T,S}
    a::T   # 1/2 support
    b::T   # a/pi^2
    c::T   # pi/a
    function LanczosKernel{S,T}() where {S,T}
        (typeof(S) == Int && S > 0 && iseven(S)) || bad_lanczos_kernel_size()
        a = T(S>>1)
        new{S,T}(a, a/T(π)^2, T(π)/a)
    end
end

struct LanczosKernelPrime{S,T} <: Kernel{T,S}
    a::T   # 1/2 support
    c::T   # pi/a
    function LanczosKernelPrime{S,T}() where {S,T}
        (typeof(S) == Int && S > 0 && iseven(S)) || bad_lanczos_kernel_size()
        a = T(S>>1)
        new{S,T}(a, T(π)/a)
    end
end

@noinline bad_lanczos_kernel_size() =
    throw(ArgumentError("Lanczos kernel size must be an even integer"))

# Outer constructors.
LanczosKernel{S}(T::Type{<:AbstractFloat} = Float64) where {S} =
    LanczosKernel{S,T}()
LanczosKernelPrime{S}(T::Type{<:AbstractFloat} = Float64) where {S} =
    LanczosKernelPrime{S,T}()

iscardinal(::Union{K,Type{K}}) where {K<:LanczosKernel} = true
iscardinal(::Union{K,Type{K}}) where {K<:LanczosKernelPrime} = false

isnormalized(::Union{K,Type{K}}) where {K<:LanczosKernel} = false
isnormalized(::Union{K,Type{K}}) where {K<:LanczosKernelPrime} = false

support(::LanczosKernel{S,T}) where {S,T} = SymmetricSupport{T,S,Open,Open}()
support(::LanczosKernelPrime{S,T}) where {S,T} = SymmetricSupport{T,S,Open,Open}()

# Expression for non-zero argument in the range (-S/2,S/2).
@inline _p(ker::LanczosKernel{S,T}, x::T) where {S,T} =
    ker.b*sin(π*x)*sin(ker.c*x)/(x*x)

# Expression for non-zero argument in the range (-S/2,S/2).
@inline function _p(ker::LanczosKernelPrime{S,T}, x::T) where {S,T}
    x1 = π*x
    s1, c1 = sin(x1), cos(x1)
    r1 = s1/x1
    x2 = ker.c*x # π*x/a
    s2, c2 = sin(x2), cos(x2)
    r2 = s2/x2
    return (c1*r2 + c2*r1 - 2*r1*r2)/x
end

(ker::LanczosKernel{S,T})(x::T) where {S,T} =
    (abs(x) ≥ ker.a ? zero(T) : x == 0 ? one(T) : _p(ker, x))

(ker::LanczosKernelPrime{S,T})(x::T) where {S,T} =
    (abs(x) ≥ ker.a ? zero(T) : x == 0 ? one(T) : _p(ker, x))

@generated function compute_weights(ker::LanczosKernel{S,T}, t::T) where {S,T}
    c = (S >> 1) # central index
    W = [Symbol("w",i) for i in 1:S] # all weights
    Expr(:block,
         Expr(:meta, :inline),
         Expr(:local, [:($w::T) for w in W]...),
         Expr(:if, :(t == zero(T)),
              Expr(:block, [:($(W[i]) = $(i == c ? 1 : 0)) for i in 1:S]...),
              Expr(:block,
                   [:($(W[i]) = _p(ker, t + $(c - i))) for i in 1:c-1]...,
                   :($(W[c]) = _p(ker, t)),
                   [:($(W[i]) = _p(ker, t - $(i - c))) for i in c+1:S]...)),
         Expr(:return, Expr(:tuple, W...)))
end

@generated function compute_weights(ker::LanczosKernelPrime{S,T},
                                    t::T) where {S,T}
    c = (S >> 1) # central index
    W = [Symbol("w",i) for i in 1:S] # all weights
    Expr(:block,
         #Expr(:meta, :inline),
         Expr(:local, [:($w::T) for w in W]...),
         [:($(W[i]) = ker(t + $(c - i))) for i in 1:c-1]...,
         :($(W[c]) = ker(t)),
         [:($(W[i]) = ker(t - $(i - c))) for i in c+1:S]...,
         Expr(:return, Expr(:tuple, W...)))
end

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

@inline frac(::Type{T}, a::Integer, b::Integer) where {T<:AbstractFloat} = (T(a)/T(b))
@inline square(x) = x*x
@inline cube(x) = x*x*x

"""
    signabs(x) -> sign(x), abs(x)

yields the sign and absolute value of `x` both with the same type as `x`.

"""
@inline signabs(x::Unsigned) = (one(x), x)
@inline signabs(x::Real) = ifelse(x > 0, (one(x), x),
                                  ifelse(x < 0, (oftype(one(x),-1), -x),
                                         (zero(x), zero(x))))

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

# Manage to call the short version of `show` for MIME output.
Base.show(io::IO, ::MIME"text/plain", ker::Kernel) = show(io, ker)

function Base.show(io::IO, ker::Kernel)
    summary(io, ker)
    print(io, '(')
    join(io, values(ker), ',')
    print(io, ')')
end

"""
    InterpolationKernels.with_eltype(T, ker)

yields an instance (resp. a type) of interpolation kernel instance (resp. type)
`ker` but with floating-point type `T`.

"""
with_eltype(::Type{T}, ::Type{K}) where {T<:AbstractFloat,K<:Kernel{T}} = K
with_eltype(::Type{T}, ker::Kernel{T}) where {T<:AbstractFloat} = ker

# The first instance below is needed to automatically do nothing when
# converting to a kernel of the same type (remmeber that convert is called to
# instanciate structure fields).
convert(::Type{K}, ker::K) where {K<:Kernel} = ker
convert(::Type{K}, ker::Kernel) where {K<:Kernel} = K(ker)

Kernel(ker::Kernel) = ker
Kernel{T}(ker::Kernel) where {T<:AbstractFloat} = with_eltype(T, ker)
Kernel{T,S}(ker::Kernel{<:Any,S}) where {T<:AbstractFloat,S} = with_eltype(T, ker)

for T in FLOATS
    @eval Base.$T(ker::Kernel) = with_eltype($T, ker)
end

"""
    InterpolationKernels.brief(ker)

yields a brief description of the kernel type or instance `ker`.

"""
brief(ker::Kernel) = brief(typeof(ker))

brief(::Type{<:BSpline{1}}) = "rectangular B-spline"
brief(::Type{<:BSplinePrime{1}}) = "derivative of rectangular B-spline"

brief(::Type{<:BSpline{2}}) = "linear B-spline"
brief(::Type{<:BSplinePrime{2}}) = "derivative of linear B-spline"

brief(::Type{<:BSpline{3}}) = "quadratic B-spline"
brief(::Type{<:BSplinePrime{3}}) = "derivative of quadratic B-spline"

brief(::Type{<:BSpline{4}}) = "cubic B-spline"
brief(::Type{<:BSplinePrime{4}}) = "derivative of cubic B-spline"

@noinline brief(::Type{<:BSpline{S}}) where {S} = "B-spline of order $S"
@noinline brief(::Type{<:BSplinePrime{S}}) where {S} =
    "derivative of B-spline of order $S"

brief(::Type{<:CubicSpline}) = "generic cubic spline"
brief(::Type{<:CubicSplinePrime}) = "derivative of generic cubic spline"

brief(::Type{<:CardinalCubicSpline}) = "cardinal cubic spline"
brief(::Type{<:CardinalCubicSplinePrime}) =
    "derivative of cardinal cubic spline"

brief(::Type{<:CatmullRomSpline}) = "Catmull & Rom cubic spline"
brief(::Type{<:CatmullRomSplinePrime}) =
    "derivative of Catmull & Rom cubic spline"

@noinline brief(::Type{<:LanczosKernel{S}}) where {S} =
    "Lanczos resampling kernel of size $S"
@noinline brief(::Type{<:LanczosKernelPrime{S}}) where {S} =
    "derivative of Lanczos resampling kernel of size $S"

# Manage to yield the derivative of (some) kernels when the notation `ker'` is
# used.
adjoint(ker::BSpline{S,T}) where {S,T} = BSplinePrime{S,T}()
adjoint(ker::CubicSpline{T}) where {T} = CubicSplinePrime{T}(ker.a, ker.b)
adjoint(ker::CardinalCubicSpline{T}) where {T} = CardinalCubicSplinePrime{T}(ker.a)
adjoint(ker::CatmullRomSpline{T}) where {T} = CatmullRomSplinePrime{T}()
adjoint(ker::LanczosKernel{S,T}) where {S,T} = LanczosKernelPrime{S,T}()

# Provide methods for all kernels.
for K in (:BSpline,             :BSplinePrime,
          :CubicSpline,         :CubicSplinePrime,
          :CardinalCubicSpline, :CardinalCubicSplinePrime,
          :CatmullRomSpline,    :CatmullRomSplinePrime,
          :LanczosKernel,       :LanczosKernelPrime)

    has_size = (K === :BSpline || K === :BSplinePrime ||
        K === :LanczosKernel || K === :LanczosKernelPrime)

    # We want that calling the kernel on a different type of real argument than
    # the floting-point type of the kernel convert the argument.
    # Unfortunately, defining:
    #
    #     (ker::$K{T})(x::Real) where {T<:AbstractFloat} = ker(convert(T,x))
    #
    # leads to ambiguities, the following is ugly but works...
    for T in FLOATS, R in (FLOATS..., :Integer)
        if R != T
            if has_size
                @eval @inline (ker::$K{S,$T})(x::$R) where {S} = ker($T(x))
            else
                @eval @inline (ker::$K{$T})(x::$R) = ker($T(x))
            end
        end
    end

    # Calling the kernel on an array.  FIXME: should be deprecated!
    if has_size
        @eval (ker::$K{S,T})(A::AbstractArray) where {S,T<:AbstractFloat} =
            map(ker, A)
    else
        @eval (ker::$K{T})(A::AbstractArray) where {T<:AbstractFloat} =
            map(ker, A)
    end

    # Change floating-point type.
    @eval $K{T}(ker::$K) where {T<:AbstractFloat} = with_eltype(T, ker)
    if has_size
        @eval begin
            $K{S,T}(ker::$K{S}) where {S,T<:AbstractFloat} = with_eltype(T, ker)
            with_eltype(::Type{T}, ::Type{<:$K{S}}) where {S,T<:AbstractFloat} = $K{S,T}
            with_eltype(::Type{T}, ker::$K{S}) where {S,T<:AbstractFloat} = $K{S,T}(values(ker)...)
            with_eltype(::Type{T}, ker::$K{T,S}) where {S,T<:AbstractFloat} = ker
        end
    else
        @eval begin
            with_eltype(::Type{T}, ::Type{<:$K}) where {T<:AbstractFloat} = $K{T}
            with_eltype(::Type{T}, ker::$K{T}) where {T<:AbstractFloat} = ker
            with_eltype(::Type{T}, ker::$K) where {T<:AbstractFloat} = $K{T}(values(ker)...)
        end
    end

    # Method `summary` prints the type of an object.
    if has_size
        @eval Base.summary(io::IO, ::$K{S,T}) where {S,T} =
            print(io, $(string(K, "{")), S, ',', T, '}')
    else
        @eval Base.summary(io::IO, ::$K{T}) where {T} =
            print(io, $(string(K, "{")), T, '}')
    end
end

end # module


text:
using Test
using MPI

MPI.Init()

comm = MPI.COMM_WORLD
MPI.Barrier(comm)
comm2 = MPI.Comm_dup(comm)
MPI.Barrier(comm2)
comm3 = MPI.Comm_dup(comm2)
MPI.Barrier(comm3)
MPI.free(comm2)
MPI.Barrier(comm3)
# Don't free comm2
MPI.free(comm3)

MPI.Finalize()
@test MPI.Finalized()


text:
# Autogenerated wrapper script for Nettle_jll for x86_64-w64-mingw32
export libhogweed, libnettle, nettle_hash

using GMP_jll
JLLWrappers.@generate_wrapper_header("Nettle")
JLLWrappers.@declare_library_product(libhogweed, "libhogweed-6.dll")
JLLWrappers.@declare_library_product(libnettle, "libnettle-8.dll")
JLLWrappers.@declare_executable_product(nettle_hash)
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll)
    JLLWrappers.@init_library_product(
        libhogweed,
        "bin\\libhogweed-6.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libnettle,
        "bin\\libnettle-8.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        nettle_hash,
        "bin\\nettle-hash.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
const config = Dict{String, Any}(
    # whether to debug AWS requests and responses
    "dbg" => false,

    # AWS credentials to connect with
    # If not set, they are taken from environment or .aws/* files
    #"id" => "",
    #"key" => "",
    #"region" => AWS.US_WEST_2,
    #"availability-zones" => ["us-west-2a", "us-west-2b"],

    # AMI and instance specifications to use for EC2 tests
    # a generic AMI to use to test ec2 operations
    # any linux AMI would do (e.g. ami-835b4efa - Ubuntu 16.04 64 bit us-west-2)
    # "ami-generic" => "ami-835b4efa",
    # the instance type to launch
    # "insttype" => "m3.medium",
    # tags to tag the instances with
    # "ownertag" => "awstest",
    # "nametag" => "awstest",
    # key to use to launch instances
    # "keyname" => "awstest",

    # AWS Account ID (required for SQS)
    #"awsAccountID" => ""
)


text:
"""
    W,P,B = vblast(H)

Find a VBLAST decomposition of `H` such that `H = pinv(W)*B*P'` or `B =
W*H*P`.  Here `P` is a permutation matrix, `B` is lower triangular with ones
on the diagonal, and `W` has orthogonal rows.

    W,P,B = vblast(H,mu)

If an SNR argument `mu` is passed in, a regularized ("MMSE") decomposition
is done, with the result that `W` will no longer have orthogonal rows and
`B` is no longer lower triangular.

# Examples
```jldoctest
julia> H= [1. 2; 3 4];W,_ = vblast(H); W
2×2 Matrix{Float64}:
 1.5  -0.5
 0.1   0.3

julia> H= BigFloat.([1.5 2; 3 4]) .+ 2im; W,_= vblast(H); W
2×2 Matrix{Complex{BigFloat}}:
      -2.0+3.0im            2.0-1.5im
 0.0779221-0.103896im  0.155844-0.103896im

```
"""
function vblast(H::Array{T,2},mu=Inf) where {T}

(K,N) = size(H);

W = zeros(T,K,K);
P = zeros(Int,N,N);
Hp = H;
k = zeros(Int,K,1);
indexH = [1:N;];
for i=1:N
    G = inv(Hp'*Hp + I/mu)*Hp'; # Do pseudoinverse
    gn = diag(G*G');
    (mx,mi) = findmin(abs.(gn))
    k[i] = indexH[mi];
    W[i,:] = G[mi,:];
    P[k[i],i]=1;
    indexH = setdiff(indexH,k[i]);
    Hp = H[:,indexH];
end
B = W*H*P  # Calculating B may not be required in some real-world scenarios
return (W,P,B)
end


text:
using BenchmarkTools

const SUITE = BenchmarkGroup()

include("bm_basis.jl")
SUITE["basis"] = basis_suite


text:
# GPU runtime library
#
# This module defines method instances that will be compiled into a target-specific image
# and will be available to the GPU compiler to call after Julia has generated code.
#
# Most functions implement, or are used to support Julia runtime functions that are expected
# by the Julia compiler to be available at run time, e.g., to dynamically allocate memory,
# box values, etc.

module Runtime

using ..GPUCompiler
using LLVM
using LLVM.Interop

## representation of a runtime method instance

struct RuntimeMethodInstance
    # either a function defined here, or a symbol to fetch a target-specific definition
    def::Union{Function,Symbol}

    return_type::Type
    types::Tuple
    name::Symbol

    # LLVM types cannot be cached, so we can't put them in the runtime method instance.
    # the actual types are constructed upon accessing them, based on a sentinel value:
    #  - nothing: construct the LLVM type based on its Julia counterparts
    #  - function: call this generator to get the type (when more control is needed)
    llvm_return_type::Union{Nothing, Function}
    llvm_types::Union{Nothing, Function}
    llvm_name::String
end

function Base.convert(::Type{LLVM.FunctionType}, rt::RuntimeMethodInstance;
                      ctx::LLVM.Context)
    types = if rt.llvm_types === nothing
        LLVMType[convert(LLVMType, typ; ctx, allow_boxed=true) for typ in rt.types]
    else
        rt.llvm_types(ctx)
    end

    return_type = if rt.llvm_return_type === nothing
        convert(LLVMType, rt.return_type; ctx, allow_boxed=true)
    else
        rt.llvm_return_type(ctx)
    end

    LLVM.FunctionType(return_type, types)
end

const methods = Dict{Symbol,RuntimeMethodInstance}()
function get(name::Symbol)
    if !haskey(methods, name)
        display(methods)
    end
    methods[name]
end

# Register a Julia function `def` as a runtime library function identified by `name`. The
# function will be compiled upon first use for argument types `types` and should return
# `return_type`. Use `Runtime.get(name)` to get a reference to this method instance.
#
# The corresponding LLVM types `llvm_types` and `llvm_return_type` will be deduced from
# their Julia counterparts. To influence that conversion, pass a callable object instead;
# this object will be evaluated at run-time and the returned value will be used instead.
#
# When generating multiple runtime functions from a single definition, make sure to specify
# different values for `name`. The LLVM function name will be deduced from that name, but
# you can always specify `llvm_name` to influence that. Never use an LLVM name that starts
# with `julia_` or the function might clash with other compiled functions.
function compile(def, return_type, types, llvm_return_type=nothing, llvm_types=nothing;
                 name=isa(def,Symbol) ? def : nameof(def), llvm_name="gpu_$name")
    meth = RuntimeMethodInstance(def,
                                 return_type, types, name,
                                 llvm_return_type, llvm_types, llvm_name)
    if haskey(methods, name)
        error("Runtime function $name has already been registered!")
    end
    methods[name] = meth

    # FIXME: if the function is a symbol, implying it will be specified by the target,
    #        we won't be able to call this function here or we'll get UndefVarErrors.
    #        work around that by generating an llvmcall stub. can we do better by
    #        using the new nonrecursive codegen to handle function lookup ourselves?
    if def isa Symbol
        args = [gensym() for typ in types]
        @eval @inline $def($(args...)) =
            ccall($("extern $llvm_name"), llvmcall, $return_type, ($(types...),), $(args...))
    end

    return
end

## exception handling

# expected functions for exception signalling
compile(:signal_exception, Nothing, ())

# expected functions for simple exception handling
compile(:report_exception, Nothing, (Ptr{Cchar},))
compile(:report_oom, Nothing, (Csize_t,))

# expected functions for verbose exception handling
compile(:report_exception_frame, Nothing, (Cint, Ptr{Cchar}, Ptr{Cchar}, Cint))
compile(:report_exception_name, Nothing, (Ptr{Cchar},))

# NOTE: no throw functions are provided here, but replaced by an LLVM pass instead
#       in order to provide some debug information without stack unwinding.

## GC

# FIXME: get rid of this and allow boxed types
T_prjlvalue(ctx) = convert(LLVMType, Any; ctx, allow_boxed=true)

function gc_pool_alloc(sz::Csize_t)
    ptr = malloc(sz)
    if ptr == C_NULL
        report_oom(sz)
        throw(OutOfMemoryError())
    end
    return unsafe_pointer_to_objref(ptr)
end

compile(gc_pool_alloc, Any, (Csize_t,), T_prjlvalue)

# expected functions for GC support
compile(:malloc, Ptr{Nothing}, (Csize_t,))

## boxing and unboxing

const tag_type = UInt
const tag_size = sizeof(tag_type)

const gc_bits = 0x3 # FIXME

# get the type tag of a type at run-time
@generated function type_tag(::Val{type_name}) where type_name
    Context() do ctx
        T_tag = convert(LLVMType, tag_type; ctx)
        T_ptag = LLVM.PointerType(T_tag)

        T_pjlvalue = convert(LLVMType, Any; ctx, allow_boxed=true)

        # create function
        llvm_f, _ = create_function(T_tag)
        mod = LLVM.parent(llvm_f)

        # this isn't really a function, but we abuse it to get the JIT to resolve the address
        typ = LLVM.Function(mod, "jl_" * String(type_name) * "_type",
                            LLVM.FunctionType(T_pjlvalue))

        # generate IR
        Builder(ctx) do builder
            entry = BasicBlock(llvm_f, "entry"; ctx)
            position!(builder, entry)

            typ_var = bitcast!(builder, typ, T_ptag)

            tag = load!(builder, typ_var)

            ret!(builder, tag)
        end

        call_function(llvm_f, tag_type)
    end
end

# we use `jl_value_ptr`, a Julia pseudo-intrinsic that can be used to box and unbox values

@inline @generated function box(val, ::Val{type_name}) where type_name
    sz = sizeof(val)
    allocsz = sz + tag_size

    # type-tags are ephemeral, so look them up at run time
    #tag = unsafe_load(convert(Ptr{tag_type}, type_name))
    tag = :( type_tag(Val(type_name)) )

    quote
        ptr = malloc($(Csize_t(allocsz)))

        # store the type tag
        ptr = convert(Ptr{tag_type}, ptr)
        Core.Intrinsics.pointerset(ptr, $tag | $gc_bits, #=index=# 1, #=align=# $tag_size)

        # store the value
        ptr = convert(Ptr{$val}, ptr+tag_size)
        Core.Intrinsics.pointerset(ptr, val, #=index=# 1, #=align=# $sz)

        unsafe_pointer_to_objref(ptr)
    end
end

@inline function unbox(obj, ::Type{T}) where T
    ptr = ccall(:jl_value_ptr, Ptr{Cvoid}, (Any,), obj)

    # load the value
    ptr = convert(Ptr{T}, ptr)
    Core.Intrinsics.pointerref(ptr, #=index=# 1, #=align=# sizeof(T))
end

# generate functions functions that exist in the Julia runtime (see julia/src/datatype.c)
for (T, t) in [Int8   => :int8,  Int16  => :int16,  Int32  => :int32,  Int64  => :int64,
               UInt8  => :uint8, UInt16 => :uint16, UInt32 => :uint32, UInt64 => :uint64,
               Bool => :bool, Float32 => :float32, Float64 => :float64]
    box_fn   = Symbol("box_$t")
    unbox_fn = Symbol("unbox_$t")
    @eval begin
        $box_fn(val)   = box($T(val), Val($(QuoteNode(t))))
        $unbox_fn(obj) = unbox(obj, $T)

        if VERSION >= v"1.8.0-DEV.600"
            compile($box_fn, Any, ($T,), T_prjlvalue; llvm_name=$"ijl_$box_fn")
            compile($unbox_fn, $T, (Any,); llvm_name=$"ijl_$unbox_fn")
        else
            compile($box_fn, Any, ($T,), T_prjlvalue; llvm_name=$"jl_$box_fn")
            compile($unbox_fn, $T, (Any,); llvm_name=$"jl_$unbox_fn")
        end
    end
end

end


text:
struct NcFiles
    u::Union{NcFile,Nothing}        # zonal velocity
    v::Union{NcFile,Nothing}        # meridional velocity
    η::Union{NcFile,Nothing}        # sea surface height anomaly
    sst::Union{NcFile,Nothing}      # tracer / sea surface temperature
    q::Union{NcFile,Nothing}        # potentival vorticity
    ζ::Union{NcFile,Nothing}        # relative vorticity
    du::Union{NcFile,Nothing}       # tendency of u [m^2/s^2]
    dv::Union{NcFile,Nothing}       # tendency of v [m^2/s^2]
    dη::Union{NcFile,Nothing}       # tendency of η [m^2/s]
    iout::Array{Integer,1}          # output index, stored in array for mutability
end

"""Generator function for "empty" NcFiles struct."""
NcFiles(x::Nothing) = NcFiles(x,x,x,x,x,x,x,x,x,[0])

"""Generator function for NcFiles struct, creating the underlying netCDF files."""
function NcFiles(feedback::Feedback,S::ModelSetup)

    if S.parameters.output

        @unpack output_vars,compression_level = S.parameters
        P = S.parameters
        @unpack x_u,x_v,x_T,x_q = S.grid
        @unpack y_u,y_v,y_T,y_q = S.grid

        @unpack run_id,runpath = feedback

        ncu = if "u" in output_vars nc_create(x_u,y_u,"u",runpath,"m/s","zonal velocity",P) else nothing end
        ncv = if "v" in output_vars nc_create(x_v,y_v,"v",runpath,"m/s","meridional velocity",P) else nothing end
        ncη = if "η" in output_vars nc_create(x_T,y_T,"eta",runpath,"m","sea surface height",P) else nothing end
        ncsst = if "sst" in output_vars nc_create(x_T,y_T,"sst",runpath,"1","sea surface temperature",P) else nothing end
        ncq = if "q" in output_vars nc_create(x_q,y_q,"q",runpath,"1/(ms)","potential vorticity",P) else nothing end
        ncζ = if "ζ" in output_vars nc_create(x_q,y_q,"relvort",runpath,"1","relative vorticity",P) else nothing end
        ncdu = if "du" in output_vars nc_create(x_u,y_u,"du",runpath,"m^2/s^2","zonal velocity tendency",P) else nothing end
        ncdv = if "dv" in output_vars nc_create(x_v,y_v,"dv",runpath,"m^2/s^2","meridional velocity tendency",P) else nothing end
        ncdη = if "dη" in output_vars nc_create(x_T,y_T,"deta",runpath,"m^2/s","sea surface height tendency",P) else nothing end

        for nc in (ncu,ncv,ncη,ncsst,ncq,ncζ,ncdu,ncdv,ncdη)
            if nc != nothing
                NetCDF.putatt(nc,"t",Dict("units"=>"s","long_name"=>"time"))
                NetCDF.putatt(nc,"x",Dict("units"=>"m","long_name"=>"zonal coordinate"))
                NetCDF.putatt(nc,"y",Dict("units"=>"m","long_name"=>"meridional coordinate"))
            end
        end

        return NcFiles(ncu,ncv,ncη,ncsst,ncq,ncζ,ncdu,ncdv,ncdη,[0])
    else
        return NcFiles(nothing)
    end
end

"""Creates a netCDF file based on grid vectors x,y the variable name, its path, unit and long_name."""
function nc_create( x::Array{T,1},
                    y::Array{T,1},
                    name::String,
                    path::String,
                    unit::String,
                    long_name::String,
                    P::Parameter) where {T<:Real}

    @unpack compression_level = P

    xdim = NcDim("x",length(x),values=x)
    ydim = NcDim("y",length(y),values=y)
    tdim = NcDim("t",0,unlimited=true)

    var = NcVar(name,[xdim,ydim,tdim],t=Float32,compress=compression_level)
    tvar = NcVar("t",tdim,t=Int32)

    nc = NetCDF.create(joinpath(path,name*".nc"),[var,tvar],mode=NC_NETCDF4)
    # add missing_value although irrelevant for ncview compatibility
    NetCDF.putatt(nc,name,Dict("units"=>unit,"long_name"=>long_name,"missing_value"=>-999999f0))
    return nc
end

function output_nc!(i::Int,
                    ncs::NcFiles,
                    Prog::PrognosticVars,
                    Diag::DiagnosticVars,
                    S::ModelSetup)

    @unpack nout = S.grid
    @unpack output = S.parameters

    if i % nout == 0 && output

        iout = ncs.iout[1]+1        # unpack and increase
        ncs.iout[1] = iout          # pack

        @unpack halo,haloη,halosstx,halossty = S.grid
        @unpack f_q,ep,dtint = S.grid
        @unpack scale,scale_inv,scale_sst = S.constants

        # CUT OFF HALOS
        # As output is before copyto!(u,u0), take u0,v0,η0
        # Tendencies calculate from the last time step, du = u_n+1-u_n etc
        # WRITING THE VARIABLES
        if ncs.u != nothing
            @views u = Float32.(scale_inv*Diag.RungeKutta.u0[halo+1:end-halo,halo+1:end-halo])
            NetCDF.putvar(ncs.u,"u",u,start=[1,1,iout],count=[-1,-1,1])
        end
        if ncs.v != nothing
            @views v = Float32.(scale_inv*Diag.RungeKutta.v0[halo+1:end-halo,halo+1:end-halo])
            NetCDF.putvar(ncs.v,"v",v,start=[1,1,iout],count=[-1,-1,1])
        end
        if ncs.η != nothing
            @views η = Float32.(Diag.RungeKutta.η0[haloη+1:end-haloη,haloη+1:end-haloη])
            NetCDF.putvar(ncs.η,"eta",η,start=[1,1,iout],count=[-1,-1,1])
        end
        if ncs.sst != nothing
            @views sst = Float32.(Prog.sst[halosstx+1:end-halosstx,halossty+1:end-halossty]/scale_sst)
            NetCDF.putvar(ncs.sst,"sst",sst,start=[1,1,iout],count=[-1,-1,1])
        end
        if ncs.q != nothing
            @views q = Float32.(scale_inv*Diag.Vorticity.q[haloη+1:end-haloη,haloη+1:end-haloη])
            NetCDF.putvar(ncs.q,"q",q,start=[1,1,iout],count=[-1,-1,1])
        end
        if ncs.ζ != nothing
            @unpack dvdx,dudy = Diag.Vorticity
            @unpack f_q = S.grid
            @views ζ = Float32.((dvdx[2:end-1,2:end-1]-dudy[2+ep:end-1,2:end-1])./abs.(f_q))
            NetCDF.putvar(ncs.ζ,"relvort",ζ,start=[1,1,iout],count=[-1,-1,1])
        end
        if ncs.du != nothing
            @views u = Float32.(scale_inv*Diag.RungeKutta.u0[halo+1:end-halo,halo+1:end-halo])
            @views du = u-Float32.(scale_inv*Prog.u[halo+1:end-halo,halo+1:end-halo])
            NetCDF.putvar(ncs.du,"du",du,start=[1,1,iout],count=[-1,-1,1])
        end
        if ncs.dv != nothing
            @views v = Float32.(scale_inv*Diag.RungeKutta.v0[halo+1:end-halo,halo+1:end-halo])
            @views dv = v-Float32.(scale_inv*Prog.v[halo+1:end-halo,halo+1:end-halo])
            NetCDF.putvar(ncs.dv,"dv",dv,start=[1,1,iout],count=[-1,-1,1])
        end
        if ncs.dη != nothing
            @views η = Float32.(Diag.RungeKutta.η0[haloη+1:end-haloη,haloη+1:end-haloη])
            @views dη = η-Float32.(Prog.η[haloη+1:end-haloη,haloη+1:end-haloη])
            NetCDF.putvar(ncs.dη,"deta",dη,start=[1,1,iout],count=[-1,-1,1])
        end

        # WRITING THE TIME
        for nc in (ncs.u,ncs.v,ncs.η,ncs.sst,ncs.q,ncs.ζ,ncs.du,ncs.dv,ncs.dη)
            if nc != nothing
                NetCDF.putvar(nc,"t",Int64[i*dtint],start=[iout])
                NetCDF.sync(nc)         # sync to view netcdf while model is still running
            end
        end
    end
end

"""Closes netCDF and progress.txt files."""
function output_close!(ncs::NcFiles,feedback::Feedback,S::ModelSetup)

    @unpack output = S.parameters

    if output
        println("All data stored.")
        write(feedback.progress_txt,"All data stored.")
        close(feedback.progress_txt)
    end
end

"""Finds the first gap in a list of integers."""
function gap(a::Array{Int,1})
    try
        return minimum([i for i in minimum(a):maximum(a) if ~(i in a)])
    catch
        return maximum(a)+1
    end
end

"""Checks output folders to determine a 4-digit run id number."""
function get_run_id_path(S::ModelSetup)

    @unpack output,outpath,get_id_mode = S.parameters

    if output
        runlist = filter(x->startswith(x,"run"),readdir(outpath))
        existing_runs = [parse(Int,id[4:end]) for id in runlist]
        if length(existing_runs) == 0           # if no runfolder exists yet
            runpath = joinpath(outpath,"run0000")
            mkdir(runpath)
            return 0,runpath
        else                                    # create next folder
            if get_id_mode == "fill"  # find the smallest gap in runfolders
                run_id = gap(existing_runs)
                runpath = joinpath(outpath,"run"*@sprintf("%04d",run_id))
                mkdir(runpath)

            elseif get_id_mode == "specific" # specify the run_id as input argument
                @unpack run_id = S.parameters
                runpath = joinpath(outpath,"run"*@sprintf("%04d",run_id))
                try # create folder if not existent
                    mkdir(runpath)
                catch # else rm folder and create new one
                    rm(runpath,recursive=true)
                    mkdir(runpath)
                end

            elseif get_id_mode == "continue" # find largest folder and count one up
                run_id = maximum(existing_runs)+1
                runpath = joinpath(outpath,"run"*@sprintf("%04d",run_id))
                mkdir(runpath)
            else
                throw(error("Order '$get_id_mode' is not valid for get_run_id_path(), choose continue or fill."))
            end
            return run_id,runpath
        end
    else
        return 0,"no runpath"
    end
end