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

using Markdown
using InteractiveUtils

# ╔═╡ afdc0ce2-a8c5-11ec-0907-1b7d0d9db74a
begin
	using Plots
	Plots.theme(:ggplot2)
	gr()
end

# ╔═╡ b4e08ac2-91c3-44c0-8a00-4b3938d5f0dd
begin
	# Math & Stats packages
	using Distributions, Random, LinearAlgebra
	rng = Random.MersenneTwister(1234)
end

# ╔═╡ 623b22e7-5c4e-4366-831c-5bc987eb8895
md"""
Version Information:
* Julia v.1.7.2
* Pluto v.0.18.4

"""

# ╔═╡ 0a5f2234-da1f-44f0-9e68-b2da7d9c4527
md"""
# NEUROSCI 4L09

## Decoding the Brain: How Bayesian Categorization Models Handle Category and Sensory Variability

**Abstract**

Brains must make inferences about the world to survive. They encode external stimuli as sensory data, then decode the sensory data to make perceptual judgements and take actions. A common perceptual judgement is categorization, the identification of the type of object stimulating the senses. Categorization is challenging, because 1) objects within the same category typically vary from one another along several physical dimensions, causing overlap with the physical features of objects from other categories, and 2) sensory responses are stochastic, such that interactions with a given object evoke variable sensory data. Both these sources of variability reduce the fidelity of the encoding process and put a greater onus on the brain to flexibly decode variable information to make useful inferences.

An optimal decoder would account for both sources of variability. However, if the variation among objects in a category or the variation of sensory data were sufficiently low, it might be reasonable for the brain to simply ignore either, or both, of these factors when making inferences about the world.

Using computer simulations, three Bayesian decoders were tested on a tactile binary categorization task: 1) an optimal decoder that accounted for category and sensory variability, 2) a decoder that accounted only for category variability, and 3) a decoder that accounted only for sensory variability. The performance results from these decoders will provide a useful reference to compare with participant data to determine which method the brain likely uses to make meaningful inferences about the world.
"""

# ╔═╡ e311bac7-d259-436b-afdb-4c32866c73b5
md"""
# Tactile Objects

--- Add picture of objecte here ---

Consider an object with two features that can be detected by sense of touch:
* Sides: 6, 7, 8, 9, or 10 sides
* Dots: 4, 5, 6, 7, or 8 mm dot spacing
"""

# ╔═╡ 0243fdb5-8109-411e-b8b0-1892e6f703d9
begin
	feat_ranges = (sides=6:10, dots=4:8)
	feat_ranges
end

# ╔═╡ 67804342-2479-453b-b551-cc9508e6cbff
begin
	"""
		TactileObject(label::Symbol, sides::Real, dots::Real)
	Tactile Object called label with two features, sides and dots

		Tactileobject(label::Symbol, feat_dist::Distribution, feat_ranges)
	Tactile Object called label constructed from bounded feature distributions

		TactileObject(category::Category, feat_ranges)
	Tactile Object in category with bounded features
	"""
	struct TactileObject
		category::Symbol
		sides::Real
		dots::Real
		# Alternative methods to construct a TactileObject
		TactileObject(label::Symbol, sides::Real, dots::Real) = new(label, sides, dots)
		TactileObject(label::Symbol, feat_dist::Distribution) = new(label, rand(feat_dist)...)
		function TactileObject(label::Symbol, feat_dist::Distribution, feat_ranges)
			sides, dots = rand(feat_dist)
			sides = clamp(sides, extrema(feat_ranges.sides)...)
			dots = clamp(dots, extrema(feat_ranges.dots)...)
			return new(label, sides, dots)
		end
		TactileObject(category::Category, feat_ranges) = TactileObject(category.name, category.dist, feat_ranges)
	end

	# Pretty print TactileObject
	Base.show(io::IO, tactobj::TactileObject) = print(io, "$(tactobj.category) with Sides: $(tactobj.sides), Dots: $(tactobj.dots)")

	"""
		get_feats(tactobj::TactileObject)

	Returns features (sides & dots) of Tactile Object as a tuple
	"""
	function get_feats(tactobj::TactileObject)
		return tactobj.sides, tactobj.dots
	end
end

# ╔═╡ b0ea37ee-8624-48ea-9be4-0b5be91a57cf
Markdown.parse("""
Example construction of unrealistic tactile objects:

	$(TactileObject(:some_object, -5, 6.3))

	$(TactileObject(:another_object, 3.5, 0.0))
""")

# ╔═╡ 6493908b-4773-4782-8616-158ad512ffcd
md"""
# Categories

There are two categories of tactile objects: 1) Elyk and 2) Noek

While they each have the same kind of features, sides and dots, the distribution of feature values differentiates the categories.
"""

# ╔═╡ 4a531fe7-d8be-42fa-b356-24c36bb2d35d
struct Category
	name::Symbol
	dist::Distribution
end

# ╔═╡ cbb058d9-a3bb-44d3-8dbc-e782dd27f172
Markdown.parse("""
## Elyks

Elyks have this distribution of feature values:
* **Mean**: ``μ_\\text{sides}`` = 7 sides, ``μ_\\text{dots}`` = 5 mm dot spacing
* **Stdev**: ``σ_\\text{sides}`` = 1.5 sides, ``σ_\\text{dots}`` = 1.5 mm dot spacing

We mathematically represent the Elyk category as a [multivariate normal distribution](https://en.wikipedia.org/wiki/Multivariate_normal_distribution):

```math
\\text{Elyk} \\sim N(\\vec{μ}_E, Σ_E)
```

where ``\\vec{μ}_{\\text{Elyk}} = [7, 5]^T`` (column vector of means)
and ``Σ_{\\text{Elyk}} = 1.5^2 I_2`` (covariance matrix)

> 	Caveats:
> 		The category model is continous, and features values exist for all real numbers
>
> 	However, in real life:
> 	* The distribution of feature values is discrete
> 	* The feature values have upper and lower bounds (6-10 sides, 4-8 mm dot spacing)
""")

# ╔═╡ a5d1b3b8-5fe5-4791-bcfa-bfc74d1db345
begin
	# Mean Sides, Dots
	elyk_μ = [7.0, 5.0]
	# Stdev Sides, Dots
	elyk_Σ = (1.5^2)I
	# Elyk Object Distribution
	elyk_dist = MvNormal(elyk_μ, elyk_Σ)

	elyk = Category(:elyk, elyk_dist)
end

# ╔═╡ 246ad164-b009-4666-8667-7d760ae09acc
begin
	# Don't set res smaller than 0.01, otherwise rendering takes too long
	res = 0.01
	sides_ax = range(extrema(feat_ranges.sides)..., step=res)
	dots_ax = range(extrema(feat_ranges.dots)..., step=res)

	elyk_pdf = [pdf(elyk_dist, [side, dot]) for side in sides_ax, dot in dots_ax]

	Plots.heatmap(sides_ax, dots_ax, elyk_pdf,
    	xlabel="Sides",
		ylabel="Dots (mm)",
    	title="Elyk Distribution",
	)
end

# ╔═╡ 70ea935b-cbb7-48c4-810d-f70e3087c4f2
# Spawn a random elyk
TactileObject(elyk, feat_ranges)

# ╔═╡ f5e57a2d-b3b9-409a-8ced-56de60f85528
Markdown.parse("""
## Noeks

Noeks have this distribution of feature values:
* **Mean**: ``μ_\\text{sides}`` = 9 sides, ``μ_\\text{dots}`` = 7 mm dot spacing
* **Stdev**: ``σ_\\text{sides}`` = 1.5 sides, ``σ_\\text{dots}`` = 1.5 mm dot spacing

We mathematically represent the Noek category as:

```math
\\text{Noek} \\sim N(\\vec{μ}_N, Σ_N)
```

where ``\\vec{μ}_{\\text{Noek}} = [9, 7]^T``
and ``Σ_{\\text{Noek}} = 1.5^2 I_2``

""")

# ╔═╡ 929bea75-2ffa-4e48-b27c-4172eb8c1e51
begin
	# Sides, Dots mean
	noek_μ = [9.0, 7.0]
	# Sides, Dots stdev
	noek_Σ = (1.5^2)I
	# Noek Object Distribution
	noek_dist = MvNormal(noek_μ, noek_Σ)

	noek = Category(:noek, noek_dist)
end

# ╔═╡ 8fcb93f6-1544-4b84-984b-aa09c1630a07
begin
	noek_pdf = [pdf(noek_dist, [side, dot]) for side in sides_ax, dot in dots_ax]

	Plots.heatmap(sides_ax, dots_ax, noek_pdf,
    	xlabel="Sides",
		ylabel="Dots (mm)",
    	title="Noek Distribution",
	)
end

# ╔═╡ 2d505656-00ca-455c-8710-dd2bf1cd6637
# Spawn a random noek
TactileObject(noek, feat_ranges)

# ╔═╡ 91b45c05-59bf-42fb-8190-aa7918d23fba
Markdown.parse("""
# Encoder

The encoder transforms the true stimulus into a sensory representation of the tactile objects

When the encoder sense an object, it can't perfectly sense (i.e. measure) the object. We'll assume this "noisy" measurement is normally distributed.

```math
\\hat{x} = \\text{Encoder}(\\vec{x}) = N(\\vec{x}, \\Sigma_\\text{encoder})
```

* ``\\vec{x}`` is the actual tactile object
* ``\\hat{x}`` is the sensation of the object created by the encoder
* ``\\Sigma_\\text{encoder}`` is the internal noise (variance) of the encoder

""")

# ╔═╡ 20a41be0-9b74-4e74-9e47-89c2e43b53a5
function encode(stimulus::TactileObject, Σ_encoder)
	sensation_dist = MvNormal([stimulus.sides, stimulus.dots], Σ_encoder)
	sides, dots = rand(sensation_dist)
	return TactileObject(:object, sides, dots)
end

