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

using Markdown
using InteractiveUtils

# ╔═╡ dcfde69e-20c7-11eb-0108-dbd2f5cb0098
begin
	using Pkg
    Pkg.activate(pwd())
	Pkg.instantiate()
end

# ╔═╡ 864115a0-20d4-11eb-16ff-71d198f9342e
begin
	using CSV
	using DataFrames
	using RDatasets
	using StatsBase
	using Plots
	using LinearAlgebra
	import Statistics: covm
end

# ╔═╡ a33c14a0-20c2-11eb-0bf8-837173d9d91c
md"""
# Assignment 1: Implement PCA
### Aadam (CS1945)
### CSE 662 - Deep Learning by Dr. Masroor Hussain

---

## Assignment Statement:
Write the code in C, Python or Matlab for $l=1$ PCA as discussed in class of Chapter no. 2. Five bonus points for $1 =< l < n$

---
"""

# ╔═╡ cf5c9eb0-20c2-11eb-05ba-2f5dd23fdaf4
md"""
# Principal Component Analysis (PCA)
PCA is a technique that is used to derive an orthogonal projection to convert a given set of observations to linearly uncorrelated variables, called **principal components**.

One of the most widely used application of PCA is **dimensionality reduction**, which helps us out, when training our models on large datasets, by minimizing noise and redundancy in the dataset. Basic idea of the dimensionality reduction is to represent an $M$-dimensional data into an $N$-dimensional subspace, where $N < M$. This would result in data loss, but our aim is to minimize it as much as possible. For that purposes, we try to find some **principal components** in the data that can represent the features as a linear combination, without loosing much information.

PCA uses **covariance matrix** to analyze:
- variance of each feature, showing if a feature is relevant or pure noise.
- linear relationship strength between pairs of features, spotting redundant features.

There are two main approaches to identify principal components:
1. Calculate the **eigenvectors** of the covariance matrix.
2. Calculate the **single value decomposition** of the covariance matrix.

Although **SVD** has higher numerical accuracy, it has lower running time, as compared to the **eigenvector decomposition** method.

For this assignment, we'll use **Eigendecomposition** method to find principal components.
"""

# ╔═╡ b6ed7820-20ce-11eb-0584-9130387e39ef
md"""
# Solution:
First, let's create a dummy dataset.
"""

# ╔═╡ 85bf8990-20d4-11eb-1a7b-e35d9e62f35a
gr();

# ╔═╡ f5d95832-20ea-11eb-2828-39ebb61c877f
md"""
### Utility functions
"""

# ╔═╡ 0024d7b0-20eb-11eb-2959-af3584eb7750
begin
	centralize(x::AbstractVector, m::AbstractVector) = (isempty(m) ? x : x - m)
	centralize(x::AbstractMatrix, m::AbstractVector) = (isempty(m) ? x : x .- m)

	decentralize(x::AbstractVector, m::AbstractVector) = (isempty(m) ? x : x + m)
	decentralize(x::AbstractMatrix, m::AbstractVector) = (isempty(m) ? x : x .+ m)

	transform(P, x, mv) = transpose(P) * centralize(x, mv)
	reconstruct(P, y, mv) = decentralize(P * y, mv)

	# extract k values/vectors
	function extract_kv(fac, ord, k)
		si = ord[1:k]
		vals = fac.values[si]
		vecs = fac.vectors[:, si]
		return (vals, vecs)
	end
end

# ╔═╡ 85976820-20d4-11eb-2afa-47a833fdccc4
md"""
## Input Data
For this example, we'll use **iris** dataset.
"""

# ╔═╡ dcce4b20-20c7-11eb-3cb1-852a46b75c31
iris = dataset("datasets", "iris") # DataFrame(rand(100, 10))

# ╔═╡ dcb4a8a2-20c7-11eb-3a1a-9fb620e84e57
size(iris)

# ╔═╡ dc975ca0-20c7-11eb-0a2b-b9e470c6e311
describe(iris)

# ╔═╡ 9241d120-20d2-11eb-0b30-25249995c00c
md"""
#### Let's create separate train and test datasets.
"""

# ╔═╡ b2a8ccc0-20d2-11eb-28eb-87dc018ae352
Xtrain = Array(iris[1:2:end,1:4])';

# ╔═╡ afb22c8e-20d3-11eb-339a-09f39705a02f
Ytrain = Array(iris[1:2:end,5]);

# ╔═╡ afe6aa10-20d3-11eb-03a4-89a46ade9feb
Xtest = Array(iris[2:2:end,1:4])';

# ╔═╡ b0004c90-20d3-11eb-2b38-452b177600d8
Ytest = Array(iris[2:2:end,5]);

# ╔═╡ b02fe810-20d3-11eb-21f4-f9939c117c12
md"""
Suppose `Xtrain` and `Xtest` are training and testing data matrix, with each observation in a column.

## Data Normalization
Now, we need to normalize the data.

"""

# ╔═╡ 399f2a20-20d4-11eb-1c3a-e76ac9707e29
mean_iris = vec(mean(Xtrain, dims=2))

# ╔═╡ c58d8a10-20e1-11eb-301a-eb715e5926d8
md"""
## Covariance Matrix Computation
Covariance matrix is an $N \times N$ symmetric matrix, where $N$ is the number of dimensions. This tells us about the relation between any two features, and the variation of data from the mean.
"""

# ╔═╡ 39b65ba2-20d4-11eb-11bd-7daa396810ac
C_iris = covm(Xtrain, mean_iris, 2);

# ╔═╡ 6493e090-20e3-11eb-0dec-0bdef8eb1c6a
DataFrame(C_iris)

# ╔═╡ c308e590-20e2-11eb-0ddb-055f8bdedb85
md"""
## Eigendecomposition of Covariance Matrix
We need to find the eigenvalues and eigenvectors of the covariance matrix in order to find the principal components of the data. These components are created in such a way that they are uncorrelated, and the first component contains the maximum amount of information possible, and the second contains the maximum remaining, and so on.
"""

# ╔═╡ 39e24da0-20d4-11eb-10e5-8320e876fb74
eg = eigen(Symmetric(C_iris));

# ╔═╡ b04fa510-20d3-11eb-0d40-21f2bf50e115
eg_values = eg.values

# ╔═╡ f010e210-20e5-11eb-1952-37e6ed84ea22
eg_vectors = eg.vectors

# ╔═╡ 35710110-20e5-11eb-0ad2-ef911c8d78bc
md"""
Let's sort the eigenvalues in reverse order, so the component which contains the most information is at the front.
"""

# ╔═╡ b1d8dc20-20de-11eb-1364-59fc5f488736
ord = sortperm(eg_values; rev=true)

# ╔═╡ 6f9d6dfe-20e5-11eb-2745-e95f930eaf88
md"""
Let's see how much information each component contains. For that, we need to find the total variance of the input, and then divide the principal component's variance with it. We'll plot the values for better visualization.
"""

# ╔═╡ b211ed80-20de-11eb-3abc-7f0b887915a5
tvar = sum(eg_values)

# ╔═╡ 9fefec80-20e6-11eb-3e16-e5f17e1d9560
eg_values[ord] ./ tvar * 100

# ╔═╡ b25a1a10-20de-11eb-3bcb-fbe70a5b974d
bar(eg_values[ord] ./ tvar * 100, title="Scree Plot")

# ╔═╡ b26da210-20de-11eb-026b-8fc4efce305c
md"""
By looking at the above graph, we can see that the first component, $PC 1$, contains more than $90 \%$ of the information, while $PC 2$ contains less than $5 \%$.

Even though many believe $PCA$ to be a dimensionality reduction algorithm, primarily it's a data transformation algorithm. It just makes our data amenable to data reduction, as shown above. Here, as most of the information is contained in the first couple of principal components, we can decide to keep only those componenets and discard the rest, resulting in a reduction in the dimensionality of the dataset.

### $l = 1$:
This means that we only take the first PC and discard the rest.
"""

# ╔═╡ d6429920-20cf-11eb-0bf1-8f74c17dc643
v, P = extract_kv(eg, ord, 1)

# ╔═╡ 72f60020-20eb-11eb-049a-d14c26e12ddb
md"""
This returns the eigen values, as well as the project matrix $P$. We can use $P$ to transform dataset. Transforming the `Xtest` using $P$ would result in:
"""

# ╔═╡ d628f6a0-20cf-11eb-11cc-01f1ca2b6f13
tr_data = transform(P, Xtest, mean_iris)

# ╔═╡ 6ae889c0-20eb-11eb-3ecb-4f13647511c3
md"""
which is a $1 \times 75$ array, instead of the original $4 \times 75$.

Let's try to reconstruct the original data from the transformed data using $P$.
"""

# ╔═╡ 1081c390-20f8-11eb-3e81-133df7444871
recon_data = reconstruct(P, tr_data, mean_iris);

# ╔═╡ 3313dc90-20f8-11eb-27eb-e9ea7e6a3805
md"""
Printing the first 5 observations from both the reconstructed data as well as the original data.
"""

# ╔═╡ d5e95590-20cf-11eb-0d5f-e3e5f2be14c0
DataFrame(recon_data)

# ╔═╡ d5b26710-20cf-11eb-0f4c-7f81055af3e3
DataFrame(Xtest[:, 1:5])

# ╔═╡ 4cd4fb30-20c3-11eb-0631-a9c8840a4cf9
md"""
As can be seen, the reconstructed data is fairly close to the original data, although not exactly the same because this is a lossy method. Even though we retained around $92\%$ information, we still lost some information by reducing the dimensions.

We can retain more information by increasing $l$, i.e. increasing the number of output dimensions.

## $1 <= l < n$
Instead of taking only $PC 1$, we will take the first three PCs.
"""

# ╔═╡ 740d6e20-20ec-11eb-1b12-1d614d3efb08
v3, P3 = extract_kv(eg, ord, 3)

# ╔═╡ 7466b1b2-20ec-11eb-1a71-e7061e6a7ad6
md"""
Again, let's first transform and reconstruct the data using `P3`.
"""

# ╔═╡ 747f1bb0-20ec-11eb-3cbf-a3a8d7a64563
tr3_data = transform(P3, Xtest, mean_iris);

# ╔═╡ 65ce7a10-20fc-11eb-1d02-f57de17db24e
size(tr3_data)

# ╔═╡ 7493b520-20ec-11eb-0e19-57cbd43d09af
recon3_data = reconstruct(P3, tr3_data, mean_iris);

# ╔═╡ 74a9ae20-20ec-11eb-1e5d-7984c0a399f3
DataFrame(recon3_data)

# ╔═╡ ff1b4e90-20f8-11eb-3924-2fa0c61677fd
DataFrame(Xtest[:, 1:5])

# ╔═╡ 053f9290-20f9-11eb-1cb6-6f5c7b26e279
md"""
We can see that using $l = 3$ gives a better approximation of the original data as compared to $l = 1$.

## Visualization
Let's visualize the first 3 principal components in a 3D plot.
"""

# ╔═╡ 9599eb10-20f9-11eb-044b-cd71da0c30a3
begin
	# group by labels, for color coding
	setosa = tr3_data[:,Ytest.=="setosa"];
	versicolor = tr3_data[:,Ytest.=="versicolor"];
	virginica = tr3_data[:,Ytest.=="virginica"];

	p = scatter(setosa[1,:],setosa[2,:],setosa[3,:],marker=:circle,linewidth=0)
	scatter!(versicolor[1,:],versicolor[2,:],versicolor[3,:],marker=:circle,linewidth=0)
	scatter!(virginica[1,:],virginica[2,:],virginica[3,:],marker=:circle,linewidth=0)
	plot!(p,xlabel="PC1",ylabel="PC2",zlabel="PC3")
end

# ╔═╡ 645934a0-20fb-11eb-1750-cf7ca28ab9e6
md"""
## Conclusion
In this notebook, we showed how we can use `PCA` to find principal components of some data, and then use those for dimensionality reduction without severly affecting and losing the information.
"""

# ╔═╡ 4d526890-20c3-11eb-082b-e912d7c9b3e2
md"""
# Resources

- [StatQuest: Principal Component Analysis (PCA), Step-by-Step](https://www.youtube.com/watch?v=FgakZw6K1QQ&t=458s&ab_channel=StatQuestwithJoshStarmer)
- [Data Analysis 6: Principal Component Analysis (PCA) - Computerphile](https://www.youtube.com/watch?v=TJdH6rPA-TI&ab_channel=Computerphile)
- [Making sense of principal component analysis, eigenvectors & eigenvalues](https://stats.stackexchange.com/questions/2691/making-sense-of-principal-component-analysis-eigenvectors-eigenvalues/140579)
- [A STEP BY STEP EXPLANATION OF PRINCIPAL COMPONENT ANALYSIS](https://builtin.com/data-science/step-step-explanation-principal-component-analysis)
- [Machine Learning — Singular Value Decomposition (SVD) & Principal Component Analysis (PCA)](https://jonathan-hui.medium.com/machine-learning-singular-value-decomposition-svd-principal-component-analysis-pca-1d45e885e491)
- [Relationship between SVD and PCA. How to use SVD to perform PCA?](https://stats.stackexchange.com/questions/134282/relationship-between-svd-and-pca-how-to-use-svd-to-perform-pca)
"""

