text:
const tex_symbols = Dict{String, TeXExpr}()
const command_to_expr = Dict{String, TeXExpr}()
const symbol_to_expr = Dict{Char, TeXExpr}()

function get_symbol_expr(symbol)
    !haskey(symbol_to_expr, symbol) && return symbol
    return copy(symbol_to_expr[symbol])
end

get_command_expr(command) = copy(command_to_expr["\\" * command])
is_supported_command(command) = haskey(command_to_expr, "\\" * command)

# Register the symbols defined for REPL completion
for (command, symbol) in latex_symbols
    symbol = first(symbol)  # Convert to Char
    expr = TeXExpr(:symbol, [symbol, command])
    tex_symbols[command] = expr
    command_to_expr[command] = expr
    symbol_to_expr[symbol] = expr
end

# Symbols missing from the REPL completion data
tex_symbols[raw"\neq"] = TeXExpr(:symbol, ['≠', raw"\neq"])

function register_command!(expr_builder, command)
    !haskey(tex_symbols, command) && return

    expr = tex_symbols[command]
    symbol = expr.args[1]
    symbol_to_expr[symbol] = expr_builder(expr)
    command_to_expr[command] = expr_builder(expr)
end

register_command!(head::Symbol, command) =
    register_command!(symbol -> TeXExpr(head, [symbol]), command)

function register_function!(expr_builder, function_name)
    func = TeXExpr(:function, [function_name])
    command_to_expr["\\" * function_name] = expr_builder(func)
end

register_function!(head::Symbol, command) =
    register_function!(func -> TeXExpr(head, [func]), command)

function register_symbol!(head::Symbol, symbol)
    symbol_to_expr[symbol] = TeXExpr(head, [symbol])
end

binary_operator_symbols = split("+ * -")
binary_operator_commands = split(raw"""
    \pm             \sqcap                   \rhd
    \mp             \sqcup                   \unlhd
    \times          \vee                     \unrhd
    \div            \wedge                   \oplus
    \ast            \setminus                \ominus
    \star           \wr                      \otimes
    \circ           \diamond                 \oslash
    \bullet         \bigtriangleup           \odot
    \cdot           \bigtriangledown         \bigcirc
    \cap            \triangleleft            \dagger
    \cup            \triangleright           \ddagger
    \uplus          \lhd                     \amalg""")

relation_symbols = split(raw"= < > :")
relation_commands = split(raw"""
    \leq        \geq        \equiv   \models
    \prec       \succ       \sim     \perp
    \preceq     \succeq     \simeq   \mid
    \ll         \gg         \asymp   \parallel
    \subset     \supset     \approx  \bowtie
    \subseteq   \supseteq   \cong    \Join
    \sqsubset   \sqsupset   \neq     \smile
    \sqsubseteq \sqsupseteq \doteq   \frown
    \in         \ni         \propto  \vdash
    \dashv      \dots       \dotplus \doteqdot""")

arrow_commands = split(raw"""
    \leftarrow              \longleftarrow           \uparrow
    \Leftarrow              \Longleftarrow           \Uparrow
    \rightarrow             \longrightarrow          \downarrow
    \Rightarrow             \Longrightarrow          \Downarrow
    \leftrightarrow         \longleftrightarrow      \updownarrow
    \Leftrightarrow         \Longleftrightarrow      \Updownarrow
    \mapsto                 \longmapsto              \nearrow
    \hookleftarrow          \hookrightarrow          \searrow
    \leftharpoonup          \rightharpoonup          \swarrow
    \leftharpoondown        \rightharpoondown        \nwarrow
    \rightleftharpoons      \leadsto""")

spaced_symbols = first.(vcat(binary_operator_symbols, relation_symbols))
register_symbol!.(:spaced, spaced_symbols)

spaced_commands = vcat(binary_operator_commands, relation_commands, arrow_commands)
register_command!.(:spaced, spaced_commands)

underover_commands = split(raw"""
    \sum \prod \coprod \bigcap \bigcup \bigsqcup \bigvee
    \bigwedge \bigodot \bigotimes \bigoplus \biguplus""")

underover_functions = split(raw"lim liminf limsup inf sup min max")

register_command!.(
    symbol -> TeXExpr(:underover, Any[symbol, nothing, nothing]),
    underover_commands)
register_function!.(
    symbol -> TeXExpr(:underover, Any[symbol, nothing, nothing]),
    underover_functions)

integral_commands = split(raw"\int \oint")

register_command!.(
    symbol -> TeXExpr(:integral, Any[symbol, nothing, nothing]),
    integral_commands)

generic_functions = split(raw"""
    arccos csc ker arcsin deg lg Pr arctan det sec arg dim
    sin cos exp sinh cosh gcd ln sup cot hom log tan
    coth tanh""")

register_function!.(identity, generic_functions)

# TODO Add to the parser what come below
punctuation_symbols = split(raw", ; . !")
punctuation_commands = split(raw"\ldotp \cdotp")

space_widths = Dict(
    raw"\,"         => 0.16667,   # 3/18 em = 3 mu
    raw"\thinspace" => 0.16667,   # 3/18 em = 3 mu
    raw"\/"         => 0.16667,   # 3/18 em = 3 mu
    raw"\>"         => 0.22222,   # 4/18 em = 4 mu
    raw"\:"         => 0.22222,   # 4/18 em = 4 mu
    raw"\;"         => 0.27778,   # 5/18 em = 5 mu
    raw"\ "         => 0.33333,   # 6/18 em = 6 mu
    raw"~"          => 0.33333,   # 6/18 em = 6 mu, nonbreakable
    raw"\enspace"   => 0.5,       # 9/18 em = 9 mu
    raw"\quad"      => 1,         # 1 em = 18 mu
    raw"\qquad"     => 2,         # 2 em = 36 mu
    raw"\!"         => -0.16667,  # -3/18 em = -3 mu
)

# Autodelim
ambi_delimiters = split(raw"""
    | \| / \backslash \uparrow \downarrow \updownarrow \Uparrow
    \Downarrow \Updownarrow . \vert \Vert \\|""")

left_delimiter = split(raw"( [ \{ < \lfloor \langle \lceil")
right_delimiter = split(raw") ] \} > \rfloor \rangle \rceil")

## Commands using a braced group as an argument
narrow_accent_map = Dict(
    raw"hat"            => raw"\circumflexaccent",
    raw"breve"          => raw"\combiningbreve",
    raw"bar"            => raw"\combiningoverline",
    raw"grave"          => raw"\combininggraveaccent",
    raw"acute"          => raw"\combiningacuteaccent",
    raw"tilde"          => raw"\combiningtilde",
    raw"dot"            => raw"\combiningdotabove",
    raw"ddot"           => raw"\combiningdiaeresis",
    raw"vec"            => raw"\combiningrightarrowabove",
       "\""             => raw"\combiningdiaeresis",
    raw"`"              => raw"\combininggraveaccent",
    raw"'"              => raw"\combiningacuteaccent",
    raw"~"              => raw"\combiningtilde",
    raw"."              => raw"\combiningdotabove",
    raw"^"              => raw"\circumflexaccent",
    raw"overrightarrow" => raw"\rightarrow",
    raw"overleftarrow"  => raw"\leftarrow",
    raw"mathring"       => raw"\circ",
)

wide_accent_commands = split(raw"\widehat \widetilde \widebar")
fontnames = split(raw"rm cal it tt sf bf default bb frak scr regular")


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

using Markdown
using InteractiveUtils

# ╔═╡ a75d8b3e-0769-11ec-0a8d-0779a2c94f44
using Distributions

# ╔═╡ ae640e10-8015-4846-972c-45d33aff4bd7
using Random

# ╔═╡ 3f22bbfc-12c6-48e1-a2ce-c2d673fdf8d1
rng = Random.seed!(42)

# ╔═╡ eeb51554-b9dd-4721-81db-3d7f862cc890
# Let's start of with the most popular, "Normal Distribution"

n = Normal()

# ╔═╡ 471b0226-a36a-447f-a72a-8f4f833f0c88
# We are creating random set of values which have a simple Normal Distribution (Bell Curve Shape), i.e. mean is zero, standard deviation is 1, skew is zero and a kurtosis of 3. In Normal Distribution, 68% of the values are within mean(+/-)std, 95% are within mean(+/-)2*std, 99.7% are within mean(+/-)3*std.

rd = rand(n, 100)

# ╔═╡ ef470561-eb8c-4c6d-b0c9-985e74f801a2
# It being a probability distribution, we can simply call functions for finding CDF and PDF. Other statistical functions can be applied similarly.

# The cumulative distribution function is the probability that a random variable will take a value less than or equal to the variable.

# The probability density function is the probability that a random variable will take a value exactly equal to itself.

# ╔═╡ da29b537-7def-4d69-aea0-02f23014e379
cdf(n, 1.5) - cdf(n, 0.5)

# ╔═╡ 025f79ce-0d3d-4354-8f16-1d1fede9c551
sum(pdf.(Uniform(0, 1), 0.8)) === pdf(Uniform(0, 1), 0.8)

# ╔═╡ 30585ae9-215e-446a-aa29-12f1c60ed082
# This package contains discrete and continuous univariate, multi-variate and matrix-variate probability distribution classes and their related functions.

# ╔═╡ 9114121f-4df4-49dd-804f-2c99033474cc
# This is a type of discrete mutlivariate.

# If you are unsure about what parameters the various distributions use, you can use `fieldnames` for the same.

fieldnames(Multinomial)

# ╔═╡ dd2223fe-20fa-4e04-80c7-a7bbb86458e5
# Let's try to estimate empirical distribution to a theoretical one.

fit(Normal, rd)

# [μ=0.0, σ=1] for the Normal distribution rd, and theoretically we got values almost equal to the assumption.

# ╔═╡ bc8e946c-2488-4177-b017-a7d455a9fc25
# DataTypes of Distributions package

# ╔═╡ eeafd017-8056-4c74-a29f-2429e3d5750c
#=

The abstract datatype of the package is `Sampleable` which comprises of objects from which we can extract samples i.e. samplers and distributions.

The type of Samplers that can be used are `VariateForm` and `ValueSupport`.

VariateForm as in if it is univariate, multi-variate or matrix-variate.

The subtypes of ValueSupport are discrete  or continuous.

ValueSupport type is used to extract basic information about the samples such as length, size, number of samples, etc. and to extract samples by using random sampler function.

A subtype of Sampleable is `Distribution` which combines both VariateForm types and ValueSupport types which in turn can be used to apply its probability density functions among other things.

=#

# ╔═╡ a6d477ff-5445-403f-9839-e2cad6152ce9
# Univariate Distributions

# ╔═╡ b9d0b2c0-914f-4f67-a839-50aa71958ebc
# NOTE: All distributions are immutable hence inplace vectorized evaluation has been deprecated

# ╔═╡ 1a2d1c6e-5cf9-48f9-b0fd-a3a608e805b8
# Continuous Univariate Distributions

# ╔═╡ e3278e81-5f16-4754-8dba-6bca7f10539d
a = Arcsine()

# ╔═╡ 5fc942b6-18b5-4133-994e-8f6df0c1fd33
params(a) # return values of parameters

# ╔═╡ 26be10a1-68af-4482-9f64-54c970b2e8d1
location(a) # get the left bound of the distribution

# ╔═╡ 9159c57b-ac3d-4bc0-bf89-ec2965953903
b = Beta()

# ╔═╡ a64f12cf-f4c2-468d-b883-5ff2ae51e605
extrema(b) # returns minimum and maximum values as a tuple

# ╔═╡ 9c7c3d1d-74b6-4767-ac3a-f30dc6190680
mean(b) # default mean is returned

# ╔═╡ 5919ad85-4b11-4a0d-95ff-d4aa32b7889f
entropy(b, 0.5) # compute entropy of the distribution w.r.t. the base value given

# ╔═╡ dcbfd625-a4c9-447f-9215-0eed31f23de6
pdfsquaredL2norm(b) # computes squared L2 norm of prob density of the distribution

# ╔═╡ ed33234c-91cb-4f92-9111-6f25fce99faf
b1 = BetaPrime()

# ╔═╡ 4c8ca066-2963-485d-982f-d338eab90c7b
rand(rng, b1) === quantile(b1, rand(rng)) # generate a random sample
# This will output false as rand(rng) is called two times. Generating a random sample is almost equivalent to quantile of the distribution

# ╔═╡ 3ec79418-2b81-430c-8cae-1abb29247c58
(α, β) = params(b1)

# ╔═╡ fd37d4e4-e0d7-4711-8d7e-2ff8507dc83b
rand(rng, Gamma(α)) / rand(rng, Gamma(β)) # this value is equal to rand(rng, b1)

# ╔═╡ e1e7a04a-684a-426b-b3bc-4e7767b2494f
b2 = Biweight()

# ╔═╡ 60fd87ba-08e2-4897-9cd8-187452151978
logpdf(b2, 0.5) # evaluate log of prob density at mass=0.5

# ╔═╡ f1875391-cc8d-49bb-a903-83a138594dd5
quantile(b2, 0.7) # computes inverse cumulative distribution function at 0.7

# ╔═╡ ba40b5f1-3dd6-4fe4-8557-3f87e95822da
c = Cauchy()

# ╔═╡ 324a3a13-87c5-4dc9-8d35-b5f7fbb93109
scale(c)

# ╔═╡ 8daba87d-59e8-44c6-af7d-6e4b0aca5a75
median(c)

# ╔═╡ f74c6cba-99f5-48b1-a6b9-b523ea55930d
c1 = Chernoff()

# ╔═╡ cd6a616f-9054-46f6-ab58-2023ef2732e2
cdf(c1, 0.5) + cdf(c1, -0.5) # evaluate tail probability; almost equivalent to 1

# ╔═╡ a1083141-cf94-408d-8181-bfe2e3daf9b3
var(c1) # evaluate variance of a basic Chernoff distribution

# ╔═╡ f7c41616-b468-4496-8f85-527d8b3fcfda
kurtosis(c1) # evaluate kurtosis of a basic Chernoff distribution

# ╔═╡ 12e16b48-bb71-4944-b915-f00f59633a28
skewness(c1) # evaluate skewness of a basic Chernoff distribution

# ╔═╡ 7e2ea0b4-68c3-46b2-bb76-43d2fcc959db
entropy(c1) # evaluate entropy of a basic Chernoff distribution

# ╔═╡ 5dc859a8-c12f-44da-8ac8-50896fb7c7e2
c2 = Chi(3) # Chi distribution with 3 degrees of freedom

# ╔═╡ 708ac013-715b-4df4-8a40-96501e093274
 dof(c2) # returns degrees of freedom

