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

using Markdown
using InteractiveUtils

# ╔═╡ 254d20da-8b47-448b-bb16-c078c86d872f
md"""
### Computing the Metabolic Demand of Transcription and Translation

In this lecture, we'll continue our discussion of transcription and translation processes. In particular, we'll discuss:

* Cell-free simulation of synthetic transcription (TX) and translation (TL) circuits
* The analysis of TX/TL systems using flux balance analysis
* Coupling TX/TL processes with metabolism in a cell-free system
"""

# ╔═╡ 8617bf13-ca94-434d-b6ea-9d982fcd2012
md"""
### Simulation of cell free synthetic circuits

One strategy to model (and analyze) gene expression programs is to build ordinary differential equations (ODEs) models and to solve these models numerically.

Let's take a look at an example that ignores metabolism:

* [Adhikari A, Vilkhovoy M, Vadhin S, Lim HE, Varner JD. Effective Biophysical Modeling of Cell-Free Transcription and Translation Processes. Front Bioeng Biotechnol. 2020;8:539081. Published 2020 Nov 26. doi:10.3389/fbioe.2020.539081](https://www.ncbi.nlm.nih.gov/labs/pmc/articles/PMC7726328/)
"""

# ╔═╡ 8aaf3475-6ea7-45ea-90f5-0c4673925258
md"""
### Accounting for the cost of transcription (TX) and translation (TL)

Nature has developed extensive control mechanisms which regulate when transcription and translation processes occur (accounted for by the $u\left(\dots\right)$ and $w\left(\dots\right)$ functions) because of their metabolic cost.

Let's take a look at using flux balance analysis to compute the metabolic demand associated with gene expression:

* [Allen TE, Palsson BØ. Sequence-based analysis of metabolic demands for protein synthesis in prokaryotes. J Theor Biol. 2003 Jan 7;220(1):1-18. doi: 10.1006/jtbi.2003.3087. PMID: 12453446.](https://pubmed.ncbi.nlm.nih.gov/12453446/)
"""

# ╔═╡ c835ff26-9df0-416c-a64f-9e04db57543b
md"""
##### Analysis of green fluorescent protein (GFP) expression

The standard flux balance analysis problem is written in `concentration` units and metabolic reaction flux, e.g., $\star$mol/gDW-hr. However, there is nothing that says we have to do that. For example, working in mass or mole units may be more convenient than working in concentration units. Let's consider the latter.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

# ╔═╡ d052ff3c-2a8d-4897-b3ce-aa08d06f8124
md"""
###### Setup: objective coefficient array
"""

# ╔═╡ d4c47ae1-d670-48f6-a2f1-61c8ff5620db
md"""
###### Setup: species bounds array
"""

# ╔═╡ 213041c0-2bfe-4c10-81da-5301b4d46b58
md"""
###### Setup: flux bounds array
"""

# ╔═╡ 981fcf59-63f4-4312-b5e6-cb74a76e5d55
begin

	# setup the parameters -
	parameters = Dict{String, Any}()

	# let's compute a literature value -
	k₁ = 100.0 			# units: 1/conc-t
	eₓ = 35.0 			# units: nt/s
	eₗ = 2.0 			# units: aa/s
	Lₓ = 922.0 			# units: nt
	Lₗ = 238.0 			# units: aa
	Kₗ = 483.0 			# units: μmol/L
	Kₓ = 0.036 			# units: μmol/L
	V = 30.0*(1/1e6) 	# units: L

	# get parameters from bionumbers -
	parameters["volume"] = V									# units: L
	parameters["RNAP_concentration"] = 0.06 					# units: μmol/L
	parameters["RNAP_elongation_rate"] = eₓ 					# units: nt/s BNID 111871
	parameters["ribsome_elongation_rate"] = eₗ 					# units: aa/s
	parameters["gene_coding_length"] = Lₓ 						# units: nt BNID 111610
	parameters["mRNA_half_life"] = 2.5*(60)						# units: s BIND 104324
	parameters["gene_concentration"] = 8.0*(V)*(1e6/1e9) 		# units: μmol/L
	parameters["saturation_constant_transcription"] = Kₓ 	# units: μmol/L
	parameters["ribosome_concentration"] = 2.3 				# units: μmol/L
	parameters["saturation_constant_translation"] = Kₗ 		# units: μmol/L

	# show -
	nothing
end

# ╔═╡ 02756d50-5881-49c8-be5a-d228b90b3912
md"""
###### Solve
"""

# ╔═╡ f12fdff4-0ce5-4d07-a278-16d4f41ede2f
md"""
__Table__: Open extent table $\dot{\epsilon}_{j}$ (units: $\star$mol/time) for the transcription and translation of green fluorescent protein (GFP).
"""

# ╔═╡ 6f852067-8c2e-45d0-baa9-73875ca49e12
md"""
__Table__: Input/output mol transfer rate $\dot{n}_{i,\star}$ (units: ⋆mol/time) for the transcription and translation of green fluorescent protein (GFP).
"""

# ╔═╡ 15d4c02f-cacd-4293-b6bc-ba4b76981caa
md"""
### Integrating TX/TL with metabolism

The Allen and Palsson study provides an approach to analyze transcription and translation processes using flux balance analysis. However, the metabolic demands are treated as exchanges.

Let's integrate the TX/TL box proposed by Allen and Palsson with a metabolic network:

