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

using Markdown
using InteractiveUtils

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

# ╔═╡ 6233804e-4b97-11ec-03c5-bbfcb0bd5450
using PlutoUI, Plots

# ╔═╡ a3562602-df40-4b5c-b191-296c26af43a4
md"""# Model of muscle force generation

This notebook works through the model of Fuglevand et al., 1993 going from motoneuron excitation to force production. Most justifcations for model components and parameter values are explained in the paper, this notebook simply aims to illustrate the different components and how they combine.

The model is split into two parts:
- the motoneuron model
- the motor unit-force model

Fuglevand et al., also include a model of surface EMG but I've not included that yet.

"""

# ╔═╡ 4e29217c-1b5a-4ac9-81fb-66b52829403d
md"""## Motoneuron model

The motoneuron combines signals from multiple sources and produces an output signal to a motor unit which in turn activates a group of muscle fibres. The input is not important in this model, only the output, so it is a simplified representation.
"""

# ╔═╡ 7410e83c-3083-4383-ae5e-eca6f7cf47ad
md""" ### Recruitment of a motoneuron

Two parameters are defined for simulation of motoneuron recruitment:

- Excitatory drive (``E(t)``): this represents the net input to the neurons trigger zone
- Recruitment threshold excitation (``\mathit{RTE}``): the minimum level of excitatory drive needed to initate a repetitive discharge in the neuron

A motoneuron remains inactive as long as its excitatory drive remains below its recuitment threshold. The first discharge of a motoneuron (becoming active) occurs when it's excitatory drive first equals or exceeds its recruitment threshold.
"""

# ╔═╡ 20737a9f-1c82-4ffb-b448-d75c45f96398
md"""#### Excitatory drive
All motoneurons are modeled to receive the same level of ``E``.
"""

# ╔═╡ 1e07451f-b0c5-4b64-bc42-e0880fbaff4b
excitation(t) = t ≤ 1.0 ? t : 1.0

# ╔═╡ 4362aebc-1df4-4a89-aabe-500b350f6db1
md"""#### Recruitment threshold excitation
Across the pool of neurons, many are given a low threshold and relatively few given a high threshold. This distribution is modelled using an exponential:

`` \mathit{RTE}_i = e^{a \cdot i}``

where ``a = ln(RR) / n``, ``RR`` is the recruitment range (i.e. ``RR=30`` is a desired 30-fold range in recruitment thresholds) and ``n`` is the number of motoneurons (``i`` goes from ``0`` to ``n-1``).

"""

# ╔═╡ 61bd56e7-9e74-41c4-86e7-ecb43c6c49bd
md"""### Rate coding

Once activated a motoneuron will repeatedly discharge.

The discharge behaviour is characterised by four parameters:
- minimum steady firing rate
- excitatory drive-firing rate relationship
- peak firing rate
- variability in interspike interval (ISI) under constant excitation

"""

# ╔═╡ 39598497-52b4-497a-8463-2c51dbb07fda
md"""#### Minimum steady firing rate
The minimum firing rate (``\mathit{MFR}``) was the same for all motoneurons and set at a value of ``8`` impulses/s.
"""

# ╔═╡ 69bed5f3-6921-4b26-832c-42617eaab6f1
md"""#### Excitatory drive-firing rate relationship
The firing rate (``\mathit{FR}``) of a motoneuron increases linearly with excitatory drive above the recruitment threshold. The gain (slope) of this relationship is assigned the same value for all motoneurons (``1.0``). Thus the ``\mathit{FR}``-``\mathit{E}`` relationship for each motoneuron ``i`` was given by:

`` \mathit{FR}_i(t) = \begin{cases} 0, & E(t) \lt \mathit{RTE}_i \\ \mathit{MFR} + g_e  \cdot  (E(t) - \mathit{RTE}_i), & E(t) \ge \mathit{RTE}_i \end{cases}``

As the ``\mathit{FR}`` is the inverse of the ``\mathit{ISI}`` and the ``\mathit{ISI}`` being the difference in time between two consecutive discharges:

`` \frac{1}{(t_{ij} - t_{ij-1})} = \begin{cases} 0, & E(t) \lt \mathit{RTE}_i \\ \mathit{MFR} + g_e  \cdot  (E(t) - \mathit{RTE}_i), & E(t) \ge \mathit{RTE}_i \end{cases}``

and therefore the time of the ``j``th discharge in motoneuron ``i`` (``t_{ij}``) can be written as:

``t_{ij} = \begin{cases} 0, & E(t) \lt \mathit{RTE}_i \\ \frac{1}{\mathit{MFR} + g_e  \cdot  (E(t) - \mathit{RTE}_i)} + t_{ij-1}, & E(t) \ge \mathit{RTE}_i \end{cases}``

The variabliity in the discharge time is added later.
"""

# ╔═╡ 21130a84-002e-478e-9d04-3e9d8746e1e1
md"""### Peak firing rate
How the peak firing rate (``\mathit{PFR}``) relates to ``\mathit{RTE}`` is unknown. The most common observation suggests that the earliest recruited motoneurons have the highest firing rates. However, some reports suggest that, over a large force range, the firing rate of high-threshold units eventually surpass that of lower threshold units. Both scenarios were tested here.

For the first scenario, the ``\mathit{PFR}`` of neuron ``i`` was scaled in inverse proportion to its ``\mathit{RTE}``:

`` \mathit{PFR}_i = \mathit{PFR}_1 - d_{PFR} \cdot \frac{\mathit{RTE}_i}{\mathit{RTE}_n}``

where ``d_{\mathit{PFR}}`` is the desired difference in peak firing rate between first and last recruited motoneurons and ``\frac{\mathit{RTE}_i}{\mathit{RTE}_n}`` is the recruitment threshold relative to the highest (and last) recruitment threshold. This makes motoneuron ``1`` have the highest firing rate and the ``n`` have the lowest.
"""

# ╔═╡ 0a8b53f0-624b-4afe-9111-3d21887ee91c
md"""### Interspike interval variability
The distribution of ISIs resembles a random process with a Gaussian probability distribution. The ISI calculated above are the mean of the distribution with the standard deviation calculated assuming a constant coefficient of variation (``\mathit{CV} = \sigma / \mu``):

`` \mathit{ISI}_{ij-1} = t_{ij} - t_{ij-1} = \mu + \sigma  \cdot  Z``

Which can be rewritten using the ``CV`` as:

`` t_{ij} = \mu + \mu \cdot CV \cdot Z + t_{ij-1}``

where ``\mu = \begin{cases} 0, & E(t) \lt \mathit{RTE}_i \\ \frac{1}{\mathit{MFR} + g_e  \cdot  (E(t) - \mathit{RTE}_i)}, & E(t) \ge \mathit{RTE}_i \end{cases}``

and ``Z`` is a Z-score drawn from a normally distributed array of values ranging from ``-3.9`` to ``3.9``.

The coefficient of variation was the same for each motoneuron and had a value of ``0.2``.
"""

# ╔═╡ c83e8528-916f-4dd0-a900-73892934743b
let
    μ = 10
    cv = 0.2
    σ = μ * cv
    zs = randn(1000) * 3.9
    isi = @. μ + σ * zs
    histogram(isi, legend = :none, xaxis = "ISI")
end

# ╔═╡ 31cedb3b-8f48-491b-a21a-046d0957dcb5
md"""### Maximum excitation
Maximum excitation (``E_\mathit{max}``) denotes the level of excitatory drive necessary to bring the last recuited motoneuron to its peak firing rate.

``E_{\mathit{max}}`` can be calculated as:

``\mathit{PFR}_n = \mathit{PFR}_1 - d_{\mathit{PFR}} \cdot \frac{\mathit{RTE}_n}{\mathit{RTE}_n} =\mathit{PFR}_1 - d_{\mathit{PFR}} ``

``\mathit{PFR}_n = \mathit{MFR} + g_e  \cdot  (E_{\mathit{max}} - \mathit{RTE}_n)``

``\mathit{PFR}_1 - d_{\mathit{PFR}} = \mathit{MFR} + g_e  \cdot  (E_{\mathit{max}} - \mathit{RTE}_n)``

``\frac{\mathit{PFR}_1 - d_{\mathit{PFR}} - \mathit{MFR}}{g_e} = (E_{\mathit{max}} - \mathit{RTE}_n)``, where ``\mathit{PFR}_1 - d_{\mathit{PFR}} - \mathit{MFR} = \mathit{PFR}_n - \mathit{MFR}``

``E_{\mathit{max}} = \mathit{RTE}_n + \frac{\mathit{PFR}_n - \mathit{MFR}}{g_e}``

Values of excitatory drive and recruitment threshold are then expressed as percentages of maximum excitation.
"""

# ╔═╡ f1079c38-4d4d-4646-9114-2fe78c8544e2
struct Motoneuron{T}
    threshold::T
    min_firing_rate::T
    peak_firing_rate::T
    ge::T
    cv::T

end

# ╔═╡ 396bb96e-9235-49b1-87f3-e753519eb86e
struct MotoneuronPool
    n::Int64
    rr::Float64
    Emax::Float64
    motoneurons::Vector{Motoneuron}

    function MotoneuronPool(n = 120, rr = 30.0, mfr = 8.0, pfr = 45.0, dpfr = 10.0, ge = 1.0, cv = 0.2)
        # recuitment thresholds
        a = log(rr) / (n - 1)
        rtes = [exp(a * i) for i = 0:n-1]

        # peak firing rates
        pfrs = [pfr - dpfr * ((rte - rtes[1]) / (rr - rtes[1])) for rte in rtes]

        # max excitation
        Emax = rtes[end] + (pfrs[end] - mfr) / ge

        # scale recruitment thresholds and gain by max excitation
        rtes /= Emax
        ge *= Emax

        # motoneurons
        mns = Motoneuron.(rtes, mfr, pfrs, ge, cv)

        return new(n, rr, Emax, mns)
    end
end

# ╔═╡ 869aa15d-777b-4930-a902-27ba09c93a8d
mean_firing_rate(m::Motoneuron, e) = e ≥ m.threshold ? min(m.peak_firing_rate, m.min_firing_rate + m.ge * (e - m.threshold)) : 0.0

# ╔═╡ 1b21b185-f528-44ec-8d31-d93b6dae69cb
md"""## Motor unit-force model

Once a motor neuron discharges it produces a motor unit twitch. This twitch can be well-approximated by:

``f_i(t) = \frac{P_i\cdot t}{T_i} \cdot  e^{\left(1 - \frac{-t}{T_i}\right)}``

where ``f_i(t)`` is the force of motor unit twitch ``i`` as a function of time, ``T_i`` is the rise (time to peak force) time to peak force of the twitch and ``P_i`` is the peak twitch force.

"""

# ╔═╡ 53e78998-8ce7-4ba5-b244-208ea8785c6a
md"""### Assignment of twitch parameters"""

