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

using Markdown
using InteractiveUtils

# ╔═╡ 2b6de7b6-832e-420f-b73c-9ebac9901547
begin
	import Distributions as di
	using Random
	using PlutoUI
end

# ╔═╡ 8b908968-8e7f-11ec-08b9-797a5b416d46
md"""
# The Basics
"There are three types of lies: lies, damned lies and statistics." -- Benjamin Disraeli (1804-1881)

This quote by Disraeli expresses a common sentiment.  Especially in the biosciences, statistics is often used in a cookie-cutter way to distinguish between significant and insignificant results, but we are all left with a bad taste in our mouth when trying to deeply think about the underlying concepts.
In this class, I would like to instead pursue the subject of data analysis using a probabilistic approach that was formulated by Bayes (published in 1763) and later refined by Laplace (1812).
# Introduction
In Science we are faced with a difficult task.  We are observing the world through experiments but what we would like to derive from this are its laws.  Unfortunately, this task is not only difficult, but it is actually an impossible task.  The reason for this is that there are many possible laws that could explain all the data that we collect.  This is in contrast to mathematical logic, where everything can be derived from a set of axioms.  The mathematician is asking: given a cause, can we work out its consequences? Scientist do the reverse: Given that certain effects have been observed, what are the underlying causes?
# Coin flip example
To illustrate the problem, we are going to start with a familiar problem: the coin flip.  Let us assume that you are visiting Las Vegas and someone invited you to a game in which you bet one dollar and after a coin flip you get paid two dollars for heads but lose your bet when tails come.  It is pretty clear from the outset that if the game uses a "fair" coin (tails:head = 1:1) then you should not win or lose any money.  To make it more interesting the person that invited you to play flips the coin 10 times and you observe that heads came 6 times and tails came 4 times.  How does this information help you to decide whether to play this game since you suspect that the coin is probably not fair (it is Las Vegas after all).

The first thing that we can do is to ask a different question: What is the probability that 6 heads out of ten resulted from a fair coin (some statistician would call the fair coin the "null hypothesis").  Here we are going to employ the Bionomial distribution which gives us what we are looking for:

$$P( X = k \mid N,p) =  \frac{N!}{(N-k)!k!}  p^k(1-p)^{N-k}$$

Notice that we wrote the Binomial distribution as a conditional probability to indicate that it represents the probability of $X=k$ given $N$ and $p$. Below we are calculating the odds for each situation:
"""

# ╔═╡ a7995c55-88d7-435c-839d-22278c0667a2
gambler = di.Binomial(100,0.5)

# ╔═╡ 0056269b-7f58-466f-acb3-f1a73f736133
with_terminal() do
	for i in 0:5:100
		println(i," ",di.pdf(gambler,i))
	end
end

# ╔═╡ 3c768871-2f62-4ae2-980b-aa9494e32a80
md"""
We just learned how to use Bayes theorem to calculate the posterior probability distribution that allows us to make a quantitative assessment of our decision whether we should accept this gambling game or not.  Such decisions are constantly made in the financial industry and Bayesian methods are commonly used in this environment.  In the financial industry the goal is to maximize profit and to minimize losses, so it is imperative to use quantitative tools.  In science, we have the problem that scientists are not in the same situation.  For scientists it is more important to get more grants, to write more papers since this is what is seen as the deliverables of their efforts.  What I would like to teach you is that by following Bayes theorem, we can improve the reliability of our research without sacrificing academic productivity.  We need to know our data though.

# Some important properties of Bayes theorem

Before we continue, I would like to point out a few important properties of Bayes theorem and its application.  The first property has to do with the way the posterior is calculated.  Let us go back to the coin flipping experiment and assume that we did an experiment with the coin in question that initially involved N=10 coin flips.  We used our prior and by combining it with the Bionomial likelihood, we were able to calculate the posterior:

$$P(p \mid \alpha,\beta, N, k) =  \frac{\Gamma(\alpha + \beta +N)}{\Gamma(\alpha+k)\Gamma(\beta + N -k)}  p^{\alpha-1+k}(1-p)^{\beta-1+N-k}$$

Lets now say that we are not quite sure about what to do and we ask the Gambler to let us flip another 10 times.  What we would do is to take the posterior that we derived, use it as our prior and then again use the Bionomial distribution as likelihood.  When looking back at the equations, it is pretty obvious that it does not matter whether you analyze the data sequentially (flip 10 times, calculate the posterior, flip 10 more times, update the posterior) or analyze it in one shot (take result of 20 flips, use prior, calculate posterior).  This is an important property of Bayes theorem.  It should not matter how you do the analysis.  It would be perfectly fine to flip the coin and update the posterior each time.  The result should be the same.  This statement is pretty obvious in our case but it holds in general. Later in this class we will see examples that are much more complicated than simple coin flips and here it sometimes makes sense to perform a complicated analysis with parts of the data and then repeat the process, rather than take all the data at once.  Also, in science it often happens that you measure some data, and you want to see whether it makes sense to persue this direction further.  That is perfectly ok.  If your inital posterior looks promising then go ahead and do some more measurements.  Bayesian analysis allows you to improve your statistics by repeating the experiment.

The next property has to do with the shape of the posterior.  We have seen last lecture that a typical outcome of your posterior can be described by the maximum of the distribution (or the 50th percentile) and its width.  In our coin flipping experiment, the maximum will tell you whether or not the gamble will be profitable for you and the width will tell you about the certainty or uncertainty that this will happen.  In the financial market you would be interested whether a stock is likely to go up or down, and what the volatility is.  So lets use our posterior to calculate these quantities.

The first thing that we want to do is to take the natural logarithm of the probability distribution.  This will not change the maximum of $P(x)$ since the transformation is monotonic.  This transformation is often using in computational approaches to probability because it turns multiplication into additions.  You will see another reason in the next paragraph.

$$logp(x) = \log(P(x | D)$$

Now we can do a Taylor expansion around $x=x_{0}$.  Since we are concerned about a peaked probability distribution we only have a constant term and a quadratic term.

$$logp(x) = logp(x_{0}) + \frac{1}{2}\frac{d^{2}logp(x)}{dx^{2}}\big | _{x_{0}}(x-x_{0})^{2}+\dots$$

Putting this back into the previous equation yields:

$$P(x|D) \propto A \exp\big[\frac{1}{2}\frac{d^{2}logp(x)}{dx^{2}}\big | _{x_{0}}(x-x_{0})^{2}\big]$$

If you look at this equation you will recognize that what we have just done is to approximate our pdf with a Gaussian distribution.  Gaussian distributions are very commonly found, and the reason for this is beyond this class, but suffice to say that in many cases posterior distributions tend to approach Gaussian distributions for large N (law of large numbers).  Let look at a Gaussian distribution:

$$P(x \mid \mu, \sigma) = \frac{1}{{\sigma \sqrt {2\pi } }}e^{{{ - \left( {x - \mu } \right)^2 } /{2\sigma ^2 }}}$$
"""

# ╔═╡ a73a746c-cfdf-434b-89ba-a77b98543a81
md"""
The maximum of this distribution can be found by finding the root of the first derivative with respect to $x$, and the standard deviation can be found by taking the second derivative.  Here is the calculation.  You can check it if you want yourself.  In particuar,

$$\sigma = \big(-\frac{d^{2}logp(x)}{dx^{2}}\big | _{x_{0}}\big)^{-\frac{1}{2}}$$
"""

# ╔═╡ 8f0d784e-7f3f-4aee-80df-a9e2f5f31157
md"""
Once we know $\mu$ and $\sigma$ we can report our posterior pdf as:

$$x = x_{0} \pm \sigma$$

Since we approximated our pdf with a Gaussian we also know that true value of x lies in within $\pm \sigma$ $67$ of the time and within $\pm 2\sigma$ $95$ of the time.

Lets get back to our example posterior and see what happens

$$P(p \mid |N,k) \propto p^{k}(1-p)^{N-k}$$ with
$$logp = c + k\log(p)+(N-k)\log(1-p)$$
"""

# ╔═╡ 5b24e37c-a5ce-40b2-97db-d850ac74523f
md"""
The derivatives are:

$$\frac{dlogp}{dp}=\frac{k}{p}-\frac{N-k}{1-p}$$
$$\frac{d^2 logp}{dp^{2}}=-\frac{k}{p^2}-\frac{N-k}{(1-p)^2}$$

The maximum can be found

$$\frac{dlogp}{dp}\big | _{p_{0}} = \frac{k}{p_{0}}-\frac{N-k}{1-p_{0}}=0$$
$$p_{0}=\frac{k}{N}$$

which is not surprising.  The second derivate is then:
$$\frac{d^2 logp}{dp^2}\big | _{p_{0}}=-\frac{N}{p_{0}(1-p_{0})}$$

$$\sigma = \sqrt{\frac{p_{0}(1-p_{0})}{N}}$$

For larger $N$, $p_{0}$ will not change much and the numerator will be essentially constant.  In this regime, the standard deviation $\mu$ is inversely proportional to the square root of the number of tries
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[compat]
Distributions = "~0.25.48"
PlutoUI = "~0.7.34"
"""

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[DensityInterface]]
deps = ["InverseFunctions", "Test"]
git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b"
uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
version = "0.4.0"

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

[[Distributions]]
deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"]
git-tree-sha1 = "38012bf3553d01255e83928eec9c998e19adfddf"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.48"

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

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

[[FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "deed294cde3de20ae0b2e0355a6c4e1c6a5ceffc"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.12.8"

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

[[OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

[[OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

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

[[PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "ee26b350276c51697c9c2d88a072b339f9f03d73"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.5"

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

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

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

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

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

[[QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "78aadffb3efd2155af139781b8a8df1ef279ea39"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.2"

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

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

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

[[Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

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

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

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

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

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

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

[[SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "8d0c8e3d0ff211d9ff4a0c2307d876c99d10bdf1"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "2.1.2"

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

[[StatsAPI]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "c3d8ba7f3fa0625b062b82853a7d5229cb728b6b"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.2.1"

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

[[StatsFuns]]
deps = ["ChainRulesCore", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "f35e1879a71cca95f4826a14cdbf0b9e253ed918"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.15"

[[SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

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

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

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

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

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

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

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

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

# ╔═╡ Cell order:
# ╠═2b6de7b6-832e-420f-b73c-9ebac9901547
# ╟─8b908968-8e7f-11ec-08b9-797a5b416d46
# ╠═a7995c55-88d7-435c-839d-22278c0667a2
# ╠═0056269b-7f58-466f-acb3-f1a73f736133
# ╟─3c768871-2f62-4ae2-980b-aa9494e32a80
# ╟─a73a746c-cfdf-434b-89ba-a77b98543a81
# ╟─8f0d784e-7f3f-4aee-80df-a9e2f5f31157
# ╟─5b24e37c-a5ce-40b2-97db-d850ac74523f
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
using BcdV3
using StaticArrays
using MimuConstants
using LinearAlgebra
using Printf
using Dates
using Random
Random.seed!(1)

function InitialPosMIMU4444BT()
    dx = 0.0063; dy = 0.0063; dz = 0.001;
    r = [-1.5*dx  0.5*dy  1*dz;   #  1
         -1.5*dx  0.5*dy -1*dz;   #  2
         -1.5*dx  1.5*dy  1*dz;   #  3
         -1.5*dx  1.5*dy -1*dz;   #  4
         -1.5*dx -1.5*dy  1*dz;   #  5
         -1.5*dx -1.5*dy -1*dz;   #  6
         -1.5*dx -0.5*dy  1*dz;   #  7
         -1.5*dx -0.5*dy -1*dz;   #  8
         -0.5*dx  0.5*dy  1*dz;   #  9
         -0.5*dx  0.5*dy -1*dz;   # 10
         -0.5*dx  1.5*dy  1*dz;   # 11
         -0.5*dx  1.5*dy -1*dz;   # 12
         -0.5*dx -1.5*dy  1*dz;   # 13
         -0.5*dx -1.5*dy -1*dz;   # 14
         -0.5*dx -0.5*dy  1*dz;   # 15
         -0.5*dx -0.5*dy -1*dz;   # 16
          0.5*dx  0.5*dy  1*dz;   # 17
          0.5*dx  0.5*dy -1*dz;   # 18
          0.5*dx  1.5*dy  1*dz;   # 19
          0.5*dx  1.5*dy -1*dz;   # 20
          0.5*dx -1.5*dy  1*dz;   # 21
          0.5*dx -1.5*dy -1*dz;   # 22
          0.5*dx -0.5*dy  1*dz;   # 23
          0.5*dx -0.5*dy -1*dz;   # 24
          1.5*dx  0.5*dy  1*dz;   # 25
          1.5*dx  0.5*dy -1*dz;   # 26
          1.5*dx  1.5*dy  1*dz;   # 27
          1.5*dx  1.5*dy -1*dz;   # 28
          1.5*dx -1.5*dy  1*dz;   # 29
          1.5*dx -1.5*dy -1*dz;   # 30
          1.5*dx -0.5*dy  1*dz;   # 31
          1.5*dx -0.5*dy -1*dz]'; # 32
    return r
end

function get_normals_icosahedron()
    f = [
        0.1876 -0.7947 0.5774
        0.6071 -0.7947 0.0000
        -0.4911 -0.7947 0.3568
        -0.4911 -0.7947 -0.3568
        0.1876 -0.7947 -0.5774
        0.9822 -0.1876 0.0000
        0.3035 -0.1876 0.9342
        -0.7946 -0.1876 0.5774
        -0.7946 -0.1876 -0.5774
        0.3035 -0.1876 -0.9342
        0.7946 0.1876 0.5774
        -0.3035 0.1876 0.9342
        -0.9822 0.1876 0.0000
        -0.3035 0.1876 -0.9342
        0.7946 0.1876 -0.5774
        0.4911 0.7947 0.3568
        -0.1876 0.7947 0.5774
        -0.6071 0.7947 0.0000
        -0.1876 0.7947 -0.5774
        0.4911 0.7947 -0.3568
    ]
    f_norm = f ./ sqrt.(sum(f.^2, dims = 2))
    return f_norm'
end

function get_callback_trace(trace, trace_p, show_every = 1)
    function callback(k::Int, dx, dlogp, log_p, r, T, b)
        trace[:f][k] = log_p
        trace[:f_abs][k] = dlogp
        trace[:f_rel][k] = abs(dlogp)/abs(log_p)

        trace[:x_abs][k] = dx

        if mod(k,show_every) == 0
            time = Dates.format(Dates.now(), "yyyy-mm-dd HH:MM:SS")
            @printf("%20s %7d logp %.16e dlogp % .8e ||dx|| %.8e\n",
                    time, k, log_p, dlogp, trace[:x_abs][k])

            # C = 0.0
            # for k = 1:8
            #     C += tr(T[k+8])
            # end
            # C /= (3*8)
            # @printf(" C_norm %.8e\n", C)
            push!(trace_p[:r], Array(r))
            push!(trace_p[:T], Array.(T))
            push!(trace_p[:b], Array.(b))
            push!(trace_p[:k], k)
        end
    end
    return callback
end

# Noise less simulations
TT = Float64
active_imus = [3,4,5,6, 27,28, 29, 30]
r_all = InitialPosMIMU4444BT()
r_m = r_all[:, active_imus]
r_ref = r_m[:,1]
r_m .-= r_ref
r_all .-= r_ref

g_m = get_normals_icosahedron()
N_orientions = size(g_m, 2)

sig_a = get_rms_noise_acc()
sig_g = get_rms_noise_gyro()
Na = size(r_m, 2)
Ng = 8

N = 40 # Dyn
N_per_orientation = Int(N/5)
Ns = ones(Int, N_orientions).*N_per_orientation
show_every = 1
N_bcd  = 3
tol_dlogp = 1.0e-9

Nt = Na + Ng
g_mag = 9.81

g = [SVector{3,TT}(g_m[:,k]) for k = 1:N_orientions] .|> x -> x/norm(x)*g_mag

b = [@SVector randn(3) for k = 1:Nt]
T = [I + 1.0e-3*@SMatrix randn(3,3) for k = 1:Nt]
T[1] = I + 1.0e-3*SMatrix{3,3, Float64}(triu(randn(3,3)))

r = SMatrix{3,Na, TT}(r_m)

w = [deg2rad(1000.0) * @SVector randn(3) for _ = 1:N]
w_dot = [deg2rad(40e3) * @SVector randn(3) for _ = 1:N]
s = [20.0 * @SVector randn(3) for _ = 1:N]

y_s_true = vcat(
    [T[k]*(-g[n]) + b[k] for k = 1:Na, n = 1:N_orientions],
    [b[k] for k = range(Na+1,length=Ng), n = 1:N_orientions],
)
y_d_true = vcat(
    [T[k]*(Ω(w[n])^2*r[:,k] + Ω(w_dot[n])*r[:,k] + s[n]) + b[k] for k = 1:Na, n = 1:N],
    [T[k]*w[n]+b[k] for k = range(Na+1,length=Ng), n = 1:N]
)

# y_s = [y_s_true[k,n] + 0.0*(k <= Na ? sig_a : sig_g)*randn(3)/sqrt(N_per_orientation) for k = 1:Nt, n = 1:N_orientions]
# y_d = [y_d_true[k,n] + 0.0*( k <= Na ? sig_a : sig_g)*randn(3) for k = 1:Nt, n = 1:N]
# y = hcat(y_s, y_d)

y = hcat(y_s_true, y_d_true)
Q_inv = [SMatrix{3,3,Float64}(I)/( k <= Na ? sig_a^2 :  sig_g^2) for k = 1:Nt]

b0 = [b_n + 1.0e-2*randn(3) for  b_n in b]
T0 = [T[k] + 1.0e-3*(k == 1 ? triu(randn(3,3)) : randn(3,3))  for k = 1:Nt]
r0_m = copy(r_m)
r0_m[:,2:end] += 1.0e-3*randn(3,Na-1)
r0 = [SVector{3}(r0_m[:,k]) for k = 1:Na]

trace = Dict(
    :f => zeros(N_bcd),
    :f_abs => zeros(N_bcd),
    :f_rel => zeros(N_bcd),
    :x_abs => zeros(N_bcd)
)
trace_p = Dict(:r => [], :T => [], :b => [], :k => [])
r_hat = copy(r0)
T_hat = copy(T0)
b_hat = copy(b0)

r_hat, T_hat, b_hat, eta = @time  bcd!(
    r_hat, T_hat, b_hat, y, Q_inv, Ns, g_mag, Val{3*(Na+Ng)}(),Val{Na}(),
    1.0e-12, 100, 1.0e-13, 200, tol_dlogp, 1.0-14, N_bcd;
    callback = get_callback_trace(
        trace, trace_p, show_every
    ), #callback
    verbose = false
);

ext: jl


text:
"""
Bridges `CP.IfThenElse` to reification.
"""
struct IfThenElse2ReificationBridge{T} <: MOIBC.AbstractBridge
    var_condition::MOI.VariableIndex
    var_true::MOI.VariableIndex
    var_false::MOI.VariableIndex
    var_condition_bin::MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}
    var_true_bin::MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}
    var_false_bin::MOI.ConstraintIndex{MOI.VariableIndex, MOI.ZeroOne}
    con_reif_condition::MOI.ConstraintIndex
    con_reif_true::MOI.ConstraintIndex
    con_reif_false::MOI.ConstraintIndex
    # Ideally, Vector{MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.Reification{<: MOI.AbstractSet}}},
    # but Julia has no notion of type erasure.
    con_if::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}
    con_else::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}}
end

function MOIBC.bridge_constraint(
    ::Type{IfThenElse2ReificationBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.IfThenElse{S1, S2, S3},
) where {T, S1, S2, S3}
    return MOIBC.bridge_constraint(
        IfThenElse2ReificationBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{IfThenElse2ReificationBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.IfThenElse{S1, S2, S3},
) where {T, S1, S2, S3}
    var_condition, var_condition_bin = MOI.add_constrained_variable(model, MOI.ZeroOne())
    var_true, var_true_bin = MOI.add_constrained_variable(model, MOI.ZeroOne())
    var_false, var_false_bin = MOI.add_constrained_variable(model, MOI.ZeroOne())

    f_scalars = MOIU.scalarize(f)
    idx_beg = 1
    idx_end = MOI.dimension(s.condition)
    con_reif_condition = MOI.add_constraint(
        model,
        MOIU.vectorize(
            [
                one(T) * var_condition,
                f_scalars[idx_beg : idx_end]...,
            ]
        ),
        CP.Reification(s.condition)
    )

    idx_beg = idx_end + 1
    idx_end = idx_end + MOI.dimension(s.true_constraint)
    con_reif_true = MOI.add_constraint(
        model,
        MOIU.vectorize(
            [
                one(T) * var_true,
                f_scalars[idx_beg : idx_end]...,
            ]
        ),
        CP.Reification(s.true_constraint)
    )

    idx_beg = idx_end + 1
    idx_end = idx_end + MOI.dimension(s.false_constraint)
    con_reif_false = MOI.add_constraint(
        model,
        MOIU.vectorize(
            [
                one(T) * var_false,
                f_scalars[idx_beg : idx_end]...,
            ]
        ),
        CP.Reification(s.false_constraint)
    )

    # (x ⟹ y) ∧ (¬x ⟹ z) is equivalent to (¬x ∨ y) ∧ (x ∨ z):
    #     [(1 - var_condition) + var_true] × [var_condition + var_false] ≥ 1
    #     [(1 - var_condition) + var_true ≥ 1] ∧ [var_condition + var_false ≥ 1]
    #     [var_true - var_condition ≥ 0] ∧ [var_condition + var_false ≥ 1]
    #     \____________ if ____________/   \____________ else ___________/
    con_if = MOI.add_constraint(
        model,
        one(T) * var_true - one(T) * var_condition,
        MOI.GreaterThan(zero(T))
    )
    con_else = MOI.add_constraint(
        model,
        one(T) * var_false + one(T) * var_condition,
        MOI.GreaterThan(one(T))
    )

    return IfThenElse2ReificationBridge(var_condition, var_true, var_false, var_condition_bin, var_true_bin, var_false_bin, con_reif_condition, con_reif_true, con_reif_false, con_if, con_else)
end

function MOI.supports_constraint(
    ::Type{IfThenElse2ReificationBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.IfThenElse{S1, S2, S3}},
) where {T, S1, S2, S3}
    return true
    # Ideally, ensure that the underlying solver supports all the needed
    # reified constraints.
end

function MOIB.added_constrained_variable_types(::Type{IfThenElse2ReificationBridge{T}}) where {T}
    return [(MOI.ZeroOne,)]
end

function MOIB.added_constraint_types(::Type{IfThenElse2ReificationBridge{T}}) where {T}
    return [
        (MOI.VectorAffineFunction{T}, CP.Reification), # TODO: how to be more precise?
        (MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T}),
    ]
end

function MOI.get(::IfThenElse2ReificationBridge{T}, ::MOI.NumberOfVariables) where {T}
    return 3
end

function MOI.get(
    ::IfThenElse2ReificationBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T}
    return 3
end

function MOI.get(
    ::IfThenElse2ReificationBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, CP.Reification,
    },
) where {T}
    return 3
end

function MOI.get(
    ::IfThenElse2ReificationBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T},
    },
) where {T}
    return 2
end

function MOI.get(
    b::IfThenElse2ReificationBridge{T},
    ::MOI.ListOfVariableIndices,
) where {T}
    return [b.var_condition, b.var_true, b.var_false]
end

function MOI.get(
    b::IfThenElse2ReificationBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VariableIndex, MOI.ZeroOne,
    },
) where {T}
    return [b.var_condition_bin, b.var_true_bin, b.var_false_bin]
end

function MOI.get(
    b::IfThenElse2ReificationBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T}, CP.Reification,
    },
) where {T}
    return [b.con_reif_condition, b.con_reif_true, b.con_reif_false]
end

function MOI.get(
    b::IfThenElse2ReificationBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, MOI.GreaterThan{T},
    },
) where {T}
    return [b.con_if, b.con_else]
end


text:
```Write a function that takes a list of strings an prints them, one per line,
 in a rectangular frame. For example the list ["Hello", "World", "in", "a", "frame"]
  gets printed as:

*********
* Hello *
* World *
* in    *
* a     *
* frame *
*********```

st = ["Hello","World","in","a","Frame"]

function string_frame(s::Vector{String})
    m_s =  maximum(length.(s))
    println(repeat("*", m_s+4))
    for i in s
        println("\r* $i",repeat(" ",abs(length(i)-m_s))," *")
    end
    println(repeat("*", m_s+4))
end

string_frame(st)

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
const doc = """Naval Fate.

Usage:
  naval_fate.jl ship new <name>...
  naval_fate.jl ship <name> move <x> <y> [--speed=<kn>]
  naval_fate.jl ship shoot <x> <y>
  naval_fate.jl mine (set|remove) <x> <y> [--moored|--drifting]
  naval_fate.jl -h | --help
  naval_fate.jl --version

Options:
  -h --help     Show this screen.
  --version     Show version.
  --speed=<kn>  Speed in knots [default: 10].
  --moored      Moored (anchored) mine.
  --drifting    Drifting mine.

"""

using DocOpt  # import docopt function

arguments = docopt(doc, version=v"2.0.0")
dump(arguments)


text:
abstract type SequenceSpec end

#  All all elements
mutable struct SequenceAll <: SequenceSpec
end

#  None  no elements
mutable struct SequenceNone <: SequenceSpec
end

#  n  elements 1 through n
struct SequenceN{T<:Integer} <: SequenceSpec
    n::T
end

#  UpTo[n] elements 1 up to at most n
struct SequenceUpToN{T<:Integer} <: SequenceSpec
    n::T
end

#  [n]  element n only
struct SequenceNOnly{T<:Integer} <: SequenceSpec
    n::T
end

#  [m,n] elements m through n
struct SequenceMN{T<:Integer,V<:Integer} <: SequenceSpec
    m::V
    n::T
end

#  [m,n] elements m through n, step s
struct SequenceMNS{T<:Integer,V<:Integer,W<:Integer} <: SequenceSpec
    m::V
    n::T
    s::W
end

seqspecstart(s::SequenceN) = 1
seqspecstart(s::SequenceUpToN) = 1
seqspecstart(s::SequenceNOnly) = s.n
seqspecstart(s::SequenceAll) = 1
seqspecstart(s::SequenceMN) = s.m
seqspecstart(s::SequenceMNS) = s.m

seqspecend(s::SequenceN,n) = s.n
seqspecend(s::SequenceUpToN,n) = min(s.n,n)
seqspecend(s::SequenceNOnly,n) = s.n
seqspecend(s::SequenceAll,n) = n
seqspecend(s::SequenceMN,n) = s.n
seqspecend(s::SequenceMNS,n) = s.n

seqspecdi(s::SequenceSpec) = 1
seqspecdi(s::SequenceMNS) = s.s

seqiter(s::SequenceSpec,n) = (seqspecstart(s),seqspecend(s,n),seqspecdi(s))

seqspecerr(x) = symerror(x, " is not a valid sequence specification")

sequencespec(n::Integer) = SequenceN(n)
sequencespec(x::Mxpr{:UpTo}) = SequenceUpToN(x[1])
sequencespec(x::Symbol) = x == :None ? SequenceNone() : x == :All ? SequenceAll() : seqspecerr(x)

sequencespec(n::Integer,nmax) = SequenceN(posnegi(nmax,n))
sequencespec(x::Mxpr{:UpTo},nmax) = SequenceUpToN(posnegi(nmax,n))
sequencespec(x::Symbol,nmax) = x == :None ? SequenceNone() : x == :All ? SequenceAll() : seqspecerr(x)  # hmm maybe we should use nmax here

function sequencespec(x::ListT)
    len = length(x)
    len == 1 && return SequenceNOnly(x[1])
    len == 2 && return SequenceMN(x[1],x[2])
    len == 3 && return SequenceMNS(x[1],x[2],x[3])
    seqspecerr(x)
end

function sequencespec(x::ListT,nmax)
    len = length(x)
    len == 1 && return SequenceNOnly(posnegi(nmax,x[1]))
    len == 2 && return SequenceMN(posnegi(nmax,x[1]),posnegi(nmax,x[2]))
    len == 3 && return SequenceMNS(posnegi(nmax,x[1]),posnegi(nmax,x[2]),x[3])
    seqspecerr(x)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 9947a66a-7fe5-11eb-0955-87881a67fbea
begin
	using PlutoUI
	using Plots

	# plotly()
end

# ╔═╡ 9872b1a8-7fe7-11eb-0b3b-cb6181515182
using SymPy

# ╔═╡ 84fc02d6-7fef-11eb-2e21-cb577c3974f0
# Pkg.add("ForwardDiff)
using ForwardDiff

# ╔═╡ e529f550-7fe3-11eb-06f8-6beaa1506e66
md"""
## AutoDiff - Exploration

Status:
 - ref. 2019-10-26 - ref. [Automatic Differentiation in 10 minutes with Julia](https://www.youtube.com/watch?v=vAp6nUMrKYg)
 - 2021-03-08 moved to `Pluto.jl`

"""

# ╔═╡ 28aa4e9c-7fe4-11eb-3d87-c5a10540c2c1
md"""
Let's start with a simple example, the computation of $\sqrt(x)$ where how autodiff works comes as both a mathematical surprise, and a computing wonder.
The example is the Babylonian algorithm, known to mankind for millenia, to compute $\sqrt(x)$:

repeat

   $~~~~t \leftarrow \frac{(t + \frac{x}{t})}{2}$

until t converges to $\sqrt(x)$

Each iteration has one add and two divides.
For illustration purposes, 10 iterations suffice
"""

# ╔═╡ 1ecd40bc-7fe6-11eb-0a2f-cb9964fc5dbf
md"""
### And now the derivative, almost by magic

In a few lines of code. No mention of $\frac{1}{2}$ over  $\sqrt(x)$. We will use the "dual number" denoted as $D$ in what follows, those where invented by the famous algebraist Clifford in 1873.
"""

# ╔═╡ 4ac7200c-7fe6-11eb-259a-8d881d2f7e3e
struct D <: Number # D is a <function, derivative> pair (a pair of floats)
    f::Tuple{Float64, Float64}
end

# ╔═╡ 54f29be2-7fe6-11eb-38c5-d163b159730a
md"""
- Sum Rule: $(x + y)' = x' + y'$
- Quotient Rule: $(\frac{x}{y})' = \frac{yx' - xy'}{y^2}$
"""

# ╔═╡ 6e7efe02-7fe6-11eb-19a7-8b72c0836c8b
begin
	import Base: +, /, -, *, convert, promote_rule

	## overload:
	+(x::D, y::D) = D(x.f .+ y.f)
	/(x::D, y::D) = D((x.f[1] / y.f[1],
			(y.f[1] * x.f[2] - x.f[1] * y.f[2]) / y.f[1]^2))
	-(x::D, y::D) = D(x.f .- y.f)
	*(x::D, y::D) = D((x.f[1] * y.f[1], (y.f[1] * x.f[2] + x.f[1] * y.f[2])))

	## convert ordinary number to Dual number, intro. 0 for derivative
	convert(::Type{D}, x::Real) = D((x, zero(x)))

	## then promote...
	promote_rule(::Type{D}, ::Type{<:Number}) = D
end

# ╔═╡ 7a9389a2-7fe5-11eb-06b8-c52b073c876a
function babylonian(x; n = 10)
    t = (1. + x) / 2.
    for i ∈ 2:n
        t = (t + x/t) / 2.
    end
    t
end

# ╔═╡ 9091553e-7fe5-11eb-2131-9b1b1097c932
begin
	α = π
	babylonian(α), √α
end

# ╔═╡ b9ccd308-7fe5-11eb-0b0c-efa652660cd5
with_terminal() do
	for α ∈ 2:3
   		println(babylonian(α), " ", √α)
	end
end

# ╔═╡ ff570916-7fe5-11eb-35ab-e915308e7bfc
## WARN: first plots require to load package which takes time

begin
	ix = 0:.01:49

	plot([x -> babylonian(x, n=i) for i ∈ 1:5],
    	ix,
    	label=["Iteration $jx" for _ ∈ 1:1, jx ∈ 1:5])

	plot!(sqrt,
    	ix,
    	c="darkslateblue",
    	label="sqrt",
    	title = "Those Babylonian really knew about √")
end

# ╔═╡ 732679da-7ff0-11eb-1145-3750276a9601
begin
	using LinearAlgebra

	n = 4
	Strang = SymTridiagonal(2 * ones(n), - ones(n-1))
end

# ╔═╡ 6e64c064-7fe6-11eb-1e6d-3dfeebace21b
md"""
The same algorithm with no rewrite at all computes properly the derivative as the check shows:
"""

# ╔═╡ 6e487e0e-7fe6-11eb-23a2-35a3a2ed7015
x₁ = 49; babylonian(D((x₁, 1))), (√x₁, .5 / √x₁)

# ╔═╡ a362fc86-7fe6-11eb-1b7a-03cd0549e790
x₂ = π; babylonian(D((x₂, 1))), (√x₂, .5 / √x₂)

# ╔═╡ a348cb0e-7fe6-11eb-1c9d-4b4278a367aa
md"""
#### It just works!

How does this work?
We will explain in a moment. Right now marvel that it does. Note we did not import any autodiff package. Everything is just basic Julia.

#### The assembler

Most folks don't read assembler, but one can see that it is short (minus comments). The shortness is a clue that suggests speed!

"""

# ╔═╡ a32df540-7fe6-11eb-073f-1b34489aa4a0
with_terminal() do
	@show @inline function babylonian_(x; n = 10)
    	 t = (1 + x) / 2.
     	for i = 2:n; t = (t + x/t) / 2. end
     	t
 	end
end

# ╔═╡ a312dfda-7fe6-11eb-0bae-37587f221888
with_terminal() do
	@code_native babylonian(D((2, 1)))
end

# ╔═╡ 37574936-7fe7-11eb-132a-2df44bf8d6bb
md"""
### Symbolically

We haven't yet explained how it works, but it may be of some value to understand that the below is mathematically equivalent, though not what the computation is doing.
Notice in the below that babylonian works on SymPy Symbols.

Note: Python and Julia are very good friends.It's not a competition! Watch how nicely we can use the same code now in SymPy.

"""

# ╔═╡ 98544e52-7fe7-11eb-1094-9382f7126ce4
with_terminal()do
	x = symbols("x")

	println("Iterations as a function of x")
	for k ∈ 1:5
    	println(simplify(babylonian(x, n=k)))
	end
end

# ╔═╡ 98399af8-7fe7-11eb-1276-13c05fa36ed9
with_terminal() do
	x = symbols("x")

	println("Derivatives as a function of x")

	for k ∈ 1:5
    	println(simplify(diff(simplify(babylonian(x, n=k)), x)))
	end
end

# ╔═╡ 981ea264-7fe7-11eb-06ef-97ee66d5619b
md"""
Let's by hand take the "derivative" of the babylonian iteration with respect to x.

Specifically $t' = \frac{dt}{dx}$. This is the old fashioned way of a human writing code.
"""

# ╔═╡ d05af69e-7feb-11eb-3874-d1b05fb6434a
function dbabylonian(x; n = 10)
    t = (1. + x) / 2.
    dt = 1. / 2.

    for i ∈ 2:n
        t = (t + x/t) / 2.
        dt = (dt + (t - x * dt) / t^2) / 2.
    end
    dt
end

# ╔═╡ d0310bf4-7feb-11eb-2d05-c9531128eb4a
md"""

Note:
  -  $t = \frac{1}{2} \times (t + \frac{x}{t})$

  - then $(\frac{dt}{dx})' = \frac{1}{2} \times (t' + \frac{x' \times t - x \times t'}{t^2})$, as $x' = 1$, we get:

$(\frac{dt}{dx})' = \frac{1}{2} \times (\frac{t' + (t - x \times t')}{t^2})$
"""

# ╔═╡ d01a58be-7feb-11eb-1f29-6d23db7914da
x= π; dbabylonian(x), .5 / √x

# ╔═╡ cffee818-7feb-11eb-079b-a7dbfdd8cdd6
md"""
What just happened?

Answer: we created an iteration by hand for t' given our iteration for t, Then we ran the iteration alongside the iteration for t.
"""

# ╔═╡ d0ca45ac-7fec-11eb-202b-87ae22b66616
babylonian(D((x, 1)))

# ╔═╡ d0addd04-7fec-11eb-0213-53ef5d729159
md"""

How did this work?

It created the same derivative iteration that we did by hand, using very general rules that are set once and need not be written by hand (and multiple dispatch).

Important: the derivative is substituted before the JIT compiler, and thus efficient compiled code is executed.
"""

# ╔═╡ d092a21e-7fec-11eb-0cc1-d1515675f521
md"""
### Dual Number Notation

Instead of $D(a, b)$ we can write: $a + b \epsilon$ where $\epsilon$ satisfies $\epsilon^2 = 0$. Some people like to recall imaginary numbers where an $i$ is introduced with $i^2 = -1$. Others like to think of how engineers just fdrop the $O(\epsilon^2)$ terms.

The four rules are:
  - 1 & 2. $(a + b\epsilon) \pm (c + d\epsilon) = (a + c) \pm (b + d)\epsilon$
  - 3.$(a + b\epsilon) \times (c + d\epsilon) = (ac) + (bc + ad)\epsilon$
  - 4.$\frac{(a + b\epsilon)}{(c + d\epsilon)} = (\frac{a}{c}) + \frac{(bc - ad)}{d^2}\epsilon$

"""

# ╔═╡ 1818fb9c-7fed-11eb-0bd3-21411f08fcdf
Base.show(io::IO, x::D) = print(io, x.f[1], " + ", x.f[2], " ϵ")

# ╔═╡ 14125e9c-7fef-11eb-3dc4-a75f84383d5d
D((1, 0))

# ╔═╡ 13e452cc-7fef-11eb-2e1d-a5e3b9936a4c
D((0, 2))^2 # should be zero!

# ╔═╡ 13ab59ea-7fef-11eb-2078-6d8ee11d1c90
D((2, 1))^2

# ╔═╡ 2886e6c4-7fef-11eb-0fff-9b3033e4f943
begin
	ϵ = D((0, 1))
	with_terminal() do
		@code_native(ϵ^2)
	end
end

# ╔═╡ 3f1f4bba-7fef-11eb-2f39-cf906aec10ae
ϵ * ϵ

# ╔═╡ 3f037186-7fef-11eb-200f-3579d0ab8e1f
ϵ^2

# ╔═╡ 3eea0d36-7fef-11eb-3ede-c7d54243194c
(1 + ϵ)^5 # note it just works (we did not train powers)

# ╔═╡ 3ece19f8-7fef-11eb-3b91-ff25fd3ca8f4
md"""
### Generalization to arbitrary roots
"""

# ╔═╡ 3eb0843a-7fef-11eb-2c67-9ba2e4d1eefc
function nth_root(x, n=2; t=1, p=10)
    for i = 1:p
        t += (x / t^(n-1) - t) / n
    end
    t
end

# ╔═╡ 5c9fbbc6-7fef-11eb-32ad-b986ca56aa9b
nth_root(2, 3), ∛2  # copied from https://www.alt-codes.net/root-symbols

# ╔═╡ 5c84a84c-7fef-11eb-3117-01f1842890f5
nth_root(2 + ϵ, 3)

# ╔═╡ 5c665cd4-7fef-11eb-2c89-e153a973b30e
nth_root(7, 12), 7^(1/12)

# ╔═╡ 5c4f64ac-7fef-11eb-1864-13f293e85fc4
let x = 2.0
	nth_root( x + ϵ), ∛x, 1/x*(2/3)/3
end

# ╔═╡ 5c33d6ce-7fef-11eb-12c1-29396395358d
md"""
### Forward Diff

Now that you understand it, you can use the official package.
"""

# ╔═╡ 84ddd280-7fef-11eb-29cb-3fc8c8547134
ForwardDiff.derivative(sqrt, 2)

# ╔═╡ 916346b8-7fef-11eb-0d53-29bd11c34510
ForwardDiff.derivative(babylonian, 2)

# ╔═╡ 99672fe6-7fef-11eb-305e-d98ce4bc8dc1
with_terminal() do
	println(@which ForwardDiff.derivative(√, 2))
end

# ╔═╡ be076e8a-7fef-11eb-0741-79beea838677
md"""
### Close Look at Convergence with big floats

the $-log10$ gives the number of correct digits. Watch the quadratic convergence right before your eyes.
"""

# ╔═╡ d5b96cac-7fef-11eb-0772-551e4d96efb2
with_terminal() do
	setprecision(3000)

	# round.(Float64.(log10.([babylonian(BigFloat(2), n=k) for k=1:10] - √BigFloat(2))), 3)
	println(round.(Float64.(log10.([babylonian(BigFloat(2), n=k) - √BigFloat(2.) for k=1:10])), sigdigits=6))
end

# ╔═╡ d59d0508-7fef-11eb-20a3-3f108b859670
struct D1{T} <: Number # D is a <function, derivative> pair
    f::Tuple{T, T}
end

# ╔═╡ d57e75ac-7fef-11eb-1cc7-db506a7e4e2f
begin
	z = D((2., 1.))
	z₁ = D1((BigFloat(2.), BigFloat(1.)))
end

# ╔═╡ fbb04502-7fef-11eb-28fd-414c2f982855
# begin
# 	#import Base: +, /, convert, promote_rule

# 	## overload:
# 	+(x::D1, y::D1) = D1(x.f .+ y.f)
# 	/(x::D1, y::D1) = D1((x.f[1] / y.f[1], (y.f[1] * x.f[2] - x.f[1] * y.f[2]) / y.f[1]^2))

# 	convert(::Type{D1{T}}, x::Real) where {T} = D1((convert(T, x), zero(T)))
# 	promote_rule(::Type{D1{T}}, ::Type{S}) where {T, S <: Number} = D1{promote_type(T, S)}
# end

# ╔═╡ 3e2da118-7ff0-11eb-36ca-e32af03b2a0c
A = randn(3, 3)

# ╔═╡ 3e0a484e-7ff0-11eb-17be-e9e7aa66fb2c
x₃= randn(3)

# ╔═╡ 3dd5c57e-7ff0-11eb-3330-cf7d02b5c632
ForwardDiff.gradient(x -> x'A*x, x₃)

# ╔═╡ Cell order:
# ╟─e529f550-7fe3-11eb-06f8-6beaa1506e66
# ╠═9947a66a-7fe5-11eb-0955-87881a67fbea
# ╟─28aa4e9c-7fe4-11eb-3d87-c5a10540c2c1
# ╠═7a9389a2-7fe5-11eb-06b8-c52b073c876a
# ╠═9091553e-7fe5-11eb-2131-9b1b1097c932
# ╠═b9ccd308-7fe5-11eb-0b0c-efa652660cd5
# ╠═ff570916-7fe5-11eb-35ab-e915308e7bfc
# ╟─1ecd40bc-7fe6-11eb-0a2f-cb9964fc5dbf
# ╠═4ac7200c-7fe6-11eb-259a-8d881d2f7e3e
# ╟─54f29be2-7fe6-11eb-38c5-d163b159730a
# ╠═6e7efe02-7fe6-11eb-19a7-8b72c0836c8b
# ╟─6e64c064-7fe6-11eb-1e6d-3dfeebace21b
# ╠═6e487e0e-7fe6-11eb-23a2-35a3a2ed7015
# ╠═a362fc86-7fe6-11eb-1b7a-03cd0549e790
# ╟─a348cb0e-7fe6-11eb-1c9d-4b4278a367aa
# ╠═a32df540-7fe6-11eb-073f-1b34489aa4a0
# ╠═a312dfda-7fe6-11eb-0bae-37587f221888
# ╟─37574936-7fe7-11eb-132a-2df44bf8d6bb
# ╠═9872b1a8-7fe7-11eb-0b3b-cb6181515182
# ╠═98544e52-7fe7-11eb-1094-9382f7126ce4
# ╠═98399af8-7fe7-11eb-1276-13c05fa36ed9
# ╟─981ea264-7fe7-11eb-06ef-97ee66d5619b
# ╠═d05af69e-7feb-11eb-3874-d1b05fb6434a
# ╟─d0310bf4-7feb-11eb-2d05-c9531128eb4a
# ╠═d01a58be-7feb-11eb-1f29-6d23db7914da
# ╟─cffee818-7feb-11eb-079b-a7dbfdd8cdd6
# ╠═d0ca45ac-7fec-11eb-202b-87ae22b66616
# ╟─d0addd04-7fec-11eb-0213-53ef5d729159
# ╟─d092a21e-7fec-11eb-0cc1-d1515675f521
# ╠═1818fb9c-7fed-11eb-0bd3-21411f08fcdf
# ╠═14125e9c-7fef-11eb-3dc4-a75f84383d5d
# ╠═13e452cc-7fef-11eb-2e1d-a5e3b9936a4c
# ╠═13ab59ea-7fef-11eb-2078-6d8ee11d1c90
# ╠═2886e6c4-7fef-11eb-0fff-9b3033e4f943
# ╠═3f1f4bba-7fef-11eb-2f39-cf906aec10ae
# ╠═3f037186-7fef-11eb-200f-3579d0ab8e1f
# ╠═3eea0d36-7fef-11eb-3ede-c7d54243194c
# ╟─3ece19f8-7fef-11eb-3b91-ff25fd3ca8f4
# ╠═3eb0843a-7fef-11eb-2c67-9ba2e4d1eefc
# ╠═5c9fbbc6-7fef-11eb-32ad-b986ca56aa9b
# ╠═5c84a84c-7fef-11eb-3117-01f1842890f5
# ╠═5c665cd4-7fef-11eb-2c89-e153a973b30e
# ╠═5c4f64ac-7fef-11eb-1864-13f293e85fc4
# ╟─5c33d6ce-7fef-11eb-12c1-29396395358d
# ╠═84fc02d6-7fef-11eb-2e21-cb577c3974f0
# ╠═84ddd280-7fef-11eb-29cb-3fc8c8547134
# ╠═916346b8-7fef-11eb-0d53-29bd11c34510
# ╠═99672fe6-7fef-11eb-305e-d98ce4bc8dc1
# ╟─be076e8a-7fef-11eb-0741-79beea838677
# ╠═d5b96cac-7fef-11eb-0772-551e4d96efb2
# ╠═d59d0508-7fef-11eb-20a3-3f108b859670
# ╠═d57e75ac-7fef-11eb-1cc7-db506a7e4e2f
# ╠═fbb04502-7fef-11eb-28fd-414c2f982855
# ╠═3e2da118-7ff0-11eb-36ca-e32af03b2a0c
# ╠═3e0a484e-7ff0-11eb-17be-e9e7aa66fb2c
# ╠═3dd5c57e-7ff0-11eb-3330-cf7d02b5c632
# ╠═732679da-7ff0-11eb-1145-3750276a9601

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_name: pascal-p/julia-notebooks


text:
function companion(A)
    N = size(A[1], 1)
    L = size(A, 1)-1
    C1 = zeros(ComplexF64, N*L, N*L)
    C2 = zeros(ComplexF64, N*L, N*L)
    res = zeros(N*L)
    C1[1:N, 1:N] .= A[1]
    for i=N+1:N*L
        C1[i,i] = 1
        C2[i,i-N] = 1
    end
    for i=1:L
        C2[1:N, N*(i-1)+1:N*i] .= -A[i+1]
    end
    F = eigen!(C1, C2)
    X = F.vectors[(L-1)*N+1:L*N, :]
    Λ = F.values

    for i=1:N*L
        X[:,i] ./= norm(X[:,i])
        R = zeros(ComplexF64, N)
        for j=1:L+1
            R += A[j] * X[:,i] .* Λ[i]^(j-1)
        end
        res[i] = norm(R) / norm(sum(A[j] .* Λ[i]^(j-1) for j=1:L+1))
    end
    return Λ, X, res
end


text:
function feedback(c::Int,ntot::Int,nbit::Int,operator)
    c += 1
    if c % 5369200 == 0
        if nbit == 16
            percent = Int(round(100c / ntot))
            print("\r\u1b[K")
            if percent == 100
                println("Sonum16: Table for ($operator) 100%.")
            else
                print("Sonum16: Table for ($operator) $percent%")
            end
        end
    end
    return c
end


text:
module Utils

include("Model.jl")
import .Model:
  model_sample

using Printf
using Images:
  Gray,
  save

using BSON:
  @save,
  @load

function img(x)
  Gray.(rotate(reshape(x, 28, 28)))
end

function gen_images(outfile, f)
  sample = hcat(img.([model_sample(f) for i = 1:10])...)
  save(outfile, sample)
end

function gen_images(
  outfile :: String,
  g,
  f,
  X :: T
) where {T <: AbstractArray}
  sample = hcat(img.([model_sample(g, f, X[:,:,:,i]).data for i = 1:size(X,4)])...)
  safe_img_save(outfile, sample)
end

function safe_img_save(x, outfile)
  try
    save(x, outfile)
  catch e
    @error @sprintf("failed write: %s", outfile) e
  end
end

function rotate(ximg :: T) where {T <: AbstractArray}
  xcopy = zeros(size(ximg)...)
  for i in range(1, stop = size(ximg, 1))
    for j in  range(1, stop = size(ximg, 2))
      xcopy[i,j] = ximg[j,i]
    end
  end
  xcopy
end

function save_model(f, g, outfile)
  model_pair = (f,g)
  try
    @save outfile model_pair
  catch e
    @error @sprintf("failed write model to: %s", outfile) e
  end
end

function load_model(infile)
  @load infile model_pair
  f,g = model_pair
  return f, g
end

end

max_stars_repo_name: aniketjivani/variational-autoencoders

max_issues_repo_name: aniketjivani/variational-autoencoders

max_forks_repo_name: aniketjivani/variational-autoencoders


text:
@inline function add_ddfp_dd(x::Tuple{T,T}, y::T) where {T<:IEEEFloat}
    hi, lo = x
    hi, lo = two_sum(hi, y, lo)
    return hi, lo
end

@inline function sub_ddfp_dd(x::Tuple{T,T}, y::T) where {T<:IEEEFloat}
    hi, lo = x
    hi, lo = two_sum(hi, lo, -y)
    return hi, lo
end

#=
@inline function mul_ddfp_dd(x::Tuple{T,T}, y::T) where {T<:IEEEFloat}
    hi, lo = x
    hihi, hilo = two_prod(y, hi)
    lohi, lolo = two_prod(y, lo)
    hi, lo = two_sum(hihi, hilo, lohi, lolo)
    return hi, lo
end
=#

# Algorithm 12 from Tight and rigourous error bounds.  relative error <= 5u²
@inline function mul_ddfp_dd(x::Tuple{T,T}, y::T) where T<:IEEEFloat
    xhi, xlo = x
    hi, lo = two_prod(xhi, y)
    t = xlo * y
    t = t + lo
    hi, lo = two_hilo_sum(hi, t)
    return hi, lo
end

@inline function dvi_ddfp_dd(x::Tuple{T,T}, y::T) where {T<:IEEEFloat}
    xhi, xlo = x
    hi = xhi / y
    uh, ul = two_prod(hi, y)
    lo = ((((xhi - uh) - ul) + xlo))/y
    hi,lo = two_hilo_sum(hi, lo)
    return hi, lo
end


text:
using Test

@testset "AbstractWrappedDataFrame" begin
    include("basic.jl")
end


text:
"""
A product map is diagonal and acts on each of the components of x separately:
`y = f(x)` becomes `y_i = f_i(x_i)`.
"""
abstract type ProductMap{T} <: CompositeLazyMap{T} end

components(m::ProductMap) = m.maps

VcatMapElement = Union{Map{<:SVector},Map{<:Number}}

ProductMap(maps::Tuple) = ProductMap(maps...)
ProductMap(maps::SVector) = ProductMap(maps...)
ProductMap(maps...) = TupleProductMap(maps...)
ProductMap(maps::VcatMapElement...) = VcatMap(maps...)
ProductMap(maps::AbstractVector) = VectorProductMap(maps)

ProductMap{T}(maps...) where {T} = _TypedProductMap(T, maps...)
_TypedProductMap(::Type{T}, maps...) where {T<:Tuple} = TupleProductMap(maps...)
_TypedProductMap(::Type{SVector{N,T}}, maps...) where {N,T} = VcatMap{N,T}(maps...)
_TypedProductMap(::Type{T}, maps...) where {T<:AbstractVector} = VectorProductMap{T}(maps...)

compatibleproductdims(d1::ProductMap, d2::ProductMap) =
	mapsize(d1) == mapsize(d2) &&
		all(map(==, map(mapsize, components(d1)), map(mapsize, components(d2))))

compatibleproduct(d1::ProductMap, d2::ProductMap) =
	compatibleproductdims(d1, d2) && compatible_domaintype(d1, d2)

isconstant(m::ProductMap) = mapreduce(isconstant, &, components(m))
islinear(m::ProductMap) = mapreduce(islinear, &, components(m))
isaffine(m::ProductMap) = mapreduce(isaffine, &, components(m))

matrix(m::ProductMap) = toexternalmatrix(m, map(matrix, components(m)))
vector(m::ProductMap) = toexternalpoint(m, map(vector, components(m)))
constant(m::ProductMap) = toexternalpoint(m, map(constant, components(m)))

jacobian(m::ProductMap, x) =
	toexternalmatrix(m, map(jacobian, components(m), tointernalpoint(m, x)))
function jacobian(m::ProductMap{T}) where {T}
	if isaffine(m)
		ConstantMap{T}(matrix(m))
	else
		ProductMap(map(jacobian, components(m)))
	end
end
# function jacdet(m::ProductMap, x)

similarmap(m::ProductMap, ::Type{T}) where {T} = ProductMap{T}(components(m))

tointernalpoint(m::ProductMap, x) = x
toexternalpoint(m::ProductMap, y) = y

applymap(m::ProductMap, x) =
	toexternalpoint(m, map(applymap, components(m), tointernalpoint(m, x)))

productmap(map1, map2) = productmap1(map1, map2)
productmap1(map1, map2) = productmap2(map1, map2)
productmap2(map1, map2) = ProductMap(map1, map2)
productmap(map1::ProductMap, map2::ProductMap) =
	ProductMap(components(map1)..., components(map2)...)
productmap1(map1::ProductMap, map2) = ProductMap(components(map1)..., map2)
productmap2(map1, map2::ProductMap) = ProductMap(map1, components(map2)...)

for op in (:inverse, :leftinverse, :rightinverse)
    @eval $op(m::ProductMap) = ProductMap(map($op, components(m)))
	@eval $op(m::ProductMap, x) = toexternalpoint(m, map($op, components(m), tointernalpoint(m, x)))
end

function composedmap(m1::ProductMap, m2::ProductMap)
	if compatibleproduct(m1, m2)
		ProductMap(map(composedmap, components(m1), components(m2)))
	else
		Composition(m1,m2)
	end
end

mapsize(m::ProductMap) = (sum(t->mapsize(t,1), components(m)), sum(t->mapsize(t,2), components(m)))

==(m1::ProductMap, m2::ProductMap) = all(map(isequal, components(m1), components(m2)))

Display.combinationsymbol(m::ProductMap) = Display.Symbol('⊗')
Display.displaystencil(m::ProductMap) = composite_displaystencil(m)
show(io::IO, mime::MIME"text/plain", m::ProductMap) = composite_show(io, mime, m)
show(io::IO, m::ProductMap) = composite_show_compact(io, m)

"""
A `VcatMap` is a product map with domain and codomain vectors
concatenated (`vcat`) into a single vector.
"""
struct VcatMap{N,T,DIM,MAPS} <: ProductMap{SVector{N,T}}
    maps    ::  MAPS
end

VcatMap(maps::Union{Tuple,Vector}) = VcatMap(maps...)
function VcatMap(maps...)
	T = numtype(maps...)
	M = sum(t->mapsize(t,1), maps)
	N = sum(t->mapsize(t,2), maps)
	# M,N = reduce((x,y) -> (x[1]+y[1],x[2]+y[2]), map(mapsize,maps))
	@assert M==N
	VcatMap{N,T}(maps...)
end

mapdim(map) = mapsize(map,2)

VcatMap{N,T}(maps::Union{Tuple,Vector}) where {N,T} = VcatMap{N,T}(maps...)
VcatMap{N,T}(maps...) where {N,T} = VcatMap{N,T}(map(Map{T}, maps)...)
function VcatMap{N,T}(maps::Map{T}...) where {N,T}
	DIM = map(mapdim,maps)
	VcatMap{N,T,DIM}(maps...)
end

VcatMap{N,T,DIM}(maps...) where {N,T,DIM} = VcatMap{N,T,DIM,typeof(maps)}(maps)

mapsize(m::VcatMap{N}) where {N} = (N,N)

tointernalpoint(m::VcatMap{N,T,DIM}, x) where {N,T,DIM} =
	convert_fromcartesian(x, Val{DIM}())
toexternalpoint(m::VcatMap{N,T,DIM}, y) where {N,T,DIM} =
	convert_tocartesian(y, Val{DIM}())

size_as_matrix(A::AbstractArray) = size(A)
size_as_matrix(A::Number) = (1,1)

# The Jacobian is block-diagonal
function toexternalmatrix(m::VcatMap{N,T}, matrices) where {N,T}
	A = zeros(T, N, N)
	l = 0
	for el in matrices
		m,n = size_as_matrix(el)
		@assert m==n
		A[l+1:l+m,l+1:l+n] .= el
		l += n
	end
	SMatrix{N,N}(A)
end

"""
A `VectorProductMap` is a product map where all components are univariate maps,
with inputs and outputs collected into a `Vector`.
"""
struct VectorProductMap{T<:AbstractVector,M} <: ProductMap{T}
    maps    ::  Vector{M}
end

VectorProductMap(maps::AbstractMap...) = VectorProductMap(maps)
VectorProductMap(maps) = VectorProductMap(collect(maps))
function VectorProductMap(maps::Vector)
	T = mapreduce(numtype, promote_type, maps)
	VectorProductMap{Vector{T}}(maps)
end

VectorProductMap{T}(maps::AbstractMap...) where {T} = VectorProductMap{T}(maps)
VectorProductMap{T}(maps) where {T} = VectorProductMap{T}(collect(maps))
function VectorProductMap{T}(maps::Vector) where {T}
	Tmaps = convert.(Map{eltype(T)}, maps)
	VectorProductMap{T,eltype(Tmaps)}(Tmaps)
end

# the Jacobian is a diagonal matrix
toexternalmatrix(m::VectorProductMap, matrices) = Diagonal(matrices)

dimension(m::VectorProductMap) = length(m.maps)

"""
A `TupleProductMap` is a product map with all components collected in a tuple.
There is no vector-valued function associated with this map.
"""
struct TupleProductMap{T,MM} <: ProductMap{T}
    maps    ::  MM
end

TupleProductMap(maps::Vector) = TupleProductMap(maps...)
TupleProductMap(maps...) = TupleProductMap(maps)
function TupleProductMap(maps::Tuple)
	T = Tuple{map(eltype, maps)...}
	TupleProductMap{T}(maps)
end

TupleProductMap{T}(maps::Vector) where {T} = TupleProductMap{T}(maps...)
TupleProductMap{T}(maps...) where {T} = TupleProductMap{T}(maps)
function TupleProductMap{T}(maps::NTuple{N,<:AbstractMap}) where {N,T <: Tuple}
	Tmaps = map((t,d) -> convert(Map{t},d), tuple(T.parameters...), maps)
	TupleProductMap{T,typeof(Tmaps)}(Tmaps)
end
TupleProductMap{T}(maps) where {T <: Tuple} = TupleProductMap{T,typeof(maps)}(maps)


text:
using Distributions, Random

#Simulates with the call gapping policy
function simCallGapping(N,τ,γ,c,tMax; seed = 1974)
    Random.seed!(seed)
    state = [0 for _ in 1:N]
    currentArm = 1
    allowedToSwitch = false
    t = 0.0
    reward = 0.0
    λ, μ = 1.0, 1/γ - 1
    rλ() = rand(Exponential(1/λ))
    rμ() = rand(Exponential(1/μ))
    nextArm(i) = (i%N)+1

    #N+1 event types in this discrete event simulation
    timeJumps = [rλ() for _ in 1:N]
    push!(timeJumps,τ) #The last one is for the call gapping

    while t<tMax
       tLast = t
       t,i = findmin(timeJumps)
       reward += state[currentArm]*(t-tLast)

       if i == N+1 #call gapping timer up
            #println("$(t),$(currentArm): gap reached")
            if state[currentArm] == 0
                #println("$(t),$(currentArm): switching from $(currentArm)")
                currentArm = nextArm(currentArm)
                allowedToSwitch = false
                timeJumps[N+1] = t+τ
                reward -= c
            else
                allowedToSwitch = true
                timeJumps[N+1] = Inf
            end
        else #arm changing state
           if state[i] == 0
              #println("$(t),$(currentArm): Arm $(i) switch to 1")
              state[i] = 1
              timeJumps[i] = t+rμ()
           else #state was 1
               #println("$(t),$(currentArm): Arm $(i) switch to 0")
               state[i] = 0
               timeJumps[i] = t+rλ()
               if allowedToSwitch && state[currentArm] == 0
                    #println("$(t),$(currentArm): switching from $(currentArm)")
                    currentArm = nextArm(currentArm)
                    allowedToSwitch = false
                    timeJumps[N+1] = t+τ
                    reward -= c
               end
           end
       end
    end
    reward/tMax
end


text:
module VortexLattice

using LinearAlgebra
using StaticArrays
using Interpolations
using WriteVTK

# value for dimensionalizing, included just for clarity in the algorithms
const RHO = 1.0

include("panel.jl")
export SurfacePanel, WakePanel, TrefftzPanel
export reflect, set_normal

include("wake.jl")
export Wake

include("geometry.jl")
export AbstractSpacing, Uniform, Sine, Cosine
export grid_to_surface_panels, wing_to_surface_panels
export lifting_line_geometry, lifting_line_geometry!
export translate, translate!, rotate, rotate!, reflect

include("reference.jl")
export Reference
export AbstractFrame, Body, Stability, Wind

include("freestream.jl")
export Freestream, trajectory_to_freestream

include("induced.jl")

include("circulation.jl")

include("system.jl")
export System
export PanelProperties, get_surface_properties

include("analyses.jl")
export steady_analysis, steady_analysis!
export unsteady_analysis, unsteady_analysis!
export spanwise_force_coefficients

include("nearfield.jl")
export body_forces
export body_forces_history
export lifting_line_coefficients, lifting_line_coefficients!

include("farfield.jl")
export far_field_drag

include("stability.jl")
export body_derivatives, stability_derivatives

include("visualization.jl")
export write_vtk

end # module


text:
# Skewness
"""
    index_skew(ts::TimeArray)

Returns the skewness of ts.
"""
function index_skew(ts::TimeArray)
    x = values(ts)
    μ = mean(x)
    σ = std(x)

    z = (x .- μ) ./ σ

    return mean(z .^ 3)
end

"""
    index_Rx(ts::TimeArray, x::Real)

Returns the relative frequency of days with precipitation above x mm.
"""
function index_Rx(ts::TimeArray, x::Real)
    return sum(values(ts .> x)) / length(ts)
end

# R10
"""
    index_R10(ts::TimeArray)

Returns the relative frequency of days with precipitation above 10 mm.
"""
function index_R10(ts::TimeArray)
    return index_Rx(ts, 10)
end

"""
    index_pXwet(ts::TimeArray, p::Real; wet_thresh::Real=1.0)

Returns the p percentile of wet days considering days with precipitation above wet_thresh (mm).
"""
function index_pXwet(ts::TimeArray, p::Real; wet_thresh::Real=1.0)
    @assert zero(p)<p<one(p) "the quantile level should be between 0 and 1."

    x = values(ts)
    x = x[x .>= wet_thresh]

    return quantile(x, p)
end

# P98Wet
"""
    index_p98wet(ts::TimeArray)

Returns the 98% percentile of wet days considering days with precipitation above 1 mm.
"""
function index_p98wet(ts::TimeArray)
    return index_pXwet(ts, 0.98)
end

# P90Wet
"""
    index_p90wet(ts::TimeArray)

Returns the 90% percentile of wet days considering days with precipitation above 1 mm.
"""
function index_p90wet(ts::TimeArray)
    return index_pXwet(ts, 0.9)
end

"""
    index_pXwetamount(ts::TimeArray, p::Real; wet_thresh=1.0)

Returns the precipitation amount from days above p percentile.
"""
function index_pXwetamount(ts::TimeArray, p::Real; wet_thresh=1.0)
    @assert zero(p)<p<one(p) "the quantile level should be between 0 and 1."

    thresh = index_pXwet(ts, p, wet_thresh=wet_thresh)

    x = values(ts)
    x = x[x .> thresh]

    return sum(x) / length(x)
end

# P98WetAmount
"""
    index_p98wetamount(ts::TimeArray)

Returns the precipitation amount from days above 98% percentile.
"""
function index_p98wetamount(ts::TimeArray)
    return index_pXwetamount(ts, 0.98)
end

"""
    index_RVxmax(ts::TimeArray, p::Real)

Returns the p-years return value of heavy precipitation.
"""
function index_RVxmax(ts::TimeArray, p::Real)

    x = collapse(ts,year,maximum)

    fm = Extremes.gevfit(values(x))
    r = Extremes.returnlevel(fm, p)

    return r.value[1]
end

# RV20_max
"""
    index_RV20max(ts::TimeArray)

Returns the 20-years return value of heavy precipitation.
"""
function index_RV20max(ts::TimeArray)
    return index_RVxmax(ts, 20.0)
end

# SDII
"""
    index_SDII(ts::TimeArray)

Returns the precipitation amount in wet days divided by the number of wet days.
"""
function index_SDII(ts::TimeArray; thresh::Real=1.0)

    x = values(ts)
    idx = x .> thresh

    return sum(x[idx]) / length(x[idx])
end


text:
module CachedArraysTests

using Test
using CellwiseValues

x = rand(10,5)

y = CachedArray(x)

@test y == x

z = y + x

a = y.array

setsize!(y,(5,3))

b = x[1:5,1:3]

@test y.array === a

@test size(y) == size(b)

@test y == b

setsize!(y,(15,30))

@test ! ( y.array === a )

@test size(y) == size(y.array)

y = CachedArray(Float64,2)

@test isa(y, AbstractArray{Float64,2})

@test isa(y, CachedArray{Float64,2,Array{Float64,2}})

@test size(y) == (0,0)

setsize!(y,(15,30))

@test size(y) == (15,30)

@test isa(CachedVector(Int),CachedVector{Int,Vector{Int}})

@test isa(CachedMatrix(Int),CachedMatrix{Int,Matrix{Int}})

end # module CachedArraysTests


text:
# functions for evaluation

"""
    genetic_evaluation!(pop::PTPopulation; method="blup", repeated=true, rel=0.5, verbose=false, updategebv=false, cg=false)

Perform genetic/genomic evaluation using `model`.
Put `repeated=true` if PE is considered in the model
This function supports the following models.

- `model="blup"`: traditional pedigree-based BLUP; stored as `ebv`
- `model="tbv"`: true breeding value (TBV) with reliability `rel` (defined as the squared correlation between TBV and predictions); stored as `gebv`
- `model="tbv_ms"`: Mendelian sampling of true breeding value (TBV) with reliability `rel` (defined as the squared correlation between TBV and predictions); stored as `gebv`

For "blup", the model does not have contemporary groups (CG); instead uses the general mean.
With `cg=true`, the group effect will be considered.

For "tbv", the prediction is `TBV + randn()*sqrt([(1-rel)/rel]*var_g)`, where `TBV` is the true breeding value, `randn()` is a standard-normal deviation, `rel` is the reliability, and `var_g` is the genetic variance defined as the initial parameter.
This formula does not account for a change of genetic variance due to selection or drift.
If you provide an array for `rel`, this function will take it as individual reliability to generate GEBV.
Note that, once GEBV is calculated, it will not be updated unless you have `updategebv=true`.
"""
function genetic_evaluation!(pop::PTPopulation; method::String="blup", repeated::Bool=true, rel::Union{Vector{Union{Missing,Float64}},Float64}=0.5, verbose::Bool=false,
   keep_file::Bool=false, updategebv::Bool=false, cg::Bool=false)
   if method=="blup"
      write_files_for_blup(pop, "phenotype.txt", "pedigree.txt", repeated=repeated)
      genetic_evaluation_blup!(pop, "phenotype.txt", "pedigree.txt", repeated, verbose=verbose, cg=cg)
      if !keep_file
         rm("phenotype.txt")
         rm("pedigree.txt")
      end
   elseif method=="tbv"
      genetic_evaluation_tbv!(pop,rel,updategebv)
   elseif method=="tbv_ms"
      genetic_evaluation_tbv_ms!(pop,rel,updategebv)
   else
      error("unknown method: $(method)")
   end
end

@inbounds function genetic_evaluation_tbv!(pop::PTPopulation, rel::Float64, updategebv::Bool)
   par = pop.par
   var_g = get_var_poly(par) + get_var_qtl(par)
   if rel>0.0
      var_error = ((1-min(rel,1.0))/min(rel,1.0)) * var_g
   else
      var_error = 100*var_g
   end
   sd_error = sqrt(var_error)
   if updategebv
      if sd_error>0
         pop.df[:,:gebv] .= pop.df[:,:tbv] .+ (randn(pop.maxAnimal).*sd_error)
      else
         pop.df[:,:gebv] .= pop.df[:,:tbv]
      end
   else
      @inbounds for i=1:pop.maxAnimal
         if ismissing(pop.gebv[i])
            pop.df[i,:gebv] = pop.df[i,:tbv] + (randn() * sd_error)
         end
      end
   end
end

function genetic_evaluation_tbv!(pop::PTPopulation, rel::Vector{Union{Missing,Float64}}, updategebv::Bool)
   par = pop.par
   var_g = get_var_poly(par) + get_var_qtl(par)
   @inbounds for i=1:pop.maxAnimal
      r = ifelse(ismissing(rel[i]),0.0,1.0)
      if r>0.0
         var_error = ((1-min(r,1.0))/min(r,1.0)) * var_g
      else
         var_error = 100*var_g
      end
      sd_error = sqrt(var_error)
      if updategebv || ismissing(pop.df[i,:gebv])
         pop.df[i,:gebv] = pop.df[i,:tbv] + (randn() * sd_error)
      end
   end
end

function genetic_evaluation_tbv!(pop::PTPopulation, updategebv::Bool)
   rel = pop.df.rel
   genetic_evaluation_tbv!(pop, rel, updategebv)
end

function genetic_evaluation_tbv_ms!(pop::PTPopulation, rel::Vector{Union{Missing,Float64}}, updategebv::Bool)
   par = pop.par
   var_g = get_var_poly(par) + get_var_qtl(par)
   @inbounds for i=1:pop.maxAnimal
      r = ifelse(ismissing(rel[i]),0.0,1.0)
      if r>0.0
         var_error = ((1-min(r,1.0))/min(r,1.0)) * var_g
      else
         var_error = 100*var_g
      end
      ms_sd_error = sqrt(var_error/2)
      if pop.df[i,:sire]==0
         tbv_s = 0.0
      else
         tbv_s = pop.df[pop.df[i,:sire],:tbv]
      end
      if pop.df[i,:dam]==0
         tbv_d = 0.0
      else
         tbv_d = pop.df[pop.df[i,:dam],:tbv]
      end
      tbv_pa = (tbv_s + tbv_d)/2
      tbv_ms = pop.df[i,:tbv] - tbv_pa
      if updategebv || ismissing(pop.df[i,:gebv])
         pop.df[i,:gebv] = tbv_ms + (randn() * ms_sd_error)
      end
   end
end

function genetic_evaluation_tbv_ms!(pop::PTPopulation, updategebv::Bool)
   rel = pop.df.rel
   genetic_evaluation_tbv_ms!(pop, rel, updategebv)
end

function genetic_evaluation_tbv_ms!(pop::PTPopulation, rel::Float64, updategebv::Bool)
   relv = copy(pop.df.rel)
   relv .= rel
   genetic_evaluation_tbv_ms!(pop, relv, updategebv)
end

"""
    rel = approximated_reliability(pop::PTPopulation; de_extra::Float64=0.0, updaterel=true)

Calculate (roughly) approximated reliability based on the number of recorded progeny and own records.
You can supply an extra daugher equivalent (DE) `de_extra` to the reliability if needed.
If `updaterel=false`, this function computes the reliability only for new individuals that have not gotten the reliability yet.
With `update=true`, it re-calculates the values for all individuals.
"""
function approximated_reliability(pop::PTPopulation; de_extra::Float64=0.0, updaterel::Bool=true)
   par = pop.par
   h2 = par.h2_poly + par.h2_qtl
   rep = par.rep
   k = (4-h2)/h2
   ndau = pop.df[:,:nrecprog]
   nrec = pop.df[:,:nrec]
   rel = zeros(pop.maxAnimal)
   @inbounds for i=1:pop.maxAnimal
      # very rough approximation; see VanRaden and Wiggans (1991) for strict methods
      # PA
      s = pop.df[i,:sire]
      d = pop.df[i,:dam]
      if s>0; rel_s=rel[s]; else; rel_s=0.0; end
      if d>0; rel_d=rel[d]; else; rel_d=0.0; end
      rel_pa = (rel_s + rel_d)/4
      de_pa = k*rel_pa/(1-rel_pa)
      # own records
      if nrec[i]>0
         rel_y = nrec[i]*h2/(1 + (nrec[i]-1)*rep)
      else
         rel_y = 0.0
      end
      de_y = k*rel_y/(1-rel_y)
      # progeny contribution
      de_dau = ndau[i]
      # total DE
      de_animal = de_pa + de_y + de_dau + de_extra
      rel_animal = de_animal/(de_animal + k)
      if updaterel || ismissing(pop.df[i,:rel])
         rel[i] = rel_animal
      else
         rel[i] = pop.df[i,:rel]
      end
   end
   return rel
end

function update_approximated_reliability!(pop::PTPopulation; de_extra::Float64=0.0, updaterel::Bool=true)
   rel = approximated_reliability(pop,de_extra=de_extra,updaterel=updaterel)
   pop.df[:,:rel] .= rel
   return nothing
end

"""
    write_files_for_blup(pop::PTPopulation, datafile, pedfile; repeated=true)

Write the data and pedigree to files.
With `repeated=false`, this function writes only the first record per animal.

Each file has the following columns:
- Phenotype file
    1. animal
    2. sire
    3. dam
    4. inbreeding (from 0.0 to 1.0)
    5. birth year
    6. sex (1=male, 2=female)
    7. true breeding value (TBV = PBV + QBV)
    8. polygenic breeding value (PBV)
    9. QTL breeding value (QBV)
    10. sequential mumber of observations
    11. mu (fixed integer 1)
    12. contemporary group code
    13. phenotype
    14. permanent environmental (PE) effect
    15. residual term

- Pedigree file
    1. animal
    2. sire
    3. dam
    4. inbupg code
    5. inbreeding (from 0.0 to 1.0)
    6. birth year
    7. sex (1=male, 2=female)
    8. number of progeny with record(s)
    9. number of own record(s)
    10. true breeding value (TBV = PBV + QBV)
    11. polygenic breeding value (PBV)
    12. QTL breeding value (QBV)
"""
function write_files_for_blup(pop::PTPopulation, datafile::String, pedfile::String; repeated::Bool=true)
   io = open(datafile, "w");
   write(io, pop, header=false, missing="0", repeated=repeated, skipmissing=true)
   close(io)

   io = open(pedfile, "w");
   write_pedigree(io, pop, header=false)
   close(io)
end

function genetic_evaluation_blup!(pop::PTPopulation, datafile::String, pedfile::String, repeated::Bool=true; algorithm::String="pcg", verbose::Bool=false, cg::Bool=false)
   (lhs,rhs) = build_mme(pop, datafile, pedfile, repeated, cg)
   if algorithm=="pcg"
      sol = similar(rhs)
      sol .= 0.0
      dlhs = diag(lhs)
      @inbounds for i=1:length(dlhs)
         if dlhs[i] ≈ 0.0
            dlhs[i] = 1.0
         end
      end
      dlhs = Diagonal(dlhs)
      #cg!(sol, lhs, rhs, Pl=dlhs, verbose=true, log=true)
      cg!(sol, lhs, rhs, Pl=dlhs, verbose=verbose)
   else
      sol = lhs \ rhs
   end
   pop.df[:,:ebv] .= sol[1:pop.maxAnimal]
end

function build_mme(pop::PTPopulation, datafile::String, pedfile::String, repeated::Bool, cg::Bool; verbose::Bool=false)
   par = pop.par
   var_a = get_var_poly(par)
   var_g = get_var_qtl(par)
   var_u = var_a + var_g
   var_pe = get_var_pe(par)
   var_e = get_var_error(par)
   var_p = par.var_p

   lambda_u = var_e/var_u
   lambda_pe = var_e/var_pe
   if repeated
      if var_pe<=0.0 && repeated
         @warn "no use of repeatability model because of repeatability <= 0"
         repeatability_model = false
      else
         repeatability_model = true
      end
   else
      repeatability_model = false
   end
   if verbose; println("repeatability model: $(repeatability_model)"); end
   if verbose; println("lambda_u: $(lambda_u)"); end
   if verbose; println("lambda_pe: $(lambda_pe)"); end

   # details
   if repeatability_model
      # y = mu + u + pe + e
      pos_y = 13
      neff = 3
      random_type = ["u","pe","fixed"]
      offset = [0, pop.maxAnimal, pop.maxAnimal+pop.maxAnimal]
      if cg
         pos_eff = [1, 1, 12]
         nlev = [pop.maxAnimal, pop.maxAnimal, pop.maxCG]
      else
         pos_eff = [1, 1, 11]
         nlev = [pop.maxAnimal, pop.maxAnimal, 1]
      end
   else
      # y = mu + u + e
      pos_y = 13
      neff = 2
      random_type = ["u","fixed"]
      offset = [0, pop.maxAnimal]
      if cg
         pos_eff = [1, 12]
         nlev = [pop.maxAnimal, pop.maxCG]
      else
         pos_eff = [1, 11]
         nlev = [pop.maxAnimal, 1]
      end
   end
   if verbose; println("neff: $(neff)"); end
   if verbose; println("nlev: $(nlev)"); end
   if verbose; println("pos_eff: $(pos_eff)"); end
   if verbose; println("offset: $(offset)"); end
   if verbose; println("random_type: $(random_type)"); end

   return build_mme_lhs_rhs(datafile,pedfile,neff,pos_eff,offset,pos_y,nlev,random_type,lambda_u,lambda_pe, verbose=verbose)
end

function build_mme_lhs_rhs(datafile::String,pedfile::String,neff::Int,pos_eff::Vector{Int},offset::Vector{Int},pos_y::Int,nlev::Vector{Int},random_type::Vector{String},lambda_u::Float64,lambda_pe::Float64; verbose::Bool=false)
   neq = sum(nlev)
   lhs = SparseMatrixDict(neq,neq)
   rhs = zeros(neq)
   build_mme_lhs_rhs_fixed!(lhs,rhs,datafile,neff,pos_eff,offset,pos_y)
   for i in 1:neff
      if random_type[i]=="u"
         if verbose; println("ainverse"); end
         build_mme_lhs_rhs_ainverse!(lhs,pedfile,offset[i],lambda_u)
      elseif random_type[i]=="pe"
         if verbose; println("diagonal"); end
         build_mme_lhs_rhs_diagonal!(lhs,offset[i],nlev[i],lambda_pe)
      end
   end
   return sparse(lhs),rhs
end

function build_mme_lhs_rhs_fixed!(lhs::SparseMatrixDict,rhs::Vector{Float64},datafile::String,neff::Int,pos_eff::Vector{Int},offset::Vector{Int},pos_y::Int)
   # design matrix
   open(datafile) do io
      while !eof(io)
         line = readline(io)
         items = split(line)
         y = tryparse(Float64,items[pos_y])
         if isnothing(y) || y==0.0
            continue
         end
         level = tryparse.(Int,items[pos_eff])
         @inbounds for i in 1:neff
            addri = offset[i] + level[i]
            for j in 1:neff
               addrj = offset[j] + level[j]
               lhs[addri,addrj] = lhs[addri,addrj] + 1.0
            end
            rhs[addri] = rhs[addri] + y
         end
      end
   end
end

function build_mme_lhs_rhs_ainverse!(lhs::SparseMatrixDict,pedfile::String,offset::Int,lambda::Float64)
   # A-inverse
   w = (1.0, -0.5, -0.5)
   open(pedfile) do io
      while !eof(io)
         line = readline(io)
         items = split(line)
         ped = tryparse.(Int,items[1:3])
         code = tryparse(Float64,items[4])/1000.0
         @inbounds for i in 1:3
            addri = offset + ped[i]
            for j in 1:3
               addrj = offset + ped[j]
               if ped[i]>0 && ped[j]>0
                  val = w[i]*w[j]*code*lambda
                  lhs[addri,addrj] = lhs[addri,addrj] + val
               end
            end
         end
      end
   end
end

function build_mme_lhs_rhs_diagonal!(lhs::SparseMatrixDict,offset::Int,nlev::Int,lambda::Float64)
   # identity matrix
   @inbounds for i in 1:nlev
      addri = offset + i
      lhs[addri,addri] = lhs[addri,addri] + lambda
   end
end


text:
function initstate(interpreter::ConcreteInterpreter,
                   domain::Domain, problem::GenericProblem)
    return initstate(interpreter, domain, problem.objtypes, problem.init)
end

function initstate(interpreter::ConcreteInterpreter,
                   domain::Domain, objtypes::AbstractDict)
   types = Set{Term}([Compound(ty, Term[o]) for (o, ty) in objtypes])
   return GenericState(types, Set{Term}(), Dict{Symbol,Any}())
end

function initstate(interpreter::ConcreteInterpreter,
                   domain::Domain, objtypes::AbstractDict,
                   fluents::AbstractVector)
    state = initstate(interpreter, domain, objtypes)
    for t in fluents
        if t.name == :(==) # Non-Boolean fluents
            @assert length(t.args) == 2 "Assignments must have two arguments."
            term, val = t.args[1], t.args[2]
            @assert !isa(term, Var) "Initial terms cannot be unbound variables."
            state[term] = evaluate(domain, state, val)
        else # Boolean fluents
            push!(state.facts, t)
        end
    end
    return state
end

function initstate(interpreter::ConcreteInterpreter,
                   domain::Domain, objtypes::AbstractDict,
                   fluents::AbstractDict)
    state = initstate(interpreter, domain, objtypes)
    for (name, val) in fluents
        set_fluent!(state, val, name)
    end
    return state
end


text:
# Make an ITensorMap representing the transfer matrix T|v> = |Tv>
function transfer_matrix(ψ::InfiniteMPS)
  N = nsites(ψ)
  ψᴴ = prime(linkinds, dag(ψ))

  # The unit cell range
  # Turn into function `eachcellindex(ψ::InfiniteMPS, cell::Integer = 1)`
  cell₁ = 1:N
  # A transfer matrix made from the 1st unit cell of
  # the infinite MPS
  # TODO: make a `Cell` Integer type and call as `ψ[Cell(1)]`
  # TODO: make a TransferMatrix wrapper that automatically
  # primes and daggers, so it can be called with:
  # T = TransferMatrix(ψ[Cell(1)])
  ψ₁ = ψ[cell₁]
  ψ₁ᴴ = ψᴴ[cell₁]
  T₀₁ = ITensorMap(
    ψ₁,
    ψ₁ᴴ;
    input_inds=unioninds(commoninds(ψ[N], ψ[N + 1]), commoninds(ψᴴ[N], ψᴴ[N + 1])),
    output_inds=unioninds(commoninds(ψ[1], ψ[0]), commoninds(ψᴴ[1], ψᴴ[0])),
  )
  return T₀₁
end

# A TransferMatrix with the dominant space projected out
struct ProjectedTransferMatrix
  T::ITensorMap
  R::ITensor
end
function (T::ProjectedTransferMatrix)(v::ITensor)
  return v - T.T(v) + (v * T.R) * δ(inds(v))
end

function operator_matrix(H::InfiniteMPO)
  N = nsites(H)

  # MPO bond dimensions
  d⃗ₕ = OffsetArray([dim(commoninds(H[n], H[n + 1])) for n in 0:N], -1)

  # MPO link indices
  l⃗ₕ = OffsetArray([commoninds(H[n], H[n + 1]) for n in 0:N], -1)

  # Collect the on-site MPO operators obtained by projecting
  # the link indices onto certain elements
  H⃡ = Vector{Matrix{ITensor}}(undef, N)
  for n in 1:N
    H⃡[n] = Matrix{ITensor}(undef, d⃗ₕ[n - 1], d⃗ₕ[n])
    for i in 1:d⃗ₕ[n - 1], j in 1:d⃗ₕ[n]
      H⃡[n][i, j] = H[n] * setelt((l⃗ₕ[n - 1] .=> i)...) * setelt((l⃗ₕ[n] .=> j)...)
      # Assume lower triangular
      if j > i
        @assert norm(H⃡[n][i, j]) ≈ 0
      end
      # Also assume nonzero diagonal entries
      # are proportional to zero
    end

    println("norm.(H⃡[n])")
    display(norm.(H⃡[n]))
  end
  return H⃡
end

# Also input C bond matrices to help compute the right fixed points
# of ψ (R ≈ C * dag(C))
function left_environment(H::InfiniteMPO, ψ::InfiniteMPS, C::InfiniteMPS)
  N = nsites(H)
  @assert N == nsites(ψ)

  # Do the 1-site case first
  @assert N == 1

  # Transfer matrix between cells
  # T(v) performs a left multiplication
  T = transpose(transfer_matrix(ψ))
  R = C[N] * prime(dag(C[N]), linkinds(ψ, (N, N + 1)))
  R = translatecell(R, -1)

  # Operator valued matrix representation of MPO
  H⃡ = operator_matrix(H)

  # MPO bond dimensions
  d⃗ₕ = OffsetArray([dim(commoninds(H[n], H[n + 1])) for n in 0:N], -1)

  # MPO link indices
  l⃗ₕ = OffsetArray([commoninds(H[n], H[n + 1]) for n in 0:N], -1)

  # MPS link indices
  l⃗ = OffsetArray([commoninds(ψ[n], ψ[n + 1]) for n in 0:N], -1)

  # Starting point of identity
  n = 1 # Assuming we have only one site
  # This is the set of left quasi fixed points
  L⃗ = Vector{ITensor}(undef, d⃗ₕ[n])
  # The last one is defined to be identity
  L⃗[d⃗ₕ[n]] = δ(l⃗[n - 1]'..., dag(l⃗[n - 1])...)

  for a in reverse(1:(d⃗ₕ[n] - 1))
    YLᵃ = emptyITensor()
    for b in (a + 1):d⃗ₕ[n]
      YLᵃ += L⃗[b] * ψ[n] * H⃡[n][b, a] * dag(ψ[n]')
    end
    YLᵃ = translatecell(YLᵃ, -1)
    Hnᵃᵃ = H⃡[n][a, a]
    λᵃ = Hnᵃᵃ[1, 1]
    δᵃ = δ(inds(Hnᵃᵃ))

    if abs(λᵃ) > 1
      error(
        "Diagonal operator of the MPO must be proportional to identity with proprtionality constant less than or equal to 1.",
      )
    end

    if Hnᵃᵃ ≠ λᵃ * δᵃ
      error("Diagonal operator of the MPO must be proportional to identity.")
    end

    if λᵃ < 1
      error("Diagonal elements proportional to identity not currently supported.")
      # Solve using C21
    elseif isone(λᵃ) # equal to identity
      println("Diagonal elements equal to identity.")
      # Solve using C25a
      # Solve:
      # (Lᵃ|[𝟏 - Tᴸ + |R)(𝟏|] = (YLᵃ| - (YLᵃ|R)(𝟏|
      b = YLᵃ - (YLᵃ * R) * δ(inds(YLᵃ))
      L⃗[a], _ = linsolve(ProjectedTransferMatrix(T, R), b)

      # Get error
      err_lhs =
        L⃗[a] - translatecell(L⃗[a] * ψ[1] * dag(prime(ψ[1], "Link")), -1) +
        L⃗[a] * R * δ(inds(L⃗[a]))
      err_rhs = YLᵃ - YLᵃ * R * δ(inds(YLᵃ))
      @show norm(err_lhs - err_rhs)
    elseif iszero(λᵃ) # is zero
      println("Diagonal element is zero")
      L⃗[a] = YLᵃ
    else
      error(
        "The diagonal operator of the MPO must be zero, identity, or proportional to identity.",
      )
    end
  end
  # TODO: make this a vector of length N
  L = emptyITensor()
  n = 1
  for a in 1:d⃗ₕ[n]
    L += translatecell(L⃗[a], 1) * setelt(l⃗ₕ[n]... => a)
  end
  return InfiniteMPS([L])
end

function right_environment(H::InfiniteMPO, ψ::InfiniteMPS, C::InfiniteMPS)
  N = nsites(H)
  @assert N == nsites(ψ)

  # Do the 1-site case first
  @assert N == 1

  # Transfer matrix between cells
  # T(v) performs a right multiplication
  T = transfer_matrix(ψ)
  L = C[N] * prime(dag(C[N]), linkinds(ψ, (N, N + 1)))

  # Operator valued matrix representation of MPO
  H⃡ = operator_matrix(H)

  # MPO bond dimensions
  d⃗ₕ = OffsetArray([dim(commoninds(H[n], H[n + 1])) for n in 0:N], -1)

  # MPO link indices
  l⃗ₕ = OffsetArray([commoninds(H[n], H[n + 1]) for n in 0:N], -1)

  # MPS link indices
  l⃗ = OffsetArray([commoninds(ψ[n], ψ[n + 1]) for n in 0:N], -1)

  # Starting point of identity
  n = 1 # Assuming we have only one site
  # This is the set of left quasi fixed points
  R⃗ = Vector{ITensor}(undef, d⃗ₕ[n])
  # The first one is defined to be identity
  R⃗[1] = δ(l⃗[n]'..., dag(l⃗[n])...)

  for a in 2:d⃗ₕ[n]
    YRᵃ = emptyITensor()
    for b in 1:(a - 1)
      YRᵃ += R⃗[b] * ψ[n] * H⃡[n][a, b] * dag(ψ[n]')
    end
    YRᵃ = translatecell(YRᵃ, 1)
    Hnᵃᵃ = H⃡[n][a, a]
    λᵃ = Hnᵃᵃ[1, 1]
    δᵃ = δ(inds(Hnᵃᵃ))

    if abs(λᵃ) > 1
      error(
        "Diagonal operator of the MPO must be proportional to identity with proprtionality constant less than or equal to 1.",
      )
    end

    if Hnᵃᵃ ≠ λᵃ * δᵃ
      error("Diagonal operator of the MPO must be proportional to identity.")
    end

    if λᵃ < 1
      error("Diagonal elements proportional to identity not currently supported.")
      # Solve using C22
    elseif isone(λᵃ) # equal to identity
      println("Diagonal elements equal to identity.")
      # Solve using C25b
      # Solve:
      # [𝟏 - Tᴿ + |𝟏)(L|]|Rᵃ) = |YRᵃ) - |𝟏)(L|YRᵃ)
      b = YRᵃ - (L * YRᵃ) * δ(inds(YRᵃ))
      R⃗[a], _ = linsolve(ProjectedTransferMatrix(T, L), b)

      # Get error
      err_lhs =
        R⃗[a] - translatecell(ψ[1] * dag(prime(ψ[1], "Link")) * R⃗[a], 1) +
        L * R⃗[a] * δ(inds(R⃗[a]))
      err_rhs = YRᵃ - L * YRᵃ * δ(inds(YRᵃ))
      @show norm(err_lhs - err_rhs)
    elseif iszero(λᵃ) # is zero
      println("Diagonal element is zero")
      L⃗[a] = YLᵃ
    else
      error(
        "The diagonal operator of the MPO must be zero, identity, or proportional to identity.",
      )
    end
  end
  # TODO: make this a vector of length N
  R = emptyITensor()
  n = 1
  for a in 1:d⃗ₕ[n]
    R += R⃗[a] * setelt(l⃗ₕ[n]... => a)
  end
  return InfiniteMPS([R])
end

vumps(H::InfiniteMPO, ψ::InfiniteMPS; kwargs...) = vumps(H, orthogonalize(ψ, :); kwargs...)

# Find the best orthogonal approximation given
# the center tensors AC and C
function ortho(AC::ITensor, C::ITensor)
  E = AC * dag(C)
  U, P = polar(E, uniqueinds(AC, C))
  l = commoninds(U, P)
  return noprime(U, l)
end

#function right_ortho(C::ITensor, AC::ITensor)
#  E = dag(C) * AC
#  U, P = polar(E, uniqueinds(C, AC))
#  l = commoninds(U, P)
#  return noprime(U, l)
#end

function vumps(H::InfiniteMPO, ψ::InfiniteCanonicalMPS; nsweeps=10)
  for sweep in 1:nsweeps
    L = left_environment(H, ψ.AL, ψ.C)
    R = right_environment(H, ψ.AR, ψ.C)

    n = 1

    # 0-site effective Hamiltonian
    H⁰ = ITensorMap([L[n], R[n]])
    vals0, vecs0, info0 = eigsolve(H⁰, ψ.C[n])
    E0 = vals0[1]
    Cⁿ = vecs0[1]
    C = InfiniteMPS([Cⁿ])

    @show E0
    @show inds(Cⁿ)

    # 1-site effective Hamiltonian
    H¹ = ITensorMap([L[n - 1], H[n], R[n]])
    vals1, vecs1, info1 = eigsolve(H¹, ψ.AL[n] * ψ.C[n]; ishermition=true)
    E1 = vals1[1]
    ACⁿ = vecs1[1]
    AC = InfiniteMPS([ACⁿ])

    @show E1
    @show inds(ACⁿ)

    ALⁿ = ortho(ACⁿ, Cⁿ)
    ψL = InfiniteMPS([ALⁿ])

    @show norm(ALⁿ * Cⁿ - ACⁿ)

    @show inds(ALⁿ)
    #ARⁿ = ortho(ACⁿ, translatecell(Cⁿ, -1))
    ARⁿ = replacetags(ALⁿ, "Left" => "Right")

    @show norm(translatecell(Cⁿ, -1) * ARⁿ - ACⁿ)
    ψR = InfiniteMPS([ARⁿ])

    ψ = InfiniteCanonicalMPS(ψL, C, ψR)
  end

  return ψ
end


text:
# This file contains the types and definitions for the
# downsampled lippmann-Schwinger equation
using Interpolations

type FastDownSampled
    # type to encapsulate the fast application of M = I + omega^2G*spadiagm(nu)
    # Fast:: FastMslow
    Fast:: FastM
    n  :: Int64
    m  :: Int64
    downsample :: Int64
    # Sampling Matrix
    S :: SparseMatrixCSC{Float64,Int64}

    function FastDownSampled(Fast, downsample)
        # Building the Sampling matrix
        index1 = 1:downsample:Fast.n
        index2 = 1:downsample:Fast.m
        Sindex = spzeros(n,m)
        for ii = 0:round(Integer,(n-1)/downsample)
            for jj = 0:round(Integer,(m-1)/downsample)
                Sindex[1+ii*downsample,1+jj*downsample ] = 1
            end
        end

        index = find(Sindex[:])
        S = speye(n*m, n*m)
        Sampling = S[index,:];
        return new( Fast, Fast.n, Fast.m, downsample,Sampling)
    end
end

import Base.*

function *(M::FastDownSampled, b::Array{Complex128,1})
    knots = (collect(1:M.downsample:M.n), collect(1:M.downsample:M.m))
    nDown = round(Integer, (M.n-1)/M.downsample +1)
    mDown = round(Integer, (M.m-1)/M.downsample +1)
    B = reshape(b, nDown,mDown);

    itp_real = interpolate(real(B), BSpline(Quadratic(Reflect())), OnGrid())
    itp_imag = interpolate(imag(B), BSpline(Quadratic(Reflect())), OnGrid())

    interpU =     itp_real[collect(M.downsample:M.n+M.downsample-1)/M.downsample,collect(M.downsample:M.m+M.downsample-1)/M.downsample ] +
              1im*itp_imag[collect(M.downsample:M.n+M.downsample-1)/M.downsample,collect(M.downsample:M.m+M.downsample-1)/M.downsample];

    u = M.Fast*interpU[:]

    # Subsampling the wavefield
    return (M.S*u)[:]
end

max_stars_repo_name: Forgotten/Lippmann-Schingwer

max_issues_repo_name: Forgotten/Lippmann-Schingwer

max_forks_repo_name: Forgotten/Lippmann-Schingwer


text:
# Autogenerated wrapper script for libpsl_jll for x86_64-apple-darwin
export libpsl, psl

using Libiconv_jll
JLLWrappers.@generate_wrapper_header("libpsl")
JLLWrappers.@declare_library_product(libpsl, "@rpath/libpsl.5.dylib")
JLLWrappers.@declare_executable_product(psl)
function __init__()
    JLLWrappers.@generate_init_header(Libiconv_jll)
    JLLWrappers.@init_library_product(
        libpsl,
        "lib/libpsl.5.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        psl,
        "bin/psl",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module TestAqua

import Aqua
import FoldsCUDA

test_all() = Aqua.test_all(
    FoldsCUDA;
    # Default `Aqua.test_all(FoldsCUDA)` does not work due to ambiguities
    # in upstream packages:
    ambiguities = false,
    # support `[extras]`:
    project_extras = false,
)

test_ambiguities() = Aqua.test_ambiguities(FoldsCUDA)

end  # module


text: Dict{String,Any}("f" => Dict{String,Any}("value" => "false","type" => "bool"),"t" => Dict{String,Any}("value" => "true","type" => "bool"))

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
#
# The functions in this file manage error messages.
# I'm trying to give you useful hints if the iteration fails. These hints
# may become move detailed/verbose/bloviated over time.
#
"""
NewtonError(resfail, iline, resnorm, toosoon, itc, maxit, armmax, printerr)
Figure out what error message to print if the iteration fails.
"""
function NewtonError(resfail, iline, resnorm, toosoon, tol, itc, maxit, armmax, printerr)
    errcode = 0
    ~toosoon || (errcode = Lottery_Winner(resnorm, tol, printerr))
    itc < maxit || (errcode = MaxitError(resnorm, maxit, printerr))
    ~iline || (errcode = LineSearchFailure(maxit, itc, armmax, printerr))
    errcode != 0 || println("Something's wrong and I don't know what!")
    #
    #
    if printerr && ~toosoon
        println("Give the history array a look to see what's happening.")
        println("  ")
    end
    return errcode
end

"""
PTCError(resnorm, maxit, delta0, toosoon, tol, printerr)
"""
function PTCError(resnorm, maxit, delta0, toosoon, tol, printerr)
    ~toosoon || (errcode = Lottery_Winner(resnorm, tol, printerr))
    #if toosoon
    #errcode = Lottery_Winner(resnorm, tol, printerr)
    #else
    if printerr && ~toosoon
        println("PTC failure; increase maxit and/or delta0")
        println("Residual norm =", "  ", resnorm)
        println("Current values: maxit  =  ", maxit, ",  delta0 = ", delta0)
        println("Give the history array a look to see what's happening.")
        println("  ")
    end
    toosoon || (errcode = 10)
    return errcode
end

"""
LineSearchFailure
"""

function LineSearchFailure(maxit, itc, armmax, printerr)
    if printerr
        println("The line search failed at iteration", " ", itc)
        println("Termination with failure")
        println("Current values: maxit  =  ", maxit, ", armmax = ", armmax)
    end
    errcode = 1
    return errcode
end

"""
MaxitError
"""
function MaxitError(resnorm, maxit, printerr)
    if printerr
        println("Maximum iterations (maxit) of ", maxit, " exceeded")
        println("Convergence failure: residual norm too large  ", resnorm)
        println("Try increasing maxit and checking your function and
                 Jacobian for bugs.")
    end
    errcode = 10
    return errcode
end
"""
Lottery_Winner(resnorm, tol)
"""
function Lottery_Winner(resnorm, tol, printerr)
    if printerr
        println("Congratulations, your initial iterate met the termination criteria.")
        println("Residual norm = ", resnorm, " Tolerance = ", tol)
        println("  ")
    end
    errcode = -1
    return errcode
end

function Krylov_Error(lmaxit, ke_report)
    if ke_report == false
        println(
            "Linear solver did not meet termination criterion at least once.
    This does not mean the nonlinear solver will fail. lmaxit= ",
            lmaxit,
        )
    end
    return true
end