text:
#=##############################################################################
# DESCRIPTION
    Lifting paneled body types definition.
# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Sep 2018
  * License   : MIT License
=###############################################################################

################################################################################
# ABSTRACT LIFTING BODY BODY TYPE
################################################################################
"""
  Implementations of AbstractBody are expected to have the following fields:
  * `U::Array{Int64,1}`                 : Indices of all panels along the upper
                                          side of the trailing edge.
  * `L::Array{Int64,1}`                 : Indices of all panels along the lower
                                          side of the trailing edge.
  * `ncellsTE::Int64`                   : Number of cells along trailing edge
  * `nnodesTE::Int64`                   : Number of nodes along trailing edge

  NOTE: U and L are assumed to have the same number of points.

  in addition to the same properties and function expected from an AbstractBody
  implementation with the following functions

  ```julia

  # Solve for distributions of potential flow
  function solve(self::BodyType, Vinfs::Array{Array{T,1},1},
                              D::Array{Array{T,1},1}, args...) where {T<:Real}
    .
    .
    .
  end

  # Outputs a vtk file with the wake
  function _savewake(self::BodyType, filename::String;
                                  len::Rtype=1.0, upper::Bool=true, optargs...)
    .
    .
    .
  end
  ```
"""
abstract type AbstractLiftingBody <: AbstractBody end

# Includes all implementations of AbstractLiftingBody
for header_name in ["liftingbody"]
  include("FLOWPanel_"*header_name*".jl")
end

# Declares implementations of AbstractLiftingBody
const LBodyTypes = Union{RigidWakeBody}

##### COMMON FUNCTIONS  ########################################################
"""
  `get_TE(self::LBodyTypes, i::Int64; upper::Bool=true)`

Returns the i-th trailing edge point.
"""
function get_TE(self::LBodyTypes, i::Int64; upper::Bool=true)
  if i>self.nnodesTE
    error("Invalid index $i; maximum is $(self.nnodesTE).")
  end

  if upper
    if i==self.nnodesTE
      return gt.get_node(self.grid, gt.get_cell(self.grid, self.U[i-1])[end])
    else
      return gt.get_node(self.grid, gt.get_cell(self.grid, self.U[i])[end-1])
    end
  else
    if i==self.nnodesTE
      return gt.get_node(self.grid, gt.get_cell(self.grid, self.L[i-1])[end-1])
    else
      return gt.get_node(self.grid, gt.get_cell(self.grid, self.L[i])[end])
    end
  end
end

"""
  `generate_loft_liftbody(bodytype::Type{LBodyTypes}, args...; optargs...)`
Generates a lofted lifting body of type `bodytype`. See documentation of
`GeometricTools.generate_loft` for a description of the arguments of this
function.
"""
function generate_loft_liftbody(bodytype::Type{LBodyTypes}, args...;
                                                  dimsplit::Int64=2, optargs...)
  # Lofts the surface geometry
  grid = gt.generate_loft(args...; optargs...)

  # Splits the quadrialateral panels into triangles
  # dimsplit = 2              # Dimension along which to split
  triang_grid = gt.GridTriangleSurface(grid, dimsplit)

  ndivs = gt.get_ndivscells(triang_grid)              # Cells in each dimension
  U = [ Base._sub2ind(ndivs, ndivs[1]-1, i) for i in 1:ndivs[2] ] # Upper LE cells
  L = [ Base._sub2ind(ndivs, 2, i) for i in 1:ndivs[2] ]          # Lower LE cells

  return bodytype(triang_grid, U, L)
end

"""
  `generate_revolution_liftbody(bodytype::Type{LBodyTypes}, args...; optargs...)`
Generates a lifting body type `bodytype` of a body of revolution. See
documentation of `GeometricTools.surface_revolution` for a description of the
arguments of this function.
"""
function generate_revolution_liftbody(bodytype::Type{LBodyTypes}, args...;
                                                  dimsplit::Int64=2,
                                                  loop_dim::Int64=2, optargs...)
  # Revolves the geometry
  grid = gt.surface_revolution(args...; loop_dim=loop_dim, optargs...)

  # Splits the quadrialateral panels into triangles
  # dimsplit = 2              # Dimension along which to split
  triang_grid = gt.GridTriangleSurface(grid, dimsplit)

  ndivs = gt.get_ndivscells(triang_grid)                 # Cells in each dimension
  U = [ Base._sub2ind(ndivs, ndivs[1]-1, i) for i in 1:ndivs[2] ] # Upper LE cells
  L = [ Base._sub2ind(ndivs, 2, i) for i in 1:ndivs[2] ]          # Lower LE cells

  return bodytype(triang_grid, U, L)
end

##### COMMON INTERNAL FUNCTIONS  ###############################################
"""
Returns the geometric matrix corresponding to only the vortex rings along a
lifting body impossing the Kutta condition at the trailing edge.
NOTE: This matrix doesn't include the wake.
"""
function _calc_G_lifting_vortexring(grid::gt.GridTriangleSurface,
                                    U::Array{Int64,1}, L::Array{Int64,1})
  return PanelSolver.G_lifting_vortexring(
                  grid.orggrid.nodes,                                  # Nodes
                  [gt.get_cell(grid, i) for i in 1:grid.ncells],       # Panels
                  U, L,                                                # TE
                  [_get_controlpoint(grid, i) for i in 1:grid.ncells], # CPs
                  [gt.get_normal(grid, i) for i in 1:grid.ncells];     # Normals
                )
end

"""
Checks correction definition of trailing edge
"""
function _checkTE(U::Array{Int64,1}, L::Array{Int64,1})
  # Same number of cells
  if size(U)!=size(L)
    return false
  end

  return true
end
##### END OF ABSTRACT LIFTING BODY #############################################


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

using Markdown
using InteractiveUtils

# ╔═╡ c3530026-5a87-4447-88f5-fc744cb48ff2
begin
	using PlutoUI
end

# ╔═╡ 2c895e14-e71a-435d-86d4-c1b101877b3c
md"""
### The End of the Beginning: Whole-Cell Metabolic Engineering

In the first lecture, we introduced metabolic engineering. Today, we'll close the loop at put together all we learned. Metabolic engineering is the practice of optimizing genetic and regulatory processes within cells to increase the cell's production of a desired small molecule or protein product of interest.

Metabolic engineers manipulate the biochemical networks cells use to convert raw materials into molecules necessary for the cell's survival. Metabolic engineering specifically seeks to:

1. Mathematically model biochemical networks, calculate the yield (product divided by substrate) of valuable products and identify parts of the network that constrain the production of these products of interest.
1. Use genetic engineering techniques to modify the biochemical network to relieve constraints limiting production. Metabolic engineers can then model the modified network to calculate the new product yield and identify new constraints (back to 1).

In this (the final lecture of Part I of the course), we'll:

* Incorporate the costs of transcription and translation into the flux balance analysis problem
* Introduce a workaround (often used in practice) to account for the metabolic cost of gene expression
* Close the loop: using a metabolic model and flux balance to optimize glycan production in _E.coli_.

"""

# ╔═╡ ce5c8e9c-7652-4c3e-98a7-98e3567dc695
md"""
### Integrating the Cost and Logic of Gene Expression With Flux Balance Analysis

The [Allen and Palsson study](https://pubmed.ncbi.nlm.nih.gov/12453446/) gave us a roadmap of integrating the cost of gene expression with the metabolic operation of the cell. Further, [Vilkhovoy et al.](https://pubmed.ncbi.nlm.nih.gov/29944340/) showed how we could do this in a cell-free system. However, how do you account for gene expression in metabolic models at the whole-genome scale?

* [Thiele I, Jamshidi N, Fleming RM, Palsson BØ. Genome-scale reconstruction of Escherichia coli's transcriptional and translational machinery: a knowledge base, its mathematical formulation, and its functional characterization. PLoS Comput Biol. 2009;5(3):e1000312. doi:10.1371/journal.pcbi.1000312](https://www.ncbi.nlm.nih.gov/labs/pmc/articles/PMC2648898/?report=classic)
* [Lerman JA, Hyduke DR, Latif H, Portnoy VA, Lewis NE, Orth JD, Schrimpe-Rutledge AC, Smith RD, Adkins JN, Zengler K, Palsson BO. In silico method for modeling metabolism and gene product expression at genome scale. Nat Commun. 2012 Jul 3;3:929. doi: 10.1038/ncomms1928. PMID: 22760628; PMCID: PMC3827721.](https://pubmed.ncbi.nlm.nih.gov/22760628/)
"""

# ╔═╡ bfd134fe-b311-4768-a53b-2c5c4c1b121a
md"""
##### Aside: A question of timescales
One open question is whether we need to solve the gene expression and flux estimation problem simultaneously. The time scale of enzyme activity (and hence metabolic flux) is much faster than gene expression. Why?

* The $k_{cat}$ for [phosphofructokinase-2 (pfk) in _E.coli_ is 9240 s$^{-1}$](https://bionumbers.hms.harvard.edu/bionumber.aspx?id=104955&ver=7&trm=kcat+pfk+in+E.+coli&org=). Thus, pfk has a characteristic time scale of $\tau\sim{k_{cat}}^{-1}$ or about 1$\times$10$^{-4}$ s.
* The $k_{cat}$ for RNAP polymerase is $k_{cat}$=$e_{x}L^{-1}$. If we assume [$e_{x}\sim{35}$ nt/s](https://bionumbers.hms.harvard.edu/bionumber.aspx?id=111871&ver=1&trm=elongation+rate+RNAP+E.coli+&org=) and average read length [$L\sim{924}$ nt](https://bionumbers.hms.harvard.edu/bionumber.aspx?id=111922&ver=3&trm=avarge+gene+length+in+E.+coli&org=), this gives $k_{cat}\sim$ 0.038 s$^{-1}$ or a timescale of $\tau\sim$26.4 s.

Thus, the time scale of metabolic flux (e.g., the timescale of enzyme activity) is approximately 5-orders of magnitude faster than gene expression. Therefore, gene expression seems constant  from the perspective of metabolic flux; hence, we treat these scales separately in traditional models.
"""

# ╔═╡ 8041cb6c-b80f-4e81-8d1d-5a2e9942d809
md"""
##### Aside: How do you account for macromolecular synthesis without the E/ME formulation?
The E/ME formulation of Palsson and coworkers is detailed and requires more time investment than perhaps you're willing to give in a metabolic engineering application. How do you account for the energy and resource cost of producing enzymes (and all the other machinery) required to make a cell?

The traditional way metabolic engineers address this question is to treat cell mass formation as just another reaction, e.g., the growth reaction $\mu$:

$$\left\{precursors\right\}~\rightarrow~cell$$

where each cellmass precursor has some stoichiometric coefficient $\sigma_{i,\mu}$ which describes how much of precursor $i$ is consumed (or produced) during cellmass formation. This formulation leads to species material balances of the form (assuming cellmass specific units, in a constant volume batch culture):

$$\frac{dx_{i}}{dt} = \sum_{j=1}^{\mathcal{R}}\sigma_{ij}\hat{r}_{j} - \left(x_{i}-\sigma_{i,\mu}\right)\mu\qquad{i=1,2,\dots,\mathcal{M}}$$

At steady state, for example in a flux balance analysis problem, these balances form the species constraints:

$$\sum_{j=1}^{\mathcal{R}}\sigma_{ij}\hat{r}_{j} - \left(x_{i}-\sigma_{i,\mu}\right)\mu = 0\qquad{i=1,2,\dots,\mathcal{M}}$$

Palsson and coworkers gave an example cell mass reaction for _E.coli_ in:

* [Orth JD, Fleming RM, Palsson BØ. Reconstruction and Use of Microbial Metabolic Networks: the Core Escherichia coli Metabolic Model as an Educational Guide. EcoSal Plus. 2010 Sep;4(1). doi: 10.1128/ecosalplus.10.2.1. PMID: 26443778.](https://pubmed.ncbi.nlm.nih.gov/26443778/)

A recent paper on sampling uncertain biomass stoichiometric coefficients:

* [Dinh et al. (2022) Quantifying the propagation of parametric uncertainty on flux balance analysis. Metabolic Engineering, 69:26-39; https://doi.org/10.1016/j.ymben.2021.10.012](https://www.sciencedirect.com/science/article/pii/S1096717621001634)

"""