# ╔═╡ 6affccfd-2e44-4692-b180-0d8c2a22e622
md"""#### Peak twitch force
Motor unit twitch force varies over a wide (100-fold) range of values, the distribution of motor unit's peak force is such that a large number produce small amounts of force whereas few produce large forces, and the force generated by a motor unit is related to its recruitment threshold. Therefore:

``P(i) = e^{b  \cdot i}``

where, similarly to before, `` b = \frac{ln(\mathit{RP})}{n}``, with ``\mathit{RP}`` being the range of peak forces and ``n`` being the number of motor units

"""

# ╔═╡ 6ea34497-caf0-4f4c-aa8a-b6143c605d81
md"""#### Contraction time
The range of motor unit contraction times is much smaller than the range of peak forces (2-5-fold), the frequency distribution is skewed towards briefer contraction times, and small force motor units have contraction times that vary over the full range of contraction time whereas intermediate to large force motor units have contraction times that are short and vary little. As such:

`` T(i) = T_L \cdot \left( \frac{1}{P_i} \right) ^{\frac{1}{c}}``

where ``T_L`` is the longest desired contraction time and ``c = \log_\mathit{RT}\mathit{RP}``, with ``\mathit{RT}`` being the desired range in contraction times.
"""

# ╔═╡ 269bdded-9e63-4930-9f84-ddbee51b90be
md""" Which gives a pool of motor units whose peak force and contraction time vary as:"""

# ╔═╡ 2f0a8a73-9ca1-43eb-90fc-9b9253b91c47
struct MotorUnit
    max_force::Float64
    contraction_time::Float64
end

# ╔═╡ e51cdb88-58b2-4884-80fb-cd4c4d76227a
struct MotorUnitPool
    n::Int64
    motorunits::Vector{MotorUnit}

    function MotorUnitPool(n = 120, rp = 100.0, rt = 3.0, tl = 0.09)
        # max forces
        b = log(rp) / (n - 1)
        ps = [exp(b * i) for i = 0:n-1]

        # contraction times
        c = log(rt, rp)
        ts = [tl * (1 / p)^(1 / c) for p in ps]

        # motor units
        mns = MotorUnit.(ps, ts)

        return new(n, mns)
    end

end

# ╔═╡ bf671903-767c-4b1f-a525-70bff744f776
twitch(mu::MotorUnit, t) = mu.max_force / mu.contraction_time * t * exp(1 - (t / mu.contraction_time))

# ╔═╡ e2036f7b-19ba-4224-85f5-a7a6289c7d3a
begin
	mup = MotorUnitPool(120, 100.0, 3.0, 0.09)
	plt_P = plot([x.max_force for x in mup.motorunits], title="forces", yaxis="force (aU)", xaxis="motor unit", legend=false)
	plt_T = plot([x.contraction_time for x in mup.motorunits], title="contraction times", yaxis="contraction time (s)", xaxis="motor unit", legend=false)
	plt_P_T = plot(plt_T[1][1][:y], plt_P[1][1][:y], st=:scatter, yaxis="force (au)", xaxis="contraction time (s)", legend=false)
	quiver!(plt_P_T[1][1][:x][[1, 60, 120]] .+ [0, 0.008, 0.01], plt_P_T[1][1][:y][[1, 60, 120]] .+ [24, 15, 0], quiver = ([0, -.008, -.01], [-24, -15, 0]), color = :black)
    annotate!(plt_P_T[1][1][:x][[1, 60, 120]] .+ [0, 0.009, 0.015], plt_P_T[1][1][:y][[1, 60, 120]] .+ [30, 18, 0], ["MU 1", "MU 60", "MU 120"], 10)
	plt_tw = plot(title="twitch forces", yaxis="force (au)", xaxis="time (s)")
	[plot!(plt_tw, t -> twitch(mup.motorunits[i], t), 0, 0.3, label="MU $i") for i in [1,40,80,120]]
end;

# ╔═╡ f6be63eb-a538-4e6f-ac9a-8db8b8183937
plot(plt_tw, title="")

# ╔═╡ 9242613e-2e19-4069-869b-2074d27e78ae
plot(plt_P, title="")

# ╔═╡ 3f832e48-076e-4fb0-9167-d6d59c8fce7b
plot(plt_T, title="")

# ╔═╡ f131b520-bd54-4185-81e4-e2ce35cfa603
plot(plt_P_T, title="")

# ╔═╡ 262e93a5-5144-4537-9ada-3d7717d808e2
md"""#### Force non-linearity

The motor unit twitch is assumed to be a linear process, this means that multiple twitches will summate and increase their force. This is seen in some extent in real motor units where the maximum force produced by a motor unit increases with firing rate. However, to stop the force becoming infinitely high a non-linear gain is imposed on the twitch force:

`` g_i(\mathit{FR}_{\mathit{norm}}) = \begin{cases} \frac{S(0.4)}{\mathit{FR}_{\mathit{norm}}}, & \mathit{FR}_{\mathit{norm}} \lt 0.4\\ \frac{S}{\mathit{FR}_{\mathit{norm}}}, & \mathit{FR}_{\mathit{norm}} \ge 0.4\end{cases}``

where ``S(\mathit{FR}_{\mathit{norm}}) = 1 - e^{-2\cdot\mathit{FR}_{\mathit{norm}}^3}``, ``\mathit{FR}_{\mathit{norm}}`` is the normalised firing rate (``T_i*\mathit{FR}``) and ``T_i`` is the contraction time for motor unit ``i``.

"""

# ╔═╡ c747175a-08ff-472b-9f7d-3c4feb82edd9
sigmoid(nfr) = 1 - exp(-2(nfr^3))

# ╔═╡ e3ce5924-656d-4453-93aa-8071182d45a7
sigmoid(mu::MotorUnit, fr) = sigmoid(mu.contraction_time * fr)

# ╔═╡ a563c138-1f12-496b-910e-cb507c635e9a
gain(nfr) = nfr > 0.4 ? sigmoid(nfr) / nfr : 1.0 * sigmoid(0.4) / 0.4

# ╔═╡ 9c734148-99bd-4e2c-b985-4735a432d811
gain(mu::MotorUnit, fr) = gain(mu.contraction_time * fr)

# ╔═╡ f3d42bf3-a473-416a-851b-ee3219797701
md"""#### Force response to a train of discharges
In motor unit ``i`` receiving a train containing ``k`` discharges the force produced by the unit is:

``F_i(t) = \sum\limits_{j=1}^k f_{ij}(t - t_{ij})``,  when ``t - t_{ij} \ge 0``

"""

# ╔═╡ 03e7dab6-b6d4-4e5a-b734-c4f707a8c959
md"""First need to transform the excitation of each motoneuron into a train of discharges:"""

# ╔═╡ 3e93f61d-a3fc-45c2-9cfe-914d03df6794
struct Stimulation{F}
    T::Float64
    excitation::F
    res::Int64 # samples/s

    Stimulation(T, excitation::F, res = 1000) where {F} = new{F}(T, excitation, res)
end

# ╔═╡ cf088732-a38b-454a-ad93-0455b03f1c91
time(s::Stimulation) = 0:1/s.res:s.T

# ╔═╡ d5aa3144-af4a-4af3-aca9-6a01ae4117f6
excitation(s::Stimulation) = s.excitation.(time(s))

# ╔═╡ 7971a856-54d7-4948-8304-e17b26060361
plot(excitation, 0, 3.0, legend=false, yaxis="excitation", xaxis="time")

# ╔═╡ eb5d2b39-1a3d-4ed9-bc9c-45aa189bec78
mean_firing_rate(s::Stimulation, mn::Motoneuron) = [mean_firing_rate(mn, e) for e in excitation(s)]

# ╔═╡ 4bc17da0-1820-41bc-822f-7e9ef4e1ec82
begin
	mnp = MotoneuronPool(120, 30.0, 8.0, 45.0, 10.0, 1.0, 0.2)
	plt_rte = plot([x.threshold for x in mnp.motoneurons], legend=false, title="recruitment thresholds (scaled)", yaxis="recuritment threshold", xaxis="motoneuron")
	plt_pfr = plot([x.peak_firing_rate for x in mnp.motoneurons], legend=false, title="peak firing rates", yaxis="peak firing rate (Hz)", xaxis="motoneuron")
	plt_efr = plot(title="excitation-firing rate relationship", yaxis="firing rate (hz)", xaxis="excitation", legend=:topleft)
	[plot!(plt_efr, e -> mean_firing_rate(mnp.motoneurons[i], e), 0, 1.0, label="MN $i") for i in [1,40,80,120]]
end;

# ╔═╡ cb3ed20d-9555-4c90-8178-39e2edd4ea0d
plot(plt_rte, title="")

# ╔═╡ 7e8dd3f2-4b50-4628-af92-9dc69da42926
plot(plt_efr, title="")

# ╔═╡ ca3fe97f-bc52-42d1-884c-3b80a573c9ab
plot(plt_pfr, title="")

# ╔═╡ e9aced3f-db1d-488f-9dca-a0359c89ff61
mean_isi(s::Stimulation, mn::Motoneuron) = 1 ./ mean_firing_rate(s, mn)

# ╔═╡ d93f8814-8cb9-41d3-b2ae-1fabfe233691
function discharges(s::Stimulation, mn::Motoneuron)
    res = s.res
    t = time(s)

    # get interspike intervals for a given excitation
    isi = mean_isi(s, mn)

    # find when motoneuron begins discharging
    i = findfirst(isi .!= Inf)
	isnothing(i) && return []
    train = [t[i]]

    # increment counter till after first discharge ended
    i += round(Int, isi[i] * res)

    # loop through until end of array
    while i < length(t)
        if isi[i] != Inf
            push!(train, t[i])
            i += round(Int, isi[i] * res)
            continue
        end
        i += 1
    end

    return train
end

# ╔═╡ 13722921-881c-4197-9c88-99b5b33cbf62
discharges(s::Stimulation, mnp::MotoneuronPool) = [discharges(s, mn) for mn in mnp.motoneurons]

# ╔═╡ 01cf60fd-6fc2-426b-bfe3-24e9340c29bb
begin
	s = Stimulation(2.0, excitation)
	trains = discharges(s, mnp)
	plt_ex = plot(time(s), excitation(s), yaxis="excitation", xaxis="time (s)", legend=false)
	function plot_train(train, yvals = (0,10); title="")
		xs = repeat(permutedims(hcat(train),(2,1)), 2)
		ys = permutedims(repeat(reshape([yvals...],1,2), size(xs, 2)...), (2,1))
		plt = plot(xs, ys, lc=:black, lw=3, legend=false, yaxis=false, yticks=[], xticks=[], title=title)

		return plt
	end
	plts = [plot_train(trains[idx], title="MU $idx") for idx in [1,30,60,90,120]]
	l = @layout [b;c;d;e;f;a{0.5h}]
	plot(plts...,plt_ex, link=:x, layout=l, titlefontsize=7, titlelocation=:left)

end

# ╔═╡ 76709458-664b-4067-8af2-2e6b70e864f9
md"""The train of discharges then determines the force of a motor unit at any instant:"""