# ╔═╡ af163e63-604a-4f4c-b57f-b0cc00378324
pdfsquaredL2norm(c2) # returns squared L2 norm of probability density function of Chi distribution

# ╔═╡ f9f6e408-30d4-4f7e-b64f-e3a8ff8e92d1
kurtosis(c2)

# ╔═╡ 98a4752e-53af-4d29-bdea-d76577fd8a5c
sampler(c2) # Samplers can frequently rely on pre-computed values to enhance efficiency. It can be equipped with this sampler technique for batch sampling.

# ╔═╡ cfdfd9a0-9017-4116-be43-fced34b58606
mean(c2)

# ╔═╡ e289e579-30e8-4292-934c-f7d7f56886e5
gradlogpdf(c2, 0.7)

# ╔═╡ c05d0cd1-db06-4818-a05c-6b4d560c7abe
logpdf(c2, 0.5) # logarithm of prob density function at 0.5

# ╔═╡ 76c4ffc3-2d4f-4f2f-b22d-5a9a2812cbdd
 c3 = Chisq(2)

# ╔═╡ f05b9262-f373-491a-b690-b3e1854bd388
params(c3)

# ╔═╡ 9197be83-f28d-462a-800e-b572b942e87e
sampler(c3)

# ╔═╡ 83de8601-722d-4ec2-8bcd-8808984b12fd
pdfsquaredL2norm(c3)

# ╔═╡ 21474254-5742-4c22-b2ad-575f4ed2f47c
logpdf(c3, 0.5)

# ╔═╡ 1d56e469-9c7c-48c3-b997-0dbd38234f83
c4 = Cosine()

# ╔═╡ fd09566a-acdf-4999-9a32-fdfd518177e8
ccdf(c4, 0.3)

# ╔═╡ ac83ee95-5f7b-4df3-9c12-1e66e472b3d6
cdf(c4, 0.7)

# ╔═╡ b0aabc61-10f1-453c-aefa-c098d614ecd6
skewness(c4)

# ╔═╡ f9295f00-f7f6-44f1-ad1d-3bab3a628c5b
mode(c4)

# ╔═╡ 754505f8-b080-4c3d-a6b2-79d69e06e771
e = Erlang()

# ╔═╡ db1654f3-99f4-4aaf-93a6-f54be2506a39
scale(e)

# ╔═╡ 3d068866-48bd-4f9c-937d-71683bbfc7b8
rate(e)

# ╔═╡ e9d696ec-5bd5-41a9-b666-4ab504f028e7
kurtosis(e)

# ╔═╡ e836d5f3-61b3-49d8-9dce-99f455884abc
var(e)

# ╔═╡ 60095213-fd06-4b3e-a131-38cbe6d74f71
mgf(e, 1.8) # evaluates moment generating function

# ╔═╡ 597d59f3-03ed-4333-b99e-dcaf536793b3
cf(e, 8.9) # evaluates characteristic function

# ╔═╡ b4fe0130-9560-4d81-b4e0-1e3dcc3a769e
e1 = Exponential()

# ╔═╡ e6052dea-e63d-4a4f-987e-a3dffe9d3bc9
suffstats(Exponential, rand(1000))

# ╔═╡ 4974f123-952c-4626-8346-9c42697590b1
fit_mle(Exponential, rand(1000))

# ╔═╡ 0f12c5c0-823f-4e14-8256-d524cb258143
f = FDist(4.5, 7)

# ╔═╡ 084c1f11-6fe7-496b-8603-071d27573e92
entropy(f)

# ╔═╡ 4f49bdff-5c44-4c2f-b6db-753a24619d76
rand(f, 50)

# ╔═╡ 997aaa92-0227-47c7-87ea-61c0fa10ccc8
g = Gamma(4.5, 7)

# ╔═╡ 28339f8a-582a-477c-ae7b-3d3e19cd50b5
fit_mle(Gamma, rand(1000))

# ╔═╡ d757ff04-4448-4bd7-8d3c-239a5e55179b
entropy(g)

# ╔═╡ 8f24b9b2-d2b7-4308-9b64-f0c7e274463f
kurtosis(g)

# ╔═╡ 3db1002f-abc3-4781-9862-2bf82a057b81
gradlogpdf(g, 0.6)

# ╔═╡ 60859f4c-a549-4851-85c8-cc1eb1a1ea9f
mgf(g, -.4)

# ╔═╡ 652fcbda-ac03-4a90-958c-a3739b08fab0
i = InverseGamma()

# ╔═╡ b4f733ed-3069-4f81-89e2-a2c069d6b25b
params(i)

# ╔═╡ 09a52cce-ca36-4d1c-8fbc-38196bffc09d
cquantile(i, .75)

# ╔═╡ 28d6e92f-cbb5-40f6-b4f1-81fdec46456d
invlogccdf(i, -.9)

# ╔═╡ f663e150-2ebe-401c-a0bb-8a910476ee02
i1 = InverseGaussian()

# ╔═╡ 5df72ac7-1f4d-4cdd-9dea-b2260b5197e0
fit_mle(InverseGaussian, rand(10000))

# ╔═╡ 0f511fed-05f6-4970-9eba-0bc58182d91f
logccdf(i1, -.75)

# ╔═╡ 444f1b65-7eab-4550-9b08-3c83a60c7f64
p = Pareto()

# ╔═╡ 268a0e38-1ec1-4fdc-a947-c65ba851e8eb
fit_mle(Pareto, rand(10000))

# ╔═╡ 0689dbf0-255e-4c30-a88b-2851914d19f3
entropy(p)

# ╔═╡ 4fa96ad7-ed04-401b-a7ac-e71c24cb73fc
pdf(p, 1.75)

# ╔═╡ a0dbf262-25c0-4d23-9b4b-1a17b20539be
l = LogNormal()

# ╔═╡ 20437f49-3bff-44a4-ba86-b3a30689450c
logccdf(l, 0.4)

# ╔═╡ b4e49482-aa94-43a8-a472-7230ece557c0
cquantile(l, 0.75)

# ╔═╡ 488a1cf3-6a06-49c7-93c9-bcf388ea5943
fit_mle(LogNormal, rand(10000))

# ╔═╡ 76f712db-0a60-4518-a06c-3bee6d2aa08e
l1 = LogitNormal()

# ╔═╡ 67147ca1-2d89-44a8-9c11-353ee6cb293c
invlogccdf(l1, 1.4)

# ╔═╡ 3a15c778-3d84-4c65-929a-22f43bca517c
fit_mle(LogitNormal, rand(1000))

# ╔═╡ 388f32bc-debf-4f54-878e-377a58cd9a66
l2 = Logistic()

# ╔═╡ b359a15e-440a-4849-935f-2ca7af3c263e
gradlogpdf(l2, 3.5)

# ╔═╡ aeb2b5ab-3214-44b2-8eee-039ef4e9684d
skewness(l2)

# ╔═╡ 78a85346-1d4c-4c20-8191-870ca37b04bd
l3 = Laplace()

# ╔═╡ 7f59cdda-709a-45e5-b6a0-53301066a6b6
mgf(l3, 7.7)

# ╔═╡ ef057b36-8339-48ff-8440-441b63e51b15
fit_mle(Laplace, rand(1000))

# ╔═╡ 5f60603b-3d67-4cd2-9f7a-fd4837db6c1d
r = Rayleigh(2)

# ╔═╡ 77411715-da83-4f24-b915-91a2a76d207c
logcdf(r, 0.4)

# ╔═╡ 87a363c1-0eb9-4b5a-b888-0054d933de7c
entropy(r)

# ╔═╡ e20b0e9b-4e81-42a6-9af9-a2dbcf0b052c
fit_mle(Rayleigh, rand(1000))

# ╔═╡ 6cc6aa0c-fbdd-42a2-afd1-7a21929699e6
s = SymTriangularDist()

# ╔═╡ 8d65ec5d-758f-4afd-8408-fef03eb3b8c3
invlogcdf(s, -0.5)

# ╔═╡ 4b926625-2d61-4326-b12c-783803b5ac99
quantile(s, 0/95)

# ╔═╡ d8d9ed48-bd64-4a0f-9dbd-42b9115d8837
t = TDist(5)

# ╔═╡ c5b30329-c34b-45c9-92b5-fd03828d2fc1
dof(t)

# ╔═╡ e92de989-f4c3-4998-9900-61dba2b79a5d
gradlogpdf(t, 1.7)

# ╔═╡ a66a69ae-a726-4cf6-94bb-bf9b88b18efb
entropy(t)

# ╔═╡ 23b2af75-73dd-4be0-ad41-760fcc533f50
t1 = TriangularDist(1, 5, 2)

# ╔═╡ 6fe6eef3-6a85-414e-9947-1a29e235b406
rand(t1)

# ╔═╡ 8ee67312-0a6a-4fa4-b7c3-519b3f146df1
logpdf(t1, 0.34)

# ╔═╡ 37e6381a-30a4-481c-8d1c-55ccb87fb353
mgf(t1, 7.4)

# ╔═╡ 1296320b-d714-4b80-87f8-1cefb1dbf868
t2 = Triweight()

# ╔═╡ a86924e6-39ed-4f58-b60c-4473f6a9bb8d
mgf(t2, 4.8)

# ╔═╡ 76151351-4f16-4d58-b2ae-c91dde660477
cf(t2, 8.98)

# ╔═╡ 8b886e2d-6eb0-43c1-9cbd-822b5e01e8ce
u = Uniform()

# ╔═╡ 1e3a631b-e4a8-412e-a4e1-3beeace38c62
location(u)

# ╔═╡ 6d232648-2dd4-402f-8b54-e79c6608afd8
kurtosis(u)

# ╔═╡ 3b3b6e04-d25c-4aa6-ad68-7466bbb6dbc3
gradlogpdf(u, 1.5) # for all values

# ╔═╡ 49e3c863-eca3-499e-828a-dca0f4039de4
v = VonMises()

# ╔═╡ 14ab4e0f-203f-4fab-9036-3696b7ef3ab8
logpdf(v, 0.9)

# ╔═╡ 866c601b-9e48-4493-ba2e-39df9864141a
cf(v, -1.6)

# ╔═╡ ff6b783b-a2a1-43fb-b38d-b0fd967841ca
sampler(v)

# ╔═╡ 245cc832-6d01-4367-ad8b-4f1b55b5390e
w = Weibull()

# ╔═╡ 820d1231-d71f-45a6-90e0-6d22d15156a5
partype(w) # shows type of parameters

# ╔═╡ e1f9c16b-bcab-4ef5-b115-77a7a9362470
cquantile(w, 0.75) # evaluates 75th percentile of the distribution

# ╔═╡ fef0ce0e-44d6-4484-80e1-45bd5c9c7e38
invlogccdf(w, 0.3)

# ╔═╡ 21578ccb-3943-48fd-b1fd-6237bdf4f80f
x = Array(rand(rng, 1000))

# ╔═╡ 0a3bb76a-3ef7-41f5-80ca-df935932da77
fit(Weibull, collect(1:10)) #log will only return a complex result if called with a complex argument, i.e. values less than 0 and the condition α > zero(α) && θ > zero(θ) is satisfied.

# ╔═╡ c17a6aa4-79a7-426c-bdee-96039b800857
fit_mle(Weibull, x) # computes maximum likelihood estimate

# ╔═╡ b88f392c-8a1a-4acf-b3cd-c31d7a4fd627
#=
NOTE: Go through https://turing.ml/dev/tutorials/00-introduction/ to get an easy understanding of Bayesian Inferential learning in Julia using Turing.jl package. Go through https://github.com/JuliaNLSolvers/Optim.jl to get a hang of optimization of univariate and multi-variate functions. We haven't covered this in the book.
=#

# ╔═╡ b3c270fb-fa77-41d4-8f3c-718f3aac6f77
# Discrete Univariate Distributions

# ╔═╡ 22f01bfb-5ea9-4f7c-975d-e906b8ae3563
Bernoulli(0.5) === Bernoulli()

# ╔═╡ ce3cc1f3-6a54-47b2-a409-f71d8dcb7afb
ber = Bernoulli(0.3)

# ╔═╡ a1a99278-7e8e-4b4e-a795-d00d5e40465e
succprob(ber) # success probability

# ╔═╡ 99786701-3a49-4433-bb05-ab5e5fa47d1e
params(ber)

# ╔═╡ ff15af23-cb50-4441-9592-1819421a6205
skewness(ber)

# ╔═╡ 64b58c58-9ce8-4e42-911d-d4f780d7c1b6
mode(ber)

# ╔═╡ 40be7660-d642-4c48-9cbe-3396106034cc
median(ber)

# ╔═╡ 5109a629-ff2a-4bf2-8e86-78f30d05a41d
logpdf(ber, 1.5)

# ╔═╡ ba7ce037-e348-48bc-8140-238e542b29c4
ccdf(ber, 0.8)

# ╔═╡ 6b058d97-d87c-442f-b938-163095b271a5
ccdf(ber, 7)

# ╔═╡ c9a2150d-6a4f-40b3-aaa3-d659b75c57d8
mgf(ber, 6.8)

# ╔═╡ a15190b3-27ac-48e2-b5bb-d36a802b38d6
bet = BetaBinomial(6, 0.5, 0.4)

# ╔═╡ 11f55be0-b133-4522-8782-0a3e1f3fde0b
params(bet)

# ╔═╡ 4199ca17-5342-46cf-80db-fd512413f2cd
ntrials(bet)

# ╔═╡ 0d6ea1fa-1741-4400-84bd-1e364c260741
kurtosis(bet)

# ╔═╡ 86b136ca-d77f-4fdf-a0c6-0d6f26d559d9
mean(bet)

# ╔═╡ 1a0b6cb1-c981-4e4b-b94f-b96bc7843e71
cdf(bet, 0.6)

# ╔═╡ 6200b397-026f-4484-9680-025bff7e73e1
quantile(bet, 0.6)

# ╔═╡ 5308d425-1862-42c6-a305-4026876576ac
bin = Binomial()

# ╔═╡ 11b450ad-d234-4964-a953-c184d9a12098
ntrials(bin)

# ╔═╡ 1bd7bb23-cae3-4314-896e-72eb0b40c514
params(bin)

# ╔═╡ 47d3d946-3d6a-475d-b9dc-1e89e47de8c0
failprob(bin)

# ╔═╡ c25e89a0-be69-4613-be6b-aabdec176dec
var(bin)

# ╔═╡ 72e46d50-99a9-4c97-90cf-ea4e682a21ae
median(bin)

# ╔═╡ 573fafea-56ed-48e8-88ec-df36f244c8f5
kurtosis(bin)

# ╔═╡ 639e791f-01e6-4481-8df9-0a1fa2482d1b
cat = Categorical(7)

# ╔═╡ 6a67111e-9224-4459-b89b-f374a08d0be6
params(cat)

