text:
"""
This file contains the code for the work by

Y. Chandak, S. Niekum, B. Castro da Silva, E. Learned-Miller, E. Brunskill, and P. S. Thomas.
Universal Off-Policy Evaluation.
In Advances in Neural Information Processing Systems, 2021

Code repo: https://github.com/yashchandak/UnO
"""

using Statistics
using Random
using Bootstrap
using BlackBoxOptim

include("./UnO_utils.jl")

#######################################################################
# Estimates for different parameters
#######################################################################??

function CDF_mean(xs, ys)
    """
    Plug-in mean estimate using a given CDF

    Args
    ----
    xs: x values
    ys: CDF estimates at those xs
    """

    ys_pdf = decumsum(ys)
    return sum(ys_pdf .* xs)
end

function CDF_variance(xs, ys)
    """
    Plug-in variance estimate using a given CDF

    Args
    ----
    xs: x values
    ys: CDF estimates at those xs
    """

    ys_pdf = decumsum(ys)
    temp_mean = sum(ys_pdf .* xs)
    return sum(ys_pdf .* ((xs .- temp_mean) .^2))
end

function CDF_entropy(xs, ys)
    """
    Plug-in mean entropy using a given CDF

    Args
    ----
    xs: x values
    ys: CDF estimates at those xs
    """
    ys_pdf = decumsum(ys)
    total = 0
    for item in ys_pdf
        if item > 0.00001
            total += item * log(item)
        end
    end
    return - total
    # return - sum(ys_pdf .* log.(ys_pdf))
end

function CDF_quantile(xs, ys, α)
    """
    Plug-in quantile estimate using a given CDF

    Args
    ----
    xs: x values
    ys: CDF estimates at those xs
    """
    id = searchsortedfirst(ys, α)
    id = min(id, length(ys))        # Overflow check because there might not be any ys > α
    return xs[id]
end

function CDF_IQR(xs, ys, α1, α2)
    """
    Plug-in inter-quantile range estimate using a given CDF

    Args
    ----
    xs: x values
    ys: CDF estimates at those xs
    """
    return CDF_quantile(xs, ys, α2) - CDF_quantile(xs, ys, α1)
end

function CDF_CVaR(xs, ys, α)
    """
    Plug-in conditional value-at-risk estimate using a given CDF

    Args
    ----
    xs: x values
    ys: CDF estimates at those xs
    """

    id = searchsortedfirst(ys, α)
    id = min(id, length(ys))     # Overflow check because there might not be any ys > α
    return (1/α) * CDF_mean(xs[1:id], clamp.(ys[1:id], 0, α))
end

#######################################################################
# CDF based bounds for different parameters
#######################################################################

function CDF_mean_bound(xs, lb_ys, ub_ys)
    """
    Bound for the mean

    Anderson (1967)
    https://statistics.stanford.edu/sites/g/files/sbiybj6031/f/AND%20ONR%2001.pdf

    Args
    ----
    xs: Sorted sample values
    lb_ys: Lower bounds CDF at those xs
    ub_ys:  Upper bound CDF at those xs
    """

    lb_mean = CDF_mean(xs, ub_ys)  # Lower bound on mean is computed using upper bound of CDF
    ub_mean = CDF_mean(xs, lb_ys)  # Upper bound on mean is computed using lower bound of CDF
    return lb_mean, ub_mean
end

function CDF_variance_bound(xs, lb_ys, ub_ys)
    """
    Bound for the variance

    Romano et. al (2002)
    https://www.tandfonline.com/doi/abs/10.1081/STA-120006065?journalCode=lsta20

    Args
    ----
    xs: Sorted sample values
    lb_ys: Lower bounds CDF at those xs
    ub_ys:  Upper bound CDF at those xs
    """

    xs_n = length(xs)
    temp_xs = fill(0.0, xs_n + 1)
    temp_ys = fill(0.0, xs_n + 1)

    ###########################
    # Lower bound for variance
    ##########################
    function min_var!(xs, lb_ys, ub_ys, temp_xs, temp_ys, x_cutpoint)
        """
        This function creates a CDF that follows the lower bound CDF
        till the cutpoint, then does a vertical jump,
        and then follows the upper bound CDF from there onwards
        """

        x_cutpoint = x_cutpoint[1]
        lower_id = searchsortedlast(xs, x_cutpoint)
        temp_xs[1:lower_id] = xs[1:lower_id]
        temp_xs[lower_id + 1] = x_cutpoint
        temp_xs[(lower_id + 2):(xs_n + 1)] = xs[(lower_id + 1):xs_n]

        temp_ys[1:lower_id] = lb_ys[1:lower_id]
        temp_ys[(lower_id + 1):(xs_n + 1)] = ub_ys[lower_id:xs_n]

        return CDF_variance(temp_xs, temp_ys)
    end

    res = bboptimize(candidate -> min_var!(xs, lb_ys, ub_ys, temp_xs, temp_ys, candidate),
                     SearchRange = (xs[1], last(xs)), NumDimensions = 1, TraceMode=:silent)
    var_lb = best_fitness(res)

    ###########################
    # Upper bound for variance
    ##########################
    function max_var!(xs, lb_ys, ub_ys, temp_xs, temp_ys, y_cutpoint)
        """
        This function creates a CDF that follows the upper bound CDF
        till the cutpoint, then does a horizontal jump,
        and then follows the lower bound CDF from there onwards
        """

        y_cutpoint = y_cutpoint[1]
        upper_id = searchsortedfirst(ub_ys, y_cutpoint)
        upper_id2 = searchsortedfirst(lb_ys, y_cutpoint)
        temp_ys[1:(upper_id - 1)] = ub_ys[1:(upper_id - 1)]
        temp_ys[upper_id:(upper_id2 -1)] .= y_cutpoint
        temp_ys[upper_id2:(xs_n)] = lb_ys[upper_id2:xs_n]

        return - CDF_variance(xs[1:xs_n], temp_ys[1:xs_n]) # Negative because using max-opt but we need min
    end

    res = bboptimize(candidate -> max_var!(xs, lb_ys, ub_ys, temp_xs, temp_ys, candidate),
                     SearchRange = (0.0, 1.0), NumDimensions = 1, TraceMode=:silent)
    var_ub = - best_fitness(res)

    return var_lb, var_ub
end

function CDF_quantile_bound(xs, lb_ys, ub_ys, α)
    """
    Bound for the quantile

    Args
    ----
    xs: Sorted sample values
    lb_ys: Lower bounds CDF at those xs
    ub_ys:  Upper bound CDF at those xs
    α: Quantile level
    """
    lb_quantile = CDF_quantile(xs, ub_ys, α)   # Lower bound on quantile is computed using upper bound of CDF
    ub_quantile = CDF_quantile(xs, lb_ys, α)   # Upper bound on quantile is computed using lower bound of CDF
    return lb_quantile, ub_quantile
end

function CDF_IQR_bound(xs, lb_ys, ub_ys, α1, α2)
    """
    Bound for the inter-quantile range

    Args
    ----
    xs: Sorted sample values
    lb_ys: Lower bounds CDF at those xs
    ub_ys:  Upper bound CDF at those xs
    α1, α2: Quantile levels
    """

    # Lower bound on IQR: lowest for α2 - maximum for α1
    lb_IQR = CDF_quantile(xs, ub_ys, α2) - CDF_quantile(xs, lb_ys, α1)
    lb_IQR = max(0, lb_IQR)

    # Upper bound on IQR: maximum for α2 - minimum for α1
    ub_IQR = CDF_quantile(xs, lb_ys, α2) - CDF_quantile(xs, ub_ys, α1)

    return lb_VaR, ub_VaR
end

function CDF_CVaR_bound(xs, lb_ys, ub_ys, α)
    """
    Bound for the mean

    # Thomas et. al (2019)
    # https://people.cs.umass.edu/~pthomas/papers/Thomas2019.pdf

    Args
    ----
    xs: Sorted sample values
    lb_ys: Lower bounds CDF at those xs
    ub_ys:  Upper bound CDF at those xs
    """

    lb_CVaR = CDF_CVaR(xs, ub_ys, α)   # Lower bound on CVaR is computed using upper bound of CDF
    ub_CVaR = CDF_CVaR(xs, lb_ys, α)   # Upper bound on CVaR is computed using lower bound of CDF
    return lb_CVaR, ub_CVaR
end

function CDF_entropy_bound(xs, lb_ys, ub_ys)
    """
    Bound for the entropy

    DeStefano et. al (2005)
    https://arxiv.org/ftp/cs/papers/0504/0504091.pdf

    The following black-box search method is naive and not ideal.
    Better way would be to actually code up the 'string-tightening'
    algorithm in the work by DeStefano et al.

    Args
    ----
    xs: Sorted sample values
    lb_ys: Lower bounds CDF at those xs
    ub_ys:  Upper bound CDF at those xs
    """

    return CDF_bbox_bound(xs, lb_ys, ub_ys, CDF_entropy)
end

function CDF_bbox_bound(xs, lb_ys, ub_ys, fn)
    """
    A generic blackbox optimization procedure to search for CDF within a band
    that optimizes for the desired parameter.

    The space of all CDFs is parameterized using "n_points".
    Increasing n_points can lower the discretization error
    but can make the optimization problem harder.

    IMP: Typically, this should be the last resort. First shot should always be to leverage
    geometric insights as discussed in Figure 3 and 6 to obtain closed for solutions.
    """

    n_points = 50
    Gmin, Gmax = xs[1], last(xs)

    # Create constraints for the BBO search from the CDF band
    lb_constraints = get_equispaced_CDF2(xs, lb_ys, Gmin, Gmax, n_points)[2]
    ub_constraints = get_equispaced_CDF2(xs, ub_ys, Gmin, Gmax, n_points)[2]

    function wrapper(candidate, presign, fn)
        """
        Projects the candidate solution in the desired space (valid CDF + bounded)
        then computes the parameter on the projected CDF

        Args
        ----
        candidate: Candidate xs and ys
        presign: Plus or negative to determine maximization or minimization
        fn: Plug-in estimator for the desired parameter
        """

        can_xs = candidate[1:n_points] .* abs(Gmax - Gmin)
        can_ys = candidate[n_points+1:2*n_points]

        candidate_xs = cumsum(can_xs)
        if sum(can_xs) > (Gmax - Gmin)
            candidate_xs = Gmin .+ (Gmax - Gmin) .* candidate_xs / sum(can_xs)
        else
            candidate_xs = Gmin .+ candidate_xs
        end

        lb_constraints = get_CDF_at_locs(xs, lb_ys, candidate_xs)
        ub_constraints = get_CDF_at_locs(xs, ub_ys, candidate_xs)
        candidate_ys = clamp.(cumsum(can_ys), lb_constraints, ub_constraints)

        return presign * fn(candidate_xs, candidate_ys)
    end

    # Ranges are set to (0, 0.1) as the values get rescaled in the wrapper

    # Lower bound for the parameter
    res = bboptimize(candidate -> wrapper(candidate, 1.0, fn),
                     SearchRange = (0.0, 0.1), NumDimensions = 2*n_points, TraceMode=:silent, method=:random_search)
    bbox_lb = best_fitness(res)

    # Upper bound for the parameter
    res = bboptimize(candidate -> wrapper(candidate, -1.0, fn), #negative fn
                     SearchRange = (0.0, 0.1), NumDimensions = 2*n_points, TraceMode=:silent, method=:random_search)
    bbox_ub = - best_fitness(res)

    return bbox_lb, bbox_ub
end

#######################################################################
# Bootstrap based bounds for different parameters
#######################################################################??

function CDF_mean_boot_bound(xs, ys, delta)
    return CDF_param_boot_bounds(xs, ys, delta, (x,y) -> CDF_mean(x,y))
end

function CDF_entropy_boot_bound(xs, ys, delta)
    return CDF_param_boot_bounds(xs, ys, delta, (x,y) -> CDF_entropy(x,y))
end

function CDF_variance_boot_bound(xs, ys, delta)
    return CDF_param_boot_bounds(xs, ys, delta, (x,y) -> CDF_variance(x,y))
end

function CDF_quantile_boot_bound(xs, ys, α, delta)
    return CDF_param_boot_bounds(xs, ys, delta, (x,y) -> CDF_quantile(x,y,α))
end

function CDF_IQR_boot_bound(xs, ys, α1, α2, delta)
    return CDF_param_boot_bounds(xs, ys, delta, (x,y) -> CDF_IQR(x,y,α1,α2))
end

function CDF_CVaR_boot_bound(xs, ys, α, delta)
    return CDF_param_boot_bounds(xs, ys, delta, (x,y) -> CDF_CVaR(x,y,α))
end

