text:
using DataFrames
using MLJ
using MagneticLocSuchowiak
using Plots

function getdata()
    DATA_PATH = "/home/sebastian/Desktop/PracaInz/data"
    CURVES = DATA_PATH * "/curves"
    LINES = DATA_PATH * "/lines"
    train_dirs = vcat(readdir(LINES, join = true), readdir(CURVES, join = true))
    TESTS = DATA_PATH * "/tests"
    test_files = readdir(TESTS, join = true)
    train_files = collect(Iterators.flatten(readdir.(train_dirs, join = true)))
    #train_files = [x for x in train_files if !contains(last(split(x, "\\")), "r")]

    ["magnetometr", "accelerometr", "orientation"]

    ALL_COLUMNS = [
        :timestep,
        :magnetometr_x,
        :magnetometr_y,
        :magnetometr_z,
        :acceleromter_x,
        :acceleromter_y,
        :acceleromter_z,
        :orientation_x,
        :orientation_y,
        :orientation_z,
        :lon,
        :lat,
        :path_id,
        :path_sample,
    ]

    TRAIN_COLUMNS = [
        :timestep,
        :magnetometr_x,
        :magnetometr_y,
        :magnetometr_z,
        :lon,
        :lat,
        :path_id,
        :path_sample,
    ]
    TEST_COLUMNS =
        [:timestep, :magnetometr_x, :magnetometr_y, :magnetometr_z,  :lon, :lat, :path_id]
    traindf = readdata(train_files, columns_to_get = TRAIN_COLUMNS, samplestocut = 3)
    testdf = readdata(test_files, columns_to_get = TEST_COLUMNS)

    coordinates = vcat(traindf[:, [:lon, :lat]], testdf[:, [:lon, :lat]])
    coordinates = coerce(coordinates, :lon => Continuous, :lat => Continuous)
    coordinatesstand = fit!(machine(Standardizer(), coordinates))
    traindf[!, [:lon, :lat]] = MLJ.transform(coordinatesstand, traindf[!, [:lon, :lat]])
    testdf[!, [:lon, :lat]] = MLJ.transform(coordinatesstand, testdf[!, [:lon, :lat]])

    # magcols = [:magnetometr_x, :magnetometr_y, :magnetometr_z]
    # mag_test = testdf[:, magcols]
    # mag_test = coerce(mag_test, :magnetometr_x => Continuous, :magnetometr_y => Continuous, :magnetometr_z => Continuous)
    # mag_test_stand = fit!(machine(Standardizer(), mag_test))
    # testdf[!, magcols] = MLJ.transform(mag_test_stand, testdf[!, magcols])
    #
    # mag_train =traindf[:, magcols]
    # mag_train = coerce(mag_train, :magnetometr_x => Continuous, :magnetometr_y => Continuous, :magnetometr_z => Continuous)
    # mag_train_stand = fit!(machine(Standardizer(), mag_train))
    # traindf[!, magcols] = MLJ.transform(mag_train_stand, traindf[!, magcols])

    return traindf, testdf, coordinatesstand
end

function msummary(lat, lon, traindf, testdf, coordinatesstand)
    println("Hyperparameters lat model: ", params(fitted_params(lat).best_model))
    println("Hyperparameters lon model: ", params(fitted_params(lon).best_model))
    println("Test set evaluation: ", testmachines(lon, lat, testdf, coordinatesstand))
    println("Train set evaluation: ", testmachines(lon, lat, traindf, coordinatesstand))
end

traindf, testdf, standardizer = getdata()
testdf = traindf[traindf[:, :path_sample].=="05", :]
traindf = traindf[.!(traindf[:, :path_sample].=="05"), :]

lon, lat = fitknn(traindf, klower = 50, kupper = 51)
lon, lat = fitforest(
    traindf,
    ntrees = [100],
    maxdepths = [50],
    minleafs = [50, 100],
)
lon, lat = fitnusvr(traindf, [0.1, 0.2, 0.3], [0.5, 1])

msummary(lat, lon, traindf, testdf, standardizer)

svmnuerror(traindf, standardizer)
svmgammaerror(traindf, standardizer)
knnkerror(traindf, standardizer)

ids = string.("tt", vcat(string.("0", 1:9), ["10", "11"]))
println("\n++++++++++++")
for id in ids
    path = testdf[testdf[:, :path_id].==id, :]
    testresult = testmachines(lon, lat, path, standardizer)
    print(id*" "*string(testresult["mean"])*"\n")
end
println("++++++++++++")

pathid = "tt03"
plotpathcomparison(pathid, testdf, lat, lon)
pathid = "l1n"
plottrain(pathid, traindf, lon, lat)

println("\n++++++++++++")
for id in ids
    path = testdf[testdf[:, :path_id].==id, :]

    realtestdf = MLJ.copy(path)
    realtestdf[!, [:lon, :lat]] =
        MLJ.inverse_transform(standardizer, realtestdf[!, [:lon, :lat]])
    meanlat = mean(realtestdf.lat)
    meanlon = mean(realtestdf.lon)
    naivepred = [(meanlon, meanlat) for i = 1:nrow(realtestdf)]
    ytest = collect(zip(realtestdf[:, :lon], realtestdf[:, :lat]))
    testresult = MagneticLocSuchowiak.evaluateresults(ytest, naivepred)

    print(id*" "*string(testresult["mean"])*"\n")
end
println("++++++++++++")

function foresterror()
    train_mean = []
    test_mean = []
    for i in 1:300
        lon, lat = fitcustommodel(traindf, RandomForestRegressor(n_trees=i))

        results_train = testmachines(lon, lat, traindf, standardizer)
        push!(train_mean, results_train["mean"])

        results_test = testmachines(lon, lat, testdf, standardizer)
        push!(test_mean, results_test["mean"])
    end
    plot(
        1:300,
        [train_mean, test_mean],
        label = ["test" "train"],
        title = "Mean Distance Error",
        ylabel = "s=[m]",
        xlabel = "n_trees",
    )
end

function svmnuerror(traindf, standardizer)
    train_mean = []
    test_mean = []
    nus = 0.1:0.1:0.9
    for i in nus
        lon, lat = fitcustommodel(traindf, NuSVR(gamma=0.1, nu=i))

        results_train = testmachines(lon, lat, traindf, standardizer)
        push!(train_mean, results_train["mean"])

        results_test = testmachines(lon, lat, testdf, standardizer)
        push!(test_mean, results_test["mean"])
    end
    display(plot(nus, [train_mean, test_mean], label=["train" "test"], title="Błąd w zależności od nu", ylabel="s=[m]", xlabel="nu"))
    savefig("svm_nu.png")
end

function svmgammaerror(traindf, standardizer)
    train_mean = []
    test_mean = []
    gammas = [0.1, 0.5, 1, 3, 5, 7, 10, 15]
    for i in gammas
        lon, lat = fitcustommodel(traindf, NuSVR(nu=0.1, gamma=i))

        results_train = testmachines(lon, lat, traindf, standardizer)
        push!(train_mean, results_train["mean"])

        results_test = testmachines(lon, lat, testdf, standardizer)
        push!(test_mean, results_test["mean"])
    end
    display(plot(gammas, [train_mean, test_mean], label=["train" "test"], title="Bląd w zależności od gamma", ylabel="s=[m]", xlabel="gamma"))
    savefig("svm_gamma.png")
end

function knnkerror(traindf, standardizer)
    train_mean = []
    test_mean = []
    ks = vcat(1:20, 21:10:200)
    for i in ks
        lon, lat = fitcustommodel(traindf, KNNRegressor(K=i))

        results_train = testmachines(lon, lat, traindf, standardizer)
        push!(train_mean, results_train["mean"])

        results_test = testmachines(lon, lat, testdf, standardizer)
        push!(test_mean, results_test["mean"])
    end
    plot(ks, [train_mean, test_mean], label=["train" "test"], title="Błąd w zależności od k", ylabel="s=[m]", xlabel="k")
    savefig("knn_k.png")
end

function plotpathcomparison(testpathfile, testdf, lat, lon)
    pathdf = testdf[testdf[:, :path_id].==testpathfile, :]
    X, y = MagneticLocSuchowiak.preparedata(pathdf)
    pred1 = MLJ.predict(lon, X)
    pred2 = MLJ.predict(lat, X)

    lonplot = plot(
        pathdf[:, :timestep],
        hcat(y[:, :lon], pred1),
        w = 3,
        title = "Longitude",
        label = ["Actual" "Predicted"],
        xlabel = "[t] = s",
        ylabel = "Standarized Logitude",
    )

    latplot = plot(
        pathdf[:, :timestep],
        hcat(y[:, :lat], pred2),
        w = 3,
        title = "Latitude",
        label = ["Actual" "Predicted"],
        xlabel = "[t] = s",
        ylabel = "Standarized Latitude",
    )

    latlonplot = plot(
        hcat(y[:, :lon], pred1),
        hcat(y[:, :lat], pred2),
        w = 3,
        title = "Path",
        label = ["Actual" "Predicted"],
        xlabel = "Standarized Logitude",
        ylabel = "Standarized Latitude",
    )

    display(lonplot)
    display(latplot)
    display(latlonplot)

    display(plot(lonplot, latplot, latlonplot, layout=grid(3,1), size=(600, 900)))
    return lonplot, latplot, latlonplot
end

function odstajace(traindf)
    sample =
        traindf[.&(traindf[:, :path_id] .== "l6n", traindf[:, :path_sample] .== "05"), :]
    plot([sample.timestep], [sample.magnetometr_x], w = 3, label = "składowa x natężenia")
    plot!(
        [sample.timestep[1:2]],
        [sample.magnetometr_x[1:2]],
        seriestype = :scatter,
        title = "",
        label = "wartości odstające",
        xlabel = "[t] = s",
        marker = (10, 5, :x),
        palette = [:red],
    )
    title!("Ścieżka l6n")
    png("Odstające")
end

function plottrain(pathid, traindf, mlon, mlat)
    singlepath = traindf[traindf[:, :path_id].==pathid, :]
    sampleslabels = unique(singlepath[:, :path_sample])
    singlesamples = [
        singlepath[singlepath[:, :path_sample].==samplelabel, :]
        for samplelabel in sampleslabels
    ]
    lons = [x[:, :lon] for x in singlesamples]
    lats = [x[:, :lat] for x in singlesamples]
    timesteps = [x[:, :timestep] for x in singlesamples]

    sampletopred = singlesamples[1]
    X, y = MagneticLocSuchowiak.preparedata(sampletopred)
    predlon = MLJ.predict(mlon, X)
    predlat = MLJ.predict(mlat, X)
    push!(lons, predlon)
    push!(lats, predlat)
    push!(timesteps, sampletopred[:, :timestep])

    labels = permutedims(vcat(sampleslabels, "prediction"))

    lonplot = plot(
        timesteps,
        lons,
        w = 2,
        title = "Longitude",
        label = labels,
        xlabel = "[t] = s",
        ylabel = "Standarized Longitude",
    )

    latplot = plot(
        timesteps,
        lats,
        w = 2,
        title = "Latitude",
        label = labels,
        xlabel = "[t] = s",
        ylabel = "Standarized Latitude",
    )

    lonlatplot = plot(
        lons,
        lats,
        title = "Path",
        label = labels,
        xlabel = "Standarized Logitude",
        ylabel = "Standarized Latitude",
    )

    display(plot(lonplot, latplot, lonlatplot, layout=grid(3,1)))
end

path1 = traindf[.&(traindf[!, :path_id] .== "l3n", traindf[!, :path_sample] .== "01"), :]
path2 = traindf[.&(traindf[!, :path_id] .== "l5n", traindf[!, :path_sample] .== "01"), :]
l = minimum([nrow(path1), nrow(path2)])
p1m = path1[1:90, :]
p2m = path2[1:90, :]
x = rms(coerce(p1m.magnetometr_x, Continuous), coerce(p2m.magnetometr_x,Continuous))
y = rms(coerce(p1m.magnetometr_y, Continuous), coerce(p2m.magnetometr_y,Continuous))
z = rms(coerce(p1m.magnetometr_z, Continuous), coerce(p2m.magnetometr_z,Continuous))
println(x)
println(y)
println(z)

createcomparisongraphs(
    path1,
    path2,
    axes = ["x", "y", "z"],
    measures = ["magnetometr"],

)

r = traindf[traindf[:, :path_id].=="l1n", :]
createfullgraph(r)
sampleslabels = unique(r[:, :path_sample])
singlesamples = [
    r[r[:, :path_sample].==samplelabel, :]
    for samplelabel in sampleslabels
]
for sample in singlesamples
    maxrow =
end

singlepath = traindf[traindf[:, :path_id].=="l2n", :]
sampleslabels = unique(singlepath[:, :path_sample])
singlesamples = [
    singlepath[singlepath[:, :path_sample].==samplelabel, :]
    for samplelabel in sampleslabels
]
orientx = [x[:, :orientation_x] for x in singlesamples]
orienty = [x[:, :orientation_y] for x in singlesamples]
orientz = [x[:, :orientation_z] for x in singlesamples]
magx = [x[:, :magnetometr_x] for x in singlesamples]
magy = [x[:, :magnetometr_y] for x in singlesamples]
magz = [x[:, :magnetometr_z] for x in singlesamples]
timesteps = [x[:, :timestep] for x in singlesamples]
idxs = [1:nrow(x) for x in singlesamples]
snip = testdf[testdf[:, :path_id].=="tt01", :]
snip = snip[1:166, :]
push!(magx, snip.magnetometr_x)
push!(magy, snip.magnetometr_y)
push!(magz, snip.magnetometr_z)
push!(idxs, 1:nrow(snip))
idxs = idxs ./ 10

Plots.plot(
    timesteps[1],
    [orientx[1], orienty[1], orientz[1]],
    label = ["orientation x" "orientation y" "orientation z"],
    xlabel = "t = [s]",
    w=2,
    layout=3,
    ylim=(-180,180)
)