# ╔═╡ 21072c0c-e840-47a4-ba7a-a783aa7a089c
probs(cat) # probability vector

# ╔═╡ 1f409f70-fee6-414b-8c58-cd71d2a8b5a3
ncategories(cat)

# ╔═╡ f25fccd1-7478-4b92-bf82-b7f3a84dd930
ccdf(cat, rand(1000))

# ╔═╡ 38ee2580-a7a8-4dd4-b3ca-281d97a65edc
ccdf(cat, [2, 5, 3, 6])

# ╔═╡ 250622de-03a0-4675-a91f-9cdb059fcca7
cdf(cat, rand(1000))

# ╔═╡ 554abb12-c1e6-4c90-aeaf-f4c8d97379a3
cdf(cat, [2, 5, 3, 6])

# ╔═╡ 6a742c0c-a6f4-45b8-b5a0-f8c9f842773b
dir = Dirac(7.8)

# ╔═╡ 7b99a92b-699e-42d3-9e0b-60cdefc9ef0f
insupport(dir, rand(100))

# ╔═╡ ab0549c7-21a3-4028-8db2-2718d25ef5a6
entropy(dir)

# ╔═╡ 8d0c5648-4ba5-42ae-ab36-25eb69b88649
var(dir)

# ╔═╡ 3b6865f9-f8ab-4f77-ad81-b654ef3de5ae
logcdf(dir, rand(100))

# ╔═╡ e9b8f83b-beae-476b-a9be-9a6c56632f06
dis = DiscreteUniform(47, 69)

# ╔═╡ cd5cc619-d400-4266-9094-a0a83f58fbf0
probval(dis)

# ╔═╡ 7c6489db-35ec-40de-b0e0-ccf93f0a90ee
span(dis)

# ╔═╡ e425cdba-1cc9-4073-8bd1-e009bf32bb04
params(dis)

# ╔═╡ 6619fb02-1afe-4350-b8cf-2ca595c3e6a0
show(dis)

# ╔═╡ 1125bd37-fb19-45df-860d-f86fbc08d823
mean(dis)

# ╔═╡ b1bdffc4-8f7c-4714-a7be-145ecf4ad3b8
skewness(dis)

# ╔═╡ 0dfe1b52-04cf-4dd9-885d-74c8e0cac799
entropy(dis)

# ╔═╡ 992a775a-d4ab-4d73-a499-459d9fa3e776
pdf(dis)

# ╔═╡ 1b235fca-23ec-46b3-ba93-33508582bef0
cdf(dis, 7)

# ╔═╡ 9fb079ba-8eea-4fca-a170-e648717af5f4
dnp = DiscreteNonParametric([2.3, 4.5, 6.5, 7.8, 9.2], [0.4, 0.3, 0.2, 0.1, 0.0])

# ╔═╡ f4b48aef-a841-4acc-a3a1-c419a135d647
params(dnp)

# ╔═╡ 364db8a0-bd2d-49a5-bffc-de56814e5888
support(dnp)

# ╔═╡ d1918575-46e4-4876-a94c-3e232b48b4f2
probs(dnp)

# ╔═╡ d997cfbb-a5dd-4728-b0a7-9a6f2ddce42f
pdf(dnp, rand(100))

# ╔═╡ 57424fc3-2161-4772-baec-a23e82e749c4
ccdf(dnp, rand(100))

# ╔═╡ efa82d81-4fdd-455d-acc0-9a2147a5de1b
maximum(dnp)

# ╔═╡ 0f6d5f54-ba8f-4ff0-ae11-a4703bc3a5ab
minimum(dnp)

# ╔═╡ 0dfc6f69-a1c7-432e-8d35-65e0b9df8d73
entropy(dnp)

# ╔═╡ a9488b85-218e-45ea-896d-d849f485e958
kurtosis(dnp)

# ╔═╡ da9be907-3b32-4bea-8c8b-94ababcd2dc5
geo = Geometric(0.8)

# ╔═╡ 13187e58-9a5e-4313-84c3-18936f0a0aef
failprob(geo)

# ╔═╡ a5d36503-0a65-4e5a-b6b2-230f459d838d
params(geo)

# ╔═╡ f25d30ae-f2ce-461d-8e2b-cba82f61c031
var(geo)

# ╔═╡ be18640e-1d7b-4c50-b3df-bfc6ebf36089
kurtosis(geo)

# ╔═╡ 7f2bd8ab-fd0d-4764-9d71-a39b495ef1e3
entropy(geo)

# ╔═╡ 55d074c4-2128-4984-bb04-260e617066d4
logccdf(geo, 7)

# ╔═╡ 276e5bb4-ae8f-4666-9abd-59a601fe96a8
invlogcdf(geo, -2.8)

# ╔═╡ 8c1c9824-9bba-4104-92f6-35f0482299f4
cquantile(geo, 2.8)

# ╔═╡ 2e5292ee-c313-41af-8ffe-766d82769ece
hy = Hypergeometric(9, 5, 7)

# ╔═╡ eac6459d-88a1-4c92-98a1-171667fb62bc
params(hy)

# ╔═╡ 388dd1f4-39fd-4d8f-a473-e37abce5e6a1
mean(hy)

# ╔═╡ 0250c6a0-1890-4011-b1e3-fbc43516b0d3
var(hy)

# ╔═╡ 965d15a2-6027-4277-9793-767301770469
mode(hy)

# ╔═╡ 666b607b-19af-44bf-9e30-52f9c0ab248b
kurtosis(hy)

# ╔═╡ 47c63f92-f0ed-494c-9e9a-14e7cc4c6e44
entropy(hy)

# ╔═╡ 6d5db40b-23ad-44db-8442-2a0f8ab6002c
neg = NegativeBinomial()

# ╔═╡ 6555ad08-6d83-4745-93d1-d016dd40e9d7
params(neg)

# ╔═╡ e227dd5b-facf-47f4-a697-bc27bed725a9
var(neg)

# ╔═╡ bc9a8b94-d96d-42d1-986b-af74d9d66c84
skewness(neg)

# ╔═╡ 87c51191-400c-4a05-b8ea-bb9fa5292d9e
logccdf(neg, 0.4)

# ╔═╡ 61842353-23ee-454b-9aae-7c5451d094b4
quantile(neg, 0.59)

# ╔═╡ a1c4bbc4-45af-4cd8-b02c-d59b94ea1cb2
po = Poisson()

# ╔═╡ a447413e-939b-4c47-84d9-00a8f0bdb197
params(po)

# ╔═╡ 74febb89-8545-4759-b4e6-95410b1a7117
rate(po)

# ╔═╡ cd1f29ed-c476-4603-b9d3-fe700cc070b0
skewness(po)

# ╔═╡ dd4fdb1b-85bf-4326-a99b-7fa4a3ac3cbe
pb = PoissonBinomial([0.5, 0.6, 0.8, 0.6, 0.3])

# ╔═╡ 3e3ad40f-4a65-469e-af7e-0410a4918d21
ntrials(pb)

# ╔═╡ 182464ed-24ef-4abf-b98a-8b0655788b06
params(pb)

# ╔═╡ b69e36a6-62b1-4dfd-8594-bbd942e64f96
mean(pb)

# ╔═╡ de4b900b-01b4-44cb-b1e8-16c580de2b74
kurtosis(pb)

# ╔═╡ 7451dbee-0353-467b-ad97-666c1c520819
entropy(pb)

# ╔═╡ 9adfcb62-50bf-4d78-8119-746c1580a92d
modes(pb)

# ╔═╡ c2f1ed3b-9b65-44db-b8ba-1cdeabb14955
quantile(pb, 0.67)

# ╔═╡ bd569a04-5aea-4081-9b0a-45690daba964
cdf(pb, rand(100))

# ╔═╡ cfb63469-12f1-47ae-82a0-660be26051d8
sk = Skellam(7.7)

# ╔═╡ c423d947-1795-4ea3-9c23-9d5558cb581f
params(sk)

# ╔═╡ 7708ce82-755d-4618-9469-f27918656474
skewness(sk)

# ╔═╡ e06e7bb2-d372-404b-897e-6bf447eab874
mean(sk)

# ╔═╡ bb7fc26f-7e86-4c45-ac77-4b16b04e266d
kurtosis(sk)

# ╔═╡ fabc0cb0-d06f-4071-8041-b225ab028984
logpdf(sk, rand(100))

# ╔═╡ e3fd3e8e-37b6-4a79-8a77-a5bece3f6af7
cdf(sk, 7)

# ╔═╡ 99aea4e7-2063-4986-ae99-83a369cf32f6
# Truncated Distributions