* [Vilkhovoy M, Horvath N, Shih CH, Wayman JA, Calhoun K, Swartz J, Varner JD. Sequence-Specific Modeling of E. coli Cell-Free Protein Synthesis. ACS Synth Biol. 2018 Aug 17;7(8):1844-1857. doi: 10.1021/acssynbio.7b00465. Epub 2018 Jul 16. PMID: 29944340.](https://pubmed.ncbi.nlm.nih.gov/29944340/)
"""

# ╔═╡ 89958d90-2a62-4306-bb52-01dd4e921c7a
md"""
### Summary and conclusions

In this lecture we:

* Introduced a continuous model of transcription (TX) and translation (TL) of a cell-free synthetic circuit
* Analyzed TX/TL processes using flux balance analysis (introduced the TX/TL box)
* Coupled TX/TL box with metabolism in a cell-free system producing a model protein

"""

# ╔═╡ 17f52bcc-ee1b-4e14-bb8b-7603e67164f4
md"""
### Next time:

* The end of the beginning. Putting everything together to solve a metabolic engineering problem
"""

# ╔═╡ 34d2f753-d74a-457a-83bc-44e08d6323ca
function ingredients(path::String)

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

# ╔═╡ 656bc222-a43f-11ec-1f09-6b332c3b672f
begin

	# load some packages -
	using PlutoUI
	using BSON
	using DataFrames
	using TOML
	using GLPK
	using PrettyTables

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

	# load the 5440/7770 code lib -
    lib = ingredients(joinpath(_PATH_TO_SRC, "Include.jl"))

	# show -
	nothing
end

# ╔═╡ c7321b9e-b741-4d2d-8fd8-f3b8954d368f
begin

	# load the GFP model file as a dictionary -
	model_file_path = joinpath(_PATH_TO_DATA,"GFP.bson")
	model = BSON.load(model_file_path, @__MODULE__)[:system]

	# show -
	nothing
end

# ╔═╡ ffaa05e2-2407-4407-8583-e46e36ffe399
model

# ╔═╡ e97e290a-d920-4797-8315-c65adf94c8d5
model["reaction_table"]

# ╔═╡ 18137921-25a7-4572-ba8c-39203031c24c
S = model["stoichiometric_matrix"]

# ╔═╡ 75394d6e-f185-4d3c-8b27-4d27c1a7c94a
(ℳ,ℛ) = size(S)

# ╔═╡ 0562d302-30dd-4c88-a9bb-5ec4100b3060
begin
	# setup the objective coefficient array -
	c_vector = model["objective_coefficient_array"]
	reaction_table = model["reaction_table"]

	# find the index of the reaction we want to maximize -
	target_reaction_id = "mRNA_EGFP_translation"
	idx_target = findfirst(x->x==target_reaction_id, reaction_table[!,:id])

	# update -
	if (isnothing(c_vector) == true)
		c_vector = zeros(ℛ)
	end
	c_vector[idx_target] = -1.0
	model["objective_coefficient_array"] = c_vector

	# show -
	nothing
end

# ╔═╡ 29a92990-1e95-49cd-aac0-7b68e1099a1d
model["species_table"][!,:symbol]

# ╔═╡ e4538ec6-678d-43dd-a01e-18a92dfe8e18
begin

	# get the default flux bounds array -
	flux_bounds_array = model["flux_bounds_array"]

	# get some stuff out of parameters -
	G = parameters["gene_concentration"]
	RNAP = parameters["RNAP_concentration"]
	RIBOSOME = parameters["ribosome_concentration"]
	Vₜ = parameters["volume"]
	t_half_life = parameters["mRNA_half_life"]

	# setup constraints on gene expression -
	Vmax_TX = RNAP*(eₓ)*(1/Lₓ)
	r_TX = Vmax_TX*(G/(Kₓ + G))
	flux_bounds_array[1,1] = (r_TX*Vₜ)
	flux_bounds_array[2,1] = (r_TX*Vₜ)
	flux_bounds_array[3,1] = (r_TX*Vₜ)

	# show -
	nothing
end

# ╔═╡ 758c8e5e-7a41-40c6-945c-bf68e93947e9
aa_map = TOML.parsefile(joinpath(_PATH_TO_DATA,"AAMap.toml"))

# ╔═╡ a1f07b91-e911-4811-b02b-64ddc487e6af
begin

	# species_bounds_array
	species_bounds_array = model["species_bounds_array"]

	# update the upper bound -
	species_bounds_array[:,2] .= Inf # everything can go out of the box at whatever rate -

	# update the lower bound -
	species_bounds_array[5:9,1] .= -1000.0

	# we need to update the bounds of all of the AA's -
	species_symbol_array = model["species_table"][!,:symbol]
	for (key, aa_symbol) ∈ aa_map

		# find index of aa_symbol -
		aa_index = findfirst(x->x==aa_symbol, species_symbol_array)

		# update the bound -
		species_bounds_array[aa_index,1] = -100.0
	end

	# show -
	nothing
end

# ╔═╡ fbd77bc8-9388-4e79-9036-bdca3f03fd4c
begin

	# estimate the optimal flux -
	result = lib.flux(S,flux_bounds_array,species_bounds_array,c_vector);

	# show nothing -
	nothing
end

# ╔═╡ 2d05bbd7-772f-4d63-b647-4d4fcb4a9cde
with_terminal() do
	println("Exit flag = $(result.exit_flag) and status flag = $(result.status_flag)")
end

# ╔═╡ ea038677-c71d-4e64-992d-c36497fe52f8
with_terminal() do

	v = result.calculated_flux_array

	# let's build a table -
	state_array = Array{Any,2}(undef, ℛ, 5)
	for i ∈ 1:ℛ
		state_array[i,1] = i
		state_array[i,2] = v[i]*(1e9/1e6)*(3600)
		state_array[i,3] = reaction_table[i,:id]
		state_array[i,4] = reaction_table[i,:forward_reaction]
		state_array[i,5] = reaction_table[i,:reverse_reaction]
	end

	# setup the header -
	header_data = (
		["i", "ϵ", "reaction name", "foward", "reverse"],
		["", "nmol/hr", "", "", ""]
	);

	pretty_table(state_array; alignment=:l, header=header_data)
end

# ╔═╡ ada9b431-b449-40c8-9a65-77a88259615e
with_terminal() do

	# get the "flux" -
	ϵ_dot = result.calculated_flux_array

	# get compute the
	n_dot_in = abs.(species_bounds_array[:,1])
	n_dot_out = n_dot_in + S*ϵ_dot

	# setup species table -
	state_array = Array{Any, 2}(undef, ℳ, 5)
	for i ∈ 1:ℳ
		state_array[i,1] = i
		state_array[i,2] = species_symbol_array[i]
		state_array[i,3] = round(n_dot_in[i]; digits=2)
		state_array[i,4] = round(n_dot_out[i]; digits=2)
		state_array[i,5] = round((n_dot_out[i] - n_dot_in[i]); digits=2)
	end

	# setup the header -
	header_data = (
		["i", "species", "nᵢ_in", "nᵢ_out", "Δn"],
		["", "", "⋆mol/time", "⋆mol/time", "⋆mol/time"]
	);

	pretty_table(state_array; header=header_data)
end

# ╔═╡ 446ef6c4-b3c5-4b17-9cec-098dac93774b
html"""
<style>
main {
    max-width: 900px;
    width: 70%;
    margin: auto;
    font-family: "Roboto, monospace";
}

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

# ╔═╡ cb4f47b8-e63a-4638-a832-38a1b38a329c
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]
BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[compat]
BSON = "~0.3.5"
DataFrames = "~1.3.2"
GLPK = "~1.0.0"
PlutoUI = "~0.7.37"
PrettyTables = "~1.3.1"
"""

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

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.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.BSON]]
git-tree-sha1 = "306bb5574b0c1c56d7e1207581516c557d105cad"
uuid = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
version = "0.3.5"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[deps.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.Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

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

[[deps.GLPK]]
deps = ["GLPK_jll", "MathOptInterface"]
git-tree-sha1 = "7971e2ce3715a873b539174137bd8c4e19ac7a8f"
uuid = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
version = "1.0.0"

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

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

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

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

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

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

[[deps.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.MathOptInterface]]
deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "JSON", "LinearAlgebra", "MutableArithmetics", "OrderedCollections", "Printf", "SparseArrays", "Test", "Unicode"]
git-tree-sha1 = "a62df301482a41cb7b1db095a4e6949ba7eb3349"
uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
version = "1.1.0"

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

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

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

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

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

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

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

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

[[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.PooledArrays]]
deps = ["DataAPI", "Future"]
git-tree-sha1 = "db3a23166af8aebf4db5ef87ac5b00d36eb771e2"
uuid = "2dfb63ee-cc39-5dd5-95bd-886bf059d720"
version = "1.4.0"

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

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

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

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

[[deps.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.SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

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

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

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

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

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

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

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

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

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

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

[[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:
# ╟─254d20da-8b47-448b-bb16-c078c86d872f
# ╟─8617bf13-ca94-434d-b6ea-9d982fcd2012
# ╟─8aaf3475-6ea7-45ea-90f5-0c4673925258
# ╟─c835ff26-9df0-416c-a64f-9e04db57543b
# ╠═c7321b9e-b741-4d2d-8fd8-f3b8954d368f
# ╠═ffaa05e2-2407-4407-8583-e46e36ffe399
# ╠═e97e290a-d920-4797-8315-c65adf94c8d5
# ╠═18137921-25a7-4572-ba8c-39203031c24c
# ╠═75394d6e-f185-4d3c-8b27-4d27c1a7c94a
# ╟─d052ff3c-2a8d-4897-b3ce-aa08d06f8124
# ╠═0562d302-30dd-4c88-a9bb-5ec4100b3060
# ╟─d4c47ae1-d670-48f6-a2f1-61c8ff5620db
# ╠═29a92990-1e95-49cd-aac0-7b68e1099a1d
# ╠═758c8e5e-7a41-40c6-945c-bf68e93947e9
# ╠═a1f07b91-e911-4811-b02b-64ddc487e6af
# ╟─213041c0-2bfe-4c10-81da-5301b4d46b58
# ╠═981fcf59-63f4-4312-b5e6-cb74a76e5d55
# ╠═e4538ec6-678d-43dd-a01e-18a92dfe8e18
# ╟─02756d50-5881-49c8-be5a-d228b90b3912
# ╠═fbd77bc8-9388-4e79-9036-bdca3f03fd4c
# ╟─2d05bbd7-772f-4d63-b647-4d4fcb4a9cde
# ╟─f12fdff4-0ce5-4d07-a278-16d4f41ede2f
# ╠═ea038677-c71d-4e64-992d-c36497fe52f8
# ╟─6f852067-8c2e-45d0-baa9-73875ca49e12
# ╠═ada9b431-b449-40c8-9a65-77a88259615e
# ╟─15d4c02f-cacd-4293-b6bc-ba4b76981caa
# ╟─89958d90-2a62-4306-bb52-01dd4e921c7a
# ╟─17f52bcc-ee1b-4e14-bb8b-7603e67164f4
# ╠═656bc222-a43f-11ec-1f09-6b332c3b672f
# ╠═34d2f753-d74a-457a-83bc-44e08d6323ca
# ╠═446ef6c4-b3c5-4b17-9cec-098dac93774b
# ╠═cb4f47b8-e63a-4638-a832-38a1b38a329c
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using PurlinLineDesign
using StructuresKit

#INPUT LAYER

mutable struct purlin_line_object

    design_code::String
    loading_direction::String
    span_geometry::Vector{Tuple{Float64, Int64, Int64}}
    spacing::Float64
    roof_slope::Float64
    cross_section_dimensions::Vector{Tuple{String, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, Float64}}
    material_properties::Vector{NTuple{4, Float64}}
    deck_details::Tuple{String, Float64, Float64, Float64, Float64, Float64}
    deck_material_properties::NTuple{4, Float64}
    frame_flange_width::Float64
    support_locations::Vector{Float64}
    bridging_locations::Vector{Float64}

    purlin_line_object() = new()

end

design_code = "AISI S100-16 ASD"

loading_direction = "gravity"   #or "uplift"

# Define the properties of each purlin segment along the line.
                #length, section_properties, material_properties
span_geometry = [(23.0*12, 1, 1),
                (2.0*12, 2, 1),
                (2.0*12, 2, 1),
                (23.0*12, 1, 1)]

# Define the purlin spacing.
spacing = 60;  #in.

# Define the roof slope.
roof_slope = 0.0;   #degrees

#Define the purlin cross-section type.
cross_section_dimensions = [("Z", 0.059, 0.91, 2.5, 8.0, 2.5, 0.91, -50.0, 0.0, 90.0, 0.0, -50.0, 3*0.059, 3*0.059, 3*0.059, 3*0.059),
                                        ("Z", 0.059, 0.91, 2.5, 8.0, 2.5, 0.91, -50.0, 0.0, 90.0, 0.0, -50.0, 3*0.059, 3*0.059, 3*0.059, 3*0.059)]

material_properties = [(29500.0, 0.30, 55.0, 70.0),
                              (29500.0, 0.30, 55.0, 70.0)]

#type="screw-fastened", thickness, fastener spacing, fastener diameter, fastener_shear_strength fastener_top_flange_location
#type="standing seam", thickness, clip spacing, clip stifness
deck_details = ("screw-fastened", 0.0179, 12.0, 0.212, 2.50, 1.25)

deck_material_properties = (29500.0, 0.30, 55.0, 70.0)

frame_flange_width = 24.0

support_locations = [0.0, 25.0*12, 50.0*12]

bridging_locations =[0.0, 50.0*12]

#MAPPING LAYER

#Create the data structure.
purlin_line = purlin_line_object()

#Map user inputs to the data structure.
purlin_line.design_code = design_code
purlin_line.loading_direction = loading_direction
purlin_line.span_geometry = span_geometry
purlin_line.spacing = spacing
purlin_line.roof_slope = roof_slope
purlin_line.cross_section_dimensions = cross_section_dimensions
purlin_line.material_properties = material_properties
purlin_line.deck_details = deck_details
purlin_line.deck_material_properties = deck_material_properties
purlin_line.frame_flange_width = frame_flange_width
purlin_line.support_locations = support_locations
purlin_line.bridging_locations = bridging_locations

#DEFINITIONS LAYER

##Calculate properties associated with each purlin line segment.

#Calculate the purlin section properties.
purlin_section_properties, purlin_free_flange_section_properties = PurlinLineDesign.calculate_purlin_section_properties(purlin_cross_section_dimensions)

plot(purlin_free_flange_section_properties[1].node_geometry[:, 1],purlin_free_flange_section_properties[1].node_geometry[:,2], seriestype = scatter)

#local buckling
#calculate along with cross-section

#bracing stiffnesses
#distortional buckling
#flexural strength
#shear strength
#web crippling strength

#ANALYSIS LAYER

#RESULTS LAYER

#stiffness provided by deck to purlin
#This needs section properties and material properties.
#Calculate for every unique segment along the span.

num_purlin_sections = size(purlin_cross_section_dimensions)[1]

#Go through the purlin cross-sections to calculate lateral and rotational stiffness.

#Define the spacing between top flange connections that restrain distortional buckling.
Lm = deck_purlin_fastener_spacing

function define_purlin_line_stiffness(purlin_line)

if deck_details[1] == "screw-fastened"

    #Define the deck to purlin screw-fastened connection spacing.
    deck_purlin_fastener_spacing = deck_details[3]

    #Define the deck to purlin screw diameter.
    deck_purlin_fastener_diameter = deck_details[4]

    #Define the nominal shear strength of the typical screw.
    Fss = deck_details[5]

    #Define the roof deck base metal thickness.
    t_roof_deck = deck_details[2]

    #Define roof deck steel elastic modulus.
    E_roof_deck = deck_material_properties[1]

    #Define roof deck steel ultimate yield stress.
    Fu_roof_deck = deck_material_properties[4]

    #Define purlin steel elastic modulus.
    E_purlin = purlin_material_properties[1][1]

    #Define purlin steel Poisson's ratio.
    μ_purlin = purlin_material_properties[1][2]

    #Define purlin steel ultimate yield stress.
    Fu_purlin =

    for i = 1:num_purlin_sections

        #Define the purlin top flange width.
        b_top = purlin_cross_section_dimensions[i][6]

        #Define purlin base metal thickness.
        t_purlin = purlin_cross_section_dimensions[i][2]

        #Define out-to-out purlin web depth.
        ho = purlin_cross_section_dimensions[i][5]

        #Define purlin top flange lip length.
        d_top = purlin_cross_section_dimensions[i][7]

        #Define purlin top flange lip angle from the horizon, in degrees.
        θ_top = purlin_cross_section_dimensions[i][11] - purlin_cross_section_dimensions[i][12]

        #Define the location from the purlin top flange pivot point to the fastener.  Assume the fastener is centered in the flange.
        c = b_top/2

        #Define the distance from the purlin web to the screw location in the top flange.  Assume it is centered in the top flange.
        deck_purlin_fastener_top_flange_location = b_top/2

        #Define the deck fastener pull-through plate stiffness.  Assume the fastener is centered between two panel ribs.
        kp = deck_pull_through_fastener_stiffness(deck_material_properties, b_top, t_roof_deck)

        #Calculate the rotational stiffness provided to the purlin by the screw-fastened connection between the deck and the purlin.  It is assumed that the deck flexural stiffness is much higher than the connection stiffness.
        kϕ = Connections.cfs_rot_screwfastened_k(b_top, c, deck_purlin_fastener_spacing, t_purlin, kp, E_purlin, CorZ)

        #Calculate the purlin distortional buckling half-wavelength.

        if deck_details[1] == "Z"
            CorZ = 1
        elseif deck_details[1] == "C"
            CorZ = 0
        end

        #Calculate top flange + lip section properties.
        Af, Jf, Ixf, Iyf, Ixyf, Cwf, xof,  hxf, hyf, yof = AISIS10016.table23131(CorZ, t_purlin, b_top, d_top, θ_top)

        #Calculate the purlin distortional buckling half-wavelength.
        Lcrd, L = AISIS10016.app23334(ho, μ_purlin, t_purlin, Ixf, xof, hxf, Cwf, Ixyf, Iyf, Lm)

        #If Lcrd is longer than the fastener spacing, then the distortional buckling will be restrained by the deck.
        if Lcrd >= Lm
            kϕ_dist = kϕ
        else
            kϕ_dist = 0.0
        end

        #Approximate the lateral stiffness provided to the top of the purlin by the screw-fastened connection between the deck and the purlin.

        #Calculate the stiffness of a single screw-fastened connection.
        Ka, ψ, α, β, Ke = Connections.cfs_trans_screwfastened_k(t_roof_deck, t_purlin, E_roof_deck, E_purlin, Fss, Fu_roof_deck, Fu_purlin, deck_purlin_fastener_diameter)

        #Convert the discrete stiffness to a distributed stiffness, divide by the fastener spacing.
        kx = Ke / deck_purlin_fastener_spacing

ext: jl

max_stars_repo_name: runtosolve/PurlinLine

max_issues_repo_name: runtosolve/PurlinLine

max_forks_repo_name: runtosolve/PurlinLine


text:
@deprecate backward(model) reverse_differentiate!(model) false
@deprecate forward(model) forward_differentiate!(model) false

@deprecate BackwardInVariablePrimal() ReverseVariablePrimal() false
@deprecate BackwardOutConstraint() ReverseConstraintFunction() false
@deprecate BackwardOutObjective() ReverseObjectiveFunction() false

@deprecate ForwardOutVariablePrimal() ForwardVariablePrimal() false
@deprecate ForwardInConstraint() ForwardConstraintFunction() false
@deprecate ForwardInObjective() ForwardObjectiveFunction() false

@deprecate QPForwBackCache(args...) QuadraticForwardReverseCache(args...) false
@deprecate ConicBackCache(args...) ConicReverseCache(args...) false

@deprecate QPDiff(args...) QuadraticDiffProblem(args...) false
@deprecate ConicDiff(args...) ConicDiffProblem(args...) false
@deprecate QPForm(args...) QuadraticProblemForm(args...) false


text:
abstract type AbstractModel end

include("objectives/obj.jl")
include("maximin_model.jl")
include("soyster_maximin_model.jl")
include("de_jong_maximin_model.jl")
include("robust_maximin_model.jl")
include("minimax_model.jl")
include("custom_model.jl")
include("ccmaximin_model.jl")
include("no_obj_model.jl")


text:
include("base/base.jl")
include("spaces/spaces.jl")

include("action_selectors/action_selectors.jl")
include("approximators/approximators.jl")
include("buffers/buffers.jl")
include("environments/environments.jl")
include("environment_models/environment_models.jl")
include("policies/policies.jl")
include("learners/learners.jl")

include("agents/agents.jl")


text:
# global mutable directory connecting benchmark names with dataframes
typealias Capture @compat Tuple{Symbol, Vector{Symbol}}

# Lensname -> Varname -> Vector of captued values
clear_captured!() = global captured = Dict{Symbol,Dict{Symbol,Vector{Any}}}()
clear_captured!()

# Convert a dict into an other one with only those keys in ks
extract{K,V}(d::Dict{K,V},ks::Vector{K}) = Dict{K,V}([k=>d[k] for k in ks])

# Capture the data and add it to global 'captured'
function capturebench!(captures::Vector{Symbol}, data::Data)
  global captured
  # From the Data object we just pull out the data values restricted
  # to the vars we want to capture
  lensname = data.lensname
  extracteddata = extract(data.data, captures)

  # Default dict behaviour, if no list then create it, otherwise append
  if haskey(captured,lensname)
    for (varname,value) in extracteddata
      if haskey(captured[lensname],varname)
        push!(captured[lensname][varname],value)
      else
        captured[lensname][varname] = Any[]
        push!(captured[lensname][varname],value)
      end
    end
  else
    captured[lensname] = [k => [v] for (k,v) in extracteddata]
  end
end

# Creates a listener for each capture and register to
# The associated data to be captured
function register_captured!(captures::Vector{Capture})
  for c in captures
    let c = c
      λ = data -> capturebench!(c[2],data)
      register!(c[1], Listener(:benchmark, λ, true, true))
    end
  end
end

# Register lenses
function setup!{C<:Capture}(captures::Vector{C})
  clear_captured!()
  register_captured!(captures)
end

# Unregister lenses and delete benchmark data
function cleanup!()
  delete_listener!(:benchmark)
  clear_captured!()
end

@doc "Quick and dirty capture:
  Evaluates `f()` and captures any values specified in `captures`.
  Also returns result of `f()`" ->
function capture{C<:Capture}(f::Function, captures::Vector{C}; exceptions = true)
  for proc in procs()
    fetch(@spawnat proc setup!(captures))
  end

  local Δt
  local Δb
  local value
  if exceptions
    try
      value, Δt, Δb = @timed(f())
    catch e
      println("Uncaught exception trickled down to capture:")
      println(e)
      value = nothing
      Δt = 0
      Δb = 0
    end
  else
    value, Δt, Δb = @timed(f())
  end
  lens(:total_time, Δt)
  local res

  # When there are multiple processors, collate all data
  res = Result()
  for proc in procs()
    res.values[proc] = remotecall_fetch(proc, ()->Lens.captured)
  end
  for proc in procs() @spawnat proc cleanup!() end
  value,res
end

# Hack for failture of type inference to detect [:a, (:a,b)] as Capture vec
capture(f::Function, captures::Vector{Any}; exceptions = true) =
  capture(f,Capture[captures...]; exceptions = exceptions)
# Convenience - if we just use a lens, assume we want the first captured var
capture(f::Function, capturename::Symbol; exceptions = true) =
  capture(f, [(capturename, [:x1])]; exceptions = exceptions)
capture(f::Function, captures::Vector{Symbol}; exceptions = true) =
  capture(f, [(capture, [:x1]) for capture in captures]; exceptions = exceptions)

macro capture(expr,captures)
   :(capture(()->$expr,$captures))
end


text:
using Base.Test

include("../src/CutShortURL.jl")

@test CutShortURL.cutshorturl("https://rahulkp220.github.io") == "http://tinyurl.com/ydgk73kx"


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Arpack"
version = v"3.5.0-3"

# Collection of sources required to build ArpackMKLBuilder
sources = [
    "https://github.com/opencollab/arpack-ng.git" =>
    "b095052372aa95d4281a645ee1e367c28255c947",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd arpack*
mkdir build
cd build
# cmake -DINTERFACE64=1 -DBUILD_SHARED_LIBS=ON -DLAPACK_FOUND=true -DLAPACK_INCLUDE_DIRS=${MKLROOT}/include -DLAPACK_LIBRARIES=mkl_rt -DBLAS_FOUND=true -DBLAS_INCLUDE_DIRS=${MKLROOT}/include -DBLAS_LIBRARIES=mkl_rt ..
CMAKE_EXE_LINKER_FLAGS="${CMAKE_EXE_LINKER_FLAGS} -Wl,-rpath -Wl,$MKLROOT/lib"
# Symbols that have float32, float64, complexf32, and complexf64 support
SDCZ_SYMBOLS="axpy copy gemv geqr2 lacpy lahqr lanhs larnv lartg \
              lascl laset scal trevc trmm trsen gbmv gbtrf gbtrs \
              gttrf gttrs pttrf pttrs"
# All symbols that have float32/float64 support (including the SDCZ_SYMBOLS above)
SD_SYMBOLS="${SDCZ_SYMBOLS} dot ger labad laev2 lamch lanst lanv2 \
            lapy2 larf larfg lasr nrm2 orm2r rot steqr swap"
# All symbols that have complexf32/complexf64 support (including the SDCZ_SYMBOLS above)
CZ_SYMBOLS="${SDCZ_SYMBOLS} dotc geru unm2r"
# Add in (s|d)*_64 symbol remappings:
for sym in ${SD_SYMBOLS}; do     SYMBOL_DEFS="${SYMBOL_DEFS} -Ds${sym}=s${sym}_64 -Dd${sym}=d${sym}_64"; done
# Add in (c|z)*_64 symbol remappings:
for sym in ${CZ_SYMBOLS}; do     SYMBOL_DEFS="${SYMBOL_DEFS} -Dc${sym}=c${sym}_64 -Dz${sym}=z${sym}_64"; done
# Add one-off symbol mappings; things that don't fit into any other bucket:
for sym in scnrm2 dznrm2 csscal zdscal dgetrf dgetrs; do     SYMBOL_DEFS="${SYMBOL_DEFS} -D${sym}=${sym}_64"; done
# Set up not only lowercase symbol remappings, but uppercase as well:
SYMBOL_DEFS="${SYMBOL_DEFS} ${SYMBOL_DEFS^^}"
FFLAGS="${FFLAGS} -fdefault-integer-8 ${SYMBOL_DEFS} -ff2c"
cmake ..  -DINTERFACE64=1 -DBUILD_SHARED_LIBS=ON -DBLAS_LIBRARIES="-L$MKLROOT/lib -lmkl_rt" -DLAPACK_LIBRARIES="-L$MKLROOT/lib -lmkl_rt" -DCMAKE_EXE_LINKER_FLAGS="${CMAKE_EXE_LINKER_FLAGS}" -DCMAKE_Fortran_FLAGS="${FFLAGS}"
make clean all
make install

# Install into ${prefix}/lib
mkdir -p ${prefix}/lib
cp lib/libarpack.so.2.0.0 ${prefix}/lib
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, libc=:glibc),
    Linux(:i686, libc=:glibc)
]
platforms = expand_gcc_versions(platforms)

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libarpack", :libarpack)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://raw.githubusercontent.com/JuliaComputing/MKL.jl/76ab0ac2b671eb4e199c1b4a31872ccb22898e4f/deps/build_MKL.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: giordano/ArpackMKLBuilder

max_issues_repo_name: giordano/ArpackMKLBuilder

max_forks_repo_name: giordano/ArpackMKLBuilder


text:
include("setup.jl")

@testset "core solver" begin
    @testset "basic example" begin
        packages = [fill("A", 2); fill("B", 4); fill("C", 4)]
        conflicts = [(1, 3), (1, 7), (3, 7), (2, 10)]
        solutions = [[1, 4, 8], [2, 3, 8], [2, 4, 7]]
        @test solutions == resolve_brute_force(packages, conflicts)
        @test solutions == resolve_brute_force(packages, Set(conflicts))
        @test solutions == resolve_core(packages, conflicts)
        @test solutions == resolve_core(packages, Set(conflicts))
        for relax = 1:5
            solutions′ = resolve_brute_force(packages, conflicts; relax)
            resolved′ = resolve_core(packages, conflicts; relax)
            @test solutions′ == resolved′
        end
        @testset "permutations" for _ = 1:100
            # permute the versions, conflicts and solutions
            packages′ = shuffle(packages)
            p = sortperm(packages′)
            conflicts′ = [(p[c[1]], p[c[2]]) for c in conflicts]
            solutions′ = [p[solution] for solution in solutions]
            foreach(sort!, solutions′)
            sort!(solutions′)
            @test solutions′ == resolve_brute_force(packages′, conflicts′)
            @test solutions′ == resolve_core(packages′, conflicts′)
        end
        @testset "unordered package names" begin
            packages = [fill(1+0im, 2); fill(2+0im, 4); fill(3+0im, 4)]
            @test solutions == resolve_brute_force(packages, conflicts)
            @test solutions == resolve_core(packages, conflicts)
        end
    end

    @testset "comprehensive tests" begin
        Random.seed!(0x53f3ce0656b85450bbb52b15fc58853f)
        count = zeros(Int, 17)
        for M = 2:5, # number of packages
            V = 1:5  # number of versions per package
            T = V^2*(M*(M-1)÷2)
            T ≤ 128 || continue
            packages = shuffle!([p for p=1:M for _=1:V])
            p = sortperm(packages)
            for C in (T ≤ 12 ? (0:2^T-1) : [randu128(3-(k%5)%3) for k=1:2^12])
                G = gen_conflicts(M, V, C)
                conflicts = Tuple{Int,Int}[(p[c[1]], p[c[2]]) for c in G]
                @assert length(conflicts) == count_ones(C % UInt128(2)^T)
                for relax = 0:M^2
                    solutions = resolve_brute_force(packages, conflicts; relax)
                    resolved = resolve_core(packages, conflicts; relax)
                    @test resolved isa Vector{Vector{Int}}
                    @test solutions == resolved
                    count[length(solutions)+1] += 1
                end
            end
        end
        # specific to chosen seed (chosen to have some large solution sets)
        @test count == [
            12693, 605486, 17731, 11964, 5792, 2073,
            777, 341, 163, 65, 23, 11, 9, 1, 0, 0, 1,
        ]
    end

    @testset "edge cases" begin
        @testset "0 pkgs" begin
            packages = []
            conflicts = Tuple{Int,Int}[]
            @test resolve_core(packages, conflicts) == [[]]
        end
        @testset "1 pkgs" begin
            for V = 1:3
                packages = fill("Only", V)
                conflicts = Tuple{Int,Int}[]
                @test resolve_core(packages, conflicts) == [[1]]
            end
        end
    end
end

@testset "resolver API" begin
    @testset "example: 0 pkgs" begin
        # also tests that pacakge and version types flow through
        for (P, V) in [(String, VersionNumber), (Real, Complex{Int})]
            versions = Dict{P,Vector{V}}()
            required = P[]
            compat = gen_compat()
            resolved = Resolver.resolve(compat, versions, required)
            @test resolved == [[]]
            @test eltype(resolved) == Vector{Pair{P,V}}
        end
    end

    @testset "example: 1 pkgs" begin
        versions = Dict("A" => 1:2)
        required = ["A"]
        compat = gen_compat()
        resolved = Resolver.resolve(compat, versions, required)
        @test resolved == [["A" => 1]]
    end

    @testset "example: 1 pkgs (VersionNumbers)" begin
        versions = Dict("A" => [v"2", v"1"])
        required = ["A"]
        compat = gen_compat()
        resolved = Resolver.resolve(compat, versions, required)
        @test resolved == [["A" => v"2"]]
    end

    @testset "example: 2 pkgs, 0 conflicts" begin
        versions = Dict(
            "A" => 1:2,
            "B" => 1:2,
        )
        required = ["A"]
        compat = gen_compat()
        resolved = Resolver.resolve(compat, versions, required)
        @test resolved == [["A" => 1]]
    end

    @testset "example: 2 pkgs, 1 conflicts" begin
        versions = Dict(
            "A" => 1:2,
            "B" => 1:2,
        )
        conflicts = Dict(
            ("A", "B") => [(1, 1)]
        )
        required = ["A"]
        compat = gen_compat(conflicts)
        resolved = Resolver.resolve(compat, versions, required)
        @test resolved == [["A" => 1]]
    end

    @testset "example: 2 pkgs, 2 deps" begin
        versions = Dict(
            "A" => 1:2,
            "B" => 1:2,
        )
        deps = Dict(
            ("A" => 1) => ["B"],
            ("A" => 2) => ["B"],
        )
        required = ["A"]
        compat = gen_compat(deps)
        resolved = Resolver.resolve(compat, versions, required)
        @test resolved == [["A" => 1, "B" => 1]]
    end

    @testset "example: 2 pkgs, 2 deps, 1 conflicts" begin
        versions = Dict(
            "A" => 1:2,
            "B" => 1:2,
        )
        deps = Dict(
            ("A" => 1) => ["B"],
            ("A" => 2) => ["B"],
        )
        conflicts = Dict(
            ("A", "B") => [(1, 1)]
        )
        required = ["A"]
        compat = gen_compat(deps, conflicts)
        resolved = Resolver.resolve(compat, versions, required)
        @test resolved == [
            ["A" => 1, "B" => 2],
            ["A" => 2, "B" => 1],
        ]
    end

    @testset "example: 2 pkgs, 2 conflicts, 2 reqs" begin
        # package type must be sortable
        # version type can be anything
        versions = Dict(
            "A" => 1:2,
            "B" => 1:2
        )
        conflicts = Dict(
            ("A", "B") => [(1, 1)],
            ("B", "A") => [(2, 2)],
        )
        required = ["B", "A"]
        compat = gen_compat(conflicts)
        resolved = Resolver.resolve(compat, versions, required)
        @test resolved == [
            ["A" => 1, "B" => 2],
            ["A" => 2, "B" => 1],
        ]
    end

    @testset "example: 2 pkgs, 2 conflicts, 2 reqs (weird types)" begin
        # package type must be sortable
        # version type can be anything
        versions = Dict(
            :A => [1+0im, 2+0im],
            :B => [1+0im, 2+0im]
        )
        conflicts = Dict(
            (:A, :B) => [(1+0im, 1+0im)],
            (:B, :A) => [(2+0im, 2+0im)],
        )
        required = [:B, :A]
        compat = gen_compat(conflicts)
        resolved = Resolver.resolve(compat, versions, required)
        @test resolved == [
            [:A => 1+0im, :B => 2+0im],
            [:A => 2+0im, :B => 1+0im],
        ]
    end

    @testset "example: 2 pkgs, 3 deps, 5 conflicts" begin
        versions = Dict(
            "A" => 1:4,
            "B" => 1:3,
        )
        deps = Dict(
            ("A" => 1) => ["B"],
            ("A" => 2) => ["B"],
            ("B" => 3) => ["A"],
        )
        conflicts = Dict(
            ("A", "B") => [
                (1, 1)
                (2, 1)
                (3, 1)
                (4, 1)
                (4, 3)
            ]
        )
        required = ["A"]
        compat = gen_compat(deps, conflicts)
        resolved = Resolver.resolve(compat, versions, required)
        @test resolved == [
            ["A" => 1, "B" => 2],
            ["A" => 3],
        ]
    end

    @testset "example: 4 pkgs, 2 deps, 2 conflicts" begin
        versions = Dict(
            "A" => 1:2,
            "B" => 1:2,
            "C" => 1:1,
            "D" => 1:1,
        )
        deps = Dict(
            ("A" => 1) => ["B", "C", "D"],
            ("A" => 2) => ["B", "C", "D"],
        )
        conflicts = Dict(
            ("A", "C") => [(1, 1)],
            ("B", "D") => [(1, 1)],
        )
        required = ["A"]
        compat = gen_compat(deps, conflicts)
        resolved = Resolver.resolve(compat, versions, required)
        @test resolved == [
            ["A" => 2, "B" => 2, "C" => 1, "D" => 1],
        ]
    end

    @testset "example: 3 pkgs, 2 deps, 2 conflicts (incomplete)" begin
        versions = Dict(
            "A" => 1:1,
            "B" => 1:1,
            "C" => 1:2,
        )
        deps = Dict(
            ("A" => 1) => ["C"],
            ("B" => 1) => ["C"],
        )
        conflicts = Dict(
            ("A", "C") => [(1, 1)],
            ("B", "C") => [(1, 2)],
        )
        required = ["A", "B"]
        compat = gen_compat(deps, conflicts)
        resolved = Resolver.resolve(compat, versions, required)
        @test resolved == [
            ["A" => 1, "C" => 2],
            ["B" => 1, "C" => 1],
        ]
    end

    @testset "example: 3 pkgs, 2 deps, 1 conflicts (incompatible)" begin
        versions = Dict(
            "A" => 1:1,
            "B" => 1:1,
            "C" => 1:1,
        )
        deps = Dict(
            ("A" => 1) => ["B"],
            ("B" => 1) => ["C"],
        )
        conflicts = Dict(
            ("A", "C") => [(1, 1)],
        )
        required = ["A"]
        compat = gen_compat(deps, conflicts)
        resolved = Resolver.resolve(compat, versions, required)
        @test resolved == [["A" => 1, "B" => 1, "C" => 1]]
    end

    @testset "example: 3 pkgs, 2 deps, 1 conflicts (incompatible)" begin
        versions = Dict(
            "A" => 1:1,
            "B" => 1:1,
            "C" => 1:1,
        )
        deps = Dict(
            ("A" => 1) => ["B", "C"],
        )
        conflicts = Dict(
            ("B", "C") => [(1, 1)],
        )
        required = ["A"]
        compat = gen_compat(deps, conflicts)
        resolved = Resolver.resolve(compat, versions, required)
        @test resolved == [["A" => 1, "B" => 1, "C" => 1]]
    end

    @testset "example: 3 pkgs, 2 deps, 3 conflicts (incompatible)" begin
        versions = Dict(
            "A" => 1:1,
            "B" => 1:1,
            "C" => 1:1,
        )
        deps = Dict(
            ("A" => 1) => ["B", "C"],
        )
        conflicts = Dict(
            ("A", "B") => [(1, 1)],
            ("A", "C") => [(1, 1)],
            ("B", "C") => [(1, 1)],
        )
        required = ["A"]
        compat = gen_compat(deps, conflicts)
        resolved = Resolver.resolve(compat, versions, required)
        @test resolved == [["A" => 1, "B" => 1, "C" => 1]]
    end

    @testset "example: 3 pkgs, 2 deps, 2 conflicts (incompatible)" begin
        versions = Dict(
            "A" => 1:2,
            "B" => 1:1,
            "C" => 1:1,
        )
        deps = Dict(
            ("A" => 1) => ["B"],
            ("A" => 2) => ["C"],
        )
        conflicts = Dict(
            ("A", "B") => [(1, 1)],
            ("A", "C") => [(2, 1)],
        )
        required = ["A"]
        compat = gen_compat(deps, conflicts)
        resolved = Resolver.resolve(compat, versions, required)
        @test resolved == [
            ["A" => 1, "B" => 1],
            ["A" => 2, "C" => 1],
        ]
    end
end


text:
@swagger """
/:
  get:
    description: Testing swagger markdown root get!
    responses:
      '200':
        description: Returns a mysterious string get.