# ╔═╡ 0067107e-ccc8-4207-97fd-f314baac72e7
function force(mu::MotorUnit, train, t)
    # train is a vector of discharge times
	isempty(train) && return 0.0

    # difference between discharge times and current time
    dts = t .- train

    # interspike intervals
    isis = diff(train)

    # firing rate = 1 / isi, initial firing rate = 0
    frs = [8, (1 ./ isis)...]
    nfrs = frs * mu.contraction_time # normalise

    # gain for given firing rate
    gs = gain.(nfrs)

    # force of each discharge
    fs = similar(train)
    fs[dts .≥ 0.0] = [twitch(mu, d_t) for d_t in dts if d_t ≥ 0.0] # only if t > d_t
    fs[dts .< 0.0] .= 0.0
    f = sum(fs .* gs) # sum across all discharges

    return f
end

# ╔═╡ 3f76fb97-febf-4c7d-abc7-814260051678
# savefig(plot(_plts[4], size=(400,400)), "muscle_model.png")

# ╔═╡ cbca5ee0-bb8b-416a-85b9-98ab26f4f95f
md"""##### See how excitation changes force for a motor unit"""

# ╔═╡ 0db994f7-31f3-4366-b17f-1c85be025d8b
md"""Excitiation:"""

# ╔═╡ 6670206e-fc4e-4463-b9e7-97575c3af300
@bind slider Slider(0.:.01:1.0, default=0.7, show_value=true)

# ╔═╡ 08a6bca6-d9e1-49ba-b9ea-96cd11e6bd9d
md"""Choose a motor unit (1-$(mup.n)):"""

# ╔═╡ 440cdecf-274e-4596-a116-1e2c5f0e6521
@bind mu_n TextField(default="60")

# ╔═╡ 713dd73c-b80b-414c-b91d-b60740d85208
md"""#### Total muscle force
The motor units are assumed to be independent of each other and therefore the total muscle force is simply the sum of each motor unit's force:

``F_M(t) = \sum \limits_{i=1}^n F_i(t)``

"""

# ╔═╡ 54103934-3af3-4e14-8360-9b4066c6c06c
struct Muscle
    n::Int64
	stim::Stimulation
    motoneurons::MotoneuronPool
    motorunits::MotorUnitPool

    function Muscle(stim; n = 120, rr = 30.0, mfr = 8.0, pfr = 45.0, dpfr = 10.0, ge = 1.0,
        cv = 0.2, rp = 100.0, rt = 3.0, tl = 0.09)
        mup = MotorUnitPool(n, rp, rt, tl)
        mnp = MotoneuronPool(n, rr, mfr, pfr, dpfr, ge, cv)

        return new(n, stim, mnp, mup)
    end

end

# ╔═╡ 52e102d0-f3ec-4983-841d-7d70d79f3219
force(mup::MotorUnitPool, trains, t) = [force(mu, train, t) for (mu,train) in zip(mup.motorunits, trains)] |> sum

# ╔═╡ 9688f3b8-e297-44f5-808e-a62466b71170
function force(m::Muscle, t)
	trains = discharges(s, m.motoneurons)
	f = force(m.motorunits, trains, t)

	return f
end

# ╔═╡ 7c8fb5b0-f71b-4482-9966-cb4e2f793a95
force(m::Muscle) = [force(m, t) for t in time(m.stim)]

# ╔═╡ 57ae3cde-16f3-4390-969c-a77530a5e5d2
begin
	plt_sig = plot(sigmoid, 0, 2.0, title="sigmoid relationship", yaxis="value", axis="normalised firing rate (T*FR)", legend=false)
	plt_gain = plot(gain, 0, 2.0, title="gain-firing rate relationship", yaxis="gain", xaxis="normalised firing rate (T*FR)", legend=false)
	plt_force1 = plot(nfr -> nfr*gain(nfr),0, 1.5, title="force-firing rate relationship", yaxis="normalised force", xaxis="normalised firing rate (T*FR)", legend=false)
	plt_force2 = plot(title="effect of firing rate on force", yaxis="force (au)", xaxis="time (s)", legend=:topleft)
	let
		tw1 = [force(mup.motorunits[120],[0.0], t) for t in 0:0.001:0.5]
		tw2 = [force(mup.motorunits[120],[0.5,0.57], t) for t in 0.5:0.001:1.0]
		plot!(plt_force2, 0:0.001:0.5, tw1, label="")
		plot!(plt_force2, 0.5:0.001:1.0, tw2, lc=1, label="")
		plot!([0.0 0.5 0.57; 0.0 0.5 0.57], [0 0 0; 5 5 5], lc=:black, lw=4, label=["discharges" "" ""])
	end
	plot(plt_sig, plt_gain, plt_force1, plt_force2, layout=(2,2), title="")
end

# ╔═╡ 2b2d0a99-01fc-424a-9792-0b0075f3bab8
let
	plts = []
	for idx in [1,40,80,120]
		_time = time(s)
		train = trains[idx]
		mu = mup.motorunits[idx]
		plt1 = plot_train(train)
		plt1[1].attr[:top_margin] = 0Plots.Measures.mm
		plt2 = plot(time(s), [force(mup.motorunits[idx], train, t) for t in _time], legend=false, title="MU $idx", titlefontsize=10, bottom_margin=0Plots.Measures.mm)
		plt = plot(plt2, plt1, layout=@layout [a{0.9h}; b])
		push!(plts, plt)
	end
	plot(plts...)
end

# ╔═╡ 227de63e-3e51-4d43-b567-dccc3ceabb2c
let
	_excitation(t) = 0.0 ≤ t ≤ 0.2 ? 0.0 : 0.2 < t ≤ 0.8 ? slider : 0.0
	s = Stimulation(1.5, _excitation)
	idx = parse(Int, mu_n)
	mn = mnp.motoneurons[idx]
	mu = mup.motorunits[idx]
	train = discharges(s, mn)
	_force = [force(mu, train, t) for t in time(s)]
	plt1 = plot(time(s), _force, ylims=(-0.5,3mu.max_force), legend=false, yaxis="force (au)", xticks=[], title="MU $idx")
	plt2 = plot_train(train)
	plt3 = plot(time(s), excitation(s), lc=:black, lw=3, grid=false, legend=false, yticks=[slider], ylims=(0,1.0), yaxis = false)
	plot(plt1, plt2, plt3, link=:x, layout=@layout [a{0.75h}; b; c])
end

# ╔═╡ c75e5bf4-c038-4e6f-9c61-d3916b233ad8
begin
	m = Muscle(s, n=120)
	m_force = force(m)
end;

# ╔═╡ 1d80cad8-1e57-42fe-b311-f778258dbf70
plot(time(s), m_force, title="muscle force", yaxis="force (au)", xaxis="time (s)", legend=false)

# ╔═╡ ac6f724f-1a7e-44d0-badf-04e6b820db81
md"""## Which model parameters determine the rate of force development (RFD) of the muscle?"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
Plots = "~1.23.6"
PlutoUI = "~0.7.19"
"""

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

# ╔═╡ Cell order:
# ╟─6233804e-4b97-11ec-03c5-bbfcb0bd5450
# ╟─a3562602-df40-4b5c-b191-296c26af43a4
# ╟─4e29217c-1b5a-4ac9-81fb-66b52829403d
# ╟─7410e83c-3083-4383-ae5e-eca6f7cf47ad
# ╟─20737a9f-1c82-4ffb-b448-d75c45f96398
# ╠═1e07451f-b0c5-4b64-bc42-e0880fbaff4b
# ╟─7971a856-54d7-4948-8304-e17b26060361
# ╟─4362aebc-1df4-4a89-aabe-500b350f6db1
# ╟─cb3ed20d-9555-4c90-8178-39e2edd4ea0d
# ╟─61bd56e7-9e74-41c4-86e7-ecb43c6c49bd
# ╟─39598497-52b4-497a-8463-2c51dbb07fda
# ╟─69bed5f3-6921-4b26-832c-42617eaab6f1
# ╟─7e8dd3f2-4b50-4628-af92-9dc69da42926
# ╟─21130a84-002e-478e-9d04-3e9d8746e1e1
# ╟─ca3fe97f-bc52-42d1-884c-3b80a573c9ab
# ╟─0a8b53f0-624b-4afe-9111-3d21887ee91c
# ╟─c83e8528-916f-4dd0-a900-73892934743b
# ╟─31cedb3b-8f48-491b-a21a-046d0957dcb5
# ╠═f1079c38-4d4d-4646-9114-2fe78c8544e2
# ╠═396bb96e-9235-49b1-87f3-e753519eb86e
# ╠═869aa15d-777b-4930-a902-27ba09c93a8d
# ╟─4bc17da0-1820-41bc-822f-7e9ef4e1ec82
# ╟─1b21b185-f528-44ec-8d31-d93b6dae69cb
# ╟─f6be63eb-a538-4e6f-ac9a-8db8b8183937
# ╟─53e78998-8ce7-4ba5-b244-208ea8785c6a
# ╟─6affccfd-2e44-4692-b180-0d8c2a22e622
# ╟─9242613e-2e19-4069-869b-2074d27e78ae
# ╟─6ea34497-caf0-4f4c-aa8a-b6143c605d81
# ╟─3f832e48-076e-4fb0-9167-d6d59c8fce7b
# ╟─269bdded-9e63-4930-9f84-ddbee51b90be
# ╟─f131b520-bd54-4185-81e4-e2ce35cfa603
# ╠═2f0a8a73-9ca1-43eb-90fc-9b9253b91c47
# ╠═e51cdb88-58b2-4884-80fb-cd4c4d76227a
# ╠═bf671903-767c-4b1f-a525-70bff744f776
# ╟─e2036f7b-19ba-4224-85f5-a7a6289c7d3a
# ╟─262e93a5-5144-4537-9ada-3d7717d808e2
# ╠═c747175a-08ff-472b-9f7d-3c4feb82edd9
# ╠═e3ce5924-656d-4453-93aa-8071182d45a7
# ╠═a563c138-1f12-496b-910e-cb507c635e9a
# ╠═9c734148-99bd-4e2c-b985-4735a432d811
# ╟─57ae3cde-16f3-4390-969c-a77530a5e5d2
# ╟─f3d42bf3-a473-416a-851b-ee3219797701
# ╟─03e7dab6-b6d4-4e5a-b734-c4f707a8c959
# ╠═3e93f61d-a3fc-45c2-9cfe-914d03df6794
# ╠═cf088732-a38b-454a-ad93-0455b03f1c91
# ╠═d5aa3144-af4a-4af3-aca9-6a01ae4117f6
# ╠═eb5d2b39-1a3d-4ed9-bc9c-45aa189bec78
# ╠═e9aced3f-db1d-488f-9dca-a0359c89ff61
# ╠═d93f8814-8cb9-41d3-b2ae-1fabfe233691
# ╠═13722921-881c-4197-9c88-99b5b33cbf62
# ╟─01cf60fd-6fc2-426b-bfe3-24e9340c29bb
# ╟─76709458-664b-4067-8af2-2e6b70e864f9
# ╠═0067107e-ccc8-4207-97fd-f314baac72e7
# ╟─2b2d0a99-01fc-424a-9792-0b0075f3bab8
# ╟─3f76fb97-febf-4c7d-abc7-814260051678
# ╟─cbca5ee0-bb8b-416a-85b9-98ab26f4f95f
# ╟─0db994f7-31f3-4366-b17f-1c85be025d8b
# ╟─6670206e-fc4e-4463-b9e7-97575c3af300
# ╟─08a6bca6-d9e1-49ba-b9ea-96cd11e6bd9d
# ╟─440cdecf-274e-4596-a116-1e2c5f0e6521
# ╟─227de63e-3e51-4d43-b567-dccc3ceabb2c
# ╟─713dd73c-b80b-414c-b91d-b60740d85208
# ╠═54103934-3af3-4e14-8360-9b4066c6c06c
# ╠═52e102d0-f3ec-4983-841d-7d70d79f3219
# ╠═9688f3b8-e297-44f5-808e-a62466b71170
# ╠═7c8fb5b0-f71b-4482-9966-cb4e2f793a95
# ╠═c75e5bf4-c038-4e6f-9c61-d3916b233ad8
# ╟─1d80cad8-1e57-42fe-b311-f778258dbf70
# ╟─ac6f724f-1a7e-44d0-badf-04e6b820db81
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