p1=plot(
    idxs,
    magx,
    w = [1 1 1 1 1 2],
    label = ["train" "train" "train" "train" "train" "test"],
)
p2=plot(
    idxs,
    magy,
    w = [1 1 1 1 1 2],
    label = ["train" "train" "train" "train" "train" "test"],
)
p3=plot(
    idxs,
    magz,
    w = [1 1 1 1 1 2],
    label = ["train" "train" "train" "train" "train" "test"],
)
plot(p1,p2,p3, layout=grid(3,1), size=(900, 1200), xlabel = "t = [s]", title = ["magnetometr_x" "magnetometr_y" "magnetometr_z"])

max_stars_repo_name: SebastianSuchowiak/MagneticLocSuchowiak

max_issues_repo_name: SebastianSuchowiak/MagneticLocSuchowiak

max_forks_repo_name: SebastianSuchowiak/MagneticLocSuchowiak


text:
# K-means algorithm

#### Interface

# C is the type of centers, an (abstract) matrix of size (d x k)
# D is the type of pairwise distance computation from samples to centers
# WC is the type of cluster weights, either Int (in the case where samples are
# unweighted) or eltype(weights) (in the case where samples are weighted).
struct KmeansResult{C<:AbstractMatrix{<:AbstractFloat},D<:Real,WC<:Real} <: ClusteringResult
    centers::C                 # cluster centers (d x k)
    assignments::Vector{Int}   # assignments (n)
    costs::Vector{D}           # cost of the assignments (n)
    counts::Vector{Int}        # number of samples assigned to each cluster (k)
    cweights::Vector{WC}       # cluster weights (k)
    totalcost::D               # total cost (i.e. objective)
    iterations::Int            # number of elapsed iterations
    converged::Bool            # whether the procedure converged
end

const _kmeans_default_init = :kmpp
const _kmeans_default_maxiter = 100
const _kmeans_default_tol = 1.0e-6
const _kmeans_default_display = :none

"""
    kmeans!(X, centers)

Update the current centers `centers` (of size `d x k` where `d` is the dimension and `k` the
number of centroids) using the samples contained in `X` (of size `d x n` where `n` is the number
of samples).
"""
function kmeans!(X::AbstractMatrix{<:Real},                # in: sample matrix (d x n)
                 centers::AbstractMatrix{<:AbstractFloat}; # in: current centers (d x k)
                 weights::Union{Nothing, AbstractVector{<:Real}}=nothing, # in: sample weights (n)
                 maxiter::Integer=_kmeans_default_maxiter, # in: maximum number of iterations
                 tol::Real=_kmeans_default_tol,            # in: tolerance of change at convergence
                 display::Symbol=_kmeans_default_display,  # in: level of display
                 distance::SemiMetric=SqEuclidean())       # in: function to compute distances
    d, n = size(X)
    dc, k = size(centers)

    d == dc || throw(DimensionMismatch("Inconsistent array dimensions for `X` and `centers`."))
    (2 <= k < n) || error("k must have 2 <= k < n.")
    if weights !== nothing
        length(weights) == n || throw(DimensionMismatch("Incorrect length of weights."))
    end

    _kmeans!(X, weights, centers, Int(maxiter), Float64(tol),
             display_level(display), distance)
end

"""
    kmeans(X, k)

K-means clustering with `k` centroids of the data contained in `X` of size `d x n` where `d` is
the dimension and `n` is the number of samples.
"""
function kmeans(X::AbstractMatrix{<:Real},                # in: sample matrix (d x n) columns = obs
                k::Integer;                               # in: number of centers
                weights::Union{Nothing, AbstractVector{<:Real}}=nothing, # in: sample weights (n)
                init::Symbol=_kmeans_default_init,        # in: initialization algorithm
                maxiter::Integer=_kmeans_default_maxiter, # in: maximum number of iterations
                tol::Real=_kmeans_default_tol,            # in: tolerance  of change at convergence
                display::Symbol=_kmeans_default_display,  # in: level of display
                distance::SemiMetric=SqEuclidean())       # in: function to calculate distance with
    d, n = size(X)
    (2 <= k < n) || throw(ArgumentError("k must be 2 <= k < n, k=$k given."))

    # initialize the centers using a type wide enough so that the updates
    # centers[i, cj] += X[i, j] * wj will occur without loss of precision through rounding
    T = float(weights === nothing ? eltype(X) : promote_type(eltype(X), eltype(weights)))
    iseeds = initseeds(init, X, k)
    centers = copyseeds!(Matrix{T}(undef, d, k), X, iseeds)

    kmeans!(X, centers;
            weights=weights, maxiter=Int(maxiter), tol=Float64(tol),
            display=display, distance=distance)
end

#### Core implementation

# core k-means skeleton
function _kmeans!(X::AbstractMatrix{<:Real},                # in: sample matrix (d x n)
                  weights::Union{Nothing, Vector{<:Real}},  # in: sample weights (n)
                  centers::AbstractMatrix{<:AbstractFloat}, # in/out: matrix of centers (d x k)
                  maxiter::Int,                             # in: maximum number of iterations
                  tol::Float64,                             # in: tolerance of change at convergence
                  displevel::Int,                           # in: the level of display
                  distance::SemiMetric)                     # in: function to calculate distance
    d, n = size(X)
    k = size(centers, 2)
    to_update = Vector{Bool}(undef, k) # whether a center needs to be updated
    unused = Vector{Int}()
    num_affected = k # number of centers to which dists need to be recomputed

    # assign containers for the vector of assignments & number of samples assigned to each cluster
    assignments = Vector{Int}(undef, n)
    counts = Vector{Int}(undef, k)

    # compute pairwise distances, preassign costs and cluster weights
    dmat = pairwise(distance, centers, X, dims=2)
    WC = (weights === nothing) ? Int : eltype(weights)
    cweights = Vector{WC}(undef, k)
    D = typeof(one(eltype(dmat)) * one(WC))
    costs = Vector{D}(undef, n)

    update_assignments!(dmat, true, assignments, costs, counts,
                        to_update, unused)
    objv = weights === nothing ? sum(costs) : dot(weights, costs)

    # main loop
    t = 0
    converged = false
    if displevel >= 2
        @printf "%7s %18s %18s | %8s \n" "Iters" "objv" "objv-change" "affected"
        println("-------------------------------------------------------------")
        @printf("%7d %18.6e\n", t, objv)
    end

    while !converged && t < maxiter
        t += 1

        # update (affected) centers
        update_centers!(X, weights, assignments, to_update, centers, cweights)

        if !isempty(unused)
            repick_unused_centers(X, costs, centers, unused, distance)
            to_update[unused] .= true
        end

        if t == 1 || num_affected > 0.75 * k
            pairwise!(dmat, distance, centers, X, dims=2)
        else
            # if only a small subset is affected, only compute for that subset
            affected_inds = findall(to_update)
            pairwise!(view(dmat, affected_inds, :), distance,
                      view(centers, :, affected_inds), X, dims=2)
        end

        # update assignments
        update_assignments!(dmat, false, assignments, costs, counts,
                            to_update, unused)
        num_affected = sum(to_update) + length(unused)

        # compute change of objective and determine convergence
        prev_objv = objv
        objv = weights === nothing ? sum(costs) : dot(weights, costs)
        objv_change = objv - prev_objv

        if objv_change > tol
            @warn("The clustering cost increased at iteration #$t")
        elseif abs(objv_change) < tol
            converged = true
        end

        # display information (if required)
        if displevel >= 2
            @printf("%7d %18.6e %18.6e | %8d\n", t, objv, objv_change, num_affected)
        end
    end

    if displevel >= 1
        if converged
            println("K-means converged with $t iterations (objv = $objv)")
        else
            println("K-means terminated without convergence after $t iterations (objv = $objv)")
        end
    end

    return KmeansResult(centers, assignments, costs, counts,
                        cweights, objv, t, converged)
end

#
#  Updates assignments, costs, and counts based on
#  an updated (squared) distance matrix
#
function update_assignments!(dmat::Matrix{<:Real},     # in:  distance matrix (k x n)
                             is_init::Bool,            # in:  whether it is the initial run
                             assignments::Vector{Int}, # out: assignment vector (n)
                             costs::Vector{<:Real},    # out: costs of the resultant assignment (n)
                             counts::Vector{Int},      # out: # samples assigned to each cluster (k)
                             to_update::Vector{Bool},  # out: whether a center needs update (k)
                             unused::Vector{Int}       # out: list of centers with no samples
                             )
    k, n = size(dmat)

    # re-initialize the counting vector
    fill!(counts, 0)

    if is_init
        fill!(to_update, true)
    else
        fill!(to_update, false)
        if !isempty(unused)
            empty!(unused)
        end
    end

    # process each sample
    @inbounds for j = 1:n
        # find the closest cluster to the i-th sample. Note that a
        # is necessarily between 1 and size(dmat, 1) === k as a result
        # and can thus be used as an index in an `inbounds` environment
        c, a = findmin(view(dmat, :, j))

        # set/update the assignment
        if is_init
            assignments[j] = a
        else  # update
            pa = assignments[j]
            if pa != a
                # if assignment changes,
                # both old and new centers need to be updated
                assignments[j] = a
                to_update[a] = true
                to_update[pa] = true
            end
        end

        # set costs and counts accordingly
        costs[j] = c
        counts[a] += 1
    end

    # look for centers that have no associated samples
    for i = 1:k
        if counts[i] == 0
            push!(unused, i)
            to_update[i] = false # this is handled using different mechanism
        end
    end
end

#
#  Update centers based on updated assignments
#
#  (specific to the case where samples are not weighted)
#
function update_centers!(X::AbstractMatrix{<:Real},        # in: sample matrix (d x n)
                         weights::Nothing,                 # in: sample weights
                         assignments::Vector{Int},         # in: assignments (n)
                         to_update::Vector{Bool},          # in: whether a center needs update (k)
                         centers::AbstractMatrix{<:AbstractFloat}, # out: updated centers (d x k)
                         cweights::Vector{Int})            # out: updated cluster weights (k)
    d, n = size(X)
    k = size(centers, 2)

    # initialize center weights
    cweights[to_update] .= 0

    # accumulate columns
    @inbounds for j in 1:n
        # skip samples assigned to a center that doesn't need to be updated
        cj = assignments[j]
        if to_update[cj]
            if cweights[cj] > 0
                for i in 1:d
                    centers[i, cj] += X[i, j]
                end
            else
                for i in 1:d
                    centers[i, cj] = X[i, j]
                end
            end
            cweights[cj] += 1
        end
    end

    # sum ==> mean
    @inbounds for j in 1:k
        if to_update[j]
            cj = cweights[j]
            for i in 1:d
                centers[i, j] /= cj
            end
        end
    end
end

#
#  Update centers based on updated assignments
#
#  (specific to the case where samples are weighted)
#
function update_centers!(X::AbstractMatrix{<:Real}, # in: sample matrix (d x n)
                         weights::Vector{W},        # in: sample weights (n)
                         assignments::Vector{Int},  # in: assignments (n)
                         to_update::Vector{Bool},   # in: whether a center needs update (k)
                         centers::AbstractMatrix{<:Real}, # out: updated centers (d x k)
                         cweights::Vector{W}        # out: updated cluster weights (k)
                         ) where W<:Real
    d, n = size(X)
    k = size(centers, 2)

    # initialize center weights
    cweights[to_update] .= 0

    # accumulate columns
    @inbounds for j in 1:n
        # skip samples with negative weights or assigned to a center
        # that doesn't need to be updated
        wj = weights[j]
        cj = assignments[j]
        if wj > 0 && to_update[cj]
            if cweights[cj] > 0
                for i in 1:d
                    centers[i, cj] += X[i, j] * wj
                end
            else
                for i in 1:d
                    centers[i, cj] = X[i, j] * wj
                end
            end
            cweights[cj] += wj
        end
    end

    # sum ==> mean
    @inbounds for j in 1:k
        if to_update[j]
            cj = cweights[j]
            for i in 1:d
                centers[i, j] /= cj
            end
        end
    end
end

#
#  Re-picks centers that get no samples assigned to them.
#
function repick_unused_centers(X::AbstractMatrix{<:Real}, # in: the sample set (d x n)
                               costs::Vector{<:Real},     # in: the current assignment costs (n)
                               centers::AbstractMatrix{<:AbstractFloat}, # out: the centers (d x k)
                               unused::Vector{Int},       # in: indices of centers to be updated
                               distance::SemiMetric)      # in: function to calculate the distance with
    # pick new centers using a scheme like kmeans++
    ds = similar(costs)
    tcosts = copy(costs)
    n = size(X, 2)

    for i in unused
        j = wsample(1:n, tcosts)
        tcosts[j] = 0
        v = view(X, :, j)
        centers[:, i] = v
        colwise!(ds, distance, v, X)
        tcosts = min(tcosts, ds)
    end
end

ext: jl


text:
export SentenceSplitter

using ..Units
using WordTokenizers

struct SentenceSplitter <: AbstractSentenceSplitter end;