# ╔═╡ 5853bdb4-ecba-4504-bb18-5cb9bafca0f4
begin
	function encoding_example()
		Σ_encoder = 1.2I
		an_elyk = TactileObject(:elyk, mean(elyk.dist)...)
		encoded_stimulus = encode(an_elyk, Σ_encoder)
		return (stimulus=an_elyk,
				noise=Σ_encoder,
				encoding=encoded_stimulus)
	end

	encoding_ex = encoding_example()

	Markdown.parse("""
	### Encoding Example

	Original stimulus:
		```$(encoding_ex.stimulus)```

	Encoder's sensation of stimulus:
		```$(encoding_ex.encoding)```

	After re-running this cell, notice that the encoder cannot consistently sense the elyk object.

	""")
end

# ╔═╡ a898528f-c913-4b8c-b814-9f0b98140979
begin
	function get_encodings(object, samples, Σ)
		return [get_feats(encode(object, Σ)) for _=1:samples]
	end

	function plot_encodings(feats, encodings, σ)
		p = Plots.scatter(encodings, markersize=2,
			title="σ=$(σ)",
			label="encoding",
			xlims=(0, 15),
			ylims=(0, 15)
		)
		Plots.scatter!(p, [feats], markersize=4, label="stimulus")
		return p
	end

	function plot_encoding_samples(subplots)
		return plot!(subplots...,
			plot_title="Encoder Sampling Distributions",
			plot_titlefontsize=12,
			xlabel="Sides", ylabel="Dots",
			legend_font_pointsize=6,
			titlefontsize=10,
			xguidefontsize=10,
			yguidefontsize=10,
			size=(500, 350)
		)
	end

	function encoder_sampler_example()
		encoder_σs = 0.5:0.5:2.0
		encoder_Σs = [(σ^2)I for σ in encoder_σs]
		samples = 1000
		object = TactileObject(elyk, feat_ranges)
		feats = get_feats(object)
		all_encodings = [get_encodings(object, samples, Σ) for Σ in encoder_Σs]
		subplots = [plot_encodings(feats, encodings, σ) for (σ, encodings) in zip(encoder_σs, all_encodings)]
		return plot_encoding_samples(subplots)
	end

	encoder_sampler_example()
end

# ╔═╡ fa662a02-4d52-4487-869a-05cb1394662a
Markdown.parse("""
### How much encoding noise is realistic?

Based on experimental data, noise values (``\\sigma``) of 1.1 to 1.6 are realistic for these tactile objects.

See the [2ifc.jl notebook]() for more details.
""")

# ╔═╡ bb6adffa-4e8a-481d-9f56-6e92192a0f1f
Markdown.parse("""
## Decoding Sensory Variability

The sensation decoder infers the true stimulus based on samples provided by the encoder.

> We'll assume that the decoder knows ``Σ_\\text{encoder}``

To deal with variability in encoder, the decoder instructs the encoder to take multiple samples of the tactile object (``\\vec{x}``). Based on those multiple samples, the decoder infers what the original stimulus most likely must have been.

Mathematically, the decoder takes in as INPUT:

```math
\\vec{X} = [\\hat{x}_1 \\hat{x}_2 \\hat{x}_3 ... \\hat{x}_N]
```

``\\vec{X}`` is a dataset of ``N`` samples of the stimulus. Each ``\\vec{x}_i`` comes from the encoder; they represent a sensation/measurement of the tactile object. We implement ``\\hat{x}_i`` as a 2-element vector: ``[\\text{ sides } \\text{ dots }]^T``

The decoder constructs a vector of hypothesized true feature values, ``\\vec{μ}``, with uniform priors.

The decoder then calculates the likelihood of the dataset of encodings for each hypothesized feature value, and uses Bayes theorem to determine the probability of each hypothesis being true:

```math
p_i = P(μ_i | \\vec{X}) = \\frac{P(\\vec{X} | μ_i)}{\\sum_i P(\\vec{X} | μ_i)}
```

For the OUTPUT, the decoder infers the true stimulus is the hypothesis with the highest probability of explaining the dataset:

```math
\\vec{x}_\\text{estimate} = \\vec{μ}[ \\text{argmax} (p_1, ..., p_i) ]
```

In summary:

```math
\\vec{x}_\\text{estimate} = \\text{sensation decoder} ( \\vec{X} )
```

""")

# ╔═╡ 5f38087b-17b0-4ef2-93b8-92c63f71aeeb
"""
	decode_sense(sensations, Σ_encoder, feat_ranges)

Given set of stochastic sensations (i.e. encodings) of a stimulus, infer the most likely true features of the stimulus
"""
function decode_sense(sensations::Vector{T}, Σ_encoder, feat_ranges) where {T<:TactileObject}
	# Hypothesize several possible true stimulus values
	step = 0.1
	sides_hyp = range(extrema(feat_ranges.sides)..., step=step)
	dots_hyp = range(extrema(feat_ranges.dots)..., step=step)
	μ_hyps = [[side, dot] for side in sides_hyp, dot in dots_hyp]

	# Assume uniform prior for μ_hyps
	# prior = 1 / length(μ_hyps)

	# Calculate likelihood of sensed stimulus given hypotheses
	# 	Assume decoder knows the variability of the encoder
	data_loglikelihood(μ, stimuli) = sum(logpdf(Distributions.MvNormal(μ, Σ_encoder), collect(get_feats(sensation))) for sensation in sensations)

	loglikelihoods = [data_loglikelihood(μ, sensations) for μ in μ_hyps]
	likelihoods = exp.(loglikelihoods .- maximum(loglikelihoods))

	# Calculate most probable true stimulus values
	posteriors = likelihoods ./ sum(likelihoods)
	estimated_features = μ_hyps[argmax(posteriors)]
	return TactileObject(:object, estimated_features...)
end

# ╔═╡ 4c0922fb-5a7d-457b-877a-da406a071c8c
begin
	function sense_decoding_example()
		an_elyk = TactileObject(elyk, feat_ranges)
		Σ_encoder = (1.4^2)I
		# Encode stimulus
		num_samples = 10
		encodings = [encode(an_elyk, Σ_encoder) for _=1:num_samples]
		# Decode sensory encoding
		sense_decoding = decode_sense(encodings, Σ_encoder, feat_ranges)
		return (stimulus=an_elyk,
				samples=num_samples,
				decoding=sense_decoding)
	end

	sensory_decoding_ex = sense_decoding_example()

	Markdown.parse("""
	### Sensory Decoder Example

	Example stimulus:
	```$(sensory_decoding_ex.stimulus)```

	Number of encoding samples:
	```$(sensory_decoding_ex.samples)```

	Decoded stimulus:
	```$(sensory_decoding_ex.decoding)```

	Notes:
	> Increasing `num_samples` improves the accuracy of the decoder (more data => better inference)
	""")
end

# ╔═╡ 3c34fb99-7e4b-4387-b3b8-b4122f60b3f2
Markdown.parse("""
## Decoding Category Variability

The category decoder infers the true category based on the features provided by the encoder and information about the category distribution.

> We'll assume the decoder has already learned the mean feature values (``μ_\\text{Elyk}`` and ``μ_\\text{Noek}``), and the variances (``Σ_\\text{Elyk}`` and ``Σ_\\text{Noek}``) of the categories.

As INPUT, this decoder receives a sensory estimate of the stimulus (``\\vec{x}_\\text{estimate}``).

Using Bayes theorem, the decoder determines the probability of the stimulus belonging to either category:

```math
\\vec{y} = P(\\text{category}_i | \\vec{x}_\\text{estimate}) = \\frac{P(\\vec{x}_\\text{estimate} | \\text{category}_i )}{\\sum_i P(\\vec{x}_\\text{estimate} | \\text{category}_i )}
```

``P(\\vec{x}_\\text{estimate} | \\text{category}_i ) = f(\\vec{x}_\\text{estimate}; \\vec{μ}_{\\text{category}_i}, Σ_{\\text{category}_i})``
where ``f`` is the [multivariate normal pdf](https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Density_function)

The category decoder OUTPUTs the probabilities that the stimulus belongs to each category, for example:

```math
\\vec{y} = [\\text{(elyk)} 0.71 \\text{ (noek)} 0.29]
```

In summary:

```math
\\vec{y} = \\text{category decoder} ( \\vec{x}_\\text{estimate} )
```

""")

# ╔═╡ e9c2b10a-c75e-438c-9860-591b5c64a55a
"""
	decode_category(stimulus, Σ_category)

Given a stimulus from a stochastic category distribution, infer the most likely true category of the stimulus
"""
function decode_category(stimulus::TactileObject, Σ_category)
	stimulus_feats = collect(get_feats(stimulus))
	# Calculate likelihoods
	loglikelihoods = [
		logpdf(Distributions.MvNormal(mean(elyk_dist), Σ_category), stimulus_feats),
		logpdf(Distributions.MvNormal(mean(noek_dist), Σ_category), stimulus_feats)
	]
	category_likelihoods = exp.(loglikelihoods .- maximum(loglikelihoods))
	# Calculate posteriors
	category_posteriors = category_likelihoods ./ sum(category_likelihoods)
	p_elyk, p_noek = category_posteriors
	return (elyk = round(p_elyk; digits=3), noek = round(p_noek; digits=3))
end

# ╔═╡ ce3da79d-4555-4693-b8c5-30cff21d8306
begin
	function example_category_decoding()
		Σ_category = cov(elyk_dist)
		# Create example stimulus & decoded stimulus
		sense_example = sense_decoding_example()
		sensed_elyk = sense_example.decoding
		# Categorize decoded stimulus
		category_decoding = decode_category(sensed_elyk, Σ_category)
		return (true_stimulus=sense_example.stimulus,
			sense_decoding=sensed_elyk,
			category_decoding=category_decoding)
	end

	category_decoding_ex = example_category_decoding()

	Markdown.parse("""
	### Category Decoder Example

	True Stimulus:
	```$(category_decoding_ex.true_stimulus)```

	Decoded stimulus:
	```$(category_decoding_ex.sense_decoding)```

	Decoded category:
	```$(category_decoding_ex.category_decoding)```

	""")
end