"""

@swagger """
/:
  put:
    description: Testing swagger markdown root put!
    responses:
      '200':
        description: Returns a mysterious string put.
"""

@swagger """
/:
  post:
    description: Testing swagger markdown root post!
    responses:
      '200':
        description: Returns a mysterious string post.
"""

@swagger """
/test:
  post:
    description: Testing swagger markdown test!
    responses:
      '200':
        description: Returns a mysterious string test.
"""

@testset "Without swagger doc defined" begin

end

@testset "Parse markdown for routes correctly" begin
    paths = SwaggerMarkdown.parse_spec()
    cleanup()

    @testset "'/' route should have GET, PUT and POST" begin
        ROOT_ROUTE = "/"
        @test haskey(paths, ROOT_ROUTE)
        r = paths[ROOT_ROUTE]
        @test haskey(r, GET)
        @test haskey(r, PUT)
        @test haskey(r, POST)

        for i in [GET, POST, PUT]
            @assert haskey(r[i], "description")
            @test r[i]["description"] == "Testing swagger markdown root $(i)!"
            @assert haskey(r[i], "responses")
            @assert haskey(r[i]["responses"], "200")
            @assert haskey(r[i]["responses"]["200"], "description")
            @test r[i]["responses"]["200"]["description"] == "Returns a mysterious string $(i)."
        end
    end

    @testset "'/test' route should have POST" begin
        TEST = "/test"
        @test haskey(paths, TEST)
        r = paths[TEST]
        @test haskey(r, POST)
        @assert haskey(r[POST], "description")
        @test r[POST]["description"] == "Testing swagger markdown test!"
        @assert haskey(r[POST], "responses")
        @assert haskey(r[POST]["responses"], "200")
        @assert haskey(r[POST]["responses"]["200"], "description")
        @test r[POST]["responses"]["200"]["description"] == "Returns a mysterious string test."
    end
end


text:
#
# x, y: Arrays of dimensions (nx,3), (ny,3)
# ifmol, ilmol: first and last indexes of array x to be considered (first and last atoms of molecule)
# jfmol, jlmol: first and last indexes of array y to be considered
# returns the the minimum distance between the points considered and
# the indexes of these points in x and y vectors
#

# If x is only a vector
"""