# ╔═╡ 2fcc54ae-5ff8-4652-8935-b7cf7baa1e97
tbin = truncated(bin, 8, 9)

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[compat]
Distributions = "~0.25.16"
"""

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

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

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

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

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "30ee06de5ff870b45c78f529a6b093b3323256a3"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "1.3.1"

[[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 = "727e463cfebd0c7b999bbf3e9e7e16f254b94193"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.34.0"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[DataAPI]]
git-tree-sha1 = "bec2532f8adb82005476c141ec23e921fc20971b"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.8.0"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

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

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

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

[[Distributions]]
deps = ["ChainRulesCore", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns"]
git-tree-sha1 = "f4efaa4b5157e0cdb8283ae0b5428bc9208436ed"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.16"

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.5"

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

[[FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "a3b7b041753094f3b17ffa9d2e2e07d8cace09cd"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.12.3"

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

[[IrrationalConstants]]
git-tree-sha1 = "f76424439413893a832026ca355fe273e93bce94"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

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

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

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

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

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

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

[[LogExpFunctions]]
deps = ["ChainRulesCore", "DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "1f5097e3bce576e1cdf6dc9f051ab8c6e196b29e"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.1"

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

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

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "2ca267b08821e86c5ef4376cffed98a46c2cb205"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.1"

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

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

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

[[PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "4dd403333bcf0909341cfe57ec115152f937d7d8"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.1"

[[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"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

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

[[QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "12fbe86da16df6679be7521dfb39fbc861e1dc7b"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.1"

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

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

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

[[Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

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

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

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

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

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

[[SpecialFunctions]]
deps = ["ChainRulesCore", "LogExpFunctions", "OpenSpecFun_jll"]
git-tree-sha1 = "a322a9493e49c5f3a10b50df3aedaf1cdb3244b7"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.6.1"

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

[[StatsAPI]]
git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.0.0"

[[StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "8cbbc098554648c84f79a463c9ff0fd277144b6c"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.10"

[[StatsFuns]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "46d7ccc7104860c38b11966dd1f72ff042f382e4"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.10"

[[SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

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

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

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

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

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

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

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╠═a75d8b3e-0769-11ec-0a8d-0779a2c94f44
# ╠═ae640e10-8015-4846-972c-45d33aff4bd7
# ╠═3f22bbfc-12c6-48e1-a2ce-c2d673fdf8d1
# ╠═eeb51554-b9dd-4721-81db-3d7f862cc890
# ╠═471b0226-a36a-447f-a72a-8f4f833f0c88
# ╠═ef470561-eb8c-4c6d-b0c9-985e74f801a2
# ╠═da29b537-7def-4d69-aea0-02f23014e379
# ╠═025f79ce-0d3d-4354-8f16-1d1fede9c551
# ╠═30585ae9-215e-446a-aa29-12f1c60ed082
# ╠═9114121f-4df4-49dd-804f-2c99033474cc
# ╠═dd2223fe-20fa-4e04-80c7-a7bbb86458e5
# ╠═bc8e946c-2488-4177-b017-a7d455a9fc25
# ╠═eeafd017-8056-4c74-a29f-2429e3d5750c
# ╠═a6d477ff-5445-403f-9839-e2cad6152ce9
# ╠═b9d0b2c0-914f-4f67-a839-50aa71958ebc
# ╠═1a2d1c6e-5cf9-48f9-b0fd-a3a608e805b8
# ╠═e3278e81-5f16-4754-8dba-6bca7f10539d
# ╠═5fc942b6-18b5-4133-994e-8f6df0c1fd33
# ╠═26be10a1-68af-4482-9f64-54c970b2e8d1
# ╠═9159c57b-ac3d-4bc0-bf89-ec2965953903
# ╠═a64f12cf-f4c2-468d-b883-5ff2ae51e605
# ╠═9c7c3d1d-74b6-4767-ac3a-f30dc6190680
# ╠═5919ad85-4b11-4a0d-95ff-d4aa32b7889f
# ╠═dcbfd625-a4c9-447f-9215-0eed31f23de6
# ╠═ed33234c-91cb-4f92-9111-6f25fce99faf
# ╠═4c8ca066-2963-485d-982f-d338eab90c7b
# ╠═3ec79418-2b81-430c-8cae-1abb29247c58
# ╠═fd37d4e4-e0d7-4711-8d7e-2ff8507dc83b
# ╠═e1e7a04a-684a-426b-b3bc-4e7767b2494f
# ╠═60fd87ba-08e2-4897-9cd8-187452151978
# ╠═f1875391-cc8d-49bb-a903-83a138594dd5
# ╠═ba40b5f1-3dd6-4fe4-8557-3f87e95822da
# ╠═324a3a13-87c5-4dc9-8d35-b5f7fbb93109
# ╠═8daba87d-59e8-44c6-af7d-6e4b0aca5a75
# ╠═f74c6cba-99f5-48b1-a6b9-b523ea55930d
# ╠═cd6a616f-9054-46f6-ab58-2023ef2732e2
# ╠═a1083141-cf94-408d-8181-bfe2e3daf9b3
# ╠═f7c41616-b468-4496-8f85-527d8b3fcfda
# ╠═12e16b48-bb71-4944-b915-f00f59633a28
# ╠═7e2ea0b4-68c3-46b2-bb76-43d2fcc959db
# ╠═5dc859a8-c12f-44da-8ac8-50896fb7c7e2
# ╠═708ac013-715b-4df4-8a40-96501e093274
# ╠═af163e63-604a-4f4c-b57f-b0cc00378324
# ╠═f9f6e408-30d4-4f7e-b64f-e3a8ff8e92d1
# ╠═98a4752e-53af-4d29-bdea-d76577fd8a5c
# ╠═cfdfd9a0-9017-4116-be43-fced34b58606
# ╠═e289e579-30e8-4292-934c-f7d7f56886e5
# ╠═c05d0cd1-db06-4818-a05c-6b4d560c7abe
# ╠═76c4ffc3-2d4f-4f2f-b22d-5a9a2812cbdd
# ╠═f05b9262-f373-491a-b690-b3e1854bd388
# ╠═9197be83-f28d-462a-800e-b572b942e87e
# ╠═83de8601-722d-4ec2-8bcd-8808984b12fd
# ╠═21474254-5742-4c22-b2ad-575f4ed2f47c
# ╠═1d56e469-9c7c-48c3-b997-0dbd38234f83
# ╠═fd09566a-acdf-4999-9a32-fdfd518177e8
# ╠═ac83ee95-5f7b-4df3-9c12-1e66e472b3d6
# ╠═b0aabc61-10f1-453c-aefa-c098d614ecd6
# ╠═f9295f00-f7f6-44f1-ad1d-3bab3a628c5b
# ╠═754505f8-b080-4c3d-a6b2-79d69e06e771
# ╠═db1654f3-99f4-4aaf-93a6-f54be2506a39
# ╠═3d068866-48bd-4f9c-937d-71683bbfc7b8
# ╠═e9d696ec-5bd5-41a9-b666-4ab504f028e7
# ╠═e836d5f3-61b3-49d8-9dce-99f455884abc
# ╠═60095213-fd06-4b3e-a131-38cbe6d74f71
# ╠═597d59f3-03ed-4333-b99e-dcaf536793b3
# ╠═b4fe0130-9560-4d81-b4e0-1e3dcc3a769e
# ╠═e6052dea-e63d-4a4f-987e-a3dffe9d3bc9
# ╠═4974f123-952c-4626-8346-9c42697590b1
# ╠═0f12c5c0-823f-4e14-8256-d524cb258143
# ╠═084c1f11-6fe7-496b-8603-071d27573e92
# ╠═4f49bdff-5c44-4c2f-b6db-753a24619d76
# ╠═997aaa92-0227-47c7-87ea-61c0fa10ccc8
# ╠═28339f8a-582a-477c-ae7b-3d3e19cd50b5
# ╠═d757ff04-4448-4bd7-8d3c-239a5e55179b
# ╠═8f24b9b2-d2b7-4308-9b64-f0c7e274463f
# ╠═3db1002f-abc3-4781-9862-2bf82a057b81
# ╠═60859f4c-a549-4851-85c8-cc1eb1a1ea9f
# ╠═652fcbda-ac03-4a90-958c-a3739b08fab0
# ╠═b4f733ed-3069-4f81-89e2-a2c069d6b25b
# ╠═09a52cce-ca36-4d1c-8fbc-38196bffc09d
# ╠═28d6e92f-cbb5-40f6-b4f1-81fdec46456d
# ╠═f663e150-2ebe-401c-a0bb-8a910476ee02
# ╠═5df72ac7-1f4d-4cdd-9dea-b2260b5197e0
# ╠═0f511fed-05f6-4970-9eba-0bc58182d91f
# ╠═444f1b65-7eab-4550-9b08-3c83a60c7f64
# ╠═268a0e38-1ec1-4fdc-a947-c65ba851e8eb
# ╠═0689dbf0-255e-4c30-a88b-2851914d19f3
# ╠═4fa96ad7-ed04-401b-a7ac-e71c24cb73fc
# ╠═a0dbf262-25c0-4d23-9b4b-1a17b20539be
# ╠═20437f49-3bff-44a4-ba86-b3a30689450c
# ╠═b4e49482-aa94-43a8-a472-7230ece557c0
# ╠═488a1cf3-6a06-49c7-93c9-bcf388ea5943
# ╠═76f712db-0a60-4518-a06c-3bee6d2aa08e
# ╠═67147ca1-2d89-44a8-9c11-353ee6cb293c
# ╠═3a15c778-3d84-4c65-929a-22f43bca517c
# ╠═388f32bc-debf-4f54-878e-377a58cd9a66
# ╠═b359a15e-440a-4849-935f-2ca7af3c263e
# ╠═aeb2b5ab-3214-44b2-8eee-039ef4e9684d
# ╠═78a85346-1d4c-4c20-8191-870ca37b04bd
# ╠═7f59cdda-709a-45e5-b6a0-53301066a6b6
# ╠═ef057b36-8339-48ff-8440-441b63e51b15
# ╠═5f60603b-3d67-4cd2-9f7a-fd4837db6c1d
# ╠═77411715-da83-4f24-b915-91a2a76d207c
# ╠═87a363c1-0eb9-4b5a-b888-0054d933de7c
# ╠═e20b0e9b-4e81-42a6-9af9-a2dbcf0b052c
# ╠═6cc6aa0c-fbdd-42a2-afd1-7a21929699e6
# ╠═8d65ec5d-758f-4afd-8408-fef03eb3b8c3
# ╠═4b926625-2d61-4326-b12c-783803b5ac99
# ╠═d8d9ed48-bd64-4a0f-9dbd-42b9115d8837
# ╠═c5b30329-c34b-45c9-92b5-fd03828d2fc1
# ╠═e92de989-f4c3-4998-9900-61dba2b79a5d
# ╠═a66a69ae-a726-4cf6-94bb-bf9b88b18efb
# ╠═23b2af75-73dd-4be0-ad41-760fcc533f50
# ╠═6fe6eef3-6a85-414e-9947-1a29e235b406
# ╠═8ee67312-0a6a-4fa4-b7c3-519b3f146df1
# ╠═37e6381a-30a4-481c-8d1c-55ccb87fb353
# ╠═1296320b-d714-4b80-87f8-1cefb1dbf868
# ╠═a86924e6-39ed-4f58-b60c-4473f6a9bb8d
# ╠═76151351-4f16-4d58-b2ae-c91dde660477
# ╠═8b886e2d-6eb0-43c1-9cbd-822b5e01e8ce
# ╠═1e3a631b-e4a8-412e-a4e1-3beeace38c62
# ╠═6d232648-2dd4-402f-8b54-e79c6608afd8
# ╠═3b3b6e04-d25c-4aa6-ad68-7466bbb6dbc3
# ╠═49e3c863-eca3-499e-828a-dca0f4039de4
# ╠═14ab4e0f-203f-4fab-9036-3696b7ef3ab8
# ╠═866c601b-9e48-4493-ba2e-39df9864141a
# ╠═ff6b783b-a2a1-43fb-b38d-b0fd967841ca
# ╠═245cc832-6d01-4367-ad8b-4f1b55b5390e
# ╠═820d1231-d71f-45a6-90e0-6d22d15156a5
# ╠═e1f9c16b-bcab-4ef5-b115-77a7a9362470
# ╠═fef0ce0e-44d6-4484-80e1-45bd5c9c7e38
# ╠═21578ccb-3943-48fd-b1fd-6237bdf4f80f
# ╠═0a3bb76a-3ef7-41f5-80ca-df935932da77
# ╠═c17a6aa4-79a7-426c-bdee-96039b800857
# ╠═b88f392c-8a1a-4acf-b3cd-c31d7a4fd627
# ╠═b3c270fb-fa77-41d4-8f3c-718f3aac6f77
# ╠═22f01bfb-5ea9-4f7c-975d-e906b8ae3563
# ╠═ce3cc1f3-6a54-47b2-a409-f71d8dcb7afb
# ╠═a1a99278-7e8e-4b4e-a795-d00d5e40465e
# ╠═99786701-3a49-4433-bb05-ab5e5fa47d1e
# ╠═ff15af23-cb50-4441-9592-1819421a6205
# ╠═64b58c58-9ce8-4e42-911d-d4f780d7c1b6
# ╠═40be7660-d642-4c48-9cbe-3396106034cc
# ╠═5109a629-ff2a-4bf2-8e86-78f30d05a41d
# ╠═ba7ce037-e348-48bc-8140-238e542b29c4
# ╠═6b058d97-d87c-442f-b938-163095b271a5
# ╠═c9a2150d-6a4f-40b3-aaa3-d659b75c57d8
# ╠═a15190b3-27ac-48e2-b5bb-d36a802b38d6
# ╠═11f55be0-b133-4522-8782-0a3e1f3fde0b
# ╠═4199ca17-5342-46cf-80db-fd512413f2cd
# ╠═0d6ea1fa-1741-4400-84bd-1e364c260741
# ╠═86b136ca-d77f-4fdf-a0c6-0d6f26d559d9
# ╠═1a0b6cb1-c981-4e4b-b94f-b96bc7843e71
# ╠═6200b397-026f-4484-9680-025bff7e73e1
# ╠═5308d425-1862-42c6-a305-4026876576ac
# ╠═11b450ad-d234-4964-a953-c184d9a12098
# ╠═1bd7bb23-cae3-4314-896e-72eb0b40c514
# ╠═47d3d946-3d6a-475d-b9dc-1e89e47de8c0
# ╠═c25e89a0-be69-4613-be6b-aabdec176dec
# ╠═72e46d50-99a9-4c97-90cf-ea4e682a21ae
# ╠═573fafea-56ed-48e8-88ec-df36f244c8f5
# ╠═639e791f-01e6-4481-8df9-0a1fa2482d1b
# ╠═6a67111e-9224-4459-b89b-f374a08d0be6
# ╠═21072c0c-e840-47a4-ba7a-a783aa7a089c
# ╠═1f409f70-fee6-414b-8c58-cd71d2a8b5a3
# ╠═f25fccd1-7478-4b92-bf82-b7f3a84dd930
# ╠═38ee2580-a7a8-4dd4-b3ca-281d97a65edc
# ╠═250622de-03a0-4675-a91f-9cdb059fcca7
# ╠═554abb12-c1e6-4c90-aeaf-f4c8d97379a3
# ╠═6a742c0c-a6f4-45b8-b5a0-f8c9f842773b
# ╠═7b99a92b-699e-42d3-9e0b-60cdefc9ef0f
# ╠═ab0549c7-21a3-4028-8db2-2718d25ef5a6
# ╠═8d0c5648-4ba5-42ae-ab36-25eb69b88649
# ╠═3b6865f9-f8ab-4f77-ad81-b654ef3de5ae
# ╠═e9b8f83b-beae-476b-a9be-9a6c56632f06
# ╠═cd5cc619-d400-4266-9094-a0a83f58fbf0
# ╠═7c6489db-35ec-40de-b0e0-ccf93f0a90ee
# ╠═e425cdba-1cc9-4073-8bd1-e009bf32bb04
# ╠═6619fb02-1afe-4350-b8cf-2ca595c3e6a0
# ╠═1125bd37-fb19-45df-860d-f86fbc08d823
# ╠═b1bdffc4-8f7c-4714-a7be-145ecf4ad3b8
# ╠═0dfe1b52-04cf-4dd9-885d-74c8e0cac799
# ╠═992a775a-d4ab-4d73-a499-459d9fa3e776
# ╠═1b235fca-23ec-46b3-ba93-33508582bef0
# ╠═9fb079ba-8eea-4fca-a170-e648717af5f4
# ╠═f4b48aef-a841-4acc-a3a1-c419a135d647
# ╠═364db8a0-bd2d-49a5-bffc-de56814e5888
# ╠═d1918575-46e4-4876-a94c-3e232b48b4f2
# ╠═d997cfbb-a5dd-4728-b0a7-9a6f2ddce42f
# ╠═57424fc3-2161-4772-baec-a23e82e749c4
# ╠═efa82d81-4fdd-455d-acc0-9a2147a5de1b
# ╠═0f6d5f54-ba8f-4ff0-ae11-a4703bc3a5ab
# ╠═0dfc6f69-a1c7-432e-8d35-65e0b9df8d73
# ╠═a9488b85-218e-45ea-896d-d849f485e958
# ╠═da9be907-3b32-4bea-8c8b-94ababcd2dc5
# ╠═13187e58-9a5e-4313-84c3-18936f0a0aef
# ╠═a5d36503-0a65-4e5a-b6b2-230f459d838d
# ╠═f25d30ae-f2ce-461d-8e2b-cba82f61c031
# ╠═be18640e-1d7b-4c50-b3df-bfc6ebf36089
# ╠═7f2bd8ab-fd0d-4764-9d71-a39b495ef1e3
# ╠═55d074c4-2128-4984-bb04-260e617066d4
# ╠═276e5bb4-ae8f-4666-9abd-59a601fe96a8
# ╠═8c1c9824-9bba-4104-92f6-35f0482299f4
# ╠═2e5292ee-c313-41af-8ffe-766d82769ece
# ╠═eac6459d-88a1-4c92-98a1-171667fb62bc
# ╠═388dd1f4-39fd-4d8f-a473-e37abce5e6a1
# ╠═0250c6a0-1890-4011-b1e3-fbc43516b0d3
# ╠═965d15a2-6027-4277-9793-767301770469
# ╠═666b607b-19af-44bf-9e30-52f9c0ab248b
# ╠═47c63f92-f0ed-494c-9e9a-14e7cc4c6e44
# ╠═6d5db40b-23ad-44db-8442-2a0f8ab6002c
# ╠═6555ad08-6d83-4745-93d1-d016dd40e9d7
# ╠═e227dd5b-facf-47f4-a697-bc27bed725a9
# ╠═bc9a8b94-d96d-42d1-986b-af74d9d66c84
# ╠═87c51191-400c-4a05-b8ea-bb9fa5292d9e
# ╠═61842353-23ee-454b-9aae-7c5451d094b4
# ╠═a1c4bbc4-45af-4cd8-b02c-d59b94ea1cb2
# ╠═a447413e-939b-4c47-84d9-00a8f0bdb197
# ╠═74febb89-8545-4759-b4e6-95410b1a7117
# ╠═cd1f29ed-c476-4603-b9d3-fe700cc070b0
# ╠═dd4fdb1b-85bf-4326-a99b-7fa4a3ac3cbe
# ╠═3e3ad40f-4a65-469e-af7e-0410a4918d21
# ╠═182464ed-24ef-4abf-b98a-8b0655788b06
# ╠═b69e36a6-62b1-4dfd-8594-bbd942e64f96
# ╠═de4b900b-01b4-44cb-b1e8-16c580de2b74
# ╠═7451dbee-0353-467b-ad97-666c1c520819
# ╠═9adfcb62-50bf-4d78-8119-746c1580a92d
# ╠═c2f1ed3b-9b65-44db-b8ba-1cdeabb14955
# ╠═bd569a04-5aea-4081-9b0a-45690daba964
# ╠═cfb63469-12f1-47ae-82a0-660be26051d8
# ╠═c423d947-1795-4ea3-9c23-9d5558cb581f
# ╠═7708ce82-755d-4618-9469-f27918656474
# ╠═e06e7bb2-d372-404b-897e-6bf447eab874
# ╠═bb7fc26f-7e86-4c45-ac77-4b16b04e266d
# ╠═fabc0cb0-d06f-4071-8041-b225ab028984
# ╠═e3fd3e8e-37b6-4a79-8a77-a5bece3f6af7
# ╠═99aea4e7-2063-4986-ae99-83a369cf32f6
# ╠═2fcc54ae-5ff8-4652-8935-b7cf7baa1e97
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

ext: jl


text:
### Main Program: parse args and call main() ###
using ArgParse

using WWURay

""" parse_cmdline
Parse the command line args to specify scene, camera, and image size """
function parse_cmdline()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--scene", "-s"
            help="scene"
            arg_type=Int
            default=1
        "--camera", "-c"
            help="camera"
            arg_type=Int
            default=1
        "height"
            help="image height"
            arg_type=Int
        "width"
            help="image width"
            arg_type=Int
        "outfile"
            help="output filename"
    end

    return parse_args(s)
end

a = parse_cmdline()

WWURay.main(a["scene"], a["camera"], a["height"], a["width"], a["outfile"])


text:
function plot_chain(bayesian_result::T,filename=nothing) where {T<:MCMCChains.Chains}
	StatsPlots.plot(bayesian_result)
end


text:
function Base.setindex!(dict::MuttsDict{K,V}, v::V, k::K) where {K,V}
    insert!(dict, k, v)
    return v
end

# d[k1, k2, ks...] is syntactic sugar for d[(k1, k2, ks...)].
setindex!(d::MuttsDict, v, k1, k2, ks...) = setindex!(d, v, tuple(k1, k2, ks...))

"""
    insert!(dict, k, v)

