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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 0a171174-2db6-4077-a4af-1da84c9f5470
begin
	using Random
	using Plots
	using Distributions
	using PlutoUI

	# Plots.jl backend
	gr()
end;

# ╔═╡ 7b93882c-9ad8-11ea-0288-0941e163f9d5
md"""
# Likelihood function

A broad brush overview.
"""

# ╔═╡ 7110fe00-83a8-4a9d-975e-d59d0f351630
md"""
## Statistical estimation

_Informal definition_: An estimator gives a guess about the true state (parameter, $\boldsymbol{\theta}$) of nature.
"""

# ╔═╡ 13f0fefb-a926-4689-9518-a5d77e0582d7
begin
	# Control and settings
	pop_mean = 5
	pop_std = 1
	n_pop_row = 20
	Random.seed!(122)

	# Population
	pop_distr = Normal(pop_mean, pop_std)
	radius = rand(pop_distr, n_pop_row^2)
	x = 1:n_pop_row
	y = 1:n_pop_row
	x_mx = x' .* ones(n_pop_row)
	y_mx = ones(n_pop_row)' .* y
end;

# ╔═╡ 26c379e4-8480-4e1e-ae8e-933eadfd69f3
md"""
Population radius distribution: ``\mathcal{N}(``$(pop_mean), $(pop_std)``^2)``
"""

# ╔═╡ 016bf804-ca81-4300-8d5b-1e063834d4d0
md"""Hit this button to resample: $(@bind do_sample Button("Sample"))"""

# ╔═╡ 6465a01e-d0e9-4670-b9af-05f7a7a6a200
 md"Change the sample size by dragging the slider: $(@bind n_sample Slider(1:1:(n_pop_row^2), default=20, show_value=true))"

# ╔═╡ 8db552a0-1a45-401c-a579-c3ffc47b27c9
begin
	do_sample
	# Control and settings
	# n_sample = 20

	# ........................
	# Sample
	# ........................
	sample_idx = sample(1:(n_pop_row^2), n_sample, replace=false)
	sample_x = x_mx[:][sample_idx]
	sample_y = y_mx[:][sample_idx]
	sample_radius = radius[sample_idx]

	# with_terminal() do
	# 	println(sample_idx)
	# end

	# ........................
	# Plotting
	# ........................
	begin
		p_pop = scatter(x_mx[:], y_mx[:], markersize=radius[:],
			label=nothing, axis=([], false), title="Population",
		)
		scatter!(sample_x, sample_y, markersize=sample_radius,
			label=nothing, markercolor="red",
		)
		p_sample = scatter(sample_x, sample_y, markersize=sample_radius,
			label=nothing, axis=([], false), markercolor="red", title="Sample",
		)
		plot(p_pop, p_sample, layout = (1, 2), legend = false, link=:all)
	end
end

# ╔═╡ 414cfce1-1bbf-4ae9-9f06-c4a4d55b7ad3
begin
	vline([mean(radius)], lw=3, label="population", yaxis=([], false), xlabel="mean")
	vline!([mean(sample_radius)], color="red", label="sample")
	plot!(xlims=(pop_mean - 3*pop_std, pop_mean + 3*pop_std))
end

# ╔═╡ a090b42b-ee45-4a63-a780-b63df6ac3d11
md"""
## Data generating process

_Informal definition_: The mathematical model that is assumed to generated the data. In statistics it should contain a probabilistic component.
"""

# ╔═╡ 7a0961e0-d332-492e-bbd2-e2127e65f1a4
md"""
For example, considering the above example:

* dataset: ``\mathbf{x}``
* to be estimated model parameter: ``\theta``

$X_i\sim\mathcal{N}(\theta, \sigma^2)$,

where $\sigma$=$(pop_std).

Usually it is more complex, for example a common linear model:

$X_i\sim\mathcal{N}(a \cdot z_i + b, \sigma)$,

where ``\theta=[a, b]``, and ``a`` and ``b`` are the slope and the intersection of the linear model respectively.

Data generating model: Steps of generating random observations from a model (the model definition is part of this!).

__Data generating process/model $\rightarrow$ likelihood function__
"""

# ╔═╡ a9f04175-2d40-4cf6-be59-b24058a4b1ce
md"""
## Likelihood function

_Informal definition_: Tells us what is the likelihood of seeing a dataset under a particular model.

* dataset: ``\mathbf{x}`` (``n`` number of observations)
* model: ``\mathcal{M}``, parametrized with ``\boldsymbol{\theta}``
* particular model: a particular set of ``\boldsymbol{\theta}`` values

$L(\theta, \mathbf{x})$

$L(\theta, \mathbf{x}) = \prod_{i}^{n} p(x_i, \theta)$

where ``p(.)`` is a probability density or mass function.
"""

# ╔═╡ 461c107c-a40e-4fb9-a4ce-9d9c387be074
md"""
## Example: [Task 2](https://docs.google.com/document/d/1vrzjcl8k1RoBaWocgTvVfLLG6LDAHiMwh_MC9PcnfsQ/edit#heading=h.fwqlqf4cfis1)
"""

# ╔═╡ afdd513d-6995-4801-9bac-59d4a2287322
begin
	x_sample = [30.1 28.7 21.2 28.6 25.8 25.1 24.2 27.3 20.0 34.8]
	sigma_pop = 5
	mean_prior_mean = 40
	mean_prior_std = 100
end;

# ╔═╡ 3f467d2d-5053-4e1e-970b-8c61152f5757
md"""
We have a vector independent observations (``n=``$(length(x_sample))):

``\mathbf{x} =`` [$(join(x_sample, ", "))]

drawn from ``\mathcal{N}(\mu, \sigma^2=``$(sigma_pop)``^2)``.

Based on this, we can write down the likelihood function:

$L(\mu, \mathbf{x}) = \prod_{i}^{n} \varphi \left( \frac {x_i - \mu}{\sigma} \right)$

where ``\varphi(.)`` is the standard normal probability density function:

$\varphi(x)={\frac {1}{{\sqrt {2\pi }}}}e^{-x^{2}}$.

Visualizing the log of the likelihood function:
"""

# ╔═╡ 2dfee269-4539-4f96-afca-4e975c2fd5f0
begin
	log_likelihood(t) = sum(logpdf(Normal(t, sigma_pop), x_sample))

	tt = range(30-2.5*sigma_pop, stop=30+2.5*sigma_pop, length=40)
	log_ll = log_likelihood.(tt)

	min_log_l = minimum(log_ll)
end;

# ╔═╡ 9a7b43ec-5d43-43c5-aa35-b776de4a48c7
md"Change the value of ``\mu`` by dragging the slider: $(@bind mu Slider((30-2.5*sigma_pop):0.2:(30+2.5*sigma_pop), default=30))"

# ╔═╡ 39e69cb1-207a-42b6-920a-1912d50edb4e
md"""
``\mu_\mathrm{selected}=``$mu MPa

``\mathrm{log} \left( L(\mu_\mathrm{selected}) \right) =``$(round(log_likelihood(mu), digits=1))
"""

# ╔═╡ 5bd1d628-beaa-47c8-98a6-f8eb36444e93
begin
	log_l_mu = log_likelihood(mu)
	begin
		plot(tt, log_ll, lw=3, xlabel="\\mu", ylabel="log-likelihood", label=nothing)
		plot!([mu, mu], [min_log_l, log_l_mu], ls=:dash, color="red", label="selected value")
		scatter!([mu], [log_l_mu], label=nothing, color="red", markersize=6)
	end
end

# ╔═╡ 684b5e7b-6de1-42af-9497-28a957aa73ad
md"""The maximum likelihood estimate is the value that maximizes the likelihood function (the same that maximizes the log-likelihood function)."""

# ╔═╡ 422e4766-7a75-4403-a189-fbf89d41a8f0
md"""
## Example: Simple linear regression

* data generating process
* constructing the likelihood function
* maximum likelihood estimate

_Task:_ estimate the slope parameter of a linear model from a sample.
"""

# ╔═╡ 5578a9d5-7d7b-4582-91e8-7dab5ba2a572
begin
	# Control and settings
	n_sample_lin = 15
	x_lims = (-2.0, 4.0)
	a_ground_truth = 1.0
	b_ground_truth = 3.0

	# Model
	linear_model(x, a) = a .* x + b_ground_truth

	eps_sigma = 0.5
	ground_truth_model(x) = linear_model(x, a_ground_truth)

	# Generate synthetic data
	x_data = rand(Uniform(x_lims[1], x_lims[2]), n_sample_lin)
	y_data = ground_truth_model.(x_data) .+ rand(Normal(0.0, eps_sigma), n_sample_lin)

	y_lims = ground_truth_model.(x_lims)
	y_range = y_lims[2] - y_lims[1]
	y_lims_lb = minimum([y_lims[1], minimum(y_data)]) - 0.2*y_range
	y_lims_ub = maximum([y_lims[2], maximum(y_data)]) + 0.2*y_range
	y_lims = (y_lims_lb, y_lims_ub)
	# with_terminal() do
	# 	println(y_lims)
	# end
end;

# ╔═╡ 534e1a70-62f4-41fa-bd03-86bb9e4bf495
md"""
Parametrization of the linear model:

$y = a \cdot x + b$

where ``b`` is known, ``b=``$(b_ground_truth), and ``a`` is to be estimated.

We have an ``n=``$n_sample_lin sample of ``x`` and ``y`` pairs:

``\mathbf{x} =`` [$(join(round.(x_data, digits=2), ", "))]

``\mathbf{y} =`` [$(join(round.(y_data, digits=2), ", "))]

This is how our data looks like:
"""

# ╔═╡ b9b80541-571c-4401-b9c9-1ae7fe95beba
begin
	scatter(x_data, y_data, xlabel="x", ylabel="y", label="observation", xlims=x_lims, ylims=y_lims, legend=:topleft)
end

# ╔═╡ 9b63c009-7984-49a6-be7d-5ba52c78a855
md"""
They do not fall on a straight line, so if we want to use a linear model it is better to formulate it with a model error. Let's use the following additive model error formulation:

$y = a \cdot x + b + \epsilon$

where ``\epsilon`` is assumed to be a realization from ``\mathcal{N}(0, \sigma_{\epsilon}=``$(eps_sigma)``^2)``. For simplicity, we assume that this distribution is fully known but in general its standard deviation is not known but also estimated. Moreover, we assume that the ``\epsilon`` values are independent between realizations.

The above equation also shows us the data generating process and we can formulate the likelihood function by recognizing that

$y_i \sim \mathcal{N}(a \cdot x_i + b, \sigma_{\epsilon}^2)$

To make the independent, identically distributed (IID) assumption explicit:

$y_i \overset{\text{IID}}{\sim} \mathcal{N}(a \cdot x_i + b, \sigma_{\epsilon}^2)$

"""

# ╔═╡ bd4dbdc8-84d0-48f8-b4fb-69805b14aaff
md"Change the value of ``a`` by dragging the slider: $(@bind a_selected Slider(-1:0.1:3, default=-0.5, show_value=true))"

# ╔═╡ 32b13ec7-b976-494a-bfff-2efedfe16ea4
begin
	# Plot
	scatter(x_data, y_data, xlabel="x", ylabel="y", label="observation", legend=:topleft, xlims=x_lims, ylims=y_lims)
	plot!([x_lims[1], x_lims[2]], linear_model.([x_lims[1], x_lims[2]], a_selected), label="selected linear model", lw=3)
end