function CDF_param_boot_bounds(xs, ys, delta, param_fn)
    """
    Bootstrap based bounds for CDF's parameters

    Args
    ----
    xs: sorted sample values
    ys: CDF values at those xs
    delta:  Acceptable failure rate
    param_fn: Plug-in function for the desired parameter
    """

    function get_param(x, y_pdf)
        """
        Helper function to re-sort data and compute the param estimate
        Resorting is required as re-sampling by bootstrap shuffles data.
        """
        order = sortperm(x)                 # Create order statistic map
        xs_boot = x[order]                  # Location of steps
        ys_boot = cumsum(y_pdf[order])      # Step heights
        return param_fn(xs_boot, ys_boot)
    end

    n = length(xs)
    idxs = collect(1:n)
    ys_pdf = decumsum(ys)

    bs = bootstrap(boot_id -> get_param(xs[boot_id], ys_pdf[boot_id]), idxs, BasicSampling(100))

    try
        # Important:
        # In the bootstrap library, due to some unknown bug a variable zalpha becomes NaN
        # Another way to avoid error is by using the following change to the BCa code IN THE BOOTSTRAP LIBRARY:
        # zalpha == zalpha || return (t0, t0, t0)  # Yash: avoid the edge case of zalpha being NaN

        bci = confint(bs, BCaConfInt(1 - delta))                # BCa Bootstrap
        _, lb_param, ub_param = bci[1]                          # return format: estimate, lower, upper
        return lb_param, ub_param

    catch e
        println("BCa bootstrap error occured. Trying percentile bootstrap...")
        bci = confint(bs, PercentileConfInt(1 - delta))         # Percentile Bootstrap
        _, lb_param, ub_param = bci[1]                          # return format: estimate, lower, upper
        return lb_param, ub_param
    end

end

max_stars_repo_name: yashchandak/UnO

max_issues_repo_name: yashchandak/UnO

max_forks_repo_name: yashchandak/UnO


text:
### A Pluto.jl notebook ###
# v0.18.4

using Markdown
using InteractiveUtils

# ╔═╡ 3e14817e-85f8-4d5b-8a8f-075d1799ed27
using AbstractTrees, PlutoUI, Random, Wordlegames

# ╔═╡ ed649770-a53a-11ec-0434-1bbbc561f8cb
# hideall
title = "Wordle games as a tree";

# ╔═╡ 021364a8-f40b-49fe-9a02-406c0044c9ed
"""
+++
title = "$title"
+++
""" |> Base.Text

# ╔═╡ c49e1040-67c8-4386-a6ce-bbbeace8bc90
md"""
# $title

As described in the previous tutorial, strategies such as maximizing the entropy or minimizing the expected pool size for the next stage can be used to select guesses automatically in Wordle or related games.

When doing so the possible games can be represented in a data structure called a [tree](https://en.wikipedia.org/wiki/Tree_(data_structure)).

Some of the terminology used with these structures is based on concepts of a family tree.

First attach some packages that will be used
"""

# ╔═╡ cec480d3-7348-4d14-9628-007893713a07
md"and create an instance of `wordle` where the guesses are chosen to maximize the entropy, which is the default criterion."

# ╔═╡ 9586bcbd-16aa-4598-8faa-8220755566e6
begin
	datadir = joinpath(pkgdir(Wordlegames), "data")
	wordle = GamePool(collect(readlines(joinpath(datadir, "Wordletargets.txt"))))
end;

# ╔═╡ ddc457c3-f207-4ed1-8a8a-7f1c4cfdaed7
md"Finally, we create a tree from the games for a random selection of 25 targets."

# ╔═╡ 52ed4682-a0d7-44cf-b878-27f1cf9cc91d
gametree25 = tree(wordle, Random.seed!(1234321), 25);

# ╔═╡ 86d0bb82-ed8f-4b03-b90c-31c55e556bae
md"""
## The AbstractTrees package

The [AbstractTrees](https://github.com/JuliaCollections/AbstractTrees.jl) package provides many methods for working with tree data structures.
One of the most useful is `print_tree` which, as the name suggests, prints the tree in a special format.
(Because the content for these tutorials is generated as [Pluto](https://github.com/fonsp/Pluto.jl.git) notebooks, we need to wrap the call to `print_tree` in `with_terminal() do ... end` to have the output displayed.
Outside of Pluto this is not necessary.
"""

# ╔═╡ 30ba17b1-9c0c-4fd6-9f40-3c3c1c986e93
with_terminal() do
	print_tree(gametree25; maxdepth=8)
end

# ╔═╡ 384abc84-3fda-4ba3-b52e-2d46052b8ce6
md"""
Each guess in a game constitutes a `"node"` in the tree.
The initial guess in any of the games is `"raise"`, which is the `"root"` node for the tree.
A node can have zero or more `"children"` which are its immediate descendents.

The nodes in this tree are each a `GameNode` struct with a `"children"` field.
"""

# ╔═╡ 0c24b1ca-b36a-46f6-9da7-e5c3e3670d29
typeof(gametree25)

# ╔═╡ 68f9375f-31af-4916-933f-fb41905819f8
fieldnames(GameNode)

# ╔═╡ e1735528-469e-4a17-8b8a-d95dca03f46c
length(gametree25.children)  # number of children of the root node

# ╔═╡ 112b0452-7ae4-49b4-bbdd-54fb00ff8bdb
md"""
The `score` field of a `GameNode` is similar to the elements of the `guesses` field of a `GamePool` object but with one important difference.
In a `GameNode` the `score` and `sc` fields are the score that will produce the guess, as opposed to the score for the guess,

That is, the first child of `"raise"` is `"pilot"` which is the next guess in a game in which `"raise"` returns a score of `"🟫🟫🟨🟫🟫"` as tiles or 9 as a decimal number.
"""

# ╔═╡ f8f1521c-fbbb-4648-af9c-0a15c7e787ff
first(gametree25.children).score

# ╔═╡ 1dd5c9fb-820f-41c6-b796-494d68479f1c
md"""
The 16 children of the root node from these 25 games are
"""

# ╔═╡ 8302e26a-53ae-4f9f-9579-26a1c839c865
[child.score.guess for child in gametree25.children]

# ╔═╡ 7f4a7f14-fe77-43fa-b2dd-37ba70cf9da3
md"""
Some of these children have many descendents.
When generating the tree the children of a node are ordered according to the size of the tree rooted at that node.

Because a subtree is exactly the same type of structure as a tree, we can print a subtree with `print_tree`.
"""

# ╔═╡ 088ca1f6-44c3-4d0d-a68b-c71881813410
with_terminal() do
	print_tree(first(gametree25.children))
end

# ╔═╡ 785fe89b-4008-437a-ac73-2bdaa26e7cf6
md"""
We see that the size of the tree rooted at `"pilot"` is 10.
"""

# ╔═╡ 8d83a193-aed2-4c1d-8bc5-614a81d16548
md"""
The "leaves" of a tree are the terminal nodes, i.e. the nodes that do not have children.
"""

# ╔═╡ 576d0904-3f3e-49e8-9ce7-ccdd8767bf1b
[leaf.score.guess for leaf in Leaves(gametree25)]

# ╔═╡ 144a9d04-f971-48d9-ace4-dd0156f591c2
md"""
It happens in this case that all of the targets that generated the tree are leaves in this tree, but that is not necessarily the case.
"""

# ╔═╡ 105a0522-729a-49ea-b0a5-62b868ba5e45
length(collect(Leaves(gametree25)))

# ╔═╡ efcb3473-6ea1-459c-9995-32b0a407a2f4
md"""
## Creating the tree structure

As a language Julia gets high marks for "composability" - the ability to adapt one package to use concepts from another package.
The use of generic functions and multiple dispatch is central to this enhanced compositibility.

All that is necessary to use many of the functions in `AbstractTrees.jl` on the trees created from a collection of games using a particular `GamePool` is to define the `GameNode` struct, the method of generating the tree, and methods for `AbstractTrees.children`, `AbstractTrees.nodetype` and `AbstractTrees.printnode`
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
Wordlegames = "1cb69566-e1cf-455f-a587-fd79a2e00f5a"

[compat]
AbstractTrees = "~0.3.4"
PlutoUI = "~0.7.38"
Wordlegames = "~0.3.0"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.8.0-beta3"
manifest_format = "2.0"
project_hash = "2c8800ff06b3096b67485eb059a224666198ba30"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "8eaf9f1b4921132a4cff3f36a1d9ba923b14a481"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.4"

[[deps.AbstractTrees]]
git-tree-sha1 = "03e0550477d86222521d254b741d470ba17ea0b5"
uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
version = "0.3.4"

[[deps.ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
version = "1.1.1"

[[deps.Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[deps.Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[deps.ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[deps.Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "96b0bc6c52df76506efc8a441c6cf1adcb1babc4"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.42.0"

[[deps.CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
version = "0.5.2+0"

[[deps.Crayons]]
git-tree-sha1 = "249fe38abf76d48563e2f4556bebd215aa317e15"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.1.1"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "ae02104e835f219b8930c7664b8012c93475c340"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.3.2"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "3daef5523dd2e769dad2365274f760ff5f282c7d"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.11"

[[deps.DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[deps.Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[deps.Downloads]]
deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
version = "1.6.0"

[[deps.FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[deps.FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[deps.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[deps.Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[deps.InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.0"

[[deps.IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "3c837543ddb02250ef42f4738347454f95079d4e"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.3"

[[deps.LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"
version = "0.6.3"

[[deps.LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
version = "7.81.0+0"

[[deps.LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[deps.LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"
version = "1.10.2+0"

[[deps.Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[deps.LinearAlgebra]]
deps = ["Libdl", "libblastrampoline_jll"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[deps.Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[deps.Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[deps.MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
version = "2.28.0+0"

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[deps.MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
version = "2022.2.1"

[[deps.NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
version = "1.2.0"

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
version = "0.3.20+0"

[[deps.OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "621f4f3b4977325b9128d5fae7a8b4829a0c2222"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.4"

[[deps.Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
version = "1.8.0"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "ColorTypes", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "670e559e5c8e191ded66fa9ea89c97f10376bb4c"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.38"

[[deps.PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "28ef6c7ce353f0b35d0df0d5930e0d072c1f5b9b"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.4.1"

[[deps.PrettyTables]]
deps = ["Crayons", "Formatting", "Markdown", "Reexport", "Tables"]
git-tree-sha1 = "dfb54c4e414caa595a1f2ed759b160f5a3ddcba5"
uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
version = "1.3.1"

[[deps.Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[deps.REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[deps.Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[deps.SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
version = "0.7.0"

[[deps.Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[deps.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[deps.Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[deps.SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[deps.SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[deps.Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[deps.TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
version = "1.0.0"

[[deps.TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[deps.Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "OrderedCollections", "TableTraits", "Test"]
git-tree-sha1 = "5ce79ce186cc678bbb5c5681ca3379d1ddae11a1"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.7.0"

[[deps.Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
version = "1.10.0"

[[deps.Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[deps.UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[deps.Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[deps.Wordlegames]]
deps = ["AbstractTrees", "DataFrames", "Random", "Tables"]
git-tree-sha1 = "4c463de78d2f3f9447b695e241eba43cc945f866"
uuid = "1cb69566-e1cf-455f-a587-fd79a2e00f5a"
version = "0.3.0"

[[deps.Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"
version = "1.2.12+1"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"
version = "5.1.0+0"

[[deps.nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"
version = "1.41.0+1"

[[deps.p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
version = "16.2.1+1"
"""

# ╔═╡ Cell order:
# ╟─ed649770-a53a-11ec-0434-1bbbc561f8cb
# ╟─021364a8-f40b-49fe-9a02-406c0044c9ed
# ╟─c49e1040-67c8-4386-a6ce-bbbeace8bc90
# ╠═3e14817e-85f8-4d5b-8a8f-075d1799ed27
# ╟─cec480d3-7348-4d14-9628-007893713a07
# ╠═9586bcbd-16aa-4598-8faa-8220755566e6
# ╟─ddc457c3-f207-4ed1-8a8a-7f1c4cfdaed7
# ╠═52ed4682-a0d7-44cf-b878-27f1cf9cc91d
# ╟─86d0bb82-ed8f-4b03-b90c-31c55e556bae
# ╠═30ba17b1-9c0c-4fd6-9f40-3c3c1c986e93
# ╟─384abc84-3fda-4ba3-b52e-2d46052b8ce6
# ╠═0c24b1ca-b36a-46f6-9da7-e5c3e3670d29
# ╠═68f9375f-31af-4916-933f-fb41905819f8
# ╠═e1735528-469e-4a17-8b8a-d95dca03f46c
# ╟─112b0452-7ae4-49b4-bbdd-54fb00ff8bdb
# ╠═f8f1521c-fbbb-4648-af9c-0a15c7e787ff
# ╟─1dd5c9fb-820f-41c6-b796-494d68479f1c
# ╠═8302e26a-53ae-4f9f-9579-26a1c839c865
# ╟─7f4a7f14-fe77-43fa-b2dd-37ba70cf9da3
# ╠═088ca1f6-44c3-4d0d-a68b-c71881813410
# ╟─785fe89b-4008-437a-ac73-2bdaa26e7cf6
# ╟─8d83a193-aed2-4c1d-8bc5-614a81d16548
# ╠═576d0904-3f3e-49e8-9ce7-ccdd8767bf1b
# ╟─144a9d04-f971-48d9-ace4-dd0156f591c2
# ╠═105a0522-729a-49ea-b0a5-62b868ba5e45
# ╟─efcb3473-6ea1-459c-9995-32b0a407a2f4
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

ext: jl

max_stars_repo_name: dmbates/WordleTutorials

max_issues_repo_name: dmbates/WordleTutorials

max_forks_repo_name: dmbates/WordleTutorials


text:
module DataUtils

using Base: StridedFastContiguousSubArray, AbstractFloat
using Flux
using Flux.Data: DataLoader
using Random, DelimitedFiles, Distributed
using ProgressMeter: Progress, next!
using StatsBase, Statistics

export get_train_validation, get_test_set

function scale_tensor(X::T) where {T<:Array{F} where F<:AbstractFloat}
    l, m, _, n = size(X)
    StatsBase.fit(ZScoreTransform, reshape(X, l, m*n))
end