# ╔═╡ dd5aed91-b294-4962-ad80-c030d85181aa
md"""
# Bayesian Decoder Models
"""

# ╔═╡ 16ee34e1-540c-422f-90d2-8d173be877c4
struct ModelParams
	feat_ranges::NamedTuple
	num_encodings::Integer
	Σ_encoder
	Σ_category
end

# ╔═╡ 84cbe1fa-696c-456a-8803-9d12494df18b
begin
	function model_debug(stimulus::TactileObject, decoded_stim, decoded_cat)
		@show stimulus
		@show decoded_stim
		@show decoded_cat
		return
	end
end

# ╔═╡ 5d3a5b68-4345-4aa5-aa06-1022043aed85
md"""
## Model 1: Inference With Sensory Variability Only
"""

# ╔═╡ c1ab6a91-721b-4356-95b8-4d4d43500ad3
"""
	model1(stimulus::TactileObject, params::ModelParams)

Decoder model that only takes into account sensory variability.

When making inferences, the model assumes Σ_category is sufficiently small to ignore.
"""
function model1(stimulus::TactileObject, params::ModelParams; debug=false)
	# Sense stimulus
	sensations = [encode(stimulus, params.Σ_encoder) for _=1:params.num_encodings]
	# Incorporate sensory variability
	estimated_stimulus = decode_sense(sensations, params.Σ_encoder, params.feat_ranges)
	# Assume sufficiently small category varibility
	assumed_Σ_category = (0.01^2)I
	estimated_category = decode_category(estimated_stimulus, assumed_Σ_category)
	debug && model_debug(stimulus, estimated_stimulus, estimated_category)
	return estimated_category
end

# ╔═╡ cda8e5c4-8789-473e-bfe2-3c084ff96ccb
md"""
## Model 2: Inference With Category Variability Only

"""

# ╔═╡ ed2894cf-cfc8-4bb1-8441-8459d101087a
"""
	model2(stimulus::TactileObject, params::ModelParams)

Decoder model that only takes into account category variability.

When making inferences, the model assumes Σ_encoder is sufficiently small to ignore.
"""
function model2(stimulus::TactileObject, params::ModelParams; debug=false)
	# Sense stimulus
	sensations = [encode(stimulus, params.Σ_encoder) for _=1:params.num_encodings]
	# Assume sufficiently small sensory varibility
	assumed_Σ_encoder = 0.001I
	estimated_stimulus = decode_sense(sensations, assumed_Σ_encoder, params.feat_ranges)
	# Incorporate category variability
	estimated_category = decode_category(estimated_stimulus, params.Σ_category)
	debug && model_debug(stimulus, estimated_stimulus, estimated_category)
	return estimated_category
end

# ╔═╡ ce6fdfe8-cd3e-469c-a384-25b564625d2d
md"""
## Model 3: Inference With Category and Sensory Variability
"""

# ╔═╡ 544642df-6a72-414a-9cd8-5f3fda9840e8
"""
	model3(stimulus::TactileObject, params::ModelParams)

Decoder model that takes into account category and sensory variability; an optimal Bayesian decoder
"""
function model3(stimulus::TactileObject, params::ModelParams; debug=false)
	# Sense stimulus
	sensations = [encode(stimulus, params.Σ_encoder) for _=1:params.num_encodings]
	# Incorporate sensory variability
	estimated_stimulus = decode_sense(sensations, params.Σ_encoder, params.feat_ranges)
	# Incorporate category variability
	estimated_category = decode_category(estimated_stimulus, params.Σ_category)
	debug && model_debug(stimulus, estimated_stimulus, estimated_category)
	return estimated_category
end

# ╔═╡ 57e2528c-008d-4d8b-8015-066da346fc5c
md"""
# Experiment: Classification Task

Given a tactile object, determine if it is an Elyk or Noek
"""

# ╔═╡ 03c92b9a-a0be-4abd-bcd1-e3b618b19a1b
struct CategoryTaskParams
	runs::Integer
end

# ╔═╡ c018a27f-4f07-4c27-9982-37e0872e3d6f
"""
	categorize(object, model, model_params, exp_params)

Compute the proportion of runs where model correctly categorizes object
"""
function categorize(stimulus::TactileObject,
	model::Function,
	model_params::ModelParams,
	exp_params::CategoryTaskParams)::Real
	# Simulate many categorization runs for given object
	correct_categorizations = 0
	for run in 1:exp_params.runs
		response = nothing
		while true
			response = model(stimulus, model_params)
			maximum(response) == 0.5 || break
		end
		if argmax(response) == stimulus.category
			correct_categorizations += 1
		end
	end
	return correct_categorizations / exp_params.runs
end

# ╔═╡ c7b2a8c7-1ec9-43a8-a9c8-22b1d21b6397
begin
	function categorization_example()
		an_elyk = TactileObject(elyk, feat_ranges)

		model_params = ModelParams(feat_ranges, 10, (1.4^2)I, cov(elyk_dist))

		exp_params = CategoryTaskParams(400)

		p_correct_1 = categorize(an_elyk, model1, model_params, exp_params)
		p_correct_2 = categorize(an_elyk, model2, model_params, exp_params)
		p_correct_3 = categorize(an_elyk, model3, model_params, exp_params)

		results = (model1=p_correct_1, model2=p_correct_2, model3=p_correct_3)
		return (stimulus=an_elyk, performance=results)
	end

	categorization_ex = categorization_example()

	Markdown.parse("""
	Example Elyk:
		```$(categorization_ex.stimulus)```

	Model performance:
		```$(categorization_ex.performance)```
	""")
end

# ╔═╡ 358c7a37-a849-4d51-b4f4-f39c61e820c9
"""
	categorization_exp(num_objects)

For each model:
* Present model with an object, either elyk or noek, to categorize.
* Run ```$(categorize)()``` method
* Repeat for desired number of objects
* Record proportion of correct categorizations
"""
function categorization_exp(num_objects=1)
	rand_elyks = [TactileObject(elyk, feat_ranges) for _=1:num_objects]
	rand_noeks = [TactileObject(noek, feat_ranges) for _=1:num_objects]

	model_params = ModelParams(feat_ranges, 10, (1.4^2)I, cov(elyk_dist))
	exp_params = CategoryTaskParams(100)

	elyk_results = Dict([
		:model1=>zeros(num_objects),
		:model2=>zeros(num_objects),
		:model3=>zeros(num_objects),
	])

	noek_results = Dict([
		:model1=>zeros(num_objects),
		:model2=>zeros(num_objects),
		:model3=>zeros(num_objects),
	])

	for (i, (rand_elyk, rand_noek)) in enumerate(zip(rand_elyks, rand_noeks))
		for model in [model1, model2, model3]
			# Elyk
			elyk_results[Symbol(model)][i] = categorize(rand_elyk, model, model_params, exp_params)
			# Noek
			noek_results[Symbol(model)][i] = categorize(rand_noek, model, model_params, exp_params)
		end
	end
	return (elyks=rand_elyks,
			elyk_results=elyk_results,
			noeks=rand_noeks,
			noek_results=noek_results)
end

# ╔═╡ 3f4d04c4-e9d2-486c-bd53-9cf4e0801b32
# Categorization Experiment Code
categorization_exp(3)

# ╔═╡ a5532278-20ff-4ad8-8be1-47de0e39ade4
# Plots
# Plot scatter of object features that each model got correct more than 50% of time

# Plot bars of model performance

# ╔═╡ d5313c20-1e37-480e-a751-3c66c306ca6d
begin
	num_objects = 20
	results = categorization_exp(num_objects)
end

# ╔═╡ 04031f0c-45c2-4c94-8bb7-84dc35427d53
elyk_performance = vcat(collect(values(results.elyk_results)))

