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

using Markdown
using InteractiveUtils

# ╔═╡ f701a8a6-3a30-4d86-9018-f201276f9369
begin

	# load some external packages
	using PlutoUI
	using DataFrames
	using BSON
	using GLPK
	using PrettyTables

	# setup my paths (where are my files?)
	_PATH_TO_ROOT = pwd()
	_PATH_TO_SRC = joinpath(_PATH_TO_ROOT,"src")
	_PATH_TO_MODEL = joinpath(_PATH_TO_ROOT,"model")

	# load the ENGRI 1120 project code library -
	include(joinpath(_PATH_TO_SRC,"Include.jl"))

	# load the model -
	MODEL = BSON.load(joinpath(_PATH_TO_MODEL,"model_v2.bson"), @__MODULE__)

	# show -
	nothing
end

# ╔═╡ 7d3b5c92-152f-416a-b82c-82682e61ab9d
md"""
### Example: Flux Balance Analysis using the ENGRI-1120 Project Networks

The FBA problem is (traditinally) a Linear Programming (LP) problem in which a linear objective function $\mathcal{O}$:

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

is minimized (or maximized) subject to constraints on the reaction bounds (written in terms of the open extent) and species bounds (species mole balances). In our design problem, these constraints take the form (shown for a single in, single out):

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

Let's set up these bounds for the design project, and run some sample calculations to get a feel for the performance of the metabolic reaction network.

"""

# ╔═╡ 762e417b-b3f3-4bc0-9c9c-f777ae273303
md"""
##### What's in the model_v2.bson file?
"""

# ╔═╡ bbfe850c-1f2c-4efb-9a9c-61bb2261370c
MODEL

# ╔═╡ a1507115-b73d-482a-b66c-734f58463e49
md"""
##### Setup the flux balance analysis calculation

Let's build the $\mathcal{M}\times\mathcal{R}$ stoichiometric array from data in the `model_v2.bson` file. Then, we can set up the reaction bounds, the species bounds, and the objective coefficient arrays (all of which are needed for the calculation).

"""

# ╔═╡ 39b2b92f-4d6b-4f8c-b840-da2fb0dc79e2
begin

	# setup the FBA calculation for the project -

	# === SELECT YOUR PRODUCT HERE ==================================================== #
	# What rate are trying to maximize? (select your product)
	# rn:R08199 = isoprene
	# rn:28235c0c-ec00-4a11-8acb-510b0f2e2687 = PGDN
	# rn:rn:R09799 = Hydrazine
	# rn:R03119 = 3G
	idx_target_rate = find_reaction_index(MODEL,:reaction_number=>"rn:28235c0c-ec00-4a11-8acb-510b0f2e2687")
	# ================================================================================= #

	# First, let's build the stoichiometric matrix from the model object -
	(cia,ria,S) = build_stoichiometric_matrix(MODEL);

	# Next, what is the size of the system? (ℳ = number of metabolites, ℛ = number of reactions)
	(ℳ,ℛ) = size(S)

	# Next, setup a default bounds array => update specific elements
	# We'll correct the directionality below -
	Vₘ = (13.7)*(3600)*(50e-9)*(1000) # units: mmol/hr
	flux_bounds = [-Vₘ*ones(ℛ,1) Vₘ*ones(ℛ,1)]

	# What is the default mol flow input array => update specific elements
	# strategy: start with nothing in both streams, add material(s) back
	n_dot_input_stream_1 = zeros(ℳ,1)	# stream 1
	n_dot_input_stream_2 = zeros(ℳ,1)	# stream 2

	# === YOU NEED TO CHANGE BELOW HERE ====================================================== #
	# Let's lookup stuff that we want/need to supply to the chip to get the reactiont to go -
	# what you feed *depends upon your product*
	compounds_that_we_need_to_supply = [
		"oxygen", "sucrose", "potassium nitrate"
	]

	# what are the amounts that we need to supply to chip (units: mmol/hr)?
	mol_flow_values = [
		10.0 ; # oxygen mmol/hr
		20.0 ; # sucrose mmol/hr
		10.0 ; # KNO3 mmol/hr
	]
	# === YOU NEED TO CHANGE ABOVE HERE ====================================================== #

	idx_supply = Array{Int64,1}()
	for compound in compounds_that_we_need_to_supply
		idx = find_compound_index(MODEL,:compound_name=>compound)
		push!(idx_supply,idx)
	end

	# supply -
	n_dot_input_stream_1[idx_supply] .= mol_flow_values

	# setup the species bounds array -
	species_bounds = [-1.0*(n_dot_input_stream_1.+n_dot_input_stream_2) 1000.0*ones(ℳ,1)]

	# Lastly, let's setup the objective function -
	c = zeros(ℛ)
	c[idx_target_rate] = -1.0

	# show -
	nothing
end

# ╔═╡ fe7efbdb-237c-4ac5-9e40-06df4264266c
md"""
#### Estimate the directionality of reaction flux

By default, we assume that all reactions are reversible. However, while this may be theoretically true we can get a better estimate of whether a reaction is likely to run backward (and hence a have lower bound $\neq$ 0) by
looking at the Gibbs energy of reaction $\Delta{G}_{r}$. In particular, the heuristic:

$$\mathrm{sgn}(\dot{\epsilon}_{j}) = -\mathrm{sgn}\left(\Delta{G}_{r,j}\right)$$

gives us a tool to assign the lower bound of $\dot{\epsilon}_{j}$ to zero, or to a non-zero value, where
$\mathrm{sgn}\left(\star\right)$ denotes the [sign function](https://en.wikipedia.org/wiki/Sign_function). If $\mathrm{sgn}(\dot{\epsilon}_{j})<0$, we have a non-zero lower bound (reaction is reversible).

"""

# ╔═╡ 537e0c28-b2d7-4d29-9592-0bb35747cd81
begin

	# update the flux bounds -> which fluxes can can backwards?
	# do determine this: sgn(v) = -1*sgn(ΔG)
	updated_flux_bounds = update_flux_bounds_directionality(MODEL,flux_bounds)

	# show -
	nothing
end

# ╔═╡ bd8f1a4a-61f8-406a-9015-c7841e0cc8d9
begin

	# compute the optimal flux -
	result = calculate_optimal_flux_distribution(S, updated_flux_bounds, species_bounds, c);

	# get the open extent vector -
	ϵ_dot = result.calculated_flux_array

	# did this converge?
	with_terminal() do

		# get exit/status information from the solver -
		exit_flag = result.exit_flag
		status_flag = result.status_flag

		# display -
		println("exit_flag = 0: $(exit_flag==0) and status_flag = 5: $(status_flag == 5)")
	end
end

# ╔═╡ d29277d6-8d05-4cc7-b38f-313c52a3f752
result

# ╔═╡ 6068c779-9367-48d6-a134-b84d3d760ae8
# compute the mol flow rate out of the device -
n_dot_output = (n_dot_input_stream_1 + n_dot_input_stream_2 + S*ϵ_dot);

# ╔═╡ 05e503a8-10fb-45ce-bd65-78c5fe11dc60
with_terminal() do

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

	# what are the reaction strings? -> we can get these from the MODEL object
	reaction_strings = MODEL[:reactions][!,:reaction_markup]

	# populate the state table -
	for reaction_index = 1:ℛ
		flux_table[reaction_index,1] = reaction_index
		flux_table[reaction_index,2] = reaction_strings[reaction_index]
		flux_table[reaction_index,3] = flux_bounds[reaction_index,1]
		flux_table[reaction_index,4] = flux_bounds[reaction_index,2]

		# clean up the display -
		tmp_value = abs(ϵ_dot[reaction_index])
		flux_table[reaction_index,5] = tmp_value < 1e-6 ? 0.0 : ϵ_dot[reaction_index]
	end

	# header row -
	flux_table_header_row = (["i","R","ϵ₁_dot LB", "ϵ₁_dot UB", "ϵᵢ_dot"],
		["","","mmol/hr", "mmol/hr", "mmol/hr"]);

	# write the table -
	pretty_table(flux_table; header=flux_table_header_row)
end

# ╔═╡ bcab8188-9820-4bbc-9abb-5d43df9d5479
with_terminal() do

	# what are the compound names and code strings? -> we can get these from the MODEL object
	compound_name_strings = MODEL[:compounds][!,:compound_name]
	compound_id_strings = MODEL[:compounds][!,:compound_id]

	# how many molecules are in the state array?
	ℳ_local = length(compound_id_strings)

	# initialize some storage -
	state_table = Array{Any,2}(undef,ℳ_local,7)

	# get the uptake array from the result -
	uptake_array = result.uptake_array

	# populate the state table -
	for compound_index = 1:ℳ_local
		state_table[compound_index,1] = compound_index
		state_table[compound_index,2] = compound_name_strings[compound_index]
		state_table[compound_index,3] = compound_id_strings[compound_index]
		state_table[compound_index,4] = n_dot_input_stream_1[compound_index]
		state_table[compound_index,5] = n_dot_input_stream_2[compound_index]

		# for display -
		tmp_value = abs(n_dot_output[compound_index])
		state_table[compound_index,6] = (tmp_value) <= 1e-6 ? 0.0 : n_dot_output[compound_index]

		# show the Δ -
		tmp_value = abs(uptake_array[compound_index])
		state_table[compound_index,7] = (tmp_value) <= 1e-6 ? 0.0 : uptake_array[compound_index]
	end

	# header row -
	state_table_header_row = (["i","name","id","n₁_dot", "n₂_dot", "n₃_dot","Δ"],
		["","","","mmol/hr", "mmol/hr", "mmol/hr", "mmol/hr"]);

	# write the table -
	pretty_table(state_table; header=state_table_header_row)
end

# ╔═╡ 5f0bd4a2-4f86-11ec-3402-0563716ffc85
html"""
<style>
main {
    max-width: 1200px;
    width: 85%;
    margin: auto;
    font-family: "Roboto, monospace";
}

a {
    color: blue;
    text-decoration: none;
}

.H1 {
    padding: 0px 30px;
}
</style>"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"

[compat]
BSON = "~0.3.4"
DataFrames = "~1.2.2"
GLPK = "~0.15.1"
PlutoUI = "~0.7.20"
PrettyTables = "~1.2.3"
"""

# ╔═╡ 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"

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

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

[[BSON]]
git-tree-sha1 = "ebcd6e22d69f21249b7b8668351ebf42d6dc87a1"
uuid = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
version = "0.3.4"

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

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

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

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

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

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

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

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

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

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