function transform_tensor!(Transform::TF, X::T) where {TF<:AbstractDataTransform, T<:Array{F} where F<:AbstractFloat}
    l, m, _, n = size(X)
    StatsBase.transform!(Transform, reshape(X, l, m*n))
end

function create_tensor(X::T) where {T<:Array{A} where A<:Array{F} where F<:AbstractFloat}
    l = length(X)
    m, n = size(X[1])
    out = Array{eltype(X[1])}(undef, l, n, 1, m)
    for i = 1:l
        for j = 1:m
            out[i, :, :, j] = X[i][j,:]
        end
    end
    return out
end

function data_prep(data_dir::String)
    files = readdir(data_dir * "/Inertial_Signals")
    @info "Processing data..."
    out = Array{Array{Float32}}(undef, length(files))
    progress = Progress(length(files))
    for i = 1:length(files)
        out[i] = readdlm(data_dir * "/Inertial_Signals/" * files[i], Float32)
        next!(progress)
    end
    return create_tensor(out)
end

function get_train_validation(X::T1, Y::T2, batch_size::Int, train_prop::F, loc, scale::Bool=true, shuffle::Bool=true) where {T1<:Array{F1} where F1<:AbstractFloat, T2<:Array{Int}, F<:AbstractFloat}
    train_prop < 1 ? nothing : error("Validation set must have entries. Please use a train_prop value less than 1.")
    train_prop > 0 ? nothing : error("Training set must have entries. Please use a train_prop value greater than 0.")
    l = length(Y)

    if shuffle
        ind = randperm(l)
    else
        ind = 1:l
    end

    T = scale_tensor(X)
    if scale
        transform_tensor!(T, X)
    end

    classes = sort!(unique(Y))

    idx = 1:Int(floor(train_prop*l))
    X_train = X[:,:,:,ind][:,:,:,idx] |> loc
    Y_train = Flux.onehotbatch(Y[ind][idx], classes) |> loc
    X_val = X[:,:,:,ind][:,:,:,last(idx)+1:l] |> loc
    Y_val = Flux.onehotbatch(Y[ind][last(idx)+1:l], classes) |> loc

    train_set = DataLoader((X_train, Y_train), batchsize=batch_size, shuffle=shuffle)

    val_set = DataLoader((X_val, Y_val), batchsize=size(Y_val,2), shuffle=shuffle)

    return train_set, val_set, T
end

function get_test_set(X, Y, T, loc, scale, shuffle)
    classes = sort!(unique(Y))
    if scale
        transform_tensor!(T, X)
    end
    X_test = X |> loc
    Y_test = Flux.onehotbatch(reshape(Y, size(Y, 1)), classes) |> loc

    test_set = DataLoader((X_test, Y_test), batchsize=size(Y, 2), shuffle=shuffle)

    return test_set
end

end


text:
#TODO: Make a constructor
#  itensor(::Vector{Pair{Block, Array}}, ::IndexSet)
#for example:
#  itensor([Block(1,1) => randn(2,2),
#           Block(2,2) => randn(3,3)], i', dag(i)))
#to set nonzero blocks.

"""
    itensor(::Array, ::QNIndexSet)

Create a block sparse ITensor from the input Array, where zeros are
dropped and nonzero blocks are determined about the zero values of
the array.
"""
function ITensors.itensor(A::Array{ElT},
                          inds::ITensors.QNIndexSet) where {ElT <: Number}
  length(A) ≠ dim(inds) && throw(DimensionMismatch("In ITensor(::Array, ::IndexSet), length of Array ($(length(A))) must match total dimension of IndexSet ($(dim(inds)))"))
  T = emptyITensor(ElT, inds)
  A = reshape(A, dims(inds))
  for vs in eachindex(T)
    Avs = A[vs]
    if !iszero(Avs)
      T[vs] = A[vs]
    end
  end
  return T
end

ITensors.itensor(A::Array{<:Number},
                 inds::ITensors.QNIndex...) =
  itensor(A, IndexSet(inds...))


text:
abstract type SignedDistanceField <: Function end

struct SDFKDTree <: SignedDistanceField
    tree::KDTree
    extrusion::Float64
end

function (sdf::SignedDistanceField)(point::AbstractPoint)
    @show "No sdf distance implemented for this type:", typeof(sdf), dump(sdf)
end

@inline function (sdf::SDFKDTree)(point::AbstractPoint)
    distance = knn(sdf.tree, point, 1)[2][1] - sdf.extrusion
end

function curve_to_kdtree_sdf(γ::FunctionCurve, extrusion::Float64, arc_length_sampling_freq::Float64)
    t_to_s, s_to_t, γs, sspan = arc_length_parameterization(γ.γ, γ.tspan)
    mins, maxs = sspan[1:2]

    γsamples = (mins:arc_length_sampling_freq:maxs)
    γsampled_points = map(γs, γsamples)

    γkdtree = KDTree(γsampled_points)
    SDFKDTree(γkdtree, extrusion)
end

struct UnionSDF{SDF1 <: SignedDistanceField, SDF2 <: SignedDistanceField} <: SignedDistanceField
    sdf1::SDF1
    sdf2::SDF2
end

@inline function (sdf::UnionSDF)(point::AbstractPoint)
    distance = min(sdf.sdf1(point), sdf.sdf2(point))
end

struct IntersectSDF{SDF1 <: SignedDistanceField, SDF2 <: SignedDistanceField} <: SignedDistanceField
    sdf1::SDF1
    sdf2::SDF2
end

@inline function (sdf::IntersectSDF)(point::AbstractPoint)
    distance = max(sdf.sdf1(point), sdf.sdf2(point))
end

struct FunctionSDF <: SignedDistanceField
    sdf::Function
end

@inline function (sdf::FunctionSDF)(point::AbstractPoint)
    distance = sdf.sdf(point)
end

function curves_to_union_sdf(γs::Vector{Function}, tspans::Vector{Tuple{Float64, Float64}},
        extrusions::Vector{Float64}, arc_length_parameterization::Float64; boundary_sdf=nothing)
    γcurves = [FunctionCurve(γ, tspan) for (γ, tspan) in zip(γs, tspans)]
    γsdfs = [curve_to_kdtree_sdf(γ, extrusion, arc_length_parameterization) for (γ, extrusion) in zip(γcurves, extrusions)]
    unionsdf = reduce((sdf1, sdf2) -> UnionSDF(sdf1, sdf2), γsdfs)
    if !(boundary_sdf isa Nothing) && boundary_sdf isa SignedDistanceField
        sdf = IntersectSDF(unionsdf, boundary_sdf)
    else
        sdf = unionsdf
    end
end

max_stars_repo_name: zoemcc/geometry

max_issues_repo_name: zoemcc/geometry

max_forks_repo_name: zoemcc/geometry


text:
# Variables utilities

"Add variable to model."
function KN_add_var(m::Model)
    ptr_int = Cint[0]
    ret = @kn_ccall(add_var, Cint, (Ptr{Cvoid}, Ptr{Cint}), m.env, ptr_int)
    _checkraise(ret)
    return ptr_int[1]
end

function KN_add_vars(m::Model, nvars::Int)
    ptr_int = zeros(Cint, nvars)
    ret = @kn_ccall(add_vars, Cint, (Ptr{Cvoid}, Cint, Ptr{Cint}), m.env, nvars, ptr_int)
    _checkraise(ret)
    return ptr_int
end

##################################################
## Upper and lower bounds utilies
##################################################

# lower bounds
function KN_set_var_lobnds(m::Model, nindex::Integer, val::Cdouble)
    ret = @kn_ccall(set_var_lobnd, Cint, (Ptr{Cvoid}, Cint, Cdouble), m.env, nindex, val)
    _checkraise(ret)
end

function KN_set_var_lobnds(m::Model, valindex::Vector{Cint}, lobnds::Vector{Cdouble})
    nvar = length(valindex)
    ret = @kn_ccall(set_var_lobnds, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}),
                    m.env, nvar, valindex, lobnds)
    _checkraise(ret)
end

function KN_set_var_lobnds(m::Model, lobnds::Vector{Cdouble})
    ret = @kn_ccall(set_var_lobnds_all, Cint, (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, lobnds)
    _checkraise(ret)
end

# upper bounds
function KN_set_var_upbnds(m::Model, nindex::Integer, val::Cdouble)
    ret = @kn_ccall(set_var_upbnd, Cint, (Ptr{Cvoid}, Cint, Cdouble), m.env, nindex, val)
    _checkraise(ret)
end

function KN_set_var_upbnds(m::Model, valindex::Vector{Cint}, upbnds::Vector{Cdouble})
    nvar = length(valindex)
    ret = @kn_ccall(set_var_upbnds, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}),
                    m.env, nvar, valindex, upbnds)
    _checkraise(ret)
end

function KN_set_var_upbnds(m::Model, upbnds::Vector{Cdouble})
    ret = @kn_ccall(set_var_upbnds_all, Cint, (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, upbnds)
    _checkraise(ret)
end

##################################################
## Fix bounds
##################################################
function KN_set_var_fxbnds(m::Model, nindex::Integer, xFxBnd::Cdouble)
    ret = @kn_ccall(set_var_fxbnd, Cint, (Ptr{Cvoid}, Cint, Cdouble),
                    m.env, nindex, xFxBnd)
    _checkraise(ret)
end

function KN_set_var_fxbnds(m::Model, xIndex::Vector{Cint}, xFxBnds::Vector{Cdouble})
    nvar = length(xIndex)
    @assert length(xFxBnds) == nvar
    ret = @kn_ccall(set_var_fxbnds, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}),
                    m.env, nvar, xIndex, xFxBnds)
    _checkraise(ret)
end

function KN_set_var_fxbnds(m::Model, xFxBnds::Vector{Cdouble})
    ret = @kn_ccall(set_var_fxbnds_all, Cint, (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, xFxBnds)
    _checkraise(ret)
end

##################################################
## Variables types
##################################################

"""
Set variable types (e.g. KN_VARTYPE_CONTINUOUS, KN_VARTYPE_BINARY,
KN_VARTYPE_INTEGER). If not set, variables are assumed to be continuous.
"""
function KN_set_var_types(m::Model, valindex::Vector{Cint}, xTypes::Vector{Cdouble})
    nvar = length(valindex)
    @assert nvar == length(xTypes)
    ret = @kn_ccall(set_var_types, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}),
                    m.env, nvar, valindex, xTypes)
    _checkraise(ret)
end

function KN_set_var_types(m::Model, xTypes::Vector{Cint})
    ret = @kn_ccall(set_var_types_all, Cint, (Ptr{Cvoid}, Ptr{Cint}),
                    m.env, xTypes)
    _checkraise(ret)
end

function KN_set_var_type(m::Model, indexVar::Integer, xType::Integer)
    ret = @kn_ccall(set_var_type, Cint, (Ptr{Cvoid}, Cint, Cint),
                    m.env, indexVar, xType)
    _checkraise(ret)
end

##################################################
## Variables properties
##################################################

"""
Specify some properties of the variables.  Currently
this API routine is only used to mark variables as linear,
but other variable properties will be added in the future.
Note: use bit-wise specification of the features:
bit value   meaning
  0     1   KN_VAR_LINEAR
default = 0 (variables are assumed to be nonlinear)

If a variable only appears linearly in the model, it can be very
helpful to mark this by enabling bit 0. This information can then
be used by Knitro to perform more extensive preprocessing. If a
variable appears nonlinearly in any constraint or the objective (or
if the user does not know) then it should not be marked as linear.
Variables are assumed to be nonlinear variables by default.
Knitro makes a local copy of all inputs, so the application may
free memory after the call.
"""
function KN_set_var_properties(m::Model, valindex::Vector{Cint}, xProperties::Vector{Cint})
    nvar = length(valindex)
    @assert nvar == length(xProperties)
    ret = @kn_ccall(set_var_properties, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cint}),
                    m.env, nvar, valindex, xProperties)
    _checkraise(ret)
end

function KN_set_var_properties(m::Model, xProperties::Vector{Cint})
    ret = @kn_ccall(set_var_properties_all, Cint, (Ptr{Cvoid}, Ptr{Cint}),
                    m.env, xProperties)
    _checkraise(ret)
end

function KN_set_var_property(m::Model, indexVar::Integer, xProperty::Integer)
    ret = @kn_ccall(set_var_property, Cint, (Ptr{Cvoid}, Cint, Cint),
                    m.env, indexVar, xProperty)
    _checkraise(ret)
end

##################################################
## Honor bounds
##################################################
"""
This API function can be used to identify which variables
should satisfy their variable bounds throughout the optimization
process (KN_HONORBNDS_ALWAYS).  The user option KN_PARAM_HONORBNDS
can be used to set ALL variables to honor their bounds.  This
routine takes precedence over the setting of KN_PARAM_HONORBNDS
and is used to customize the settings for individual variables.
Knitro makes a local copy of all inputs, so the application may
free memory after the call.
"""
function KN_set_var_honorbnds(m::Model, nindex::Integer, xHonorBound::Cint)
    ret = @kn_ccall(set_var_honorbnd, Cint,
                    (Ptr{Cvoid}, Cint, Cint),
                    m.env, nindex, xHonorBound)
    _checkraise(ret)
end

function KN_set_var_honorbnds(m::Model, valindex::Vector{Cint}, xHonorBound::Vector{Cint})
    nvar = length(valindex)
    @assert length(xHonorBound) == nvar
    ret = @kn_ccall(set_var_honorbnds, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cint}),
                    m.env, nvar, valindex, xHonorBound)
    _checkraise(ret)
end