function (splitter::SentenceSplitter)(text)
  WordTokenizers.split_sentences(text) .|> StringDocument
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
"""
    CalculusWithJulia

A package to accompany [notes](calculuswithjulia.github.io) on using Julia with calculus.

This package does two things:

* It loads a few other packages making it easier to use (and install) the functionality provided by them and

* It defines a handful of functions for convenience. The exported ones
are `e`, `unzip`, `rangeclamp` `tangent`, `secant`, `D` (and the prime
notation), `divergence`, `gradient`, `curl`, and `∇`.

## Packages loaded by `CalculusWithJulia`

* The `SpecialFunctions` is loaded giving access to a few special functions used in these notes, e.g., `airyai`, `gamma`

* The `ForwardDiff` package is loaded giving access to its  `derivative`,  `gradient`, `jacobian`, and `hessian` functions for finding automatic derivatives of functions. In addition, this package defines `'` (for functions) to return a derivative (which commits [type piracy](https://docs.julialang.org/en/v1/manual/style-guide/index.html#Avoid-type-piracy-1)), `∇` to find the gradient (`∇(f)`), the divergence (`∇⋅F`). and the curl (`∇×F`), along with `divergence` and `curl`.

* The `LinearAlgebra` package is loaded for access to several of its functions for working with vectors `norm`, `cdot` (`⋅`), `cross` (`×`), `det`.

* The `PlotUtils` package is loaded so that its `adapted_grid` function is available.

## Packages with extra features added when loaded

The `Julia` package `Requires` allows for additional code to be run when another package is loaded. The following packages have additional code to load:

* `SymPy`: for symbolic math.

* `Plots`: the `Plots` package provides a plotting interface.

Several plot recipes are provided to ease the creation of plots in the notes.
`plotif`, `trimplot`, and `signchart` are used for plotting univariate functions;
`plot_polar` and `plot_parametric` are used to plot curves in 2 or 3 dimensions;
`plot_parametric` also makes the plotting og parameterically defined surfaces easier;
`vectorfieldplot` and `vectorfieldplot3d` can be used to plot vector fields; and
`arrow` is a simplified interface to `quiver` that also indicates 3D vectors.

* `MDBM`: the  `MDBM` package is useful for finding contours in ``2`` or ``3`` dimensions.

## Other packages with recurring roles in the notes

* `Roots` is used to find zeros of univariate functions

* `QuadGK` and `HCubature` are used for numeric integration

* `ImplicitEquations` is used to make 2-dimesional plots of implicit equations.

"""
module CalculusWithJulia

using Reexport
using RecipesBase
import PlotUtils
import ForwardDiff
export ForwardDiff

@reexport using LinearAlgebra
@reexport using SpecialFunctions

using EllipsisNotation
import EllipsisNotation: Ellipsis
export ..
import IntervalSets
import IntervalSets: ClosedInterval

include("multidimensional.jl")
include("derivatives.jl")
include("integration.jl")
include("plot-utils.jl")
include("plot-recipes.jl")
include("WeaveSupport/WeaveSupport.jl")

using Requires

function __init__()
    @require SymPy="24249f21-da20-56a4-8eb1-6a02cf4ae2e6" include("sympy.jl")
    @require Polynomials="f27b6e38-b328-58d1-80ce-0feddd5e7a45" include("polynomials.jl")
    @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("plots.jl")
    @require MDBM="ea0cff06-48de-41e3-bd0e-d3c1feffd247" include("implicit_equation.jl")
end

e = exp(1)
export e

export unzip, rangeclamp
export tangent, secant, D
export riemann
export divergence, gradient, curl, ∇

end # module


text:
using ReservoirComputing, DifferentialEquations, Statistics, Random

u0 = [1.0,0.0,0.0]
tspan = (0.0,1000.0)
datasize = 100000
tsteps = range(tspan[1], tspan[2], length = datasize)

function lorenz(du,u,p,t)
    p = [10.0,28.0,8/3]
    du[1] = p[1]*(u[2]-u[1])
    du[2] = u[1]*(p[2]-u[3]) - u[2]
    du[3] = u[1]*u[2] - p[3]*u[3]
end

function prior_model_data_generator(u0, tspan, tsteps, model = lorenz)
    prob = ODEProblem(lorenz, u0, tspan)
    sol = Array(solve(prob, saveat = tsteps))
    return sol
end

train_len = 10000

ode_prob = ODEProblem(lorenz, u0, tspan)
ode_sol = solve(ode_prob, saveat = tsteps)
ode_data =Array(ode_sol)
input_data = ode_data[:, 1:train_len]
target_data = ode_data[:, 2:train_len+1]

test_data = ode_data[:, train_len+1:end][:, 1:1000]
predict_len = size(test_data, 2)
tspan_train = (tspan[1], ode_sol.t[train_len])

hybrid = Hybrid(prior_model_data_generator, u0, tspan_train, train_len)

Random.seed!(77)
esn = ESN(input_data,
    reservoir = RandSparseReservoir(300),
    variation = hybrid)

output_layer = train(esn, target_data, StandardRidge(0.3))

output = esn(Generative(predict_len), output_layer)

@test mean(abs.(test_data[1:100] .- output[1:100])) ./ mean(abs.(test_data[1:100])) < 0.11

max_stars_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_issues_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_forks_repo_name: MartinuzziFrancesco/EchoStateNetwork


text:
function load_league(league :: String)
    mkpath("form/$league")
    df = CSV.File(download("https://projects.fivethirtyeight.com/soccer-api/club/spi_matches.csv")) |> DataFrame
    return @linq df |> where(:league .== league)
end

max_stars_repo_name: shalinkpatel/football-form

max_issues_repo_name: shalinkpatel/football-form

max_forks_repo_name: shalinkpatel/football-form


text:
function compute_RMSE(U, V, Y, r, d1, d2, rows_t, vals_t, cols_t)
	res = 0.0
	n = 0.0
	for i = 1:d1
		tmp = nzrange(Y, i)
		d2_bar = rows_t[tmp];
		vals_d2_bar = vals_t[tmp];
		ui = U[:, i]
		len = size(d2_bar)[1]
		for j = 1:len
			J = d2_bar[j];
			vj = V[:, J]
			res += (vals_d2_bar[j] - dot(ui,vj))^2
			n += 1.0
		end
	end
	return (res / n)^0.5
end

function compute_RMSE_train(U, V, X, r, d1, d2, rows, vals, cols)
	res = 0.0
	n = 0.0
	for i = 1:d1
		tmp = nzrange(X, i)
		d2_bar = rows[tmp];
		vals_d2_bar = vals[tmp];
		ui = U[:, i]
		len = size(d2_bar)[1]
		for j = 1:len
			J = d2_bar[j];
			vj = V[:, J]
			res += (vals_d2_bar[j] - dot(ui,vj))^2
			n += 1.0
		end
	end
	return (res / n)^0.5
end

function objective(U, V, X, d1, lambda, rows, vals)
	res = 0.0
	res = lambda * (vecnorm(U) ^ 2 + vecnorm(V) ^ 2)
	for i in 1:d1
		tmp = nzrange(X, i)
		d2_bar = rows[tmp];
		vals_d2_bar = vals[tmp];
		len = size(d2_bar)[1];
		for k in 1:len
			j = d2_bar[k]
			res += (vals_d2_bar[k] - dot(U[:,i], V[:,j]))^2
		end
	end
	return res
end

function compute_pairwise_error_ndcg(U, V, Y, r, d1, d2, rows_t, vals_t, cols_t, ndcg_k)
	sum_error = 0.; ndcg_sum = 0.;
	for i = 1:d1
		tmp = nzrange(Y, i)
		d2_bar = rows_t[tmp];
		vals_d2_bar = vals_t[tmp];
		ui = U[:, i]
		len = size(d2_bar)[1]
		score = zeros(len)
		for j = 1:len
			J = d2_bar[j];
			vj = V[:, J]
			score[j] = dot(ui,vj)
		end
		error_this = 0
		n_comps_this = 0
		for j in 1:(len - 1)
			jval = vals_d2_bar[j]
			for k in (j + 1):len
				kval = vals_d2_bar[k]
				if score[j] >= score[k] && jval < kval
					error_this += 1
				end
				if score[j] <= score[k] && jval > kval
					error_this += 1
				end
				n_comps_this += 1
			end
		end
		sum_error += error_this / n_comps_this

		p1 = sortperm(score, rev = true)
		p1 = p1[1:ndcg_k]
		M1 = vals_d2_bar[p1]
		p2 = sortperm(vals_d2_bar, rev = true)
		p2 = p2[1:ndcg_k]
		M2 = vals_d2_bar[p2]
		dcg = 0.; dcg_max = 0.
		for k = 1:ndcg_k
			dcg += (2 ^ M1[k] - 1) / log2(k + 1)
			dcg_max += (2 ^ M2[k] - 1) / log2(k + 1)
		end
		ndcg_sum += dcg / dcg_max
	end
	return sum_error / d1, ndcg_sum / d1
end

function update(U, V, X, r, d1, d2, lambda, rows, vals, stepsize, cols, threshold_user, threshold_item)
	l = rand(1:length(rows))
	j = rows[l]
	i = cols[l]

	# SSE for user side
	if rand(1)[1] > threshold_user
		i = rand(1:d1)
	end

    	# SSE for item side
	if rand(1)[1] > threshold_item
		j = rand(1:d2)
	end

	eij = vals[l] - dot(U[:,i], V[:,j])
	ui = U[:,i] + stepsize * (eij * V[:,j] - lambda * U[:,i])
	vj = V[:,j] + stepsize * (eij * U[:,i] - lambda * V[:,j])
	for k in 1:r
		U[k, i] = ui[k]
	end
	for k in 1:r
		V[k, j] = vj[k]
	end
	return U, V
end

function main(train, test, r, lambda, threshold_user, threshold_item)
	X = readdlm(train, ',' , Int64);
	x = vec(X[:,1]);
	y = vec(X[:,2]);
	v = vec(X[:,3]);
	Y = readdlm(test, ',' , Int64);
	xx = vec(Y[:,1]);
	yy = vec(Y[:,2]);
	vv = vec(Y[:,3]);
	# userid; movieid
	n = max(maximum(x), maximum(xx)); msize = max(maximum(y), maximum(yy));
	println("Training dataset ", train, " and test dataset ", test, " are loaded. \n There are ", n, " users and ", msize, " items in the dataset.")
	X = sparse(x, y, v, n, msize); # userid by movieid
	Y = sparse(xx, yy, vv, n, msize);
	# julia column major
	# now moveid by userid
	X = X';
	Y = Y';

	# too large to debug the algorithm, subset a small set: 500 by 750
	rows = rowvals(X);
	vals = nonzeros(X);
	cols = zeros(Int, size(vals)[1]);

	d2, d1 = size(X);
	cc = 0;
	for i = 1:d1
		tmp = nzrange(X, i);
		nowlen = size(tmp)[1];
		for j = 1:nowlen
			cc += 1
			cols[cc] = i
		end
	end

	rows_t = rowvals(Y);
	vals_t = nonzeros(Y);
	cols_t = zeros(Int, size(vals_t)[1]);
	cc = 0;
	for i = 1:d1
		tmp = nzrange(Y, i);
		nowlen = size(tmp)[1];
		for j = 1:nowlen
			cc += 1
			cols_t[cc] = i
		end
	end

	# initialize U, V
	U = randn(r, d1);
	V = randn(r, d2);

	stepsize = 0.01
	totaltime = 0.00000;
	println("iter time objective_function train_RMSE test_RMSE")

	nowobj = objective(U, V, X, d1, lambda, rows, vals)
	rmse = compute_RMSE(U, V, Y, r, d1, d2, rows_t, vals_t, cols_t)
	rmse_tr = compute_RMSE_train(U, V, X, r, d1, d2, rows, vals, cols)
	println("[", 0, ", ", totaltime, ", ", nowobj, ", ", rmse_tr, ", ", rmse, "],")

	for iter in 1:3000000000
		tic();

		U, V = update(U, V, X, r, d1, d2, lambda, rows, vals, stepsize, cols, threshold_user, threshold_item)

		totaltime += toq();

	 	if iter % 10000000 == 0
	 		nowobj = objective(U, V, X, d1, lambda, rows, vals)
	 		rmse = compute_RMSE(U, V, Y, r, d1, d2, rows_t, vals_t, cols_t)
	 		rmse_tr = compute_RMSE_train(U, V, X, r, d1, d2, rows, vals, cols)
			println("[", iter, ", ", totaltime, ", ", nowobj, ", ", rmse_tr, ", ", rmse, "],")
		end
	end
	return V, U
end


text:
########################################################################
# Server
########################################################################

abstract type AbstractServer end

include("fedProxServer.jl")
include("ScaffoldServer.jl")
include("fedDCDServer.jl")


text:
function get_slider(;
    timestep=0.01,
    gravity=-9.81,
    spring=0.0,
    damper=0.0,
    T=Float64)

    # Parameters
    joint_axis = [0.0; 0.0; 1.0]
    len = 1.0
    width, depth = 0.1, 0.1
    child_vertex = [0.0; 0.0; len / 2.0] # joint connection point

    # Links
    origin = Origin{T}()
    pbody = Box(width, depth, len, len)

    # Constraints
    joint_between_origin_and_pbody = JointConstraint(Prismatic(origin, pbody, joint_axis;
        child_vertex=child_vertex,
        spring=spring,
        damper=damper))

    bodies = [pbody]
    joints = [joint_between_origin_and_pbody]

    mech = Mechanism(origin, bodies, joints,
        gravity=gravity,
        timestep=timestep,
        spring=spring,
        damper=damper)

    return mech
end

function initialize_slider!(mechanism::Mechanism{T};
    position=0.0) where T

    body = mechanism.bodies[1]
    joint = mechanism.joints[1]
    child_vertex = joint.translational.vertices[2]
    set_maximal_configurations!(mechanism.origin, body,
        child_vertex=child_vertex - [0, 0, position])
end

function get_nslider(;
    timestep=0.01,
    gravity=-9.81,
    spring=0.0,
    damper=0.0,
    num_bodies=5,
    T=Float64)

    # Parameters
    ex = [0.0; 0.0; 1.0]
    h = 1.0
    r = 0.05
    vert11 = [0.0; r; 0.0]
    vert12 = -vert11

    # Links
    origin = Origin{T}()
    bodies = [Cylinder(r, h, h, color=RGBA(1.0, 0.0, 0.0)) for i = 1:num_bodies]

    # Constraints
    jointb1 = JointConstraint(Prismatic(origin, bodies[1], ex;
        child_vertex = 0 * vert11))

    if num_bodies > 1
        joints = [
            jointb1;
            [JointConstraint(Prismatic(bodies[i - 1], bodies[i], ex;
                parent_vertex=vert12,
                child_vertex=vert11,
                spring=spring,
                damper=damper)) for i = 2:num_bodies]
            ]
    else
        joints = [jointb1]
    end

    mech = Mechanism(origin, bodies, joints,
        gravity=gravity,
        timestep=timestep)

    return mech