max_stars_repo_name: TomRottier/muscle-model

max_issues_repo_name: TomRottier/muscle-model

max_forks_repo_name: TomRottier/muscle-model


text:
module PFsimulator
using Robotlib
using DataFrames
# using Debug
using DSP
using Plots
import Plots.plot
using PFstateEstimator
import PyPlot
import Convex
import SCS
using GLM, DataFrames

include("particlePlot.jl")
export SimParams, SimData, generate_simdata, randse3, rand_traj, perturb_traj, traj2frames1D, traj2frames3D, ctraj_intrap, simulate_force, simulate_sensor, plot_tracking, run_tracking, get_search_traj, plot_errors, plot_meas, plot_trace, get_dataframe

sq(T) = squeeze(T[1:3,4,:],2)'
rms(x::Matrix) = sqrt(mean(x.^2,1))[:]
rms(x::Vector) = sqrt(mean(x.^2))

type SimParams
    h       # Sample time
    T       # Time steps
    σt      # std for translation error in FK
    σR      # std for rotation error in FK
    σq      # std for resolver/encoder error in joints
    σm      # Sensor noise
    σf      # Force noise
    force_ref
    Kreact  # Reactive force, f_react = Kreact * f
    Kcomp   # Compliance matrix
    T_TF_M  # Tool flange to sensor
end

type SimData{T}
    force::AbstractMatrix{T}
    meas::AbstractArray{T,3}
    traj_act::AbstractArray{T,3}
    traj_nom::AbstractArray{T,3}
    traj_meas::AbstractArray{T,3}
    search_traj::Vector{Matrix{T}}
    traj_type::Symbol
end

type SimResult{T}
    name::AbstractString
    x̂::AbstractArray{T,3}
    se::AbstractMatrix{T}
    et::AbstractVector{T}
    er::T
    simdata::SimData
    simparams::SimParams
    pfparams::PFparams
    trace::AbstractVector{PFtrace}
end
SimResult(name, x̂, se, et, er, simdata, simparams, pfparams) = SimResult(name, x̂, se, et, er, simdata, simparams, pfparams,PFtrace[])

# """
# `generate_simdata(traj_nom, p::SimParams)`
# """

function generate_simdata(traj_nom, p::SimParams; optimize = false)
    force       = simulate_force(traj_nom, p)
    traj_act, traj_meas = perturb_traj(traj_nom, force, p)
    search_traj = get_search_traj(traj_nom, p.T_TF_M, optimize)
    meas        = simulate_sensor(traj_nom, traj_act, search_traj, p)
    singular_axis = svd(squeeze(traj_nom[1:3,4,:],2))[1][:,3]
    traj_type   = singular_axis ⋅ [0,0,1] ≈ 1 ? :xy : :yz

    return SimData(force, meas, traj_act, traj_nom, traj_meas, search_traj,traj_type)
end

function randse3(σt,σR)
    R = expω(skew(σR*randn(3)))
    t = σt*randn(3)
    return [R t; 0 0 0 1]
end

function ctraj_intrap(N,T_start, T_end)
    T_out = zeros(4,4,N) # Initialize output trajectory
    xi_start = twistcoords(logT(T_start)) # Convert to start wrench
    xi_end = twistcoords(logT(T_end)) # Convert to end wrench
    xi_traj = zeros(6,N) # Initialize wrench trajectory
    for i = 1:6
        xi_traj[i,:] = linspace(xi_start[i], xi_end[i], N)
    end
    for n = 1:N
        T_out[:,:,n] = expm(skew4(xi_traj[:,n]))
    end
    return T_out
end

"""
`traj2frames1D(traj_in, increment)`
Returns frames along traj_in, with y along seam, z pointing down towards
the work object, and x is chosen so that the frames are orthonormal,
right oriented. T_out orientations are expressed in base frame.
"""
function traj2frames1D(traj_in, increment)
    N = length(traj_in)
    T_out = zeros(4,4,N) # Initialize output frames
    # differentiate:
    dx = centralDiff(traj_in)
    dy = convert(Float64, increment)
    for n = 1:N
        ŷ = [dx[n]; dy; 0.0] / norm([dx[n]; dy; 0.0]) # Points along seam
        ẑ = [0, 0, -1.0] # Always pointing down
        x̂ = ŷ × ẑ # Guarantees orthonormal, right oriented frames
        T_out[1:3, 1:3, n] = [x̂ ŷ ẑ]
        T_out[1:3, 4, n] = [traj_in[n], dy*(n-1), 0]
    end
    T_out[4,4,:] = 1
    return T_out
end