# ╔═╡ e7b50234-9c24-438f-aff8-c065c534f741
begin
	Plots.plot(1:num_objects, elyk_performance)
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[compat]
Distributions = "~0.25.52"
Plots = "~1.27.1"
"""

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

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "af92965fb30777147966f58acb05da51c5616b5f"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.3"

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

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

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

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

[[deps.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 = "4b859a208b2397a7a623a03449e4636bdb17bcf2"
uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
version = "1.16.1+1"

[[deps.Calculus]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad"
uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
version = "0.5.1"

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

[[deps.ChangesOfVariables]]
deps = ["ChainRulesCore", "LinearAlgebra", "Test"]
git-tree-sha1 = "bf98fa45a0a4cee295de98d4c1462be26345b9a1"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.2"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.DualNumbers]]
deps = ["Calculus", "NaNMath", "SpecialFunctions"]
git-tree-sha1 = "90b158083179a6ccbce2c7eb1446d5bf9d7ae571"
uuid = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74"
version = "0.6.7"

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.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 = "a6c850d77ad5118ad3be4bd188919ce97fffac47"
uuid = "d2c73de3-f751-5644-a686-071e5b155ba9"
version = "0.64.0+0"

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

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

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

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

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

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

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

[[deps.HypergeometricFunctions]]
deps = ["DualNumbers", "LinearAlgebra", "SpecialFunctions", "Test"]
git-tree-sha1 = "65e4589030ef3c44d3b90bdc5aac462b4bb05567"
uuid = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
version = "0.3.8"

[[deps.IniFile]]
git-tree-sha1 = "f550e6e32074c939295eb5ea6de31849ac2c9625"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.1"

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

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "91b5dcf362c5add98049e6c29ee756910b03051d"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.3"

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

[[deps.IterTools]]
git-tree-sha1 = "fa6287a4469f5e048d763df38279ee729fbd44e5"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.4.0"

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

[[deps.JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "abc9885a7ca2052a736a600f7fa66209f96506e1"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.4.1"

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

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

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

[[deps.LERC_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "bf36f528eec6634efc60d7ec062008f171071434"
uuid = "88015f11-f218-50d7-93a8-a6af411a945d"
version = "3.0.0+1"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.NaNMath]]
git-tree-sha1 = "737a5957f387b17e74d4ad2f440eb330b39a62c5"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "1.0.0"

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

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

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

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

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

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

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

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

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

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

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "85b5da0fa43588c75bb1ff986493443f821c70b7"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.3"

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

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

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

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

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

[[deps.Preferences]]
deps = ["TOML"]
git-tree-sha1 = "d3538e7f8a790dc8903519090857ef8e1283eecd"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.5"

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

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

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

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

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

[[deps.RecipesBase]]
git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.2.1"

[[deps.RecipesPipeline]]
deps = ["Dates", "NaNMath", "PlotUtils", "RecipesBase"]
git-tree-sha1 = "995a812c6f7edea7527bb570f0ac39d0fb15663c"
uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c"
version = "0.5.1"

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

[[deps.RelocatableFolders]]
deps = ["SHA", "Scratch"]
git-tree-sha1 = "cdbd3b1338c72ce29d9584fdbe9e9b70eeb5adca"
uuid = "05181044-ff0b-4ac5-8273-598c1e38db00"
version = "0.1.3"

[[deps.Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "838a3a4188e2ded87a4f9f184b4b0d78a1e91cb7"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.3.0"

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

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

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

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

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

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

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

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

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

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

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "5ba658aeecaaf96923dce0da9e703bd1fe7666f9"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.1.4"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "6976fab022fea2ffea3d945159317556e5dad87c"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.4.2"

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

[[deps.StatsAPI]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "c3d8ba7f3fa0625b062b82853a7d5229cb728b6b"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.2.1"

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.Unzip]]
git-tree-sha1 = "34db80951901073501137bdbc3d5a8e7bbd06670"
uuid = "41fe7b60-77ed-43a1-b4f0-825fd5a5650d"
version = "0.1.2"

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

[[deps.Wayland_protocols_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4528479aa01ee1b3b4cd0e6faef0e04cf16466da"
uuid = "2381bf8a-dfd0-557d-9999-79630e7b1b91"
version = "1.25.0+0"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.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:
# ╠═afdc0ce2-a8c5-11ec-0907-1b7d0d9db74a
# ╠═b4e08ac2-91c3-44c0-8a00-4b3938d5f0dd
# ╠═623b22e7-5c4e-4366-831c-5bc987eb8895
# ╟─0a5f2234-da1f-44f0-9e68-b2da7d9c4527
# ╟─e311bac7-d259-436b-afdb-4c32866c73b5
# ╠═0243fdb5-8109-411e-b8b0-1892e6f703d9
# ╠═67804342-2479-453b-b551-cc9508e6cbff
# ╠═b0ea37ee-8624-48ea-9be4-0b5be91a57cf
# ╟─6493908b-4773-4782-8616-158ad512ffcd
# ╠═4a531fe7-d8be-42fa-b356-24c36bb2d35d
# ╟─cbb058d9-a3bb-44d3-8dbc-e782dd27f172
# ╠═a5d1b3b8-5fe5-4791-bcfa-bfc74d1db345
# ╠═246ad164-b009-4666-8667-7d760ae09acc
# ╠═70ea935b-cbb7-48c4-810d-f70e3087c4f2
# ╟─f5e57a2d-b3b9-409a-8ced-56de60f85528
# ╠═929bea75-2ffa-4e48-b27c-4172eb8c1e51
# ╠═8fcb93f6-1544-4b84-984b-aa09c1630a07
# ╠═2d505656-00ca-455c-8710-dd2bf1cd6637
# ╟─91b45c05-59bf-42fb-8190-aa7918d23fba
# ╠═20a41be0-9b74-4e74-9e47-89c2e43b53a5
# ╠═5853bdb4-ecba-4504-bb18-5cb9bafca0f4
# ╠═a898528f-c913-4b8c-b814-9f0b98140979
# ╟─fa662a02-4d52-4487-869a-05cb1394662a
# ╟─bb6adffa-4e8a-481d-9f56-6e92192a0f1f
# ╠═5f38087b-17b0-4ef2-93b8-92c63f71aeeb
# ╠═4c0922fb-5a7d-457b-877a-da406a071c8c
# ╟─3c34fb99-7e4b-4387-b3b8-b4122f60b3f2
# ╠═e9c2b10a-c75e-438c-9860-591b5c64a55a
# ╠═ce3da79d-4555-4693-b8c5-30cff21d8306
# ╠═dd5aed91-b294-4962-ad80-c030d85181aa
# ╠═16ee34e1-540c-422f-90d2-8d173be877c4
# ╠═84cbe1fa-696c-456a-8803-9d12494df18b
# ╟─5d3a5b68-4345-4aa5-aa06-1022043aed85
# ╠═c1ab6a91-721b-4356-95b8-4d4d43500ad3
# ╟─cda8e5c4-8789-473e-bfe2-3c084ff96ccb
# ╠═ed2894cf-cfc8-4bb1-8441-8459d101087a
# ╠═ce6fdfe8-cd3e-469c-a384-25b564625d2d
# ╠═544642df-6a72-414a-9cd8-5f3fda9840e8
# ╟─57e2528c-008d-4d8b-8015-066da346fc5c
# ╠═03c92b9a-a0be-4abd-bcd1-e3b618b19a1b
# ╠═c018a27f-4f07-4c27-9982-37e0872e3d6f
# ╠═c7b2a8c7-1ec9-43a8-a9c8-22b1d21b6397
# ╠═358c7a37-a849-4d51-b4f4-f39c61e820c9
# ╠═3f4d04c4-e9d2-486c-bd53-9cf4e0801b32
# ╠═a5532278-20ff-4ad8-8be1-47de0e39ade4
# ╠═d5313c20-1e37-480e-a751-3c66c306ca6d
# ╠═04031f0c-45c2-4c94-8bb7-84dc35427d53
# ╠═e7b50234-9c24-438f-aff8-c065c534f741
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
"""
    IntegDomainModule

Module to manage integration domains.

An integration domain consists of the finite elements that approximate the
geometry, the function to supply the "missing" (other) dimension, indication
whether or not the integration domain represents an axially symmetric
situation, and integration rule used to evaluate integrals over the domain.
"""
module IntegDomainModule

__precompile__(true)

using ..FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import ..FESetModule: AbstractFESet, AbstractFESet0Manifold, AbstractFESet1Manifold, AbstractFESet2Manifold, AbstractFESet3Manifold, Jacobian, bfun, bfundpar
import ..FENodeSetModule: FENodeSet
import ..IntegRuleModule: AbstractIntegRule

"""
    IntegDomain{S<:AbstractFESet, F<:Function}

Integration domain.

- `T` = type of finite element set.  The type of the FE set will be dependent
  upon the operations required. For instance, for interior (volume) integrals
  such as body load or the stiffness hexahedral H8 may be used, whereas for
  boundary  (surface) integrals quadrilateral Q4 would be needed.
- `F` = type of function to return the "other" dimension.

An integration domain consists of the finite elements that approximate the
geometry, the function to supply the "missing" (other) dimension, indication
whether or not the integration domain represents an axially symmetric
situation, and integration rule used to evaluate integrals over the domain.
"""
mutable struct IntegDomain{S<:AbstractFESet, F<:Function}
    fes::S # finite element set object
    integration_rule::AbstractIntegRule  # integration rule object
    otherdimension::F # function to compute the "other" dimension (thickness, or cross-sectional area)
    axisymmetric::Bool
end

"""
    IntegDomain(fes::S, integration_rule::AbstractIntegRule) where {S<:AbstractFESet}

Construct with the default orientation matrix (identity), and the other
dimension  being the default 1.0.
"""
function  IntegDomain(fes::S, integration_rule::AbstractIntegRule) where {S<:AbstractFESet}
    return IntegDomain(fes, integration_rule, otherdimensionunity, false)
end

"""
    IntegDomain(fes::S, integration_rule::AbstractIntegRule,
      otherdimension::FFlt) where {S<:AbstractFESet}

Construct with the default orientation matrix (identity), and constant other
dimension.
"""
function  IntegDomain(fes::S, integration_rule::AbstractIntegRule, otherdimension::FFlt) where {S<:AbstractFESet}
    function otherdimensionfu(loc::FFltMat,
        conn::CC, N::FFltMat)::FFlt where {CC}
        return otherdimension::FFlt
    end
    return IntegDomain(fes, integration_rule, otherdimensionfu, false)
end

"""
    IntegDomain(fes::S, integration_rule::AbstractIntegRule,
      axisymmetric::Bool) where {S<:AbstractFESet}

Construct with the default orientation matrix (identity), for axially
symmetric models. The other dimension is the default unity (1.0).

This will probably be called when `axisymmetric = true`, since the default is
`axisymmetric = false`.
"""
function IntegDomain(fes::S, integration_rule::AbstractIntegRule, axisymmetric::Bool) where {S<:AbstractFESet}
    return IntegDomain(fes, integration_rule, otherdimensionunity, axisymmetric)
end

"""
    IntegDomain(fes::S, integration_rule::AbstractIntegRule, axisymmetric::Bool,
      otherdimension::FFlt) where {S<:AbstractFESet}

Construct for axially symmetric models. The other dimension is given as a number.
"""
function IntegDomain(fes::S, integration_rule::AbstractIntegRule, axisymmetric::Bool, otherdimension::FFlt) where {S<:AbstractFESet}
    function otherdimensionfu(loc::FFltMat, conn::CC, N::FFltMat)::FFlt where {CC}
        return otherdimension::FFlt
    end
    return IntegDomain(fes, integration_rule, otherdimensionfu, axisymmetric)
end

"""
    otherdimensionunity(loc::FFltMat, conn::CC, N::FFltMat)::FFlt
              where {CC}

Evaluate the other dimension: default is 1.0.
"""
function otherdimensionunity(loc::FFltMat, conn::CC, N::FFltMat)::FFlt where {CC}
    return 1.0
end

"""
    Jacobianpoint(self::IntegDomain{T}, J::FFltMat,
                loc::FFltMat, conn::CC,
                N::FFltMat)::FFlt where {T<:AbstractFESet0Manifold, CC}

Evaluate the point Jacobian.