# ╔═╡ Cell order:
# ╟─a33c14a0-20c2-11eb-0bf8-837173d9d91c
# ╟─cf5c9eb0-20c2-11eb-05ba-2f5dd23fdaf4
# ╟─dcfde69e-20c7-11eb-0108-dbd2f5cb0098
# ╟─b6ed7820-20ce-11eb-0584-9130387e39ef
# ╠═864115a0-20d4-11eb-16ff-71d198f9342e
# ╟─85bf8990-20d4-11eb-1a7b-e35d9e62f35a
# ╟─f5d95832-20ea-11eb-2828-39ebb61c877f
# ╠═0024d7b0-20eb-11eb-2959-af3584eb7750
# ╟─85976820-20d4-11eb-2afa-47a833fdccc4
# ╠═dcce4b20-20c7-11eb-3cb1-852a46b75c31
# ╠═dcb4a8a2-20c7-11eb-3a1a-9fb620e84e57
# ╠═dc975ca0-20c7-11eb-0a2b-b9e470c6e311
# ╟─9241d120-20d2-11eb-0b30-25249995c00c
# ╠═b2a8ccc0-20d2-11eb-28eb-87dc018ae352
# ╠═afb22c8e-20d3-11eb-339a-09f39705a02f
# ╠═afe6aa10-20d3-11eb-03a4-89a46ade9feb
# ╠═b0004c90-20d3-11eb-2b38-452b177600d8
# ╟─b02fe810-20d3-11eb-21f4-f9939c117c12
# ╠═399f2a20-20d4-11eb-1c3a-e76ac9707e29
# ╟─c58d8a10-20e1-11eb-301a-eb715e5926d8
# ╠═39b65ba2-20d4-11eb-11bd-7daa396810ac
# ╟─6493e090-20e3-11eb-0dec-0bdef8eb1c6a
# ╟─c308e590-20e2-11eb-0ddb-055f8bdedb85
# ╠═39e24da0-20d4-11eb-10e5-8320e876fb74
# ╠═b04fa510-20d3-11eb-0d40-21f2bf50e115
# ╠═f010e210-20e5-11eb-1952-37e6ed84ea22
# ╟─35710110-20e5-11eb-0ad2-ef911c8d78bc
# ╠═b1d8dc20-20de-11eb-1364-59fc5f488736
# ╟─6f9d6dfe-20e5-11eb-2745-e95f930eaf88
# ╠═b211ed80-20de-11eb-3abc-7f0b887915a5
# ╟─9fefec80-20e6-11eb-3e16-e5f17e1d9560
# ╠═b25a1a10-20de-11eb-3bcb-fbe70a5b974d
# ╟─b26da210-20de-11eb-026b-8fc4efce305c
# ╠═d6429920-20cf-11eb-0bf1-8f74c17dc643
# ╟─72f60020-20eb-11eb-049a-d14c26e12ddb
# ╠═d628f6a0-20cf-11eb-11cc-01f1ca2b6f13
# ╟─6ae889c0-20eb-11eb-3ecb-4f13647511c3
# ╠═1081c390-20f8-11eb-3e81-133df7444871
# ╟─3313dc90-20f8-11eb-27eb-e9ea7e6a3805
# ╠═d5e95590-20cf-11eb-0d5f-e3e5f2be14c0
# ╠═d5b26710-20cf-11eb-0f4c-7f81055af3e3
# ╟─4cd4fb30-20c3-11eb-0631-a9c8840a4cf9
# ╠═740d6e20-20ec-11eb-1b12-1d614d3efb08
# ╟─7466b1b2-20ec-11eb-1a71-e7061e6a7ad6
# ╠═747f1bb0-20ec-11eb-3cbf-a3a8d7a64563
# ╠═65ce7a10-20fc-11eb-1d02-f57de17db24e
# ╠═7493b520-20ec-11eb-0e19-57cbd43d09af
# ╠═74a9ae20-20ec-11eb-1e5d-7984c0a399f3
# ╠═ff1b4e90-20f8-11eb-3924-2fa0c61677fd
# ╟─053f9290-20f9-11eb-1cb6-6f5c7b26e279
# ╟─9599eb10-20f9-11eb-044b-cd71da0c30a3
# ╟─645934a0-20fb-11eb-1750-cf7ca28ab9e6
# ╟─4d526890-20c3-11eb-082b-e912d7c9b3e2


text:
"""
    KrylovKit

A Julia package collecting a number of Krylov-based algorithms for linear problems,
singular value and eigenvalue problems and the application of functions of linear maps or
operators to vectors.

KrylovKit accepts general functions or callable objects as linear maps, and general Julia
objects with vector like behavior as vectors.

The high level interface of KrylovKit is provided by the following functions:
*   [`linsolve`](@ref): solve linear systems
*   [`eigsolve`](@ref): find a few eigenvalues and corresponding eigenvectors
*   [`geneigsolve`](@ref): find a few generalized eigenvalues and corresponding vectors
*   [`svdsolve`](@ref): find a few singular values and corresponding left and right
    singular vectors
*   [`exponentiate`](@ref): apply the exponential of a linear map to a vector
"""
module KrylovKit

using LinearAlgebra
using Printf
const IndexRange = AbstractRange{Int}

export linsolve, eigsolve, geneigsolve, svdsolve, schursolve, exponentiate, expintegrator
export orthogonalize, orthogonalize!, orthonormalize, orthonormalize!
export basis, rayleighquotient, residual, normres, rayleighextension
export initialize, initialize!, expand!, shrink!
export ClassicalGramSchmidt, ClassicalGramSchmidt2, ClassicalGramSchmidtIR
export ModifiedGramSchmidt, ModifiedGramSchmidt2, ModifiedGramSchmidtIR
export LanczosIterator, ArnoldiIterator, GKLIterator
export CG, GMRES, BiCGStab, Lanczos, Arnoldi, GKL, GolubYe
export KrylovDefaults, ClosestTo, EigSorter
export RecursiveVec, InnerProductVec

include("algorithms.jl")

# Structures to store a list of basis vectors
"""
    abstract type Basis{T} end

An abstract type to collect specific types for representing a basis of vectors of type `T`.

Implementations of `Basis{T}` behave in many ways like `Vector{T}` and should have a
`length`, can be indexed (`getindex` and `setindex!`), iterated over (`iterate`), and
support resizing (`resize!`, `pop!`, `push!`, `empty!`, `sizehint!`).

The type `T` denotes the type of the elements stored in an `Basis{T}` and can be any custom type that has vector like behavior (as defined in the docs of [`KrylovKit`](@ref)).

See [`OrthonormalBasis`](@ref).
"""
abstract type Basis{T} end

include("orthonormal.jl")

# Dense linear algebra structures and functions used in the algorithms below
include("dense/givens.jl")
include("dense/linalg.jl")
include("dense/packedhessenberg.jl")
include("dense/reflector.jl")

# Simple coordinate basis vector, i.e. a vector of all zeros and a single one on position `k`:
"""
    SimpleBasisVector(m, k)

Construct a simple struct `SimpleBasisVector <: AbstractVector{Bool}` representing a
coordinate basis vector of length `m` in the direction of `k`, i.e. for `e_k =
SimpleBasisVector(m, k)` we have `length(e_k) = m` and `e_k[i] = (i == k)`.
"""
struct SimpleBasisVector <: AbstractVector{Bool}
    m::Int
    k::Int
end
Base.size(e::SimpleBasisVector) = (e.m,)
@inline function Base.getindex(e::SimpleBasisVector, i)
    @boundscheck Base.checkbounds(e, i)
    return e.k == i
end

# Krylov factorizations and their iterators, the central objects for writing algorithms in KrylovKit
abstract type KrylovFactorization{T,S} end
abstract type KrylovIterator{F,T} end

"""
        basis(fact::KrylovFactorization)

Return the list of basis vectors of a [`KrylovFactorization`](@ref), which span the Krylov
subspace. The return type is a subtype of `Basis{T}`, where `T` represents the type of the
vectors used by the problem.
"""
function basis end

"""
    rayleighquotient(fact::KrylovFactorization)

Return the Rayleigh quotient of a [`KrylovFactorization`](@ref), i.e. the reduced matrix
within the basis of the Krylov subspace. The return type is a subtype of
`AbstractMatrix{<:Number}`, typically some structured matrix type.
"""
function rayleighquotient end

"""
    residual(fact::KrylovFactorization)

Return the residual of a [`KrylovFactorization`](@ref). The return type is some vector of
the same type as used in the problem. See also [`normres(F)`](@ref) for its norm, which
typically has been computed already.
"""
function residual end

"""
    normres(fact::KrylovFactorization)

Return the norm of the residual of a [`KrylovFactorization`](@ref). As this has typically
already been computed, it is cheaper than (but otherwise equivalent to) `norm(residual(F))`.
"""
function normres end

"""
    rayleighextension(fact::KrylovFactorization)

Return the vector `b` appearing in the definition of a [`KrylovFactorization`](@ref).
"""
function rayleighextension end

"""
    shrink!(fact::KrylovFactorization, k)

Shrink an existing Krylov factorization `fact` down to have length `k`. Does nothing if
`length(fact)<=k`.
"""
function shrink! end

"""
    expand!(iter::KrylovIterator, fact::KrylovFactorization)

Expand the Krylov factorization `fact` by one using the linear map and parameters in `iter`.
"""
function expand! end

"""
    initialize!(iter::KrylovIterator, fact::KrylovFactorization)

Initialize a length 1 Krylov factorization corresponding to `iter` in the already existing
factorization `fact`, thereby destroying all the information it currently holds.
"""
function initialize! end

"""
    initialize(iter::KrylovIterator)

Initialize a length 1 Krylov factorization corresponding to `iter`.
"""
function initialize end

# iteration for destructuring into components
Base.iterate(F::KrylovFactorization) = (basis(F), Val(:rayleighquotient))
Base.iterate(F::KrylovFactorization, ::Val{:rayleighquotient}) = (rayleighquotient(F), Val(:residual))
Base.iterate(F::KrylovFactorization, ::Val{:residual}) = (residual(F), Val(:normres))
Base.iterate(F::KrylovFactorization, ::Val{:normres}) = (normres(F), Val(:rayleighextension))
Base.iterate(F::KrylovFactorization, ::Val{:rayleighextension}) = (rayleighextension(F), Val(:done))
Base.iterate(F::KrylovFactorization, ::Val{:done}) = nothing

# some often used tools
function checkposdef(z)
    r = checkhermitian(z)
    r > 0 ||
        error("operator does not appear to be positive definite: diagonal element $z")
    return r
end
function checkhermitian(z, n = abs(z))
    imag(z) <= sqrt(max(eps(n),eps(one(n)))) ||
        error("operator does not appear to be hermitian: diagonal element $z")
    return real(z)
end

include("krylov/lanczos.jl")
include("krylov/arnoldi.jl")
include("krylov/gkl.jl")

"""
    abstract type KrylovFactorization{T,S<:Number}
    mutable struct LanczosFactorization{T,S<:Real}    <: KrylovFactorization{T,S}
    mutable struct ArnoldiFactorization{T,S<:Number}  <: KrylovFactorization{T,S}

Structures to store a Krylov factorization of a linear map `A` of the form
```julia
    A * V = V * B + r * b'.
```
For a given Krylov factorization `fact` of length `k = length(fact)`, the basis `A is
obtained via [`basis(fact)`](@ref basis) and is an instance of some subtype of
[`Basis{T}`](@ref Basis), with also `length(V) == k` and where `T` denotes the type of
vector like objects used in the problem. The Rayleigh quotient `B` is obtained as
[`rayleighquotient(fact)`](@ref) and `typeof(B)` is some subtype of `AbstractMatrix{S}`
with `size(B) == (k,k)`, typically a structured matrix. The residual `r` is obtained as
[`residual(fact)`](@ref) and is of type `T`. One can also query [`normres(fact)`](@ref) to
obtain `norm(r)`, the norm of the residual. The vector `b` has no dedicated name and often
takes a default form (see below). It should be a subtype of `AbstractVector` of length `k`
and can be obtained as [`rayleighextension(fact)`](@ref) (by lack of a better dedicated
name).

In particular, `LanczosFactorization` stores a Lanczos factorization of a real symmetric or
complex hermitian linear map and has `V::OrthonormalBasis{T}` and
`B::SymTridiagonal{S<:Real}`. `ArnoldiFactorization` stores an Arnoldi factorization of a
general linear map and has `V::OrthonormalBasis{T}` and
[`B::PackedHessenberg{S<:Number}`](@ref PackedHessenberg). In both
cases, `b` takes the default value ``e_k``, i.e. the unit vector of all zeros and a one in
the last entry, which is represented using [`SimpleBasisVector`](@ref).

A Krylov factorization `fact` can be destructured as `V, B, r, nr, b = fact` with `nr = norm(r)`.

`LanczosFactorization` and `ArnoldiFactorization` are mutable because they can
[`expand!`](@ref) or [`shrink!`](@ref). See also [`KrylovIterator`](@ref) (and in
particular [`LanczosIterator`](@ref) and [`ArnoldiIterator`](@ref)) for iterators that
construct progressively expanding Krylov factorizations of a given linear map and a
starting vector.
"""
KrylovFactorization, LanczosFactorization, ArnoldiFactorization

"""
    abstract type KrylovIterator{F,T}
    struct LanczosIterator{F,T,O<:Orthogonalizer} <: KrylovIterator{F,T}
    struct ArnoldiIterator{F,T,O<:Orthogonalizer} <: KrylovIterator{F,T}

    LanczosIterator(f, v₀, [orth::Orthogonalizer = KrylovDefaults.orth], keepvecs::Bool = true)
    ArnoldiIterator(f, v₀, [orth::Orthogonalizer = KrylovDefaults.orth])

Iterators that take a linear map of type `F` and an initial vector of type `T` and generate
an expanding `KrylovFactorization` thereof.

In particular, for a real symmetric or complex hermitian linear map `f`, `LanczosIterator`
uses the [Lanczos iteration](https://en.wikipedia.org/wiki/Lanczos_algorithm) scheme to
build a successively expanding `LanczosFactorization`. While `f` cannot be tested to be
symmetric or hermitian directly when the linear map is encoded as a general callable object
or function, it is tested whether the imaginary part of `dot(v, f(v))` is sufficiently
small to be neglected.

Similarly, for a general linear map `f`, `ArnoldiIterator` iterates over progressively
expanding `ArnoldiFactorizations` using the [Arnoldi
iteration](https://en.wikipedia.org/wiki/Arnoldi_iteration).

The optional argument `orth` specifies which [`Orthogonalizer`](@ref) to be used. The
default value in [`KrylovDefaults`](@ref) is to use [`ModifiedGramSchmidtIR`](@ref), which
possibly uses reorthogonalization steps. For `LanczosIterator`, one can use to discard the
old vectors that span the Krylov subspace by setting the final argument `keepvecs` to
`false`. This, however, is only possible if an `orth` algorithm is used that does not rely
on reorthogonalization, such as `ClassicalGramSchmidt()` or `ModifiedGramSchmidt()`. In
that case, the iterator strictly uses the Lanczos three-term recurrence relation.

When iterating over an instance of `KrylovIterator`, the values being generated are subtypes
of [`KrylovFactorization`](@ref), which can be immediately destructured into a
[`basis`](@ref), [`rayleighquotient`](@ref), [`residual`](@ref), [`normres`](@ref) and
[`rayleighextension`](@ref), for example as
```julia
for (V,B,r,nr,b) in ArnoldiIterator(f, v₀)
    # do something
    nr < tol && break # a typical stopping criterion
end
```
Note, however, that if `keepvecs=false` in `LanczosIterator`, the basis `V` cannot be
extracted. Since the iterators don't know the dimension of the underlying vector space of
objects of type `T`, they keep expanding the Krylov subspace until `normres` falls below
machine precision `eps` for the given `eltype(T)`.

The internal state of `LanczosIterator` and `ArnoldiIterator` is the same as the return
value, i.e. the corresponding `LanczosFactorization` or `ArnoldiFactorization`. However, as
Julia's Base iteration interface (using `Base.iterate`) requires that the state is not
mutated, a `deepcopy` is produced upon every next iteration step.

Instead, you can also mutate the `KrylovFactorization` in place, using the following
interface, e.g. for the same example above
```julia
iterator = ArnoldiIterator(f, v₀)
factorization = initialize(iterator)
while normres(factorization) > tol
    expand!(iterator, factorization)
    V,B,r,nr,b = factorization
    # do something
end
```
Here, [`initialize(::KrylovIterator)`](@ref) produces the first Krylov factorization of
length 1, and `expand!(::KrylovIterator, ::KrylovFactorization)`(@ref) expands the
factorization in place. See also
[`initialize!(::KrylovIterator,::KrylovFactorization)`](@ref) to initialize in an already
existing factorization (most information will be discarded) and
[`shrink!(::KrylovFactorization, k)`](@ref) to shrink an existing factorization down to
length `k`.
"""
KrylovIterator, LanczosIterator, ArnoldiIterator

# A general structure to pass on convergence information
"""
    struct ConvergenceInfo{S,T}
        converged::Int
        residual::T
        normres::S
        numiter::Int
        numops::Int
    end

Used to return information about the solution found by the iterative method.
*   `converged`: the number of solutions that have converged according to an appropriate error
    measure and requested tolerance for the problem. Its value can be zero or one for
    [`linsolve`](@ref) and [`exponentiate`](@ref), or any integer `>= 0` for
    [`eigsolve`](@ref), [`schursolve`](@ref) or [`svdsolve`]().
*   `residual:` the (list of) residual(s) for the problem, or `nothing` for problems without
    the concept of a residual (i.e. `exponentiate`). This is a single vector (of the same
    type as the type of vectors used in the problem) for `linsolve`, or a `Vector` of such
    vectors for `eigsolve`, `schursolve` or `svdsolve`.
*   `normres`: the norm of the residual(s) (in the previous field) or the value of any other
    error measure that is appropriate for the problem. This is a `Real` for `linsolve` and
    `exponentiate`, and a `Vector{<:Real}` for `eigsolve`, `schursolve` and `svdsolve`. The
    number of values in `normres` that are smaller than a predefined tolerance corresponds
    to the number `converged` of solutions that have converged.
*   `numiter`: the number of iterations (sometimes called restarts) used by the algorithm.
*   `numops`: the number of times the linear map or operator was applied
"""
struct ConvergenceInfo{S,T}
    converged::Int # how many vectors have converged, 0 or 1 for linear systems, exponentiate, any integer for eigenvalue problems
    residual::T
    normres::S
    numiter::Int
    numops::Int
end

function Base.show(io::IO, info::ConvergenceInfo)
    print(io, "ConvergenceInfo: ")
    info.converged == 0 && print(io, "no converged values ")
    info.converged == 1 && print(io, "one converged value ")
    info.converged > 1 && print(io, "$(info.converged) converged values ")
    println(io, "after $(info.numiter) iterations and $(info.numops) applications of the linear map;")
    println(io, "norms of residuals are given by $((info.normres...,)).")
end

# eigsolve en schursolve
include("eigsolve/eigsolve.jl")
include("eigsolve/lanczos.jl")
include("eigsolve/arnoldi.jl")
include("eigsolve/geneigsolve.jl")
include("eigsolve/golubye.jl")
include("eigsolve/svdsolve.jl")

# linsolve
include("linsolve/linsolve.jl")
include("linsolve/cg.jl")
include("linsolve/gmres.jl")
include("linsolve/bicgstab.jl")

# exponentiate
include("matrixfun/exponentiate.jl")
include("matrixfun/expintegrator.jl")

apply(A::AbstractMatrix, x::AbstractVector) = A*x
apply(f, x) = f(x)

apply!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, A, x)
apply!(y, f, x) = copyto!(y, f(x))

include("recursivevec.jl")
include("innerproductvec.jl")

end

ext: jl


text:
function gallery()

s01 = raw"""
                  _oo0oo_
                 o8888888o
                 88i . i88
                 (| -_- |)
                 0|  =  |0
               ___|i---i|___
             .i |||     ||| i.
            | |||||  :  ||||| |
           | _||||| -:- |||||- |
          |   | |||  -  ||| |   |
          | |_|  i|---|i  |_| |
          |  .-|__  i-i  ___|-. |
        ___i. .i  |--.--|  i. .i___
     .ii i<  i.___|_<|>_|___.i >i ii.
    | | :  i- |i.;i| _ |i;.i| - i : | |
    |  | i_.   |_ __| |__ _|   .-i |  |
=====i-.____i.___ |_____|___.-i___.-i=====
                  i=---=i
"""

s02 = raw"""
   ┌─┐       ┌─┐
┌──┘ ┴───────┘ ┴──┐
│                 │
│       ───       │
│  ─┬┘       └┬─  │
│                 │
│       ─┴─       │
│                 │
└───┐         ┌───┘
    │         │
    │         │
    │         │
    │         └──────────────┐
    │                        │
    │                        ├─┐
    │                        ┌─┘
    │                        │
    └─┐  ┐  ┌───────┬──┐  ┌──┘
      │ ─┤ ─┤       │ ─┤ ─┤
      └──┴──┘       └──┴──┘
"""

s03 = raw"""
                                  .  .
                  -.            |_|||        |||  ~~~~~~::::... |~
               ___-==_       _-~o~  ||    |||  ||            _|~~-
       __---~~~.==~|||=_    -_--~|_-~|-   |||   ||        _|~
   _-~~     .=~    |  ||-_    i-~7  |-   |  ||    |      |
 .~       .~       |   || -_    |  |-   |   ||      |   |
|  ____  |         |     || ~-_|  ||- _|   .||       | |
|~~    ~~|--~~~~--_ |     ~==-|   | |~--===~~        .|
         i         ~-|      ||    |-~|~~       __--~~
                     |-~~-_| |    |   ~|_   _-~            ||
                          |  |     |__   ||~                |__
                      _--~ _| | .-~~____--~-|                  ~~==.
                     ((->|~   i.|||i -_|    ~~-| ,              . _||
                                -_     ~|      ~~---l__i__i__i--~~_|
                                _-~-__   ~)  |--______________--~~
                              ||.-~~~-~_--~- |-------~~~~~~~~
                                     ||.-~~~--|
"""

s04 = raw"""
                      iiiii.
 iiiiiiii.           Bd    Bo
4B       iB          Ba      B
iB        iB        JB       BF
 ib        Bk       B>       B
  Bk        Br     JB       dB
  iB         B     Bs       B~
   iB        sB    BE       B
    B         BL   Bi      BF
     B.       4B   B      B*sss*b
     iB        B.  B     B      BF
      iB       RB  BF     Bi      B
       Bk      ?B u*     dF      .B
       ^B.      Bi     zB      uBe
        BBb             BE.dW@eBi    ?B
         BB           .oBB dBc    ?F
          B      .dBBi . zoB>   BBr .uF
          BL .uB*i      B&Bk   .BdBF
           Bi            ii^iBPiBP9B
          JP              .oBu:BP B
          B          ..ueBi      ii  Bi
         dB          BF              B
         B     ....udE             4B
          BB    aaaai Br            @B
           ^BL         B            BF
             RN        4N           B
              *Bb                  dB
               Bk                 BF
               Bb                BF
                 Bii               BF
                 sB                B
                  BL               B
                  iB               B
                   B               B
"""

s05 = raw"""
                  ___====-_  _-====___
            _--^^^#####||      ||#####^^^--_
         _-^##########|| (    ) ||##########^-_
        -############||  ||^^||  ||############-
      _|############||   (@::@)   ||############|_
     |#############((     ||||     ))#############|
    -###############||    (oo)    ||###############-
   -#################||  | VV |  ||#################-
  -###################|||      |||###################-
 _#||##########||######(   ||   )######||##########||#_
 || |#||#||#|||  |#||##|  |  |  |##||#|  |||#||#||#| ||
 i  ||  V  V  i   V  |#|| |  | ||#|  V   i  V  V  ||  i
    i   i  i      i   | | |  | | |   i      i  i   i
                     (  | |  | |  )
                    __| | |  | | |__
                   (vvv(VVV)(VVV)vvv)
"""

s06 = raw"""
              ,----------------,              ,---------,
         ,-----------------------,          ,i        ,i|
       ,i                      ,i|        ,i        ,i  |
      +-----------------------+  |      ,i        ,i    |
      |  .-----------------.  |  |     +---------+      |
      |  |                 |  |  |     | -==----i|      |
      |  |  I LOVE DOS!    |  |  |     |         |      |
      |  |  Bad command or |  |  ||----|i---=    |      |
      |  |  C:|>_          |  |  |   ,||==== ooo |      ;
      |  |                 |  |  |  || |((((     |    ,i
      |  i-----------------i  |,i .;i| |((((     |  ,i
      +-----------------------+  ;;  | |         |,i
         |_)______________(_|  ||i   | +---------+
    ___________________________|___  i,
   |  oooooooooooooooo  .o.  oooo |,   |,i-----------
  | ==ooooooooooooooo==.o.  ooo= ||   ,i|--{)B     ,i
 |_==__==========__==_ooo__ooo=_|i   |___________,i
"""

s07 = raw"""
             ..                     .''.
            c  '.          _,,    .'    '
           ,     ,   ___  7###. .'       '
          .       `-'   ''|####;#'. .##, .
          ,               '####.##,######.
         .                 '  ' '' ,####.
        '                           '##'.
       '                                .
__,,--.--                                '
      '                                  '
  .--|-        &&                        -'-'''-.._
     |         '&              &&         '
    __,-                       '&        -'---.
  ''  ,                                   ;
       .             kk.                 -.-,_
        -_           ''KK                ;
          ''--,_                        ,
                ''--,__           __,,--'
                       ''--..,--''     Pr59
"""

    return Items(s01, s02, s03, s04, s05, s06, s07)
end

function highlight(content::String)
    panel = Panel(
        content,
        justify          = :left,
        style            = "gold1",
        title            = "🌺 Welcome to Art Gallery 🌺",
        title_style      = "bold cyan centre",
        title_justify    = :center,
        subtitle         = string(today()),
        subtitle_style   = "bold cyan",
        subtitle_justify = :right,
        width            = 79
    )

    return print(panel)
end


text:
# ---
# title: 1674. Minimum Moves to Make Array Complementary
# id: problem1674
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Greedy
# link: <https://leetcode.com/problems/minimum-moves-to-make-array-complementary/description/>
# hidden: true
# ---
#
# You are given an integer array `nums` of **even** length `n` and an integer
# `limit`. In one move, you can replace any integer from `nums` with another
# integer between `1` and `limit`, inclusive.
#
# The array `nums` is **complementary** if for all indices `i` ( **0-indexed**
# ), `nums[i] + nums[n - 1 - i]` equals the same number. For example, the array
# `[1,2,3,4]` is complementary because for all indices `i`, `nums[i] + nums[n -
# 1 - i] = 5`.
#
# Return the _**minimum** number of moves required to make _`nums`
# _**complementary**_.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,2,4,3], limit = 4
#     Output: 1
#     Explanation: In 1 move, you can change nums to [1,2, _2_ ,3] (underlined elements are changed).
#     nums[0] + nums[3] = 1 + 3 = 4.
#     nums[1] + nums[2] = 2 + 2 = 4.
#     nums[2] + nums[1] = 2 + 2 = 4.
#     nums[3] + nums[0] = 3 + 1 = 4.
#     Therefore, nums[i] + nums[n-1-i] = 4 for every i, so nums is complementary.
#
#
# **Example 2:**
#
#
#
#     Input: nums = [1,2,2,1], limit = 2
#     Output: 2
#     Explanation: In 2 moves, you can change nums to [ _2_ ,2,2, _2_ ]. You cannot change any number to 3 since 3 > limit.
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1,2,1,2], limit = 2
#     Output: 0
#     Explanation: nums is already complementary.
#
#
#
#
# **Constraints:**
#
#   * `n == nums.length`
#   * `2 <= n <= 105`
#   * `1 <= nums[i] <= limit <= 105`
#   * `n` is even.
#
#
## @lc code=start
using LeetCode

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


text:
"""
 Gets the current slack of the constraint, for feasible solution it's always positive.
 For double sided inequalities it's the least slack that is given.
"""
function getslack(constraint::ConstraintRef)
  i = constraint.idx
  row_val = MathProgBase.getconstrsolution(internalmodel(constraint.m))[i]
  row_lb = MathProgBase.getconstrLB(internalmodel(constraint.m))[i]
  row_ub = MathProgBase.getconstrUB(internalmodel(constraint.m))[i]
  row_basis = MathProgBase.getbasis(internalmodel(constraint.m))[2][i]
  if row_basis != :Basic
    return 0.0
  end
  return min(row_ub - row_val, row_val - row_lb)
end

"""
  Computes the perturbation range for costs and rhs, in which the current basis remains optimal.
"""
function getoptranges(model::Model)
  col_basis, row_basis = MathProgBase.getbasis(internalmodel(model))
  cost_int = zeros(model.numCols,2)  # intervals for non-basic costs
  cost_int[col_basis .== :NonbasicAtLower, 1] = -model.redCosts[col_basis .== :NonbasicAtLower]
  cost_int[col_basis .== :NonbasicAtLower, 2] .= Inf

  cost_int[col_basis .== :NonbasicAtUpper, 1] .= -Inf
  cost_int[col_basis .== :NonbasicAtUpper, 2] = model.redCosts[col_basis .== :NonbasicAtUpper]

  if model.objSense == :Max
    cost_int[:,1], cost_int[:,2] = cost_int[:,2], cost_int[:,1]
    cost_int[isinf.(cost_int)] *= -1.0;
  end
  # intervals for basic cost costcoefficient
  #c_N - c_B B^-1 N >= 0 => c_N - pi N >= dc_B B^-1 N
  #
  A = MathProgBase.getconstrmatrix(internalmodel(model))
  Ident = spdiagm(0=>ones(A.m))
  B = hcat(A[:,col_basis .== :Basic], Ident[:,row_basis .== :Basic])
  N = hcat(A[:,col_basis .== :NonbasicAtLower], A[:,col_basis .== :NonbasicAtUpper],
    Ident[:,row_basis .== :NonbasicAtLower], Ident[:,row_basis .== :NonbasicAtUpper])
  y = model.linconstrDuals

  c_red = vcat(model.redCosts[col_basis .== :NonbasicAtLower],
   model.redCosts[col_basis .== :NonbasicAtUpper],
   -y[row_basis .== :NonbasicAtLower], -y[row_basis .== :NonbasicAtUpper])

  cost_int[col_basis .== :Basic, 1] .= -Inf
  cost_int[col_basis .== :Basic, 2] .= Inf
  for i = 1:sum(col_basis .== :Basic)
    j = (1:model.numCols)[col_basis .== :Basic][i]
    dcb = zeros(A.m)
    dcb[i] = 1.0
    A_red = ((B'\dcb)'*N)'
    neg_ind = ((A_red .< -1e-6) .& (c_red .> 1e-8)) .|  ((A_red .> 1e-6) .& (c_red .< -1e-8))
    cost_int[j, 1] = maximum([cost_int[j, 1]; c_red[neg_ind]./A_red[neg_ind]])
    pos_ind = ((A_red .> 1e-6) .& (c_red .> 1e-8)) .|  ((A_red .< -1e-6) .& (c_red .< -1e-8))
    cost_int[j, 2] = minimum([cost_int[j, 2]; c_red[pos_ind]./A_red[pos_ind]])
  end
  #Intervals for rhs: U >= B^-1(b+td) >= L  =>  t B^-1d >= L-B^-1b && t B^-1d <= U -B^-1b
  B = hcat(A[:,col_basis .== :Basic], -Ident[:,row_basis .== :Basic])
  x_B = model.colVal[col_basis .== :Basic]
  x_lb = model.colLower[col_basis .== :Basic]
  x_ub = model.colUpper[col_basis .== :Basic]
  row_val = MathProgBase.getconstrsolution(internalmodel(model))[row_basis .== :Basic]
  row_lb = MathProgBase.getconstrLB(internalmodel(model))[row_basis .== :Basic]
  row_ub = MathProgBase.getconstrUB(internalmodel(model))[row_basis .== :Basic]
  x_lb_d = [x_lb - x_B; row_lb - row_val]
  x_ub_d = [x_ub - x_B; row_ub - row_val]

  rhs_int = zeros(A.m,2)
  rhs_int[:,1] .= -Inf
  rhs_int[:,2] .= Inf
  rhs_int[row_basis .== :Basic, 1] = row_val - row_ub
  rhs_int[row_basis .== :Basic, 2] = row_val - row_lb
  for i = (1:A.m)[row_basis .!= :Basic]
    ei = zeros(A.m)
    ei[i] = 1.0
    db = B\ei
    neg_ind = db .< -1e-7
    rhs_int[i, 1] = maximum([rhs_int[i, 1]; x_ub_d[neg_ind]./db[neg_ind]])
    rhs_int[i, 2] = minimum([rhs_int[i, 2]; x_lb_d[neg_ind]./db[neg_ind]])
    pos_ind = db .> 1e-7
    rhs_int[i, 1] = maximum([rhs_int[i, 1]; x_lb_d[pos_ind]./db[pos_ind]])
    rhs_int[i, 2] = minimum([rhs_int[i, 2]; x_ub_d[pos_ind]./db[pos_ind]])
  end
  return cost_int, rhs_int;
end
"""
  Computes the perturbation range for var, in which the current basis remains optimal.
"""
function getoptrange(var::Variable)
  dc, db = getoptranges(var.m)
  return dc[var.col,:]
end

"""
  Computes the perturbation range for constraint, in which the current basis remains optimal.
"""
function getoptrange(constraint::ConstraintRef)
  dc, db = getoptranges(constraint.m)
  return db[constraint.idx,:]
end

max_stars_repo_name: EdvinAblad/JuMPLPUtils

max_issues_repo_name: EdvinAblad/JuMPLPUtils

max_forks_repo_name: EdvinAblad/JuMPLPUtils


text:
function bat_fit(hist, fit_function, prior; nsamples=10^5, nchains=4, params_first=true, norm_expected=true)
    likelihood = let h = hist, f = fit_function
        # Histogram counts for each bin as an array:
        observed_counts = h.weights

        # Histogram binning:
        bin_edges = h.edges[1]
        bin_edges_left = bin_edges[1:end-1]
        bin_edges_right = bin_edges[2:end]
        bin_widths = bin_edges_right - bin_edges_left
        bin_centers = midpoints(bin_edges)

        params -> begin
            if params_first
                expected_counts = bin_widths .* f(params, bin_centers)
            else
                expected_counts = bin_widths .* f(bin_centers, params)
            end
            if norm_expected
                expected_counts ./= sum(expected_counts)
                expected_counts .*= sum(observed_counts)
            end
            ll_value = sum(logpdf.(Poisson.(expected_counts), observed_counts))
            return LogDVal(ll_value)
        end
    end
    parshapes = varshape(prior)
    posterior = PosteriorDensity(likelihood, prior);
    return bat_sample(posterior, MCMCSampling(mcalg = MetropolisHastings(), nsteps = nsamples, nchains = nchains)).result;
end

#
# Wrapper for fits that might need several iterations
#
function bat_fit_wrapper(histogram, model, prior, nsamples_list; ch_str="XX", params_first=false, norm_expected=false)
    parameters, errors = [], []
    run = 1
    for ns in nsamples_list
        @info("Ch" * ch_str * " | try " * string(run) * " of " * string(length(nsamples_list)) * ": " * string(ns) * " samples")
        samples = try bat_fit(histogram, model, prior; nsamples=ns, nchains=4, params_first=params_first, norm_expected=norm_expected);
        catch
            IJulia.clear_output(true)
            run += 1
            continue
        end
        IJulia.clear_output(true)
        @info("Success: Ch" * ch_str * " | " * string(ns) * " samples")
        parameters = mode(samples)
        errors = std(samples)
        break
    end
    return [parameters, errors]
end

#
# Peak model with low E tail
#
function peak(E, par)
    n = try par[1] catch; par[1][1] end
    σ = try par[2] catch; par[2][1] end
    µ = try par[3] catch; par[3][1] end

    bkg_l = try par[4] catch; par[4][1] end
    bkg_r = try par[5] catch; par[5][1] end
    b = (bkg_r - bkg_l) / (maximum(E) - minimum(E))
    a = bkg_l - b * minimum(E)

    c = try par[6] catch; par[6][1] end
    d = try par[7] catch; par[7][1] end
    β = try par[8] catch; par[8][1] end

    d = b > 0 ? -d : d
    correction = d < 0 ? -1*minimum(d/2 .* erfc.((E .- µ) ./ (sqrt(2)*σ))) : 0.0
    result = @. n / sqrt(2*pi*σ^2) * exp(-(E - µ)^2 / (2*σ^2)) + # Gaussian
        a + b*E + # Linear bkg
        d/2 * erfc((E - µ) / (sqrt(2)*σ)) + correction +  # Low-E step
        c/(2*β) * exp((E-µ) / β + σ^2 / (2*β^2)) * erfc((E - µ) / (sqrt(2)*σ) + σ / (sqrt(2)*β)) # Low-E tail
end
function peak_twisted(par, E)
    n = try par[1] catch; par[1][1] end
    σ = try par[2] catch; par[2][1] end
    µ = try par[3] catch; par[3][1] end

    bkg_l = try par[4] catch; par[4][1] end
    bkg_r = try par[5] catch; par[5][1] end
    b = (bkg_r - bkg_l) / (maximum(E) - minimum(E))
    a = bkg_l - b * minimum(E)

    c = try par[6] catch; par[6][1] end
    d = try par[7] catch; par[7][1] end
    β = try par[8] catch; par[8][1] end

    d = b > 0 ? -d : d
    correction = d < 0 ? -1*minimum(d/2 .* erfc.((E .- µ) ./ (sqrt(2)*σ))) : 0.0
    result = @. n / sqrt(2*pi*σ^2) * exp(-(E - µ)^2 / (2*σ^2)) + # Gaussian
        a + b*E + # Linear bkg
        d/2 * erfc((E - µ) / (sqrt(2)*σ)) + correction +  # Low-E step
        c/(2*β) * exp((E-µ) / β + σ^2 / (2*β^2)) * erfc((E - µ) / (sqrt(2)*σ) + σ / (sqrt(2)*β)) # Low-E tail
end

function broad_peak(E, par)
    n = try par[1] catch; par[1][1] end
    σ = try par[2] catch; par[2][1] end
    µ = try par[3] catch; par[3][1] end

    bkg_l = try par[4] catch; par[4][1] end
    bkg_r = try par[5] catch; par[5][1] end
    b = (bkg_r - bkg_l) / (maximum(E) - minimum(E))
    a = bkg_l - b * minimum(E)

    c = try par[6] catch; par[6][1] end
    d = try par[7] catch; par[7][1] end
    β = try par[8] catch; par[8][1] end

    n2 = try par[9]  catch; par[9][1] end
    σ2 = try par[10] catch; par[10][1] end

    d = b > 0 ? -d : d
    correction = d < 0 ? -1*minimum(d/2 .* erfc.((E .- µ) ./ (sqrt(2)*σ))) : 0.0
    result = @. n / sqrt(2*pi*σ^2) * exp(-(E - µ)^2 / (2*σ^2)) + # Gaussian
        n2 / sqrt(2*pi*σ2^2) * exp(-(E - µ)^2 / (2*σ2^2)) + # Gaussian 2
        a + b*E + # Linear bkg
        d/2 * erfc((E - µ) / (sqrt(2)*σ)) .+ correction +  # Low-E step
        c/(2*β) * exp((E-µ) / β + σ^2 / (2*β^2)) * erfc((E - µ) / (sqrt(2)*σ) + σ / (sqrt(2)*β)) # Low-E tail
end

function tail_model(x, par)
    scale = try par[1] catch; par[1][1] end
    σ     = try par[2] catch; par[2][1] end
    μ     = try par[3] catch; par[3][1] end
    C     = try par[4] catch; par[4][1] end
    l     = try par[5] catch; par[5][1] end
    e     = try par[6] catch; par[6][1] end
    t     = try par[7] catch; par[7][1] end
    d     = 0
    return @. scale * exp(-0.5 * ((x - μ)^2) / (σ^2)) / (sqrt(2 * π * σ^2)) + C * (exp(e*(x-l)) + d) / (exp((x-l)*t) +1)
end

function peak_model(x, par)
    scale = par[1]
    σ     = par[2]
    μ     = par[3]

    #step = x <= µ ? par[4] : par[5]
    if length(x) > 1
        step = []
        for i in x
            push!(step, i <= µ ? par[4] : par[5])
        end
        #mwa_step = isodd(round(length(x)/4)) ? Int.(round(length(x)/4)) : Int.(round(length(x)/4)-1)
        s = 4 * par[2] / 0.1
        mwa_step = isodd(Int.(round(s))) ? Int.(round(s)) : Int.(round(s)) + 1
        step = mwa(Array{Float64,1}(step),mwa_step)
    else
        step = x <= µ ? par[4] : par[5]
    end
    return scale .* exp.(-0.5 * ((x .- μ).^2) / (σ^2)) ./ (sqrt(2 * π * σ^2)) .+ step
end
function peak_model2(par, x)
    scale = par.scale
    σ     = par.σ
    μ     = par.µ
    cp0   = par.cp0
    #m0    = par.m0
    cp1   = par.cp1
    #m1    = par.m1
#     s     = par.s

    if length(x) > 1
        step = []
        for i in x
            push!(step, i <= µ ? cp0 : cp1)
        end
        s = 4 * σ / 0.1
#         mwa_step = isodd(round(length(x)/4)) ? Int.(round(length(x)/4)) : Int.(round(length(x)/4)-1)
        mwa_step = isodd(Int.(round(s))) ? Int.(round(s)) : Int.(round(s)) + 1
        step = mwa(Array{Float64,1}(step),mwa_step)
    else
        step = x <= µ ? cp0 : cp1
    end
    return @. scale * exp(-0.5 * ((x - μ)^2) / (σ^2)) / (sqrt(2 * π * σ^2)) + step
end

function model(x, par)
    scale = par[1]
    σ     = par[2]
    μ     = par[3]
    cp0   = par[4]
    m1    = par[5]
    return @. scale * exp(-0.5 * ((x - μ)^2) / (σ^2)) / (sqrt(2 * π * σ^2)) + m1 * x + cp0
end
function model2(par, x)
    scale = par.scale
    σ     = par.σ
    μ     = par.µ
    cp0   = par.cp0
    m1    = par.m1
    #m2    = par[9]
    return @. scale * exp(-0.5 * ((x - μ)^2) / (σ^2)) / (sqrt(2 * π * σ^2)) + m1 * x + cp0
end
function double_model(x, par)
    scale = par[1]
    σ     = par[2]
    μ     = par[3]
    scale2 = par[4]
    σ2     = par[5]
    μ2     = par[6]
    cp0   = par[7]
    m1    = par[8]
    #m2    = par[9]
    return @. scale * exp(-0.5 * ((x - μ)^2) / (σ^2)) / (sqrt(2 * π * σ^2)) + scale2 * exp(-0.5 * ((x - μ2)^2) / (σ2^2)) / (sqrt(2 * π * σ2^2)) + m1 * x + cp0# + m2 * x
end
function double_model2(par,x)
    scale = par.scale
    σ     = par.σ
    μ     = par.µ
    scale2 = par.scale2
    σ2     = par.σ2
    μ2     = par.µ2
    cp0   = par.cp0
    m1    = par.m1
    #m2    = par[9]
    return @. scale * exp(-0.5 * ((x - μ)^2) / (σ^2)) / (sqrt(2 * π * σ^2)) + scale2 * exp(-0.5 * ((x - μ2)^2) / (σ2^2)) / (sqrt(2 * π * σ2^2)) + m1 * x + cp0# + m2 * x
end
# function tail_model(x, par)
#     scale = par[1]
#     σ     = par[2]
#     μ     = par[3]
#     C     = par[4]
#     l     = par[5]
#     e     = par[6]
#     t     = par[7]
#     d     = 0
#     return @. scale * exp(-0.5 * ((x - μ)^2) / (σ^2)) / (sqrt(2 * π * σ^2)) + C * (exp(e*(x-l)) + d) / (exp((x-l)/t) +1)
# end
# function tail_model2(par, x)
#     scale = par.scale
#     σ     = par.σ
#     μ     = par.µ
#     C     = par.C
#     l     = par.l
#     e     = par.e
#     t     = par.t
#     d     = 0
#     return @. scale * exp(-0.5 * ((x - μ)^2) / (σ^2)) / (sqrt(2 * π * σ^2)) + C * (exp(e*(x-l)) + d) / (exp((x-l)/t) +1)
# end
@. linmodel(x, p) = p[1]*x + p[2]
@. hypmodel_simple(x, p) = sqrt(p[1]+p[2]/x)
@. hypmodel(x, p) = sqrt(p[1]+p[2]/(x^2))
@. poly3model(x, p) = p[1]*x^3 + p[2]*x^2 + p[3]*x + p[4]
@. atanmodel(x, p) = p[1] * atan(p[2]*x + p[3]) + p[4]
@. sqrt_fct(x,p) = sqrt( p[1] + p[2] * x + p[3] * x * x )

function fit_bat(hist, fit_function, prior, fit_ranges; nsamples=1*10^5, nchains=4, second_ff=fit_function)

    likelihood = let h = hist, f = fit_function
        # Histogram counts for each bin as an array:
        observed_counts = h.weights

        # Histogram binning:
        bin_edges = h.edges[1]
        bin_edges_left = bin_edges[1:end-1]
        bin_edges_right = bin_edges[2:end]
        bin_widths = bin_edges_right - bin_edges_left
        bin_centers = (bin_edges_right + bin_edges_left) / 2

        params -> begin
            # Log-likelihood for a single bin:
            function bin_log_likelihood(i)
                # Simple mid-point rule integration of fit function `f` over bin:
                expected_counts = bin_widths[i] * f(params, bin_centers[i])
                logpdf(Poisson(expected_counts), observed_counts[i])
            end

            # Sum log-likelihood over bins:
            idxs = eachindex(observed_counts)
            ll_value = bin_log_likelihood(idxs[1])
            for i in idxs[2:end]
                ll_value += bin_log_likelihood(i)
            end

            # Wrap `ll_value` in `LogDVal` so BAT knows it's a log density-value.
            return LogDVal(ll_value)
        end
    end
    parshapes = varshape(prior)
    posterior = PosteriorDensity(likelihood, prior);

    pretunesamples = 10000
    max_ncycles = 30
    T = Float64
    tuning = BAT.AdaptiveMetropolisTuning(
        r = 1.0,
        λ = 0.5,
        α = 0.05..0.9,
        β = 1.5,
        c = 1e-4..1e2
    )
    convergence = BAT.BrooksGelmanConvergence(
        threshold = T(5),
        corrected = false
    )
    init = BAT.MCMCInitStrategy(
        init_tries_per_chain = 8..128,
        max_nsamples_init = pretunesamples,
        max_nsteps_init = pretunesamples * 10,
        max_time_init = Inf
    )
    burnin = BAT.MCMCBurninStrategy(
        max_nsamples_per_cycle = pretunesamples,
        max_nsteps_per_cycle = pretunesamples * 10,
        max_time_per_cycle = Inf,
        max_ncycles = max_ncycles
    )
    parshapes = varshape(prior)
    posterior = PosteriorDensity(likelihood, prior);
    nsamples = 1*10^4
    nchains = 4

    samples = bat_sample(Philox4x((123, 456)), posterior, (nsamples, nchains), MetropolisHastings(),
        max_nsteps = 10 * nsamples,
        max_time = Inf,
        tuning = tuning,
        init = init,
        burnin = burnin,
        convergence = convergence,
    ).result;
    params = mode(samples)[1]
    error = std(samples)[1]
    # println("Truth: $p0")
    println("Mode: $(mode(samples))")
    println("Stddev: $(std(samples))")

    fitf = RadiationSpectra.FitFunction{Float64}( second_ff, 1, length(params));
    # set_parameter_bounds!(fitf, bounds)
    RadiationSpectra.set_fitranges!(fitf, (fit_ranges, ) )

    par = []
    err = []
    for x in eachindex(params)
        push!(par, params[x])
        push!(err, error[x])
    end
#     params = [params.scale[1], params.σ[1], params.µ[1], params.scale2[1], params.σ2[1], params.µ2[1], params.cp0[1], params.m1[1]]
#     error = [error.scale[1], error.σ[1], error.µ[1], error.scale2[1], error.σ2[1], error.µ2[1], error.cp0[1], error.m1[1]]
    fitf.fitted_parameters = par

    return par, err, fitf
end


text:
"""
Inverse operation of scatter
"""
function gather(input::AbstractArray{T,N}, index::AbstractArray{<:Integer,N}, dims::Integer) where {T,N}
    @assert dims <= N "Specified dimensions must lower or equal to the rank of input matrix."
    out = similar(index, T)
    @inbounds for x = CartesianIndices(out)
        tup = collect(Tuple(x))
        tup[dims] = index[x]
        view(out, x) .= view(input, tup...)
    end
    return out
end

function gather(input::Matrix{T}, index::Array{Int}) where T
    out = Array{T}(undef, size(input,1), size(index)...)
    @inbounds for ind = CartesianIndices(index)
        view(out, :, ind) .= view(input, :, index[ind])
    end
    return out
end

gather(input::Fill{T,2,<:Any}, index::Array{Int}) where T = gather(Matrix(input), index)

function gather_indices(X::Array{T}) where T
    Y = DefaultDict{T,Vector{CartesianIndex}}(CartesianIndex[])
    @inbounds for (ind, val) = pairs(X)
        push!(Y[val], ind)
    end
    Y
end


text:
module SCIPCInterface

include("scip_defines.jl")
include("scip_enums.jl")
include("scip_typedefs.jl")
include("scip_structs.jl")
include("scip_functions.jl")
include("scip_constructors.jl")

end


text:
"""
    create_hypercube{T<:Real}(origin::Vector{T}, edgelength::T)::HyperRectVolume

*AHMI-internal, not part of stable public API.*

creates a hypercube shaped spatial volume
"""
function create_hypercube(
    origin::Vector{T},
    edgelength::T)::HyperRectVolume{T} where {T<:AbstractFloat}

    dim = length(origin)
    lo = zeros(T, dim)
    hi = zeros(T, dim)

    _setcubeboundaries!(lo, hi, origin, edgelength)

    HyperRectVolume(lo, hi)
end

"""
    create_hypercube!{T<:Real}(origin::Vector{T}, edgelength::T)::HyperRectVolume

*AHMI-internal, not part of stable public API.*

resizes a hypercube shaped spatial volume
"""
function modify_hypercube!(
    rect::HyperRectVolume{T},
    neworigin::Vector{T},
    newedgelength::T) where {T<:AbstractFloat}

    _setcubeboundaries!(rect.lo, rect.hi, neworigin, newedgelength)

    nothing
end

@inline function _setcubeboundaries!(
    lo::Vector{T},
    hi::Vector{T},
    origin::Vector{T},
    edgelength::T) where {T<:AbstractFloat}

    for i = 1:length(lo)
        lo[i] = origin[i] - edgelength * 0.5
        hi[i] = origin[i] + edgelength * 0.5
    end

    nothing
end

"""
    find_hypercube_centers(dataset::DataSet{T, I}, whiteningresult::WhiteningResult, settings::HMISettings)::Vector{I}

*AHMI-internal, not part of stable public API.*

finds possible starting points for the hyperrectangle creation
"""
function find_hypercube_centers(
    dataset::DataSet{T, I},
    whiteningresult::WhiteningResult{T},
    settings::HMISettings)::Bool where {T<:AbstractFloat, I<:Integer}

    sortLogProb = sortperm(dataset.logprob, rev = true)

    NMax = min(settings.max_startingIDs, round(I, sqrt(dataset.N * settings.max_startingIDs_fraction)))
    NConsidered = round(I, sqrt(dataset.N) * settings.max_startingIDs_fraction)
    @debug "Considered seed samples $NConsidered"

    discardedsamples = falses(dataset.N)

    testlength = find_density_test_cube_edgelength(dataset.data[:, sortLogProb[1]], dataset, round(I, sqrt(dataset.N)))
    @debug "Edge length of global mode cube: $testlength"

    maxprob = dataset.logprob[sortLogProb[1]]
    startingsamples = zeros(I, NMax)
    cntr = 0
    @showprogress for n::I in sortLogProb[1:NConsidered]
        if discardedsamples[n]
            continue
        end
        if cntr == NMax || dataset.logprob[n] < maxprob - log(whiteningresult.targetprobfactor)
            break
        end

        cntr += 1
        startingsamples[cntr] = n

        cubevol = create_hypercube(dataset.data[:, n], testlength)
        cube = IntegrationVolume(dataset, cubevol, true)

        discardedsamples[cube.pointcloud.pointIDs] .= true
    end
    resize!(startingsamples, cntr)

    success = true
    if cntr < settings.warning_minstartingids
        stop = floor(I, length(sortLogProb) * 0.2)
        step = floor(I, stop / settings.warning_minstartingids)
        if step == 0 step = 1 end
        startingsamples = sortLogProb[1:step:stop]
        success = false
        @warn "Returned minimum number of starting points: $(settings.warning_minstartingids)"
    end

    @debug "Selected Starting Samples: $cntr out of $(dataset.N) points"
    dataset.startingIDs = startingsamples

    success
end

function find_density_test_cube_edgelength(
    mode::Vector{T},
    dataset::DataSet{T, I},
    points::I = 100)::T where {T<:AbstractFloat, I<:Integer}

    find_density_test_cube(mode, dataset, points)[1]
end

function find_density_test_cube(
    mode::Vector{T},
    dataset::DataSet{T, I},
    points::I)::Tuple{T, IntegrationVolume{T, I}} where {T<:AbstractFloat, I<:Integer}

    P = dataset.P

    l::T = 1.0
    tol::T = 1.0
    mult::T = 2.0^(1.0 / P)
    last_change = 0

    rect = create_hypercube(mode, l)
    intvol = IntegrationVolume(dataset, rect, false)
    pt = intvol.pointcloud.points

    iterations = 0
    while pt < points / tol || pt > points * tol
        iterations += 1
        tol += 0.001 * iterations
        if pt > points
            l /= mult
            mult = last_change == -1 ? mult^2.0 : mult^0.5
            last_change = -1
        else
            l *= mult
            mult = last_change == 1 ? mult^2.0 : mult^0.5
            last_change = 1
        end

        modify_hypercube!(rect, mode, l)
        modify_integrationvolume!(intvol, dataset, rect, false)
        pt = intvol.pointcloud.points
    end

    ## @debug "Tolerance Test Cube: Iterations $iterations\tPoints: $(intvol.pointcloud.points)\ttarget Points: $points"

    l, intvol
end

function create_initialhypercube(
    mode::Array{T, 1},
    dataset::DataSet{T, I},
    targetprobfactor::T) where {T<:AbstractFloat, I<:Integer}

    edgelength::T = 1.0

    cube = create_hypercube(mode, edgelength)
    vol = IntegrationVolume(dataset, cube, true)

    while vol.pointcloud.points > 0.01 * dataset.N
        edgelength *= 0.5^(1/dataset.P)
        modify_hypercube!(cube, mode, edgelength)
        modify_integrationvolume!(vol, dataset, cube, true)
    end

    tol = 1.0
    step = 0.7
    direction = 0
    PtsIncrease = 0.0

    it = 0
    while vol.pointcloud.probfactor < targetprobfactor / tol || vol.pointcloud.probfactor > targetprobfactor
        tol += 0.01 * it
        it += 1

        if vol.pointcloud.probfactor > targetprobfactor
            #decrease side length
            edgelength *= step

            step = adjuststepsize!(step,direction == -1)
            direction = -1
        else
            #increase side length
            edgelength /= step

            step = adjuststepsize!(step, direction == 1)
            direction = 1
        end
        PtsIncrease = vol.pointcloud.points
        modify_hypercube!(cube, mode, edgelength)
        modify_integrationvolume!(vol, dataset, cube, true)

        PtsIncrease = vol.pointcloud.points / PtsIncrease

        if vol.pointcloud.points > 0.01 * dataset.N && vol.pointcloud.probfactor < targetprobfactor
            break
        end
    end

    return cube, vol
end

@enum Edge LowerEdge UpperEdge
@enum Adaption IncreaseVolume DecreaseVolume NoChange Init
function modify_edge!(
    dataset::DataSet{T, I},
    change_mod::T,
    adaption::Adaption,
    edge::Edge,
    spatialvolume::HyperRectVolume{T},
    searchvol::HyperRectVolume{T},
    vol::IntegrationVolume{T, I, HyperRectVolume{T}},
    newvol::IntegrationVolume{T, I, HyperRectVolume{T}},
    dim::I) where {T<:AbstractFloat, I<:Integer}

    margin = spatialvolume.hi[dim] - spatialvolume.lo[dim]
    @assert margin != 0

    buffer = (edge == UpperEdge) ? spatialvolume.hi[dim] : spatialvolume.lo[dim]
    if edge == UpperEdge
        if adaption == IncreaseVolume
            spatialvolume.hi[dim] += margin * change_mod
        elseif adaption == DecreaseVolume
            spatialvolume.hi[dim] -= margin * change_mod
        else
            @error "No edge modification possible: Volume change not specified"
        end
    else
        if adaption == IncreaseVolume
            spatialvolume.lo[dim] -= margin * change_mod
        elseif adaption == DecreaseVolume
            spatialvolume.lo[dim] += margin * change_mod
        else
            @error "No edge modification possible: Volume change not specified"
        end
    end
    prevpts = vol.pointcloud.points
    resize_integrationvol!(newvol, vol, dataset, dim, spatialvolume, false, searchvol)
    pts_change = newvol.pointcloud.points / prevpts
    return pts_change, buffer, margin
end
function accept_edge_modification!(
    vol::IntegrationVolume{T, I, HyperRectVolume{T}},
    newvol::IntegrationVolume{T, I, HyperRectVolume{T}}) where {T<:AbstractFloat, I<:Integer}

    copy!(vol, newvol)
end
function reject_edge_modification!(
    vol::IntegrationVolume{T, I, HyperRectVolume{T}},
    newvol::IntegrationVolume{T, I, HyperRectVolume{T}},
    spatialvolume::HyperRectVolume{T},
    buffer::T,
    edge::Edge,
    dim::I) where {T<:AbstractFloat, I<:Integer}

    newvol.pointcloud.points = vol.pointcloud.points
    if edge == UpperEdge
        spatialvolume.hi[dim] = buffer
    elseif edge == LowerEdge
        spatialvolume.lo[dim] = buffer
    else
        @error "No Edge modification to reject"
    end
end

function adapt_dimension!(
    dataset::DataSet{T, I},
    vol::IntegrationVolume{T, I, HyperRectVolume{T}},
    newvol::IntegrationVolume{T, I, HyperRectVolume{T}},
    spatialvolume::HyperRectVolume{T},
    searchvol::HyperRectVolume{T},
    edge::Edge,
    dim::I,
    increase::T,
    decrease::T,
    targetprobfactor::T,
    max_iterations_per_dimension::I) where {T<:AbstractFloat, I<:Integer}

    ptsTolInc::T = dataset.tolerance
    ptsTolDec::T = dataset.tolerance * 1.1

    last_change = Init
    current_iteration = 0
    local change = false

    while last_change != NoChange &&
        vol.pointcloud.probfactor > 1.0 &&
        current_iteration < max_iterations_per_dimension

        current_iteration += 1
        pts_change, buffer, margin = modify_edge!(dataset, increase, IncreaseVolume, edge, spatialvolume, searchvol, vol, newvol, dim)

        if newvol.pointcloud.probfactor < targetprobfactor && pts_change > (1.0 + increase / ptsTolInc) && last_change != DecreaseVolume
            accept_edge_modification!(vol, newvol)
            change = true
            last_change = IncreaseVolume
        else
            reject_edge_modification!(vol, newvol, spatialvolume, buffer, edge, dim)
            pts_change, buffer, margin = modify_edge!(dataset, decrease, DecreaseVolume, edge, spatialvolume, searchvol, vol, newvol, dim)

            if pts_change > (1.0 - decrease / ptsTolDec) && last_change != IncreaseVolume
                accept_edge_modification!(vol, newvol)
                change = true
                last_change = DecreaseVolume
            else
                reject_edge_modification!(vol, newvol, spatialvolume, buffer, edge, dim)
                last_change = NoChange
            end
        end
    end

    change, current_iteration
end

function adapt_hyperrectangle!(
    dataset::DataSet{T, I},
    vol::IntegrationVolume{T, I, HyperRectVolume{T}},
    newvol::IntegrationVolume{T, I, HyperRectVolume{T}},
    spatialvolume::HyperRectVolume{T},     #spatial volume
    searchvol::HyperRectVolume{T}, #adaptively proposed changes to the spatial volume
    increase::T,
    decrease::T,
    dim::I,
    targetprobfactor::T,
    max_iterations_per_dimension::I) where {T<:AbstractFloat, I<:Integer}

    change1, iterations1 = adapt_dimension!(dataset, vol, newvol, spatialvolume, searchvol,
        LowerEdge, dim, increase, decrease, targetprobfactor, max_iterations_per_dimension)

    change2, iterations2 = adapt_dimension!(dataset, vol, newvol, spatialvolume, searchvol,
        UpperEdge, dim, increase, decrease, targetprobfactor, max_iterations_per_dimension)

    change1 && change2, iterations1 + iterations2
end

"""
    create_hyperrectangle(...)

*AHMI-internal, not part of stable public API.*

This function creates a hyper-rectangle around each starting sample.
It starts by building a hyper-cube  and subsequently adapts each face individually,
thus turning the hyper-cube into a hyper-rectangle.
The faces are adjusted in a way to match the shape of the distribution as best as possible.
"""
function create_hyperrectangle(
    id::I,
    dataset::DataSet{T, I},
    targetprobfactor::T,
    settings::HMISettings) where {T<:AbstractFloat, I<:Integer}

    mode = dataset.data[:, id]
    initialcube, vol = create_initialhypercube(mode, dataset, targetprobfactor)

    local change = true
    volbuffer::IntegrationVolume{T, I} = deepcopy(vol)

    spatialvolume::HyperRectVolume{T} = deepcopy(vol.spatialvolume)
    searchvol::HyperRectVolume{T} = deepcopy(spatialvolume)

    increase_default = settings.rect_increase
    increase = increase_default
    decrease = 1.0 - 1.0 / (1.0 + increase)

    min_points = 5
    max_iterations_per_dimension = 20
    total_iterations = 0

    while change && vol.pointcloud.probfactor > 1.0

        if vol.pointcloud.points * increase < min_points
            increase *= ceil(I, min_points / (vol.pointcloud.points * increase))
            decrease = 1.0 - 1.0 / (1.0 + increase)
        elseif increase > increase_default && vol.pointcloud.points * increase_default > 2 * min_points
            increase = increase_default
            decrease = 1.0 - 1.0 / (1.0 + increase)
        end

        for p::I = 1:dataset.P
            change, iterations = adapt_hyperrectangle!(dataset, vol, volbuffer, spatialvolume, searchvol,
                increase, decrease, p, targetprobfactor, max_iterations_per_dimension)
            total_iterations += iterations
        end
    end

    res = search(dataset, vol.spatialvolume, true)
    resize!(vol.pointcloud.pointIDs, res.points)
    copyto!(vol.pointcloud.pointIDs, res.pointIDs)
    vol.pointcloud.points = res.points
    vol.pointcloud.maxLogProb = res.maxLogProb
    vol.pointcloud.minLogProb = res.minLogProb
    vol.pointcloud.probfactor = exp(vol.pointcloud.maxLogProb - vol.pointcloud.minLogProb)
    vol.pointcloud.probweightfactor = exp(vol.pointcloud.maxWeightProb - vol.pointcloud.minWeightProb)

    vol, initialcube, total_iterations
end

@inline function adjuststepsize!(Step, Increase::Bool)
    if Increase
        return Step * 0.5
    else
        return Step * 2.0
    end
end


text:
function compute_cholesky_diff(collection::Mongoc.Collection)
    choleskytimes = getcholeskytime_all(collection)
    for d in collection
        if haskey(d["features"], "solver") && haskey(d["features"]["solver"], "solving_time")
            instance_name = d["_id"]["instance_name"]
            oid = d["_id"]["oid"]
            time = d["features"]["solver"]["solving_time"]
            typeof(time) <: AbstractString && (time = parse(Float64, time))
            time_diff = time - choleskytimes[instance_name]
            setfeatures_decomposition!(collection,
                                       instance_name,
                                       oid,
                                       Mongoc.BSON("solver.cholesky_diff" => time_diff))
        end
    end
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
using Yao.EasyBuild, Yao
using Test
using LinearAlgebra
using YaoBlocks.ConstGate

single_swap_test_circuit(ϕ::Real) = hadamard_test_circuit(SWAP, ϕ)
single_swap_test(reg, ϕ::Real) = hadamard_test(SWAP, reg, ϕ)

@testset "state overlap" begin
    reg1 = rand_state(3) |> focus!(1,2)
    rho1 = reg1 |> ρ
    reg2 = rand_state(3) |> focus!(1,2)
    rho2 = reg2 |> ρ
    reg3 = rand_state(3) |> focus!(1,2)
    rho3 = reg3 |> ρ
    desired = tr(Matrix(rho1)*Matrix(rho2))
    c = swap_test_circuit(2, 2, 0)
    res = expect(put(5, 1=>Z), join(join(reg2, reg1), zero_state(1)) |> c) |> tr
    @test desired ≈ res
    desired = tr(Matrix(rho1)*Matrix(rho2)*Matrix(rho3)) |> real
    c = swap_test_circuit(2, 3, 0)
    res = expect(put(7, 1=>Z), reduce(join, [reg3, reg2, reg1, zero_state(1)]) |> c) |> tr |> real
    @test desired ≈ res
    desired = tr(Matrix(rho1)*Matrix(rho2)*Matrix(rho3)) |> imag
    c = swap_test_circuit(2, 3, -π/2)
    res = expect(put(7, 1=>Z), reduce(join, [reg3, reg2, reg1, zero_state(1)]) |> c) |> tr |> real
    @test desired ≈ res
end

@testset "hadamard test" begin
    nbit = 4
    U = put(nbit, 2=>Rx(0.2))
    reg = rand_state(nbit)

    @test hadamard_test(U, reg, 0.0) ≈ real(expect(U, reg))
    @test hadamard_test(U, reg, -π/2) ≈ imag(expect(U, reg))

    reg = zero_state(2) |> EasyBuild.singlet_block()
    @test single_swap_test(reg, 0) ≈ -1

    reg = zero_state(2)
    @test single_swap_test(reg, 0) ≈ 1
    reg = product_state(2, 0b11)
    @test single_swap_test(reg, 0) ≈ 1
end


text:
# # default copy is to loop through
# # override this for most operators.
function default_BlockBandedMatrix(S::Operator)
    ret = BlockBandedMatrix(Zeros, S)

    @inbounds for J=blockaxes(ret,2), K=blockcolrange(ret,Int(J))
        ret[K,J] = view(S,K,J)
    end
    ret
end
#

# diagblockshift gives the shift for the diagonal block of an operator
# that is, trace an operator down the diagonal.  What blocks correspond to the
# block diagonal?
# this is used to determine how many blocks to pad in the QR decomposition, as
# every lower block gets added to the upper daigonal

diagblockshift(a,b) = error("Developer: Not implemented for blocklengths $a, $b")

function diagblockshift(a::AbstractRange, b::AbstractRange)
    @assert step(a) == step(b)
    first(b)-first(a)
end
diagblockshift(a::AbstractUnitRange, b::AbstractUnitRange) = first(b)-first(a)

#TODO: generalize
function diagblockshift(a::AbstractFill{Int},b::AbstractFill{Int})
    @assert getindex_value(a) == getindex_value(b)
    0
end

diagblockshift(a::AbstractFill{Int},b::Vcat{Int,1,<:Tuple{V1,<:AbstractFill{Int}}}) where {V1 <: AbstractVector{Int}} =
    max(0,-diagblockshift(b,a))

function diagblockshift(a::Vcat{Int,1,<:Tuple{V1,<:AbstractFill{Int}}},b::AbstractFill{Int}) where V1 <: AbstractVector{Int}
    @assert getindex_value(a.args[end]) == getindex_value(b)
    isempty(a.args[1]) && return diagblockshift(a.args[2],b)
    a1, b1 = a[1],b[1]
    a1 == b1 && return diagblockshift(Vcat(a.args[1][2:end],a.args[2]),b)
    a1 >  b1 && length(a.args[1]) == 1 && return 0
    a1 >  b1 && return max(0,-1+diagblockshift(flatten(([a1-b1;a.args[1][2:end]],a.args[2]),b)))
    a1 <  b1 && length(a.args[1]) == 1 && return 1
    # a1 <  b1 &&
    return 1+diagblockshift(Vcat(a.args[1][2:end],a.args[2]),Vcat([b1-a1],b))
end

function diagblockshift(a::Vcat{Int,1,<:Tuple{V1,<:AbstractFill{Int}}},
                        b::Vcat{Int,1,<:Tuple{V2,<:AbstractFill{Int}}}) where {V1 <: AbstractVector{Int},V2 <: AbstractVector{Int}}
    isempty(a.args[1]) && return diagblockshift(a.args[2],b)
    isempty(b.args[1]) && return diagblockshift(a,b.args[2])
    a1, b1 = a[1],b[1]
    a1 == b1 && return diagblockshift(Vcat(a.args[1][2:end],a.args[2]),Vcat(b.args[1][2:end],b.args[2]))
    a1 >  b1 && return max(0,-1+diagblockshift(Vcat([a1-b1;a.args[1][2:end]],a.args[2]),
                                               Vcat(b.args[1][2:end],b.args[2])))
    # a1 <  b1 &&
    return 1+diagblockshift(Vcat(a.args[1][2:end],a.args[2]),Vcat([b1-a1;b.args[1][2:end]],b.args[2]))
end

diagblockshift(op::Operator) =
    diagblockshift(blocklengths(domainspace(op)),blocklengths(rangespace(op)))

function CachedOperator(::Type{BlockBandedMatrix},op::Operator;padding::Bool=false)
    l,u=blockbandwidths(op)
    padding && (u+=l+diagblockshift(op))
    data=BlockBandedMatrix{eltype(op)}(undef,
                                        Vector{Int}(), Vector{Int}(),
                                        (l,u))
    CachedOperator(op,data,(0,0),domainspace(op),rangespace(op),(-l,u),padding)
end

## Grow cached operator
#
function resizedata!(B::CachedOperator{T,BlockBandedMatrix{T}},::Colon,col::Integer) where {T<:Number}
    if col > size(B,2)
        throw(ArgumentError("Cannot resize beyound size of operator"))
    end

    if col > B.datasize[2]
        if B.datasize[2] == 0
            datablocksize = Block(0)
        else
            datablocksize = block(domainspace(B),B.datasize[2])
            bs = blockstop(domainspace(B),datablocksize)
            if bs ≠ B.datasize[2]
                error("Internal Error: $(B.datasize) is not lined up with the block $datablocksize as the last column doesn't end at $bs")
            end
        end

        l,u = blockbandwidths(B.data)
        J = block(domainspace(B),col)
        col = blockstop(domainspace(B),J)  # pad to end of block

        rows = blocklengths(rangespace(B.op))[1:Int(J)+l]
        cols = blocklengths(domainspace(B.op))[1:Int(J)]

        b_size = BlockBandedSizes(Vector{Int}(rows), Vector{Int}(cols), l, u)

        resize!(B.data.data, bb_numentries(b_size))
        B.data = _BlockBandedMatrix(B.data.data, b_size)

        JR = datablocksize+1:J
        KR=blockcolstart(B.data,first(JR)):blockcolstop(B.data,last(JR))
        copyto!(view(B.data,KR,JR), view(B.op,KR,JR))

        B.datasize = (blockstop(rangespace(B),last(KR)),col)
    end

    B
end

function resizedata!(B::CachedOperator{T,BlockBandedMatrix{T}},n::Integer,m::Integer) where {T<:Number}
    N = block(rangespace(B),n)
    m̃ = blockstart(domainspace(B),N)
    resizedata!(B,:,max(m,m̃))
end

## QR
# we use a RaggedMatrix to represent the growing lengths of the
# householder reflections
QROperator(R::CachedOperator{T,BlockBandedMatrix{T}}) where {T} =
    QROperator(R,RaggedMatrix{T}(undef, 0, Int[]),0)

# function resizedata!(QR::QROperator{CachedOperator{T,BlockBandedMatrix{T},
#                                           MM,DS,RS,BI}},
#                      ::Colon, col) where {T,MM,DS,RS,BI}
#     if col ≤ QR.ncols
#         return QR
#     end
#
#     MO=QR.R_cache
#     W=QR.H
#
#     R=MO.data
#     cs=colstop(MO,col)
#
#     if cs ≥ MO.datasize[1]
#         resizedata!(MO,cs+100,:)  # add 100 rows
#         R=MO.data
#     end
#
#     if col > size(W,2)
#         m=size(W,2)
#         resize!(W.cols,2col+1)
#
#         for j=m+1:2col
#             cs=colstop(MO,j)
#             W.cols[j+1]=W.cols[j] + cs-j+1
#             W.m=max(W.m,cs-j+1)
#         end
#
#         resize!(W.data,W.cols[end]-1)
#     end
#
#     for k=QR.ncols+1:col
#         cs = colstop(R,k)
#         W[1:cs-k+1,k] = view(R,k:cs,k) # diagonal and below
#         wp=view(W,1:cs-k+1,k)
#         W[1,k]+= flipsign(norm(wp),W[1,k])
#         normalize!(wp)
#
#         # scale banded entries
#         for j=k:rowstop(R,k)
#             v=view(R,k:cs,j)
#             dt=dot(wp,v)
#             LinearAlgebra.axpy!(-2*dt,wp,v)
#         end
#
#         # scale banded/filled entries
#         for j=rowstop(R,k)+1:rowstop(R,cs)
#             csrt=colstart(R,j)
#             v=view(R,csrt:cs,j)  # shift down each time
#             wp2=view(wp,csrt-k+1:cs-k+1)
#             dt=dot(wp2,v)
#             LinearAlgebra.axpy!(-2*dt,wp2,v)
#         end
#     end
#     QR.ncols=col
#     QR
# end

# always resize by column
resizedata!(QR::QROperator{CachedOperator{T,BlockBandedMatrix{T},
                               MM,DS,RS,BI}},
                     ::Colon, col::Int) where {T,MM,DS,RS,BI} =
    resizedata!(QR, :, block(domainspace(QR.R_cache),col))

function resizedata!(QR::QROperator{CachedOperator{T,BlockBandedMatrix{T},
                               MM,DS,RS,BI}},
                     ::Colon, COL::Block) where {T<:BlasFloat,MM,DS,RS,BI}
     MO = QR.R_cache
     W = QR.H
     R = MO.data

     l,u = blockbandwidths(R)

     ds = domainspace(QR)
     col = blockstop(ds, COL)  # last column

     if col ≤ QR.ncols
         return QR
     end

     J_start = Int(block(ds, QR.ncols+1)) # first block-column not factorized
     @assert blockstart(ds, J_start) == QR.ncols+1  # we want to make sure we haven't partially factorized a column

     # last block, convoluted def to match blockbandedmatrix
     J_col = Int(COL)
     K_end = Int(blockcolstop(MO, J_col))  # last row block in last column
     J_end = Int(blockrowstop(MO, K_end))  # QR will affect up to this column
     j_end = blockstop(domainspace(MO), J_end)  # we need to resize up this column
     sz = sizeof(T)

     if j_end ≥ MO.datasize[2]
         # add columns up to column rs, which is last column affected by QR
         resizedata!(MO, :, j_end)
         R = MO.data
     end

     if col > size(W,2)
         # resize Householder matrix
         m = size(W,2)
         resize!(W.cols, 2col+1)

         for j=m+1:2col
             cs = blockstop(rangespace(MO), blockcolstop(MO, block(domainspace(MO),j)))
             W.cols[j+1]=W.cols[j] + cs-j+1
             W.m=max(W.m,cs-j+1)
         end

         resize!(W.data, W.cols[end]-1)
     end

     w = pointer(W.data)
     r = pointer(R.data)

     bs = R.block_sizes

     for j =QR.ncols+1 : col   # first column of block
         bi = findblockindex.(bs.axes, (j, j)) # converts from global indices to block indices
         K1, J1 = Int.(block.(bi))  # this is the diagonal block corresponding to j
         κ, ξ = blockindex.(bi)

         st = bs.block_strides[J1]  # the stride of the matrix
         shft = bs.block_starts[K1,J1]-1 + st*(ξ-1) + κ-1 # the index of the pointer to the j, j entry

         K_CS = Int(blockcolstop(R,J1)) # last row in J-th blockcolumn
         k_end = last(bs.axes[1][Block(K_CS)])

         w_j = W.cols[j]  # the data index  for the j-th column of W
         wp = w+sz*(w_j-1)          # j-th column of W

         M = k_end - j + 1 # the number of entries we are diagonalizing. we know the stride tells us the total number of rows

         BLAS.blascopy!(M, r+sz*shft, 1, wp, 1) # copy the column into W

         # we need to scale the first entry and then normalize
         W.data[w_j] += flipsign(BLAS.nrm2(M,wp,1), W.data[w_j])
         normalize!(M, wp)

         # scale rest of columns in first block
         # for ξ_2 = 2:

         for ξ_2 = ξ:length(bs.axes[2][Block(J1)])
             # we now apply I-2v*v' in place
             r_sh = r+sz*(shft + st*(ξ_2-ξ)) # the pointer the (j,ξ_2)-th entry
             dt = BandedMatrices.dot(M, wp, 1, r_sh, 1)
             BLAS.axpy!(M, -2*dt, wp, 1, r_sh ,1)
         end

         for J = J1+1:min(K1+u,J_end)
             st = bs.block_strides[J]
             shft = bs.block_starts[K1,J] + κ-2 # the index of the pointer to the j, j entry
             for ξ_2 = axes(bs.axes[2][Block(J)],1)
                 # we now apply I-2v*v' in place
                 # r_sh = r+sz*(shft + st*(ξ_2-1)) # the pointer the (j,ξ_2)-th entry

                 # TODO: remove these debugging statement
                 # @assert w_j-1 + M ≤ length(W.data)
                 # @assert shft + st*(ξ_2-1) + M ≤ length(R.data)
                 # @assert 0 ≤ w_j-1
                 # if ! (0 ≤ shft + st*(ξ_2-1))
                 #     @show shft, st, ξ_2, l, u
                 #     @show κ, bs.block_starts[K1,J]
                 #     @show K1, J
                 #     @show MO.op
                 # end
                 # dt = dot(M, wp, 1, r_sh, 1)
                 # BLAS.axpy!(M, -2*dt, wp, 1, r_sh ,1)

                 dt = dot(view(W.data, w_j:w_j+M-1) ,
                                    view(R.data, shft + st*(ξ_2-1) +1:shft + st*(ξ_2-1) +M))
                 BLAS.axpy!(-2*dt, view(W.data, w_j:w_j+M-1) ,
                                    view(R.data, shft + st*(ξ_2-1) +1:shft + st*(ξ_2-1) +M))
             end
         end
     end

    QR.ncols=col
    QR
end


text:
using Nemo, GibbsTypePriors, JLD

grid_k = [collect(1:25:1000); 1000]

Vnk_numerical_accuracy = GibbsTypePriors.Vnk_NGG.(1000, grid_k, 1.2, 0.2) |> x -> accuracy_bits.(x)
save("test/graphical_tests/saves_for_graphical_tests/accuracy_Vnk_sigma_02_1000.jld", "Vnk_numerical_accuracy", Vnk_numerical_accuracy)

max_stars_repo_name: konkam/GibbsTypePriors

max_issues_repo_name: konkam/GibbsTypePriors

max_forks_repo_name: konkam/GibbsTypePriors


text:
using Test, Logging, Random
using DecisionProgramming

@info "Testing ChandeNode"
@test isa(ChanceNode("A", [], ["a", "b"]), ChanceNode)
@test isa(ChanceNode("B", [], ["x", "y", "z"]), ChanceNode)
@test_throws MethodError ChanceNode(1, [], ["x", "y", "z"])
@test_throws MethodError ChanceNode("B", [1], ["y", "z"])
@test_throws MethodError ChanceNode("B", ["A"], [1, "y", "z"])
@test_throws MethodError ChanceNode("B", ["A"])

@info "Testing DecisionNode"
@test isa(DecisionNode("D", [], ["x", "y"]), DecisionNode)
@test isa(DecisionNode("E", ["C"], ["x", "y"]), DecisionNode)
@test_throws MethodError DecisionNode(1, [], ["x", "y", "z"])
@test_throws MethodError DecisionNode("D", [1], ["y", "z"])
@test_throws MethodError DecisionNode("D", ["A"], [1, "y", "z"])
@test_throws MethodError DecisionNode("D", ["A"])

@info "Testing ValueNode"
@test isa(ValueNode("V", []), ValueNode)
@test isa(ValueNode("V", ["E", "D"]), ValueNode)
@test_throws MethodError ValueNode(1, [])
@test_throws MethodError ValueNode("V", [2])

@info "Testing State"
@test isa(States(State[1, 2, 3]), States)
@test_throws DomainError States(State[0, 1])

@info "Testing paths"
@test vec(collect(paths(States(State[2, 3])))) == [(1, 1), (2, 1), (1, 2), (2, 2), (1, 3), (2, 3)]
@test vec(collect(paths(States(State[2, 3]), Dict(Node(1)=>State(2))))) == [(2, 1), (2, 2), (2, 3)]
@test vec(collect(paths(States(State[2, 3]), FixedPath(Dict(Node(1)=>State(2)))))) == [(2, 1), (2, 2), (2, 3)]

@info "Testing Probabilities"
@test isa(Probabilities(Node(1), [0.4 0.6; 0.3 0.7]), Probabilities)
@test isa(Probabilities(Node(1), [0.0, 0.4, 0.6]), Probabilities)
@test_throws DomainError Probabilities(Node(1), [1.1, 0.1])

@info "Testing DefaultPathProbability"
P = DefaultPathProbability(
    [Node(1), Node(2)],
    [Node[], [Node(1)]],
    [Probabilities(Node(1), [0.4, 0.6]), Probabilities(Node(2), [0.3 0.7; 0.9 0.1])]
)
@test isa(P, DefaultPathProbability)
@test P((State(1), State(2))) == 0.4 * 0.7

@info "Testing Utilities"
@test isa(Utilities(Node(1), Utility[-1.1, 0.0, 2.7]), Utilities)
@test isa(Utilities(Node(1), Utility[-1.1 0.0; 2.7 7.0]), Utilities)

@info "Testing DefaultPathUtility"
U = DefaultPathUtility(
    [Node[2], Node[1, 2]],
    [Utilities(Node(3), Utility[1.0, 1.4]), Utilities(Node(4), Utility[1.0 1.5; 0.6 3.4])]
)
@test isa(U, DefaultPathUtility)
@test U((State(2), State(1))) == Utility(1.0 + 0.6)

@info "Testing InfluenceDiagram"
diagram = InfluenceDiagram()
@test isa(diagram, InfluenceDiagram)
push!(diagram.Nodes, ChanceNode("A", [], ["a", "b"]))
@test isa(diagram, InfluenceDiagram)

@info "Testing add_node! and validate_node"
diagram = InfluenceDiagram()
add_node!(diagram, ChanceNode("A", [], ["a", "b"]))
@test isa(diagram.Nodes[1], ChanceNode)
@test_throws DomainError add_node!(diagram, ChanceNode("A", [], ["a", "b"]))
@test_throws DomainError add_node!(diagram, ChanceNode("C", ["B", "B"], ["a", "b"]))
@test_throws DomainError add_node!(diagram, ChanceNode("C", ["C", "B"], ["a", "b"]))
@test_throws DomainError add_node!(diagram, ChanceNode("C", [], ["a"]))
@test_throws DomainError add_node!(diagram, DecisionNode("A", [], ["a", "b"]))
@test_throws DomainError add_node!(diagram, DecisionNode("C", ["B", "B"], ["a", "b"]))
@test_throws DomainError add_node!(diagram, DecisionNode("C", ["C", "B"], ["a", "b"]))
@test_throws DomainError add_node!(diagram, DecisionNode("C", [], ["a"]))
@test_throws DomainError add_node!(diagram, ValueNode("A", []))
@test_throws DomainError add_node!(diagram, ValueNode("C", ["B", "B"]))
@test_throws DomainError add_node!(diagram, ValueNode("C", ["C", "B"]))
add_node!(diagram, ChanceNode("C", ["A"], ["a", "b"]))
add_node!(diagram, DecisionNode("D", ["A"], ["c", "d"]))
add_node!(diagram, ValueNode("V", ["A"]))
@test length(diagram.Nodes) == 4
@test isa(diagram.Nodes[3], DecisionNode)
@test isa(diagram.Nodes[4], ValueNode)

@info "Testing generate_arcs!"
diagram = InfluenceDiagram()
@test_throws DomainError generate_arcs!(diagram)
add_node!(diagram, ChanceNode("A", [], ["a", "b"]))
add_node!(diagram, ChanceNode("C", ["A"], ["a", "b", "c"]))
add_node!(diagram, ValueNode("V", ["A", "C"]))
generate_arcs!(diagram)
@test diagram.Names == ["A", "C", "V"]
@test diagram.I_j == [[], Node[1], Node[1, 2]]
@test diagram.States == [["a", "b"], ["a", "b", "c"]]
@test diagram.S == [State(2), State(3)]
@test diagram.C == Node[1, 2]
@test diagram.D == Node[]
@test diagram.V == Node[3]
@test diagram.X == Probabilities[]
@test diagram.Y == Utilities[]

#Non-existent node B
diagram = InfluenceDiagram()
add_node!(diagram, ChanceNode("A", [], ["a", "b"]))
add_node!(diagram, ChanceNode("C", ["B"], ["a", "b", "c"]))
add_node!(diagram, ValueNode("V", ["A", "C"]))
@test_throws DomainError generate_arcs!(diagram)

#Cylic
diagram = InfluenceDiagram()
add_node!(diagram, ChanceNode("A", ["R"], ["a", "b"]))
add_node!(diagram, ChanceNode("R", ["C", "A"], ["a", "b", "c"]))
add_node!(diagram, ChanceNode("C", ["A"], ["a", "b", "c"]))
add_node!(diagram, ValueNode("V", ["A", "C"]))
@test_throws DomainError generate_arcs!(diagram)

#Value node in I_j
diagram = InfluenceDiagram()
add_node!(diagram, ChanceNode("A", ["R"], ["a", "b"]))
add_node!(diagram, ChanceNode("R", ["C", "A"], ["a", "b", "c"]))
add_node!(diagram, ChanceNode("C", ["A", "V"], ["a", "b", "c"]))
add_node!(diagram, ValueNode("V", ["A"]))
@test_throws DomainError generate_arcs!(diagram)

@info "Testing ProbabilityMatrix"
diagram = InfluenceDiagram()
add_node!(diagram, ChanceNode("A", [], ["a", "b"]))
add_node!(diagram, ChanceNode("B", ["A"], ["a", "b", "c"]))
add_node!(diagram, DecisionNode("D", ["A"], ["a", "b", "c"]))
generate_arcs!(diagram)
@test ProbabilityMatrix(diagram, "A") == zeros(2)
@test ProbabilityMatrix(diagram, "B") == zeros(2, 3)
@test_throws DomainError ProbabilityMatrix(diagram, "C")
@test_throws DomainError ProbabilityMatrix(diagram, "D")
X_A = ProbabilityMatrix(diagram, "A")
X_A["a"] = 0.2
@test X_A  == [0.2, 0]
X_A["b"] = 0.9
@test X_A  == [0.2, 0.9]
@test_throws DomainError add_probabilities!(diagram, "A", X_A)
X_A["b"] = 0.8
@test add_probabilities!(diagram, "A", X_A) == [[0.2, 0.8]]
@test_throws DomainError add_probabilities!(diagram, "A", X_A)

@info "Testing UtilityMatrix"
diagram = InfluenceDiagram()
add_node!(diagram, ChanceNode("A", [], ["a", "b"]))
add_node!(diagram, DecisionNode("D", ["A"], ["a", "b", "c"]))
add_node!(diagram, ValueNode("V", ["A", "D"]))
generate_arcs!(diagram)
@test UtilityMatrix(diagram, "V") == fill(Inf, (2, 3))
@test_throws DomainError UtilityMatrix(diagram, "C")
@test_throws DomainError UtilityMatrix(diagram, "D")
Y_V = UtilityMatrix(diagram, "V")
@test_throws DomainError add_utilities!(diagram, "V", Y_V)
Y_V["a", :] = [1, 2, 3]
Y_V["b", "c"] = 4
Y_V["b", "a"] = 5
Y_V["b", "b"] = 6
@test Y_V  == [1 2 3; 5 6 4]
add_utilities!(diagram, "V", Y_V)
@test diagram.Y == [[1 2 3; 5 6 4]]
@test_throws DomainError add_utilities!(diagram, "V", Y_V)

@info "Testing generate_diagram!"
diagram = InfluenceDiagram()
add_node!(diagram, ChanceNode("A", [], ["a", "b"]))
add_node!(diagram, DecisionNode("D", ["A"], ["a", "b", "c"]))
add_node!(diagram, ValueNode("V", ["A", "D"]))
generate_arcs!(diagram)
add_utilities!(diagram, "V", [-1 2 3; 5 6 4])
add_probabilities!(diagram, "A", [0.2, 0.8])
generate_diagram!(diagram)
@test diagram.translation == Utility(0)

@info "Testing positive and negative path utility translations"
generate_diagram!(diagram, positive_path_utility=true)
@test diagram.translation == Utility(2)
@test all(diagram.U(s, diagram.translation) > 0 for s in paths(diagram.S))
generate_diagram!(diagram, negative_path_utility=true)
@test diagram.translation == Utility(-7)
@test all(diagram.U(s, diagram.translation) < 0 for s in paths(diagram.S))


text:
module Countries

# Data obtained from https://datahub.io/core/country-codes
# Consider downloading data in build.jl?

using DelimitedFiles, Reexport
@reexport using Currencies

const (data,headers) = readdlm(joinpath(@__DIR__,"data","country-codes.csv"),',',header=true)

struct Country{T} end

const (nrow,ncol) = size(data)
for i in 1:nrow
    country = data[i,7]
    if length(country) == 2
        @eval begin
            $(Symbol(country)) = Country{Symbol($(country))}()
            name(::Country{Symbol($(country))}) = $(data[i,3])
            code(::Country{Symbol($(country))}) = $(data[i,9])
            capital(::Country{Symbol($(country))}) = $(data[i,29])
            continent(::Country{Symbol($(country))}) = $(data[i,30])
            isdeveloping(::Country{Symbol($(country))}) = isequal($(data[i,32]),"Developing")
            region(::Country{Symbol($(country))}) = $(data[i,50])
            subregion(::Country{Symbol($(country))}) = $(data[i,53])
            Base.show(io::Base.IO,::Country{Symbol($(country))}) = print(io,$(country))
        end
        ccys = Vector{Currencies.Currency}()
        currencies = split(data[i,10],",")
        for currency in currencies
            if isdefined(Currencies,Symbol(currency))
                @eval ccy = Currencies.$(Symbol(currency))
                push!(ccys,ccy)
            end
        end
        @eval currencies(::Country{Symbol($(country))}) = $(ccys)
    end
end

end # module


text:
#
#----------------------------------------------------------------------
export bdfExtK, bdfExtK!
#----------------------------------------------------------------------
""" k'th Order Backward Difference Formula with Extrapolation for integrating

 du = f(t) + g(t)           (1)
 dt

 where f(t) is rapidly changing, and g(t) is considerably slower.
 For example, we treat diffusion implicity, and convection explicityly
 while time-stepping Navier-Stokes as diffusion timescales are an order
 of magnitude faster (pressure is also treated implicitly but that is
 another story) This allows for larger time-steps based on the intertial term.

 Implementation:

 input: t: array of timesteps [t(i+1), t(i), ..., t(i-k+1)] where
           k is the order of the extrapolation

 output: a: extrapolation coefficients for f(t) (k)
         b: interpolation coefficients for du/dt (k+1)

 (1) ⟺   u' - f(t) = g(t)

 where du/dt(t=i+1) = u(i+1)b(1) + ... + u(i-k+1)*b(k+1),
           f(t=i+1) = f(i+1)
           g(t=i+1) = g(i)a(1) + ... + g(i-k+1)*a(k)

"""
function bdfExtK(t;k=3)
    t  = unique(t)
    kk = length(t) - 1
    t1 = t[1]
    t0 = t[2:end]

    a = interpMat(t1,t0)
    b = derivMat(t)[1,:]'

    if(kk < k)
        global a = hcat(a,zeros(1,k-kk))
        global b = hcat(b,zeros(1,k-kk))
    else
        a = a[1:k]
        b = b[1:k+1]
    end

    if(kk==0) # steady state
        global a[1] = 1.0
    end

    return reshape(a,k), reshape(b,k+1)
end
#--------------------------------------#
function bdfExtK!(a::Vector
                 ,b::Vector
                 ,t::Vector)

    k = length(t) - 1
    aa,bb = bdfExtK(t;k=k)

    a .= aa
    b .= bb

    return
end
#----------------------------------------------------------------------
export TimeStepper
#----------------------------------------------------------------------
struct TimeStepper{T,U}
    time::Vector{T}
    bdfA::Vector{T}
    bdfB::Vector{T}

    istep::Array{U,1} # current step
    fstep::Array{U,1} # final step

    dt::Array{T,1} # time step size
    Ti::Array{T,1} # start time
    Tf::Array{T,1} # end time

end
#--------------------------------------#
function TimeStepper(Ti,Tf,dt,k=3) # add flag for variable dt

    time = [Ti] .* ones(k+1)
    bdfA,bdfB = bdfExtK(time)

    istep = [0]
    fstep = [0] # steady

    if(dt != 0)
        global fstep = [Int64(ceil(Tf/dt))]
    end

    return TimeStepper{Float64,Int64}(time,bdfA,bdfB
                                     ,istep,fstep
                                     ,[dt],[Ti],[Tf])
end
#--------------------------------------#
function updateHist!(tstep::TimeStepper)

    @unpack time, istep, dt, bdfA, bdfB = tstep

    updateHist!(time)

    istep  .+= 1
    time[1] += dt[1]
    bdfExtK!(bdfA,bdfB,time)

    return
end
#----------------------------------------------------------------------


text:
"""
Generate a Ricker Wavelet. Reference:
Frequencies of the Ricker wavelet, Yanghua Wang, GEOPHYSICS, VOL. 80, NO. 2
Bandwidth is roughly 1.2 * fpeak for ricker

# Keyword Arguments

* `fqdom::Float64`: dominant frequency
* `tgrid`: time-domain grid
* `tpeak::Float64=tgrid[1]+1.5/fqdom`: the peak of the ricker in time (has a default)
* `tpeak_rand_flag` : randomly place peak of ricker
"""
function ricker(fqdom::Float64,
		tgrid::StepRangeLen;
		tpeak::Float64=tgrid[1]+1.5/fqdom, # using approximate half width of ricker
		attrib::AbstractString="",
		trim_tol::Float64=0.0,
		maxamp::Float64=1.0,
		tpeak_rand_flag=false,
		maxamp_rand_flag=false,
		phase_rand_flag=false,
		)
	(tpeak < tgrid[1]+1.5/fqdom) && error("cannot output Ricker for given tgrid and tpeak")
	(tpeak > tgrid[end]-1.5/fqdom) && error("cannot output Ricker for given tgrid and tpeak")

	if(tpeak_rand_flag)
		tpeak=rand(Uniform(tgrid[1]+1.5/fqdom,tgrid[end]-1.5/fqdom))
	end
	if(maxamp_rand_flag)
		maxamp=rand(Uniform(-1,1))
	end
	isapprox(fqdom,0.0) && error("dominant frequency cannot be zero")

	#! some constants
	pf = (π*π)*(fqdom^2.0)
	nt = length(tgrid)
	δt = step(tgrid)

	# a vector is odd number of samples (nt + 1 corresponds to time zero)
	wav = zero(tgrid);
	# k = (1 - 2* pf * t^2) * Exp[-pf *t^2]
	# Simplify[D[k,t]]
	# FortranForm[Simplify[D[k,t]]]
	if(occursin("[DIFF]",attrib))
			# ricker after a time derivative
			for it = 1:nt
				tsquare = (tgrid[it]-tpeak) * (tgrid[it]-tpeak)
				t       = -1.0 * (tgrid[it]-tpeak)
				wav[it] = (2.0 * pf * t * (-3.0 + 2.0 * pf * tsquare)) * exp(-1.0 * pf * tsquare)
			end
	else
	#! ricker wavelet
		for it = 1:nt
			tsquare = (tgrid[it]-tpeak) * (tgrid[it]-tpeak)
			wav[it] = (1.0 - 2.0 * pf * tsquare) * exp(-1.0e0 * pf * tsquare) * maxamp
		end
	end

	isapprox(maximum(abs.(wav)),0.0) && warn("wavelet is zeros")

	# apply random phase to the ricker wavelet
	if(phase_rand_flag)
		wav=apply_rand_phase(wav)
	end

	if(trim_tol != 0.0)
		return wav[abs(wav).>=trim_tol]
	else
		return wav
	end
end

"""
ormbsy wavelet

* `tperc::Float64=0.0` : the wavelet is tapered in time using this percentage value
"""
function  ormsby(
		fqdom,
		tgrid::StepRangeLen;
		fracbandwidth::Float64=1.2, # using same value as Ricker
		f1::Float64=fqdom-fracbandwidth*0.5*fqdom,
		f4::Float64=fqdom+fracbandwidth*0.5*fqdom,
		f2::Float64=0.25*f4+0.75*f1,
		f3::Float64=0.25*f1+0.75*f4,
		tpeak::Float64=tgrid[1]+1.5/(0.25*(f1+f4)), # using approximate half width
		trim_tol::Float64=0.0,
		tperc::Float64=0.0,
		tpeak_rand_flag=false,
		maxamp_rand_flag=false,
		phase_rand_flag=false,
		)

	# some constants
	A43 = (pi*f4)^2 / (pi*f4 - pi*f3);
	A34 = (pi*f3)^2 / (pi*f4 - pi*f3);
	A21 = (pi*f2)^2 / (pi*f2 - pi*f1);
	A12 = (pi*f1)^2 / (pi*f2 - pi*f1);

	wav = zeros(length(tgrid));

	if(tpeak_rand_flag)
		tpeak=rand(Uniform(tgrid[1]+1.5/(0.25*(f1+f4)),tgrid[end]-1.5/(0.25*(f1+f4))))
	end

	# ormsby wavelet
	for it = 1: length(tgrid)
		t = tgrid[it]-tpeak
		S4 = (sinc(pi*f4*t))^2
		S3 = (sinc(pi*f3*t))^2
		S2 = (sinc(pi*f2*t))^2
		S1 = (sinc(pi*f1*t))^2

		wav[it] =  ((A43 * S4 - A34 * S3) - (A21 * S2 - A12 * S1))
	end

	isapprox(maximum(abs.(wav)),0.0) && warn("wavelet is zeros")

	# normalize
	wav /= maximum(abs.(wav));

	if(maxamp_rand_flag)
		wav *= rand(Uniform(-1,1))
	end

	# apply random phase to the ricker wavelet
	if(phase_rand_flag)
		wav=apply_rand_phase(wav)
	end

	if(trim_tol != 0.0)
		wav = wav[abs(wav).>=trim_tol]
	end

	wav = wav .* Utils.taper(ones(length(wav)),tperc)
	return wav

end

function apply_rand_phase(wav)
	nt=length(wav)
	W=FFTW.rfft(wav);
	θ=rand(Uniform(-Float64(pi),Float64(pi)))
	for i in eachindex(W)
		W[i]=W[i]*complex(cos(θ),sin(θ))
	end
	wav=FFTW.irfft(W,nt)
end

"""
Return a ricker wavelet for a given input `Medium`.

# Arguments
* `mod::Medium` : Medium
* `nλ::Int64=10` : number of wavelengths (P-wave) in the medium
* `tmaxfrac::Float64=1.0` : by default the maximum modelling time is computed using the average velocity and the diagonal distance of the medium,
use this fraction to increase of reduce the maximum time

# Keyword Arguments
* all the keywords arguments of the `ricker` method can be used.
"""
function ricker(mod::Medium, nλ::Int=10, tmaxfrac::Real=1.0, epsilon=inv(sqrt(ndims(mod))); args... )
	@assert(!iszero(mod))
	fqdom, tgrid = get_fqdom_tgrid(mod, nλ, tmaxfrac, epsilon)
	wav=ricker(fqdom, tgrid; args...)
	return wav, tgrid
end

"""
Same as ricker, but return ormsby...
"""
function ormsby(mod::Medium, nλ::Int=10, tmaxfrac::Real=1.0, epsilon=inv(sqrt(ndims(mod))); args... )
	@assert(!iszero(mod))
	fqdom, tgrid = get_fqdom_tgrid(mod, nλ, tmaxfrac, epsilon)
	wav=ormsby(fqdom, tgrid; args...)
	return wav, tgrid
end

"""
Return dominant source frequency, and its temporal grid for a finite-difference simulation, for given number of wavelengths `nλ` in the medium.
The model has `nλ` wavelengths, and the maximum modeling time is determined by `tmaxfrac`.
`epsilon` is Courant number.
"""
function get_fqdom_tgrid( mod::Medium, nλ::Int, tmaxfrac::Real, epsilon)

	# maximum distance (diagnol) the wave travels
	d = sqrt(sum([(m[1]-m[end])^2 for m in mod.mgrid]))

	# dominant wavelength using mod dimensions
	λdom=d*inv(real(nλ))

	# average P velocity
	vavg=mod.ref[:vp]

	fqdom = vavg/λdom

	# use two-way maximum distance to get tmax
	tmax=2.0*d*inv(vavg)*tmaxfrac

	# choose sampling interval to obey max freq of source wavelet
	δmin = minimum(step.(mod.mgrid))
	vmax=try
		sqrt(mod.bounds[:vp][2]^2 + mod.bounds[:vs][2]^2) # see Virieux (1986)
	catch
		mod.bounds[:vp][2]
	end
	δt=epsilon*δmin/vmax

	# check if δt is reasonable
	#(δt > 0.1/fqdom) : error("decrease spatial sampling or nλ")

	tgrid=range(0.0, stop=tmax, step=δt)

	return fqdom, tgrid
end