# TODO: test more thoroughly:
"""
´traj2frames3D(traj_in, fix = "z")´
Returns frames along traj_in (3*N matrix), with y along seam, z pointing down towards
the work object, and x is chosen so that the frames are orthonormal,
right oriented. T_out orientations are expressed in base frame. Fix is the direction in
which the seam does not vary.
"""
function traj2frames3D(traj_in, fix = "z")

    if size(traj_in,1) != 3
        traj_in = traj_in'
        if size(traj_in,1) != 3
            println("Input argument traj_in should be a 3*N matrix")
        end
    end

    N = size(traj_in,2)
    T_out = zeros(4,4,N) # Initialize output frames
    # differentiate to get seam direction:
    traj_dir = centralDiff(traj_in')'
    # Normalize to get unit direction vectors
    for n = 1:N
        traj_dir[:,n] /= norm(traj_dir[:,n])
    end
    for n = 1:N
        ŷ = vec(traj_dir[:,n]) # Points along seam
        if fix == "z"
            ẑ = [0, 0, -1.0]
            x̂ = ŷ × ẑ
        else
            x̂ = [-1., 0, 0]
            ẑ = x̂ × ŷ
        end
        T_out[1:3, 1:3, n] = [x̂ ŷ ẑ]
        T_out[1:3, 4, n] = traj_in[:,n]
    end
    T_out[4,4,:] = 1
    return T_out
end

function perturb_traj(traj_nom, force, p::SimParams)
    # Simulate low-frequency noise:
    fkin_noise          = similar(traj_nom)
    fkin_noise_twist    = zeros(6,size(traj_nom,3))
    for i = 1:size(traj_nom,3)
        fkin_noise[:,:,i]       = randse3(p.σt,p.σR)
        # fkin_noise[:,:,i] = randse3(1,1)
        fkin_noise_twist[:,i]   = twistcoords(logm(fkin_noise[:,:,i]))
    end

    #low-pass filter noise:
    L = min(500,p.T)
    fkin_noise_twist_filt = DSP.filtfilt(ones(L),[(L)], fkin_noise_twist')'
    fkin_noise_twist_filt = min(fkin_noise_twist_filt,1) # TODO: link this to dFK or set a max_error in simparams
    fkin_noise_twist_filt = max(fkin_noise_twist_filt,-1)

    # Convert fkin noise back to Ti:
    for i = 1:1:size(fkin_noise_twist_filt,2)
        fkin_noise[:,:,i] = expm(skew4(fkin_noise_twist_filt[:,i]))
    end

    # Addressing noise affecting resolvers/encoders:
    eq                  = p.σq*randn(6,size(traj_nom,3))
    xi                  = DH2twistsPOE(DH7600())
    q_center_of_traj    = ikinePOE(xi, traj_nom[:,:,round(Int,size(traj_nom,3)/2)], ones(6,1))[1] #ones seems to work as initial guess, but not zeros
    J                   = jacobianPOE(q_center_of_traj,xi)[1]
    J[1:3,:]           *=1000
    J[4:6,:]           /=1000
    e_C_twist           = J*eq # The covariance matrix of e_C is then J*σq*J'

    # TODO: perturb due to force, simulate deflections due to forces DONE!!
    force_perturb_twist =  p.Kcomp*force
    traj_meas           = similar(traj_nom)
    traj_act            = similar(traj_nom)
    force_perturb       = similar(traj_nom)
    e_C                 = similar(traj_nom)
    # Perturb traj_act and traj_meas:
    for i = 1:size(traj_nom,3)
        e_C[:,:,i]            = expm(skew4(e_C_twist[:,i]))
        force_perturb[:,:,i]  = expm(skew4(force_perturb_twist[:,i]))
        traj_meas[:,:,i]      = traj_nom[:,:,i]*e_C[:,:,i] # Noise affecting traj_meas
        traj_act[:,:,i]       = traj_nom[:,:,i]*fkin_noise[:,:,i]
        traj_act[1:3,1:3,i]   = traj_act[1:3,1:3,i]*force_perturb[1:3,1:3,i]
        traj_act[1:3,4,i]    += force_perturb[1:3,4,i]
        #resolver/encoder noise, and force perturbations, affecting traj_act
    end
    return traj_act, traj_meas
end

function perturb_simparams(si::SimParams, f)
    # force_ref
    s    = deepcopy(si)
    s.σt = rand(logspace(log10(s.σt/f),log10(s.σt*f),1000000))
    s.σR = rand(logspace(log10(s.σR/f),log10(s.σR*f),1000000))
    s.σq = rand(logspace(log10(s.σq/f),log10(s.σq*f),1000000))
    s.σm = Float64[rand(logspace(log10(s.σm[i]/f),log10(s.σm[i]*f),1000000)) for i = eachindex(s.σm)]
    s.σf = rand(logspace(log10(s.σf/f),log10(s.σf*f),1000000))
    s.force_ref = Float64[rand(logspace(log10(s.force_ref[i]/f),log10(s.force_ref[i]*f),1000000)) for i = eachindex(s.force_ref)]
    # s.T_TF_M[2,3,:] = Float64[rand(logspace(log10(s.T_TF_M[2,3,i]/f),log10(s.T_TF_M[2,3,i]*f),1000000)) for i = eachindex(s.T_TF_M[2,3,:])]
    return s
end

function perturb_pfparams!(p::PFparams, f)
    p.σW    = Float64[rand(logspace(log10(p.σW[i]/f),log10(p.σW[i]*f),1000000)) for i = eachindex(p.σW)]
    p.σV    = Float64[rand(logspace(log10(p.σV[i]/f),log10(p.σV[i]*f),1000000)) for i = eachindex(p.σV)]
    p.σVFK  = Float64[rand(logspace(log10(p.σVFK[i]/f),log10(p.σVFK[i]*f),1000000)) for i = eachindex(p.σVFK)]
end

function simulate_force(traj_nom, p::SimParams)
    return p.force_ref .+ p.σf.*randn(6,p.T) # + p.Kreact*p.force_ref
end

# TODO: Skriv vettiga test för var funktion i denna filen!!!
function simulate_sensor(traj_nom, traj_act, search_traj, p::SimParams)
    T_TF_M      = p.T_TF_M
    Nsensors    = size(T_TF_M,3)
    meas        = zeros(3,p.T,Nsensors)
    for j = 1:Nsensors
        meas[:,:,j] = PFstateEstimator.distToSeam([0,0], search_traj[j], traj_act, T_TF_M[:,:,j], tries = 5)[2]'
        for i = 1:p.T
            meas[3,i,j] = xyθ(traj_nom[:,:,i]*T_TF_M[:,:,j], traj_act[:,:,i]*T_TF_M[:,:,j]) # TODO: minus sometimes works better, find out the convention
        end
    end
    nZ = meas .!= 0
    meas += (p.σm.*randn(3,p.T,Nsensors)).*nZ
    return meas
end

function run_tracking(name, simdata, simparams, pfparams, seed=1; savetrace = true)
    srand(seed)
    T       = simparams.T
    x̂       = zeros(4,4,T)
    se      = zeros(3,T)
    meas    = permutedims(simdata.meas,[1,3,2])
    tic()
    state = init_pf(simdata.traj_meas[:,:,1], pfparams.Npart*pfparams.NpartInit, pfparams.σW, meas[:,:,1], simparams.T_TF_M)
    if savetrace
        trace   = Array(PFtrace,T)
        for i = 1:T
            x̂[:,:,i], se[:,i], state, trace[i] = PFstateEstimator.pfStateEstimator(state, meas[:,:,i], simdata.force[:,i], simdata.traj_meas[:,:,i], simdata.traj_nom, simdata.search_traj, simparams.T_TF_M, pfparams, savetrace = true)
        end
    else
        for i = 1:T
            x̂[:,:,i], se[:,i], state           = PFstateEstimator.pfStateEstimator(state, meas[:,:,i], simdata.force[:,i], simdata.traj_meas[:,:,i], simdata.traj_nom, simdata.search_traj, simparams.T_TF_M, pfparams, savetrace = false)
        end
    end
    t = toq()
    et, er = tracking_error(x̂, simdata)
    println(name, " took ",t, " seconds. Error ",round(et,3), "   ", round(er,3))

    savetrace ? SimResult(name, x̂, se, et, er, simdata, simparams, pfparams, trace) : SimResult(name, x̂, se, et, er, simdata, simparams, pfparams)

end

function plot_tracking(simresult)
    x̂ = simresult.x̂
    se = simresult.se
    time = linspace(0,simresult.simparams.h*simresult.simparams.T,simresult.simparams.T)
    p1 = subplot(time,sq(x̂),c=:blue,lab="x̂",nc=1, title=simresult.name)
    subplot!(time,sq(simresult.simdata.traj_act), c=:red, lab="traj_act")
    subplot!(time,sq(simresult.simdata.traj_nom),c=:green,lab="traj_nom")
    subplot!(time,sq(x̂)+2se',c=:blue,lab="+2σx̂", l=:dash)
    subplot!(time,sq(x̂)-2se',c=:blue,lab="-2σx̂", l=:dash)

    p = subplot(time,traj2quat(x̂)',c=:blue,lab="x̂",nc=1, title=simresult.name*" rotation")
    subplot!(time,traj2quat(simresult.simdata.traj_act)', c=:red, lab="traj_act")
    subplot!(time,traj2quat(simresult.simdata.traj_nom)',c=:green,lab="traj_nom")

    return p1
end

plot_angle(simresult) = plot(Rangle(simresult.x̂,simresult.simdata.traj_act,true),nc=1, title=simresult.name, lab="Angle (degrees)")

function plot_meas(dim,simparams,simdata,base="M")
    if dim == 0
        dim = 1:size(simparams.T_TF_M,3)
    end
    for d in dim
        T_TF_M = simparams.T_TF_M[:,:,d]
        if base == "M"
            plot(simdata.meas[:,:,d]')
        elseif base == "RB"
            meas = zeros(3,simparams.T)
            for i = 1:simparams.T
                meas[:,i] = (trinv(simdata.traj_nom[:,:,i]*T_TF_M[:,:,d])*[simdata.meas[1:2,i,d];0;1])[1:3]
            end
            plot(meas[1:3,:]', lab="Sensor $d")
        end
    end
end

function plot_meas(dim, simresult::SimResult, base="M")
    plot_meas(dim,simresult.simparams,simresult.simdata,base)
    title!(simresult.name)
end

function get_dataframe(simresultlists...; onehot=false)
    # σW::Vector{T}       # noise in state update
    # σV::Vector{T}       # Noise in sensor reading
    # σVFK::Vector{T}     # Noise in FK
    # σt      # std for translation error in FK
    # σR      # std for rotation error in FK
    # σq      # std for resolver/encoder error in joints
    # σm      # Sensor noise
    # σf      # Force noise
    # force_ref
    frame   = DataFrame()
    Nexp    = length(simresultlists)
    MC      = length(simresultlists[1])
    et      = [simresultlists[i][mc].et for i in 1:Nexp, mc in 1:MC]
    et0     = [norm(simresultlists[i][mc].simdata.traj_act[1:3,4,1]-simresultlists[i][mc].simdata.traj_nom[1:3,4,1]) for i in 1:Nexp, mc in 1:MC]
    names   = [simresultlists[i][mc].name for i in 1:Nexp, mc in 1:MC]
    traj_type = [string(simresultlists[i][mc].simdata.traj_type) for i in 1:Nexp, mc in 1:MC]
    ex      = Float64[et[i,mc][1] for i in 1:Nexp, mc in 1:MC]
    ey      = Float64[et[i,mc][2] for i in 1:Nexp, mc in 1:MC]
    ez      = Float64[et[i,mc][3] for i in 1:Nexp, mc in 1:MC]
    er      = Float64[simresultlists[i][mc].er for i in 1:Nexp, mc in 1:MC]
    et0     = Float64[et0[i,mc] for i in 1:Nexp, mc in 1:MC]

    Npart   = [simresultlists[i][mc].pfparams.Npart for i in 1:Nexp, mc in 1:MC]
    σW1      = [simresultlists[i][mc].pfparams.σW[1] for i in 1:Nexp, mc in 1:MC]
    σW2      = [simresultlists[i][mc].pfparams.σW[2] for i in 1:Nexp, mc in 1:MC]
    σW3      = [simresultlists[i][mc].pfparams.σW[3] for i in 1:Nexp, mc in 1:MC]
    σV1      = [simresultlists[i][mc].pfparams.σV[1] for i in 1:Nexp, mc in 1:MC]
    σV2      = [simresultlists[i][mc].pfparams.σV[2] for i in 1:Nexp, mc in 1:MC]
    σV3      = [simresultlists[i][mc].pfparams.σV[3] for i in 1:Nexp, mc in 1:MC]
    σVFK1    = [simresultlists[i][mc].pfparams.σVFK[1] for i in 1:Nexp, mc in 1:MC]
    σVFK2    = [simresultlists[i][mc].pfparams.σVFK[2] for i in 1:Nexp, mc in 1:MC]
    σVFK3    = [simresultlists[i][mc].pfparams.σVFK[3] for i in 1:Nexp, mc in 1:MC]
    σVFK4    = [simresultlists[i][mc].pfparams.σVFK[4] for i in 1:Nexp, mc in 1:MC]
    σVFK5    = [simresultlists[i][mc].pfparams.σVFK[5] for i in 1:Nexp, mc in 1:MC]

    σt    = [simresultlists[i][mc].simparams.σt for i in 1:Nexp, mc in 1:MC]
    σR    = [simresultlists[i][mc].simparams.σR for i in 1:Nexp, mc in 1:MC]
    σq    = [simresultlists[i][mc].simparams.σq for i in 1:Nexp, mc in 1:MC]
    σm1    = [simresultlists[i][mc].simparams.σm[1] for i in 1:Nexp, mc in 1:MC]
    σm2    = [simresultlists[i][mc].simparams.σm[2] for i in 1:Nexp, mc in 1:MC]
    σm3    = [simresultlists[i][mc].simparams.σm[3] for i in 1:Nexp, mc in 1:MC]
    σf    = [simresultlists[i][mc].simparams.σf for i in 1:Nexp, mc in 1:MC]
    force_ref1    = [simresultlists[i][mc].simparams.force_ref[1] for i in 1:Nexp, mc in 1:MC]
    force_ref2    = [simresultlists[i][mc].simparams.force_ref[2] for i in 1:Nexp, mc in 1:MC]
    force_ref3    = [simresultlists[i][mc].simparams.force_ref[3] for i in 1:Nexp, mc in 1:MC]
    force_ref4    = [simresultlists[i][mc].simparams.force_ref[4] for i in 1:Nexp, mc in 1:MC]
    force_ref5    = [simresultlists[i][mc].simparams.force_ref[5] for i in 1:Nexp, mc in 1:MC]
    force_ref6    = [simresultlists[i][mc].simparams.force_ref[6] for i in 1:Nexp, mc in 1:MC]

    sens0   = [all(simresultlists[i][mc].simdata.meas .== 0) for i in 1:Nexp, mc in 1:MC]
    sens1   = [!sens0[i,mc] && size(simresultlists[i][mc].simparams.T_TF_M,3) == 1 for i in 1:Nexp, mc in 1:MC]
    sens2   = [!sens0[i,mc] && size(simresultlists[i][mc].simparams.T_TF_M,3) == 2 for i in 1:Nexp, mc in 1:MC]
    sens    = map(Int,sens1) + 2map(Int,sens2)

    frame[:names]       = names[:]
    frame[:ex]          = ex[:]
    frame[:ey]          = ey[:]
    frame[:ez]          = ez[:]
    frame[:er]          = er[:]
    frame[:et0]         = et0[:]

    frame[:Npart]       = Npart[:]
    frame[:logNpart]    = map(log10,Npart[:])
    frame[:sens]        = sens[:]
    frame[:sens0]       = sens0[:]
    frame[:sens1]       = sens1[:]
    frame[:sens2]       = sens2[:]
    frame[:traj_type]   = traj_type[:]
    frame[:σW1]         = σW1[:]
    frame[:σW2]         = σW2[:]
    frame[:σW3]         = σW3[:]
    frame[:σV1]         = σV1[:]
    frame[:σV2]         = σV2[:]
    frame[:σV3]         = σV3[:]
    frame[:σVFK1]       = σVFK1[:]
    frame[:σVFK2]       = σVFK2[:]
    frame[:σVFK3]       = σVFK3[:]
    frame[:σVFK4]       = σVFK4[:]
    frame[:σVFK5]       = σVFK5[:]
    frame[:σt]          = σt[:]
    frame[:σR]          = σR[:]
    frame[:σq]          = σq[:]
    frame[:σm1]         = σm1[:]
    frame[:σm2]         = σm2[:]
    frame[:σm3]         = σm3[:]
    frame[:σf]          = σf[:]
    frame[:force_r1]  = force_ref1[:]
    frame[:force_r2]  = force_ref2[:]
    frame[:force_r3]  = force_ref3[:]
    frame[:force_r4]  = force_ref4[:]
    frame[:force_r5]  = force_ref5[:]
    frame[:force_r6]  = force_ref6[:]

    pool!(frame, [:names])
    pool!(frame, [:traj_type])
    return frame
end

function plot_errors(simresultlists...)
    Nexp    = length(simresultlists)
    MC      = length(simresultlists[1])
    et      = [simresultlists[i][mc].et for i in 1:Nexp, mc in 1:MC]
    names   = [simresultlists[i][1].name for i in 1:Nexp]
    ex      = Float64[et[i,mc][1] for i in 1:Nexp, mc in 1:MC]
    ey      = Float64[et[i,mc][2] for i in 1:Nexp, mc in 1:MC]
    ez      = Float64[et[i,mc][3] for i in 1:Nexp, mc in 1:MC]
    er      = Float64[simresultlists[i][mc].er for i in 1:Nexp, mc in 1:MC]
    PyPlot.figure()
    PyPlot.subplot(2,2,1)
    PyPlot.boxplot(ex',labels=names), PyPlot.title("Error x")
    PyPlot.subplot(2,2,2)
    PyPlot.boxplot(ey',labels=names), PyPlot.title("Error y")
    PyPlot.subplot(2,2,3)
    PyPlot.boxplot(ez',labels=names), PyPlot.title("Error z")
    PyPlot.subplot(2,2,4)
    PyPlot.boxplot(er',labels=names), PyPlot.title("Error rot")
    # PyPlot.show()

    # subplot([mx my mz mr], t=:bar, title=["Mean error X" "Mean error y" "Mean error z" "Mean error r"], xlabel=names)
    # subplot([max may maz mar], t=:bar, title=["Max error X" "Max error y" "Max error z" "Max error r"], xlabel=names)
    nothing
end

function plot_resample(res)
    N = res.pfparams.Npart
    plot([res.trace[i].Neff for i = 1:res.simparams.T],lab="Effective number of particles", title=res.name)
    scatter!(N*[res.trace[i].resample for i = 1:res.simparams.T], c=:green, lab="Resample")
    scatter!(1.2N*[res.trace[i].reset for i = 1:res.simparams.T]-0.1N, c=:red, m=:rect, lab="Reset")
    scatter!(1.4N*[res.trace[i].runMeasUpdate for i = 1:res.simparams.T]-0.2N, c=:magenta, m=:circle, lab="Measurement update")
end

function plot_trace(res; xIndices = 1:3, yIndices = 1:3)
    T       = res.simparams.T
    N       = res.pfparams.Npart
    y       = res.simdata.meas[:,:,1]
    xreal   = [sq(res.simdata.traj_act)'; R2rpy(res.simdata.traj_act,deg=true)]
    pdata   = Void
    T_TF_M  = res.simparams.T_TF_M[:,:,1]
    for i = 1:T
        ti      = res.trace[i]
        x       = [sq(ti.x)'; R2rpy(ti.x,deg=true)]
        w       = ti.w
        a       = ti.resample_i
        yhat    = PFstateEstimator.distToSeam([0,0], res.simdata.search_traj[1], ti.x, T_TF_M, tries = 20)[2]'
        for j = 1:N
            yhat[3,j] = xyθ(res.simdata.traj_nom[:,:,i]*T_TF_M, ti.x[:,:,j]*T_TF_M)
        end
        pdata = plotPoints(x, w, y, yhat, N, a, i, xreal, 0, 0, pdata, density = true, xIndices = xIndices, yIndices = yIndices, leftOnly = true)
    end

end

function Plots.plot(res::SimResult, args...)
    pyplot()
    for arg in args
        arg == :meas && plot_meas(0,res)
        arg == :tracking && plot_tracking(res)
        arg == :resample && plot_resample(res)
        arg == :angle && plot_angle(res)
        #TODO: nom, act, meas, x̂
    end
    if !isempty([:nom, :act, :traj_meas] ∩ args)
        subplot(n=3,nc=1)
        for arg in args
            arg == :nom && subplot!(sq(res.simdata.traj_nom), c=:green, lab="traj_nom")
            arg == :act && subplot!(sq(res.simdata.traj_act), c=:red, lab="traj_act")
            arg == :traj_meas && subplot!(sq(res.simdata.traj_meas), c=:magenta, lab="traj_meas")
        end
    end
end

function tracking_error(x̂, simdata)
    et = rms(sq(x̂) - sq(simdata.traj_act))
    er = rms(Rangle(x̂,simdata.traj_act,true))
    return et,er
end

# """
# This function calculates an efficient search trajectory.
# """

function get_search_traj(traj_nom,T_TF_M, optimize = false)
    Nsensors = size(T_TF_M,3)
    search_traj = deepcopy(traj_nom)
    search_traj = cat(4,search_traj,search_traj)
    for ns = 1:Nsensors
        for i = 1:size(search_traj,3)
            search_traj[:,:,i,ns] = search_traj[:,:,i,ns]*T_TF_M[:,:,ns]
        end
    end
    search_traj = permutedims(squeeze(search_traj[1:3,4,:,:],2),[2 1 3])
    rettraj = Array(Array{Float64,2},Nsensors)
    if !optimize
        for ns = 1:Nsensors
            D = [conv(diff(diff(search_traj[:,1,ns])),ones(3)) conv(diff(diff(search_traj[:,2,ns])),ones(3)) conv(diff(diff(search_traj[:,3,ns])),ones(3))]
            inds = any(abs(D) .> 1e-3,2)
            rettraj[ns] = search_traj[:,:,ns]
        end
        return rettraj
    end
    for ns = 1:Nsensors
        # return search_traj
        N = size(search_traj,1)
        Tv = Convex.Variable(size(search_traj,1,2))
        diffs = Tv-search_traj[:,:,ns]
        f = Convex.sumsquares(diffs)

        for i = 1:N-2
            Di = abs(Tv[i,:] -2Tv[i+1,:] + Tv[i+2,:])
            # D[i,:] = [Di[1],Di[2],Di[3]]
            # f += 100*max(Di[1],Di[2],Di[3])
            f += 1*sum(Di)
        end
        p = Convex.minimize(f,[Convex.norm(diffs,Inf) < 0.05])
        Convex.solve!(p,SCS.SCSSolver(verbose=false, max_iters=100000, normalize=false, eps = 1e-5))
        if p.status != :Optimal
            warn("Optimization did not converge in first step, using full search_traj, status: ", p.status)
            rettraj[ns] = search_traj[:,:,ns]
            continue
        end
        T = Tv.value
        D = diff(diff(T))
        if false
            Plots.subplot(search_traj[:,:,ns],layout=[1,1,1])
            Plots.subplot!(T,c=:red)
            Plots.subplot(D,layout=[1,1,1])
        end

        diffs = Tv-search_traj[:,:,ns]
        f = Convex.sumsquares(diffs)

        c = Convex.norm(diffs,Inf) < 0.05 # Constrain the approximation error to be no more than this
        nulls = falses(N)
        for i = 2:N-1
            if abs(D[i-1]) < 0.0002
                f += 100*sum(Tv[i-1,:] -2Tv[i,:] + Tv[i+1,:])
                nulls[i] = true
            end
        end
        p = Convex.minimize(f,c)
        Convex.solve!(p,SCS.SCSSolver(verbose=false, max_iters=150000, normalize=false, eps = 1e-5))
        if p.status != :Optimal
            warn("Optimization did not converge in second step, using full search_traj, status: ", p.status)
            rettraj[ns] = search_traj[:,:,ns]
            continue
        end
        T = Tv.value
        rettraj[ns] = T[!nulls,:]

        if false
            Plots.subplot(search_traj[:,:,ns],layout=[1,1,1],c=:blue,lab="search_traj")
            Plots.subplot!(T,c=:red,lab="optresult")
            # Plots.subplot(D,layout=[1,1,1])

            Plots.subplot!(find(!nulls),rettraj[ns],c=:green,lab="rettraj")
            gui()
        end
    end
    return rettraj
end

vars(vs, _) = vs
vars(vs, s::Symbol) = isdefined(s) ? vs : push!(vs, s)
function vars(vs, e::Expr)
    for arg in e.args
        vars(vs, arg)
    end
    vs
end
extractvars(ex::Expr) = vars(Set{Symbol}(), ex)

# """
# Write formulas on the form
#
# `f = response + (var1 .== val1) - (var1 .== val2) ~ rhs`
# """

function dglm(f::Formula,  df::DataFrame, args...)
    l = f.lhs.args
    # Assert operator is minus
    @assert l[1] == :-
    response    = l[2].args[2]
    diffsym = symbol("diff_",response)

    cond1       = l[2].args[3].args
    cond1symb   = cond1[1]
    comparison1 = cond1[2]
    value1      = cond1[3]

    cond2       = l[3].args
    cond2symb   = cond2[1]
    comparison2 = cond2[2]
    value2      = cond2[3]
    # Determine the indices
    ind1 = eval(Expr(:comparison, df[cond1symb], comparison1, value1))
    ind2 = eval(Expr(:comparison, df[cond2symb], comparison2, value2))
    @assert (sum(ind1) == sum(ind2)) "The two conditions did not evaluate to true for the same number of elements"
    diff = df[response][ind1] - df[response][ind2]
    dfd = DataFrame()
    dfd[diffsym] = diff
    V = extractvars(f.rhs)
    for sym = names(df)
        # Only add the relevant symbols to the new dataframe
        if sym ∉ V && sym ∉ f.rhs.args # Super weird bug here. Sometimes the first condition evaluates to false, but the if-statement is still entered!!
            continue
        end
        # @bp all(df[sym][ind1] .!= df[sym][ind2])
        @assert all(df[sym][ind1] .== df[sym][ind2]) "The entries considered for modeling does not have equal values in the two selected datasets"
        dfd[sym] = df[sym][ind1]
    end
    glm($(diffsym) ~ $(f.rhs), dfd, args...)

end

using PyCall
@pyimport matplotlib2tikz
function savetikz(path, fig = PyPlot.gcf(), extra=[""])
    matplotlib2tikz.save(path,fig, figureheight = "\\figureheight", figurewidth = "\\figurewidth", extra = pybuiltin("set")(extra))
end

end

ext: jl


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    boundingbox(geometry)

Axis-aligned bounding box of the `geometry`.
"""
boundingbox(geometry) = boundingbox(vertices(geometry))

# fallback implementation treats geometry as a set of points
function boundingbox(points::AbstractVector{Point{Dim,T}}) where {Dim,T}
  xmin = MVector(ntuple(i->typemax(T), Dim))
  xmax = MVector(ntuple(i->typemin(T), Dim))
  for p in points
    x = coordinates(p)
    @. xmin = min(x, xmin)
    @. xmax = max(x, xmax)
  end
  Box(Point(xmin), Point(xmax))
end

# -----------
# PRIMITIVES
# -----------

boundingbox(p::Primitive) = boundingbox(boundary(p))

boundingbox(b::Box) = b

function boundingbox(s::Sphere{Dim,T}) where {Dim,T}
  c = center(s)
  r = radius(s)
  r⃗ = Vec(ntuple(i->r, Dim))
  Box(c - r⃗, c + r⃗)
end


text:
using SparseArrays

@testset "Matrix" begin
  R = FlintZZ
  M = SMatSpace(R, 3, 3)

  @test R == @inferred base_ring(M)

  A = identity_matrix(FlintZZ, 3)
  Asparse = sparse_matrix(A)

  @test R == @inferred base_ring(Asparse)
  @test M == @inferred parent(Asparse)

  @test 3 == @inferred nrows(Asparse)
  @test 3 == @inferred ncols(Asparse)
  @test 3 == @inferred nnz(Asparse)
  @test (3, 3) == @inferred size(Asparse)

  B = identity_matrix(FlintZZ, 3)
  Bsparse = sparse_matrix(B)

  Csparse = sparse_matrix(FlintZZ)
  @test base_ring(Csparse) == FlintZZ
  @test nrows(Csparse) == 0
  @test ncols(Csparse) == 0

  A = fmpz[1 2; 0 0; 3 4]

  Dsparse = sparse_matrix(A)
  @test base_ring(Dsparse) == FlintZZ
  @test nrows(Dsparse) == 3
  @test ncols(Dsparse) == 2
  @test nnz(Dsparse) == 4

  Esparse = sparse_matrix(FlintZZ)
  @test base_ring(Esparse) == FlintZZ
  @test nrows(Esparse) == 0
  @test ncols(Esparse) == 0
  @test nnz(Esparse) == 0

  # Hash

  @test hash(Asparse) == @inferred hash(Bsparse)

  # Equality

  @test @inferred ==(sparse_matrix(FlintZZ, [0 0 0; 1 1 1]), sparse_matrix(FlintZZ, [0 0 0; 1 1 1]))

  B = zero_matrix(FlintZZ, 2, 3)
  Bsparse = sparse_matrix(B)
  @test nrows(Bsparse) == 2
  @test ncols(Bsparse) == 3
  @test nnz(Bsparse) == 0
  @test Asparse !== Bsparse

  # Sparsity

  @test isapprox(2/3, @inferred sparsity(Asparse))

  # Density

  @test isapprox(1/3, @inferred density(Asparse))

  # Copy

  Csparse = @inferred copy(Asparse)
  @test Csparse == Asparse

  # Entry access

  A = identity_matrix(FlintZZ, 3)
  Asparse = sparse_matrix(A)
  for i in 1:3
    for j in 1:3
      b = @inferred A[i, j]
      if i == j
        @test b == one(FlintZZ)
      else
        @test iszero(b)
      end
    end
  end

  # Row access

  D = sparse_matrix(FlintZZ, [1 5 3; 0 0 0; 0 1 0])
  @test sparse_row(FlintZZ, [1, 2, 3], fmpz[1, 5, 3]) == @inferred D[1]
  @test sparse_row(FlintZZ) == @inferred D[2]
  @test sparse_row(FlintZZ, [2], fmpz[1]) == @inferred D[3]

  D[1] = sparse_row(FlintZZ, [1], fmpz[1])
  D[2] = sparse_row(FlintZZ, [2], fmpz[1])
  D[3] = sparse_row(FlintZZ, [3], fmpz[1])
  @test D == sparse_matrix(identity_matrix(FlintZZ, 3))

  # Modular reduction

  D = sparse_matrix(FlintZZ, [1 5 3; 5 5 5; -4 1 1])
  D = mod_sym!(D, fmpz(5))
  @test nrows(D) == 3
  @test ncols(D) == 3
  @test D.rows[1] == sparse_row(FlintZZ, [1, 3], fmpz[1, -2])
  @test D.rows[2] == sparse_row(FlintZZ)
  @test D.rows[3] == sparse_row(FlintZZ, [1, 2, 3], fmpz[1, 1, 1])

  # Random row

  D = sparse_matrix(FlintZZ, [1 5 3; 5 5 5; -4 1 1])
  r = @inferred rand_row(D)
  @test any(isequal(r), D.rows)

  # Change of ring

  R = ResidueRing(FlintZZ, 5)
  D = sparse_matrix(FlintZZ, [1 5 3; 5 5 5; -4 1 1])
  D_R = @inferred change_base_ring(R, D)
  @test D_R == sparse_matrix(R, map(R, [1 0 3; 0 0 0; 1 1 1]))

  # Transpose

  D = sparse_matrix(FlintZZ, [1 5 3; 0 0 0; 0 1 0])
  Dt = @inferred transpose(D)
  @test Dt == sparse_matrix(FlintZZ, [1 0 0; 5 0 1; 3 0 0])
  @test Dt == D'

  # Iterator interface

  D = sparse_matrix(FlintZZ, [1 5 3; 0 0 0; 0 1 0])
  @test D.rows == [ r for r in D]
  @test D.rows == @inferred collect(D)

  # Multiplications

  D = sparse_matrix(FlintZZ, [1 5 3; 0 0 0; 0 1 0])
  v = fmpz[1, 2, 3]
  w = @inferred mul(D, v)
  @test w == fmpz[20, 0, 2]
  w = @inferred mul(D, view(v, 1:3))
  @test w == fmpz[20, 0, 2]

  v = fmpz[1 2 3; 0 0 4; 0 0 0]
  w = @inferred mul(D, v)
  @test w == fmpz[1 2 23; 0 0 0; 0 0 4]
  v = fmpz[1 1 1; 1 2 3; 0 0 4; 0 0 0]
  w = @inferred mul(D, view(v, 2:4, :))
  @test w == fmpz[1 2 23; 0 0 0; 0 0 4]

  v = matrix(FlintZZ, fmpz[1 2 3; 0 0 4; 0 0 0])
  w = @inferred mul(D, v)
  @test w == matrix(FlintZZ, fmpz[1 2 23; 0 0 0; 0 0 4])

  v = sparse_row(FlintZZ, [2], fmpz[1])
  w = @inferred mul(v, D)
  @test w == sparse_row(FlintZZ)

  # Addition

  D = sparse_matrix(FlintZZ, [1 5 3; 0 0 0; 0 1 0])
  E = sparse_matrix(FlintZZ, [0 0 0; 1 0 1; 0 0 1])
  DplusE = @inferred D + E
  @test DplusE == sparse_matrix(FlintZZ, [1 5 3; 1 0 1; 0 1 1])
  @test D + E == E + D

  # Subtraction

  D = sparse_matrix(FlintZZ, [1 5 3; 0 1 0; 0 1 0])
  E = sparse_matrix(FlintZZ, [1 5 0; 1 0 1; 0 0 0])
  DminusE = @inferred D - E
  @test DminusE == sparse_matrix(FlintZZ, [0 0 3; -1 1 -1; 0 1 0])

  # Scalar multiplication

  D = sparse_matrix(FlintZZ, [1 5 3; 0 0 0; 0 1 0])
  E = @inferred 0 * D
  @test E == zero_matrix(SMat, FlintZZ, 3)
  E = @inferred BigInt(2) * D
  @test E == sparse_matrix(FlintZZ, [2 10 6; 0 0 0; 0 2 0])
  E = @inferred fmpz(2) * D
  @test E == sparse_matrix(FlintZZ, [2 10 6; 0 0 0; 0 2 0])

  R = ResidueRing(FlintZZ, 6)
  D = sparse_matrix(R, [1 2 2; 0 0 1; 2 2 2])
  E = @inferred fmpz(3) * D
  @test E == sparse_matrix(R, [3 0 0; 0 0 3; 0 0 0])
  E = @inferred Int(3) * D
  @test E == sparse_matrix(R, [3 0 0; 0 0 3; 0 0 0])
  E = @inferred R(3) * D
  @test E == sparse_matrix(R, [3 0 0; 0 0 3; 0 0 0])

  # Submatrix

  D = sparse_matrix(FlintZZ, [1 5 3; 0 0 0; 0 1 0])
  E = @inferred sub(D, 2:3, 2:3)
  @test base_ring(E) == FlintZZ
  @test nrows(E) == 2
  @test ncols(E) == 2
  @test E == sparse_matrix(FlintZZ, [0 0; 1 0])

  # Vertical concatenation
  D = sparse_matrix(FlintZZ, [1 5 3; 0 0 0; 0 1 0])
  E = @inferred vcat(D, D)
  @test base_ring(E) == FlintZZ
  @test nrows(E) == 6
  @test ncols(E) == 3
  @test E == sparse_matrix(FlintZZ, [1 5 3; 0 0 0; 0 1 0;
                                     1 5 3; 0 0 0; 0 1 0])
  # Horizontal concatenation

  D = sparse_matrix(FlintZZ, [1 5 3; 0 0 0; 0 1 0])
  E = @inferred hcat(D, D)
  @test base_ring(E) == FlintZZ
  @test nrows(E) == 3
  @test ncols(E) == 6
  @test E == sparse_matrix(FlintZZ, [1 5 3 1 5 3; 0 0 0 0 0 0; 0 1 0 0 1 0])

  # Hadamard bound

  D = sparse_matrix(FlintZZ, [1 5 3; 0 1 0; 0 1 0])
  h = @inferred hadamard_bound2(D)
  @test h == fmpz(35)

  # Maximum/minium
  D = sparse_matrix(FlintZZ, [1 5 3; 0 -10 0; 0 1 0])
  b = @inferred maximum(abs, D)
  @test b == fmpz(10)
  D = sparse_matrix(FlintQQ, [1 2 fmpq(9, 4); 0 -10 0; 0 1 0])
  b = @inferred maximum(D)
  @test b == fmpq(9, 4)
  b = @inferred minimum(D)
  @test b == fmpq(-10)

  D = sparse_matrix(FlintZZ, [0 2 0; 0 0 1; 0 0 0])
  @test @inferred isupper_triangular(D)
  D = sparse_matrix(FlintZZ, [0 0 2; 0 0 1; 0 0 0])
  @test !isupper_triangular(D)
  D = sparse_matrix(FlintZZ, [0 0 0; 0 0 0; 0 0 0])
  @test @inferred isupper_triangular(D)

  # Zero and identity matrix

  D = @inferred identity_matrix(SMat, FlintZZ, 3)
  @test D == sparse_matrix(FlintZZ, [1 0 0; 0 1 0; 0 0 1]);
  D = @inferred zero_matrix(SMat, FlintZZ, 3)
  @test D == sparse_matrix(FlintZZ, [0 0 0; 0 0 0; 0 0 0]);
  D = @inferred zero_matrix(SMat, FlintZZ, 4, 3)
  @test D == sparse_matrix(FlintZZ, [0 0 0; 0 0 0; 0 0 0; 0 0 0]);

  # Concatenation syntax

  D = sparse_matrix(FlintZZ, [1 5 3; 0 -10 0; 0 1 0])
  E = [D D D]
  @test E == hcat(hcat(D, D), D)
  E = [D; D; D]
  @test E == vcat(vcat(D, D), D)

  # Is one

  D = sparse_matrix(FlintZZ, [1 5 3; 0 -10 0; 0 1 0])
  b = @inferred isone(D)
  @test !b
  D = sparse_matrix(FlintZZ, [1 0 0; 0 1 0; 0 0 1])
  b = @inferred isone(D)
  @test b
  D = sparse_matrix(FlintZZ, [1 0 0; 0 0 1; 0 1 0])
  b = @inferred isone(D)
  @test !b
  D = sparse_matrix(FlintZZ, [1 0 0; 0 1 0; 0 0 1; 0 0 0])
  b = @inferred isone(D)
  @test !b

  D = sparse_matrix(FlintZZ, [0 0 0; 0 0 0])
  b = @inferred iszero(D)
  @test b
  D = sparse_matrix(FlintZZ, [0 0 0; 0 0 1])
  b = @inferred iszero(D)
  @test !b
  D = sparse_matrix(FlintZZ, [0 0 1; 0 0 0])
  b = @inferred iszero(D)
  @test !b

  # Conversion to julia types
  D = sparse_matrix(FlintZZ, [1 5 3; 0 -10 0; 0 1 0])
  E = SparseArrays.sparse(D)
  @test Matrix(E) == fmpz[1 5 3; 0 -10 0; 0 1 0]
  @test Array(E) == fmpz[1 5 3; 0 -10 0; 0 1 0]
end


text:
export Residual, get_volume_residual_function, get_surface_residual_function

using FastGaussQuadrature

"""
Function to obtain Jacobian from `Finel` struct
"""
function from_finel_jac!(Jac::Matrix{Float64}, fin::Finel, x...)
    deriv_order=zeros(Int64, length(x))
    for j=1:length(x)
        deriv_order[j]=1

        f=get_interpfun(fin.domain.interp; deriv_order=deriv_order)

        for (i, coefs) in enumerate(fin._coord_linpars)
            Jac[i, j]=f(x..., coefs...)
        end

        deriv_order[j]=0
    end
end

"""
Function to return polynomial coefficients from array notation
"""
function get_polycoefs(vars::Tuple{Vector{Float64}, InterpolationFunction}...)
    [
        regress(finel_get_regr(fin, intp), vals) for (vals, intp) in vars
    ]
end

"""
Definition for volume residual obtention configurations struct

* `variables`: vector of tuples. First element indicates a symbol for a variable; the second, an `InterpolationFunction`
    struct to serve as its interpolator. The third should be a vector of tuples, each with a symbol and a vecor of ints
    indicating derivative order in each axis

    Example:
        ```[
            (
                :u, interp1, [
                    (:du!dx, [1, 0]), (du!dy, [0, 1])
                ]
            ),
            (
                :v, interp2, []
            )
        ]```
* `ex`: expression evaluating the residuals.

    Examples:
        * `:(u*du!dx+v-1.0)`
        * `begin w=sqrt(u); w*v-du!dx end`
* `extra_args`: symbols for any additional arguments to be passed
"""
mutable struct Residual
    variables::Vector{
        Tuple{
            Symbol, InterpolationFunction, Vector{
                Tuple{Symbol, Vector{Int64}}
            }
        }
    }
    ex::Expr
    _nvr::Int64
    _weights::Vector{Float64}
    _abscissas::Vector{Float64}
    extra_args::Vector{Symbol}
end

"""
Constructor for Residual taking only public values

* `variables`: vector of tuples. First element indicates a symbol for a variable; the second, an `InterpolationFunction`
    struct to serve as its interpolator. The third should be a vector of tuples, each with a symbol and a vecor of ints
    indicating derivative order in each axis

    Example:
        ```[
            (
                :u, interp1, [
                    (:du!dx, [1, 0]), (du!dy, [0, 1])
                ]
            ),
            (
                :v, interp2, []
            )
        ]```
* `ex`: expression evaluating the residuals.

    Examples:
        * `:(u*du!dx+v-1.0)`
        * `begin w=sqrt(u); w*v-du!dx end`

* `nquad`: Gaussian quadrature order (Gauss Legendre)
* `extra_args`: symbols for any extra arguments to be passed for residual computation
"""
function Residual(
    variables,
    ex::Expr,
    nquad::Int64;
    extra_args=[]
)
    _nvr=length(variables)

    for (_, _, vv) in variables
        _nvr+=length(vv)
    end

    x, w=gausslegendre(nquad)
    x=(1.0.+x)./2
    w./=2

    Residual(
        variables,
        ex,
        _nvr,
        w,
        x,
        extra_args
    )
end

"""
Function to obtain function for punctual evaluation volume residual

* `resd`: an `Residual` instance

* return: a function (in the format `resfun(v1, dv1!dx, dv2!dy, v2...)` evaluating the residual given other variables)
"""
function get_punctual_resfunc(resd::Residual)
    pronames=Vector{Symbol}()

    for v in resd.variables
        push!(pronames, v[1])

        for ds in v[3]
            push!(pronames, ds[1])
        end
    end

    q=quote
        bfun=($(pronames...), $(resd.extra_args...)) -> ( $(resd.ex) )
    end

    try
        eval(q)

        return bfun
    catch
        throw(
            error(
                "get_punctual_resfunc:ERROR:unable to compile punctual residual function for expression $(resd.ex)"
            )
        )
    end
end

"""
Transform vectors of values into vectors of absolute coordinate polynomial coefficients
"""
function vals_to_polys(fin::Finel, resd::Residual, vals::Vector{Vector{Float64}})
    return [
        regress(finel_get_regr(fin, vintp[2]), vs) for (vs, vintp) in zip(vals, resd.variables)
    ]
end

"""
Obtains variables of interest given absolute coordinates and polynomial coefficients
"""
function polys_to_vars(resd::Residual, pols::Vector{Vector{Float64}}, pos::Vector{Float64})
    vs=Vector{Float64}(undef, resd._nvr)

    nv=1
    for (v, p) in zip(resd.variables, pols)
        f=get_interpfun(v[2])

        vs[nv]=f(pos..., p...)
        nv+=1

        for dv in v[3]
            f=get_interpfun(v[2]; deriv_order=dv[2])

            vs[nv]=f(pos..., p...)
            nv+=1
        end
    end

    vs
end

"""
Function to produce local coordinates given percentages in the space between restrictions
"""
function loc_coordinate_inrestr!(dmn::Domain, percs::Vector{Float64}, coords::Vector{Float64}; exempt::Union{Int64, Nothing}=nothing)
    w=1.0

    for (i, p) in enumerate(percs)
        br, fr=dmn._restriction_funcs[i][1](coords...), dmn._restriction_funcs[i][2](coords...)

        if i!=exempt
            w*=(fr-br)
        end

        coords[i]=(1.0-p)*br+p*fr
    end

    w
end

"""
Get volume evaluation of residual
"""
function get_volume_residual_function(resd::Residual)
    locfunc=get_punctual_resfunc(resd)
    nquad=length(resd._abscissas)

    volres=(vvals::Vector{Vector{Float64}}, fin::Finel, args...) -> begin
        naux=length(fin.domain.interp.coords)

        pols=vals_to_polys(fin, resd, vvals)
        rngs=[1:nquad for i=1:naux]
        locpercs=Vector{Float64}(undef, naux)
        locpos=Vector{Float64}(undef, naux)

        nd=length(fin._coord_linpars)

        Jac=Matrix{Float64}(undef, nd, naux)
        pos=Vector{Float64}(undef, nd)

        R=0.0

        for cinds in Iterators.product(rngs...)
            w=1.0

            for (i, c) in enumerate(cinds)
                w*=resd._weights[c]

                locpercs[i]=resd._abscissas[c]
            end

            w*=loc_coordinate_inrestr!(fin.domain, locpercs, locpos)

            get_abs_coords!(pos, fin, locpos)

            from_finel_jac!(Jac, fin, locpos...)

            w*=orthonormal_det!(Jac)

            R+=locfunc(polys_to_vars(resd, pols, pos)..., args...)*w
        end

        R
    end

    volres
end

"""
Get surface evaluation of residual
"""
function get_surface_residual_function(resd::Residual) # havent changed anything yet
    locfunc=get_punctual_resfunc(resd)
    nquad=length(resd._abscissas)

    surfres=(vvals::Vector{Vector{Float64}}, fin::Finel, args...) -> begin
        naux=length(fin.domain.interp.coords)

        pols=vals_to_polys(fin, resd, vvals)
        rngs=[1:nquad for i=1:(naux-1)]
        locpercs=Vector{Float64}(undef, naux)
        locpos=Vector{Float64}(undef, naux)

        nd=length(fin._coord_linpars)

        Jac=Matrix{Float64}(undef, nd, naux)
        pos=Vector{Float64}(undef, nd)

        R=0.0

        for dim=1:length(fin.domain.interp.coords)
            for nb=1:2
                for cinds in Iterators.product(rngs...)
                    w=1.0

                    j=1
                    for i=1:length(fin.domain.interp.coords)
                        if i!=dim
                            c=cinds[j]

                            w*=resd._weights[c]

                            locpercs[i]=resd._abscissas[c]

                            j+=1
                        else
                            locpercs[i]=(
                                nb==1 ?
                                    0.0 : 1.0
                            )
                        end
                    end

                    w*=loc_coordinate_inrestr!(fin.domain, locpercs, locpos; exempt=dim)

                    get_abs_coords!(pos, fin, locpos)

                    from_finel_jac!(Jac, fin, locpos...)

                    restr_grad=fin.domain._restriction_derivatives[dim][nb](locpos...)

                    locvals=locfunc(polys_to_vars(resd, pols, pos)..., args...)
                    shat=normal_volspace!(Jac, dim, restr_grad, nb==1)

                    R+=dot(locvals,
                        shat)*w
                end
            end
        end

        R
    end

    surfres
end