end

function initialize_nslider!(mechanism::Mechanism{T};
    position=0.2,
    relative_position=0.0) where T

    pbody = mechanism.bodies[1]

    # set position and velocities
    set_maximal_configurations!(mechanism.origin, pbody,
        parent_vertex=[0.0, 0.0, position])

    # set relative positions
    previd = pbody.id
    for body in mechanism.bodies[2:end]
        set_maximal_configurations!(get_body(mechanism, previd), body,
            parent_vertex=[0.0, -0.1, relative_position])
        previd = body.id
    end

    zero_velocity!(mechanism)
end


text:
module Caffe

using ProtoBuf
import ..Merlin
using ..Merlin: Var

include("caffe_pb.jl")

#function __init__()
#end

"""
Load Caffe model.

## 👉 Example
```julia
using Merlin
using Merlin.Caffe

path = "<xxx.caffemodel>"
model = Caffe.load(path)
```
"""
function load(path)
  np = open(path) do io
    readproto(io, NetParameter())
  end

  if length(np.layer) > 0
    layers = np.layer
    dict = Dict(
      "Convolution" => conv,
      "Data" => data,
      "Dropout" => dropout,
      "InnerProduct" => inner_product,
      "Pooling" => pooling,
      "ReLU" => relu,
      "Softmax" => softmax,
      "SoftmaxWithLoss" => softmax_loss,
    )
  else
    layers = np.layers
    t = V1LayerParameter_LayerType
    dict = Dict(
      t.CONVOLUTION => conv,
      t.DATA => data,
      t.DROPOUT => dropout,
      t.INNER_PRODUCT => inner_product,
      t.POOLING => pooling,
      t.RELU => relu,
      t.SOFTMAX => softmax,
      t.SOFTMAX_LOSS => softmax_loss,
    )
  end

  x = Var(:x)
  names = Dict()
  for l in layers
    names[l.name] = l
    f = dict[l._type]
    #x = f(x)
  end
  x
  #Merlin.@graph x
end

function conv(layer)
  blobs = layer.blobs
  w = blobs[1].data
  b = length(blobs) > 1 ? blobs[2].data : nothing
  num = blobs[1].num > 0 ? Int(blobs[1].num) : Int(blobs[1].shape.dim[1])
  channels = blobs[1].channels > 0 ? Int(blobs[1].channels) : Int(blobs[1].shape.dim[2])
  param = layer.convolution_param
  kernel, stride, pad = kernel_size(param), stride_size(param), pad_size(param)

  # group
  #=
  ngroups = param.group
  n_in = channels * ngroups
  n_out = num
  for i = 1:ngroups
    in_slice = slice((i-1) * n_in, i * n_in)
    out_slice = slice((i-1) * n_out, i * n_out)
    w = func.W.data[out_slice, in_slice]
  end
  p.group == 1 && (w = reshape(w, kernel[1], kernel[2], channels, num))
  =#

  Merlin.Conv(Var(w), Var(b), stride, pad)
end

function data(layer)
end

function dropout(layer)
  ratio = Float64(layer.dropout_param.dropout_ratio)
  Merlin.Dropout(ratio)
end

function inner_product(layer)
  param = layer.inner_product_param
  if param.axis != 1
    throw("Axis: $(axis) in inner-product is unsupported.")
  end

  blobs = layer.blobs
  b = blobs[1]
  if blob.height > 0
    height, width = Int(b.height), Int(b.width)
  elseif length(b.shape.dim) == 2
    height, width = Int(b.shape.dim[1]), Int(b.shape.dim[2])
  elseif length(b.shape.dim) == 4
    height, width = Int(b.shape.dim[3]), Int(b.shape.dim[4])
  else
    throw("Dimensions: $(b.shape.dim) in inner-product is unsupported.")
  end

  w = reshape(blobs[1].data, height, width)
  b = param.bias_term ? vec(blobs[2].data) : nothing
  Merlin.Linear(Var(w), Var(b))
end

function pooling(layer)
  param = layer.convolution_param
  kernel, stride, pad = kernel_size(param), stride_size(param), pad_size(param)
  types = PoolingParameter_PoolMethod
  if p.pool == types.MAX
    mode = :max
  elseif p.pool == types.AVE
    mode = :ave
  else
    warn("Mode: $(mode) in Pooling is unsupported. Set to max.")
    mode = :max
  end
  Merlin.Pooling(mode, kernel, stride, pad)
end

function relu(layer)
  slope = Float64(layer.relu_param.negative_slope)
  if slope != 0.0
    warn("Leaky ReLU is unsupported. Slope is set to 0.")
  end
  Merlin.relu
end

function softmax(layer)
  axis = Int(layer.softmax_param.axis)
  Merlin.Softmax(axis)
end

function softmax_loss(layer)
  axis = Int(layer.softmax_param.axis)
  Merlin.SoftmaxCrossEntropy(axis)
end

########## Internal ##########

function kernel_size(param)
  param.kernel_h > 0 && return Int(param.kernel_h), Int(param.kernel_w)
  s = Int(param.kernel_size[1])
  s, s
end

function stride_size(param)
  param.stride_h > 0 && return Int(param.stride_h), Int(param.stride_w)
  length(param.stride) == 0 && return 1, 1
  s = Int(param.stride[1])
  s, s
end

function pad_size(param)
  param.pad_h > 0 && return Int(param.pad_h), Int(param.pad_w)
  length(param.pad) == 0 && return 0, 0
  s = Int(param.pad[1])
  s, s
end

end


text:
@testset "WeightedSetMosaic" begin
    using OptEnrichedSetCover: WeightedSetMosaic

    @testset "empty" begin
        sm = SetMosaic(Set{Symbol}[])
        wsm = assignweights(sm, [Dict{Int, Float64}()])
        @test wsm isa WeightedSetMosaic
        @test originalmosaic(wsm) == sm
        @test nelements(wsm) == 0
        @test nsets(wsm) == 0
        @test nexperiments(wsm) == 1
    end

    @testset "empty but with elements" begin
        sm = SetMosaic(Set{Symbol}[], Set([:a, :b]))
        wsm = assignweights(sm, [Dict{Int, Float64}()])
        @test wsm isa WeightedSetMosaic
        @test originalmosaic(wsm) == sm
        @test nelements(wsm) == 2
        @test nsets(wsm) == 0
        @test nexperiments(wsm) == 1

        # named
        nwsm = assignweights(sm, Dict(:X => Dict{Int, Float64}()))
        @test originalmosaic(nwsm) == sm
        @test nelements(nwsm) == 2
        @test nsets(nwsm) == 0
        @test nexperiments(nwsm) == 1
    end

    @testset "[:a :b] [:c :d] [:a :b :c :d]" begin
        sm = SetMosaic([Set([:a, :b]), Set([:c, :d]), Set([:a, :b, :c, :d])])

        wsm = assignweights(sm, [Dict(1 => -5.0, 3 => -2.0)])

        @test nelements(wsm) == 4
        @test nexperiments(wsm) == 1
        @test nsets(wsm) == 2
        @test setweight(wsm, 1, 1) == -5
        @test setweight(wsm, 2, 1) === missing
        @test setweight(wsm, 3, 1) == -2
        @test_throws KeyError setweight(wsm, 1, 2)
        @test_throws KeyError setweight(wsm, 4, 1)

        wsm_copy = copy(wsm)
        @test nelements(wsm_copy) == nelements(wsm)
        @test nexperiments(wsm_copy) == nexperiments(wsm)
        @test wsm_copy.original === wsm.original
        @test wsm_copy.setXexp_weights !== wsm.setXexp_weights
        @test wsm_copy.setXexp_weights == wsm.setXexp_weights

        # ignore weight of nonexisting set
        wsm2 = assignweights(sm, [Dict(4 => -8.0)])
        @test nelements(wsm2) == 4
        @test nsets(wsm2) == 0
        @test_throws KeyError setweight(wsm2, 4, 1)

        # mask with maweightsx_min_weight, [:a :b :c :d] is excluded
        wsm3 = assignweights(sm, [Dict(1 => -5.0, 3 => -2.0)], max_min_weight=-2.5)
        @test nelements(wsm3) == 4
        @test nsets(wsm3) == 1
        @test setweight(wsm3, 1, 1) == -5
        @test setweight(wsm3, 2, 1) === missing
        @test setweight(wsm3, 3, 1) === missing

        wsm4 = assignweights(sm, [Dict(1 => -5.0, 3 => -2.0)], max_setsize=3)
        @test nelements(wsm4) == 4
        @test nsets(wsm4) == 1
        @test setweight(wsm4, 1, 1) == -5
        @test setweight(wsm4, 2, 1) === missing
        @test setweight(wsm4, 3, 1) === missing
    end

    @testset "A=[:a :b] B=[:c :d] C=[:a :b :c :d]" begin
        sm = SetMosaic(Dict(:A=>Set([:a, :b]), :B=>Set([:c, :d]), :C=>Set([:a, :b, :c, :d])))
        wsm = assignweights(sm, [Dict(:A => -5.0, :C => -2.0)])

        @test nelements(wsm) == 4
        @test nsets(wsm) == 2
        @test nexperiments(wsm) == 1
        @test setweight(wsm, :A, 1) == -5
        @test setweight(wsm, :B, 1) === missing
        @test setweight(wsm, :C, 1) == -2
    end

    @testset "multimask: A=[:a :b] B=[:c :d] C=[:a :b :c :d], mask=[[:a :b] [:a :b :c] [:a :b :d]]" begin
        sm = SetMosaic(Dict(:A=>Set([:a, :b]), :B=>Set([:c, :d]), :C=>Set([:a, :b, :c, :d])))
        setweights = [Dict(:A => -5.0,             :C => -2.0),
                      Dict(:A => -6.0, :B => -1.0, :C => -8.0),
                      Dict(            :B => -2.0, :C => -3.0)]
        wsm = assignweights(sm, setweights)
        @test nelements(wsm) == 4
        @test nsets(wsm) == 3
        @test nexperiments(wsm) == 3
        @test_throws KeyError setweight(wsm, :A, 4)
        @test_throws KeyError setweight(wsm, :D, 1)
        @test setweight(wsm, :A, 1) == -5
        @test setweight(wsm, :A, 2) == -6
        @test setweight(wsm, :A, 3) === missing
        @test setweight(wsm, :B, 1) === missing
        @test setweight(wsm, :B, 2) == -1
        @test setweight(wsm, :B, 3) == -2
        @test setweight(wsm, :C, 1) == -2
        @test setweight(wsm, :C, 2) == -8
        @test setweight(wsm, :C, 3) == -3

        wsm2 = assignweights(sm, setweights, max_min_weight=-3)
        @test nsets(wsm2) == 2
        @test setweight(wsm2, :A, 1) == -5
        @test setweight(wsm2, :A, 2) == -6
        @test setweight(wsm2, :A, 3) === missing
        @test setweight(wsm2, :B, 1) === missing
        @test setweight(wsm2, :B, 2) === missing
        @test setweight(wsm2, :B, 3) === missing
        @test setweight(wsm2, :C, 1) == -2
        @test setweight(wsm2, :C, 2) == -8
        @test setweight(wsm2, :C, 3) == -3

        nsetweights = Dict(:X => setweights[1],
                           :Y => setweights[2],
                           :Z => setweights[3])
        nwsm2 = assignweights(sm, nsetweights, max_weight=-2, max_setsize=3)
        @test nsets(nwsm2) == 2 # C filtered out because it's too big
        @test nexperiments(nwsm2) == 3
        @test_throws KeyError setweight(nwsm2, :A, :U)
        @test_throws KeyError setweight(nwsm2, :D, :X)
        @test setweight(nwsm2, :A, :X) == -5
        @test setweight(nwsm2, :A, :Y) == -6
        @test setweight(nwsm2, :A, :Z) === missing
        @test setweight(nwsm2, :B, :X) === missing
        @test setweight(nwsm2, :B, :Y) === missing
        @test setweight(nwsm2, :B, :Z) == -2
        @test setweight(nwsm2, :C, :X) === missing
        @test setweight(nwsm2, :C, :Y) === missing
        @test setweight(nwsm2, :C, :Z) === missing
    end
end


text:
function build_incomebalance!(m, jm)
    # Add income balance constraints
    for c in m._demands
        level_name = get_name(c.consumer)

        ex6a = :(
            JuMP.@NLexpression(
                $jm,
                +($((:($(swap_our_param_with_jump_param(jm, en.quantity)) *
                $(get_jump_variable_for_commodity(jm, en.commodity))) for en in c.endowments)...)) -
                $(jm[level_name])
            )
        )

        ex6b = eval(swap_our_param_with_jump_param(jm, ex6a))
        Complementarity.add_complementarity(jm, jm[level_name], ex6b, string("F_", level_name))
        push!(m._nlexpressions, ex6b)
    end
end


text:
# Save current mesh with some context information as an HDF5 file.
function save_mesh_file(mesh::TreeMesh, output_directory, timestep=0)
  save_mesh_file(mesh, output_directory, timestep, mpi_parallel(mesh))
end