[[DataFrames]]
deps = ["Compat", "DataAPI", "Future", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrettyTables", "Printf", "REPL", "Reexport", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
git-tree-sha1 = "d785f42445b63fc86caa08bb9a9351008be9b765"
uuid = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
version = "1.2.2"

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

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

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

[[Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[GLPK]]
deps = ["BinaryProvider", "CEnum", "GLPK_jll", "Libdl", "MathOptInterface"]
git-tree-sha1 = "8c9b0ce6d476800c7a7e8cc0db5a60e241cf107a"
uuid = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
version = "0.15.1"

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

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

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

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

[[InvertedIndices]]
git-tree-sha1 = "bee5f1ef5bf65df56bdd2e40447590b272a5471f"
uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
version = "1.1.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"

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

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

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

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

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

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

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

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

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

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

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

[[MutableArithmetics]]
deps = ["LinearAlgebra", "SparseArrays", "Test"]
git-tree-sha1 = "7bb6853d9afec54019c1397c6eb610b9b9a19525"
uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
version = "0.3.1"

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

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

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

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

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

[[PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "db3a23166af8aebf4db5ef87ac5b00d36eb771e2"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.4.0"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

# ╔═╡ Cell order:
# ╟─7d3b5c92-152f-416a-b82c-82682e61ab9d
# ╟─762e417b-b3f3-4bc0-9c9c-f777ae273303
# ╠═bbfe850c-1f2c-4efb-9a9c-61bb2261370c
# ╟─a1507115-b73d-482a-b66c-734f58463e49
# ╠═39b2b92f-4d6b-4f8c-b840-da2fb0dc79e2
# ╟─fe7efbdb-237c-4ac5-9e40-06df4264266c
# ╠═537e0c28-b2d7-4d29-9592-0bb35747cd81
# ╠═bd8f1a4a-61f8-406a-9015-c7841e0cc8d9
# ╠═d29277d6-8d05-4cc7-b38f-313c52a3f752
# ╠═6068c779-9367-48d6-a134-b84d3d760ae8
# ╠═05e503a8-10fb-45ce-bd65-78c5fe11dc60
# ╠═bcab8188-9820-4bbc-9abb-5d43df9d5479
# ╠═f701a8a6-3a30-4d86-9018-f201276f9369
# ╟─5f0bd4a2-4f86-11ec-3402-0563716ffc85
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using ManifoldsBase
using ManifoldsBase:
    @manifold_element_forwards, @manifold_vector_forwards, @default_manifold_fallbacks
using ManifoldsBase: DefaultManifold, AbstractNumbers, RealNumbers, ComplexNumbers
import ManifoldsBase:
    number_eltype,
    check_point,
    distance,
    embed!,
    exp!,
    inner,
    isapprox,
    log!,
    parallel_transport_to!,
    retract!,
    inverse_retract!
import Base: angle, convert
using LinearAlgebra
using DoubleFloats
using ForwardDiff
using ReverseDiff
using StaticArrays
using Test

struct CustomDefinedRetraction <: ManifoldsBase.AbstractRetractionMethod end
struct CustomUndefinedRetraction <: ManifoldsBase.AbstractRetractionMethod end
struct CustomDefinedInverseRetraction <: ManifoldsBase.AbstractInverseRetractionMethod end

struct DefaultPoint{T} <: AbstractManifoldPoint
    value::T
end
DefaultPoint(v::T) where {T} = DefaultPoint{T}(v)
convert(::Type{DefaultPoint{T}}, v::T) where {T} = DefaultPoint(v)
Base.size(p::DefaultPoint) = size(p.value)
Base.eltype(v::DefaultPoint) = eltype(v.value)

struct DefaultTVector{T} <: TVector
    value::T
end
DefaultTVector(v::T) where {T} = DefaultTVector{T}(v)
Base.size(X::DefaultTVector) = size(X.value)
Base.eltype(X::DefaultTVector) = eltype(X.value)

ManifoldsBase.@manifold_element_forwards DefaultPoint value
ManifoldsBase.@manifold_vector_forwards DefaultTVector value
ManifoldsBase.@default_manifold_fallbacks ManifoldsBase.DefaultManifold DefaultPoint DefaultTVector value value

function ManifoldsBase._injectivity_radius(::DefaultManifold, ::CustomDefinedRetraction)
    return 10.0
end
function ManifoldsBase._retract(M::DefaultManifold, p, X, ::CustomDefinedRetraction)
    return retract_custom(M, p, X)
end
function retract_custom(::DefaultManifold, p::DefaultPoint, X::DefaultTVector)
    return DefaultPoint(p.value + X.value)
end
function ManifoldsBase._inverse_retract(
    M::DefaultManifold,
    p,
    q,
    ::CustomDefinedInverseRetraction,
)
    return inverse_retract_custom(M, p, q)
end
function inverse_retract_custom(::DefaultManifold, p::DefaultPoint, q::DefaultPoint)
    return DefaultTVector(q.value - p.value)
end
struct MatrixVectorTransport{T} <: AbstractVector{T}
    m::Matrix{T}
end
# dummy retractions, inverse retracions for fallback tests - mutating should be enough
ManifoldsBase.retract_polar!(::DefaultManifold, q, p, X) = (q .= p .+ X)
ManifoldsBase.retract_project!(::DefaultManifold, q, p, X) = (q .= p .+ X)
ManifoldsBase.retract_qr!(::DefaultManifold, q, p, X) = (q .= p .+ X)
function ManifoldsBase.retract_exp_ode!(
    ::DefaultManifold,
    q,
    p,
    X,
    m::AbstractRetractionMethod,
    B::ManifoldsBase.AbstractBasis,
)
    return (q .= p .+ X)
end
ManifoldsBase.retract_pade!(::DefaultManifold, q, p, X, i) = (q .= p .+ X)
ManifoldsBase.retract_softmax!(::DefaultManifold, q, p, X) = (q .= p .+ X)
ManifoldsBase.get_embedding(M::DefaultManifold) = M # dummy embedding
ManifoldsBase.inverse_retract_polar!(::DefaultManifold, Y, p, q) = (Y .= q .- p)
ManifoldsBase.inverse_retract_project!(::DefaultManifold, Y, p, q) = (Y .= q .- p)
ManifoldsBase.inverse_retract_qr!(::DefaultManifold, Y, p, q) = (Y .= q .- p)
ManifoldsBase.inverse_retract_softmax!(::DefaultManifold, Y, p, q) = (Y .= q .- p)
function ManifoldsBase.inverse_retract_nlsolve!(
    ::DefaultManifold,
    Y,
    p,
    q,
    m::NLSolveInverseRetraction,
)
    return (Y .= q .- p)
end
function ManifoldsBase.vector_transport_along_project!(
    ::DefaultManifold,
    Y,
    p,
    X,
    c::AbstractVector,
)
    return (Y .= X)
end

Base.getindex(x::MatrixVectorTransport, i) = x.m[:, i]

Base.size(x::MatrixVectorTransport) = (size(x.m, 2),)

@testset "Testing Default (Euclidean)" begin
    M = ManifoldsBase.DefaultManifold(3)
    types = [
        Vector{Float64},
        SizedVector{3,Float64},
        MVector{3,Float64},
        Vector{Float32},
        SizedVector{3,Float32},
        MVector{3,Float32},
        Vector{Double64},
        MVector{3,Double64},
        SizedVector{3,Double64},
        DefaultPoint{Vector{Float64}},
    ]

    @test repr(M) == "DefaultManifold(3; field = ℝ)"
    @test isa(manifold_dimension(M), Integer)
    @test manifold_dimension(M) ≥ 0
    @test base_manifold(M) == M
    @test number_system(M) == ManifoldsBase.ℝ
    @test ManifoldsBase.representation_size(M) == (3,)

    p = zeros(3)
    m = PolarRetraction()
    @test injectivity_radius(M) == Inf
    @test injectivity_radius(M, p) == Inf
    @test injectivity_radius(M, m) == Inf
    @test injectivity_radius(M, p, m) == Inf
    @test default_retraction_method(M) == ExponentialRetraction()
    @test default_inverse_retraction_method(M) == LogarithmicInverseRetraction()

    rm = ManifoldsBase.ExponentialRetraction()
    irm = ManifoldsBase.LogarithmicInverseRetraction()

    # Representation sizes not equal
    @test ManifoldsBase.check_size(M, zeros(3, 3)) isa DomainError
    @test ManifoldsBase.check_size(M, zeros(3), zeros(3, 3)) isa DomainError
    rm2 = CustomDefinedRetraction()
    rm3 = CustomUndefinedRetraction()

    for T in types
        @testset "Type $T" begin
            pts = convert.(Ref(T), [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]])

            @test injectivity_radius(M, pts[1]) == Inf
            @test injectivity_radius(M, pts[1], rm) == Inf
            @test injectivity_radius(M, rm) == Inf
            @test injectivity_radius(M, rm2) == 10
            @test injectivity_radius(M, pts[1], rm2) == 10

            tv1 = log(M, pts[1], pts[2])

            for pt in pts
                @test is_point(M, pt)
            end
            @test is_vector(M, pts[1], tv1; atol = eps(eltype(pts[1])))

            tv2 = log(M, pts[2], pts[1])
            @test isapprox(M, pts[2], exp(M, pts[1], tv1))
            @test isapprox(M, pts[1], exp(M, pts[1], tv1, 0))
            @test isapprox(M, pts[2], exp(M, pts[1], tv1, 1))
            @test isapprox(M, pts[1], exp(M, pts[2], tv2))
            @test is_point(M, retract(M, pts[1], tv1))
            @test isapprox(M, pts[1], retract(M, pts[1], tv1, 0))

            @test is_point(M, retract(M, pts[1], tv1, rm))
            @test isapprox(M, pts[1], retract(M, pts[1], tv1, 0, rm))

            new_pt = exp(M, pts[1], tv1)
            retract!(M, new_pt, pts[1], tv1)
            @test is_point(M, new_pt)
            for p in pts
                @test isapprox(M, p, zero_vector(M, p), log(M, p, p); atol = eps(eltype(p)))
                @test isapprox(
                    M,
                    p,
                    zero_vector(M, p),
                    inverse_retract(M, p, p);
                    atol = eps(eltype(p)),
                )
                @test isapprox(
                    M,
                    p,
                    zero_vector(M, p),
                    inverse_retract(M, p, p, irm);
                    atol = eps(eltype(p)),
                )
            end
            zero_vector!(M, tv1, pts[1])
            @test isapprox(M, pts[1], tv1, zero_vector(M, pts[1]))
            log!(M, tv1, pts[1], pts[2])
            @test norm(M, pts[1], tv1) ≈ sqrt(inner(M, pts[1], tv1, tv1))

            @test isapprox(M, exp(M, pts[1], tv1, 1), pts[2])
            @test isapprox(M, exp(M, pts[1], tv1, 0), pts[1])

            @test distance(M, pts[1], pts[2]) ≈ norm(M, pts[1], tv1)

            @test mid_point(M, pts[1], pts[2]) == convert(T, [0.5, 0.5, 0.0])
            midp = allocate(pts[1])
            @test mid_point!(M, midp, pts[1], pts[2]) === midp
            @test midp == convert(T, [0.5, 0.5, 0.0])

            q = copy(M, pts[1])
            Ts = [0.0, 1.0 / 2, 1.0]
            @testset "Geodesic interface test" begin
                @test isapprox(M, geodesic(M, pts[1], tv1)(0.0), pts[1])
                @test isapprox(M, geodesic(M, pts[1], tv1)(1.0), pts[2])
                g! = geodesic!(M, pts[1], tv1)
                g!(q, 0.0)
                isapprox(M, q, pts[1])
                g!(q, 0.0)
                isapprox(M, q, pts[2])
                geodesic!(M, q, pts[1], tv1, 1.0 / 2)
                isapprox(M, q, midp)
                @test isapprox(M, geodesic(M, pts[1], tv1, 1.0), pts[2])
                @test isapprox(M, geodesic(M, pts[1], tv1, 1.0 / 2), midp)
                @test isapprox(M, shortest_geodesic(M, pts[1], pts[2])(0.0), pts[1])
                @test isapprox(M, shortest_geodesic(M, pts[1], pts[2])(1.0), pts[2])
                sg! = shortest_geodesic!(M, pts[1], pts[2])
                sg!(q, 0.0)
                isapprox(M, q, pts[1])
                sg!(q, 1.0)
                isapprox(M, q, pts[2])
                @test isapprox(M, shortest_geodesic(M, pts[1], pts[2], 0.0), pts[1])
                @test isapprox(M, shortest_geodesic(M, pts[1], pts[2], 1.0), pts[2])
                shortest_geodesic!(M, q, pts[1], pts[2], 0.5)
                isapprox(M, q, midp)
                @test all(
                    isapprox.(Ref(M), geodesic(M, pts[1], tv1, Ts), [pts[1], midp, pts[2]]),
                )
                @test all(
                    isapprox.(
                        Ref(M),
                        shortest_geodesic(M, pts[1], pts[2], Ts),
                        [pts[1], midp, pts[2]],
                    ),
                )
                Q = [copy(M, q), copy(M, q), copy(M, q)]
                geodesic!(M, Q, pts[1], tv1, Ts)
                @test all(isapprox.(Ref(M), Q, [pts[1], midp, pts[2]]))
                shortest_geodesic!(M, Q, pts[1], pts[2], Ts)
                @test all(isapprox.(Ref(M), Q, [pts[1], midp, pts[2]]))
            end

            @testset "basic linear algebra in tangent space" begin
                @test isapprox(
                    M,
                    pts[1],
                    0 * tv1,
                    zero_vector(M, pts[1]);
                    atol = eps(eltype(pts[1])),
                )
                @test isapprox(M, pts[1], 2 * tv1, tv1 + tv1)
                @test isapprox(M, pts[1], 0 * tv1, tv1 - tv1)
                @test isapprox(M, pts[1], (-1) * tv1, -tv1)
            end

            @testset "Hat and vee in the tangent space" begin
                X = log(M, pts[1], pts[2])
                a = vee(M, pts[1], X)
                b = similar(a)
                vee!(M, b, pts[1], X)
                Y = hat(M, pts[1], a)
                Z = similar(Y)
                hat!(M, Z, pts[1], a)
                @test a == b
                @test X == Y
                @test Z == X
                @test a == ((T <: DefaultPoint) ? vec(X.value) : vec(X))
            end

            @testset "broadcasted linear algebra in tangent space" begin
                @test isapprox(M, pts[1], 3 * tv1, 2 .* tv1 .+ tv1)
                @test isapprox(M, pts[1], -tv1, tv1 .- 2 .* tv1)
                @test isapprox(M, pts[1], -tv1, .-tv1)
                v = similar(tv1)
                v .= 2 .* tv1 .+ tv1
                @test isapprox(M, pts[1], v, 3 * tv1)
            end

            @testset "project test" begin
                # point
                @test isapprox(M, pts[1], project(M, pts[1]))
                pt = similar(pts[1])
                project!(M, pt, pts[1])
                @test isapprox(M, pt, pts[1])

                @test isapprox(M, pts[1], embed(M, pts[1]))
                pt = similar(pts[1])
                embed!(M, pt, pts[1])
                @test isapprox(M, pt, pts[1])

                # tangents
                @test isapprox(M, pts[1], tv1, project(M, pts[1], tv1))
                tv = similar(tv1)
                project!(M, tv, pts[1], tv1)
                @test isapprox(M, pts[1], tv, tv1)

                @test isapprox(M, pts[1], tv1, embed(M, pts[1], tv1))
                tv = similar(tv1)
                embed!(M, tv, pts[1], tv1)
                @test isapprox(M, pts[1], tv, tv1)
            end

            @testset "vector transport" begin
                # test constructor
                @test default_vector_transport_method(M) == ParallelTransport()
                X1 = log(M, pts[1], pts[2])
                X2 = log(M, pts[1], pts[3])
                X1t1 = vector_transport_to(M, pts[1], X1, pts[3])
                X1t2 = zero(X1t1)
                vector_transport_to!(M, X1t2, pts[1], X1, X2, ProjectionTransport())
                X1t3 = vector_transport_direction(M, pts[1], X1, X2)
                @test ManifoldsBase.is_vector(M, pts[3], X1t1)
                @test ManifoldsBase.is_vector(M, pts[3], X1t3)
                @test isapprox(M, pts[3], X1t1, X1t3)
                # along a `Vector` of points
                c = [pts[1], pts[1]]
                X1t4 = vector_transport_along(M, pts[1], X1, c)
                @test isapprox(M, pts[1], X1, X1t4)
                X1t5 = allocate(X1)
                vector_transport_along!(M, X1t5, pts[1], X1, c)
                @test isapprox(M, pts[1], X1, X1t5)
                # transport along more than one interims point
                @test vector_transport_along(M, pts[1], X1, pts[2:3]) == X1
                X1t6 = allocate(X1)
                vector_transport_along!(M, X1t6, pts[1], X1, pts[2:3])
                @test isapprox(M, pts[1], X1, X1t6)
                # along a custom type of points
                if T <: DefaultPoint
                    S = eltype(pts[1].value)
                    mat = reshape(pts[1].value, length(pts[1].value), 1)
                else
                    S = eltype(pts[1])
                    mat = reshape(pts[1], length(pts[1]), 1)
                end
                c2 = MatrixVectorTransport{S}(mat)
                X1t4c2 = vector_transport_along(M, pts[1], X1, c2)
                @test isapprox(M, pts[1], X1, X1t4c2)
                X1t5c2 = allocate(X1)
                vector_transport_along!(M, X1t5c2, pts[1], X1, c2)
                @test isapprox(M, pts[1], X1, X1t5c2)
                # On Euclidean Space Schild & Pole are identity
                @test vector_transport_to(
                    M,
                    pts[1],
                    X2,
                    pts[2],
                    SchildsLadderTransport(),
                ) == X2
                @test vector_transport_to(M, pts[1], X2, pts[2], PoleLadderTransport()) ==
                      X2
                @test vector_transport_to(
                    M,
                    pts[1],
                    X2,
                    pts[2],
                    ScaledVectorTransport(ParallelTransport()),
                ) == X2

                # along is also the identity
                c = [
                    mid_point(M, pts[1], pts[2]),
                    pts[2],
                    mid_point(M, pts[2], pts[3]),
                    pts[3],
                ]
                Xtmp = allocate(X2)
                @test vector_transport_along(M, pts[1], X2, c, SchildsLadderTransport()) ==
                      X2
                @test vector_transport_along!(
                    M,
                    Xtmp,
                    pts[1],
                    X2,
                    c,
                    SchildsLadderTransport(),
                ) == X2
                @test vector_transport_along(M, pts[1], X2, c, PoleLadderTransport()) == X2
                @test vector_transport_along!(
                    M,
                    Xtmp,
                    pts[1],
                    X2,
                    c,
                    PoleLadderTransport(),
                ) == X2
                @test vector_transport_along(M, pts[1], X2, c, ParallelTransport()) == X2
                @test vector_transport_along!(
                    M,
                    Xtmp,
                    pts[1],
                    X2,
                    c,
                    ParallelTransport(),
                ) == X2
                @test ManifoldsBase._vector_transport_along!(
                    M,
                    Xtmp,
                    pts[1],
                    X2,
                    c,
                    ParallelTransport(),
                ) == X2
                @test vector_transport_along(M, pts[1], X2, c, ProjectionTransport()) == X2
                @test vector_transport_along!(
                    M,
                    Xtmp,
                    pts[1],
                    X2,
                    c,
                    ProjectionTransport(),
                ) == X2
                # check mutating ones with defaults
                p = allocate(pts[1])
                ManifoldsBase.pole_ladder!(M, p, pts[1], pts[2], pts[3])
                # -log_p3 p == log_p1 p2
                @test isapprox(M, pts[3], -log(M, pts[3], p), log(M, pts[1], pts[2]))
                ManifoldsBase.schilds_ladder!(M, p, pts[1], pts[2], pts[3])
                @test isapprox(M, pts[3], log(M, pts[3], p), log(M, pts[1], pts[2]))

                @test repr(ParallelTransport()) == "ParallelTransport()"
                @test repr(ScaledVectorTransport(ParallelTransport())) ==
                      "ScaledVectorTransport(ParallelTransport())"
            end

            @testset "ForwardDiff support" begin
                exp_f(t) = distance(M, pts[1], exp(M, pts[1], t * tv1))
                d12 = distance(M, pts[1], pts[2])
                for t in 0.1:0.1:0.9
                    @test d12 ≈ ForwardDiff.derivative(exp_f, t)
                end

                retract_f(t) = distance(M, pts[1], retract(M, pts[1], t * tv1))
                for t in 0.1:0.1:0.9
                    @test ForwardDiff.derivative(retract_f, t) ≥ 0
                end
            end

            isa(pts[1], Union{Vector,SizedVector}) && @testset "ReverseDiff support" begin
                exp_f(t) = distance(M, pts[1], exp(M, pts[1], t[1] * tv1))
                d12 = distance(M, pts[1], pts[2])
                for t in 0.1:0.1:0.9
                    @test d12 ≈ ReverseDiff.gradient(exp_f, [t])[1]
                end

                retract_f(t) = distance(M, pts[1], retract(M, pts[1], t[1] * tv1))
                for t in 0.1:0.1:0.9
                    @test ReverseDiff.gradient(retract_f, [t])[1] ≥ 0
                end
            end
        end
    end

    @testset "mid_point on 0-index arrays" begin
        M = ManifoldsBase.DefaultManifold(1)
        p1 = fill(0.0)
        p2 = fill(1.0)
        @test isapprox(M, fill(0.5), mid_point(M, p1, p2))
    end

    @testset "Retraction" begin
        a = NLSolveInverseRetraction(ExponentialRetraction())
        @test a.retraction isa ExponentialRetraction

    end

    @testset "copy of points and vectors" begin
        M = ManifoldsBase.DefaultManifold(2)
        for (p, X) in (
            ([2.0, 3.0], [4.0, 5.0]),
            (DefaultPoint([2.0, 3.0]), DefaultTVector([4.0, 5.0])),
        )
            q = similar(p)
            copyto!(M, q, p)
            @test p == q
            r = copy(M, p)
            @test r == p
            Y = similar(X)
            copyto!(M, Y, p, X)
            @test Y == X
            Z = copy(M, p, X)
            @test Z == X
        end

        p1 = DefaultPoint([2.0, 3.0])
        p2 = copy(p1)
        @test (p1 == p2) && (p1 !== p2)
    end
    @testset "further vector and point automatic forwards" begin
        M = ManifoldsBase.DefaultManifold(3)
        p = DefaultPoint([1.0, 0.0, 0.0])
        q = DefaultPoint([0.0, 0.0, 0.0])
        X = DefaultTVector([0.0, 1.0, 0.0])
        Y = DefaultTVector([1.0, 0.0, 0.0])
        @test angle(M, p, X, Y) ≈ π / 2
        @test inverse_retract(M, p, q, LogarithmicInverseRetraction()) == -Y
        @test retract(M, q, Y, CustomDefinedRetraction()) == p
        @test retract(M, q, Y, ExponentialRetraction()) == p
        # rest not implemented - so they also fall back even onto mutating
        Z = similar(Y)
        r = similar(p)
        # test passthrough using the dummy implementations
        for retr in [PolarRetraction, ProjectionRetraction, QRRetraction, SoftmaxRetraction]
            @test retract(M, q, Y, retr()) == DefaultPoint(q.value + Y.value)
            @test retract!(M, r, q, Y, retr()) == DefaultPoint(q.value + Y.value)
        end
        @test retract(
            M,
            q,
            Y,
            ODEExponentialRetraction(PolarRetraction(), DefaultBasis()),
        ) == DefaultPoint(q.value + Y.value)
        @test retract!(
            M,
            r,
            q,
            Y,
            ODEExponentialRetraction(PolarRetraction(), DefaultBasis()),
        ) == DefaultPoint(q.value + Y.value)
        @test retract(M, q, Y, PadeRetraction(2)) == DefaultPoint(q.value + Y.value)
        @test retract!(M, r, q, Y, PadeRetraction(2)) == DefaultPoint(q.value + Y.value)
        @test retract!(M, r, q, Y, EmbeddedRetraction(ExponentialRetraction())) ==
              DefaultPoint(q.value + Y.value)
        @test retract(M, q, Y, EmbeddedRetraction(ExponentialRetraction())) ==
              DefaultPoint(q.value + Y.value)
        p2 = allocate(p, eltype(p.value), size(p.value))
        @test size(p2.value) == size(p.value)
        X2 = allocate(X, eltype(X.value), size(X.value))
        @test size(X2.value) == size(X.value)
        X3 = ManifoldsBase.allocate_result(M, log, p, q)
        @test log!(M, X3, p, q) == log(M, p, q)
        @test X3 == log(M, p, q)
        @test log!(M, X3, p, q) == log(M, p, q)
        @test X3 == log(M, p, q)
        @test inverse_retract(M, p, q, CustomDefinedInverseRetraction()) == -Y
        X4 = ManifoldsBase.allocate_result(M, inverse_retract, p, q)
        @test inverse_retract!(M, X4, p, q) == inverse_retract(M, p, q)
        @test X4 == inverse_retract(M, p, q)
        # rest not implemented but check passthrough
        for r in [
            PolarInverseRetraction,
            ProjectionInverseRetraction,
            QRInverseRetraction,
            SoftmaxInverseRetraction,
        ]
            @test inverse_retract(M, q, p, r()) == DefaultTVector(p.value - q.value)
            @test inverse_retract!(M, Z, q, p, r()) == DefaultTVector(p.value - q.value)
        end
        @test inverse_retract(
            M,
            q,
            p,
            EmbeddedInverseRetraction(LogarithmicInverseRetraction()),
        ) == DefaultTVector(p.value - q.value)
        @test inverse_retract(M, q, p, NLSolveInverseRetraction(ExponentialRetraction())) ==
              DefaultTVector(p.value - q.value)
        @test inverse_retract!(
            M,
            Z,
            q,
            p,
            EmbeddedInverseRetraction(LogarithmicInverseRetraction()),
        ) == DefaultTVector(p.value - q.value)
        @test inverse_retract!(
            M,
            Z,
            q,
            p,
            NLSolveInverseRetraction(ExponentialRetraction()),
        ) == DefaultTVector(p.value - q.value)
        c = ManifoldsBase.allocate_coordinates(M, p, Float64, manifold_dimension(M))
        @test c isa Vector
        @test length(c) == 3
        @test 2.0 \ X == DefaultTVector(2.0 \ X.value)
        @test X + Y == DefaultTVector(X.value + Y.value)
        @test +X == X
        @test (Y .= X) === Y
        # vector transport pass through
        @test vector_transport_to(M, p, X, q, ProjectionTransport()) == X
        @test vector_transport_direction(M, p, X, X, ProjectionTransport()) == X
        @test vector_transport_to!(M, Y, p, X, q, ProjectionTransport()) == X
        @test vector_transport_direction!(M, Y, p, X, X, ProjectionTransport()) == X
        @test vector_transport_along(M, p, X, [p], ProjectionTransport()) == X
        @test vector_transport_along!(M, Z, p, X, [p], ProjectionTransport()) == X
        @test vector_transport_to(M, p, X, :q, ProjectionTransport()) == X
        @test parallel_transport_to(M, p, X, q) == X
        @test parallel_transport_direction(M, p, X, X) == X
        @test parallel_transport_along(M, p, X, []) == X
        @test parallel_transport_to!(M, Y, p, X, q) == X
        @test parallel_transport_direction!(M, Y, p, X, X) == X
        @test parallel_transport_along!(M, Y, p, X, []) == X
    end
    @testset "DefaultManifold  and ONB" begin
        M = ManifoldsBase.DefaultManifold(3)
        p = [1.0, 0.0, 0.0]
        CB = get_basis(M, p, DefaultOrthonormalBasis())
        @test CB.data == [[1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]]
    end
    @testset "Show methods" begin
        @test repr(CayleyRetraction()) == "CayleyRetraction()"
        @test repr(PadeRetraction(2)) == "PadeRetraction(2)"
    end
    @testset "Further TestArrayRepresentation" begin
        M = ManifoldsBase.DefaultManifold(3)
        p = [1.0, 0.0, 0.0]
        X = [1.0, 0.0, 0.0]
        @test is_point(M, p, true)
        @test is_vector(M, p, X, true)
        pF = [1.0, 0.0]
        XF = [0.0, 0.0]
        m = ExponentialRetraction()
        @test_throws DomainError is_point(M, pF, true)
        @test_throws DomainError is_vector(M, p, XF, true)
        @test_throws DomainError is_vector(M, pF, XF, true; check_point = true)
        @test injectivity_radius(M) == Inf
        @test injectivity_radius(M, p) == Inf
        @test injectivity_radius(M, p, m) == Inf
        @test injectivity_radius(M, m) == Inf
    end
end

ext: jl


text:
function >>(left::Expr, right::Expr)
    quote
        $left
        $right
    end
end

function mapcatexpr(f, xs)
    body = quote end
    for x in xs
        next = f(x)
        body = body >> next
    end
    prettify(body)
end

function brackets(i)
    var -> :($var[$i])
end

function condname(var)
    Symbol("cond__" * string(var))
end

function footer(observe)
    quote
        loglikelihood += current.loglikelihood
        if loglikelihood == -Inf
            $observe
            return traceroot
        end
        previous = current
        current = current.next
    end
end

function varbody(model, dist, var, access)
    invar = access(var)
    random = !haskey(model.conditionings, var)
    value = random ? :(rand($dist)) : access(condname(var))

    quote
        if current isa Leaf
            $invar = $value
            current = extend!(previous, $(quot(var)), $dist, $invar, $random)

            # TODO: we could move this if out of the runtime to avoid checking this
            # on each node
            if traceroot isa Leaf
                traceroot = current
            end
        else
            $invar = current.value
        end
    end
end


text:
# SIF portion

sif_fname = "/home/psomkuti/s5p_sif_2018-2019_gap_filled.h5"
spectra_fname = "/home/psomkuti/jLRS/data/photosystem_spectra.txt"

ps_spectra_wl = Float32[]
ps_spectra_phiI = Float32[]
ps_spectra_phiII = Float32[]
ps_spectra_raw = readlines(spectra_fname)

for line in ps_spectra_raw
    a, b, c = (p -> parse(Float32, p)).(split(line, ","))
    push!(ps_spectra_wl, a)
    push!(ps_spectra_phiI, b)
    push!(ps_spectra_phiII, c)
end

function create_TROPOMI_SIF_from_locbounds(lon_min,
                                           lat_min,
                                           lon_max,
                                           lat_max)

    # SIF data array: month, lon, lat
    sif_h5 = h5open(sif_fname, "r")

    lon_bounds = sif_h5["lon"][:]
    lat_bounds = sif_h5["lat"][:]

    idx_lon_min = searchsortedfirst(lon_bounds, lon_min) - 1
    idx_lon_max = searchsortedfirst(lon_bounds, lon_max)

    idx_lat_min = searchsortedfirst(lat_bounds, lat_min) - 1
    idx_lat_max = searchsortedfirst(lat_bounds, lat_max)

    # Load SIF subset and switch dimensions so they're
    # time, lon, lat ->
    # lat, lon, time
    sif_subset = permutedims(
        sif_h5["sif"][:,idx_lon_min:idx_lon_max,idx_lat_min:idx_lat_max],
        [3,2,1]
    )

    dates = [DateTime(1999, i, 15) for i in 1:12]

    return TROPOMISIFData(
        sif_subset,
        lon_bounds[idx_lon_min:idx_lon_max],
        lat_bounds[idx_lat_min:idx_lat_max],
        dates,
        NaN
    )

end

function calculate_tropomi_sif(loctime::GeolocationTime,
                               sza::Real,
                               sif_data::TROPOMISIFData)

    idx_x = searchsortedfirst(sif_data.lon_bounds, loctime.loc.lon) - 1
    idx_y = searchsortedfirst(sif_data.lat_bounds, loctime.loc.lat) - 1

    newtime = DateTime(
        Dates.year(sif_data.time_bounds[1]),
        Dates.month(loctime.time),
        Dates.day(loctime.time),
        Dates.hour(loctime.time),
        Dates.minute(loctime.time),
        Dates.second(loctime.time)
    )

    idx_t_left = searchsortedfirst(sif_data.time_bounds, newtime) - 1

    if idx_t_left == 0
        idx_t_left = 12
    end

    if idx_t_left == 12
        idx_t_right = 1
    else
        idx_t_right = idx_t_left + 1
    end

    t_fac = (newtime - sif_data.time_bounds[idx_t_left]) /
        (sif_data.time_bounds[idx_t_right] - sif_data.time_bounds[idx_t_left])

    sif_val = (1.0 - t_fac) * sif_data.data[idx_y, idx_x, idx_t_left] +
        t_fac * sif_data.data[idx_y, idx_x, idx_t_right]

    # Keep in mind: this is SIF at 755nm
    return sif_val * cos(deg2rad(sza))

end

function calculate_fAPAR(NDVI::Real)

    linear_fAPAR = 1.35 * NDVI - 0.32

    if linear_fAPAR < 0.0
        linear_fAPAR = 0.0
    end

    if linear_fAPAR > 1.0
        linear_fAPAR = 1.0
    end

    return linear_fAPAR

end

######
###### J. Johnson's SIF model
###### https://github.com/jenjohnson/johnson-berry-2021-pres
######

function photosynthetic_yields(
    Qin, # input PAR, umol PPFD / m2 / s - you want this to be an array
    Tin, # input leaf temp, C
    Cin, # input mesophyll CO2, ubar
    Oin;  # input atmospheric O2, mbar
    Abs=0.85, # this is fAPAR essentially (?), mol PPFD absorbed by the leaf, PPFD incident / mol absorbed
    alpha_opt="static",
    beta = 0.52, # mol PPFD absorbed by PSII, PPFD / mol absorbed
    CB6F = (350.0 / 300.0) / 1e6, # Cyt b6f density, mol sites per m2
    RUB = (100.0 / 3.6) / 1e6, # Rubisco density, mol sites per m2
    Rds = 0.01, # dark respiration scalar
    Kf = 0.05e9, # Rate constant for fluorescence at PSII and PSI, 1/s
    Kd = 0.55e9, # Rate constant for constitutive heat loss at PSII and PSI, 1/s
    Kp1 = 14.5e9, # Rate constant for photochemistry at PSI, 1/s
    Kp2 = 4.5e9,  # Rate constant for photochemistry at PSII, 1/s
    Kn1 = 14.5e9, # Rate constant for regulated heat loss at PSI, 1/s
    Ku2 = 0.0e9,  # Rate constant for exciton sharing PSII, 1/s
    kq = 300.0, # Cyt b6f kcat for PQH2, mol mol / e- / sites / s
    nl = 0.75,  # ATP per e- in linear flow, ATP/e-
    nc = 1.00,  # ATP per e- in cyclic flow, ATP/e-
    kc = 3.6,   # Rubisco kcat for CO2, mol CO2 mol / sites / s
    ko = 3.6 * 0.27, # Rubisco kcat for O2, mol O2 mol / sites / s
    Kc = 260.0 / 1e6,  # Rubisco Km for CO2, bar
    Ko = 179000.0 / 1e6, # Rubisco Km for O2, bar
    theta1 = 1., # curvature parameter for Aj/Ac transition
    eps1 = 0., # PSI transfer function
    eps2 = 1. # PSII transfer function
)

    @assert (alpha_opt == "static") | (alpha_opt == "dynamic") "Bad parameter!"

    # Inputs
    Q = Qin ./ 1e6 # PAR in mol PPFD / m2 /s
    T = Tin
    C = Cin / 1e6 # mesophyll CO2 partial pressure, bar
    O = Oin / 1e3 # atmospheric O2 partial pressure, O2

    # Calculations

    Vqmax = CB6F * kq           # Maximum Cyt b6f activity, mol e-1 m-2 s-1
    Vcmax = RUB * kc            # Maximum Rubisco activity, mol CO2 m-2 s-1
    Rd = Vcmax * Rds            # Mitochondrial respiration, mol CO2 m-2 s-1
    S = (kc / Kc) * (Ko ./ ko)   # Rubisco specificity for CO2/O2, dimensionless
    gammas = O / (2 * S)        # CO2 compensation point in the absence of Rd, bar
    eta = (1.0 - (nl / nc) + (3.0 + 7.0 * gammas / C) / ((4.0 + 8.0 * gammas / C) * nc)) # PS I/II ETR
    phi1P_max = Kp1 / (Kp1 + Kd + Kf) # Maximum photochemical yield PS I

    if alpha_opt == "static"
        a2 = Abs .* beta
        a1 = Abs .- a2
    elseif alpha_opt == "dynamic"
        solve_xcs(Abs, CB6F, Kd, Kf, Kp2, Ku2, Q, eta, kq, phi1P_max) =
            @. (-sqrt((Kd + Kf + Ku2) * (CB6F^2 * Kd^3 * kq^2 * phi1P_max^2 + CB6F^2 * Kf^3 * kq^2 * phi1P_max^2 + CB6F^2 * Kd * Kp2^2 * eta^2 * kq^2 + CB6F^2 * Kf * Kp2^2 * eta^2 * kq^2 + CB6F^2 * Kp2^2 * Ku2 * eta^2 * kq^2 + CB6F^2 * Kd * Kf^2 * kq^2 * phi1P_max^2 * 3.0 + CB6F^2 * Kd^2 * Kf * kq^2 * phi1P_max^2 * 3.0 + CB6F^2 * Kd * Kp2^2 * kq^2 * phi1P_max^2 + CB6F^2 * Kd^2 * Kp2 * kq^2 * phi1P_max^2 * 2.0 + CB6F^2 * Kf * Kp2^2 * kq^2 * phi1P_max^2 + CB6F^2 * Kf^2 * Kp2 * kq^2 * phi1P_max^2 * 2.0 + CB6F^2 * Kd^2 * Ku2 * kq^2 * phi1P_max^2 + CB6F^2 * Kf^2 * Ku2 * kq^2 * phi1P_max^2 + CB6F^2 * Kp2^2 * Ku2 * kq^2 * phi1P_max^2 + Abs^2 * Kd * Kp2^2 * Q^2 * eta^2 * phi1P_max^2 + Abs^2 * Kf * Kp2^2 * Q^2 * eta^2 * phi1P_max^2 + Abs^2 * Kp2^2 * Ku2 * Q^2 * eta^2 * phi1P_max^2 + CB6F^2 * Kd * Kf * Kp2 * kq^2 * phi1P_max^2 * 4.0 + CB6F^2 * Kd * Kf * Ku2 * kq^2 * phi1P_max^2 * 2.0 + CB6F^2 * Kd * Kp2 * Ku2 * kq^2 * phi1P_max^2 * 2.0 + CB6F^2 * Kf * Kp2 * Ku2 * kq^2 * phi1P_max^2 * 2.0 + CB6F^2 * Kd * Kp2^2 * eta * kq^2 * phi1P_max * 2.0 + CB6F^2 * Kd^2 * Kp2 * eta * kq^2 * phi1P_max * 2.0 + CB6F^2 * Kf * Kp2^2 * eta * kq^2 * phi1P_max * 2.0 + CB6F^2 * Kf^2 * Kp2 * eta * kq^2 * phi1P_max * 2.0 + CB6F^2 * Kp2^2 * Ku2 * eta * kq^2 * phi1P_max * 2.0 + CB6F^2 * Kd * Kf * Kp2 * eta * kq^2 * phi1P_max * 4.0 + CB6F^2 * Kd * Kp2 * Ku2 * eta * kq^2 * phi1P_max * 2.0 + CB6F^2 * Kf * Kp2 * Ku2 * eta * kq^2 * phi1P_max * 2.0 + Abs * CB6F * Kd * Kp2^2 * Q * eta * kq * phi1P_max^2 * 2.0 + Abs * CB6F * Kd * Kp2^2 * Q * eta^2 * kq * phi1P_max * 2.0 + Abs * CB6F * Kd^2 * Kp2 * Q * eta * kq * phi1P_max^2 * 2.0 + Abs * CB6F * Kf * Kp2^2 * Q * eta * kq * phi1P_max^2 * 2.0 + Abs * CB6F * Kf * Kp2^2 * Q * eta^2 * kq * phi1P_max * 2.0 + Abs * CB6F * Kf^2 * Kp2 * Q * eta * kq * phi1P_max^2 * 2.0 - Abs * CB6F * Kp2^2 * Ku2 * Q * eta * kq * phi1P_max^2 * 2.0 + Abs * CB6F * Kp2^2 * Ku2 * Q * eta^2 * kq * phi1P_max * 2.0 + Abs * CB6F * Kd * Kf * Kp2 * Q * eta * kq * phi1P_max^2 * 4.0 + Abs * CB6F * Kd * Kp2 * Ku2 * Q * eta * kq * phi1P_max^2 * 2.0 + Abs * CB6F * Kf * Kp2 * Ku2 * Q * eta * kq * phi1P_max^2 * 2.0)) + CB6F * Kd^2 * kq * phi1P_max + CB6F * Kf^2 * kq * phi1P_max + Abs * Kd^2 * Q * phi1P_max^2 * 2.0 + Abs * Kf^2 * Q * phi1P_max^2  *  2.0 + CB6F * Kd * Kp2 * eta * kq + CB6F * Kf * Kp2 * eta * kq + CB6F * Kp2 * Ku2 * eta * kq + CB6F * Kd * Kf * kq * phi1P_max * 2.0 + CB6F * Kd * Kp2 * kq * phi1P_max + CB6F * Kf * Kp2 * kq * phi1P_max + CB6F * Kd * Ku2 * kq * phi1P_max + CB6F * Kf * Ku2 * kq * phi1P_max + CB6F * Kp2 * Ku2 * kq * phi1P_max + Abs * Kd * Kf * Q * phi1P_max^2 * 4.0 + Abs * Kd * Kp2 * Q * phi1P_max^2 * 2.0 + Abs * Kf * Kp2 * Q * phi1P_max^2  *  2.0 + Abs * Kd * Ku2 * Q * phi1P_max^2 * 2.0 + Abs * Kf * Ku2 * Q * phi1P_max^2 * 2.0 + Abs * Kd * Kp2 * Q * eta * phi1P_max + Abs * Kf * Kp2 * Q * eta * phi1P_max + Abs * Kp2 * Ku2 * Q * eta * phi1P_max) ./ (Q * phi1P_max * (Kd^2 * phi1P_max + Kf^2 * phi1P_max + Kd * Kp2 * eta + Kf * Kp2 * eta + Kp2 * Ku2 * eta + Kd * Kf * phi1P_max * 2.0 + Kd * Kp2 * phi1P_max + Kf * Kp2 * phi1P_max + Kd * Ku2 * phi1P_max + Kf * Ku2 * phi1P_max) * 2.0)
            #(-sqrt((Kd+Kf+Ku2)*(CB6F^2*Kd^3*kq^2*phi1P_max^2+CB6F^2*Kf^3*kq^2*phi1P_max^2+CB6F^2*Kd*Kp2^2*eta^2*kq^2+CB6F^2*Kf*Kp2^2*eta^2*kq^2+CB6F^2*Kp2^2*Ku2*eta^2*kq^2+CB6F^2*Kd*Kf^2*kq^2*phi1P_max^2*3.0+CB6F^2*Kd^2*Kf*kq^2*phi1P_max^2*3.0+CB6F^2*Kd*Kp2^2*kq^2*phi1P_max^2+CB6F^2*Kd^2*Kp2*kq^2*phi1P_max^2*2.0+CB6F^2*Kf*Kp2^2*kq^2*phi1P_max^2+CB6F^2*Kf^2*Kp2*kq^2*phi1P_max^2*2.0+CB6F^2*Kd^2*Ku2*kq^2*phi1P_max^2+CB6F^2*Kf^2*Ku2*kq^2*phi1P_max^2+CB6F^2*Kp2^2*Ku2*kq^2*phi1P_max^2+Abs^2*Kd*Kp2^2*Q^2*eta^2*phi1P_max^2+Abs^2*Kf*Kp2^2*Q^2*eta^2*phi1P_max^2+Abs^2*Kp2^2*Ku2*Q^2*eta^2*phi1P_max^2+CB6F^2*Kd*Kf*Kp2*kq^2*phi1P_max^2*4.0+CB6F^2*Kd*Kf*Ku2*kq^2*phi1P_max^2*2.0+CB6F^2*Kd*Kp2*Ku2*kq^2*phi1P_max^2*2.0+CB6F^2*Kf*Kp2*Ku2*kq^2*phi1P_max^2*2.0+CB6F^2*Kd*Kp2^2*eta*kq^2*phi1P_max*2.0+CB6F^2*Kd^2*Kp2*eta*kq^2*phi1P_max*2.0+CB6F^2*Kf*Kp2^2*eta*kq^2*phi1P_max*2.0+CB6F^2*Kf^2*Kp2*eta*kq^2*phi1P_max*2.0+CB6F^2*Kp2^2*Ku2*eta*kq^2*phi1P_max*2.0+CB6F^2*Kd*Kf*Kp2*eta*kq^2*phi1P_max*4.0+CB6F^2*Kd*Kp2*Ku2*eta*kq^2*phi1P_max*2.0+CB6F^2*Kf*Kp2*Ku2*eta*kq^2*phi1P_max*2.0+Abs*CB6F*Kd*Kp2^2*Q*eta*kq*phi1P_max^2*2.0+Abs*CB6F*Kd*Kp2^2*Q*eta^2*kq*phi1P_max*2.0+Abs*CB6F*Kd^2*Kp2*Q*eta*kq*phi1P_max^2*2.0+Abs*CB6F*Kf*Kp2^2*Q*eta*kq*phi1P_max^2*2.0+Abs*CB6F*Kf*Kp2^2*Q*eta^2*kq*phi1P_max*2.0+Abs*CB6F*Kf^2*Kp2*Q*eta*kq*phi1P_max^2*2.0-Abs*CB6F*Kp2^2*Ku2*Q*eta*kq*phi1P_max^2*2.0+Abs*CB6F*Kp2^2*Ku2*Q*eta^2*kq*phi1P_max*2.0+Abs*CB6F*Kd*Kf*Kp2*Q*eta*kq*phi1P_max^2*4.0+Abs*CB6F*Kd*Kp2*Ku2*Q*eta*kq*phi1P_max^2*2.0+Abs*CB6F*Kf*Kp2*Ku2*Q*eta*kq*phi1P_max^2*2.0))+CB6F*Kd^2*kq*phi1P_max+CB6F*Kf^2*kq*phi1P_max+Abs*Kd^2*Q*phi1P_max^2*2.0+Abs*Kf^2*Q*phi1P_max^2*2.0+CB6F*Kd*Kp2*eta*kq+CB6F*Kf*Kp2*eta*kq+CB6F*Kp2*Ku2*eta*kq+CB6F*Kd*Kf*kq*phi1P_max*2.0+CB6F*Kd*Kp2*kq*phi1P_max+CB6F*Kf*Kp2*kq*phi1P_max+CB6F*Kd*Ku2*kq*phi1P_max+CB6F*Kf*Ku2*kq*phi1P_max+CB6F*Kp2*Ku2*kq*phi1P_max+Abs*Kd*Kf*Q*phi1P_max^2*4.0+Abs*Kd*Kp2*Q*phi1P_max^2*2.0+Abs*Kf*Kp2*Q*phi1P_max^2*2.0+Abs*Kd*Ku2*Q*phi1P_max^2*2.0+Abs*Kf*Ku2*Q*phi1P_max^2*2.0+Abs*Kd*Kp2*Q*eta*phi1P_max+Abs*Kf*Kp2*Q*eta*phi1P_max+Abs*Kp2*Ku2*Q*eta*phi1P_max)/(Q*phi1P_max*(Kd^2*phi1P_max+Kf^2*phi1P_max+Kd*Kp2*eta+Kf*Kp2*eta+Kp2*Ku2*eta+Kd*Kf*phi1P_max*2.0+Kd*Kp2*phi1P_max+Kf*Kp2*phi1P_max+Kd*Ku2*phi1P_max+Kf*Ku2*phi1P_max)*2.0)

        a2 = solve_xcs(Abs, CB6F, Kd, Kf, Kp2, Ku2, Q, eta, kq, phi1P_max)
        a1 = Abs .- a2
    end

    # Calculate limiting rates for gas-exchange and electron transport

    # Expressions for potential Rubisco-limited rates (_c)
    #   N.B., see Eqns. 32-33
    JP700_j = @. (Q * Vqmax) / (Q + Vqmax / (a1 * phi1P_max))
    JP680_j = @. JP700_j / eta
    Vc_j = @. JP680_j / (4.0 * (1.0 + 2.0 * gammas / C))
    Vo_j = @. Vc_j * 2.0 * gammas / C
    Ag_j = @. Vc_j - Vo_j / 2.0

    # Expressions for potential Rubisco-limited rates (_c)
    Vc_c = @. C * Vcmax / (C + Kc * (1.0 + O / Ko))
    Vo_c = @. Vc_c * 2.0 * gammas / C
    Ag_c = @. Vc_c - Vo_c / 2.0
    JP680_c = @. Ag_c * 4.0 * (1.0 + 2.0 * gammas / C) / (1.0 - gammas / C)
    # This quantity here is a scalar, but needs to be repeated to be of the
    # same size as JP680_c.
    JP700_c = repeat([JP680_c * eta], length(JP700_j))

    # Define anonymous function for quadratic to smooth transitions
    #   N.B., this returns an array with the two roots of the quadratic

    tr(l1, l2, th) = [
        ((l1 + l2) + sqrt((l1 + l2)^2 - 4.0 * th * l1 * l2)) / (2.0 * th),
        ((l1 + l2) - sqrt((l1 + l2)^2 - 4.0 * th * l1 * l2)) / (2.0 * th)
    ]

    # Select minimum PS1 ETR
    JP700_a = minimum.(tr.(JP700_j, JP700_c, Ref(theta1)))
    JP680_a = minimum.(tr.(JP680_j, JP680_c, Ref(theta1)))

    # Select minimum Ag_a
    Ag_a  = (
        minimum.(tr.(Ag_j, Ag_c, Ref(theta1))) .* (C > gammas) +
        maximum.(tr.(Ag_j, Ag_c, Ref(theta1))) .* (C <= gammas)
    )
    An_a = Ag_a .- Rd

    if alpha_opt == "dynamic"
        # Derive a2/a1 at light saturation point and update
        # N.B., this represents dynamic optimization of a2/a1 under
        # limiting light, and then under saturating light holds a2/a1 at
        # the values attained at the light saturation point.

        #I = findmin(hcat(JP700_j, JP700_c), dims=2)[2]
        I = (p -> p.I[2]).(findmin(hcat(JP700_j, JP700_c), dims=2)[2])
        which = findall(diff(I, dims=1) .== 1)[end]

        a2_new = repeat([a2[which.I[1]]], length(a2))
        a2_old = copy(a2)

        a2_update = zeros(length(a2))
        a2_update[a2_old .> a2_new] = a2_old[a2_old .> a2_new]
        a2_update[a2_old .<= a2_new] = a2_new[a2_old .<= a2_new]

        #println(a2_old .> a2_new)
        #println(a2_old .<= a2_new)

        a2 = copy(a2_update)
        a1 = Abs .- a2
    end

    CB6F_a = @. JP700_j / kq         # Eqns. 21, 30a, 34
    phi1P_a = @. JP700_a / (Q * a1)  # Eqn. 20
    q1_a = @. phi1P_a / phi1P_max    # Eqn. 19a
    phi2P_a = @. JP680_a / (Q * a2)  # Eqn. 26
    q2_a = @. 1.0 - CB6F_a / CB6F    # Eqns. 28 and 34

    # N.B., rearrange Eqn. 25a to solve for Kn2_a

    Kn2_a = @. (
        sqrt(Kp2^2 * phi2P_a^2 - 2.0 * Kp2^2 * phi2P_a * q2_a +
             Kp2^2 * q2_a.^2 - 4.0 * Kp2 * Ku2 * phi2P_a^2 * q2_a .+
             2.0 * Kp2 * Ku2 * phi2P_a^2 + 2.0 * Kp2 * Ku2 * phi2P_a * q2_a .+
             Ku2.^2 * phi2P_a^2)
        - Kp2 * phi2P_a + Ku2 * phi2P_a + Kp2 * q2_a) / (2.0 * phi2P_a) - Kf - Ku2 - Kd

    # Photosystem II (Eqns. 23a-23e and 25a-25d)
    phi2p_a = @. (q2_a) * Kp2 / (Kp2 + Kn2_a + Kd + Kf + Ku2)
    phi2n_a = @. (q2_a) * Kn2_a / (Kp2 + Kn2_a + Kd + Kf + Ku2) +
        (1 - q2_a) * Kn2_a / (Kn2_a + Kd + Kf + Ku2)
    phi2d_a = @. (q2_a) * Kd / (Kp2 + Kn2_a + Kd + Kf + Ku2) +
        (1 - q2_a) * Kd / (Kn2_a + Kd + Kf + Ku2)
    phi2f_a = @. (q2_a) * Kf / (Kp2 + Kn2_a + Kd + Kf + Ku2) +
        (1 - q2_a) * Kf ./ (Kn2_a + Kd + Kf + Ku2)
    phi2u_a = @. (q2_a) * Ku2 / (Kp2 + Kn2_a + Kd + Kf + Ku2) +
        (1 - q2_a) * Ku2 / (Kn2_a + Kd + Kf + Ku2)

    phi2P_a = @. phi2p_a / (1.0 - phi2u_a)
    phi2N_a = @. phi2n_a / (1.0 - phi2u_a)
    phi2D_a = @. phi2d_a / (1.0 - phi2u_a)
    phi2F_a = @. phi2f_a / (1.0 - phi2u_a)

    # For Photosystem I (Eqns. 19a-19d)
    phi1P_a = @. q1_a * Kp1 / (Kp1 + Kd + Kf)
    phi1N_a = @. (1.0 - q1_a) * Kn1 / (Kn1 + Kd + Kf)
    phi1D_a = @. q1_a * Kd / (Kp1 + Kd + Kf) + (1.0 - q1_a) * Kd / (Kn1 + Kd + Kf)
    phi1F_a = @. q1_a * Kf / (Kp1 + Kd + Kf) + (1.0 - q1_a) * Kf / (Kn1 + Kd + Kf)

    # PAM measured fluorescence levels (Eqns. 38-42)
    #   N.B., hardcoding of a2(1) for dark-adapted value
    ### this means that we MUST supply a wavelength array with a low radiance
    ### as its first value - for the case of the dynamic model!!
    Fm_a = @. a2[1] * Kf / (Kd + Kf) * eps2 + a1[1] * Kf / (Kn1 + Kd + Kf) * eps1
    Fo_a = @. a2[1] * Kf / (Kp2 + Kd + Kf) * eps2 + a1[1] * Kf / (Kp1 + Kd + Kf) * eps1
    Fmp_a = @. a2 * Kf / (Kn2_a + Kd + Kf) * eps2 + a1 * Kf / (Kn1 + Kd + Kf) * eps1
    Fop_a = @. a2 * Kf / (Kp2 + Kn2_a + Kd + Kf) * eps2 + a1 * Kf / (Kp1 + Kd + Kf) * eps1
    Fs_a = @. a2 * phi2F_a * eps2 + a1 * phi1F_a * eps1

    # PAM indices
    PAM1_a = @. 1.0 - Fs_a / Fmp_a # PhiP
    PAM2_a = @. Fs_a * (1.0 / Fmp_a - 1.0 / Fm_a) # PhiN
    PAM3_a = @. Fs_a / Fm_a # PhiD + PhiF
    PAM4_a = @. Q * 0.85 / 2 * PAM1_a # ETR
    PAM5_a = @. (Fmp_a - Fs_a) / (Fmp_a - Fop_a) # qP
    PAM6_a = @. (Fmp_a - Fs_a) * Fop_a / ((Fmp_a - Fop_a) * Fs_a) # qL
    PAM7_a = @. PAM4_a / (1.0 - PAM5_a) # kPuddle
    PAM8_a = @. PAM4_a / (1.0 - PAM6_a) # kLake
    PAM9_a = @. Fm_a / Fmp_a - 1.0 # NPQ

    return PAM1_a, PAM2_a, PAM3_a, PAM4_a, PAM5_a, PAM6_a, PAM7_a, PAM8_a, PAM9_a, Fs_a, An_a, a1 .* phi1F_a, a2 .* phi2F_a

end

function model_fluorescence(
    Q::Number, # We want this to be a number, not a vector
    T,
    C,
    O,
    wavelength, # this should be in nm
    NDVI
)

    fAPAR_chl = calculate_fAPAR(NDVI)

    if fAPAR_chl <= 1e-3
        return 0.0
    end

    fAPAR_abs = 0.85
    fAPAR_int = fAPAR_chl / fAPAR_abs
    fESC = 1.0

    p1, p2 = photosynthetic_yields(
        Q, # Q in umol PPFD m^-2 s^-1, INCIDENT PAR
        T, # leaf temp in C
        C, # CO2 pp in ubar
        O, # O2 pp in mbar
        Abs=fAPAR_abs * fAPAR_chl, # Total leaf absorbance to PAR [mol / mol]
        alpha_opt="static"
    )[[12, 13]]

    if p1 isa Vector
        p1 = p1[1]
    end

    if p2 isa Vector
        p2 = p2[1]
    end

    # Find where we need to sample the spectra
    wl_idx = argmin(abs.(ps_spectra_wl .- wavelength))
    hcNA = 1.0 / 8.359347229111778 # kg m^3 / mol / s2

    # Q: incoming radiance in PPFD mol / m2 / s
    # outgoing radiance at wavelength [PPFD umol/m2/s/nm]
    # (this is radiating out in to the full space)
    # (phiI, phiII is in units of 1/nm, p1 and p2 is in [1])
    SIF_PPFD = fAPAR_chl * Q * (p1 * ps_spectra_phiI[wl_idx] + p2 * ps_spectra_phiII[wl_idx])

    # outgoing radiance at wavelength [PPFD mol/m2/s/nm]
    SIF_PPFD *= 1e-3
    # outgoing radiance at wavelength [W/m2]
    SIF_PPFD *= hcNA * 1e9
    # outgoing radiance at wavelength [W/m2/um]
    SIF_PPFD /= 1e0 * wavelength
    # outgoing radiance at wavelength [W/m2/um/sr]
    SIF_PPFD /= 4 * pi
    # this final quantity should be roughly peaking around 1-3 W/m2/sr/um

end

function create_SIF_scene(
    instrument::String,
    mode::String,
    loctime::GeolocationTime,
    vza::Real,
    uncertainty_function::Function,
    vnp_sd::VNPData
    )

    # solar azimuth (saa) at this point unused!
    this_sza, this_saa = calculate_solar_angles(loctime)

    # Obtain irradiance information (downwelling radiance at surface), already SZA-corrected
    # (unlike in L2 algorithms, where irradiance needs to be multiplied by
    #  mu0 to account for normal component)
    this_irradiance = calculate_BOA_irradiance(this_sza, 757.0)

    # Calculate PPFD for this scene
    this_PPFD = calculate_PPFD(this_sza)

    # black sky albedos from BRDFs?
    this_nir = calculate_reflectance(loctime, this_sza, "M7", vnp_sd)
    this_vis = calculate_reflectance(loctime, this_sza, "M5", vnp_sd)

    # These wavelengths are for VIIRS M5 and M7,
    # factor of 1000 takes us from W/m2/nm/sr to W/m2/um/sr
    this_nir_radiance = this_nir * calculate_BOA_irradiance(this_sza, 865.0) * 1000 / pi
    this_vis_radiance = this_vis * calculate_BOA_irradiance(this_sza, 672.0) * 1000 / pi

    this_ndvi = (this_nir - this_vis) / (this_nir + this_vis)

    this_nirv = this_ndvi * this_nir
    this_nirv_radiance = this_nirv * calculate_BOA_irradiance(this_sza, 865.0) * 1000 / pi

    # Reflectance at ~757 nm is roughly between VIIRS bands M5 and M7
    refl_M5 = calculate_reflectance(loctime, this_sza, "M5", vnp_sd)
    refl_M7 = calculate_reflectance(loctime, this_sza, "M7", vnp_sd)

    this_reflectance = 0.5 * (refl_M5 + refl_M7)

    # irradiance is in /nm, but we want /um
    this_TOA_radiance = this_irradiance * this_reflectance * 1000 / pi

    # This is in W/m2/sr/um
    this_sif = model_fluorescence(
        this_PPFD,
        25.0,
        200.0,
        209.0,
        757.0, # wavelength in nm
        this_ndvi
    )

    # The function needed to calculate the uncertainty must be supplied
    this_sif_ucert = uncertainty_function(this_TOA_radiance)

    # Scale SIF and uncertainty values up for 740nm
    # Note that the native calculation is still done for ~755nm,
    # however the resultws are scaled up to be comparable to e.g. TROPOMI

    this_sif *= 1.56
    this_sif_ucert *= 1.56

    # Noisified SIF
    if this_sif_ucert > 0.
        # Apply global seed variable to sampling
        # Random.seed!(jLRS_seed)
        this_measured_sif = rand(Normal(this_sif, this_sif_ucert))
    else
        this_measured_sif = this_sif
    end

    # At the moment no cloud or aerosol data,
    # could add ISCCP sampler in here
    this_od = 0.0

    this_scene = Scene(
        instrument,
        mode, # sampling mode comes from the instrument
        loctime, # location time comes from the instrument
        this_measured_sif,
        this_sif,
        this_sif_ucert,
        this_sza, # SZA comes from calculcations (via loctime)
        vza, # viewing zenith comes from the instrument,
        this_nir, # NIR reflectance
        this_vis, # VIS reflectance
        this_ndvi, # NDVI from VIIRS
        this_nirv, # NIRv calculated from NDVI * NIR
        this_nirv_radiance, # NIRv * L0(868nm)
        this_irradiance, # Irradiance at the surface and some ref. wl
        this_PPFD, # Integrated irradiance at PAR wavelengths 400nm to 700nm
        this_reflectance, # Reflectance comes from BRDF sampling and SZA
        this_od # optical depth may come from ISCCP one day
    )

    return this_scene

end


text:
cross(o::Vector{Float64}, a::Vector{Float64}, b::Vector{Float64}) = (a[1]-o[1])*(b[2]-o[2]) - (a[2]-o[2])*(b[1]-o[1])

function convex(points::Vector{Vector{Float64}})
    lower = Vector{Float64}[]
    upper = Vector{Float64}[]
    lt(a,b) = a[1] == b[1] ? a[2] < b[2] : a[1] < b[1]
    sort!(points, lt=lt)
    for i=1:length(points)
        while length(lower) >= 2 && cross(lower[end-2], lower[end-1], points[i]) <= 0
            pop!(lower)
        end
        push!(lower, points[i])
    end
    for i=length(points):-1:1
        while length(upper) >= 2 && cross(upper[end-2], upper[end-1], points[i]) <= 0
            pop!(upper)
        end
        push!(upper, points[i])
    end
    pop!(upper)
    pop!(lower)
    coords = [lower, upper]
    push!(coords, coords[1])
    GeoJSON.Feature(GeoJSON.Polygon(Vector{Vector{Float64}}[coords]))
end


text:
"""
    crossover!(model, de, group)

Performs crossover step for each particle pt in the chain

# Arguments

- `model`: model containing a likelihood function with data and priors
- `de`: differential evolution object
- `group`: a group of particles
"""
function crossover!(model, de, group)
    for pt in group
        crossover!(model, de, group, pt)
    end
    return nothing
end

"""
    crossover!(model, de, group, pt::Particle)

Performs crossover step for for a given particle pt

# Arguments

- `model`: model containing a likelihood function with data and priors
- `de`: differential evolution object
- `pt`: a target particle to be updated
"""
function crossover!(model, de, group, pt::Particle)
    if rand() ≤ de.θsnooker
        # generate the proposal
        proposal,pz = snooker_update!(de, pt, group)
        log_adj = adjust_loglike(pt, proposal, pz)
        # compute the weight of the proposal: prior loglikelihood + data loglikelihood
        de.evaluate_fitness!(de, model, proposal)
        # accept proposal according to Metropolis-Hastings rule
        de.update_particle!(de, pt, proposal, log_adj)
    else
        # generate the proposal
        proposal = de.generate_proposal(de, pt, group)
        # compute the weight of the proposal: prior loglikelihood + data loglikelihood
        de.evaluate_fitness!(de, model, proposal)
        # accept proposal according to Metropolis-Hastings rule
        de.update_particle!(de, pt, proposal)
    end
end

"""
    crossover!(model, de, group)

Performs crossover step for each particle pt in the chain

# Arguments

- `model`: model containing a likelihood function with data and priors
- `de`: differential evolution object
- `group`: a group of particles
- `block`: a vector of boolean values indicating which parameters to update
"""
function crossover!(model, de, group, block)
    for pt in group
        crossover!(model, de, group, pt, block)
    end
    return nothing
end

"""
    crossover!(model, de, group, pt::Particle, block)

Performs crossover step for a given particle pt.

# Arguments

- `model`: model containing a likelihood function with data and priors
- `de`: differential evolution object
- `group`: a group of particles
- `block`: a vector of boolean values indicating which parameters to update
"""
function crossover!(model, de, group, pt::Particle, block)
    if rand() ≤ de.θsnooker
        # generate the proposal
        proposal,pz = snooker_update!(de, pt, group)
        reset!(proposal, pt, block)
        log_adj = adjust_loglike(pt, proposal, pz)
        # compute the weight of the proposal: prior loglikelihood + data loglikelihood
        de.evaluate_fitness!(de, model, proposal)
        # accept proposal according to Metropolis-Hastings rule
        de.update_particle!(de, pt, proposal, log_adj)
    else
        # generate the proposal
        proposal = de.generate_proposal(de, pt, group)
        reset!(proposal, pt, block)
        # compute the weight of the proposal: prior loglikelihood + data loglikelihood
        de.evaluate_fitness!(de, model, proposal)
        # accept proposal according to Metropolis-Hastings rule
        de.update_particle!(de, pt, proposal)
    end
end

"""
    resample(de, group, n, replace)

Sample a random particle from previously accepted values for snooker update.

# Arguments

- `de`: differential evolution object
- `group`: a group of particles
- `n`: number of particles to sample
- `replace`: sample with replacement if true
"""
function resample(de, group, n, replace)
    P′ = Vector{eltype(group)}(undef,n)
    mx_idx = de.iter - 1
    idx = sample_indices(de.samples, mx_idx, n; replace)
    for i in 1:n
        P′[i] = Particle(;Θ=de.samples[idx[i][1],:,idx[i][2]])
    end
    return P′
end

sample_indices(x, ub, n; replace) = @views sample(CartesianIndices(x[1:ub,1,:]), n; replace)

"""
    sample(de::DE, group_diff, n, replace)

Sample a random particle.

# Arguments

- `de`: differential evolution object
- `group`: a group of particles
- `n`: number of particles to sample
- `replace`: sample with replacement if true
"""
function sample(de::DE, group_diff, n, replace)
    return sample(group_diff, n; replace)
end

"""
    random_gamma(de, Pt, group)

Generate proposal according to θ' = θt + γ1(θm − θn) + γ2(θb − θt) + b.
γ2=0 after burnin

# Arguments

- `de`: differential evolution object
- `Pt`: current particle
- `group`: a group of particles
"""
function random_gamma(de, Pt, group)
    # select the base particle θb
    Pb = select_base(group)
    # group without Pb
    group_diff = setdiff(group, [Pt])
    # sample particles for θm and θn
    Pm,Pn = de.sample(de, group_diff, 2, false)
    # sample gamma weights
    γ₁ = rand(Uniform(.5, 1))
    # set γ₂ = 0 after burnin
    γ₂ = de.iter > de.burnin ? 0.0 : rand(Uniform(.5, 1))
    # sample noise for each parameter
    b = Uniform(-de.ϵ, de.ϵ)
    # compute proposal value
    Θp = Pt + γ₁ * (Pm - Pn) + γ₂ * (Pb - Pt) + b
    # reset each parameter to previous value with probability (1-κ)
    recombination!(de, Pt, Θp)
    return Θp
end

"""
    fixed_gamma(de, Pt, group)

Generate proposal according to θ' = θt + γ(θm − θn) + b
where γ = 2.38.

# Arguments

- `de`: differential evolution object
- `Pt`: current particle
- `group`: a group of particles
"""
function fixed_gamma(de, Pt, group)
    group_diff = setdiff(group, [Pt])
    # sample particles for θm and θn
    Pm,Pn = de.sample(de, group_diff, 2, false)
    # sample gamma weights
    γ = 2.38
    # sample noise for each parameter
    b = Uniform(-de.ϵ, de.ϵ)
    # compute proposal value
    Θp = Pt + γ * (Pm - Pn) + b
    # reset each parameter to previous value with probability (1-κ)
    recombination!(de, Pt, Θp)
    return Θp
end

"""
    variable_gamma(de, Pt, group)

Generate proposal according to θ' = θt + γ(θm − θn) + b
where γ = 2.38/√(2d) where d is the number of parameters

# Arguments

- `de`: differential evolution object
- `Pt`: current particle
- `group`: a group of particles
"""
function variable_gamma(de, Pt, group)
    Np = sum(length.(Pt.Θ))
    group_diff = setdiff(group, [Pt])
    # sample particles for θm and θn
    Pm,Pn = de.sample(de, group_diff, 2, false)
    γ = 2.38 / sqrt(2 * Np)
    # sample noise for each parameter
    b = Uniform(-de.ϵ, de.ϵ)
    # compute proposal value
    Θp = Pt + γ * (Pm - Pn) + b
    # reset each parameter to previous value with probability (1-κ)
    recombination!(de, Pt, Θp)
    return Θp
end

"""
    snooker_update!(de, Pt, group)

Performs snooker update during crossover

# Arguments

- `de`: differential evolution object
- `Pt`: current particle
- `group`: a group of particles
"""
function snooker_update!(de, Pt, group)
    # sample 3 particles without replacement
    Pz,Pm,Pn = de.sample(de, group, 3, false)
    # difference vector for old value and random particle z
    Pd = Pt - Pz
    # project Pm on to Pd
    Pr1 = project(Pm, Pd)
    # project Pn on to Pd
    Pr2 = project(Pn, Pd)
    # random γ value
    γ = rand(Uniform(1.2, 2.2))
    # sample noise for each parameter
    b = Uniform(-de.ϵ, de.ϵ)
    # compute proposal value
    Θp = Pt + γ * (Pr1 - Pr2) + b
    # reset each parameter to previous value with probability (1-κ)
    recombination!(de, Pt, Θp)
    return Θp,Pz
end

"""
    adjust_loglike(Pt, proposal, Pz)

The adjusted log likelihood component for a snooker update.

- `Pt`: the target particle from iteration n - 1
- `proposal`: the proposal particle on iteration n
- `Pz`: the particle formed by the projection of particles m and n
"""
function adjust_loglike(Pt, proposal, Pz)
    Np = sum(length.(Pt.Θ))
    adj1 = norm(proposal - Pz)^(Np - 1)
    adj2 = norm(Pt - Pz)^(Np - 1)
    return log(adj1 / adj2)
end

"""
    select_base(group)

Selects base particle θb with probability proportional to weight.

- `group`: a group of particles
"""
function select_base(group)
    w = map(x -> x.weight, group)
    θ = exp.(w) / sum(exp.(w))
    p = sample(group, Weights(θ))
    return p
end

"""
    recombination!(de, pt::Particle, pp::Particle)

Resets parameters of proposal to previous value with probability
(1-κ).

- `de`: differential evolution object
- `pt`: current partical
- `pp`: proposal particle
"""
function recombination!(de, pt::Particle, pp::Particle)
    de.κ == 1.0 ? (return nothing) : nothing
    N = length(pt.Θ)
    for i in 1:N
        if isa(pp.Θ[i], Array)
            recombination!(de, pt.Θ[i], pp.Θ[i])
        else
            pp.Θ[i] = rand() <= (1 - de.κ) ?  pt.Θ[i] : pp.Θ[i]
        end
    end
    return nothing
end

# Handles elements within an array
function recombination!(de, Θt, Θp)
    N = length(Θt)
    for i in 1:N
        Θp[i] = rand() <= (1 - de.κ) ? Θt[i] : Θp[i]
    end
    return nothing
end

"""
    reset!(p1::Particle, p2::Particle, idx)

During block updates, all parameters are updated. This function resets
parameters in the proposal `p1` with values from the previous particle.

# Arguments

`p1`: proposal particle which will have values reset
`p2`: previous particle
`idx`: boolean vector indicating which values are updated in block. False
values are reset.
"""
reset!(p1::Particle, p2::Particle, idx) = reset!(p1.Θ, p2.Θ, idx)

function reset!(Θ1, Θ2, idx)
    for i in 1:length(Θ1)
        reset!(Θ1, Θ2, idx[i], i)
    end
    return nothing
end

function reset!(Θ1, Θ2, idx::Bool, i)
    !idx ? Θ1[i] = Θ2[i] : nothing
    return nothing
end

function reset!(Θ1, Θ2, idx::Array{Bool,N}, i) where {N}
    return reset!(Θ1[i], Θ2[i], idx)
end


text:
using Multicombinations
using Base.Test

# multichoose

@test multichoose(3,2) == 6
@test multichoose(3,0) == 1
@test multichoose(0,3) == 0
@test multichoose(0,0) == 1
@test multichoose(-1,-1) == 0

# multicombinations

@test collect(multicombinations("abc",3)) == ["aaa","aab","aac","abb","abc","acc","bbb","bbc","bcc","ccc"]
@test collect(multicombinations("åäö",3)) == ["ååå","ååä","ååö","åää","åäö","åöö","äää","ääö","äöö","ööö"]
@test collect(multicombinations("abc",0)) == collect(combinations("abc",0))
@test collect(multicombinations("",1)) == collect(combinations("",1))
@test collect(multicombinations("",0)) == collect(combinations("",0))
@test collect(multicombinations("",-1)) == collect(combinations("",-1))

# integersums

@test collect(integersums(3,3)) == {[3,0,0],[2,1,0],[2,0,1],[1,2,0],[1,1,1],[1,0,2],[0,3,0],[0,2,1],[0,1,2],[0,0,3]}
@test collect(integersums(3,0)) == {[0,0,0]}
@test collect(integersums(0,3)) == {}
@test collect(integersums(0,0)) == {[]}
@test collect(integersums(-1,-1)) == {}

# is2mc

@test map(is2mc,integersums(3,3)) == collect(multicombinations(1:3,3))
@test map(is2mc,integersums(3,0)) == collect(multicombinations(1:3,0))
@test map(is2mc,integersums(0,3)) == collect(multicombinations([],3))
@test map(is2mc,integersums(0,0)) == collect(multicombinations([],0))

# length

@test length(collect(multicombinations("abc",3))) == length(multicombinations("abc",3))
@test length(collect(multicombinations("åäö",3))) == length(multicombinations("åäö",3))
@test length(collect(integersums(3,3))) == length(integersums(3,3))


text:
using Revise
using caramel
using Test
using JuMP

function inventory_feasibility(V_::Vector)
    starting = sum(V(i).START for i in V_)
    demands = sum(sum(V(i).d) for i in V_)
    minim = sum(V(i).MIN for i in V_)

    return starting - minim >= demands
end

function EM_check(E_::Vector, M_::Vector)
    #unique md in edge list == all m in M
    from_E = unique(md.(E.(E_)))
    from_M = M_

    return sort(from_E) == sort(from_M)
end

function EV_check(E_::Vector, V_::Vector)
    #unique vtx in ori and dst of edge_list == all in V
    from_E = unique(union(ori.(E.(E_)), dst.(E.(E_))))
    from_V = V_

    return sort(from_E) == sort(from_V)
end

@testset "caramel.jl" begin
    include("raw_test.jl")
end


text:
module PLY

export plyopen, plyclose, plygetvertexproperties, plygetvertices, plygetvertexnormals, PLY_CHAR, PLY_SHORT, PLY_INT, PLY_UCHAR, PLY_USHORT, PLY_UINT, PLY_FLOAT, PLY_DOUBLE

const PLY_CHAR = 1
const PLY_SHORT = 2
const PLY_INT = 3
const PLY_UCHAR = 4
const PLY_USHORT = 5
const PLY_UINT = 6
const PLY_FLOAT = 7
const PLY_DOUBLE = 8

function plyopen(fn::ASCIIString)
	nelems = zeros(Int32,1)
	elem_names = Array{ASCIIString,1}[]
	push!(elem_names, ASCIIString[])
	push!(elem_names, ASCIIString[])
	push!(elem_names[1], "testtesttest")
	push!(elem_names[2], "testtesttest")

	file_type = zeros(Int32,1)
	version = [0.0f0]
	ply = ccall( (:ply_open_for_reading, "libplyfile"),
	            Ptr{Void},
							(Ptr{Uint8}, Ptr{Int32}, Ptr{Ptr{Ptr{Uint8}}}, Ptr{Int32}, Ptr{Float32}),
							fn,nelems,elem_names,file_type,version)
	return ply,nelems,elem_names,file_type,version
end

function plygetvertexproperties(ply::Ptr{Void})
	nelems = zeros(Int32,1)
	nprops = zeros(Int32,1)
	datatype = zeros(Int32,1)

	props = ccall( (:ply_get_vertex_properties, "libplyfile"),
	               Ptr{Void},
								 (Ptr{Void}, Ptr{Int32}, Ptr{Int32}, Ptr{Int32}),
								 ply, nelems, nprops, datatype )

	return nelems[1],nprops[1],datatype[1]
end

function plygetvertices(ply::Ptr{Void}, nelems, datatype, data::Array{Float32,2})
	nelems = convert(Int32, nelems)
	datatype = convert(Int32, datatype)

	ccall( (:ply_get_vertices, "libplyfile"),
	       Void,
				 (Ptr{Void}, Int32, Int32, Ptr{Float32}),
				 ply, nelems, datatype, pointer(data) )
end

function plygetvertexnormals(ply::Ptr{Void}, nelems, datatype, data::Array{Float32,2})
	nelems = convert(Int32, nelems)
	datatype = convert(Int32, datatype)

	ccall( (:ply_get_vertex_normals, "libplyfile"),
	       Void,
				 (Ptr{Void}, Int32, Int32, Ptr{Float32}),
				 ply, nelems, datatype, pointer(data) )
end

function plyclose(handle::Ptr{Void})
	ccall( (:ply_close, "libplyfile"),
	        Void,
					(Ptr{Void},),
					handle)
end

end # module PLY


text:
# # Practical error bounds for the forces
#
# This is a simple example showing how to compute error estimates for the forces
# on a ``{\rm TiO}_2`` molecule, from which we can either compute asymptotically
# valid error bounds or increase the precision on the computation of the forces.
#
# The strategy we follow is described with more details in [^CDKL2021] and we
# will use in comments the density matrices framework. We will also needs
# operators and functions from
# [`src/scf/newton.jl`](https://dftk.org/blob/master/src/scf/newton.jl).
#
# [^CDKL2021]:
#     E. Cancès, G. Dusson, G. Kemlin, and A. Levitt
#     *Practical error bounds for properties in plane-wave electronic structure
#     calculations* Preprint, 2021. [arXiv](https://arxiv.org/abs/2111.01470)
using DFTK
using LinearAlgebra
using ForwardDiff
using LinearMaps
using IterativeSolvers

# ## Setup
# We setup manually the ``{\rm TiO}_2`` configuration from
# [Materials Project](https://materialsproject.org/materials/mp-2657/).
Ti = ElementPsp(:Ti, psp=load_psp("hgh/lda/ti-q4.hgh"))
O = ElementPsp(:O, psp=load_psp("hgh/lda/o-q6.hgh"))
atoms = [Ti => [[0.5, 0.5, 0.5],
                [0.0, 0.0, 0.0]],
         O  => [[0.19542, 0.80458, 0.5],
                [0.80458, 0.19542, 0.5],
                [0.30458, 0.30458, 0.0],
                [0.69542, 0.69542, 0.0]]]
lattice = [[8.79341  0.0      0.0];
           [0.0      8.79341  0.0];
           [0.0      0.0      5.61098]];

# We apply a small displacement to one of the ``\rm Ti`` atoms to get nonzero
# forces.
el, pos = atoms[1]
atoms[1] = Pair(el, pos .+ [[-0.22, 0.28, 0.35] / 10, [0, 0, 0]]);

# We build a model with one k-point only, not too high `Ecut_ref` and small
# tolerance to limit computational time. These parameters can be increased for
# more precise results.
model = model_LDA(lattice, atoms)
kgrid = [1, 1, 1]
Ecut_ref = 35
basis_ref = PlaneWaveBasis(model; Ecut=Ecut_ref, kgrid)
tol = 1e-8;

# ## Computations

# We compute the reference solution ``P_*`` from which we will compute the
# references forces.
scfres_ref = self_consistent_field(basis_ref, tol=tol, callback=info->nothing)
ψ_ref, _ = DFTK.select_occupied_orbitals(basis_ref, scfres_ref.ψ,
                                         scfres_ref.occupation);

# We compute a variational approximation of the reference solution with
# smaller `Ecut`. `ψr`, `ρr` and `Er` are the quantities computed with `Ecut`
# and then extended to the reference grid.
#
# !!! note "Choice of convergence parameters"
#     Be careful to choose `Ecut` not too close to `Ecut_ref`.
#     Note also that the current choice `Ecut_ref = 35` is such that the
#     reference solution is not converged and `Ecut = 15` is such that the
#     asymptotic regime (crucial to validate the approach) is barely established.
Ecut = 15
basis = PlaneWaveBasis(model; Ecut=Ecut, kgrid)
scfres = self_consistent_field(basis, tol=tol, callback=info->nothing)
ψr = DFTK.transfer_blochwave(scfres.ψ, basis, basis_ref)
ρr = compute_density(basis_ref, ψr, scfres.occupation)
Er, hamr = energy_hamiltonian(basis_ref, ψr, scfres.occupation; ρ=ρr);

# We then compute several quantities that we need to evaluate the error bounds.

# - Compute the residual ``R(P)``, and remove the virtual orbitals, as required
#   in [`src/scf/newton.jl`](https://github.com/JuliaMolSim/DFTK.jl/blob/fedc720dab2d194b30d468501acd0f04bd4dd3d6/src/scf/newton.jl#L121).
res = DFTK.compute_projected_gradient(basis_ref, ψr, scfres.occupation)
res, occ = DFTK.select_occupied_orbitals(basis_ref, res, scfres.occupation)
ψr, _ = DFTK.select_occupied_orbitals(basis_ref, ψr, scfres.occupation);

# - Compute the error ``P-P_*`` on the associated orbitals ``ϕ-ψ`` after aligning
#   them: this is done by solving ``\min |ϕ - ψU|`` for ``U`` unitary matrix of
#   size ``N\times N`` (``N`` being the number of electrons) whose solution is
#   ``U = S(S^*S)^{-1/2}`` where ``S`` is the overlap matrix ``ψ^*ϕ``.
function compute_error(basis, ϕ, ψ)
    map(zip(ϕ, ψ)) do (ϕk, ψk)
        S = ψk'ϕk
        U = S*(S'S)^(-1/2)
        ϕk - ψk*U
    end
end
err = compute_error(basis_ref, ψr, ψ_ref);

# - Compute ``{\bm M}^{-1}R(P)`` with ``{\bm M}^{-1}`` defined in [^CDKL2021]:
P = [PreconditionerTPA(basis_ref, kpt) for kpt in basis_ref.kpoints]
map(zip(P, ψr)) do (Pk, ψk)
    DFTK.precondprep!(Pk, ψk)
end
function apply_M(φk, Pk, δφnk, n)
    DFTK.proj_tangent_kpt!(δφnk, φk)
    δφnk = sqrt.(Pk.mean_kin[n] .+ Pk.kin) .* δφnk
    DFTK.proj_tangent_kpt!(δφnk, φk)
    δφnk = sqrt.(Pk.mean_kin[n] .+ Pk.kin) .* δφnk
    DFTK.proj_tangent_kpt!(δφnk, φk)
end
function apply_inv_M(φk, Pk, δφnk, n)
    DFTK.proj_tangent_kpt!(δφnk, φk)
    op(x) = apply_M(φk, Pk, x, n)
    function f_ldiv!(x, y)
        x .= DFTK.proj_tangent_kpt(y, φk)
        x ./= (Pk.mean_kin[n] .+ Pk.kin)
        DFTK.proj_tangent_kpt!(x, φk)
    end
    J = LinearMap{eltype(φk)}(op, size(δφnk, 1))
    δφnk = cg(J, δφnk, Pl=DFTK.FunctionPreconditioner(f_ldiv!),
              verbose=false, reltol=0, abstol=1e-15)
    DFTK.proj_tangent_kpt!(δφnk, φk)
end
function apply_metric(φ, P, δφ, A::Function)
    map(enumerate(δφ)) do (ik, δφk)
        Aδφk = similar(δφk)
        φk = φ[ik]
        for n = 1:size(δφk,2)
            Aδφk[:,n] = A(φk, P[ik], δφk[:,n], n)
        end
        Aδφk
    end
end
Mres = apply_metric(ψr, P, res, apply_inv_M);

# We can now compute the modified residual ``R_{\rm Schur}(P)`` using a Schur
# complement to approximate the error on low-frequencies[^CDKL2021]:
#
# ```math
# \begin{bmatrix}
# (\bm \Omega + \bm K)_{11} & (\bm \Omega + \bm K)_{12} \\
# 0 & {\bm M}_{22}
# \end{bmatrix}
# \begin{bmatrix}
# P_{1} - P_{*1} \\ P_{2}-P_{*2}
# \end{bmatrix} =
# \begin{bmatrix}
# R_{1} \\ R_{2}
# \end{bmatrix}.
# ```

# - Compute the projection of the residual onto the high and low frequencies:
resLF = DFTK.transfer_blochwave(res, basis_ref, basis)
resHF = res - DFTK.transfer_blochwave(resLF, basis, basis_ref);

# - Compute ``{\bm M}^{-1}_{22}R_2(P)``:
e2 = apply_metric(ψr, P, resHF, apply_inv_M);

# - Compute the right hand side of the Schur system:
## Rayleigh coefficients needed for `apply_Ω`
Λ = map(enumerate(ψr)) do (ik, ψk)
    Hk = hamr.blocks[ik]
    Hψk = Hk * ψk
    ψk'Hψk
end
ΩpKe2 = DFTK.apply_Ω(e2, ψr, hamr, Λ) .+ DFTK.apply_K(basis_ref, e2, ψr, ρr, occ)
ΩpKe2 = DFTK.transfer_blochwave(ΩpKe2, basis_ref, basis)
rhs = resLF - ΩpKe2;

# - Solve the Schur system to compute ``R_{\rm Schur}(P)``: this is the most
#   costly step, but inverting ``\bm{\Omega} + \bm{K}`` on the small space has
#   the same cost than the full SCF cycle on the small grid.
ψ, _ = DFTK.select_occupied_orbitals(basis, scfres.ψ, scfres.occupation)
e1 = DFTK.solve_ΩplusK(basis, ψ, rhs, occ; tol_cg=tol)
e1 = DFTK.transfer_blochwave(e1, basis, basis_ref)
res_schur = e1 + Mres;

# ## Error estimates

# We start with different estimations of the forces on the first ``\rm Ti`` atom,
# in direction 1.
# - Forces computed with the reference solution:
f_ref = compute_forces(scfres_ref)
## [1][1][1] means that we look at the forces on the first atom type (Ti), the
## first atom of this type, in the first direction.
println("F(P_*) = $(f_ref[1][1][1])")

# - Forces computed with the variational approximation:
f = compute_forces(scfres)
println("F(P) = $(f[1][1][1])")

# We then try to improve ``F(P)`` using the first order linearization:
#
# ```math
# F(P) = F(P_*) + {\rm d}F(P)\cdot(P-P_*).
# ```

# To this end, we use the `ForwardDiff.jl` package to compute ``{\rm d}F(P)``
# using automatic differentiation.
function df(basis, occupation, ψ, δψ, ρ)
    δρ = DFTK.compute_δρ(basis, ψ, δψ, occupation)
    ForwardDiff.derivative(ε -> compute_forces(basis, ψ.+ε.*δψ, occupation; ρ=ρ+ε.*δρ), 0)
end;

# - Computation of the forces by a linearization argument if we have access to
#   the actual error ``P-P_*``:
df_err = df(basis_ref, occ, ψr, DFTK.proj_tangent(err, ψr), ρr)
println("F(P) - df(P)⋅(P-P_*) = $(f[1][1][1]-df_err[1][1][1])")

# - Computation of the forces by a linearization argument when replacing the
#   error ``P-P_*`` by the modified residual:
df_schur = df(basis_ref, occ, ψr, res_schur, ρr)
println("F(P) - df(P)⋅Rschur(P) = $(f[1][1][1]-df_schur[1][1][1])")

# Then, we estimate the error on the forces made by the different computations
# above:

# - Relative error on the forces with no post-processing:
println("|F(P) - F(P_*)| / |F(P_*)| = $(norm(f-f_ref)/norm(f_ref))")

# - Relative error made by the linearization ``F(P) - {\rm d}F(P)⋅(P-P_*)`` if
#   we had access to the actual error ``P-P_*``, which is not the case in
#   practice (we are aiming at reaching this precision):
println("|F(P) - dF(P)⋅(P-P_*) - F(P_*)| / |F(P_*)| = $(norm(f-df_err-f_ref)/norm(f_ref))")

# - Relative error made by replacing ``P-P_*`` by the modified residual
#   ``R_{\rm Schur}(P)`` (computable in practice) in the linearization (note
#   how closer we are to the previous one):
println("|F(P) - dF(P)⋅Rschur(P) - F(P_*)| / |F(P_*)| = $(norm(f-df_schur-f_ref)/norm(f_ref))")