# ╔═╡ 06c3eafa-10b3-4a47-9a06-39f6abd29db3
md"""
### Using flux balance analysis to improve metabolic performance
Ultimately a metabolic engineer uses a model to improve the performance of a production host or cell-free system. There is a huge variety of design approaches to improve metabolic function; let's look at a recent one from Cornell:

* [Wayman JA, Glasscock C, Mansell TJ, DeLisa MP, Varner JD. Improving designer glycan production in Escherichia coli through model-guided metabolic engineering. Metab Eng Commun. 2019;9:e00088. Published 2019 Mar 29. doi:10.1016/j.mec.2019.e00088](https://www.ncbi.nlm.nih.gov/labs/pmc/articles/PMC6454127/)

"""

# ╔═╡ b15b540b-8f13-4f21-b2da-510d34143b33
md"""
### Summary and conclusions

Today, we:

* Incorporated the costs of transcription and translation into the flux balance analysis problem
* Introduced a workaround (often used in practice) to account for the metabolic cost of gene expression
* Closed the loop: using a metabolic model and flux balance to optimize glycan production in _E.coli_.
"""

# ╔═╡ a818128b-d642-49e0-8c7c-db36b37ea882
md"""
### Alas, if we just had some more time, it would have been cool to talk about:

##### So many interesting tools:

* [Lewis NE, Nagarajan H, Palsson BO. Constraining the metabolic genotype-phenotype relationship using a phylogeny of in silico methods. Nat Rev Microbiol. 2012 Feb 27;10(4):291-305. doi: 10.1038/nrmicro2737. PMID: 22367118; PMCID: PMC3536058.](https://pubmed.ncbi.nlm.nih.gov/22367118/)

##### Better constraints give better flux estimates:

* [Buescher JM, Antoniewicz MR, Boros LG, et al. A roadmap for interpreting (13)C metabolite labeling patterns from cells. Curr Opin Biotechnol. 2015;34:189-201. doi:10.1016/j.copbio.2015.02.003](https://www.ncbi.nlm.nih.gov/labs/pmc/articles/PMC4552607/)

##### Flux balance analysis prediction of mutant response is not correct:

* [Segrè D, Vitkup D, Church GM. Analysis of optimality in natural and perturbed metabolic networks. Proc Natl Acad Sci U S A. 2002 Nov 12;99(23):15112-7. doi: 10.1073/pnas.232349399. Epub 2002 Nov 1. PMID: 12415116; PMCID: PMC137552.](https://pubmed.ncbi.nlm.nih.gov/12415116/)

##### Molecular crowding plays a role in metabolism (non-representative random sample):

* [Beg QK, Vazquez A, Ernst J, de Menezes MA, Bar-Joseph Z, Barabási AL, Oltvai ZN. Intracellular crowding defines the mode and sequence of substrate uptake by Escherichia coli and constrains its metabolic activity. Proc Natl Acad Sci U S A. 2007 Jul 31;104(31):12663-8. doi: 10.1073/pnas.0609845104. Epub 2007 Jul 24. PMID: 17652176; PMCID: PMC1937523.](https://pubmed.ncbi.nlm.nih.gov/17652176/)'
* [Conrado RJ, Mansell TJ, Varner JD, DeLisa MP. Stochastic reaction-diffusion simulation of enzyme compartmentalization reveals improved catalytic efficiency for a synthetic metabolic pathway. Metab Eng. 2007 Jul;9(4):355-63. doi: 10.1016/j.ymben.2007.05.002. Epub 2007 May 26. PMID: 17601761.](https://pubmed.ncbi.nlm.nih.gov/17601761/)
* [Klumpp S, Scott M, Pedersen S, Hwa T. Molecular crowding limits translation and cell growth. Proc Natl Acad Sci U S A. 2013 Oct 15;110(42):16754-9. doi: 10.1073/pnas.1310377110. Epub 2013 Sep 30. PMID: 24082144; PMCID: PMC3801028.](https://pubmed.ncbi.nlm.nih.gov/22760628/)

"""

# ╔═╡ 84d8a1c4-483c-4cc2-bdd0-7c88d26fd319
TableOfContents(title="📚 Lecture Outline", indent=true, depth=5, aside=true)