function KN_set_var_honorbnds(m::Model, xHonorBound::Vector{Cint})
    ret = @kn_ccall(set_var_honorbnds_all, Cint, (Ptr{Cvoid}, Ptr{Cint}),
                    m.env, xHonorBound)
    _checkraise(ret)
end

##################################################
## Naming variables
##################################################
"""
Set names for model components passed in by the user/modeling
language so that Knitro can internally print out these names.
Knitro makes a local copy of all inputs, so the application may
free memory after the call.
"""
function KN_set_var_names(m::Model, nindex::Integer, name::String)
    ret = @kn_ccall(set_var_name, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cchar}),
                    m.env, nindex, name)
    _checkraise(ret)
end

function KN_set_var_names(m::Model, varIndex::Vector{Cint}, names::Vector{String})
    nvar = length(varIndex)
    ret = @kn_ccall(set_var_names, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Ptr{Char}}),
                    m.env, nvar, varIndex, names)
    _checkraise(ret)
end

function KN_set_var_names(m::Model, names::Vector{String})
    ret = @kn_ccall(set_var_names_all, Cint, (Ptr{Cvoid}, Ptr{Ptr{Cchar}}),
                    m.env, names)
    _checkraise(ret)
end

##################################################
## Initial values
##################################################
# primal values
"""
Set initial values for primal variables.  If not set, variables
may be initialized as 0 or initialized by Knitro based on some
initialization strategy (perhaps determined by a user option).
"""
function KN_set_var_primal_init_values(m::Model, xinitval::Vector{Cdouble})
    ret = @kn_ccall(set_var_primal_init_values_all, Cint,
                    (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, xinitval)
    _checkraise(ret)
end

function KN_set_var_primal_init_values(m::Model, indx::Integer, xinitval::Cdouble)
    ret = @kn_ccall(set_var_primal_init_value, Cint,
                    (Ptr{Cvoid}, Cint, Cdouble), m.env, indx, xinitval)
    _checkraise(ret)
end

# dual values
function KN_set_var_dual_init_values(m::Model, xinitval::Vector{Cdouble})
    ret = @kn_ccall(set_var_dual_init_values_all, Cint,
                    (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, xinitval)
    _checkraise(ret)
end

function KN_set_var_dual_init_values(m::Model, indx::Integer, xinitval::Cdouble)
    ret = @kn_ccall(set_var_dual_init_value, Cint,
                    (Ptr{Cvoid}, Cint, Cdouble), m.env, indx, xinitval)
    _checkraise(ret)
end

##################################################
## Scalings
##################################################
"""
Set an array of variable scaling and centering values to
perform a linear scaling
  x[i] = xScaleFactors[i] * xScaled[i] + xScaleCenters[i]
for each variable. These scaling factors should try to
represent the "typical" values of the "x" variables so that the
scaled variables ("xScaled") used internally by Knitro are close
to one.  The values for xScaleFactors should be positive.
If a non-positive value is specified, that variable will not
be scaled.
"""
function KN_set_var_scalings(m::Model, nindex::Integer,
                             xScaleFactors::Cdouble, xScaleCenters::Cdouble)
    ret = @kn_ccall(set_var_scaling, Cint,
                    (Ptr{Cvoid}, Cint, Cdouble, Cdouble),
                    m.env, nindex, xScaleFactors, xScaleCenters)
    _checkraise(ret)
end
KN_set_var_scalings(m::Model, nindex::Integer, xScaleFactors::Cdouble) =
    KN_set_var_scalings(m, nindex, xScaleFactors, 0.)

function KN_set_var_scalings(m::Model, valindex::Vector{Cint},
                             xScaleFactors::Vector{Cdouble}, xScaleCenters::Vector{Cdouble})
    nvar = length(valindex)
    @assert nvar = length(xScaleFactors) == length(xScaleCenters)
    ret = @kn_ccall(set_var_scalings, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}),
                    m.env, nvar, valindex, xScaleFactors, xScaleCenters)
    _checkraise(ret)
end
KN_set_var_scalings(m::Model, xIndex::Vector{Cint}, xScaleFactors::Vector{Cdouble}) =
    KN_set_var_scalings(m, xIndex, xScaleFactors, zeros(Cdouble, length(xScaleFactors)))

function KN_set_var_scalings(m::Model,
                             xScaleFactors::Vector{Cdouble}, xScaleCenters::Vector{Cdouble})
    @assert length(xScaleFactors) == length(xScaleCenters)
    ret = @kn_ccall(set_var_scalings_all, Cint, (Ptr{Cvoid}, Ptr{Cdouble}, Ptr{Cdouble}),
                    m.env, xScaleFactors, xScaleCenters)
    _checkraise(ret)
end
KN_set_var_scalings(m::Model, xScaleFactors::Vector{Cdouble}) =
    KN_set_var_scalings(m, xScaleFactors, zeros(Cdouble, length(xScaleFactors)))

##################################################
## Feasibility tolerance
##################################################
"""
Set custom absolute feasibility tolerances to use for the
termination tests.
The user options KN_PARAM_FEASTOL/KN_PARAM_FEASTOLABS define
a single tolerance that is applied equally to every constraint
and variable.  This API function allows the user to specify
separate feasibility termination tolerances for each constraint
and variable.  Values specified through this function will override
the value determined by KN_PARAM_FEASTOL/KN_PARAM_FEASTOLABS. The
tolerances should be positive values.  If a non-positive value is
specified, that constraint or variable will use the standard tolerances
based on  KN_PARAM_FEASTOL/KN_PARAM_FEASTOLABS.
The variables are considered to be satisfied when
    x[i] - xUpBnds[i] <= xFeasTols[i]  for all i=1..n, and
    xLoBnds[i] - x[i] <= xFeasTols[i]  for all i=1..n
The regular constraints are considered to be satisfied when
    c[i] - cUpBnds[i] <= cFeasTols[i]  for all i=1..m, and
    cLoBnds[i] - c[i] <= cFeasTols[i]  for all i=1..m
The complementarity constraints are considered to be satisfied when
    min(x1_i, x2_i) <= ccFeasTols[i]  for all i=1..ncc,
where x1 and x2 are the arrays of complementary pairs.
Knitro makes a local copy of all inputs, so the application
may free memory after the call.
"""
function KN_set_var_feastols(m::Model, nindex::Integer, xFeasTol::Cdouble)
    ret = @kn_ccall(set_var_feastol, Cint, (Ptr{Cvoid}, Cint, Cdouble),
                    m.env, nindex, xFeasTol)
    _checkraise(ret)
end

function KN_set_var_feastols(m::Model, xIndex::Vector{Cint}, xFeasTols::Vector{Cdouble})
    nvar = length(xIndex)
    @assert length(xFeasTols) == nvar
    ret = @kn_ccall(set_var_feastols, Cint,
                    (Ptr{Cvoid}, Cint, Ptr{Cint}, Ptr{Cdouble}),
                    m.env, nvar, xIndex, xFeasTols)
    _checkraise(ret)
end

function KN_set_var_feastols(m::Model, xFeasTols::Vector{Cdouble})
    ret = @kn_ccall(set_var_feastols_all, Cint, (Ptr{Cvoid}, Ptr{Cdouble}),
                    m.env, xFeasTols)
    _checkraise(ret)
end


text:
primitive type DiseaseState 8 end

DiseaseState(x::UInt8) = reinterpret(DiseaseState, x)
UInt8(x::DiseaseState) = reinterpret(UInt8, x)

# `DieaseState` const representations
const State_S = DiseaseState(0b0001)
const State_E = DiseaseState(0b0010)
const State_I = DiseaseState(0b0100)
const State_R = DiseaseState(0b1000)

const DiseaseStates = Vector{DiseaseState}

const _DiseaseStateVector = [State_S, State_E, State_I, State_R]
const _DiseaseStateCharVector = ['S', 'E', 'I', 'R']

# `DiseaseState` conversion to/from `Char`
function Base.convert(::Type{Char}, x::DiseaseState)
  return _DiseaseStateCharVector[_DiseaseStateVector .== Ref(x)][1]
end

function Base.show(io::IO, x::DiseaseState)
  return print(io, convert(Char, x))
end

Base.convert(::Type{DiseaseStates}, x::Type{SEIR}) = [State_S, State_E, State_I, State_R]
Base.convert(::Type{DiseaseStates}, x::Type{SEI}) = [State_S, State_E, State_I]
Base.convert(::Type{DiseaseStates}, x::Type{SIR}) = [State_S, State_I, State_R]
Base.convert(::Type{DiseaseStates}, x::Type{SI}) = [State_S, State_I]

Base.convert(::Type{Tuple}, x::Type{SEIR}) = (State_S, State_E, State_I, State_R)
Base.convert(::Type{Tuple}, x::Type{SEI}) = (State_S, State_E, State_I)
Base.convert(::Type{Tuple}, x::Type{SIR}) = (State_S, State_I, State_R)
Base.convert(::Type{Tuple}, x::Type{SI}) = (State_S, State_I)

Base.iterate(::Type{S}, v...) where {S <: DiseaseStateSequence} = iterate(convert(Tuple, S), v...)


text:
insts = OrderedDict()
insts["minimal"] = [
    ((2, 2, true),),
    ((2, 2, false),),
    ]
insts["fast"] = [
    ((2, 3, true),),
    ((2, 3, false),),
    ((5, 6, true),),
    ((5, 6, false),),
    ((10, 20, true),),
    ((10, 20, false),),
    ((20, 40, true),),
    ((20, 40, false),),
    ]
insts["various"] = [
    ((15, 20, true),),
    ((15, 20, false),),
    ((30, 40, true),),
    ((30, 40, false),),
    ((60, 80, true),),
    ((60, 80, false),),
    ((30, 160, true),),
    ((30, 160, false),),
    ]
return (MatrixQuadraticJuMP, insts)


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## rand!(::BitArray) && bitrand

function rand!(rng::AbstractRNG, B::BitArray)
    isempty(B) && return B
    Bc = B.chunks
    rand!(rng, Bc)
    Bc[end] &= Base._msk_end(B)
    return B
end

"""
    bitrand([rng=GLOBAL_RNG], [dims...])

Generate a `BitArray` of random boolean values.

# Examples
```jldoctest
julia> rng = MersenneTwister(1234);

julia> bitrand(rng, 10)
10-element BitArray{1}:
  true
  true
  true
 false
  true
 false
 false
  true
 false
  true
```
"""
bitrand(r::AbstractRNG, dims::Dims)   = rand!(r, BitArray(uninitialized, dims))
bitrand(r::AbstractRNG, dims::Integer...) = rand!(r, BitArray(uninitialized, convert(Dims, dims)))

bitrand(dims::Dims)   = rand!(BitArray(uninitialized, dims))
bitrand(dims::Integer...) = rand!(BitArray(uninitialized, convert(Dims, dims)))

## randstring (often useful for temporary filenames/dirnames)

"""
    randstring([rng=GLOBAL_RNG], [chars], [len=8])

Create a random string of length `len`, consisting of characters from
`chars`, which defaults to the set of upper- and lower-case letters
and the digits 0-9. The optional `rng` argument specifies a random
number generator, see [Random Numbers](@ref).

# Examples
```jldoctest
julia> srand(0); randstring()
"c03rgKi1"

julia> randstring(MersenneTwister(0), 'a':'z', 6)
"wijzek"

julia> randstring("ACGT")
"TATCGGTC"
```

!!! note
    `chars` can be any collection of characters, of type `Char` or
    `UInt8` (more efficient), provided [`rand`](@ref) can randomly
    pick characters from it.
"""
function randstring end

let b = UInt8['0':'9';'A':'Z';'a':'z']
    global randstring
    randstring(r::AbstractRNG, chars=b, n::Integer=8) = String(rand(r, chars, n))
    randstring(r::AbstractRNG, n::Integer) = randstring(r, b, n)
    randstring(chars=b, n::Integer=8) = randstring(GLOBAL_RNG, chars, n)
    randstring(n::Integer) = randstring(GLOBAL_RNG, b, n)
end

## randsubseq & randsubseq!

# Fill S (resized as needed) with a random subsequence of A, where
# each element of A is included in S with independent probability p.
# (Note that this is different from the problem of finding a random
#  size-m subset of A where m is fixed!)
function randsubseq!(r::AbstractRNG, S::AbstractArray, A::AbstractArray, p::Real)
    0 <= p <= 1 || throw(ArgumentError("probability $p not in [0,1]"))
    n = length(A)
    p == 1 && return copyto!(resize!(S, n), A)
    empty!(S)
    p == 0 && return S
    nexpected = p * length(A)
    sizehint!(S, round(Int,nexpected + 5*sqrt(nexpected)))
    if p > 0.15 # empirical threshold for trivial O(n) algorithm to be better
        for i = 1:n
            rand(r) <= p && push!(S, A[i])
        end
    else
        # Skip through A, in order, from each element i to the next element i+s
        # included in S. The probability that the next included element is
        # s==k (k > 0) is (1-p)^(k-1) * p, and hence the probability (CDF) that
        # s is in {1,...,k} is 1-(1-p)^k = F(k).   Thus, we can draw the skip s
        # from this probability distribution via the discrete inverse-transform
        # method: s = ceil(F^{-1}(u)) where u = rand(), which is simply
        # s = ceil(log(rand()) / log1p(-p)).
        # -log(rand()) is an exponential variate, so can use randexp().
        L = -1 / log1p(-p) # L > 0
        i = 0
        while true
            s = randexp(r) * L
            s >= n - i && return S # compare before ceil to avoid overflow
            push!(S, A[i += ceil(Int,s)])
        end
        # [This algorithm is similar in spirit to, but much simpler than,
        #  the one by Vitter for a related problem in "Faster methods for
        #  random sampling," Comm. ACM Magazine 7, 703-718 (1984).]
    end
    return S
