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

    The main object defined in this file is: BiBlock.
    It represents a single block and is the simplest composite unit
    that can be used for practical settings of inference and smoothing.
    Additionally, it contains the option to sample with a preconditioned
    Crank-Nicolson step.

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

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

        DEFINITION

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

"""
    mutable struct BiBlock{L,TGP,TGPl,TW,TWn,TX}
        b::Block{L,TGP,TGPl,TW,TWn,TX}
        b°::Block{L,TGP,TGPl,TW,TWn,TX}
        ρ::Float64
        accpt_history::Vector{Bool}
    end

Composite unit that allows for sampling of a single block. It provides two
`Block`s: one proposal `b°`, one accepted `b` that can be used for smoothing or
inference problems. `ρ` is a memory parameter of the preconditioned
Crank-Nicolson scheme and `accpt_history` stores the history of accept/reject
decisions (useful for MCMC).

    function BiBlock(
        sp::SamplingPair,
        range::UnitRange{Int64},
        ρ=0.0,
        last_block=false,
        ll_hist_len=0
    )

Base constructor.
"""
mutable struct BiBlock{L,TGP,TGPl,TW,TWn,TX}
    b::Block{L,TGP,TGPl,TW,TWn,TX}
    b°::Block{L,TGP,TGPl,TW,TWn,TX}
    ρ::Float64
    accpt_history::Vector{Bool}

    function BiBlock(
            sp::SamplingPair{TGP,TGPl,TW,TWn,TX},
            range::UnitRange{Int64},
            ρ=0.0,
            last_block=false,
            ll_hist_len=0
        ) where {TGP,TGPl,TW,TWn,TX}
        b = Block(sp.u, range, last_block, ll_hist_len)
        b° = Block(sp.u°, range, last_block, ll_hist_len)
        new{last_block,TGP,TGPl,TW,TWn,TX}(
            b, b°, ρ,
            Vector{Bool}(undef, ll_hist_len)
        )
    end
end

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

        IMPUTATION

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

"""
    draw_proposal_path!(bb::BiBlock)

Sample a proposal path, compute log-likelihood along the way. Assumes
`bb.b.XX[1].x[1]` is a starting point. Uses preconditioned Crank-Nicolson scheme
with memory parameter set as `bb.ρ`.
"""
function draw_proposal_path!(bb::BiBlock) end

function draw_proposal_path!(bb::BiBlock{false})
    success, bb.b°.ll = _draw_proposal_path!(bb)
    success || return false # anything?
    y1 = bb.b°.XX[end-1].x[end]
    success, ll°_last = _draw_proposal_path_last_segment!(bb, y1)
    bb.b°.ll += ll°_last
    success
end

function draw_proposal_path!(bb::BiBlock{true})
    success, bb.b°.ll = _draw_proposal_path!(bb)
    success
end

function _draw_proposal_path!(bb::BiBlock)
    rand!(
        bb.b.PP, bb.b°.XX, bb.b°.WW, bb.b.WW, bb.ρ, _LL, bb.b.XX[1].x[1];
        Wnr=bb.b°.Wnr
    )
end

function _draw_proposal_path_last_segment!(bb::BiBlock, y1)
    rand!(
        bb.b.P_last[1], bb.b°.XX[end], bb.b°.WW[end], bb.b.WW[end], bb.ρ, _LL,
        y1; Wnr=bb.b°.Wnr
    )
end

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

        ACCEPT/REJECT DECISION

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

"""
    accept_reject_proposal_path!(bb::BiBlock, mcmciter)

Accept/reject decision of the Metropolis-Hastings algorithm for the step of path
imputation.
"""
function accept_reject_proposal_path!(bb::BiBlock, mcmciter)
    accepted = rand(Exponential(1.0)) > -(bb.b°.ll - bb.b.ll)
    accepted && swap_paths!(bb)
    set_accepted!(bb, mcmciter, accepted)
    save_ll!(bb, mcmciter)
    accepted && swap_ll!(bb)
end

"""
    set_accepted!(bb::BiBlock, i::Int, v)

Commit the accept/reject decision `v` to acceptance history of BiBlock `b` at
the position `i`.
"""
set_accepted!(bb::BiBlock, i::Int, v) = (bb.accpt_history[i] = v)

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

        SWAPS

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

"""
    swap_paths!(bb::BiBlock)

Swap `XX` and `WW` containers between proposal-acceptance pair.
"""
function swap_paths!(bb::BiBlock)
    swap_XX!(bb)
    swap_WW!(bb)
end

"""
    swap_XX!(bb::BiBlock)

Swap `XX` containers between proposal-acceptance pair.
"""
function swap_XX!(bb::BiBlock)
    for i in eachindex(bb.b.XX)
        bb.b.XX[i], bb.b°.XX[i] = bb.b°.XX[i], bb.b.XX[i]
    end
end

"""
    swap_WW!(bb::BiBlock)

Swap `WW` containers between proposal-acceptance pair.
"""
function swap_WW!(bb::BiBlock)
    for i in eachindex(bb.b.WW)
        bb.b.WW[i], bb.b°.WW[i] = bb.b°.WW[i], bb.b.WW[i]
    end
end

"""
    swap_PP!(bb::BiBlock)

Swap `PP` containers (including PP_last) between proposal-acceptance pair.
"""
function swap_PP!(bb::BiBlock) end

function swap_PP!(bb::BiBlock{true})
    _swap_PP!(bb)
end

function swap_PP!(bb::BiBlock{false})
    _swap_PP!(bb)
    bb.b.P_last[1], bb.b°.P_last[1] = bb.b°.P_last[1], bb.b.P_last[1]
    bb.b.P_excl[1], bb.b°.P_excl[1] = bb.b°.P_excl[1], bb.b.P_excl[1]
    for i in eachindex(bb.b.Pb_excl)
        bb.b.Pb_excl[i], bb.b°.Pb_excl[i] = bb.b°.Pb_excl[i], bb.b.Pb_excl[i]
    end
end

function _swap_PP!(bb::BiBlock)
    for i in eachindex(bb.b.PP)
        bb.b.PP[i], bb.b°.PP[i] = bb.b°.PP[i], bb.b.PP[i]
    end
end

"""
    swap_ll!(bb::BiBlock)

Swap `ll` containers between proposal-acceptance pair.
"""
function swap_ll!(bb::BiBlock)
    bb.b.ll, bb.b°.ll = bb.b°.ll, bb.b.ll
end

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

        UTILITY

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

"""
    ll_of_accepted(bb::BiBlock, i)

Return the log-likelihood of the path that was accepted at the `i`th iteration.
"""
function ll_of_accepted(bb::BiBlock, i)
    return (bb.accpt_history[i] ? bb.b°.ll_history[i] : bb.b.ll_history[i])
end

"""
    accpt_rate(bb::BiBlock, range)

Compute the acceptance rate over the `range` of MCMC accept/reject history.
"""
accpt_rate(bb::BiBlock, range) = sum(bb.accpt_history[range])/length(range)

"""
    loglikhd!(b::BiBlock)

Compute the log-likelihood for the accepted block, evaluated at a sampled path
and store the result in an internal field `ll`.
"""
loglikhd!(bb::BiBlock) = loglikhd!(bb.b)

"""
    loglikhd°!(b::BiBlock)

Compute the log-likelihood for the proposal block, evaluated at a sampled path
and store the result in an internal field `ll`.
"""
loglikhd°!(bb::BiBlock) = loglikhd!(bb.b°)

"""
    save_ll!(bb::BiBlock, i::Int)

Commit the current proposal and accepted log-likelihood fields `ll` to history,
at index `i`.
"""
function save_ll!(bb::BiBlock, i::Int)
    save_ll!(bb.b, i)
    save_ll!(bb.b°, i)
end

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

        SETTING ARTIFICIAL OBSERVATIONS

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

"""
    GP.set_obs!(bb::BiBlock)

Freeze an artificial observation at the terminal point of the block. For a
terminal block nothing is done.
"""
function GP.set_obs!(bb::BiBlock) end

function GP.set_obs!(bb::BiBlock{false})
    set_obs!(bb.b.P_last[1], bb.b.XX[end].x[end])
    set_obs!(bb.b°.P_last[1], bb.b.XX[end].x[end])
end

GP.set_obs!(bb::BiBlock{true}) = nothing

"""
    GP.recompute_guiding_term!(bb::BiBlock)

Recompute the guiding terms of both the proposal and the accepted laws.
"""
function GP.recompute_guiding_term!(bb::BiBlock)
    recompute_guiding_term!(bb.b)
    recompute_guiding_term!(bb.b°)
end

"""
    find_W_for_X!(bb::BiBlock)

Find the Wiener process `bb.b.WW` that reconstructs path `bb.b.XX` under the
accepted law `bb.b.PP` (possibly including `bb.b.P_last`).
"""
find_W_for_X!(bb::BiBlock) = find_W_for_X!(bb.b)

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

        SETTING PARAMETERS

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

"""
    is_critical_update(bb::BiBlock, pnames)

Verify whether the update characterized by a list of parameter names stored in
`pnames` is `critical` in a sense of prompting for recomputation of the guiding
term.
"""
GP.is_critical_update(bb::BiBlock, pnames) = GP.is_critical_update(
    bb.b.PP, pnames.θ_local_aux, pnames.θ_local_obs
)

"""
    set_proposal_law!(
        bb::BiBlock,
        θ°,
        pnames,
        critical_change=GP.is_critical_update(bb, pnames);
        skip=0
    )

Set the parameters in `bb.b°.PP` and `bb.b°.P_last` to `θ°` and make sure that
all other parameters are shared with `bb.b.PP` and `bb.b.P_last`. Recompute the
guiding term if needed, and then, compute the proposal trajectory `bb.b°.XX` for
the proposal point `θ°`.
"""
function set_proposal_law!(
        bb::BiBlock,
        θ°,
        pnames,
        critical_change=GP.is_critical_update(bb, pnames);
        skip=0
    )
    critical_change = DD.set_parameters!(bb, θ°, pnames, critical_change)
    critical_change && recompute_guiding_term!(bb.b°)
    recompute_path!(bb.b°, bb.b.WW; skip=skip)
end

"""
    DD.set_parameters!(
        bb::BiBlock,
        θ°,
        pnames,
        critical_change=is_critical_update(bb, pnames)
    )

Set the parameters in `bb.b°.PP` and `bb.b°.P_last` to `θ°` and make sure that
all other parameters are shared with `bb.b.PP` and `bb.b.P_last`.
"""
function DD.set_parameters!(
        bb::BiBlock,
        θ°,
        pnames,
        critical_change=GP.is_critical_update(bb, pnames)
    )
    GP.equalize_obs_params!(bb) && (critical_change = true)
    GP.equalize_law_params!(bb, pnames) && (critical_change = true)
    DD.set_parameters!(bb.b°.PP, θ°, pnames.PP)
    DD.set_parameters!(bb.b°.P_last, θ°, pnames.P_last)
    DD.set_parameters!(bb.b°.P_excl, θ°, pnames.P_excl)
    DD.set_parameters!(bb.b°.Pb_excl, θ°, pnames.Pb_excl)
    critical_change
end

function DD.set_parameters!(PP::AbstractArray{<:GuidProp}, θ°, pnames)
    DD.set_parameters!(PP, θ°, pnames.updt, pnames.updt_aux, pnames.updt_obs)
end

"""
    GP.equalize_obs_params!(bb::BiBlock)

Make sure that in the corresponding pairs of `GuidProp` structs of both `bb.b`
and `bb.b°`, the parameters `θ` of the `obs` fields are the same. If not, then
set the ones in `bb.b°` to be the same as the ones in `bb.b`.
"""
function GP.equalize_obs_params!(bb::BiBlock)
    # equalization of observations that does not matter for this update, but
    # will once the blocks switch...
    GP.equalize_obs_params!(bb.b.P_excl, bb.b°.P_excl)

    # no equalization on P_last or Pb_excl

    # and finally equalization that matters currently, may yield critical change
    GP.equalize_obs_params!(bb.b.PP, bb.b°.PP)
end

"""
    GP.equalize_law_params!(bb::BiBlock, pnames)

Make sure that in the corresponding pairs of `GuidProp` structs of both `bb.b`
and `bb.b°`, the `variable` parameters of each law are the same. If not, then
set the ones in `bb.b°` to be the same as the ones in `bb.b`.
"""
function GP.equalize_law_params!(bb::BiBlock, pnames) end

function GP.equalize_law_params!(bb::BiBlock{true}, pnames)
    _eql_PP!(bb.b.PP, bb.b°.PP, pnames.PP.var, pnames.PP.var_aux)
end

function GP.equalize_law_params!(bb::BiBlock{false}, pnames)
    critical_change = _eql_PP!(
        bb.b.PP, bb.b°.PP, pnames.PP.var, pnames.PP.var_aux
    )

    _eql_PP!(
        bb.b.P_last, bb.b°.P_last, pnames.P_last.var, pnames.P_last.var_aux
    ) && (critical_change = true)

    # this is just not needed now, but will be for the next blocks
    _eql_PP!(
        bb.b.P_excl, bb.b°.P_excl, pnames.P_excl.var, pnames.P_excl.var_aux
    )
    _eql_PP!(
        bb.b.Pb_excl, bb.b°.Pb_excl, pnames.Pb_excl.var, pnames.Pb_excl.var_aux
    )

    critical_change
end

"""
    _eql_PP!(PP, PP°, var_p_names, var_p_aux_names)

Go through all `GuidProp` structs in  `PP` and `PP°` and make sure that the
parameters listed in `var_p_names` and `var_p_aux_names` agree. Equalize them if
they do not. `var_p_names` should list all parameter names from the target law
and `var_p_aux_names` should list all parameters from the auxiliary law.
"""
function _eql_PP!(PP, PP°, var_p_names, var_p_aux_names)
    if !DD.same_entries(PP, PP°, var_p_names)
        GP.equalize_law_params!(
            PP, PP°, var_p_names, var_p_aux_names
        ) && return true
    end
    false
end


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

using Markdown
using InteractiveUtils

# ╔═╡ c7e41465-41ad-4f11-baec-bfaae9edbf59
begin
	using PlutoUI
	PlutoUI.TableOfContents(aside=true)
end

# ╔═╡ d0233ece-1f69-11eb-3da1-43b8ddc9fbe1
begin
	# Necessary packages
	using LinearAlgebra
	using Random
	using Statistics
	using Clustering
	using Plots
end

# ╔═╡ 00d72ad7-ea80-4954-8c75-590d4e9ccdc4
md"""
# K-means Algorithm

Data clustering is one of the main mathematical applications variety of algorithms have been developed to tackle the problem. K-means is one of the basic algorithms for data clustering.

__Prerequisites__

The reader should be familiar with basic linear algebra.

__Competences__

The reader should be able to recognise applications where K-means algorithm can be efficiently used and use it.

__Credits.__ The notebook was initially derived from M.Sc. Thesis of Ivančica Mirošević.
"""

