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

using Markdown
using InteractiveUtils

# ╔═╡ 1f62d332-8303-11eb-31e8-4b34a42c4b9f
begin
  using Pkg
    Pkg.activate("MLJ_env", shared=true)

    using DataFrames
    using PlutoUI
    using Test
    using Printf
    using Plots
    using LinearAlgebra
    using Random
    using Statistics
	using StatsBase
	using Distributions
	using Images
	using ImageIO
end

# ╔═╡ f7ad6d8c-8467-11eb-081d-2ba27071eec5
using Colors #, ImageMetadata

# ╔═╡ d67867cc-8302-11eb-3d40-d5a50515cdb9
md"""
## C04w04: Fitting Gaussian Mixture Models with EM

In this assignment we will
  - implement the EM algorithm for a Gaussian mixture model
  - apply your implementation to cluster images
  - explore clustering results and interpret the output of the EM algorithm
"""

# ╔═╡ 53e17eba-8303-11eb-2471-eb8521bf88d2
md"""
### Implementing the EM algorithm for Gaussian mixture models

In this section, we will implement the EM algorithm taking the following steps:
  - Provide a log likelihood function for this model.
  - Implement the EM algorithm.
  - Create some synthetic data.
  - Visualize the progress of the parameters during the course of running EM.
  - Visualize the convergence of the model.
"""

# ╔═╡ 68c8f63a-8303-11eb-0ae3-c3f39e9acf8d
md"""
#### Log likelihood

We provide a function to calculate log likelihood for mixture of Gaussians. The log likelihood quantifies the probability of observing a given set of data under a particular setting of the parameters in our model.

We will use this to assess the convergence of our EM algorithm; specifically, we will keep looping through EM update steps until the log likehood ceases to increase at a certain rate.
"""

# ╔═╡ 84b9158c-8303-11eb-2ce4-c109ce3ff51a
function log_sum_exp(z)
    """
    Compute log(Σᵢ exp(zᵢ)) for some array z.
    """
    max_z = maximum(z)
    max_z + log(sum(exp.(z .- max_z)))
end

# ╔═╡ c01ba57e-8303-11eb-34e7-dfb911bac2f2
function log_likelihood(data, weights, μ, σ)
    """
    Compute the log likelihood of the data for a Gaussian mixture model
	with the given parameters.
    """
    nclusters, ndim = length(μ), size(data[1])[1]
	DT = eltype(data[1])
    llh = zero(DT)

    for d ∈ data
        z = zeros(DT, nclusters)
        for k ∈ 1:nclusters
            ## Compute: (x - μ)ᵀ ×  Σ⁻¹ × (x - μ)
            δ = reshape(d .- μ[k], :, 1) # make sure 2-dim!
            exp_term = dot(δ', inv(σ[k]) * δ)

            ## Compute log likelihood contribution for this data point
			## and this cluster
            z[k] += log(weights[k])
            z[k] -= 1/2. * (ndim * log(2. * π) + log(det(σ[k])) + exp_term)
    	end
		## Increment llh contribution of this data point across all clusters
		llh += log_sum_exp(z)
  	end
    llh
end

# ╔═╡ 835e2f8e-8444-11eb-3f3c-d534adaf4bd6
function approx(actᵥ::T, expᵥ::T; ϵ=1e-8) where T <: Real
	abs(actᵥ - expᵥ) ≤ ϵ
end

# ╔═╡ b0655cea-8445-11eb-1fbc-9350beaed253
function approx(actᵥ::AbstractArray{T}, expᵥ::AbstractArray{T};
		ϵ=1e-8) where T <: Real
	all(t -> approx(t[1], t[2]; ϵ), zip(actᵥ, expᵥ))
end

# ╔═╡ 2d4fb59a-8449-11eb-3f02-a5ece0576a9b
function approx(actᵥ::AbstractArray{T}, expᵥ::AbstractArray{T};
		ϵ=1e-8) where {T <: AbstractArray{T1} where T1 <: Real}
	all(t -> approx(t[1], t[2]; ϵ), zip(actᵥ, expᵥ))
end

# ╔═╡ 580675c0-830a-11eb-16e4-b963a391b138
md"""
#### E-step: assign cluster responsibilities, given current parameters

The first step in the EM algorithm is to compute cluster responsibilities.

Let $r_{ik}$ denote the responsibility of cluster $k$ for data point $i$.
Note that cluster responsibilities are fractional parts. *Cluster responsibilities for a single data point $i$ should sum to 1*.

$$r_{i1} + r_{i2} + \ldots + r_{iK} = 1$$

To figure how much a cluster is responsible for a given data point, we compute the likelihood of the data point under the  particular cluster assignment, multiplied by the weight of the cluster.

For data point $i$ and cluster $k$, this quantity is

$$r_{ik} \propto \pi_k N(x_i | \mu_k, \Sigma_k)$$

where $N(x_i | \mu_k, \Sigma_k)$ is the Gaussian distribution for cluster $k$ (with mean $\mu_k$ and covariance $\Sigma_k$).

We used $\propto$ because the quantity $N(x_i | \mu_k, \Sigma_k)$ is not yet the responsibility we want.
 To ensure that all responsibilities over each data point add up to 1, we add the normalization constant in the denominator:

$$r_{ik} = \frac{\pi_k N(x_i | \mu_k, \Sigma_k)}{\sum_{k=1}^{K} \pi_k N(x_i | \mu_k, \Sigma_k)}$$

Complete the following function that computes $r_{ik}$ for all data points $i$ and clusters $k$.
"""

# ╔═╡ 4d69bcec-830c-11eb-0c80-314c8cded6db
function compute_responsibilities(data, weights, μ, σ)
    """
    E-step fo EM Algorithm:
    - compute responsibilities, given the current parameters
    """
	@assert length(data) > 0

    ndata, nclusters = length(data), length(μ)
    resp = zeros(eltype(data[1]), ndata, nclusters)

    ## Update resp. matrix so that resp[i, k] is the responsibility of
	## cluster k for data point i.
    for i ∈ 1:ndata
        for k ∈ 1:nclusters
			d = MvNormal(μ[k], σ[k])
			resp[i, k] = weights[k] * pdf(d, data[i])
		end
	end
    ## Add up responsibilities over each data point and normalize
    row_sums = reshape(sum(resp, dims=2), :, 1)
    resp ./= row_sums
end

# ╔═╡ 31835154-830d-11eb-09fd-196046609868
begin
  # const ϵ = 1e-8
  resp₀ = compute_responsibilities([[1., 2.], [-1., -2.]],
	[0.3, 0.7],                         # weights
    [[0., 0.], [1., 1.]],               # μ
	[[1.5 0.; 0. 2.5], [1. 1.; 1. 2.]]  # σ
  )

  @test size(resp₀) == (2, 2)
  @test approx(resp₀, [0.10512733 0.89487267; 0.46468164 0.53531836])
end

# ╔═╡ 430e58a8-8310-11eb-2a53-ed2608f43122
md"""
#### M-step: Update parameters, given current cluster responsibilities

Once the cluster responsibilities are computed, we update the parameters (*weights*, *means*, and *covariances*) associated with the clusters.

**Computing soft counts**. Before updating the parameters, we first compute what is known as "soft counts".

The soft count of a cluster is the sum of all cluster responsibilities for that cluster:

$$N^{\text{soft}}_k = r_{1k} + r_{2k} + \ldots + r_{Nk} = \sum_{i=1}^{N} r_{ik}$$

where we loop over data points. Note that, unlike k-means, we must loop over every single data point in the dataset. This is because all clusters are represented in all data points, to a varying degree.

"""

# ╔═╡ 42ef055c-8310-11eb-1ae5-d5ae8bf22c1f
"""
  Compute the total responsibility assigned to each cluster, which will be useful when
  implementing M-steps below. In the lectures this is called N^soft
"""
compute_soft_counts(resp) = sum(resp; dims=1)

# ╔═╡ 93988ed8-8310-11eb-15f4-554b537eacd6
md"""
**Updating weights.** The cluster weights show us how much each cluster is represented over all data points.

The weight of cluster $k$ is given by the ratio of the soft count
$N^{\text{soft}}_{k}$ to the total number of data points $N$:

$$\hat{\pi}_k = \frac{N^{\text{soft}}_{k}}{N}$$

Notice that $N$ is equal to the sum over the soft counts $N^{\text{soft}}_{k}$ of all clusters.

Complete the following function:
"""

# ╔═╡ 937da226-8310-11eb-1966-6bcc8ff3cb47
function compute_weights(counts)
    nclusters, n = length(counts), sum(counts)
    weights = zeros(eltype(counts), nclusters)

	for k ∈ 1:nclusters
        ## Update the weight for cluster k using the M-step update rule for
		## the cluster weight: π̂ₖ
        weights[k] = counts[k] / n
	end
    weights
end

# ╔═╡ 93653664-8310-11eb-23f2-d5236b8bc0ed
begin
	resp₁ = compute_responsibilities([[1., 2.], [-1.,-2.], [0., 0.]],
		[0.3, 0.7],                         # weights
		[[0., 0.], [1., 1.]],               # μ
		[[1.5  0.; 0. 2.5], [1. 1.; 1. 2.]] # σ
  	)
	counts₁ = compute_soft_counts(resp₁)
	weights₁ = compute_weights(counts₁)

	@test approx(weights₁, [0.27904865942515705 0.720951340574843])
end

# ╔═╡ 93490854-8310-11eb-3de5-5b24918d8840
with_terminal() do
  println(" - counts: $(counts₁)")
  println(" - weigths: $(weights₁)")
  println("Checkpoint passed!")
end

# ╔═╡ 42d1c078-8310-11eb-2262-bb5ae92a2b27
md"""
**Updating means**. The mean of each cluster is set to the [weighted average](https://en.wikipedia.org/wiki/Weighted_arithmetic_mean) of all data points, weighted by the cluster responsibilities:

$$\hat{\mu}_k = \frac{1}{N_k^{\text{soft}}} \sum_{i=1}^N r_{ik}x_i$$

Complete the following function:
"""

# ╔═╡ 1b27764e-8315-11eb-3854-f9e5b591744a
function compute_means(data, resp, counts)
    nclusters, ndata = length(counts), length(data)
    means = [zeros(eltype(data[1]), length(data[1])) for _ ∈ 1:nclusters]

    for k in 1:nclusters
        ## Update means for cluster k using the M-step update rule
		## for the mean variables.
        ## This will assign the variable means[k] to be our estimate for μ̂ₖ.
        weighted_sum = sum(
          [resp[i, k] * data[i] for i ∈ 1:ndata]
        )
        means[k] = weighted_sum / counts[k]
	end
    means
end

# ╔═╡ 34b1bb5e-8317-11eb-0260-5f1d3132a7ba
begin
	data_tmp = [[1., 2.], [-1.,-2.]]
	resp₂ = compute_responsibilities(data_tmp,
		[0.3, 0.7],                         # weights
		[[0., 0.], [1., 1.]],               # μ
		[[1.5  0.; 0. 2.5], [1. 1.; 1. 2.]] # σ
	)
	counts₂ = compute_soft_counts(resp₂)
	means₂ = compute_means(data_tmp, resp₂, counts₂)
	exp_means₂ = [[-0.6310085, -1.262017], [0.25140299, 0.50280599]]
	@test approx(means₂[1], exp_means₂[1])
	@test approx(means₂[2], exp_means₂[2])
end

# ╔═╡ 1ad0cbe4-8315-11eb-341c-dbce64acb374
with_terminal() do
  println(" - counts: $(counts₂)")
  println(" - means: $(means₂)")
  println("Checkpoint passed!")
end

# ╔═╡ 1a99af8a-8315-11eb-0f74-eb97741a40ca
md"""

**Updating covariances**.  The covariance of each cluster is set to the weighted average of all [outer products](https://people.duke.edu/~ccc14/sta-663/LinearAlgebraReview.html), weighted by the cluster responsibilities:

$$\hat{\Sigma}_k = \frac{1}{N^{\text{soft}}_k}\sum_{i=1}^N r_{ik} (x_i - \hat{\mu}_k)(x_i - \hat{\mu}_k)^T$$

The "outer product" in this context refers to the matrix product

$$(x_i - \hat{\mu}_k)(x_i - \hat{\mu}_k)^T.$$
Letting $(x_i - \hat{\mu}_k)$ to be $d \times 1$ column vector, this product is a $d \times d$ matrix. Taking the weighted average of all outer products gives us the covariance matrix, which is also $d \times d$.

Reminder:
The inner product (of two column vectors $1 \times n$) denoted as $<v, w> = v^{t}w$ is a scalar, whereas the outer product of the same two column vectors denoted as $v \otimes w = vw^{t}$ is a matrix.

Complete the following function:

"""

# ╔═╡ 986e362c-8317-11eb-0eef-cd1de11cd4b2
function compute_covariances(data, resp, counts, μ)
    nclusters, ndim, ndata = length(counts), size(data[1])[1], length(data)
    covs = [zeros(eltype(data[1]), ndim, ndim) for _ ∈ 1:nclusters]

    for k ∈ 1:nclusters
        ## Update covariances for cluster k using the M-step update rule
		## for covariance variables.
        ## This will assign the variable covariances[k] to be the estimate for Σ̂ₖ.
        w_sum = zeros(eltype(data[1]), ndim, ndim)
        for i ∈ 1:ndata
            w_sum += resp[i, k] * ((data[i] - μ[k]) * (data[i] - μ[k])')
		end
		covs[k] = (1. / counts[k]) * w_sum
	end
    covs
end

# ╔═╡ 57c76562-8373-11eb-0e43-256ec189f399
begin
	## deps on prev test cell for resp₂, counts₂, means₂
	covs₂ = compute_covariances(data_tmp, resp₂, counts₂, means₂)

	exp_covs₂ = [
		[0.60182827 1.20365655; 1.20365655 2.4073131],
		[0.93679654 1.87359307; 1.87359307 3.74718614]
	]

	@test approx(covs₂[1], exp_covs₂[1])
	@test approx(covs₂[2], exp_covs₂[2])
end

# ╔═╡ a12bc08e-8371-11eb-29f0-e14f2c17c1d1
with_terminal() do
  println(" - covariances: $(covs₂)")
  println("Checkpoint passed!")
end

# ╔═╡ c59d7e92-8373-11eb-09d1-55814d7e9f9e
md"""
#### The EM algorithm

Now let us write a function that takes initial parameter estimates and runs EM.
"""

# ╔═╡ c5808488-8373-11eb-3d9a-51c955524390
function em(data, init_μ, init_σ, init_weights;
    maxiter=1000, ϵ=1e-4)
    ## Make copies of initial parameters, which we will update during each iteration
    μ, σ = copy(init_μ), copy(init_σ)
    weights = copy(init_weights)

    ## Infer dimensions of dataset and the number of clusters
    ndata, ndim, nclusters = length(data), size(data[1])[1], length(μ)

    ## Initialize some useful variables
    resp = zeros(eltype(data[1]), ndata, nclusters)
    llh = log_likelihood(data, weights, μ, σ)
    llh_trace = [llh]
	n_iter = 1
    for it ∈ 1:maxiter
		n_iter += 1
        it % 5 == 0 && println("Iteration $(it) / $(maxiter)")

        ## E-step: compute responsibilities
        resp = compute_responsibilities(data, weights, μ, σ)

        ## M-step
        ## Calc. the total resp. assigned to each cluster, useful when
        ## implementing M-steps below. In the lectures this is called N^{soft}
        counts = compute_soft_counts(resp)

        ## update weight of cluster k using M-step update rule for cluster weight π̂ₖ
        weights = compute_weights(counts)

        ## update means of cluster k using the M-step update rule for the mean vars.
        ## assign the variable μ[k] to be our estimate for μ̂ₖ.
        μ = compute_means(data, resp, counts)

        ## update covs ofcluster k using the M-step update rule for cov vars.
        ## assign the variable σ[k] to be the estimate for Σ̂ₖ.
        σ = compute_covariances(data, resp, counts, μ)

        ## Compute the loglikelihood at this iteration
        llh_latest = log_likelihood(data, weights, μ, σ)
        push!(llh_trace, llh_latest)

        ## Check for convergence in log-likelihood and store
        (llh_latest - llh) < ϵ && llh_latest > -Base.Inf && break
        llh = llh_latest
  	end
  	Dict{Symbol, Any}(
    	:weights => weights,
    	:means => μ,
    	:covs => σ,
    	:loglik => llh_trace,
    	:resp => resp
  	)
end

# ╔═╡ c56791bc-8373-11eb-00e3-230f82206892
md"""
#### Testing the implementation on the simulated data

To help us develop and test our implementation, we will generate some observations from a mixture of Gaussians and then run our EM algorithm to discover the mixture components.
We'll begin with a function to generate the data, and a quick plot to visualize its output for a 2-dimensional mixture of three Gaussians.

Now we will create a function to generate data from a mixture of Gaussians model.
"""

# ╔═╡ a790272e-8375-11eb-2a40-df02f3d1ca24
function generate_mog_data(ndata, μ, σ, weights; seed=42)
    """
    Creates a list of data points
    """
    nclusters = length(weights)
    data = []

	Random.seed!(seed)
	rng = MersenneTwister(seed)

    for i ∈ 1:ndata
        ## Use sample (from StatsBase) and weights to pick a:
		## 1 ≤ cluster_id < nclusters
		k = sample(rng, 1:nclusters, Weights(weights), 1; replace=false)[1]

		## Use np.random.multivariate_normal to create data from this cluster
		x = rand(MvNormal(μ[k], σ[k]))  ## ????

		# np.random.multivariate_normal(means[k], covariances[k])
        push!(data, x)
	end
	data
end

# ╔═╡ 00261792-8441-11eb-3969-55f797a05e23

# ╔═╡ a740cde6-8375-11eb-0f2a-efc1f23f8f38
begin
  init_μ = [
	[5., 0.], # mean of cluster 1
	[1., 1.], # mean of cluster 2
	[0., 5.]  # mean of cluster 3
  ]

  init_σ =[
	[.5 0.; 0 .5],       # covariance of cluster 1
	[.92 .38; .38 .91],  # covariance of cluster 2
	[.5 0.; 0 .5]]       # covariance of cluster 3

  init_weights = [1/4., 1/2., 1/4.]  # weights of each cluster

  data = generate_mog_data(100, init_μ, init_σ, init_weights)
end

# ╔═╡ a72ae382-8375-11eb-15fc-c750a071dccf
scatter(map(t -> t[1], data), map(t -> t[2], data), color=[:lightblue], lw=3)

# ╔═╡ b3e522d6-844c-11eb-2aa5-750b3ee4c212
md"""
**Checkpoint** for *em* function
"""

# ╔═╡ 9aee7716-83d7-11eb-3f3c-d534adaf4bd6
 begin
 	m_data = [[0.31845368, 5.03182099], [1.45029339, 2.15717573],
 [1.02645431, 4.76042375], [ 5.23493625, -0.81138851], [ 4.65710743, -0.08532351],
 [0.43746553, 4.93778385], [-0.12756373,  5.14296293], [4.18200738, 0.24819247],
 [ 1.65646672, -0.87991369], [-0.69508006,  5.03848976], [6.01348014, 1.87722666],
 [ 0.20099208, -0.93213203], [-0.38268222,  0.26365374], [ 5.03394554, -1.66740986],
 [ 1.53889119, -0.12731715], [0.81448617, 1.63543709], [1.7219243 , 0.28460715],
 [4.90757786, 0.56716088], [0.60512066, 4.87555125], [ 1.30867865, -0.58711045],
 [0.19570654, 5.73836612], [3.97078181, 0.39328234], [0.47190777, 1.47637309],
 [1.32157017, 0.56486368], [3.89988459, 0.92433204], [0.11205488, 6.4072012 ],
 [0.40717315, 1.67134485], [0.56231532, 5.01349864],
 [1.82113515, 0.78892223], [-0.77578232,  4.62650898], [ 0.45032469, -1.01099583],
 [2.25024431, 0.64832921], [1.47157789, 1.68387047], [-0.0420993 ,  0.50923158],
 [1.18348882, 0.56978228], [5.52201954, 0.46015202], [0.69012115, 1.40085221],
 [0.02126064, 4.76228496], [5.32547107, 1.46856369], [2.58605452, 2.11944448],
 [0.182161  , 4.44577711], [ 0.11836846, -0.33361633], [0.66083532, 1.03879141],
 [1.91659334, 0.63739381], [0.28487706, 4.5665366 ], [0.4764496 , 5.72581327],
 [3.10453287, 2.19997423], [-0.85376105,  4.17389402], [0.08542074, 4.38781614],
 [0.21418408, 5.60131432], [-0.29932464,  4.62574683], [2.25835904, 1.1068076 ],
 [ 4.77358263, -0.06002466], [0.73093444, 1.40220082], [0.39250487, 0.57133562],
 [0.8949198, 1.9300788], [-0.25767001,  4.11851619], [0.55042767, 0.21762904],
 [-0.34493741,  4.57098776], [4.71632515, 0.63293438], [1.67820052, 1.29111918],
 [ 1.4862169 , -0.07691204], [5.46625127, 0.01184258], [-1.33846048, -0.20047555],
 [0.01458455, 5.5685916 ], [1.49605116, 4.82940436], [ 5.75963176, -0.13066323],
 [1.70380591, 1.31376129], [ 4.76328545, -0.33344549], [ 0.8945073 , -1.13036243],
 [1.61279831, 0.66792909], [-0.20028537,  5.0349852 ], [-0.79147695,  4.72403614],
 [1.33337808, 1.427949  ], [1.81545643, 2.37052768], [1.02894019, 3.21262401],
 [-0.1914201 ,  4.03020414], [1.45589226, 0.85329139], [2.23603443, 0.87710097],
 [0.32608726, 0.09240739], [-0.18786611,  5.85099597], [4.29931854, 1.35428737],
 [0.49849313, 0.68444377], [-0.8336478 ,  5.81745508], [-0.57795091, -1.6674263 ],
 [1.83292993, 0.7963827 ],[1.01736806, 0.8455941 ], [0.48493064, 0.45675926],
 [-0.47669162,  3.91542167],[0.96494312, 0.65663175], [1.04999151, 0.38350349],
 [1.17458921, 0.73066238],[0.04638356, 4.86196943], [1.66742479, 1.54181187],
 [0.94372389, 1.32117047],[ 5.63522868, -0.0138849 ],[0.52699679, 5.37301739],
 [4.80426637, 0.8017483], [-0.00321541,  0.2083691],[1.38949165, 0.95020545]
	];

 	m_μs = [[0.19570654, 5.73836612], [6.01348014, 1.87722666],
		[0.52699679, 5.37301739]];

	m_σs =[[ 3.41361848  -1.91940532; -1.91940532  4.68785076],
		[ 3.41361848 -1.91940532; -1.91940532  4.68785076],
		[ 3.41361848 -1.91940532; -1.91940532  4.68785076]];

	m_weights = [0.3333333333333333, 0.3333333333333333, 0.3333333333333333];

	m_results = em(m_data, m_μs, m_σs, m_weights);
	:done
end

# ╔═╡ 4b50e4b0-8444-11eb-3969-55f797a05e23
begin
	@test approx(m_results[:means], Array{Float64,1}[
			Float64[0.021382851931258, 4.947729003491],
			Float64[4.942392352282, 0.3136531084758],
			Float64[1.0818112539095, 0.73903507897]])

	@test approx(m_results[:weights], Float64[0.3007102300634,
			0.17993710074063, 0.5193526691959])

	@test approx(m_results[:covs], Array{Float64,2}[
			[0.2932613984949 0.05048454589297; 0.05048454589297 0.352815373349],
			[0.3556437022527 -0.014948748777286; -0.014948748777286 0.6669502500144],
			[0.6711499192647 0.33058964530721; 0.33058964530721 0.9042972427651]])

	@test approx(m_results[:weights], [0.3007102300634,
			0.17993710074063, 0.5193526691959])

	@test approx(m_results[:loglik][1:5], [-541.3161249736, -372.1355927523,
			-366.9935696838, -365.6599199088, -364.3348647749, -362.87960462340,
			-361.3750403053])
end

# ╔═╡ a7135b54-8375-11eb-2d63-392c4fa94aea
md"""
**Now we will fit a mixture of Gaussians** to this data using our implementation of the EM algorithm. As with k-means, it is important to ask how we obtain an initial configuration of mixing weights and component parameters.

In this simple case, we'll take three random points to be the initial cluster means, use the empirical covariance of the data to be the initial covariance in each cluster (a clear overestimate), and set the initial mixing weights to be uniform across clusters.
"""

# ╔═╡ a6fc6796-8375-11eb-02af-01cee06a06dd
begin
	seed₂ = 2010
	Random.seed!(seed₂)
	rng = MersenneTwister(seed₂)

	## Initialization of parameters
	const K = 3
	chosen = sample(rng, 1:length(data), K; replace=false)
  	#        np.random.choice(length(data), K, replace=False)
  	@assert length(chosen) == K

  	init_μ₁ = [data[x] for x in chosen];
  	init_σ₁ = [cov(data) for _ ∈ 1:K]; # cov(data, rowvar=0) => col. is an observ.
  	init_w = [1.0 / K for _ ∈ 1:K];

  	## Run EM
  	results = em(data, init_μ₁, init_σ₁, init_w);
	:done
end

# ╔═╡ a6b33828-8375-11eb-027a-130ef7a14128
md"""
**Note**. Like k-means, EM is prone to converging to a local optimum. In practice, you may want to run EM multiple times with different random initialization. We have omitted multiple restarts to keep the assignment reasonably short.

Our algorithm returns a dictionary with five elements:
  - 'loglik': a record of the log likelihood at each iteration
  - 'resp': the final responsibility matrix
  - 'means': a list of K means
  - 'covs': a list of K covariance matrices
  - 'weights': the weights corresponding to each model component

"""

# ╔═╡ 4db52910-844d-11eb-3969-55f797a05e23
results[:means]

# ╔═╡ 6f4476bc-844d-11eb-1fbc-9350beaed253
results[:weights]

# ╔═╡ 5ab06896-844d-11eb-3f3c-d534adaf4bd6
md"""
*For the following quiz questions, please round your answer to the nearest thousandth (3 decimals)*

**Quiz Question**: What is the weight that EM assigns to the first component after running the above codeblock?
"""

# ╔═╡ a69ff718-8375-11eb-3965-89ecd2088d97
round(results[:weights][1]; digits=3)

# ╔═╡ a687f546-8375-11eb-3435-1fa822f71ff1
md"""
**Quiz Question**: Using the same set of results, obtain the mean that EM assigns the second component. What is the mean in the first dimension?
"""

# ╔═╡ a66ede62-8375-11eb-3746-758cb7efc8a8
## Answer:   the second component / first dim
round(results[:means][2][1], digits=3)

# ╔═╡ a657ed74-8375-11eb-1d97-1bd417fa42ee
md"""
**Quiz Question**: Using the same set of results, obtain the covariance that EM assigns the third component. What is the variance in the first dimension?
"""

# ╔═╡ a63e00e4-8375-11eb-10db-897f94a0dc32
# Answer:   covs third component / first dim
results[:covs][3, 1]

# ╔═╡ c54c0456-8373-11eb-1cb6-11bf225bc6b0
md"""
### Plot progress of parameters

One useful feature of testing our implementation on low-dimensional simulated data is that we can easily visualize the results.

We will use the following `plot_contours` function to visualize the Gaussian components over the data at three different points in the algorithm's execution:

  1. At initialization (using initial_mu, initial_cov, and initial_weights)
  2. After running the algorithm to completion
  3. After just 12 iterations (using parameters estimates returned when setting `maxiter=12`)
"""

# ╔═╡ 9ba36afc-83d7-11eb-32d4-4d4d0df107d4
function bivariate_normal(X, Y;
    σ_x=1.0, σ_y=1.0, μ_x=0.0, μ_y=0.0, σ_xy=0.0)
    Xμ = X .- μ_x
    Yμ = Y .- μ_y

    ρ = σ_xy / (σ_x * σ_y)
    z = Xμ.^2 ./ σ_x^2 + Yμ.^2 ./ σ_y^2 - 2 * ρ * Xμ .* Yμ ./ (σ_x * σ_y)

  	den = 2. * π * σ_x * σ_y * √(1 - ρ^2)
	exp.(-z / (2. * (1. - ρ^2))) / den
end

# ╔═╡ cd5eb61e-83d9-11eb-1df3-9389091527e8
function meshgrid(x, y)
  @assert length(x) == length(y)
  n = length(x)

  vx, vy = copy(x), copy(y)
  for _ ∈ 1:n-1
    vx = hcat(vx, x)
    vy = hcat(vy, y)
  end
  vx, vy = reshape(vx, :, n)', reshape(vy, :, n)'
  Matrix(vx), Matrix(vy)
end

# ╔═╡ 5a708ed0-83d8-11eb-147c-abf3f80f8c92
function plot_contours(data, means, covs, title)
    δ = 0.025
    k = length(means)
    X = Y = range(-3.0; stop=9.0, step=δ)

	col = [:green, :red, :indigo]
	parms = []
	scatter([_x[1] for _x ∈ data], [_y[2] for _y ∈ data], legend=false, title=title)

	for i ∈ 1:K
		μ, σ = means[i], covs[i]
		σ_x, σ_y = √(σ[1, 1]), √(σ[2, 2])
		σ_xy = σ[1, 2] / (σ_x * σ_y)

		push!(parms, (σ_x=σ_x, σ_y=σ_y, μ_x=μ[1], μ_y=μ[2], σ_xy=σ_xy))
	end

	contour!(X, Y, (x, y) -> bivariate_normal(x, y; parms[1]...), color=col[1])
  	contour!(X, Y, (x, y) -> bivariate_normal(x, y; parms[2]...), color=col[2])
  	contour!(X, Y, (x, y) -> bivariate_normal(x, y; parms[3]...), color=col[3])
end

# ╔═╡ 9b4a0a5e-83d7-11eb-3f02-a5ece0576a9b
plot_contours(data, init_μ₁, init_σ₁, "Initial clusters")

# ╔═╡ e9ae52fe-845e-11eb-15c1-65ebc1f4df8e
begin
	# Parameters after running EM to convergence
	plot_contours(data, results[:means], results[:covs], "Final clusters")
end

# ╔═╡ 3ab0e342-845f-11eb-321d-61972f4e6bd4
md"""
Fill in the following code block to visualize the set of parameters we get after running EM for 12 iterations.
"""

# ╔═╡ 4321d874-845f-11eb-107c-c32e7c640947
begin
	res12 = em(data, init_μ₁, init_σ₁, init_w);
	plot_contours(data, res12[:means], res12[:covs], "Clusters after 12 iterations")
end

# ╔═╡ ffd6d3a2-845f-11eb-1308-39662d5cf06e
md"""
**Quiz Question**: Plot the loglikelihood that is observed at each iteration. Is the loglikelihood plot monotonically increasing, monotonically decreasing, or neither [multiple choice]?
"""

# ╔═╡ a3ac0552-845f-11eb-2283-6defb781f7c7
plot(1:length(results[:loglik]), results[:loglik], lw=3,
	xlabel="Iteration",
	ylabel="Log-likelihood",
	legend=false)

# ╔═╡ 125d76fa-8460-11eb-1d49-57e4eea41e89
md"""
 - The loglikelihood is monotonically increasing
"""

# ╔═╡ c72c03ba-8469-11eb-2fc0-47dab53b1ce9
md"""
## Fitting a Gaussian mixture model for image data¶
"""

# ╔═╡ 075f9b50-8464-11eb-35b3-c945b1e49ccf
begin
	const IMG_DIR = "../../ML_UW_Spec/data/images"
	images = []

	for dir_entry ∈ cd(readdir, IMG_DIR)
		println("Found $(dir_entry)")
		sdir = join([IMG_DIR, dir_entry], "/")

		for img_file ∈ readdir(sdir)
			push!(images, load(join([sdir, img_file], "/")))
		end
	end
end

# ╔═╡ a27f740e-8465-11eb-171f-3ff85b79e29b
length(images)

# ╔═╡ 77bd32d4-8466-11eb-2713-6b11add92ba3
mosaicview(images[1:12]...; nrow=3, ncol=4)

# ╔═╡ fb990ba8-8466-11eb-13be-7d3f33c0576c
mosaicview(images[400:411]...; nrow=3, ncol=4)

# ╔═╡ 36448e26-8467-11eb-12eb-23daf58ea752
mosaicview(images[700:711]...; nrow=3, ncol=4)

# ╔═╡ d1b3f548-8466-11eb-304e-119b7c99959e
mosaicview(images[1000:1011]...; nrow=3, ncol=4)

# ╔═╡ 8c8995d0-8467-11eb-3efb-d1a7961fdefa
img_data = [
	[mean(red.(img)), mean(green.(img)), mean(blue.(img))] for img ∈ images
]

# ╔═╡ bcc9335c-8469-11eb-1492-d310a071b09a
md"""
We need to come up with initial estimates for the mixture weights and component parameters. Let's take three images to be our initial cluster centers, and let's initialize the covariance matrix of each cluster to be diagonal with each element equal to the sample variance from the full data. As in our test on simulated data, we'll start by assuming each mixture component has equal weight.

This may take a few minutes to run.
"""

# ╔═╡ 7b6cc3d8-846e-11eb-315c-5544ae999cc8
var_col(img_data, i::Integer) = var(x[i] for x ∈ img_data)

# ╔═╡ daffa734-8469-11eb-35ea-376de2c6e1e3
begin
	seed₃ = 2010
	Random.seed!(seed₃)
	rng₃ = MersenneTwister(seed₃)

	# Initalize parameters
	const KC = 4
	const NC = 3  # 3 colors R, G, B
	n_images = length(images)

	im_init_μs = [img_data[x] for x ∈ sample(rng₃, 1:n_images, KC; replace=false)]
	@assert length(im_init_μs) == KC

	id_mat = Matrix{Float64}(I, NC, NC);
	im_cov = id_mat .* vec([var_col(img_data, c) for c in 1:NC]);

	im_init_σs = [im_cov for _ ∈ 1:KC];     # [cov, cov, cov, cov]
	im_init_ws = [1.0 / KC for _ ∈ 1:KC]; # [1/4., 1/4., 1/4., 1/4.]

	## Run our EM algorithm on the image data using the above initializations.
	## This should converge in about 125 iterations
	res_img = em(img_data, im_init_μs, im_init_σs, im_init_ws);
	:done
end

# ╔═╡ 7fb68656-8471-11eb-3969-55f797a05e23
function _plot(data_x, data_y)
  plot(data_x, data_y, lw=4, xlabel="Iteration", ylabel="Log-likelihood'")
end

# ╔═╡ ab34bece-8471-11eb-3f3c-d534adaf4bd6
begin
	llh_ = res_img[:loglik]
	_plot(1:length(llh_), llh_)
end

# ╔═╡ cccb8bd0-8471-11eb-1fbc-9350beaed253
md"""
The log likelihood increases so quickly on the first few iterations that we can barely see the plotted line. Let's plot the log likelihood after the first three iterations to get a clearer view of what's going on:
"""

# ╔═╡ d63fd518-8471-11eb-147c-abf3f80f8c92
_plot(3:length(llh_), llh_[3:end])

# ╔═╡ d624eb68-8471-11eb-32d4-4d4d0df107d4
md"""
#### Evaluating uncertainty

Next we'll explore the evolution of cluster assignment and uncertainty. Remember that the EM algorithm represents uncertainty about the cluster assignment of each data point through the responsibility matrix. Rather than making a 'hard' assignment of each data point to a single cluster, the algorithm computes the responsibility of each cluster for each data point, where the responsibility corresponds to our certainty that the observation came from that cluster.

We can track the evolution of the responsibilities across iterations to see how these 'soft' cluster assignments change as the algorithm fits the Gaussian mixture model to the data; one good way to do this is to plot the data and color each point according to its cluster responsibilities. Our data are three-dimensional, which can make visualization difficult, so to make things easier we will plot the data using only two dimensions, taking just the [R G], [G B] or [R B] values instead of the full [R G B] measurement for each observation.

"""

# ╔═╡ 37bd479e-8472-11eb-3106-e3d88fb54de4
# TODO...

# ╔═╡ 37a19df0-8472-11eb-2662-6157d71dd17f
md"""
#### Interpreting each cluster

Let's dig into the clusters obtained from our EM implementation. Recall that our goal in this section is to cluster images based on their RGB values. We can evaluate the quality of our clustering by taking a look at a few images that 'belong' to each cluster. We hope to find that the clusters discovered by our EM algorithm correspond to different image categories - in this case, we know that our images came from four categories ('cloudy sky', 'rivers', 'sunsets', and 'trees and forests'), so we would expect to find that each component of our fitted mixture model roughly corresponds to one of these categories.

If we want to examine some example images from each cluster, we first need to consider how we can determine cluster assignments of the images from our algorithm output. This was easy with k-means - every data point had a 'hard' assignment to a single cluster, and all we had to do was find the cluster center closest to the data point of interest. Here, our clusters are described by probability distributions (specifically, Gaussians) rather than single points, and our model maintains some uncertainty about the cluster assignment of each observation.

One way to phrase the question of cluster assignment for mixture models is as follows: how do we calculate the distance of a point from a distribution? Note that simple Euclidean distance might not be appropriate since (non-scaled) Euclidean distance doesn't take direction into account.  For example, if a Gaussian mixture component is very stretched in one direction but narrow in another, then a data point one unit away along the 'stretched' dimension has much higher probability (and so would be thought of as closer) than a data point one unit away along the 'narrow' dimension.

In fact, the correct distance metric to use in this case is known as [Mahalanobis distance](https://en.wikipedia.org/wiki/Mahalanobis_distance). For a Gaussian distribution, this distance is proportional to the square root of the negative log likelihood. This makes sense intuitively - reducing the Mahalanobis distance of an observation from a cluster is equivalent to increasing that observation's probability according to the Gaussian that is used to represent the cluster. This also means that we can find the cluster assignment of an observation by taking the Gaussian component for which that observation scores highest. We'll use this fact to find the top examples that are 'closest' to each cluster.

__Quiz Question:__ Calculate the likelihood (score) of the first image in our data set (`images[0]`) under each Gaussian component through a call to `multivariate_normal.pdf`.  Given these values, what cluster assignment should we make for this image? Hint: don't forget to use the cluster weights.
"""

# ╔═╡ 3787dc3a-8472-11eb-129f-b5e34bda1328
res_img[:weights]

# ╔═╡ 37702d1a-8472-11eb-1df3-9389091527e8
res_img[:means]

# ╔═╡ d606647c-8471-11eb-3f02-a5ece0576a9b
# d = MvNormal(μ[k], σ[k])
#			resp[i, k] = weights[k] * pdf(d, data[i])

# ╔═╡ d92f92e4-8472-11eb-17cf-eb0a1c444ee9
begin
	proba_img = Float64[]
	for k in 1:KC
		d = MvNormal(res_img[:means][k], res_img[:covs][k])
		push!(proba_img, res_img[:weights][k] * pdf(d, img_data[1]))
	end
	proba_img
end

# ╔═╡ d9132afa-8472-11eb-38c5-39d44a1e57f7
md"""
Now we calculate cluster assignments for the entire image dataset using the result of running EM for 60 iterations above.
"""

# ╔═╡ d8f3fd56-8472-11eb-3856-99f30a89c6d3
begin
	p_weights = res_img[:weights]
	p_means = res_img[:means]
	p_covs = res_img[:covs]
	rgb = img_data

	_N = length(img_data)   ## number of images
	_K = length(p_means)    ## number of clusters

	assignments = zeros(Float64, _N)
	probs = zeros(Float64, _N)

	for i ∈ 1:_N
    	## Compute the score of data point i under each Gaussian component:
    	p = [
			p_weights[k] * pdf(MvNormal(p_means[k], p_covs[k]), rgb[i]) for k ∈ 1:_K
    	]

    	## Compute assign. of each data point to a given cluster based on the scores:
    	assignments[i] = argmax(p)
    	## For data point i, store the corresponding score under this cluster assign:
    	probs[i] = maximum(p)
	end

	df_assign = DataFrame(:cluster_num => Int.(assignments), :probs => probs, :image => images);
	:done
end

# ╔═╡ 9a3f1a20-8475-11eb-1fbc-9350beaed253
last(df_assign, 5)

# ╔═╡ 9a2134e2-8475-11eb-3f3c-d534adaf4bd6
md"""
We'll use the 'assignment' DataFrame to find the top images from each cluster by sorting the datapoints within each cluster by their score under that cluster (stored in probs).

Create a function that returns the top 5 images assigned to a given category in our data.
"""

# ╔═╡ 9a00824c-8475-11eb-3969-55f797a05e23
function get_top_images(df, cluster; k=5)
    images_in_cluster = df[df.cluster_num .== cluster, :]
	top_images = sort(images_in_cluster, [:probs], rev=true)[1:k, :]
    top_images.image
end

# ╔═╡ 17ac70be-8476-11eb-3f02-a5ece0576a9b
mosaicview(get_top_images(df_assign, 1)...; nrow=2, ncol=3)

# ╔═╡ 8b254648-8477-11eb-1df3-9389091527e8
mosaicview(get_top_images(df_assign, 2)...; nrow=2, ncol=3)

# ╔═╡ 8b068190-8477-11eb-147c-abf3f80f8c92
mosaicview(get_top_images(df_assign, 3)...; nrow=2, ncol=3)

# ╔═╡ 8ae806ca-8477-11eb-32d4-4d4d0df107d4
mosaicview(get_top_images(df_assign, 4)...; nrow=2, ncol=3)

# ╔═╡ a5dc1b88-8477-11eb-129f-b5e34bda1328
md"""
These look pretty good! Our algorithm seems to have done a good job overall at 'discovering' the four categories that from which our image data was drawn. It seems to have had the most difficulty in distinguishing between rivers and cloudy skies, probably due to the similar color profiles of images in these categories; if we wanted to achieve better performance on distinguishing between these categories, we might need a richer representation of our data than simply the average [R G B] values for each image.
"""

# ╔═╡ Cell order:
# ╟─d67867cc-8302-11eb-3d40-d5a50515cdb9
# ╠═1f62d332-8303-11eb-31e8-4b34a42c4b9f
# ╟─53e17eba-8303-11eb-2471-eb8521bf88d2
# ╟─68c8f63a-8303-11eb-0ae3-c3f39e9acf8d
# ╠═84b9158c-8303-11eb-2ce4-c109ce3ff51a
# ╠═c01ba57e-8303-11eb-34e7-dfb911bac2f2
# ╠═835e2f8e-8444-11eb-3f3c-d534adaf4bd6
# ╠═b0655cea-8445-11eb-1fbc-9350beaed253
# ╠═2d4fb59a-8449-11eb-3f02-a5ece0576a9b
# ╟─580675c0-830a-11eb-16e4-b963a391b138
# ╠═4d69bcec-830c-11eb-0c80-314c8cded6db
# ╠═31835154-830d-11eb-09fd-196046609868
# ╟─430e58a8-8310-11eb-2a53-ed2608f43122
# ╟─42ef055c-8310-11eb-1ae5-d5ae8bf22c1f
# ╟─93988ed8-8310-11eb-15f4-554b537eacd6
# ╠═937da226-8310-11eb-1966-6bcc8ff3cb47
# ╠═93653664-8310-11eb-23f2-d5236b8bc0ed
# ╠═93490854-8310-11eb-3de5-5b24918d8840
# ╟─42d1c078-8310-11eb-2262-bb5ae92a2b27
# ╠═1b27764e-8315-11eb-3854-f9e5b591744a
# ╠═34b1bb5e-8317-11eb-0260-5f1d3132a7ba
# ╠═1ad0cbe4-8315-11eb-341c-dbce64acb374
# ╟─1a99af8a-8315-11eb-0f74-eb97741a40ca
# ╠═986e362c-8317-11eb-0eef-cd1de11cd4b2
# ╠═57c76562-8373-11eb-0e43-256ec189f399
# ╠═a12bc08e-8371-11eb-29f0-e14f2c17c1d1
# ╟─c59d7e92-8373-11eb-09d1-55814d7e9f9e
# ╠═c5808488-8373-11eb-3d9a-51c955524390
# ╟─c56791bc-8373-11eb-00e3-230f82206892
# ╠═a790272e-8375-11eb-2a40-df02f3d1ca24
# ╠═00261792-8441-11eb-3969-55f797a05e23
# ╠═a740cde6-8375-11eb-0f2a-efc1f23f8f38
# ╠═a72ae382-8375-11eb-15fc-c750a071dccf
# ╟─b3e522d6-844c-11eb-2aa5-750b3ee4c212
# ╠═9aee7716-83d7-11eb-3f3c-d534adaf4bd6
# ╠═4b50e4b0-8444-11eb-3969-55f797a05e23
# ╟─a7135b54-8375-11eb-2d63-392c4fa94aea
# ╠═a6fc6796-8375-11eb-02af-01cee06a06dd
# ╟─a6b33828-8375-11eb-027a-130ef7a14128
# ╠═4db52910-844d-11eb-3969-55f797a05e23
# ╠═6f4476bc-844d-11eb-1fbc-9350beaed253
# ╟─5ab06896-844d-11eb-3f3c-d534adaf4bd6
# ╠═a69ff718-8375-11eb-3965-89ecd2088d97
# ╟─a687f546-8375-11eb-3435-1fa822f71ff1
# ╠═a66ede62-8375-11eb-3746-758cb7efc8a8
# ╟─a657ed74-8375-11eb-1d97-1bd417fa42ee
# ╠═a63e00e4-8375-11eb-10db-897f94a0dc32
# ╟─c54c0456-8373-11eb-1cb6-11bf225bc6b0
# ╠═9ba36afc-83d7-11eb-32d4-4d4d0df107d4
# ╠═cd5eb61e-83d9-11eb-1df3-9389091527e8
# ╠═5a708ed0-83d8-11eb-147c-abf3f80f8c92
# ╠═9b4a0a5e-83d7-11eb-3f02-a5ece0576a9b
# ╠═e9ae52fe-845e-11eb-15c1-65ebc1f4df8e
# ╟─3ab0e342-845f-11eb-321d-61972f4e6bd4
# ╠═4321d874-845f-11eb-107c-c32e7c640947
# ╟─ffd6d3a2-845f-11eb-1308-39662d5cf06e
# ╠═a3ac0552-845f-11eb-2283-6defb781f7c7
# ╟─125d76fa-8460-11eb-1d49-57e4eea41e89
# ╟─c72c03ba-8469-11eb-2fc0-47dab53b1ce9
# ╠═075f9b50-8464-11eb-35b3-c945b1e49ccf
# ╠═a27f740e-8465-11eb-171f-3ff85b79e29b
# ╠═77bd32d4-8466-11eb-2713-6b11add92ba3
# ╠═fb990ba8-8466-11eb-13be-7d3f33c0576c
# ╠═36448e26-8467-11eb-12eb-23daf58ea752
# ╠═d1b3f548-8466-11eb-304e-119b7c99959e
# ╠═f7ad6d8c-8467-11eb-081d-2ba27071eec5
# ╠═8c8995d0-8467-11eb-3efb-d1a7961fdefa
# ╟─bcc9335c-8469-11eb-1492-d310a071b09a
# ╠═7b6cc3d8-846e-11eb-315c-5544ae999cc8
# ╠═daffa734-8469-11eb-35ea-376de2c6e1e3
# ╠═7fb68656-8471-11eb-3969-55f797a05e23
# ╠═ab34bece-8471-11eb-3f3c-d534adaf4bd6
# ╟─cccb8bd0-8471-11eb-1fbc-9350beaed253
# ╠═d63fd518-8471-11eb-147c-abf3f80f8c92
# ╟─d624eb68-8471-11eb-32d4-4d4d0df107d4
# ╠═37bd479e-8472-11eb-3106-e3d88fb54de4
# ╟─37a19df0-8472-11eb-2662-6157d71dd17f
# ╠═3787dc3a-8472-11eb-129f-b5e34bda1328
# ╠═37702d1a-8472-11eb-1df3-9389091527e8
# ╠═d606647c-8471-11eb-3f02-a5ece0576a9b
# ╠═d92f92e4-8472-11eb-17cf-eb0a1c444ee9
# ╟─d9132afa-8472-11eb-38c5-39d44a1e57f7
# ╠═d8f3fd56-8472-11eb-3856-99f30a89c6d3
# ╠═9a3f1a20-8475-11eb-1fbc-9350beaed253
# ╟─9a2134e2-8475-11eb-3f3c-d534adaf4bd6
# ╠═9a00824c-8475-11eb-3969-55f797a05e23
# ╠═17ac70be-8476-11eb-3f02-a5ece0576a9b
# ╠═8b254648-8477-11eb-1df3-9389091527e8
# ╠═8b068190-8477-11eb-147c-abf3f80f8c92
# ╠═8ae806ca-8477-11eb-32d4-4d4d0df107d4
# ╟─a5dc1b88-8477-11eb-129f-b5e34bda1328


text:
module TensorPolynomialBasesTests

using TensorPolynomialBases
using Test

@testset "UtilsTests" begin
  include("UtilsTests.jl")
end

@testset "InterfacesTests" begin
  include("InterfacesTests.jl")
end

@testset "MonomialBasesTests" begin
  include("MonomialBasesTests.jl")
end

@testset "QGradMonomialBasesTests" begin
  include("QGradMonomialBasesTests.jl")
end

@testset "QCurlGradMonomialBasesTests" begin
  include("QCurlGradMonomialBasesTests.jl")
end

end # module

ext: jl


text:
meanstd(x) = (μ = mean(x); Ref(μ .± std(x, mean = μ)))

function descriptive_stats(df)
    d = stack(df, [:turning_point, :center_of_search], variable_name = :point_type, value_name = :point_xy)
    gd = groupby(d, [:point_type, :group, :nest_coverage])
    g = copy(combine(gd, :point_xy => meanstd => :point_xy, nrow))
    g.id = repeat(1:nrow(g)÷2, outer = 2)
    d = unstack(g, [:group, :nest_coverage, :nrow], :point_type, :point_xy)
    sort!(d, [:nest_coverage, :group])
    select!(d, [1,2,4,5,3])
    open("table1.txt", "w") do io
        pretty_table(io, d, ["Group" "Burrow coverage" "Turning point" "Gravity center" "n";
                               "" "" "μ ± σ" "μ ± σ" ""],
                     hlines = [1,7],
                     alignment = [:l, :l, :c, :c, :r],
                     formatters = (v,i,j) -> 3 ≤ j ≤ 4  ? myformat(v) : v
                    )
    end
end

myformat(_::Missing) = "-"
myformat(x::Measurement) = string(round(Int, x.val), "±", round(Int, x.err))
myformat(xs::AbstractVector{String}) = string("(", xs[1], ",", xs[2], ")")
myformat(xs::AbstractVector{Float64}) = string("(", round(xs[1], digits = 2), ",", round(xs[2], digits = 2), ")")
myformat(xs::AbstractVector{Measurement{T}}) where {T <: Real} = myformat(myformat.(xs))

########### speed

speed(d) = (norm(p2.xy - p1.xy)/(p2.t - p1.t) for (p1, p2) in zip(d, lag(d, -1, default = d[end])) if p1.xy ≠ p2.xy)
function foo(d)
    length(d) < 10 && return missing
    s = speed(d)
    μ = mean(s)
    μ ± std(s, mean = μ)
end
function speeds!(df)
    df.homing_speed = [foo(track.rawcoords[1:track.tp]) for track in df.track]
    df.search_speed = [foo(track.rawcoords[track.tp:end]) for track in df.track]
    v = combine(groupby(df, :group), :homing_speed => mean ∘ skipmissing, :search_speed => mean ∘ skipmissing)
    println("speeds:")
    println(v)
    μ = mean(skipmissing(vcat(df.homing_speed, df.search_speed)))
    @printf "mean speed is %s cm/s\n" μ
end

################### displaced closed nest stats

roundsig(x) = x ≤ 0.05 ? @sprintf("%.2g", x) : "ns"

HypothesisTests.pvalue(t::StatsModels.TableRegressionModel) = GLM.coeftable(t).cols[4][2]

function displaced_stats(df)
    gdf = groupby(df, [:group, :nest_coverage])
    g = gdf[[(group = group, nest_coverage = "closed") for group in ("none", "right","left","towards", "away")]]

    n = 10^6
    tbls = map((:turning_point, :center_of_search)) do point
        data = combine(g, point => (x -> x .- Ref(mean(x))) => :centered)
        x1 = [r.centered for r in eachrow(data) if r.group == "none"]
        x2 = [r.centered for r in eachrow(data) if r.group ≠ "none"]
        tx = ApproximatePermutationTest(first.(x1), first.(x2), var, n)
        ty = ApproximatePermutationTest(last.(x1), last.(x2), var, n)
        @printf "%i permutations sampled at random from %i possible permutations\n" n factorial(big(length(x1)+length(x2)))
        σ = pvalue.([tx, ty], tail = :left)
        nσ = length(x1) + length(x2)
        @printf "Is the %s variance of the none group, (%i, %i) significantly smaller than the variancve of the displaced groups, (%i, %i)? P = (%.2g, %.2g) (n = %i)\n" replace(string(point), '_' => ' ') var(x1)... var(x2)... σ... nσ
        data = combine(g, point => (x -> first.(x)) => :x, point => (x -> last.(x)) => :y, :group => (x -> first.(intended.(x))) => :ix, :group => (x -> last.(intended.(x))) => :iy)
        tx = lm(@formula(ix ~ x), data)
        ty = lm(@formula(iy ~ y), data)
        nμ = nrow(data)
        @printf "Is the effect (%.2f, %.2f) of the displacement on the %s significant? P = (%.2g, %.2g) (n = %i)\n" coef(tx)[2] coef(ty)[2] replace(string(point), '_' => ' ') pvalue(tx) pvalue(ty) nμ
        combine(g, [:group, point] => ((g, x) -> myformat(roundsig.(pvalue.([
                                                                             OneSampleTTest(first.(x), first(intended(g[1]))),
                                                                             OneSampleTTest(last.(x),   last(intended(g[1])))
                                                                            ])))) => point, nrow => :n)
    end
    tbl = hcat(tbls[1][!, Not(All(:n, :nest_coverage))], tbls[2][!, Not(All(:group, :nest_coverage))])
    m = Matrix(tbl)
    writedlm("table2.csv", m, ',')
end

############### transfer far stats

onesample(xy, g) = myformat([pvalue(OneSampleTTest(f.(xy), f(intended(g[1])))) for f in (first, last)])

function transfer_stats(df)
    gdf = groupby(df, :group)
    groups = ("far", "back")
    g = gdf[[(group = g,) for g in groups]]
    tbls = map((:turning_point, :center_of_search)) do point
        combine(g, [point, :group] => onesample => point, nrow)
    end
    tbl = hcat(tbls[1][:, Not(:nrow)], tbls[2][:, Not(:group)])
    m = Matrix(tbl)
    writedlm("table3.csv", m, ',')
end


text:
include("qwsearch/qwsearch.jl")
include("qwevolution/qwevolution.jl")


text:
module SmithNormalForm

using LinearAlgebra
using SparseArrays
using Base.CoreLogging

import Base: show, summary
import LinearAlgebra: diagm, diag

export snf, smith, Smith

include("bezout.jl")
include("snf.jl")

# ---------------------------------------------------------------------------------------- #

struct Smith{P,Q<:AbstractMatrix{P},V<:AbstractVector{P}} <: Factorization{P}
    S::Q
    Sinv::Q
    T::Q
    Tinv::Q
    SNF::V
    Smith{P,Q,V}(S::AbstractMatrix{P}, Sinv::AbstractMatrix{P},
                 T::AbstractMatrix{P}, Tinv::AbstractMatrix{P},
                 D::AbstractVector{P}) where {P,Q,V} = new(S, Sinv, T, Tinv, D)
end
Smith(S::AbstractMatrix{P}, T::AbstractMatrix{P}, SNF::AbstractVector{P}) where {P} =
    Smith{P,typeof(S),typeof(SNF)}(S, similar(S, 0, 0), T, similar(T, 0, 0), SNF)

# ---------------------------------------------------------------------------------------- #

"""
    smith(X::AbstractMatrix{P}; inverse::Bool=true) --> Smith{P,Q,V}

Return a Smith normal form of an integer matrix `X` as a `Smith` structure (of element type
`P`, matrix type `Q`, and invariant factor type `V`).

The Smith normal form is well-defined for any matrix ``m×n`` matrix `X` with elements in a
principal domain (PID; e.g., integers) and provides a decomposition of `X` into ``m×m``,
``m×n``, and `S`, `Λ`, and `T` as `X = SΛT`, where `Λ` is a diagonal matrix with entries
("invariant factors") Λᵢ ≥ Λᵢ₊₁ ≥ 0 with nonzero entries divisible in the sense Λᵢ | Λᵢ₊₁.
The invariant factors can be obtained from [`diag(::Smith)`](@ref).

`S` and `T` are invertible matrices; if the keyword argument `inverse` is true (default),
the inverse matrices are computed and returned as part of the `Smith` factorization.
"""
function smith(X::AbstractMatrix{P}; inverse::Bool=true) where {P}
    S, T, D, Sinv, Tinv = snf(X, inverse=inverse)
    SNF = diag(D)
    return Smith{P, typeof(X), typeof(SNF)}(S, Sinv, T, Tinv, SNF)
end

# ---------------------------------------------------------------------------------------- #

"""
    diagm(F::Smith) --> AbstractMatrix

Return the Smith normal form diagonal matrix `Λ` from a factorization `F`.
"""
function diagm(F::Smith{P}) where P
    rows = size(F.S, 1)
    cols = size(F.T, 1)
    D    = issparse(F.SNF) ? spzeros(P, rows, cols) : zeros(P, rows, cols)
    for (i,Λᵢ) in enumerate(diag(F))
        D[i,i] = Λᵢ
    end
    return D
end

"""
    diag(F::Smith) --> AbstractVector

Return the invariant factors (or, equivalently, the elementary divisors) of a Smith normal
form `F`.
"""
diag(F::Smith) = F.SNF

summary(io::IO, F::Smith{P,Q,V}) where {P,Q,V} = print(io, "Smith{", P, ",", Q, ",", V, "}")
function show(io::IO, ::MIME"text/plain", F::Smith)
    summary(io, F)
    println(io, " with:")
    Base.print_matrix(io, diagm(F), " Λ = "); println(io)
    Base.print_matrix(io, F.S, " S = "); println(io)
    Base.print_matrix(io, F.T, " T = ")
    return nothing
end

end # module


text:
using FilterLogging
using Test

@testset "FilterLogging.jl" begin
    # Write your tests here.
end

max_stars_repo_name: benelsen/FilterLogging

max_issues_repo_name: benelsen/FilterLogging

max_forks_repo_name: benelsen/FilterLogging


text:
#!/usr/bin/env julia

# Make sure we have all the CUDA gubbins set up.
using Pkg
Pkg.activate(".")
Pkg.instantiate()
using CUDA

# Actual CUDA kernel.
function pichunk!(threads, n, result)
  rank=threadIdx().x
  lower=floor((rank-1) * n/threads) + 1
  upper=floor((rank) * n/threads)
  step = 1.0/n
  for i=lower:upper
    x = (i + 0.5) * step
    result[rank] += 4.0/(1.0 + x*x)
  end
  return
end

# Sort out CUDA thread setup.
function _picalc(nthr, slices)
  result = CUDA.fill(0.0, nthr)
  @cuda threads=nthr pichunk!(nthr, slices, result)
  synchronize()
  s = sum(result)/slices
  return(s)
end

# Run the benchmark.
function picalc(numsteps, nthr)

  println("Calculating PI using:")
  println("  ", numsteps, " slices")
  println("  ", nthr, " CUDA threads(s)")
  start = time()

  mypi = _picalc(nthr, numsteps)

  elapsed = time() - start

  println("Obtained value of PI: ", mypi)
  println("Time taken: ", round(elapsed, digits=3), " seconds")

end

# Parse arguments.
numsteps = if length(ARGS) > 0
  parse(Int, ARGS[1])
else
  1_000_000_000
end

threads = if length(ARGS) > 1
  parse(Int, ARGS[2])
else
  64
end

# Warm things up
print("  Warming up...")
warms = time()
dummy = _picalc(1,10)
warmt = time() - warms
println("done. [", round(warmt, digits=3), "s]\n")

# Run the full example
picalc(numsteps,threads)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule grep_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("grep")
JLLWrappers.@generate_main_file("grep", UUID("47013893-54eb-5c9c-83b4-9a24496ece36"))
end  # module grep_jll


text:
# How-to guides
push!(LOAD_PATH,"../src/")

ENV["GKSwstype"]="100"

using Documenter
using Plots
using Literate
using WaterWaves1D

# generate example

DOC_OUTPUT  = joinpath(@__DIR__, "src", "examples")
NB_OUTPUT   = joinpath(@__DIR__, "src", "notebooks")

 #Literate.markdown(joinpath(@__DIR__, "..",  "examples","ManyModelsforWaterWaves.jl"), DOC_OUTPUT, documenter=true)
 #Literate.notebook(joinpath(@__DIR__, "..",  "notebooks","Example.ipynb"), NB_OUTPUT, execute=false)
 #Literate.markdown(joinpath(@__DIR__, "..",  "notebooks","FullDispersion.ipynb"), DOC_OUTPUT, documenter=true)
 #Literate.notebook(joinpath(@__DIR__, "..",  "notebooks","FullDispersion.ipynb"), NB_OUTPUT, execute=false)

makedocs(modules=[WaterWaves1D],
         doctest = false,
         authors = "Vincent Duchene",
         format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
         sitename = "WaterWaves1D.jl",
         #linkcheck = true,
         #checkdocs = :all,
         pages = [
            "Home" => "index.md",
            "Background" => "background.md",
            "Problems" => "problems.md",
            "How-to..." => "how-to.md",
            "Example" => "example.md",
            "Plot recipes" => "plot_recipes.md",
            "Library" => "library.md"
                   ])

#                  "Animation"     => "examples/animation.md",
#                 "Example"       => "examples/two_problems.md",

deploydocs(
    repo   = "github.com/WaterWavesModels/WaterWaves1D.jl.git"
 )


text:
using ArgParse,CairoMakie,Distributions

CairoMakie.activate!(type="pdf")

function parse_argument()
end

function main()
    graph_name ="exponential_dist_mass"

    xs::Array{Float64} = 0:0.1:5.0

    graph_params = [Dict("params"=>Dict("λ"=>1.0), "linestyle"=>nothing, "label"=>"Ex(1.0)"),
                    Dict("params"=>Dict("λ"=>1.5), "linestyle"=>:dash,   "label"=>"Ex(1.5)"),
                    Dict("params"=>Dict("λ"=>0.5), "linestyle"=>:dot,    "label"=>"Ex(0.5)")]

    fig = Figure(resolution = (800, 600))

    axis = Axis(fig[1,1], title = "exponential distribution mass function")

    for params in graph_params
        λ = params["params"]["λ"]
        ys = [pdf(Exponential(λ), x) for x in xs]
        lines!(fig[1,1], xs, ys, linestyle=params["linestyle"], label=params["label"])
    end

    axislegend(axis)

    save("$graph_name.pdf", fig)
end

main()

max_stars_repo_name: seisyuu-hantatsushi/julia

max_issues_repo_name: seisyuu-hantatsushi/julia

max_forks_repo_name: seisyuu-hantatsushi/julia


text:
#Backtracking Sequential Coloring algorithm
include("custom_graph.jl")

function bsc_color(G::CGraph)
    V = num_vertices(G)
    F = zeros(Int64, V)
    freeColors = [Vector{Int64}() for _ in 1:V] #set of free colors for each vertex
    U = zeros(Int64, 0) #stores set of free colors
    #F = zeros(Int64, 0) #stores final coloring of vertices

    sorted_vertices = order_by_degree(G)
    start = 1
    optColorNumber = V + 1
    x = sorted_vertices[1]
    colors[0] = 0
    push!(U, 1)
    freeColors[x] = U

    while (start >= 1)
        back = false
        for i = 1:V
            if i > start
                x = find_uncolored_vertex(sorted_vertices, F)
                U = free_colors(x,F,G,optColorNumber)
                sort(U)
            end
            if length(U) > 0
                k = U[1]
                F[x] = k
                deleteat!(U,1)
                freeColors[x] = copy(U)
                l = colors[i-1]
                colors[i] = max(k,l)
            else
                start = i - 1
                back = true
                break
            end
        end
        if back
            if start >= 1
                x = A[start]
                F[x] = 0 #uncolor x
                U = freeColors[x]
            end
        else
            F_opt = F
            optColorNumber = colors[V-1]
            i = least_index(sorted_vertices,optColorNumber,G)
            start = i - 1
            if start < 1
                break #leave the while loop
            end
            uncolor_all(F, sorted_vertices, start, G)
            for i = 0:start
                x  = A[i]
                U = freeColors[x]
                U = remove_colors(U, optColorNumber)
                freeColors[x] = copy(U)
            end
        end
    end
    return F_opt
end

"""
    vertex_degree(G,z)

Find the degree of the vertex z which belongs to the graph G.
"""
function degree(G::CGraph,z::Int64)
    return length(neighbors(G,z))
end

function sorted_vertices(G::CGraph)
    V = length(vertices(G))
    marked = zeros(Int64,V)
    sv = zeros(Int64,0)
    max_degree = -1
    max_degree_vertex = -1
    for i = 1:V
        max_degree = -1
        max_degree_vertex = -1
        for j = 1:V
            if j != i
                if degree(G,j) > max_degree && marked[j] == 0
                    max_degree = degree(G,j)
                    max_degree_vertex = j
                end
            end
        end
        push!(sv,max_degree_vertex)
        marked[max_degree_vertex] = 1
    end
    return sv
end

#find uncolored vertex of maximal degree of saturation
function find_uncolored_vertex(sv::Array{Int64,1}, G::CGraph)
    colors = zeros(Int64,0)
    max_colors = -1
    max_color_index = -1
    for i = 1:length(vertices(G))
        if F[i] != 0
            for j in neighbors(G,i)
                if F[j] != 0 && F[j] in colors == false
                    push!(colors, F[j])
                end
            end
            if length(colors) > max_colors
                max_colors = length(colors)
                max_color_index = i
            end
        end
        colors = zeros(Int64,0)
    end
    for i = 1:length(vertices(G))
        if A[i] == max_color_index
            return i
        end
    end

end

#set of free colors of x, which are < optColorNumber
function free_colors(x::Int64, F::Array{Int64,1}, G::CGraph, max_color::Int64)
    colors = zeros(Int64,0)
    for color in 1:max_color
        present = true
        for y in neighbors(G,x)
            if F[y] == color
                present = false
                break
            end
        end
        if present
            push!(colors,color)
        end
    end
    return colors
end

#least index with F(A[i]) = optColorNumber
function least_index(A::Array{Int64, 1}, F::Array{Int64,1}, optColorNumber::Int64, G::CGraph)
    for i = 1:length(G.vertices)
        if F[A[i]] == optColorNumber
            return i
        end
    end
end

#uncolor all vertices A[i] with i >= start
function uncolor_all(F::Array{Int64,1}, A::Array{Int64,1}, start::Int64, G::CGraph)
    for i = start:length(G.vertices)
        F[A[i]] = 0
    end
end

#remove from U all colors >= optColorNumber
function remove_colors(U::Array{Int64,1}, optColorNumber::Int64)
    modified_U = zeros(Int64,0)
    for i = 1:length(U)
        if U[i] < optColorNumber
            push!(mmodified_U, U[i])
        end
    end
    return modified_U
end


text:
# Autogenerated wrapper script for libLAS_jll for i686-linux-gnu-cxx03
export liblas, liblas_c

using boost_jll
using GDAL_jll
using PROJ_jll
using libgeotiff_jll
JLLWrappers.@generate_wrapper_header("libLAS")
JLLWrappers.@declare_library_product(liblas, "liblas.so.3")
JLLWrappers.@declare_library_product(liblas_c, "liblas_c.so.3")
function __init__()
    JLLWrappers.@generate_init_header(boost_jll, GDAL_jll, PROJ_jll, libgeotiff_jll)
    JLLWrappers.@init_library_product(
        liblas,
        "lib/liblas.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        liblas_c,
        "lib/liblas_c.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    BlockWeighter(side)

A weighting method that assigns weights to points in spatial object
based on blocks of given `side`. The number (n) of points inside a
block determines the weights (1/n) of these points.
"""
struct BlockWeighter{T} <: AbstractWeighter
  side::T
end

function weight(object::AbstractSpatialObject, weighter::BlockWeighter)
  p = partition(object, BlockPartitioner(weighter.side))

  weights = Vector{Float64}(undef, npoints(object))
  for s in subsets(p)
    n = length(s)
    weights[s] .= 1/n
  end

  SpatialWeights(domain(object), weights)
end


text:
function getparam(d...;init=Knet.xavier)
	et = ctx.dtype
	atype = ctx.device == :gpu ? Knet.KnetArray{et} : Array{et}
	Knet.Param(atype(init(d...)))
end

const _layernames = Dict{String,Int}()

function get_layername(layername; kwargs...)
	info = Dict(kwargs)
	if haskey(info, :name)
		info[:name]
	else
		if ! haskey(_layernames, layername)
			_layernames[layername] = 1
		end
		id = _layernames[layername]
		_layernames[layername] += 1
		"$(layername)_$(id)"
	end
end

abstract type LazyLayer <: Layer end

## Generic code, X can be a Tensor or a Tuple
function (layer::Layer)(X)
	@debug "Calling $(typeof(layer))" layer
    call(layer, X)
end

## Generic code, X can be a Tensor or a Tuple
function (layer::LazyLayer)(X)
	@assert hasproperty(layer, :built) "LazyLayer $layer without built property detected"
	@assert isa(layer.built, Bool) "LazyLayer $layer built property not a boolean"

    if ! layer.built
		@debug "Building $(typeof(layer))" layer
        build(layer, size(X)[1:end-1])
		layer.built = true
    end
	@debug "Calling $(typeof(layer))" layer
    call(layer, X)
end

"""
Regular densely-connected NN layer.

Dense implements the operation: output = activation(dot(input, weight) + bias) where
activation is the element-wise activation function passed as the activation argument,
weight is a weights matrix created by the layer, and bias is a bias vector created
by the layer (only applicable if use_bias is true).

# Usage

```julia
layer = Dense(10, relu)
layer = Dense(100, use_bias=false)
```

"""
mutable struct Dense <:LazyLayer
    units::Int
	activation::Function
	use_bias::Bool
	name::String
	built::Bool
	init::NamedTuple
	params::NamedTuple

	function Dense(units::Int, activation=identity; use_bias=true,
		initw = Knet.xavier, initb = zeros, kwargs...)
		@assert units > 0 "Units of a Dense layer should be > 0"
		name = get_layername("dense"; kwargs...)
	    new(units, activation, use_bias, name, false,
		(w=initw, b=initb), (w=nothing, b=nothing))
	end
end

function call(layer::Dense, X::Tensor)
	X = Knet.mat(X) # Flatten if required

	w,b = layer.params
	if layer.use_bias
    	layer.activation.(w*X .+ b)
	else
		layer.activation.(w*X)
	end
end

function build(layer::Dense, shape::Shape)
	# nInput = length(shape) > 1 ? *(shape...) : shape[1]
	nInput = *(shape...)

    w = getparam(layer.units, nInput, init=layer.init.w)
	b = nothing
    if layer.use_bias
		b = getparam(layer.units, init=layer.init.b)
	end
	layer.params = (w=w,b=b)
end

"""
Flattening Layer. Photon by default already has flattening funcitonality
build into the Dense layer, so you won't need to include a separate Flatten
layer before a Dense layer.
"""
mutable struct Flatten <: Layer
	dims
	Flatten(dims=nothing) = new(dims)
end

function call(layer::Flatten, X::Tensor)
	layer.dims === nothing ? Knet.mat(X) : Knet.mat(X, dims=layer.dims)
end

## Activation
mutable struct Activation <:Layer
	activation
end

function call(c::Activation, X::Tensor)
	c.activation.(X)
end

"""
Dropout layer with optional the rate (between 0 and 1) of dropout. If
no rate is specified, 0.5 (so 50%) will be used.
"""
struct Dropout <: Layer
	rate::Float64

	function Dropout(rate=0.5)
		@assert rate <= 1.0
		@assert rate >= 0.0
		new(rate)
	end
end

function call(d::Dropout, X::Tensor)
	Knet.dropout(X, d.rate)
end

"""
Batch normalization layer (Ioffe and Szegedy, 2014). Normalizes the input at each batch,
i.e. applies a transformation that maintains the mean activation close to 0
and the activation standard deviation close to 1.

Finally, if activation is not nothing, it is applied to the outputs as well.
"""
struct BatchNorm <: Layer

  	moments
	activation

	function BatchNorm(activation=identity)
		new(Knet.bnmoments(), activation)
	end
end

function call(bn::BatchNorm, X::Tensor)
	bn.activation.(Knet.batchnorm(X, bn.moments))
end

"""
Beginning of allowing for a single model instance to run on multiple devices
(expiremental)
"""
struct ContextSwitch <: Layer

  	device::Symbol
	deviceId::Int
	dtype::Type

	function ContextSwitch(;device=ctx.device,deviceId=ctx.deviceId,dtype=ctx.dtype)
		new(device,deviceId,dtype)
	end
end

function call(c::ContextSwitch, X::Tensor)
	setContext(device=c.device, deviceId=c.deviceId, dtype=c.dtype)

	if c.device == :cpu
		X = convert(Array{c.dtype},X)
	elseif c.device == :gpu
		X = convert(Knet.KnetArray{c.dtype},X)
	end
	X
end

@debug "Loaded Core modules"


text:
# This file was generated, do not modify it. # hide
acc = (TN + TP) / length(y)
prec = TP /  (TP + FP)
rec  = TP / (TP + FN)
@show round(acc, sigdigits=3)
@show round(prec, sigdigits=3)
@show round(rec, sigdigits=3)

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
using DataStructures
if VERSION >= v"0.5.0-dev+7720"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

@testset "CircularBuffer" begin
    cb = CircularBuffer{Int}(5)
    @test length(cb) == 0
    @test capacity(cb) == 5
    # throws ArgumentError on v0.4 and BoundsError on v0.5 (diverged at 0.5.0-dev+5230)
    if VERSION >= v"0.5.0-dev+5230"
        @test_throws BoundsError first(cb)
    else
        @test_throws ArgumentError first(cb)
    end
    @test isempty(cb) == true
    @test isfull(cb) == false

    push!(cb, 1)
    @test length(cb) == 1
    @test capacity(cb) == 5
    @test isfull(cb) == false

    append!(cb, 2:8)
    @test length(cb) == capacity(cb)
    @test size(cb) == (length(cb),)
    @test isempty(cb) == false
    @test isfull(cb) == true
    @test convert(Array, cb) == Int[4,5,6,7,8]
    @test cb[1] == 4
    @test cb[2] == 5
    @test cb[3] == 6
    @test cb[4] == 7
    @test cb[5] == 8
    @test_throws BoundsError cb[6]
    @test_throws BoundsError cb[3:6]
    @test cb[3:4] == Int[6,7]
    @test cb[[1,5]] == Int[4,8]

    cb[3] = 999
    @test convert(Array, cb) == Int[4,5,999,7,8]

    # Test unshift
    cb = CircularBuffer{Int}(5)  # New, empty one for full test coverage
    for i in -5:5
        unshift!(cb, i)
    end
    arr = convert(Array, cb)
    @test arr == Int[5, 4, 3, 2, 1]
    for (idx, n) in enumerate(5:1)
        @test arr[idx] == n
    end
end


text:
include("file_blacklist.jl")
include("auxiliary_functions.jl")

function get_folders(path)
    folder_list = []
    try
        file_list = readdir(path)
        filtered_list = filter(folder_and_not_blacklisted, file_list)
        folder_list = []

        for file in filtered_list
            path_to_file = string(path, "/",file)
            subfolders = get_folders(path_to_file)

            push!(folder_list, path_to_file)
            folder_list = cat(folder_list, subfolders, dims=1)
        end
    catch e
        @info e
    end

    return folder_list
end

function generate_toc_pages(folder_list)

    for folder in folder_list

        folder_content = []
        try
            folder_content = sort(filter(f -> !folder_blacklisted(f) && !file_blacklisted(f), readdir(folder)), by=x->lowercase(x))
        catch e
            @info e
        end

        start_content = ""
        try
            if(isfile(string(folder, "/index_content.md")))
                start_content = open(f->read(f, String), string(folder, "/index_content.md"))
            end
        catch e
            @info e
        end

        title_raw = last(split(folder, "/"))
        title = apply_formatting(replace(title_raw, "-" => " "))

        html_head_content = "\n\t{{ insert head.html }}\n\t
\n\t\t
$title
\n\t\t\t
\n"
        html_end_content = "\t\t\t\t
\n\t\t\t
\n\t\t
\n\t\t{{ insert page_foot.html }}\n\t\n"
        toc_content = ""

        for file in folder_content
            file_edited = apply_formatting(replace(replace(file, ".md" => ""), "-" => " "))
            if contains(file, ".md")
                toc_content = string(toc_content, "\t\t\t\t ""))\">\n\t\t\t\t
$file_edited

\n\t\t\t\t\t ""))/$(replace(file, ".md" => "")).jpg\">\n\t\t\t\t\t

\n\t\t\t\t\t

\n\t\t\t\t\t

\n\t\t\t\t
\n\t\t\t\t")
            end
        end

        toc_write = string(html_head_content, toc_content, html_end_content)
        open(string(folder, "/index.html"), "w") do io
            write(io, toc_write)
        end

        # toc_content = "@def title = \"$title\" \n@def tags = [\"toc\"] \n# $title\n$start_content\n@@toc-wrapper\n"

        # for md in folder_content
        #     if contains(md, ".md")

        #       toc_content =  string(toc_content, "\t@@image-var-1 [![](/assets/toc-previews/$(replace(folder, "./" => ""))/$(replace(md, ".md" => "")).jpg)](", replace(md, ".md" => ""), ")\n\t\t@@circle-1@@\n\t\t@@circle-2@@\n\t\t@@circle-3@@\n\t@@\n") #\n ### [", apply_formatting(replace(md, ".md" => "")), "]($(replace(md, ".md" => "")))
        #     end
        # end

        # toc_content =  string(toc_content, "@@")

        # open(string(folder, "/index.md"), "w") do io
        #     write(io, toc_content)
        # end
    end
end

toc_page_folders  = get_folders(".")
generate_toc_pages(toc_page_folders)


text:
using Gumbo, HTTP, AbstractTrees

export get_problem_title

get_problem_title(pn) = begin
    p = HTTP.get("https://projecteuler.net/problem=$pn").body |> String |> parsehtml

    for elem in PreOrderDFS(p.root)
        if elem isa HTMLElement
            if tag(elem) == :h2
                return elem.children[1].text
            end
        end
    end
end


text:
file_labels = Dict(
    "free_convection" => "Free convection",
    "strong_wind" => "Strong wind",
    "strong_wind_no_coriolis" => "Strong wind, no rotation",
    "weak_wind_strong_cooling" => "Weak wind, strong cooling",
    "strong_wind_weak_cooling" => "Strong wind, weak cooling",
    "strong_wind_weak_heating" => "Strong wind, weak heating"
)

zs = Dict(
    "u" => 𝒟 -> 𝒟.u.z, "v" => 𝒟 -> 𝒟.v.z, "T"  => 𝒟 -> 𝒟.T.z,
    "uw" => 𝒟 -> 𝒟.uw.z, "vw" => 𝒟 -> 𝒟.vw.z, "wT" => 𝒟 -> 𝒟.wT.z
)

scaling_factor = Dict(
     "u" => 1,
     "v" => 1,
     "T" => 1,
    "uw" => 1e4,
    "vw" => 1e4,
    "wT" => 1e4,
    "T" => 1
)

x_labels = Dict(
     "u" => "U (m/s)",
     "v" => "V (m/s)",
     "T" => "T (m/s)",
    "uw" => "U'W' x 10⁴ (m²/s²)",
    "vw" => "V'W' x 10⁴ (m²/s²)",
    "wT" => "W'T' x 10⁴ (C⋅m/s)",
    "T" => "T (C)"
)

# titles = Dict(
#     "uw" => "Zonal momentum flux, U'W'",
#     "vw" => "Meridional momentum flux, V'W'",
#     "wT" => "Temperature flux, W'T'",
#     "T" => "Temperature, T",
# )

function animate_prediction(xs, name, 𝒟, test_file; filename=name, legend_labels=["" for i in 1:length(xs)], directory="Output")
    filepath = pwd() * "/" * directory * "/"
    isdir(dirname(filepath)) || mkpath(filepath)

    anim = @animate for n in 1:size(xs[1],2)
        x_max = maximum([maximum(x) for x in xs]).*scaling_factor[name]
        x_min = minimum([minimum(x) for x in xs]).*scaling_factor[name]

        fig = plot(xlim=(x_min, x_max), legend=:bottom, size=(400,400), xlabel=x_labels[name], ylabel="Depth (m)")
        for i in reverse(1:length(xs))
            plot!(fig, xs[i][:,n].*scaling_factor[name], zs[name](𝒟), label=legend_labels[i], title=file_labels[test_file]*", $(round(𝒟.t[n]/86400, digits=1)) days", linewidth=4, la=0.5, palette=:Set1_3)
        end

    end

    gif(anim, pwd() * "/$(directory)/$(filename).gif", fps=20)
end


text:
using DandelionSlack
using Base.Test
using FactCheck

include("helpers_test.jl")
include("webapi_test.jl")
include("json_deserialize_test.jl")
include("methods_test.jl")
include("types_test.jl")
include("newtype_test.jl")
include("events_test.jl")
include("rtm_test.jl")
include("throttling_test.jl")


text:
# stores Table as Julia code, precisely as NamedTuple of Arrays.

_convert_to_named_tuple(table) = (; pairs(eachcol(table))...) # e.g. for DataFrames
_convert_to_named_tuple(table:: T) where {T <: NamedTuple{<: AbstractArray}} = table # Do not change anything for NamedTuple of arrays

function _convert_to_named_tuple(table:: T) where {T <: AbstractArray{<: NamedTuple}}
    column_names = keys(first(table))
    data = [getproperty.(table, col) for col in column_names]
    return NamedTuple{column_names}(data)
end

function write_table!(::TableIOInterface.JuliaFormat, filename:: AbstractString, tablename:: AbstractString, table)
    _checktable(table)
    _checktablename(tablename)
    table_as_tuple = _convert_to_named_tuple(table)
    open(filename, "a") do file
        println(file, tablename, " = ", table_as_tuple)
    end
    nothing
end


text:
module Miletus

using Dates, LinearAlgebra
using Dates: days
using Statistics: mean

export CurrencyUnit, CurrencyQuantity, USD, GBP,
 value, process, ivol, fit, fit_ivol,
 SingleStock, CoreModel, CoreForwardModel, GeomBMModel,
 montecarlo, MonteCarloModel, LeastSquaresMonteCarlo

import Optim
import StatsBase: fit
import StatsFuns: normcdf, norminvcdf, normpdf, invsqrt2, sqrthalfπ, invsqrt2π, logistic, logit
import Base.Math: @horner
import ForwardDiff: Dual

extract_derivative(x::Dual) = ForwardDiff.partials(x,1)
extract_derivative(x::Real) = zero(x)

import Reexport: @reexport
@reexport using BusinessDays

const Φ = normcdf
const invΦ = norminvcdf
const ϕ = normpdf

include("daycounts.jl")
include("termstructure.jl")
using .TermStructure

include("currency.jl")
using .Currency

include("utils/math.jl")
include("utils/black.jl")
include("utils/ivol.jl")

include("observables.jl")
include("contracts.jl")
include("process.jl")

include("models/abstractmodel.jl")
include("models/core.jl")
include("models/coreforward.jl")
include("models/geombm.jl")
include("models/binomial.jl")
include("models/binomial_ivol.jl")

include("models/montecarlo.jl")
include("models/lsmc.jl")

include("print.jl")

include("greeks.jl")
export vega, delta, rho, greek

end # module


text:
"""
"""
function pair_arrays(a::AbstractArray,b::AbstractArray)
  ArrayPair(a,b)
end

struct ArrayPair{T,N,A,B} <:AbstractArray{T,N}
  a::A
  b::B

  function ArrayPair(a::AbstractArray,b::AbstractArray)
    @assert length(a) == length(b)
    A = typeof(a); B = typeof(b)
    Ta = eltype(a); Tb = eltype(b)
    T = Tuple{Ta,Tb}
    N = ndims(a)
    new{T,N,A,B}(a,b)
  end

end

Base.size(p::ArrayPair) = size(p.a)

Base.IndexStyle(::Type{ArrayPair{T,N,A,B}}) where {T,N,A,B} = IndexStyle(A)

@inline Base.getindex(p::ArrayPair,i::Integer) = (p.a[i], p.b[i])

@inline Base.getindex(p::ArrayPair{T,N},i::Vararg{Integer,N}) where {T,N} = (p.a[i...], p.b[i...])

uses_hash(::Type{<:ArrayPair}) = Val{true}()

function array_cache(hash,p::ArrayPair)
  ca = array_cache(hash,p.a)
  cb = array_cache(hash,p.b)
  (ca,cb)
end

@inline function getindex!(cache,p::ArrayPair,i)
  ca, cb = cache
  ai = getindex!(ca,p.a,i)
  bi = getindex!(cb,p.b,i)
  (ai,bi)
end

"""
"""
function unpair_arrays(pair::AbstractArray{<:Tuple})
  p = testitem(pair)
  N = length(p)
  @assert N > 0
  _unpair_arrays(pair,Val{1}(),Val{N}())
end

function unpair_arrays(pair::ArrayPair)
  pair.a, pair.b
end

function  _unpair_arrays(pair,::Val{n},::Val{N}) where {n,N}
  pn = UnpairedArray(Val{n}(),pair)
  ps = _unpair_arrays(pair,Val{n+1}(),Val{N}())
  (pn, ps...)
end

function  _unpair_arrays(pair,::Val{N},::Val{N}) where N
  pn = UnpairedArray(Val{N}(),pair)
  (pn, )
end

struct UnpairedArray{I,T,N,P<:AbstractArray} <: AbstractArray{T,N}
  pair::P
  function UnpairedArray(::Val{I},pair::AbstractArray{<:Tuple,N}) where {I,N}
    p = testitem(pair)
    T = typeof(p[I])
    P = typeof(pair)
    new{I,T,N,P}(pair)
  end
end

Base.size(a::UnpairedArray) = size(a.pair)

Base.IndexStyle(::Type{UnpairedArray{I,T,N,P}}) where {T,N,I,P} = IndexStyle(P)

@inline function Base.getindex(a::UnpairedArray{I},i::Integer) where I
  a.pair[i][I]
end

@inline function Base.getindex(a::UnpairedArray{I,T,N},i::Vararg{Integer,N}) where {I,T,N}
  a.pair[i...][I]
end

uses_hash(::Type{<:UnpairedArray}) = Val{true}()

function array_cache(hash,p::UnpairedArray)
  array_cache(hash,p.pair)
end

@inline function getindex!(cache,a::UnpairedArray{I},i) where I
  p = getindex!(cache,a.pair,i)
  p[I]
end


text:
value(x::Type{Tracker.TrackedReal{T}}) where T = T
value(x::Type{Tracker.TrackedArray{T,N,A}}) where {T,N,A} = Array{T,N}
value(x::Tracker.TrackedReal)  = x.data
value(x::Tracker.TrackedArray) = x.data

@inline fastpow(x::Tracker.TrackedReal, y::Tracker.TrackedReal) = x^y
@inline Base.any(f::Function,x::Tracker.TrackedArray) = any(f,Tracker.data(x))

# Support adaptive with non-tracked time
@inline function ODE_DEFAULT_NORM(u::Tracker.TrackedArray,t) where {N}
  sqrt(sum(abs2,value(u)) / length(u))
end
@inline function ODE_DEFAULT_NORM(u::AbstractArray{<:Tracker.TrackedReal,N},t) where {N}
  sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
end
@inline function ODE_DEFAULT_NORM(u::Array{<:Tracker.TrackedReal,N},t) where {N}
  sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
end
@inline ODE_DEFAULT_NORM(u::Tracker.TrackedReal,t) = abs(value(u))

# Support TrackedReal time, don't drop tracking on the adaptivity there
@inline function ODE_DEFAULT_NORM(u::Tracker.TrackedArray,t::Tracker.TrackedReal) where {N}
  sqrt(sum(abs2,u) / length(u))
end
@inline function ODE_DEFAULT_NORM(u::AbstractArray{<:Tracker.TrackedReal,N},t::Tracker.TrackedReal) where {N}
  sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip(u,Iterators.repeated(t))) / length(u))
end
@inline function ODE_DEFAULT_NORM(u::Array{<:Tracker.TrackedReal,N},t::Tracker.TrackedReal) where {N}
  sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip(u,Iterators.repeated(t))) / length(u))
end
@inline ODE_DEFAULT_NORM(u::Tracker.TrackedReal,t::Tracker.TrackedReal) = abs(u)

function solve_up(prob::DiffEqBase.DEProblem,sensealg::Union{AbstractSensitivityAlgorithm,Nothing},u0::Tracker.TrackedArray,p::Tracker.TrackedArray,args...;kwargs...)
  Tracker.track(solve_up,prob,sensealg,u0,p,args...;kwargs...)
end

function solve_up(prob::DiffEqBase.DEProblem,sensealg::Union{AbstractSensitivityAlgorithm,Nothing},u0::Tracker.TrackedArray,p,args...;kwargs...)
  Tracker.track(solve_up,prob,sensealg,u0,p,args...;kwargs...)
end

function solve_up(prob::DiffEqBase.DEProblem,sensealg::Union{AbstractSensitivityAlgorithm,Nothing},u0,p::Tracker.TrackedArray,args...;kwargs...)
  Tracker.track(solve_up,prob,sensealg,u0,p,args...;kwargs...)
end

Tracker.@grad function solve_up(prob,sensealg::Union{Nothing,AbstractSensitivityAlgorithm},
                                u0,p,args...;
                                kwargs...)
  _solve_adjoint(prob,sensealg,Tracker.data(u0),Tracker.data(p),args...;kwargs...)
end


text:
# tuple
(kernel::GateKernel)(qubits::Int...) = QGate(kernel)(qubits...)
(gate::VarGate)(qubit::Int...) = (gate, [qubit...])
(block::QGateBlock)(qubit::Int...) = (block,[qubit...])
getindex(kernel::GateKernel, param...) = VarGate(kernel, [param...])
getindex(gate::VarGate, params...) = (gate.param = [params...]; return gate)
(gate::ConstGate)(qubit::Int...) = (gate, [qubit...])

#array
(kernel::GateKernel)(qubits::Vector{Int}) = QGate(kernel)(qubits)
(gate::VarGate)(qubit::Vector{Int}) = (gate, qubit)
(block::QGateBlock)(qubit::Vector{Int}) = (block, qubit)
getindex(kernel::GateKernel, param::Vector) = VarGate(kernel, param)
getindex(gate::VarGate, params::Vector) = (gate.param = params; return gate)
(gate::ConstGate)(qubit::Vector{Int}) = (gate, qubit)

apply!(state::QState, tup:: Tuple{Union{QGate, QGateBlock}, Vector{Int}} ; kwargs...) = apply!(state, tup[1], tup[2]; kwargs...)

function measure!(state::QState; kwargs...)
	qubits = get(kwargs, :qubits, nothing)
	qubit = get(kwargs, :qubit, nothing)

	shots = get(kwargs, :shots, 1024)
	if qubits ==nothing
		qubit == nothing && (throw(ArgumentError("Qubit(s) to be measured not specified")))
		return measure!(state, qubits, shots; kwargs...)
	else
		qubit != nothing && (throw(ArgumentError("Both qubits and qubit are specified")))
		return measure!(state, qubit, shots; kwargs...)
	end
end

randomQCircuit(;size::Int, depth::Int = rand(1:10)) = QCircuit(randomQGateBlock(size,depth), size)


text:
# 1D quadrature

integrate(poly::Poly, a, b) = (pp = polyint(poly); fapply(pp,b) - fapply(pp,a))
integrate(x::Number, a, b) = x*(b-a)

const gp = Base.gauss(Float64, 7)
function integrate(f::Function,a, b)
    gx, gw = gp
    @assert a == -1 && b == 1

    result = gw[1]*f(gx[1])
    for i in 2:length(gw)
        result += gw[i]*f(gx[i])
    end
    result
end

do_quad_ref(f::Poly,p::DG1D) = integrate(f,-1,1)

# ElemJacobian
elemJ(x::DG1D,k) = (1//2)*(x.mesh[k,:r] - x.mesh[k,:l])

# 2D quadrature

# Get quadrature weights (from lecture notes)
function quadp(p)
    if p == 1
        return (
            [Vertex2(1//3,1//3)],[0.5])
    elseif p == 2
        return (
            [Vertex2(2//3,1//6),Vertex2(1//6,2//3),
             Vertex2(1//6,1//6)],
             [1//6,1//6,1//6])
    elseif p == 3
        return (
            [Vertex2(0.1550510257, 0.1785587282),
             Vertex2(0.6449489742, 0.0750311102),
             Vertex2(0.1550510257, 0.6663902460),
             Vertex2(0.6449489742, 0.2800199154)],
             [0.1590206908,0.0909793091,
              0.1590206908,0.0909793091])
    elseif p == 4
        return (
            [Vertex2(0.4459484909, 0.4459484909),
             Vertex2(0.1081030181, 0.4459484909),
             Vertex2(0.4459484909, 0.1081030181),
             Vertex2(0.0915762135, 0.0915762135),
             Vertex2(0.8168475729, 0.0915762135),
             Vertex2(0.0915762135, 0.8168475729)],
             [0.1116907948,0.1116907948,0.1116907948,
              0.0549758718,0.0549758718,0.0549758718])
    else
        error("Unimplemented")
    end
end

function do_quad(f,porder)
    result = zero(typeof(f(zero(Vertex2))))
    for (x,w) in zip(quadp(porder)...)
        result += w*f(x)
    end
    result
end

do_quad_ref(f::Function,p::Galerkin2D) = do_quad(f,clamp(porder(p)*2,1,4))
do_quad_ref(f::Function,p::Galerkin2D,edge::Edge) = do_quad_ref(p.mesh,edge,f,porder(p)*2)

do_quad(c::Cell2D, f, p) = det(Ak(c))*do_quad(x->f(𝜒⁻¹(c,x)),p)
function elemJ(p::Galerkin2D,c)
    d = det(Ak(c))
    d < 0 && error("Misoriented cell")
    d
end

#point = 𝜒⁻¹(edge,x)
#@show p

macro do_quad_ref2d_Δ1(fargs, var, block)
    esc(quote
        args = $fargs
        p = args[1]
        edge = args[2]
        points = args[3]
        gx, gw = points
        # Manual loop peeling
        x,w = gx[1], gw[1]
        $var = x
        result = w*begin
            $block
        end
        for i in 2:length(gw)
            x,w = gx[i],gw[i]
            $var = x
            result += w*begin
                $block
            end
        end
        (norm(edge)/2)*result
    end)
end

function do_quad_ref(m::Meshes.Mesh, edge::Edge, f, porder; points=Base.gauss(Float64,porder))
    gx, gw = points
    result = 0.0
    for i in 1:length(gw)
        x,w = gx[i],gw[i]
        result += w*f(x)
    end
    (norm(edge)/2)*result
end


text:
using RoME, Distributions, TransformUtils
# dev test script for Pose2 and Pose2Point2BearingRange using bss sampler
using RoMEPlotting

const TU = TransformUtils

# drive in a circle (copy hexagonal increase)

# start with an empty factor graph object
N = 100
fg = initfg()

# Add the first pose :x0
addVariable!(fg, :x0, Pose2)

# Add at a fixed location PriorPose2 to pin :x0 to a starting location
addFactor!(fg, [:x0], PriorPose2(MvNormal(zeros(3), 0.01*eye(3))))

# Drive around in a hexagon
n_sides = 16
side_length = 3.0
for i in 0:(n_sides-1)
  psym = Symbol("x$i")
  nsym = Symbol("x$(i+1)")
  addVariable!(fg, nsym, Pose2)
  pp = Pose2Pose2(MvNormal([side_length;0;2pi/n_sides], diagm([0.01;0.01;0.001].^2)))
  addFactor!(fg, [psym;nsym], pp )
end

# two acoustic beacons

addVariable!(fg, :l1, Point2)
addVariable!(fg, :l2, Point2)

pL1 = IIF.Prior(MvNormal([50.0; 50.0], 0.1^2*eye(2)))
pL2 = IIF.Prior(MvNormal([-50.0; 50.0], 0.1^2*eye(2)))
addFactor!(fg,[:l1],pL1)
addFactor!(fg,[:l2],pL2)

# doautoinit!(fg, :l1, N=N)
# doautoinit!(fg, :l2, N=N)
ensureAllInitialized!(fg)

# build probabilities from parametric Gaussian but use bss sampler

# calculate the range bearings from each pose to

GTposes = KDE.getKDEMax.(getBelief.(fg, ls(fg)[1]))
GTlandm = KDE.getKDEMax.(getBelief.(fg, ls(fg)[2]))

brDict = Dict{Symbol, Dict{Symbol,Any}}()

for poseidx in 1:length(GTposes), lmidx in 1:length(GTlandm)
    gtpose = GTposes[poseidx]
    landm = GTlandm[lmidx]

    px, py = gtpose[1], gtpose[2]
    lx, ly = landm[1], landm[2]
    delta = [landm[1]-gtpose[1];landm[2]-gtpose[2]]
    d = norm(delta)
    wyaw = atan2(delta[2], delta[1])
    # must still incorporate pose yaw
    ryaw = TU.wrapRad(wyaw - gtpose[3])

    # generate range measurement samplable belief
    r_rng = collect(linspace(0, 2*d, 128)) # range vector from 0 to two times the distance
    n_rng = Normal(d, 0.1)
    w_rng = pdf(n_rng, r_rng) # weights
    # TODO: check for weight renormalization
    bss_rng = IIF.AliasingScalarSampler(r_rng, w_rng)

    r_bear = collect(linspace(-pi, pi, 512)) # range vector from 0 to t
    n_bear = Normal(ryaw, 0.1)
    w_bear = pdf(n_bear, r_bear)
    bss_bear = IIF.AliasingScalarSampler(r_bear, w_bear)

    posesym = ls(fg)[1][poseidx]
    lmsym = ls(fg)[2][lmidx]
    if !haskey(brDict, posesym)
      brDict[posesym] = Dict{Symbol, Any}()
    end
    brDict[posesym][lmsym] = Pose2Point2BearingRange(bss_bear, bss_rng)
end

for posesym in [:x1], lmsym in [:l1;:l2]
  addFactor!(fg, [posesym;lmsym], brDict[posesym][lmsym])
end

drawGraph(fg,engine="dot")

# t = string(now())
# savejld(fg, file="presolve_$t.jld")
solveTree!(fg)
# savejld(fg, file="postsolve_$t.jld")

using Gadfly
# show that it all works
# poses
pl = drawPoses(fg, spscale=1.0)

# landmarks
pl = drawPosesLandms(fg, spscale=2.0)
Gadfly.draw(Gadfly.PDF("/tmp/test3.pdf", 20cm, 10cm),pl)

stuff = plotLocalProduct(fg, :x1)

# then use this to validate the Sandshark pipeline

# pvt: make this an example or test somewhere?
# loosen


text:
using NearestNeighbors, PDMats

export AnalogForecasting

"""
    AnalogForecasting(k, xt, catalog)

parameters of the analog forecasting method

- k            : number of analogs
- neighborhood : global analogs
- catalog      : catalog with analogs and successors
- regression   : (:locally_constant, :increment, :local_linear)
- sampling     : (:gaussian, :multinomial)

"""
struct AnalogForecasting <: AbstractForecasting

    kdt::KDTree
    k::Int64 # number of analogs
    neighborhood::BitArray{2}
    catalog::Catalog
    regression::Symbol
    sampling::Symbol

    function AnalogForecasting(
        k::Int64,
        xt::TimeSeries,
        catalog::Catalog;
        regression = :local_linear,
        sampling = :gaussian,
        leafsize = 50,
    )

        neighborhood = trues((xt.nv, xt.nv)) # global analogs

        kdt = KDTree(catalog.analogs, leafsize = leafsize)

        new(kdt, k, neighborhood, catalog, regression, sampling)

    end

    function AnalogForecasting(
        k::Int64,
        xt::TimeSeries,
        catalog::Catalog,
        neighborhood::BitArray{2},
        regression::Symbol,
        sampling::Symbol,
        leafsize = 50,
    )

        kdt = KDTree(catalog.analogs, leafsize = leafsize)

        new(kdt, k, neighborhood, catalog, regression, sampling)

    end

end

"""
    Apply the analog method on catalog of historical data
    to generate forecasts.
"""
function (forecasting::AnalogForecasting)(x::Array{Float64,2})

    nv, np = size(x)
    xf = zeros(Float64, (nv, np))
    xf_mean = zeros(Float64, (nv, np))
    ivar = [1]
    condition = true

    # global analog forecasting
    index_knn, dist_knn = knn(forecasting.kdt, x, forecasting.k)

    # parameter of normalization for the kernels
    λ = median(Iterators.flatten(dist_knn))
    @assert λ > 0.0
    weights = [ones(forecasting.k) for i = 1:np]

    while condition

        if all(forecasting.neighborhood)
            ivar_neighboor = 1:nv
            ivar = 1:nv
            condition = false
        else
            ivar_neighboor = findall(vec(forecasting.neighborhood[ivar, :]))
        end

        # compute weights
        for ip = 1:np
            for k = 1:forecasting.k
                weights[ip][k] = exp(-dist_knn[ip][k]^2 / λ^2)
            end
            s = sum(weights[ip])
            for k = 1:forecasting.k
                weights[ip][k] /= s
            end

            if any(isnan.(weights[ip]))
                @show ip
                @show λ
                @show dist_knn[ip]
                @show weights[ip]
                throw(" Some nan values in weights ")
            end

        end

        # initialization
        xf_tmp = zeros(Float64, (last(ivar), forecasting.k))

        if forecasting.regression == :local_linear
            local_linear = LocalLinear(forecasting.k, ivar, ivar_neighboor)
        end

        for ip = 1:np
            if forecasting.regression == :locally_constant

                # compute the analog forecasts
                xf_tmp[ivar, :] .= forecasting.catalog.successors[ivar, index_knn[ip]]
                # weighted mean and covariance
                xf_mean[ivar, ip] = sum(xf_tmp[ivar, :] .* weights[ip]', dims = 2)
                Exf = xf_tmp[ivar, :] .- xf_mean[ivar, ip]
                cov_xf = Symmetric(
                    1.0 ./ (1.0 .- sum(weights[ip] .^ 2)) .* (Exf .* weights[ip]') * Exf',
                )

            elseif forecasting.regression == :increment # method "locally-incremental"
                # compute the analog forecasts
                xf_tmp[ivar, :] .= (
                    x[ivar, ip] .+ forecasting.catalog.successors[ivar, index_knn[ip]] .-
                    forecasting.catalog.analogs[ivar, index_knn[ip]]
                )

                # weighted mean and covariance
                xf_mean[ivar, ip] = sum(xf_tmp[ivar, :] .* weights[ip]', dims = 2)
                Exf = xf_tmp[ivar, :] .- xf_mean[ivar, ip]
                cov_xf = Symmetric(
                    1.0 ./ (1.0 .- sum(weights[ip] .^ 2)) .* (Exf .* weights[ip]') * Exf',
                )

            elseif forecasting.regression == :local_linear

                # define analogs, successors and weights

                # compute centered weighted mean and weighted covariance
                cov_xf = compute(
                    local_linear,
                    x,
                    xf_tmp,
                    xf_mean,
                    ip,
                    forecasting.catalog.analogs[ivar_neighboor, index_knn[ip]],
                    forecasting.catalog.successors[ivar, index_knn[ip]],
                    weights[ip],
                )

                if cov_xf == 0 # error in pinv back to locally constant

                    xf_tmp[ivar, :] .= forecasting.catalog.successors[ivar, index_knn[ip]]
                    # weighted mean and covariance
                    xf_mean[ivar, ip] = sum(xf_tmp[ivar, :] .* weights[ip]', dims = 2)
                    Exf = xf_tmp[ivar, :] .- xf_mean[ivar, ip]
                    cov_xf = Symmetric(
                        1.0 ./ (1.0 .- sum(weights[ip] .^ 2)) .* (Exf .* weights[ip]') *
                        Exf',
                    )
                else

                    # constant weights for local linear
                    weights[ip] .= 1.0 / forecasting.k
                end

            else

                @error "regression: locally_constant, :increment, :local_linear "
            end

            # Gaussian sampling
            if forecasting.sampling == :gaussian

                ϵ = 1e-8
                while !isposdef(cov_xf)
                    cov_xf .= ensure_pos_sym(cov_xf, ϵ = ϵ)
                    ϵ *= 10
                end

                # random sampling from the multivariate Gaussian distribution
                d = MvNormal(xf_mean[ivar, ip], cov_xf)
                xf[ivar, ip] .= rand!(d, xf[ivar, ip])

                # Multinomial sampling
            elseif forecasting.sampling == :multinomial

                # random sampling from the multinomial distribution
                # of the weights
                # this speedup is due to Peter Acklam
                cumprob = cumsum(weights[ip])
                R = rand()
                M = 1::Int64
                N = length(cumprob)
                for i = 1:N-1
                    M += R > cumprob[i]
                end
                igood = M
                xf[ivar, ip] .= xf_tmp[ivar, igood]

            else

                @error " sampling = :gaussian or :multinomial"

            end

        end

        if all(ivar .== [nv]) || length(ivar) == nv
            condition = false
        else
            ivar .+= 1
        end

    end

    xf, xf_mean

end


text:
module Evolve

import Laplacian

# Return the x values associated with a given choice of grid divisions and grid
# width
function X(L :: Real, N :: Integer, xmin :: Real)
    return (L/N)*[0:N]+xmin
end
X(L :: Real, N :: Integer) = X(L, N, 0)

# Return a time evolution matrix A for a given potential energy vector and box width.
function A(V :: Vector, L :: Real)

    ϵ = L / (length(V) - 1)

    AA  = SymTridiagonal(0.25(1+ϵ^2*V)im+0.5,-0.125im*ones(length(V)-1))
    AA.dv[1] = AA.dv[end] = 1
    AA.ev[1] = AA.ev[end] = zero(typeof(AA[1,1]))
    return AA

end
A(V :: Function, L :: Real, N :: Integer, xmin :: Real) = A(V(X(L,N,xmin)),L)
A(V :: Function, L :: Real, N :: Integer) = A(V, L, N, 0)

# Return a reverse time evolution matrix B for a given potential energy vector
# and box width.
function B(V :: Vector, L :: Real)

    ϵ = L / (length(V) - 1)

    AA  = SymTridiagonal(0.25(1+ϵ^2*V)im-0.5,-0.125im*ones(length(V)-1))
    AA.dv[1] = AA.dv[end] = 1
    AA.ev[1] = AA.ev[end] = zero(typeof(AA[1,1]))
    return AA

end
B(V :: Function, L :: Real, N :: Integer, xmin :: Real) = B(V(X(L,N,xmin)),L)
B(V :: Function, L :: Real, N :: Integer) = B(V, L, N, 0)

# Return a constructor for an update matrix A with the specified potential energy
# function
function A(V :: Function)
    return (L :: Real, N :: Integer, xmin :: Real) -> A(V, L, N, xmin)
end

# Return the Hamiltonian for a given potential energy vector and box width.
function H(V :: Vector, L :: Real)

    ϵ = L / (length(V) - 1)

    AA = SymTridiagonal((ϵ^-2)-V, -0.5(ϵ^-2)ones(length(V)-1))
    AA.dv[1] = AA.dv[end] = 1
    AA.ev[1] = AA.ev[end] = zero(typeof(AA[1,1]))
    return AA

end
H(V :: Function, L :: Real, N :: Integer, xmin :: Real) = H(V(X(L,N,xmin)),L)
H(V :: Function, L :: Real, N :: Integer) = H(V, L, N, 0)

# Method for setting ϕ₀
function Φ(k₀ :: Real, σ :: Real, x₀ :: Real)

    return (L :: Real, N :: Integer, xmin :: Real) -> begin
        xs = (L/N) .* [0:N] + xmin
        swp = map((x) -> exp(im*k₀*x - (x-x₀)^2 / (2σ^2)), xs)
        swp *= (1 / (π*σ^2))^(1/4)          # normalize
        swp[1] = swp[end] = 0               # enforce ∞ square well boundary
        return swp
    end

end
Φ(k₀ :: Real, σ :: Real) = Φ(k₀, σ, 0)          # start in the center

# Function for evolving a given state vector through n additional steps
# (first step not in result). Keep running until done or until test fails.
function evolve(A :: AbstractArray, ϕ :: Vector, n :: Integer, test :: Function)

    ϕs = zeros(Complex, length(ϕ), n)

    # Solve for χ using Thomas algorithm
    χ  = tommy(A, ϕ)
    ϕs[:,1] = χ - ϕ

    for j in 1:(n-1)
        if test(ϕs[:,j])
            ϕs = ϕs[:,1:j]
            break
        end
        χ = tommy(A, ϕs[:,j])
        ϕs[:,j+1] = χ - ϕs[:,j]
    end

    return ϕs

end
evolve(A :: AbstractArray, ϕ :: Vector, n :: Integer) = evolve(A, ϕ, n, x->true)

# Method for evolving a state vector by passing functions
evolve(A :: Function, ϕ :: Function, n :: Integer, L :: Real, N :: Integer, xmin :: Real) = evolve(A(L,N,xmin), ϕ(L,N,xmin), n)

# Solve Ax = b where A is tridiagonal using the Thomas algorithm
# based on MATLAB code: http://www.mathworks.com/matlabcentral/fileexchange/1359-fast-tridiagonal-system-solver/content/thomas.m
function tommy(a :: Vector, b :: Vector, c :: Vector, d :: Vector)

    n = length(a)
    m = zeros(typeof(a[1]), n)
    l = zeros(typeof(a[1]), n-1)
    y = zeros(typeof(a[1]), n)
    x = zeros(typeof(a[1]), n)

    m[1] = a[1]
    y[1] = d[1]

    # Perform LU factorization and forward substitution to solve Ly = d for y
    for i in 2:n
        l[i-1] = c[i-1] / m[i-1]
        m[i] = a[i] - l[i-1] * b[i-1]

        y[i] = d[i] - l[i-1] * y[i-1]
    end

    x[n] = y[n] / m[n]

    # Begin backward substitution to solve Ux = y for x
    for i in (n-1):-1:1
        x[i] = (y[i] - b[i] * x[i+1]) / m[i]
    end

    return x

end

tommy(A :: Matrix, d :: Vector) = tommy(diag(A), diag(A, -1), diag(A, +1), d)
tommy(A :: Tridiagonal, d :: Vector) = tommy(A.d, A.du, A.dl, d)
tommy(A :: SymTridiagonal, d :: Vector) = tommy(A.dv, A.ev, A.ev, d)

# Create a pocket function centered at zero
function pocket(L :: Real, h :: Real , b :: Real, s :: Real, w :: Real, xmin :: Real, x :: Real)
    if x < (xmin + 0.5L - 0.5w - s) || x > (xmin + 0.5L + 0.5w +s)
        V = 0
    elseif x < (xmin + 0.5L - 0.5w) || x > (xmin + 0.5L + 0.5w)
        V = h
    else
        V = b
    end
    return V
end

end


text:
using PyPlot
using Logging

include("psmsdst.jl")
greedygen = include("solvers/softgreedy.jl")

function analyze()
    problemname = "PSMSDSTdata/PSMSDST6.txt"

    problem = parsepsmsdst(problemname)

    hardnesses = 5:0.25/2.0:8
    costs = Vector()

    for hardness in hardnesses
        solver,name = greedygen(attempts=0,hardness=hardness)

        solution = solver(problem)

        starttimes, finishtimes, cost = evaluatesolution(problem,solution)

        @info "Looking for best hardness" hardness cost

        push!(costs,cost)
    end

    plt.plot(hardnesses,costs)
    plt.show()
end

if abspath(PROGRAM_FILE) == @__FILE__
    analyze()
end

max_stars_repo_name: mkualquiera/huepl

max_issues_repo_name: mkualquiera/huepl

max_forks_repo_name: mkualquiera/huepl


text:
# This gets run by setup.sh

using Pkg

Pkg.activate(@__DIR__)
Pkg.add(name="ExponentialUtilities", version=ARGS[1])

using ExponentialUtilities


text:
using Revise, ApproxOperator

elements, nodes = importmsh("./msh/patchtest.msh")
nₚ = length(nodes[:x])
nₑ = length(elements["Ω"])

set𝓖!(elements["Ω"],:TriGI3)
# set𝓖!(elements["Γᵗ₁"],:SegGI2)
# set𝓖!(elements["Γᵗ₂"],:SegGI2)
set𝓖!(elements["Γᵍ"],:SegGI2)
elements["Γᵍ"] = Element{:Tri3}(elements["Ω"],elements["Γᵍ"])

E = 1.0
ν = 0.0
prescribe!(elements["Ω"],:b₁,(x,y,z)->0.0)
prescribe!(elements["Ω"],:b₂,(x,y,z)->0.0)
# prescribe!(elements["Γᵗ₁"],:t₁,(x,y,z)->E/(1-ν))
# prescribe!(elements["Γᵗ₁"],:t₂,(x,y,z)->E/(1+ν))
# prescribe!(elements["Γᵗ₂"],:t₁,(x,y,z)->E/(1+ν))
# prescribe!(elements["Γᵗ₂"],:t₂,(x,y,z)->E/(1-ν))
prescribe!(elements["Γᵍ"],:g₁,(x,y,z)->1.0+x+y)
prescribe!(elements["Γᵍ"],:g₂,(x,y,z)->1.0+x+y)
prescribe!(elements["Γᵍ"],:n₁₁,(x,y,z)->1.0)
prescribe!(elements["Γᵍ"],:n₁₂,(x,y,z)->0.0)
prescribe!(elements["Γᵍ"],:n₂₂,(x,y,z)->1.0)

coefficient = (:E=>E,:ν=>ν,:α=>1e7)
ops = [Operator(:∫∫εᵢⱼσᵢⱼdxdy,coefficient...),
       Operator(:∫∫vᵢbᵢdxdy,coefficient...),
       Operator(:∫vᵢtᵢds,coefficient...),
       Operator(:∫σᵢⱼnⱼgᵢds,coefficient...),
       Operator(:∫vᵢgᵢds,coefficient...),
       Operator(:Hₑ_PlaneStress,coefficient...)]

k = zeros(2*nₚ,2*nₚ)
f = zeros(2*nₚ)

ops[1](elements["Ω"],k)
# ops[3](elements["Γᵗ₁"],f)
# ops[3](elements["Γᵗ₂"],f)
ops[4](elements["Γᵍ"],k,f)
# ops[5](elements["Γᵍ"],k,f)

d = k\f

d₁ = zeros(nₚ)
d₂ = zeros(nₚ)
d₁ .= d[1:2:2*nₚ-1]
d₂ .= d[2:2:2*nₚ]
push!(nodes,:d₁=>d₁,:d₂=>d₂)
prescribe!(elements["Ω"],:u,(x,y,z)->1.0+x+y)
prescribe!(elements["Ω"],:v,(x,y,z)->1.0+x+y)
prescribe!(elements["Ω"],:∂u∂x,(x,y,z)->1.0)
prescribe!(elements["Ω"],:∂v∂x,(x,y,z)->1.0)
prescribe!(elements["Ω"],:∂u∂y,(x,y,z)->1.0)
prescribe!(elements["Ω"],:∂v∂y,(x,y,z)->1.0)
h1,l2 = ops[6](elements["Ω"])

max_stars_repo_name: jmwjc/ApproxOperator

max_issues_repo_name: jmwjc/ApproxOperator

max_forks_repo_name: jmwjc/ApproxOperator