end

"""
    randsubseq!(S, A, p)

Like [`randsubseq`](@ref), but the results are stored in `S`
(which is resized as needed).
"""
randsubseq!(S::AbstractArray, A::AbstractArray, p::Real) = randsubseq!(GLOBAL_RNG, S, A, p)

randsubseq(r::AbstractRNG, A::AbstractArray{T}, p::Real) where {T} =
    randsubseq!(r, T[], A, p)

"""
    randsubseq(A, p) -> Vector

Return a vector consisting of a random subsequence of the given array `A`, where each
element of `A` is included (in order) with independent probability `p`. (Complexity is
linear in `p*length(A)`, so this function is efficient even if `p` is small and `A` is
large.) Technically, this process is known as "Bernoulli sampling" of `A`.
"""
randsubseq(A::AbstractArray, p::Real) = randsubseq(GLOBAL_RNG, A, p)

## rand Less Than Masked 52 bits (helper function)

"Return a sampler generating a random `Int` (masked with `mask`) in ``[0, n)``, when `n <= 2^52`."
ltm52(n::Int, mask::Int=nextpow2(n)-1) = LessThan(n-1, Masked(mask, UInt52Raw(Int)))

## shuffle & shuffle!

"""
    shuffle!([rng=GLOBAL_RNG,] v::AbstractArray)

In-place version of [`shuffle`](@ref): randomly permute `v` in-place,
optionally supplying the random-number generator `rng`.

# Examples
```jldoctest
julia> rng = MersenneTwister(1234);

julia> shuffle!(rng, Vector(1:16))
16-element Array{Int64,1}:
  2
 15
  5
 14
  1
  9
 10
  6
 11
  3
 16
  7
  4
 12
  8
 13
```
"""
function shuffle!(r::AbstractRNG, a::AbstractArray)
    n = length(a)
    @assert n <= Int64(2)^52
    mask = nextpow2(n) - 1
    for i = n:-1:2
        (mask >> 1) == i && (mask >>= 1)
        j = 1 + rand(r, ltm52(i, mask))
        a[i], a[j] = a[j], a[i]
    end
    return a
end

shuffle!(a::AbstractArray) = shuffle!(GLOBAL_RNG, a)

"""
    shuffle([rng=GLOBAL_RNG,] v::AbstractArray)

Return a randomly permuted copy of `v`. The optional `rng` argument specifies a random
number generator (see [Random Numbers](@ref)).
To permute `v` in-place, see [`shuffle!`](@ref). To obtain randomly permuted
indices, see [`randperm`](@ref).

# Examples
```jldoctest
julia> rng = MersenneTwister(1234);

julia> shuffle(rng, Vector(1:10))
10-element Array{Int64,1}:
  6
  1
 10
  2
  3
  9
  5
  7
  4
  8
```
"""
shuffle(r::AbstractRNG, a::AbstractArray) = shuffle!(r, copymutable(a))
shuffle(a::AbstractArray) = shuffle(GLOBAL_RNG, a)

## randperm & randperm!

"""
    randperm([rng=GLOBAL_RNG,] n::Integer)

Construct a random permutation of length `n`. The optional `rng`
argument specifies a random number generator (see [Random Numbers](@ref)).
To randomly permute an arbitrary vector, see [`shuffle`](@ref)
or [`shuffle!`](@ref).

# Examples
```jldoctest
julia> randperm(MersenneTwister(1234), 4)
4-element Array{Int64,1}:
 2
 1
 4
 3
```
"""
randperm(r::AbstractRNG, n::Integer) = randperm!(r, Vector{Int}(uninitialized, n))
randperm(n::Integer) = randperm(GLOBAL_RNG, n)

"""
    randperm!([rng=GLOBAL_RNG,] A::Array{<:Integer})

Construct in `A` a random permutation of length `length(A)`. The
optional `rng` argument specifies a random number generator (see
[Random Numbers](@ref)). To randomly permute an arbitrary vector, see
[`shuffle`](@ref) or [`shuffle!`](@ref).

# Examples
```jldoctest
julia> randperm!(MersenneTwister(1234), Vector{Int}(uninitialized, 4))
4-element Array{Int64,1}:
 2
 1
 4
 3
```
"""
function randperm!(r::AbstractRNG, a::Array{<:Integer})
    n = length(a)
    @assert n <= Int64(2)^52
    n == 0 && return a
    a[1] = 1
    mask = 3
    @inbounds for i = 2:n
        j = 1 + rand(r, ltm52(i, mask))
        if i != j # a[i] is uninitialized (and could be #undef)
            a[i] = a[j]
        end
        a[j] = i
        i == 1+mask && (mask = 2mask + 1)
    end
    return a
end

randperm!(a::Array{<:Integer}) = randperm!(GLOBAL_RNG, a)

## randcycle & randcycle!

"""
    randcycle([rng=GLOBAL_RNG,] n::Integer)

Construct a random cyclic permutation of length `n`. The optional `rng`
argument specifies a random number generator, see [Random Numbers](@ref).

# Examples
```jldoctest
julia> randcycle(MersenneTwister(1234), 6)
6-element Array{Int64,1}:
 3
 5
 4
 6
 1
 2
```
"""
randcycle(r::AbstractRNG, n::Integer) = randcycle!(r, Vector{Int}(uninitialized, n))
randcycle(n::Integer) = randcycle(GLOBAL_RNG, n)

"""
    randcycle!([rng=GLOBAL_RNG,] A::Array{<:Integer})

Construct in `A` a random cyclic permutation of length `length(A)`.
The optional `rng` argument specifies a random number generator, see
[Random Numbers](@ref).

# Examples
```jldoctest
julia> randcycle!(MersenneTwister(1234), Vector{Int}(uninitialized, 6))
6-element Array{Int64,1}:
 3
 5
 4
 6
 1
 2
```
"""
function randcycle!(r::AbstractRNG, a::Array{<:Integer})
    n = length(a)
    n == 0 && return a
    @assert n <= Int64(2)^52
    a[1] = 1
    mask = 3
    @inbounds for i = 2:n
        j = 1 + rand(r, ltm52(i-1, mask))
        a[i] = a[j]
        a[j] = i
        i == 1+mask && (mask = 2mask + 1)
    end
    return a
end

randcycle!(a::Array{<:Integer}) = randcycle!(GLOBAL_RNG, a)

## random UUID generation

import Base: UUID, uuid_version

"""
    uuid1([rng::AbstractRNG=GLOBAL_RNG]) -> UUID

Generates a version 1 (time-based) universally unique identifier (UUID), as specified
by RFC 4122. Note that the Node ID is randomly generated (does not identify the host)
according to section 4.5 of the RFC.

# Examples
```jldoctest
julia> rng = MersenneTwister(1234);

julia> Random.uuid1(rng)
2cc938da-5937-11e7-196e-0f4ef71aa64b
```
"""
function uuid1(rng::AbstractRNG=GLOBAL_RNG)
    u = rand(rng, UInt128)

    # mask off clock sequence and node
    u &= 0x00000000000000003fffffffffffffff

    # set the unicast/multicast bit and version
    u |= 0x00000000000010000000010000000000

    # 0x01b21dd213814000 is the number of 100 nanosecond intervals
    # between the UUID epoch and Unix epoch
    timestamp = round(UInt64, time() * 1e7) + 0x01b21dd213814000
    ts_low = timestamp & typemax(UInt32)
    ts_mid = (timestamp >> 32) & typemax(UInt16)
    ts_hi = (timestamp >> 48) & 0x0fff

    u |= UInt128(ts_low) << 96
    u |= UInt128(ts_mid) << 80
    u |= UInt128(ts_hi) << 64

    UUID(u)
end

"""
    uuid4([rng::AbstractRNG=GLOBAL_RNG]) -> UUID

Generates a version 4 (random or pseudo-random) universally unique identifier (UUID),
as specified by RFC 4122.

# Examples
```jldoctest
julia> rng = MersenneTwister(1234);

julia> Random.uuid4(rng)
82015f10-44cc-4827-996e-0f4ef71aa64b
```
"""
function uuid4(rng::AbstractRNG=GLOBAL_RNG)
    u = rand(rng, UInt128)
    u &= 0xffffffffffff0fff3fffffffffffffff
    u |= 0x00000000000040008000000000000000
    UUID(u)
end

max_stars_repo_name: ashaywalke/julia

max_issues_repo_name: ashaywalke/julia

max_forks_repo_name: ashaywalke/julia


text:
# ---
# title: 1157. Online Majority Element In Subarray
# id: problem1157
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Array, Binary Search, Segment Tree
# link: <https://leetcode.com/problems/online-majority-element-in-subarray/description/>
# hidden: true
# ---
#
# Implementing the class `MajorityChecker`, which has the following API:
#
#   * `MajorityChecker(int[] arr)` constructs an instance of MajorityChecker with the given array `arr`;
#   * `int query(int left, int right, int threshold)` has arguments such that:
#     * `0 <= left <= right < arr.length` representing a subarray of `arr`;
#     * `2 * threshold > right - left + 1`, ie. the threshold is always a strict majority of the length of the subarray
#
# Each `query(...)` returns the element in `arr[left], arr[left+1], ...,
# arr[right]` that occurs at least `threshold` times, or `-1` if no such element
# exists.
#
#
#
# **Example:**
#
#
#
#     MajorityChecker majorityChecker = new MajorityChecker([1,1,2,2,1,1]);
#     majorityChecker.query(0,5,4); // returns 1
#     majorityChecker.query(0,3,3); // returns -1
#     majorityChecker.query(2,3,2); // returns 2
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 20000`
#   * `1 <= arr[i] <= 20000`
#   * For each query, `0 <= left <= right < len(arr)`
#   * For each query, `2 * threshold > right - left + 1`
#   * The number of queries is at most `10000`
#
#
## @lc code=start
using LeetCode

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


text:
struct Foo
    t::Vector{DateTime}
    x::Vector{Float64}
end

ts = DateTime(2018,01,13):Minute(1):DateTime(2018,01,14)
x = randn(length(ts))

f = Foo(ts,x)

DaylightPlots.tlims(f::Foo) = DaylightPlots.tlims(f.t,f.x)

@recipe f(f::Foo) = (f.t,f.x)

daylight(f)


text:
"""
    fix_age(x::AbstractString)::Int

Return age after parsing `x`.

```jldoctest
julia> fix_age = Codex.Questionnaires.fix_age;

julia> fix_age("27-jun-93")
26

julia> fix_age("1993")
26