- `J` = Jacobian matrix
- `loc` = location of the quadrature point in physical coordinates,
- `conn` = connectivity of the element,
- `N` = matrix of basis function values at the quadrature point.
"""
function Jacobianpoint(self::IntegDomain{T}, J::FFltMat, loc::FFltMat, conn::CC, N::FFltMat)::FFlt where {T<:AbstractFESet0Manifold, CC}
    return Jacobian(self.fes, J)::FFlt
end

"""
    Jacobiancurve(self::IntegDomain{T}, J::FFltMat,
        loc::FFltMat, conn::CC,
        N::FFltMat)::FFlt where {T<:AbstractFESet0Manifold, CC}

Evaluate the curve Jacobian.

- `J` = Jacobian matrix
- `loc` = location of the quadrature point in physical coordinates,
- `conn` = connectivity of the element,
- `N` = matrix of basis function values at the quadrature point.
"""
function Jacobiancurve(self::IntegDomain{T}, J::FFltMat, loc::FFltMat, conn::CC, N::FFltMat)::FFlt where {T<:AbstractFESet0Manifold, CC}
    Jac = Jacobianpoint(self, J, loc, conn, N)
    if self.axisymmetric
        return Jac*2*pi*loc[1];
    else
        return Jac*self.otherdimension(loc, conn,  N)
    end
end

"""
    Jacobiansurface(self::IntegDomain{T}, J::FFltMat,
                loc::FFltMat, conn::CC,
                N::FFltMat)::FFlt where {T<:AbstractFESet0Manifold, CC}

Evaluate the surface Jacobian.

For the zero-dimensional cell, the surface Jacobian is (i) the product of the
point Jacobian and the other dimension (units of length squared); or,  when
used as axially symmetric (ii) the product of the point Jacobian and the
circumference of the circle through the point `loc` times the other dimension
(units of length).

- `J` = Jacobian matrix
- `loc` = location of the quadrature point in physical coordinates,
- `conn` = connectivity of the element,
- `N` = matrix of basis function values at the quadrature point.
"""
function Jacobiansurface(self::IntegDomain{T}, J::FFltMat, loc::FFltMat, conn::CC, N::FFltMat)::FFlt where {T<:AbstractFESet0Manifold, CC}
    Jac = Jacobianpoint(self, J, loc, conn, N)::FFlt
    if self.axisymmetric
        return Jac*2*pi*loc[1]*self.otherdimension(loc, conn,  N);
    else
        return Jac*self.otherdimension(loc, conn,  N)
    end
end

"""
    Jacobianvolume(self::IntegDomain{T}, J::FFltMat,
          loc::FFltMat, conn::CC,
          N::FFltMat)::FFlt where {T<:AbstractFESet0Manifold, CC}

Evaluate the volume Jacobian.

For the zero-dimensional cell, the volume Jacobian is (i) the product of the
point Jacobian and the other dimension (units of length cubed); or,  when used
as axially symmetric (ii) the product of the point Jacobian and the
circumference of the circle through the point `loc` and the other dimension
(units of length squared).

- `J` = Jacobian matrix
- `loc` = location of the quadrature point in physical coordinates,
- `conn` = connectivity of the element,
- `N` = matrix of basis function values at the quadrature point.
"""
function Jacobianvolume(self::IntegDomain{T}, J::FFltMat, loc::FFltMat, conn::CC, N::FFltMat)::FFlt where {T<:AbstractFESet0Manifold, CC}
    Jac = Jacobianpoint(self, J, loc, conn, N)
    if self.axisymmetric
        return Jac*2*pi*loc[1]*self.otherdimension(loc, conn,  N);
    else
        return Jac*self.otherdimension(loc, conn,  N)
    end
end

"""
    Jacobianmdim(self::IntegDomain{T}, J::FFltMat,
      loc::FFltMat, conn::CC,
      N::FFltMat, m::FInt)::FFlt where {T<:AbstractFESet0Manifold, CC}

Evaluate the manifold Jacobian for an m-dimensional manifold.

For an 0-dimensional finite element,  the manifold Jacobian is for
- m=0: +1
- m=1: `Jacobiancurve`
- m=2: `Jacobiansurface`
- m=3: `Jacobianvolume`
"""
function Jacobianmdim(self::IntegDomain{T}, J::FFltMat, loc::FFltMat, conn::CC, N::FFltMat, m::FInt)::FFlt where {T<:AbstractFESet0Manifold, CC}
    @assert (m >= 0)  && (m <= 3)
    if (m==3)
        return Jacobianvolume(self, J, loc, conn, N)
    elseif (m==2)
        return Jacobiansurface(self, J, loc, conn, N)
    elseif (m==1)
        return Jacobiancurve(self, J, loc, conn, N)
    else # (m==0)
        return Jacobianpoint(self, J, loc, conn, N)
    end
end

"""
    Jacobiancurve(self::IntegDomain{T}, J::FFltMat,
              loc::FFltMat, conn::CC,
              N::FFltMat)::FFlt where {T<:AbstractFESet1Manifold, CC}

Evaluate the curve Jacobian.

- `J` = Jacobian matrix
- `loc` = location of the quadrature point in physical coordinates,
- `conn` = connectivity of the element,
- `N` = matrix of basis function values at the quadrature point.
"""
function Jacobiancurve(self::IntegDomain{T}, J::FFltMat, loc::FFltMat, conn::CC, N::FFltMat)::FFlt where {T<:AbstractFESet1Manifold, CC}
    return Jacobian(self.fes, J)::FFlt
end

"""
    Jacobiansurface(self::IntegDomain{T}, J::FFltMat,
                loc::FFltMat, conn::CC,
                N::FFltMat)::FFlt where {T<:AbstractFESet1Manifold, CC}

Evaluate the surface Jacobian.

For the one-dimensional cell,  the surface Jacobian is (i) the product of the
curve Jacobian and the other dimension (units of length); or,  when used as
axially symmetric (ii) the product of the curve Jacobian and the circumference
of the circle through the point `loc`.

- `J` = Jacobian matrix
- `loc` = location of the quadrature point in physical coordinates,
- `conn` = connectivity of the element,
- `N` = matrix of basis function values at the quadrature point.
"""
function Jacobiansurface(self::IntegDomain{T}, J::FFltMat, loc::FFltMat, conn::CC, N::FFltMat)::FFlt where {T<:AbstractFESet1Manifold, CC}
    Jac = Jacobiancurve(self, J, loc, conn, N)
    if self.axisymmetric
        return Jac*2*pi*loc[1];
    else
        return Jac*self.otherdimension(loc, conn,  N)
    end
end

"""
    Jacobianvolume(self::IntegDomain{T}, J::FFltMat,
                loc::FFltMat, conn::CC,
                N::FFltMat)::FFlt where {T<:AbstractFESet1Manifold, CC}

Evaluate the volume Jacobian.

For the one-dimensional cell,  the volume Jacobian is (i) the product of the
curve Jacobian and the other dimension (units of length squared); or,  when
used as axially symmetric (ii) the product of the curve Jacobian and the
circumference of the circle through the point `loc` and the other dimension
(units of length).

- `J` = Jacobian matrix
- `loc` = location of the quadrature point in physical coordinates,
- `conn` = connectivity of the element,
- `N` = matrix of basis function values at the quadrature point.
"""
function Jacobianvolume(self::IntegDomain{T}, J::FFltMat, loc::FFltMat, conn::CC, N::FFltMat)::FFlt where {T<:AbstractFESet1Manifold, CC}
    Jac = Jacobiancurve(self, J, loc, conn, N)
    if self.axisymmetric
        return Jac*2*pi*loc[1]*self.otherdimension(loc, conn,  N);
    else
        return Jac*self.otherdimension(loc, conn,  N)
    end
end

"""
    Jacobianmdim(self::IntegDomain{T}, J::FFltMat,
    loc::FFltMat, conn::CC,
    N::FFltMat, m::FInt)::FFlt where {T<:AbstractFESet1Manifold, CC}

Evaluate the manifold Jacobian for an m-dimensional manifold.

For an 1-dimensional finite element,  the manifold Jacobian is for
- m=1: `Jacobiancurve`
- m=2: `Jacobiansurface`
- m=3: `Jacobianvolume`
"""
function Jacobianmdim(self::IntegDomain{T}, J::FFltMat, loc::FFltMat, conn::CC, N::FFltMat, m::FInt)::FFlt where {T<:AbstractFESet1Manifold, CC}
    @assert (m >= 1) && (m <= 3)
    if (m==3)
        return Jacobianvolume(self, J, loc, conn, N)
    elseif (m==2)
        return Jacobiansurface(self, J, loc, conn, N)
    else # (m==1)
        return Jacobiancurve(self, J, loc, conn, N)
    end
end

"""
    Jacobiansurface(self::IntegDomain{T}, J::FFltMat,
                loc::FFltMat, conn::CC,
                 N::FFltMat)::FFlt where {T<:AbstractFESet2Manifold, CC}

Evaluate the surface Jacobian.

- `J` = Jacobian matrix
- `loc` = location of the quadrature point in physical coordinates,
- `conn` = connectivity of the element,
- `N` = matrix of basis function values at the quadrature point.
"""
function Jacobiansurface(self::IntegDomain{T}, J::FFltMat, loc::FFltMat, conn::CC, N::FFltMat)::FFlt where {T<:AbstractFESet2Manifold, CC}
    return Jacobian(self.fes, J)
end

"""
    Jacobianvolume(self::IntegDomain{T}, J::FFltMat,
                loc::FFltMat, conn::CC,
                N::FFltMat)::FFlt where {T<:AbstractFESet2Manifold, CC}

Evaluate the volume Jacobian.

For the two-dimensional cell,  the volume Jacobian is (i) the product of the
surface Jacobian and the other dimension (units of length); or,  when used as
axially symmetric (ii) the product of the surface Jacobian and the
circumference of the circle through the point `loc` (units of length).