Inserts the `(k,v)` pair in the dictionary, replacing any previous entry for `k`.
Returns +1 if a new entry was inserted, or 0 if a previous entry was replaced.
"""
function Base.insert!(dict::MuttsDict{K,V}, k::K, v::V) where {K,V}
    @dassert1 is_mutable(dict)
    _insert!(dict, hash(k), k, v)
end

function _insert!(dict::MuttsDict{K,V}, hashcode::UInt64, k::K, v::V) where {K,V}
    root = dict.root
    if isa(root, Leaf{K,V})
        # Branch specialized for isa(root,Leaf{K,V})
        (!is_mutable(root)) && error("insert!($(typeof(dict)), ..): dict is immutable")
        (new_root,size_change) = _insert!(root, hashcode, k, v)
        dict.root = new_root
    else
        # Branch specialized for isa(root,INode{K,V})
        (!is_mutable(root)) && error("insert!($(typeof(dict)), ..): dict is immutable")
        config = dict.config
        height = length(config.inode_capacity)
        pathmask = UInt64(0)
        pathhash = UInt64(0)
        size_change = _insert!(root, config, height, hashcode, k, v,
            pathmask, pathhash)
    end
    @dassert1(size_change <= 1)
    dict.n += size_change

    # Check if we need to change configurations
    if dict.n == dict.config.next_config_size
        new_config = _config_for_size(dict.n)
        @dassert1 new_config.next_config_size > dict.n
        if length(new_config.inode_capacity) > length(dict.config.inode_capacity)
            # Increase depth
            @dassert1 length(new_config.inode_capacity) == 1+length(dict.config.inode_capacity)
            mark_immutable!(dict.root)
            dict.root = INode{K,V}(new_config.inode_capacity[1], dict.root)
        elseif isa(dict.root, INode{K,V})
            # Possible resizing of root inode
            height = length(config.inode_capacity)
            dict.root = _get_mutable_version(dict.root, new_config, height)
        end
        dict.config = new_config
    end
    return size_change
end

function _insert!(inode::INode{K,V},
    config::NodeConfig,
    height::Int,
    hashcode::UInt64,
    k::K,
    v::V,
    pathmask=UInt64(0),
    pathhash=UInt64(0)
) where {K,V}
    @dassert1 inode.is_mutable
    mask = reinterpret(UInt64, length(inode.entries)-1)
    shift = reinterpret(UInt64, height+1)*UInt64(8)
    slot = 1 + ((hashcode >> shift) & mask)
    pathmask = pathmask | (mask << shift)
    pathhash = pathhash | ((slot-1) << shift)
    entries = inode.entries
    if isa(entries, Vector{Leaf{K,V}})
        @dassert1 height==1
        leaf = @inbounds entries[slot]
        leaf = _get_mutable_version(leaf, pathmask, pathhash)
        (leaf,size_change) = _insert!(leaf, hashcode, k, v)
        @inbounds entries[slot] = leaf
        return size_change
    else
        child_inode = @inbounds entries[slot]
        child_inode = _get_mutable_version(child_inode, config, height-1)
        size_change = _insert!(child_inode, config, height-1, hashcode, k, v, pathmask, pathhash)
        @inbounds inode.entries[slot] = child_inode
        return size_change
    end
end

# Returns (new_leaf,size_change)
# Note: In Julia 1.4 returning a `(mutable_object,Int)` tuple still allocates,
# so we're allocating on every insert. This is fixed in newer Julia versions.
function _insert!(leaf::Leaf{K,V}, hashcode::UInt64, k::K, v::V) where {K,V}
    @dassert1 is_mutable(leaf)
    num_slots = length(leaf.entries)
    max_probe2 = min(max_probe-1, num_slots-1)
    for i=0:max_probe2
        index = _hashindex(hashcode, i, num_slots)
        if _slot_occupied(leaf, index, num_slots)
            # If the key is here replace the entry, otherwise keep probing
            (ki,vi) = @inbounds leaf.entries[index]
            if (k === ki) || isequal(k, ki)
                _put_entry(leaf, index, k, v)
                return (leaf,0)
            end
        else
            # Found an empty slot - place the entry here
            _put_entry(leaf, index, k, v)
            return (leaf,1)
        end
    end

    # Reached maximum probe length, resize.  Rather then modifying this leaf in place,
    # we create a new leaf with greater capacity and copy the entries there. This is
    # simpler, and handles the case where because of rotten luck we exceed max_probe
    # inserting to the new leaf, and have to enlarge the leaf a second time (or more).
    new_capacity = _next_leaf_table_size(num_slots)
    new_leaf = Leaf{K,V}(new_capacity)

    for i=1:length(leaf.entries)
        slot = reinterpret(UInt64,i)
        if _slot_occupied(leaf, slot, num_slots)
            (ki,vi) = @inbounds leaf.entries[i]
            h = hash(ki)
            (new_leaf,_) = _insert!(new_leaf, h, ki, vi)
        end
    end
    return _insert!(new_leaf, hashcode, k, v)
end


text:
#Generate 2D data from 4 spherical (σ=0.25) gaussian clusters centered at [-1,1], [1,-1], [-1,1] and [1,1]
"""
    generate_4_clusters(npts_per_cluster; σ=0.25)
    generate_4_clusters(::DataType, npts_per_cluster; σ=0.25)

Generate 2D data from 4 spherical gaussian clusters, with standard deviation `σ`, centered at [-1,1], [1,-1], [-1,1] and [1,1]. `DataType` can be `Float32` or `Float64`, default is `Float64`.
"""
function generate_4_clusters(T::DataType, npts_per_cluster; σ = 0.25)
    c = convert(Matrix{T},[1 1 -1 -1; 1 -1 1 -1])
    data = Matrix{T}(undef,2,4*npts_per_cluster)
    for i ∈ axes(c,2)
        rpts = randn(T,(2,npts_per_cluster))
        sindx = (i-1)*npts_per_cluster+1
        eindx = i*npts_per_cluster
        data[:,sindx:eindx] = σ*rpts .+ c[:,i]
    end
    return data
end

generate_4_clusters(npts_per_cluster; σ = 0.25) = generate_4_clusters(Float64, npts_per_cluster; σ)

#Generate 2D data arranged in a circle (default radius, r =0.5) around centres at [-1,1], [1,-1], [-1,1] and [1,1]
"""
    generate_4_circle_clusters(npts_per_cluster = 8; r = 0.5)
    generate_4_circle_clusters(::DataType, npts_per_cluster = 8; r = 0.5)

Generate 2D data arranged in a circle (default radius, r =0.5) around centres at [-1,1], [1,-1], [-1,1] and [1,1]. `DataType` can be `Float32` or `Float64`, default is `Float64`.
"""
function generate_4_circle_clusters(T::DataType, npts_per_cluster::Int = 8; r = 0.5)
    c = convert(Matrix{T},[1 1 -1 -1; 1 -1 1 -1])
    data = Matrix{T}(undef,2,4*npts_per_cluster)
    for i ∈ axes(c,2)
        sindx = (i-1)*npts_per_cluster+1
        eindx = i*npts_per_cluster
        pts = generate_circle_cluster(T, npts_per_cluster)
        data[:,sindx:eindx] = r*pts .+ c[:,i]
    end
    return data
end

generate_4_circle_clusters(npts_per_cluster::Int = 8; r = 0.5) = generate_4_circle_clusters(Float64, npts_per_cluster; r)

#Generate 2D data arranged in a circle of radius one around [0,0]
"""
    generate_circle_cluster(npts = 8)
    generate_circle_cluster(::DataType, npts = 8)

Generate 2D data arranged in a circle of radius one around [0,0]. `DataType` can be `Float32` or `Float64`, default is `Float64`.
"""
function generate_circle_cluster(T::DataType, npts::Int = 8)
    Δθ = 2π/npts
    θ = 0:Δθ:2π-Δθ
    x = Matrix{T}(undef,2,npts)
    for i ∈ eachindex(θ)
        x[1,i] = cos(θ[i])
        x[2,i] = sin(θ[i])
    end
    return x
end

generate_circle_cluster(npts::Int = 8) = generate_circle_cluster(Float64, npts)


text:
# `backward_density` functions return gates storing derivatives

# The logic for gradient computation here is analogous to `gradients.jl`
# when interpreting the coefficient vector of DensityMatrix as statevector.

backward_density(g::XGate, ::AbstractMatrix) = g
backward_density(g::YGate, ::AbstractMatrix) = g
backward_density(g::ZGate, ::AbstractMatrix) = g

backward_density(g::HadamardGate, ::AbstractMatrix) = g

backward_density(g::SGate, ::AbstractMatrix) = g
backward_density(g::TGate, ::AbstractMatrix) = g

backward_density(g::SdagGate, ::AbstractMatrix) = g
backward_density(g::TdagGate, ::AbstractMatrix) = g

function backward_density(g::RxGate, Δ::AbstractMatrix)
    c = cos(g.θ[])
    s = sin(g.θ[])
    RxGate(-s*(Δ[3, 3] + Δ[4, 4]) + c*(Δ[4, 3] - Δ[3, 4]))
end

function backward_density_mixed_add(g::RxGate, Δ::AbstractMatrix)
    c = cos(0.5*g.θ[])
    s = sin(0.5*g.θ[])
    RxGate(0.5*(-s*tr(Δ) + c*(Δ[4, 3] - Δ[3, 4])))
end

function backward_density_mixed_sub(g::RxGate, Δ::AbstractMatrix)
    c = cos(0.5*g.θ[])
    RxGate(0.5*c*(Δ[1, 2] + Δ[2, 1]))
end

function backward_density(g::RyGate, Δ::AbstractMatrix)
    c = cos(g.θ[])
    s = sin(g.θ[])
    RyGate(-s*(Δ[2, 2] + Δ[4, 4]) + c*(Δ[2, 4] - Δ[4, 2]))
end

function backward_density_mixed_add(g::RyGate, Δ::AbstractMatrix)
    c = cos(0.5*g.θ[])
    s = sin(0.5*g.θ[])
    RyGate(0.5*(-s*tr(Δ) + c*(Δ[2, 4] - Δ[4, 2])))
end

function backward_density_mixed_sub(g::RyGate, Δ::AbstractMatrix)
    c = cos(0.5*g.θ[])
    RyGate(0.5*c*(Δ[1, 3] + Δ[3, 1]))
end

function backward_density(g::RzGate, Δ::AbstractMatrix)
    c = cos(g.θ[])
    s = sin(g.θ[])
    RzGate(-s*(Δ[2, 2] + Δ[3, 3]) + c*(Δ[3, 2] - Δ[2, 3]))
end

function backward_density_mixed_add(g::RzGate, Δ::AbstractMatrix)
    c = cos(0.5*g.θ[])
    s = sin(0.5*g.θ[])
    RzGate(0.5*(-s*tr(Δ) + c*(Δ[3, 2] - Δ[2, 3])))
end

function backward_density_mixed_sub(g::RzGate, Δ::AbstractMatrix)
    c = cos(0.5*g.θ[])
    RzGate(0.5*c*(Δ[1, 4] + Δ[4, 1]))
end

function backward_density(g::RotationGate, Δ::AbstractMatrix)
    θ = norm(g.nθ)
    if θ == 0
        RotationGate([Δ[4, 3] - Δ[3, 4], Δ[2, 4] - Δ[4, 2], Δ[3, 2] - Δ[2, 3]])
    else
        n = g.nθ/θ
        c = cos(θ)
        s = sin(θ)
        # cross-product matrix
        K = [0 -n[3] n[2]; n[3] 0 -n[1]; -n[2] n[1] 0]
        # -K^2
        K2 = I - reshape(kron(n, n), 3, 3)
        dRθ = c*K - s*K2
        dn = K2 / θ
        dRn1 = [2*(1-c)*n[1] (1-c)*n[2] (1-c)*n[3]; (1-c)*n[2] 0 -s; (1-c)*n[3] s 0]
        dRn2 = [0 (1-c)*n[1] s; (1-c)*n[1] 2*(1-c)*n[2] (1-c)*n[3]; -s (1-c)*n[3] 0]
        dRn3 = [0 -s (1-c)*n[1]; s 0 (1-c)*n[2]; (1-c)*n[1] (1-c)*n[2] 2*(1-c)*n[3]]
        RotationGate([sum((n[i]*dRθ + (dRn1*dn[1, i] + dRn2*dn[2, i] + dRn3*dn[3, i])) .* Δ[2:4, 2:4]) for i in 1:3])
    end
end

function backward_density_mixed_add(g::RotationGate, Δ::AbstractMatrix)
    dsn = [Δ[4, 3] - Δ[3, 4], Δ[2, 4] - Δ[4, 2], Δ[3, 2] - Δ[2, 3]]
    θ = norm(g.nθ)
    if θ == 0
        return RotationGate(0.5 * dsn)
    else
        n = g.nθ/θ
        c = cos(0.5*θ)
        s = sin(0.5*θ)
        dnθ = -0.5*s * tr(Δ) * n + (0.5*c - s/θ) * dot(n, dsn) * n + (s/θ) * dsn
        return RotationGate(dnθ)
    end
end

function backward_density_mixed_sub(g::RotationGate, Δ::AbstractMatrix)
    dsn = [Δ[1, 2] + Δ[2, 1], Δ[1, 3] + Δ[3, 1], Δ[1, 4] + Δ[4, 1]]
    θ = norm(g.nθ)
    if θ == 0
        return RotationGate(0.5 * dsn)
    else
        n = g.nθ/θ
        c = cos(0.5*θ)
        s = sin(0.5*θ)
        dnθ = (0.5*c - s/θ) * dot(n, dsn) * n + (s/θ) * dsn
        return RotationGate(dnθ)
    end
end

function backward_density(g::PhaseShiftGate, Δ::AbstractMatrix)
    # agrees with backward pass for Rz(θ) since global prefactor cancels
    c = cos(g.ϕ[])
    s = sin(g.ϕ[])
    PhaseShiftGate(sum([-s -c; c -s] .* Δ[2:3, 2:3]))
end

function backward_density_mixed_add(g::PhaseShiftGate, Δ::AbstractMatrix)
    c = cos(g.ϕ[])
    s = sin(g.ϕ[])
    PhaseShiftGate(0.5*(s*((Δ[1, 4] + Δ[4, 1]) - tr(Δ)) + c*(Δ[3, 2] - Δ[2, 3])))
end

function backward_density_mixed_sub(g::PhaseShiftGate, Δ::AbstractMatrix)
    c = cos(g.ϕ[])
    s = sin(g.ϕ[])
    PhaseShiftGate(0.5*(c*((Δ[1, 4] + Δ[4, 1]) - tr(Δ)) + s*(Δ[2, 3] - Δ[3, 2])))
end

backward_density(g::SwapGate, Δ::AbstractMatrix) = g

function backward_density(g::EntanglementXXGate, Δ::AbstractMatrix)
    c = cos(g.θ[])
    s = sin(g.θ[])
    EntanglementXXGate(-s*(Δ[3, 3] + Δ[8, 8] + Δ[4, 4] + Δ[7, 7] + Δ[9,  9] + Δ[14, 14] + Δ[10, 10] + Δ[13, 13])
                       +c*(Δ[8, 3] - Δ[3, 8] + Δ[4, 7] - Δ[7, 4] + Δ[14, 9] - Δ[ 9, 14] + Δ[13, 10] - Δ[10, 13]))
end

function backward_density_mixed_add(g::EntanglementXXGate, Δ::AbstractMatrix)
    c = cos(0.5*g.θ[])
    s = sin(0.5*g.θ[])
    EntanglementXXGate(0.5*(-s*tr(Δ) + c*(Δ[4, 7] - Δ[7, 4] + Δ[8, 3] - Δ[3, 8] + Δ[14, 9] - Δ[9, 14] + Δ[13, 10] - Δ[10, 13])))
end

function backward_density_mixed_sub(g::EntanglementXXGate, Δ::AbstractMatrix)
    c = cos(0.5*g.θ[])
    s = sin(0.5*g.θ[])
    EntanglementXXGate(0.5*c*(Δ[6, 1] + Δ[1, 6] + Δ[5, 2] + Δ[2, 5] + Δ[15, 12] + Δ[12, 15] - Δ[16, 11] - Δ[11, 16]))
end

function backward_density(g::EntanglementYYGate, Δ::AbstractMatrix)
    c = cos(g.θ[])
    s = sin(g.θ[])
    EntanglementYYGate(-s*(Δ[2,  2] + Δ[12, 12] + Δ[10, 10] + Δ[4,  4] + Δ[5,  5] + Δ[15, 15] + Δ[7,  7] + Δ[13, 13])
                       +c*(Δ[2, 12] - Δ[12,  2] + Δ[10,  4] - Δ[4, 10] + Δ[5, 15] - Δ[15,  5] + Δ[7, 13] - Δ[13,  7]))
end

function backward_density_mixed_add(g::EntanglementYYGate, Δ::AbstractMatrix)
    c = cos(0.5*g.θ[])
    s = sin(0.5*g.θ[])
    EntanglementYYGate(0.5*(-s*tr(Δ) + c*(Δ[10, 4] - Δ[4, 10] + Δ[2, 12] - Δ[12, 2] + Δ[5, 15] - Δ[15, 5] + Δ[7, 13] - Δ[13, 7])))
end

function backward_density_mixed_sub(g::EntanglementYYGate, Δ::AbstractMatrix)
    c = cos(0.5*g.θ[])
    s = sin(0.5*g.θ[])
    EntanglementYYGate(0.5*c*(Δ[11, 1] + Δ[1, 11] + Δ[9, 3] + Δ[3, 9] + Δ[14, 8] + Δ[8, 14] - Δ[16, 6] - Δ[6, 16]))
end

function backward_density(g::EntanglementZZGate, Δ::AbstractMatrix)
    c = cos(g.θ[])
    s = sin(g.θ[])
    EntanglementZZGate(-s*(Δ[3,  3] + Δ[14, 14] + Δ[15, 15] + Δ[2,  2] + Δ[12, 12] + Δ[5,  5] + Δ[9, 9] + Δ[8, 8])
                       +c*(Δ[3, 14] - Δ[14,  3] + Δ[15,  2] - Δ[2, 15] + Δ[12,  5] - Δ[5, 12] + Δ[9, 8] - Δ[8, 9]))
end

function backward_density_mixed_add(g::EntanglementZZGate, Δ::AbstractMatrix)
    c = cos(0.5*g.θ[])
    s = sin(0.5*g.θ[])
    EntanglementZZGate(0.5*(-s*tr(Δ) + c*(Δ[15, 2] - Δ[2, 15] + Δ[3, 14] - Δ[14, 3] + Δ[12, 5] - Δ[5, 12] + Δ[9, 8] - Δ[8, 9])))
end

function backward_density_mixed_sub(g::EntanglementZZGate, Δ::AbstractMatrix)
    c = cos(0.5*g.θ[])
    s = sin(0.5*g.θ[])
    EntanglementZZGate(0.5*c*(Δ[16, 1] + Δ[1, 16] + Δ[13, 4] + Δ[4, 13] + Δ[10, 7] + Δ[7, 10] - Δ[6, 11] - Δ[11, 6]))
end

# shortcuts for non-parametric gates

backward_density(g::ControlledGate{XGate}, ::AbstractMatrix) = g
backward_density(g::ControlledGate{YGate}, ::AbstractMatrix) = g
backward_density(g::ControlledGate{ZGate}, ::AbstractMatrix) = g

backward_density(g::ControlledGate{HadamardGate}, ::AbstractMatrix) = g

backward_density(g::ControlledGate{SGate}, ::AbstractMatrix) = g
backward_density(g::ControlledGate{TGate}, ::AbstractMatrix) = g
backward_density(g::ControlledGate{SdagGate}, ::AbstractMatrix) = g
backward_density(g::ControlledGate{TdagGate}, ::AbstractMatrix) = g

backward_density(g::ControlledGate{SwapGate}, ::AbstractMatrix) = g

function backward_density(g::ControlledGate{G}, Δ::AbstractMatrix) where {G}

    # number of target wires
    T = target_wires(g)

    # number of control wires
    C = control_wires(g)

    # for a single control qubit:
    # controlled-U = |0><0| x I + |1><1| x U = I + |1><1| x (U - I)
    # (generalizes to several control qubits)

    # mixed term |1><1| x (U - I) ρ + ρ |1><1| x (U† - I)
    Δadd = zeros(eltype(Δ), 4^T, 4^T)
    Δsub = zeros(eltype(Δ), 4^T, 4^T)
    eo = binary_digits(C, 0)
    for p in 0:2^C-1
        binary_digits!(eo, p)
        Δr = reshape(Δ, 4^T, 4^C, 4^T, 4^C)
        for j in 1:C
            # use another variable for type stability
            Δt = reshape(Δr, 4^T, 4, 4^(C-j), 4^T, 4, 4^(C-j))
            if eo[j] == 0
                Δr = 0.5 * (Base.view(Δt, :, 1, :, :, 1, :) + Base.view(Δt, :, 2, :, :, 2, :) + Base.view(Δt, :, 3, :, :, 3, :) + Base.view(Δt, :, 4, :, :, 4, :) - Base.view(Δt, :, 1, :, :, 4, :) - Base.view(Δt, :, 4, :, :, 1, :))
            else
                Δr = 0.5 * (Base.view(Δt, :, 3, :, :, 2, :) - Base.view(Δt, :, 2, :, :, 3, :))
            end
        end
        nodd::Int = sum(eo)
        # sign factor from even powers of i
        signfac = 1 - 2*(((nodd + 1) ÷ 2) % 2)
        if iseven(nodd)
            Δadd += (2*signfac) * reshape(Δr, 4^T, 4^T)
        else
            Δsub += (2*signfac) * reshape(Δr, 4^T, 4^T)
        end
    end
    # backward_density_mixed_add(g.U, Δadd) delayed until later
    dU = backward_density_mixed_sub(g.U, Δsub)

    # conjugation by |1><1| x (U - I)
    Δr = reshape(Δ, 4^T, 4^C, 4^T, 4^C)
    for j in 1:C
        # use another variable for type stability
        Δt = reshape(Δr, 4^T, 4, 4^(C-j), 4^T, 4, 4^(C-j))
        Δr = 0.5 * (Base.view(Δt, :, 1, :, :, 1, :) + Base.view(Δt, :, 4, :, :, 4, :) - Base.view(Δt, :, 1, :, :, 4, :) - Base.view(Δt, :, 4, :, :, 1, :))
    end
    dU += backward_density(g.U, reshape(Δr, 4^T, 4^T))
    Δadd -= 2 * reshape(Δr, 4^T, 4^T)
    dU += backward_density_mixed_add(g.U, Δadd)

    return ControlledGate{G}(dU, g.M)
end

backward_density(g::MatrixGate, ::AbstractMatrix) = g

# shortcuts for non-parametric gates

backward_density(cg::CircuitGate{M,XGate}, ::DensityMatrix, ::DensityMatrix) where {M} = cg
backward_density(cg::CircuitGate{M,YGate}, ::DensityMatrix, ::DensityMatrix) where {M} = cg
backward_density(cg::CircuitGate{M,ZGate}, ::DensityMatrix, ::DensityMatrix) where {M} = cg

backward_density(cg::CircuitGate{M,HadamardGate}, ::DensityMatrix, ::DensityMatrix) where {M} = cg

backward_density(cg::CircuitGate{M,SGate}, ::DensityMatrix, ::DensityMatrix) where {M} = cg
backward_density(cg::CircuitGate{M,TGate}, ::DensityMatrix, ::DensityMatrix) where {M} = cg
backward_density(cg::CircuitGate{M,SdagGate}, ::DensityMatrix, ::DensityMatrix) where {M} = cg
backward_density(cg::CircuitGate{M,TdagGate}, ::DensityMatrix, ::DensityMatrix) where {M} = cg

backward_density(cg::CircuitGate{M,SwapGate}, ::DensityMatrix, ::DensityMatrix) where {M} = cg

function backward_density(cg::CircuitGate{M,G}, ρ::DensityMatrix, Δ::DensityMatrix) where {M,G}
    @assert ρ.N == Δ.N
    CircuitGate{M,G}(cg.iwire, backward_density(cg.gate, rdm(ρ.N, cg.iwire, Δ.v, ρ.v, 4)))
end

function backward_density(m::Moment, ρ::DensityMatrix, Δ::DensityMatrix)
    dgates = CircuitGate[]
    for cg in reverse(m.gates)
        Udag = Base.adjoint(cg)
        ρ = apply(ρ, Udag)
        # backward step of quantum state
        pushfirst!(dgates, backward_density(cg, ρ, Δ))
        Δ = apply(Δ, Udag)
    end
    return Moment(dgates), ρ, Δ
end

function backward_density(moments::Vector{Moment}, ρ::DensityMatrix, Δ::DensityMatrix)
    dmoments = Moment[]
    for moment in reverse(moments)
        # backward step of quantum state
        (m, ρ, Δ) = backward_density(moment, ρ, Δ)
        pushfirst!(dmoments, m)
    end
    return dmoments, Δ
end

"""
    gradients(c::Circuit{N}, ρ::DensityMatrix, Δ::AbstractVector{<:Real}) where {N}

Perform a backward pass to compute gradients of a (fictitious) cost function with
respect to the circuit parameters of `c` and input density matrix `ρ`. `Δ` contains
the cost function gradients with respect to the circuit outputs (measurement operator averages).
The gradients with respect to the circuit parameters are returned in a duplicate circuit;
the overall return value is the tuple (dc::Circuit{N}, dρ::DensityMatrix).
"""
function gradients(c::Circuit{N}, ρ::DensityMatrix, Δ::AbstractVector{<:Real}) where {N}
    @assert(ρ.N == N)
    # length of circuit output vector must match gradient vector
    @assert(length(Δ) == length(c.meas))
    # forward pass through unitary gates
    ρ = apply(ρ, c.moments)
    # gradient of cost function with respect to ρ
    ρbar = density_from_matrix(0.5^N * sum([Δ[i] * sparse_matrix(c.meas[i], N) for i in 1:length(Δ)]))
    # backward pass through unitary gates
    dcgc, ρbar = backward_density(c.moments, ρ, ρbar)
    dmeas = MeasurementOperator.([Δ[i] * matrix(ρ) for i in 1:length(Δ)], (Tuple(1:N),))
    return Circuit{N}(dcgc, dmeas), ρbar
end


text:
"""
	struct OptArgs

Storing parameters for `SplitGraph.runopt` function.

### General
- `γ::Real = 2`
- `itmax::Int64 = 15`: maximal number of iterations of naive MCCs / SA cycles.
- `likelihood_sort::Bool = true`: sort equivalent configurations using likelihood test
  based on branch length.
- `resolve::Bool = true`: try to resolve trees while finding MCCs.
- `seq_lengths`: lengths of sequences that trees were built from.
  Used in likelihood calculations.
  This is initialized from other input arguments, and defaults to sequences of length one.
### Simulated annealing
- `Md::Real = 10`:  number of SA iterations (per temperature) for a tree of `n` leaves is
  `ceil(Int, n/Md)`.