Returns the the minimum distance between the points considered.

"""
function minimumdistance(x::T, y::Vector{T}) where {T}
    jatom = 0
    dmin = +Inf
    for (j, yj) in pairs(y)
        d = distance(x, yj)
        if d < dmin
            jatom = j
            dmin = d
        end
    end
    return dmin, 1, jatom
end

# If both are vectors of vectors

function minimumdistance(x::Vector{T}, y::Vector{T}) where {T}
    iatom = 0
    jatom = 0
    dmin = +Inf
    for (i, xi) in pairs(x)
        for (j, yj) in pairs(y)
            d = distance(xi, yj)
            if d < dmin
                iatom = i
                jatom = j
                dmin = d
            end
        end
    end
    return dmin, iatom, jatom
end

# Function that returns the distance of a reference atom as well, to be used for
# computation of the volume shell by Monte-Caro integration

function minimumdistance(x::Vector{T}, y::Vector{T}, jrefatom::Int) where {T}
    iatom = 0
    jatom = 0
    drefatom = +Inf
    dmin = +Inf
    for (i, xi) in pairs(x)
        for (j, yj) in pairs(y)
            d = distance(xi, yj)
            # Minimum distance of any solvent atom to the solute
            if d < dmin
                iatom = i
                jatom = j
                dmin = d
            end
            # Minimum distance of the reference atom to the solute
            if j == jrefatom
                if d < drefatom
                    drefatom = d
                end
            end
        end
    end
    return dmin, iatom, jatom, drefatom
end

#
# With periodic boundary conditions
#

function minimumdistance(x::Vector{T}, y::Vector{T}, sides::T) where {T}
    iatom = 0
    jatom = 0
    dmin = +Inf
    for (i, xi) in pairs(x)
        for (j, yj) in pairs(y)
            d = distance(xi, yj, sides)
            if d < dmin
                iatom = i
                jatom = j
                dmin = d
            end
        end
    end
    return dmin, iatom, jatom
end

# Function that returns the distance of a reference atom as well, to be used for
# computation of the volume shell by Monte-Caro integration

function minimumdistance(x::Vector{T}, y::Vector{T}, jrefatom::Int, sides::T) where {T}
    iatom = 0
    jatom = 0
    drefatom = +Inf
    dmin = +Inf
    for (i, xi) in pairs(x)
        for (j, yj) in pairs(y)
            d = distance(xi, yj, sides)
            # Minimum distance of any solvent atom to the solute
            if d < dmin
                iatom = i
                jatom = j
                dmin = d
            end
            # Minimum distance of the reference atom to the solute
            if j == jrefatom
                if d < drefatom
                    drefatom = d
                end
            end
        end
    end
    return dmin, iatom, jatom, drefatom
end


text:
module NaturalSort
export natural

function natural(x::AbstractString, y::AbstractString)
    statex = start(x)
    statey = start(y)
    donex = false
    doney = false
    while !(donex = done(x, statex)) & !(doney = done(y, statey))
        cx, statex = next(x, statex)
        cy, statey = next(y, statey)
        if isnumber(cx) && isnumber(cy)
            # Skip leading zeros
            while cx == '0' && !(donex = done(x, statex))
                cx, statex = next(x, statex)
            end
            while cy == '0' && !(doney = done(y, statey))
                cy, statey = next(y, statey)
            end

            # Begin comparing numbers
            diff = false
            lt = false
            while true
                isnumx = isnumber(cx)
                isnumy = isnumber(cy)
                if isnumx && isnumy
                    if !diff && cx != cy
                        # Keep track of how numbers differ, in case the lengths match
                        diff = true
                        lt = cx < cy
                    end

                    donex = done(x, statex)
                    doney = done(y, statey)
                    if donex || doney
                        if donex && !doney && isnumber(next(y, statey)[1])
                            # Number in y is longer than number in x
                            return true
                        elseif !donex && doney && isnumber(next(x, statex)[1])
                            # Number in x is longer than number in y
                            return false
                        end
                        # Both numbers ended and same length
                        return diff ? lt : donex && !doney
                    end

                    cx, statex = next(x, statex)
                    cy, statey = next(y, statey)
                elseif isnumx
                    # Number in x is longer than number in y
                    return false
                elseif isnumy
                    # Number in y is longer than number in x
                    return true
                elseif diff
                    # Numbers were same length but different
                    return lt
                else
                    # Numbers were the same
                    break
                end
            end
        end
        if cx != cy
            return cx < cy
        end
    end

    return donex && !doney
end
end # module


text:
"""
    AbstractEddyViscosityClosure <: AbstractTurbulenceClosure