# ╔═╡ cf56ff3d-d9fd-4222-997f-750896ac0440
begin
	# Log-likelihood
	function log_like_fun(t)
		ll = 0
		for ii = 1:n_sample_lin
			ll += logpdf(Normal(linear_model.(x_data[ii], t), eps_sigma), y_data[ii])
		end
		return ll
	end

	tl = range(-1.0, stop=3.0, length=40)
	log_like = log_like_fun.(tl)

	min_log_like = minimum(log_like)
	log_like_a = log_like_fun(a_selected)

	# Plot
	plot(tl, log_like, lw=3, label=nothing, xlabel="a", ylabel="log-likelihood")
	plot!([a_selected, a_selected], [min_log_like, log_like_a], ls=:dash, color="red", label="selected model")
	scatter!([a_selected], [log_like_a], label=nothing, color="red", markersize=6)
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[compat]
Distributions = "~0.25.32"
Plots = "~1.23.6"
PlutoUI = "~0.7.20"
"""

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

[[AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "0bc60e3006ad95b4bb7497698dd7c6d649b9bc06"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.1"

[[Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.1"

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

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

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

[[Bzip2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
version = "1.0.8+0"

[[Cairo_jll]]
deps = ["Artifacts", "Bzip2_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "f2202b55d816427cd385a9a4f3ffb226bee80f99"
uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
version = "1.16.1+0"

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

[[ChangesOfVariables]]
deps = ["LinearAlgebra", "Test"]
git-tree-sha1 = "9a1d594397670492219635b35a3d830b04730d62"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.1"

[[ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random"]
git-tree-sha1 = "a851fec56cb73cfdf43762999ec72eff5b86882a"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.15.0"

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

[[Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

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

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

[[Contour]]
deps = ["StaticArrays"]
git-tree-sha1 = "9f02045d934dc030edad45944ea80dbd1f0ebea7"
uuid = "d38c429a-6771-53c6-b99e-75d170b6e991"
version = "0.5.7"

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

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

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

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

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

[[DensityInterface]]
deps = ["InverseFunctions", "Test"]
git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b"
uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
version = "0.4.0"

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

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

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

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

[[EarCut_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "3f3a2501fa7236e9b911e0f7a588c657e822bb6d"
uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
version = "2.2.3+0"

[[Expat_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b3bfd02e98aedfa5cf885665493c5598c350cd2f"
uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
version = "2.2.10+0"

[[FFMPEG]]
deps = ["FFMPEG_jll"]
git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8"
uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
version = "0.4.1"

[[FFMPEG_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "Pkg", "Zlib_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
git-tree-sha1 = "d8a578692e3077ac998b50c0217dfd67f21d1e5f"
uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
version = "4.4.0+0"

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

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

[[Fontconfig_jll]]
deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "21efd19106a55620a188615da6d3d06cd7f6ee03"
uuid = "a3f928ae-7b40-5064-980b-68af3947d34b"
version = "2.13.93+0"

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

[[FreeType2_jll]]
deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "87eb71354d8ec1a96d4a7636bd57a7347dde3ef9"
uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
version = "2.10.4+0"

[[FriBidi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "aa31987c2ba8704e23c6c8ba8a4f769d5d7e4f91"
uuid = "559328eb-81f9-559d-9380-de523a88c83c"
version = "1.0.10+0"

[[GLFW_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libglvnd_jll", "Pkg", "Xorg_libXcursor_jll", "Xorg_libXi_jll", "Xorg_libXinerama_jll", "Xorg_libXrandr_jll"]
git-tree-sha1 = "0c603255764a1fa0b61752d2bec14cfbd18f7fe8"
uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89"
version = "3.3.5+1"

[[GR]]
deps = ["Base64", "DelimitedFiles", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Printf", "Random", "Serialization", "Sockets", "Test", "UUIDs"]
git-tree-sha1 = "30f2b340c2fff8410d89bfcdc9c0a6dd661ac5f7"
uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
version = "0.62.1"

[[GR_jll]]
deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Pkg", "Qt5Base_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "fd75fa3a2080109a2c0ec9864a6e14c60cca3866"
uuid = "d2c73de3-f751-5644-a686-071e5b155ba9"
version = "0.62.0+0"

[[GeometryBasics]]
deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
git-tree-sha1 = "58bcdf5ebc057b085e58d95c138725628dd7453c"
uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
version = "0.4.1"

[[Gettext_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046"
uuid = "78b55507-aeef-58d4-861c-77aaff3498b1"
version = "0.21.0+0"

[[Glib_jll]]
deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "74ef6288d071f58033d54fd6708d4bc23a8b8972"
uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
version = "2.68.3+1"

[[Graphite2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "344bf40dcab1073aca04aa0df4fb092f920e4011"
uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472"
version = "1.3.14+0"

[[Grisu]]
git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe"
version = "1.0.2"

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "0fa77022fe4b511826b39c894c90daf5fce3334a"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.17"

[[HarfBuzz_jll]]
deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "Graphite2_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg"]
git-tree-sha1 = "129acf094d168394e80ee1dc4bc06ec835e510a3"
uuid = "2e76f6c2-a576-52d4-95c1-20adfe4de566"
version = "2.8.1+1"

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

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

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

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

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

[[InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "a7254c0acd8e62f1ac75ad24d5db43f5f19f3c65"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.2"

[[IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

[[IterTools]]
git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.3.0"

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

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

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[JpegTurbo_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d735490ac75c5cb9f1b00d8b5509c11984dc6943"
uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
version = "2.1.0+0"

[[LAME_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c"
uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d"
version = "3.100.1+0"

[[LZO_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6"
uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
version = "2.10.1+0"

[[LaTeXStrings]]
git-tree-sha1 = "f2355693d6778a178ade15952b7ac47a4ff97996"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.3.0"

[[Latexify]]
deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"]
git-tree-sha1 = "a8f4f279b6fa3c3c4f1adadd78a621b13a506bce"
uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
version = "0.15.9"

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

[[Libffi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "0b4a5d71f3e5200a7dff793393e09dfc2d874290"
uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490"
version = "3.2.2+1"

[[Libgcrypt_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll", "Pkg"]
git-tree-sha1 = "64613c82a59c120435c067c2b809fc61cf5166ae"
uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4"
version = "1.8.7+0"

[[Libglvnd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll", "Xorg_libXext_jll"]
git-tree-sha1 = "7739f837d6447403596a75d19ed01fd08d6f56bf"
uuid = "7e76a0d4-f3c7-5321-8279-8d96eeed0f29"
version = "1.3.0+3"

[[Libgpg_error_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c333716e46366857753e273ce6a69ee0945a6db9"
uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"
version = "1.42.0+0"

[[Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.1+1"

[[Libmount_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "9c30530bf0effd46e15e0fdcf2b8636e78cbbd73"
uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
version = "2.35.0+0"

[[Libtiff_jll]]
deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"]
git-tree-sha1 = "340e257aada13f95f98ee352d316c3bed37c8ab9"
uuid = "89763e89-9b03-5906-acba-b20f662cd828"
version = "4.3.0+0"

[[Libuuid_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7f3efec06033682db852f8b3bc3c1d2b0a0ab066"
uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
version = "2.36.0+0"

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

[[LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "be9eef9f9d78cecb6f262f3c10da151a6c5ab827"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.5"

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

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "3d3e902b31198a27340d0bf00d6ac452866021cf"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.9"

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

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

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

[[Measures]]
git-tree-sha1 = "e498ddeee6f9fdb4551ce855a46f54dbd900245f"
uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
version = "0.3.1"

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

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

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

[[NaNMath]]
git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.5"

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

[[Ogg_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7937eda4681660b4d6aeeecc2f7e1c81c8ee4e2f"
uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
version = "1.3.5+0"

[[OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[OpenSSL_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "15003dcb7d8db3c6c857fda14891a539a8f2705a"
uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
version = "1.1.10+0"

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

[[Opus_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "51a08fb14ec28da2ec7a927c4337e4332c2a4720"
uuid = "91d4177d-7536-5919-b921-800302f37372"
version = "1.3.2+0"

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

[[PCRE_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b2a7af664e098055a7529ad1a900ded962bca488"
uuid = "2f80f16e-611a-54ab-bc61-aa92de5b98fc"
version = "8.44.0+0"

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

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "ae4bbcadb2906ccc085cf52ac286dc1377dceccc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.1.2"

[[Pixman_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b4f5d02549a10e20780a24fce72bea96b6329e29"
uuid = "30392449-352a-5448-841d-b1acce4e97dc"
version = "0.40.1+0"

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

[[PlotThemes]]
deps = ["PlotUtils", "Requires", "Statistics"]
git-tree-sha1 = "a3a964ce9dc7898193536002a6dd892b1b5a6f1d"
uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a"
version = "2.0.1"

[[PlotUtils]]
deps = ["ColorSchemes", "Colors", "Dates", "Printf", "Random", "Reexport", "Statistics"]
git-tree-sha1 = "b084324b4af5a438cd63619fd006614b3b20b87b"
uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
version = "1.0.15"

[[Plots]]
deps = ["Base64", "Contour", "Dates", "Downloads", "FFMPEG", "FixedPointNumbers", "GR", "GeometryBasics", "JSON", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs", "UnicodeFun"]
git-tree-sha1 = "0d185e8c33401084cab546a756b387b15f76720c"
uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
version = "1.23.6"

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

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

[[Qt5Base_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "xkbcommon_jll"]
git-tree-sha1 = "ad368663a5e20dbb8d6dc2fddeefe4dae0781ae8"
uuid = "ea2cea3b-5b76-57ae-a6ef-0a8af62496e1"
version = "5.15.3+0"

[[QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "78aadffb3efd2155af139781b8a8df1ef279ea39"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.2"

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

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

[[RecipesBase]]
git-tree-sha1 = "a4425fe1cde746e278fa895cc69e3113cb2614f6"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.2.0"

[[RecipesPipeline]]
deps = ["Dates", "NaNMath", "PlotUtils", "RecipesBase"]
git-tree-sha1 = "7ad0dfa8d03b7bcf8c597f59f5292801730c55b8"
uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c"
version = "0.4.1"

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

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

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

[[Scratch]]
deps = ["Dates"]
git-tree-sha1 = "0b4b7f1393cff97c33891da2a0bf69c6ed241fda"
uuid = "6c6a2e73-6563-6170-7368-637461726353"
version = "1.1.0"

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

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

[[Showoff]]
deps = ["Dates", "Grisu"]
git-tree-sha1 = "91eddf657aca81df9ae6ceb20b959ae5653ad1de"
uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f"
version = "1.0.3"

[[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", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "f0bccf98e16759818ffc5d97ac3ebf87eb950150"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.8.1"

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3c76dde64d03699e074ac02eb2e8ba8254d428da"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.13"

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

[[StatsAPI]]
git-tree-sha1 = "0f2aa8e32d511f758a2ce49208181f7733a0936a"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.1.0"

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

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

[[StructArrays]]
deps = ["Adapt", "DataAPI", "StaticArrays", "Tables"]
git-tree-sha1 = "2ce41e0d042c60ecd131e9fb7154a3bfadbf50d3"
uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
version = "0.6.3"

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

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

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

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

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

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

[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

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

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

[[UnicodeFun]]
deps = ["REPL"]
git-tree-sha1 = "53915e50200959667e78a92a418594b428dffddf"
uuid = "1cfade01-22cf-5700-b092-accc4b62d6e1"
version = "0.4.1"

[[Wayland_jll]]
deps = ["Artifacts", "Expat_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "3e61f0b86f90dacb0bc0e73a0c5a83f6a8636e23"
uuid = "a2964d1f-97da-50d4-b82a-358c7fce9d89"
version = "1.19.0+0"

[[Wayland_protocols_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "66d72dc6fcc86352f01676e8f0f698562e60510f"
uuid = "2381bf8a-dfd0-557d-9999-79630e7b1b91"
version = "1.23.0+0"

[[XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.12+0"

[[XSLT_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"]
git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a"
uuid = "aed1982a-8fda-507f-9586-7b0439959a61"
version = "1.1.34+0"

[[Xorg_libX11_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
git-tree-sha1 = "5be649d550f3f4b95308bf0183b82e2582876527"
uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"
version = "1.6.9+4"

[[Xorg_libXau_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4e490d5c960c314f33885790ed410ff3a94ce67e"
uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec"
version = "1.0.9+4"

[[Xorg_libXcursor_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXfixes_jll", "Xorg_libXrender_jll"]
git-tree-sha1 = "12e0eb3bc634fa2080c1c37fccf56f7c22989afd"
uuid = "935fb764-8cf2-53bf-bb30-45bb1f8bf724"
version = "1.2.0+4"

[[Xorg_libXdmcp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fe47bd2247248125c428978740e18a681372dd4"
uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05"
version = "1.1.3+4"

[[Xorg_libXext_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "b7c0aa8c376b31e4852b360222848637f481f8c3"
uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
version = "1.3.4+4"

[[Xorg_libXfixes_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "0e0dc7431e7a0587559f9294aeec269471c991a4"
uuid = "d091e8ba-531a-589c-9de9-94069b037ed8"
version = "5.0.3+4"

[[Xorg_libXi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXfixes_jll"]
git-tree-sha1 = "89b52bc2160aadc84d707093930ef0bffa641246"
uuid = "a51aa0fd-4e3c-5386-b890-e753decda492"
version = "1.7.10+4"

[[Xorg_libXinerama_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll"]
git-tree-sha1 = "26be8b1c342929259317d8b9f7b53bf2bb73b123"
uuid = "d1454406-59df-5ea1-beac-c340f2130bc3"
version = "1.1.4+4"

[[Xorg_libXrandr_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll"]
git-tree-sha1 = "34cea83cb726fb58f325887bf0612c6b3fb17631"
uuid = "ec84b674-ba8e-5d96-8ba1-2a689ba10484"
version = "1.5.2+4"

[[Xorg_libXrender_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "19560f30fd49f4d4efbe7002a1037f8c43d43b96"
uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa"
version = "0.9.10+4"

[[Xorg_libpthread_stubs_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "6783737e45d3c59a4a4c4091f5f88cdcf0908cbb"
uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74"
version = "0.1.0+3"

[[Xorg_libxcb_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
git-tree-sha1 = "daf17f441228e7a3833846cd048892861cff16d6"
uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b"
version = "1.13.0+3"

[[Xorg_libxkbfile_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "926af861744212db0eb001d9e40b5d16292080b2"
uuid = "cc61e674-0454-545c-8b26-ed2c68acab7a"
version = "1.1.0+4"

[[Xorg_xcb_util_image_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "0fab0a40349ba1cba2c1da699243396ff8e94b97"
uuid = "12413925-8142-5f55-bb0e-6d7ca50bb09b"
version = "0.4.0+1"

[[Xorg_xcb_util_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll"]
git-tree-sha1 = "e7fd7b2881fa2eaa72717420894d3938177862d1"
uuid = "2def613f-5ad1-5310-b15b-b15d46f528f5"
version = "0.4.0+1"

[[Xorg_xcb_util_keysyms_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "d1151e2c45a544f32441a567d1690e701ec89b00"
uuid = "975044d2-76e6-5fbe-bf08-97ce7c6574c7"
version = "0.4.0+1"

[[Xorg_xcb_util_renderutil_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "dfd7a8f38d4613b6a575253b3174dd991ca6183e"
uuid = "0d47668e-0667-5a69-a72c-f761630bfb7e"
version = "0.3.9+1"

[[Xorg_xcb_util_wm_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "e78d10aab01a4a154142c5006ed44fd9e8e31b67"
uuid = "c22f9ab0-d5fe-5066-847c-f4bb1cd4e361"
version = "0.4.1+1"

[[Xorg_xkbcomp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxkbfile_jll"]
git-tree-sha1 = "4bcbf660f6c2e714f87e960a171b119d06ee163b"
uuid = "35661453-b289-5fab-8a00-3d9160c6a3a4"
version = "1.4.2+4"

[[Xorg_xkeyboard_config_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xkbcomp_jll"]
git-tree-sha1 = "5c8424f8a67c3f2209646d4425f3d415fee5931d"
uuid = "33bec58e-1273-512f-9401-5d533626f822"
version = "2.27.0+4"

[[Xorg_xtrans_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "79c31e7844f6ecf779705fbc12146eb190b7d845"
uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10"
version = "1.4.0+3"

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

[[Zstd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "cc4bf3fdde8b7e3e9fa0351bdeedba1cf3b7f6e6"
uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
version = "1.5.0+0"

[[libass_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "HarfBuzz_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "5982a94fcba20f02f42ace44b9894ee2b140fe47"
uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0"
version = "0.15.1+0"

[[libfdk_aac_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "daacc84a041563f965be61859a36e17c4e4fcd55"
uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280"
version = "2.0.2+0"

[[libpng_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c"
uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
version = "1.6.38+0"

[[libvorbis_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"]
git-tree-sha1 = "c45f4e40e7aafe9d086379e5578947ec8b95a8fb"
uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a"
version = "1.3.7+0"

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

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

[[x264_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fea590b89e6ec504593146bf8b988b2c00922b2"
uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a"
version = "2021.5.5+0"

[[x265_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "ee567a171cce03570d77ad3a43e90218e38937a9"
uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76"
version = "3.5.0+0"

[[xkbcommon_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll", "Wayland_protocols_jll", "Xorg_libxcb_jll", "Xorg_xkeyboard_config_jll"]
git-tree-sha1 = "ece2350174195bb31de1a63bea3a41ae1aa593b6"
uuid = "d8fb68d0-12a3-5cfd-a85a-d49703b185fd"
version = "0.9.1+5"
"""