- `J` = Jacobian matrix
- `loc` = location of the quadrature point in physical coordinates,
- `conn` = connectivity of the element,
- `N` = matrix of basis function values at the quadrature point.
"""
function Jacobianvolume(self::IntegDomain{T}, J::FFltMat, loc::FFltMat, conn::CC, N::FFltMat)::FFlt where {T<:AbstractFESet2Manifold, CC}
    Jac = Jacobiansurface(self, J, loc, conn, N)::FFlt
    if self.axisymmetric
        return Jac*2*pi*loc[1];
    else
        return Jac*self.otherdimension(loc, conn,  N)
    end
end

"""
    Jacobianmdim(self::IntegDomain{T}, J::FFltMat,
                loc::FFltMat, conn::CC,
                N::FFltMat, m::FInt)::FFlt where {T<:AbstractFESet2Manifold, CC}

Evaluate the manifold Jacobian for an m-dimensional manifold.

For an 2-dimensional finite element,  the manifold Jacobian is for
- m=2: `Jacobiansurface`
- m=3: `Jacobianvolume`
"""
function Jacobianmdim(self::IntegDomain{T}, J::FFltMat, loc::FFltMat, conn::CC, N::FFltMat, m::FInt)::FFlt where {T<:AbstractFESet2Manifold, CC}
    @assert (m >= 2) && (m <= 3)
    if (m==3)
        return Jacobianvolume(self, J, loc, conn, N)
    else # (m==2)
        return Jacobiansurface(self, J, loc, conn, N)
    end
end

"""
    Jacobianvolume(self::IntegDomain{T}, J::FFltMat,
                loc::FFltMat, conn::CC,
                N::FFltMat)::FFlt where {T<:AbstractFESet3Manifold, CC}
Evaluate the volume Jacobian.

- `J` = Jacobian matrix
- `loc` = location of the quadrature point in physical coordinates,
- `conn` = connectivity of the element,
- `N` = matrix of basis function values at the quadrature point.
"""
function Jacobianvolume(self::IntegDomain{T}, J::FFltMat, loc::FFltMat, conn::CC, N::FFltMat)::FFlt where {T<:AbstractFESet3Manifold, CC}
    return Jacobian(self.fes, J)::FFlt
end

"""
    Jacobianmdim{T<:AbstractFESet3Manifold}(self::IntegDomain, J::FFltMat,
                loc::FFltMat, conn::FIntMat, N::FFltMat, m::FInt)

Evaluate the manifold Jacobian for an m-dimensional manifold.

For an 3-dimensional cell,  the manifold Jacobian is
- m=3: `Jacobianvolume`
"""
function Jacobianmdim(self::IntegDomain{T}, J::FFltMat, loc::FFltMat, conn::CC, N::FFltMat, m::FInt)::FFlt where {T<:AbstractFESet3Manifold, CC}
    @assert (m == 3)
    return Jacobianvolume(self, J, loc, conn, N)
end

"""
    integrationdata(self::IntegDomain)

Calculate the data needed for  numerical quadrature for the integration rule
stored by the integration domain.
"""
function  integrationdata(self::IntegDomain)
    return integrationdata(self, self.integration_rule)
end

"""
    integrationdata(self::IntegDomain, integration_rule::T) where {T<:AbstractIntegRule}

Calculate the data needed for numerical quadrature.

For given integration domain, compute the quantities needed for numerical
integration.

# Return
`npts`, `Ns`, `gradNparams`, `w`, `pc` = number of quadrature points, arrays of
basis function values at the quadrature points,  arrays of gradients of basis
functions  with respect  to the parametric coordinates, array of weights and
array of locations of the quadrature points.
"""
function  integrationdata(self::IntegDomain, integration_rule::T) where {T<:AbstractIntegRule}
    pc::FFltMat = integration_rule.param_coords;
    w::FFltMat  =  integration_rule.weights ;
    npts::FInt = integration_rule.npts;
    # Precompute basis f. values + basis f. gradients wrt parametric coor
    Ns = FFltMat[];
    gradNparams = FFltMat[];
    for j=1:npts
        push!(Ns, bfun(self.fes,vec(pc[j,:])));
        push!(gradNparams, bfundpar(self.fes,vec(pc[j,:])));
    end
    return npts, reshape(Ns, 1,npts), reshape(gradNparams, 1,npts), w, pc
end

end

ext: jl


text:
module AltMusicTrigger

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Alt Music" => Ahorn.EntityPlacement(
        Maple.AltMusicTrigger,
        "rectangle"
    )
)

function Ahorn.editingOptions(trigger::Maple.AltMusicTrigger)
    return Dict{String, Any}(
        "track" => sort(collect(keys(Maple.Songs.songs)))
    )
end

end


text:
using Pkg; Pkg.activate(@__DIR__)
using MuJoCo
mj_activate("/home/taylor/.mujoco/bin/mjkey.txt") # set location to MuJoCo key path

using LyceumMuJoCo, LyceumMuJoCoViz
using FiniteDiff

using IterativeLQR
using LinearAlgebra
using Random

# ## load MuJoCo model
path = joinpath(@__DIR__, "../../../env/hopper/deps/hopper.xml")

include("mujoco_model.jl")
hopper = MuJoCoModel(path)
sim = LyceumMuJoCo.MJSim(hopper.m, hopper.d)

# @show hopper.m.body_mass
# @show hopper.m.body_inertia[:, :]
# length(hopper.d.ctrl)

# LyceumMuJoCoViz.visualize(sim) #, trajectories=[states]) # ctrl + LEFT (access trajectory mode)

T = 100
states = Array(undef, statespace(sim), T)
sim.d.qpos #.=
sim.d.qvel #.= zeros(7)
sim.d.ctrl
for t = 1:T
    # sim.d.qpos .= [0.0; 0.0; 0.55; 0.5]
    # sim.d.qvel .= [0.0; 0.0; 0.0; 0.0]
    sim.d.ctrl .= 1.0 * randn(3)#[0.0; 0.0; 1.0 * hopper.m.body_mass[2] * 9.81]#[5.0; 10.0; 0.0]
    LyceumMuJoCo.step!(sim)
    states[:, t] .= getstate(sim)
end
visualize(sim, trajectories=[states])
# 2.0 * hopper.m.body_mass[2] * 9.81

# ## horizon
T = 201

# ## hopper
nx = hopper.nx
nu = hopper.nu

# ## model
dyn = Dynamics(
    (y, x, u, w) -> f!(y, hopper, x, u),
    (dx, x, u, w) -> fx!(dx, hopper, x, u),
    (du, x, u, w) -> fu!(du, hopper, x, u),
    nx, nx, nu)

model = [dyn for t = 1:T-1]

# ## initial conditions
x1 = [0.0; 1.21; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0]
xM = [0.25; 2.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0]#[0.5; 0.5; 1.05; 0.0; 0.0; 0.0; 0.5; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0]
xT = [0.5; 1.21; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0]#[0.5; 0.5; 0.55; 0.0; 0.0; 0.0; 0.5; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0; 0.0]

# # rollout
ū = [0.01 * rand(3) for t = 1:T-1]
x̄ = rollout(model, x1, ū)

# ## objective
obj1 = (x, u, w) -> transpose(x - xM) * Diagonal([10.0 * ones(6); 1.0e-2 * ones(6)]) * (x - xM) + transpose(u) * Diagonal([1.0; 1.0; 1.0]) * u
obj2 = (x, u, w) -> transpose(x - xT) * Diagonal([10.0 * ones(6); 1.0e-2 * ones(6)]) * (x - xT) + transpose(u) * Diagonal([1.0; 1.0; 1.0]) * u
objT = (x, u, w) -> transpose(x - xT) * Diagonal([ones(6); ones(6)]) * (x - xT)

ct1 = IterativeLQR.Cost(obj1, nx, nu)
ct2 = IterativeLQR.Cost(obj2, nx, nu)
cT = IterativeLQR.Cost(objT, nx, 0)
obj = [[ct1 for t = 1:100]..., [ct2 for t = 1:100]..., cT]

# ## constraints
goal(x, u, w) = (x - xT)[collect(1:6)]

cont = Constraint()
conT = Constraint(goal, nx, 0)
cons = [[cont for t = 1:T-1]..., conT]

# ## problem
prob = solver(model, obj, cons,
    opts=Options(
        verbose=true,
        max_al_iter=5))
initialize_controls!(prob, ū)
initialize_states!(prob, x̄)

# ## solve
solve!(prob)

# ## solution
x_sol, u_sol = get_trajectory(prob)
@show IterativeLQR.eval_obj(prob.m_data.obj.costs, prob.m_data.x, prob.m_data.u, prob.m_data.w)
@show prob.s_data.iter[1]
@show norm(goal(prob.m_data.x[T], zeros(0), zeros(0)), Inf)

# ##
# using Plot
# plot(hcat(x_sol...)')
# plot(hcat(u_sol..., u_sol[end])', linetype=:steppost)

# ## visualize
x_vis = [[x̄[1] for t = 1:10]..., x̄..., [x̄[end] for t = 1:10]...]
x_vis = [[x_sol[1] for t = 1:100]..., x_sol..., [x_sol[end] for t = 1:100]...]
states = Array(undef, statespace(sim), length(x_vis))
for t = 1:length(x_vis)
    sim.d.qpos .= x_vis[t][1:6]
    sim.d.qvel .= x_vis[t][6 .+ (1:6)]
    # sim.d.ctrl .=
    # LyceumMuJoCo.step!(sim)
    states[:, t] .= getstate(sim)
end
visualize(sim, trajectories=[states])


text:
using PMDPs, PMDPs.LP
using JuMP, GLPK

@testset "LP" begin

    pp = simple_pp()
    traces = [simple_short_trace(pp), empty_trace(pp),]
    exhaust = simple_exhaust_capacity_trace(pp)
    diff = simple_trace_different_utility_revenue_optim_allocs(pp)
    mg = PMDPs.PMDPg(pp)

    # revenue
    r, u, alloc, action_seq, requests = LP.MILP_hindsight_pricing(mg, traces[1])
    @test @ntuple(r,u,alloc,action_seq) == (r=30.,u=3, alloc=[1.,1.,1.],action_seq=[10.,10.,10.])

    r, u, alloc, action_seq, requests = LP.MILP_hindsight_pricing(mg, traces[2])
    @test @ntuple(r,u,alloc,action_seq) == (r=0.,u=0, alloc=[],action_seq=[])

    r, u, alloc, action_seq, requests = LP.MILP_hindsight_pricing(mg, exhaust)
    @test @ntuple(r,u,alloc,action_seq) == (r=55.,u=6, alloc=[1., 1., 1., 0., 0.,1., 1.],
                action_seq = [10.0, 10.0, 10.0, PMDPs.REJECT_ACTION, PMDPs.REJECT_ACTION, 15.0, 10.0])

    r, u, alloc, action_seq, requests = LP.MILP_hindsight_pricing(mg, diff)
    @test @ntuple(r,u,alloc,action_seq) == (r=30.,u=4, alloc=[1., 1., 1., 0., 0.],
                action_seq = [10.0, 10., 10., PMDPs.REJECT_ACTION,  PMDPs.REJECT_ACTION])

    # utilization
    obj_val, u, alloc, action_seq, requests = LP.MILP_hindsight_pricing(mg, traces[1]; objective=:utilization)
    @test @ntuple(obj_val,u,alloc,action_seq) == (obj_val=3.,u=3, alloc=[1.,1.,1.],action_seq=[10.,10.,10.])

    obj_val, u, alloc, action_seq, requests = LP.MILP_hindsight_pricing(mg, traces[2]; objective=:utilization)
    @test @ntuple(obj_val,u,alloc,action_seq) == (obj_val=0.,u=0, alloc=[],action_seq=[])

    obj_val, u, alloc, action_seq, requests = LP.MILP_hindsight_pricing(mg, exhaust; objective=:utilization)
    @test @ntuple(obj_val,u,alloc,action_seq) == (obj_val=6.,u=6, alloc=[1., 1., 1., 0., 0.,1., 1.],
                action_seq = [10.0, 10.0, 10.0, PMDPs.REJECT_ACTION, PMDPs.REJECT_ACTION, 15.0, 10.0])

    obj_val, u, alloc, action_seq, requests = LP.MILP_hindsight_pricing(mg, diff; objective=:utilization)
    @test @ntuple(obj_val,u,alloc,action_seq) == (obj_val=5., u=5, alloc=[1., 1., 0., 1., 0.,],
                action_seq = [10., 10., PMDPs.REJECT_ACTION, 10., PMDPs.REJECT_ACTION])
end


text:
@testset "sparse(I, J, V, m, n, +)" begin
    m, n = 3, 5
    s, t, w = [1,2], [2,3], [0.5,0.5]

    test_rrule(sparse, s, t, w, m, n, +)
end

@testset "SparseMatrixCSC(A)" begin
    A = rand(5, 3)
    test_rrule(SparseMatrixCSC, A)
    test_rrule(SparseMatrixCSC{Float32,Int}, A, rtol=1e-4)
end

@testset "SparseVector(v)" begin
    v = rand(5)
    test_rrule(SparseVector, v)
    test_rrule(SparseVector{Float32}, Float32.(v), rtol=1e-4)
end

@testset "findnz" begin
    A = sprand(5, 5, 0.5)
    dA = similar(A)
    rand!(dA.nzval)
    I, J, V = findnz(A)
    V̄ = rand!(similar(V))
    test_rrule(findnz, A ⊢ dA, output_tangent=(zeros(length(I)), zeros(length(J)), V̄))

    v = sprand(5, 0.5)
    dv = similar(v)
    rand!(dv.nzval)
    I, V = findnz(v)
    V̄ = rand!(similar(V))
    test_rrule(findnz, v ⊢ dv, output_tangent=(zeros(length(I)), V̄))
end


text:
export BrickletAmbientLightV2IlluminanceCallbackThreshold
struct BrickletAmbientLightV2IlluminanceCallbackThreshold
    option::Char
    min::Integer
    max::Integer
end

export BrickletAmbientLightV2Configuration
struct BrickletAmbientLightV2Configuration
    illuminance_range::Integer
    integration_time::Integer
end

export BrickletAmbientLightV2Identity
struct BrickletAmbientLightV2Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletAmbientLightV2
"""
Measures ambient light up to 64000lux
"""
mutable struct BrickletAmbientLightV2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletAmbientLightV2(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_ambient_light_v2")
        deviceInternal = package.BrickletAmbientLightV2(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_illuminance
"""
    $(SIGNATURES)