function save_mesh_file(mesh::TreeMesh, output_directory, timestep,
                        mpi_parallel::Val{false})
  # Create output directory (if it does not exist)
  mkpath(output_directory)

  # Determine file name based on existence of meaningful time step
  if timestep > 0
    filename = joinpath(output_directory, @sprintf("mesh_%06d.h5", timestep))
  else
    filename = joinpath(output_directory, "mesh.h5")
  end

  # Open file (clobber existing content)
  h5open(filename, "w") do file
    # Add context information as attributes
    n_cells = length(mesh.tree)
    attributes(file)["ndims"] = ndims(mesh)
    attributes(file)["n_cells"] = n_cells
    attributes(file)["n_leaf_cells"] = count_leaf_cells(mesh.tree)
    attributes(file)["minimum_level"] = minimum_level(mesh.tree)
    attributes(file)["maximum_level"] = maximum_level(mesh.tree)
    attributes(file)["center_level_0"] = mesh.tree.center_level_0
    attributes(file)["length_level_0"] = mesh.tree.length_level_0
    attributes(file)["periodicity"] = collect(mesh.tree.periodicity)

    # Add tree data
    file["parent_ids"] = @view mesh.tree.parent_ids[1:n_cells]
    file["child_ids"] = @view mesh.tree.child_ids[:, 1:n_cells]
    file["neighbor_ids"] = @view mesh.tree.neighbor_ids[:, 1:n_cells]
    file["levels"] = @view mesh.tree.levels[1:n_cells]
    file["coordinates"] = @view mesh.tree.coordinates[:, 1:n_cells]
  end

  return filename
end

# Save current mesh with some context information as an HDF5 file.
function save_mesh_file(mesh::TreeMesh, output_directory, timestep,
                        mpi_parallel::Val{true})
  # Create output directory (if it does not exist)
  mpi_isroot() && mkpath(output_directory)

  # Determine file name based on existence of meaningful time step
  if timestep >= 0
    filename = joinpath(output_directory, @sprintf("mesh_%06d.h5", timestep))
  else
    filename = joinpath(output_directory, "mesh.h5")
  end

  # Since the mesh is replicated on all ranks, only save from MPI root
  if !mpi_isroot()
    return filename
  end

  # Open file (clobber existing content)
  h5open(filename, "w") do file
    # Add context information as attributes
    n_cells = length(mesh.tree)
    attributes(file)["ndims"] = ndims(mesh)
    attributes(file)["n_cells"] = n_cells
    attributes(file)["n_leaf_cells"] = count_leaf_cells(mesh.tree)
    attributes(file)["minimum_level"] = minimum_level(mesh.tree)
    attributes(file)["maximum_level"] = maximum_level(mesh.tree)
    attributes(file)["center_level_0"] = mesh.tree.center_level_0
    attributes(file)["length_level_0"] = mesh.tree.length_level_0
    attributes(file)["periodicity"] = collect(mesh.tree.periodicity)

    # Add tree data
    file["parent_ids"] = @view mesh.tree.parent_ids[1:n_cells]
    file["child_ids"] = @view mesh.tree.child_ids[:, 1:n_cells]
    file["neighbor_ids"] = @view mesh.tree.neighbor_ids[:, 1:n_cells]
    file["levels"] = @view mesh.tree.levels[1:n_cells]
    file["coordinates"] = @view mesh.tree.coordinates[:, 1:n_cells]
  end

  return filename
end

"""
    load_mesh(restart_file::AbstractString; n_cells_max)

Load the mesh from the `restart_file`.
"""
function load_mesh(restart_file::AbstractString; n_cells_max)
  if mpi_isparallel()
    load_mesh_parallel(restart_file; n_cells_max=n_cells_max)
  else
    load_mesh_serial(restart_file; n_cells_max=n_cells_max)
  end
end

function load_mesh_serial(restart_file::AbstractString; n_cells_max)
  ndims_ = h5open(restart_file, "r") do file
    read(attributes(file)["ndims"])
  end

  mesh = TreeMesh(SerialTree{ndims_}, n_cells_max)
  load_mesh!(mesh, restart_file)
end

function load_mesh!(mesh::SerialTreeMesh, restart_file::AbstractString)
  # Determine mesh filename
  filename = get_restart_mesh_filename(restart_file, mpi_parallel(mesh))
  mesh.current_filename = filename
  mesh.unsaved_changes = false

  # Read mesh file
  h5open(filename, "r") do file
    # Set domain information
    mesh.tree.center_level_0 = read(attributes(file)["center_level_0"])
    mesh.tree.length_level_0 = read(attributes(file)["length_level_0"])
    mesh.tree.periodicity    = Tuple(read(attributes(file)["periodicity"]))

    # Set length
    n_cells = read(attributes(file)["n_cells"])
    resize!(mesh.tree, n_cells)

    # Read in data
    mesh.tree.parent_ids[1:n_cells] = read(file["parent_ids"])
    mesh.tree.child_ids[:, 1:n_cells] = read(file["child_ids"])
    mesh.tree.neighbor_ids[:, 1:n_cells] = read(file["neighbor_ids"])
    mesh.tree.levels[1:n_cells] = read(file["levels"])
    mesh.tree.coordinates[:, 1:n_cells] = read(file["coordinates"])
  end

  return mesh
end

function load_mesh_parallel(restart_file::AbstractString; n_cells_max)
  if mpi_isroot()
    ndims_ = h5open(restart_file, "r") do file
      read(attributes(file)["ndims"])
    end
    MPI.Bcast!(Ref(ndims_), mpi_root(), mpi_comm())
  else
    ndims_ = MPI.Bcast!(Ref(0), mpi_root(), mpi_comm())[]
  end

  mesh = TreeMesh(ParallelTree{ndims_}, n_cells_max)
  load_mesh!(mesh, restart_file)
end

function load_mesh!(mesh::ParallelTreeMesh, restart_file::AbstractString)
  # Determine mesh filename
  filename = get_restart_mesh_filename(restart_file, mpi_parallel(mesh))
  mesh.current_filename = filename
  mesh.unsaved_changes = false

  if mpi_isroot()
    h5open(filename, "r") do file
      # Set domain information
      mesh.tree.center_level_0 = read(attributes(file)["center_level_0"])
      mesh.tree.length_level_0 = read(attributes(file)["length_level_0"])
      mesh.tree.periodicity    = Tuple(read(attributes(file)["periodicity"]))
      MPI.Bcast!(collect(mesh.tree.center_level_0), mpi_root(), mpi_comm())
      MPI.Bcast!(collect(mesh.tree.length_level_0), mpi_root(), mpi_comm())
      MPI.Bcast!(collect(mesh.tree.periodicity),    mpi_root(), mpi_comm())

      # Set length
      n_cells = read(attributes(file)["n_cells"])
      MPI.Bcast!(Ref(n_cells), mpi_root(), mpi_comm())
      resize!(mesh.tree, n_cells)

      # Read in data
      mesh.tree.parent_ids[1:n_cells] = read(file["parent_ids"])
      mesh.tree.child_ids[:, 1:n_cells] = read(file["child_ids"])
      mesh.tree.neighbor_ids[:, 1:n_cells] = read(file["neighbor_ids"])
      mesh.tree.levels[1:n_cells] = read(file["levels"])
      mesh.tree.coordinates[:, 1:n_cells] = read(file["coordinates"])
      @views MPI.Bcast!(mesh.tree.parent_ids[1:n_cells],      mpi_root(), mpi_comm())
      @views MPI.Bcast!(mesh.tree.child_ids[:, 1:n_cells],    mpi_root(), mpi_comm())
      @views MPI.Bcast!(mesh.tree.neighbor_ids[:, 1:n_cells], mpi_root(), mpi_comm())
      @views MPI.Bcast!(mesh.tree.levels[1:n_cells],          mpi_root(), mpi_comm())
      @views MPI.Bcast!(mesh.tree.coordinates[:, 1:n_cells],  mpi_root(), mpi_comm())
    end
  else # non-root ranks
    # Set domain information
    mesh.tree.center_level_0 = MPI.Bcast!(collect(mesh.tree.center_level_0), mpi_root(), mpi_comm())
    mesh.tree.length_level_0 = MPI.Bcast!(collect(mesh.tree.length_level_0), mpi_root(), mpi_comm())[1]
    mesh.tree.periodicity    = Tuple(MPI.Bcast!(collect(mesh.tree.periodicity),    mpi_root(), mpi_comm()))

    # Set length
    n_cells = MPI.Bcast!(Ref(0), mpi_root(), mpi_comm())[]
    resize!(mesh.tree, n_cells)

    # Read in data
    @views MPI.Bcast!(mesh.tree.parent_ids[1:n_cells],      mpi_root(), mpi_comm())
    @views MPI.Bcast!(mesh.tree.child_ids[:, 1:n_cells],    mpi_root(), mpi_comm())
    @views MPI.Bcast!(mesh.tree.neighbor_ids[:, 1:n_cells], mpi_root(), mpi_comm())
    @views MPI.Bcast!(mesh.tree.levels[1:n_cells],          mpi_root(), mpi_comm())
    @views MPI.Bcast!(mesh.tree.coordinates[:, 1:n_cells],  mpi_root(), mpi_comm())
  end

  # Partition mesh
  partition!(mesh)

  return mesh
end


text:
module TestCachingOptimizer

using Test

import MathOptInterface

const MOI = MathOptInterface
const MOIU = MOI.Utilities

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

###
### Helper structs
###

struct DummyModelAttribute <: MOI.AbstractModelAttribute end

struct DummyEvaluator <: MOI.AbstractNLPEvaluator end

struct DummyVariableAttribute <: MOI.AbstractVariableAttribute end

struct DummyConstraintAttribute <: MOI.AbstractConstraintAttribute end

mutable struct NoFreeVariables <: MOI.AbstractOptimizer
    inner::MOIU.Model{Float64}
    NoFreeVariables() = new(MOIU.Model{Float64}())
end

MOI.is_empty(model::NoFreeVariables) = MOI.is_empty(model.inner)

MOI.empty!(model::NoFreeVariables) = MOI.empty!(model.inner)

function MOI.get(model::NoFreeVariables, attr::MOI.AnyAttribute, idx::Vector)
    return MOI.get(model.inner, attr, idx)
end

function MOI.get(model::NoFreeVariables, attr::MOI.AnyAttribute, args...)
    return MOI.get(model.inner, attr, args...)
end

function MOI.supports_add_constrained_variables(
    ::NoFreeVariables,
    ::Type{MOI.Reals},
)
    return false
end

function MOI.supports_add_constrained_variable(
    ::NoFreeVariables,
    ::Type{<:MOI.AbstractScalarSet},
)
    return true
end

function MOI.add_constrained_variable(
    model::NoFreeVariables,
    set::MOI.AbstractScalarSet,
)
    return MOI.add_constrained_variable(model.inner, set)
end

function MOI.supports_add_constrained_variables(
    ::NoFreeVariables,
    ::Type{<:MOI.AbstractVectorSet},
)
    return true
end

function MOI.add_constrained_variables(
    model::NoFreeVariables,
    set::MOI.AbstractVectorSet,
)
    return MOI.add_constrained_variables(model.inner, set)
end

MOI.supports_incremental_interface(::NoFreeVariables) = true

function MOI.copy_to(dest::NoFreeVariables, src::MOI.ModelLike; kwargs...)
    return MOI.Utilities.default_copy_to(dest, src; kwargs...)
end

###
### The tests
###

function test_MOI_Test()
    # It seems like these loops might take a while. But the first one takes
    # _forever_, and then the rest are really quick (like 140s vs 0.4s).
    for state in
        (MOIU.NO_OPTIMIZER, MOIU.EMPTY_OPTIMIZER, MOIU.ATTACHED_OPTIMIZER)
        for mode in (MOIU.MANUAL, MOIU.AUTOMATIC)
            model = MOIU.CachingOptimizer(MOIU.Model{Float64}(), mode)
            if state != MOIU.NO_OPTIMIZER
                optimizer = MOIU.MockOptimizer(
                    MOIU.Model{Float64}(),
                    supports_names = false,
                )
                MOIU.reset_optimizer(model, optimizer)
                if state == MOIU.ATTACHED_OPTIMIZER
                    MOIU.attach_optimizer(model)
                end
            end
            @test MOIU.state(model) == state
            @test MOIU.mode(model) == mode
            MOI.Test.runtests(
                model,
                MOI.Test.Config(exclude = Any[MOI.optimize!]),
                exclude = [
                    "test_attribute_SolverName",
                    "test_attribute_SolverVersion",
                ],
            )
        end
    end
    return
end

# !!! warning
#     This is some type piracy! To enable CachingOptimizer to pass some MOI.Test
#     functions with MockOptimizer, we overload `setup_test` to setup the inner
#     mock optimizer.
#
#     This is pretty fragile! It requires the inner optimizer to be attached,
#     amongst other things. It's used by `test_compute_conflict` below, but not
#     by test_MOI_Test above (test_MOI_Test has `exclude = Any[MOI.optimize!]`).
function MOI.Test.setup_test(
    f::Any,
    model::MOI.Utilities.CachingOptimizer{
        MOI.Utilities.MockOptimizer{
            MOI.Utilities.UniversalFallback{MOI.Utilities.Model{Float64}},
        },
        MOI.Utilities.Model{Float64},
    },
    config::MOI.Test.Config,
)
    MOI.Test.setup_test(f, model.optimizer, config)
    return
end

function test_compute_conflict()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    model = MOI.Utilities.CachingOptimizer(MOI.Utilities.Model{Float64}(), mock)
    MOI.Test.runtests(
        model,
        MOI.Test.Config(),
        include = ["test_solve_conflict"],
    )
    return
end

"""
Without an optimizer attached (i.e., `MOI.state(model) == NO_OPTIMIZER`) we
need throw nice errors for attributes that are based on the optimizer. For
`AbstractModelAttribute`s that `is_set_by_optimize` returns `true` for, we
overload `TerminationStatus`, `PrimalStatus`, or `DualStatus` to return
sane default values. Otherwise we throw a nice error.
"""
function test_default_attributes()
    model = MOIU.CachingOptimizer(MOIU.Model{Float64}(), MOIU.MANUAL)
    @test MOIU.state(model) == MOIU.NO_OPTIMIZER

    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.NO_SOLUTION
    @test MOI.get(model, MOI.DualStatus()) == MOI.NO_SOLUTION
    x = MOI.add_variables(model, 2)
    attr = MOI.VariablePrimal()
    exception = ErrorException(
        "Cannot query $(attr) from caching optimizer because no optimizer" *
        " is attached.",
    )
    @test_throws exception MOI.get(model, MOI.VariablePrimal(), x[1])
    @test_throws exception MOI.get(model, MOI.VariablePrimal(), x)
    for attr in (
        MOI.SolverName(),
        MOI.Silent(),
        MOI.TimeLimitSec(),
        MOI.NumberOfThreads(),
        MOI.ResultCount(),
    )
        @test_throws(
            ErrorException(
                "Cannot query $(attr) from caching optimizer because no " *
                "optimizer is attached.",
            ),
            MOI.get(model, attr),
        )
    end
    return
end

function test_copyable_solver_attributes()
    cache = MOIU.UniversalFallback(MOIU.Model{Float64}())
    cached = MOIU.CachingOptimizer(cache, MOIU.MANUAL)
    MOI.set(cached, MOI.Silent(), true)
    MOI.set(cached, MOI.TimeLimitSec(), 0.0)
    MOI.set(cached, MOI.NumberOfThreads(), 1)
    mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()))
    MOIU.reset_optimizer(cached, mock)
    @test MOI.get(mock, MOI.Silent())
    @test MOI.get(cached, MOI.Silent())
    @test MOI.get(mock, MOI.TimeLimitSec()) == 0.0
    @test MOI.get(cached, MOI.TimeLimitSec()) == 0.0
    @test MOI.get(mock, MOI.NumberOfThreads()) == 1
    @test MOI.get(cached, MOI.NumberOfThreads()) == 1
    MOI.set(cached, MOI.Silent(), false)
    MOI.set(cached, MOI.TimeLimitSec(), 1.0)
    MOI.set(cached, MOI.NumberOfThreads(), 2)
    @test !MOI.get(mock, MOI.Silent())
    @test !MOI.get(cached, MOI.Silent())
    @test MOI.get(mock, MOI.TimeLimitSec()) ≈ 1.0
    @test MOI.get(cached, MOI.TimeLimitSec()) ≈ 1.0
    @test MOI.get(mock, MOI.NumberOfThreads()) == 2
    @test MOI.get(cached, MOI.NumberOfThreads()) == 2
    mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()))
    MOIU.reset_optimizer(cached, mock)
    @test !MOI.get(mock, MOI.Silent())
    @test !MOI.get(cached, MOI.Silent())
    @test MOI.get(mock, MOI.TimeLimitSec()) ≈ 1.0
    @test MOI.get(cached, MOI.TimeLimitSec()) ≈ 1.0
    @test MOI.get(mock, MOI.NumberOfThreads()) == 2
    @test MOI.get(cached, MOI.NumberOfThreads()) == 2
    MOI.set(cached, MOI.Silent(), true)
    MOI.set(cached, MOI.TimeLimitSec(), 0.0)
    MOI.set(cached, MOI.NumberOfThreads(), 1)
    @test MOI.get(mock, MOI.Silent())
    @test MOI.get(cached, MOI.Silent())
    @test MOI.get(mock, MOI.TimeLimitSec()) == 0.0
    @test MOI.get(cached, MOI.TimeLimitSec()) == 0.0
    @test MOI.get(mock, MOI.NumberOfThreads()) == 1
    @test MOI.get(cached, MOI.NumberOfThreads()) == 1
    mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()))
    MOIU.reset_optimizer(cached, mock)
    @test MOI.get(mock, MOI.Silent())
    @test MOI.get(cached, MOI.Silent())
    @test MOI.get(mock, MOI.TimeLimitSec()) == 0.0
    @test MOI.get(cached, MOI.TimeLimitSec()) == 0.0
    @test MOI.get(mock, MOI.NumberOfThreads()) == 1
    @test MOI.get(cached, MOI.NumberOfThreads()) == 1
    return
end

function test_mapping_of_variables()
    mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()))
    model = MOIU.CachingOptimizer(
        MOIU.UniversalFallback(MOIU.Model{Float64}()),
        mock,
    )
    MOIU.attach_optimizer(model)
    x = MOI.add_variable(model)
    y = first(MOI.get(mock, MOI.ListOfVariableIndices()))
    @test x != y # Otherwise, these tests will trivially pass
    @test !MOI.is_valid(model, y)
    @test !MOI.is_valid(mock, x)
    cx = MOI.add_constraint(model, x, MOI.GreaterThan(1.0))
    cy = first(
        MOI.get(
            mock,
            MOI.ListOfConstraintIndices{
                MOI.VariableIndex,
                MOI.GreaterThan{Float64},
            }(),
        ),
    )
    @test !MOI.is_valid(model, cy)
    @test !MOI.is_valid(mock, cx)
    @test MOI.get(mock, MOI.ConstraintFunction(), cy) == y

    c2x = MOI.add_constraint(model, 2.0x, MOI.GreaterThan(1.0))
    c2y = first(
        MOI.get(
            mock,
            MOI.ListOfConstraintIndices{
                MOI.ScalarAffineFunction{Float64},
                MOI.GreaterThan{Float64},
            }(),
        ),
    )
    @test !MOI.is_valid(model, c2y)
    @test !MOI.is_valid(mock, c2x)
    @test MOI.get(model, MOI.ConstraintFunction(), c2x) ≈ 2.0x
    @test MOI.get(model, MOI.ConstraintSet(), c2x) == MOI.GreaterThan(1.0)
    @test MOI.get(mock, MOI.ConstraintFunction(), c2y) ≈ 2.0y
    @test MOI.get(mock, MOI.ConstraintSet(), c2y) == MOI.GreaterThan(1.0)

    MOI.set(model, MOI.ConstraintSet(), c2x, MOI.GreaterThan(2.0))
    @test MOI.get(model, MOI.ConstraintFunction(), c2x) ≈ 2.0x
    @test MOI.get(model, MOI.ConstraintSet(), c2x) == MOI.GreaterThan(2.0)
    @test MOI.get(mock, MOI.ConstraintFunction(), c2y) ≈ 2.0y
    @test MOI.get(mock, MOI.ConstraintSet(), c2y) == MOI.GreaterThan(2.0)

    MOI.set(model, MOI.ConstraintFunction(), c2x, 3.0x)
    @test MOI.get(model, MOI.ConstraintFunction(), c2x) ≈ 3.0x
    @test MOI.get(model, MOI.ConstraintSet(), c2x) == MOI.GreaterThan(2.0)
    @test MOI.get(mock, MOI.ConstraintFunction(), c2y) ≈ 3.0y
    @test MOI.get(mock, MOI.ConstraintSet(), c2y) == MOI.GreaterThan(2.0)

    MOI.set(model, MOI.ConstraintSet(), c2x, MOI.GreaterThan(4.0))
    @test MOI.get(model, MOI.ConstraintFunction(), c2x) ≈ 3.0x
    @test MOI.get(model, MOI.ConstraintSet(), c2x) == MOI.GreaterThan(4.0)
    @test MOI.get(mock, MOI.ConstraintFunction(), c2y) ≈ 3.0y
    @test MOI.get(mock, MOI.ConstraintSet(), c2y) == MOI.GreaterThan(4.0)

    MOI.set(model, MOI.ObjectiveFunction{typeof(x)}(), x)
    @test MOI.get(mock, MOI.ObjectiveFunction{typeof(y)}()) == y

    MOI.set(model, MOI.ObjectiveFunction{typeof(2.0x)}(), 2.0x)
    @test MOI.get(mock, MOI.ObjectiveFunction{typeof(2.0y)}()) ≈ 2.0y

    MOI.set(model, DummyModelAttribute(), 2.0x + 1.0)
    @test MOI.get(model, DummyModelAttribute()) ≈ 2.0x + 1.0
    @test MOI.get(mock, DummyModelAttribute()) ≈ 2.0y + 1.0

    for (attr, cache_index, optimizer_index) in
        [(DummyVariableAttribute(), x, y), (DummyConstraintAttribute(), cx, cy)]
        MOI.set(model, attr, cache_index, 1.0x)
        @test MOI.get(model, attr, cache_index) ≈ 1.0x
        @test MOI.get(mock, attr, optimizer_index) ≈ 1.0y

        MOI.set(model, attr, [cache_index], [3.0x])
        @test MOI.get(model, attr, [cache_index])[1] ≈ 3.0x
        @test MOI.get(mock, attr, [optimizer_index])[1] ≈ 3.0y
    end

    @testset "RawSolver" begin
        MOI.get(model, MOI.RawSolver()) === mock
    end

    @testset "HeuristicCallback" begin
        attr = MOI.HeuristicCallback()
        f(callback_data) = nothing
        MOI.set(model, attr, f)
        @test MOI.get(model, attr) === f
    end

    @testset "CallbackVariablePrimal" begin
        attr = MOI.CallbackVariablePrimal(nothing)
        err =
            ErrorException("No mock callback primal is set for variable `$y`.")
        @test_throws err MOI.get(model, attr, x)
        MOI.set(mock, attr, y, 1.0)
        @test_throws MOI.InvalidIndex(x) MOI.get(mock, attr, x)
        @test MOI.get(mock, attr, y) == 1.0
        @test MOI.get(model, attr, x) == 1.0
    end

    @testset "LazyConstraint" begin
        sub = MOI.LazyConstraint(nothing)
        @test MOI.supports(model, sub)
        MOI.submit(model, sub, 2.0x, MOI.GreaterThan(1.0))
        @test mock.submitted[sub][1][1] ≈ 2.0y
        @test mock.submitted[sub][1][2] == MOI.GreaterThan(1.0)
    end
    @testset "HeuristicSolution" begin
        sub = MOI.HeuristicSolution(nothing)
        @test MOI.supports(model, sub)
        MOI.submit(model, sub, [x], [1.0])
        @test mock.submitted[sub][1][1] == [y]
        @test mock.submitted[sub][1][2] == [1.0]
    end

    nlp_data =
        MOI.NLPBlockData([MOI.NLPBoundsPair(1.0, 2.0)], DummyEvaluator(), false)
    MOI.set(model, MOI.NLPBlock(), nlp_data)
    for nlp_data in
        [MOI.get(model, MOI.NLPBlock()), MOI.get(mock, MOI.NLPBlock())]
        @test nlp_data.constraint_bounds == [MOI.NLPBoundsPair(1.0, 2.0)]
        @test nlp_data.evaluator isa DummyEvaluator
        @test nlp_data.has_objective == false
    end
    return
end

function test_CachingOptimizer_MANUAL_mode()
    m = MOIU.CachingOptimizer(MOIU.Model{Float64}(), MOIU.MANUAL)
    @test MOIU.state(m) == MOIU.NO_OPTIMIZER

    s = MOIU.MockOptimizer(MOIU.Model{Float64}(), supports_names = false)
    @test MOI.is_empty(s)
    MOIU.reset_optimizer(m, s)
    @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER

    v = MOI.add_variable(m)
    x = MOI.add_variables(m, 2)
    saf = MOI.ScalarAffineFunction(
        MOI.ScalarAffineTerm.([1.0, 2.0, 3.0], [v; x]),
        0.0,
    )
    @test MOI.supports(m, MOI.ObjectiveFunction{typeof(saf)}())
    MOI.set(m, MOI.ObjectiveFunction{typeof(saf)}(), saf)
    @test MOI.get(
        m,
        MOIU.AttributeFromModelCache(MOI.ObjectiveFunction{typeof(saf)}()),
    ) ≈ saf
    @test MOI.get(m, MOI.ObjectiveFunction{typeof(saf)}()) ≈ saf

    @test_throws AssertionError MOI.optimize!(m)

    MOIU.attach_optimizer(m)
    @test MOIU.state(m) == MOIU.ATTACHED_OPTIMIZER
    @test MOI.get(
        m,
        MOIU.AttributeFromOptimizer(MOI.ObjectiveFunction{typeof(saf)}()),
    ) ≈ saf

    @test MOI.supports(m, MOI.ObjectiveSense())
    MOI.set(m, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    @test MOI.get(m, MOI.ObjectiveSense()) == MOI.MAX_SENSE
    @test MOI.get(m, MOIU.AttributeFromModelCache(MOI.ObjectiveSense())) ==
          MOI.MAX_SENSE
    @test MOI.get(m, MOIU.AttributeFromOptimizer(MOI.ObjectiveSense())) ==
          MOI.MAX_SENSE

    @test MOI.supports(
        m,
        MOIU.AttributeFromOptimizer(MOIU.MockModelAttribute()),
    )
    MOI.set(m, MOIU.AttributeFromOptimizer(MOIU.MockModelAttribute()), 10)
    @test MOI.get(m, MOIU.AttributeFromOptimizer(MOIU.MockModelAttribute())) ==
          10

    MOI.set(m, MOIU.AttributeFromOptimizer(MOI.ResultCount()), 1)
    @test MOI.supports(
        m,
        MOIU.AttributeFromOptimizer(MOI.VariablePrimal()),
        typeof(v),
    )
    MOI.set(m, MOIU.AttributeFromOptimizer(MOI.VariablePrimal()), v, 3.0)

    MOI.optimize!(m)

    @test MOI.get(m, MOI.VariablePrimal(), v) == 3.0
    @test MOI.get(m, MOI.VariablePrimal(), [v]) == [3.0]
    @test MOI.get(m, MOIU.AttributeFromOptimizer(MOI.VariablePrimal()), v) ==
          3.0

    @test MOI.supports_constraint(
        m.model_cache,
        MOI.VariableIndex,
        MOI.LessThan{Float64},
    )
    @test MOI.supports_constraint(
        m.optimizer.inner_model,
        MOI.VariableIndex,
        MOI.LessThan{Float64},
    )
    @test MOI.supports_constraint(
        m.optimizer,
        MOI.VariableIndex,
        MOI.LessThan{Float64},
    )
    @test MOI.supports_constraint(m, MOI.VariableIndex, MOI.LessThan{Float64})
    lb = MOI.add_constraint(m, v, MOI.LessThan(10.0))
    MOI.set(m, MOI.ConstraintSet(), lb, MOI.LessThan(11.0))
    @test MOI.get(m, MOI.ConstraintSet(), lb) == MOI.LessThan(11.0)
    @test MOI.get(m, MOI.ConstraintFunction(), lb) == v

    MOIU.drop_optimizer(m)
    @test MOIU.state(m) == MOIU.NO_OPTIMIZER

    MOI.set(m, MOI.ConstraintSet(), lb, MOI.LessThan(12.0))
    @test MOI.get(m, MOI.ConstraintSet(), lb) == MOI.LessThan(12.0)

    MOI.delete(m, x[2])
    @test_throws MOI.InvalidIndex{typeof(x[2])} MOI.delete(m, x[2])
    @test !MOI.is_valid(m, x[2])

    # TODO: test more constraint modifications

    @test sprint(show, m) == MOI.Utilities.replace_acronym("""
    $(MOIU.CachingOptimizer{MOI.AbstractOptimizer,MOIU.Model{Float64}})
    in state NO_OPTIMIZER
    in mode MANUAL
    with model cache $(MOIU.Model{Float64})
    with optimizer nothing""")

    MOI.empty!(s)
    MOIU.reset_optimizer(m, s)
    @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER
    MOIU.attach_optimizer(m)
    @test MOIU.state(m) == MOIU.ATTACHED_OPTIMIZER
    MOI.empty!(m)
    @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER
end

function test_CachingOptimizer_AUTOMATIC_mode()
    m = MOIU.CachingOptimizer(MOIU.Model{Float64}(), MOIU.AUTOMATIC)
    @test MOIU.state(m) == MOIU.NO_OPTIMIZER

    v = MOI.add_variable(m)
    @test MOI.supports(m, MOI.VariableName(), typeof(v))
    MOI.set(m, MOI.VariableName(), v, "v")
    @test MOI.get(m, MOI.VariableName(), v) == "v"

    s = MOIU.MockOptimizer(MOIU.Model{Float64}(), supports_names = false)
    @test MOI.is_empty(s)
    MOIU.reset_optimizer(m, s)
    @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER

    saf = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, v)], 0.0)
    @test MOI.supports(m, MOI.ObjectiveFunction{typeof(saf)}())
    MOI.set(m, MOI.ObjectiveFunction{typeof(saf)}(), saf)
    @test MOI.get(
        m,
        MOIU.AttributeFromModelCache(MOI.ObjectiveFunction{typeof(saf)}()),
    ) ≈ saf

    MOI.optimize!(m)
    @test MOIU.state(m) == MOIU.ATTACHED_OPTIMIZER
    @test MOI.get(m, MOIU.AttributeFromOptimizer(MOI.TerminationStatus())) ==
          MOI.OPTIMIZE_NOT_CALLED

    @test MOI.get(m, MOI.VariableName(), v) == "v"
    @test MOI.get(m, MOIU.AttributeFromModelCache(MOI.VariableName()), v) == "v"
    # The caching optimizer does not copy names
    @test MOI.get(m, MOIU.AttributeFromOptimizer(MOI.VariableName()), v) == ""

    @test MOI.supports(m, MOI.ObjectiveSense())
    MOI.set(m, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    @test MOI.get(m, MOIU.AttributeFromModelCache(MOI.ObjectiveSense())) ==
          MOI.MAX_SENSE
    @test MOI.get(m, MOIU.AttributeFromOptimizer(MOI.ObjectiveSense())) ==
          MOI.MAX_SENSE

    @test MOI.supports(
        m,
        MOIU.AttributeFromOptimizer(MOIU.MockModelAttribute()),
    )
    MOI.set(m, MOIU.AttributeFromOptimizer(MOIU.MockModelAttribute()), 10)
    @test MOI.get(m, MOIU.AttributeFromOptimizer(MOIU.MockModelAttribute())) ==
          10

    MOI.set(m, MOIU.AttributeFromOptimizer(MOI.ResultCount()), 1)
    @test MOI.supports(
        m,
        MOIU.AttributeFromOptimizer(MOI.VariablePrimal()),
        typeof(v),
    )
    MOI.set(m, MOIU.AttributeFromOptimizer(MOI.VariablePrimal()), v, 3.0)

    MOI.optimize!(m)

    @test MOI.get(m, MOIU.AttributeFromOptimizer(MOI.VariablePrimal()), v) ==
          3.0

    @testset "Modify not allowed" begin
        s.modify_allowed = false
        MOI.modify(
            m,
            MOI.ObjectiveFunction{typeof(saf)}(),
            MOI.ScalarConstantChange(1.0),
        )
        s.modify_allowed = true
        @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER
        MOIU.attach_optimizer(m)
        @test MOIU.state(m) == MOIU.ATTACHED_OPTIMIZER
        ci = MOI.add_constraint(m, saf, MOI.EqualTo(0.0))
        s.modify_allowed = false
        MOI.modify(m, ci, MOI.ScalarCoefficientChange(v, 1.0))
        s.modify_allowed = true
        @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER
        MOIU.attach_optimizer(m)
        @test MOIU.state(m) == MOIU.ATTACHED_OPTIMIZER
    end

    # Simulate that constraints cannot be added
    @testset "Add constraint not allowed" begin
        s.add_con_allowed = false
        MOI.add_constraint(
            m,
            MOI.VectorOfVariables([v]),
            MOI.SecondOrderCone(1),
        )
        s.add_con_allowed = true
        @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER
        MOI.empty!(m)
        @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER
    end

    @testset "Add variable not allowed" begin
        s.add_var_allowed = false # Simulate optimizer that cannot add variables incrementally
        MOI.add_variable(m)
        @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER
        @test_throws MOI.AddVariableNotAllowed MOIU.attach_optimizer(m)
        @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER

        s.add_var_allowed = true
        MOIU.attach_optimizer(m)
        @test MOIU.state(m) == MOIU.ATTACHED_OPTIMIZER
        s.add_var_allowed = false
        MOI.add_variables(m, 2)
        s.add_var_allowed = true
        @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER
    end

    @testset "Delete not allowed" begin
        vi = MOI.add_variable(m)
        s.delete_allowed = false # Simulate optimizer that cannot delete variable
        MOI.delete(m, vi)
        s.delete_allowed = true
        @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER
        MOIU.attach_optimizer(m)
        @test MOIU.state(m) == MOIU.ATTACHED_OPTIMIZER

        vi = MOI.add_variable(m)
        ci = MOI.add_constraint(m, vi, MOI.EqualTo(0.0))
        s.delete_allowed = false # Simulate optimizer that cannot delete constraint
        MOI.delete(m, ci)
        s.delete_allowed = true
        @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER
        MOIU.attach_optimizer(m)
        @test MOIU.state(m) == MOIU.ATTACHED_OPTIMIZER
    end

    @test sprint(show, m) == MOI.Utilities.replace_acronym("""
    $(MOIU.CachingOptimizer{MOI.AbstractOptimizer,MOIU.Model{Float64}})
    in state ATTACHED_OPTIMIZER
    in mode AUTOMATIC
    with model cache $(MOIU.Model{Float64})
    with optimizer $(MOIU.MockOptimizer{MOIU.Model{Float64}})""")

    MOI.empty!(m)
    @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER
    return
end

function test_empty_model_and_optimizer()
    s = MOIU.MockOptimizer(MOIU.Model{Float64}(), supports_names = false)
    model = MOIU.Model{Float64}()
    m = MOIU.CachingOptimizer(model, s)
    @test m isa MOIU.CachingOptimizer{typeof(s),typeof(model)}
    @test MOI.is_empty(m)
    @test MOIU.state(m) == MOIU.EMPTY_OPTIMIZER
    @test MOIU.mode(m) == MOIU.AUTOMATIC
    @test MOI.get(m, MOI.SolverName()) == "Mock"
    @test sprint(show, m) == MOI.Utilities.replace_acronym("""
    $(MOIU.CachingOptimizer{MOIU.MockOptimizer{MOIU.Model{Float64}},MOIU.Model{Float64}})
    in state EMPTY_OPTIMIZER
    in mode AUTOMATIC
    with model cache $(MOIU.Model{Float64})
    with optimizer $(MOIU.MockOptimizer{MOIU.Model{Float64}})""")
end

function test_empty_model_nonempty_optimizer()
    s = MOIU.MockOptimizer(MOIU.Model{Float64}(), supports_names = false)
    MOI.add_variable(s)
    model = MOIU.Model{Float64}()
    @test MOI.is_empty(model)
    @test !MOI.is_empty(s)
    @test_throws AssertionError MOIU.CachingOptimizer(model, s)
    return
end

function test_nonempty_model()
    s = MOIU.MockOptimizer(MOIU.Model{Float64}(), supports_names = false)
    model = MOIU.Model{Float64}()
    MOI.add_variable(model)
    @test !MOI.is_empty(model)
    @test MOI.is_empty(s)
    @test MOIU.CachingOptimizer(model, s) isa MOIU.CachingOptimizer
    return
end

function _constrained_variables_test(model)
    @test !MOI.supports_add_constrained_variables(model, MOI.Reals)
    @test MOI.supports_add_constrained_variable(model, MOI.ZeroOne)
    @test !MOI.supports_constraint(model, MOI.VariableIndex, MOI.ZeroOne)
    @test MOI.supports_add_constrained_variables(model, MOI.Nonnegatives)
    @test !MOI.supports_constraint(
        model,
        MOI.VectorOfVariables,
        MOI.Nonnegatives,
    )
    scalar_set = MOI.ZeroOne()
    x, cx = MOI.add_constrained_variable(model, scalar_set)
    vector_set = MOI.Nonnegatives(2)
    y, cy = MOI.add_constrained_variables(model, vector_set)
    constraint_types = Set([
        (MOI.VariableIndex, MOI.ZeroOne),
        (MOI.VectorOfVariables, MOI.Nonnegatives),
    ])
    @test Set(MOI.get(model.model_cache, MOI.ListOfConstraintTypesPresent())) ==
          constraint_types
    if MOIU.state(model) == MOIU.EMPTY_OPTIMIZER
        MOIU.attach_optimizer(model)
    end
    @test Set(MOI.get(model.optimizer, MOI.ListOfConstraintTypesPresent())) ==
          constraint_types
    return
end

function test_constrained_variables()
    cache = NoFreeVariables()
    optimizer = NoFreeVariables()
    model = MOIU.CachingOptimizer(cache, optimizer)
    _constrained_variables_test(model)
    return
end

function test_constrained_variables_AUTOMATIC()
    cache = NoFreeVariables()
    optimizer = NoFreeVariables()
    model = MOIU.CachingOptimizer(cache, MOIU.AUTOMATIC)
    MOIU.reset_optimizer(model, optimizer)
    _constrained_variables_test(model)
    return
end

struct Issue1220 <: MOI.AbstractOptimizer
    optimizer_attributes::Dict{Any,Any}
    Issue1220() = new(Dict{Any,Any}())
end

MOI.is_empty(model::Issue1220) = isempty(model.optimizer_attributes)

function MOI.get(model::Issue1220, ::MOI.ListOfOptimizerAttributesSet)
    return collect(keys(model.optimizer_attributes))
end

MOI.supports(::Issue1220, ::MOI.AbstractOptimizerAttribute) = true

MOI.supports(::Issue1220, ::MOI.NumberOfThreads) = false

function MOI.get(model::Issue1220, attr::MOI.AbstractOptimizerAttribute)
    return model.optimizer_attributes[attr]
end

function MOI.set(model::Issue1220, attr::MOI.AbstractOptimizerAttribute, value)
    model.optimizer_attributes[attr] = value
    return
end

function test_issue1220_dont_pass_raw_parameter()
    model = MOIU.CachingOptimizer(Issue1220(), Issue1220())
    MOI.set(model, MOI.Silent(), true)
    MOI.set(model, MOI.RawOptimizerAttribute("foo"), "bar")
    MOI.set(model, MOI.NumberOfThreads(), 1)
    MOIU.reset_optimizer(model, Issue1220())
    @test MOI.get(model, MOI.Silent()) == true
    @test_throws KeyError MOI.get(model, MOI.RawOptimizerAttribute("foo"))
    @test_throws KeyError MOI.get(model, MOI.NumberOfThreads())
    return
end

function test_status_codes()
    optimizer = MOI.Utilities.MockOptimizer(MOI.Utilities.Model{Float64}())
    model = MOI.Utilities.CachingOptimizer(
        MOI.Utilities.Model{Float64}(),
        optimizer,
    )
    MOI.Utilities.attach_optimizer(model)
    MOI.set(optimizer, MOI.TerminationStatus(), MOI.OPTIMAL)
    MOI.set(optimizer, MOI.PrimalStatus(), MOI.FEASIBLE_POINT)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    return
end

mutable struct CopyToAndOptimizer <: MOI.AbstractOptimizer
    is_dirty::Bool
end

function MOI.empty!(x::CopyToAndOptimizer)
    x.is_dirty = false
    return
end

MOI.is_empty(x::CopyToAndOptimizer) = !x.is_dirty

function MOI.optimize!(x::CopyToAndOptimizer, ::MOI.ModelLike)
    x.is_dirty = true
    return MOI.Utilities.IndexMap(), false
end

function test_copy_to_and_optimize!()
    optimizer = CopyToAndOptimizer(false)
    model = MOI.Utilities.CachingOptimizer(
        MOI.Utilities.Model{Float64}(),
        optimizer,
    )
    MOI.optimize!(model)
    @test MOI.Utilities.state(model) == MOI.Utilities.EMPTY_OPTIMIZER
    @test !MOI.is_empty(optimizer)
    MOI.empty!(model)
    @test MOI.is_empty(optimizer)
    @test MOI.Utilities.state(model) == MOI.Utilities.EMPTY_OPTIMIZER
    return
end

end  # module

TestCachingOptimizer.runtests()


text:
module PkgPractice

export plusTwo

plusTwo(x) = return x+2

end # module

max_stars_repo_name: mcmory/PkgPractice

max_issues_repo_name: mcmory/PkgPractice

max_forks_repo_name: mcmory/PkgPractice


text:
@test [combinations([])...] == []
@test [combinations(['a', 'b', 'c'])...] == Any[['a'],['b'],['c'],['a','b'],['a','c'],['b','c'],['a','b','c']]

@test [combinations("abc",3)...] == Any[['a','b','c']]
@test [combinations("abc",2)...] == Any[['a','b'],['a','c'],['b','c']]
@test [combinations("abc",1)...] == Any[['a'],['b'],['c']]
@test [combinations("abc",0)...] == Any[[]]
@test [combinations("abc",-1)...] == []

@test filter(x->iseven(x[1]),[combinations([1,2,3],2)...]) == Any[[2,3]]

# multiset_combinations
@test [multiset_combinations("aabc", 5)...] == Any[]
@test [multiset_combinations("aabc", 2)...] == Any[['a','a'],['a','b'],['a','c'],['b','c']]
@test [multiset_combinations("aabc", 1)...] == Any[['a'],['b'],['c']]
@test [multiset_combinations("aabc", 0)...] == Any[Char[]]
@test [multiset_combinations("aabc", -1)...] == Any[]
@test [multiset_combinations("", 1)...] == Any[]
@test [multiset_combinations("", 0)...] == Any[Char[]]
@test [multiset_combinations("", -1)...] == Any[]

# with_replacement_combinations
@test [with_replacement_combinations("abc", 2)...] == Any[['a','a'],['a','b'],['a','c'],
                                                              ['b','b'],['b','c'],['c','c']]
@test [with_replacement_combinations("abc", 1)...] == Any[['a'],['b'],['c']]
@test [with_replacement_combinations("abc", 0)...] == Any[Char[]]
@test [with_replacement_combinations("abc", -1)...] == Any[]
@test [with_replacement_combinations("", 1)...] == Any[]
@test [with_replacement_combinations("", 0)...] == Any[Char[]]
@test [with_replacement_combinations("", -1)...] == Any[]

#cool-lex iterator
@test_throws DomainError [CoolLexCombinations(-1, 1)...]
@test_throws DomainError [CoolLexCombinations(5, 0)...]
@test [CoolLexCombinations(4,2)...] == Vector[[1,2], [2,3], [1,3], [2,4], [3,4], [1,4]]
@test isa(iterate(CoolLexCombinations(1000, 20))[2], Combinatorics.CoolLexIterState{BigInt})

# Power set
@test collect(powerset([])) == Any[[]]
@test collect(powerset(['a', 'b', 'c'])) == Any[[],['a'],['b'],['c'],['a','b'],['a','c'],['b','c'],['a','b','c']]
@test collect(powerset(['a', 'b', 'c'], 1)) == Any[['a'],['b'],['c'],['a','b'],['a','c'],['b','c'],['a','b','c']]
@test collect(powerset(['a', 'b', 'c'], 1, 2)) == Any[['a'],['b'],['c'],['a','b'],['a','c'],['b','c']]


text:
mutable struct Limiter
  Δt::Float64
  t₀::Float64
end

function Limiter(Δt::Number)
  return Limiter(Float64(Δt), time())
end

function (lim::Limiter)()
  t = time()
  diff = lim.Δt - (t - lim.t₀)

  if diff > 0.0
    sleep(diff)
  else
    yield()
  end

  lim.t₀ = t
end

max_stars_repo_name: tlqmj/phasespace

max_issues_repo_name: tlqmj/phasespace

max_forks_repo_name: tlqmj/phasespace


text:
using Base: ReinterpretArray

import StaticArrays
using StaticArrays: SVector

const Vector2 = SVector{2, Float64}
const Vector3 = SVector{3, Float64}
const MosiVector = Union{Vector2, Vector3}

const Vector2s = Vector{Vector2}
const Vector3s = Vector{Vector3}
const AbstractVector2s = AbstractVector{Vector2}
const AbstractVector3s = AbstractVector{Vector3}

flatten(x::Vector{T}) where T <: MosiVector = reinterpret(Float64, x)
unflatten(::Type{T}, x::Vector{Float64}) where T <: MosiVector = reinterpret(T, x)
unflatten(x::ReinterpretArray{Float64, 1, T, Vector{T}}) where T <: MosiVector = x.parent

to_matrix(v::Vector2s) = reshape(flatten(v), 2, length(v))
to_matrix(v::Vector3s) = reshape(flatten(v), 3, length(v))

ndims(::Type{Vector2}) = 2
ndims(::Type{Vector3}) = 3
ndims(::T) where T <: MosiVector = ndims(T)


text:
"""
	firstdifference(obj::AbstractDataFrame, PID::Symbol, TID::Symbol;
		gap::Integer = 1)

	This function returns the first-difference transformation of the dataframe.
	It uses a panel ID to perform the transformation panel wise and a temporal
	ID which is used to determine the order. The gap determines if a valid
	observation may skip time periods. The output has missing values when an
	observation is not valid. The first observation for the temporal ID is kept
	as the original value.
"""
function firstdifference(obj::AbstractDataFrame,
						 PID::Symbol,
						 TID::Symbol)
	varlist = names(obj)
	obj = copy(obj)
	allowmissing!(obj)
	sort!(obj, cols = [PID, TID])
	categorical!(obj)
	Step = by(obj, PID) do subdf
        frequency(subdf[TID])
    end
	Step = filter(elem -> elem ≥ zero(typeof(elem)), Step[:x1])
	@assert !isempty(Step) "No positive step."
	Step = minimum(Step)
	Gaps = gaps(obj, Step, PID = PID, TID = TID)
	Gaps = Gaps[:Valid]
	categorical = setdiff(names(obj)[broadcast(<:,
						typeof.(obj.columns),
						AbstractCategoricalVector)],
						[PID, TID])
	todiff = obj[:, union([PID],
					setdiff(names(obj), union([TID], categorical)))]
	todiff = by(todiff, PID) do subdf
		if size(subdf, 1) == 1
			output = subdf[:,2:end]
		else
			output = subdf[:,2:end]
			for (name, col) ∈ eachcol(output)
				output[:,name] = vcat(missing, diff(col))
			end
		end
		return output
	end
	output = hcat(obj[[PID, TID]], todiff, obj[categorical])
	output = output[Gaps,varlist]
	return output
end


text:
a = parse(Int, readline())
b, c = parse.(split(readline())) # parse.(Int, split(readline()))
s = chomp(readline())
println(string(a + b + c) * " " * s)

max_stars_repo_name: minoki/my-atcoder-solutions

max_issues_repo_name: minoki/my-atcoder-solutions

max_forks_repo_name: minoki/my-atcoder-solutions


text:
using Test

@testset "ApplicationBuilder.jl" begin

# TODO: Make the tests work on Windows and Linux!!! :'(
@static if Sys.isapple()

@testset "Test ApplicationBuilder (by compiling examples/*.jl)" begin
    include("ApplicationBuilder.jl")
end
@testset "Command-line interface (compiling examples/*.jl)" begin
    include("build_app-cli.jl")
end

else  # Windows and Linux

@testset "bundle.jl" begin
    include("bundle.jl")
end

end
end


text:
using Test
using LogicCircuits
using ProbabilisticCircuits

@testset "Circuit saver test" begin
    mktempdir() do tmp

        circuit, vtree = load_struct_prob_circuit(
                            zoo_psdd_file("little_4var.psdd"), zoo_vtree_file("little_4var.vtree"))

        # load, save, and load as .psdd
        save_circuit("$tmp/temp.psdd", circuit, vtree)
        save(vtree, "$tmp/temp.vtree");
        load_struct_prob_circuit("$tmp/temp.psdd", "$tmp/temp.vtree")

        # save and load as .sdd
        save_circuit("$tmp/temp.sdd", circuit, vtree)
        save(vtree, "$tmp/temp.vtree")

    end

    #TODO add some actual @test statements
    @test true
end


text:
export AugmentDataLabel

function AugmentDataLabel(data::AbstractArray{T,N}, label, image_weights,flip_dims,permute_dims) where {T,N}

  for i in flip_dims
    if rand() > 0.5
      data          = reverse(data,dims=i)
      if isempty(label) == false
        label         = reverse(label,dims=i)
      end
      if isempty(image_weights) == false
        image_weights = reverse(image_weights,dims=i)
      end
    end
  end

  if isempty(permute_dims)==false
    if rand() > 0.5
      permute_dims_augment       = shuffle(permute_dims)
      permutation                = range(1,stop=N,step=1); permutation = convert(Vector{Any},permutation)
      permutation[permute_dims] .=  permute_dims_augment

      data          = permutedims(data,permutation)
      if isempty(label) == false
        label         = permutedims(label,permutation)
      end
      if isempty(image_weights) == false
        image_weights = permutedims(image_weights,permutation)
      end
    end
  end

  return data, label, image_weights
end


text:
# Example 1 of
# Guglielmi, Protasov,
# "Computing Lyapunov exponents of switching systems."
# The Lyapunov Exponent was unknown
# found by SOS-24
# It is between
# (-0.7163555099121912 given by [(2,1.0786514),(1,2.467595)])
# -0.7157393916973624 given by [(2,1.0576859),(1,2.492785)] found by SOS-24
# -0.7148034620151038 given by [(2,1.02),(1,2.539)] improved from SOS-24 by hand
# and
# -0.7120361328125 given by SOS-24

@testset "Example 1" begin
    s = ContinuousSwitchedSystem([logm([1 1; -1 1]/3),logm([1 1; -1 0]/3)])
    #smp = DiscretePeriodicSwitching(s, [1, 2])
    const expected = [-0.6483688354492188, # 1
                      -0.6867463399790471, # 2
                      -0.703767519177436,  # 3
                      -0.7071168233766547, # 4
                      -0.7073838521852485, # 5
                      -0.7096574117458425, # 6
                      -0.710710268191155,  # 7
                      -0.7107560445583425, # 8
                      -0.7113358785427175, # 9
                      -0.7117890677652648, # 10
                      -0.7118256033710013, # 11
                      -0.7120087088397513] # 12
    for factory in sdp_factories
        for d in 1:12
            tol = (d < 10 ? (ismosek(factory) ? 1e-5 : 1e-4) : (ismosek(factory) ? 1e-4 : 1e-3))
            lb, ub = soslyapb(s, d, factory=factory, tol=tol)
            @test lb == -Inf
            if ismosek(factory)
                @test isapprox(ub, expected[d])
            else
                @test isapprox(ub, expected[d], rtol=2*tol)
            end
#               for k in 1:100
#                   psw = sosbuildsequence(s, d, p_0=:Random, niter=20)
#                   @test psw !== nothing
#                   @show psw.period
#                   @show psw.growthrate
#               end
        end
    end
    #@test getsmp(s) == smp
    #@test isapprox(s.lb, 3.917384715148, rtol=1e-12)
    if isempty(sdp_factories)
        @test isapprox(s.ub, Inf)
    else
        @test isapprox(s.ub, expected[12], rtol=1e-3)
        #@test_approx_eq(s.ub, -0.7107086181640625)
    end
end


text:
module ASCII
using Dates, Mmap, Printf, SeisIO, SeisIO.FastIO, SeisIO.Formats

# imports
include("ASCII/imports.jl")
include("ASCII/GeoCSV.jl")
include("ASCII/SLIST.jl")

# exports
export  formats,
        read_geocsv_file!,
        read_geocsv_slist!,
        read_geocsv_tspair!,
        read_slist!
end


text:
# This is not run by runtests.jl. Only intended to be checked manually.
using Plots
const test_funcs = [
    (x -> 2, 0, 1),
    (x -> x, -1, 1),
    (x -> 5x, -1, 1),
    (x -> 1 / x, 0, 1),
    (x -> 1 / x, -1, 1),
    (x -> tan(x), 0, 4),
    (x -> -1 / x, 0, 1),
    (x -> 1 / abs(x), -1, 1),
    (x -> 1e6x, -1, 1),
    (x -> 1e50x, -1, 1),
    (x -> log(1 + sin(cos(x))), -6, 6),
    (x -> sin(x^3) + cos(x^3), 0, 6.28),
    (x -> sin(x), -5, 200),
    (x -> sin(1 / x), -2, 2),
    (x -> sin(1 / x), 0, 2),
    (x -> sin(x^4), -4, 4),
    (x -> sin(300x), -4, 4),
    (x -> 1 + x^2 + 0.0125log(abs(1 - 3(x - 1))), -2, 2),
    (x -> sin(exp(x)), -6, 6),
    (x -> 1 / sin(x), -10, 10),
    (x -> sin(x) / x, -6, 6),
    (x -> tan(x^3 - x + 1) + 1 / (x + 3exp(x)), -2, 2),
]
##
plots = [plot(func...) for func in test_funcs]
pitchfork = plot(x -> sqrt(x), 0, 10)
plot!(pitchfork, x -> -sqrt(x), 0, 10)
plot!(pitchfork, x -> 0, -10, 0)
plot!(pitchfork, x -> 0, 0, 10, linestyle = :dash)
push!(plots, pitchfork)
display.(plots)