Abstract supertype for turbulence closures that are defined by an isotropic viscosity
and isotropic diffusivities.
"""
abstract type AbstractEddyViscosityClosure <: AbstractTurbulenceClosure end

#####
##### Viscosity-stress products
#####

"""
    ν_σᶜᶜᶜ(i, j, k, grid, ν, σᶜᶜᶜ, u, v, w)

Multiply the array `ν` located at `ᶜᶜᶜ` by a function

    `σᶜᶜᶜ(i, j, k, grid, u, v, w)`

at index `i, j, k` and location `ᶜᶜᶜ`.
"""
@inline ν_σᶜᶜᶜ(i, j, k, grid, ν, σᶜᶜᶜ, u, v, w) =
    @inbounds ν[i, j, k] * σᶜᶜᶜ(i, j, k, grid, u, v, w)

"""
    ν_σᶠᶠᶜ(i, j, k, grid, ν, σᶠᶠᶜ, u, v, w)

Multiply the array `ν` located at `ᶜᶜᶜ` by a function

    `σᶠᶠᶜ(i, j, k, grid, u, v, w)`

at index `i, j, k` and location `ᶠᶠᶜ`.
"""
@inline ν_σᶠᶠᶜ(i, j, k, grid, ν, σᶠᶠᶜ, u, v, w) =
    @inbounds ℑxyᶠᶠᵃ(i, j, k, grid, ν) * σᶠᶠᶜ(i, j, k, grid, u, v, w)

# These functions are analogous to the two above, but for different locations:
@inline ν_σᶠᶜᶠ(i, j, k, grid, ν, σᶠᶜᶠ, u, v, w) =
    @inbounds ℑxzᶠᵃᶠ(i, j, k, grid, ν) * σᶠᶜᶠ(i, j, k, grid, u, v, w)

@inline ν_σᶜᶠᶠ(i, j, k, grid, ν, σᶜᶠᶠ, u, v, w) =
    @inbounds ℑyzᵃᶠᶠ(i, j, k, grid, ν) * σᶜᶠᶠ(i, j, k, grid, u, v, w)

#####
##### Stress divergences
#####

viscous_flux_ux(i, j, k, grid, clock, closure::AbstractEddyViscosityClosure, U, diffusivities) = 2 * ν_σᶜᶜᶜ(i, j, k, grid, diffusivities.νₑ, Σ₁₁, U.u, U.v, U.w)
viscous_flux_uy(i, j, k, grid, clock, closure::AbstractEddyViscosityClosure, U, diffusivities) = 2 * ν_σᶠᶠᶜ(i, j, k, grid, diffusivities.νₑ, Σ₁₂, U.u, U.v, U.w)
viscous_flux_uz(i, j, k, grid, clock, closure::AbstractEddyViscosityClosure, U, diffusivities) = 2 * ν_σᶠᶜᶠ(i, j, k, grid, diffusivities.νₑ, Σ₁₃, U.u, U.v, U.w)

viscous_flux_vx(i, j, k, grid, clock, closure::AbstractEddyViscosityClosure, U, diffusivities) = 2 * ν_σᶠᶠᶜ(i, j, k, grid, diffusivities.νₑ, Σ₂₁, U.u, U.v, U.w)
viscous_flux_vy(i, j, k, grid, clock, closure::AbstractEddyViscosityClosure, U, diffusivities) = 2 * ν_σᶜᶜᶜ(i, j, k, grid, diffusivities.νₑ, Σ₂₂, U.u, U.v, U.w)
viscous_flux_vz(i, j, k, grid, clock, closure::AbstractEddyViscosityClosure, U, diffusivities) = 2 * ν_σᶜᶠᶠ(i, j, k, grid, diffusivities.νₑ, Σ₂₃, U.u, U.v, U.w)

viscous_flux_wx(i, j, k, grid, clock, closure::AbstractEddyViscosityClosure, U, diffusivities) = 2 * ν_σᶠᶜᶠ(i, j, k, grid, diffusivities.νₑ, Σ₃₁, U.u, U.v, U.w)
viscous_flux_wy(i, j, k, grid, clock, closure::AbstractEddyViscosityClosure, U, diffusivities) = 2 * ν_σᶜᶠᶠ(i, j, k, grid, diffusivities.νₑ, Σ₃₂, U.u, U.v, U.w)
viscous_flux_wz(i, j, k, grid, clock, closure::AbstractEddyViscosityClosure, U, diffusivities) = 2 * ν_σᶜᶜᶜ(i, j, k, grid, diffusivities.νₑ, Σ₃₃, U.u, U.v, U.w)


text:
# Morphological open
""" For dim != 2, this function still creates a 2d circle.
"""
function kernelcircle(r, dim=2)
    n = 2 * ceil(Int, r) + 1
    A = centered(ones(Int32, n, n, ones(Int, dim-2)...))
    r2 = r^2
    for I in CartesianIndices(A)
        x, y = Tuple(I)
        if x^2 + y^2 > r2
            A[I] = 0
        end
    end
    A
end

function dilation(A::AbstractArray{T, D}, r, alg=ImageFiltering.Algorithm.FIR()) where {D, T}
    k = kernelcircle(r, D)
    imfilter(A, k, alg) .>= 0.5
end
function erosion(A::AbstractArray{T, D}, r, alg=ImageFiltering.Algorithm.FIR()) where {D, T}
    k = kernelcircle(r, D)
    imfilter(A, k, alg) .>= sum(k) - 0.5
end
function closing(A, r)
    A1 = dilation(A, r)
    erosion(A1, r)
end
function opening(A::AbstractArray{T, D}, r, alg) where {D, T}
    k = kernelcircle(r, D)
    A1 = imfilter(A, k, alg) .>= sum(k) - 0.5
    imfilter(A1, k, alg) .>= 0.5
end
function opening(A, r)
    opening(A, r, ImageFiltering.Algorithm.FIR())
end

function internalgradient(A::AbstractArray{T, D}) where {D, T}
    Bool.(A - erosion(A, 1))
end
function externalgradient(A::AbstractArray{T, D}) where {D, T}
    Bool.(dilation(A, 1) - A)
end
function gradient(A::AbstractArray{T, D}) where {D, T}
    Bool.(dilation(A, 1) - erosion(A, 1))
end


text:
using Distributions

const DISTRIBUTIONS = [Degenerate, Gaussian, Bernoulli]

# The role of these tests is to ensure that the
# the general functionality always works and
# doesn't result in NaN or Inf weights.
# Not that the RBM performance is still good.
if :integration in TEST_GROUPS
    @testset "RBM Integration" begin
        n_vis = 100
        n_hid = 10

        for T in [Float32, Float64]
            rbm = RBM(T, Gaussian, Bernoulli, n_vis, n_hid)
            Boltzmann.test(rbm; debug=true)
        end
    end
end

if :benchmark in TEST_GROUPS
    @testset "RBM Benchmark" begin
        n_vis = 100
        n_hid = 10

        results = DataFrame()
        for T in [Float32, Float64]
            rbm = RBM(T, Gaussian, Bernoulli, n_vis, n_hid)
            df = Boltzmann.benchmark(rbm; debug=true)
            df[:Type] = fill(T, size(df, 1))
            results = vcat(results, df)
        end

        # For now just print the output, but we may want load a CSV containing
        # existing master benchmarks and compare the two here.
        println(results)
    end
end


text:
module Resolver

const SetOrVector{T} = Union{AbstractSet{T}, AbstractVector{T}}

function resolve(
    compat   :: Function, # (p₁ => v₁, p₂ => v₂) -> Bool
    versions :: AbstractDict{P, <:AbstractVector{V}},
    required :: SetOrVector{P},
) where {P,V}
    vertices, conflicts = vertices_and_conflicts(
        compat,
        versions,
        required,
    )
    packages = P[p for (p, v) in vertices]
    dummies = Bool[isnothing(v) for (p, v) in vertices]
    let resolved
        M = length(versions)
        for relax = 0:M*(M-1)÷2
            solutions = resolve_core(packages, conflicts; dummies, relax)
            # sort solutions
            for S in solutions
                sort!(S, by = v -> packages[v])
            end
            sort!(solutions)
            # map back to version identifiers
            resolved = Vector{Pair{P,V}}[
                [vertices[v] for v in S if !dummies[v]]
                    for S in solutions
            ]
            # only keep the most satisfying solutions
            sat = [sum(p in required for (p, v) in S; init=0) for S in resolved]
            resolved = resolved[sat .≥ maximum(sat; init = 1-isempty(versions))]
            !isempty(resolved) && break
        end
        resolved
    end
end

function vertices_and_conflicts(
    compat   :: Function, # (p₁ => v₁, p₂ => v₂) -> Bool
    versions :: AbstractDict{P, <:AbstractVector{V}},
    required :: SetOrVector{P},
) where {P,V}
    # check compat callback function signature
    hasmethod(compat, Tuple{Pair{P,V}, Pair{P,V}}) &&
    hasmethod(compat, Tuple{Pair{P,V}, Pair{P,Nothing}}) ||
        throw(ArgumentError("compat: callback takes two package-version pairs"))

    # check package versions data structure
    let seen = Set{V}()
        for (p, vers) in versions
            isempty(vers) &&
                throw(ArgumentError("versions: package $p: no versions"))
            for v in vers
                v in seen &&
                    throw(ArgumentError("versions: package $p: duplicate version $v"))
                push!(seen, v)
            end
            empty!(seen)
        end
    end

    # check required packages set
    for p in required
        p in keys(versions) ||
            throw(ArgumentError("required: package $p: not in versions"))
    end

    # co-compute reachable versions and conflicts between them
    package_names = sort!(collect(keys(versions)))
    reachable = Pair{P,Int}[p => Int(p in required) for p in package_names]
    conflicts = Set{NTuple{2,Int}}()
    conflicted = Set(required)
    while true
        clean = true
        for (i₁, (p₁, k₁)) in enumerate(reachable),
            (i₂, (p₂, k₂)) in enumerate(reachable)
            p₁ < p₂ || continue
            v₁ = get(versions[p₁], k₁, nothing)
            v₂ = get(versions[p₂], k₂, nothing)
            (v₁ === nothing || compat(p₁ => v₁, p₂ => v₂)) &&
            (v₂ === nothing || compat(p₂ => v₂, p₁ => v₁)) && continue
            push!(conflicts, minmax(i₁, i₂))
            for (p, k) in (p₁ => k₁ + 1, p₂ => k₂ + 1)
                if k ≤ length(versions[p]) + (p ∈ required) && (p => k) ∉ reachable
                    push!(reachable, p => k)
                    clean = false
                end
                push!(conflicted, p)
            end
        end
        clean && break
    end
    vertices = Pair{P,Union{V,Nothing}}[
        p => get(versions[p], k, nothing)
            for (p, k) in reachable if p in conflicted
    ]

    return vertices, conflicts
end

function resolve_core(
    packages  :: AbstractVector,
    conflicts :: SetOrVector{<:NTuple{2,Integer}};
    dummies   :: Vector{Bool} = fill(false, length(packages)),
    relax     :: Integer = 0,
)
    # check conflicts
    for (v₁, v₂) in conflicts
        v₁ in keys(packages) ||
            throw(ArgumentError("conflicts: invalid version index: $v₁"))
        v₂ in keys(packages) ||
            throw(ArgumentError("conflicts: invalid version index: $v₂"))
        v₁ != v₂ ||
            throw(ArgumentError("conflicts: package $(packages[v₁]): self-conflict $v₁"))
        packages[v₁] != packages[v₂] ||
            throw(ArgumentError("conflicts: package $(packages[v₁]): conflict between $v₁, $v₂"))
    end

    # conflict count matrix
    N = length(packages)
    X = zeros(UInt32, N, N)
    for (v₁, p₁) in enumerate(packages),
        (v₂, p₂) in enumerate(packages)
        x = (v₁, v₂) ∈ conflicts || (v₂, v₁) ∈ conflicts
        d = x && (dummies[v₁] || dummies[v₂])
        X[v₁, v₂] = X[v₂, v₁] = p₁ == p₂ ? relax + 2 : d ? relax + 1 : x
    end

    # deduplicate nodes by conflict sets
    keep = UInt32[]
    seen = Set{Vector{UInt32}}()
    for (v, p) in enumerate(packages)
        X[:, v] in seen && continue
        push!(seen, X[:, v])
        push!(keep, v)
    end
    X = X[keep, keep]
    packages = packages[keep]

    # package indices
    P = zeros(UInt32, length(packages))
    let indices = Dict{eltype(packages), UInt32}()
        for (v, p) in enumerate(packages)
            P[v] = get!(indices, p, length(indices) + 1)
        end
    end

    # find all optimal solutions
    solutions = optimal_solutions(P, X, Int(relax))

    # translate back to original version indices
    Vector{Int}[Int[keep[v] for v in S] for S in solutions]
end

# The core solver is based on naive k-clique listing algorithm in Chiba &
# Nishizeki 1985, "Arboricity and subgraph listing algorithms". It is also
# similar to the well-known Bron-Kerbosch maximal clique listing algorithm, but
# doesn't require keeping track of the set of excluded nodes because the clique
# size is known in advance. This version lists only optimal solutions in the
# following sense: we have a total ordering on the versions of each package and
# we must choose some version of each package (possibly the special "no version"
# version); this induces a partial ordering on solutions and the solutions that
# are listed are those that are optimal by this partial ordering. These are the
# conflict-free solutions that a user might want to choose based on version
# preferences. If there is a conflict between the latest versions of packages A
# and B, they might want the latest version of A and an older version of B or
# the latest version of B and an older version of A. Whereas the total number of
# feasible conflict-free solutions is exponential in the number of conflict
# edges, the number of optimal solutions tends to be fairly small.
#
# The pivoting optimization that can be used to speed up Bron-Kerbosch cannot be
# used here because if we both limit vertices considered at each recursion level
# to the neighborhood of some node and require choosing a node that dominates a
# specific solution in the set of already-found solutions, then we miss valid,
# optimal solutions. However, the domination requirement acts as a form of
# pivoting anyway: after finding an optimal solution, we only consider versions
# that improve on that solution for some package. Those versions must also
# conflict with that solution since otherwise it wouldn't be optimal. So rather
# than pivoting around a node, we are pivoting around an entire solution.

function optimal_solutions(P::Vector{T}, X::Matrix{T}, c::Int = 0) where {T<:Integer}
    # number of packages & versions
    M = maximum(P, init=0)
    L = maximum(X, init=1) + M
    N = size(X, 1)

    # conflict count vector, level vector, solution vector, solutions set
    C = zeros(T, N)
    S = zeros(T, M)
    solutions = typeof(S)[]
    M == 0 && return push!(solutions, S)

    function search!(c::Int, r::Int = 1, d::Int = 1)
        for j = 1:N
            # check subgraph inclusion at this level
            c′ = c - C[j]
            c′ ≥ 0 || continue
            # record version choice
            S[r] = j
            # advance dominance frontier (if necessary & possible)
            l = length(solutions)
            d′ = d
            if d ≤ l
                # check for advancement
                j < solutions[d][P[j]] || continue
                # advance past already-dominated solutions
                while true
                    d′ += 1
                    d′ ≤ l || break
                    S′ = solutions[d′]
                    any(k < S′[P[k]] for r′ = 1:r for k = S[r′]) || break
                end
            end
            # check for a complete solution
            if r == M
                if d′ > l # it's optimal, save it!
                    push!(solutions, sort(S, by = k -> P[k]))
                end
                break
            end
            # restrict next subgraph to neighbors without too many conflicts
            for k = 1:N
                C[k] += X[k, j]
            end
            # recursion
            search!(c′, r + 1, d′)
            # restore the conflict count & level vectors
            for k = 1:N
                C[k] -= X[k, j]
            end
            # exclude vertex from future iterations
            C[j] += r*L
        end
        for j = 1:N
            dr = divrem(C[j], L)
            if dr[1] ≥ r
                C[j] = dr[2]
            end
        end
    end
    search!(c)

    foreach(sort!, solutions)
    return sort!(solutions)
end

end # module


text:
function countRGB(s)
    numR = numG = numB = 0
    for c in s
        if c == 'R'
            numR += 1
        elseif c == 'G'
            numG += 1
        elseif c == 'B'
            numB += 1
        end
    end
    (numR, numG, numB)
end

function numequaldistance(s)
    ret = 0
    for i = 1:length(s)-2
        for d = 1:(length(s)-i)÷2
            j = i+d
            k = j+d
            if s[i] != s[j] && s[j] != s[k] && s[k] != s[i]
                ret += 1
            end
        end
    end
    ret
end

function solve()
    readline()
    s = readline()
    numR, numG, numB = countRGB(s)
    ned = numequaldistance(s)
    numR*numG*numB - ned
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
using Dates

function main()
    dtstr = "March 7 2009 7:30pm" # Base.Dates doesn't handle "EST"
    cleandtstr = replace(dtstr, r"(am|pm)"i, "")
    dtformat = dateformat"U dd yyyy HH:MM"
    dtime = parse(DateTime, cleandtstr, dtformat) +
        Hour(12 * contains(dtstr, r"pm"i)) # add 12h for the pm
    println(Dates.format(dtime + Hour(12), dtformat))
end

main()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# kolmogorov_smirnov.jl
# Kolmogorov–Smirnov
#
# Copyright (C) 2014   Christoph Sawade
#
# Permission is hereby granted, free of charge, to any person obtaining
# a copy of this software and associated documentation files (the
# "Software"), to deal in the Software without restriction, including
# without limitation the rights to use, copy, modify, merge, publish,
# distribute, sublicense, and/or sell copies of the Software, and to
# permit persons to whom the Software is furnished to do so, subject to
# the following conditions:
#
# The above copyright notice and this permission notice shall be
# included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
# LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
# OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
# WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

export
    ExactOneSampleKSTest,
    ApproximateOneSampleKSTest, ApproximateTwoSampleKSTest

abstract type KSTest <: HypothesisTest end
abstract type ApproximateKSTest <: KSTest end
abstract type ExactKSTest <: KSTest end

population_param_of_interest(x::KSTest) = ("Supremum of CDF differences", 0.0, x.δ) # parameter of interest: name, value under h0, point estimate
default_tail(test::KSTest) = :both

## ONE SAMPLE KS-TEST

# compute supremum of differences between target and empirical cdf before and after the jump of the empirical cdf.
function ksstats(x::AbstractVector{T}, d::UnivariateDistribution) where T<:Real
    n = length(x)
    cdfs = cdf.(Ref(d), sort(x))
    δp = maximum((1:n) / n - cdfs)
    δn = -minimum((0:n-1) / n - cdfs)
    δ = max(δn, δp)
    (n, δ, δp, δn)
end

### EXACT KOLMOGOROV SMIRNOV TEST

struct ExactOneSampleKSTest <: ExactKSTest
    n::Int      # number of observations
    δ::Float64  # supremum of CDF differences
    δp::Float64 # supremum of the positive CDF differences
    δn::Float64 # supremum of the negative CDF differences
end

"""
    ExactOneSampleKSTest(x::AbstractVector{<:Real}, d::UnivariateDistribution)

