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

using Markdown
using InteractiveUtils

# ╔═╡ 1c630e40-628c-11ec-1793-05688dabf6c9
begin
	using Pkg; Pkg.activate(@__DIR__); Pkg.instantiate()
	Pkg.precompile()

	using PlutoUI
	using Statistics
	using StatsBase
	using RDatasets
	using Plots
	using StatsPlots
	using KernelDensity
	using Distributions
	using LinearAlgebra
	using HypothesisTests
	using PyCall
	using MLBase
end;

# ╔═╡ 7242a33b-ec37-48d6-afca-8c8387a45a83
PlutoUI.TableOfContents(aside=true, indent=true, depth=3)

# ╔═╡ dd2aa8b3-9da3-418d-bf31-74574483360d
html"""

"""

# ╔═╡ abc828f1-fd14-4329-8b47-f4135300f937
md"""
# 📊 Statistics

Having a solid understanding of statistics in data science allows us to understand our data better, and allows us to create a quantifiable evaluation of any future conclusions.

[YouTube Link](https://www.youtube.com/watch?v=AAGxWEJ_eWk)
"""

# ╔═╡ a8aae3d4-198c-482f-a972-05be31be39c1
md"""
## Load dataset

In this notebook, we will use eruption data on the faithful geyser (Yellowstone, USA). The data will contain wait times between every consecutive times the geyser goes off (in minutes) and the length of the eruptions (also in minutes).

Let's get the data first using the [RDatasets.jl](https://github.com/JuliaStats/RDatasets.jl) package...

**Note:** the dataset will be loaded as a `DataFrame`.
"""

# ╔═╡ 63420679-9b77-4412-ba1c-fc471da6b1e9
data = dataset("datasets", "faithful")

# ╔═╡ 9d56b2cc-c4c6-4ce2-8d9c-33174fe5bc33
md"""
We can use the `describe` function to get some stats from our data...
"""

# ╔═╡ 67ac4a1d-eb15-473f-a548-07425fd627d8
describe(data)

# ╔═╡ 5ee886cc-2730-40cf-8497-507dcdc08c5f
md"""
Let's call each column (feature) separately...
"""

# ╔═╡ 110fda28-1b77-44d6-8485-a4ec6e012011
begin
	eruption = data[!, :Eruptions]
	# eruption = data.Eruptions

	waittime = data[!, :Waiting]
	# waittime = data.Waiting
end;

# ╔═╡ 6e093fbb-1e3a-4131-b3c2-b5891e734fa7
md"""
Now, we can plot both variables in the same scatter plot...
"""

# ╔═╡ ed9da9c2-53e4-45d6-9ec1-02c8bf54c1a5
begin
	scatter(eruption, label="Eruptions")
	scatter!(waittime, label="Wait Time", xlabel="Index", ylabel="Time (minutes)")
end

# ╔═╡ 13137e30-cfde-459a-adb4-1ddc07e626e8
md"""
## Statistics plots

As you can see, this doesn't tell us much about the data... Let's try some statistical plots.
"""

# ╔═╡ ff9d63e7-808a-4465-a68d-659a222afbae
md"""
### Boxplot
"""

# ╔═╡ c3672771-248a-47f9-a019-2de03cdd79f4
boxplot(["Eruption Length"], eruption, legend=false,
		size=(200,400), whisker_width=1, ylabel="Time (minutes)")

# ╔═╡ 46b55947-428a-483f-a00f-a232e5329fbd
md"""
Statistical plots such as a boxplot (and a violin plot), can provide a much better understanding of the data. Here, we immediately see that the median time of each eruption is about 4 minutes.
"""

# ╔═╡ e0933f5c-50b5-47cf-bae8-7a45a219fefb
md"""
### Histogram

The next plot we will see is a histogram plot. Check its arguments below...

- `x`: AbstractVector of values to be binned.
- `bins`: Integer, NTuple{2,Integer}, AbstractVector or Symbol. Default is :auto (the Freedman-Diaconis rule). For histogram-types, defines the approximate number of bins to aim for, or the auto-binning algorithm to use (`:sturges`, `:sqrt`, `:rice`, `:scott` or `:fd`). For fine-grained control pass a Vector of break values, e.g. `range(minimum(x), stop = maximum(x), length = 25)`.
- `weights`: Vector of weights for the values in x, for weighted bin counts.
- `normalize`: Bool or Symbol. Histogram normalization mode. Possible values are: `false`/`:none` (no normalization, default), `true`/`:pdf` (normalize to a discrete Probability Density Function, where the total area of the bins is 1), `:probability` (bin heights sum to 1) and `:density` (the area of each bin, rather than the height, is equal to the counts - useful for uneven bin sizes).
- `bar_position`: Symbol. Choose from `:overlay` (default), `:stack`. (warning: May not be implemented fully).
- `bar_width`: nothing or Number. Width of bars in data coordinates. When nothing, chooses based on x (or y when orientation = `:h`).
- `bar_edges`: Bool. Align bars to edges (true), or centers (the default).
- `orientation`: Symbol. Horizontal or vertical orientation for bar types. Values `:h`, `:hor`, `:horizontal` correspond to horizontal (sideways, anchored to y-axis), and `:v`, `:vert`, and `:vertical` correspond to vertical (the default).

Let's plot eruption length data...
"""

# ╔═╡ 943c8862-b909-4750-96da-7bda62c5fc8c
histogram(eruption, bins=:sqrt, label="Eruption", color=:orange,
		  xlabel="Eruption Length (minutes)", ylabel="Frequency")

# ╔═╡ ed480ae6-3302-4d0e-b83e-53a9c395d5ce
md"""
### Kernel density estimates

Next, we will see how we can fit a kernel density estimation function to our data. We will make use of the [KernelDensity.jl](https://github.com/JuliaStats/KernelDensity.jl) package.

With Kernel density estimates, we can see how things are changing...
"""

# ╔═╡ ed58e41f-3f73-40c1-92ec-5f2845c62b97
ρ = kde(eruption)

# ╔═╡ 952506dc-6cd8-4ecd-ace8-b9d914560289
md"""
If we want the histogram and the kernel density graph to be aligned we need to remember that the "density contribution" of every point added to one of these histograms is `density * nb of elements * bin width`. We do this, since KDE are densities, while the y-axis is showing the frequency.

**Note:** Read more about kernel density estimates [here](https://en.wikipedia.org/wiki/Kernel_density_estimation).
"""

# ╔═╡ 3ea44987-7ba7-4df6-a82a-c11a98ea784d
# the range of density data (x-axis)
ρ.x

# ╔═╡ 31274ae3-20da-4f52-9153-6f6b4d4472cb
# density data (y-axis)
ρ.density

# ╔═╡ 53f3e80a-2c88-405f-9131-21d07b132389
# density * nb of elements * bin width
scaled_density = ρ.density .* length(eruption) .* 0.2

# ╔═╡ ccfe5310-cc85-457c-bfbc-9671062b9088
begin
	histogram(eruption, bins=:sqrt, label="Eruption", color=:lightblue,
		  	  xlabel="Eruption Length (minutes)", ylabel="Frequency")

	plot!(ρ.x, scaled_density, linewidth=1.5, color=:black, label="KDE fit")
end

# ╔═╡ cbf471e5-523c-4550-835b-5c45b9648d8c
md"""
## Probability distributions

First, we will take a look at one probablity distribution, namely the Normal Distribution and verify that it generates a bell curve.
"""

# ╔═╡ 914412a4-69d4-4d9f-b575-03cfa03a7873
# generating random normal data
normaldata = randn(100_000)

# ╔═╡ 86c9cf10-2ade-45de-8649-a8dd5b2ae586
begin
	# generating the KDE
	ρₙ = kde(normaldata)

	# scaling density estimates
	scaled_ρₙ = ρₙ.density .* length(normaldata) * 0.1
end

# ╔═╡ e137810d-9e3c-43c8-905c-4efc2213a130
begin
	# plotting histogram
	histogram(normaldata, color=:lightblue, label="Z ~ N(0,1)",
			  xlabel="Z", ylabel="Frequency")

	# plotting KDE
	plot!(ρₙ.x, scaled_ρₙ, linewidth=3, color=:black, label="KDE fit")
end

# ╔═╡ f6430d83-a92b-4c03-a07e-8a6146cc9e1c
md"""
Another way to generate the same plot is via using the [Distributions.jl](https://github.com/JuliaStats/Distributions.jl) package and choosing the probability distribution you want, and then drawing random numbers from it.
"""

# ╔═╡ 4ffd9a7a-c62c-40a4-8ec1-4e2351a2a66a
md"""
### Normal distribution

As an example, we will use `N = Normal()` below to draw random numbers from a [Normal distribution](https://en.wikipedia.org/wiki/Normal_distribution). The default $\mu$ and $\sigma$ values are 0 and 1 (standard Normal distribution).
"""

# ╔═╡ 978e2f4c-b411-459c-8024-4f6205484f63
begin
	# drawing numbers ~N(0,1)
	N = Normal()
	gaussdata = rand(N, 100_000)

	# calculating KDE
	rhoN = kde(gaussdata)
	scaled_rhoN = rhoN.density .* length(gaussdata) .* 0.1
end

# ╔═╡ d1483a5e-0b7e-4317-adb0-8f024fad21eb
begin
	# plotting histogram
	histogram(gaussdata, color=:lightblue, label="Z ~ N(0,1)",
			  xlabel="Z", ylabel="Frequency")

	# plotting KDE
	plot!(rhoN.x, scaled_rhoN, linewidth=3, color=:black, label="KDE fit")
end

# ╔═╡ 62470113-95dc-4b92-9e1d-3118efa0bd08
md"""
### Binomial distribution

Now, let's draw random numbers from a [Binomial distribution](https://en.wikipedia.org/wiki/Binomial_distribution).
"""

# ╔═╡ 7749b6b5-cf01-4812-b1ce-6fe3204b0bcf
begin
	# drawing numbers ~B with n=40 trials
	B = Binomial(40)
	binomdata = rand(B, 100_000)

	# calculating KDE
	rhoB = kde(binomdata)
	scaled_rhoB = rhoB.density .* length(binomdata) .* 0.5
end

# ╔═╡ ff68a550-9cc6-416e-96b8-07de164e5fc4
begin
	# plotting histogram
	histogram(binomdata, color=:lightblue, label="Z ~ B(40,0.5)",
			  xlabel="Z", ylabel="Frequency")

	# plotting KDE
	plot!(rhoB.x, scaled_rhoB, linewidth=3, color=:black, label="KDE fit")
end

# ╔═╡ 693de4f8-88b1-476f-9319-64edafb049af
md"""
### Fit data to a distribution

Next, we will try to fit a given set of numbers to a parametric distribution. Thus, we are going to use the StatsBase.jl `fit` function.
"""

# ╔═╡ 7f9949d9-22e2-4d4e-92ec-6e52e01d10aa
begin
	# creating 1000 random values
	randomdata = rand(1000)
	# fitting random values to a Standard Normal distribution
	fitted_distrib = fit(Normal, randomdata)
	# drawing 1000 values from the new fitted distribution
	new_data = rand(fitted_distrib, 1000)
end

# ╔═╡ 98301aae-61ca-46ab-936a-fbb310854778
# note that `fitted_distrib` is a `Distribution`
fitted_distrib

# ╔═╡ f270c785-3a7a-4e43-9f21-b9a6e199383b
begin
	# plotting new data (Normal distribution)
	histogram(new_data, nbins=20, fillalpha=0.3, label="Normal Data")

	# plotting original data
	histogram!(randomdata, nbins=20, linecolor=:red, fillalpha=0.3,
			   label="Original Data")
end

# ╔═╡ 2d5e5740-a3c6-4448-b79f-f231a0ac9d96
md"""
We can try to fit eruptions data to a Normal distribution, although eruptions distribution does not seem to be Normal.
"""

# ╔═╡ ef651442-e846-4a5e-9f3b-0da63afeca16
begin
	erup_distrib = fit(Normal, eruption)
	normal_erup = rand(erup_distrib, 1000)
end

# ╔═╡ a5a54536-4b60-4a9c-b423-2804cce37382
erup_distrib

# ╔═╡ 25b15bda-9bb4-4ed9-90e4-c1f937051314
begin
	# plotting fitted data (Normal distribution)
	histogram(normal_erup, nbins=20, fillalpha=0.3, label="Normal Data")

	# plotting original data
	histogram!(eruption, nbins=20, linecolor=:red, fillalpha=0.3,
			   label="Original Data")
end

# ╔═╡ b82db811-d97d-4efd-81bc-a48780b23dcb
md"""
## Hypothesis testing

Next, we will perform hypothesis testing using the [HypothesisTests.jl](https://github.com/JuliaStats/HypothesisTests.jl) package.

Let's try to use the `OneSampleTTest` function to perform a one sample t-test of random sampled values.

The one-sample t-test is a statistical hypothesis test used to determine whether an unknown population mean is different from a specific value. In this case, the default value under $H_0$ is 0.
"""

# ╔═╡ af54608d-209f-4d16-9573-b4b8634d2acf
# generating 1000 random normal values
hypdata = randn(1000)

# ╔═╡ 867dab59-1c8d-44d9-9add-9f3999f1b86c
# performing one sample t-test
OneSampleTTest(hypdata)

# ╔═╡ 47c7f6df-0696-4b1e-b6c2-7eb8e782eb82
md"""
**Note:** as we can see above, we cannot reject the null hypothesis $H_0$. This was expected, since we drew sample values from a Standard Normal distribution, with $\mu=0$.

Now we can perform a one sample t-test of eruption data. We will probably be able to reject $H_0$, since the eruption length mean is close to 3.5 minutes...
"""

# ╔═╡ 6555c2a2-227f-419b-9b8b-d563767f3ea5
OneSampleTTest(eruption)

# ╔═╡ 78585d9d-5696-41dd-80f0-abab1585b008
md"""
## Correlation statistics

First of all, let's see visually if there is a correlation between `eruption` and `waittime` variables using a scatter plot...
"""

# ╔═╡ 650a1af5-b413-4794-bf51-77fa34cfb820
scatter(eruption, waittime, legend=false, color=:orange,
		xlabel="Eruption Length (minutes)", markersize=3,
		ylabel = "Time between Eruptions (minutes)")

# ╔═╡ 80db913d-faa8-43d5-9794-e17b02823cc5
md"""
These two variables may have a strong positive correlation. We can calculate both `Pearson` and `Spearman` correlation coefficients to check that.
"""

# ╔═╡ fcb84ac7-3bca-4284-8939-a55bd78a02fa
begin
	spearman_r = corspearman(eruption, waittime)
	"r(Spearman) = $spearman_r"
end

# ╔═╡ 960b8739-10d1-490d-ac49-e51527386161
begin
	pearson_r = cor(eruption, waittime)
	"r(Pearson) = $pearson_r"
end

# ╔═╡ 5bd50a69-4edb-435a-9263-e0b5d47851aa
md"""
Interesting! This means that the next time you visit Yellowstone National part ot see the faithful geysser and you have to wait for too long for it to go off, you will likely get a longer eruption!

Currently we are using the p-value of Spearman and Pearson correlation coefficients from Python. But you can follow the formula [here](https://stackoverflow.com/questions/53345724/how-to-use-julia-to-compute-the-pearson-correlation-coefficient-with-p-value) to implement your own.

Let's import the `scipy.stats` Python module.

**Note:** sometimes there are some issues getting Python and Julia to communicate as desired. We can explicity add `scipy.stats` using the [Conda.jl](https://github.com/JuliaPy/Conda.jl) package:

```julia
julia> using Conda
julia> Conda.add("scipy")
```
"""

# ╔═╡ d1ef6cd2-fda1-47f6-af12-8492ea70a539
scipy_stats = pyimport("scipy.stats");

# ╔═╡ c5f77715-02cb-437a-a45f-3f3a4dfba67b
begin
	rₛ, pvalueₛ = scipy_stats.spearmanr(eruption, waittime)
	"r(Spearman) = $(rₛ) || p-value = $(pvalueₛ)"
end

# ╔═╡ 615e7c40-bafe-475f-959d-29c5f5e845e9
begin
	rₚ, pvalueₚ = scipy_stats.pearsonr(eruption, waittime)
	"r(Pearson) = $(rₚ) || p-value = $(pvalueₚ)"
end

# ╔═╡ 32bc5a24-638e-43bc-a5c1-bd0a85db7849
md"""
## Confusion matrix & AUC

Finally, we will cover basic tools you will need such as confusion matrix and AUC scores. We use the [MLBase.jl](https://github.com/JuliaStats/MLBase.jl) package for that.

First, let's calculate a confusion matrix, given a target vector `y₁` and a prediction vector `ŷ₁`...
"""

# ╔═╡ 838a8a32-95f6-48c0-b87d-28a66ef431c5
begin
	# generating target and prediction vectors
	y₁ = [1, 1, 1, 1, 1, 1, 1, 2]
	yhat₁ = [1, 1, 2, 2, 1, 1, 1, 1]
end;

# ╔═╡ 68d7cd6c-bd76-406b-a3d7-4528ee9671fd
# compute default confusion matrix
C = confusmat(2, y₁, yhat₁)

# ╔═╡ 2bcf8647-368d-4b74-8e1a-4a4b1a577fdc
md"""
We may normalize cell values per class...
"""

# ╔═╡ 86ad3850-f631-42b2-8d19-bf20cba89c6b
C ./ sum(C, dims=2)

# ╔═╡ 38a2bb94-1d8f-4e1a-a07f-3bc1332301fd
md"""
Also, we may compute the accuracy (correct rate) from confusion matrix...

```math
accuracy = \frac{TP}{TP+TN+FP+FN}
```
"""

# ╔═╡ ea2c50c6-9404-4ba4-8fe7-4d1e1b3246ff
acc = sum(diag(C)) / length(y₁)

# ╔═╡ 70cbc1f9-171e-4664-87d1-6e73f0bc841c
md"""
Instead of calculating accuracy manually, we could use the `correctrate` function...
"""

# ╔═╡ 188c0219-8a54-4f28-95fd-9f2f525c12be
acc == correctrate(y₁, yhat₁)

# ╔═╡ 063d30ad-0ac4-4d50-b66b-9b06bf155fd8
md"""
Now, we can use the `roc` function to calculate ROC scores...
"""

# ╔═╡ f3cca2ef-17ec-44e0-996d-8a0cb2f2fcb1
begin
	# generating target and prediction vectors
	y₂ = [1, 1, 1, 1, 1, 1, 1, 0]
	yhat₂ = [1, 1, 0, 0, 1, 1, 1, 1]
end;

# ╔═╡ 9af05d0f-5daf-48f4-9949-79db70ae4bbc
ROC = MLBase.roc(y₂, yhat₂)

# ╔═╡ 92f98fe2-9154-4029-8282-45dee3f11d31
md"""
We may also use both `recall` and `precision` functions to compute these metrics...
"""

# ╔═╡ 0cfa5643-b65b-4374-a227-8f036ce32bbf
begin
	rec = recall(ROC)
	"Recall = $rec"
end

# ╔═╡ 93ea2a29-bac3-4629-bac5-5d03e8e45d80
begin
	prec = precision(ROC)
	"Precision = $prec"
end

# ╔═╡ Cell order:
# ╟─1c630e40-628c-11ec-1793-05688dabf6c9
# ╟─7242a33b-ec37-48d6-afca-8c8387a45a83
# ╟─dd2aa8b3-9da3-418d-bf31-74574483360d
# ╟─abc828f1-fd14-4329-8b47-f4135300f937
# ╟─a8aae3d4-198c-482f-a972-05be31be39c1
# ╠═63420679-9b77-4412-ba1c-fc471da6b1e9
# ╟─9d56b2cc-c4c6-4ce2-8d9c-33174fe5bc33
# ╠═67ac4a1d-eb15-473f-a548-07425fd627d8
# ╟─5ee886cc-2730-40cf-8497-507dcdc08c5f
# ╠═110fda28-1b77-44d6-8485-a4ec6e012011
# ╟─6e093fbb-1e3a-4131-b3c2-b5891e734fa7
# ╠═ed9da9c2-53e4-45d6-9ec1-02c8bf54c1a5
# ╟─13137e30-cfde-459a-adb4-1ddc07e626e8
# ╟─ff9d63e7-808a-4465-a68d-659a222afbae
# ╠═c3672771-248a-47f9-a019-2de03cdd79f4
# ╟─46b55947-428a-483f-a00f-a232e5329fbd
# ╟─e0933f5c-50b5-47cf-bae8-7a45a219fefb
# ╠═943c8862-b909-4750-96da-7bda62c5fc8c
# ╟─ed480ae6-3302-4d0e-b83e-53a9c395d5ce
# ╠═ed58e41f-3f73-40c1-92ec-5f2845c62b97
# ╟─952506dc-6cd8-4ecd-ace8-b9d914560289
# ╠═3ea44987-7ba7-4df6-a82a-c11a98ea784d
# ╠═31274ae3-20da-4f52-9153-6f6b4d4472cb
# ╠═53f3e80a-2c88-405f-9131-21d07b132389
# ╠═ccfe5310-cc85-457c-bfbc-9671062b9088
# ╟─cbf471e5-523c-4550-835b-5c45b9648d8c
# ╠═914412a4-69d4-4d9f-b575-03cfa03a7873
# ╠═86c9cf10-2ade-45de-8649-a8dd5b2ae586
# ╠═e137810d-9e3c-43c8-905c-4efc2213a130
# ╟─f6430d83-a92b-4c03-a07e-8a6146cc9e1c
# ╟─4ffd9a7a-c62c-40a4-8ec1-4e2351a2a66a
# ╠═978e2f4c-b411-459c-8024-4f6205484f63
# ╠═d1483a5e-0b7e-4317-adb0-8f024fad21eb
# ╟─62470113-95dc-4b92-9e1d-3118efa0bd08
# ╠═7749b6b5-cf01-4812-b1ce-6fe3204b0bcf
# ╠═ff68a550-9cc6-416e-96b8-07de164e5fc4
# ╟─693de4f8-88b1-476f-9319-64edafb049af
# ╠═7f9949d9-22e2-4d4e-92ec-6e52e01d10aa
# ╠═98301aae-61ca-46ab-936a-fbb310854778
# ╠═f270c785-3a7a-4e43-9f21-b9a6e199383b
# ╟─2d5e5740-a3c6-4448-b79f-f231a0ac9d96
# ╠═ef651442-e846-4a5e-9f3b-0da63afeca16
# ╠═a5a54536-4b60-4a9c-b423-2804cce37382
# ╠═25b15bda-9bb4-4ed9-90e4-c1f937051314
# ╟─b82db811-d97d-4efd-81bc-a48780b23dcb
# ╠═af54608d-209f-4d16-9573-b4b8634d2acf
# ╠═867dab59-1c8d-44d9-9add-9f3999f1b86c
# ╟─47c7f6df-0696-4b1e-b6c2-7eb8e782eb82
# ╠═6555c2a2-227f-419b-9b8b-d563767f3ea5
# ╟─78585d9d-5696-41dd-80f0-abab1585b008
# ╠═650a1af5-b413-4794-bf51-77fa34cfb820
# ╟─80db913d-faa8-43d5-9794-e17b02823cc5
# ╠═fcb84ac7-3bca-4284-8939-a55bd78a02fa
# ╠═960b8739-10d1-490d-ac49-e51527386161
# ╟─5bd50a69-4edb-435a-9263-e0b5d47851aa
# ╠═d1ef6cd2-fda1-47f6-af12-8492ea70a539
# ╠═c5f77715-02cb-437a-a45f-3f3a4dfba67b
# ╠═615e7c40-bafe-475f-959d-29c5f5e845e9
# ╟─32bc5a24-638e-43bc-a5c1-bd0a85db7849
# ╠═838a8a32-95f6-48c0-b87d-28a66ef431c5
# ╠═68d7cd6c-bd76-406b-a3d7-4528ee9671fd
# ╟─2bcf8647-368d-4b74-8e1a-4a4b1a577fdc
# ╠═86ad3850-f631-42b2-8d19-bf20cba89c6b
# ╟─38a2bb94-1d8f-4e1a-a07f-3bc1332301fd
# ╠═ea2c50c6-9404-4ba4-8fe7-4d1e1b3246ff
# ╟─70cbc1f9-171e-4664-87d1-6e73f0bc841c
# ╠═188c0219-8a54-4f28-95fd-9f2f525c12be
# ╟─063d30ad-0ac4-4d50-b66b-9b06bf155fd8
# ╠═f3cca2ef-17ec-44e0-996d-8a0cb2f2fcb1
# ╠═9af05d0f-5daf-48f4-9949-79db70ae4bbc
# ╟─92f98fe2-9154-4029-8282-45dee3f11d31
# ╠═0cfa5643-b65b-4374-a227-8f036ce32bbf
# ╠═93ea2a29-bac3-4629-bac5-5d03e8e45d80

max_stars_repo_name: fnaghetini/DataScience

max_issues_repo_name: fnaghetini/DataScience

max_forks_repo_name: fnaghetini/DataScience


text:
export Sphere1

"""
$(TYPEDEF)

Sphere1 is a S1 mechanization of one Circular rotation, with `theta in [-pi,pi)`.
"""
struct Sphere1 <: InferenceVariable
  dims::Int
  labels::Vector{String}
  manifolds::Tuple{Symbol}
  Sphere1(;labels::Vector{<:AbstractString}=String[]) = new(1, labels, (:Circular,))
end

ext: jl


text:
export FacElem, FacElemMon

export transform

################################################################################
#
#  Multiplicative representation
#
################################################################################

# abstract nonsense

const FacElemMonDict = ObjectIdDict()

function call{T}(x::FacElemMon{T})
  z = FacElem{T}()
  z.parent = x
  return z
end

function FacElem{B}(base::Array{B, 1}, exp::Array{fmpz, 1})

  length(base) == 0 && error("Array must not be empty")

  length(base) != length(exp) && error("not the same length")

  z = FacElem{B}()

  for i in 1:length(base)
    if exp[i] == 0
      continue
    end
    z.fac[base[i]] = exp[i]
  end

  z.parent = FacElemMon{typeof(base[1])}(parent(base[1]))
  return z
end

parent(x::FacElem) = x.parent

base_ring(x::FacElemMon) = x.base_ring

base_ring(x::FacElem) = base_ring(parent(x))

base(x::FacElem) = keys(x.fac)

function deepcopy{B}(x::FacElem{B})
  z = FacElem{B}()
  z.fac = _deepcopy(x.fac)
  if isdefined(x, :parent)
    z.parent = x.parent
  end
  return z
end

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, x::FacElemMon)
  print(io, "Factored elements over $(x.base_ring)")
end

function show(io::IO, x::FacElem)
  print(io, "Factored element with data\n$(x.fac)")
end

################################################################################
#
#  Inverse
#
################################################################################

function inv(x::FacElem{nf_elem})
  y = deepcopy(x)
  for a in base(y)
    y.fac[a] = -y.fac[a]
  end
  return y
end

################################################################################
#
#  Exponentiation
#
################################################################################

function pow!{T <: Union{fmpz, Integer}}(z::FacElem, x::FacElem, y::T)
  z.fac = _deepcopy(x.fac)
  for a in base(x)
    # this should be inplac
    z.fac[a] = y*x.fac[a]
  end
end

function ^(x::FacElem, y::fmpz)
  z = parent(x)()
  if y == 0
    return z
  else
    z.fac = _deepcopy(x.fac)
    for a in base(x)
      # this should be inplac
      z.fac[a] = y*x.fac[a]
    end
    return z
  end
end

^(x::FacElem, y::Integer) = ^(x, fmpz(y))

################################################################################
#
#  Multiplication
#
################################################################################

function mul!{B}(z::FacElem{B}, x::FacElem{B}, y::FacElem{B})
  z.fac = _deepcopy(x.fac)
  for a in base(y)
    if haskey(x.fac, a)
      z.fac[a] = z.fac[a] + y.fac[a]
    else
      z.fac[a] = y.fac[a]
    end
  end
end

function *{B}(x::FacElem{B}, y::FacElem{B})
  z = deepcopy(x)
  for a in base(y)
    if haskey(x.fac, a)
      z.fac[a] = z.fac[a] + y.fac[a]
    else
      z.fac[a] = y.fac[a]
    end
  end
  return z
end

function *{B}(x::FacElem{B}, y::B)
  z = deepcopy(x)
  if haskey(x.fac, y)
    z.fac[y] = z.fac[y] + 1
  else
    z.fac[y] = 1
  end
  return z
end

function div{B}(x::FacElem{B}, y::FacElem{B})
  z = deepcopy(x)
  for a in base(y)
    if haskey(x.fac, a)
      z.fac[a] = z.fac[a] - y.fac[a]
    else
      z.fac[a] = -y.fac[a]
    end
  end
  return z
end

################################################################################
#
#  Apply transformation
#
################################################################################

function _transform{T}(x::Array{FacElem{T}, 1}, y::fmpz_mat)
  length(x) != rows(y) &&
              error("Length of array must be number of rows of matrix")

  z = Array(FacElem{T}, cols(y))

  t = parent(x[1])()

  for i in 1:cols(y)
    z[i] = x[1]^y[1,i]
    for j in 2:rows(y)
      if y[j, i] == 0
        continue
      end
      pow!(t, x[j], y[j, i])
      mul!(z[i], z[i], t)
    end
  end
  return z
end

function transform{T}(x::Array{FacElem{T}, 1}, y::fmpz_mat)
  return _transform(x, y)
end

################################################################################
#
#  Evaluate
#
################################################################################

doc"""
***
  evaluate{T}(x::FacElem{T}) -> T

> Expands or evaluates the factored element, i.e. actually computes the
> value.
"""
function evaluate{T}(x::FacElem{T})
  function ev(d)#d::Dict{T, fmpz})
    z = one(base_ring(x))
    if length(d)==0
      return z
    elseif length(d)==1
      x = first(d)
      return x[1]^x[2]
    end
    b = similar(d)
    for (k,v) in d
      if v>-10 && v<10
        z *= k^v
      else
        r = isodd(v) ? 1 :0
        vv = div(v-r, 2)
        if vv!=0
          b[k] = vv
        end
        if r!=0
          z*= k
        end
      end
    end
    return ev(b)^2*z
  end

  return ev(x.fac)
end

doc"""
***
  naive_evaluate{T}(x::FacElem{T}) -> T

> Expands or evaluates the factored element, i.e. actually computes the
> value. Uses the obvious naive algorithm. Faster for input in finite rings.
"""
function naive_evaluate{T}(x::FacElem{T})
  z = one(base_ring(x))
  d = x.fac
  if length(d)==0
    return z
  end
  for (k,v) in d
    z *= k^v
  end
  return z
end

################################################################################
#
#  Auxillary deep copy functions
#
################################################################################

function _deepcopy{K, V}(x::Dict{K, V})
  z = typeof(x)()
  for a in keys(x)
    z[a] = deepcopy(x[a])
  end
  return z
end

################################################################################
#
#  Parent object overloading
#
################################################################################

function call{T}(F::FacElemMon{T}, a::T)
  z = F()
  z.fac[a] = fmpz(1)
  return z
end


text:
n = 2
m = 1

A = [1.0 1.0; 0.0 1.0]
B = [0.0; 1.0]

Q = [1.0 0.0; 0.0 1.0]
R = [1.0]

x_init = [1.0; 0.0]

nz = 2 * n + m + 2 * n

function unpack(z)
	x2 = z[1:n]
	u1 = z[n .+ (1:m)]
	x1 = z[n + m .+ (1:n)]
	y2 = z[n + m + n .+ (1:n)]
	y1 = z[n + m + n + n .+ (1:n)]

	return x2, u1, x1, y2, y1
end

z0 = rand(nz)
unpack(z0)
ones(1)' * R * ones(1)
function lagrangian(z, x_init)
	x2, u1, x1, y2, y1 = unpack(z)

	L = 0.0

	# objective
	L += transpose(x2) * Q * x2
	L += (transpose(u1) * R * u1)[1]
	L += transpose(x1) * Q * x1

	# constraints
	L += transpose(y2) * (A * x1 + B * u1[1] - x2)
	L += transpose(y1) * (x_init - x1)

	return L
end

@variables z_sym[1:nz]
@variables x_sym[1:n]

L = lagrangian(z_sym, x_sym)
L = simplify.(L)

dL = Symbolics.gradient(L, z_sym)
ddL = Symbolics.hessian(L, z_sym)
ddL_x = Symbolics.jacobian(dL, x_sym)

L_grad = eval(Symbolics.build_function(dL, z_sym, x_sym)[1])
L_hess = eval(Symbolics.build_function(ddL, z_sym, x_sym)[1])
Lx_hess = eval(Symbolics.build_function(ddL_x, z_sym, x_sym)[1])

function solve(z0, x_init)
	z = copy(z0)

	∇L = L_grad(z, x_init)

	println()
	for i = 1:10
		r_norm = norm(∇L)
		println("iter ($i)")
		println("	grad norm: $(r_norm)")
		r_norm < 1.0e-8 && break

		∇²L = L_hess(z, x_init)

		Δ = ∇²L \ ∇L

		α = 1.0
		ẑ = z - α * Δ

		∇L̂ = L_grad(ẑ, x_init)

		iter = 1
		while norm(∇L̂) > norm(∇L)
			α *= 0.5
			ẑ = z - α * Δ
			∇L̂ = L_grad(ẑ, x_init)
			iter += 1

			if iter > 100
				break
			end
		end

		z = ẑ
		∇L = ∇L̂
	end

	return z, ∇L, L_hess(z, x_init)
end

z_sol, ∇L_sol, ∇²L_sol = solve(z0, x_init)
K_direct = (∇²L_sol \ Lx_hess(z_sol, x_init))[n .+ (1:m), :]

# LQR solution
K, P = tvlqr([A], [B], [Q, Q], [R[1,1]])

@assert norm(K_direct - K[1]) < 1.0e-8

max_issues_repo_name: thowell/DirectMotionPlanning


text:
"""
    CurrentSensor(; name)

Sensor to measure the current through a branch in amperes.
"""
function CurrentSensor(; name)
    @named p = Pin()
    @named n = Pin()
    @variables i(t)=1.0
    eqs = [
        p.v ~ n.v
        i ~ p.i
        i ~ -n.i
    ]
    ODESystem(eqs, t, [i], [], systems=[p, n]; name=name)
end

"""
    PotentialSensor(; name)

Sensor to measure the potential in volt.
"""
function PotentialSensor(; name)
    @named p = Pin()
    @variables phi(t)=1.0
    eqs = [
        p.i ~ 0
        phi ~ p.v
    ]
    ODESystem(eqs, t, [phi], [], systems=[p]; name=name)
end

"""
    VoltageSensor(; name)

Sensor to measure the potential difference between two pins in volt.
"""
function VoltageSensor(; name)
    @named p = Pin()
    @named n = Pin()
    @variables v(t)=1.0
    eqs = [
        p.i ~ 0
        n.i ~ 0
        v ~ p.v - n.v
    ]
    ODESystem(eqs, t, [v], []; systems=[p, n], name=name)
end

"""
    PowerSensor(; name)

Sensor to measure the power in watt.
"""
function PowerSensor(; name)
    @named pc = Pin()
    @named nc = Pin()
    @named pv = Pin()
    @named nv = Pin()
    @named voltage_sensor = VoltageSensor()
    @named current_sensor = CurrentSensor()
    @variables power(t)=1.0
    eqs = [
        connect(voltage_sensor.p, pv)
        connect(voltage_sensor.n, nv)
        connect(current_sensor.p, pc)
        connect(current_sensor.n, nc)
        power ~ current_sensor.i * voltage_sensor.v
    ]
    ODESystem(eqs, t, [power], []; systems=[pc, nc, pv, nv, voltage_sensor, current_sensor], name=name)
end

"""
    MultiSensor(; name)

Sensor to measure current [A], voltage [V] and power [W].
"""
function MultiSensor(; name)
    @named pc = Pin()
    @named nc = Pin()
    @named pv = Pin()
    @named nv = Pin()
    @named voltage_sensor = VoltageSensor()
    @named current_sensor = CurrentSensor()
    sts = @variables begin
        i(t)=1.0
        v(t)=1.0
    end
    eqs = [
        connect(voltage_sensor.p, pv)
        connect(voltage_sensor.n, nv)
        connect(current_sensor.p, pc)
        connect(current_sensor.n, nc)
        i ~ current_sensor.i
        v ~ voltage_sensor.v
    ]
    ODESystem(eqs, t, sts, []; systems=[pc, nc, pv, nv, voltage_sensor, current_sensor], name=name)
end


text:
"""
    KMarkovUpdater

Updater that stores the k most recent observations as the belief.

Example:

```julia
up = KMarkovUpdater(5)
s0 = initialstate(pomdp, rng)
initial_observation = gen(DDNNode(:o), pomdp, s0, rng)
initial_obs_vec = fill(initial_observation, 5)
hr = HistoryRecorder(rng=rng, max_steps=100)
hist = simulate(hr, pomdp, policy, up, initial_obs_vec, s0)
```
"""
struct KMarkovUpdater <: Updater
    k::Int
end

# this is a new type because `history(::AbstractVector)` and `currentobs(::AbstractVector)` are too broad
struct PreviousObservations{K,O} <: AbstractVector{O}
    storage::NTuple{K,O}
end

PreviousObservations(v::AbstractVector) = PreviousObservations(tuple(v...))

Base.size(::PreviousObservations{K}) where {K} = (K,)
Base.getindex(p::PreviousObservations, i::Int) = p.storage[i]
POMDPs.history(p::PreviousObservations) = collect((o=obs,) for obs in p)
POMDPs.currentobs(p::PreviousObservations) = p[end]