julia> fix_age("23 jaar")
23
```
"""
function fix_age(x::AbstractString)::Int
    x = string(x)
    if endswith(x, "jaar")
        x = x[1:end-4]
    end
    value = tryparse(Int, x)
    if !isnothing(value) && 18 < value && value < 65
        return value
    elseif !isnothing(value) && 65 < value
        return 2019 - value
    elseif x == "27-jun-93"
        return 2019 - 1993
    else
        throw(MissingException("No known fix for the age $x"))
    end
end

fix_age(x::Int) = x

"""
    unify_demographics(df)

Returns an simplified and unified DataFrame which is the same for 2018, 2019 and 2020.
It may throw out some data which we don't need at the time of writing.
"""
function unify_demographics(df)
    df.age = fix_age.(df.v1)
    select!(df, :id, :completed_at, :age)
    df
end


text:
using DataInterpolations, Test

@testset "Interface" begin include("interface.jl") end
@testset "Interpolation Tests" begin include("interpolation_tests.jl") end
@testset "Derivative Tests" begin include("derivative_tests.jl") end
@testset "Integral Tests" begin include("integral_tests.jl") end
@testset "Online Tests" begin include("online_tests.jl") end
@testset "Regularization Smoothing" begin include("regularization.jl") end


text:
macro debugcb(ex)
    seen = false
    ex = MacroTools.prewalk(ex) do x
        if !seen && MacroTools.@capture(MacroTools.shortdef(x), f_(xs__) = body_)
            seen = true
            def = MacroTools.splitdef(x)
            debug_ex = :($Base.@info $f)
            push!(debug_ex.args, map(MacroTools.namify, def[:args])...)
            # debug_ex = :(println($f, $(map(MacroTools.namify, def[:args])...)))
            pushfirst!(def[:body].args, debug_ex)
            return MacroTools.combinedef(def)
        else
            return x
        end
    end
    return esc(ex)
end


text:
using CSV
using DataFrames
using Plots
include("/Users/anguslewis/Documents/SFFMProject/DFQExamples.jl/plot_utils.jl")

@halfwidth_plot_defaults()

plotlyjs()
function p()
ks_data = CSV.read((@__DIR__)*"/../data/meshs_ks_func_count_2.csv",DataFrame)

plot()
linestyles_vec = [:solid,:dash,:dashdot]
markerstyles_vec = [:cross,:diamond,:circle]
colours = [1;2;4]
for (c,col) in enumerate(names(ks_data))
    plot!(log10.(1:2:21),ks_data[:,col],
        label=col,linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],color=colours[c],
        linewidth=2)
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="KS error")
# plot!(legend=(0.7,0.9))
plot!()
@add_lines!(ks_data,("Unif","QBDRAP"),@__DIR__)
savefig((@__DIR__)*"/meshs_ks_error_formatted.pdf")
end
p()

function p()
l2_pdf_data = CSV.read((@__DIR__)*"/../data/meshs_l2_pdf_func_count_2.csv",DataFrame)
colours = [1;2;4]

plot()
linestyles_vec = [:solid,:dash,:dashdot]
markerstyles_vec = [:cross,:diamond,:circle]
for (c,col) in enumerate(names(l2_pdf_data))
    plot!(log10.(1:2:21),l2_pdf_data[:,col],
        label=col,linestyle=linestyles_vec[c],
        marker=markerstyles_vec[c],color=colours[c],
        linewidth=2)
end
plot!(xlabel="Dimension")
plot!(ylabel="Error"); error_ticks!(plot!())
plot!(title="L² error - PDF")
# plot!(legend=(0.7,0.9))
plot!()
@add_lines!(l2_pdf_data,("Unif","QBDRAP"),@__DIR__)
savefig((@__DIR__)*"/meshs_l2_pdf_error_formatted.pdf")
end
p()

function p()
@fullwidth_plot_defaults()

os = 1:2:7
plot(layout=(1,length(os)))
linestyles_vec = [:solid,:dash,:dashdot]
markerstyles_vec = [:cross,:diamond,:circle]
xs = range(eps(),1-eps(),length=200)
for (c,o) in enumerate(os)
    model = BoundedFluidQueue(fill(0.0,1,1),[1.0],1.0)

    mesh = DGMesh(0.0:1.0,o)
    dq = DiscretisedFluidQueue(model,mesh)
    d0 = SFMDistribution((x,i)->2*Float64(x<=0.5),dq,TrapezoidRule;fun_evals=10_001)
    plot!(xs,cdf(d0).(xs,1),subplot=c,label=(c==1 ? "DG" : false),
        linestyle=linestyles_vec[1],
        # marker=markerstyles_vec[1],
        linewidth=2)

    mesh = DGMesh(0.0:(1.0./o):1.0,1)
    dq = DiscretisedFluidQueue(model,mesh)
    d0 = SFMDistribution((x,i)->2*Float64(x<=0.5),dq,TrapezoidRule;fun_evals=10_001)
    plot!(xs,cdf(d0).(xs,1),subplot=c,label=(c==1 ? "unif." : false),
        linestyle=linestyles_vec[2],
        # marker=markerstyles_vec[2],
        linewidth=2)

    mesh = FRAPMesh(0.0:1.0,o)
    dq = DiscretisedFluidQueue(model,mesh)
    d0 = SFMDistribution((x,i)->2*Float64(x<=0.5),dq,TrapezoidRule;fun_evals=10_001)
    plot!(xs,cdf(d0).(xs,1),subplot=c,label=(c==1 ? "QBD-RAP" : false),
        linestyle=linestyles_vec[3],
        # marker=markerstyles_vec[3],
        linewidth=2)

    plot!(xs,x->2*x*(x<0.5) + Float64(x>=0.5),subplot=c,label=(c==1 ? "truth" : false),
        linestyle=:dot,
        # marker=markerstyles_vec[3],
        color=:grey,
        linewidth=2)

    plot!(subplot=c, yticks=(c==1 ? (0:0.2:1.2) : false),
        ylims = (-0.1,1.1), xticks = 0.25:0.5:0.75,
        ylabel=(c==1 ? "probability" : false), grid=false,
        title = string("Dim. ",o),
        legend=(0.93,0.4))
end
plot!()
savefig((@__DIR__)*"/cdfs_formatted.pdf")
end
p()

function p()
os = 1:2:7
plot(layout=(1,length(os)))
linestyles_vec = [:solid,:dash,:dashdot]
markerstyles_vec = [:cross,:diamond,:circle]
xs = range(eps(),1-eps(),length=200)
for (c,o) in enumerate(os)
    model = BoundedFluidQueue(fill(0.0,1,1),[1.0],1.0)

    mesh = DGMesh(0.0:1.0,o)
    dq = DiscretisedFluidQueue(model,mesh)
    d0 = SFMDistribution((x,i)->2*Float64(x<=0.5),dq,TrapezoidRule;fun_evals=10_001)
    plot!(xs,pdf(d0).(xs,1),subplot=c,label=(c==1 ? "DG" : false),
        linestyle=linestyles_vec[1],
        # marker=markerstyles_vec[1],
        linewidth=2)

    mesh = DGMesh(0.0:(1.0./o):1.0,1)
    dq = DiscretisedFluidQueue(model,mesh)
    d0 = SFMDistribution((x,i)->2*Float64(x<=0.5),dq,TrapezoidRule;fun_evals=10_001)
    plot!(xs,pdf(d0).(xs,1),subplot=c,label=(c==1 ? "unif." : false),
        linestyle=linestyles_vec[2],
        # marker=markerstyles_vec[2],
        linewidth=2)

    mesh = FRAPMesh(0.0:1.0,o)
    dq = DiscretisedFluidQueue(model,mesh)
    d0 = SFMDistribution((x,i)->2*Float64(x<=0.5),dq,TrapezoidRule;fun_evals=10_001)
    plot!(xs,pdf(d0).(xs,1),subplot=c,label=(c==1 ? "QBD-RAP" : false),
        linestyle=linestyles_vec[3],
        # marker=markerstyles_vec[3],
        linewidth=2)

    plot!(xs,x->2*Float64(x<=0.5),subplot=c,label=(c==1 ? "truth" : false),
        linestyle=:dot,
        # marker=markerstyles_vec[3],
        color=:grey,
        linewidth=2)

    plot!(subplot=c, yticks=(c==1 ? (-2:2:6) : false),
        ylims = (-3,7.2), xticks = 0.25:0.5:0.75,
        ylabel=(c==1 ? "probability" : false), grid=false,
        title = string("Dim. ",o),
        legend=(0.1,0.8))
end
plot!()
savefig((@__DIR__)*"/pdfs_formatted.pdf")
end
p()


text:
module TestMeasures

using MLJBase, Test
import Distributions
using CategoricalArrays
using Statistics
import LossFunctions
using StableRNGs
using OrderedCollections: LittleDict

rng  = StableRNGs.StableRNG(123)

@testset "aggregation" begin
    v = rand(5)
    @test aggregate(v, mae) ≈ mean(v)
    @test aggregate(v, TruePositive()) ≈ sum(v)
    @test aggregate(v, rms) ≈ sqrt(mean(v.^2))
    λ = rand()
    @test aggregate(λ, rms) === λ
    @test aggregate(aggregate(v, l2), l2) == aggregate(v, l2)
    m = LittleDict([0, 1, 2, 3, 4], v)
    @test aggregate(m, MTPR()) == mean(v)
end

@testset "metadata" begin
    measures()
    measures(m -> m.target_scitype <: AbstractVector{<:Finite} &&
             m.supports_weights)
    info(rms)
    @test true
end

@testset "coverage" begin
    # just checking that the traits work not that they're correct
    @test orientation(BrierScore()) == :score
    @test orientation(auc) == :score
    @test orientation(rms) == :loss

    @test reports_each_observation(auc) == false
    @test is_feature_dependent(auc) == false

    @test MLJBase.distribution_type(auc) == MLJBase.UnivariateFinite
end

@testset "MLJBase.value" begin
    yhat = randn(rng,5)
    X = (weight=randn(rng,5), x1 = randn(rng,5))
    y = randn(rng,5)
    w = randn(rng,5)

    @test MLJBase.value(mae, yhat, nothing, y, nothing) ≈ mae(yhat, y)
    @test MLJBase.value(mae, yhat, nothing, y, w) ≈ mae(yhat, y, w)

    spooky(yhat, y) = abs.(yhat - y) |> mean
    @test MLJBase.value(spooky, yhat, nothing, y, nothing) ≈ mae(yhat, y)

    cool(yhat, y, w) = abs.(yhat - y) .* w |> mean
    MLJBase.supports_weights(::Type{typeof(cool)}) = true
    @test MLJBase.value(cool, yhat, nothing, y, w) ≈ mae(yhat, y, w)

    funky(yhat, X, y) = X.weight .* abs.(yhat - y) |> mean
    MLJBase.is_feature_dependent(::Type{typeof(funky)}) = true
    @test MLJBase.value(funky, yhat, X, y, nothing) ≈ mae(yhat, y, X.weight)

    weird(yhat, X, y, w) = w .* X.weight .* abs.(yhat - y) |> mean
    MLJBase.is_feature_dependent(::Type{typeof(weird)}) = true
    MLJBase.supports_weights(::Type{typeof(weird)}) = true
    @test MLJBase.value(weird, yhat, X, y, w) ≈ mae(yhat, y, X.weight .* w)
end

mutable struct DRegressor <: Deterministic end
MLJBase.target_scitype(::Type{<:DRegressor}) =
    AbstractVector{<:Continuous}

mutable struct D2Regressor <: Deterministic end
MLJBase.target_scitype(::Type{<:D2Regressor}) =
    AbstractVector{Continuous}

mutable struct DClassifier <: Deterministic end
MLJBase.target_scitype(::Type{<:DClassifier}) =
    AbstractVector{<:Finite}

mutable struct PClassifier <: Probabilistic end
MLJBase.target_scitype(::Type{<:PClassifier}) =
    AbstractVector{<:Finite}

mutable struct PRegressor <: Probabilistic end
MLJBase.target_scitype(::Type{<:PRegressor}) =
    AbstractVector{<:Continuous}

mutable struct PCountRegressor <: Probabilistic end
MLJBase.target_scitype(::Type{<:PCountRegressor}) =
    AbstractVector{<:Count}

@testset "default_measure" begin
    @test MLJBase.default_measure(DRegressor()) == rms
    @test MLJBase.default_measure(D2Regressor()) == rms
    @test MLJBase.default_measure(DClassifier()) == misclassification_rate
    @test MLJBase.default_measure(PClassifier()) == log_loss

    @test MLJBase.default_measure(DRegressor) == rms
    @test MLJBase.default_measure(D2Regressor) == rms
    @test MLJBase.default_measure(DClassifier) == misclassification_rate
    @test MLJBase.default_measure(PClassifier) == log_loss

    @test MLJBase.default_measure(PRegressor) == log_loss
    @test MLJBase.default_measure(PCountRegressor) == log_loss
end

include("confusion_matrix.jl")
include("roc.jl")
include("continuous.jl")
include("finite.jl")
include("probabilistic.jl")
include("loss_functions_interface.jl")

@testset "show method for measures" begin
    io = IOBuffer()
    for meta in measures()
        m = eval(Meta.parse("$(meta.name)()"))
        show(io, MIME("text/plain"), m)
    end
end

@testset "missing and NaN values in aggregation" begin
    v =[1, 2, missing, 5, NaN]
    @test MLJBase.Sum()(v) == 8
    @test MLJBase.RootMeanSquare()(v) ≈ sqrt((1 + 4 + 25)/3)
    @test_throws MLJBase.ERR_NOTHING_LEFT_TO_AGGREGATE MLJBase.Mean()(Float32[])
end

end
true


text:
#### EDMF Functions

bound(x, x_bounds) = min(max(x, x_bounds[1]), x_bounds[2])

inside_bounds(x, x_bounds) = x > x_bounds[1] && x < x_bounds[2]

stable(obukhov_length::FT) where {FT} = obukhov_length > 0
unstable(obukhov_length::FT) where {FT} = obukhov_length < 0

"""
    StabilityDependentParam

A parameter that has two values:
 - `stable` for stable conditions
 - `unstable` for unstable conditions
"""
struct StabilityDependentParam{FT}
    stable::FT
    unstable::FT
end

function (sdp::StabilityDependentParam{FT})(obukhov_length::FT) where {FT}
    return unstable(obukhov_length) ? sdp.unstable : sdp.stable
end

function compute_tendencies_en_O2!(
    grid::Grid{FT},
    q_tendencies,
    q,
    aux_O2,
    cv,
    aux,
    tri_diag,
    params,
) where {FT}
    gm, en, ud, sd, al = allcombinations(q_tendencies)
    k_1 = first_interior(grid, Zmin())
    @inbounds for k in over_elems_real(grid)
        q_tendencies[cv, k, en] =
            aux_O2[cv][:press, k] +
            aux_O2[cv][:buoy, k] +
            aux_O2[cv][:shear, k] +
            aux_O2[cv][:entr_gain, k]
        # aux_O2[cv][:rain_src, k]
    end
    q_tendencies[cv, k_1, en] = FT(0)

    construct_tridiag_diffusion_O2!(grid, q, aux, params, tri_diag[cv])
    k_1 = first_interior(grid, Zmin())
    @inbounds for k in over_elems(grid)
        tri_diag[cv][:f, k] =
            aux[:ρ_0, k] * q[:a, k, en] * q[cv, k, en] * (1 / params[:Δt][1]) +
            q_tendencies[cv, k, en]
    end
    tri_diag[cv][:f, k_1] =
        aux[:ρ_0, k_1] *
        q[:a, k_1, en] *
        q[cv, k_1, en] *
        (1 / params[:Δt][1]) + q[cv, k_1, en]
end

function compute_tendencies_gm_scalars!(
    grid,
    q_tendencies,
    q,
    aux,
    params,
    tri_diag,
)
    gm, en, ud, sd, al = allcombinations(q)
    @unpack SurfaceModel = params
    k_1 = first_interior(grid, Zmin())
    Δzi = grid.Δzi
    α_1 = aux[:α_0, k_1]
    ae_1 = q[:a, k_1, en]
    Δt = params[:Δt][1]
    @inbounds for k in over_elems(grid)
        q_tendencies[:q_tot, k, gm] += aux[:mf_tend_q_tot, k]
        q_tendencies[:θ_liq, k, gm] += aux[:mf_tend_θ_liq, k]
    end

    q_tendencies[:q_tot, k_1, gm] += SurfaceModel.ρq_tot_flux * Δzi * α_1 / ae_1
    q_tendencies[:θ_liq, k_1, gm] += SurfaceModel.ρθ_liq_flux * Δzi * α_1 / ae_1

    @inbounds for k in over_elems_real(grid)
        tri_diag[:q_tot][:ρaK, k] =
            q[:a, k, en] * aux[:K_h, k, gm] * aux[:ρ_0, k]
        tri_diag[:θ_liq][:ρaK, k] =
            q[:a, k, en] * aux[:K_h, k, gm] * aux[:ρ_0, k]
    end
    construct_tridiag_diffusion_O1!(grid, q, aux, Δt, tri_diag[:q_tot])
    construct_tridiag_diffusion_O1!(grid, q, aux, Δt, tri_diag[:θ_liq])

    @inbounds for k in over_elems(grid)
        tri_diag[:q_tot][:f, k] =
            q[:q_tot, k, gm] + Δt * q_tendencies[:q_tot, k, gm]
    end

    @inbounds for k in over_elems(grid)
        tri_diag[:θ_liq][:f, k] =
            q[:θ_liq, k, gm] + Δt * q_tendencies[:θ_liq, k, gm]
    end
end

function compute_tendencies_ud!(grid, q_tendencies, q, aux, params)
    gm, en, ud, sd, al = allcombinations(q)
    @inbounds for i in ud
        @inbounds for k in over_elems_real(grid)
            a_k = q[:a, k, i]
            α_0_kp = aux[:α_0, k]
            w_env = q[:w, k, en]
            ρ_k = aux[:ρ_0, k]
            w_i = q[:w, k, i]
            ε_model = aux[:ε_model, k, i]
            δ_model = aux[:δ_model, k, i]
            θ_liq_env = q[:θ_liq, k, en]
            q_tot_env = q[:q_tot, k, en]
            B_k = aux[:buoy, k, i]
            ρa_k = ρ_k * a_k
            ρaw_k = ρa_k * w_i

            a_cut = q[:a, Cut(k), i]
            θ_liq_cut = q[:θ_liq, Cut(k), i]
            q_tot_cut = q[:q_tot, Cut(k), i]
            w_cut = q[:w, Cut(k), i]
            ρ_cut = aux[:ρ_0, Cut(k)]

            ρaw_cut = ρ_cut .* a_cut .* w_cut
            ρawθ_liq_cut = ρaw_cut .* θ_liq_cut
            ρawq_tot_cut = ρaw_cut .* q_tot_cut
            ρaww_cut = ρ_cut .* a_cut .* w_cut .* w_cut

            tendencies = 0
            adv = -α_0_kp * advect(ρaw_cut, w_cut, grid)
            tendencies += adv
            ε_term = a_k * w_i * ε_model
            tendencies += ε_term
            δ_term = -a_k * w_i * δ_model
            tendencies += δ_term
            q_tendencies[:a, k, i] = tendencies

            adv = -advect(ρaww_cut, w_cut, grid)
            exch = ρaw_k * (-δ_model * w_i + ε_model * w_env)
            buoy = ρa_k * B_k
            nh_press = aux[:nh_press, k, i]

            tendencies = (adv + exch + buoy + nh_press)
            q_tendencies[:w, k, i] = tendencies

            tendencies_θ_liq = 0
            tendencies_q_tot = 0

            tendencies_θ_liq += -advect(ρawθ_liq_cut, w_cut, grid)
            tendencies_q_tot += -advect(ρawq_tot_cut, w_cut, grid)

            tendencies_θ_liq +=
                ρaw_k * (ε_model * θ_liq_env - δ_model * θ_liq_cut[2])
            tendencies_q_tot +=
                ρaw_k * (ε_model * q_tot_env - δ_model * q_tot_cut[2])

            q_tendencies[:θ_liq, k, i] = tendencies_θ_liq
            q_tendencies[:q_tot, k, i] = tendencies_q_tot
        end
    end
end

function saturation_adjustment_sd!(grid, q, aux, params)
    gm, en, ud, sd, al = allcombinations(q)
    @unpack param_set = params
    @inbounds for i in sd
        @inbounds for k in over_elems_real(grid)
            ts = ActiveThermoState(param_set, q, aux, k, i)
            q_liq = PhasePartition(ts).liq
            T = air_temperature(ts)
            aux[:T, k, i] = T
            aux[:q_liq, k, i] = q_liq
        end
    end
end

include(joinpath("Models", "entr_detr.jl"))
include(joinpath("Models", "buoyancy.jl"))
include(joinpath("Models", "mixing_length.jl"))
include(joinpath("Models", "microphysics.jl"))
include(joinpath("Models", "surface.jl"))
include(joinpath("Models", "pressure.jl"))
include(joinpath("Models", "eddy_diffusivity.jl"))

function top_of_updraft(grid::Grid, q::StateVec, params)
    @unpack w_bounds, a_bounds = params
    gm, en, ud, sd, al = allcombinations(q)

    z_star_a = zeros(length(ud))
    z_star_w = zeros(length(ud))
    k_2 = first_interior(grid, Zmax())
    @inbounds for i in ud
        z_star_a[i] = min([
            !inside_bounds(q[:a, k, i], a_bounds) ? grid.zc[k] :
                grid.zc[k_2 + 1] for k in over_elems_real(grid)
        ]...)
        z_star_w[i] = min([
            !inside_bounds(q[:w, k, i], w_bounds) ? grid.zc[k] :
                grid.zc[k_2 + 1] for k in over_elems_real(grid)
        ]...)
    end
    return z_star_a, z_star_w
end

function filter_scalars!(grid, q, aux, params)
    gm, en, ud, sd, al = allcombinations(q)
    z_star_a, z_star_w = top_of_updraft(grid, q, params)
    @inbounds for i in ud
        @inbounds for k in over_elems_real(grid)
            aux[:HVSD_a, k, i] = 1 - heaviside(grid.zc[k] - z_star_a[i], 1)
            aux[:HVSD_w, k, i] = 1 - heaviside(grid.zc[k] - z_star_w[i], 1)
        end
    end

    @inbounds for i in ud
        @inbounds for k in over_elems_real(grid)[2:end]
            q[:w, k, i] =
                bound(q[:w, k, i] * aux[:HVSD_w, k, i], params[:w_bounds])
            q[:a, k, i] =
                bound(q[:a, k, i] * aux[:HVSD_w, k, i], params[:a_bounds])

            weight = aux[:HVSD_w, k, i]
            q[:θ_liq, k, i] =
                weight * q[:θ_liq, k, i] + (1 - weight) * q[:θ_liq, k, gm]
            q[:q_tot, k, i] =
                weight * q[:q_tot, k, i] + (1 - weight) * q[:q_tot, k, gm]
        end
    end
end

function compute_cv_gm!(grid, q, ϕ, ψ, cv, tke_factor)
    gm, en, ud, sd, al = allcombinations(q)
    @inbounds for k in over_elems(grid)
        Δϕ = q[ϕ, k, en] - q[ϕ, k, gm]
        Δψ = q[ψ, k, en] - q[ψ, k, gm]
        q[cv, k, en] =
            tke_factor * q[:a, k, en] * Δϕ * Δψ + q[:a, k, en] * q[cv, k, en]
        @inbounds for i in ud
            Δϕ = q[ϕ, k, i] - q[ϕ, k, gm]
            Δψ = q[ψ, k, i] - q[ψ, k, gm]
            q[cv, k, en] += tke_factor * q[:a, k, i] * Δϕ * Δψ
        end
    end
end

function compute_mf_gm!(grid, q, aux)
    gm, en, ud, sd, al = allcombinations(q)
    domain_c = over_elems_real(grid)

    @inbounds for i in ud
        @inbounds for k in domain_c
            aux[:mf_aux, k, i] =
                (q[:w, k, i] - q[:w, k, en]) * aux[:ρ_0, k] * q[:a, k, i]
        end
    end

    @inbounds for k in domain_c
        aux[:mf_θ_liq, k] = sum([
            aux[:mf_aux, k, i] * (q[:θ_liq, k, i] - q[:θ_liq, k, en])
            for i in ud
        ])
        aux[:mf_q_tot, k] = sum([
            aux[:mf_aux, k, i] * (q[:q_tot, k, i] - q[:q_tot, k, en])
            for i in ud
        ])
    end

    @inbounds for k in domain_c
        aux[:mf_tend_θ_liq, k] =
            -aux[:α_0, k] * grad(aux[:mf_θ_liq, Cut(k)], grid)
        aux[:mf_tend_q_tot, k] =
            -aux[:α_0, k] * grad(aux[:mf_q_tot, Cut(k)], grid)
    end
end

function compute_cv_shear!(grid::Grid{FT}, q, aux, aux_O2, ϕ, ψ, cv) where {FT}
    gm, en, ud, sd, al = allcombinations(q)
    is_tke = cv == :tke
    tke_factor = is_tke ? FT(0.5) : 1
    grad_u = 0
    grad_v = 0
    @inbounds for k in over_elems_real(grid)
        if is_tke
            grad_u = ∇_dn(q[:u, Cut(k), gm], grid)
            grad_v = ∇_dn(q[:v, Cut(k), gm], grid)
            grad_ϕ = ∇_dn(q[ϕ, Cut(k), en], grid)
            grad_ψ = ∇_dn(q[ψ, Cut(k), en], grid)
        else
            grad_ϕ = grad(q[ϕ, Cut(k), en], grid)
            grad_ψ = grad(q[ψ, Cut(k), en], grid)
        end
        ρaK = aux[:ρ_0, k] * q[:a, k, en] * aux[:K_h, k, gm]
        aux_O2[cv][:shear, k] =
            tke_factor * 2 * ρaK * (grad_ϕ * grad_ψ + grad_u^2 + grad_v^2)
    end
end

function compute_cv_interdomain_src!(
    grid::Grid{FT},
    q,
    aux,
    aux_O2,
    ϕ,
    ψ,
    cv,
    tke_factor,
) where {FT}
    gm, en, ud, sd, al = allcombinations(q)
    @inbounds for k in over_elems(grid)
        aux_O2[cv][:interdomain, k] = FT(0)
        @inbounds for i in ud
            Δϕ = q[ϕ, k, i] - q[ϕ, k, en]
            Δψ = q[ψ, k, i] - q[ψ, k, en]
            aux_O2[cv][:interdomain, k] +=
                tke_factor * q[:a, k, i] * (1 - q[:a, k, i]) * Δϕ * Δψ
        end
    end
end

function compute_cv_env!(
    grid::Grid{FT},
    q,
    aux,
    aux_O2,
    ϕ,
    ψ,
    cv,
    tke_factor,
) where {FT}
    gm, en, ud, sd, al = allcombinations(q)
    @inbounds for k in over_elems(grid)
        if q[:a, k, en] > 0
            @inbounds for i in sd
                Δϕ = q[ϕ, k, i] - q[ϕ, k, gm]
                Δψ = q[ψ, k, i] - q[ψ, k, gm]
                q[cv, k, en] -= tke_factor * q[:a, k, i] * Δϕ * Δψ
            end
            q[cv, k, en] = q[cv, k, en] / q[:a, k, en]
        else
            q[cv, k, en] = FT(0)
        end
    end
end

function cleanup_covariance!(grid::Grid{FT}, q) where {FT}
    gm, en, ud, sd, al = allcombinations(q)
    @inbounds for k in over_elems_real(grid)
        q[:tke, k, en] = q[:tke, k, en] < eps(FT) ? FT(0) : q[:tke, k, en]
    end
end

function construct_tridiag_diffusion_O1!(
    grid::Grid{FT},
    q,
    aux,
    Δt,
    tri_diag,
) where {FT}
    gm, en, ud, sd, al = allcombinations(aux)
    k_1 = first_interior(grid, Zmin())
    k_2 = first_interior(grid, Zmax())
    dzi = grid.Δzi
    @inbounds for k in over_elems_real(grid)
        ρaK_cut = tri_diag[:ρaK, Cut(k)]
        X = aux[:ρ_0, k] * q[:a, k, en] / Δt
        Z = ρaK_cut[1] * dzi * dzi
        Y = ρaK_cut[2] * dzi * dzi
        if k == k_1
            Z = FT(0)
        elseif k == k_2
            Y = FT(0)
        end
        tri_diag[:a, k] = -Z / X
        tri_diag[:b, k] = 1 + Y / X + Z / X
        tri_diag[:c, k] = -Y / X
    end
end

function construct_tridiag_diffusion_O2!(
    grid::Grid{FT},
    q,
    aux,
    params,
    tri_diag,
) where {FT}
    gm, en, ud, sd, al = allcombinations(q)
    Δzi = grid.Δzi
    Δzi2 = grid.Δzi2
    dti = 1 / params[:Δt][1]
    k_1 = first_interior(grid, Zmin())
    k_2 = first_interior(grid, Zmax())

    @inbounds for k in over_elems_real(grid)
        ρ_0_cut = aux[:ρ_0, Cut(k)]
        ae_cut = q[:a, Cut(k), en]
        w_cut = q[:w, Cut(k), en]
        ρa_K = q[:a, Cut(k), en] .* aux[:K_h, Cut(k), gm] .* aux[:ρ_0, Cut(k)]

        D_env = sum([
            ρ_0_cut[2] * q[:a, k, i] * q[:w, k, i] * aux[:ε_model, k, i]
            for i in ud
        ])

        l_mix = max(aux[:l_mix, k, gm], FT(1))
        tke_env = max(q[:tke, k, en], FT(0))

        tri_diag[:a, k] = (-ρa_K[2] * Δzi2)
        tri_diag[:b, k] = (
            ρ_0_cut[2] * ae_cut[2] * dti -
            ρ_0_cut[2] * ae_cut[2] * w_cut[2] * Δzi +
            ρa_K[3] * Δzi2 +
            ρa_K[2] * Δzi2 +
            D_env +
            ρ_0_cut[2] * ae_cut[2] * params[:tke_diss_coeff] * sqrt(tke_env) / l_mix
        )
        tri_diag[:c, k] =
            (ρ_0_cut[3] * ae_cut[3] * w_cut[3] * Δzi - ρa_K[3] * Δzi2)
    end

    tri_diag[:a, k_1] = FT(0)
    tri_diag[:b, k_1] = FT(1)
    tri_diag[:c, k_1] = FT(0)

    tri_diag[:b, k_2] += tri_diag[:c, k_2]
    tri_diag[:c, k_2] = FT(0)
end


text:
# Initialize results

function calculateResults(LDDEP::AbstractLDDEProblem, method::DiscretizationMethod, DiscretizationLength::Real; n_steps::Int64=nStepOfLength(DiscretizationLength, method.Δt), calculate_additive::Bool=false)
    result = Result(LDDEP, method, DiscretizationLength,n_steps=n_steps, calculate_additive=calculate_additive)
    calculateResults!(result)
    return (result)
end
function calculateResults!(result::AbstractResult)
    calculateProportionalResult!(result)
    calculateDelayResult!(result)
    if result.calculate_additive
        calculateAdditiveResult!(result)
    end
    # return (result)
end

function calculateProportionalResult!(rst::AbstractResult)
    registerProportionalResult!(rst,
        rst.method(rst.A_avgs, rst))
end
function calculateDelayResult!(rst::AbstractResult)
    registerDelayResult!(rst,
        rst.method.(rst.LDDEP.Bs, Ref(rst)))
end
function calculateAdditiveResult!(rst::AbstractResult)
    registerAdditiveResult!(rst,
        rst.method(rst.LDDEP.c, rst))
end
function registerProportionalResult!(rst::AbstractResult, submxs::Vector{<:SubMX})
    rst.subMXs[1] .= submxs
end

function registerDelayResult!(rst::AbstractResult, submxs::Vector{<:Vector{<:SubMX}})
    for (i,smx) in enumerate(submxs)
        rst.subMXs[i+1] .= smx
    end
end

function registerAdditiveResult!(rst::AbstractResult, subvs::Vector{<:SubV})
    rst.subVs[:] .= subvs;
end
###############################################################################
############################## Discrete Mapping ###############################
###############################################################################
# For every time step
function DiscreteMapping(LDDEP::AbstractLDDEProblem, method::DiscretizationMethod, DiscretizationLength::Real; n_steps::Int64=nStepOfLength(DiscretizationLength, method.Δt), calculate_additive::Bool=false)
    result = calculateResults(LDDEP, method, DiscretizationLength,n_steps = n_steps, calculate_additive = calculate_additive);
    DiscreteMapping(DiscreteMappingSteps(result)...)
end

function DiscreteMappingSteps(rst::AbstractResult{d}) where d
    mappingMXs = [spdiagm(-d => ones(rst.n-d)) for i in 1:rst.n_steps]
    mappingVs = [spzeros(rst.n) for i in rst.subVs]
    for smx_t in rst.subMXs
        addSubmatrixToResult!.(mappingMXs,smx_t)
    end
    addSubvectorToResults!.(mappingVs, rst.subVs)

    (rst.ts, mappingMXs, mappingVs)
end

# If not all state space values are needed
function DiscreteMapping(LDDEP::AbstractLDDEProblem, method::DiscretizationMethod, DiscretizationLength::Real, idxs::AbstractVector{<:Integer}; n_steps::Int64=nStepOfLength(DiscretizationLength, method.Δt), calculate_additive::Bool=false)
    result = calculateResults(LDDEP, method, DiscretizationLength,n_steps = n_steps, calculate_additive = calculate_additive);
    DiscreteMapping(DiscreteMappingSteps(result, idxs)...)
end
function DiscreteMappingSteps(rst::AbstractResult{d},idxs::AbstractVector{<:Integer}) where d
    mappingMXs = [spdiagm(-d => ones(rst.n-d)) for i in 1:rst.n_steps]
    mappingVs = [spzeros(rst.n) for i in rst.subVs]
    for smx_t in rst.subMXs
        addSubmatrixToResult!.(mappingMXs,smx_t)
    end
    addSubvectorToResults!.(mappingVs, rst.subVs)

    (rst.ts, getindex.(mappingMXs,Ref(idxs),Ref(idxs)), getindex.(mappingVs,Ref(idxs)))
end

###############################################################################
############################# Calculate Functions #############################
###############################################################################
function DiscreteMapping_1step(LDDEP::AbstractLDDEProblem, method::DiscretizationMethod, DiscretizationLength::Real; n_steps::Int64=nStepOfLength(DiscretizationLength, method.Δt), calculate_additive::Bool=false)
    mp=DiscreteMapping(LDDEP, method, DiscretizationLength,n_steps = n_steps, calculate_additive = calculate_additive);
    result = initializeOneStepMapping(mp);
    result.ts .= [mp.ts[1],mp.ts[end]];
    result.mappingMXs[1] = prodl(mp.mappingMXs);
    if calculate_additive
        result.mappingVs[1]=reduce_additive(mp.mappingMXs, mp.mappingVs);
    end
    return result
end

# If not all state space values are needed
function DiscreteMapping_1step(LDDEP::AbstractLDDEProblem, method::DiscretizationMethod, DiscretizationLength::Real, idxs::AbstractVector{<:Integer}; n_steps::Int64=nStepOfLength(DiscretizationLength, method.Δt), calculate_additive::Bool=false)
    mp=DiscreteMapping(LDDEP, method, DiscretizationLength,n_steps = n_steps, calculate_additive = calculate_additive);
    result = initializeOneStepMapping(mp);
    result.ts .= [mp.ts[1],mp.ts[end]];
    result.mappingMXs[1] = prodl(mp.mappingMXs,idxs);
    if calculate_additive
        result.mappingVs[1]=reduce_additive(mp.mappingMXs, mp.mappingVs, idxs);
    end
    return result
end

# Some attributes of the mappings
function fixPointOfMapping(dm::DiscreteMapping)
    (I - prodl(dm.mappingMXs)) \ Vector(reduce_additive(dm.mappingMXs, dm.mappingVs))
end
function fixPointOfMapping(dm::DiscreteMapping, idxs::AbstractVector{<:Integer})
    (I - prodl(dm.mappingMXs,idxs)) \ Vector(reduce_additive(dm.mappingMXs, dm.mappingVs, idxs))
end

function spectralRadiusOfMapping(dm::DiscreteMapping; args...)
    abs(eigs(prodl(dm.mappingMXs); args...)[1][1])
end
function spectralRadiusOfMapping(dm::DiscreteMapping, idxs::AbstractVector{<:Integer}; args...)
    abs(eigs(prodl(dm.mappingMXs, idxs); args...)[1][1])
end


text:
@testset "regular float 64" begin
    X = 3 * rand(100, 2)
    Δ = dist(X)
    Y = wda_smacof(Δ, anchors = X)

    X = 3 * rand(300, 2)
    Δ = dist(X)
    # @time Y = smacof(Δ, anchors = X)
    # @show mse(Y, X)
    @time Y = wda_smacof(Δ, anchors = X)
    @show mse(Y, X)

    # for i in 1:2
    #     X = 3 * rand(500, 2)
    #     Δ = dist(X)

    #     @time Y = wda_smacof(Δ, anchors = X)
    #     @show mse(Y, X)
    # end
end


text:
using SimpleDirectMediaLayer
const SDL2 = SimpleDirectMediaLayer
using Test

SDL2_pkg_dir = joinpath(@__DIR__, "..","..")

@testset "Image" begin
# Test that you can init
@test 0 == SDL2.Init(UInt32(SDL2.INIT_VIDEO))

@test SDL2.IMG_Load(joinpath(SDL2_pkg_dir,"assets","cat.bmp")) != Ptr{Nothing}(C_NULL)
@test SDL2.IMG_Load(joinpath(SDL2_pkg_dir,"assets","cat.png")) != Ptr{Nothing}(C_NULL)
@test SDL2.IMG_Load(joinpath(SDL2_pkg_dir,"assets","cat.jpg")) != Ptr{Nothing}(C_NULL)

SDL2.Quit()
end

max_stars_repo_name: Kyjor/julia-game

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: Kyjor/julia-game

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: Kyjor/julia-game

max_forks_repo_licenses: MIT Unlicense


text:
struct PSOOptions{L,M}

    # PSO Options similar to MATLABs PSO implemtation options
    NDims::Int
    Display::String
    DisplayInterval::Int
    FunctionTolerance::Real
    FunValCheck::String
    InertiaRange::SVector{2}
    InitialSwarmMatrix::SizedMatrix{M,L}
    InitialSwarmSpan::SVector{M}
    LowerBounds::SVector{M}
    UpperBounds::SVector{M}
    MaxIterations::Int
    MaxStallIterations::Int
    MaxStallTime::AbstractFloat
    MaxTime::AbstractFloat
    MinNeighborsFraction::AbstractFloat
    ObjectiveLimit::AbstractFloat
    SelfAdjustmentWeight::AbstractFloat
    SocialAdjustmentWeight::AbstractFloat
    SwarmSize::Int
    UseParallel::Bool
    HybridFcn::Bool
    HybridOptimizer::Union{Nothing, Optim.AbstractOptimizer}
    HybridOptimizerOpts::Union{Nothing, Optim.Options}
    HybridAutoDiff::Symbol

    # CreationFcn # Impliment later to allow for different (potentially user defined) swarm initialization functions
    # OutputFcn   # Impliment later to allow user defined output function
    # PlotFcn     # Impliment later to allow user defined plotting function

end

function PSOOptions(NDims; Display = "iter", DisplayInterval = 1, FunctionTolerance = 1.0e-6, FunValCheck = "off",
                    InertiaRange = SVector{2}([0.1, 1.1]), InitialSwarmMatrix = nothing, InitialSwarmSpan = nothing,
                    LowerBounds = nothing, UpperBounds = nothing, MaxIterations = -1, MaxStallIterations = 20,
                    MaxStallTime = Inf64, MaxTime = Inf, MinNeighborsFraction = 0.25, ObjectiveLimit = -Inf,
                    SelfAdjustmentWeight = 1.49, SocialAdjustmentWeight = 1.49, SwarmSize = -1, UseParallel = false,
                    HybridFcn = false, HybridOptimizer = nothing, HybridOptimizerOpts = nothing, HybridAutoDiff = :finite)

    if MaxIterations == -1
        MaxIterations = 200*NDims
    end
    if SwarmSize == -1
        SwarmSize = min(100, 10*NDims)
    end
    if InitialSwarmMatrix === nothing
        InitialSwarmMatrix = SizedMatrix{NDims, 2}(fill(NaN, (NDims, 2)))
    end
    if InitialSwarmSpan === nothing
        InitialSwarmSpan = SVector{NDims}(fill(2000, (NDims, 1)))
    end
    if LowerBounds === nothing
        LowerBounds = SVector{NDims}(fill(-Inf, (NDims, 1)))
    end
    if UpperBounds === nothing
        UpperBounds = SVector{NDims}(fill(Inf, (NDims, 1)))
    end
    if UseParallel == true && Threads.nthreads() < 2
        throw(ArgumentError("More than one thread must be active to use parallel computing."))
    end
    if HybridFcn && HybridOptimizer === nothing
        HybridOptimizer = NelderMead()
    end
    if HybridOptimizer !== nothing
        HybridFcn = true
    end
    if HybridFcn && HybridOptimizerOpts === nothing
        HybridOptimizerOpts = Optim.Options()
    end
    if InertiaRange[1] > InertiaRange[2]
        throw(ArgumentError("InertiaRange[2] must be greater than InertiaRange[1]!"))
    end

    M = NDims
    L = size(InitialSwarmMatrix)[2]
    PSOOptions{L,M}(NDims, Display, DisplayInterval, FunctionTolerance, FunValCheck, InertiaRange, InitialSwarmMatrix,
                    InitialSwarmSpan, LowerBounds, UpperBounds, MaxIterations, MaxStallIterations, MaxStallTime, MaxTime,
                    MinNeighborsFraction, ObjectiveLimit, SelfAdjustmentWeight, SocialAdjustmentWeight, SwarmSize, UseParallel,
                    HybridFcn, HybridOptimizer, HybridOptimizerOpts, HybridAutoDiff)
end

max_stars_repo_name: GrantHecht/JuSwarm

max_forks_repo_name: GrantHecht/JuSwarm


text:
function minkowski_distance(set1::Array, set2::Array, p::Int64)

    @assert length(set1) == length(set2) "Sets don't have the same length."

    sum = zero(eltype(set1))

    for i in 1:length(set1)
        sum += abs(set2[i] - set1[i])^p
    end

    return sum^(1/p)
end