# ╔═╡ Cell order:
# ╟─7b93882c-9ad8-11ea-0288-0941e163f9d5
# ╟─7110fe00-83a8-4a9d-975e-d59d0f351630
# ╟─0a171174-2db6-4077-a4af-1da84c9f5470
# ╟─13f0fefb-a926-4689-9518-a5d77e0582d7
# ╟─26c379e4-8480-4e1e-ae8e-933eadfd69f3
# ╟─016bf804-ca81-4300-8d5b-1e063834d4d0
# ╟─6465a01e-d0e9-4670-b9af-05f7a7a6a200
# ╟─8db552a0-1a45-401c-a579-c3ffc47b27c9
# ╟─414cfce1-1bbf-4ae9-9f06-c4a4d55b7ad3
# ╟─a090b42b-ee45-4a63-a780-b63df6ac3d11
# ╟─7a0961e0-d332-492e-bbd2-e2127e65f1a4
# ╟─a9f04175-2d40-4cf6-be59-b24058a4b1ce
# ╟─461c107c-a40e-4fb9-a4ce-9d9c387be074
# ╟─afdd513d-6995-4801-9bac-59d4a2287322
# ╟─3f467d2d-5053-4e1e-970b-8c61152f5757
# ╟─2dfee269-4539-4f96-afca-4e975c2fd5f0
# ╟─9a7b43ec-5d43-43c5-aa35-b776de4a48c7
# ╟─39e69cb1-207a-42b6-920a-1912d50edb4e
# ╟─5bd1d628-beaa-47c8-98a6-f8eb36444e93
# ╟─684b5e7b-6de1-42af-9497-28a957aa73ad
# ╟─422e4766-7a75-4403-a189-fbf89d41a8f0
# ╟─5578a9d5-7d7b-4582-91e8-7dab5ba2a572
# ╟─534e1a70-62f4-41fa-bd03-86bb9e4bf495
# ╟─b9b80541-571c-4401-b9c9-1ae7fe95beba
# ╟─9b63c009-7984-49a6-be7d-5ba52c78a855
# ╟─bd4dbdc8-84d0-48f8-b4fb-69805b14aaff
# ╟─32b13ec7-b976-494a-bfff-2efedfe16ea4
# ╟─cf56ff3d-d9fd-4222-997f-750896ac0440
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
abstract type DiagonalizationType end

function ed_size end
function spectral_size end
function diagonalize! end

struct Full <: DiagonalizationType end
Base.:(==)(::Full, ::Full) = true

ed_size(::Full) = 1
spectral_size(::Full, model) = basissize(model.basis)

function diagonalize!(callback, H, ::Full)
    callback(1, eigen!(Hermitian(Matrix(H)))...)
end

abstract type SparseMethod end

struct ShiftInvertARPACK <: SparseMethod end
struct ShiftInvert <: SparseMethod end
struct POLFED <: SparseMethod end

struct Sparse{M<:SparseMethod} <: DiagonalizationType
    method::M
    σ::Vector{Float64}
    count::Int64
    normalize::Bool
end
#Base.:(==)(s1::Sparse, s2::Sparse) = s1.σ == s2.σ && s1.count == s2.count

spectral_size(sp::Sparse, model) = sp.count
ed_size(sp::Sparse) = length(sp.σ)

function diagonalize!(callback, H, sp::Sparse{ShiftInvertARPACK})
    if sp.normalize
        emin, emax = eigsolve(H, 1, :SR)[1][1], eigsolve(H, 1, :LR)[1][1]
        H .= (H - emin*I) ./ (emax-emin)
    end
    for (i, σ) in enumerate(sp.σ)
        (d,v,nconv,niter,nmult,resid) = eigs(H; nev=sp.count, sigma=σ, check=1, tol=1e-12)
        logmsg("$nconv converged after $niter iterations.")
        callback(i, LinearAlgebra.sorteig!(d,v)...)
    end
end

function diagonalize!(callback, H, sp::Sparse{ShiftInvert})
    if sp.normalize
        emin, emax = eigsolve(H, 1, :SR)[1][1], eigsolve(H, 1, :LR)[1][1]
        H .= (H - emin*I) ./ (emax-emin)
    end

    L = size(H, 1)
    for (i, σ) in enumerate(sp.σ)
        iH = let F = factorize(H-σ*I)
            #LinearMap{eltype(H)}(x -> F \ x, L, ismutating=false, issymmetric=true)
            x -> F \ x
        end
        evals, evecs, info = eigsolve(iH, L, sp.count; ishermitian=true, krylovdim=max(20,2*sp.count+1))
        @. evals = σ + 1 / evals
        logmsg("Shift-invert: $(info.converged) converged after $(info.numiter) iterations.")
        callback(i, LinearAlgebra.sorteig!(evals, hcat(evecs...))...)
    end
end

function diagonalize!(callback, H, sp::Sparse{POLFED})
    emin, emax = eigsolve(H, 1, :SR)[1][1], eigsolve(H, 1, :LR)[1][1]

    if sp.normalize
        H .= (H - emin*I) ./ (emax-emin)
        emin, emax = 0.0, 1.0
    end

    L = size(H, 1)
    Δϵ = (emax-emin)/size(H, 1) # average level spacing

    for (i, σ) in enumerate(sp.σ)
        # this order estimate needs to be improved!
        order = order_estimate(σ, Δϵ, (emin,emax))
        δH = delta(H, σ, (emin, emax); order)
        _, vecs, info = eigsolve(δH, L, sp.count; ishermitian=true, krylovdim=max(20,2*sp.count+1))
        logmsg("POLFED with order $order: $(info.converged) converged after $(info.numiter) iterations.")
        # compute eigenvals from vectors:
        vals = dot.(vecs, Ref(H), vecs) ./ dot.(vecs, vecs)

        callback(i, LinearAlgebra.sorteig!(vals, hcat(vecs...)))
    end