Perform a one-sample exact Kolmogorov–Smirnov test of the null hypothesis that the data in
vector `x` comes from the distribution `d` against the alternative hypothesis that the
sample is not drawn from `d`.

Implements: [`pvalue`](@ref)
"""
function ExactOneSampleKSTest(x::AbstractVector{T}, d::UnivariateDistribution) where T<:Real
    if length(x) > length(unique(x))
        @warn("This test is inaccurate with ties")
    end

    ExactOneSampleKSTest(ksstats(x, d)...)
end

testname(::ExactOneSampleKSTest) = "Exact one sample Kolmogorov-Smirnov test"

function show_params(io::IO, x::ExactOneSampleKSTest, ident="")
    println(io, ident, "number of observations:   $(x.n)")
end

function pvalue(x::ExactKSTest; tail=:both)
    if tail == :left
        pvalue(KSOneSided(x.n), x.δn; tail=:right)
    elseif tail == :right
        pvalue(KSOneSided(x.n), x.δp; tail=:right)
    elseif tail == :both
        pvalue(KSDist(x.n), x.δ; tail=:right)
    else
        throw(ArgumentError("tail=$(tail) is invalid"))
    end
end

### APPROXIMATE KOLMOGOROV SMIRNOV TEST

struct ApproximateOneSampleKSTest <: ApproximateKSTest
    n::Int      # number of observations
    δ::Float64  # supremum of CDF differences
    δp::Float64 # supremum of the positive CDF differences
    δn::Float64 # suproemum of the negative CDF differences
end

"""
    ApproximateOneSampleKSTest(x::AbstractVector{<:Real}, d::UnivariateDistribution)

Perform an asymptotic one-sample Kolmogorov–Smirnov test of the null hypothesis that the
data in vector `x` comes from the distribution `d` against the alternative hypothesis
that the sample is not drawn from `d`.

Implements: [`pvalue`](@ref)
"""
function ApproximateOneSampleKSTest(x::AbstractVector{T}, d::UnivariateDistribution) where T<:Real
    if length(x) > length(unique(x))
        @warn("This test is inaccurate with ties")
    end

    ApproximateOneSampleKSTest(ksstats(x, d)...)
end

testname(::ApproximateOneSampleKSTest) = "Approximate one sample Kolmogorov-Smirnov test"

function show_params(io::IO, x::ApproximateOneSampleKSTest, ident="")
    println(io, ident, "number of observations:   $(x.n)")
    println(io, ident, "KS-statistic:             $(sqrt(x.n)*x.δ)")
end

# one-sided: http://www.encyclopediaofmath.org/index.php/Kolmogorov-Smirnov_test
function pvalue(x::ApproximateOneSampleKSTest; tail=:both)
    if tail == :left
        exp(-2*x.n*x.δn^2)
    elseif tail == :right
        exp(-2*x.n*x.δp^2)
    elseif tail == :both
        pvalue(Kolmogorov(), sqrt(x.n)*x.δ; tail=:right)
    else
        throw(ArgumentError("tail=$(tail) is invalid"))
    end
end

## TWO SAMPLE KS-TEST

### APPROXIMATE KOLMOGOROV SMIRNOV TEST

struct ApproximateTwoSampleKSTest <: ApproximateKSTest
    n_x::Int    # number of observations
    n_y::Int    # number of observations
    δ::Float64  # supremum of CDF differences
    δp::Float64 # supremum of the positive CDF differences
    δn::Float64 # suproemum of the negative CDF differences
end

"""
    ApproximateTwoSampleKSTest(x::AbstractVector{<:Real}, y::AbstractVector{<:Real})

Perform an asymptotic two-sample Kolmogorov–Smirnov-test of the null hypothesis that `x`
and `y` are drawn from the same distribution against the alternative hypothesis that they
come from different distributions.

Implements: [`pvalue`](@ref)

# External links

  * [Approximation of one-sided test (Encyclopedia of Mathematics)
    ](https://www.encyclopediaofmath.org/index.php/Kolmogorov-Smirnov_test)
"""
function ApproximateTwoSampleKSTest(x::AbstractVector{T}, y::AbstractVector{S}) where {T<:Real, S<:Real}
    n_x, n_y = length(x), length(y)
    if n_x+n_y > length(unique([x; y]))
        @warn("This test is inaccurate with ties")
    end

    ApproximateTwoSampleKSTest(ksstats(x, y)...)
end

testname(::ApproximateTwoSampleKSTest) = "Approximate two sample Kolmogorov-Smirnov test"

function show_params(io::IO, x::ApproximateTwoSampleKSTest, ident="")
    n = x.n_x*x.n_y/(x.n_x+x.n_y)
    println(io, ident, "number of observations:   [$(x.n_x),$(x.n_y)]")
    println(io, ident, "KS-statistic:              $(sqrt(n)*x.δ)")
end

function pvalue(x::ApproximateTwoSampleKSTest; tail=:both)
    n = x.n_x*x.n_y/(x.n_x+x.n_y)
    if tail == :left
        exp(-2*n*x.δn^2)
    elseif tail == :right
        exp(-2*n*x.δp^2)
    elseif tail == :both
        pvalue(Kolmogorov(), sqrt(n)*x.δ; tail=:right)
    else
        throw(ArgumentError("tail=$(tail) is invalid"))
    end
end

# compute supremum of differences between empirical cdfs.
function ksstats(x::AbstractVector{T}, y::AbstractVector{S}) where {T<:Real, S<:Real}
    n_x, n_y = length(x), length(y)
    sort_idx = sortperm([x; y])
    pdf_diffs = [ones(n_x)/n_x; -ones(n_y)/n_y][sort_idx]
    cdf_diffs = cumsum(pdf_diffs)
    δp = maximum(cdf_diffs)
    δn = -minimum(cdf_diffs)
    δ = max(δp, δn)
    (n_x, n_y, δ, δp, δn)
end


text:
BroadcastStyle(::Type{<:judiMultiSourceVector}) = ArrayStyle{judiMultiSourceVector}()

function similar(bc::Broadcast.Broadcasted{ArrayStyle{judiMultiSourceVector}}, ::Type{ElType}) where ElType
    # Scan the inputs for the ArrayAndChar:
    A = find_msv(bc)
    return similar(A, ElType)
end

"`A = find_msv(As)` returns the first judiMultiSourceVector among the arguments."
find_msv(bc::Base.Broadcast.Broadcasted) = find_msv(bc.args)
find_msv(args::Tuple) = find_msv(find_msv(args[1]), Base.tail(args))
find_msv(x) = x
find_msv(::Tuple{}) = nothing
find_msv(a::judiMultiSourceVector, rest) = a
find_msv(::Any, rest) = find_msv(rest)

extrude(x::judiMultiSourceVector) = extrude(x.data)

# Add broadcasting by hand due to the per source indexing
for func ∈ [:lmul!, :rmul!, :rdiv!, :ldiv!]
    @eval begin
        $func(ms::judiMultiSourceVector, x::Number) = $func(ms.data, x)
        $func(x::Number, ms::judiMultiSourceVector) = $func(x, ms.data)
    end
end

# Broadcasted custom type
check_compat() = true
get_src(ms::judiMultiSourceVector, j) = ms.data[j]
get_src(v::Vector{<:Array}, j) = v[j]
get_src(v::Array{T, N}, j) where {T<:Number, N} = v
get_src(n::Number, j) = n

struct MultiSource <: Base.AbstractBroadcasted
    m1
    m2
    op
end

function materialize(bc::MultiSource)
    m1, m2 = materialize(bc.m1), materialize(bc.m2)
    ms = similar(find_msv((m1, m2)))
    check_compat(m1, m2)
    for i=1:ms.nsrc
        ms.data[i] = materialize(broadcasted(bc.op, get_src(m1, i), get_src(m2, i)))
    end
    ms
end

function materialize!(ms::judiMultiSourceVector, bc::MultiSource)
    m1, m2 = materialize(bc.m1), materialize(bc.m2)
    check_compat(ms, m1)
    check_compat(ms, m2)
    for i=1:ms.nsrc
        broadcast!(bc.op, ms.data[i], get_src(m1, i), get_src(m2, i))
    end
    nothing
end

for op ∈ [:+, :-, :*, :/]
    # Two multi source vectors
    @eval begin
        $(op)(ms1::judiMultiSourceVector, ms2::judiMultiSourceVector) = materialize(broadcasted($op, ms1, ms2))
    end
    # External types and broadcasting
    for LT in [Number, Real, Vector{<:Array}, Array{<:Number, <:Integer}]
        # +/*/... julia type vs multi source
        @eval $(op)(ms1::$(LT), ms2::judiMultiSourceVector) = materialize(broadcasted($op, ms1, ms2))
        @eval $(op)(ms1::judiMultiSourceVector, ms2::$(LT)) = materialize(broadcasted($op, ms1, ms2))
        # broadcasted julia type vs multi source vector or broadcast
        for MS in [judiMultiSourceVector, MultiSource]
            @eval broadcasted(::typeof($op), ms1::$(LT), ms2::$(MS)) = MultiSource(ms1, ms2, $op)
            @eval broadcasted(::typeof($op), ms1::$(MS), ms2::$(LT)) = MultiSource(ms1, ms2, $op)
        end
    end
    # multi source with multi source
    @eval broadcasted(::typeof($op), ms1::judiMultiSourceVector, ms2::MultiSource) = MultiSource(ms1, ms2, $op)
    @eval broadcasted(::typeof($op), ms1::MultiSource, ms2::judiMultiSourceVector) = MultiSource(ms1, ms2, $op)
    @eval broadcasted(::typeof($op), ms1::judiMultiSourceVector, ms2::judiMultiSourceVector) = MultiSource(ms1, ms2, $op)
    @eval broadcasted(::typeof($op), ms1::MultiSource, ms2::MultiSource) = MultiSource(ms1, ms2, $op)
end


text:
@inline as_to_rad(x) = deg2rad(x / 3600.0)
@inline as_to_rad_mod(x) = deg2rad(rem(x, 1296000.0) / 3600.0)
@inline rad_to_as(x) = rad2deg(x) * 3600.0

@inline μas_to_rad(x) = deg2rad(x / 3600e6)
@inline rad_to_μas(x) = rad2deg(x) * 3600e6

@inline jc(ep) = ((julian1(ep) - J2000) + julian2(ep)) / DAYS_PER_CENTURY

function fundamental_arguments(tdb_jc)

    # Mean Anomaly of the Moon (l, F1)
    M_m  = as_to_rad( @evalpoly(tdb_jc,   485_868.249_036, 1_717_915_923.217_8,  31.879_2,  0.051_635, -0.000_244_70) )

    # Mean Anomaly of the Sun (l', F2)
    M_s  = as_to_rad( @evalpoly(tdb_jc, 1_287_104.793_048,   129_596_581.048_1,  -0.553_2,  0.000_136, -0.000_011_49) )

    # Mean Argument of Latitude of the Moon (F, F3)
    u_sm = as_to_rad( @evalpoly(tdb_jc,   335_779.526_232, 1_739_527_262.847_8, -12.751_2, -0.001_037,  0.000_004_17) )

    # Mean Elongation of the Moon from the Sun (D, F4)
    D_s  = as_to_rad( @evalpoly(tdb_jc, 1_072_260.703_692, 1_602_961_601.209_0,  -6.370_6,  0.006_593, -0.000_031_69) )

    # Mean Longitude of the ascending Node of the Moon (Ω, F5)
    Ω_m  = as_to_rad( @evalpoly(tdb_jc,   450_160.398_036,    -6_962_890.543_1,   7.472_2,  0.007_702, -0.000_059_39) )

    # Mean Heliocentric Longitudes of the Planets
    # Souchay et al. (1999)
    λ_M1 = @evalpoly(tdb_jc, 4.402_608_842, 2608.790_314_157_4) # Mercury
    λ_M2 = @evalpoly(tdb_jc, 3.176_146_697, 1021.328_554_621_1) # Venus
    λ_M3 = @evalpoly(tdb_jc, 1.753_470_314,  628.307_584_999_1) # Earth
    λ_M4 = @evalpoly(tdb_jc, 6.203_480_913,  334.061_242_670_0) # Mars
    λ_M5 = @evalpoly(tdb_jc, 0.599_546_497,   52.969_096_264_1) # Jupiter
    λ_M6 = @evalpoly(tdb_jc, 0.874_016_757,   21.329_910_496_0) # Saturn
    λ_M7 = @evalpoly(tdb_jc, 5.481_293_872,    7.478_159_856_7) # Uranus
    λ_M8 = @evalpoly(tdb_jc, 5.311_886_287,    3.813_303_563_8) # Neptune

    # General Precession
    # Kinoshita and Souchay (1990)
    p_λ  = @evalpoly(tdb_jc, 0.0, 0.024_381_75, 0.000_005_386_91)

    return SVector(M_m, M_s, u_sm, D_s, Ω_m, λ_M1, λ_M2, λ_M3, λ_M4, λ_M5, λ_M6, λ_M7, λ_M8, p_λ)
end


text:
using Makie

cmp = :pu_or
# sheen = 32.0f0
sheen = 16.0f0
camera_translation = (0.0, 0.0, 0.0)
camera_rotation = (0.0, 0.2, 0.0)

# xlims = [-0.1,  1.1, -0.1,  1.1, -0.1,  1.1, -0.1,  1.1]
# ylims = [-0.1,  1.1,  1.1, -0.1,  1.1,  1.1, -0.1, -0.1]
# zlims = [ 0.6,  0.6,  0.6,  0.6, -0.6, -0.6, -0.6, -0.6]

lim = FRect3D((-0.1, -0.1, -0.6), (1.2, 1.2, 1.2))

ϕ01(x,y) = 0.5 * Float32(0 < x < 1) * Float32(0 <  y <  1)

ϕ11(x,y) = 0.5  * Float32(0 < x < 1/2) * Float32(0 < y < 1/2)
ϕ12(x,y) = 0.5  * Float32(1/2 < x < 1) * Float32(0 < y < 1/2)
ϕ13(x,y) = 0.5  * Float32(0 < x < 1/2) * Float32(1/2 < y < 1)
ϕ14(x,y) = 0.5  * Float32(1/2 < x < 1) * Float32(1/2 < y < 1)

ϕ21(x,y) = 0.5  * Float32(2/4 < x < 3/4) * Float32(2/4 < y < 3/4)
ϕ22(x,y) = 0.5  * Float32(3/4 < x < 4/4) * Float32(2/4 < y < 3/4)
ϕ23(x,y) = 0.5  * Float32(2/4 < x < 3/4) * Float32(3/4 < y < 4/4)
ϕ24(x,y) = 0.5  * Float32(3/4 < x < 4/4) * Float32(3/4 < y < 4/4)

vx = -0.01: 0.0103 : 1.01
vy = -0.01: 0.0103 : 1.01

plot1 = Scene();
z = Float32[ϕ01(x,y) for x in vx, y in vy];
surface!(plot1, vx, vy, z, colormap=cmp, shininess=sheen, show_axis=false, shading=false, colorrange=(-0.5, 0.5), limits=lim);
translate_cam!(plot1, camera_translation);
rotate_cam!(plot1, camera_rotation);

plot2 = Scene();
z = Float32[ϕ11(x,y) + ϕ12(x,y) - ϕ13(x,y) - ϕ14(x,y) for x in vx, y in vy];
surface!(plot2, vx, vy, z, colormap=cmp, show_axis=false, shading=false, colorrange=(-0.5, 0.5), limits=lim);
translate_cam!(plot2, camera_translation);
rotate_cam!(plot2, camera_rotation);

plot3 = Scene();
z = Float32[ϕ11(x,y) - ϕ12(x,y) + ϕ13(x,y) - ϕ14(x,y) for x in vx, y in vy];
surface!(plot3, vx, vy, z, colormap=cmp, show_axis=false, colorrange=(-0.5, 0.5), limits=lim);
translate_cam!(plot3, camera_translation);
rotate_cam!(plot3, camera_rotation);

plot4 = Scene();
z = Float32[ϕ11(x,y) - ϕ12(x,y) - ϕ13(x,y) + ϕ14(x,y) for x in vx, y in vy];
surface!(plot4, vx, vy, z, colormap=cmp, show_axis=false, colorrange=(-0.5,0.5), limits=lim);
translate_cam!(plot4, camera_translation);
rotate_cam!(plot4, camera_rotation);

plot5 = Scene();
z = Float32[ϕ21(x,y) + ϕ22(x,y) - ϕ23(x,y) - ϕ24(x,y) for x in vx, y in vy];
surface!(plot5, vx, vy, z, colormap=cmp, show_axis=false, colorrange=(-0.5, 0.5), limits=lim);
translate_cam!(plot5, camera_translation);
rotate_cam!(plot5, camera_rotation);

plot6 = Scene();
z = Float32[ϕ21(x,y) - ϕ22(x,y) + ϕ23(x,y) - ϕ24(x,y) for x in vx, y in vy];
surface!(plot6, vx, vy, z, colormap=cmp, show_axis=false, colorrange=(-0.5, 0.5), limits=lim);
translate_cam!(plot6, camera_translation);
rotate_cam!(plot6, camera_rotation);

plot7 = Scene();
z = Float32[ϕ21(x,y) - ϕ22(x,y) - ϕ23(x,y) + ϕ24(x,y) for x in vx, y in vy];
surface!(plot7, vx, vy, z, colormap=cmp, show_axis=false, colorrange=(-0.5, 0.5), limits=lim);
translate_cam!(plot7, camera_translation);
rotate_cam!(plot7, camera_rotation);

out_scene = vbox(plot1, plot2, plot3, plot4, plot5, plot6, plot7, parent=Scene(resolution=(2400, 400), clear=false))
save("./out/plots/basis_functions.png", out_scene)


text:
using MPIMeasurements
using Base.Test
using Unitful
using Compat

# define Grid
shp = [3,3,3]
fov = [3.0,3.0,3.0]Unitful.mm
ctr = [0,0,0]Unitful.mm
positions = RegularGridPositions(shp,fov,ctr)

# create Robot
dR = DummyRobot()

res = performTour!(dR, hallSensorRegularScanner, positions, DummyMeasObj())
#move back to park position after measurement has finished
movePark(robot)


text:
using Test, DocFree

out = plusTwo(3)

@test out == 5


text:
#### thermodynamics

export new_thermo_state, recover_thermo_state

"""
    new_thermo_state(atmos::AtmosModel, state::Vars, aux::Vars)

Create a new thermodynamic state, based on the `state`, and _not_
the `aux` state.

!!! note
    This method calls the iterative saturation adjustment
    procedure for EquilMoist models.
"""
new_thermo_state(atmos::AtmosModel, state::Vars, aux::Vars) =
    new_thermo_state(atmos, atmos.moisture, state, aux)

function new_thermo_state(
    atmos::AtmosModel,
    moist::DryModel,
    state::Vars,
    aux::Vars,
)
    e_int = internal_energy(atmos, state, aux)
    return PhaseDry(atmos.param_set, e_int, state.ρ)
end

function new_thermo_state(
    atmos::AtmosModel,
    moist::EquilMoist,
    state::Vars,
    aux::Vars,
)
    e_int = internal_energy(atmos, state, aux)
    return PhaseEquil(
        atmos.param_set,
        e_int,
        state.ρ,
        state.moisture.ρq_tot / state.ρ,
        moist.maxiter,
        moist.tolerance,
    )
end

function new_thermo_state(
    atmos::AtmosModel,
    moist::NonEquilMoist,
    state::Vars,
    aux::Vars,
)
    e_int = internal_energy(atmos, state, aux)
    q = PhasePartition(
        state.moisture.ρq_tot / state.ρ,
        state.moisture.ρq_liq / state.ρ,
        state.moisture.ρq_ice / state.ρ,
    )

    return PhaseNonEquil{eltype(state), typeof(atmos.param_set)}(
        atmos.param_set,
        e_int,
        state.ρ,
        q,
    )
end

"""
    recover_thermo_state(atmos::AtmosModel, state::Vars, aux::Vars)

Recover the thermodynamic state, based on the `state` _and_
the `aux` state, based on the existing temperature.

!!! note
    This method does _not_ call the iterative saturation adjustment
    procedure.
"""
recover_thermo_state(atmos::AtmosModel, state::Vars, aux::Vars) =
    recover_thermo_state(atmos, atmos.moisture, state, aux)

function recover_thermo_state(
    atmos::AtmosModel,
    moist::DryModel,
    state::Vars,
    aux::Vars,
)
    e_int = internal_energy(atmos, state, aux)
    return PhaseDry(atmos.param_set, e_int, state.ρ)
end

function recover_thermo_state(
    atmos::AtmosModel,
    moist::EquilMoist,
    state::Vars,
    aux::Vars,
)
    e_int = internal_energy(atmos, state, aux)
    return PhaseEquil{eltype(state), typeof(atmos.param_set)}(
        atmos.param_set,
        e_int,
        state.ρ,
        state.moisture.ρq_tot / state.ρ,
        aux.moisture.temperature,
    )
end

function recover_thermo_state(
    atmos::AtmosModel,
    moist::NonEquilMoist,
    state::Vars,
    aux::Vars,
)
    e_int = internal_energy(atmos, state, aux)
    q = PhasePartition(
        state.moisture.ρq_tot / state.ρ,
        state.moisture.ρq_liq / state.ρ,
        state.moisture.ρq_ice / state.ρ,
    )

    return PhaseNonEquil{eltype(state), typeof(atmos.param_set)}(
        atmos.param_set,
        e_int,
        state.ρ,
        q,
    )
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA