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

using Markdown
using InteractiveUtils

# ╔═╡ 7ed39fbc-8c1a-4d10-a6b7-7608e6e87395
md"""
### Estimation of flux through a metabolic network

Our overall objective as a metabolic engineer is to maximize the `performance` of a metabolic network, e.g., develop a system that produces a desired product with the highest possible yield, or at the maximum possible rate, etc.

Toward this goal, we developed material balances that describe the dynamic and steady-state `concentration` of chemical species (metabolites) in a metabolic network. In particular, the `concentration` of metabolites (written in cell mass-specific units, e.g., $\star$mol/gDW) for a metabolic network with $\mathcal{R}$ reactions and $\mathcal{M}$ metabolites operating in batch culture can be written (in matrix-vector form) as:

$$\frac{d\mathbf{x}}{dt} = \mathbf{S}\mathbf{v} - \mu\mathbf{x}$$

where $\mathbf{x}$ denbotes the $\mathcal{M}\times{1}$ column vector of metabolite concentrations, $\mathbf{S}$ denotes the $\mathcal{M}\times\mathcal{R}$ stoichiometricx matrix, $\mathbf{v}$ denotes the $\mathcal{R}\times{1}$ column vector of (net)reaction rates, and $\mu\mathbf{x}$ denotes the diliution due to growth term. Note: the form of the dilution terms depends upon the concentration basis.

We also introduced some tools that we might use to achieve our engineering objective. In particular, we introduced some tools to explore the structure and properties of the stoichiometric matrix $\mathbf{S}$, and Flux Balance Analysis (FBA), a tool to compute optimal flux through a metabolic network.

Today, we will:

* Compute [extreme pathways](https://pubmed.ncbi.nlm.nih.gov/10716907/) using the [expa](https://pubmed.ncbi.nlm.nih.gov/15613397/) algorithm
* Compute optimal `flux` through a metabolic network using [flux balance analysis (FBA)](https://www.ncbi.nlm.nih.gov/labs/pmc/articles/PMC3108565/)

"""

# ╔═╡ a55ad25a-e4ff-4af7-9684-c1b02fec900d
md"""
### Extreme pathways
The set of extreme pathways, a generating set for all possible steady-state flux distributions in a metabolic network, can be computed from the stoichiometric matrix $\mathbf{S}$. In particular, any steady-state flux through a metabolic network can be written as:

$$\mathbf{v} = \sum_{i=1}^{\mathcal{P}}\alpha_{i}\mathbf{P}_{i}$$

where $\mathcal{P}$ denotes the number of convex basis pathways, and $\mathbf{P}_{i}$ denotes the $\mathcal{R}\times{1}$ extreme pathway basis vector. The term(s) $\alpha_{i}\geq{0}$ denotes the weight of extreme pathway $i$.

For more information on extreme pathways, check out:
* [Lecture 12: Pathways. Systems Biology: Constraint-based Reconstruction and Analysis, Bernhard Ø. Palsson, Cambridge University Press, 2015](https://www.youtube.com/watch?v=5o46XASlNDw)
"""

# ╔═╡ dc712cba-1f18-46ec-a7f1-b25b0aca1e0a
md"""
##### Computing extreme pathways using expa
We've implemented the expa algorithm from Palsson and coworkers:

* [Bell SL, Palsson BØ. Expa: a program for calculating extreme pathways in biochemical reaction networks. Bioinformatics. 2005 Apr 15;21(8):1739-40. doi: 10.1093/bioinformatics/bti228. Epub 2004 Dec 21. PMID: 15613397.](https://pubmed.ncbi.nlm.nih.gov/15613397/)

Let's use expa to explore the toy network taken from (note: we've renumbered the reactions):

* [Schilling CH, Letscher D, Palsson BO. Theory for the systemic definition of metabolic pathways and their use in interpreting metabolic function from a pathway-oriented perspective. J Theor Biol. 2000 Apr 7;203(3):229-48. doi: 10.1006/jtbi.2000.1073. PMID: 10716907.](https://pubmed.ncbi.nlm.nih.gov/10716907/)
"""

# ╔═╡ 2b9ca1ce-9080-4c42-91a1-e5005ddc4ee5
pidx = 7

# ╔═╡ 4b32c7ae-397f-43c4-ac43-d3107b4de9c4
md"""
__Table 1__: Extreme pathways table for the toy network of [Schilling CH, Letscher D, Palsson BO. Theory for the systemic definition of metabolic pathways and their use in interpreting metabolic function from a pathway-oriented perspective. J Theor Biol. 2000 Apr 7;203(3):229-48. doi: 10.1006/jtbi.2000.1073. PMID: 10716907.](https://pubmed.ncbi.nlm.nih.gov/10716907/)
"""

# ╔═╡ 009df867-8c26-422a-8adf-a83da0b667bd
md"""
### Flux Balance Analysis (FBA)
Flux balance analysis (FBA) is a mathematical modeling and analysis approach that estimates the _intracellular_ reaction rate (metabolic flux) of carbon and energy throughout a metabolic network (units: $\star$mol/gDW-time or $\star$mol/L-time for cell-free networks). FBA arguably the most widely used tool in this space. However, there are alternatives to FBA, such as metabolic flux analysis (MFA), but these alternatives vary more in the solution approach than the structure of the estimation problem.

Let's look at the following reference to understand better the different components of a flux balance analysis problem:

* [Orth, J., Thiele, I. & Palsson, B. What is flux balance analysis?. Nat Biotechnol 28, 245–248 (2010). https://doi.org/10.1038/nbt.1614](https://www.ncbi.nlm.nih.gov/labs/pmc/articles/PMC3108565/)

Computational tools for working with flux balance analysis models:

* [Heirendt, Laurent et al. "Creation and analysis of biochemical constraint-based models using the COBRA Toolbox v.3.0." Nature protocols vol. 14,3 (2019): 639-702. doi:10.1038/s41596-018-0098-2](https://pubmed.ncbi.nlm.nih.gov/30787451/)

"""

# ╔═╡ 3f98b575-eb2c-40b3-ac6f-43c55e4ec9d3
md"""
##### Flux balance analysis problem structure
The FBA problem is typically encoded as a linear programming (LP) problem of the form:

$$\max_{v}\sum_{i=1}^{\mathcal{R}}c_{i}v_{i}$$

subject to the constraints:

$$\begin{eqnarray}
\mathbf{S}\mathbf{v} &=& \mathbf{0}\\
\mathcal{L}\leq&\mathbf{v}&\leq\mathcal{U}\\
\dots &=& \dots
\end{eqnarray}$$

where $\mathbf{S}$ denotes the stoichiometric matrix, $c_{i}$ denote the objective coefficients,
$\mathbf{v}$ denotes the metabolic flux (the unknown that we are trying to estimate), and
$\mathcal{L}$ ($\mathcal{U}$) denote the permissible lower (upper) bounds on the _unkown_ metabolic flux. The first set of constraints enforces conservation of mass, while the second imparts thermodynamic and kinetic information into the calculation. Finally, there are potentially other types of constraints (both linear and nonlinear) used in this type of problem (we will not cover these here, but these additional constraints may be important in specific applications that we see later).
"""

# ╔═╡ 72e5c813-2545-4d7c-9437-dc3f98f4e6e7
md"""
##### Toy FBA example
Let's consider the `toy` network above. Assume this network is operating in cells growing in a continuous culture at a steady-state, where we use cell mass specific units.  We'll use the [GNU Linear Programming Kit (GLPK)](https://www.gnu.org/software/glpk/) and the [GLPK.jl](https://github.com/jump-dev/GLPK.jl) package to solve the linear programming problem.
"""

# ╔═╡ 6dd988f9-28d3-47cb-8b32-a415132b388b
md"""
__Table 2__: Flux table for maximizing toy network. The [GNU Linear Programming Kit solver](https://www.gnu.org/software/glpk/) produced an estimate of metabolic flux ($\star$mol/gDW-hr) given the problem constraints.
"""

# ╔═╡ a40ef7cf-3e46-4f4b-9af9-5a0c88d15ac9
md"""
### Some fun with constraints and the FBA problem formulation
The standard flux balance analysis problem is written in `concentration` units and metabolic reaction flux, e.g., $\star$mol/gDW-hr. However, there is nothing that says we have to do that. For example, it may be more convenient to work in mass or mole units instead of working in concentration units. Let's consider the latter.

We know from our material and energy balance class that the open mole balance around component $i$ in the _logical_ control volume is given by:

$$\sum_{s=1}^{\mathcal{S}}d_{s}\dot{n}_{is} + \sum_{j=1}^{\mathcal{R}}\sigma_{ij}\dot{\epsilon}_{j} = 0\qquad{i=1,2,\dots,\mathcal{M}}$$

The first term is the rate of `transport` into and from the control volume (units: mol $i$/time) from $\mathcal{S}$ possible `streams`; the `transport` terms can be physical or logical where $d_{s} = 1$ if the stream $s$ enters control volume, $d_{s}=-1$ is stream $s$ exits the control volume. The second summation denotes the reaction terms, where $\sigma_{ij}$ denote the stoichiometric coefficient describing the connection between metabolite $i$ and reaction $j$ and $\dot{\epsilon}_{j}$ denote the open extent of reaction (units: mol/time).

Let's suppose that we have a single logical stream entering (s=1) and exiting (s=2) the control volume. In this case, the open mole balance becomes:

$$\dot{n}_{i,2} = \dot{n}_{i,1} + \sum_{j=1}^{\mathcal{R}}\sigma_{ij}\dot{\epsilon}_{j}\qquad{i=1,2,\dots,\mathcal{M}}$$

These balances can be used as constraints to find the optimal open extent of reaction. In particular, we know that we actually pass $\alpha\leq{S\dot{\epsilon}}\leq\beta$ to the solve. Thus, because $\dot{n}_{i,2}\geq{0}$, the FBA problem is subject to the mol constraints:

$$\dot{n}_{i,1} + \sum_{j=1}^{\mathcal{R}}\sigma_{ij}\dot{\epsilon}_{j}\geq{0}\qquad\forall{i}$$

In other words, when searching for the optimal set of $\dot{\epsilon}_{j}$ we have to select values that give physically realistic answers (we can't have a negative mol flow rate). Next, the $\dot{\epsilon}_{j}$ terms (just flux) are bounded from above and below:

$$\mathcal{L}_{j}\leq\dot{\epsilon}_{j}\leq\mathcal{U}_{j}\qquad{j=1,2\dots,\mathcal{R}}$$

where the $\mathcal{L}_{j}$ and $\mathcal{U}_{j}$ denote the lower and upper bounds that $\dot{\epsilon}_{j}$ can take, remember that the open extents $\dot{\epsilon}_{j}$ are just reaction rates times the volume. Thus, the lower and upper bounds describe the permissible range we expect the rate _could_ obtain.

Putting everything together gives a slightly different problem formulation to compute the mol/time flux through a reaction network. An objective:

$$\mathcal{O} = \sum_{j=1}^{\mathcal{R}}c_{j}\dot{\epsilon}_{j}$$

is minimized (or maximized) subject to a collection of linear constraints:

$$\begin{eqnarray}
\sum_{j=1}^{\mathcal{R}}\sigma_{ij}\dot{\epsilon}_{j}&\geq&{-\dot{n}_{i,1}}\qquad\forall{i}\\
\dot{n}_{i,2}&\geq&{0}\qquad\forall{i}\\
\mathcal{L}_{j}&\leq\dot{\epsilon}_{j}\leq&\mathcal{U}_{j}\qquad{j=1,2\dots,\mathcal{R}}
\end{eqnarray}$$

"""

# ╔═╡ a7dd770b-22b3-4561-a85b-6071e7142c6d
md"""
##### Example: flux balance analysis with mole flow rates and reaction open extent
Let's suppose the objective was to maximize the production of $P_{x}$. Let's set up a flux balance analysis problem to solve this problem.
"""

# ╔═╡ 0fd5a620-682e-48d8-b1e6-ba110018f9e3
md"""
__Table 3__: Flux table for maximizing $P$ production. The [GNU Linear Programming Kit solver](https://www.gnu.org/software/glpk/) produced an estimate of the (open) reaction extents $\dot{\epsilon}_{j}~\forall{j}$ given the problem constraints.
"""

# ╔═╡ bc52eb4b-9761-459c-8a6e-c457f13d7335
md"""
__Table 4__: State table for maximizing $P$ production. The flux balance analysis problem produces estimates for the optimal value of the (open) reaction extents $\dot{\epsilon}_{j}~\forall{j}$ (units: $\star$mol/hr). We can then use the steady-state species mol balances to compute the output composition (stream 1 input, stream 2 is an output)
"""

# ╔═╡ b620232d-6c8f-478c-93e5-9b4743518ff6
md"""
### Summary and conclusions
In this lecture we:

* Computed [extreme pathways](https://pubmed.ncbi.nlm.nih.gov/10716907/) using the [expa](https://pubmed.ncbi.nlm.nih.gov/15613397/) algorithm
* Computed the optimal `flux` through a metabolic network using [flux balance analysis (FBA)](https://www.ncbi.nlm.nih.gov/labs/pmc/articles/PMC3108565/)

"""

# ╔═╡ ebb0d7cf-fd92-45be-ab98-5bb4de2c2ba2
md"""
### Next time:

* Compare a dynamic experiment with a flux balance analysis computation: [Varma A, Palsson BO. Stoichiometric flux balance models quantitatively predict growth and metabolic by-product secretion in wild-type Escherichia coli W3110. Appl Environ Microbiol. 1994 Oct;60(10):3724-31. doi: 10.1128/aem.60.10.3724-3731.1994. PMID: 7986045; PMCID: PMC201879](https://pubmed.ncbi.nlm.nih.gov/7986045/)

* Flux problem basics: how do we compute the bounds of a reaction?

"""

# ╔═╡ 0e46d5ca-197f-47ca-a075-7281e3087460
TableOfContents(title="📚 Table of Contents", indent=true, depth=5, aside=true)

# ╔═╡ 907af081-4a01-4a51-bda9-3a34295f9208
function ingredients(path::String)

	# this is from the Julia source code (evalfile in base/loading.jl)
	# but with the modification that it returns the module instead of the last object
	name = Symbol("lib")
	m = Module(name)
	Core.eval(m,
        Expr(:toplevel,
             :(eval(x) = $(Expr(:core, :eval))($name, x)),
             :(include(x) = $(Expr(:top, :include))($name, x)),
             :(include(mapexpr::Function, x) = $(Expr(:top, :include))(mapexpr, $name, x)),
             :(include($path))))
	m
end

# ╔═╡ c51f8430-8147-11ec-3501-498c9bf67b57
begin

	# import some packages -
	using PlutoUI
	using LinearAlgebra
	using RowEchelon
	using IterativeSolvers
	using Combinatorics
	using Plots
	using GLPK
	using PrettyTables

	# setup paths -
	const _PATH_TO_NOTEBOOK = pwd()
	const _PATH_TO_DATA = joinpath(_PATH_TO_NOTEBOOK,"data")
	const _PATH_TO_FIGS = joinpath(_PATH_TO_NOTEBOOK,"figs")
	const _PATH_TO_SRC = joinpath(_PATH_TO_NOTEBOOK,"src")

	# load the class lib as LectureLib -
	lib = ingredients(joinpath(_PATH_TO_SRC, "Include.jl"));

	# return -
	nothing
end

# ╔═╡ 306e390f-acbe-4b2a-8e4f-3571003359ad
PlutoUI.LocalResource(joinpath(_PATH_TO_FIGS,"Fig-expa-ToyNetwork.png"))

# ╔═╡ be4b1854-1a14-4452-ad7e-22d614740a10
begin

	# Setup a collection of reaction strings -
	reaction_array = Array{String,1}()

	# Setup a collection of reaction strings -
	reaction_array = Array{String,1}()

	# encode the reactions -
	# internal reactions -
	push!(reaction_array,"v₁,A,B,false")
	push!(reaction_array,"v₂,B,C,true")
	push!(reaction_array,"v₃,C,D,true")
	push!(reaction_array,"v₄,C,E,false")

	# exchange reactions -
	push!(reaction_array,"b₁,∅,A,true")
	push!(reaction_array,"b₂,∅,B,true")
	push!(reaction_array,"b₃,D,∅,true")
	push!(reaction_array,"b₄,E,∅,true")

	# compute the stoichiometric matrix -
	# the optional expand arguement = should we split reversible reactions? (default: false)
	(S, species_array, reaction_name_array) = lib.build_stoichiometric_matrix(reaction_array;
		expand=true);

	# show -
	nothing
end

# ╔═╡ 8965f69d-3014-46b3-816d-7d6f7fd57adf
(ℳ,ℛ) = size(S)

# ╔═╡ d0e0e084-c706-45cf-b40f-47ca2e5ce088
S

# ╔═╡ 8be1e489-7382-4315-8c4c-111abdead290
begin

	# compute the extreme pathways Tableu -
	PM = lib.expa(S)

	# P constaints the extreme pathways (rows) and 𝒩 is the "balanced" array (should be all zeros) -
	P = PM[:,1:ℛ]
	𝒩 = PM[:,(ℛ+1):end]

	# show -
	nothing
end

# ╔═╡ a655928e-a275-4bb1-bca3-39640804ed63
size(P)

# ╔═╡ 9108a028-a5fb-428e-9fd8-2c7f0d80f4af
rank(P)

# ╔═╡ 0a325629-0a14-426f-a99e-9bccdc2e0dfa
begin

	# find non-zero index for reactions -
	idx_non_zero = findall(x->x!=0.0,P[pidx,:]);

	# which reactions are being used?
	reaction_list_expa = reaction_name_array[idx_non_zero]

	# generate the reaction strings -
	reaction_display_string_array = lib.generate_reaction_strings(reaction_array; expand=true)

	# make a table -
	expa_pathway_table = Array{Any,2}(undef,length(idx_non_zero),2)

	for i ∈ 1:length(reaction_list_expa)
		expa_pathway_table[i,1] = reaction_list_expa[i]
		expa_pathway_table[i,2] = reaction_display_string_array[idx_non_zero[i]]
	end

	expa_header_row = (["name", "reaction string"])
	pretty_table(HTML, expa_pathway_table; header=expa_header_row, tf=tf_html_minimalist)

end

# ╔═╡ 45222d77-906b-4393-b2f9-63452300c28e
begin

	# species bounds array -
	# Sv = 0
	# Solver by default takes L <= Sv <= U
	species_bounds_array = zeros(ℳ,2);

	# specify the flux_bounds -
	flux_bounds_array = [

		# v₁ (ℒ,𝒰)
		0.0 100.0 		; # 1 v₁ units: mmol/gDW-hr

		# v₂
		0.0 100.0 		; # 2 Fv₂ units: mmol/gDW-L
		0.0 100.0 		; # 3 Rv₂ units: mmol/gDW-L

		# v₃
		0.0 100.0 		; # 4 Fv₃ units: mmol/gDW-L
		0.0 100.0 		; # 5 Rv₃ units: mmol/gDW-L

		# v₄
		0.0 100.0 		; # 6 v₄ units: mmol/gDW-L

		# b₁
		0.0 100.0 		; # 7 Fb₁ units: mmol/gDW-L
		0.0 0.0 		; # 8 Rb₁ units: mmol/gDW-L

		# b₂
		0.0 100.0 		; # 9 Fb₂ units: mmol/gDW-L
		0.0 100.0 		; # 10 Rb₂ units: mmol/gDW-L

		# b₃
		0.0 100.0 		; # 11 Fb₃ units: mmol/gDW-L
		0.0 0.0 		; # 12 Rb₃ units: mmol/gDW-L

		# b₄
		0.0 100.0 		; # 13 Fb₄ units: mmol/gDW-L
		0.0 0.0 		; # 14 Rb₄ units: mmol/gDW-L
	];

	# specify the objective -
	c = zeros(ℛ)
	c[13] = -1 # negative: default is minimize

	# compute the optimal flux distribution -
	result = lib.flux(S, flux_bounds_array,species_bounds_array,c);

	# show -
	nothing
end

# ╔═╡ 87436aa6-289f-48f7-93ed-865360f27d92
result

# ╔═╡ 39f3f612-114e-4f71-be75-ec427607e577
begin
	# initialize -
	state_table = Array{Any,2}(undef,ℛ,3)

	# compute the reaction string array -
	reaction_string_array = lib.generate_reaction_strings(reaction_array; expand=true)

	# what is v?
	v = result[2]

	# build state table?
	for i ∈ 1:ℛ
		state_table[i,1] = i
		state_table[i,2] = reaction_string_array[i]
		state_table[i,3] = v[i]
	end
	header_row = (["index", "reaction string", "optimal flux"], ["","","mmol/gDW-L"])
	pretty_table(HTML, state_table; header=header_row, tf=tf_html_minimalist)
end

# ╔═╡ 3a091b76-76cd-479c-b382-d5f2e3e75c59
PlutoUI.LocalResource(joinpath(_PATH_TO_FIGS,"Fig-ToyNetwork-CBT_v2.png"))

# ╔═╡ 1f31384d-13f8-4f31-9e23-9f6ad18eda23
begin

	# setup flow rate in
	# units: *mol/time
	n_dot_in = [
		10.0 	; # 1 A₁
		10.0 	; # 2 A₂
		0.0 	; # 3 B
		0.0 	; # 4 P
		0.0 	; # 5 C
		0.0 	; # 6 x
		0.0 	; # 7 y
	]

	# Setup STM -
	# notice: the STM does NOT have the exchange fluxes b*
	S₂ = [

		# r₁ r₂ r₃
		-1.0 0.0 0.0 ; # 1 A₁
		0.0 0.0 -1.0 ; # 2 A₂
		1.0 -1.0 0.0 ; # 3 B
		0.0 1.0 0.0  ; # 4 P
		0.0 0.0 1.0  ; # 5 C
		-1.0 0.0 1.0 ; # 6 x
		1.0 0.0 -1.0 ; # 7 y
	];

	# set the flux bounds array -
	# assume: all reactions are *irreversible* -
	flux_bounds_array₂ = [

		# ℒ 𝒰
		0.0 1000.0 	; # 1 r₁
		0.0 1000.0  	; # 2 r₂
		0.0 20.0 	; # 3 r₃
	];

	# set the species bounds array -
	# units: *mol/time
	species_bounds_array₂ = [

		# ℒ lower     𝒰 upper
		-n_dot_in[1] 1000.0 				; # 1 A₁
		-n_dot_in[2] 1000.0 				; # 2 A₂
		-n_dot_in[3] 1000.0 				; # 3 B
		-n_dot_in[4] 1000.0 				; # 4 P
		-n_dot_in[5] 1000.0 				; # 5 C
		-n_dot_in[6] 1000.0 				; # 6 x
		-n_dot_in[7] 1000.0 				; # 7 y
	];

	# max P formation -
	(ℳ₂,ℛ₂) = size(S₂)
	obj_vector₂ = zeros(ℛ₂)
	obj_vector₂[2] = -1.0

	# compute the optimal flux -
	result₂= lib.flux(S₂, flux_bounds_array₂, species_bounds_array₂, obj_vector₂)

	# initialize some storage -
	flux_table = Array{Any,2}(undef,ℛ₂,4)

	# reactions -
	reaction_strings = [
		"A₁ + x --> B + y" 	; # 1 r₁
		"B --> P" 			; # 2 r₂
		"A₂ + y --> C + x" 	; # 3 r₃
	]

	# populate the state table -
	for reaction_index ∈ 1:ℛ₂
		flux_table[reaction_index,1] = reaction_strings[reaction_index]
		flux_table[reaction_index,2] = result₂.calculated_flux_array[reaction_index]
		flux_table[reaction_index,3] = flux_bounds_array₂[reaction_index,1]
		flux_table[reaction_index,4] = flux_bounds_array₂[reaction_index,2]
	end

	# header row -
	flux_table_header_row = (["Reaction","ϵᵢ_dot", "ϵ₁_dot LB", "ϵ₁_dot UB"],
		["","mol/time", "mol/time", "mol/time"]);

	# write the table -
	pretty_table(HTML, flux_table; header=flux_table_header_row, tf=tf_html_minimalist)
end

# ╔═╡ c1368f87-bd93-448f-ad10-b70040f5c854
begin

	# setup: state table -
	ϵ = result₂.calculated_flux_array
	n_dot_out = n_dot_in + S₂*ϵ

	state_symbol_array = [
		"A₁", "A₂", "B", "P", "C", "x","y"
	]

	# build state table -
	state_table₂ = Array{Any,2}(undef,ℳ₂,5)
	for state_index ∈ 1:ℳ₂

		state_table₂[state_index,1] = state_index
		state_table₂[state_index,2] = state_symbol_array[state_index]
		state_table₂[state_index,3] = n_dot_in[state_index]
		state_table₂[state_index,4] = n_dot_out[state_index]
		state_table₂[state_index,5] = n_dot_out[state_index] - n_dot_in[state_index]
	end

	# header row -
	state_table_header_row = (["index i","state","nᵢ_dot_in", "nᵢ_dot_out", "Δnᵢ_dot"],
		["","", "mol/time", "mol/time", "mol/time"]);

	# write the table -
	pretty_table(HTML, state_table₂; header=state_table_header_row, tf=tf_html_minimalist)
end

# ╔═╡ d201abce-202c-44f6-98a3-67e47c2a99f4
html"""
<style>
main {
    max-width: 900px;
    width: 75%;
    margin: auto;
    font-family: "Roboto, monospace";
}

a {
    color: blue;
    text-decoration: none;
}
</style>"""

# ╔═╡ 647315d2-e0d7-4a55-bab8-e8ff9d784b97
html"""
<script>
	// initialize -
	var section = 0;
	var subsection = 0;
	var subsubsection = 0;
	var headers = document.querySelectorAll('h3, h5, h6');

	// main loop -
	for (var i=0; i < headers.length; i++) {

		var header = headers[i];
	    var text = header.innerText;
	    var original = header.getAttribute("text-original");
	    if (original === null) {

			// Save original header text
	        header.setAttribute("text-original", text);
	    } else {

			// Replace with original text before adding section number
	        text = header.getAttribute("text-original");
	    }

	    var numbering = "";
	    switch (header.tagName) {
	        case 'H3':
	            section += 1;
	            numbering = section + ".";
	            subsection = 0;
	            break;
	        case 'H5':
	            subsection += 1;
	            numbering = section + "." + subsection;
	            break;
			case 'H6':
	            subsubsection += 1;
	            numbering = section + "." + subsection + "." + subsubsection;
	            break;
	    }
		// update the header text
		header.innerText = numbering + " " + text;
	};
</script>"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
RowEchelon = "af85af4c-bcd5-5d23-b03a-a909639aa875"

[compat]
Combinatorics = "~1.0.2"
GLPK = "~0.15.3"
IterativeSolvers = "~0.9.2"
Plots = "~1.25.7"
PlutoUI = "~0.7.32"
PrettyTables = "~1.3.1"
RowEchelon = "~0.2.1"
"""

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

julia_version = "1.7.1"
manifest_format = "2.0"

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

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

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

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

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

[[deps.BenchmarkTools]]
deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
git-tree-sha1 = "940001114a0147b6e4d10624276d56d531dd9b49"
uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
version = "1.2.2"

[[deps.BinaryProvider]]
deps = ["Libdl", "Logging", "SHA"]
git-tree-sha1 = "ecdec412a9abc8db54c0efc5548c64dfce072058"
uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
version = "0.5.10"

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

[[deps.CEnum]]
git-tree-sha1 = "215a9aa4a1f23fbd05b92769fdd62559488d70e9"
uuid = "fa961155-64e5-5f13-b03f-caf6b980ea82"
version = "0.4.1"

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

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

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

[[deps.CodecBzip2]]
deps = ["Bzip2_jll", "Libdl", "TranscodingStreams"]
git-tree-sha1 = "2e62a725210ce3c3c2e1a3080190e7ca491f18d7"
uuid = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd"
version = "0.7.2"

[[deps.CodecZlib]]
deps = ["TranscodingStreams", "Zlib_jll"]
git-tree-sha1 = "ded953804d019afa9a3f98981d99b33e3db7b6da"
uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
version = "0.7.0"

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

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

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

[[deps.Combinatorics]]
git-tree-sha1 = "08c8b6831dc00bfea825826be0bc8336fc369860"
uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
version = "1.0.2"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.GLPK]]
deps = ["BinaryProvider", "CEnum", "GLPK_jll", "Libdl", "MathOptInterface"]
git-tree-sha1 = "6f4e9754ee93e2b2ff40c0b0a6b4cdffd289190d"
uuid = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
version = "0.15.3"

[[deps.GLPK_jll]]
deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "fe68622f32828aa92275895fdb324a85894a5b1b"
uuid = "e8aa6df9-e6ca-548a-97ff-1f85fc5b8b98"
version = "5.0.1+0"

[[deps.GMP_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "781609d7-10c4-51f6-84f2-b8444358ff6d"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.IterativeSolvers]]
deps = ["LinearAlgebra", "Printf", "Random", "RecipesBase", "SparseArrays"]
git-tree-sha1 = "1169632f425f79429f245113b775a0e3d121457c"
uuid = "42fd0dbc-a981-5370-80f2-aaf504508153"
version = "0.9.2"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.MathOptInterface]]
deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "JSON", "LinearAlgebra", "MutableArithmetics", "OrderedCollections", "Printf", "SparseArrays", "Test", "Unicode"]
git-tree-sha1 = "38062215a56109442d4ec25a0a2970cbfef55bbc"
uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
version = "0.10.7"

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

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

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

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

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

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

[[deps.MutableArithmetics]]
deps = ["LinearAlgebra", "SparseArrays", "Test"]
git-tree-sha1 = "73deac2cbae0820f43971fad6c08f6c4f2784ff2"
uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
version = "0.3.2"

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

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

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

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

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

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

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

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

[[deps.Parsers]]
deps = ["Dates"]
git-tree-sha1 = "0b5cfbb704034b5b4c1869e36634438a047df065"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.2.1"

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

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

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

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

[[deps.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", "Unzip"]
git-tree-sha1 = "7e4920a7d4323b8ffc3db184580598450bde8a8e"
uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
version = "1.25.7"

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

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

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

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

[[deps.Profile]]
deps = ["Printf"]
uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"

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

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

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

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

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

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

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

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

[[deps.RowEchelon]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f479526c4f6efcbf01e7a8f4223d62cfe801c974"
uuid = "af85af4c-bcd5-5d23-b03a-a909639aa875"
version = "0.2.1"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.TranscodingStreams]]
deps = ["Random", "Test"]
git-tree-sha1 = "216b95ea110b5972db65aa90f88d8d89dcb8851c"
uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
version = "0.9.6"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

# ╔═╡ Cell order:
# ╟─7ed39fbc-8c1a-4d10-a6b7-7608e6e87395
# ╟─a55ad25a-e4ff-4af7-9684-c1b02fec900d
# ╟─dc712cba-1f18-46ec-a7f1-b25b0aca1e0a
# ╟─306e390f-acbe-4b2a-8e4f-3571003359ad
# ╠═be4b1854-1a14-4452-ad7e-22d614740a10
# ╠═8965f69d-3014-46b3-816d-7d6f7fd57adf
# ╠═d0e0e084-c706-45cf-b40f-47ca2e5ce088
# ╠═8be1e489-7382-4315-8c4c-111abdead290
# ╠═a655928e-a275-4bb1-bca3-39640804ed63
# ╠═9108a028-a5fb-428e-9fd8-2c7f0d80f4af
# ╠═2b9ca1ce-9080-4c42-91a1-e5005ddc4ee5
# ╟─4b32c7ae-397f-43c4-ac43-d3107b4de9c4
# ╟─0a325629-0a14-426f-a99e-9bccdc2e0dfa
# ╟─009df867-8c26-422a-8adf-a83da0b667bd
# ╟─3f98b575-eb2c-40b3-ac6f-43c55e4ec9d3
# ╟─72e5c813-2545-4d7c-9437-dc3f98f4e6e7
# ╠═45222d77-906b-4393-b2f9-63452300c28e
# ╠═87436aa6-289f-48f7-93ed-865360f27d92
# ╟─6dd988f9-28d3-47cb-8b32-a415132b388b
# ╟─39f3f612-114e-4f71-be75-ec427607e577
# ╟─a40ef7cf-3e46-4f4b-9af9-5a0c88d15ac9
# ╟─a7dd770b-22b3-4561-a85b-6071e7142c6d
# ╟─3a091b76-76cd-479c-b382-d5f2e3e75c59
# ╟─0fd5a620-682e-48d8-b1e6-ba110018f9e3
# ╠═1f31384d-13f8-4f31-9e23-9f6ad18eda23
# ╟─bc52eb4b-9761-459c-8a6e-c457f13d7335
# ╟─c1368f87-bd93-448f-ad10-b70040f5c854
# ╟─b620232d-6c8f-478c-93e5-9b4743518ff6
# ╟─ebb0d7cf-fd92-45be-ab98-5bb4de2c2ba2
# ╠═0e46d5ca-197f-47ca-a075-7281e3087460
# ╠═c51f8430-8147-11ec-3501-498c9bf67b57
# ╠═907af081-4a01-4a51-bda9-3a34295f9208
# ╟─d201abce-202c-44f6-98a3-67e47c2a99f4
# ╠═647315d2-e0d7-4a55-bab8-e8ff9d784b97
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
# This file is a part of GERDADeepLearning.jl, licensed under the MIT License (MIT).

using HDF5, Base.Threads, MultiThreadingTools

type EventFormatException <: Exception
  msg::String
end

function create_libroot(h5file, libname)
  libroot = g_create(h5file, libname)
  attrs(libroot)["type"] = "EventLibrary"
  attrs(libroot)["version"] = "1.0"
  return libroot
end

"""
The keylist must be iterable and implement listname.
The entries must support the string method.
"""
function create_extendible_hdf5_files(output_dir, keylists, detector_names, sample_size, chunk_size, label_keys)

  sample_size = Int64(sample_size) # Prevents HDF5 error
  chunk_size = Int64(chunk_size)

  h5files = [h5open(joinpath(output_dir, "$dname.h5"), "w") for dname in detector_names]

  label_arrays = Dict[]

  for (i,h5file) in enumerate(h5files)
    libroot = create_libroot(h5file, detector_names[i])

    detector_labels = Dict()

    waveforms = d_create(libroot, "waveforms", Float32, ((sample_size,chunk_size), (sample_size,-1)), "chunk", (sample_size,chunk_size)) # requires Int64 sizes
    set_dims!(waveforms, (sample_size, 0))
    detector_labels[:waveforms] = waveforms

    labelsg = g_create(libroot, "labels")
    push!(label_arrays, detector_labels)
    for (label_key, dtype) in label_keys
      labeld = d_create(labelsg, "$label_key", dtype, ((chunk_size,), (-1,)), "chunk", (chunk_size,))
      set_dims!(labeld, (0,))
      detector_labels[label_key] = labeld
    end

    properties = g_create(libroot, "properties")
    attrs(properties)["name"] = detector_names[i]
    attrs(properties)["sampling_rate"] = 100e6
    attrs(properties)["detector_id"] = i
    attrs(properties)["detector_name"] = detector_names[i]
    attrs(properties)["waveform_type"] = "raw"
    for (i,keylist) in enumerate(keylists)
      klg = g_create(properties, "Keylist$i")
      attrs(klg)["name"] = listname(keylist)
      attrs(klg)["entries"] = [string(key) for key in keylist]
    end
  end

  return h5files, label_arrays
end

function lazy_read_all(dir::AbstractString)
  files = readdir(dir)
  files = filter(f->endswith(f, ".h5"), files)
  names = [file[1:end-3] for file in files]
  libs = [lazy_read_library(joinpath(dir, files[i]), names[i]) for i in 1:length(files)]
  result = DLData(libs)
  result.dir = dir
  return result
end

function lazy_read_library(h5_filepath, libname)
  init = lib -> _initialize_from_file(lib, h5_filepath, libname)
  lib =  EventLibrary(init)

  ifile = h5open(h5_filepath, "r")
  libroot = ifile[libname]

  if read(attrs(libroot)["type"]) != "EventLibrary"
    throw(EventFormatException("Not a valid EventLibrary: $(read(attrs(libroot)["type"]))"))
  end
  if read(attrs(libroot)["version"]) != "1.0"
    throw(EventFormatException("Unknown EventLibrary version: $(read(attrs(libroot)["version"]))"))
  end

  # Read properties
  props = libroot["properties"]
  for key in names(attrs(props))
    value = read(attrs(props)[key])
    lib.prop[Symbol(key)] = value
  end
  lib.prop[:eventcount] = size(libroot["waveforms"], 2)
  # TODO keylists

  # Initialize label names
  labels = libroot["labels"]
  for key in names(labels)
    lib.labels[Symbol(key)] = zeros(Float32, 0)
  end

  close(ifile)

  return lib
end

function _initialize_from_file(lib::EventLibrary, h5_filepath, libname)
  h5open(h5_filepath, "r") do ifile
    libroot = ifile[libname]

    # Read waveforms
    try
      lib.waveforms = read(libroot["waveforms"])
    catch err
      info("Illegal waveform data for lib $(lib[:name]).")
      lib.waveforms = zeros(Float32, 1, 0)
    end

    # Read labels
    labels = libroot["labels"]
    for key in names(labels)
      data = read(labels[key])
      lib.labels[Symbol(key)] = data
    end
  end
end

export write_all
function write_all_sequentially(data::DLData, dir::AbstractString, uninitialize::Bool)
  isdir(dir) || mkdir(dir)
  for lib in data
    filepath = joinpath(dir, lib[:name]*".h5")
    write_lib(lib, filepath, uninitialize)
  end
end

export write_all
function write_all_multithreaded(data::DLData, dir::AbstractString, uninitialize::Bool)
  isdir(dir) || mkdir(dir)
  @everythread begin
    for i in threadpartition(1:length(data))
      lib = data.entries[i]
        filepath = joinpath(dir, lib[:name]*".h5")
        write_lib(lib, filepath, uninitialize)
    end
  end
end

export write_lib
function write_lib(lib::EventLibrary, filepath::AbstractString, uninitialize::Bool)
  h5open(filepath, "w") do h5file
    libroot = create_libroot(h5file, lib[:name])

    write(libroot, "waveforms", waveforms(lib))

    labelsg = g_create(libroot, "labels")
    for (key,value) in lib.labels
      write(labelsg, string(key), value)
    end

    properties = g_create(libroot, "properties")
    for(key,value) in lib.prop
      try
        attrs(properties)[string(key)] = value
      catch err
        threadsafe_info("Cannot write property $key: $err")
      end
    end
  end

  if uninitialize
    dispose(lib)
    lib.initialization_function = lib2 -> _initialize_from_file(lib2, filepath, lib[:name])
  end
end

function _str_to_sym_dict(dict)
  result = Dict{Symbol, Any}()
  for (key,value) in dict
    result[Symbol(key)] = value
  end
  return result
end

function segh5_to_hdf5(formattype::AbstractString, keylists::Vector{Vector{AbstractString}}, output_dir::AbstractString, verbosity::Integer)
  all_files = vcat(keylists...)

  sample_size = _segh5_get_waveform_length(all_files[1])
  verbosity > 2 && info("Determined number of samples per waveform = $sample_size.")

  label_keys = [:keylist=>Int32, :filenum=>Int32, :E=>Float32, :E1=>Float32, :E2=>Float32, :E3=>Float32, :E4=>Float32]

  # Create HDF5 extendible storage
  h5files, h5_arrays = create_extendible_hdf5_files(output_dir, [], [formattype], sample_size, 256, label_keys)
  out_keylist = h5_arrays[1][:keylist]
  out_filenum = h5_arrays[1][:filenum]
  out_E = h5_arrays[1][:E]
  out_E1 = h5_arrays[1][:E1]
  out_E2 = h5_arrays[1][:E2]
  out_E3 = h5_arrays[1][:E3]
  out_E4 = h5_arrays[1][:E4]
  out_waveforms = h5_arrays[1][:waveforms]

  try
      file_count = length(all_files)
      file_i = 0
      for(kl_i, filelist) in enumerate(keylists)
          for file in filelist
            if !isfile(file)
              verbosity >= 1 && info("Skipping entry because file does not exist: $file")
            else
              verbosity >= 2 && info("Reading file $file_i / $file_count...")
              verbosity >= 3 && info("File: $file")

              h5in = h5open(file)

              # println(size(out_waveforms))
              # println(size(h5in["DAQ_Data"]["daq_pulses"]))
              # flush(STDOUT)

              f_energies = read(h5in["Processed_data"]["energies"])
              f_pulses = h5in["DAQ_Data"]["daq_pulses"][1,:,:][1,:,:]

              extend_h5_array(out_keylist, fill(kl_i, size(f_energies,2)))
              extend_h5_array(out_filenum, fill(file_i, size(f_energies,2)))

              extend_h5_array(out_E, f_energies[1,:])
              extend_h5_array(out_E1, f_energies[2,:])
              extend_h5_array(out_E2, f_energies[3,:])
              extend_h5_array(out_E3, f_energies[4,:])
              extend_h5_array(out_E4, f_energies[5,:])

              extend_h5_matrix(out_waveforms, f_pulses)

              out_E
            end
            file_i += 1
          end
      end
  catch exc
      #verbosity >= 1 && info("Fatal error during conversion to HDF5. Closing HDF5 files and rethrowing error.")
      for h5file in h5files close(h5file) end
      rethrow(exc)
  end

verbosity >= 3 && info("Finished conversion without fatal errors. Closing HDF5 files...")

for h5file in h5files close(h5file) end

verbosity >= 3 && info("All HDF5 files closed.")
end

function _segh5_get_waveform_length(filepath)
  file = h5open(filepath)
  return size(file["DAQ_Data"]["daq_pulses"])[2] # #channels, #samples, #events
end

function extend_h5_array(array::HDF5Dataset, newdata::Vector)
  set_dims!(array, (length(array)+length(newdata),) )
  array[(end-length(newdata)+1) : end] = newdata
end

function extend_h5_matrix(matrix::HDF5Dataset, newdata::Matrix)
  set_dims!(matrix, (size(matrix,1), size(matrix,2)+size(newdata,2)))
  matrix[:, (end-size(newdata, 2)+1) : end] = newdata
end

ext: jl


text:
module FileIO
# This module handles all the io of files, and parsing.
using DelimitedFiles, Dates, LinearAlgebra, UnitfulAtomic
using UnitfulAtomic.Unitful: ustrip, uconvert
using UnitfulAtomic: bohr
using ..DFControl
using ..Utils
using ..Calculations
using ..Structures
using ..Jobs

include("qe.jl")
include("wannier.jl")
include("job.jl")

function parse_file(filename::AbstractString, parse_funcs::Vector{<:Pair{String}};out = Dict{Symbol,Any}(),
                    extra_parse_funcs::Vector{<:Pair} = Pair{String,Function}[])

    open(filename, "r") do f
        lc = 0
        while !eof(f)
            line = strip(readline(f))
            lc += 1
            if isempty(line)
                continue
            end
            for pf in (parse_funcs, extra_parse_funcs)
                func = getfirst(x -> occursin(x[1], line), pf)
                func === nothing && continue
                try
                    func[2](out, line, f)
                catch
                    @warn "File corruption or parsing error detected executing parse function \n$(func[2]) in file $filename at line $lc: \"$line\".\nTrying to continue smoothly."
                end
            end
        end
    end
    return out
end

function cif2structure(cif_file::String; structure_name = "NoName")
    tmpdir = dirname(cif_file)
    tmpfile = joinpath(tmpdir, "tmp.in")
    @assert splitext(cif_file)[2] == ".cif" error("Please specify a valid cif calculation file")
    run(`$pythonpath $cif2cellpath $cif_file --no-reduce -p quantum-espresso -o $tmpfile`)

    bla, structure = qe_read_calculation(tmpfile; structure_name = structure_name)
    rm(tmpfile)
    return structure
end

Base.length(::Type{<:Real}) = 1

function Base.parse(::Type{NamedTuple{names,types}},
                    spl::Vector{<:AbstractString}) where {names,types}
    @assert sum(length.(types.parameters)) == length(spl)
    tbuf = Vector{Union{types.parameters...}}(undef, length(types.parameters))
    counter = 1
    for (i, typ) in enumerate(types.parameters)
        if typ <: AbstractVector
            l = length(typ)
            tbuf[i] = parse(typ, spl[counter:counter+l-1])
            counter += l
        else
            tbuf[i] = parse(typ, spl[counter])
            counter += 1
        end
    end
    pstring = (tbuf...,)
    return NamedTuple{names}(pstring)
end
Base.parse(::Type{T}, s::AbstractString) where {T<:NamedTuple} = parse(T, split(s))

function Base.parse(::Type{Point{N,T}}, spl::Vector{<:AbstractString}) where {N,T}
    @assert N == length(spl)
    return Point{N,T}(parse.(T, spl))
end
Base.parse(::Type{T}, s::AbstractString) where {T<:Point} = parse(T, split(s))

function write_xsf(filename::AbstractString, structure::Structure)
    open(filename, "w") do f
        write(f, "CRYSTAL\n")
        c = ustrip.(structure.cell')
        write(f, "PRIMVEC\n")
        write(f, "$(c[1,1]) $(c[1,2]) $(c[1,3])\n")
        write(f, "$(c[2,1]) $(c[2,2]) $(c[2,3])\n")
        write(f, "$(c[3,1]) $(c[3,2]) $(c[3,3])\n")
        write(f, "CONVVEC\n")
        write(f, "$(c[1,1]) $(c[1,2]) $(c[1,3])\n")
        write(f, "$(c[2,1]) $(c[2,2]) $(c[2,3])\n")
        write(f, "$(c[3,1]) $(c[3,2]) $(c[3,3])\n")
        write(f, "PRIMCOORD\n")
        write(f, "$(length(structure.atoms)) 1\n")
        for at in structure.atoms
            n = at.element.symbol
            p = ustrip.(at.position_cart)
            write(f, "$n $(p[1]) $(p[2]) $(p[3])\n")
        end
    end
end

end


text:
#precompute a few numbers
const 𝐬𝐪𝛑 = √π
const 𝐨𝐬𝐪𝛑𝐥𝐧2 = 1/sqrt(π/log(2.0))
const 𝐬𝐪𝐥𝐧2 = sqrt(log(2.0))
const 𝐜₂ = 100.0*𝐡*𝐜/𝐤

#-------------------------------------------------------------------------------
# wavenumber truncation of line shapes

cutline(ν, νl, Δνcut)::Bool = abs(ν - νl) > Δνcut ? true : false

function includedlines(ν::Real,
                       νl::AbstractVector{<:Real},
                       Δνcut::Real)::Vector{Int64}
    findall(x -> !cutline(ν, x, Δνcut), νl)
end

function includedlines(ν::AbstractVector{<:Real},
                       νl::AbstractVector{<:Real},
                       Δνcut::Real)::Vector{Int64}
    findall((νl .> minimum(ν) - Δνcut) .& (νl .< maximum(ν) + Δνcut))
end

#-------------------------------------------------------------------------------
# Chebyshev polynomial fit for Qref/Q

function chebyQrefQ(T::Real, n::Int64, a::Vector{Float64})
    #check the temperature range
    @assert TMIN <= T <= TMAX "temperature outside of Qref/Q interpolation range [$TMIN, $TMAX]"
    #map T to [-1,1]
    τ::Float64 = 2*(T - TMIN)/(TMAX - TMIN) - 1
    #values of first two chebys at τ
    c₁::Float64 = 1.0
    c₂::Float64 = τ
    #value of expansion after first two terms
    @inbounds y = a[1] + a[2]*c₂
    for k = 3:n
        #next cheby value
        c₃::Float64 = 2*τ*c₂ - c₁
        #contribute to expansion
        @inbounds y += a[k]*c₃
        #swap values
        c₁ = c₂
        c₂ = c₃
    end
    #return the inverse, Qref/Q
    return 1.0/y
end

#-------------------------------------------------------------------------------
# special strategy for line profiles from sorted vectors of wavenumbers

function surf!(σ::AbstractVector,
               f::F,
               ν::AbstractVector,
               νl::AbstractVector,
               Δνcut::Real,
               A::Vararg{Vector{Float64},N}) where {F<:Function, N}
    @assert all(diff(ν) .> 0) "wavenumber vectors must be sorted in ascending order"
    L = length(νl)
    j₁ = 1 #tracking index to avoid searching from beginning every time
    for i ∈ eachindex(ν)
        #temporary stack value
        σᵢ::Float64 = 0.0
        #find the first line that isn't cut off
        j = j₁
        @inbounds while (j <= L) && cutline(ν[i], νl[j], Δνcut)
            j += 1
        end
        #only proceed if there is a line to include
        if j <= L
            #update the starting index for the search
            j₁ = j
            #evaluate line profiles until one gets cut off, then move on
            @inbounds while (j <= L) && !cutline(ν[i], νl[j], Δνcut)
                #let block is required for good performance
                args = let k = j
                    @inbounds ntuple(n->A[n][k], N)
                end
                @inbounds σᵢ += f(ν[i], νl[j], args...)
                j += 1
            end
        end
        #set the array value
        σ[i] = σᵢ
    end
end

#-------------------------------------------------------------------------------
# temperature scaling of line intensity

export scaleintensity

"""
    scaleintensity(S, νl, Epp, M, I, T)

Compute the [temperature scaling for line intensity](https://hitran.org/docs/definitions-and-units/#mjx-eqn-eqn-intensity-temperature-dependence).

# Arguments
* `S`: spectal line intensity at 296 K [cm``^{-1}``/(molecule``\\cdot``cm``^{-2}``)]
* `νl`: wavenumber of line [cm``^{-1}``]
* `Epp`: lower-state energy of transition [cm``^{-1}``]
* `M`: [HITRAN molecular identification number](https://hitran.org/docs/molec-meta)
* `I`: [HITRAN local isotopologue number](https://hitran.org/docs/iso-meta/)
* `T`: temperature [K]
"""
function scaleintensity(S, νl, Epp, M::Int16, I::Int16, T)
    #arguments to exp
    a = -𝐜₂*Epp
    b = -𝐜₂*νl
    #numerator and denominator
    n = exp(a/T)*(1 - exp(b/T))
    d = exp(a/𝐓ᵣ)*(1 - exp(b/𝐓ᵣ))
    #check if there is an approximating function
    if MOLPARAM[M].hascheb[I]
        QrefQ = chebyQrefQ(T, MOLPARAM[M].ncheb[I], MOLPARAM[M].cheb[I])
    else
        throw("no interpolating polynomial available to compute Qref/Q for isotopologue $I of $(MOLPARAM[M].name) ($(MOLPARAM[M].formula))")
        #QrefQ = (𝐓ᵣ/T)^1.5
    end
    #shifted line intensity
    S*QrefQ*(n/d)
end

function scaleintensity(sl::SpectralLines, idx::Vector{Int64}, T)::Vector{Float64}
    Sₛ = zeros(Float64, length(idx))
    for i ∈ eachindex(idx)
        j = idx[i]
        @inbounds Sₛ[i] = scaleintensity(sl.S[j], sl.ν[j], sl.Epp[j], sl.M[j], sl.I[j], T)
    end
    return Sₛ
end

#-------------------------------------------------------------------------------
# doppler broadening

export αdoppler, fdoppler, doppler, doppler!

"""
    αdoppler(νl, μ, T)

Compute doppler (gaussian) broadening coefficient from line wavenumber `νl` [cm``^{-1}``], gas molar mass `μ` [kg/mole], and temperature `T` [K].
"""
αdoppler(νl, μ, T) = (νl/𝐜)*sqrt(2.0*𝐑*T/μ)

function αdoppler(sl::SpectralLines, i::Vector{Int64}, T)::Vector{Float64}
    αdoppler.(view(sl.ν,i), view(sl.μ,i), T)
end

"""
    fdoppler(ν, νl, α)

Evaluate doppler (gaussian) profile

# Arguments
* `ν`: profile evaluation wavenumber [cm``^{-1}``]
* `νl`: wavenumber of absorption line [cm``^{-1}``]
* `α`: doppler (gaussian) broadening coefficient
"""
fdoppler(ν, νl, α) = exp(-(ν - νl)^2/α^2)/(α*𝐬𝐪𝛑)

"""
    doppler(ν, νl, S, α)

Evaluate doppler (gaussian) absoption cross-section [cm``^2``/molecule]

# Arguments
* `ν`: profile evaluation wavenumber [cm``^{-1}``]
* `νl`: wavenumber of absorption line [cm``^{-1}``]
* `S`: line absoption intensity [cm``^{-1}``/(molecule``\\cdot``cm``^{-2}``)]
* `α`: doppler (gaussian) broadening coefficient
"""
doppler(ν, νl, S, α) = S*fdoppler(ν, νl, α)

"""
    doppler(ν, sl, T, P, Pₚ, Δνcut=25)

Evaluate a single doppler (gaussian) absoption cross-section [cm``^2``/molecule]. Temperature scaling and doppler profiles are evaluated along the way.

# Arguments
* `ν`: wavenumber indicating where to evaluate [cm``^{-1}``]
* `sl`: [`SpectralLines`](@ref)
* `T`: temperature [K]
* `P`: air pressure [Pa]
* `Pₚ`: partial pressure [Pa]
* `Δνcut`: profile truncation distance [cm``^{-1}``]
"""
function doppler(ν, sl::SpectralLines, T, P, Pₚ, Δνcut=25.0)
    i = includedlines(ν, sl.ν, Δνcut)
    S = scaleintensity(sl, i, T)
    α = αdoppler(sl, i, T)
    sum(doppler.(ν, view(sl.ν,i), S, α))
end

"""
    doppler!(σ, ν, sl, T, P, Pₚ, Δνcut=25)

Identical to [`doppler`](@ref), but fills the vector of cross-sections (`σ`) in-place.
"""
function doppler!(σ::AbstractVector,
                  ν::AbstractVector,
                  sl::SpectralLines,
                  T,
                  P,
                  Pₚ,
                  Δνcut=25.0)
    i = includedlines(ν, sl.ν, Δνcut)
    S = scaleintensity(sl, i, T)
    α = αdoppler(sl, i, T)
    surf!(σ, doppler, ν, view(sl.ν,i), Δνcut, S, α)
end

"""
    doppler(ν, sl, T, P, Pₚ, Δνcut=25)

Compute a vector of doppler (gaussian) absorption cross-sections [cm``^2``/molecule] from a [`SpectralLines`](@ref) object. Temperature scaling and doppler profiles are evaluated along the way.

# Arguments
* `ν`: vector of wavenumbers indicating where to evaluate [cm``^{-1}``]
* `sl`: [`SpectralLines`](@ref)
* `T`: temperature [K]
* `P`: air pressure [Pa]
* `Pₚ`: partial pressure [Pa]
* `Δνcut`: profile truncation distance [cm``^{-1}``]
"""
function doppler(ν::AbstractVector,
                 sl::SpectralLines,
                 T,
                 P,
                 Pₚ,
                 Δνcut=25.0)::Vector{Float64}
    σ = zeros(Float64, length(ν))
    doppler!(σ, ν, sl, T, P, Pₚ, Δνcut)
    return σ
end

#-------------------------------------------------------------------------------
# pressure broadening

export γlorentz, florentz, lorentz, lorentz!

"""
    γlorentz(γa, γs, na, T, P, Pₚ)

Compute lorentzian broadening coefficient

# Arguments
* `γa`: air-broadened half width at half maximum (HWHM) [cm``^{-1}``/atm] at 296 K and 1 atm
* `γs`: self-broadened half width at half maximum (HWHM) [cm``^{-1}``/atm] at 296 K and 1 atm
* `na`: coefficient of temperature dependence of air-broadened half width
* `T`: temperature [K]
* `P`: air pressure [Pa]
* `Pₚ`: partial pressure [Pa]
"""
function γlorentz(γa, γs, na, T, P, Pₚ)
    ((𝐓ᵣ/T)^na)*(γa*(P - Pₚ) + γs*Pₚ)/𝐀
end

function γlorentz(sl::SpectralLines, i::Vector{Int64}, T, P, Pₚ)::Vector{Float64}
    γlorentz.(view(sl.γa,i), view(sl.γs,i), view(sl.na,i), T, P, Pₚ)
end

"""
    florentz(ν, νl, γ)

Evaluate lorentz profile

# Arguments
* `ν`: profile evaluation wavenumber [cm``^{-1}``]
* `νl`: wavenumber of absorption line [cm``^{-1}``]
* `γ`: lorentzian broadening coefficient
"""
florentz(ν, νl, γ) = γ/(π*((ν - νl)*(ν - νl) + γ*γ))

"""
    lorentz(ν, νl, S, γ)

Evaluate lorentzian absoption cross-section [cm``^2``/molecule]

# Arguments
* `ν`: profile evaluation wavenumber [cm``^{-1}``]
* `νl`: wavenumber of absorption line [cm``^{-1}``]
* `S`: line absoption intensity [cm``^{-1}``/(molecule``\\cdot``cm``^{-2}``)]
* `γ`: lorentzian broadening coefficient
"""
lorentz(ν, νl, S, γ) = S*florentz(ν, νl, γ)

"""
    lorentz(ν, sl, T, P, Pₚ, Δνcut=25)

Compute a single lorentzian absorption cross-sections [cm``^2``/molecule] from a [`SpectralLines`](@ref) object. Temperature scaling and lorentzian profiles are evaluated along the way.

# Arguments
* `ν`: single wavenumber indicating where to evaluate [cm``^{-1}``]
* `sl`: [`SpectralLines`](@ref)
* `T`: temperature [K]
* `P`: air pressure [Pa]
* `Pₚ`: partial pressure [Pa]
* `Δνcut`: profile truncation distance [cm``^{-1}``]
"""
function lorentz(ν, sl::SpectralLines, T, P, Pₚ, Δνcut=25.0)
    i = includedlines(ν, sl.ν, Δνcut)
    S = scaleintensity(sl, i, T)
    γ = γlorentz(sl, i, T, P, Pₚ)
    sum(lorentz.(ν, view(sl.ν,i), S, γ))
end

"""
    lorentz!(σ, ν, sl, T, P, Pₚ, Δνcut=25)

Identical to [`lorentz`](@ref), fills the vector of cross-sections (`σ`) in-place.
"""
function lorentz!(σ::AbstractVector,
                  ν::AbstractVector,
                  sl::SpectralLines,
                  T,
                  P,
                  Pₚ,
                  Δνcut=25.0)
    i = includedlines(ν, sl.ν, Δνcut)
    S = scaleintensity(sl, i, T)
    γ = γlorentz(sl, i, T, P, Pₚ)
    surf!(σ, lorentz, ν, view(sl.ν,i), Δνcut, S, γ)
end

"""
    lorentz(ν, sl, T, P, Pₚ, Δνcut=25)

Compute a vector of lorentzian absorption cross-sections [cm``^2``/molecule] from a [`SpectralLines`](@ref) object. Temperature scaling and lorentzian profiles are evaluated along the way.

# Arguments
* `ν`: vector of wavenumbers indicating where to evaluate [cm``^{-1}``]
* `sl`: [`SpectralLines`](@ref)
* `T`: temperature [K]
* `P`: air pressure [Pa]
* `Pₚ`: partial pressure [Pa]
* `Δνcut`: profile truncation distance [cm``^{-1}``]
"""
function lorentz(ν::AbstractVector,
                 sl::SpectralLines,
                 T,
                 P,
                 Pₚ,
                 Δνcut=25.0)::Vector{Float64}
    σ = zeros(Float64, length(ν))
    lorentz!(σ, ν, sl, T, P, Pₚ, Δνcut)
    return σ
end

#-------------------------------------------------------------------------------
# voigt profile

export fvoigt, voigt, voigt!

"""
    fvoigt(ν, νl, α, γ)

Evaluate Voigt profile

# Arguments
* `ν`: profile evaluation wavenumber [cm``^{-1}``]
* `νl`: wavenumber of absorption line [cm``^{-1}``]
* `α`: doppler (gaussian) broadening coefficient
* `γ`: lorentzian broadening coefficient
"""
function fvoigt(ν, νl, α, γ)
    #inverse of the doppler parameter
    β = 1/α
    #factor for real and complex parts of Faddeeva args, avoiding β division
    d = 𝐬𝐪𝐥𝐧2*β
    #arguments to Faddeeva function
    x = (ν - νl)*d
    y = γ*d
    #evaluate real part of Faddeeva function
    f = faddeyeva(x,y)
    #final calculation, avoiding α division by using β again
    𝐨𝐬𝐪𝛑𝐥𝐧2*β*f
end

"""
    voigt(ν, νl, S, α, γ)

Evaluate Voigt absoption cross-section [cm``^2``/molecule]

# Arguments
* `ν`: profile evaluation wavenumber [cm``^{-1}``]
* `νl`: wavenumber of absorption line [cm``^{-1}``]
* `S`: line absoption intensity [cm``^{-1}``/(molecule``\\cdot``cm``^{-2}``)]
* `α`: doppler (gaussian) broadening coefficient
* `γ`: lorentzian broadening coefficient
"""
voigt(ν, νl, S, α, γ) = S*fvoigt(ν, νl, α, γ)

"""
    voigt(ν, sl::SpectralLines, T, P, Pₚ, Δνcut=25)

Evaluate Voigt absorption cross-section at a single wavenumber.
"""
function voigt(ν, sl::SpectralLines, T, P, Pₚ, Δνcut=25.0)
    i = includedlines(ν, sl.ν, Δνcut)
    S = scaleintensity(sl, i, T)
    α = αdoppler(sl, i, T)
    γ = γlorentz(sl, i, T, P, Pₚ)
    sum(voigt.(ν, view(sl.ν,i), S, α, γ))
end

"""
    voigt!(σ, ν, sl, T, P, Pₚ, Δνcut=25)

Identical to [`voigt`](@ref), but fills the vector of cross-sections (`σ`) in-place.
"""
function voigt!(σ::AbstractVector,
                ν::AbstractVector,
                sl::SpectralLines,
                T,
                P,
                Pₚ,
                Δνcut=25.0)
    i = includedlines(ν, sl.ν, Δνcut)
    S = scaleintensity(sl, i, T)
    α = αdoppler(sl, i, T)
    γ = γlorentz(sl, i, T, P, Pₚ)
    surf!(σ, voigt, ν, view(sl.ν,i), Δνcut, S, α, γ)
end

"""
    voigt(ν, sl, T, P, Pₚ, Δνcut=25)

Compute a vector of Voigt absorption cross-sections [cm``^2``/molecule] from a [`SpectralLines`](@ref) object. Temperature scaling and Voigt profiles are evaluated along the way.

# Arguments
* `ν`: vector of wavenumbers indicating where to evaluate [cm``^{-1}``]
* `sl`: [`SpectralLines`](@ref)
* `T`: temperature [K]
* `P`: air pressure [Pa]
* `Pₚ`: partial pressure [Pa]
* `Δνcut`: profile truncation distance [cm``^{-1}``]
"""
function voigt(ν::AbstractVector,
               sl::SpectralLines,
               T,
               P,
               Pₚ,
               Δνcut=25.0)::Vector{Float64}
    σ = zeros(Float64, length(ν))
    voigt!(σ, ν, sl, T, P, Pₚ, Δνcut)
    return σ
end

#-------------------------------------------------------------------------------
# sublorentzian profile for CO2
# Perrin, M. Y., and J. M. Hartman. “Temperature-Dependent Measurements and Modeling of Absorption by CO2-N2 Mixtures in the Far Line-Wings of the 4.3 Μm CO2 Band.” Journal of Quantitative Spectroscopy and Radiative Transfer, vol. 42, no. 4, 1989, pp. 311–17.

export ΧPHCO2, PHCO2, PHCO2!

"""
    ΧPHCO2(ν, νl, T)

Compute the `Χ` (Chi) factor for sub-lorentzian CO2 line profiles, as in
* [Perrin, M. Y., and J. M. Hartman. “Temperature-Dependent Measurements and Modeling of Absorption by CO2-N2 Mixtures in the Far Line-Wings of the 4.3 Μm CO2 Band.” Journal of Quantitative Spectroscopy and Radiative Transfer, vol. 42, no. 4, 1989, pp. 311–17.](https://www.sciencedirect.com/science/article/abs/pii/0022407389900770)

# Arguments
* `ν`: profile evaluation wavenumber [cm``^{-1}``]
* `νl`: wavenumber of absorption line [cm``^{-1}``]
* `T`: temperature [K]
"""
function ΧPHCO2(ν, νl, T)
    Δν = abs(ν - νl)
    if Δν < 3.0
        return 1.0
    end
    B1 = 0.0888 - 0.16*exp(-0.0041*T)
    if Δν < 30.0
        return exp(-B1*(Δν - 3.0))
    end
    B2 = 0.0526*exp(-0.00152*T)
    if Δν < 120.0
        return exp(-B1*27.0 - B2*(Δν - 30.0))
    end
    return exp(-B1*27.0 - B2*90.0 - 0.0232*(Δν - 120.0))
end

"""
    PHCO2(ν, νl, S, α)

Evaluate Perrin & Hartman sub-lorentzian absoption cross-section [cm``^2``/molecule] for CO2

# Arguments
* `ν`: profile evaluation wavenumber [cm``^{-1}``]
* `νl`: wavenumber of absorption line [cm``^{-1}``]
* `T`: temperature [K]
* `S`: line absoption intensity [cm``^{-1}``/(molecule``\\cdot``cm``^{-2}``)]
* `α`: doppler (gaussian) broadening coefficient
* `γ`: lorentzian broadening coefficient
"""
function PHCO2(ν, νl, T, S, α, γ)
    Χ = ΧPHCO2(ν, νl, T)
    voigt(ν, νl, S, α, Χ*γ)
end

"""
    PHCO2(ν, sl, T, P, Pₚ, Δνcut=500)

Compute a single Perrin & Hartman sub-lorentzian CO2 absorption cross-sections [cm``^2``/molecule] from a [`SpectralLines`](@ref) object. Temperature scaling and profiles are evaluated along the way.

# Arguments
* `ν`: single wavenumber indicating where to evaluate [cm``^{-1}``]
* `sl`: [`SpectralLines`](@ref)
* `T`: temperature [K]
* `P`: air pressure [Pa]
* `Pₚ`: partial pressure [Pa]
* `Δνcut`: profile truncation distance [cm``^{-1}``]
"""
function PHCO2(ν, sl::SpectralLines, T, P, Pₚ, Δνcut=500.0)
    i = includedlines(ν, sl.ν, Δνcut)
    S = scaleintensity(sl, i, T)
    α = αdoppler(sl, i, T)
    γ = γlorentz(sl, i, T, P, Pₚ)
    sum(PHCO2.(ν, view(sl.ν,i), T, S, α, γ))
end

"""
    PHCO2!(σ, ν, sl, T, P, Pₚ, Δνcut=500)

Identical to [`PHCO2`](@ref), but fills the vector of cross-sections (`σ`) in-place.
"""
function PHCO2!(σ::AbstractVector,
                ν::AbstractVector,
                sl::SpectralLines,
                T,
                P,
                Pₚ,
                Δνcut=500.0)
    i = includedlines(ν, sl.ν, Δνcut)
    S = scaleintensity(sl, i, T)
    α = αdoppler(sl, i, T)
    γ = γlorentz(sl, i, T, P, Pₚ)
    f(ν, νl, S, α, γ) = PHCO2(ν, νl, T, S, α, γ) #shove T into the function
    surf!(σ, f, ν, view(sl.ν,i), Δνcut, S, α, γ)
end

"""
    PHCO2(ν, sl, T, P, Pₚ, Δνcut=500)

Compute a vector of Perrin & Hartman sub-lorentzian CO2 absorption cross-sections [cm``^2``/molecule] from a [`SpectralLines`](@ref) object. Temperature scaling and profiles are evaluated along the way.

# Arguments
* `ν`: vector of wavenumbers indicating where to evaluate [cm``^{-1}``]
* `sl`: [`SpectralLines`](@ref)
* `T`: temperature [K]
* `P`: air pressure [Pa]
* `Pₚ`: partial pressure [Pa]
* `Δνcut`: profile truncation distance [cm``^{-1}``]
"""
function PHCO2(ν::AbstractVector,
               sl::SpectralLines,
               T,
               P,
               Pₚ,
               Δνcut=500.0)::Vector{Float64}
    σ = zeros(Float64, length(ν))
    PHCO2!(σ, ν, sl, T, P, Pₚ, Δνcut)
    return σ
end