Returns the illuminance of the ambient light sensor. The measurement range goes
up to about 100000lux, but above 64000lux the precision starts to drop.

.. versionchanged:: 2.0.2\$nbsp;(Plugin)
  An illuminance of 0lux indicates an error condition where the sensor cannot
  perform a reasonable measurement. This can happen with very dim or very bright
  light conditions. In bright light conditions this might indicate that the sensor
  is saturated and the configuration should be modified (:func:`Set Configuration`)
  to better match the light conditions.

If you want to get the illuminance periodically, it is recommended to use the
:cb:`Illuminance` callback and set the period with
:func:`Set Illuminance Callback Period`.
"""
function get_illuminance(device::BrickletAmbientLightV2)
    return device.deviceInternal.get_illuminance()
end

export set_illuminance_callback_period
"""
    $(SIGNATURES)

Sets the period with which the :cb:`Illuminance` callback is triggered
periodically. A value of 0 turns the callback off.

The :cb:`Illuminance` callback is only triggered if the illuminance has changed
since the last triggering.
"""
function set_illuminance_callback_period(device::BrickletAmbientLightV2, period)
    device.deviceInternal.set_illuminance_callback_period(period)
end

export get_illuminance_callback_period
"""
    $(SIGNATURES)

Returns the period as set by :func:`Set Illuminance Callback Period`.
"""
function get_illuminance_callback_period(device::BrickletAmbientLightV2)
    return device.deviceInternal.get_illuminance_callback_period()
end

export set_illuminance_callback_threshold
"""
    $(SIGNATURES)

Sets the thresholds for the :cb:`Illuminance Reached` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Callback is turned off"
 "'o'",    "Callback is triggered when the illuminance is *outside* the min and max values"
 "'i'",    "Callback is triggered when the illuminance is *inside* the min and max values"
 "'<'",    "Callback is triggered when the illuminance is smaller than the min value (max is ignored)"
 "'>'",    "Callback is triggered when the illuminance is greater than the min value (max is ignored)"
"""
function set_illuminance_callback_threshold(device::BrickletAmbientLightV2, option, min, max)
    device.deviceInternal.set_illuminance_callback_threshold(option, min, max)
end

export get_illuminance_callback_threshold
"""
    $(SIGNATURES)

Returns the threshold as set by :func:`Set Illuminance Callback Threshold`.
"""
function get_illuminance_callback_threshold(device::BrickletAmbientLightV2)
    return device.deviceInternal.get_illuminance_callback_threshold()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the period with which the threshold callbacks

* :cb:`Illuminance Reached`,

are triggered, if the thresholds

* :func:`Set Illuminance Callback Threshold`,

keep being reached.
"""
function set_debounce_period(device::BrickletAmbientLightV2, debounce)
    device.deviceInternal.set_debounce_period(debounce)
end

export get_debounce_period
"""
    $(SIGNATURES)

Returns the debounce period as set by :func:`Set Debounce Period`.
"""
function get_debounce_period(device::BrickletAmbientLightV2)
    return device.deviceInternal.get_debounce_period()
end

export set_configuration
"""
    $(SIGNATURES)

Sets the configuration. It is possible to configure an illuminance range
between 0-600lux and 0-64000lux and an integration time between 50ms and 400ms.

.. versionadded:: 2.0.2\$nbsp;(Plugin)
  The unlimited illuminance range allows to measure up to about 100000lux, but
  above 64000lux the precision starts to drop.

A smaller illuminance range increases the resolution of the data. A longer
integration time will result in less noise on the data.

.. versionchanged:: 2.0.2\$nbsp;(Plugin)
  If the actual measure illuminance is out-of-range then the current illuminance
  range maximum +0.01lux is reported by :func:`Get Illuminance` and the
  :cb:`Illuminance` callback. For example, 800001 for the 0-8000lux range.

.. versionchanged:: 2.0.2\$nbsp;(Plugin)
  With a long integration time the sensor might be saturated before the measured
  value reaches the maximum of the selected illuminance range. In this case 0lux
  is reported by :func:`Get Illuminance` and the :cb:`Illuminance` callback.

If the measurement is out-of-range or the sensor is saturated then you should
configure the next higher illuminance range. If the highest range is already
in use, then start to reduce the integration time.
"""
function set_configuration(device::BrickletAmbientLightV2, illuminance_range, integration_time)
    device.deviceInternal.set_configuration(illuminance_range, integration_time)
end