function initialize_belief(bu::KMarkovUpdater, obs_vec::AbstractVector)
    if length(obs_vec) != bu.k
        error("KMarkovUpdater: The length of the initial observation vector
               does not match the number of observation to stack\n"*throw_example(bu))
    end
    return PreviousObservations(obs_vec)
end

function update(bu::KMarkovUpdater, old_b::AbstractVector{O}, action, obs) where {O}
    if !isa(obs, O)
        error("KMarkovUpdater: Observation did not match previous observation type.\n"*throw_example(bu))
    end
    return PreviousObservations((old_b[end-(bu.k-2):end]..., obs))
end

function initialize_belief(bu::KMarkovUpdater, obs_vec)
    error("KMarkovUpdater: To initialize the belief, pass in a vector of observation.\n"*throw_example(bu))
end

function throw_example(bu::KMarkovUpdater)
    example = """
    Did you forget to pass the initial observation to the simulator?
    Example:
    ```julia
    up = KMarkovUpdater(5)
    s0 = initialstate(pomdp, rng)
    initial_observation = gen(DDNNode(:o), pomdp, s0, rng)
    initial_obs_vec = fill(initial_observation, 5)
    hr = HistoryRecorder(rng=rng, max_steps=100)
    hist = simulate(hr, pomdp, policy, up, initial_obs_vec, s0)
    ```
    """
    return example
end


text:
using Distributions
using Base.Test

# Core testing procedure

function test_mixture(g::UnivariateMixture, n::Int, ns::Int)
    X = zeros(n)
    for i = 1:10
        X[i] = rand(g)
    end

    cs = components(g)
    pr = probs(g)
    @assert length(cs) == length(pr)

    # mean
    mu = 0.0
    for k = 1:length(cs)
        mu += pr[k] * mean(cs[k])
    end
    @test_approx_eq mean(g) mu

    # ground-truth
    p0 = zeros(n)
    for i = 1:n
        p0_i = 0.0
        x_i = X[i]
        for k = 1:length(cs)
            p0_i += pr[k] * pdf(cs[k], x_i)
        end
        p0[i] = p0_i
    end
    lp0 = log(p0)

    for i = 1:n
        @test_approx_eq pdf(g, X[i]) p0[i]
        @test_approx_eq logpdf(g, X[i]) lp0[i]
    end

    p_e = pdf(g, X)
    lp_e = logpdf(g, X)
    @assert isa(p_e, Vector{Float64}) && length(p_e) == n
    @assert isa(lp_e, Vector{Float64}) && length(lp_e) == n

    @test_approx_eq p_e p0
    @test_approx_eq lp_e lp0

    # sampling
    Xs = rand(g, ns)
    @test isa(Xs, Vector{Float64})
    @test length(Xs) == ns
    @test_approx_eq_eps mean(Xs) mean(g) 0.01
end

function test_mixture(g::MultivariateMixture, n::Int, ns::Int)
    X = zeros(length(g), n)
    for i = 1:10
        X[:,i] = rand(g)
    end

    cs = components(g)
    pr = probs(g)
    @assert length(cs) == length(pr)

    # mean
    mu = 0.0
    for k = 1:length(cs)
        mu += pr[k] * mean(cs[k])
    end
    @test_approx_eq mean(g) mu

    # ground-truth
    p0 = zeros(n)
    for i = 1:n
        p0_i = 0.0
        x_i = X[:,i]
        for k = 1:length(cs)
            p0_i += pr[k] * pdf(cs[k], x_i)
        end
        p0[i] = p0_i
    end
    lp0 = log(p0)

    for i = 1:n
        @test_approx_eq pdf(g, X[:,i]) p0[i]
        @test_approx_eq logpdf(g, X[:,i]) lp0[i]
    end

    p_e = pdf(g, X)
    lp_e = logpdf(g, X)
    @assert isa(p_e, Vector{Float64}) && length(p_e) == n
    @assert isa(lp_e, Vector{Float64}) && length(lp_e) == n

    @test_approx_eq p_e p0
    @test_approx_eq lp_e lp0

    # sampling
    Xs = rand(g, ns)
    @test isa(Xs, Matrix{Float64})
    @test size(Xs) == (length(g), ns)
    @test_approx_eq_eps vec(mean(Xs, 2)) mean(g) 0.01
end

# Tests

println("    testing UnivariateMixture")

g_u = MixtureModel(Normal, [(0.0, 1.0), (2.0, 1.0), (-4.0, 1.5)], [0.2, 0.5, 0.3])
@test isa(g_u, MixtureModel{Univariate, Continuous, Normal})
@test length(components(g_u)) == 3
test_mixture(g_u, 1000, 10^6)

println("    testing MultivariateMixture")
g_m = MixtureModel(
    IsoNormal[ MvNormal([0.0, 0.0], 1.0),
               MvNormal([0.2, 1.0], 1.0),
               MvNormal([-0.5, -3.0], 1.6) ],
    [0.2, 0.5, 0.3])
@test isa(g_m, MixtureModel{Multivariate, Continuous, IsoNormal})
@test length(components(g_m)) == 3
@test length(g_m) == 2
test_mixture(g_m, 1000, 10^6)


text:
# ---
# title: 1477. Find Two Non-overlapping Sub-arrays Each With Target Sum
# id: problem1477
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/find-two-non-overlapping-sub-arrays-each-with-target-sum/description/>
# hidden: true
# ---
#
# Given an array of integers `arr` and an integer `target`.
#
# You have to find **two non-overlapping sub-arrays** of `arr` each with sum
# equal `target`. There can be multiple answers so you have to find an answer
# where the sum of the lengths of the two sub-arrays is **minimum**.
#
# Return _the minimum sum of the lengths_ of the two required sub-arrays, or
# return _**-1**_ if you cannot find such two sub-arrays.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [3,2,2,4,3], target = 3
#     Output: 2
#     Explanation: Only two sub-arrays have sum = 3 ([3] and [3]). The sum of their lengths is 2.
#
#
# **Example 2:**
#
#
#
#     Input: arr = [7,3,4,7], target = 7
#     Output: 2
#     Explanation: Although we have three non-overlapping sub-arrays of sum = 7 ([7], [3,4] and [7]), but we will choose the first and third sub-arrays as the sum of their lengths is 2.
#
#
# **Example 3:**
#
#
#
#     Input: arr = [4,3,2,6,2,3,4], target = 6
#     Output: -1
#     Explanation: We have only one sub-array of sum = 6.
#
#
# **Example 4:**
#
#
#
#     Input: arr = [5,5,4,4,5], target = 3
#     Output: -1
#     Explanation: We cannot find a sub-array of sum = 3.
#
#
# **Example 5:**
#
#
#
#     Input: arr = [3,1,1,1,5,1,2,1], target = 3
#     Output: 3
#     Explanation: Note that sub-arrays [1,2] and [2,1] cannot be an answer because they overlap.
#
#
#
#
# **Constraints:**
#
#   * `1 <= arr.length <= 10^5`
#   * `1 <= arr[i] <= 1000`
#   * `1 <= target <= 10^8`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module NEPTypes
    using ..NEPCore
    using SparseArrays
    using LinearAlgebra
    using InteractiveUtils

    # Specializalized NEPs
    export ProjectableNEP
    export DEP
    export PEP
    export REP
    export SPMF_NEP

    export AbstractSPMF
    export SumNEP, SPMFSumNEP, GenericSumNEP
    export Proj_NEP
    export Proj_SPMF_NEP
    export DerSPMF

    export create_proj_NEP
    export get_Av
    export get_fv

    export interpolate
    export set_projectmatrices!
    export expand_projectmatrices!

    export compute_rf

    # We overload these
    import ..NEPCore.compute_Mder
    import ..NEPCore.compute_Mlincomb
    import ..NEPCore.compute_Mlincomb!
    import ..NEPCore.compute_MM
    import ..NEPCore.compute_resnorm
    import Base.size
    import SparseArrays.issparse

    # This is workaround for Julia bug https://github.com/JuliaLang/julia/issues/29224
    # This can be removed if the fix is backported in julia 1.0.1
    import Base.*
    function (*)(A::SubArray{<:Complex},B::Matrix{<:Real})
        return copy(A)*B;
    end

    "Returns the greatest type of all array elements."
    promote_typeof(A::AbstractArray) = isempty(A) ? eltype(A) : mapreduce(typeof, promote_type, A)

    "Returns the greatest element type of all array elements."
    promote_eltype(A::AbstractArray{<:AbstractArray}) = isempty(A) ? eltype(eltype(A)) : mapreduce(eltype, promote_type, A)

    #
"""
    abstract ProjectableNEP <: NEP

A ProjectableNEP is a NEP which can be projected, i.e., one can construct the problem ``W'*M(λ)Vw=0`` with the [`Proj_NEP`](@ref).
A NEP which is of this must have the function [`create_proj_NEP(orgnep::ProjectableNEP)`](@ref) implemented.
This function must return a `Proj_NEP`.

 See also: [`set_projectmatrices!`](@ref).

# Example:
```julia-repl
julia> nep=nep_gallery("dep0");
julia> typeof(nep)
DEP{Float64,Array{Float64,2}}
julia> isa(nep,ProjectableNEP)
true
julia> projnep=create_proj_NEP(nep);
julia> e1 = Matrix(1.0*I,size(nep,1),1);
julia> set_projectmatrices!(projnep,e1,e1);
julia> compute_Mder(nep,3.0)[1,1]
-2.942777908030041
julia> compute_Mder(projnep,3.0)
1×1 Array{Complex{Float64},2}:
 -2.942777908030041 + 0.0im
```
"""
    abstract type ProjectableNEP <: NEP end

    #######################################################
    ### Sum of products matrices and functions

"""
    abstract  AbstractSPMF <: ProjectableNEP

An AbstractSPMF is an abstract class representing NEPs which can be represented
as a sum of products of matrices and functions ``M(λ)=Σ_i A_i f_i(λ)``,
where i = 0,1,2,..., all of the matrices are of size ``n×n`` and ``f_i`` are functions.

Any AbstractSPMF has to have implementations of [`get_Av()`](@ref) and [`get_fv()`](@ref) which return the
functions and matrices.
"""
    abstract  type AbstractSPMF{T} <: ProjectableNEP end # See issue #17

"""
    get_Av(nep::AbstractSPMF)

Returns an array of matrices ``A_i`` in the AbstractSPMF: ``M(λ)=Σ_i A_i f_i(λ)``
"""
    function get_Av(nep::AbstractSPMF) # Dummy function which enforces that you have to implement
        error("You need to implement get_Av for all AbstractSPMFs")
    end
"""
    get_Av(nep::AbstractSPMF)

Returns an Array of functions (that can be evaluated both as scalar and matrix functions) ``f_i`` in the AbstractSPMF: ``M(λ)=Σ_i A_i f_i(λ)``
"""
    function get_fv(nep::AbstractSPMF) # Dummy function which enforces that you have to implement
        error("You need to implement get_fv for all AbstractSPMFs")
    end

"""
    struct SPMF_NEP{T<:AbstractMatrix,Ftype}  <: AbstractSPMF{T}
    function SPMF_NEP(AA, fii [,check_consistency=true] [,Schur_fact = false]
                      [,align_sparsity_patterns = false] [,Ftype=ComplexF64])

An `SPMF_NEP` is a `NEP` defined by a *S*um of *P*roducts of
*M*atrices and *F*unctions,
i.e.,
```math
M(λ)=∑_i A_i f_i(λ).
```
All of the matrices ``A_0,...`` are of size ``n×n``
and ``f_i`` are a functions. The  functions ``f_i`` must be defined
for matrices in the standard matrix function sense.
The constructor creates a `SPMF_NEP` consisting
of matrices `AA` and functions `fii`.

# Parameters

* `AA` is a `Vector` of matrices. The matrices have to be of the same type. If you need a NEP with different types you can use [`SumNEP`](@ref) to construct a sum of two `SPMF_NEP`.

* `fii` is a `Vector` of functions. Each function takes one parameter `S`. The functions must be available both as a scalar valid function and a matrix function. If `S` is a square matrix, `fii[k](S)` musst also be a square matrix. If `S` is a scalar `fii[k](S)` is a scalar.

* `check_consistency` (default `true`) determines if we should initiate by running tests to verify that the `fii` satisfies the conditions that every function is valid both for matrices and scalars. This is done by using `@code_typed` and the functions need to be type-stable in that sense.

* `align_sparsity_patterns` (default `false`) has effect only for sparse matrices (`SparseMatrixCSC`). If `align_sparsity_patterns=true` the `SparseMatrixCSC` matrices will be replaced by equivalent `SparseMatrixCSC` matrices where the `colptr` and `rowval` are identical. This increases the speed of some functions, e.g., `compute_Mder`. If `align_sparsity_patterns=true` the matrices in the NEP should be considered read only. If the sparsity patterns are completely or mostly distinct, it may be more efficient to set this flag to false.

* `Ftype` (default `ComplexF64`) determines an underlying type of the functions. The output of any function should be "smaller" than the promoted type of the input and `Ftype`. More precisely, if `F=fii[k]`, then the type logic is as follows `eltype(F(λ))=promote_type(eltype(λ),Ftype)`.

* `Schur_fact` (default `false`) determines if the `compute_MM` function should triangularize the matrix before carrying out the computation. This can be faster for large matrices.

# Example
```julia-repl
julia> A0=[1 3; 4 5]; A1=[3 4; 5 6];
julia> id_op=S -> one(S) # Note: We use one(S) to be valid both for matrices and scalars
julia> exp_op=S -> exp(S)
julia> nep=SPMF_NEP([A0,A1],[id_op,exp_op]);
julia> compute_Mder(nep,1)-(A0+A1*exp(1))
2×2 Array{Float64,2}:
 0.0  0.0
 0.0  0.0
```

"""
    struct SPMF_NEP{T<:AbstractMatrix,Ftype}  <: AbstractSPMF{T}
    # Logic behind Ftype:
    #  The eltype(F(λ))=promote_type(eltype(λ),Ftype)
        n::Int
        A::Vector{T}   # Array of Array of matrices
        fi::Vector{Function}  # Array of functions
        Schur_factorize_before::Bool # Tells if you want to do the Schur-factorization
        sparsity_patterns_aligned:: Bool
    end

    # Alias for SPMFs which contain sparse matrices
    SPMF_NEP_sparse{T<:AbstractSparseMatrix,Ftype}=SPMF_NEP{T,Ftype}

    function SPMF_NEP(AA::Vector{<:AbstractMatrix}, fii::Vector{<:Function};
                       Schur_fact = false,
                       check_consistency=true, Ftype=ComplexF64,
                       align_sparsity_patterns=false)

         matrices_are_sparse = issparse.(AA)
         if !(all(matrices_are_sparse) || all(.!matrices_are_sparse))
             error("Mixing sparse and dense matrices is not allowed in SPMF_NEP. Either use a consistent " *
                "format, or split your problem into a sparse and a dense part and use SumNEP.")
         end

         T=Float64;
         if (check_consistency)
             for t=1:length(fii)
                 # Scalar leads to scalars:
                 s=one(T);
                 ci=@code_typed(fii[t](s)) # ci[end] gives the return type
                 if !(ci[end] <: Number)
                     @warn "It seems you have not provided valid matrix-functions for defining SPMF_NEP. The functions fii should return a scalar if evaluated in a scalar and a matrix if evaluated in a matrix. If you want to disable to input checking, set check_consistency=false in SPMF_NEP."
                 end
                 S=ones(T,2,2);
                 ci=@code_typed(fii[t](S))
                 if !(ci[end] <: Matrix)
                     @warn "It seems you have not provided valid matrix-functions for defining SPMF_NEP. The functions fii should return a scalar if evaluated in a scalar and a matrix if evaluated in a matrix. If you want to disable to input checking, set check_consistency=false in SPMF_NEP."
                 end
             end
         end

         if (size(AA,1)==0)
             return SPMF_NEP(0); # Create empty SPMF_NEP.
         end
         n=size(AA[1],1);

         if(size(AA,1) != size(fii,1))
             error("Inconsistency: Number of supplied matrices = ", length(AA), " but the number of supplied functions are = ", length(fii))
         end
         for i = 2:size(AA,1)
             if size(AA[i]) != size(AA[1])
                 error("The dimensions of the matrices mismatch: size(AA[1])",
                       size(AA[1]),"!=",size(AA[i]),"=size(AA[",i,"])")
             end
         end

         if !(eltype(AA) <: SparseMatrixCSC)
             # Dense
             this=SPMF_NEP{promote_typeof(AA),Ftype}(n,AA,fii,Schur_fact,false);
         else
             # Sparse: Potentially do the joint sparsity pattern trick.
             if (!align_sparsity_patterns)
                 # No aligning, just create it
                 this=SPMF_NEP{promote_typeof(AA),Ftype}(n,AA,fii,Schur_fact,false);
             else
                 TT = eltype(AA[1])
                 As=form_aligned_sparsity_patterns(AA,TT)
                 this=SPMF_NEP{promote_typeof(As),Ftype}(n,As,fii,Schur_fact,true);
             end
         end
         return this
    end
    function SPMF_NEP(n) # Create an empty NEP of size n x n
        Z=zeros(n,n)
        return SPMF_NEP{AbstractMatrix,Complex}(n,Vector{Matrix}(),Vector{Function}(),false);
    end

""" Return a vector of sparse matrices which are the same as AA but also have the same sparsity pattern. """
    function form_aligned_sparsity_patterns(AA::Vector{<:SparseMatrixCSC},T)

        # Create a Zero matrix with the correct sparsity pattern.
        Zero = LinearAlgebra.fillstored!(copy(AA[1]), 1)
        for i = 2:size(AA,1)
            Zero .+= LinearAlgebra.fillstored!(copy(AA[i]), 1)
        end
        Zero = T.(Zero)
        Zero.nzval[:] .= T(0) # Set the nonzero elements to zero (keep pattern)

        # Create an vector of sparse matrices. Set the elements corresponding to
        # AA[i] in each of them, such that As[i]=AA[i] but with different sparsity patterns.
        As = Vector{SparseMatrixCSC{T,Int}}()
        @inbounds for A in AA
            # Create a new zero matrix. Note all S-matrices have the same colptr, so
            # modification of As[i] will change all. The As[i] matrices should therefore
            # be viewed as read-only objects.
            S = SparseMatrixCSC(Zero.m, Zero.n,
                                Zero.colptr, Zero.rowval,
                                fill(zero(T),size(Zero.nzval,1)))
            # Set the nonzero elements of A in S
            for col = 1:size(A, 2)
                for j in nzrange(A, col)
                    S[A.rowval[j], col] = A.nzval[j]
                end
            end
            push!(As, S)
        end
        return As;
    end

    function compute_MM(nep::SPMF_NEP{T,Ftype},S,V) where {T,Ftype}

        AA=get_Av(nep);
        ff=get_fv(nep);
        m=size(ff,1); n=size(nep,1); p=size(S,1);

        # Type logic including Ftype
        FStype=promote_type(eltype(S),Ftype) # eltype of f(S)
        T0=promote_type(eltype(V),eltype(AA[1]),FStype) # Output type

        # Always return a dense matrix
        Z=zeros(T0,n,p);

        # Precompute schur factorization (temporarily disabled)
        #if (nep.Schur_factorize_before) #Optimize if allowed to factorize before
        #    (TT, Q, ) = schur(S)  # Currently not used
        #end

        VFi=Matrix{promote_type(FStype,eltype(V))}(undef,n,p); # Type of V*f(S)
        local Fi=Matrix{FStype}(undef,p,p);
        for i=1:m
            ## Compute Fi=f_i(S) in an optimized way
            if (isdiag(S)) # optimize if S is diagonal
                Sd=diag(S);
                local Fid::Vector{FStype}
                if (norm(Sd .- Sd[1])==0) # Optimize further if S is a multiple of identity
                    Fid=fill(ff[i](Sd[1]),p)
                else  # Diagonal but not constant
                    Fid0=Vector{Number}(undef,p)
                    for j=1:p
                        Fid0[j]=ff[i](Sd[j])
                    end
                    Fid=Vector{FStype}(Fid0);
                end
                Fi .= Diagonal(Fid)
            else
                Fi .= ff[i](S);
            end
            mul!(VFi,V,Fi);
            Z[:,:] .+= AA[i]*VFi;
        end
        return Z
    end

""" Internal helper function: returns (TZ,x) where x is a vector of function evaluation of nep.fi[i](λ). TZ is greatest of the types in vector x."""
    function compute_Mder_fi_and_output_type(nep::AbstractSPMF,λ::Number)
        ff=get_fv(nep);
        AA=get_Av(nep)
        x = [f(λ) for f in ff]

        # figure out the return type, as the greatest type of all input
        Tx = promote_typeof(x)
        TA = promote_eltype(AA) # Greatest type of all A-matrices
        TZ = promote_type(TA,Tx)  # output type
        return (TZ,x);
    end

""" Internal helper function: Computes a linear combination with coeffs x of the A-matrices in an
SPMF. The result will be stored in an  AbstractMatrix with eltype T.  """
    function compute_linear_combination_SPMF(nep::AbstractSPMF,::Type{T},x::Vector) where {T<:Number}
        AA=get_Av(nep);
        # Fall back version just sum it up, e.g., for dense matrices
        Z=mapreduce(i->AA[i]*x[i],+,1:length(AA))
        #@assert eltype(Z) <: T
    end
    # Specialized function sparse matrices, where we can exploit sparsity
    function compute_linear_combination_SPMF(nep::SPMF_NEP_sparse,::Type{T},x::Vector) where {T<:Number}
        local Z::SparseMatrixCSC{T,Int}
        if (nep.sparsity_patterns_aligned)
             # Sparsity patterns are aligned, so just change
             # the value entries in the sparse matrix object
             Z = SparseMatrixCSC(nep.A[1].m, nep.A[1].n,
                                 copy(nep.A[1].colptr), copy(nep.A[1].rowval),
                                 copy(convert.(T, nep.A[1].nzval .* x[1])))
             for k = 2:length(nep.A)
                 Z.nzval .+= nep.A[k].nzval .* x[k]
             end
         else
             AA=nep.A;
             # No alignment of sparsity pattern. Naive summing.
             Z=mapreduce(i->AA[i]*x[i],+,1:length(AA))
         end
         return Z
    end

    function compute_Mder(nep::AbstractSPMF,λ::Number)
        TZ,x = compute_Mder_fi_and_output_type(nep,λ)
        # Try to do "smart" summing
        Z = compute_linear_combination_SPMF(nep,TZ,x);
        return Z
     end

    # For higher derivatives
    function compute_Mder(nep::SPMF_NEP{T,Ftype},λ::Number,i::Integer) where {T,Ftype}
        if (i==0)
            return compute_Mder(nep,λ);
        else
            n=size(nep,1);

            # Jordan matrix trick
            k=i+1;
       	    S=diagm(0 => fill(λ,k), -1 => (1:k-1))
            TS=eltype(S);

            # Type logic
            Fλtype=promote_type(TS,Ftype);
            TT=promote_type(Fλtype,eltype(nep.A[1])); # Return type

            # Compute the derivatives of the individual functions with the jordan matrix trick
            fk=Vector{TT}(undef,size(nep.A,1))
            for j=1:size(nep.A,1)
                fk[j] = nep.fi[j](S)[end,1]; # Contains the derivative
            end

            # Try to do "smart" summing
            return compute_linear_combination_SPMF(nep,TT,fk);
        end
    end

    ###########################################################
    # Delay eigenvalue problems - DEP
    #

    """
    type DEP <: AbstractSPMF
    function DEP(AA::Vector{AbstractMatrix} [,tauv::Vector=[0,1.0]])

A `DEP` (Delay Eigenvalue problem) is a problem defined by
the sum
```math
M(λ)=-λI + Σ_i A_i exp(-τ_i λ)
```
where all of the matrices are of size ``n×n``. This type of
NEP describes the stability of time-delay systems.

The construction takes the system matrices ``A_i``, and `tauv` is a vector of the values  ``τ_i``.

# Example:
```julia-repl
julia> A0=randn(3,3); A1=randn(3,3);
julia> tauv=[0,0.2] # Vector with delays
julia> dep=DEP([A0,A1],tauv)
julia> λ=3.0;
julia> M1=compute_Mder(dep,λ)
julia> M2=-λ*I+A0+A1*exp(-tauv[2]*λ)
julia> norm(M1-M2)
0.0
```
"""
    struct DEP{Z<:Real, T<:AbstractMatrix} <: AbstractSPMF{T}
        n::Int
        A::Array{T,1}     # An array of matrices (full or sparse matrices)
        tauv::Vector{Z}   # the delays (which are always real)
    end
    function DEP(AA::Vector{T},tauv::Vector=[0,1.0]) where {T<:AbstractMatrix}
        n=size(AA[1],1)
        if (!isreal(tauv))
            error("Incorrect construction of DEP. The delays need to be real.")
        end

        # Note: this enforces that eltype(tauv)=real(eltype(A[1]))
        tauvconv::Vector{real(eltype(AA[1]))}=Vector{real(eltype(AA[1]))}(tauv);

        this=DEP{eltype(tauvconv),T}(n,AA,tauvconv);
        return this;
    end

# Compute the ith derivative of a DEP
    function compute_Mder(nep::DEP,λ::Number,i::Integer=0)
        local M, J
        # T is eltype(nep.A[1]) unless λ complex, then T is complex(eltype(nep.A[1]))
        # T can be determined compile time, since DEP parametric type
        T=isa(λ,Complex) ? complex(eltype(nep.A[1])) : eltype(nep.A[1]);

        TM=promote_type(eltype(λ),T); # output type
        if (issparse(nep.A[1])) # Can be determined compiled time since DEP parametric type
            M=spzeros(TM,nep.n,nep.n)
            J=sparse(TM(1)I, nep.n, nep.n)
        else
            M=zeros(TM,nep.n,nep.n)
            J=Matrix{TM}(I, nep.n, nep.n)
        end
        if i==0; M=-λ*J;  end
        if i==1; M=-one(TM)*J; end
        for j=1:size(nep.A,1)
            a=exp(-nep.tauv[j]*λ)*(-nep.tauv[j])^i;
            M += nep.A[j]*a
        end
        return M
    end

#  Computes the sum ``Σ_i M_i V f_i(S)`` for a DEP
    function compute_MM(nep::DEP,S,V)

        T1=promote_type(promote_type(eltype(S),eltype(V)),eltype(nep.A[1]));
        T=promote_type(T1,eltype(nep.tauv));

        Z::Matrix{T}=-V*S;
        for j=1:size(nep.A,1)
            Z+=nep.A[j]*V*exp(Matrix(-nep.tauv[j]*S))
        end
        return Z
    end

    #  Fetch the Av's, since they are not explicitly stored in DEPs
    function get_Av(nep::DEP)
        local J
        if issparse(nep)
            J = sparse(eltype(nep.A[1])(1)I, nep.n, nep.n)
        else
            J = Matrix{eltype(nep.A[1])}(I, nep.n, nep.n)
        end
        return [J, nep.A...]
    end
    #  Fetch the Fv's, since they are not explicitly stored in DEPs
    function get_fv(nep::DEP)
        fv = Array{Function,1}(undef, size(nep.A,1)+1)
        # First function is -λ
        fv[1] = S -> -S

        # The other functions are exp(-tauv[j]*S)
        for i=1:size(nep.A,1)
            if nep.tauv[i] == 0 # Zero delay means constant term
                fv[i+1] = S -> one(S)
            else
                # introducing tau is needed to define type stable functions
                tau=nep.tauv[i] #delay of the i-th term
                fv[i+1] = S -> exp(-tau*S)
            end
        end

        return fv
    end

    include("types_poly.jl");

    ###########################################################
    # Rational eigenvalue problem - REP
"""
    function REP(A,roots,poles)

A `REP`-call creates a rational eigenvalue problem. The `REP` is defined by the
sum ``Σ_i A_i s_i(λ)/q_i(λ)``, where i = 0,1,2,..., all of the
matrices are of size ``n×n`` and ``s_i`` and ``q_i`` are polynomials.
The constructor takes the roots and poles as input of polynomials with
normalized highest coefficient. The NEP is defined as
```math
-λI+A_0+A_1\\frac{p(λ)}{q(λ)}
```
where `p` has the roots `roots` and `q` has the roots `poles`.

# Example
```julia-repl
julia> A0=[1 2; 3 4]; A1=[3 4; 5 6];
julia> nep=REP([A0,A1],[1,3], [4,5,6]);
julia> compute_Mder(nep,3)
2×2 Array{Float64,2}:
 Inf  Inf
 Inf  Inf
julia> (λ,x)=quasinewton(nep,v=[1;0])
(-0.3689603779201249 + 0.0im, Complex{Float64}[-2.51824+0.0im, 1.71283+0.0im])
julia> -λ*x+A0*x+A1*x*(λ-1)*(λ-3)/((λ-4)*(λ-5)*(λ-6))
2-element Array{Complex{Float64},1}:
 -2.5055998942313806e-13 + 0.0im
   1.318944953254686e-13 + 0.0im
```
"""
    function REP(AA::Vector,roots::Vector,poles::Vector)

        n=size(AA[1],1)
        A=[one(AA[1]),AA[1],AA[2]];
        roots_float=float.(roots);
        poles_float=float.(poles);
        nep=SPMF_NEP(A,
                     [S-> -S,S->one(S),
                     S-> root_eval(S,poles_float)\root_eval(S,roots_float)]);

        return nep;
    end

    function root_eval(S,a::Vector)
        F=mapreduce(i-> S-a[i]*I, (S1,S2) -> S1*S2, 1:size(a,1))
        return F;
    end

    #######################################################
    ### Represents a projected NEP
"""
    abstract type Proj_NEP <: NEP

`Proj_NEP` represents a projected `NEP`. The projection is defined
as the NEP
```math
N(λ)=W^HM(λ)V
```
where ``M(λ)`` is a base NEP and `W` and `V` rectangular matrices representing
a basis of the projection spaces.
Instances are created with `create_proj_NEP`. See [`create_proj_NEP`](@ref)
for examples.

Any `Proj_NEP` needs to implement two functions to manipulate the projection:

* [`set_projectmatrices!`](@ref): Set matrices `W` and `V`
* [`expand_projectmatrices!`](@ref): Effectively expand the matrices `W` and `V` with one column.

"""
    abstract type Proj_NEP <: NEP end

"""
    pnep=create_proj_NEP(orgnep::ProjectableNEP[,maxsize [,T]])

Create a NEP representing a projected problem ``N(λ)=W^HM(λ)V``,
 where the base NEP is represented by `orgnep`.
The optional parameter `maxsize::Int` determines how large the projected
problem can be and `T` is the Number type used for the projection matrices
(default `ComplexF64`).
These are needed for memory preallocation reasons.
Use [`set_projectmatrices!`](@ref) and [`expand_projectmatrices!`](@ref)
 to specify projection matrices ``V`` and ``W``.

# Example:
The following example illustrates that a projection
of a `NEP` is also a `NEP` and we can for instance
call `compute_Mder` on it:
```julia-repl
julia> nep=nep_gallery("pep0")
julia> V=Matrix(1.0*I,size(nep,1),2);
julia> W=Matrix(1.0*I,size(nep,1),2);
julia> pnep=create_proj_NEP(nep);
julia> set_projectmatrices!(pnep,W,V);
julia> compute_Mder(pnep,3.0)
2×2 Array{Complex{Float64},2}:
  6.08082+0.0im  -5.47481+0.0im
 0.986559+0.0im  -6.98165+0.0im
julia> W'*compute_Mder(nep,3.0)*V  # Gives the same result
2×2 Array{Float64,2}:
 6.08082   -5.47481
 0.986559  -6.98165
```
If you know that you will only use real projection matrices, you
can specify this in at the creation:
```julia-repl
julia> pnep=create_proj_NEP(nep,2,Float64);
julia> set_projectmatrices!(pnep,W,V);
julia> compute_Mder(pnep,3.0)
2×2 Array{Float64,2}:
 6.08082   -5.47481
 0.986559  -6.98165
```
"""
    function create_proj_NEP(orgnep::ProjectableNEP)
         error("Not implemented. All ProjectableNEP have to implement create_proj_NEP.")
    end
    function create_proj_NEP(orgnep::AbstractSPMF,
                             maxsize::Int=min(size(orgnep,1),201),
                             T::Type{<:Number}=ComplexF64)
        return Proj_SPMF_NEP(orgnep,maxsize,T);
    end

"""
    struct Proj_SPMF_NEP <: Proj_NEP

This type represents the (generic) way to project NEPs which are
`AbstractSPMF`. See examples in [`create_proj_NEP`](@ref).
"""
    mutable struct Proj_SPMF_NEP  <: Proj_NEP
        orgnep::AbstractSPMF
        nep_proj::SPMF_NEP; # An instance of the projected NEP
        orgnep_Av::Vector
        orgnep_fv::Vector
        projnep_B_mem::Vector # A vector of matrices
        function Proj_SPMF_NEP(nep::AbstractSPMF,maxsize::Int,
                               subspace_eltype=ComplexF64)
            this = new(nep)

            this.orgnep_Av = get_Av(nep)
            if     (size(this.orgnep_Av,1) != 1) && (size(this.orgnep_Av,2) == 1) # Stored as column vector - do nothing
            elseif (size(this.orgnep_Av,1) == 1) && (size(this.orgnep_Av,2) == 1) # It is a single entry - do nothing
            elseif (size(this.orgnep_Av,1) == 1) && (size(this.orgnep_Av,2) != 1) # Stored as a row-vector
                this.orgnep_Av = vec(this.orgnep_Av)
            else
                error("The given array should be a vector but is of size ", size(this.orgnep_Av), ".")
            end

            this.orgnep_fv = get_fv(nep)
            if     (size(this.orgnep_fv,1) != 1) && (size(this.orgnep_fv,2) == 1) # Stored as column vector - do nothing
            elseif (size(this.orgnep_fv,1) == 1) && (size(this.orgnep_fv,2) == 1) # It is a single entry - do nothing
            elseif (size(this.orgnep_fv,1) == 1) && (size(this.orgnep_fv,2) != 1) # Stored as a row-vector
                this.orgnep_fv = vec(this.orgnep_fv)
            else
                error("The given array should be a vector but is of size ", size(this.orgnep_fv), ".")
            end

            # compute greatest eltype of Av:
            Av_eltype = mapreduce(eltype,promote_type,this.orgnep_Av);

            Bmat_eltype=promote_type(subspace_eltype,Av_eltype)

            # Construct memory for the projected problem matrices
            this.projnep_B_mem=Vector{Matrix{Bmat_eltype}}(undef,size(this.orgnep_fv,1));
            for k=1:size(this.orgnep_fv,1)
                this.projnep_B_mem[k]=zeros(Bmat_eltype,maxsize,maxsize);
            end

            # Reset the projectmatrices
            set_projectmatrices!(this,zeros(size(this.orgnep,1),0),zeros(size(this.orgnep,1),0))

            return this
        end
    end

"""
    set_projectmatrices!(pnep::Proj_NEP,W,V)
Set the projection matrices for the NEP to `W` and `V`, i.e.,
corresponding the NEP: ``N(λ)=W^HM(λ)V``. See also [`create_proj_NEP`](@ref).

# Example
This illustrates if `W` and `V` are vectors of ones, the projected problem
becomes the sum of the rows and columns of the original NEP.
```julia-repl
julia> nep=nep_gallery("pep0")
julia> pnep=create_proj_NEP(nep);
julia> V=ones(200,1);  W=ones(200,1);
julia> set_projectmatrices!(pnep,W,V);
julia> compute_Mder(pnep,0)
1×1 Array{Complex{Float64},2}:
 48.948104019482756 + 0.0im
julia> sum(compute_Mder(nep,0),dims=[1,2])
1×1 Array{Float64,2}:
 48.948104019482955
```

"""
    function set_projectmatrices!(nep::Proj_SPMF_NEP,W,V)
        ## Sets the left and right projected basis and computes
        ## the underlying projected NEP
        m=size(nep.orgnep_Av,1);
        k=size(V,2);
        @assert(k <= size(nep.projnep_B_mem[1],1)) # Don't go outside the prealloc memory
        # For over all i: Compute the expanded matrices
        WT=copy(W')
        B = map(i -> begin
                # Compute the projecte matrix
                nep.projnep_B_mem[i][1:k,1:k]=WT*nep.orgnep_Av[i]*V;
                view(nep.projnep_B_mem[i],1:k,1:k);
                end, 1:m)
        # Keep the sequence of functions for SPMFs
        nep.nep_proj=SPMF_NEP(B,nep.orgnep_fv,check_consistency=false)
    end

"""
    expand_projectmatrices!(nep::Proj_SPMF_NEP, Wnew, Vnew)

The projected NEP is updated by adding the last column of `Wnew` and `Vnew`
to the basis. Note that `Wnew` and `Vnew` contain also the "old" basis vectors.
See also [`create_proj_NEP`](@ref)

# Example:

In the following example you see that the expanded projected problem
has one row and column more, and the leading subblock is the same
as the smaller projected NEP.
```julia-repl
julia> nep=nep_gallery("pep0"); n=size(nep,1);
julia> V=Matrix(1.0*I,n,2); W=Matrix(1.0*I,n,2);
julia> pnep=create_proj_NEP(nep);
julia> set_projectmatrices!(pnep,W,V);
julia> compute_Mder(pnep,0)
2×2 Array{Complex{Float64},2}:
 0.679107+0.0im   -0.50376+0.0im
 0.828413+0.0im  0.0646768+0.0im
julia> Vnew=[V ones(n)]
julia> Wnew=[W ones(n)]
julia> expand_projectmatrices!(pnep,Wnew,Vnew);
julia> compute_Mder(pnep,0)
3×3 Array{Complex{Float64},2}:
 0.679107+0.0im   -0.50376+0.0im  -12.1418+0.0im
 0.828413+0.0im  0.0646768+0.0im   16.3126+0.0im
 -17.1619+0.0im   -10.1628+0.0im   48.9481+0.0im
```

"""
    function expand_projectmatrices!(nep::Proj_SPMF_NEP,Wnew::AbstractMatrix,Vnew::AbstractMatrix)
        w=Wnew[:,end];
        v=Vnew[:,end];
        Av=nep.orgnep_Av;
        k=size(Vnew, 2)-1;
        m = size(nep.orgnep_Av,1);
        @assert(k+1 <= size(nep.projnep_B_mem[1],1)) # Don't go outside the prealloc memory
        WnewT=copy(Wnew[:,1:k]');
        # For over all i: Compute expanded part of matrices:
        B = map(i -> begin
                # Expand the B-matrices
                nep.projnep_B_mem[i][1:k,k+1]=WnewT*Av[i]*v;
                nep.projnep_B_mem[i][k+1,1:(k+1)]=w'*Av[i]*view(Vnew,:,1:(k+1));
                view(nep.projnep_B_mem[i],1:(k+1),1:(k+1));
                end, 1:m)
        # Keep the sequence of functions for SPMFs
        nep.nep_proj=SPMF_NEP(B,nep.orgnep_fv,check_consistency=false)
    end

    # Use delagation to the nep_proj
    compute_MM(nep::Union{Proj_SPMF_NEP},par...)=compute_MM(nep.nep_proj,par...)
    # Use MM to compute Mlincomb for SPMFs
    compute_Mlincomb(nep::Proj_SPMF_NEP,λ::Number,V::AbstractVecOrMat)=
             compute_Mlincomb(nep.nep_proj,λ,V)
    compute_Mlincomb(nep::Proj_SPMF_NEP,λ::Number,V::AbstractVecOrMat,a::Vector)=
             compute_Mlincomb(nep.nep_proj,λ,V,a)
    compute_Mder(nep::Union{Proj_SPMF_NEP},λ::Number)=compute_Mder(nep.nep_proj,λ,0)
    compute_Mder(nep::Union{Proj_SPMF_NEP},λ::Number,i::Integer)=compute_Mder(nep.nep_proj,λ,i)

    function size(nep::Proj_NEP,dim)
        n = size(nep.nep_proj,2);
        return n
    end
    function size(nep::Proj_NEP)
        n = size(nep.nep_proj,2);
        return (n,n)
    end

    function get_Av(nep::Proj_SPMF_NEP)
        return get_Av(nep.nep_proj);
    end

    function get_fv(nep::Proj_SPMF_NEP)
        return nep.orgnep_fv
    end

    function issparse(nep::Proj_NEP)
        return false; # A projected NEP is (essentially) never sparse
    end

"""
    struct GenericSumNEP{NEP1<:NEP,NEP2<:NEP}  <: NEP

See also: [`SumNEP`](@ref), [`SPMFSumNEP`](@ref)
"""
    struct GenericSumNEP{NEP1<:NEP,NEP2<:NEP}  <: NEP
        nep1::NEP1
        nep2::NEP2
    end

"""
    struct SPMFSumNEP{NEP1<:AbstractSPMF,NEP2<:AbstractSPMF}  <: AbstractSPMF{AbstractMatrix}

See also: [`SumNEP`](@ref), [`GenericSumNEP`](@ref)
"""
    struct SPMFSumNEP{NEP1<:AbstractSPMF,NEP2<:AbstractSPMF}  <: AbstractSPMF{AbstractMatrix}
        nep1::NEP1
        nep2::NEP2
    end

    AnySumNEP=Union{GenericSumNEP,SPMFSumNEP};
    # Creator functions for SumNEP
    function SumNEP(nep1::AbstractSPMF,nep2::AbstractSPMF)
        return SPMFSumNEP(nep1,nep2);
    end
    """
    SumNEP{nep1::NEP,nep2::NEP}
    SumNEP{nep1::AbstractSPMF,nep2::AbstractSPMF}

`SumNEP` is a function creating an object that corresponds to a sum of two NEPs,
i.e., if nep is created by `SumNEP` it is defined by
```math
M(λ)=M_1(λ)+M_2(λ)
```
where ``M_1`` and ``M_2`` are defined by `nep1` and `nep2`.

# Example:
```julia-repl
julia> nep1=DEP([ones(3,3),randn(3,3)])
julia> nep2=PEP([ones(3,3),randn(3,3),randn(3,3)])
julia> sumnep=SumNEP(nep1,nep2);
julia> s=3.0;
julia> M=compute_Mder(sumnep,s);
3×3 Array{Float64,2}:
  8.54014     6.71897   7.12007
 -0.943908  -13.0795   -0.621659
  6.03155    -7.26726  -6.42828
julia> M1=compute_Mder(nep1,s);
julia> M2=compute_Mder(nep2,s);
julia> M1+M2  # Same as M
3×3 Array{Float64,2}:
  8.54014     6.71897   7.12007
 -0.943908  -13.0795   -0.621659
  6.03155    -7.26726  -6.42828
```

See also: [`SPMFSumNEP`](@ref), [`GenericSumNEP`](@ref)

"""
    function SumNEP(nep1::NEP,nep2::NEP)
        return GenericSumNEP(nep1,nep2);
    end

    # Delegate all interface functions
    size(snep::AnySumNEP)=size(snep.nep1)
    size(snep::AnySumNEP,d)=size(snep.nep1,d)
    compute_Mlincomb(nep::AnySumNEP, λ::Number, V::AbstractVecOrMat) =
        (compute_Mlincomb(nep.nep1, λ, V)+compute_Mlincomb(nep.nep2,λ,V))
    compute_Mder(nep::AnySumNEP, λ::Number,i::Int = 0) =
        (compute_Mder(nep.nep1,λ,i)+compute_Mder(nep.nep2,λ,i))
    compute_MM(nep::AnySumNEP, S::AbstractMatrix,V::AbstractMatrix) =
        (compute_MM(nep.nep1,S,V)+compute_MM(nep.nep2,S,V))

    # For SPMFSumNEP, also delegate the get_Av() and get_fv()
    get_Av(nep::SPMFSumNEP) = [get_Av(nep.nep1); get_Av(nep.nep2)]
    get_fv(nep::SPMFSumNEP) = [get_fv(nep.nep1); get_fv(nep.nep2)]

   #######################################################
   ### Functions in common for many NEPs in NEPTypes

   #
"""
    size(nep::Union{DEP,PEP,SPMF_NEP})
    size(nep::Union{DEP,PEP,SPMF_NEP},dim)

Overloads the size functions for NEPs storing size in `nep.n`
"""
    function size(nep::Union{DEP,PEP,SPMF_NEP})
        return (nep.n,nep.n)
    end
    function size(nep::Union{DEP,PEP,SPMF_NEP},dim)
        return nep.n
    end

"""
    issparse(nep)
Returns true/false if the NEP is sparse (if `compute_Mder()` returns sparse)
"""
    function issparse(nep::Union{DEP,PEP,SPMF_NEP})
        return issparse(nep.A[1])
    end

    function get_Av(nep::Union{SPMF_NEP,PEP})
        return nep.A;
    end
    function get_fv(nep::SPMF_NEP)
        return nep.fi;
    end

    include("nep_deflation.jl")
    include("low_rank_nep.jl")
    include("errmeasure.jl")

    # structure exploitation for DEP
    function compute_Mlincomb(nep::DEP,λ::Number,V::AbstractVecOrMat,
                              a::Vector=ones(eltype(V),size(V,2)))
        n=size(V,1); k=size(V,2);
        # Type logic
        TT=promote_type(eltype(V),typeof(λ),eltype(nep.A[1]),eltype(nep.tauv),eltype(a))

        # initialize variables
        z=zeros(TT,n); Vw = Vector{TT}(undef, n); AVw = Vector{TT}(undef, n)
        # with a direct computation one can see that
        # z=-λV[:,1]-V[:,2]+\sum_{j=1}^{length(nep.tauv)} nep.Av[j+1] (V w)
        # where w is the vector with the scaled delays
        for j=1:length(nep.tauv)
            w=Array{TT,1}(exp(-λ*nep.tauv[j])*(-nep.tauv[j]) .^(0:k-1))
            mul!(Vw, V, a.*w)
            mul!(AVw, nep.A[j], Vw)
            z[:] .+= AVw;
        end

        # distinguis the case V is a vector and V is a matrix
        # fix with the proper derivative count
        if (V isa AbstractVector)
            z .-= a[1]*λ*V
        elseif k==1
            z .-= a[1]*λ*V[:]
        else
            z .+= muladd(-λ*a[1],V[:,1],-a[2]*V[:,2])
        end
        return z
    end

    compute_Mlincomb!(nep::DEP,λ::Number,V::AbstractVecOrMat, a::Vector=ones(size(V,2)))=compute_Mlincomb(nep,λ,V, a)

    function compute_Mlincomb!(nep::SPMF_NEP{T,Ftype},
                               λ::Number,
                               V::AbstractVecOrMat,
                               a::Vector=ones(size(V,2))) where {T,Ftype}

        local n,k;
        n=size(V,1);
        k=size(V,2);

    	# we need to assume that the elements of a are different than zero.
    	V[:,findall(x->x==0,a)] .= 0
    	a[findall(x->x==0,a)] .= 1
        local S,TS;
        if (V isa AbstractVector)
            #Vector means just compute matrix vector
            S=λ
            TS = eltype(λ)
        else
            # V matrix means compute linear combination of derivatives. Use
            # scaling trick
            TS= promote_type(typeof(λ),eltype(a));
       	    S=diagm(0 => fill(λ,k), -1 => (a[2:k]./a[1:k-1]).*(1:k-1))
        end

        # Type logic
        Fλtype=promote_type(TS,Ftype);
        TT=promote_type(Fλtype,eltype(V)); # Return type

        z=zeros(TT,n)
        for i=1:size(nep.A,1)
            # Get the function value if V is a vector,
            # otherwise get a vector of scaled derivatives
            Fi1=(V isa AbstractVector) ? nep.fi[i](S) : nep.fi[i](S)[:,1]

            VFi1=V*Fi1
            z .+= nep.A[i]*VFi1
        end

    	return a[1]*reshape(z,n);
    end

    compute_Mlincomb(nep::SPMF_NEP,λ::Number,V::AbstractVecOrMat, a::Vector=ones(size(V,2)))=compute_Mlincomb!(nep,λ,copy(V), copy(a))

    function compute_Mlincomb(
                        nep::PEP,
                        λ::Number,
                        V::AbstractVecOrMat,
                        a::Vector=ones(eltype(V),size(V,2)))

        # Type logic
        TT=promote_type(typeof(λ),eltype(V),eltype(a),eltype(nep.A[1]));

        n=size(nep,1)
        z=zeros(TT,n)
        d=length(nep.A)-1
        k=min(size(V,2),d+1)
        tmp = Vector{TT}(undef, n)

        if iszero(λ)
            for j=0:k-1
                mul!(tmp, nep.A[j+1], view(V,:,j+1))
                z .+= a[j+1] .* factorial(j) .* tmp
            end
        else
            for j=0:k-1
                for i=j:d
                    mul!(tmp, nep.A[i+1], view(V,:,j+1))
                    z .+= a[j+1] .* λ^(i-j) .* (factorial(i)/factorial(i-j)) .* tmp
                end
            end
        end
        return z
    end

    ########## THE FOLLOWING TYPE DerSPMF AND FUNCTION WILL BE DOCUMENTED AND THE CODE BETTER ORGANIZED ##########
    """
        struct DerSPMF{T<:AbstractMatrix,FDtype,TT<:Number} <: AbstractSPMF{T}

    A DerSPMF is a representation of a NEP defined by a Sum of Products of Matrices and Functions [`SPMF_NEP`](@ref). This format makes the execution of [`compute_Mlincomb`](@ref) for the specified `σ` more efficient.
    """
    struct DerSPMF{T<:AbstractMatrix,TT<:Number,FDtype} <: AbstractSPMF{T}
        spmf::AbstractSPMF{T}
        σ::TT
        fD::Matrix{FDtype}
    end

    # implement all the functions
    function size(nep::DerSPMF)
        return size(nep.spmf)
    end
    function size(nep::DerSPMF,dim)
        return size(nep.spmf,dim)
    end
    function get_fv(nep::DerSPMF)
        return get_fv(nep.spmf)
    end
    function get_Av(nep::DerSPMF)
        return get_Av(nep.spmf)
    end

"""
    newspmf=DerSPMF(spmf,σ,m)

Creates a `DerSPMF` representing the NEP `spmf` where the first `m` derivatives of the functions `f_i` in the number
`σ` are precomputed. This will in general speed up the execution of [`compute_Mlincomb`](@ref) for the
selected shift `σ`.

# Parameters

* `spmf` is the original `AbstractSPMF`.

* `σ::Number` specifies where the derivatives will be precomputed.

# Example
```julia-repl
julia> A0=[1 3; 4 5]; A1=[3 4; 5 6];
julia> id_op=S -> one(S) # Note: We use one(S) to be valid both for matrices and scalars
julia> exp_op=S -> exp(S)
julia> nep=SPMF_NEP([A0,A1],[id_op,exp_op]);
julia> m=5 # Precompute 5 derivatives
julia> σ=3.3
julia> dnep=DerSPMF(nep,σ,m)
julia> V=rand(2,m)
julia> z=compute_Mlincomb(dnep,σ,V)  # Same as below ...
2-element Array{Float64,1}:
39.47327945131233
64.31391434063991
julia> z=compute_Mlincomb(nep,σ,V)  # .. but more efficient
2-element Array{Complex{Float64},1}:
 381.008192939787 + 0.0im
 601.379183090249 + 0.0im
```
"""
    function DerSPMF(spmf::AbstractSPMF,σ::Number,m::Int)
          # Compute DD-matrix from get_fv(spmf)

          Av=get_Av(spmf)
          fv=get_fv(spmf)

          TT=promote_type(typeof(σ),eltype(Av[1]))

          # test if the functions fv introduce a super type
          SS=diagm(0=> σ*ones(TT,2m+2),  -1 => (1:2m+1))
          for t=1:length(fv)
              ci=@code_typed(fv[t](SS[1])); TT=promote_type(eltype(ci[end]),TT);
          end

          p=length(fv)
          # matrix for the computation of derivatives

          fD=Matrix{TT}(undef, 2*m+2,p)
          for t=1:p fD[:,t]=fv[t](SS)[:,1] end
          return DerSPMF(spmf,σ,fD);
    end

    function compute_Mlincomb(
                        nep::DerSPMF{T,FDtype},
                        λ::Number,
                        V::AbstractVecOrMat,
                        a::Vector=ones(size(V,2))) where {T,FDtype}

        if λ!=nep.σ
            return compute_Mlincomb(nep.spmf,λ,V,a)
        else

            local n,k,p
            p=size(nep.fD,2)

            if (V isa AbstractVector)
                k=1; n=length(V)
            else
                n,k=size(V)
            end

            Av=get_Av(nep)
            # Type logic
            TT=promote_type(eltype(V),typeof(λ),eltype(Av[1]),eltype(a))
            z=zeros(TT,n)
            VafD=V*(a.*view(nep.fD,1:k,:));
            @inbounds for j=1:p
                z .+= Av[j]*(view(VafD,:,j))
            end
            return z
        end
    end

    include("nep_type_helpers.jl");

end  # End Module


text:
using Dates

mutable struct Logger
    out::String
    levels::Dict
    log::Function
    function Logger(levels::Dict; out = pwd() * "logs/log.txt")
        log(level::Int64, message::String) = _log(level, message, levels, out)
        log(message::String) = _log(1, message, levels, out)
        log(http::HTTP.Stream, message::String) = _log(http, message)
        new(out, levels, log)
    end
    function Logger(; out = pwd() * "/logs/log.txt")
        if contains(out, "src")
            out = pwd() * "../logs/log.txt"
        end
        levels = Dict(1 => Crayon(foreground = :light_cyan),
        2 => Crayon(foreground = :light_yellow),
        3 => Crayon(foreground = :yellow, bold = true),
        4 => Crayon(foreground = :red, bold = true))
        Logger(levels; out = out)
    end
end
# print(, "In red. ", Crayon(bold = true), "Red and bold")
function _log(level::Int64, message::String, levels::Dict, out::String)
    time = now()
    if level > 1
        open(out, "w") do o
            seekend(o)
            try
                write(o, "[" * string(time) * "]: $message\n")
            catch
                try
                    touch(out)
                    write(o, "[" * string(time) * "]: $message"\n)
                catch
                    throw(ArgumentError("Cannot access logs."))
                end
            end
        end
    end
    println(Crayon(foreground = :light_gray, bold = true), "[", levels[level],
     string(time), Crayon(foreground = :light_gray, bold = true), "]: ",
     message)
end

function _log(http::HTTP.Stream, message::String)
    write(http, "<script>console.log('" * message * "');</script>")
end


text:
module Fields

export
    Face, Cell,
    AbstractField, AbstractLocatedField,
    Field, CellField, FaceFieldX, FaceFieldY, FaceFieldZ,
    interior, interiorparent,
    xnode, ynode, znode, location,
    set!,
    VelocityFields, TracerFields, tracernames, PressureFields, Tendencies

include("field_utils.jl")
include("field.jl")
include("set!.jl")
include("field_tuples.jl")
include("show_fields.jl")

end


text:
module GalsimBenchmark
using DataFrames
import FITSIO
import StaticArrays
import WCS

import Celeste: Model, DeterministicVI, ParallelRun, Infer
import Celeste.Model: CatalogEntry
import Celeste.ParallelRun: one_node_single_infer, one_node_joint_infer

const GALSIM_BENCHMARK_DIR = joinpath(Pkg.dir("Celeste"), "benchmark", "galsim")
const FILENAME = "output/galsim_test_images.fits"

function make_psf(psf_sigma_px)
    alphaBar = [1.; 0.]
    xiBar = [0.; 0.]
    tauBar = [psf_sigma_px^2 0.; 0. psf_sigma_px^2]
    [
        Model.PsfComponent(
            alphaBar[k],
            StaticArrays.SVector{2, Float64}(xiBar),
            StaticArrays.SMatrix{2, 2, Float64, 4}(tauBar)
        )
        for k in 1:2
    ]
end

immutable FitsExtension
    pixels::Matrix{Float32}
    header::FITSIO.FITSHeader
end

function read_fits(filename; read_sdss_psf=false)
    println("Reading $filename...")
    @assert isfile(filename)
    fits = FITSIO.FITS(filename)
    println("Found $(length(fits)) extensions.")

    extensions::Vector{FitsExtension} = []
    for extension in fits
        push!(extensions, FitsExtension(read(extension), FITSIO.read_header(extension)))
    end

    # assume WCS same for each extension
    wcs = WCS.from_header(FITSIO.read_header(fits[1], String))[1]

    close(fits)
    extensions, wcs
end

function make_images(band_pixels, psf, wcs, epsilon, iota)
    # assume dimensions equal for all images
    height, width = size(band_pixels[1])
    [
        Model.Image(
            height,
            width,
            band_pixels[band],
            band,
            wcs,
            psf,
            0, # SDSS run
            0, # SDSS camcol
            0, # SDSS field
            fill(epsilon, height, width),
            fill(iota, height),
            Model.RawPSF(Array(Float64, 0, 0), 0, 0, Array(Float64, 0, 0, 0)),
        )
        for band in 1:5
    ]
end

function typical_band_relative_intensities(is_star::Bool)
    source_type_index = is_star ? 1 : 2
    prior_parameters::Model.PriorParams = Model.load_prior()
    # Band relative intensities are a mixture of lognormals. Which mixture component has the most
    # weight?
    dominant_component = indmax(prior_parameters.k[:, source_type_index])
    # What are the most typical log relative intensities for that component?
    inter_band_ratios = exp(
        prior_parameters.c_mean[:, dominant_component, source_type_index]
        - diag(prior_parameters.c_cov[:, :, dominant_component, source_type_index])
    )
    Float64[
        1 / inter_band_ratios[2] / inter_band_ratios[1],
        1 / inter_band_ratios[2],
        1,
        inter_band_ratios[3],
        inter_band_ratios[3] * inter_band_ratios[4],
    ]
end

function typical_reference_brightness(is_star::Bool)
    source_type_index = is_star ? 1 : 2
    prior_parameters::Model.PriorParams = Model.load_prior()
    exp(
        prior_parameters.r_mean[source_type_index]
        - prior_parameters.r_var[source_type_index]
    )
end

# Since we're considering elliptical galaxy shapes, angle is only meaningful up to rotations of 180
# deg. This finds an equivalent angle in [0, 180) deg.
function canonical_angle(params)
    angle_radians = params[Model.ids.e_angle]
    while angle_radians < 0
        angle_radians += pi
    end
    while angle_radians > pi
        angle_radians -= pi
    end
    angle_radians
end

const BENCHMARK_PARAMETER_LABELS = String[
    "X center (world coords)",
    "Y center (world coords)",
    "Minor/major axis ratio",
    "Angle (degrees)",
    "Half-light radius (pixels)",
    "Brightness (nMgy)",
    "Color band 1-2 ratio",
    "Color band 2-3 ratio",
    "Color band 3-4 ratio",
    "Color band 4-5 ratio",
    "Probability of galaxy",
]

function get_field(header::FITSIO.FITSHeader, key::String)
    if haskey(header, key)
        header[key]
    else
        NA
    end
end

function actual_values(ids, star_galaxy_index, params)
    Float64[
        params[ids.u[1]],
        params[ids.u[2]],
        params[ids.e_axis],
        canonical_angle(params) * 180 / pi,
        params[ids.e_scale] * sqrt(params[ids.e_axis]),
        exp(params[ids.r1[star_galaxy_index]]),
        exp(params[ids.c1[1, star_galaxy_index]]),
        exp(params[ids.c1[2, star_galaxy_index]]),
        exp(params[ids.c1[3, star_galaxy_index]]),
        exp(params[ids.c1[4, star_galaxy_index]]),
        params[ids.a[2]],
    ]
end

function get_expected_dataframe(header)
    DataFrame(
        label=fill(header["CL_DESCR"], length(BENCHMARK_PARAMETER_LABELS)),
        field=BENCHMARK_PARAMETER_LABELS,
        expected=Any[
            get_field(header, "CL_X1"),
            get_field(header, "CL_Y1"),
            get_field(header, "CL_RTIO1"),
            get_field(header, "CL_ANGL1"),
            get_field(header, "CL_RADP1"),
            get_field(header, "CL_FLUX1"),
            get_field(header, "CL_C12_1"),
            get_field(header, "CL_C23_1"),
            get_field(header, "CL_C34_1"),
            get_field(header, "CL_C45_1"),
            header["CL_TYPE1"] == "star" ? 0 : 1,
        ])
end

function benchmark_comparison_data(single_infer_params, joint_infer_params, header)
    ids = Model.ids
    star_galaxy_index = header["CL_TYPE1"] == "star" ? 1 : 2
    comparison_dataframe = get_expected_dataframe(header)
    comparison_dataframe[:single_infer_actual] =
        actual_values(ids, star_galaxy_index, single_infer_params)
    comparison_dataframe[:joint_infer_actual] =
        actual_values(ids, star_galaxy_index, joint_infer_params)
    comparison_dataframe
end

function assert_counts_match_expected_flux(band_pixels::Vector{Matrix{Float32}},
                                           header::FITSIO.FITSHeader,
                                           iota::Float64)
    if !header["CL_NOISE"]
        expected_flux_nmgy = prod(size(band_pixels[3])) * header["CL_SKY"]
        for source_index in 1:header["CL_NSRC"]
            expected_flux_nmgy += header[string("CL_FLUX", source_index)]
        end
        expected_flux_counts = expected_flux_nmgy * iota
        @assert abs(sum(band_pixels[3]) - expected_flux_counts) / expected_flux_counts < 1e-3
    end
end

function make_catalog_entry(x_position_world_coords, y_position_world_coords)
    CatalogEntry(
        [x_position_world_coords, y_position_world_coords],
        false, # is_star
        # sample_star_fluxes
        typical_band_relative_intensities(true) .* typical_reference_brightness(true),
        # sample_galaxy_fluxes
        typical_band_relative_intensities(false) .* typical_reference_brightness(false),
        0.1, # gal_frac_dev
        0.7, # gal_ab
        pi / 4, # gal_angle
        4., # gal_scale
        "sample", # objid
        0, # thing_id
    )
end

function make_catalog(header::FITSIO.FITSHeader)
    catalog = CatalogEntry[]
    num_sources = header["CL_NSRC"]
    for source_index in 1:num_sources
        if num_sources == 1
            initial_position = [0.005335 for i in 1:2] # center of image
        else
            initial_position = [
                header[string("CL_X", source_index)],
                header[string("CL_Y", source_index)],
            ]
        end
        push!(catalog, make_catalog_entry(initial_position[1], initial_position[2]))
    end
    catalog
end

# Returns a data frame with one row for each test case and  parameter name, with columns
# * label: test case name
# * field: human-readable parameter name
# * expected: "ground truth" expected value
# * single_infer_actual: inferred MAP value or posterior median from single-source inference
# * joint_infer_actual: ditto, for multi-source joint inference
function main(; test_case_names=String[], print_fn=println,
              infer_source_callback=DeterministicVI.infer_source)
    all_benchmark_data = []
    extensions, wcs = read_fits(joinpath(GALSIM_BENCHMARK_DIR, FILENAME))
    @assert length(extensions) % 5 == 0 # one extension per band for each test case

    for test_case_index in 1:div(length(extensions), 5)
        first_band_index = (test_case_index - 1) * 5 + 1
        header = extensions[first_band_index].header
        this_test_case_name = header["CL_DESCR"]
        if !isempty(test_case_names) && !in(this_test_case_name, test_case_names)
            continue
        end
        println("Running test case '$this_test_case_name'")
        iota = header["CL_IOTA"]
        psf = make_psf(header["CL_SIGMA"])
        n_sources = header["CL_NSRC"]

        band_pixels = [
            extensions[index].pixels for index in first_band_index:(first_band_index + 4)
        ]
        assert_counts_match_expected_flux(band_pixels, header, iota)

        images = make_images(band_pixels, psf, wcs, header["CL_SKY"], iota)
        catalog = make_catalog(header)

        # we're only scoring one object per image
        target_sources = [1,]
        neighbor_map = Infer.find_neighbors(target_sources, catalog, images)

        target_sources_joint = collect(1:n_sources)
        neighbor_map_joint = Infer.find_neighbors(target_sources_joint, catalog, images)

        single_results = one_node_single_infer(
            catalog, target_sources, neighbor_map, images;
            infer_source_callback=infer_source_callback)
        joint_results = one_node_joint_infer(catalog, target_sources_joint,
                                             neighbor_map_joint, images)

        benchmark_data = benchmark_comparison_data(single_results[1].vs,
                                                   joint_results[1].vs,
                                                   header)
        print_fn(repr(benchmark_data))
        push!(all_benchmark_data, benchmark_data)
    end

    full_data = vcat(all_benchmark_data...)
    print_fn(repr(full_data[!isna(full_data[:expected]), :]))
    full_data
end

end # module GalsimBenchmark


text:
#using MAT
using WAV
#using Debug
#using Profile
using NumericExtensions
using HDF5

require("lock.jl")

#@iprofile begin
function ftest(data::Array{Float32,2}, tapers::Array{Float32,2}, p::Float64)

  (NC,C)=size(data)
  (NK,K)=size(tapers);
  N=NC

  Kodd = [1:2:K]
  Keven = [2:2:K]

  tapers2 = broadcast(*, tapers, ones(Float32,1,1,1))  # f x K x 1
  data = broadcast(*, data, ones(Float32,1,1,1))       # f x C x 1
  data = permutedims(data,(1,3,2))                     # f x 1 x C
  data = broadcast(*, data, tapers2)                   # f x K x C

  # use slow rfft for now b/c fast plan_rfft is buggy and in flux
  J = rfft(data,1)                                     # f x K x C

  Jp = J[:,Kodd,:]                               # f x K x C
  H0 = sum(tapers2[:,Kodd],1)                    # 1 x K
  H0sq = sumsq(H0)                               # 1
  JpH0 = squeeze(sum(broadcast(*, Jp, H0),2),2)  # f x C
  A = JpH0./H0sq                                 # f x C
  Kp = size(Jp,2)
  Ap = broadcast(*, A, ones(Float32,1,1,1))      # f x C x 1
  Ap = permutedims(Ap,(1,3,2))                   # f x 1 x C
  Jhat = broadcast(*, Ap, H0)                    # f x K x C
  num=(K-1)*abs2(A).*H0sq
  den=squeeze(sum(sqrdiff(Jp,Jhat),2),2)+squeeze(sum(abs2(J[:,Keven,:]),2),2)
  Fval=num./den

  return Fval

end
#end #profile

#@iprofile begin
function brown_puckette(x,k,fs)
  nfft = length(x)
  X=rfft(x,1)
  Xh0 = 0.5*(X[k]-0.5*X[k+1]-0.5*X[k-1])
  Xh1 = 0.5*exp(im*2*pi*(k-1)/nfft) *
      (X[k] - 0.5*exp(im*2*pi/nfft)*X[k+1] - 0.5*exp(-1im*2*pi/nfft)*X[k-1])
  phi0 = atan2(imag(Xh0),real(Xh0))
  phi1 = atan2(imag(Xh1),real(Xh1))
  if (phi1-phi0)<0
    phi1 = phi1+2*pi
  end
  freq = (phi1-phi0)*fs/(2*pi)

  period = fs/freq
  last = int(floor(period * floor(length(x)/period)))
  if last>0
    real_part = mean(x[1:last] .* cos([1:last].*(2*pi/period)))
    imag_part = mean(x[1:last] .* sin([1:last].*(2*pi/period)))
    amp = 2*abs(real_part + im*imag_part)
  else
    amp = NaN
  end
  return freq,amp
end
#end #profile

#@debug function do_it(params)
function do_it(params)
  FILEPATH = params[1]
  FILENAMES = params[2]
  FILENAME = params[3]
  FILETYPE = params[4]
  FILELEN_TIC = params[5]
  NCHANNELS = params[6]
  t = params[7]
  NWINDOWS_PER_WORKER = params[8]
  NFFT = params[9]
  PVAL = params[10]
  FS = params[11]
  tapers = params[12]
  sig = params[13]

  NSAMPLES = (NFFT>>1)*(NWINDOWS_PER_WORKER+1);
  if FILETYPE=="ch"
    for i = 1:NCHANNELS
      filei=string(FILEPATH,"/",FILENAMES[i])
      fid = open(filei,"r")
      seek(fid, t*4)
      tmp2 = read(fid, Float32, min(NSAMPLES, FILELEN_TIC-t))
      if(i==1)  dd=similar(tmp2, (size(tmp2,1), NCHANNELS))  end
      dd[:,i]=tmp2;
      close(fid)
    end
  else #if FILETYPE=="wav"
    dd, fs, nbits, extra =
        wavread("$FILENAME.wav"; subrange=(t+1):min(t+NSAMPLES, FILELEN_TIC))
    dd=float32(dd)
  end
  if (size(dd,1) < NSAMPLES)
    dd = [dd, zeros(Float32, (NSAMPLES-size(dd,1),NCHANNELS))]
  end

  idx={}
  for j = 1:NWINDOWS_PER_WORKER
    ddd=dd[[1:NFFT].+(j-1)*(NFFT>>1),:]
    F = ftest(ddd, tapers, PVAL)
#@bp
    for l = 1:NCHANNELS
      tmp = 1 .+ find(F[2:end-1,l] .> sig)
      for m = 1:length(tmp)
        freq, amp = brown_puckette(ddd[:,l],tmp[m],FS)
        push!(idx, [t/(NFFT>>1)+float64(j), freq, amp, float64(l)])
      end
    end
  end
  @printf("done with %.1f sec\n",t/FS);
  return [x[y] for x in idx, y=1:4]
end

max_stars_repo_name: JaneliaSciComp/Ax

max_issues_repo_name: JaneliaSciComp/Ax

max_forks_repo_name: JaneliaSciComp/Ax


text:
module HelpersTests

using Test

@testset "Macros" begin include("MacrosTests.jl") end

@testset "HelperFunctions" begin include("HelperFunctionsTests.jl") end

@testset "GridapTypes" begin include("GridapTypesTests.jl") end

end # module