end

ext: jl


text:
module ElectromagneticFields

    using Documenter
    using LinearAlgebra

    export MagneticEquilibrium

    include("equilibrium.jl")

    export AnalyticEquilibrium, AnalyticPerturbation, ZeroPerturbation

    export CartesianEquilibrium, CartesianPerturbation

    export load_equilibrium, periodicity

    include("analytic/analytic_equilibrium.jl")
    include("analytic/cartesian_equilibrium.jl")

    export ABC, AxisymmetricTokamakCartesian, AxisymmetricTokamakCylindrical,
           AxisymmetricTokamakToroidal, EzCosZ,
           Solovev, SolovevXpoint, SolovevQuadratic, SymmetricQuadratic, ThetaPinch
    export SolovevITER, SolovevNSTX, SolovevXpointITER, SolovevXpointNSTX
    export @abc_equilibrium,
           @axisymmetric_tokamak_equilibrium_cartesian,
           @axisymmetric_tokamak_equilibrium_cylindrical,
           @axisymmetric_tokamak_equilibrium_toroidal,
           @ezcosz_perturbation,
           @solovev_equilibrium,
           @solovev_xpoint_equilibrium,
           @solovev_equilibrium_quadratic,
           @symmetric_quadratic_equilibrium,
           @theta_pinch_equilibrium

    include("analytic/abc.jl")
    include("analytic/axisymmetric_tokamak_cartesian.jl")
    include("analytic/axisymmetric_tokamak_cylindrical.jl")
    include("analytic/axisymmetric_tokamak_toroidal.jl")
    include("analytic/ezcosz.jl")
    include("analytic/solovev_common.jl")
    include("analytic/solovev_psi.jl")
    include("analytic/solovev.jl")
    include("analytic/solovev_quadratic.jl")
    include("analytic/solovev_xpoint.jl")
    include("analytic/symmetric_quadratic.jl")
    include("analytic/theta_pinch.jl")

end


text:
# Receta para graficar función de inflación con datos

@recipe function plot(inflfn::InflationFunction, data::CountryStructure)

    # Computar trayectoria de inflación
    traj_infl = inflfn(data)

    # Etiquetas para EnsembleFunction
    _label = inflfn isa EnsembleFunction ?
        reshape(measure_name(inflfn), 1, :) : measure_name(inflfn)

    label --> _label
    legend --> :topright

    infl_dates(data), traj_infl
end


text:
# Sat 29 Dec 2012 11:48:33 PM EST
#
# NeHe Tut 12 - Display Lists
#
# Q - quit
# Up/Down - increase/decrease rotation of each cube about it's own x-axis
# Left/Right - increase/decrease rotation of each cube about it's own y-axis

# load necessary GL/SDL routines

using OpenGL
@OpenGL.version "1.0"
@OpenGL.load
using SDL

# initialize variables

bpp       = 16
wintitle  = "NeHe Tut 12"
icontitle = "NeHe Tut 12"
width     = 640
height    = 480

T0        = 0
Frames    = 0

xrot      = 0
yrot      = 0

boxcol    = [1.0 0.0 0.0;
             1.0 0.5 0.0;
             1.0 1.0 0.0;
             0.0 1.0 0.0;
             0.0 1.0 1.0]

topcol    = [0.5 0.0  0.0;
             0.5 0.25 0.0;
             0.5 0.5  0.0;
             0.0 0.5  0.0;
             0.0 0.5  0.5]

keystate_checked   = false
lastkeycheckTime   = 0
key_repeatinterval = 75 #ms

# open SDL window with an OpenGL context

SDL_Init(SDL_INIT_VIDEO)
#videoInfo = SDL_GetVideoInfo()
videoFlags = (SDL_OPENGL | SDL_GL_DOUBLEBUFFER | SDL_HWPALETTE | SDL_RESIZABLE)
#if videoInfo.hw_available
    videoFlags = (videoFlags | SDL_HWSURFACE)
#else
    #videoFlags = (videoFlags | SDL_SWSURFACE)
#end
#if videoInfo.blit_hw
    videoFlags = (videoFlags | SDL_HWACCEL)
#end
SDL_GL_SetAttribute(SDL_GL_DOUBLEBUFFER, 1)
SDL_SetVideoMode(width, height, bpp, videoFlags)
SDL_WM_SetCaption(wintitle, icontitle)

glViewport(0, 0, width, height)
glClearColor(0.0, 0.0, 0.0, 0.5)
glClearDepth(1.0)
glDepthFunc(GL_LEQUAL)
glEnable(GL_DEPTH_TEST)
glShadeModel(GL_SMOOTH)
glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST)

#enable simple lighting

glEnable(GL_LIGHT0)
glEnable(GL_LIGHTING)
glEnable(GL_COLOR_MATERIAL)

glMatrixMode(GL_PROJECTION)
glLoadIdentity()

gluPerspective(45.0,width/height,0.1,100.0)

glMatrixMode(GL_MODELVIEW)

# build the display lists

box = glGenLists(2)

glNewList(box, GL_COMPILE)
    glBegin(GL_QUADS)
        # Bottom Face
        glTexCoord(0.0, 1.0)
        glVertex(-1.0, -1.0, -1.0)
        glTexCoord(1.0, 1.0)
        glVertex(1.0, -1.0, -1.0)
        glTexCoord(1.0, 0.0)
        glVertex(1.0, -1.0,  1.0)
        glTexCoord(0.0, 0.0)
        glVertex(-1.0, -1.0,  1.0)

        # Front Face
        glTexCoord(1.0, 0.0)
        glVertex(-1.0, -1.0,  1.0)
        glTexCoord(0.0, 0.0)
        glVertex(1.0, -1.0,  1.0)
        glTexCoord(0.0, 1.0)
        glVertex(1.0,  1.0,  1.0)
        glTexCoord(1.0, 1.0)
        glVertex(-1.0,  1.0,  1.0)

        # Back Face
        glTexCoord(0.0, 0.0)
        glVertex(-1.0, -1.0, -1.0)
        glTexCoord(0.0, 1.0)
        glVertex(-1.0,  1.0, -1.0)
        glTexCoord(1.0, 1.0)
        glVertex(1.0,  1.0, -1.0)
        glTexCoord(1.0, 0.0)
        glVertex(1.0, -1.0, -1.0)

        # Right Face
        glTexCoord(0.0, 0.0)
        glVertex(1.0, -1.0, -1.0)
        glTexCoord(0.0, 1.0)
        glVertex(1.0,  1.0, -1.0)
        glTexCoord(1.0, 1.0)
        glVertex(1.0,  1.0,  1.0)
        glTexCoord(1.0, 0.0)
        glVertex(1.0, -1.0,  1.0)

        # Left Face
        glTexCoord(1.0, 0.0)
        glVertex(-1.0, -1.0, -1.0)
        glTexCoord(0.0, 0.0)
        glVertex(-1.0, -1.0,  1.0)
        glTexCoord(0.0, 1.0)
        glVertex(-1.0,  1.0,  1.0)
        glTexCoord(1.0, 1.0)
        glVertex(-1.0,  1.0, -1.0)
    glEnd()
glEndList()

top = uint32(box+1)

glNewList(top, GL_COMPILE)
    glBegin(GL_QUADS)
        # Top Face
        glTexCoord(1.0, 1.0)
        glVertex(-1.0, 1.0, -1.0)
        glTexCoord(1.0, 0.0)
        glVertex(-1.0, 1.0,  1.0)
        glTexCoord(0.0, 0.0)
        glVertex(1.0, 1.0,  1.0)
        glTexCoord(0.0, 1.0)
        glVertex(1.0, 1.0, -1.0)
    glEnd()
glEndList()

# load textures from images

tex   = Array(Uint32,1) # generating 1 texture

img, w, h = glimread(expanduser("~/.julia/SDL/Examples/NeHe/tut12/cube.bmp"))

glGenTextures(1,tex)
glBindTexture(GL_TEXTURE_2D,tex[1])
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR_MIPMAP_NEAREST)
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR)
glTexImage2D(GL_TEXTURE_2D, 0, 3, w, h, 0, GL_RGB, GL_UNSIGNED_BYTE, img)
gluBuild2DMipmaps(GL_TEXTURE_2D, 3, w, h, GL_RGB, GL_UNSIGNED_BYTE, img)

# enable texture mapping

glEnable(GL_TEXTURE_2D)

# main drawing loop

while true
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

    glBindTexture(GL_TEXTURE_2D, tex[1])

    for yloop = 1:5
	      for xloop = 1:yloop
            glLoadIdentity()

            glTranslate(1.4+2.8xloop-1.4yloop, ((6.0-yloop)*2.4)-7.0, -20.0)

            glRotate(45.0-(2.0yloop)+xrot, 1.0, 0.0, 0.0)
            glRotate(45.0+yrot, 0.0, 1.0, 0.0)

            glColor(boxcol[yloop,:])
            glCallList(box)

            glColor(topcol[yloop,:])
            glCallList(top)
        end
    end

    SDL_GL_SwapBuffers()

    Frames +=1
    t = SDL_GetTicks()
    if (t - T0) >= 5000
        seconds = (t - T0)/1000
        fps = Frames/seconds
        println("$Frames frames in $seconds seconds = $fps")
        T0 = t
        Frames = 0
    end

    SDL_PumpEvents()
    if SDL_GetTicks() - lastkeycheckTime >= key_repeatinterval
        keystate         = SDL_GetKeyState()
        keystate_checked = true
        lastkeycheckTime = SDL_GetTicks()
    end

    # Sampling rates for event processing are so fast that a single key press
    # lasts through several iterations of this loop.  This means that one press
    # can be seen as 50 or more presses by the SDL event system.  To correct
    # for this, we only check keypresses every 75ms.

    if keystate_checked == true
        if keystate[SDLK_q] == true
            break
        elseif keystate[SDLK_LEFT] == true
            yrot -=0.2
        elseif keystate[SDLK_RIGHT] == true
            yrot +=0.2
        elseif keystate[SDLK_UP] == true
            xrot -=0.2
        elseif keystate[SDLK_DOWN] == true
            xrot +=0.2
        end
        keystate_checked = false
    end
end


text:
#=
array_initialization:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-10
=#

module ArrayInit
    # methods
    function main()
        println("Not init array = ", Array{Float32}(undef, 3, 2))
        println("Init 0 array = ", zeros(Float32, 3, 2))
        println("Init 1 array = ", ones(Float32, 3, 2))
        println("Init uniform distribution = ", rand(Float32, 3, 2))
        println("Init normal distribution = ", randn(Float32, 3, 2))
        println("Init fill 8 = ", fill(8, 3, 2))
        println("Init zeros fill 10 = ", fill!(zeros(Int32, 4, 4), 10))
        println("Similar not init = ", similar(fill(8, 4, 4)))
        println("Similar init float32 = ", similar(fill(8, 4, 4), Float32))
        println("Similar init float32 3x2 = ", similar(fill(8, 4, 4), Float32, 3, 2))
        println("Reshape 3x2 -> 2x3: ", reshape(zeros(Float32, 3, 2), 2, 3))
        println("List Comprehension [x for x in 1:10] = ", [x for x in 1:10])
        println("List Comprehension [10a + b for a in 0:3, b in 0:4] = ", [10a + b for a in 0:3, b in 0:4])
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .ArrayInit
    ArrayInit.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
"""
    ginv(sys; type = "1-2", mindeg = false, fast = true, atol = 0, atol1 = atol, atol2 = atol, rtol,
              offset = sqrt(ϵ)) -> (sysinv, info)

Compute for a descriptor system `sys = (A-λE,B,C,D)` with the transfer function matrix `G(λ)`
a generalized inverse system `sysinv = (Ai-λEi,Bi,Ci,Di)` with the transfer function matrix `Gi(λ)`
such that two or more of the following _Moore-Penrose conditions_ are satisfied:

         (1) G(λ)*Gi(λ)*G(λ) = G(λ);
         (2) Gi(λ)*G(λ)*Gi(λ) = Gi(λ);
         (3) G(λ)*Gi(λ) = (G(λ)*Gi(λ))';
         (4) Gi(λ)*G(λ) = (Gi(λ)*G(λ))'.

The desired type of the computed generalized inverse can be specified using the keyword parameter
`type` as follows:

      "1-2"     - for a generalized inverse which satisfies conditions (1) and (2) (default);
      "1-2-3"   - for a generalized inverse which satisfies conditions (1), (2) and (3);
      "1-2-4"   - for a generalized inverse which satisfies conditions (1), (2) and (4);
      "1-2-3-4" - for the Moore-Penrose pseudoinverse, which satisfies all conditions (1)-(4).

The vector `poles` specified as a keyword argument, can be used to specify the desired eigenvalues
alternatively to or jointly with enforcing a desired stability degree `sdeg` of the poles of the
computed generalized inverse.

The keyword argument `mindeg` can be used to specify the option to determine a minimum order
generalized inverse, if `mindeg = true`, or a particular generalized inverse which has
possibly non-minimal order, if `mindeg = false` (default).

To assess the presence of zeros on the boundary of the stability domain `Cs`, a boundary offset  `β`
can be specified via the keyword parameter `offset = β`.
Accordingly, for a continuous-time setting,
the boundary of `Cs` contains the complex numbers with real parts within the interval `[-β,β]`,
while for a discrete-time setting, the boundary of `Cs` contains
the complex numbers with moduli within the interval `[1-β,1+β]`.
The default value used for `β` is `sqrt(ϵ)`, where `ϵ` is the working machine precision.

The returned named tuple `info` has the components `info.nrank`, `info.nfp`, `info.fnorm` and `info.tcond`,
where:
`info.nrank` is the normal rank of `G(λ)`,
`info.nfp` is the number of freely assignable poles of the inverse `Gi(λ)`,
`info.fnorm` is the maximum of norms of employed feedback gains (also for pole assignment) and
`info.tcond` is the maximum of condition numbers of employed non-orthogonal transformations (see below).

The rank determinations in the performed reductions
 are based on rank revealing QR-decompositions with column pivoting,
if `fast = true`, or the more reliable SVD-decompositions, if `fast = false`.
The computation of a minimum order inverse is performed by solving suitable minimum dynamic cover
problems. These computations involve using
non-orthogonal transformations whose maximal condition number is returned in `info.tcond`, in conjunction with
using feedback gains (also for pole assignment), whose maximal norms are returned in `info.fnorm`.
High values of these quantities indicate a potential loss of numerical stability of computations.

The keyword arguments `atol1`, `atol2` and `rtol`, specify, respectively,
the absolute tolerance for the nonzero elements of `A`, `B`, `C` and `D`,
the absolute tolerance for the nonzero elements of `E`,
and the relative tolerance for the nonzero elements of `A`, `E`, `B`, `C` and `D`.
The default relative tolerance is `n*ϵ`, where `ϵ` is the working machine epsilon
and `n` is the maximum dimension of state, input and output vectors of the system `sys`.
The keyword argument `atol` can be used to simultaneously set `atol1 = atol` and `atol2 = atol`.

_Method:_  The methods proposed in [1] are employed in conjunction with
full rank factorizations computed using the approach of [2].

_References:_

[1] A. Varga. Computing generalized inverse systems using matrix pencil methods.
    Int. J. of Applied Mathematics and Computer Science, vol. 11, pp. 1055-1068, 2001.

[2] Varga, A. A note on computing range space bases of rational matrices.
       arXiv:1707.0048, [https://arxiv.org/abs/1707.00489](https://arxiv.org/abs/1707.00489), 2017.

"""
function ginv(sys::DescriptorStateSpace{T}; type::String = "1-2", fast::Bool = true, mindeg::Bool = false,
              poles::Union{AbstractVector,Missing} = missing, sdeg::Union{Real,Missing} = missing,
              atol::Real = zero(real(T)), atol1::Real = atol, atol2::Real = atol,
              rtol::Real = (max(sys.nx,sys.nu,sys.ny)*eps(real(float(one(T)))))*iszero(min(atol1,atol2)),
              offset::Real = sqrt(eps(float(real(T))))) where T
   p, m = size(sys)
   Ts = sys.Ts;

   T1 = T <: BlasFloat ? T : promote_type(Float64,T)
   ONE = one(T1)

   disc = !iszero(Ts)
   sdeg_nomissing = !ismissing(sdeg)
   poles_nomissing = !ismissing(poles)

   stabilize = sdeg_nomissing || poles_nomissing

   # set default values of sdeg if poles = missing
   # sdeg_nomissing || (sdeg = disc ? real(T1)(0.95) : real(T1)(-0.05))

   # sort desired eigenvalues
   if poles_nomissing
      tempc = poles[imag.(poles) .> 0]
      if !isempty(tempc)
         tempc1 = conj(poles[imag.(poles) .< 0])
         isequal(tempc[sortperm(real(tempc))],tempc1[sortperm(real(tempc1))]) ||
                 error("poles must be a self-conjugated complex vector")
      end
      # check that all eigenvalues are inside of the stability region
      sdeg_nomissing && ( ((disc && any(abs.(poles) .> sdeg) )  || (!disc && any(real.(poles) .> sdeg)))  &&
            error("The elements of poles must lie in the stability region of interest") )
   end

   nrank = gnrank(sys; atol1, atol2, rtol);
   nfp = 0;
   tcond = 1;
   fnorm = 0;
   if p == m && m == nrank
      # compute a standard inverse for an invertible system
      sysinv = inv(sys);
      info = (nrank = nrank, nfp = nfp, tcond = 1, fnorm = 0)
      return sysinv, info
   end
   disc = !iszero(sys.Ts)
   if type == "1-2"
      # compute an (1,2)-inverse
      if min(p,m) == nrank
         if m == nrank
            # compute a left inverse
            sysinv, info1 = glsol([sys; I], m; atol1, atol2, rtol, fast, poles, sdeg, mindeg)
            info1.nrank != nrank && (@warn "Inconsistent rank evaluations: check tolerances")
            mindeg || (nfp = info1.nl)
         else
            # compute a right inverse
            sysinv, info1 = grsol([sys I], p; atol1, atol2, rtol, fast, poles, sdeg, mindeg)
            info1.nrank != nrank && (@warn "Inconsistent rank evaluations: check tolerances")
            mindeg || (nfp = info1.nr)
         end
         tcond = info1.tcond;
         fnorm = info1.fnorm;
      else
         if mindeg
            # compute a full rank factorization G = U*V
            U, V, info1 = grange(sys; atol1, atol2, rtol, offset, fast)
            info1.nrank != nrank && (@warn "Inconsistent rank evaluations: check tolerances")
            # compute a left inverse of least order of U
            UL, info1, = glsol([U; I], nrank; atol1, atol2, rtol, fast, poles, sdeg, mindeg)
            tcond = info1.tcond;
            fnorm = info1.fnorm;
            # compute a right inverse of least order of V
            VR, info1, = grsol([V I], nrank; atol1, atol2, rtol, fast, poles, sdeg, mindeg)
            tcond = max(tcond, info1.tcond);
            fnorm = max(fnorm, info1.fnorm);
            sysinv = gminreal(VR*UL; atol1, atol2, rtol)
         else
            # compute a generalized 1-2 inverse using the Kronecker form
            # of the system pencil for a system with no full rank TFM
            a, e, b, c, d = dssdata(T1,sys)
            n = size(a,1);
            A, E, Q, Z, νr, μr, νi, nf, νl, μl = sklf(a, e, b, c, d; atol1, atol2, rtol, fast)
            B = Q[n+1:end,:]'; C =  -Z[n+1:end,:]; D = zeros(T1,m,p);
            nreg = nf+sum(νi); nr1 = sum(νr); nl = sum(μl); ninv = nr1+nreg+nl;
            mr = n+m-ninv; # pr = n+p-ninv;
            if stabilize
               # make spurious zeros stable
               if nr1 > 0
                  i1 = 1:nr1; j1 = 1:mr; j2 = mr+1:mr+nr1;
                  F, = saloc(A[i1,j2], E[i1,j2], A[i1,j1]; evals = poles, sdeg, disc,
                             atol1, atol2, atol3 = atol1, rtol, sepinf = false, fast)
                  mul!(view(A,i1,j2), view(A,i1,j1), F, ONE, ONE)
                  mul!(view(C,:,j2), view(C,:,j1), F, ONE, ONE)
                  fnorm = norm(F)
               end
               if nl > 0
                  i1 = nr1+nreg+1:nr1+nreg+nl; j1=n+m-nl+1:n+m; i2 = ninv+1:n+p;
                  K, = saloc(A[i1,j1]', E[i1,j1]', A[i2,j1]'; evals = poles, sdeg, disc,
                             atol1, atol2, atol3 = atol1, rtol, sepinf = false, fast); K = copy(K')
                  mul!(view(A,i1,j1), K, view(A,i2,j1), ONE, ONE)
                  mul!(view(B,i1,:), K, view(B,i2,:), ONE, ONE)
                  fnorm = max(fnorm,norm(K))
               end
            end
            nfp = nr1+nl;
            # select a square system
            ic = mr+1:n+m; ir = 1:ninv;
            sysinv = gminreal(dss(view(A,ir,ic), view(E,ir,ic), view(B,ir,:), view(C,:,ic), D; Ts);
                              atol1, atol2, rtol)
         end
      end
   elseif type == "1-2-3"
      # compute an (1,2,3)-inverse
      if p == nrank
         # for full row rank, any (1,2)-inverse is an (1,2,3)-inverse
         sysinv, info1 = grsol([sys I], p; atol1, atol2, rtol, fast, poles, sdeg, mindeg)
         mindeg || (nfp = info1.nr)
      else
         # compute the full-rank factorization G = U*G1 with U inner
         U, G1, info1 = grange(sys; inner = true, atol1, atol2, rtol, offset, fast)
         info1.nrank != nrank && (@warn "Inconsistent rank evaluations: check tolerances")
         G2, info1 = grsol([ G1 I],nrank; atol1, atol2, rtol, fast, poles, sdeg, mindeg)
         mindeg || (nfp = info1.nr)
         sysinv = gminreal(G2*U'; atol1, atol2, rtol)
      end
      tcond = info1.tcond;
      fnorm = info1.fnorm;
   elseif type == "1-2-4"
      # compute an (1,2,4)-inverse
      if m == nrank
         # for full column rank, any (1,2)-inverse is an (1,2,4)-inverse
         sysinv, info1 = glsol([sys; I], m; atol1, atol2, rtol, fast, poles, sdeg, mindeg)
         mindeg || (nfp = info1.nl)
      else
         # compute the full-rank factorization G = G1*V with V coinner
         V, G1, info1 = gcrange(sys; coinner = true, atol1, atol2, rtol, offset, fast)
         info1.nrank != nrank && (@warn "Inconsistent rank evaluations: check tolerances")
         G2, info1 = glsol([G1; I], nrank; atol1, atol2, rtol, fast, poles, sdeg, mindeg)
         mindeg || (nfp = info1.nl)
         sysinv = gminreal(V'*G2; atol1, atol2, rtol)
      end
      tcond = info1.tcond;
      fnorm = info1.fnorm;
   elseif type == "1-2-3-4"
      # compute an (1,2,3,4)-inverse
      if m == nrank
         # for full column rank, any (1,2,3)-inverse is an (1,2,3,4)-inverse
         # compute the full-rank factorization G = U*G2 with U inner and
         # G2 square and invertible
         U, G2, = grange(sys; inner = true, atol1, atol2, rtol, offset, fast)
         sysinv = gminreal(G2\U'; atol1, atol2, rtol)
      elseif p == nrank
         # for full row rank, any (1,2,4)-inverse is an (1,2,3,4)-inverse
         # compute the full-rank factorization G = G2*V with V coinner and
         # G2 square and invertible
         V, G2, = gcrange(sys; coinner = true, atol1, atol2, rtol, offset, fast)
         sysinv = gminreal(V'/G2; atol1, atol2, rtol)
      else
         # compute the full-rank factorization G = U*G1 with U inner
         U, G1, info1 = grange(sys; inner = true, atol1, atol2, rtol, offset, fast)
         info1.nrank != nrank && (@warn "Inconsistent rank evaluations: check tolerances")
         # compute the full-rank factorization G1 = G2*V with V coinner
         V, G2, = gcrange(G1; coinner = true, atol1, atol2, rtol, offset, fast)
         # compute the Moore-Penrose pseudo-inverse
         sysinv = gminreal(V'*(G2\U'); atol1, atol2, rtol)
      end
   else
      error("no such computable inverse type")
   end
   info = (nrank = nrank, nfp = nfp, tcond = tcond, fnorm = fnorm)
   return sysinv, info
   #  end GINV
end


text:
abstract type AbstractJacobiWeight{T} <: Weight{T} end

axes(::AbstractJacobiWeight{T}) where T = (Inclusion(ChebyshevInterval{T}()),)

==(w::AbstractJacobiWeight, v::AbstractJacobiWeight) = w.a == v.a && w.b == v.b

broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(*), w::AbstractJacobiWeight, v::AbstractJacobiWeight) =
    JacobiWeight(w.a + v.a, w.b + v.b)

broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(sqrt), w::AbstractJacobiWeight) =
    JacobiWeight(w.a/2, w.b/2)

broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(Base.literal_pow), ::Base.RefValue{typeof(^)}, w::AbstractJacobiWeight, ::Base.RefValue{Val{k}}) where k =
    JacobiWeight(k * w.a, k * w.b)

struct JacobiWeight{T} <: AbstractJacobiWeight{T}
    a::T
    b::T
    JacobiWeight{T}(a, b) where T = new{T}(convert(T,a), convert(T,b))
end

JacobiWeight(a::V, b::T) where {T,V} = JacobiWeight{promote_type(T,V)}(a,b)
jacobiweight(a,b, d::AbstractInterval{T}) where T = JacobiWeight(a,b)[affine(d,ChebyshevInterval{T}())]

==(A::JacobiWeight, B::JacobiWeight) = A.b == B.b && A.a == B.a

function getindex(w::JacobiWeight, x::Number)
    x ∈ axes(w,1) || throw(BoundsError())
    (1-x)^w.a * (1+x)^w.b
end

sum(P::JacobiWeight) = jacobimoment(P.a, P.b)

struct LegendreWeight{T} <: AbstractJacobiWeight{T} end
LegendreWeight() = LegendreWeight{Float64}()
legendreweight(d::AbstractInterval{T}) where T = LegendreWeight{float(T)}()[affine(d,ChebyshevInterval{T}())]

function getindex(w::LegendreWeight{T}, x::Number) where T
    x ∈ axes(w,1) || throw(BoundsError())
    one(T)
end

getproperty(w::LegendreWeight{T}, ::Symbol) where T = zero(T)

sum(::LegendreWeight{T}) where T = 2one(T)

_weighted(::LegendreWeight, P) = P

broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(*), ::LegendreWeight{T}, ::LegendreWeight{V}) where {T,V} =
    LegendreWeight{promote_type(T,V)}()

broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(sqrt), w::LegendreWeight{T}) where T = w

broadcasted(::LazyQuasiArrayStyle{1}, ::typeof(Base.literal_pow), ::Base.RefValue{typeof(^)}, w::LegendreWeight, ::Base.RefValue{Val{k}}) where k = w

# support auto-basis determination

singularities(a::AbstractAffineQuasiVector) = singularities(a.x)
singularitiesbroadcast(_, L::LegendreWeight) = L # Assume we stay smooth
singularitiesbroadcast(::typeof(exp), L::LegendreWeight) = L
singularitiesbroadcast(::typeof(Base.literal_pow), ::typeof(^), L::LegendreWeight, ::Val) = L
for op in (:+, :-,:*)
    @eval begin
        singularitiesbroadcast(::typeof($op), ::LegendreWeight{T}, ::LegendreWeight{V}) where {T,V} = LegendreWeight{promote_type(T,V)}()
        singularitiesbroadcast(::typeof($op), L::LegendreWeight, ::NoSingularities) = L
        singularitiesbroadcast(::typeof($op), ::NoSingularities, L::LegendreWeight) = L
    end
end
singularitiesbroadcast(::typeof(/), ::NoSingularities, L::LegendreWeight) = L # can't find roots
singularitiesbroadcast(::typeof(*), x::AbstractJacobiWeight, y::AbstractJacobiWeight) = JacobiWeight(x.a+y.a, x.b+y.b)

_parent(::NoSingularities) = NoSingularities()
_parent(a) = parent(a)
_parentindices(a::NoSingularities, b...) = _parentindices(b...)
_parentindices(a, b...) = parentindices(a)
singularitiesbroadcast(F::Function, G::Function, V::SubQuasiArray, K) = singularitiesbroadcast(F, G, parent(V), K)[parentindices(V)...]
singularitiesbroadcast(F, V::Union{NoSingularities,SubQuasiArray}...) = singularitiesbroadcast(F, map(_parent,V)...)[_parentindices(V...)...]

singularitiesbroadcast(::typeof(*), ::LegendreWeight, b::JacobiWeight) = b
singularitiesbroadcast(::typeof(*), a::JacobiWeight, ::LegendreWeight) = a

abstract type AbstractJacobi{T} <: OrthogonalPolynomial{T} end

singularities(::AbstractJacobi{T}) where T = LegendreWeight{T}()
singularities(::Inclusion{T,<:AbstractInterval}) where T = LegendreWeight{T}()
singularities(d::Inclusion{T,<:Interval}) where T = LegendreWeight{T}()[affine(d,ChebyshevInterval{T}())]

struct Legendre{T} <: AbstractJacobi{T} end
Legendre() = Legendre{Float64}()

legendre() = Legendre()
legendre(d::AbstractInterval{T}) where T = Legendre{float(T)}()[affine(d,ChebyshevInterval{T}()), :]

==(::Legendre, ::Legendre) = true

OrthogonalPolynomial(w::LegendreWeight{T}) where {T} = Legendre{T}()
orthogonalityweight(::Legendre{T}) where T = LegendreWeight{T}()

function qr(P::Legendre)
    Q = Normalized(P)
    QuasiQR(Q, Diagonal(Q.scaling))
end

struct Jacobi{T} <: AbstractJacobi{T}
    a::T
    b::T
    Jacobi{T}(a, b) where T = new{T}(convert(T,a), convert(T,b))
end

Jacobi(a::V, b::T) where {T,V} = Jacobi{float(promote_type(T,V))}(a, b)

jacobi(a,b) = Jacobi(a,b)
jacobi(a,b, d::AbstractInterval{T}) where T = Jacobi(a,b)[affine(d,ChebyshevInterval{T}()), :]

Jacobi(P::Legendre{T}) where T = Jacobi(zero(T), zero(T))

OrthogonalPolynomial(w::JacobiWeight) = Jacobi(w.a, w.b)
orthogonalityweight(P::Jacobi) = JacobiWeight(P.a, P.b)

const WeightedJacobi{T} = WeightedBasis{T,<:JacobiWeight,<:Jacobi}

WeightedJacobi(a,b) = JacobiWeight(a,b) .* Jacobi(a,b)
WeightedJacobi{T}(a,b) where T = JacobiWeight{T}(a,b) .* Jacobi{T}(a,b)

axes(::AbstractJacobi{T}) where T = (Inclusion{T}(ChebyshevInterval{real(T)}()), OneTo(∞))
==(P::Jacobi, Q::Jacobi) = P.a == Q.a && P.b == Q.b
==(P::Legendre, Q::Jacobi) = Jacobi(P) == Q
==(P::Jacobi, Q::Legendre) = P == Jacobi(Q)
==(A::WeightedJacobi, B::WeightedJacobi) = A.args == B.args
==(A::WeightedJacobi, B::Jacobi{T}) where T = A == JacobiWeight(zero(T),zero(T)).*B
==(A::WeightedJacobi, B::Legendre) = A == Jacobi(B)
==(A::Jacobi{T}, B::WeightedJacobi) where T = JacobiWeight(zero(T),zero(T)).*A == B
==(A::Legendre, B::WeightedJacobi) = Jacobi(A) == B

###
# transforms
###

function grid(Pn::SubQuasiArray{T,2,<:AbstractJacobi,<:Tuple{<:Inclusion,<:AbstractUnitRange}}) where T
    kr,jr = parentindices(Pn)
    ChebyshevGrid{1,T}(maximum(jr))
end

function ldiv(::Legendre{V}, f::AbstractQuasiVector) where V
    T = ChebyshevT{V}()
    [cheb2leg(paddeddata(T \ f)); zeros(V,∞)]
end

function ldiv(P::Jacobi{V}, f::AbstractQuasiVector) where V
    T = ChebyshevT{V}()
    [cheb2jac(paddeddata(T \ f), P.a, P.b); zeros(V,∞)]
end

########
# Mass Matrix
#########

legendre_massmatrix(::Type{T}) where T = Diagonal(convert(T,2) ./ (2(0:∞) .+ 1))

function legendre_massmatrix(Ac, B)
    A = parent(Ac)
    P = Legendre{eltype(B)}()
    (P\A)'*legendre_massmatrix(eltype(P))*(P\B)
end

@simplify *(Ac::QuasiAdjoint{<:Any,<:Legendre}, B::Legendre) = legendre_massmatrix(Ac, B)
@simplify *(Ac::QuasiAdjoint{<:Any,<:AbstractJacobi}, B::AbstractJacobi) = legendre_massmatrix(Ac,B)

# 2^{a + b + 1} {\Gamma(n+a+1) \Gamma(n+b+1) \over (2n+a+b+1) \Gamma(n+a+b+1) n!}.

function jacobi_massmatrix(a, b)
    n = 0:∞
    Diagonal(2^(a+b+1) .* (exp.(loggamma.(n .+ (a+1)) .+ loggamma.(n .+ (b+1)) .- loggamma.(n .+ (a+b+1)) .- loggamma.(n .+ 1)) ./ (2n .+ (a+b+1))))
end

@simplify function *(wAc::QuasiAdjoint{<:Any,<:WeightedBasis{<:Any,<:AbstractJacobiWeight}}, wB::WeightedBasis{<:Any,<:AbstractJacobiWeight})
    w,A = arguments(parent(wAc))
    v,B = arguments(wB)
    A'*((w .* v) .* B)
end
@simplify function *(Ac::QuasiAdjoint{<:Any,<:AbstractJacobi}, wB::WeightedBasis{<:Any,<:AbstractJacobiWeight,<:AbstractJacobi})
    A = parent(Ac)
    w,B = arguments(wB)
    P = Jacobi(w.a, w.b)
    (P\A)' * jacobi_massmatrix(w.a, w.b) * (P \ B)
end

########
# Jacobi Matrix
########

jacobimatrix(::Legendre{T}) where T = _BandedMatrix(Vcat(((zero(T):∞)./(1:2:∞))', Zeros{T}(1,∞), ((one(T):∞)./(1:2:∞))'), ∞, 1,1)

# These return vectors A[k], B[k], C[k] are from DLMF. Cause of MikaelSlevinsky we need an extra entry in C ... for now.
function recurrencecoefficients(::Legendre{T}) where T
    n = zero(T):∞
    ((2n .+ 1) ./ (n .+ 1), Zeros{T}(∞), n ./ (n .+ 1))
end