# ╔═╡ 235e1901-eed3-4934-9214-ffdd82cf21d6
md"""
## Definitions

__Data clustering problem__ is the following: partition the given set of $m$ objects of the same type into $k$ subsets according to some criterion. Additional request may be to find the optimal $k$.

__K-means clustering problem__ is the following: partition the set  $X=\{x_{1},x_{2},\cdots ,x_{m}\}$ , where $x_{i}\in\mathbb{R}^{n}$, into $k$ _clusters_ $\pi=\{C_{1},C_{2},...,C_{k}\}$ such that

$$
J(\pi)=\sum_{i=1}^{k}\sum_{x\in
C_{i}}\| x-c_{i}\|_{2}^{2} \to \min$$

over all possible partitions. Here $c_{i}=\displaystyle\frac{1}{|C_{i}|}\sum_{x\in C_{i}} x$ is the mean of points in $C_i$ and $|C_i|$ is the cardinality of $C_i$.
"""

# ╔═╡ 26640b90-1f69-11eb-029e-4dc3a825f639
md"
## K-means clustering algorithm

1. __Initialization__: Choose initial set of $k$ means $\{c_1,\ldots,c_k\}$ (for example, by choosing randomly $k$ points from $X$).
2. __Assignment step__: Assign each point $x$ to one nearest mean $c_i$.
3. __Update step__: Compute the new means.
4. __Iteration__: Repeat Steps 2 and 3 until the assignment no longer changes.
"

# ╔═╡ 48914840-1f69-11eb-30a9-3943767bc261
md"
## First Variation clustering algorithm

A __first variation__ of a partition $\pi=\{C_1,\ldots,C_k\}$ is
a partition $\pi^{\prime}=\{C_{1}^{\prime},\cdots ,C_{k}^{\prime }\}$
obtained by moving a single point $x$ from a cluster  $C_{i}$ to a cluster $C_{j}$. Notice that $\pi$ is a first variation of itself.

A __next partition__ of the partition $\pi$ is a partition
$\mathop{\mathrm{next}}(\pi)=\mathop{\mathrm{arg min}}\limits_{\pi^{\prime}} J(\pi^{\prime})$.

We have the following algorithm:

1. Choose initial partition $\pi$.
2. Compute $\mathop{\mathrm{next}}(\pi)$
3. If $J(\mathop{\mathrm{next}}(\pi))<J(\pi)$, set $\pi=\mathop{\mathrm{next}}(\pi)$ and go to Step 2
4. Stop.
"

# ╔═╡ 2b76953c-edf2-41cf-9d75-7efb2d6a58f1
md"""
## Facts

1. The k-means clustering problem is NP-hard.

2. In the k-means algorithm, $J(\pi)$ decreases in every iteration.

3. K-means algorithm can converge to a local minimum.

4. Each iteration of the k-means algorithm requires $O(mnk)$ operations.

4. K-means algorithm is implemented in the function `kmeans()` in the package [Clustering.jl](https://github.com/JuliaStats/Clustering.jl).

5.  $J(\pi)=\mathop{\mathrm{trace}}(S_W)$, where

$$
S_{W}=\sum\limits_{i=1}^k\sum\limits_{x\in C_{i}}
(x-c_i)(x-c_{i})^{T}
=\sum_{i=1}^k\frac{1}{2|C_{i}|}\sum_{x\in C_{i}}\sum_{y \in C_{i}}
(x-y)(x-y)^{T}.$$

Let $c$ denote the mean of $X$. Then $S_W=S_T-S_B$, where

$$
\begin{aligned}
S_{T}&=\sum_{x\in X}(x-c)(x-c)^{T} =
\frac{1}{2m}\sum_{i=1}^m\sum_{j=1}^m
(x_{i}-x_{j})(x_{i}-x_{j})^{T}, \\
S_{B}&=\sum_{i=1}^k|C_{i}|(c_{i}-c)(c_{i}-c)^{T} =
\frac{1}{2m}\sum_{i=1}^k\sum_{j=1}^k|C_{i}||C_{j}|
(c_{i}-c_{j})(c_{i}-c_{j})^{T}.
\end{aligned}$$

6. In order to try to avoid convergence to local minima, the k-means algorithm can be enhanced with first variation by adding the following steps:
    1. Compute $\mathop{\mathrm{next}}(\pi)$.
    2. If $J(\mathop{\mathrm{next}}(\pi))<J(\pi)$, set $\pi=\mathop{\mathrm{next}}(\pi)$ and go to Step 2.

"""

# ╔═╡ 86d55833-ee5a-4f73-b0bf-e8b6e6d65617
function myKmeans(X::Vector{T}, k::Int) where T
    # X is Array of Arrays
    m,n=length(X),length(X[1])
    C=Vector{Int}(undef,m)
    # Choose random k means among X
    c=X[randperm(m)[1:k]]
    # This is just to start the while loop
    cnew=copy(c)
    cnew[1]=cnew[1].+(1.0,1.0)
    # Loop
    iterations=0
    while cnew!=c
        iterations+=1
        cnew=copy(c)
        # Assignment
        for i=1:m
            C[i]=findmin([norm(X[i].-c[j]) for j=1:k])[2]
        end
        # Update
        for j=1:k
          c[j]=(mean([x[1] for x in X[C.==j]]),mean([x[2] for x in X[C.==j]]))
        end
    end
    C,c,iterations
end

# ╔═╡ e15fa4cf-930a-4065-b839-c2ec944afcfa
md"""

## Examples

### Random clusters

We generate $k$ random clusters around points with integer coordinates.
"""

# ╔═╡ bacdc83d-1acd-46df-8d06-3e0279978a33
begin
	# Generate points as Tuple()
	k=5
	Random.seed!(1235)
	centers= [Tuple(rand(-5:5,2)) for i=1:k]
	# Number of points in cluster
	sizes=rand(10:50,k)
	csizes=cumsum(sizes)
	# X is array of arrays
	X=Vector{Tuple{Float64,Float64}}(undef,sum(sizes))
	X[1:csizes[1]]=[centers[1].+Tuple((rand(2).-0.5)/2) for i=1:sizes[1]]
	for j=2:k
		X[csizes[j-1]+1:csizes[j]]=[centers[j].+Tuple((rand(2).-0.5)/2) for i=1:sizes[j]]
	end
	centers, sizes
end

# ╔═╡ 99aba610-1f75-11eb-1b5e-2f2ee211ed7f
X

# ╔═╡ fceab134-0921-459b-8d91-952a44abb07b
begin
	# Plot
	scatter(X,label="Points",title="Points")
	scatter!(centers,markershape = :hexagon, ms = 6,label="Centers")
end

# ╔═╡ f559b870-608a-426e-a23b-b66ef6d4b47f
# Plot the solution
function plotKmeansresult(C::Vector,c::Vector,X::Vector)
    scatter()
    # Clusters
    for j=1:k
        scatter!(X[findall(C.==j)],label="Cluster $j")
    end
    # Means
    scatter!(c,markershape=:hexagon,ms=6,color=:red,label="Centers")
	plot!(title="Computed clusters")
end

# ╔═╡ ba84fe97-4dd0-47d3-b980-f5e5a8de6fa3
X

# ╔═╡ 310a4957-cc58-424b-845c-3ebaf7db77f3
md"""
__What happens?__

We see that the algorithm, although simple, for this example
converges to a local minimum.

Let us try the function `kmeans()` from the package `Clustering.jl`.
The inputs to `kmeans()` are:

* a matrix whose columns are points,
* number of clusters we are looking for, and,
* optionally, the method to compute initial means.

If we choose `init=:rand`, the results are similar. If we choose
`init=:kmpp`, which is the default, the results are better, but convergence to a local minimum is still possible.

__Run the clustering several times!__

```
seeding_algorithm(s::Symbol) =
    s == :rand ? RandSeedAlg() :
    s == :kmpp ? KmppAlg() :
    s == :kmcen ? KmCentralityAlg() :
    error("Unknown seeding algorithm $s")
```
"""

# ╔═╡ f9af8a75-48fd-4e82-b43e-61867bead6f9
methods(kmeans)

# ╔═╡ 7e869699-8eb1-4e38-905d-a9448d037841
begin
	Xₘ=transpose([[x[1] for x in X] [x[2] for x in X]])
	output=kmeans(Matrix(Xₘ),k,init=:kmpp)
end

# ╔═╡ 4f8152ee-ce86-4df6-81d2-be2df33df980
fieldnames(KmeansResult)

# ╔═╡ 95fb2e66-e811-47f4-a772-8270d0c5d83b
output.centers

# ╔═╡ 26b7bac0-1f72-11eb-16fc-c1d74920c462
output.assignments

# ╔═╡ 7bb0e9ce-72dd-4623-ab3d-a6d18fb9dd15
# We need to modify the plotting function
function plotKmeansresult(out::KmeansResult,X::AbstractArray)
    k=size(out.centers,2)
    scatter(aspect_ratio=1)
    # Clusters
    for j=1:k
        scatter!(X[1,findall(out.assignments.==j)], X[2,findall(out.assignments.==j)], label="Cluster $j")
    end
    # Means
    scatter!(out.centers[1,:], out.centers[2,:], markershape=:hexagon,ms=6,color=:red,label="Centers")
end

# ╔═╡ 8e5b89ce-e6b5-4473-bd98-a2285331418c
begin
	C,c,iterations=myKmeans(X,k)
	plotKmeansresult(C,c,X)
end

# ╔═╡ 406bb1b2-a8bb-4cd2-95c3-3975f99e3f32
begin
	out=kmeans(Xₘ,k,init=:kmpp)
	plotKmeansresult(out,Xₘ)
end

# ╔═╡ 2445e75a-77f9-4c27-b400-44d8b49fc03b
md"""
### Concentric rings

The k-means algorithm works well if clusters can be separated by hyperplanes. In this example it is not the case.
"""

# ╔═╡ 21bdc615-76e8-46cc-ac86-310d28f6cd25
begin
	# Number of rings, try also k=3
	k₁=2
	# Center
	Random.seed!(5361)
	center=[rand(-5:5);rand(-5:5)]
	# Radii
	radii=randperm(10)[1:k₁]
	# Number of points in circles
	sizes₁=rand(1000:2000,k₁)
	center,radii,sizes₁
end

# ╔═╡ 9a4060b7-ba5a-401f-90f2-b9780ed93ddb
begin
	# Generate points
	X₁=Array{Float64}(undef,2,sum(sizes₁))
	csizes₁=cumsum(sizes₁)
	# Random angles
	ϕ=2*π*rand(sum(sizes₁))
	for i=1:csizes₁[1]
		X₁[:,i]=center+radii[1]*[cos(ϕ[i]);sin(ϕ[i])] + (rand(2).-0.5)/50
	end
	for j=2:k₁
		for i=csizes₁[j-1]+1:csizes₁[j]
			X₁[:,i]=center+radii[j]*[cos(ϕ[i]);sin(ϕ[i])] + (rand(2).-0.5)/50
		end
	end
	scatter(X₁[1,:],X₁[2,:],title="Concentric Rings", aspect_ratio=1,label="Points")
end

# ╔═╡ 0be24dd5-f3b0-4daa-aad8-5dbcdfca9313
begin
	out₁=kmeans(X₁,k₁,init=:rand)
	plotKmeansresult(out₁,X₁)
end

# ╔═╡ 1637ce15-f55a-4d66-8a00-45a555bf5d61

# ╔═╡ Cell order:
# ╟─c7e41465-41ad-4f11-baec-bfaae9edbf59
# ╟─00d72ad7-ea80-4954-8c75-590d4e9ccdc4
# ╟─235e1901-eed3-4934-9214-ffdd82cf21d6
# ╟─26640b90-1f69-11eb-029e-4dc3a825f639
# ╟─48914840-1f69-11eb-30a9-3943767bc261
# ╟─2b76953c-edf2-41cf-9d75-7efb2d6a58f1
# ╠═d0233ece-1f69-11eb-3da1-43b8ddc9fbe1
# ╠═86d55833-ee5a-4f73-b0bf-e8b6e6d65617
# ╟─e15fa4cf-930a-4065-b839-c2ec944afcfa
# ╠═bacdc83d-1acd-46df-8d06-3e0279978a33
# ╠═99aba610-1f75-11eb-1b5e-2f2ee211ed7f
# ╠═fceab134-0921-459b-8d91-952a44abb07b
# ╠═f559b870-608a-426e-a23b-b66ef6d4b47f
# ╠═8e5b89ce-e6b5-4473-bd98-a2285331418c
# ╠═ba84fe97-4dd0-47d3-b980-f5e5a8de6fa3
# ╟─310a4957-cc58-424b-845c-3ebaf7db77f3
# ╠═f9af8a75-48fd-4e82-b43e-61867bead6f9
# ╠═7e869699-8eb1-4e38-905d-a9448d037841
# ╠═4f8152ee-ce86-4df6-81d2-be2df33df980
# ╠═95fb2e66-e811-47f4-a772-8270d0c5d83b
# ╠═26b7bac0-1f72-11eb-16fc-c1d74920c462
# ╠═7bb0e9ce-72dd-4623-ab3d-a6d18fb9dd15
# ╠═406bb1b2-a8bb-4cd2-95c3-3975f99e3f32
# ╟─2445e75a-77f9-4c27-b400-44d8b49fc03b
# ╠═21bdc615-76e8-46cc-ac86-310d28f6cd25
# ╠═9a4060b7-ba5a-401f-90f2-b9780ed93ddb
# ╠═0be24dd5-f3b0-4daa-aad8-5dbcdfca9313
# ╠═1637ce15-f55a-4d66-8a00-45a555bf5d61

ext: jl

max_stars_repo_path: src/Module C - Applications/L09 K-means Algorithm.jl

max_stars_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_issues_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course

max_forks_repo_name: ivanslapnicar/GIAN-Applied-NLA-Course


text:
#(When threading works, not yet):
#       To get the benefit of threading start julia with multiple threads:
#       JULIA_NUM_THREADS=x julia (where you choose x)

#A multiprocess write method for saving large image sequences in NRRD format
#This function is useful when retrieving a stack from img takes a lot of time due to lazy computation
function multiproc_write(fname::AbstractString, img::AbstractArray{T}, out_type=T) where {T}
    mmapa = prep_nrrd_write(fname, img, out_type)
    multiproc_write!(mmapa, img)
end

function prep_nrrd_write(fname::AbstractString, img::AbstractArray{T,N}, out_type=T) where {T,N}
    bname, ext = splitext(fname)
    if !isempty(ext) && !in(ext, [".nhdr", ".nrrd"])
        error("Must write in NRRD format (.nhdr or .nrrd extension)")
    end
    s = open(bname*".raw", "w+")
    mmapa = Mmap.mmap(s, Array{out_type, N}, size(img); shared=true)
    header = NRRD.headerinfo(out_type, axes(img))
    header["datafile"] = bname*".raw"
    open(fname, "w") do io
        write(io, magic(format"NRRD"))
        NRRD.write_header(io, "0004", header)
    end
    return mmapa
end

function multiproc_write!(mmapA::AbstractArray{T,3}, img::AbstractArray{T2,3}) where {T, T2}
    @assert size(mmapA) == size(img)
    mmapA[:,:,:] = img[:,:,:]
end

#Multithreading would be nice, but currently this doesn't seem to work well, perhaps because IO operations and external C library calls are restricted to one thread (see julia parallel docs as off march 2018)
function multiproc_write!(mmapA::AbstractArray{T,4}, img::AbstractArray{T2,4}) where {T, T2}
    @assert size(mmapA) == size(img)
    np = nworkers()
    if np == 0
        error("No workers available.  Add them with addprocs() (also make sure code and the input image are loaded on the worker processes")
    end
    print("Defining image and output array on other processes...\n")
    if np > 1
        @eval @everywhere img = $(img)
    end
    #nthreads = Base.Threads.nthreads()
    #print("Writing image with $(nthreads) threads\n")
    print("Writing image with $(np) procs\n")
    cur_i = 1
    while cur_i<=size(img,4)
        i_start = cur_i
        cur_i = min(size(img,4), i_start + np-1)
        idxs = [i_start:cur_i...]
        rrefs = []
        ws = workers()
#    Threads.@threads for i=1:size(img,4)
        for i=1:length(idxs)
            push!(rrefs, remotecall(getindex, ws[i], img, :,:,:,idxs[i]))
        end
        for i=1:length(idxs)
            mmapA[:,:,:,idxs[i]] = fetch(rrefs[i])
        end
        print("Wrote $cur_i stacks out of $(size(img,4)) so far.\n")
        cur_i+=1
    end
    print("Finished\n")
end


text:
function ellipse_traj(x_start, x_goal, z, T)
	dist = x_goal - x_start
	a = 0.5 * dist
	b = z
	z̄ = 0.0
	# x = range(x_start, stop = x_goal, length = T)
	x = circular_projection_range(x_start, stop = x_goal, length = T)
	z = sqrt.(max.(0.0, (b^2) * (1.0 .- ((x .- (x_start + a)).^2.0) / (a^2.0))))
	return x, z
end

function circular_projection_range(start; stop=1.0, length=10)
	dist = stop - start
	θr = range(π, stop=0, length=length)
	r = start .+ dist * ((1 .+ cos.(θr))./2)
	return r
end

function initial_configuration(θ1, θ2, θ3, θ4;
	l_torso = 0.367,
    l_thigh = 0.2,
    l_leg = 0.2,
    w_torso = 0.267)

	q_body = zeros(6)
	q_body[1] = 0.0

	# position
	q_body[3] = l_thigh * cos(θ1) + l_leg * cos(θ2 - θ1)

	# orientation
	mrp = MRP(RotZ(0.0))
	q_body[4:6] = [mrp.x; mrp.y; mrp.z]

	# feet positions (in body frame)
	q_f1 = [0.5 * l_torso - l_thigh * sin(θ1) + l_leg * sin(θ2 - θ1); -0.5 * w_torso; 0.0]
	q_f2 = [0.5 * l_torso - l_thigh * sin(θ3) + l_leg * sin(θ4 - θ3); 0.5 * w_torso; 0.0]
	q_f3 = [-0.5 * l_torso - l_thigh * sin(θ3) + l_leg * sin(θ4 - θ3); -0.5 * w_torso; 0.0]
	q_f4 = [-0.5 * l_torso - l_thigh * sin(θ1) + l_leg * sin(θ2 - θ1); 0.5 * w_torso; 0.0]

	q_body[1] = 0.5 * (0.5 * (q_f1[1] + q_f2[1]) + 0.5 * (q_f3[1] + q_f4[1]))

	return q_body, q_f1, q_f2, q_f3, q_f4
end

function generate_gait(θ1, θ2, θ3, θ4, T, Tm, h;
    l_torso = 0.367,
    l_thigh = 0.2,
    l_leg = 0.2,
    w_torso = 0.267)

	q_body, q_f1, q_f2, q_f3, q_f4 = initial_configuration(
        θ1, θ2, θ3, θ4,
        l_torso = l_torso,
        l_thigh = l_thigh,
        l_leg = l_leg,
        w_torso = w_torso)

	# visualize!(vis, model, [q_body], [q_f1], [q_f2], [q_f3], [q_f4])

	# feet positions
	pf1 = copy(q_f1)
	pf2 = copy(q_f2)

	pf3 = copy(q_f3)
	pf4 = copy(q_f4)

	strd = 2 * (pf1 - pf2)[1]

	# q_shift1 = zeros(model.nq)
	# q_shift1[1] = 0.5 * strd
	# q_shift1[10] = strd
	# q_shift1[13] = strd

	body_shift = [0.5 * strd; 0.0; 0.0; 0.0; 0.0; 0.0]
	foot_shift = [1.0 * strd; 0.0; 0.0]

	qM_body = q_body + body_shift
	qM_f1 = copy(q_f1)
	qM_f2 = q_f2 + foot_shift
	qM_f3 = q_f3 + foot_shift
	qM_f4 = copy(q_f4)

	qT_body = qM_body + body_shift
	qT_f1 = copy(qM_f1) + foot_shift
	qT_f2 = copy(qM_f2)
	qT_f3 = copy(qM_f3)
	qT_f4 = copy(qM_f4) + foot_shift

	q_body_ref = linear_interpolation(q_body, qT_body, T)#[q_body, linear_interpolation(q_body, qM_body, Tm - 1)..., linear_interpolation(qM_body, qT_body, Tm - 1)...]
	q_f1_ref = [q_f1, linear_interpolation(q_f1, qM_f1, Tm - 1)..., linear_interpolation(qM_f1, qT_f1, Tm - 1)...]
	q_f2_ref = [q_f2, linear_interpolation(q_f2, qM_f2, Tm - 1)..., linear_interpolation(qM_f2, qT_f2, Tm - 1)...]
	q_f3_ref = [q_f3, linear_interpolation(q_f3, qM_f3, Tm - 1)..., linear_interpolation(qM_f3, qT_f3, Tm - 1)...]
	q_f4_ref = [q_f4, linear_interpolation(q_f4, qM_f4, Tm - 1)..., linear_interpolation(qM_f4, qT_f4, Tm - 1)...]

	# visualize!(vis, model, q_body_ref, q_f1_ref, q_f2_ref, q_f3_ref, q_f4_ref,
	#     Δt = h)

	T_fix = 5
	zh = 0.05
	xf1_el, zf1_el = ellipse_traj(pf1[1], pf1[1] + strd, zh, Tm - T_fix)
	xf1 = [[pf1[1] for t = 1:Tm + T_fix]..., xf1_el[2:end]...]
	zf1 = [[pf1[3] for t = 1:Tm + T_fix]..., zf1_el[2:end]...]
	pf1_ref = [[xf1[t]; pf1[2];  zf1[t]] for t = 1:T]
	f1_contact = [t < (Tm + T_fix) ? 1 : (t == T ? 1 : 0) for t = 1:T]

	xf4_el, zf4_el = ellipse_traj(pf4[1], pf4[1] + strd, zh, Tm - T_fix)
	xf4 = [[pf4[1] for t = 1:Tm + T_fix]..., xf4_el[2:end]...]
	zf4 = [[pf4[3] for t = 1:Tm + T_fix]..., zf4_el[2:end]...]
	pf4_ref = [[xf4[t]; pf4[2]; zf4[t]] for t = 1:T]
	f4_contact = [t < (Tm + T_fix) ? 1 : (t == T ? 1 : 0) for t = 1:T]

	xf2_el, zf2_el = ellipse_traj(pf2[1], pf2[1] + strd, zh, Tm - T_fix)
	xf2 = [[xf2_el[1] for t = 1:T_fix]..., xf2_el..., [xf2_el[end] for t = 1:Tm-1 + T_fix]...]
	zf2 = [[zf2_el[1] for t = 1:T_fix]..., zf2_el..., [zf2_el[end] for t = 1:Tm-1 + T_fix]...]
	pf2_ref = [[xf2[t]; pf2[2]; zf2[t]] for t = 1:T]
	f2_contact = [t <= T_fix ? 1 : (t > Tm ? 1 : 0) for t = 1:T]

	xf3_el, zf3_el = ellipse_traj(pf3[1], pf3[1] + strd, zh, Tm - T_fix)
	xf3 = [[xf3_el[1] for t = 1:T_fix]..., xf3_el..., [xf3_el[end] for t = 1:Tm-1]...]
	zf3 = [[zf3_el[1] for t = 1:T_fix]..., zf3_el..., [zf3_el[end] for t = 1:Tm-1]...]
	pf3_ref = [[xf3[t]; pf3[2]; zf3[t]] for t = 1:T]
	f3_contact = [t <= T_fix ? 1 : (t > Tm ? 1 : 0) for t = 1:T]

	tr = range(0, stop = tf, length = T)
	# plot(tr, hcat(pf1_ref...)', labels = "")
	# plot!(tr, hcat(pf4_ref...)', labels = "")
	#
	# plot!(tr, hcat(pf2_ref...)', labels = "")
	# plot!(tr, hcat(pf3_ref...)', labels = "")

	# visualize!(vis, model, q_body_ref,
	#     pf1_ref, pf2_ref, pf3_ref, pf4_ref,
	#     Δt = h)

	# p_fr = plot(tr, f1_contact, linetype = :steppost, color = :black, label = "FR", xaxis = false, x_ticks = false, yaxis = false, yticks = false)
	# p_rl = plot(tr, f4_contact, linetype = :steppost, color = :black, label = "RL", xaxis = false, x_ticks = false, yaxis = false, yticks = false)
	#
	# p_fl = plot(tr, f2_contact, linetype = :steppost, color = :black, label = "FL", xaxis = false, x_ticks = false, yaxis = false, yticks = false)
	# p_rr = plot(tr, f3_contact, linetype = :steppost, color = :black, label = "RR", yaxis = false, yticks = false, xaxis = "time (s)")
	#
	# plot(p_fr, p_rl, p_fl, p_rr, layout = (4, 1))

	contact_modes = [[f1_contact[t];
	                  f2_contact[t];
	                  f3_contact[t];
	                  f4_contact[t]] for t = 1:T]

	contact_modes_2 = [contact_modes..., contact_modes[2:end]...]

	plot(hcat(contact_modes_2...)', linetype = :steppost, label = "")

	T2 = 2 * T - 1
	tf2 = 2 * tf
	tr2 = range(0, stop = tf2, length = T2)
	pf1_ref_2 = [pf1_ref..., [p + foot_shift for p in pf1_ref[2:end]]...]
	pf2_ref_2 = [pf2_ref..., [p + foot_shift for p in pf2_ref[2:end]]...]
	pf3_ref_2 = [pf3_ref..., [p + foot_shift for p in pf3_ref[2:end]]...]
	pf4_ref_2 = [pf4_ref..., [p + foot_shift for p in pf4_ref[2:end]]...]
	qT_body_2 = q_body + 4 * body_shift
	q_body_ref_2 = linear_interpolation(q_body, qT_body_2, T2)#[q_body, linear_interpolation(q_body, qM_body, Tm - 1)..., linear_interpolation(qM_body, qT_body, Tm - 1)...]

	# plot(tr2, hcat(pf1_ref_2...)', labels = "")
	# plot!(tr2, hcat(pf2_ref_2...)', labels = "")
	# plot!(tr2, hcat(pf3_ref_2...)', labels = "")
	# plot!(tr2, hcat(pf4_ref_2...)', labels = "")

	return q_body_ref_2, pf1_ref_2, pf2_ref_2, pf3_ref_2, pf4_ref_2, contact_modes_2
	# return q_body_ref, pf1_ref, pf2_ref, pf3_ref, pf4_ref, contact_modes

end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using GridInterpolations
using LocalFunctionApproximation
using LocalApproximationValueIteration
using POMDPs
using POMDPModels
using POMDPModelTools
using StaticArrays
using JLD2, FileIO
using Random
using Logging
using Distributions
using CMSSPs

# scale_file = ARGS[1]
# simtime_file = ARGS[2]
# cost_file = ARGS[3]
# policy_name = ARGS[4]
# poly_or_exp = ARGS[5]

scale_file = "./paramsets/scale-1.toml"
simtime_file = "./paramsets/simtime-1.toml"
cost_file = "./paramsets/cost-1.toml"
policy_name = "dreamr-uf-params111"
poly_or_exp = "poly"
num_iters = 50

rng = MersenneTwister(10)

params = parse_params(scale_file=scale_file, simtime_file=simtime_file, cost_file=cost_file)

uf_mdp = get_uf_mdp(MultiRotorUAVState, MultiRotorUAVAction, params)

if poly_or_exp == "poly"

    xy_spacing = polyspace_symmetric(params.scale_params.XY_LIM, params.scale_params.XY_AXISVALS)
    xydot_spacing = polyspace_symmetric(params.scale_params.XYDOT_LIM, params.scale_params.XYDOT_AXISVALS)

elseif poly_or_exp == "exp"

    xy_spacing = log2space_symmetric(params.scale_params.XY_LIM, params.scale_params.XY_AXISVALS)
    xydot_spacing = log2space_symmetric(params.scale_params.XYDOT_LIM, params.scale_params.XYDOT_AXISVALS)

end

state_grid = RectangleGrid(xy_spacing, xy_spacing, xydot_spacing, xydot_spacing)
@show length(vertices(state_grid))

lfa = LocalGIFunctionApproximator(state_grid)

uf_policy = infinite_horizon_VI_localapprox(uf_mdp, lfa, num_iters, true, params.scale_params.MC_GENERATIVE_NUMSAMPLES, rng)

policy_fn = string(policy_name,".jld2")
save_localapproxvi_policy_to_jld2(policy_fn, uf_policy, uf_mdp)


text:
using PyPlot, Distributions

using TransitionalMCMC

# Prior Bounds
lb, ub  = -5, 5

# Prior Density and sampler
logprior(x) = logpdf(Uniform(lb, ub), x[1]) + logpdf(Uniform(lb, ub), x[2])
priorRnd(Nsamples) = rand(Uniform(lb, ub), Nsamples, 2)

# Log Likelihood
function logLik(x)
    return -1 * ((x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2)
end

Nsamples = 200

samps, acc = tmcmc(logLik, logprior, priorRnd, Nsamples, 5, 2)

plt.scatter(samps[:,1], samps[:,2])


text:
module GalacticNLopt

using Reexport, GalacticOptim, GalacticOptim.SciMLBase
@reexport using NLopt

(f::NLopt.Algorithm)() = f

function __map_optimizer_args!(prob::OptimizationProblem, opt::NLopt.Opt;
    callback=nothing,
    maxiters::Union{Number,Nothing}=nothing,
    maxtime::Union{Number,Nothing}=nothing,
    abstol::Union{Number,Nothing}=nothing,
    reltol::Union{Number,Nothing}=nothing,
    local_method::Union{NLopt.Algorithm,NLopt.Opt,Nothing}=nothing,
    local_maxiters::Union{Number,Nothing}=nothing,
    local_maxtime::Union{Number,Nothing}=nothing,
    local_options::Union{NamedTuple,Nothing}=nothing,
    kwargs...)

    if local_method !== nothing
        if isa(local_method, NLopt.Opt)
            if ndims(local_method) != length(prob.u0)
                error("Passed local NLopt.Opt optimization dimension does not match OptimizationProblem dimension.")
            end
            local_meth = local_method
        else
            local_meth = NLopt.Opt(local_method, length(prob.u0))
        end

        if !isnothing(local_options)
            for j in Dict(pairs(local_options))
                eval(Meta.parse("NLopt." * string(j.first) * "!"))(local_meth, j.second)
            end
        end

        if !(isnothing(local_maxiters))
            NLopt.maxeval!(local_meth, local_maxiters)
        end

        if !(isnothing(local_maxtime))
            NLopt.maxtime!(local_meth, local_maxtime)
        end

        NLopt.local_optimizer!(opt, local_meth)
    end

    # add optimiser options from kwargs
    for j in kwargs
        eval(Meta.parse("NLopt." * string(j.first) * "!"))(opt, j.second)
    end

    if prob.ub !== nothing
        NLopt.upper_bounds!(opt, prob.ub)
    end

    if prob.lb !== nothing
        NLopt.lower_bounds!(opt, prob.lb)
    end

    if !(isnothing(maxiters))
        NLopt.maxeval!(opt, maxiters)
    end

    if !(isnothing(maxtime))
        NLopt.maxtime!(opt, maxtime)
    end

    if !isnothing(abstol)
        NLopt.ftol_abs!(opt, abstol)
    end
    if !isnothing(reltol)
        NLopt.ftol_rel!(opt, reltol)
    end

    return nothing
end

function SciMLBase.__solve(prob::OptimizationProblem, opt::Union{NLopt.Algorithm,NLopt.Opt};
    maxiters::Union{Number,Nothing}=nothing,
    maxtime::Union{Number,Nothing}=nothing,
    local_method::Union{NLopt.Algorithm,NLopt.Opt,Nothing}=nothing,
    local_maxiters::Union{Number,Nothing}=nothing,
    local_maxtime::Union{Number,Nothing}=nothing,
    local_options::Union{NamedTuple,Nothing}=nothing,
    abstol::Union{Number,Nothing}=nothing,
    reltol::Union{Number,Nothing}=nothing,
    progress=false,
    callback=(args...) -> (false),
    kwargs...)
    local x

    maxiters = GalacticOptim._check_and_convert_maxiters(maxiters)
    maxtime = GalacticOptim._check_and_convert_maxtime(maxtime)
    local_maxiters = GalacticOptim._check_and_convert_maxiters(local_maxiters)
    local_maxtime = GalacticOptim._check_and_convert_maxtime(local_maxtime)

    f = GalacticOptim.instantiate_function(prob.f, prob.u0, prob.f.adtype, prob.p)

    _loss = function (θ)
        x = f.f(θ, prob.p)
        callback(θ, x...)
        return x[1]
    end

    fg! = function (θ, G)
        if length(G) > 0
            f.grad(G, θ)
        end

        return _loss(θ)
    end

    if isa(opt, NLopt.Opt)
        if ndims(opt) != length(prob.u0)
            error("Passed NLopt.Opt optimization dimension does not match OptimizationProblem dimension.")
        end
        opt_setup = opt
    else
        opt_setup = NLopt.Opt(opt, length(prob.u0))
    end

    prob.sense === GalacticOptim.MaxSense ? NLopt.max_objective!(opt_setup, fg!) : NLopt.min_objective!(opt_setup, fg!)

    __map_optimizer_args!(prob, opt_setup, maxiters=maxiters, maxtime=maxtime, abstol=abstol, reltol=reltol, local_method=local_method, local_maxiters=local_maxiters, local_options=local_options; kwargs...)

    t0 = time()
    (minf, minx, ret) = NLopt.optimize(opt_setup, prob.u0)
    t1 = time()

    SciMLBase.build_solution(prob, opt, minx, minf; original=opt_setup, retcode=ret)
end

end


text:
using UAParser, YAML, Test

# a helper macro for creating tests
macro testparseval(obj::Symbol, valname::String, outputname::Symbol)
    valsymb = Symbol(valname)
    esc(quote
        if $obj[$valname] == nothing
            ismissing($outputname.$valsymb)
        else
            $obj[$valname] == $outputname.$valsymb
        end
    end)
end

#Since parser is known not to be 100% compatible with source Python parser (see README), testing accuracy
#Accuracy stats derived from test files for v0.6 of package
#Any future changes to parsing code will be evaluated against this standard

@testset "parse_device" begin
    #Test 1: Validation of parsedevice
    test_device = YAML.load(open(joinpath(dirname(@__FILE__), "data", "test_device.yaml")))

    pass = 0
    fail = 0
    for test_case in test_device["test_cases"]
        old_fail = fail
        str = test_case["user_agent_string"]
        dev = parsedevice(str)
        test_case["family"] == dev.family ? pass += 1 : fail += 1
        fail == old_fail || @info """When parsing device for "$str", expect $test_case, but got $dev"""
    end

    @test pass/(pass + fail) >= 0.999

end

@testset "parse_os" begin
    #Test 2: Validation of parseos
    test_os = YAML.load(open(joinpath(dirname(@__FILE__), "data", "test_os.yaml")))

    pass = 0
    fail = 0

    for test_case in test_os["test_cases"]
        old_fail = fail
        str = test_case["user_agent_string"]
        os = parseos(str)
        test_case["family"] == os.family ? pass += 1 : fail += 1
        (@testparseval test_case "major" os) ? pass += 1 : fail += 1
        (@testparseval test_case "minor" os) ? pass += 1 : fail += 1
        (@testparseval test_case "patch_minor" os) ? pass += 1 : fail += 1
        fail == old_fail || @info """When parsing os for "$str", expect $test_case, but got $os"""
    end

    @test pass/(pass + fail) >= 0.999

end

@testset "parse_ua" begin
    #Test 4: Validation of parseuseragent
    test_ua = YAML.load(open(joinpath(dirname(@__FILE__), "data", "test_ua.yaml")))

    pass = 0
    fail = 0

    for test_case in test_ua["test_cases"]
        old_fail = fail
        str = test_case["user_agent_string"]
        ua = parseuseragent(str)
        test_case["family"] == ua.family ? pass += 1 : fail += 1
        (@testparseval test_case "major" ua) ? pass += 1 : fail += 1
        (@testparseval test_case "minor" ua) ? pass += 1 : fail += 1
        (@testparseval test_case "patch" ua) ? pass += 1 : fail += 1
        fail == old_fail || @info """When parsing ua for "$str", expect $test_case, but got $ua"""
    end

    @test pass/(pass + fail) >= 0.999

end


text:
# #############################################################################
# #############################################################################
#PSS clustering methods

#This file contains functions to quantify and analyze the degree of clustering of points on a Poincare Surface of Section, which would indicate the presence of a nearby periodic mode.

#Finding the order of periodicity of a series of bounces from their positions on the PSS is a constrained clustering problem, where the data is already ordered such that the clusters are sampled periodically, so points from the same cluster are separated by uniform intervals equal to the number of clusters.

#   clustcount,clustsepmean = findcluster(bncpts::Array{Float64,2};threshhold::Float64=0.1,samplesize::Int64=30,maxtestinterval::Int64=0)
#   This function estimates the number of clusters in the sequence <bncpts> of (theta,chi) bouncepoint locations on the Poincare Surface of Section, but averaging the pairwise distances of <samplesize> bouncepoints separated by a fixed interval in the sequence, and finding the highest interval not exceeding <maxtestinterval> (no limit if set to 0) for which there is a sharp drop in average distance to a proportion below <threshhold>.
#   The number <clustcount> of clusters and the associated mean separation <clustsepmean> of points within each cluster are returned. If no clear indication of clustering is found (no drop in mean intra-cluster separation below threshhold), then <clustcount> is set to 1.

# #############################################################################
# #############################################################################
#Initiate

#Dependencies
#require("util.jl")
#require("boundary.jl")
#require("refractive_index.jl")
#require("ray_solver.jl")

# #############################################################################
# #############################################################################
#Distance definitions on the Poincare Surface of Section

#Convert incidence angle chi to value that grows approximately linearly with
#the periodicity number of whispering gallery modes in a circular cavity. Returns
#the number of sides of a regular polygon-shaped mode if given the appropriate
#incidence angle (half the interior angle).
function chi2whisporder(chi::Float64)
    return 2*pi/(pi-2*chi)::Float64
end

#Metric on the PSS
#Treats (theta,chi) pairs as positions on a half-infinite cylinder embedded in R^3,
#theta being the angular position and chi2whisporder(chi) being the vertical position.
function pssmetric(theta1::Float64,theta2::Float64,chi1::Float64,chi2::Float64)
    return hypot(2*sin(0.5*(theta1-theta2)),
                 chi2whisporder(chi1)-chi2whisporder(chi2))::Float64
end

#Distance function between 2 points on the PSS
#Column vectors
function pssdist{N1,N2}(p1::Array{Float64,N1},p2::Array{Float64,N2})
    return pssmetric(p1[1],p2[1],p1[2],p2[2])::Float64
end

# #############################################################################
# #############################################################################
#Algorithm to estimate number of clusters

#This algorithm relies on the following properties of clusters in the specific context of k-periodic bounce points:

#   - points an interval of k apart in the bounce sequence are separated by a very
#     small PSS distance.
#   - due to the fractal nature of the PSS, the clusters of bounce points themselves
#     appear in clusters of a divisor of k.

#This means that one can progressively test increasing intervals until a sharp drop in PSS distance indicates a periodic cluster, and subsequently increment the test intervals by the number of clusters, since k should be some multiple of it. The details of the algorithm are as follows:

#   1) Distance between pairs of points a fixed interval apart in the bounce sequence
#      is estimated by the mean of such pairwise distances from all such pairs
#      available in the bounce data. This approximates the distance between their
#      corresponding clusters, so a sudden low value indicates that each pair belongs
#      to the same cluster.
#   2) The cluster count is initially estimated to be 1
#   3) The cluster distance for various intervals of bounce points is computed, with
#      the interval incremented by the cluster count, since that should be a divisor
#      of the "actual" order of periodicity.
#   4) If a sharp drop is detected below a threshhold proportion of the
#      typical cluster distances, then it is likely that there is a clustering
#      present. The cluster count is now estimated to be the current test interval.
#   5) Repeat (3) and (4) until the test intervals is too high relative to the total
#      number of bounce points for meaningful clusters to form. If the cluster count
#      remains 1, then it is inconclusive whether there is any clustering present at
#      all.

#Mean PSS distance of pairs of bouncepoints a specified interval apart
function periodic_pssdist(bncpts::Array{Float64,2},order::Int64,samplesize::Int64=0)
    const maxsamplesize::Int64 = size(bncpts,1)-order
    if samplesize == 0 || samplesize > maxsamplesize
        samplesize = maxsamplesize
    end
    return mean([pssdist(bncpts[i,:],bncpts[i+order,:]) for i=1:samplesize])
end

#Estimate the number of clusters of bouncepoints
function findcluster(bncpts::Array{Float64,2};threshhold::Float64=0.1,samplesize::Int64=30,maxtestinterval::Int64=0)
    #Invalid maximum test interval means no limit
    if maxtestinterval <= 0
        maxtestinterval = size(bncpts,1) - samplesize - 1
    end

    #Initiate
    clustcount::Int64 = 1
    testinterval::Int64 = 2
    clustsepsum::Float64 = periodic_pssdist(bncpts,1,samplesize)
    clustsepcount::Int64 = 1

    while(testinterval <= maxtestinterval)
        #Get mean PSS distance for test interval
        testdist::Float64 = periodic_pssdist(bncpts,testinterval,samplesize)

        #Test against threshhold
        if testdist < threshhold*clustsepsum/clustsepcount
            #Sharp drop detected
            #cluster count likely to be multiple of test interval
            clustcount = testinterval
            #Reset intra-cluster distance data
            clustsepsum = testdist
            clustsepcount = 1
        else
            #No significant drop detected
            #Collect intra-cluster separation data
            clustsepsum += testdist
            clustsepcount += 1
        end

        #Increment test interval by the current cluster count estimate
        testinterval += clustcount
    end

    #Return best estimate of number of clusters, and mean separation of points
    #within each cluster
    return clustcount, clustsepsum/clustsepcount
end


text:
#NetwAnalysis: Circuit element tools
#-------------------------------------------------------------------------------

#==Helper functions
===============================================================================#
#convenient alias... place in CircuitAnalysis?
mag(x) = abs(x)
mag2(x) = abs2(x)

#==Main functions
===============================================================================#

#Compute k stability factor:
function kstab(S::SParameters{2})
	(s11, s12, s21, s22) = mx2elem(S)
	m2det = mag2(s11*s22 - s21*s12)
	return (1 - mag2(s11) - mag2(s22) + m2det)/(2*mag(s12*s21))
end

#Compute Mason's gain (Unilateral gain):
function Ugain(S::SParameters{2})
	k = kstab(S)
	s21_s12 = S[2,1] / S[1,2]
	return mag2(s21_s12 - 1) / (2*k*mag(s21_s12) - 2*real(s21_s12))
end

#Compute Maximum Available Gain:
function MAG(S::SParameters{2})
	k = kstab(S)
	kadj = max(abs(k), 1) #Degenerates MAG into MSG if k<1
	return mag(S[2,1] / S[1,2]) * (kadj - sqrt(kadj*kadj - 1))
end

#Compute Maximum Stable Gain:
MSG(S::SParameters{2}) = mag(S[2,1] / S[1,2])

#==Access gains using a single exported function
===============================================================================#
gain(::DS{:U}, np::Network) = Ugain(np)
gain(::DS{:MA}, np::Network) = MAG(np)
gain(::DS{:MS}, np::Network) = MSG(np)
gain(s::Symbol, np::Network) = gain(DS(s), np)

#==Auto-convert provided networks
===============================================================================#

fnlist = Symbol[:kstab, :MSG, :MAG, :Ugain]
for fn in fnlist; @eval begin #CODEGEN------------------------------------------

$fn(np::Network) = $fn(Network(:S, np))

end; end #CODEGEN---------------------------------------------------------------

#Last line


text:
using Test
using MPI
using Pkg

Pkg.add(url = "https://github.com/clima/climatemachine.jl", rev = "master")
using ClimateMachine

# To test coupling
#using ClimaCoupler.Coupling

# To create meshes (borrowed from Ocean for now!)
using ClimateMachine.Ocean.Domains

# To setup some callbacks
using ClimateMachine.GenericCallbacks

# To invoke timestepper
using ClimateMachine.ODESolvers

Pkg.add(url = "https://github.com/sandreza/QuickVizExample.jl.git", rev = "main")
using QuickVizExample
using GLMakie

ClimateMachine.init()

# Use toy balance law for now
include("CplTestingBL.jl")

# TODO: Use RecordState to produce plots every X steps,
# create time slider for components
mutable struct RecordState
    atmos::Vector{Any}
    ocean::Vector{Any}
end
RecordState() = RecordState(Any[], Any[])
function (rs::RecordState)()
    println("***record callback***")
    push!(rs.atmos, deepcopy(mA.state))
    push!(rs.ocean, deepcopy(mO.state))
end
GenericCallbacks.init!(f::RecordState, solver, Q, param, t) = nothing
function GenericCallbacks.call!(f::RecordState, solver, Q, param, t)
    f()
    return nothing
end
GenericCallbacks.fini!(f::RecordState, solver, Q, param, t) = nothing

rs = RecordState()
recordcallback = EveryXSimulationSteps(rs, 10)

## Visualization
"""
    coupledviz(modellist, plottime)

Create a 3d volume plots for coupled model components.

Assumes Atmos, Ocean, Land ordering in the model list.
"""
function coupledviz(modellist, plottime = 0.0, np_def = 100; title = "")
    states = Array{Any, 1}(nothing, length(modellist))
    statenames = similar(states)
    names = ["θAtmos, t=$plottime", "θOcean, t=$plottime", "θLand, t=$plottime"]
    for (i, model) in enumerate(modellist)
        gh = GridHelper(model.grid)
        x, y, z = coordinates(model.grid)
        xr = range(minimum(x), maximum(x), length = np_def)
        yr = range(minimum(y), maximum(y), length = np_def)
        zr = range(minimum(z), maximum(z), length = np_def)

        # interpolate data
        fld = zeros(np_def, np_def, np_def)
        ϕ = ScalarField(copy(x), gh)
        ϕ .= view(model.state, :, 1, :)
        fld[:, :, :] .= view(ϕ(xr, yr, zr), :, :, :)
        states[i] = fld
        statenames[i] = names[i]
    end
    volumeslice(states, statenames = statenames, title = title)
    return states
end

#include("simple_2testcomp.jl")
models = map(m -> m.component_model, simulation.component_list)
coupledviz(models, 10.0, 100; title = "Test Plot")

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
module LandRegistryPricePaidDB

using CSV
using DataFrames
using Dates
using Pipe
using Serialization
using Index1024

LR = LandRegistryPricePaidDB
export LR

const LANDREGDIR = "/home/matt/wren/UkGeoData/landreg/"
const DATADIR = joinpath(LANDREGDIR, "complete")
const NEWESTDIR = joinpath(LANDREGDIR, "newest_per_address")
const TDATADIR = joinpath(LANDREGDIR, "test")
const mask = 0x00ffffffffffffff

export csv_files, postcode_to_UInt64, create_kvs, load_postcode_price_index, load_postcode_csv, prices_for_postcode

#==
    see https://www.gov.uk/guidance/about-the-price-paid-data#explanations-of-column-headers-in-the-ppd
    TUID
        Transaction unique identifier
        A reference number which is generated automatically recording each published sale.
        The number is unique and will change each time a sale is recorded.
    Property Type
        D = Detached, S = Semi-Detached, T = Terraced, F = Flats/Maisonettes, O = Other
    old / new
        Y = a newly built property, N = an established residential building
    duration
        F = Freehold, L= Leasehold etc.
    paon
        Primary Addressable Object Name (e.g. house number "159")
    saon
        Secondary Addressable Object Name (e.g. "Flat 7")
    ppd category
        A = Standard Price Paid
        B = Additional Price Paid
=##
const header = ["TUID", "price", "date", "postcode", "propertyType", "oldOrNew", "duration", "paon", "saon", "street", "locality", "townCity", "district", "country", "PPDcategory", "recordStatus"]

csv_files(datadir=DATADIR) = filter(f->endswith(f, ".csv"), readdir(datadir))

function read_csvs(datadir=DATADIR)
    files = csv_files(datadir)
    dfs = Vector{DataFrame}(undef, length(files))

    int(t) =  ismissing(t) ? missing : parse(Int64, t)
    date(t) = ismissing(t) ? missing : Date(split(t, " ")[1], "yyyy-mm-dd")

    Threads.@threads for i in 1:length(files)
        dfs[i] = @pipe CSV.read(files[i], DataFrame; header_t, types=String) |>
                    filter!(row->row.propertyType != "O", _) |>
                    transform!(_, [
                        :date => ByRow(t->date(t)) => :date_t,
                        :price => ByRow(t->int(t)) => :price_t,
                        [:postcode, :paon, :saon, :street] => ByRow((a,b,c,d)->hash((a,b,c,d))) => :address_hash
                    ]) |>
                    select!(_, Not([:TUID, :date, :price])) |>
                    rename!(_, [:price_t=>:price, :date_t=>:date])
    end
    dfs
end

function consolidate(dfs)
    df = empty(DataFrame([
        :address_hash=>[zero(UInt64),zero(UInt64)],
        :date => [now(), now()],
        :price => [0, missing],
        :postcode => ["", missing],
        :propertyType => ["", missing],
        :oldOrNew => ["", missing],
        :duration => ["", missing],
        :paon =>  ["", missing],
        :saon  => ["", missing],
        :street => ["", missing],
        :locality => ["", missing],
        :townCity => ["", missing],
        :district => ["", missing],
        :country => ["", missing],
        :PPDcategory => ["", missing],
        :recordStatus => ["", missing]
        ]
        ))
    reduce((df, d)->append!(df, d), dfs, init=df)
end

function create_df(datadir=DATADIR, df_file="LandReg.df")
    @pipe read_csvs(datadir) |> consolidate(_) |> save(joinpath(datadir, df_file), _)
end

function save(io::IO, object)
    serialize(io, object)
    object
end

function save(memofile::AbstractString, object)
    try
        open(memofile, "w+") do io
           return save(io, object)
        end
    catch e
        println(stderr, e)
    end
    object
end

function load(memofile="")
    if memofile == ""
        memofile = joinpath(DATADIR, "LandReg.df")
    end
    if filesize(memofile) > 0
        open(memofile, "r") do io
            return deserialize(io)
        end
    end
end

#==
grep address_hash landreg_hashed.csv > landreg_sorted_uniq.csv
grep -v address_hash landreg_hashed.csv |
    sort -t ',' -k 1,2 |
    awk ' BEGIN { FS=","; previd=""; prev=""}  $1 != previd { print prev; prev=$0; previd=$1; next}  { prev=$0 } ' | grep -v "^$" >> landreg_sorted_uniq.csv

some have no postscodes

awk ' BEGIN { FS="," } $5 == "O" { next } $4 != "" { gsub(/-.*/, "", $2); print $4","$2","$3 }' landreg_sorted_uniq.csv | sort > landreg_pc_year_price.csv
==#

function postcode_to_UInt64(pc)
    if ismissing(pc)
        return 0
    end

    m = match(r"([A-Z]+)([0-9]+([A-Z]+)?) ?([0-9]+)([A-Z]+)", replace(pc, " "=>""))
    if m == nothing || m[1] === nothing || m[2] === nothing || m[4] === nothing || m[5] === nothing
        return 0
    end
    reduce((a,c) -> UInt64(a) << 8 + UInt8(c), collect(lpad(m[1], 2) * lpad(m[2], 2) * m[4] * m[5]), init=0)

end

function UInt64_to_postcode(u)
    if u == 0
        return ""
    end
    cs = Char[]
    while u > 0
        push!(cs, Char(u & 0xff))
        u >>= 8
    end
    part(cs) = replace(String(reverse(cs)), " "=>"")
    "$(part(cs[4:end])) $(part(cs[1:3]))"
end

function count_lines(io, pcode)
    lines = 0
    pos = position(io)
    while (line = readline(io)) !== nothing
        lines += 1
        newcode = split(line, ",", limit=2)[1]
        if newcode != pcode
            return newcode, lines, pos
        end
        pos = position(io)
    end
    return "", lines, pos
end

#==
LR.build_index_file(joinpath(LR.LANDREGDIR, "Postcode_Year_Price.index"), LR.create_kvs(joinpath(LR.LANDREGDIR, "landreg_pc_year_price.csv")))
==#

create_kvs(fname) = open(create_kvs, fname)

function create_kvs(io::IO)
    kvs = Dict{UInt64, DataAux}()
    readline(io)
    pos = position(io)
    pcode, lines, nextpos = count_lines(io, "postcode")
    while pcode != ""
        newcode, lines, nextpos = count_lines(io, pcode)
        kvs[postcode_to_UInt64(pcode)] = (data=pos, aux=lines)
        pcode = newcode
        pos = nextpos
    end
    kvs
end

function csv(io::IO, offset)
    buff = IOBuffer()
    write(buff, readline(io; keep=true))
    seek(io, offset)
    write(buff, readline(io; keep=true))
    seekstart(buff)
    CSV.File(buff)
end

prices_for_postcode(idx, pcode, csvfile) = open(csvfile) do io prices_for_postcode(idx, pcode, io) end

function prices_for_postcode(idx, pcode, csvio::IO)
    (offset, lines) = get(idx, postcode_to_UInt64(pcode), (0,0))
    if lines > 0
        seek(csvio, offset)
        return CSV.File(csvio; header=["postcode", "year", "price"], limit=lines)
    end
end

load_postcode_price_index(idxfile=joinpath(LR.LANDREGDIR, "Postcode_Year_Price.index")) = open_index(idxfile)

load_postcode_csv(csvfile=joinpath(LR.LANDREGDIR, "landreg_pc_year_price.csv")) = open(csvfile, "r")

#####
end


text:
@testset "FileDFG Tests" begin

    if typeof(dfg) <: CloudGraphsDFG
        @warn "TEST: Nuking all data for user '$(dfg.userId)', robot '$(dfg.robotId)'!"
        clearRobot!!(dfg)
    end

    # Same graph as iifInterfaceTests.jl
    numNodes = 10

    #change ready and backendset for x7,x8 for improved tests on x7x8f1
    verts = map(n -> addVariable!(dfg, Symbol("x$n"), ContinuousScalar, labels = [:POSE]), 1:numNodes)
    #TODO fix this to use accessors
    verts[7].ready = 1
    # verts[7].backendset = 0
    verts[8].ready = 0
    verts[8].backendset = 1
    #call update to set it on cloud
    updateVariable!(dfg, verts[7])
    updateVariable!(dfg, verts[8])

    facts = map(n -> addFactor!(dfg, [verts[n], verts[n+1]], LinearConditional(Normal(50.0,2.0))), 1:(numNodes-1))

    # Save and load the graph to test.
    saveFolder = "/tmp/fileDFG"
    saveDFG(dfg, saveFolder)

    copyDfg = DistributedFactorGraphs._getDuplicatedEmptyDFG(dfg)
    retDFG = loadDFG(saveFolder, Main, copyDfg)
    @test symdiff(ls(dfg), ls(retDFG)) == []
    @test symdiff(lsf(dfg), lsf(retDFG)) == []
    for var in ls(dfg)
        @test getVariable(dfg, var) == getVariable(retDFG, var)
    end
    for fact in lsf(dfg)
        @test getFactor(dfg, fact) == getFactor(retDFG, fact)
    end
end


text:
module Proportions
export get_proportion,get_proportion_round, get_proportion_round_add1251,
get_proportion_2d,get_proportion_round_2d,get_proportion_round_add1251_2d

# Include functions
include("onedproportion.jl")
include("twodproportion.jl")

#moduleend
end

max_stars_repo_name: hillelawaskar/Proportions

max_issues_repo_name: hillelawaskar/Proportions


text:
lines = readlines("input.txt")

function main()
    Σ = 0
    pairs = Dict('(' => ')', '[' => ']', '{' => '}', '<' => '>')
    revpairs = Dict(v => k for (k,v) = pairs)
    penalties = Dict(')' => 3, ']' => 57, '}' => 1197, '>' => 25137)
    for line = lines
        stack = []
        for char = collect(line)
            if char in ['(', '[', '{', '<']
                pushfirst!(stack, char)
            else
                if popfirst!(stack) != revpairs[char]
                    Σ += penalties[char]
                    break
                end
            end
        end
    end
    println(Σ)
end

main()


text:
function bubble_sort(list)
	for i = 1:length(list)-1
		for j = 2:length(list)
			if list[j-1] > list[j]
				temp = list[j-1]
				list[j-1] = list[j]
				list[j] = temp
			end
		end
	end
end

list = [
  5009, 2151, 8661, 1568, 4724, 1836, 2525, 1148, 2732, 7560, 1241, 6165, 4884, 4969, 5597, 8714, 3405, 8339, 3865, 9327, 1715, 9840, 1510, 3439, 8560, 4305, 7551, 6657, 5914, 3477, 5012, 3195, 6682, 2223, 7915, 3774, 3661, 1978, 9151, 9218, 7005, 7841, 5598, 3243, 5739, 6202, 8544, 4125, 1292, 6439, 8142, 5216, 4742, 1770, 2835, 1107, 2395, 4264, 9777, 8225, 6213, 1476, 6360, 4129,
  2934, 2096, 8196, 4801, 5428, 4726, 4039, 8155, 5456, 5659, 8303, 6499, 9642, 8508, 4532, 5507, 4716, 6904, 5426, 2308, 3498, 3024, 5507, 9505, 4118, 3516, 8666, 4180, 3233, 6458, 1911, 8371, 8963, 1528, 8271, 5883, 3049, 7134, 2911, 8029, 6763, 7728, 5247, 4264, 1962, 7601, 6742, 8249, 6821, 2077, 4663, 1664, 4230, 4691, 3209, 3541, 3591, 2658, 6742, 7425, 4076, 2686, 1541, 5146,
  5748, 3864, 6798, 1271, 1425, 3397, 1563, 1061, 7045, 1513, 5603, 2916, 1266, 8285, 1104, 8972, 2128, 1341, 9219, 4220, 5169, 4018, 3813, 6854, 6370, 4869, 6557, 9212, 6008, 6972, 4270, 7935, 5306, 4786, 8753, 2026, 4632, 5832, 9073, 1443, 9108, 8582, 8756, 3463, 1800, 7496, 5173, 9386, 5563, 3185, 8928, 4269, 5573, 2582, 6458, 3929, 9681, 1712, 6896, 7071, 5454, 7436, 7988, 5447,
  7240, 1178, 6305, 7924, 4490, 4369, 2167, 8256, 4912, 5328, 4248, 7956, 3761, 7512, 2226, 5158, 3550, 3560, 8528, 5047, 9852, 7706, 3680, 1383, 7105, 1758, 6901, 3760, 2802, 9076, 4453, 3891, 8116, 3584, 2451, 9618, 4305, 5652, 3436, 6745, 2781, 8321, 6263, 1160, 2958, 6136, 6815, 6243, 5317, 8490, 7126, 7363, 1189, 2053, 9585, 4347, 4348, 5647, 2023, 5997, 6711, 3800, 6947, 7592,
  6487, 6078, 7438, 9167, 2777, 7846, 3734, 9722, 5617, 2146, 6478, 4362, 3098, 4868, 9664, 6738, 9049, 1335, 6317, 6340, 6320, 9501, 5063, 3754, 9778, 8433, 5022, 9815, 6200, 1264, 1225, 3022, 9629, 1226, 7355, 6835, 6039, 1769, 1521, 1791, 4664, 5021, 6551, 5974, 8965, 3744, 9947, 2021, 8078, 8416, 9012, 2247, 4668, 7669, 9829, 3502, 5916, 2863, 8092, 4444, 1060, 4544, 7445, 2001,
  6879, 8363, 1614, 4333, 7278, 9626, 9495, 3004, 5264, 1533, 4470, 9047, 5376, 4266, 8675, 7852, 5263, 1498, 6156, 2414, 2668, 1845, 4496, 4280, 5520, 5030, 2657, 3230, 4236, 1515, 1472, 1021, 5326, 9487, 7451, 1948, 8489, 9370, 5715, 6409, 9294, 6592, 7757, 2563, 4819, 6693, 3877, 6089, 8105, 5421, 9836, 5996, 4591, 7729, 9265, 2870, 9557, 9112, 7338, 9183, 3188, 1707, 9989, 6662,
  3639, 9954, 9069, 9811, 7399, 9438, 4727, 5635, 1195, 3364, 3363, 2955, 4434, 9714, 6741, 9966, 7645, 4200, 4537, 5415, 4812, 9243, 9634, 4923, 3973, 4350, 9201, 2803, 3565, 9247, 4841, 5948, 2650, 8830, 2123, 6390, 3294, 3552, 7908, 2299, 7975, 7761, 5386, 8510, 1509, 8219, 8939, 4472, 2697, 5975, 8058, 6845, 4766, 9177, 1653, 1195, 7519, 8131, 6381, 5597, 1843, 6771, 1038, 9968,
  1490, 5557, 2887, 7211, 4390, 6466, 2863, 1686, 7556, 1131, 5934, 4681, 8201, 3625, 4260, 1390, 9363, 4507, 6617, 8289, 9770, 1818, 2352, 7445, 5670, 4084, 9163, 4937, 4531, 6063, 9499, 4452, 7420, 3233, 8327, 8441, 3263, 1479, 8056, 2541, 2000, 5324, 3008, 4025, 1463, 7831, 5282, 2875, 5845, 8560, 1364, 6460, 9908, 8141, 4026, 4785, 5790, 5867, 8689, 4558, 3768, 8161, 5428, 3962,
  1197, 7043, 3046, 7815, 1962, 3073, 8538, 9715, 3138, 7867, 5347, 2604, 1731, 1585, 3086, 4123, 5915, 2153, 4333, 4263, 6088, 4441, 9907, 1278, 6415, 1165, 3609, 2221, 1065, 4935, 4833, 2076, 8166, 6050, 9628, 4355, 5029, 3449, 7211, 9794, 9740, 8661, 9795, 1249, 7096, 8842, 5829, 9564, 4547, 9474, 2782, 6382, 1605, 4199, 3202, 8850, 6417, 3742, 1419, 1629, 2390, 2721, 1319, 6203,
  7930, 1890, 8241, 1657, 5306, 2010, 6724, 7182, 7726, 6614, 3670, 3293, 7823, 2371, 3989, 6220, 8667, 1866, 2431, 2659, 5004, 8134, 8341, 8474, 2779, 5232, 7782, 8510, 9923, 9270, 7241, 6655, 9106, 4995, 8647, 8427, 1383, 5768, 4488, 4688, 9548, 3881, 8993, 4979, 9088, 1088, 5289, 7383, 6305, 1602, 4895, 2045, 8231, 7592, 1804, 1629, 7694, 5881, 4805, 6182, 7377, 8979, 6146, 1698,
  9268, 8164, 8937, 7348, 6559, 4423, 8099, 5189, 8676, 9755, 1602, 8921, 9494, 7532, 7139, 7976, 7825, 3690, 2935, 4634, 8055, 6601, 7828, 4864, 8446, 9283, 8175, 2216, 1245, 1453, 1561, 8150, 2453, 7828, 9967, 2076, 9556, 5403, 9028, 1511, 8237, 9480, 7756, 6223, 9596, 3636, 7228, 4518, 2362, 7231, 6579, 7099, 3550, 9001, 2797, 5155, 3839, 8427, 6822, 5968, 5604, 6020, 2778, 8037,
  8317, 4659, 3224, 6910, 8579, 3135, 4045, 6795, 7163, 1647, 7260, 1448, 5960, 3977, 8424, 9793, 6135, 9838, 9960, 6119, 3882, 1012, 6361, 6252, 1337, 2729, 6221, 9994, 1669, 9593, 5606, 9172, 5372, 8750, 4153, 4803, 7772, 9880, 8232, 7063, 8577, 5745, 8282, 5241, 8589, 5330, 3159, 9008, 5618, 7419, 5230, 4797, 1080, 6773, 8960, 5415, 1530, 8668, 1689, 4227, 3476, 4840, 7961, 6386,
  6103, 9897, 1836, 4360, 8739, 2633, 7004, 3316, 1244, 7784, 5915, 7486, 4192, 8802, 4255, 3870, 3138, 5355, 9492, 2442, 2879, 9844, 2740, 6191, 3751, 6196, 1063, 3931, 4938, 5843, 7759, 8737, 8649, 7810, 4539, 2693, 1785, 6162, 1362, 1031, 7434, 8358, 1197, 3578, 5971, 4152, 6080, 8857, 7466, 5311, 9790, 6831, 5997, 2067, 1090, 3042, 3803, 7063, 3936, 1470, 1923, 7904, 7103, 2185,
  7411, 3783, 5014, 9079, 6087, 6186, 9424, 8085, 6179, 1155, 3628, 9364, 7631, 7227, 5364, 6403, 5924, 7492, 5250, 8952, 4945, 4785, 8873, 4576, 4690, 3316, 1712, 9549, 4447, 4666, 1695, 3888, 1134, 8975, 4472, 6537, 9300, 1945, 4858, 7132, 4891, 3020, 9108, 3508, 6604, 3523, 5549, 4694, 5863, 5102, 3746, 4961, 9445, 7743, 1767, 5509, 6355, 6185, 2852, 8361, 1587, 8298, 3645, 7156,
  3902, 4901, 2803, 9304, 1002, 4738, 3049, 2348, 3856, 6261, 7433, 2570, 3660, 5426, 6520, 6581, 6248, 9069, 3419, 1800, 4622, 6090, 1425, 3479, 2757, 2066, 6076, 7316, 3049, 8435, 6625, 3616, 4401, 7867, 5556, 6180, 3632, 9458, 5429, 9861, 4853, 7684, 3424, 8337, 3896, 5021, 8799, 6578, 2809, 1558, 7297, 4779, 2088, 1007, 7426, 8640, 4547, 8349, 4391, 3040, 5975, 9359, 4605, 8286,
  3712, 6009, 3001, 3510, 9618, 6217, 8746, 6130, 2280, 3199, 8769, 9672, 1707, 8449, 7622, 6023, 8071, 1330, 3861, 2919, 6879, 5535, 6461, 3720, 6053, 2198, 8350, 2781, 2516, 4499, 3310, 6990, 6501, 9443, 3765, 4101, 2659, 5805, 5480, 2686, 7898, 6521, 4077, 1965, 4149, 5976, 2511, 6315, 2804, 1448, 1773, 5527, 6001, 3619, 2434, 7911, 2392, 6081, 1074, 8116, 2505, 9600, 3390, 3385,
  5965, 4980, 7226, 4918, 6358, 4069, 3210, 5225, 9473, 9291, 8450, 3246, 9434, 9228, 9587, 6687, 8184, 9159, 7564, 5839, 5097, 1322, 9682, 5790, 9584, 7184, 9055, 1653, 4620, 5953, 6896, 2568, 2885, 2277, 4575, 5942, 9809, 3897, 4615, 5330, 6766, 7800, 3781, 1114, 5276, 2384, 3770, 9676, 4318, 6325, 1656, 7545, 7658, 4878, 7466, 8522, 7339, 8346, 7116, 3016, 6282, 2002, 1457, 8464,
  3866, 5586, 1569, 6358, 3823, 7483, 2508, 6827, 7501, 5763, 1514, 5719, 4113, 7288, 6731, 6898, 7319, 6252, 2693, 9824, 4701, 5593, 9554, 3137, 4698, 4231, 7582, 3348, 8853, 2808, 5346, 3277, 1826, 9368, 3490, 6833, 4128, 6983, 9501, 9007, 6902, 7465, 3427, 4947, 4079, 9338, 5386, 6420, 1844, 8005, 1255, 9742, 9907, 9580, 7944, 8139, 6791, 9127, 8062, 6518, 9382, 8554, 8766, 1567,
  9515, 9517, 5179, 9148, 3786, 6371, 8665, 6044, 7829, 7588, 2862, 4081, 1310, 7413, 3923, 7623, 3825, 2154, 2745, 4512, 8764, 7786, 3270, 2532, 2533, 7132, 2706, 2542, 8408, 2390, 9616, 6182, 7923, 4179, 9654, 2457, 7637, 6976, 5101, 6710, 8616, 6050, 9542, 7995, 2446, 5053, 1486, 8387
]

bubble_sort(list)

max_stars_repo_name: bichanna/slap

max_issues_repo_name: adiyenuubarii/slap

max_forks_repo_name: adiyenuubarii/slap


text:
module MyAstroUtils

using Reexport
@reexport using FITSIO
@reexport using DataFrames
@reexport using SortMerge

include("misc.jl")
include("coords.jl")
include("fits.jl")
include("DB.jl")

end # module


text:
# This code demonstrates the use of SolveDSGE for solving linear rational expectations models.
# The demonstration model is a simple forward-looking new Keynesian model with monetary policy
# conducted using a Taylor rule.

# The model is as follows:

# pi_{t} = beta*E_{t}pi_{t+1} + kappa*y_{t} + u_{t}
#  y_{t} = E_{t}y_{t+1} - mu*(R_{t}-E_{t}pi_{t+1} + v_{t}
#  R_{t} = phi_pi*E_{t}pi_{t+1} + phi_y*y_{t}

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

# Bring the SolveDSGE module into global scope

using SolveDSGE

# Assign values to model parameters

const beta = 0.99
const mu = 0.5
const phi_pi = 1.5
const rho_u = 0.8
const rho_v = 0.8
const kappa = 0.1
const phi_y = 0.5
const sigma = [1.0 0.0; 0.0 1.0]

const cutoff = 1.0
const tol    = 1e-8

# Solve the model using Blanchard and Kahn (1980)

# Construct the model matrices

b = [1.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0; 0.0 0.0 beta 0.0; 0.0 0.0 mu*(1.0-phi_pi) 1.0]
a = [rho_u 0.0 0.0 0.0; 0.0 rho_v 0.0 0.0; -1.0 0.0 1.0 -kappa; 0.0 -1.0 0.0 (1.0+mu*phi_y)]
c = [1.0 0.0; 0.0 1.0; 0.0 0.0; 0.0 0.0]

c = b\c
a = b\a

nx = 2
ny = 2

# Put the model into Blanchard_Kahn_Form type

m_bk = Blanchard_Kahn_Form(nx,ny,a,c,sigma)

# Do some basic checking to determine whether the model matrices are conformable

assessment = check_model_form(m_bk)
println(assessment)

# Solve the model by solving an eigenvalue problem

soln_bk_1 = solve_re(m_bk,cutoff)
responses_bk_1 = impulses(soln_bk_1,5,1)

# Solve the model by solving an algebraic Riccati equation

soln_bk_2 = solve_re(m_bk,cutoff,tol)
responses_bk_2 = impulses(soln_bk_2,5,1)

# Solve the model using Klein (2000)

# Construct the model matrices

b = [1.0 0.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0 0.0; 0.0 0.0 beta 0.0 0.0; 0.0 0.0 mu 1.0 0.0; 0.0 0.0 phi_pi 0.0 0.0]
a = [rho_u 0.0 0.0 0.0 0.0; 0.0 rho_v 0.0 0.0 0.0; -1.0 0.0 1.0 -kappa 0.0; 0.0 -1.0 0.0 1.0 mu; 0.0 0.0 0.0 -phi_y 1.0]
c = [1.0 0.0; 0.0 1.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]

nx = 2
ny = 3

# Put the model into Klein_Form type

m_k = Klein_Form(nx,ny,a,b,c,sigma)

# Do some basic checking to determine whether the model matrices are conformable

assessment = check_model_form(m_k)
println(assessment)

# Solve the model by solving a generalized eigenvalue problem

soln_k = solve_re(m_k,cutoff)
responses_k = impulses(soln_k,5,1)

# Solve the model using Sims (2001)

# Construct the model matrices

gamma0 = [1.0 0.0 0.0 0.0 0.0 0.0 0.0;
          0.0 1.0 0.0 0.0 0.0 0.0 0.0;
          0.0 0.0 1.0 0.0 0.0 0.0 0.0;
          0.0 0.0 0.0 1.0 0.0 0.0 0.0;
          0.0 0.0 0.0 -phi_y 1.0 -phi_pi 0.0;
          -1.0 0.0 1.0 -kappa 0.0 -beta 0.0;
          0.0 -1.0 0.0 1.0 mu -mu -1.0]

gamma1 = [rho_v 0.0 0.0 0.0 0.0 0.0 0.0;
          0.0 rho_u 0.0 0.0 0.0 0.0 0.0;
          0.0 0.0 0.0 0.0 0.0 1.0 0.0;
          0.0 0.0 0.0 0.0 0.0 0.0 1.0;
          0.0 0.0 0.0 0.0 0.0 0.0 0.0;
          0.0 0.0 0.0 0.0 0.0 0.0 0.0;
          0.0 0.0 0.0 0.0 0.0 0.0 0.0]

c = zeros(7)

psi = [1.0 0.0;
       0.0 1.0;
       0.0 0.0;
       0.0 0.0;
       0.0 0.0;
       0.0 0.0;
       0.0 0.0]

pi = [0.0 0.0;
      0.0 0.0
      1.0 0.0
      0.0 1.0;
      0.0 0.0;
      0.0 0.0;
      0.0 0.0]

# Put the model into Sims_Form type

m_sims = Sims_Form(gamma0,gamma1,c,psi,pi,sigma)

# Do some basic checking to determine whether the model matrices are conformable

assessment = check_model_form(m_sims)
println(assessment)

# Solve the model by solving a generalized eigenvalue problem

soln_sims = solve_re(m_sims,cutoff)

# Solve the model in structural form

# Construct the model matrices

a0 = [1.0 0.0 0.0 0.0 0.0; 0.0 1.0 0.0 0.0 0.0; -1.0 0.0 1.0 -kappa 0.0; 0.0 -1.0 0.0 1.0 mu; 0.0 0.0 0.0 -phi_y 1.0]
a1 = [rho_u 0.0 0.0 0.0 0.0; 0.0 rho_v 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0]
a2 = [0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 beta 0.0 0.0; 0.0 0.0 mu 1.0 0.0; 0.0 0.0 phi_pi 0.0 0.0]
a4 = [1.0 0.0; 0.0 1.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]

# Put the model in Structural_Form type

m_bp = Binder_Pesaran_Form(a0,a1,a2,a4,sigma)

# Do some basic checking to determine whether the model matrices are conformable

assessment = check_model_form(m_bp)
println(assessment)

# Solve the by by solving a generalized eigenvalue problem

soln_bp_1 = solve_re(m_bp,cutoff)
responses_bp_1 = impulses(soln_bp_1,5,1)

# Solve the model using Binder and Pesaran's (1995) brute force method

soln_bp_2 = solve_re(m_bp,cutoff,tol)
responses_bp_2 = impulses(soln_bp_2,5,1)


text:
using SpecialFunctions
using DiffRules: DiffRules, @define_diffrule, diffrule, diffrules, hasdiffrule

# Hand code the identity because it's really fundamental. It doesn't need to generate a new
# node on the computational graph since it does nothing, but it is useful to have it's
# gradient implemented for use in higher-order functions.
import Base.identity
@explicit_intercepts identity Tuple{Any}
@inline ∇(::typeof(identity), ::Type{Arg{1}}, p, y, ȳ, x) = ȳ
@inline ∇(::typeof(identity), ::Type{Arg{1}}, x::Real) = one(x)

# Ignore functions that have complex ranges. This may change when Nabla supports complex
# numbers.
ignored_fs = [(:SpecialFunctions, :hankelh1),
              (:SpecialFunctions, :hankelh2),
              (:Base, :log1p),
              (:Base, :rem2pi),
              (:Base, :mod),
              (:Base, :atan),
              (:Base, :rem)]

unary_sensitivities, binary_sensitivities = [], []

for (package, f, arity) in diffrules()
    (package == :NaNMath || (package, f) in ignored_fs) && continue

    @eval import $package: $f
    if arity == 1
        push!(unary_sensitivities, (package, f))
        ∂f∂x = diffrule(package, f, :x)
        @eval @explicit_intercepts $f Tuple{∇Scalar}
        @eval @inline ∇(::typeof($f), ::Type{Arg{1}}, p, y, ȳ, x::∇Scalar) = ȳ * $∂f∂x
        @eval @inline ∇(::typeof($f), ::Type{Arg{1}}, x::∇Scalar) = $∂f∂x
    elseif arity == 2
        push!(binary_sensitivities, (package, f))
        ∂f∂x, ∂f∂y = diffrule(package, f, :x, :y)
        @eval @explicit_intercepts $f Tuple{∇Scalar, ∇Scalar}
        @eval ∇(::typeof($f), ::Type{Arg{1}}, p, z, z̄, x::∇Scalar, y::∇Scalar) = z̄ * $∂f∂x
        @eval ∇(::typeof($f), ::Type{Arg{2}}, p, z, z̄, x::∇Scalar, y::∇Scalar) = z̄ * $∂f∂y
    else
        error("Cannot implement sensitivity for $package.$f: arity $arity not supported.")
    end
end

# Add method to resolve exponentiation ambiguity.
^(n::Node{<:Real}, p::Integer) = invoke(^, Tuple{Node{<:Real}, Real}, n, p)

import Base: float
@explicit_intercepts float Tuple{∇ArrayOrScalar}
∇(::typeof(float), ::Type{Arg{1}}, p, y, ȳ, x) = float(ȳ)


text:
using ViscosityDrop

### Setup simulation and run
filename = "example"
dir = "."
Ny, Nz = 64, 64
do_all(filename, dir, Ny=Ny, Nz=Nz, yLength=2, zLength=2,
                      ν=1e-1,
                      z_position=0.6, width=0.2, amplitude=-1e-1,
                      Δt=0.002, stop_time=120.0, iteration_interval=1000,
                      schedule_interval=0.5, output_fields = :all)

### Open simulation
using JLD2
file = jldopen("./example.jld2")

Ly = file["grid/Ly"]
Lz = file["grid/Lz"]
yb = range(-Ly/2,Ly/2,length=Ny)
zb = range(-Lz,0,length=Nz)

iterations = parse.(Int, keys(file["timeseries/t"]))

### Make animations
using Printf
using Plots

# Buoyancy
anim1 = @animate for iter in iterations
    t = file["timeseries/t/$iter"]*7.4e-2
    b_snapshot = file["timeseries/b/$iter"][1, :, :]

    contourf(yb, zb, b_snapshot',color=:inferno,
             title="Buoyancy, t=$(@sprintf("%.2f", t)) ms",
             aspect_ratio=:equal)
end
gif(anim1, "buoyancy.gif", fps = 15)

# Vorticity
anim2 = @animate for iter in iterations
    t = file["timeseries/t/$iter"]*7.4e-2
    ω_snapshot = file["timeseries/ω/$iter"][1, 1:Ny, 1:Nz]

    contourf(yb, zb, ω_snapshot',color=:balance,
             title="Vorticity, t=$(@sprintf("%.2f", t)) ms",
             aspect_ratio=:equal)
end
gif(anim2, "vorticity.gif", fps = 15)


text:
using DataStructures

#####
##### Abstract type
#####

abstract type DynamicalSystemModel end

#####
##### SINDY model object and methods
#####
mutable struct sindyModel <: DynamicalSystemModel
    ts::Vector
    # After fitting
    A::Matrix
    # For augmenting data
    library::OrderedDict{Function,Any}
    # For printing equations
    variable_names
    # For initial training and retraining
    optimizer::SparseSolver
end

# Convenience method with a default BAD optimizer
sindyModel(ts, A, lib, var) =
    sindyModel(ts, A, lib, var, denseSolver())
# Initializer for creating templates
z = zeros(1,1)
sindyModel(lib, var, opt::SparseSolver) =
    sindyModel([0], z, lib, var, opt)

(m::sindyModel)(X) = m.A*augment_data(m, X)
(m::sindyModel)(X, t) = m.A*augment_data(m, X) # To match syntax below
(m::sindyModel)(u::AbstractArray,p,t) = m(u, t) # OrdinaryDiffEq syntax

# Functions for both uncontrolled and controlled
augment_data(m::DynamicalSystemModel, X) =
    calc_augmented_data(X, m.library)

#####
##### SINDYc model object and methods
#####
mutable struct sindycModel <: DynamicalSystemModel
    ts::Vector
    # After fitting
    A::Matrix
    B::Matrix
    # Interpolated Control signal
    U::Matrix
    U_func::Function
    # For augmenting data
    library::OrderedDict{Function,Any}
    # For printing equations
    variable_names
    # For initial training and retraining
    optimizer::SparseSolver
end

# Convenience method with a default BAD optimizer
sindycModel(ts, A, B, U, Uf, lib, var) =
    sindycModel(ts, A, B, U, Uf, lib, var, denseSolver())
# Initializer for creating templates
z = zeros(1,1)
sindycModel(lib, var, opt::SparseSolver) =
    sindycModel([0], z, z, z, ()->0, lib, var, opt)

(m::sindycModel)(X) = m.A*augment_data(m, X) .+ m.B*m.U
(m::sindycModel)(X, t) = m.A*augment_data(m, X) .+ m.B*m.U_func(t)
(m::sindycModel)(u::AbstractArray,p,t) = m(u, t) # OrdinaryDiffEq syntax

intrinsic_dynamics(m::sindycModel, X) = m.A*augment_data(m, X)
control_signal(m::sindycModel, t) = m.B*m.U_func(t)
control_signal(m::sindycModel) = m.B*m.U

# Basic integration functions
function simulate_model(m::DynamicalSystemModel, u0)
    condition(u,t,integrator) = any(abs.(u).>1e4)
    cb = DiscreteCallback(condition, terminate!)
    prob = ODEProblem(m, u0, (m.ts[1], m.ts[end]), [0], callback=cb)
    return Array(solve(prob, Tsit5(), saveat=m.ts));
end

###
### Export
###
export sindyModel, sindycModel,
    augment_data, intrinsic_dynamics, control_signal,
    simulate_model


text:
using AutoFeatures
using Test

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


text:
using Distributions, Test

@testset "truncated DiscreteUniform" begin
    # just test equivalence of truncation results
    bounds = [(1, 10), (-3, 7), (-5, -2)]
    @testset "lower=$lower, upper=$upper" for (lower, upper) in bounds
        d = DiscreteUniform(lower, upper)
        @test truncated(d, -Inf, Inf) == d
        @test truncated(d, nothing, nothing) === d
        @test truncated(d, lower - 0.1, Inf) == d
        @test truncated(d, lower - 0.1, nothing) == d
        @test truncated(d, -Inf, upper + 0.1) == d
        @test truncated(d, nothing, upper + 0.1) == d
        @test truncated(d, lower + 0.3, Inf) == DiscreteUniform(lower + 1, upper)
        @test truncated(d, lower + 0.3, nothing) == DiscreteUniform(lower + 1, upper)
        @test truncated(d, -Inf, upper - 0.5) == DiscreteUniform(lower, upper - 1)
        @test truncated(d, nothing, upper - 0.5) == DiscreteUniform(lower, upper - 1)
        @test truncated(d, lower + 1.5, upper - 1) == DiscreteUniform(lower + 2, upper - 1)
    end
end


text:
using Gadfly

n = 10
plot(x=rand(n), y=rand(n), color=8.^rand(n),
     Scale.color_log2(minvalue=1, maxvalue=8))


text:
const ArrayLike = Union{AbstractArray, Tuple}
const StringLike = Union{AbstractString, Symbol}

function cycle(v::AbstractVector, i::Int)
    ax = axes(v, 1)
    return v[first(ax) + mod(i - first(ax), length(ax))]
end

"""
    iscontinuous(v::AbstractArray)

Determine whether `v` should be treated as a continuous or categorical vector.
"""
iscontinuous(::AbstractArray) = false
iscontinuous(::AbstractArray{<:Number}) = true

for sym in [:hidexdecorations!, :hideydecorations!, :hidedecorations!]
    @eval function $sym(ae::AxisEntries; kwargs...)
        axis = Axis(ae)
        isnothing(axis) || $sym(axis; kwargs...)
    end
end

for sym in [:linkxaxes!, :linkyaxes!, :linkaxes!]
    @eval function $sym(ae::AxisEntries, aes::AxisEntries...)
        axs = filter(!isnothing, map(Axis, (ae, aes...)))
        $sym(axs...)
    end
end

function hideinnerdecorations!(aes::Matrix{AxisEntries})
    options = (label=true, ticks=true, minorticks=true, grid=false, minorgrid=false)
    foreach(ae -> hidexdecorations!(ae; options...), aes[1:end-1, :])
    foreach(ae -> hideydecorations!(ae; options...), aes[:, 2:end])
end

function deleteemptyaxes!(aes::Matrix{AxisEntries})
    for (i, ae) in enumerate(aes)
        if isempty(ae.entries)
            axis = Axis(ae)
            if !isnothing(axis)
                delete!(axis)
                aes[i] = AxisEntries(nothing, ae.entries, ae.labels, ae.scales)
            end
        end
    end
end

extend_extrema((l1, u1), (l2, u2)) = min(l1, l2), max(u1, u2)

function assert_equal(a, b)
    @assert a == b
    return a
end

function unnest(arr::AbstractArray{<:AbstractArray})
    inner_size = mapreduce(size, assert_equal, arr)
    outer_size = size(arr)
    flattened = reduce(vcat, map(vec, vec(arr)))
    return reshape(flattened, inner_size..., outer_size...)
end

unnest(arr::NTuple{<:Any, <:AbstractArray}) = unnest(collect(arr))

adjust_index(rg1, rg2, idx::Integer) = idx in rg2 ? idx : only(rg2)
adjust_index(rg1, rg2, idxs::AbstractArray) = map(idx -> adjust_index(rg1, rg2, idx), idxs)
adjust_index(rg1, rg2, ::Colon) = rg1 == rg2 ? Colon() : fill(only(rg2), length(rg1))

maybewrap(x::ArrayLike) = x
maybewrap(x) = fill(x)

unwrap(x) = x
unwrap(x::AbstractArray{<:Any, 0}) = x[]


text:
# Code for factor analysis model. We write the code for it seperately because it has a very special structure.

## Define the one shot function FactorAnalysisModel

struct ProblemFactorAnalysisModel #{A <: AbstractMatrix{ <: Real }, I <: Integer, R <: Real, V <: AbstractVector{ <:Real }}
    Σ #::A # This is part of the data, which is a positive semidefinite matrix
    r #::I # represents rank of the decision variable matrix X
    M #::R # upper bound on the operator-2 norm of decision variable matrix X
    Z0 #::A # one initial state provided by the user (for the matrix part)
    z0 #::V # one initial state provided by the user (the vector part)
end

## This structure contains all the necessary information to describe a state

mutable struct StateFactorAnalysisModel #{R <: Real, A <: AbstractMatrix{ <: Real }, V <: AbstractVector{ <: Real}, I <: Integer} # Done
    X #::A # the first iterate matrix part
    x #::V # the first iterate vector part
    Y #::A # the second iterate matrix part
    y #::V # the second iterate vector part
    Z #::A # the third iterate matrix part
    z #::V # the third iterate vector part
    i #::I # iteration number
    fxd_pnt_gap #::R # fixed point gap for state
    fsblt_gap #::R # how feasible the current point is, this is measured by || x- Π_C(x) ||
    μ #::R # current value of μ
    γ #::R # current value of γ
end

## this function constructs the very first state for the problem and settings for the factor analysis model
function StateFactorAnalysisModel(problem::ProblemFactorAnalysisModel, settings::Settings)

    # this function constructs the very first state for the problem and settings
    Z0 = copy(problem.Z0)
    z0 = copy(problem.z0)
    X0 = zero(Z0)
    x0 = zero(z0)
    Y0 = zero(Z0)
    y0 = zero(z0)
    i = copy(settings.n_iter_max) # iteration number where best fixed point gap is reached
    fxd_pnt_gap = copy(settings.big_M)
    fsblt_gap = copy(settings.big_M)
    μ = copy(settings.μ_max)
    γ = γ_from_μ(μ, settings)

    return StateFactorAnalysisModel(X0, x0, Y0, y0, Z0, z0, i, fxd_pnt_gap, fsblt_gap, μ, γ)

end

## structure that describes the information to intialize our algorithm:

mutable struct InitInfoFactorAnalysisModel #{ A <: AbstractMatrix{ <: Real }, V <: AbstractVector{ <: Real}, R <: Real}

    # this will have Z0, z0, γ, μ
    Z0 #::A # initial condition (matrix part) to start an inner iteration
    z0 #::V # initial condition (vector part) to start an inner iteration
    μ #::R # μ required start an inner iteration
    γ #::R # γ required to start an inner iteration

end

## Of course, we will need a constructor for the InitInfoFactorAnalysisModel type, given the user input. The user inputs are problem::problem, settings::settings.

function InitInfoFactorAnalysisModel(problem::ProblemFactorAnalysisModel, settings::Settings)

    ## constructs the initial NExOS information
    Z0 = copy(problem.Z0)
    z0 = copy(problem.z0) # need to understand if copy is necessary
    μ = copy(settings.μ_max)
    γ = γ_from_μ(μ, settings)

    return InitInfoFactorAnalysisModel(Z0, z0, μ, γ)

end

## update the state for factor analysis model

function update_state_fam!(state::StateFactorAnalysisModel, init_info::InitInfoFactorAnalysisModel, problem::ProblemFactorAnalysisModel, settings::Settings)

    # extract information from the state and init info
    X = state.X
    x = state.x
    Y = state.Y
    y = state.y
    Z = init_info.Z0
    z = init_info.z0 # TODO: Need to write the InitInfoFactorAnalysisModel
    β = settings.β
    γ = init_info.γ
    μ = init_info.μ

    # create best point variables
    best_point_X = X
    best_point_x = x
    best_point_Y = y
    best_point_y = y
    best_point_Z = Z
    best_point_z = z
    best_fxd_pnt_gap = settings.big_M
    best_fsblt_gap = settings.big_M

    i = 1
    while i <= settings.n_iter_max
        X, x, Y, y, Z, z = inner_iteration_fam(X, x, Y, y, Z, z, β, γ, μ, problem) # Done: Need to write inner iteration factor analysis model: inner_iteration_fam
        crnt_fxd_pnt_gap = norm(x-y, Inf)+norm(X-Y, Inf)
        Y_fsbl = Π_exact(RankSet(problem.M, problem.r), Y)
        crnt_fsblt_gap = norm(Y-Y_fsbl, Inf) # the vector part y is not changed in the second step of the factor analysis inner iteration, so we do not need norm(y-y_fsbl, Inf) part

        # update the best points so far if we have lower objective value
        if crnt_fxd_pnt_gap <= best_fxd_pnt_gap
            best_point_X = X
            best_point_x = x
            best_point_Y = Y
            best_point_y = y
            best_point_Z = Z
            best_point_z = z
            best_fxd_pnt_gap = crnt_fxd_pnt_gap
            best_fsblt_gap = crnt_fsblt_gap
        end # if

        # display important information
        if settings.verbose == true
            if mod(i, settings.freq) == 0
                @info "μ = $μ | iteration = $i | log fixed point gap = $(log10(crnt_fxd_pnt_gap)) | log feasibility gap = $(log10(crnt_fsblt_gap))"
            end # if
        end # if

        # termination condition check
        if crnt_fxd_pnt_gap <= settings.tol
            if settings.verbose == true
                @info "μ = $μ | log fixed point gap reached $(log10(settings.tol))"
            end
            break
        end #if

        # increment the count by 1
        i = i+1

    end # while

    ## update the state
    state.X = best_point_X
    state.x = best_point_x
    state.Y = best_point_Y
    state.y = best_point_y
    state.Z = best_point_Z
    state.z = best_point_z
    state.i = i
    state.fxd_pnt_gap = best_fxd_pnt_gap
    state.fsblt_gap = best_fsblt_gap
    state.μ = μ
    state.γ = γ

    # Print information about the best point
    if settings.verbose == true
        @info "information about the best state found for μ = $(state.μ)"
        @info "μ = $(state.μ) | log fixed point gap = $(log10(state.fxd_pnt_gap)) | log feasibility gap = $(log10(state.fsblt_gap)) | inner iterations = $(state.i)"
    end

    return state
end

# inner iteration function for FactorAnalysisModel
function inner_iteration_fam(X, x, Y, y, Z, z, β, γ, μ, problem::ProblemFactorAnalysisModel)
    #(X::A, x::V, Y::A, y::V, Z::A, z::V, β::R, γ::R, μ::R, problem::ProblemFactorAnalysisModel) where { A <: AbstractMatrix{ <: Real }, V <: AbstractVector{ <: Real}, R <: Real}
    # old iteration: x = prox_NExOS(problem.f, β, γ, z)
    X, x = prox_NExOS_fam(problem.Σ, problem.M, γ, Z, z)
    Y = Π_NExOS(RankSet(problem.M, problem.r), β, γ, μ, 2*X - Z)
    y = max.(2*x - z,0)
    Z = Z + Y - X
    z = z + y - x
    # println("size of the new matrices = ", size(Z), size(z))
    return X, vec(x), Y, vec(y), Z, vec(z)
end

# Proximal operator funciton evaluation using Convex, this is much slower than using JuMP, so I am commenting it for now.

# function prox_NExOS_fam(Σ, M, γ, X, d) #(Σ::A, M::R, γ::R, X::A, d::V) where {R <: Real, A <: AbstractMatrix{R}, V <:  AbstractVector{R}} # For now M is not used, may use it in a future version
#
#   # This functions takes the input data Σ, γ, X, d and evaluates
#   # the proximal operator of the function f at the point (X,D)
#
#   # Data extraction
#   # ---------------
#   n = length(d) # dimension of the problem
#   # println("*****************************")
#   # println(size(d))
#   # println("the value of d is = ", d)
#   # println("the type of d is", typeof(d))
#   D = LinearAlgebra.diagm(d) # creates the diagonal matrix D that embed
#
#   # Create the variables
#   #  --------------------
#   X_tl = Convex.Semidefinite(n) # Here Semidefinite(n) encodes that
#   # X_tl ≡ ̃X is a positive semidefinite matrix
#   d_tl = Convex.Variable(n) # d_tl ≡ ̃d
#   D_tl = diagm(d_tl) # Create the diagonal matrix ̃D from ̃d
#
#   # Create terms of the objective function, which we write down
#   #  in three parts
#   #  ----------------------------------------------------------
#   t1 = square(norm2(Σ - X_tl - D_tl)) # norm2 computes Frobenius norm in Convex.jl
#   t2 = square(norm2(X-X_tl))
#   t3 = square(norm2(D-D_tl))
#
#   # Create objective
#   # ----------------
#   objective = t1 + (1/(2*γ))*(t2 + t3) # the objective to be minimized
#
#   # create the problem instance
#   # ---------------------------
#   convex_problem = Convex.minimize(objective, [d_tl >= 0, Σ - D_tl  in :SDP])
#
#   # set the solver
#   # --------------
#   convex_solver = () -> SCS.Optimizer(verbose=false)#COSMO.Optimizer(verbose=false)
#
#   # solve the problem
#   # -----------------
#   Convex.solve!(convex_problem, convex_solver)
#
#   # get the optimal solution
#   # ------------------------
#   X_sol = X_tl.value
#   d_sol = d_tl.value
#   # println("d_sol = ", d_sol)
#
#   # return the output
#   return X_sol, vec(d_sol)
#
# end

## put everything in a function
# uses JuMP
function prox_NExOS_fam(Σ, M, γ, X, d; X_tl_sv = nothing, d_tl_sv = nothing, warm_start = false)

	# This functions takes the input data Σ, γ, X, d and evaluates the proximal operator of the function f at the point (X,d)

	n = length(d)

	# prox_model = JuMP.Model(optimizer_with_attributes(COSMO.Optimizer, "verbose" => false))

	prox_model = JuMP.Model(optimizer_with_attributes(Mosek.Optimizer))

	@variables( prox_model,
	begin
		d_tl[1:n] >= 0
		X_tl[1:n, 1:n], PSD
	end
	)

	if warm_start == true
		println("warm start enabled")
		# Set warm-start
		set_start_value.(X_tl, X_tl_sv) # Warm start
		set_start_value.(d_tl, d_tl_sv) # Warm start
	    println("norm difference is = ", norm(start_value.(X_tl) - X_tl_sv))
	end

    t_1 = vec(Σ - X_tl - diagm(d_tl))
	t_2 = vec(X_tl-X)
	t_3 = vec(diagm(d_tl)-diagm(d))
	obj = t_1'*t_1 + ((1/(2*γ))*(t_2'*t_2 + t_3'*t_3))

	@objective(prox_model, Min, obj)

	@constraints(prox_model, begin
		Symmetric(Σ - diagm(d_tl)) in PSDCone()
	end)

	set_silent(prox_model)

	JuMP.optimize!(prox_model)

	# obj_val = JuMP.objective_value(prox_model)
	X_sol = JuMP.value.(X_tl)
	d_sol = JuMP.value.(d_tl)

	return X_sol, d_sol

end

function update_init_info_experimental_fam!(state::StateFactorAnalysisModel, init_info::InitInfoFactorAnalysisModel, problem::ProblemFactorAnalysisModel, settings::Settings)

    init_info.μ = init_info.μ*settings.μ_mult_fact
    init_info.γ = γ_from_μ(init_info.μ, settings)
    init_info.z0 = state.z
    init_info.Z0 = state.Z

    return init_info

end


text:
using LinearAlgebra, DiffEqOperators, Random, Test
################################################################################
# Test BoundaryPaddedArray up to 5 dimensions
################################################################################

for dimensionality in 2:5
    for dim in 1:dimensionality
        sizes = rand(4:10, dimensionality)
        A = rand(sizes...)
        lower = Array(selectdim(A, dim, 1))
        upper = Array(selectdim(A, dim, size(A)[dim]))

        Apad = DiffEqOperators.BoundaryPaddedArray{Float64, dim, dimensionality, dimensionality-1, typeof(A), typeof(lower)}(lower, upper, selectdim(A, dim, 2:(size(A)[dim]-1)))

        @test A == Array(Apad) #test Concretization of BoundaryPaddedMatrix

        for I in CartesianIndices(A)  #test getindex for all indicies of Apad
            @test A[I] == Apad[I]
        end
    end
end

################################################################################
# Test ComposedBoundaryPaddedMatrix
################################################################################

n = 5
m = 7
A = rand(n,m)
A[1,1] = A[end,1] = A[1,end] = A[end,end] = 0.0

lower = Vector[A[1,2:(end-1)], A[2:(end-1),1]]
upper = Vector[A[end,2:(end-1)], A[2:(end-1),end]]

Apad = DiffEqOperators.ComposedBoundaryPaddedMatrix{Float64, typeof(A), typeof(lower[1])}(lower, upper, A[2:(end-1), 2:(end-1)])

@test A == Array(Apad) #test Concretization of BoundaryPaddedMatrix

for i in 1:n, j in 1:m #test getindex for all indicies of Apad
    @test A[i,j] == Apad[i,j]
end

################################################################################
# Test ComposedBoundaryPaddedTensor{3}
################################################################################

n = 4
m = 5
o = 7
A = rand(n,m,o)
A[1,1,:] = A[end,1, :] = A[1,end, :] = A[end,end, :] = zeros(o)
A[1,:,1] = A[end, :, 1] = A[1,:,end] = A[end,:,end] = zeros(m)
A[:,1,1] = A[:,end,1] = A[:,1,end] = A[:,end,end] = zeros(n)

lower = Matrix[A[1,2:(end-1), 2:(end-1)], A[2:(end-1),1, 2:(end-1)], A[2:(end-1), 2:(end-1), 1]]
upper = Matrix[A[end, 2:(end-1), 2:(end-1)], A[2:(end-1), end, 2:(end-1)], A[2:(end-1), 2:(end-1), end]]

Apad = DiffEqOperators.ComposedBoundaryPadded3Tensor{Float64, typeof(A), typeof(lower[1])}(lower, upper, A[2:(end-1), 2:(end-1), 2:(end-1)])

@test A == Array(Apad) #test Concretization of BoundaryPaddedMatrix

for I in CartesianIndices(A) #test getindex for all indicies of Apad
    @test A[I] == Apad[I]
end


text:
# Ref FluxML/Flux.jl#1209
@testset "BPTT" begin
  seq = [rand(Float32, (2,1)) for i = 1:3]
  for r ∈ [RNN,]
    rnn = r(2,3)
    Flux.reset!(rnn)
    grads_seq = gradient(Flux.params(rnn)) do
      sum(rnn.(seq)[3])
    end
    Flux.reset!(rnn);
    bptt = gradient(Wh->sum(tanh.(rnn.cell.Wi * seq[3] + Wh *
                                  tanh.(rnn.cell.Wi * seq[2] + Wh *
                                        tanh.(rnn.cell.Wi * seq[1] +
                                            Wh * rnn.cell.state
                                        + rnn.cell.b)
                                  + rnn.cell.b)
                            + rnn.cell.b)),
                    rnn.cell.Wh)
    @test grads_seq[rnn.cell.Wh] ≈ bptt[1]
  end
end


text:
# function black_formula(optType::OptionType, strike::Float64, forward::Float64, stdDev::Float64,
#                       discount::Float64 = 1.0, displacement::Float64 = 0.0)
#   # TODO check parameters
#   if stdDev == 0.0
#     return max((forward - strike) * value(optType), 0.0) * discount
#   end
#
#   forward = forward + displacement
#   strike = strike + displacement
#
#   if strike == 0.0
#     return isa(optType, Call) ? forward * discount : 0.0
#   end
#
#   d1 = log(forward / strike) / stdDev + 0.5 * stdDev
#   d2 = d1 - stdDev
#   phi = Normal()
#   nd1 = phi(value(optionType) * d1)
#   nd2 = phi(value(optionType) * d2)
#   result = discount * optionType * (forward * nd1 - strike * nd2)
#   result >= 0.0 || error("negative value")
#
#   return result
# end

function black_scholes_theta(process::AbstractBlackScholesProcess, val::Float64, delta_::Float64, gamma_::Float64)
  u = state_variable(process).value
  r = zero_rate(process.riskFreeRate, 0.0, ContinuousCompounding()).rate
  q = zero_rate(process.dividendYield, 0.0, ContinuousCompounding()).rate
  v = local_vol(process.localVolatility, 0.0, u)

  return r * val - (r - q) * u * delta_ - 0.5 * v * v * u * u * gamma_
end

function black_formula_std_dev_derivative(strike::Float64, forward::Float64, stdDev::Float64, expiry::Float64, disc::Float64, displacement::Float64)
  # TOOD check params
  stdDev >= 0.0 || error("stdDev must be non-negative")
  disc > 0.0 || error("discount must be positive")

  forward = forward + displacement
  strike = strike + displacement

  if stdDev == 0.0 || strike == 0.0
    return 0.0
  end

  d1 = log(forward/strike) / stdDev + 0.5 * stdDev

  return disc * forward * distribution_derivative(Normal(), d1)
end

black_formula_vol_derivative(strike::Float64, forward::Float64, stdDev::Float64, expiry::Float64, disc::Float64, displacement::Float64) =
                            black_formula_std_dev_derivative(strike, forward, stdDev, expiry, disc, displacement) * sqrt(expiry)


text:
module RareEventABC

# package code goes here

using Distributions, DataFrames, ProgressMeter
import StatsFuns.logsumexp
import Distributions: logpdf, std, quantile, rand, length, _rand!, _logpdf

export doSMC, SMCoptions, doPMMH       ##Algorithms
export plotEpidemic                    ##Application specific
export ess, MvIndependent              ##Miscellaneous

include("delta.jl")
include("mvIndependent.jl")
include("misc.jl")
include("slicesample.jl")
include("smc.jl")
include("iidnormal.jl")
include("epidemics.jl")
include("pmmh.jl")

end # module


text:
using JetPack, Jets, Test

n1,n2 = 101,201

@testset "JopMix, linearity test, T=$T, mix=$mix" for T in (Float32,Float64), mix in ((4,3), (0,3), (3,0))
    A = JopMix(JetSpace(T,n1,n2), mix)
    lhs,rhs = linearity_test(A)
    @test lhs ≈ rhs
end

@testset "JopMix, dot product test, T=$T, mix=$mix" for T in (Float32,Float64), mix in ((4,3), (0,3), (3,0))
    A = JopMix(JetSpace(T,n1,n2), mix)
    lhs,rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs
end