# ╔═╡ 7988263a-a5f5-11ec-3836-150105d04264
html"""
<style>
main {
    max-width: 900px;
    width: 70%;
    margin: auto;
    font-family: "Roboto, monospace";
}

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

# ╔═╡ dee24b55-a5c5-4a47-9e1c-88ea0cf1cd61
html"""
<script>
	// initialize -
	var section = 0;
	var subsection = 0;
	var subsubsection = 0;
	var headers = document.querySelectorAll('h3, h5, h6');

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

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

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

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

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

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
PlutoUI = "~0.7.37"
"""

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

julia_version = "1.7.2"
manifest_format = "2.0"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

# ╔═╡ Cell order:
# ╟─2c895e14-e71a-435d-86d4-c1b101877b3c
# ╟─ce5c8e9c-7652-4c3e-98a7-98e3567dc695
# ╟─bfd134fe-b311-4768-a53b-2c5c4c1b121a
# ╟─8041cb6c-b80f-4e81-8d1d-5a2e9942d809
# ╟─06c3eafa-10b3-4a47-9a06-39f6abd29db3
# ╟─b15b540b-8f13-4f21-b2da-510d34143b33
# ╟─a818128b-d642-49e0-8c7c-db36b37ea882
# ╠═84d8a1c4-483c-4cc2-bdd0-7c88d26fd319
# ╠═c3530026-5a87-4447-88f5-fc744cb48ff2
# ╟─7988263a-a5f5-11ec-3836-150105d04264
# ╟─dee24b55-a5c5-4a47-9e1c-88ea0cf1cd61
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002

ext: jl


text:
using Dates

struct MinimaxNode
    snakes::Vector{Snake}
    food::Food
    height::Int
    width::Int
    depth::Int
    action::Union{Direction, Nothing}
    winner::Union{Int, Nothing}
end

function iterative_deepening(board::Boardstate, simulation_time::Float64, height::Int, width::Int, you_index::Int)
    simulation_time = Millisecond(floor(Int,simulation_time))
    snakes = board.snakes
    food = board.food
    root = MinimaxNode(snakes, food, height, width, 0, nothing, nothing)
    action_nodes = possible_successors(root, you_index)
    best_action = nothing
    depth = 1

    start = now()
    last_time = Millisecond(0)
    values = Tuple{Float64, MinimaxNode}[]
    while now() - start + last_time < simulation_time
        start_run = now()
        empty!(values)

        for a_node in action_nodes
            if now() - start < simulation_time
                score = minimax(a_node, depth, -Inf, Inf, (you_index % 2) + 1, you_index)
                push!(values, (score, a_node))
            end
        end
        depth += 1
        if length(values) == length(action_nodes)
            sort!(values, by= x -> x[1], rev=true)
            best_action = values[1][2].action
            action_nodes = [v[2] for v in values]
        end
        last_time = now() - start_run
    end
    # println("$depth")
    # println([(v[1], v[2].action) for v in values])
    return best_action
end

function minimax(node::MinimaxNode, depth::Int, alpha::Float64, beta::Float64, player::Int, you_index::Int)
    if depth == 0 || node.winner !== nothing
        return score(node, you_index) # metric or win/lose
    end

    next_player = (player % 2) + 1
    if player == 1
        value = -Inf
        for successor in possible_successors(node, player)
            value = max(value, minimax(successor, depth - 1, alpha, beta, next_player, you_index))
            if value >= beta
                break # beta cutoff
            end
            alpha = max(alpha, value)
        end
        return value
    else
        value = Inf
        for successor in possible_successors(node, player)
            value = min(value, minimax(successor, depth - 1, alpha, beta, next_player, you_index))
            if value <= alpha
                break # alpha cutoff
            end
            beta = min(beta, value)
        end
        return value
    end
end

function possible_successors(node::MinimaxNode, player::Int)
    successors = MinimaxNode[]
    snek = node.snakes[player]
    for a in possible_actions(snek)
        new_snek, new_food = deepcopy(snek), deepcopy(node.food)
        move_snake!(new_snek, new_food, a)
        snakes = deepcopy(node.snakes)
        snakes[player] = new_snek
        winner = check_winner(snakes, player, node.height, node.width)
        n = MinimaxNode(snakes, new_food, node.height, node.width, node.depth + 1, a, winner)
        push!(successors, n)
    end

    return successors
end

function check_winner(snakes, player, height, width)
    opponent = (player % 2) + 1
    snek = snakes[player]
    opponent_snake = snakes[opponent]
    # starve
    if snek.health == 0
        return opponent
    end

    if get_head(snek) in opponent_snake.body
        # head2head
        if get_head(snek) == get_head(opponent_snake)
            # RAM
            if length(snek.body) > length(opponent_snake.body)
                return player
            end
            # draw
            if length(snek.body) == length(opponent_snake.body)
                return -1
            end
        end

        # lost (head2head) collision with opponent_snake
        return opponent
    end
    # collision with self
    if get_head(snek) in get_body_no_head(snek)
        return opponent
    end
    if !check_bounds(get_head(snek), height, width)
        return opponent
    end
end

function score(node::MinimaxNode, player::Int)
    discount = 0.99 ^ node.depth
    h_factor = 2
    if length(node.snakes[player].body) >= 10
        l_factor = -10
    else
        l_factor = 10
    end
    f_factor = 1
    if node.winner == player
        return 200 * discount
    end
    if node.winner == (player % 2) + 1
        return -200 * discount
    end
    # draw
    if node.winner == -1
        return -50
    end
    snek = node.snakes[player]
    health_diff = 2 * snek.health - sum([s.health for s in node.snakes])
    len_diff = 2 * length(node.snakes[player].body) - sum([length(s.body) for s in node.snakes])
    flood = floodfill(node.snakes, node.width, node.height, player)

    return h_factor * health_diff + l_factor * len_diff + f_factor * flood * discount
end

max_stars_repo_name: rhotertj/solide-schlange

max_issues_repo_name: rhotertj/solide-schlange

max_forks_repo_name: rhotertj/solide-schlange


text:
using Plots

function plotsegment!(pl1, pl2, leafid, segs, col)
    i = leafid
    atroot = false
    while !atroot
        gg = segs[i]
        t = gg[1]
        u = gg[2]
        plot!(pl1, t, getindex.(u,1),color=col)
        plot!(pl2, t, getindex.(u,2),color=col)
        ipar = tree.Par[i]
        i = ipar
        if i==1
            atroot=true
        end
    end
end

function plotmessage!(pl1, pl2, leafid, messages, col)
    i = leafid
    atroot = false
    while !atroot
        m = messages[i]
        t = m.ts
        u = [u[1] for u in m.soldis] # mean value
        plot!(pl1, t, getindex.(u,1),color=col)
        plot!(pl2, t, getindex.(u,2),color=col)
        ipar = tree.Par[i]
        i = ipar
        if i==1
            atroot=true
        end
    end
end

## trace plots of pars
fig1 = plot(0:iters, getindex.(first.(θs), 1), color=:blue,  title="θ₁")
fig2 = plot(0:iters, getindex.(first.(θs), 2), color=:blue, xlabel="iterations", title="θ₂")
plot!(fig1, 0:iters, fill(θ0[1][1], iters+1), lw=2.0, color=:green)
plot!(fig2, 0:iters, fill(θ0[1][2], iters+1), lw=2.0, color=:green)
pl = plot(fig1, fig2, layout = (2, 1), legend = false)
savefig(pl, "./Figures/trace_theta.png")

## trace plots of pars
fig1 = plot(0:iters, getindex.(last.(θs), 1), color=:blue, title="σ₁")
fig2 = plot(0:iters, getindex.(last.(θs), 2), color=:blue, xlabel="iterations", title="σ₂")
plot!(fig1, 0:iters, fill(θ0[2][1], iters+1), lw=2.0, color=:green)
plot!(fig2, 0:iters, fill(θ0[2][2], iters+1), lw=2.0, color=:green)
pl = plot(fig1, fig2, layout = (2, 1), legend = false)
savefig(pl, "./Figures/trace_sigma.png")

## forward simulated paths
cols = repeat([:blue, :red, :magenta, :orange],2)
pl1 = scatter(tree.T[tree.lids], getindex.(Xd,1)[tree.lids],color=:black,legend=false, title="x₁")
pl2 = scatter(tree.T[tree.lids], getindex.(Xd,2)[tree.lids],color=:black,legend=false, xlabel="time", title="x₂")
for i in 2:tree.n
    global gg
    gg = segs[i]
    t = gg[1]
    u = gg[2]
    col = sample(cols)
    plot!(pl1, t, getindex.(u,1),color=col)
    plot!(pl2, t, getindex.(u,2),color=col)
end
pl = plot(pl1, pl2, layout=(2,1), legend=false)
display(pl)
savefig(pl, "./Figures/dataforward.png")

## just the observations
cols = repeat([:blue, :red, :magenta, :orange],2)
pl1 = scatter(tree.T[tree.lids], getindex.(Xd,1)[tree.lids],color=:black,legend=false, title="x₁")
pl2 = scatter(tree.T[tree.lids], getindex.(Xd,2)[tree.lids],color=:black,legend=false,  xlabel="time",  title="x₂")
for i in 2:tree.n
    global gg
    gg = segs[i]
    t = gg[1]
    u = gg[2]
    col = :white
    plot!(pl1, t, 0*t,color=col)
    plot!(pl2, t, 0*t,color=col)
end
pl = plot(pl1, pl2, layout=(2,1), legend=false)
display(pl)
savefig(pl, "./Figures/observations.png")

## some segments of guided process at final iteration
leaves = Xd[tree.lids]
tleaves = tree.T[tree.lids]
pl1 = scatter(tleaves, getindex.(leaves,1) ,color=:black, legend=false, title="x₁")
pl2 = scatter(tleaves, getindex.(leaves,2) ,color=:black, legend=false,
 xlabel="time",  title="x₂")

cols = repeat([:blue, :red, :green, :magenta, :orange],2)
for i ∈ 0:3
    plotsegment!(pl1, pl2, tree.lids[end-i], guidedsegs, cols[i+1])
end
pl = plot(pl1, pl2, layout = (2, 1), legend = false)
display(pl)
savefig(pl, "./Figures/guidedforward.png")

## plot expected tree at final iteration
leaves = Xd[tree.lids]
tleaves = tree.T[tree.lids]
pl1 = scatter(tleaves, getindex.(leaves,1) ,color=:black, legend=false, title="x₁");
pl2 = scatter(tleaves, getindex.(leaves,2) ,color=:black, legend=false, title="x₂",  xlabel="time");

θpred = ([mean(getindex.(θs1,1)), mean(getindex.(θs1,2))], [mean(getindex.(θs2,1)), mean(getindex.(θs2,2))] )
Z = [myinnov(forwardguiding_input[4][i].ts, 𝕏, 0) for i ∈ 2:tree.n]
Xᵒ, guidedsegsᵒ, llᵒ, 𝐋ᵒ = fwguidtree!(forwardguiding_input[1],
                                     forwardguiding_input[2],
                                     forwardguiding_input[3],
                                     forwardguiding_input[4],
                                     tree, f, g, θpred, Z, EM(false); apply_time_change=false)

for i ∈ 0:3
  plotsegment!(pl1, pl2, tree.lids[end-i], guidedsegsᵒ, cols[i+1])
end
pl = plot(pl1, pl2, layout = (2, 1), legend = false)
display(pl)
savefig(pl, "./Figures/expected_guidedforward.png")

# plot messages
# cols = repeat([:blue, :red, :green, :magenta, :orange],2)
# for i ∈ 0:3
#     plotmessage!(pl1, pl2, tree.lids[end-i], messages, cols[i+1])
# end
# pl = plot(pl1, pl2, layout = (2, 1), legend = false)
# display(pl)
# savefig(pl, "./Figures/messages.png")


text:
using MultipleTesting
using StatsBase
using Distributions
using CategoricalArrays
using IndependentHypothesisWeighting

Ps = rand(BetaUniformMixtureModel(0.7, 0.3),10000)

sum(adjust(Ps, BenjaminiHochberg()) .<= 0.1)

ihw_opt = IHW(folds=2, weight_learner=IndependentHypothesisWeighting.UnitWeightLearner(), α=0.1)
ihw_fit = fit(ihw_opt, Ps, rand(10000))

sum(adjust(ihw_fit) .<= 0.1)

Xs = CategoricalVector(sample(1:2, 10000))
Ps1 = rand(BetaUniformMixtureModel(0.7, 0.2),10000)
Ps2 = rand(Uniform(),10000)
Ps = Ps1 .* (Xs.==1) .+ Ps2 .* (Xs.==2)

ihw_opt = IHW(folds=2, weight_learner=GBHLearner(Storey(0.5)), α=0.1)
ihw_fit = fit(ihw_opt, Ps, Xs)

sum(adjust(ihw_fit) .<= 0.1)

sum(adjust(Ps, BenjaminiHochberg()) .<= 0.1)

ihw_grenander = IHW(folds=2, weight_learner=GrenanderLearner(), α=0.1)
ihw_grenander_fit = fit(ihw_grenander, Ps, Xs)

unique(ihw_grenander_fit.weights)
sum(adjust(ihw_grenander_fit) .<= 0.1)

ihw_discr = IHW(folds=2, weight_learner = DiscretizationWeightLearner(GrenanderLearner(), QuantileSlicingDiscretizer(2)))
Xs_cont = rand(10000) .+ convert(Vector{Float64}, Xs)

ihw_grenander_fit_discr = fit(ihw_discr, Ps, Xs_cont)

# now compare weights above to implementation based on JuMP

using JuMP
using Clp
using LinearAlgebra
using MLDataPattern

for fold_idx in [1;2]
	train_idx, test_idx = ihw_grenander_fit.kfolds[fold_idx]

	Ps_train = Ps[train_idx]
	Xs_train = Xs[train_idx]

	Xs_test = Xs[test_idx]
	m1_test = sum(Xs_test .== 1)
	m2_test = sum(Xs_test .== 2)
	ms_test = [m1_test; m2_test]
	Ps_train_1 = Ps_train[Xs_train .== 1]
	Ps_train_2 = Ps_train[Xs_train .== 2]

	gren_fit_1 = fit(IndependentHypothesisWeighting.Grenander(), Ps_train_1)
	gren_fit_2 = fit(IndependentHypothesisWeighting.Grenander(), Ps_train_2)

	mymodel = Model(Clp.Optimizer)

	@variable(mymodel, ts[1:2] >= 0)
	@variable(mymodel, Fs[1:2] >= 0)

	offsets_1 = gren_fit_1.Fs[1:(end-1)] .- gren_fit_1.fs .* gren_fit_1.locs[1:(end-1)]
	@constraint(mymodel,  Fs[1] .<= offsets_1 .+ gren_fit_1.fs .* ts[1])

	offsets_2 = gren_fit_2.Fs[1:(end-1)] .- gren_fit_2.fs .* gren_fit_2.locs[1:(end-1)]
	@constraint(mymodel,  Fs[2] .<= offsets_2 .+ gren_fit_2.fs .* ts[2])

	@constraint(mymodel, dot(ms_test, ts .- 0.1 .* Fs) <= 0)

	@objective(mymodel, Max, dot(ms_test, Fs))

	optimize!(mymodel)

	solver_ts = JuMP.value.(ts)

	solver_ws = solver_ts * sum(ms_test) ./ sum(ms_test .* solver_ts)
	@test sort(collect(ihw_grenander_fit.weighting_fits[fold_idx])) ≈ sort(solver_ws)
end


text:
module ExportWebAssembly

export irgen, write_bitcode, @extern

using LLVM

const jlctx = Ref{LLVM.Context}()

function __init__()
    jlctx[] = LLVM.Context(convert(LLVM.API.LLVMContextRef,
                                   cglobal(:jl_LLVMContext, Nothing)))
end

function module_setup(mod::LLVM.Module)
    triple!(mod, "wasm32-unknown-unknown-wasm")
    datalayout!(mod, "e-m:e-p:32:32-i64:64-n32:64-S128")
end

# make function names safe
safe_fn(fn::String) = replace(fn, r"[^aA-zZ0-9_]", "_")
safe_fn(f::Core.Function) = safe_fn(String(typeof(f).name.mt.name))
safe_fn(f::LLVM.Function) = safe_fn(LLVM.name(f))

# needed?
function raise_exception(insblock::BasicBlock, ex::Value)
end

function irgen(@nospecialize(func), @nospecialize(tt); optimize = true)
    # get the method instance
    isa(func, Core.Builtin) && error("function is not a generic function")
    world = typemax(UInt)
    meth = which(func, tt)
    sig_tt = Tuple{typeof(func), tt.parameters...}
    (ti, env) = ccall(:jl_type_intersection_with_env, Any,
                      (Any, Any), sig_tt, meth.sig)::Core.SimpleVector
    meth = Base.func_for_method_checked(meth, ti)
    linfo = ccall(:jl_specializations_get_linfo, Ref{Core.MethodInstance},
                  (Any, Any, Any, UInt), meth, ti, env, world)

    # set-up the compiler interface
    function hook_module_setup(ref::Ptr{Cvoid})
        ref = convert(LLVM.API.LLVMModuleRef, ref)
        module_setup(LLVM.Module(ref))
    end
    function hook_raise_exception(insblock::Ptr{Cvoid}, ex::Ptr{Cvoid})
        insblock = convert(LLVM.API.LLVMValueRef, insblock)
        ex = convert(LLVM.API.LLVMValueRef, ex)
        raise_exception(BasicBlock(insblock), Value(ex))
    end
    dependencies = Vector{LLVM.Module}()
    function hook_module_activation(ref::Ptr{Cvoid})
        ref = convert(LLVM.API.LLVMModuleRef, ref)
        push!(dependencies, LLVM.Module(ref))
    end
    params = Base.CodegenParams(cached=false,
                                track_allocations=false,
                                code_coverage=false,
                                static_alloc=false,
                                prefer_specsig=true,
                                module_setup=hook_module_setup,
                                module_activation=hook_module_activation,
                                raise_exception=hook_raise_exception)

    # get the code
    mod = let
        ref = ccall(:jl_get_llvmf_defn, LLVM.API.LLVMValueRef,
                    (Any, UInt, Bool, Bool, Base.CodegenParams),
                    linfo, world, #=wrapper=#false, #=optimize=#false, params)
        if ref == C_NULL
            error(jlctx[], "the Julia compiler could not generate LLVM IR")
        end

        llvmf = LLVM.Function(ref)
        LLVM.parent(llvmf)
    end

    # the main module should contain a single jfptr_ function definition,
    # e.g. jlcall_kernel_vadd_62977
    @show isdeclaration(first(functions(mod)))
    # @show length(functions(mod))
    # definitions = filter(f->!isdeclaration(f), functions(mod))
    definitions = filter(f->!isdeclaration(f), collect(functions(mod)))
    # definitions = collect(functions(mod))
    wrapper = let
        fs = collect(filter(f->startswith(LLVM.name(f), "jfptr_"), definitions))
        @assert length(fs) == 1
        fs[1]
    end

    # the jlcall wrapper function should point us to the actual entry-point,
    # e.g. julia_kernel_vadd_62984
    entry_tag = let
        m = match(r"jfptr_(.+)_\d+", LLVM.name(wrapper))
        @assert m != nothing
        m.captures[1]
    end
    unsafe_delete!(mod, wrapper)
    entry = let
        re = Regex("julia_$(entry_tag)_\\d+")
        llvmcall_re = Regex("julia_$(entry_tag)_\\d+u\\d+")
        fs = collect(filter(f->occursin(re, LLVM.name(f)) &&
                               !occursin(llvmcall_re, LLVM.name(f)), definitions))
        if length(fs) != 1
            compiler_error(f, tt, cap, "could not find single entry-point";
                           entry=>entry_tag, available=>[LLVM.name.(definitions)])
        end
        fs[1]
    end

    LLVM.name!(entry, string(nameof(func)))

    # link in dependent modules
    link!.(Ref(mod), dependencies)

    # clean up incompatibilities
    for llvmf in functions(mod)
        # only occurs in debug builds
        delete!(function_attributes(llvmf), EnumAttribute("sspreq", 0, jlctx[]))

        # make function names safe for ptxas
        # (LLVM ought to do this, see eg. D17738 and D19126), but fails
        # TODO: fix all globals?
        llvmfn = LLVM.name(llvmf)
        if !isdeclaration(llvmf)
            llvmfn′ = safe_fn(llvmf)
            if llvmfn != llvmfn′
                LLVM.name!(llvmf, llvmfn′)
            end
        end
    end

    # tm = machine(jlctx[].cap, triple(mod))
    triple = "wasm32-unknown-unknown-wasm"
    triple = "x86_64-unknown-linux-gnu"
    tm = TargetMachine(Target(triple), triple)

    # GPU code is _very_ sensitive to register pressure and local memory usage,
    # so forcibly inline every function definition into the entry point
    # and internalize all other functions (enabling ABI-breaking optimizations).
    # FIXME: this is too coarse. use a proper inliner tuned for GPUs
    optimize && ModulePassManager() do pm
        no_inline = EnumAttribute("noinline", 0, jlctx[])
        always_inline = EnumAttribute("alwaysinline", 0, jlctx[])
        for f in filter(f->f!=entry && !isdeclaration(f), collect(functions(mod)))
            attrs = function_attributes(f)
            if !(no_inline in collect(attrs))
                push!(attrs, always_inline)
            end
            linkage!(f, LLVM.API.LLVMInternalLinkage)
        end
        always_inliner!(pm)
        run!(pm, mod)
    end

    optimize && ModulePassManager() do pm
        # add_library_info!(pm, triple(mod))
        add_transform_info!(pm, tm)
        ccall(:jl_add_optimization_passes, Cvoid,
              (LLVM.API.LLVMPassManagerRef, Cint),
              LLVM.ref(pm), Base.JLOptions().opt_level)

        # CUDAnative's JIT internalizes non-inlined child functions, making it possible
        # to rewrite them (whereas the Julia JIT caches those functions exactly);
        # this opens up some more optimization opportunities
        dead_arg_elimination!(pm)   # parent doesn't use return value --> ret void

        global_optimizer!(pm)
        global_dce!(pm)
        strip_dead_prototypes!(pm)

        run!(pm, mod)
    end
    return mod
end

jltype(x::Symbol) = eval(x)
llvmtype(x) = LLVMType(ccall(:julia_type_to_llvm, LLVM.API.LLVMTypeRef, (Any, Bool), x, false))
llvmtype(x::Symbol) = llvmtype(jltype(x))

"""
    @extern(fun, returntype, argtypes, args...)

Creates a call to an external function meant to be defined at link time.

The symbol `fun` is the external function with a return type `returntype`
and arguments given as a tuple type given by `argtypes`.

Note that the types must be base types handled by LLVM (mainly bits types).
Julia boxed types cannot be used. This includes arrays.

Note that arguments are not autoconverted to the proper type like `ccall`.
"""
macro extern(funname, returntype, argtypes, args...)
    funname = funname.value
    llvmreturntype = llvmtype(returntype)
    llvmargtypes = llvmtype.(tuple(argtypes.args[2:end]...))
    llvmargs = join(("$a %$(idx-1)" for (idx, a) in enumerate(llvmargtypes)), ", ")
    llvmargtypes = join(llvmargtypes, ", ")
    declarationstr = "declare $llvmreturntype @$funname($llvmargtypes)"
    runstr = """
        %ret = call $llvmreturntype @$funname($llvmargs)
        ret $llvmreturntype %ret
        """
    esc(Expr(:call, Base.llvmcall, (declarationstr, runstr), returntype, Tuple{jltype.(argtypes.args[2:end])...}, args...))
end

module Utils
    include("../deps/deps.jl")
end
include("io.jl")

end # module


text:
using OptEnrichedSetCover
using Test, LinearAlgebra, DataFrames

const OESC = OptEnrichedSetCover

# run the tests from jlfile
macro testfile(jlfile)
    quote
        @testset "\"$($jlfile)\" tests" begin
            include($jlfile)
        end
    end
end

@testfile "test_set_score.jl"
@testfile "test_sparse_mask_matrix.jl"
@testfile "test_mosaic.jl"
@testfile "test_masked_mosaic.jl"
@testfile "test_weighted_mosaic.jl"
@testfile "test_multiobj_cover_problem.jl"
@testfile "test_cover_enumerator.jl"
@testfile "test_parallel.jl"


text:
# empirical data for Germany after 9.3. scaled to 50K HHs
#emptotinf = [4.711538462,	5.837339744,	7.644230769,	8.325320513,	14.72355769,	18.36939103,	23.21714744,	29.13461538,	37.08733974,	49.38701923,	61.37820513,	79.51923077,	88.99439103,	99.65144231,	116.4102564,	132.1554487,	149.53125,	176.0336538,	203.8100962,	231.1498397,	248.7780449,	267.96875,	287.6923077,	311.9871795,	339.7195513,	365.2203526,	384.9839744,	401.1338141,	414.1586538,	431.3421474,	453.9102564,	473.4815705, 489.4671474,	500.4326923,	512.2355769,	521.1217949,	526.2780449,	539.8758013,	551.6746795,	566.494391,	574.286859,	581.6666667,	589.2027244,	594.1145833,	603.5576923,	613.4975962,	620.9895833,	627.0552885,	632.0913462,	636.0496795,	640.6730769,	647.1915064,	653.0809295,	657.3597756,	660.9254808,	663.7179487,	665.6730769,	669.0985577,	673.7259615,	678.8060897,	683.4455128]

# empirical data smoothed R0
empR0sm = [3.111478515,	2.587691482,	2.400553023,	2.187280465,	2.057180151,	2.457779377,	2.850747778,	3.483690962,	3.663591612,	3.438431489,	3.13360927,	2.916868444,	2.380304109,	2.257667713,	2.051074872,	1.769917342,	1.518917287,	1.358632478,	1.352162677,	1.439888262,	1.558839418,	1.566312254,	1.454341178,	1.30002569,	1.107143519,	0.990820486,	0.982068753,	1.010676322,	1.011066964,	0.963784487,	0.893908707,	0.824248406,	0.80602619,	0.851298179,	0.884557006,	0.874564867,	0.810013341,	0.721855519,	0.648114235,	0.669354718,	0.798497082,	0.912641321,	0.988772116,	0.975282103,	0.906741323,	0.798235484,	0.755734188,	0.803951625,	0.877166325,	0.896926093,	0.86843844,	0.821292575,	0.729887915,	0.723910481,	0.787991881,	0.828291584,	0.826018956,	0.789068446,	0.726492331,	0.703305688,	0.821437328]

# empirical data Germany after 9.3. scaled to 100k HHs in percent
empcas = [0.0000025,	0.0000025,	0.00000375,	0.00000375,	0.00000875,	0.00001125,	0.00001375,	0.00002125,	0.00003,	0.000035,	0.000055,	0.00008375,	0.000105,	0.0001175,	0.00015375,	0.00019625,	0.0002575,	0.00033375,	0.0004275,	0.00054125,	0.00066625,	0.00080625,	0.00096875,	0.00115,	0.00138375,	0.00159375,	0.001805,	0.00198,	0.0022625,	0.00252,	0.00293625,	0.00325875,	0.00345875,	0.00342,	0.0037775,	0.0039925,	0.0041175,	0.004755,	0.005065,	0.00544,	0.00557375,	0.0057325,	0.0060775,	0.00629125,	0.00659875,	0.00696875,	0.0072,	0.00734625,	0.00747,	0.0076575,	0.0078925,	0.00808375,	0.00827875,	0.00842,	0.008515,	0.0085825,	0.00874125,	0.00874125,	0.00909375,	0.00924,	0.0093875]

# empirical data for Germany after 9.3. scaled to 100K HHs
emptotinf = [9.423076923,	11.67467949,	15.28846154,	16.65064103,	29.44711538,	36.73878205,	46.43429487,	58.26923077,	74.17467949,	98.77403846,	122.7564103,	159.0384615,	177.9887821,	199.3028846,	232.8205128,	264.3108974,	299.0625,	352.0673077,	407.6201923,	462.2996795,	497.5560897,	535.9375,	575.3846154,	623.974359,	679.4391026,	730.4407051,	769.9679487,	802.2676282,	828.3173077,	862.6842949,	907.8205128,	946.963141,	978.9342949,	1000.865385,	1024.471154,	1042.24359,	1052.55609,	1079.751603,	1103.349359,	1132.988782,	1148.573718,	1163.333333,	1178.405449,	1188.229167,	1207.115385,	1226.995192,	1241.979167,	1254.110577,	1264.182692,	1272.099359,	1281.346154,	1294.383013,	1306.161859,	1314.719551,	1321.850962,	1327.435897,	1331.346154,	1338.197115,	1347.451923, 1357.6121794, 1366.8910256]


text:
# Inserting
# ---------

# Split a leaf into two, returning (leftnode, rightnode)
function split!(left::LeafNode{K, V, B}) where {K, V, B}
    right = LeafNode{K, V, B}()
    right.right = left.right
    right.left = left
    left.right = right
    if right.right !== nothing
        right.right = notnothing(right.right)
        right.right.left = right
    end
    right.parent = left.parent

    m = length(left)
    resize!(right, m - div(m, 2))
    copyto!(right.entries, 1, left.entries, div(m, 2) + 1, length(right))
    copyto!(right.keys, 1, left.keys, div(m, 2) + 1, length(right))
    resize!(left, div(m, 2))

    left.maxend = last(left.keys[1])
    for entry in left
        left.maxend = max(left.maxend, last(entry))
    end

    right.maxend = last(right.keys[1])
    for entry in right
        right.maxend = max(right.maxend, last(entry))
    end

    return left, right
end

# Split an internal node in two, returning (leftnode, rightnode)
function split!(left::InternalNode{K, V, B}) where {K, V, B}
    right = InternalNode{K, V, B}()
    right.right = left.right
    right.left = left
    left.right = right
    if right.right !== nothing
        right.right = notnothing(right.right)
        right.right.left = right
    end
    right.parent = left.parent

    m = length(left)

    resize!(right.children, m - div(m, 2))
    resize!(right.maxends, m - div(m, 2))
    resize!(right.keys, m - div(m, 2) - 1)

    copyto!(right.children, 1, left.children, div(m, 2)+1, length(right.children))
    copyto!(right.maxends, 1, left.maxends, div(m, 2)+1, length(right.maxends))
    copyto!(right.keys, 1, left.keys, div(m, 2)+1, length(right.keys))

    resize!(left.children, div(m, 2))
    resize!(left.maxends, div(m, 2))
    resize!(left.keys, div(m, 2) - 1)

    left.maxend = left.children[1].maxend
    for child in left.children
        left.maxend = max(left.maxend, child.maxend)
    end

    right.maxend = right.children[1].maxend
    for child in right.children
        right.maxend = max(right.maxend, child.maxend)
    end

    for child in right.children
        child.parent = right
    end

    return left, right
end

# Find the maximum Segment end point is a subtree
function nodemaxend(t::InternalNode{K, V, B}) where {K, V, B}
    maxend = zero(K)
    for i in 1:length(t.children)
        maxend = max(maxend, t.children[i].maxend)
    end
    return maxend
end

function nodemaxend(t::LeafNode{K, V, B}) where {K, V, B}
    maxend = zero(K)
    for i in 1:length(t)
        maxend = max(maxend, last(t.keys[i]))
    end
    return maxend
end

sameparent(u::Node{K, V, B}, v::Node{K, V, B}) where {K, V, B} = (u.parent == v.parent)

# Find the first leaf node in the tree
function firstleaf(t::SegmentBTree)
    return firstleaf(t.root)
end

function firstleaf(t::InternalNode)
    return firstleaf(t.children[1])
end

function firstleaf(t::LeafNode)
    return t
end

# Find the root node
function root(t::Node)
    while t.parent !== nothing
        t = notnothing(t.parent)
    end
    return t
end

# This version of nextleafkey is used in nextintersection at aggressively avoid
# allocation, and dereferecing nullables more than once.
macro nextleafkey(leaf, nleaf, i)
    quote
        if $(esc(i)) < length($(esc(leaf)))
            $(esc(i)) += 1
        else
            $(esc(i)) = 1
            $(esc(nleaf)) = $(esc(leaf)).right
        end
    end
end

# Find the minimum Segment in a subtree
function minkey(t::InternalNode)
    return minkey(t.children[1])
end

function minkey(t::LeafNode{K, V, B}) where {K, V, B}
    return t.keys[1]
end

function Base.push!(t::SegmentBTree{K, V, B}, entry::V,
                    unique_key::Bool=false, update::Bool=true) where {K, V, B}
    return _push!(t, t.root, entry, unique_key, update)
end

function _push!(t::SegmentBTree{K, V, B},
                node::InternalNode{K, V, B},
                entry::V, unique_key::Bool, update::Bool) where {K, V, B}
    i = findidx(node, entry) # key index
    j = i <= length(node) - 1 && entry < node.keys[i] ? i : i + 1 # child index
    return _push!(t, node.children[j], entry, unique_key, update)
end

function _push!(t::SegmentBTree{K, V, B},
                node::LeafNode{K, V, B}, entry::V,
                unique_key::Bool, update::Bool) where {K, V, B}
    i = max(1, findidx(node, entry))

    # key exists and we are replacing it
    if i <= length(node) && unique_key &&
        first(node.keys[i]) == first(entry) && last(node.keys[i]) == last(entry)
        if update
            node.entries[i] = entry
            return entry
        else
            return node.entries[i]
        end
    end

    insert!(node, i, entry)
    if length(node) == 1 || last(entry) > node.maxend
        node.maxend = last(entry)
    end
    t.n += 1

    # split when full
    if length(node) == B
        maxend = node.maxend
        leftleaf, rightleaf = split!(node)
        median = rightleaf.keys[1]

        # travel back up the tree setting maxend values and splitting as needed
        parent = node.parent
        leftnode = leftleaf
        child = leftleaf
        rightnode = rightleaf
        hassplit = true
        while parent !== nothing
            p = notnothing(parent)
            if hassplit
                i = findidx(p, entry) # key index
                j = i <= length(p) - 1 && entry < p.keys[i] ? i : i + 1 # child index
                p.maxends[j] = p.children[j].maxend
                insert!(p.children, j + 1, rightnode)
                insert!(p.maxends, j + 1, p.children[j+1].maxend)
                p.maxend = max(p.maxend, maxend)
                insert!(p.keys, j, median)

                # split when full
                if length(p) == B
                    leftnode, rightnode = split!(p)
                    median = minkey(rightnode)
                    maxend = max(leftnode.maxend, rightnode.maxend)
                else
                    hassplit = false
                end
            else
                p.maxend = max(p.maxend, last(entry))
                ifind = findfirst(isequal(child), p.children)
                if ifind === nothing
                    ifind = 0
                end
                p.maxends[ifind] = child.maxend
            end
            child = p
            parent = p.parent
        end

        if hassplit
            t.root = InternalNode{K, V, B}()
            push!(t.root.keys, median)
            push!(t.root.children, leftnode)
            push!(t.root.children, rightnode)
            t.root.maxend = maxend
            push!(t.root.maxends, leftnode.maxend)
            push!(t.root.maxends, rightnode.maxend)
            leftnode.parent = t.root
            rightnode.parent = t.root
        end
    else
        # travel back up the tree setting maxend values
        parent = node.parent
        child = node
        while parent !== nothing
            p = notnothing(parent)
            p.maxend = max(p.maxend, last(entry))
            ifind = findfirst(isequal(child), p.children)
            if ifind === nothing
                ifind = 0
            end
            p.maxends[ifind] = child.maxend
            parent = p.parent
            child = p
        end
    end

    return entry
end


text:
# 2D SLAM with velocity states

"""
$(TYPEDEF)
"""
mutable struct DynPose2VelocityPrior{T1,T2} <: IncrementalInference.FunctorSingleton where {T1 <: IIF.SamplableBelief,T2 <: IIF.SamplableBelief}
  Zpose::T1
  Zvel::T2
  DynPose2VelocityPrior{T1,T2}() where {T1 <: IIF.SamplableBelief,T2 <: IIF.SamplableBelief} = new{T1,T2}()
  DynPose2VelocityPrior{T1,T2}(z1::T1,z2::T2) where {T1 <: IIF.SamplableBelief,T2 <: IIF.SamplableBelief} = new{T1,T2}(z1,z2)
end
DynPose2VelocityPrior(z1::T1,z2::T2) where {T1 <: IIF.SamplableBelief, T2 <: IIF.SamplableBelief} = DynPose2VelocityPrior{T1,T2}(z1,z2)
getSample(dp2v::DynPose2VelocityPrior{T1,T2}, N::Int=1) where {T1 <: IIF.SamplableBelief, T2 <: IIF.SamplableBelief} = ([rand(dp2v.Zpose,N);rand(dp2v.Zvel,N)], )

"""
$(TYPEDEF)
"""
mutable struct VelPose2VelPose2{T1,T2} <: IncrementalInference.FunctorPairwiseMinimize where {T1 <: IIF.SamplableBelief, T2 <: IIF.SamplableBelief}
  Zpose::Pose2Pose2{T1} #Zpose::T1
  Zvel::T2
  reuseres::Vector{Vector{Float64}}
  VelPose2VelPose2{T1,T2}() where {T1 <: IIF.SamplableBelief, T2 <: IIF.SamplableBelief} = new{T1,T2}()
  VelPose2VelPose2{T1,T2}(z1::T1, z2::T2) where {T1 <: IIF.SamplableBelief, T2 <: IIF.SamplableBelief} = new{T1,T2}(Pose2Pose2(z1),z2,[zeros(3) for i in 1:Threads.nthreads()])
end
VelPose2VelPose2(z1::T1, z2::T2) where {T1 <: IIF.SamplableBelief, T2 <: IIF.SamplableBelief} = VelPose2VelPose2{T1,T2}(z1, z2)
getSample(vp2vp2::VelPose2VelPose2, N::Int=1) = ([rand(vp2vp2.Zpose.z,N);rand(vp2vp2.Zvel,N)], )
function (vp2vp2::VelPose2VelPose2{T1,T2})(
                res::Array{Float64},
                userdata::FactorMetadata,
                idx::Int,
                meas::Tuple,
                Xi::Array{Float64,2},
                Xj::Array{Float64,2}  ) where {T1 <: IIF.SamplableBelief, T2 <: IIF.SamplableBelief}
  #
  z = meas[1][:,idx]
  wxi, wxj = Xi[:,idx], Xj[:,idx]
  # @show z, wxi, wxj
  dt = (userdata.variableuserdata[2].ut - userdata.variableuserdata[1].ut)*1e-6   # roughly the intended use of userdata
  # fill!(vp2vp2.reuseres[Threads.threadid()], 0.0)
  vp2vp2.Zpose(vp2vp2.reuseres[Threads.threadid()], userdata, idx, meas, Xi, Xj)
  wDXij = (wxj[4:5]-wxi[4:5])
  bDXij = TransformUtils.R(-wxi[3])*wDXij
  # calculate the residual
  res[1] = sum((vp2vp2.reuseres[Threads.threadid()]).^2)
  res[1] += sum((z[4:5] - bDXij).^2)
  dx = se2vee(SE2(wxi[1:3]) \ SE2(wxj[1:3]))
  res[1] += sum((dx[1:2]/dt - 0.5*(wxj[4:5]+wxi[4:5])).^2)  # first order integration
  res[1]
end

# import RoME: VelPose2VelPose2, IIF.SamplableBelief
# import IncrementalInference: FactorMetadata

"""
$(TYPEDEF)
"""
mutable struct DynPose2Pose2{T} <: IncrementalInference.FunctorPairwise where {T <: IIF.SamplableBelief}
  Zpose::Pose2Pose2{T} #Zpose::T1
  # reuseres::Vector{Float64}
  partial::Tuple{Int,Int,Int}
  DynPose2Pose2{T}() where {T <: IIF.SamplableBelief} = new{T}()
  DynPose2Pose2{T}(z1::T) where {T <: IIF.SamplableBelief} = new{T}(Pose2Pose2(z1), (1,2,3))
end
DynPose2Pose2(z1::T) where {T <: IIF.SamplableBelief} = DynPose2Pose2{T}(z1)
getSample(vp2vp2::DynPose2Pose2, N::Int=1) = (rand(vp2vp2.Zpose.z,N), )
function (vp2vp2::DynPose2Pose2{T})(
                res::Array{Float64},
                userdata::FactorMetadata,
                idx::Int,
                meas::Tuple,
                Xi::Array{Float64,2},
                Xj::Array{Float64,2}  ) where {T <: IIF.SamplableBelief}
  #
  vp2vp2.Zpose(res, userdata, idx, meas, Xi, Xj)
  nothing
end

## Compare functions

function compare(a::MvNormal, b::MvNormal; tol::Float64=1e-10)::Bool
  TP = true
  TP = TP && norm(a.μ - b.μ)<tol
  TP = TP && sum(norm.(a.Σ.mat - b.Σ.mat))<tol
  return TP
end

function compare(a::DynPose2VelocityPrior, b::DynPose2VelocityPrior)::Bool
  RoME.compare(a.Zpose, b.Zpose) && RoME.compare(a.Zvel, b.Zvel)
end

function compare(a::VelPose2VelPose2, b::VelPose2VelPose2; tol::Float64=1e-10)::Bool
  TP = true
  TP = TP && RoME.compare(a.Zpose, b.Zpose)
  TP = TP && RoME.compare(a.Zvel, b.Zvel)
  TP = TP && norm(a.reuseres - b.reuseres) < tol
  return TP
end

## Packing types

"""
$(TYPEDEF)
"""
mutable struct PackedDynPose2VelocityPrior <: IncrementalInference.PackedInferenceType
  strpose::AbstractString
  strvel::AbstractString
  PackedDynPose2VelocityPrior() = new()
  PackedDynPose2VelocityPrior(z1::AS, z2::AS) where {AS <: AbstractString} = new(z1, z2)
end

function convert(::Type{PackedDynPose2VelocityPrior}, d::DynPose2VelocityPrior)
  return PackedDynPose2VelocityPrior(string(d.Zpose),string(d.Zvel))
end
function convert(::Type{DynPose2VelocityPrior}, d::PackedDynPose2VelocityPrior)
  posedistr = extractdistribution(d.strpose)
  veldistr = extractdistribution(d.strvel)
  return DynPose2VelocityPrior(posedistr, veldistr)
end

"""
$(TYPEDEF)
"""
mutable struct PackedVelPose2VelPose2 <: IncrementalInference.PackedInferenceType
  strpose::AbstractString
  strvel::AbstractString
  PackedVelPose2VelPose2() = new()
  PackedVelPose2VelPose2(z1::AS, z2::AS) where {AS <: AbstractString} = new(z1, z2)
end

function convert(::Type{PackedVelPose2VelPose2}, d::VelPose2VelPose2)
  return PackedVelPose2VelPose2(string(d.Zpose.z),string(d.Zvel))
end
function convert(::Type{VelPose2VelPose2}, d::PackedVelPose2VelPose2)
  posedistr = extractdistribution(d.strpose)
  veldistr = extractdistribution(d.strvel)
  return VelPose2VelPose2(posedistr, veldistr)
end

"""
$(TYPEDEF)
"""
mutable struct PackedDynPose2Pose2 <: IncrementalInference.PackedInferenceType
  strpose::AbstractString
  PackedDynPose2Pose2() = new()
  PackedDynPose2Pose2(z1::AS) where {AS <: AbstractString} = new(z1)
end

function convert(::Type{PackedDynPose2Pose2}, d::DynPose2Pose2)
  return PackedDynPose2Pose2(string(d.Zpose.z))
end
function convert(::Type{DynPose2Pose2}, d::PackedDynPose2Pose2)
  posedistr = extractdistribution(d.strpose)
  return DynPose2Pose2(posedistr)
end

#


text:
using Pkg

Pkg.activate(".")

import Accounts: Domain, API

using .Domain, .API

donald_email = create(EMAIL, "donald@duckcity.com")
donald_work = create(WORK,
  """
  Donalds Hardware Store
  attn. Donald Duck
  1190 Seven Seas Dr
  FL 32830 Lake Buena Vista
  USA
  """
)

addresses = [donald_email, donald_work]

donald = create("Donald Duck", addresses)

println(donald)


text:
using Documenter
using BenchmarkTools
using InvertedIndices
using FunctionIndices
using Latexify

makedocs(;
    sitename = "FunctionIndices",
    pagas = ["index.md", "performance.md", "references.md"],
)

deploydocs(; repo = "github.com/wangl-cc/FunctionIndices.jl", push_preview = true)


text:
# condition number

LinearAlgebra.cond(A::LinearOperator, p::Real=2) = cond(coefficients(A), p)

# transpose, adjoint

LinearAlgebra.transpose(A::LinearOperator) = transpose(coefficients(A))

LinearAlgebra.adjoint(A::LinearOperator) = adjoint(coefficients(A))

# factorize

LinearAlgebra.factorize(A::LinearOperator) = factorize(coefficients(A))

# lu

LinearAlgebra.lu(A::LinearOperator; kwargs...) =
    LinearOperator(domain(A), codomain(A), lu(coefficients(A); kwargs...))

LinearAlgebra.lu!(A::LinearOperator; kwargs...) =
    LinearOperator(domain(A), codomain(A), lu!(coefficients(A); kwargs...))

# eigen

LinearAlgebra.eigvals(A::LinearOperator; kwargs...) = eigvals(coefficients(A); kwargs...)

function LinearAlgebra.eigvecs(A::LinearOperator; kwargs...)
    Ξ = eigvecs(coefficients(A); kwargs...)
    Ξ_ = Vector{Sequence{typeof(domain(A)),Vector{eltype(Ξ)}}}(undef, size(Ξ, 2))
    @inbounds for i ∈ axes(Ξ, 2)
        Ξ_[i] = Sequence(domain(A), Ξ[:,i])
    end
    return Ξ_
end

function LinearAlgebra.eigen(A::LinearOperator; kwargs...)
    Λ, Ξ = eigen(coefficients(A); kwargs...)
    Ξ_ = Vector{Sequence{typeof(domain(A)),Vector{eltype(Ξ)}}}(undef, size(Ξ, 2))
    @inbounds for i ∈ axes(Ξ, 2)
        Ξ_[i] = Sequence(domain(A), Ξ[:,i])
    end
    return Λ, Ξ_
end

# kernel

function LinearAlgebra.nullspace(A::LinearOperator; kwargs...)
    Ξ = nullspace(coefficients(A); kwargs...)
    Ξ_ = Vector{Sequence{typeof(domain(A)),Vector{eltype(Ξ)}}}(undef, size(Ξ, 2))
    @inbounds for i ∈ axes(Ξ, 2)
        Ξ_[i] = Sequence(domain(A), Ξ[:,i])
    end
    return Ξ_
end


text:
module ImageDrawTest

using Test
using ImageDraw, ImageCore, ColorVectorSpace

tests = [
    "core.jl",
    "line2d.jl",
    "ellipse2d.jl",
    "circle2d.jl",
    "paths.jl",
    "cross.jl",
    "rectangle.jl",
    "polygon2d.jl"
]

for t in tests
    @testset "$t" begin
        include(t)
    end
end

end


text:
module MathOptInterfaceOSQP

include("modcaches.jl")
using .ModificationCaches

using SparseArrays
using MathOptInterface
using MathOptInterface.Utilities
using LinearAlgebra: rmul!

export Optimizer, OSQPSettings, OSQPModel

const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities
const CI = MOI.ConstraintIndex
const VI = MOI.VariableIndex

const SparseTriplets = Tuple{Vector{Int}, Vector{Int}, Vector{<:Any}}

const Affine = MOI.ScalarAffineFunction{Float64}
const Quadratic = MOI.ScalarQuadraticFunction{Float64}
const VectorAffine = MOI.VectorAffineFunction{Float64}

const Interval = MOI.Interval{Float64}
const LessThan = MOI.LessThan{Float64}
const GreaterThan = MOI.GreaterThan{Float64}
const EqualTo = MOI.EqualTo{Float64}
const IntervalConvertible = Union{Interval, LessThan, GreaterThan, EqualTo}

const Zeros = MOI.Zeros
const Nonnegatives = MOI.Nonnegatives
const Nonpositives = MOI.Nonpositives
const SupportedVectorSets = Union{Zeros, Nonnegatives, Nonpositives}

import OSQP

lower(::Zeros, i::Int) = 0.0
lower(::Nonnegatives, i::Int) = 0.0
lower(::Nonpositives, i::Int) = -Inf
upper(::Zeros, i::Int) = 0.0
upper(::Nonnegatives, i::Int) = Inf
upper(::Nonpositives, i::Int) = 0.0

# TODO: just use ∈ on 0.7 (allocates on 0.6):
function _contains(haystack, needle)
    for x in haystack
        x == needle && return true
    end
    false
end

mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::OSQP.Model
    hasresults::Bool
    results::OSQP.Results
    is_empty::Bool
    silent::Bool
    settings::Dict{Symbol, Any} # need to store these, because they should be preserved if empty! is called
    sense::MOI.OptimizationSense
    objconstant::Float64
    constrconstant::Vector{Float64}
    modcache::ProblemModificationCache{Float64}
    warmstartcache::WarmStartCache{Float64}
    rowranges::Dict{Int, UnitRange{Int}}

    function Optimizer(; kwargs...)
        inner = OSQP.Model()
        hasresults = false
        results = OSQP.Results()
        is_empty = true
        sense = MOI.MIN_SENSE
        objconstant = 0.
        constrconstant = Float64[]
        modcache = ProblemModificationCache{Float64}()
        warmstartcache = WarmStartCache{Float64}()
        rowranges = Dict{Int, UnitRange{Int}}()
        optimizer = new(inner, hasresults, results, is_empty, false,
                        Dict{Symbol, Any}(:verbose => true), sense, objconstant,
                        constrconstant, modcache, warmstartcache, rowranges)
        for (key, value) in kwargs
            MOI.set(optimizer, MOI.RawParameter(key), value)
        end
        return optimizer
    end
end

MOI.get(::Optimizer, ::MOI.SolverName) = "OSQP"

MOI.supports(::Optimizer, ::MOI.Silent) = true
function MOI.set(optimizer::Optimizer, ::MOI.Silent, value::Bool)
    optimizer.silent = value
    if !MOI.is_empty(optimizer)
        if optimizer.silent
            OSQP.update_settings!(optimizer.inner; :verbose => false)
        else
            OSQP.update_settings!(optimizer.inner; :verbose => optimizer.settings[:verbose])
        end
    end
end
MOI.get(optimizer::Optimizer, ::MOI.Silent) = optimizer.silent

MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, limit::Real)
    MOI.set(model, OSQPSettings.TimeLimit(), limit)
    return
end

function MOI.set(model::Optimizer, attr::MOI.TimeLimitSec, ::Nothing)
    delete!(model.settings, :time_limit)
    OSQP.update_settings!(model.inner, time_limit=0.0)
    return
end

function MOI.get(model::Optimizer, ::MOI.TimeLimitSec)
    return get(model.settings, :time_limit, nothing)
end

hasresults(optimizer::Optimizer) = optimizer.hasresults

function MOI.empty!(optimizer::Optimizer)
    optimizer.inner = OSQP.Model()
    optimizer.hasresults = false
    optimizer.results = OSQP.Results()
    optimizer.is_empty = true
    optimizer.sense = MOI.MIN_SENSE # model parameter, so needs to be reset
    optimizer.objconstant = 0.
    optimizer.constrconstant = Float64[]
    optimizer.modcache = ProblemModificationCache{Float64}()
    optimizer.warmstartcache = WarmStartCache{Float64}()
    empty!(optimizer.rowranges)
    optimizer
end

MOI.is_empty(optimizer::Optimizer) = optimizer.is_empty

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike; copy_names=false)
    copy_names && error("Copying names is not supported.")
    MOI.empty!(dest)
    idxmap = MOIU.IndexMap(dest, src)
    assign_constraint_row_ranges!(dest.rowranges, idxmap, src)
    dest.sense, P, q, dest.objconstant = processobjective(src, idxmap)
    A, l, u, dest.constrconstant = processconstraints(src, idxmap, dest.rowranges)
    settings = copy(dest.settings)
    if dest.silent
        settings[:verbose] = false
    end
    OSQP.setup!(dest.inner; P = P, q = q, A = A, l = l, u = u, settings...)
    dest.modcache = ProblemModificationCache(P, q, A, l, u)
    dest.warmstartcache = WarmStartCache{Float64}(size(A, 2), size(A, 1))
    processprimalstart!(dest.warmstartcache.x, src, idxmap)
    processdualstart!(dest.warmstartcache.y, src, idxmap, dest.rowranges)
    dest.is_empty = false
    idxmap
end

"""
Set up index map from `src` variables and constraints to `dest` variables and constraints.
"""
function MOIU.IndexMap(dest::Optimizer, src::MOI.ModelLike)
    idxmap = MOIU.IndexMap()
    vis_src = MOI.get(src, MOI.ListOfVariableIndices())
    for i in eachindex(vis_src)
        idxmap[vis_src[i]] = VI(i)
    end
    i = 0
    for (F, S) in MOI.get(src, MOI.ListOfConstraints())
        MOI.supports_constraint(dest, F, S) || throw(MOI.UnsupportedConstraint{F, S}())
        cis_src = MOI.get(src, MOI.ListOfConstraintIndices{F, S}())
        for ci in cis_src
            i += 1
            idxmap[ci] = CI{F, S}(i)
        end
    end
    idxmap
end

function assign_constraint_row_ranges!(rowranges::Dict{Int, UnitRange{Int}}, idxmap::MOIU.IndexMap, src::MOI.ModelLike)
    startrow = 1
    for (F, S) in MOI.get(src, MOI.ListOfConstraints())
        cis_src = MOI.get(src, MOI.ListOfConstraintIndices{F, S}())
        for ci_src in cis_src
            set = MOI.get(src, MOI.ConstraintSet(), ci_src)
            ci_dest = idxmap[ci_src]
            endrow = startrow + MOI.dimension(set) - 1
            rowranges[ci_dest.value] = startrow : endrow
            startrow = endrow + 1
        end
    end
end

function constraint_rows(rowranges::Dict{Int, UnitRange{Int}}, ci::CI{<:Any, <:MOI.AbstractScalarSet})
    rowrange = rowranges[ci.value]
    length(rowrange) == 1 || error()
    first(rowrange)
end
constraint_rows(rowranges::Dict{Int, UnitRange{Int}}, ci::CI{<:Any, <:MOI.AbstractVectorSet}) = rowranges[ci.value]
constraint_rows(optimizer::Optimizer, ci::CI) = constraint_rows(optimizer.rowranges, ci)

"""
Return objective sense, as well as matrix `P`, vector `q`, and scalar `c` such that objective function is `1/2 x' P x + q' x + c`.
"""
function processobjective(src::MOI.ModelLike, idxmap)
    sense = MOI.get(src, MOI.ObjectiveSense())
    n = MOI.get(src, MOI.NumberOfVariables())
    q = zeros(n)
    if sense != MOI.FEASIBILITY_SENSE
        function_type = MOI.get(src, MOI.ObjectiveFunctionType())
        if function_type == MOI.ScalarAffineFunction{Float64}
            faffine = MOI.get(src, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}())
            P = spzeros(n, n)
            processlinearterms!(q, faffine.terms, idxmap)
            c = faffine.constant
        elseif function_type == MOI.ScalarQuadraticFunction{Float64}
            fquadratic = MOI.get(src, MOI.ObjectiveFunction{MOI.ScalarQuadraticFunction{Float64}}())
            I = [Int(idxmap[term.variable_index_1].value) for term in fquadratic.quadratic_terms]
            J = [Int(idxmap[term.variable_index_2].value) for term in fquadratic.quadratic_terms]
            V = [term.coefficient for term in fquadratic.quadratic_terms]
            symmetrize!(I, J, V)
            P = sparse(I, J, V, n, n)
            processlinearterms!(q, fquadratic.affine_terms, idxmap)
            c = fquadratic.constant
        else
            throw(MOI.UnsupportedAttribute(MOI.ObjectiveFunction{function_type}()))
        end
        sense == MOI.MAX_SENSE && (rmul!(P, -1); rmul!(q, -1); c = -c)
    else
        P = spzeros(n, n)
        q = zeros(n)
        c = 0.
    end
    sense, P, q, c
end

function processlinearterms!(q, terms::Vector{<:MOI.ScalarAffineTerm}, idxmapfun::Function = identity)
    # This is currently needed to avoid depwarns. TODO: make this nice again:
    if q isa VectorModificationCache
        q[:] = 0
    else
        q .= 0
    end
    for term in terms
        var = term.variable_index
        coeff = term.coefficient
        q[idxmapfun(var).value] += coeff
    end
end

function processlinearterms!(q, terms::Vector{<:MOI.ScalarAffineTerm}, idxmap::MOIU.IndexMap)
    processlinearterms!(q, terms, var -> idxmap[var])
end

function symmetrize!(I::Vector{Int}, J::Vector{Int}, V::Vector)
    n = length(V)
    (length(I) == length(J) == n) || error()
    for i = 1 : n
        if I[i] != J[i]
            push!(I, J[i])
            push!(J, I[i])
            push!(V, V[i])
        end
    end
end

function processconstraints(src::MOI.ModelLike, idxmap, rowranges::Dict{Int, UnitRange{Int}})
    m = mapreduce(length, +, values(rowranges), init=0)
    l = Vector{Float64}(undef, m)
    u = Vector{Float64}(undef, m)
    constant = Vector{Float64}(undef, m)
    bounds = (l, u)
    I = Int[]
    J = Int[]
    V = Float64[]
    for (F, S) in MOI.get(src, MOI.ListOfConstraints())
        processconstraints!((I, J, V), bounds, constant, src, idxmap, rowranges, F, S)
    end
    l .-= constant
    u .-= constant
    n = MOI.get(src, MOI.NumberOfVariables())
    A = sparse(I, J, V, m, n)
    (A, l, u, constant)
end

function processconstraints!(triplets::SparseTriplets, bounds::Tuple{<:Vector, <:Vector}, constant::Vector{Float64},
        src::MOI.ModelLike, idxmap, rowranges::Dict{Int, UnitRange{Int}},
        F::Type{<:MOI.AbstractFunction}, S::Type{<:MOI.AbstractSet})
    cis_src = MOI.get(src, MOI.ListOfConstraintIndices{F, S}())
    for ci in cis_src
        s = MOI.get(src, MOI.ConstraintSet(), ci)
        f = MOI.get(src, MOI.ConstraintFunction(), ci)
        rows = constraint_rows(rowranges, idxmap[ci])
        processconstant!(constant, rows, f)
        processlinearpart!(triplets, f, rows, idxmap)
        processconstraintset!(bounds, rows, s)
    end
    nothing
end

function processconstant!(c::Vector{Float64}, row::Int, f::Affine)
    c[row] = MOI.constant(f, Float64)
    nothing
end

function processconstant!(c::Vector{Float64}, rows::UnitRange{Int}, f::VectorAffine)
    for (i, row) in enumerate(rows)
        c[row] = f.constants[i]
    end
end

function processlinearpart!(triplets::SparseTriplets, f::MOI.ScalarAffineFunction, row::Int, idxmap)
    (I, J, V) = triplets
    for term in f.terms
        var = term.variable_index
        coeff = term.coefficient
        col = idxmap[var].value
        push!(I, row)
        push!(J, col)
        push!(V, coeff)
    end
end

function processlinearpart!(triplets::SparseTriplets, f::MOI.VectorAffineFunction, rows::UnitRange{Int}, idxmap)
    (I, J, V) = triplets
    for term in f.terms
        row = rows[term.output_index]
        var = term.scalar_term.variable_index
        coeff = term.scalar_term.coefficient
        col = idxmap[var].value
        push!(I, row)
        push!(J, col)
        push!(V, coeff)
    end
end

function processconstraintset!(bounds::Tuple{<:Vector, <:Vector}, row::Int, s::IntervalConvertible)
    processconstraintset!(bounds, row, MOI.Interval(s))
end

function processconstraintset!(bounds::Tuple{<:Vector, <:Vector}, row::Int, interval::Interval)
    l, u = bounds
    l[row] = interval.lower
    u[row] = interval.upper
    nothing
end

function processconstraintset!(bounds::Tuple{<:Vector, <:Vector}, rows::UnitRange{Int}, s::S) where {S<:SupportedVectorSets}
    l, u = bounds
    for (i, row) in enumerate(rows)
        l[row] = lower(s, i)
        u[row] = upper(s, i)
    end
end

function processprimalstart!(x, src::MOI.ModelLike, idxmap)
    has_primal_start = false
    for attr in MOI.get(src, MOI.ListOfVariableAttributesSet())
        if attr isa MOI.VariablePrimalStart
            has_primal_start = true
        end
    end
    if has_primal_start
        vis_src = MOI.get(src, MOI.ListOfVariableIndices())
        for vi in vis_src
            value = MOI.get(src, MOI.VariablePrimalStart(), vi)
            if value != nothing
                x[idxmap[vi].value] = value
            end
        end
    end
end

function processdualstart!(y, src::MOI.ModelLike, idxmap, rowranges::Dict{Int, UnitRange{Int}})
    for (F, S) in MOI.get(src, MOI.ListOfConstraints())
        has_dual_start = false
        for attr in MOI.get(src, MOI.ListOfConstraintAttributesSet{F, S}())
            if attr isa MOI.ConstraintDualStart
                has_dual_start = true
            end
        end
        if has_dual_start
            cis_src = MOI.get(src, MOI.ListOfConstraintIndices{F, S}())
            for ci in cis_src
                rows = constraint_rows(rowranges, idxmap[ci])
                dual = MOI.get(src, MOI.ConstraintDualStart(), ci)
                if dual != nothing
                    for (i, row) in enumerate(rows)
                        y[row] = -dual[i] # opposite dual convention
                    end
                end
            end
        end
    end
end

## Standard optimizer attributes:
MOI.get(optimizer::Optimizer, ::MOI.ObjectiveSense) = optimizer.sense
function MOI.get(optimizer::Optimizer, a::MOI.NumberOfVariables)
    OSQP.dimensions(optimizer.inner)[1]
end

function MOI.get(optimizer::Optimizer, a::MOI.ListOfVariableIndices)
    [VI(i) for i = 1 : MOI.get(optimizer, MOI.NumberOfVariables())] # TODO: support for UnitRange would be nice
end

## Solver-specific optimizer attributes:
module OSQPSettings

using MathOptInterface
using OSQP

export OSQPAttribute, isupdatable

abstract type OSQPAttribute <: MathOptInterface.AbstractOptimizerAttribute end

# TODO: just use ∈ on 0.7 (allocates on 0.6):
function _contains(haystack, needle)
    for x in haystack
        x == needle && return true
    end
    false
end

for setting in fieldnames(OSQP.Settings)
    Attribute = Symbol(mapreduce(uppercasefirst, *, split(String(setting), '_'))) # to camelcase
    @eval begin
        export $Attribute
        struct $Attribute <: OSQPAttribute end
        Base.Symbol(::$Attribute) = $(QuoteNode(setting))
        isupdatable(::$Attribute) = $(_contains(OSQP.UPDATABLE_SETTINGS, setting))
    end
end
end # module

using .OSQPSettings

_symbol(param::MOI.RawParameter) = Symbol(param.name)
_symbol(a::OSQPAttribute) = Symbol(a)
OSQPSettings.isupdatable(param::MOI.RawParameter) = _contains(OSQP.UPDATABLE_SETTINGS, _symbol(param))
function MOI.set(optimizer::Optimizer, a::Union{OSQPAttribute, MOI.RawParameter}, value)
    (isupdatable(a) || MOI.is_empty(optimizer)) || throw(MOI.SetAttributeNotAllowed(a))
    setting = _symbol(a)
    optimizer.settings[setting] = value
    if !MOI.is_empty(optimizer)
        OSQP.update_settings!(optimizer.inner; setting => value)
    end
end

function MOI.get(optimizer::Optimizer, a::Union{OSQPAttribute, MOI.RawParameter})
    return optimizer.settings[_symbol(a)]
end

## Optimizer methods:
function MOI.optimize!(optimizer::Optimizer)
    processupdates!(optimizer.inner, optimizer.modcache)
    processupdates!(optimizer.inner, optimizer.warmstartcache)
    OSQP.solve!(optimizer.inner, optimizer.results)
    optimizer.hasresults = true
    # Copy previous solution into warm start cache without setting the dirty bit:
    copyto!(optimizer.warmstartcache.x.data, optimizer.results.x)
    copyto!(optimizer.warmstartcache.y.data, optimizer.results.y)
    nothing
end

## Optimizer attributes:
MOI.get(optimizer::Optimizer, ::MOI.RawSolver) = optimizer.inner
MOI.get(optimizer::Optimizer, ::MOI.ResultCount) = optimizer.hasresults ? 1 : 0

MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{Quadratic}) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true

function MOI.set(optimizer::Optimizer, a::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}, obj::MOI.ScalarAffineFunction{Float64})
    MOI.is_empty(optimizer) && throw(MOI.SetAttributeNotAllowed(a))
    optimizer.modcache.P[:] = 0
    processlinearterms!(optimizer.modcache.q, obj.terms)
    optimizer.objconstant = MOI.constant(obj)
    nothing
end

function MOI.set(optimizer::Optimizer, a::MOI.ObjectiveFunction{Quadratic}, obj::Quadratic)
    MOI.is_empty(optimizer) && throw(MOI.SetAttributeNotAllowed(a))
    cache = optimizer.modcache
    cache.P[:] = 0
    for term in obj.quadratic_terms
        row = term.variable_index_1.value
        col = term.variable_index_2.value
        coeff = term.coefficient
        row > col && ((row, col) = (col, row)) # upper triangle only
        cache.P[row, col] += coeff
    end
    processlinearterms!(optimizer.modcache.q, obj.affine_terms)
    optimizer.objconstant = MOI.constant(obj)
    nothing
end

function MOI.get(optimizer::Optimizer, a::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(optimizer, a)
    rawobj = optimizer.results.info.obj_val + optimizer.objconstant
    ifelse(optimizer.sense == MOI.MAX_SENSE, -rawobj, rawobj)
end

error_not_solved() = error("Problem is unsolved.")
function check_has_results(optimizer::Optimizer)
    if !hasresults(optimizer)
        error_not_solved()
    end
end

# Since these aren't explicitly returned by OSQP, I feel like it would be better to have a fallback method compute these:
function MOI.get(optimizer::Optimizer, a::MOI.SolveTime)
    check_has_results(optimizer)
    return optimizer.results.info.run_time
end

function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString)
    return string(optimizer.results.info.status)
end

function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus)
    hasresults(optimizer) || return MOI.OPTIMIZE_NOT_CALLED
    osqpstatus = optimizer.results.info.status
    if osqpstatus == :Unsolved
        return MOI.OPTIMIZE_NOT_CALLED
    elseif osqpstatus == :Interrupted
        return MOI.INTERRUPTED
    elseif osqpstatus == :Dual_infeasible
        return MOI.DUAL_INFEASIBLE
    elseif osqpstatus == :Primal_infeasible
        return MOI.INFEASIBLE
    elseif osqpstatus == :Max_iter_reached
        return MOI.ITERATION_LIMIT
    elseif osqpstatus == :Solved
        return MOI.OPTIMAL
    elseif osqpstatus == :Solved_inaccurate
        return MOI.ALMOST_OPTIMAL
    elseif osqpstatus == :Primal_infeasible_inaccurate
        return MOI.ALMOST_INFEASIBLE
    else
        @assert osqpstatus == :Non_convex
        return MOI.INVALID_MODEL
    end
end

function MOI.get(optimizer::Optimizer, a::MOI.PrimalStatus)
    if a.N > MOI.get(optimizer, MOI.ResultCount())
        return MOI.NO_SOLUTION
    end
    osqpstatus = optimizer.results.info.status
    if osqpstatus == :Unsolved
        return MOI.NO_SOLUTION
    elseif osqpstatus == :Primal_infeasible
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif osqpstatus == :Solved
        return MOI.FEASIBLE_POINT
    elseif osqpstatus == :Primal_infeasible_inaccurate
        return MOI.NEARLY_INFEASIBILITY_CERTIFICATE
    elseif osqpstatus == :Dual_infeasible
        return MOI.INFEASIBILITY_CERTIFICATE
    else # :Interrupted, :Max_iter_reached, :Solved_inaccurate, :Non_convex (TODO: good idea? use OSQP.SOLUTION_PRESENT?)
        return MOI.NO_SOLUTION
    end
end

function MOI.get(optimizer::Optimizer, a::MOI.DualStatus)
    if a.N > MOI.get(optimizer, MOI.ResultCount())
        return MOI.NO_SOLUTION
    end
    osqpstatus = optimizer.results.info.status
    if osqpstatus == :Unsolved
        return MOI.NO_SOLUTION
    elseif osqpstatus == :Dual_infeasible
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif osqpstatus == :Primal_infeasible
        return MOI.INFEASIBILITY_CERTIFICATE
    elseif osqpstatus == :Primal_infeasible_inaccurate
        return MOI.NEARLY_INFEASIBILITY_CERTIFICATE
    elseif osqpstatus == :Solved
        return MOI.FEASIBLE_POINT
    else # :Interrupted, :Max_iter_reached, :Solved_inaccurate, :Non_convex (TODO: good idea? use OSQP.SOLUTION_PRESENT?)
        return MOI.NO_SOLUTION
    end
end

## Variables:
function MOI.is_valid(optimizer::Optimizer, vi::VI)
    vi.value ∈ 1 : MOI.get(optimizer, MOI.NumberOfVariables())
end

## Variable attributes:
function MOI.get(optimizer::Optimizer, a::MOI.VariablePrimal, vi::VI)
    MOI.check_result_index_bounds(optimizer, a)
    x = ifelse(_contains(OSQP.SOLUTION_PRESENT, optimizer.results.info.status), optimizer.results.x, optimizer.results.dual_inf_cert)
    x[vi.value]
end

function MOI.set(optimizer::Optimizer, a::MOI.VariablePrimalStart, vi::VI, value)
    MOI.is_empty(optimizer) && throw(MOI.SetAttributeNotAllowed(a))
    optimizer.warmstartcache.x[vi.value] = value
end

## Constraints:
function MOI.is_valid(optimizer::Optimizer, ci::CI)
    MOI.is_empty(optimizer) && return false
    ci.value ∈ keys(optimizer.rowranges)
end

function MOI.set(optimizer::Optimizer, a::MOI.ConstraintDualStart, ci::CI, value)
    MOI.is_empty(optimizer) && throw(MOI.SetAttributeNotAllowed(a))
    rows = constraint_rows(optimizer, ci)
    for (i, row) in enumerate(rows)
        optimizer.warmstartcache.y[row] = -value[i] # opposite dual convention
    end
    nothing
end

# function modification:
function MOI.set(optimizer::Optimizer, attr::MOI.ConstraintFunction, ci::CI{Affine, <:IntervalConvertible}, f::Affine)
    MOI.is_valid(optimizer, ci) || throw(MOI.InvalidIndex(ci))
    row = constraint_rows(optimizer, ci)
    optimizer.modcache.A[row, :] = 0
    for term in f.terms
        col = term.variable_index.value
        coeff = term.coefficient
        optimizer.modcache.A[row, col] += coeff
    end
    Δconstant = optimizer.constrconstant[row] - f.constant
    optimizer.constrconstant[row] = f.constant
    optimizer.modcache.l[row] += Δconstant
    optimizer.modcache.u[row] += Δconstant
    nothing
end

function MOI.set(optimizer::Optimizer, attr::MOI.ConstraintFunction, ci::CI{VectorAffine, <:SupportedVectorSets}, f::VectorAffine)
    MOI.is_valid(optimizer, ci) || throw(MOI.InvalidIndex(ci))
    rows = constraint_rows(optimizer, ci)
    for row in rows
        optimizer.modcache.A[row, :] = 0
    end
    for term in f.terms
        row = rows[term.output_index]
        col = term.scalar_term.variable_index.value
        coeff = term.scalar_term.coefficient
        optimizer.modcache.A[row, col] += coeff
    end
    for (i, row) in enumerate(rows)
        Δconstant = optimizer.constrconstant[row] - f.constants[i]
        optimizer.constrconstant[row] = f.constants[i]
        optimizer.modcache.l[row] += Δconstant
        optimizer.modcache.u[row] += Δconstant
    end
end

# set modification:
function MOI.set(optimizer::Optimizer, attr::MOI.ConstraintSet, ci::CI{Affine, S}, s::S) where {S <: IntervalConvertible}
    MOI.is_valid(optimizer, ci) || throw(MOI.InvalidIndex(ci))
    interval = S <: Interval ? s : MOI.Interval(s)
    row = constraint_rows(optimizer, ci)
    constant = optimizer.constrconstant[row]
    optimizer.modcache.l[row] = interval.lower - constant
    optimizer.modcache.u[row] = interval.upper - constant
    nothing
end

function MOI.set(optimizer::Optimizer,  attr::MOI.ConstraintSet, ci::CI{VectorAffine, S}, s::S) where {S <: SupportedVectorSets}
    MOI.is_valid(optimizer, ci) || throw(MOI.InvalidIndex(ci))
    rows = constraint_rows(optimizer, ci)
    for (i, row) in enumerate(rows)
        constant = optimizer.constrconstant[row]
        optimizer.modcache.l[row] = lower(s, i) - constant
        optimizer.modcache.u[row] = upper(s, i) - constant
    end
    nothing
end

# partial function modification:
function MOI.modify(optimizer::Optimizer, ci::CI{Affine, <:IntervalConvertible}, change::MOI.ScalarCoefficientChange)
    MOI.is_valid(optimizer, ci) || throw(MOI.InvalidIndex(ci))
    row = constraint_rows(optimizer, ci)
    optimizer.modcache.A[row, change.variable.value] = change.new_coefficient
    nothing
end

# TODO: MultirowChange?

MOI.supports_constraint(optimizer::Optimizer, ::Type{Affine}, ::Type{<:IntervalConvertible}) = true
MOI.supports_constraint(optimizer::Optimizer, ::Type{VectorAffine}, ::Type{<:SupportedVectorSets}) = true

## Constraint attributes:
function MOI.get(optimizer::Optimizer, a::MOI.ConstraintDual, ci::CI)
    MOI.check_result_index_bounds(optimizer, a)
    y = ifelse(_contains(OSQP.SOLUTION_PRESENT, optimizer.results.info.status), optimizer.results.y, optimizer.results.prim_inf_cert)
    rows = constraint_rows(optimizer, ci)
    -y[rows]
end

# Objective modification
function MOI.modify(optimizer::Optimizer, attr::MOI.ObjectiveFunction, change::MOI.ScalarConstantChange)
    MOI.is_empty(optimizer) && throw(MOI.ModifyObjectiveNotAllowed(change))
    optimizer.objconstant = change.new_constant
end

function MOI.modify(optimizer::Optimizer, attr::MOI.ObjectiveFunction, change::MOI.ScalarCoefficientChange)
    MOI.is_empty(optimizer) && throw(MOI.ModifyObjectiveNotAllowed(change))
    optimizer.modcache.q[change.variable.value] = change.new_coefficient
end

# There is currently no ScalarQuadraticCoefficientChange.

MOIU.@model(OSQPModel, # modelname
    (), # scalarsets
    (MOI.Interval, MOI.LessThan, MOI.GreaterThan, MOI.EqualTo), # typedscalarsets
    (MOI.Zeros, MOI.Nonnegatives, MOI.Nonpositives), # vectorsets
    (), # typedvectorsets
    (), # scalarfunctions
    (MOI.ScalarAffineFunction, MOI.ScalarQuadraticFunction), # typedscalarfunctions
    (), # vectorfunctions
    (MOI.VectorAffineFunction,) # typedvectorfunctions
)

end # module


text:
#==========================================================================================#

# CROSS PRODUCTS

crossprod(x::AbstractMatrix)                          = Symmetric(x' * x)
crossprod(x::AbstractMatrix, w::UnitWeights)          = Symmetric(x' * x)
crossprod(x::AbstractMatrix, w::AbstractVector)       = Symmetric(x' * (Diagonal(w) * x))
crossprod(x::AbstractMatrix, w::AbstractSparseMatrix) = Symmetric(x' * Matrix(w * x))
crossprod(x::AbstractMatrix, w::AbstractMatrix)       = Symmetric(x' * (w * x))

# SCALING OF SYMMETRIC MATRICES

function LinearAlgebra.lmul!(s::Number, X::Symmetric)
    n = size(X, 1)
    if X.uplo == "U"
        @inbounds @simd for i in 1:n
            X.data[i, i] *= s
            for j in 1:(i - 1)
                X.data[j, i] *= s
                X.data[i, j] = X.data[j, i]
            end
        end
    else
        @inbounds @simd for i in 1:n
            X.data[i, i] *= s
            for j in (i + 1):n
                X.data[j, i] *= s
                X.data[i, j] = X.data[j, i]
            end
        end
    end
    X
end

#==========================================================================================#

# DISTANCE ON THE GLOBE

havd(x::Float64) = 0.5 * (1.0 - cosd(x))

function geodistance(y0::Float64, x0::Float64, y1::Float64, x1::Float64)
    h = havd(y1 - y0) + havd(x1 - x0) * cosd(y0) * cosd(y1)
    return 12742.0 * asin(min(1.0, sqrt(h)))
end

#==========================================================================================#

# KERNELS

bartlett(x::Real)        = bartlett(float(x))
bartlett(x::Float64)     = (abs(x) >= 1.0) ? 0.0 : (1.0 - abs(x))
truncated(x::Real)       = truncated(float(x))
truncated(x::Float64)    = (abs(x) >= 1.0) ? 0.0 : 1.0
tukeyhanning(x::Real)    = tukeyhanning(float(x))
tukeyhanning(x::Float64) = (abs(x) >= 1.0) ? 0.0 : 0.5 * (1.0 + cospi(x))

parzen(x::Real) = parzen(float(x))

function parzen(x::Float64)
    if abs(x) <= 0.5
        return 1.0 - 6.0 * abs2(x) + 6 * abs(x)^3
    elseif abs(x) <= 1.0
        return 2.0 * (1.0 - abs(x))^3
    else
        return 0.0
    end
end

gallant(x::Real)    = parzen(float(x))
gallant(x::Float64) = parzen(x)

#==========================================================================================#

# FIRST AND SECOND STAGES

first_stage(obj::TwoStageModel)  = obj.first_stage
second_stage(obj::TwoStageModel) = obj.second_stage

switch_stage(obj::OneStageModel) = obj
switch_stage(obj::TwoStageModel) = obj.second_stage
switch_stage(obj::ParEstimate)   = obj

#==========================================================================================#

# RETRIEVAL

getvector(obj::ParModel, x::Symbol) = getvector(switch_stage(obj).sample, x)
getmatrix(obj::ParModel, args...)   = getmatrix(switch_stage(obj).sample, args...)
getnames(obj::ParModel, args...)    = getnames(switch_stage(obj).sample, args...)
getcorr(obj::ParModel)              = getcorr(switch_stage(obj).sample)
getnonmissing(obj::AnyModel)        = getnonmissing(switch_stage(obj).sample)
getweights(obj::AnyModel)           = getweights(switch_stage(obj).sample)

#==========================================================================================#

# REWEIGHTING

reweight(w::UnitWeights, v::ProbabilityWeights)     = v
reweight(w::AbstractWeights, v::ProbabilityWeights) = ProbabilityWeights(w .* v)

#==========================================================================================#

# FORMATTER

function formatter(X::Array{<: Real}, d::Int)

    table = format.("{:.$(d)f}", X)

    for c = 1:size(table, 2)
        l = maximum(length.(table[:, c]))
        table[:, c] .= format.("{:>$(l)s}", table[:, c])
    end

    return table
end

max_stars_repo_name: lbittarello/Microeconometrics

max_issues_repo_name: lbittarello/Microeconometrics

max_forks_repo_name: lbittarello/Microeconometrics


text:
module AmbientForcing

    using OrdinaryDiffEq: solve, ODEFunction, ODEProblem, Tsit5
    using ForwardDiff: jacobian
    using LinearAlgebra: nullspace, I
    include("AmbientForcingODE.jl")
    include("RandInitVectors.jl")
    export ambient_forcing
    export random_force
    export idx_exclusive
    export constraint_equations
end