function jacobimatrix(J::Jacobi)
    b,a = J.b,J.a
    n = 0:∞
    B = Vcat(2 / (a+b+2),  2 .* (n .+ 2) .* (n .+ (a+b+2)) ./ ((2n .+ (a+b+3)) .* (2n .+ (a+b+4))))
    A = Vcat((b-a) / (a+b+2), (b^2-a^2) ./ ((2n .+ (a+b+2)) .* (2n .+ (a+b+4))))
    C = 2 .* (n .+ a) .* (n .+ b) ./ ((2n .+ (a+b)) .* (2n .+ (a+b+1)))

    _BandedMatrix(Vcat(C',A',B'), ∞, 1,1)
end

function recurrencecoefficients(P::Jacobi)
    n = 0:∞
    ñ = 1:∞
    a,b = P.a,P.b
    A = Vcat((a+b+2)/2, (2ñ .+ (a+b+1)) .* (2ñ .+ (a+b+2)) ./ ((2*(ñ .+ 1)) .* (ñ .+ (a+b+1))))
    # n = 0 is special to avoid divide-by-zero
    B = Vcat((a-b)/2, (a^2 - b^2) * (2ñ .+ (a + b+1)) ./ ((2*(ñ .+ 1)) .* (ñ .+ (a+b+1)) .* (2ñ .+ (a+b))))
    C = ((n .+ a) .* (n .+ b) .* (2n .+ (a+b+2))) ./ (ñ .* (n .+ (a +b+1)) .* (2n .+ (a+b)))
    (A,B,C)
end

@simplify *(X::Identity, P::Legendre) = ApplyQuasiMatrix(*, P, P\(X*P))

##########
# Conversion
##########

\(A::Jacobi, B::Legendre) = A\Jacobi(B)
\(A::Legendre, B::Jacobi) = Jacobi(A)\B

function \(A::Jacobi, B::Jacobi)
    T = promote_type(eltype(A), eltype(B))
    a,b = B.a,B.b
    if A.a ≈ a && A.b ≈ b
        Eye{T}(∞)
    elseif isone(-a-b) && A.a == a && A.b == b+1
        _BandedMatrix(Vcat((((0:∞) .+ a)./((1:2:∞) .+ (a+b)))',
                            Vcat(1,((2:∞) .+ (a+b))./((3:2:∞) .+ (a+b)))'), ∞, 0,1)
    elseif isone(-a-b) && A.a == a+1 && A.b == b
        _BandedMatrix(Vcat((-((0:∞) .+ b)./((1:2:∞) .+ (a+b)))',
                            Vcat(1,((2:∞) .+ (a+b))./((3:2:∞) .+ (a+b)))'), ∞, 0,1)
    elseif A.a == a && A.b == b+1
        _BandedMatrix(Vcat((((0:∞) .+ a)./((1:2:∞) .+ (a+b)))',
                            (((1:∞) .+ (a+b))./((1:2:∞) .+ (a+b)))'), ∞, 0,1)
    elseif A.a == a+1 && A.b == b
        _BandedMatrix(Vcat((-((0:∞) .+ b)./((1:2:∞) .+ (a+b)))',
                            (((1:∞) .+ (a+b))./((1:2:∞) .+ (a+b)))'), ∞, 0,1)
    elseif A.a ≥ a+1
        J = Jacobi(a+1,b)
        (A \ J) * (J \ B)
    elseif A.b ≥ b+1
        J = Jacobi(a,b+1)
        (A \ J) * (J \ B)
    else
        error("not implemented for $A and $B")
    end
end

function \(A::Jacobi, w_B::WeightedJacobi)
    a,b = A.a,A.b
    (JacobiWeight(zero(a),zero(b)) .* A) \ w_B
end

function \(w_A::WeightedJacobi, B::Jacobi)
    a,b = B.a,B.b
    w_A \ (JacobiWeight(zero(a),zero(b)) .* B)
end

function broadcastbasis(::typeof(+), w_A::WeightedJacobi, w_B::WeightedJacobi)
    wA,A = w_A.args
    wB,B = w_B.args

    w = JacobiWeight(min(wA.a,wB.a), min(wA.b,wB.b))
    P = Jacobi(max(A.a,B.a + w.a - wB.a), max(A.b,B.b + w.b - wB.b))
    w .* P
end

function \(w_A::WeightedJacobi, w_B::WeightedJacobi)
    wA,A = w_A.args
    wB,B = w_B.args

    if wA == wB
        A \ B
    elseif B.a == A.a && B.b == A.b+1 && wB.b == wA.b+1 && wB.a == wA.a
        _BandedMatrix(Vcat((((2:2:∞) .+ 2A.b)./((2:2:∞) .+ (A.a+A.b)))', ((2:2:∞)./((2:2:∞) .+ (A.a+A.b)))'), ∞, 1,0)
    elseif B.a == A.a+1 && B.b == A.b && wB.b == wA.b && wB.a == wA.a+1
        _BandedMatrix(Vcat((((2:2:∞) .+ 2A.a)./((2:2:∞) .+ (A.a+A.b)))', -((2:2:∞)./((2:2:∞) .+ (A.a+A.b)))'), ∞, 1,0)
    elseif wB.a ≥ wA.a+1
        J = JacobiWeight(wB.a-1,wB.b) .* Jacobi(B.a-1,B.b)
        (w_A\J) * (J\w_B)
    elseif wB.b ≥ wA.b+1
        J = JacobiWeight(wB.a,wB.b-1) .* Jacobi(B.a,B.b-1)
        (w_A\J) * (J\w_B)
    else
        error("not implemented for $A and $wB")
    end
end

\(A::Legendre, wB::WeightedJacobi) = Jacobi(A) \ wB

##########
# Derivatives
##########

# Jacobi(a+1,b+1)\(D*Jacobi(a,b))
@simplify function *(D::Derivative{<:Any,<:AbstractInterval}, S::Jacobi)
    A = _BandedMatrix((((1:∞) .+ (S.a + S.b))/2)', ∞, -1,1)
    ApplyQuasiMatrix(*, Jacobi(S.a+1,S.b+1), A)
end

# Jacobi(a-1,b-1)\ (D*w*Jacobi(a,b))
@simplify function *(D::Derivative{<:Any,<:AbstractInterval}, WS::WeightedJacobi)
    w,S = WS.args
    a,b = S.a, S.b
    if w.a == 0 && w.b == 0
        D*S
    elseif iszero(w.a) && w.b == b #L_6
        A = _BandedMatrix((b:∞)', ∞, 0,0)
        ApplyQuasiMatrix(*, JacobiWeight(w.a,b-1) .* Jacobi(a+1,b-1), A)
    elseif iszero(w.b) && w.a == a #L_6^t
        A = _BandedMatrix(-(a:∞)', ∞, 0,0)
        ApplyQuasiMatrix(*, JacobiWeight(a-1,w.b) .* Jacobi(a-1,b+1), A)
    elseif w.a == a && w.b == b # L_1^t
        A = _BandedMatrix((-2*(1:∞))', ∞, 1,-1)
        ApplyQuasiMatrix(*, JacobiWeight(a-1,b-1) .* Jacobi(a-1, b-1), A)
    elseif iszero(w.a)
        W = (JacobiWeight(w.a, b-1) .* Jacobi(a+1, b-1)) \ (D * (JacobiWeight(w.a,b) .* S))
        J = Jacobi(a+1,b) # range Jacobi
        C1 = J \ Jacobi(a+1, b-1)
        C2 = J \ Jacobi(a,b)
        ApplyQuasiMatrix(*, JacobiWeight(w.a,w.b-1) .* J, (w.b-b) * C2 + C1 * W)
    elseif iszero(w.b)
        W = (JacobiWeight(a-1, w.b) .* Jacobi(a-1, b+1)) \ (D * (JacobiWeight(a,w.b) .* S))
        J = Jacobi(a,b+1) # range Jacobi
        C1 = J \ Jacobi(a-1, b+1)
        C2 = J \ Jacobi(a,b)
        ApplyQuasiMatrix(*, JacobiWeight(w.a-1,w.b) .* J, -(w.a-a) * C2 + C1 * W)
    else
        error("Not implemented")
    end
end

function \(L::Legendre, WS::WeightedBasis{Bool,JacobiWeight{Bool},Jacobi{Bool}})
    w,S = WS.args
    if w.b && w.a
        @assert S.b && S.a
        _BandedMatrix(Vcat(((2:2:∞)./(3:2:∞))', Zeros(1,∞), (-(2:2:∞)./(3:2:∞))'), ∞, 2,0)
    elseif w.b && !w.a
        @assert S.b && !S.a
        _BandedMatrix(Ones{eltype(L)}(2,∞), ∞, 1,0)
    elseif !w.b && w.a
        @assert !S.b && S.a
        _BandedMatrix(Vcat(Ones{eltype(L)}(1,∞),-Ones{eltype(L)}(1,∞)), ∞, 1,0)
    else
        error("Not implemented")
    end
end

@simplify function *(St::QuasiAdjoint{Bool,Jacobi{Bool}}, WS::WeightedBasis{Int,JacobiWeight{Int},Jacobi{Bool}})
    w = parent(W)
    (w.b == 2 && S.b && w.a == 2 && S.a && parent(St) == S) || throw(ArgumentError())
    W_sqrt = Diagonal(JacobiWeight(true,true))
    L = Legendre()
    A = PInv(L)*W_sqrt*S
    A'*(L'L)*A
end

###
# Splines
###

function \(A::Legendre, B::HeavisideSpline)
    @assert B.points == -1:2:1
    Vcat(1, Zeros(∞,1))
end

###
# sum
###

function _sum(P::Legendre{T}, dims) where T
    @assert dims == 1
    Hcat(convert(T, 2), Zeros{T}(1,∞))
end

_sum(p::SubQuasiArray{T,1,Legendre{T},<:Tuple{Inclusion,Int}}, ::Colon) where T = parentindices(p)[2] == 1 ? convert(T, 2) : zero(T)


text:
function [X]=Cart(DG,P1,P2,TypeE,Param)

if strcmp(TypeE,'Y')
  ksi=DG.xwY;
  n=DG.OrdPolyY+1;
  X=zeros(n,2);
  X(:,1)=.5*(P1(1)*(1-ksi)+...
    P2(1)*(1+ksi));
  X(:,2)=.5*(P1(2)*(1-ksi)+...
    P2(2)*(1+ksi));
else
  ksi=DG.xwX;
  n=DG.OrdPolyX+1;
  X=zeros(n,2);
  X(:,1)=.5*(P1(1)*(1-ksi)+...
    P2(1)*(1+ksi));
  X(:,2)=.5*(P1(2)*(1-ksi)+...
    P2(2)*(1+ksi));
  X(:,1)=.5*(P1(1)*(1-ksi)+...
    P2(1)*(1+ksi));
  hSX1=hS(P1(1),Param);
  hSX2=hS(P2(1),Param);
  if abs(hSX2-hSX1)>0
    for i=1:n
      X(i,2)=P1(2)+(P2(2)-P1(2))/(hSX2-hSX1)*(hS(X(i,1),Param)-hSX1);
    end
  else
    X(:,2)=.5*(P1(2)*(1-ksi)+...
      P2(2)*(1+ksi));
  end
end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule USRPHardwareDriver_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("USRPHardwareDriver")
JLLWrappers.@generate_main_file("USRPHardwareDriver", UUID("e339f254-40e1-517e-839f-dc67f1dbd84e"))
end  # module USRPHardwareDriver_jll


text:
"""
Package: Forecast

    boxcox(x::Vector)
    boxcox(x::Vector, λ::Float64)
    boxcox(x::Vector, λ::Vector)

Compute a Box-Cox power transformation given λ or (λ1,λ2) for data containing negative values, or compute an optimal power transformation if no λ or (λ1,λ2) is provided.

```math
x(\\lambda) =
\\begin{cases}
 \\dfrac{x_i^\\lambda - 1}{\\lambda} & \\text{if } \\lambda \\neq 0 \\\\
 \\ln x_i & \\text{if } \\lambda = 0
\\end{cases}
```
for negative values

```math
x(\\boldsymbol{\\lambda}) =
\\begin{cases}
 \\dfrac{(x_i + \\lambda_2)^{\\lambda_1} - 1}{\\lambda_1} & \\text{if } \\lambda_1 \\neq 0 \\\\
 \\ln (x_i + \\lambda_2) & \\text{if } \\lambda_1 = 0
\\end{cases}
```

# Arguments
- `x`: Vector to be transformed.
- `λ`: Exponent/s for the tranformation

# Returns
A vector with a boxcox tarnsformation for `x` or a Dict with :x boxcox tranformed and the optimal :λ

# Reference
Box, G. E. P. and Cox, D. R. (1964). An analysis of transformations, Journal of the Royal Statistical Society, Series B, 26, 211-252. A

# Examples
```julia-repl
julia> x = rand(100)
julia> bc = boxcox(x)
julia> iboxcox(bc[:x],bc[:λ]) ≈ x

julia> x = rand(100) .- 0.5
julia> bc = boxcox(x)
julia> iboxcox(bc[:x],bc[:λ]) ≈ x
```
"""
function boxcox(x::AbstractVector{T}) where T<:Real

    negative = any(x -> x < 0, x)

    if negative
        λ2 = -(minimum(x) - T(1))
        x = x .+ λ2
    end

    f(λ) = -pvalue(JarqueBeraTest(boxcox(x,λ)))
    λ = Optim.minimizer(optimize(f,-3.0,3.0,Brent()))

    return Dict(:x => boxcox(x,λ), :λ => negative ? (λ,λ2) : λ)

end

boxcox(x::AbstractVector{<:Real}, λ::Real)  = λ == 0.0 ? log.(x) : (x .^ λ .- 1.0) ./ λ
boxcox(x::AbstractVector{<:Real},
       λ::Tuple{Real,Real}) =
           λ[1] == 0.0 ? log.(x .+ λ[2]) : ( (x .+ λ[2]) .^ λ[1] .- 1.0) ./ λ[1]

"""
Package: Forecast

    iboxcox(x::Vector, λ::Float64)
    iboxcox(x::Vector, λ::Vector)

Compute the inverse transformation of a Box-Cox power transformation given λ.

# Arguments
- `x`: Vector with a boxcox tranformation to be inverted.
- `λ`: Exponent for the inverse tranformation.

# Returns
A vector with witht the inverse transformation of x given λ.

# Reference
Box, G. E. P. and Cox, D. R. (1964). An analysis of transformations, Journal of the Royal Statistical Society, Series B, 26, 211-252. A

# Examples
```julia-repl
julia> x = rand(100)
julia> bc = boxcox(x)
julia> iboxcox(bc[:x],bc[:λ]) ≈ x

julia> x = rand(100) .- 0.5
julia> bc = boxcox(x)
julia> iboxcox(bc[:x],bc[:λ]) ≈ x
```
"""
iboxcox(x::AbstractVector{<:Real}, λ::Real) = λ == 0.0 ? exp.(x) : (1.0 .+ λ .* x ) .^ (1/λ)
iboxcox(x::AbstractVector{<:Real},
        λ::Tuple{Real,Real}) =
    λ[1] == 0.0 ? exp.(x) .- λ[2] : (1.0 .+ λ[1] .* x ) .^ (1/λ[1]) .- λ[2]


text:
module QCurlGradMonomialBasesTests

using Test
using Gridap.TensorValues
using Gridap.Fields
using Gridap.Polynomials

xi = Point(2,3)
np = 5
x = fill(xi,np)

order = 0
D = 2
T = Float64
V = VectorValue{D,T}
G = gradient_type(V,xi)
b = QCurlGradMonomialBasis{D}(T,order)

@test num_terms(b) == 4
@test get_order(b) == 0

xi = Point(2,3,5)
np = 5
x = fill(xi,np)

order = 0
D = 3
T = Float64
V = VectorValue{D,T}
G = gradient_type(V,xi)
b = QCurlGradMonomialBasis{D}(T,order)

v = V[
  (1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 0.0, 1.0),
  (2.0, 0.0, 0.0), (0.0, 3.0, 0.0), (0.0, 0.0, 5.0)]

g = G[
  (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
  (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
  (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
  (1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
  (0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0),
  (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0)]

bx = repeat(permutedims(v),np)
∇bx = repeat(permutedims(g),np)
test_field(b,x,bx,grad=∇bx)

xi = Point(2,3)
np = 5
x = fill(xi,np)

order = 1
D = 2
T = Float64
V = VectorValue{D,T}
G = gradient_type(V,xi)
b = QCurlGradMonomialBasis{D}(T,order)

v = V[
  (1.0, 0.0), (0.0, 1.0), (2.0, 0.0), (0.0, 3.0),
  (4.0, 0.0), (0.0, 9.0), (3.0, 0.0), (0.0, 2.0),
  (6.0, 0.0), (0.0, 6.0), (12.0, 0.0), (0.0, 18.0)]

g = G[
  (0.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 0.0),
  (1.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 1.0),
  (4.0, 0.0, 0.0, 0.0), (0.0, 0.0, 0.0, 6.0),
  (0.0, 1.0, 0.0, 0.0), (0.0, 0.0, 1.0, 0.0),
  (3.0, 2.0, 0.0, 0.0), (0.0, 0.0, 3.0, 2.0),
  (12.0, 4.0, 0.0, 0.0),(0.0, 0.0, 9.0, 12.0)]

bx = repeat(permutedims(v),np)
∇bx = repeat(permutedims(g),np)
test_field(b,x,bx,grad=∇bx)

end # module


text:
export randomscopes!

"""
    randomscopes!

Draw scopes from prior distribution.
"""
function randomscopes!(spn::AbstractRegionGraphNode, x::AbstractMatrix{<:Real})
    (N,D) = size(x)

    _randomscopes!(spn, collect(1:D))
end

function _randomscopes!(n::RegionGraphNode, dims::AbstractVector{Int})
    _randomscopes!.(children(n), Ref(dims))
    setscope!(n, dims)
end

function _randomscopes!(n::PartitionGraphNode, dims::AbstractVector{Int})
    K = length(children(n))
    α = postparams(n.prior, DirichletSufficientStats(K))
    prior = α / sum(α)

    z = map(d -> randomindex(prior), 1:length(dims))
    for (i,child) in enumerate(children(n))
        _randomscopes!(child, dims[z .== i])
    end
    setscope!(n, dims)
end

function _randomscopes!(n::FactorizedDistributionGraphNode, dims::AbstractVector{Int})
    setscope!(n, dims)
end

function _randomscopes!(n::FactorizedMixtureGraphNode, dims::AbstractVector{Int})
    setscope!(n, dims)
end

max_stars_repo_name: crinders/BayesianSumProductNetworks

max_issues_repo_name: crinders/BayesianSumProductNetworks

max_forks_repo_name: crinders/BayesianSumProductNetworks


text:
# The API server serving functions from srvrfn.jl.
# Runs server in blocking mode when invoked directly with "--runsrvr" argument.
# Call `run_srvr` to start server otherwise.
using JuliaWebAPI
using HttpCommon
using Logging
using Compat
using ZMQ

include("srvrfn.jl")

const SRVR_ADDR = "tcp://127.0.0.1:9999"
const JSON_RESP_HDRS = Dict{String,String}("Content-Type" => "application/json; charset=utf-8")
const BINARY_RESP_HDRS = Dict{String,String}("Content-Type" => "application/octet-stream")

function run_srvr(fmt, tport, async=false)
    Logging.configure(level=INFO, filename="apisrvr_test.log")
    Logging.info("queue is at $SRVR_ADDR")

    api = APIResponder(tport, fmt)
    Logging.info("responding with: $api")

    register(api, testfn1; resp_json=true, resp_headers=JSON_RESP_HDRS)
    register(api, testfn2)
    register(api, testbinary; resp_headers=BINARY_RESP_HDRS)
    register(api, testArray)
    register(api, testFile; resp_json=true, resp_headers=JSON_RESP_HDRS)

    process(api; async=async)
end

function test_preproc(req::Request, res::Response)
    for (n,v) in req.headers
        if lowercase(n) == "juliawebapi"
            res.status = parse(Int, v)
            return false
        end
    end
    JuliaWebAPI.default_preproc(req, res)
end

function run_httprpcsrvr(fmt, tport, async=false)
    run_srvr(fmt, tport, true)
    apiclnt = APIInvoker(ZMQTransport(SRVR_ADDR, REQ, false), fmt)
    if async
        @async run_http(apiclnt, 8888, test_preproc)
    else
        run_http(apiclnt, 8888, test_preproc)
    end
end

# run in blocking mode if invoked with run flag
!isempty(ARGS) && (ARGS[1] == "--runsrvr") && run_srvr(JuliaWebAPI.JSONMsgFormat(), JuliaWebAPI.ZMQTransport(SRVR_ADDR, ZMQ.REP, true), false)

# run http rpc server if invoked with flag
!isempty(ARGS) && (ARGS[1] == "--runhttprpcsrvr") && run_httprpcsrvr(JuliaWebAPI.JSONMsgFormat(), JuliaWebAPI.ZMQTransport(SRVR_ADDR, ZMQ.REP, true), false)

max_stars_repo_name: tanmaykm/JuliaWebAPIPlugins

max_issues_repo_name: tanmaykm/JuliaWebAPIPlugins

max_forks_repo_name: tanmaykm/JuliaWebAPIPlugins


text:
using Compat
using Compat.Test
using AccurateArithmetic

@test true