- `cooling_schedule = :geometric`: type of cooling schedule `(:geometric, :linear)`
- `Tmin::Float64 = 1e-3`: minimal temperature of SA.
- `Tmax::Float64 = 1`: maximal temperature of SA.
- `αT::Float64 = 0.95`: ratio between terms in the geometric cooling.
- `dT::Float64 = 1e-2`: temperature step in the linear cooling.
### Verbosity
- `verbose::Bool=false`: first level of verbosity
- `vv::Bool = false`: second level of verbosity
"""
@with_kw struct OptArgs
	γ::Real = 2
	itmax::Int64 = 15
	likelihood_sort::Bool = true
	resolve::Bool = true
	seq_lengths::Vector{Int} = [1, 1]
	# For the annealing
	Md::Real = 5
	Tmin::Float64 = 1e-3
	Tmax::Float64 = 1.; @assert Tmax > Tmin
	αT::Float64 = 0.95
	dT::Float64 = 1e-2
	cooling_schedule = :geometric
	Trange = get_cooling_schedule(; Tmin, Tmax, dT, αT, type=cooling_schedule)
	sa_rep::Int64 = 1
	# Verbosity
	verbose::Bool = false
	vv::Bool = false
end

function get_cooling_schedule(;
	Tmax=1., Tmin = 1e-3, dT = 1e-2, αT = 0.95, type=:geometric,
)
	if type == :geometric
		return get_geometric_cooling_schedule(Tmin, Tmax, αT)
	elseif type == :linear
		return get_linear_cooling_schedule(Tmin, Tmax, dT)
	else
		error("Unknown `cooling_schedule` field: $(type). See `?OptArgs` for allowed values.")
	end
end

function get_geometric_cooling_schedule(Tmin, Tmax, αT)
	n = ceil(Int, (log(Tmin) - log(Tmax)) / log(αT))
	return [αT^i * Tmax for i in 0:n]
end

get_linear_cooling_schedule(Tmin, Tmax, dT) = return reverse(Tmin:dT:Tmax)

max_stars_repo_name: PierreBarrat/RecombTools

max_issues_repo_name: PierreBarrat/RecombTools

max_forks_repo_name: PierreBarrat/RecombTools


text:
function max_violation_dynamics(prob::Problem{T,Discrete})::T where T <: AbstractFloat
    n = prob.model.n; m = prob.model.m
    max_viol = norm(prob.X[1] - prob.x0,Inf)
    X̄ = zeros(n)
    for k = 1:prob.N-1
        X̄ .= 0
        evaluate!(X̄, prob.model, prob.X[k], prob.U[k][1:m], get_dt(prob,prob.U[k]))
        max_viol = max(max_viol,norm(X̄ - prob.X[k+1],Inf))
    end
    return max_viol
end

# assume dircol solve for continuous problems and defaults to implicit rk3 integration
function max_violation_dynamics(prob::Problem{T,Continuous})::T where T <:AbstractFloat
    n,m,N = size(prob)
    X̄ = zeros(prob.model.n)

    fVal = [zeros(prob.model.n) for k = 1:N]
    fValm = [zeros(prob.model.n) for k = 1:N]
    Xm = [zeros(prob.model.n) for k = 1:N]
    Um = [zeros(prob.model.m) for k = 1:N-1]
    dt = get_dt_traj(prob,prob.U)

    # Calculate midpoints
    for k = 1:N
        evaluate!(fVal[k], prob.model, prob.X[k], prob.U[k])
    end
    for k = 1:N-1
        Xm[k] = (prob.X[k] + prob.X[k+1])/2 + dt[k]/8*(fVal[k] - fVal[k+1])
        Um[k] = (prob.U[k][1:m] + prob.U[k+1][1:m])*0.5
        evaluate!(fValm[k], prob.model, Xm[k], Um[k])
    end

    max_viol = norm(prob.X[1] - prob.x0,Inf)

    for k = 1:N-1
        mv = norm(-prob.X[k+1] + prob.X[k] + dt[k]*(fVal[k] + 4*fValm[k] + fVal[k+1])/6,Inf)
        max_viol = max(max_viol,mv)
    end

    return max_viol
end

function max_violation_direct(prob::Problem)
    max(max_violation(prob),max_violation_dynamics(prob))
end

""" $(SIGNATURES)
Get the row and column lists of a sparse matrix, with ordered elements
"""
function get_rc(A::SparseMatrixCSC)
    row,col,inds = findnz(A)
    v = sortperm(inds)
    row[v],col[v]
end

function convertInf!(A::VecOrMat{Float64},infbnd=1.1e20)
    infs = isinf.(A)
    A[infs] = sign.(A[infs])*infbnd
    return nothing
end

"""Number of collocation constraints"""
num_colloc(prob::Problem)::Int = (prob.N-1)*prob.model.n

#######################
#   COST FUNCTIONS    #
#######################

"Generate state midpoint according to quadrature rule"
function gen_xm_cubic(prob::Problem)
    ẋ = zeros(prob.model.n); ẏ = zeros(prob.model.n)

    function xm(y,x,v,u,h)
        prob.model.f(ẋ,x,u)
        prob.model.f(ẏ,y,v)

        0.5*(y+x) + h/8*(ẋ - ẏ)
    end
end

function gen_stage_cost(prob::Problem)
    xm = gen_xm_cubic(prob)
    obj = prob.obj

    function cost(X,U,H)
        N = length(X)
        J = 0.0
        for k = 1:N-1
            Xm = xm(X[k+1],X[k],U[k+1],U[k],H[k])
            Um = 0.5*(U[k] + U[k+1])
            J += H[k]/6*(stage_cost(obj[k],X[k],U[k]) + 4*stage_cost(obj[k],Xm,Um) + stage_cost(obj[k],X[k+1],U[k+1]))
        end
        J += stage_cost(obj[N],X[N])
        return J
    end
end

function gen_stage_cost_gradient(prob::Problem)
    n = prob.model.n; m = prob.model.m; N = prob.N
    function fc(z)
        ż = zeros(eltype(z),n)
        prob.model.f(ż,z[1:n],z[n .+ (1:m)])
        return ż
    end

    function fc(x,u)
        ẋ = zero(x)
        prob.model.f(ẋ,x,u)
        return ẋ
    end

    ∇fc(z) = ForwardDiff.jacobian(fc,z)
    ∇fc(x,u) = ∇fc([x;u])
    dfcdx(x,u) = ∇fc(x,u)[:,1:n]
    dfcdu(x,u) = ∇fc(x,u)[:,n .+ (1:m)]

    xm(y,x,v,u,h) = 0.5*(y + x) + h/8*(fc(x,u) - fc(y,v))
    dxmdy(y,x,v,u,h) = 0.5*I - h/8*dfcdx(y,v)
    dxmdx(y,x,v,u,h) = 0.5*I + h/8*dfcdx(x,u)
    dxmdv(y,x,v,u,h) = -h/8*dfcdu(y,v)
    dxmdu(y,x,v,u,h) = h/8*dfcdu(x,u)

    dℓdx(obj,x,u) = obj.Q*x + obj.q + obj.H'*u
    dℓdu(obj,x,u) = obj.R*u + obj.r + obj.H*x

    dgdx(obj,y,x,v,u,h) = h/6*(dℓdx(obj,x,u) + 4*dxmdx(y,x,v,u,h)'*dℓdx(obj,xm(y,x,v,u,h),0.5*(u+v)))
    dgdy(obj,y,x,v,u,h) = h/6*(4.0*dxmdy(y,x,v,u,h)'*dℓdx(obj,xm(y,x,v,u,h),0.5*(u+v))+ dℓdx(obj,y,v))
    dgdu(obj,y,x,v,u,h) = h/6*(dℓdu(obj,x,u) + 4*(dxmdu(y,x,v,u,h)'*dℓdx(obj,xm(y,x,v,u,h),0.5*(u+v)) + 0.5*dℓdu(obj,xm(y,x,v,u,h),0.5*(u+v))))
    dgdv(obj,y,x,v,u,h) = h/6*(4*(dxmdv(y,x,v,u,h)'*dℓdx(obj,xm(y,x,v,u,h),0.5*(u+v)) + 0.5*dℓdu(obj,xm(y,x,v,u,h),0.5*(u+v))) + dℓdu(obj,y,v))

    nn = 2*(n+m)
    _tmp_ = zeros(n)

    function _cost_grad!(∇g,X,U,H)
        shift = 0
        ∇g .= 0 # set all to zero, for additive inplace operations
        for k = 1:N-1
            obj = prob.obj[k]
            x = X[k]; y = X[k+1]; u = U[k]; v = U[k+1]; h = H[k]

            _xm = xm(y,x,v,u,h)
            _um = 0.5*(u+v)
            _dℓdx = dℓdx(obj,x,u)
            _dℓdu = dℓdu(obj,x,u)
            _dℓdy = dℓdx(obj,y,v)
            _dℓdv = dℓdu(obj,y,v)
            _dℓmdx = dℓdx(obj,_xm,_um)
            _dℓmdu = dℓdu(obj,_xm,_um)
            _dxmdx = dxmdx(y,x,v,u,h)
            _dxmdy = dxmdy(y,x,v,u,h)
            _dxmdu = dxmdu(y,x,v,u,h)
            _dxmdv = dxmdv(y,x,v,u,h)

            _∇g = view(∇g,shift .+ (1:nn))
            _∇g[1:n] += h/6*(_dℓdx + 4*_dxmdx'*_dℓmdx)
            _∇g[n .+ (1:m)] += h/6*(_dℓdu + 4*(_dxmdu'*_dℓmdx + 0.5*_dℓmdu))
            _∇g[(n+m) .+ (1:n)] += h/6*(4.0*_dxmdy'*_dℓmdx + dℓdx(obj,y,v))
            _∇g[(2*n+m) .+ (1:m)] += h/6*(4*(_dxmdv'*_dℓmdx + 0.5*_dℓmdu) + dℓdu(obj,y,v))
            shift += (n+m)
        end

        gradient!(_tmp_, prob.obj[N], X[N])
        ∇g[(N-1)*(n+m) .+ (1:n)] += _tmp_
    end
end

cost(prob::Problem, solver::DIRCOLSolver) = cost(prob, solver.Z)
cost(prob::Problem, Z::Primals) = cost(prob.obj, Z.X, Z.U, get_dt_traj(prob))

##############################
#   COST FUNCTION GRADIENT   #
##############################

function cost_gradient!(grad_f, prob::Problem, X::AbstractVectorTrajectory, U::AbstractVectorTrajectory, H::Vector)
    n,m,N = size(prob)
    grad = reshape(grad_f, n+m, N)
    part = (x=1:n, u=n .+ (1:m))
    dt = get_dt_traj(prob,U)
    for k = 1:N-1
        grad_k = PartedVector(view(grad,:,k), part)
        gradient!(grad_k, prob.obj[k], X[k], U[k])
        grad_k .*= dt[k]
    end
    grad_k = PartedVector(view(grad,1:n,N), part)
    gradient!(grad_k, prob.obj[N], X[N])
    return nothing
end

############################
#    DYNAMICS FUNCTIONS    #
############################

# function traj_points!(prob::Problem, solver::DIRCOLSolver{T,HermiteSimpson}, X, U) where T
#     n,m,N = size(prob)
#     dt = prob.dt
#     Xm = solver.X_
#     fVal = solver.fVal
#     X,U = Z.X, Z.U
#     for k = 1:N-1
#         Xm[k] = (X[k] + X[k+1])/2 + dt/8*(fVal[k] - fVal[k+1])
#     end
#     return Xm
# end
#
# function TrajectoryOptimization.dynamics!(prob::Problem{T,Continuous}, solver::DirectSolver, X, U) where T<:AbstractFloat
#     for k = 1:prob.N
#         evaluate!(solver.fVal[k], prob.model, X[k], U[k])
#     end
# end

#####################
#    CONSTRAINTS    #
#####################

function update_constraints!(g, prob::Problem, solver::DIRCOLSolver, X, U) where T
    n,m,N = size(prob)
    p = solver.p
    pcum = [0; cumsum(p)]
    for k = 1:N
        if p[k] > 0
            k == N ? part = :terminal : part = :stage
            part_c = solver.C[k].parts
            inds = pcum[k] .+ (1:p[k])
            if k == N
                evaluate!(PartedArray(view(g, inds), part_c), prob.constraints[k], X[k])
            else
                evaluate!(PartedArray(view(g, inds), part_c), prob.constraints[k], X[k], U[k])
            end
        end
    end
end

# function constraint_jacobians(prob::Problem, solver::DirectSolver, X::AbstractVectorTrajectory, U::AbstractVectorTrajectory)
#     n,m,N = size(prob)
#     for k = 1:N
#         if k == prob.N
#             jacobian!(solver.∇C[k], prob.constraints[k], X[k])
#         else
#             jacobian!(solver.∇C[k], prob.constraints[k], X[k], U[k])
#         end
#     end
# end

function partition_constraint_jacobian(jac::AbstractMatrix, prob::Problem)
    n,m,N = size(prob)
    p_colloc = (N-1)*n
    jac_colloc = view(jac, 1:p_colloc, :)
    jac_custom = view(jac, p_colloc+1:size(jac,1), :)
    return jac_colloc, jac_custom
end

function constraint_jacobian!(jac, prob::Problem, solver::DirectSolver, X, U)
    n,m,N = size(prob)
    p_colloc = num_colloc(prob)
    p = solver.p
    off = 0
    for k = 1:N
        part_c = solver.∇C[k].parts
        if k == N
            n_blk = p[k]*n
            block = PartedArray(reshape(view(jac, off .+ (1:n_blk)), p[k], n), part_c)
            jacobian!(block, prob.constraints[k], X[k])
        else
            n_blk = p[k]*(n+m)
            block = PartedArray(reshape(view(jac, off .+ (1:n_blk)), p[k], n+m), part_c)
            jacobian!(block, prob.constraints[k], X[k], U[k])
        end
        off += n_blk
    end
end

function constraint_jacobian_sparsity!(jac::AbstractMatrix, prob::Problem)
    n,m,N = size(prob)
    p_colloc = num_colloc(prob)
    nG_colloc = p_colloc*2(n+m)
    jac_colloc, jac_custom = partition_constraint_jacobian(jac, prob)
    collocation_constraint_jacobian_sparsity!(jac_colloc, prob)

    p = num_constraints(prob)
    off1 = 0
    off2 = 0
    off = nG_colloc
    for k = 1:N
        k == N ? b2 = n : b2 = n+m
        n_blk = p[k]*b2
        block = view(jac_custom, off1 .+ (1:p[k]), (off2 .+ (1:b2)))
        block .= reshape(off .+ (1:n_blk), p[k], b2)
        off1 += p[k]
        off2 += b2
        off += n_blk
    end
end

#################################
#    COLLOCATION CONSTRAINTS    #
#################################

function collocation_constraints!(g, prob::Problem, solver::DIRCOLSolver{T,HermiteSimpson}, X, U) where T
    n,m,N = size(prob)
    dt = prob.dt
    fVal = solver.fVal  #[zero(X[1]) for k = 1:N]
    Xm = solver.X_  #[zero(X[1]) for k = 1:N-1]
    g_colloc = reshape(g, n, N-1)

    # Calculate midpoints
    for k = 1:N
        evaluate!(fVal[k], prob.model, X[k], U[k])
    end
    for k = 1:N-1
        Xm[k] = (X[k] + X[k+1])/2 + dt/8*(fVal[k] - fVal[k+1])
    end
    fValm = copy(fVal[1])
    for k = 1:N-1
        Um = (U[k] + U[k+1])*0.5
        evaluate!(fValm, prob.model, Xm[k], Um)
        g_colloc[:,k] = -X[k+1] + X[k] + dt*(fVal[k] + 4*fValm + fVal[k+1])/6
    end
end

# Calculate jacobian
function calc_block!(vals::PartedMatrix, F1,F2,Fm,dt)
    n,m = size(F1.xu)
    In = Diagonal(I, n)
    Im = Diagonal(I, m)
    vals.x1 .= dt/6*(F1.xx + 4Fm.xx*( dt/8*F1.xx + In/2)) + In
    vals.u1 .= dt/6*(F1.xu + 4Fm.xx*( dt/8*F1.xu) + 4Fm.xu*(Im/2))
    vals.x2 .= dt/6*(F2.xx + 4Fm.xx*(-dt/8*F2.xx + In/2)) - In
    vals.u2 .= dt/6*(F2.xu + 4Fm.xx*(-dt/8*F2.xu) + 4Fm.xu*(Im/2))
    return nothing
end

function collocation_constraint_jacobian!(jac, prob::Problem, solver::DIRCOLSolver{T,HermiteSimpson}, X, U) where T
    n,m,N = size(prob)
    dt = prob.dt

    # Compute dynamics jacobians
    F = solver.∇F
    for k = 1:N
        jacobian!(F[k], prob.model, X[k], U[k])
    end

    # Calculate midpoints
    fVal = solver.fVal  # [zeros(n) for k = 1:N]
    Xm = solver.X_  # [zeros(n) for k = 1:N-1]
    for k = 1:N
        evaluate!(fVal[k], prob.model, X[k], U[k])
    end
    for k = 1:N-1
        Xm[k] = (X[k] + X[k+1])/2 + dt/8*(fVal[k] - fVal[k+1])
    end

    # Collocation jacobians
    Fm = PartedMatrix(prob.model)
    n_blk = 2(n+m)n
    off = 0
    In = Matrix(I,n,n)
    Im = Matrix(I,m,m)
    part = create_partition2((n,),(n,m,n,m), Val((:x1,:u1,:x2,:u2)))
    for k = 1:N-1
        block = PartedArray(reshape(view(jac, off .+ (1:n_blk)), n, 2(n+m)), part)
        Um = (U[k] + U[k+1])/2
        jacobian!(Fm, prob.model, Xm[k], Um)
        calc_block!(block, F[k], F[k+1], Fm, dt)
        off += n_blk
    end
end

function collocation_constraint_jacobian_sparsity!(jac::AbstractMatrix, prob::Problem)
    n,m,N = size(prob)
    n_blk = 2(n+m)n

    blk = 1:n_blk
    off1 = 0
    off2 = 0
    for k = 1:N-1
        block = view(jac, off1 .+ (1:n), off2 .+ (1:2(n+m)))
        block .= reshape((k-1)*n_blk .+ blk, n, 2(n+m))
        off1 += n
        off2 += n+m
    end
end

# generate DIRCOL functions
function gen_dircol_functions(prob::Problem{T}, solver::DIRCOLSolver) where T
    n,m,N = size(prob)
    NN = N*(n+m)
    p_colloc = num_colloc(prob)
    p = num_constraints(prob)
    P = p_colloc + sum(p)
    dt = prob.dt

    part_f = create_partition2(prob.model)
    part_z = create_partition(n,m,N,N)
    pcum = cumsum(p)

    jac_structure = spzeros(P, NN)
    constraint_jacobian_sparsity!(jac_structure, prob)
    r,c = get_rc(jac_structure)

    function eval_f(Z)
        X,U = unpack(Z,part_z)
        cost(prob.obj, X, U, get_dt_traj(prob))
    end

    function eval_grad_f(Z, grad_f)
        X,U = unpack(Z, part_z)
        cost_gradient!(grad_f, prob, X, U, get_dt_traj(prob))
    end

    function eval_g(Z, g)
        X,U = unpack(Z,part_z)
        g_colloc = view(g,1:p_colloc)
        g_custom = view(g,p_colloc+1:length(g))

        collocation_constraints!(g_colloc, prob, solver, X, U)
        update_constraints!(g_custom, prob, solver, X, U)
    end

    function eval_jac_g(Z, mode, rows, cols, vals)
        if mode == :Structure
            copyto!(rows,r)
            copyto!(cols,c)
        else
            X,U = unpack(Z, part_z)

            nG_colloc = p_colloc * 2(n+m)
            jac_colloc = view(vals, 1:nG_colloc)
            collocation_constraint_jacobian!(jac_colloc, prob, solver, X, U)

            # General constraint jacobians
            jac_custom = view(vals, nG_colloc+1:length(vals))
            constraint_jacobian!(jac_custom, prob, solver, X, U)
        end

        return nothing
    end
    return eval_f, eval_g, eval_grad_f, eval_jac_g
end

function remove_bounds!(prob::Problem)
    n,m,N = size(prob)
    bounds = [BoundConstraint(n,m) for k = 1:prob.N]

    # Initial Time step
    if :bound ∈ labels(prob.constraints[1])
        bnd_init = remove_bounds!(prob.constraints[1])[1]
    else
        bnd_init = bounds[1]
    end
    bounds[1] = BoundConstraint(n,m, x_min=prob.x0, u_min=bnd_init.u_min,
                                     x_max=prob.x0, u_max=bnd_init.u_max)

    # All time steps
    for k = 2:prob.N
        bnd = remove_bounds!(prob.constraints[k])
        if !isempty(bnd)
            bounds[k] = bnd[1]::BoundConstraint
        end
    end
if :goal ∈ labels(prob.constraints[N])
        goal = pop!(prob.constraints[N])
        xf = zeros(n)
        evaluate!(xf, goal, zero(xf))
        term_bound = BoundConstraint(n,m, x_min=-xf, u_min=bounds[N-1].u_min,
                                          x_max=-xf, u_max=bounds[N-1].u_max)
        bounds[N] = term_bound::BoundConstraint
    end
    # Terminal time step
    #TODO handle control at Nth U differently

    if :goal ∈ labels(prob.constraints[N])
        goal = pop!(prob.constraints[N])
        xf = zeros(n)
        evaluate!(xf, goal, zero(xf))
        term_bound = BoundConstraint(n,m, x_min=-xf, u_min=bounds[N-1].u_min,
                                          x_max=-xf, u_max=bounds[N-1].u_max)
        bounds[N] = term_bound::BoundConstraint
    end
    return bounds
end

function get_bounds(prob::Problem, bounds::Vector{<:BoundConstraint})
    n,m,N = size(prob)
    p_colloc = num_colloc(prob)
    Z = Primals(prob, true)

    Z.equal ? uN = N : uN = N-1
    x_U = [zeros(n) for k = 1:N]
    x_L = [zeros(n) for k = 1:N]
    u_U = [zeros(m) for k = 1:uN]
    u_L = [zeros(m) for k = 1:uN]
    for k = 1:uN
        x_U[k] = bounds[k].x_max
        x_L[k] = bounds[k].x_min
        u_U[k] = bounds[k].u_max
        u_L[k] = bounds[k].u_min
    end
    #TODO handle control at Nth U differently
    if Z.equal
        u_U[N] = bounds[N-1].u_max
        u_L[N] = bounds[N-1].u_min
    end
    if !Z.equal
        x_U = bounds[N].x_max
        x_L = bounds[N].x_min
    end
    z_U = Primals(x_U,u_U)
    z_L = Primals(x_L,u_L)

    # Constraints
    p = num_constraints(prob)
    g_U = [PartedVector(prob.constraints[k]) for k = 1:N-1]
    g_L = [PartedVector(prob.constraints[k]) for k = 1:N-1]
    push!(g_U, PartedVector(prob.constraints[N], :terminal))
    push!(g_L, PartedVector(prob.constraints[N], :terminal))
    for k = 1:N
        if p[k] > 0
            g_L[k].inequality .= -Inf
        end
    end
    g_U = vcat(zeros(p_colloc), g_U...)
    g_L = vcat(zeros(p_colloc), g_L...)

    convertInf!(z_U.Z)
    convertInf!(z_L.Z)
    convertInf!(g_U)
    convertInf!(g_L)
    return z_U.Z, z_L.Z, g_U, g_L
end

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
export Portfolio

function Portfolio(cbmo::Float64,mu::Vector,
                   Q::Any,u::Vector{Float64},
                   kappa::Int64,z0::Vector{Float64}
                   ;relax::Function=NL,
                   t::Float64=0.0,
		   n::Int64=5,
		   tol::Float64=1e-6)

	c1=z->sum(z[1:n]);lc1=1.0;uc1=1.0;
	c2=z->sum(z[n+1:2*n]);lc2=n-kappa;uc2=Inf;
	c3=z->z[1:n].*z[n+1:2*n];lc3=-Inf*ones(n);uc3=zeros(n);
	#c3=z->sum(z[1:n].*z[n+1:2*n]);lc3=-Inf;uc3=zeros(1);
	c3=z->relax(z[1:n],z[n+1:2*n],t=t);lc3=-Inf*ones(n);uc3=zeros(n);
	cons=z->vcat(c1(z),c2(z),c3(z));
	fobj=z->cbmo*sqrt(z[1:n]'*Q*z[1:n])-z[1:n]'*mu;
	lvar=zeros(2*n);uvar=vcat(u,ones(n));
	mp=ADNLPModel(fobj, z0, lvar=lvar, uvar=uvar, c=cons, lcon=vcat(lc1,lc2,lc3), ucon=vcat(uc1,uc2,uc3))

 return mp
end

function Portfolio(cbmo::Float64,mu::Vector,
                   Q::Any,u::Vector{Float64},
                   kappa::Int64,z0::Vector{Float64},
		   tol::Float64=1e-6)
	n=size(Q,1)

	c1=z->sum(z[1:n]);lc1=1.0;uc1=1.0;
	c2=z->sum(z[n+1:2*n]);lc2=n-kappa;uc2=Inf;
	c3=z->z[1:n].*z[n+1:2*n];lc3=-Inf*ones(n);uc3=zeros(n);
	cons=z->vcat(c1(z),c2(z),c3(z));
	fobj=z->cbmo*sqrt(z[1:n]'*Q*z[1:n])-z[1:n]'*mu;
	lvar=zeros(2*n);uvar=vcat(u,ones(n));
	mp=ADNLPModel(fobj, z0, lvar=lvar, uvar=uvar, c=cons, lcon=vcat(lc1,lc2,lc3), ucon=vcat(uc1,uc2,uc3))

	#ex1=JuMP.Model()
	#JuMP.@variable(ex1,x[1:2],start=1.0)
	#JuMP.@NLobjective(ex1,Min,x[1]-x[2])
	#JuMP.@constraint(ex1,1-x[2]>=0)
	#ex1=MathProgNLPModel(ex1)

	#cG=z->z[1:n]
	#G=ADNLPModel(x->(), z0, c=cG, lcon=zeros(n))
	#cH=z->z[n+1:2*n]
	#H=ADNLPModel(x->(), z0, c=cH, lcon=zeros(n))

	G=JuMP.Model()
	JuMP.@variable(G,x[1:2*n],start=1.0)
	JuMP.@constraint(G,x[1:n].>=0)
	JuMP.@NLobjective(G,Min,0.0)
	G=MathProgNLPModel(G)
	H=JuMP.Model()
	JuMP.@variable(H,x[1:2*n],start=1.0)
	JuMP.@constraint(H,x[n+1:2*n].>=0)
	JuMP.@NLobjective(H,Min,0.0)
	H=MathProgNLPModel(H)

 return mp,G,H
end


text:
#####
##### Horizontal Laplacians
#####

@inline function ∇²hᶜᶜᵃ(i, j, k, grid, c)
    return 1/Vᵃᵃᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_∂xᶠᵃᵃ, c) +
                                    δyᵃᶜᵃ(i, j, k, grid, Ay_∂yᵃᶠᵃ, c))
end

@inline function ∇²hᶠᶜᵃ(i, j, k, grid, u)
    return 1/Vᵃᵃᶜ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, Ax_∂xᶜᵃᵃ, u) +
                                    δyᵃᶜᵃ(i, j, k, grid, Ay_∂yᵃᶠᵃ, u))
end

@inline function ∇²hᶜᶠᵃ(i, j, k, grid, u)
    return 1/Vᵃᵃᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_∂xᶠᵃᵃ, u) +
                                    δyᵃᶠᵃ(i, j, k, grid, Ay_∂yᵃᶜᵃ, u))
end

@inline ∇²hᶜᶜᵃ(i, j, k, grid, F::FU, args...) where FU <: Function = ∂²xᶜᵃᵃ(i, j, k, grid, F, args...) + ∂²yᵃᶜᵃ(i, j, k, grid, F, args...)
@inline ∇²hᶠᶜᵃ(i, j, k, grid, F::FU, args...) where FU <: Function = ∂²xᶠᵃᵃ(i, j, k, grid, F, args...) + ∂²yᵃᶜᵃ(i, j, k, grid, F, args...)
@inline ∇²hᶜᶠᵃ(i, j, k, grid, F::FU, args...) where FU <: Function = ∂²xᶜᵃᵃ(i, j, k, grid, F, args...) + ∂²yᵃᶠᵃ(i, j, k, grid, F, args...)

#####
##### 3D Laplacian
#####

"""
    ∇²(i, j, k, grid, c)

Calculates the Laplacian of c via

    1/V * [δxᶜᵃᵃ(Ax * ∂xᶠᵃᵃ(c)) + δyᵃᶜᵃ(Ay * ∂yᵃᶠᵃ(c)) + δzᵃᵃᶜ(Az * ∂zᵃᵃᶠ(c))]

which will end up at the location `ccc`.
"""
@inline function ∇²(i, j, k, grid, c)
    return 1/Vᵃᵃᶠ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, Ax_∂xᶜᵃᵃ, c) +
                                    δyᵃᶠᵃ(i, j, k, grid, Ay_∂yᵃᶜᵃ, c) +
                                    δzᵃᵃᶠ(i, j, k, grid, Az_∂zᵃᵃᶜ, c))
end

#####
##### Horizontal biharmonic operators
#####

@inline ∇⁴hᶜᶜᵃ(i, j, k, grid, c) = ∇²hᶜᶜᵃ(i, j, k, grid, ∇²hᶜᶜᵃ, c)
@inline ∇⁴hᶠᶜᵃ(i, j, k, grid, c) = ∇²hᶠᶜᵃ(i, j, k, grid, ∇²hᶠᶜᵃ, c)
@inline ∇⁴hᶜᶠᵃ(i, j, k, grid, c) = ∇²hᶜᶠᵃ(i, j, k, grid, ∇²hᶜᶠᵃ, c)