export get_configuration
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Configuration`.
"""
function get_configuration(device::BrickletAmbientLightV2)
    return device.deviceInternal.get_configuration()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletAmbientLightV2)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletAmbientLightV2, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
using Documenter
using JDF

makedocs(
    sitename = "JDF",
    format = Documenter.HTML(),
    modules = [JDF]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "<repository url>"
)=#

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
module DC
export DCPort, gen_dc_updates
using ...Network
using ...Signals

# struct DCPort{T <: AbstractFloat}
#     address::Address
#     current::T # value of the current from the port now.
#     stack::Vector{TimedDC{T}} # accepted instructions of pairs of start-time and current.
# end

# function gen_dc_updates(dt::T, ports::Vector{DCPort{T}}) where{T <: AbstractFloat}
#     reduce(ports; init=(zero(T), Vector{DCPort{T}}(undef, 0))) do acc, port
#         keep, take = take_due_signals(dt, port.stack)
#         #println("take: $(take).")
#         last_t_dc = reduce((acc2, x2) -> acc2.t <= x2.t ? x2 : acc2,
#                            take;
#                            init=TimedDC(- dt, port.current))
#         #println("last_t_dc: $(last_t_dc)")
#         (acc[1] + last_t_dc.current, [acc[2]..., DCPort(port.address, last_t_dc.current, keep)])
#     end
# end # return (i_syn, updates.ports_dc) for DC acceptor.

end # module end


text:
"""
    TypeO(var_1::TypeN)

Lorem ipsum dolor sit amet.
"""
struct TypeO

  var_1::TypeN

end


text:
include("examples/example_adjoint_nonuniform.jl");

q₀ = [ 0.60 * x for x in qₙ];
S = 1000;
β = 0.001;
w = 0.000051;                                                               # Априорный параметр апроксимации дельта-функции
q_final, Js, Qs = minimize(q₀, u₀, ulₘ, urₘ,
                           Xₙ, N, Tₘ, M,
                           ε, f1_data, f2_data,
                           S = S, β = β, w = w,
                           create_mesh=mshfrm,
                           showProgress = true);

frames = [1:20; 21:div(S,50):S];    # Первые двадцать без пропусков
frames = collect(1:div(S, 50):S);   # С пропусками, чтобы всего было 50 кадров
make_minimzation_gif(Js, Qs, qₙ, Xₙ, name = "Minimization_nonuniform.gif",
                     frames_to_write = frames, β = β)


text:
function getswimag(data, options)
    combined_mag = combine_echoes_swi(data.mag, data.TEs, options.mag_combine)
    savenii(combined_mag, "combined_mag", options.writesteps, data.header)
    swimag = sensitivity_correction(combined_mag, data, options)
    if options.mag_softplus !== false
        savenii(swimag, "sensitivity_corrected_mag", options.writesteps, data.header)
        swimag = softplus_scaling(swimag, options.mag_softplus)
    end
    savenii(swimag, "swimag", options.writesteps, data.header)
    return swimag
end

function sensitivity_correction(combined_mag, data, options)
    sensitivity = options.mag_sens
    if isnothing(sensitivity)
        sensitivity = getsensitivity(data.mag[:,:,:,1], getpixdim(data))
    end
    savenii(sensitivity, "sensitivity", options.writesteps, data.header)
    return combined_mag ./ sensitivity
end

#getpixdim(data) = data.header.pixdim[2:(1+ndims(data.mag))]
getpixdim(data) =  data.header.step

function combine_echoes_swi(mag, TEs, type)
    if ndims(mag) == 3 # only one echo
        return copy(mag)
    elseif type == :SNR
        return RSS(mag)
    elseif type == :average
        return dropdims(sum(mag; dims=4); dims=4)
    elseif type == :last
        return mag[:,:,:,end]
    elseif typeof(type) <: Pair
        type, para = type
        if type == :CNR
            if length(para) == 2
                (w1, w2) = para
                field = :B7T
            else
                (w1, w2, field) = para
            end
            weighting = calculate_cnr_weighting(TEs, w1, w2; field=field)
            return combine_weighted(mag, weighting)
        elseif type == :SE
            TE_SE = para
            return simulate_single_echo_mag(mag, TEs, TE_SE)
        elseif type == :closest
            TE_SE = para
            eco = findmin(abs.(TEs .- TE_SE))[2]
            return mag[:,:,:,eco]
        elseif type == :echo
            return mag[:,:,:,para]
        elseif type == :average
            return sum(mag[:,:,:,para]; dims=4)
        end
    end
    throw("ERROR: $type not defined for combination of echoes!")
end

function calculate_cnr_weighting(TEs, w1::Tuple, w2::Tuple; unusedargs...)
    S(TE, w) = w[1] * exp(-TE / w[2])
    w(TE, w1, w2) = S(TE, w1) - S(TE, w2)
    return w.(TEs, Ref(w1), Ref(w2))
end

function calculate_cnr_weighting(TEs, w1, w2; field=:B7T)
    T2s, factor = gettissue_easy(field)
    S(TE, tissue) = factor[tissue] * exp(-TE / T2s[tissue])
    w(TE, w1, w2) = S(TE, w1) - S(TE, w2)
    return w.(TEs, w1, w2)
end

function combine_weighted(mag, w)
    combined = zeros(size(mag)[1:3])
    for ieco in 1:size(mag, 4)
        combined .+= mag[:,:,:,ieco] .* w[ieco]
    end
    return combined ./ sum(w)
end

function softplus_scaling(mag, para)
    q = estimatequantile(mag, 0.8)
    if para === true
        return softplus.(mag, q/2)
    elseif para isa Number
        return softplus.(mag, para * q)
    elseif para isa Tuple
        return softplus.(mag, para[1] * q, para[2])
    else
        error("wrong input for softplus: got $para")
    end
end

function softplus(val, offset, factor=2)
    f = factor / offset
    sp(x) = log(1 + exp(f * (x - offset))) / f
    return sp(val) - sp(0)
end


text:
# Function to get the features for each grid
function getfeatures(grid)
  r = map((x) -> x.r, grid)
  g = map((x) -> x.g, grid)
  b = map((x) -> x.b, grid)
  return vec([
    mean(r), mean(g), mean(b),
    median(r), median(g), median(b),
    var(r), var(g), var(b)
  ])
end

# Take a guess at every element in the other image

img2 = load("./img/test.png");
gsize = 10;

# We do k-means but with a moving window
L = unique(labels)
cons = Dict([x => zeros(Float64, size(img2)) for x in L])
cons["STD"] = zeros(Float64, size(img2))
for i in 1:2:(size(img2,1)-gsize)
  for j in 1:2:(size(img2,2)-gsize)
    sp = img2[i:(i+gsize-1),j:(j+gsize-1)]
    sf = getfeatures(sp)
    cl = classify(1, [sf], training_set, labels, k=3)
    cons[cl][i:(i+gsize-1),j:(j+gsize-1)] .+= 1
    cons["STD"][i:(i+gsize-1),j:(j+gsize-1)] .+= 1
  end
end

for l in L
  cons[l] = cons[l]./cons["STD"]
  cons[l][find(cons[l].==0)] .= 0.0
end

C = fill(0, size(img2))
for i in eachindex(C)
  if cons["STD"][i] > 0
    u = [cons[l][i] for l in L]
    m = find(x -> x.== maximum(u), u)[1]
    C[i] = m
  end
end

bg = readpng("./img/test.png")
w = bg.width
h = bg.height

begin
  Drawing(w, h)
  origin()
  placeimage(bg, -w/2, -h/2)
  for i in 1:size(C,1)
    for j in 1:size(C,2)
      xinit = (i-1)
      yinit = (j-1)
      setopacity(0.4)
      sethue(colorant"#999")
      if C[i,j] > 0
        if L[C[i,j]] == "W"
          sethue("blue")
        end
        if L[C[i,j]] == "F"
          sethue("green")
        end
        if L[C[i,j]] == "U"
          sethue("red")
        end
        if L[C[i,j]] == "N"
          sethue("purple")
        end
        rect(yinit-w/2, xinit-h/2, 1, 1, :fill)
      end
    end
  end
  finish()
  preview()
end

max_stars_repo_name: tpoisot/WaterFinder

max_issues_repo_name: tpoisot/WaterFinder

max_forks_repo_name: tpoisot/WaterFinder


text:
import QuasiGeostrophy: compute
struct Wrapper{T, S} <: AbstractExpression
    data::T
    meta_data::S
end
struct WrapperMetaData{T}
    io_name::T
end

function Base.show(io::IO, field::Wrapper{T, S}) where {T <: Char, S}
    color = 230
    printstyled(io, field.data, color = color)
end
function Base.show(io::IO, field::Wrapper{T, S}) where {T, S <: WrapperMetaData}
    color = 230
    printstyled(io, field.meta_data.io_name, color = color)
end

compute(a::Wrapper) = a.data

macro wrapper(expr)
    rewritten_expr = _wrapper(expr)
    return rewritten_expr
end

function _wrapper(expr::Expr)
    symb = expr.args[1]
    val  = expr.args[2]
    if expr.head != :(=)
        println( "@wrapper macro not in proper form")
        println( "must be ")
        println( "@wrapper a=1 b=2 c=3")
        return error()
    end
    string_symb = String(symb)
    new_expr = :($(esc(symb)) =  Wrapper($val, WrapperMetaData($string_symb)))
    return new_expr
end

macro wrapper(exprs...)
    rewritten_exprs = [_wrapper(expr) for expr in exprs]
    return Expr(:block, rewritten_exprs...)
end

function rhs(qg_system)
    qg_system[1].lhs.data.data .= evaluate(qg_system[1].rhs)
    qg_system[2].lhs.data.data .= evaluate(qg_system[2].rhs)
    rhs¹ = evaluate(qg_system[3].rhs)
    rhs² = evaluate(qg_system[4].rhs)
    return [rhs¹, rhs²]
end
function evolve_system(qg_system, Δt; filter = 1)
    data = [qg_system[3].lhs.operand.data.data, qg_system[4].lhs.operand.data.data]
    olddata = copy(data)
    stage1 = rhs(qg_system)
    data .+= Δt .* stage1
    stage2 = rhs(qg_system)
    data .= olddata .+ Δt/2 .* (stage1 + stage2)
    qg_system[3].lhs.operand.data.data .= filter * data[1]
    qg_system[4].lhs.operand.data.data .= filter * data[2]
    return nothing
end
#
function cfl(ψ¹, ψ², ∂x, ∂y)
    F = ψ¹.data.metadata.transform.backward
    x,y = ψ¹.data.metadata.grid.grid
    Δx = x[2] - x[1]
    Δy = y[2] - y[1]
    u1 = F * evaluate(∂y(ψ¹))
    v1 = F * evaluate(∂x(ψ¹))
    u2 = F * evaluate(∂y(ψ²))
    v2 = F * evaluate(∂x(ψ²))
    u1 = maximum(abs.(u1))
    v1 = maximum(abs.(v1))
    u2 = maximum(abs.(u2))
    v2 = maximum(abs.(v2))
    return minimum([Δx / u1, Δx / u2, Δy / v1, Δy / v2]) # gives a Δt
end

# building some cutoff filters
function cutoff(k, cutoff)
    return k > cutoff ? 0.0 : 1.0
end

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: upiterbarg/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy