text:
# ------------------------------------------------------------------------------------------
# # Introduction to DataFrames
# **[Bogumił Kamiński](http://bogumilkaminski.pl/about/), Apr 21, 2018**
# ------------------------------------------------------------------------------------------

using DataFrames
using BenchmarkTools

# ------------------------------------------------------------------------------------------
# ## Performance tips
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Access by column number is faster than by name
# ------------------------------------------------------------------------------------------

x = DataFrame(rand(5, 1000))
@btime x[500];
@btime x[:x500];

# ------------------------------------------------------------------------------------------
# ### When working with data `DataFrame` use barrier functions or type annotation
# ------------------------------------------------------------------------------------------

function f_bad() # this function will be slow
    srand(1); x = DataFrame(rand(1000000,2))
    y, z = x[1], x[2]
    p = 0.0
    for i in 1:nrow(x)
        p += y[i]*z[i]
    end
    p
end

@btime f_bad();

@code_warntype f_bad() # the reason is that Julia does not know the types of columns in `DataFrame`

# solution 1 is to use barrier function (it should be possible to use it in almost any code)
function f_inner(y,z)
   p = 0.0
   for i in 1:length(y)
       p += y[i]*z[i]
   end
   p
end

function f_barrier() # extract the work to an inner function
    srand(1); x = DataFrame(rand(1000000,2))
    f_inner(x[1], x[2])
end

function f_inbuilt() # or use inbuilt function if possible
    srand(1); x = DataFrame(rand(1000000,2))
    dot(x[1], x[2])
end

@btime f_barrier();
@btime f_inbuilt();

# solution 2 is to provide the types of extracted columns
# it is simpler but there are cases in which you will not know these types
function f_typed()
    srand(1); x = DataFrame(rand(1000000,2))
    y::Vector{Float64}, z::Vector{Float64} = x[1], x[2]
    p = 0.0
    for i in 1:nrow(x)
        p += y[i]*z[i]
    end
    p
end

@btime f_typed();

# ------------------------------------------------------------------------------------------
# ### Consider using delayed `DataFrame` creation technique
# ------------------------------------------------------------------------------------------

function f1()
    x = DataFrame(Float64, 10^4, 100) # we work with DataFrame directly
    for c in 1:ncol(x)
        d = x[c]
        for r in 1:nrow(x)
            d[r] = rand()
        end
    end
    x
end

function f2()
    x = Vector{Any}(100)
    for c in 1:length(x)
        d = Vector{Float64}(10^4)
        for r in 1:length(d)
            d[r] = rand()
        end
        x[c] = d
    end
    DataFrame(x) # we delay creation of DataFrame after we have our job done
end

@btime f1();
@btime f2();

# ------------------------------------------------------------------------------------------
# ### You can add rows to a `DataFrame` in place and it is fast
# ------------------------------------------------------------------------------------------

x = DataFrame(rand(10^6, 5))
y = DataFrame(transpose(1.0:5.0))
z = [1.0:5.0;]

@btime vcat($x, $y); # creates a new DataFrame - slow
@btime append!($x, $y); # in place - fast

x = DataFrame(rand(10^6, 5)) # reset to the same starting point
@btime push!($x, $z); # add a single row in place - fastest

# ------------------------------------------------------------------------------------------
# ### Allowing `missing` as well as `categorical` slows down computations
# ------------------------------------------------------------------------------------------

using StatsBase

function test(data) # uses countmap function to test performance
    println(eltype(data))
    x = rand(data, 10^6)
    y = categorical(x)
    println(" raw:")
    @btime countmap($x)
    println(" categorical:")
    @btime countmap($y)
    nothing
end

test(1:10)
test([randstring() for i in 1:10])
test(allowmissing(1:10))
test(allowmissing([randstring() for i in 1:10]))

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
include("matcher.jl")
f = open(input)

#Given a lexem,
function token(lexem,line::Int, col::Int)
	#TODO MUDAR INDEXAÇÃOO (USAR REGISTRO)
	tkn = Token(string(lexem),line,col,"",0)
	#tkn = Dict("lexem"=>string(lexem),"line"=>line, "col"=>col, "categ_nom" => " ","categ_num" => " ")
	matchedcateg =  matchlexem(lexem)

	if (matchedcateg != false)
		tkn.categ_num = matchedcateg
		tkn.categ_nom = tks_names[matchedcateg]
	else
       	tkn.categ_num = Int(ID)
       	tkn.categ_nom = tks_names[Int(ID)]
    end
	return tkn
end

function producer()

	ch = Channel{Any}(1)
	#It says if the head of the reader is actually inside a string. In this mode blank spaces are not ignored
	stringmode = false

	@async begin
		open(input) do file
			l = 0
		 for ln in eachline(file)  #=1:length(lines)=#
		 	l+=1
			comment_mode = false
			chunks = split(ln," ")
			lineq = Queue(Token)
			lexembuff = ""
			col = 0
			chunk = -999
			for cki = 1:length(chunks) #chunk in chunks
				chunk = chunks[cki]

				if(length(chunk) > 0)
					if (chunk[1] == '#')
						break
					end
					i=1
					#for i=1:length(chunk)
					while(i<=length(chunk))
					if (chunk[i] == '"')
						stringmode = !stringmode
					end

						if contains(==,chars_numbers,chunk[i])
							lexembuff = string(lexembuff,chunk[i])
						else

							if length(lexembuff) > 0
								enqueue!(lineq,token(lexembuff,l,col+=1))
								lexembuff = ""
							end
#TODO FAZER UMA FUNÇÃO PARA ELSE-IF
							if (contains(==,separators,chunk[i]))
								enqueue!(lineq,token(string(chunk[i]),l,col+=1))
							elseif chunk[i] == '>'
									if chunk[i+1] == '='
										enqueue!(lineq,token(">=",l,col+=1))
										i+=1
									else
										enqueue!(lineq,token(">",l,col+=1))
									end

							elseif chunk[i] == '<'
									if chunk[i+1] == '='
										enqueue!(lineq,token("<=",l,col+=1))
										i+=1
									else
										enqueue!(lineq,token("<",l,col+=1))
									end

							elseif chunk[i] == '='
									if chunk[i+1] == '='
										enqueue!(lineq,token("==",l,col+=1))
										i+=1
									else
										enqueue!(lineq,token("=",l,col+=1))
									end
							elseif chunk[i] == ':'
									if chunk[i+1] == ':'
										enqueue!(lineq,token("::",l,col+=1))
										i+=1
									end
							end

						end
						i+=1
					end
				end
				if (length(lexembuff) > 0 && !stringmode)
								enqueue!(lineq,token(lexembuff,l,col+=1))
								lexembuff = ""
				elseif stringmode
					lexembuff = string(lexembuff," ")
				end

				if (l == countlines(input) && cki == length(chunks))
					#eoftkn = Dict("lexem"=>"\\EOF","line"=>L, "col"=>col+1, "categ_nom" => "EOF","categ_num" => " ")
					#enqueue!(lineq,token("\\eof",l,col+=1))
				end
			end
			#The pop! command for Queue is not working...lol
			lineq_vector = collect(lineq)
			for i in lineq_vector
				put!(ch,i)
			end
			put!(ch,"eol")
		end
	end

	close(ch)
end
	return ch
end

token(tkn::String) = token(string(tkn),0,0)

ch = producer()

#	@printf("%.3f",x)

g = open(input,"r")
function nextToken()
	tkn = -1
	lex = -1
	col = -1
	nom = -1
	num = -1
	try
		tkn = take!(ch)
		#todo colocar o print no programa de teste e tirar do nextotken

		#println("[",line,", ",col,"] (",categ_num,", ",categ_nom,") {",lexem,"}")

	catch
		#info("End of file reached")
		tkn = Token("\\eof",0,0,"EOF",Int(EOF))
	end
	if typeof(tkn) == Token
		lex = tkn.lexem
		col = tkn.column
		nom = tkn.categ_nom
		num = tkn.categ_num
		#println("\(\"$lex\",$col,$nom,$num\)")
	else
		#info("eol")
		return false
	end
	if tkn.column == 1
		header = @sprintf("%4s  %s\n",lpad(tkn.line,4),readline(g))
		write(tree_file,header)
		println()
		print(header)
	end
	return tkn;
end

#include("sinthatic.jl")

ext: jl

max_stars_repo_name: NaelsonDouglas/namelesscompiler

max_issues_repo_name: NaelsonDouglas/namelesscompiler

max_forks_repo_name: NaelsonDouglas/namelesscompiler


text:
using Test
using PETSc, MPI, LinearAlgebra, SparseArrays

m,n = 20,20
x = randn(n)
V = PETSc.VecSeq(x)

@test norm(x) ≈ norm(V) rtol=10eps()

S = sprand(m,n,0.1) + I
M = PETSc.MatSeqAIJ(S)

@test norm(S) ≈ norm(M) rtol=10eps()

w = M*x
@test w ≈ S*x

ksp = PETSc.KSP(M; ksp_rtol=1e-8, pc_type="jacobi", ksp_monitor=true)
#PETSc.settolerances!(ksp; rtol=1e-8)

@test PETSc.gettype(ksp) == "gmres" # default

pc = PETSc.PC(ksp)
@test PETSc.gettype(pc) == "jacobi"

y = ksp \ w
@test S*y ≈ w  rtol=1e-8

w = M'*x
@test w ≈ S'*x

y = ksp' \ w
@test S'*y ≈ w rtol=1e-8

f!(y,x) = y .= 2 .*x

M = PETSc.MatShell{Float64}(f!,10,10)

x = rand(10)

@test M*x ≈ 2x
@test PETSc.KSP(M) \ x ≈ x/2


text:
"
    _ftIdx(ft, grp)

Finds the index of grp (string) in ft.labels which corresponds to ft.mat.
For Index i for the grp  returned by this function ft[i, :, :] returns the
2D array [[TP, FP], [FN, TN]] for that group.
"
function _ftIdx(ft::FairTensor, grp)
    idx = findfirst(x->x==string.(grp), ft.labels)
    if idx==nothing throw(ArgumentError("$grp not found in the fairness tensor")) end
    return idx
end

# Helper function for calculating TruePositive, FalseNegative, etc.
# If grp is :, it calculates combined value for all groups, else the specified group
_calcmetric(ft::FairTensor, grp, inds...) = typeof(grp)==Colon ? sum(ft[:, inds...]) : sum(ft[_ftIdx(ft, grp), inds...])

(::TruePositive)(ft::FairTensor; grp=:) = _calcmetric(ft, grp, 1, 1)
(::FalsePositive)(ft::FairTensor; grp=:) = _calcmetric(ft, grp, 1, 2)
(::FalseNegative)(ft::FairTensor; grp=:) = _calcmetric(ft, grp, 2, 1)
(::TrueNegative)(ft::FairTensor; grp=:) = _calcmetric(ft, grp, 2, 2)

# The functions true_positive, false_negative, etc are instances of TruePositive, FalseNegative, etc.
# So on using true_positive(fair_tensor) will use the above defined functions
(::TPR)(ft::FairTensor; grp=:) = 1/(1+false_negative(ft; grp=grp)/(1e-15 + true_positive(ft; grp=grp)))
(::TNR)(ft::FairTensor; grp=:) = 1/(1+false_positive(ft; grp=grp)/(1e-15 + true_negative(ft; grp=grp)))
(::FPR)(ft::FairTensor; grp=:) = 1-true_negative_rate(ft; grp=grp)
(::FNR)(ft::FairTensor; grp=:) = 1-true_positive_rate(ft; grp=grp)

(::FDR)(ft::FairTensor; grp=:) = 1/(1+false_positive(ft; grp=grp)/(1e-15 + true_positive(ft; grp=grp)))
(::Precision)(ft::FairTensor; grp=:) = 1 - false_discovery_rate(ft; grp=grp)
(::NPV)(ft::FairTensor; grp=:) =  1/(1+false_negative(ft; grp=grp)/(1e-15 + true_negative(ft; grp=grp)))

"""
    disparity(M, ft; refGrp=nothing, func=/)

Computes disparity for fairness tensor `ft` with respect to an array of metrics `M`.

For any class A and a reference Group B, `disparity = func(metric(A), metric(B))`.
By default `func` is `/` .

A dataframe is returned with disparity values for all combinations of metrics and classes.
It contains a column named labels for the classes and has a column for disparity of each metric in M.
The column names are metric names appended with `_disparity`.

## Keywords

* `refGroup=nothing` : The reference group
* `func=/` : The function used to evaluate disparity. This function should take 2 arguments.
The second argument shall correspond to reference group.

Please note that division by 0 will result in NaN
"""
function disparity(M::Vector{<:Measure}, ft::FairTensor{C}; refGrp=nothing, func=/) where C
    refGrp!==nothing || throw(ArgumentError("Value of reference group needs to be provided"))
    refGrpIdx = _ftIdx(ft, refGrp)
    df = DataFrame(labels=ft.labels)
    for m in M
        colName = string(m) * "_disparity"
        colDisparity = Symbol(colName)
        # TODO : _ftIdx is repeatedly called internally. Instead the value can be stored and reused.
        baseVal = m(ft; grp=refGrp)
        arr = zeros(Float64, C)
        for i in 1:C
            arr[i] = func(m(ft; grp=ft.labels[i]), baseVal)
        end
        df[:, colDisparity] = arr
    end
    return df
end


text:
import ..UncertainValues: UncertainScalarPopulation, UncertainScalarKDE

"""
    RawValues

Indicates that instead of summarising each bin, vectors of raw values should
be returned for a binned resampling.
"""
struct RawValues end

const BINREPR = Union{UncertainScalarKDE, UncertainScalarPopulation, RawValues}

"""
    BinnedResampling(left_bin_edges, n::Int; bin_repr = UncertainScalarKDE)
    BinnedResampling(UncertainScalarKDE, left_bin_edges, n::Int)
    BinnedResampling(UncertainScalarPopulation, left_bin_edges, n::Int)
    BinnedResampling(RawValues, left_bin_edges, n::Int)

Indicates that binned resampling should be performed.

## Fields

- `left_bin_edges`. The left edgepoints of the bins. Either a range or some
    custom type which implements `minimum` and `step` methods.
- `n`. The number of draws. Each point in the dataset is sampled `n` times.
    If there are `m` points in the dataset, then the total number of draws
    is `n*m`.
- `bin_repr`. A type of uncertain value indicating how each bin should be
    summarised (`UncertainScalarKDE` for kernel density estimated distributions
    in each bin, `UncertainScalarPopulation` to represent values in each bin
    as an equiprobable population) or not summarise but return raw values
    falling in each bin (`RawValues`).

## Examples

```julia
using UncertainData

# Resample on a grid from 0 to 200 in steps of 20
grid = 0:20:200

# The number of samples per point in the dataset
n_draws = 10000

# Create the resampling scheme. Use kernel density estimates to distribution
# in each bin.
resampling = BinnedResampling(grid, n_draws, bin_repr = UncertainScalarKDE)

# Represent each bin as an equiprobably population
resampling = BinnedResampling(grid, n_draws, bin_repr = UncertainScalarPopulation)

# Keep raw values for each bin (essentially the same as UncertainScalarPopulation,
# but avoids storing an additional vector of weights for the population members).
resampling = BinnedResampling(grid, n_draws, bin_repr = RawValues)
```
"""
Base.@kwdef struct BinnedResampling{R, B} <: AbstractBinnedUncertainValueResampling where {R <: BINREPR, B}
    bin_repr::Type{R} = UncertainScalarKDE
    left_bin_edges::B
    n::Int
end

BinnedResampling(left_bin_edges, n::Int; bin_repr = UncertainScalarKDE) =
    BinnedResampling(bin_repr, left_bin_edges, n)

function Base.show(io::IO, b::BinnedResampling{R, B}) where {R, B}
    T = typeof(b)
    println(io, "$(T.name){bin_repr: $R, left_bin_edges: $B, n=$(b.n)}")
end

"""
    BinnedWeightedResampling(left_bin_edges, weights, n::Int; bin_repr = UncertainScalarKDE)
    BinnedWeightedResampling(UncertainScalarKDE, left_bin_edges, weights, n::Int)
    BinnedWeightedResampling(UncertainScalarPopulation, left_bin_edges, weights, n::Int)
    BinnedWeightedResampling(RawValues, left_bin_edges, weights, n::Int)

Indicates that binned resampling should be performed, but weighting each
point in the dataset differently.

## Fields

- `left_bin_edges`. The left edgepoints of the bins. Either a range or some
    custom type which implements `minimum` and `step` methods.
- `weights`. The relative probability weights assigned to each point.
- `n`. The total number of draws. These are distributed among the
    points of the dataset according to `weights`.
- `bin_repr`. A type of uncertain value indicating how each bin should be
    summarised (`UncertainScalarKDE` for kernel density estimated distributions
    in each bin, `UncertainScalarPopulation` to represent values in each bin
    as an equiprobable population) or not summarise but return raw values
    falling in each bin (`RawValues`).

## Examples

```julia
using UncertainData, StatsBase

# Resample on a grid from 0 to 200 in steps of 20
grid = 0:20:200

# Assume our dataset has 50 points. We'll assign random weights to them.
wts = Weights(rand(50))

# The total number of draws (on average 1000000/50 = 20000 draws per point
# if weights are equal)
n_draws = 10000000

# Create the resampling scheme. Use kernel density estimates to distribution
# in each bin.
resampling = BinnedWeightedResampling(grid, wts, n_draws, bin_repr = UncertainScalarKDE)

# Represent each bin as an equiprobably population
resampling = BinnedWeightedResampling(grid, wts, n_draws, bin_repr = UncertainScalarPopulation)

# Keep raw values for each bin (essentially the same as UncertainScalarPopulation,
# but avoids storing an additional vector of weights for the population members).
resampling = BinnedWeightedResampling(grid, wts n_draws, bin_repr = RawValues)
```
"""
Base.@kwdef struct BinnedWeightedResampling{R, B, W} <: AbstractBinnedUncertainValueResampling where {R <: BINREPR, B, W}
    bin_repr::Type{R} = UncertainScalarKDE
    left_bin_edges::B
    weights::W
    n::Int
end

BinnedWeightedResampling(left_bin_edges, weights, n::Int; bin_repr = UncertainScalarKDE) =
    BinnedWeightedResampling(bin_repr, left_bin_edges, weights, n)

function Base.show(io::IO, b::BinnedWeightedResampling{R, B, W}) where {R, B, W}
    T = typeof(b)
    println(io, "$(T.name){bin_repr: $R, left_bin_edges: $B, weights: $W, n=$(b.n)}")
end

"""
    BinnedMeanResampling

Binned resampling where each bin is summarised using
the mean of all draws falling in that bin.

## Fields

- `left_bin_edges`. The left edgepoints of the bins. Either a range or some
    custom type which implements `minimum` and `step` methods.
- `n`. The number of draws. Each point in the dataset is sampled `n` times.
    If there are `m` points in the dataset, then the total number of draws
    is `n*m`.

## Examples

```julia
using UncertainData

# Resample on a grid from 0 to 200 in steps of 20
grid = 0:20:200

# The number of samples per point in the dataset
n_draws = 10000

# Create the resampling scheme
resampling = BinnedMeanResampling(grid, n_draws)
```
"""
struct BinnedMeanResampling{B} <: AbstractBinnedSummarisedResampling
    left_bin_edges::B
    n::Int
end

function Base.show(io::IO, b::BinnedMeanResampling{B}) where {B}
    T = typeof(b)
    println(io, "$(T.name){left_bin_edges=$(b.left_bin_edges), n=$(b.n)}")
end

"""
    BinnedMeanWeightedResampling

Binned resampling where each bin is summarised using the mean of all draws
falling in that bin. Points in the dataset are sampled with probabilities
according to `weights`.

## Fields

- `left_bin_edges`. The left edgepoints of the bins. Either a range or some
    custom type which implements `minimum` and `step` methods.
- `weights`. The relative probability weights assigned to each point.
- `n`. The total number of draws. These are distributed among the
    points of the dataset according to `weights`.

## Examples

```julia
using UncertainData, StatsBase

# Resample on a grid from 0 to 200 in steps of 20
grid = 0:20:200

# Assume our dataset has 50 points. We'll assign random weights to them.
wts = Weights(rand(50))

# The total number of draws (on average 1000000/50 = 20000 draws per point
# if weights are equal)
n_draws = 10000000

# Create the resampling scheme
resampling = BinnedMeanWeightedResampling(grid, wts, n_draws)
```
"""
struct BinnedMeanWeightedResampling{B} <: AbstractBinnedSummarisedResampling
    left_bin_edges::B
    weights
    n::Int
end

export
BinnedResampling,
BinnedWeightedResampling,
BinnedMeanResampling,
BinnedMeanWeightedResampling,
RawValues


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

using Markdown
using InteractiveUtils

# ╔═╡ 09dffe9d-11a6-43eb-9e16-53c46943bf7c
begin
  using TikzPictures
  #using LaTeXStrings
end

# ╔═╡ f0d7de38-1ed7-11ec-074a-97be510efaed
md"""
### Definitions
- **_acyclic_**
  - An undirected graph is **_acyclic_** if it does not contain any simple cycle of length ``k \ge 3``
    - The reason for defining in this way is to exclude trivial cycles like ``\langle v_0, v_1, v_0\rangle`` and ``\langle v_0, v_1, v_0, v_1, v_0\rangle``. In other words, we admit as cycle those who have at least three distinct vertices, e.g. $(HTML("<br>")) $(TikzPicture(L"\node at (0,0) [name=v0,shape=circle,draw=white,fill=gray!50,thick] {0};\node at (1,0) [name=v1,shape=circle,draw=white,fill=gray!50] {1};\node  at (1,-0.7) [name=v2,shape=circle,draw=white,fill=gray!50,thick] {2};\draw (v0) to (v1);\draw (v1) to (v2);\draw (v2) to (v0);"))
  - A directed graph is **_acyclic_** if it simply does not contain any cycle (i.e. no self-loops, no cycles of any length).
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
TikzPictures = "37f6aa50-8035-52d0-81c2-5a1d08754b2d"

[compat]
TikzPictures = "~3.4.1"
"""

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

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

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

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

[[Bzip2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "19a35467a82e236ff51bc17a3a44b69ef35185a2"
uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
version = "1.0.8+0"

[[Cairo_jll]]
deps = ["Artifacts", "Bzip2_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "f2202b55d816427cd385a9a4f3ffb226bee80f99"
uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
version = "1.16.1+0"

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

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

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

[[Expat_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b3bfd02e98aedfa5cf885665493c5598c350cd2f"
uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
version = "2.2.10+0"

[[Fontconfig_jll]]
deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "21efd19106a55620a188615da6d3d06cd7f6ee03"
uuid = "a3f928ae-7b40-5064-980b-68af3947d34b"
version = "2.13.93+0"

[[FreeType2_jll]]
deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "87eb71354d8ec1a96d4a7636bd57a7347dde3ef9"
uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
version = "2.10.4+0"

[[Gettext_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046"
uuid = "78b55507-aeef-58d4-861c-77aaff3498b1"
version = "0.21.0+0"

[[Glib_jll]]
deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "7bf67e9a481712b3dbe9cb3dac852dc4b1162e02"
uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
version = "2.68.3+0"

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

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

[[JpegTurbo_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d735490ac75c5cb9f1b00d8b5509c11984dc6943"
uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
version = "2.1.0+0"

[[LZO_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6"
uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
version = "2.10.1+0"

[[LaTeXStrings]]
git-tree-sha1 = "c7f1c695e06c01b95a67f0cd1d34994f3e7db104"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.2.1"

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

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

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

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

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

[[Libffi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "761a393aeccd6aa92ec3515e428c26bf99575b3b"
uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490"
version = "3.2.2+0"

[[Libgcrypt_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll", "Pkg"]
git-tree-sha1 = "64613c82a59c120435c067c2b809fc61cf5166ae"
uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4"
version = "1.8.7+0"

[[Libgpg_error_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c333716e46366857753e273ce6a69ee0945a6db9"
uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"
version = "1.42.0+0"

[[Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.1+1"

[[Libmount_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "9c30530bf0effd46e15e0fdcf2b8636e78cbbd73"
uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
version = "2.35.0+0"

[[Libtiff_jll]]
deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"]
git-tree-sha1 = "340e257aada13f95f98ee352d316c3bed37c8ab9"
uuid = "89763e89-9b03-5906-acba-b20f662cd828"
version = "4.3.0+0"

[[Libuuid_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7f3efec06033682db852f8b3bc3c1d2b0a0ab066"
uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
version = "2.36.0+0"

[[LittleCMS_jll]]
deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pkg"]
git-tree-sha1 = "110897e7db2d6836be22c18bffd9422218ee6284"
uuid = "d3a379c0-f9a3-5b72-a4c0-6bf4d2e8af0f"
version = "2.12.0+0"

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

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

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

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

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

[[OpenJpeg_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libtiff_jll", "LittleCMS_jll", "Pkg", "libpng_jll"]
git-tree-sha1 = "76374b6e7f632c130e78100b166e5a48464256f8"
uuid = "643b3616-a352-519d-856d-80112ee9badc"
version = "2.4.0+0"

[[PCRE_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b2a7af664e098055a7529ad1a900ded962bca488"
uuid = "2f80f16e-611a-54ab-bc61-aa92de5b98fc"
version = "8.44.0+0"

[[Pixman_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b4f5d02549a10e20780a24fce72bea96b6329e29"
uuid = "30392449-352a-5448-841d-b1acce4e97dc"
version = "0.40.1+0"

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

[[Poppler_jll]]
deps = ["Artifacts", "Cairo_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "OpenJpeg_jll", "Pkg", "libpng_jll"]
git-tree-sha1 = "e11443687ac151ac6ef6699eb75f964bed8e1faa"
uuid = "9c32591e-4766-534b-9725-b71a8799265b"
version = "0.87.0+2"

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

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

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

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

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

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

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

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

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

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

[[Tectonic]]
deps = ["Pkg"]
git-tree-sha1 = "acf12eccb390a78653ee805cd527898f01f78a85"
uuid = "9ac5f52a-99c6-489f-af81-462ef484790f"
version = "0.6.1"

[[TikzPictures]]
deps = ["LaTeXStrings", "Poppler_jll", "Requires", "Tectonic"]
git-tree-sha1 = "a08671c0979063a437378f6410bb75a465f3cd1c"
uuid = "37f6aa50-8035-52d0-81c2-5a1d08754b2d"
version = "3.4.1"

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

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

[[XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.12+0"

[[XSLT_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"]
git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a"
uuid = "aed1982a-8fda-507f-9586-7b0439959a61"
version = "1.1.34+0"

[[Xorg_libX11_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
git-tree-sha1 = "5be649d550f3f4b95308bf0183b82e2582876527"
uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"
version = "1.6.9+4"

[[Xorg_libXau_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4e490d5c960c314f33885790ed410ff3a94ce67e"
uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec"
version = "1.0.9+4"

[[Xorg_libXdmcp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fe47bd2247248125c428978740e18a681372dd4"
uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05"
version = "1.1.3+4"

[[Xorg_libXext_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "b7c0aa8c376b31e4852b360222848637f481f8c3"
uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
version = "1.3.4+4"

[[Xorg_libXrender_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "19560f30fd49f4d4efbe7002a1037f8c43d43b96"
uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa"
version = "0.9.10+4"

[[Xorg_libpthread_stubs_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "6783737e45d3c59a4a4c4091f5f88cdcf0908cbb"
uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74"
version = "0.1.0+3"

[[Xorg_libxcb_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
git-tree-sha1 = "daf17f441228e7a3833846cd048892861cff16d6"
uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b"
version = "1.13.0+3"

[[Xorg_xtrans_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "79c31e7844f6ecf779705fbc12146eb190b7d845"
uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10"
version = "1.4.0+3"

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

[[Zstd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "cc4bf3fdde8b7e3e9fa0351bdeedba1cf3b7f6e6"
uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
version = "1.5.0+0"

[[libpng_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c"
uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
version = "1.6.38+0"

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

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

# ╔═╡ Cell order:
# ╠═09dffe9d-11a6-43eb-9e16-53c46943bf7c
# ╟─f0d7de38-1ed7-11ec-074a-97be510efaed
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
export shear, shear!

"""
    shear(arr, Δ, shear_dir_dim=1, shear_dim=2)

Shears an array by the amount of `Δ` pixels via an FFT approach. `Δ` is the relative shift between the
top and bottom row shifted with respect to each other.
`shear_dir_dim` decides the direction of the shear and `shear_dim` is the
second dimension where the shear happens.
There is also `shear!` available.

For complex arrays we use `fft`, for real array we use `rfft`.
"""
function shear(arr::AbstractArray, Δ, shear_dir_dim=1, shear_dim=2; fix_nyquist=false)
    return shear!(copy(arr), Δ, shear_dir_dim, shear_dim, fix_nyquist=fix_nyquist)
end

"""
    shear!(arr, Δ, shear_dir_dim=1, shear_dim=2)

For more details see `shear.`

For complex arrays we can completely avoid large memory allocations.
For real arrays, we need at least allocate on array in the fourier space.
"""
function shear!(arr::AbstractArray{<:Complex, N}, Δ, shear_dir_dim=1, shear_dim=2; fix_nyquist=false) where N
    fft!(arr, shear_dir_dim)

    # stores the maximum amount of shift
    shift = reshape(fftfreq(size(arr, shear_dir_dim)), NDTools.select_sizes(arr, shear_dir_dim))

    apply_shift_strength!(arr, arr, shift, shear_dir_dim, shear_dim, Δ, fix_nyquist)

    # go back to real space
    ifft!(arr, shear_dir_dim)
    return arr
end

function shear!(arr::AbstractArray{<:Real, N}, Δ, shear_dir_dim=1, shear_dim=2; fix_nyquist=false) where N
    p = plan_rfft(arr, shear_dir_dim)
    arr_ft = p * arr

    # stores the maximum amount of shift
    shift = reshape(rfftfreq(size(arr, shear_dir_dim)), NDTools.select_sizes(arr_ft, shear_dir_dim))

    apply_shift_strength!(arr_ft, arr, shift, shear_dir_dim, shear_dim, Δ, fix_nyquist)
    # go back to real space

    # overwrites arr in-place
    ldiv!(arr, p, arr_ft)
    return arr
end

function apply_shift_strength!(arr, arr_orig, shift, shear_dir_dim, shear_dim, Δ, fix_nyquist=false)
    #applies the strength to each slice
    shift_strength = reshape(fftpos(1, size(arr, shear_dim), CenterFirst) .- 0.5, NDTools.select_sizes(arr, shear_dim))

    # do the exp multiplication in place
    e = cispi.(2 .* Δ .* shift .* shift_strength)
    # for even arrays we need to fix real property of highest frequency
    if iseven(size(arr_orig, shear_dir_dim))
        inds = NDTools.slice_indices(axes(e), shear_dir_dim, fft_center(size(arr_orig, shear_dir_dim)))
        r = real.(view(e, inds...))
        if fix_nyquist
            inv_r = 1 ./ r
            inv_r = map(x -> (isinf(x) ? 0 : x), inv_r)
            e[inds...] .= inv_r
        else
            e[inds...] .= r
        end
    end
    arr .*= e
    return arr
end


text:
# load JSON configuration files for channel and instrument parameters
#
# Copyright 2017 Raytheon BBN Technologies

module config
import ..yaml
#import .yaml.load_file_relative

cfg_folder = joinpath(Pkg.dir("QGL"), "cfg")
cfg_path = joinpath(cfg_folder, "cfg_path.txt")
# for simplicity here use single config file

if isdir(cfg_folder) && isfile(cfg_path)
    f = open(cfg_path)
    cfg_yaml = readstring(f)
    close(f)

else
	println("Please provide path to yaml settings file:")
	 cfg_yaml = chomp(readline())

	if !isdir(cfg_folder)
		mkdir(cfg_folder)
	end
	open(cfg_path, "w") do f
        write(f, cfg_yaml)
    end
end

cfg = yaml.load_file_relative(cfg_yaml)
get_qubit_params() = cfg["qubits"]
get_marker_params() = cfg["markers"]
get_edge_params() = cfg["edges"]
get_instrument_params() = cfg["instruments"]
sequence_files_path = cfg["config"]["AWGDir"]

end


text:
export echo, echo!

function echo(args...)
    options, pos_args = parse_args(args)
    _echo(options, pos_args)
end

function _echo(options, args)
    @show options
    @show args
    args
end

echo!(args...) = print(join(echo(args...), " "))


text:
# File for computing Surface-Plasmon-Polariton dispersion relation

" Parallel wavevector for SPP mode on a single interface "
kspp(eps1,eps2,w) = w/c0*sqrt(eps1*eps2/(eps1+eps2))

" Dispersion relation for SPP modes on a film "
function kk(kpar,eps,w)
    return sqrt(kpar^2 - eps*(w/c0)^2)
end

function kspp(kpar,a,eps1,eps2,eps3,w)
    k1 = kk(kpar,eps1,w)
    k2 = kk(kpar,eps2,w)
    k3 = kk(kpar,eps3,w)

    f1 = (k1/eps1 + k2/eps2)/(k1/eps1-k2/eps2)
    f2 = (k1/eps1 + k3/eps3)/(k1/eps1-k3/eps3)

    return exp(-4.0*k1*a)-f1*f2
end

function solve_kpar(a,eps1,eps2,eps3,w)
    f(kpar) = kssp(kpar,a,eps1,eps2,eps3,w)
    return fzero(f,w/c0,order = 16)
end

function dispersion_relation(a,material1 :: OptProp, material2 :: OptProp, material3 :: OptProp )
    wv = collect(linspace(1.0e13,1.0e15,1000))
    eps1 = permittivity.([material1],wv)
    eps2 = permittivity.([material2],wv)
    eps3 = permittivity.([material3],wv)
    kpar = zeros(Float64,1000)

    for i=1:1000
        kpar[i] = solve_kpar(a,eps1[i],eps2[i],eps3[i],wv[i])
    end

    return kpar
end


text:
const SigBitsForFloat64 = 80

elliptick(modulus::Float64) = Float64(elliptick(ArbFloat(modulus, SigBitsForFloat64)))
elliptice(modulus::Float64) = Float64(elliptice(ArbFloat(modulus, SigBitsForFloat64)))
ellipticf(modulus::Float64) = Float64(ellipticf(ArbFloat(modulus, SigBitsForFloat64)))
ellipticpi(nu::Float64, modulus::Float64) = Float64(ellipticpi(ArbFloat(nu, SigBitsForFloat64), ArbFloat(modulus, SigBitsForFloat64)))

elliptic_k(modulus::Float64) = Float64(elliptic_k(ArbFloat(modulus, SigBitsForFloat64)))
elliptic_e(modulus::Float64) = Float64(elliptic_e(ArbFloat(modulus, SigBitsForFloat64)))
elliptic_f(modulus::Float64) = Float64(elliptic_f(ArbFloat(modulus, SigBitsForFloat64)))
elliptic_pi(nu::Float64, modulus::Float64) = Float64(elliptic_pi(ArbFloat(nu, SigBitsForFloat64), ArbFloat(modulus, SigBitsForFloat64)))

ellipticrf(x::Float64, y::Float64, z::Float64) =
    Float64(ellipticrf(ArbFloat(x, SigBitsForFloat64), ArbFloat(y, SigBitsForFloat64), ArbFloat(z, SigBitsForFloat64)))
ellipticrg(x::Float64, y::Float64, z::Float64) =
    Float64(ellipticrg(ArbFloat(x, SigBitsForFloat64), ArbFloat(y, SigBitsForFloat64), ArbFloat(z, SigBitsForFloat64)))
ellipticrj(x::Float64, y::Float64, z::Float64) =
    Float64(ellipticrj(ArbFloat(x, SigBitsForFloat64), ArbFloat(y, SigBitsForFloat64), ArbFloat(z, SigBitsForFloat64)))


text:
using SCF
using Documenter

makedocs(;
    modules=[SCF],
    authors="Stefanos Carlström <stefanos.carlstrom@gmail.com>",
    repo="https://github.com/JuliaAtoms/SCF.jl/blob/{commit}{path}#{line}",
    sitename="SCF.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="juliaatoms.org",
        assets=String["assets/latex.js"],
        mathengine = Documenter.MathJax()
    ),
    pages=[
        "Home" => "index.md",
        "Quantum Systems" => "quantum_systems.md",
        "Fock Operators" => "fock_operators.md",
        "Self-Consistent Iteration" => "self_consistent_iteration.md"
    ],
    doctest=false,
)

deploydocs(;
    repo="github.com/JuliaAtoms/SCF.jl",
)


text:
function init_constraint!(
    com::CS.CoM,
    constraint::ReifiedConstraint,
    fct::Union{MOI.VectorOfVariables,VAF{T}},
    set::RS;
    active = true
) where {A,T<:Real,RS<:ReifiedSet{A}}
    inner_constraint = constraint.inner_constraint
    complement_constraint = constraint.complement_constraint

    variables = com.search_space
    rei_vidx = constraint.indices[1]
    rei_var = variables[rei_vidx]

    feasible = init_constraint!(
        com,
        inner_constraint,
        inner_constraint.fct,
        inner_constraint.set;
    )
    # map the bounds to the indicator constraint
    constraint.bound_rhs = inner_constraint.bound_rhs
    # the reified variable can't be activated if inner constraint is infeasible
    if !feasible && active
        !rm!(com, rei_var, Int(constraint.activate_on)) && return false
    end

    if complement_constraint !== nothing
        feasible = init_constraint!(
            com,
            complement_constraint,
            complement_constraint.fct,
            complement_constraint.set;
        )
    end
    # still feasible
    return true
end

function prune_constraint!(
    com::CS.CoM,
    constraint::ReifiedConstraint,
    fct::Union{MOI.VectorOfVariables,VAF{T}},
    set::RS;
    logs = true,
) where {A,T<:Real,RS<:ReifiedSet{A}}
    # 1. if the inner constraint is solved then the reified variable can be set to activate_on
    # 2. if the inner constraint is infeasible the reified variable can be set to !activate_on
    # 3. if the reified constraint is active then prune can be called for the inner constraint
    # 4. if the reified constraint is fixed to inactive one can complement prune

    variables = com.search_space
    rei_vidx = constraint.indices[1]
    inner_constraint = constraint.inner_constraint
    complement_constraint = constraint.complement_constraint
    activate_on = Int(constraint.activate_on)
    activate_off = activate_on == 1 ? 0 : 1

    # 3
    if issetto(variables[rei_vidx], activate_on)
        !activate_inner!(com, constraint) && return false
        return prune_constraint!(
            com,
            inner_constraint,
            inner_constraint.fct,
            inner_constraint.set,
        )
    # 4
    elseif issetto(variables[rei_vidx], activate_off) && complement_constraint !== nothing
        !activate_complement_inner!(com, constraint) && return false
        return prune_constraint!(
            com,
            complement_constraint,
            complement_constraint.fct,
            complement_constraint.set,
        )
    # 1
    elseif is_constraint_solved(
        com,
        inner_constraint,
        inner_constraint.fct,
        inner_constraint.set,
    )
        !fix!(com, variables[rei_vidx], activate_on) && return false
    # 2
    elseif is_constraint_violated(
        com,
        inner_constraint,
        inner_constraint.fct,
        inner_constraint.set,
    )
        !fix!(com, variables[rei_vidx], activate_off) && return false
    end
    return true
end

function still_feasible(
    com::CS.CoM,
    constraint::ReifiedConstraint,
    fct::Union{MOI.VectorOfVariables,VAF{T}},
    set::RS,
    vidx::Int,
    val::Int,
) where {A,T<:Real,RS<:ReifiedSet{A}}
    inner_constraint = constraint.inner_constraint
    variables = com.search_space
    activate_on = Int(constraint.activate_on)
    activate_off = activate_on == 1 ? 0 : 1
    rei_vidx = constraint.indices[1]

    # if currently activated check if inner constraint is feasible
    if (vidx == rei_vidx && val == activate_on) || issetto(variables[rei_vidx], activate_on)
        # check if already violated
        violated = is_constraint_violated(
            com,
            inner_constraint,
            inner_constraint.fct,
            inner_constraint.set,
        )
        violated && return false
        return still_feasible(
            com,
            inner_constraint,
            inner_constraint.fct,
            inner_constraint.set,
            vidx,
            val,
        )
    end
    # if inner constraint can't be activated it shouldn't be solved
    if (vidx == rei_vidx && val == activate_off) || issetto(variables[rei_vidx], activate_off)
        if all(i == vidx || isfixed(com.search_space[i]) for i in inner_constraint.indices)
            values = [
                i == vidx ? val : value(com.search_space[i])
                for i in inner_constraint.indices
            ]
            return !is_constraint_solved(
                inner_constraint,
                inner_constraint.fct,
                inner_constraint.set,
                values,
            )
        end
    end
    return true
end

function is_constraint_solved(
    constraint::ReifiedConstraint,
    fct::Union{MOI.VectorOfVariables,VAF{T}},
    set::RS,
    values::Vector{Int},
) where {A,T<:Real,RS<:ReifiedSet{A}}
    activate_on = Int(constraint.activate_on)
    inner_constraint = constraint.inner_constraint
    return is_constraint_solved(
        inner_constraint,
        inner_constraint.fct,
        inner_constraint.set,
        values[2:end],
    ) == (values[1] == activate_on)
end

"""
    is_constraint_violated(
        com::CoM,
        constraint::ReifiedConstraint,
        fct::Union{MOI.VectorOfVariables,VAF{T}},
        set::RS
    )  where {A,T<:Real,RS<:ReifiedSet{A}}

Checks if the constraint is violated as it is currently set. This can happen inside an
inactive reified or indicator constraint.
"""
function is_constraint_violated(
    com::CoM,
    constraint::ReifiedConstraint,
    fct::Union{MOI.VectorOfVariables,VAF{T}},
    set::RS,
) where {A,T<:Real,RS<:ReifiedSet{A}}
    if all(isfixed(var) for var in com.search_space[constraint.indices])
        return !is_constraint_solved(
            constraint,
            fct,
            set,
            [CS.value(var) for var in com.search_space[constraint.indices]],
        )
    end

    reified_vidx = constraint.indices[1]
    reified_var = com.search_space[reified_vidx]
    if isfixed(reified_var) && CS.value(reified_var) == Int(constraint.activate_on)
        inner_constraint = constraint.inner_constraint
        return is_constraint_violated(
            com,
            inner_constraint,
            inner_constraint.fct,
            inner_constraint.set,
        )
    end
    return false
end

function changed!(com::CS.CoM, constraint::ReifiedConstraint, fct, set)
    inner_constraint = constraint.inner_constraint
    changed!(com, inner_constraint, inner_constraint.fct, inner_constraint.set)
    if constraint.complement_constraint !== nothing
        complement_constraint = constraint.complement_constraint
        changed!(com, complement_constraint, complement_constraint.fct, complement_constraint.set)
    end
end


text:
struct BoolOr <: PushBinaryOp{Bool} end
execute_op(i::BoolOr, a::Bool, b::Bool) = (a||b)
opname(::Type{BoolOr}) = "or"

struct Eq{I} <: PushInstructionA2tB1{I,Bool} end
execute_op(i::Eq{I}, a::I, b::I) where {I} = (b==a)
oldopname(::Type{<:Eq}) = "="

struct Lt{I} <: PushInstructionA2tB1{I,Bool} end
execute_op(i::Lt{I}, a::I, b::I) where {I} = (b<a)
oldopname(::Type{<:Lt}) = "<"

struct Gt{I} <: PushInstructionA2tB1{I,Bool} end
execute_op(i::Gt{I}, a::I, b::I) where {I} = (b>a)
oldopname(::Type{<:Gt}) = ">"


text:
function AIC(x,y)
    N, k = size(x)
    @assert N == length(y)
    return N*log(RSS(x,y)/N) + 2k
end

function BIC(x,y)
    N,k = size(x)
    @assert N == length(y)
    return N*log(RSS(x,y)/N) + k*log(N)
end

function AIC_min(x,y)
    _,p = size(x)
    aic_min = Array{Float64}(undef,p)
    aic_min_com = Array{Any}(undef,p)
    for i in 1:p
        combi = combinations(1:p,i) |> collect
        aic_min[i],aic_index = [AIC(x[:,c],y) for c in combi] |> findmin
        aic_min_com[i] = combi[aic_index]
    end
    return aic_min, aic_min_com
end

function BIC_min(x,y)
    _,p = size(x)
    bic_min = Array{Float64}(undef,p)
    bic_min_com = Array{Any}(undef,p)
    for i in 1:p
        combi = combinations(1:p,i) |> collect
        bic_min[i],bic_index = [BIC(x[:,c],y) for c in combi] |> findmin
        bic_min_com[i] = combi[bic_index]
    end
    return bic_min, bic_min_com
end

function AR2(x,y)
    N,k = size(x)
    @assert N == length(y)
    return 1 - RSS(x,y)/(N-k-1)/TSS(y)*(N-1)
end

function AR2_max(x,y)
    _,p = size(x)
    ar2_max = Array{Float64}(undef,p)
    ar2_max_com = Array{Any}(undef,p)
    for i in 1:p
        combi = combinations(1:p,i) |> collect
        ar2_max[i],ar2_index = [AR2(x[:,c],y) for c in combi] |> findmax
        ar2_max_com[i] = combi[ar2_index]
    end
    return ar2_max, ar2_max_com
end


text:
module TestConstructor

using Test
using DiffCheckNumbers

constant = DiffCheck{Float64, ()}(1.0, NamedTuple{(), NTuple{0, Float64}}(( )))
x = DiffCheck{Float64, (:x, )}(1.0, NamedTuple{(:x, ), NTuple{1, Float64}}((1.0, )))
xy = DiffCheck{Float64, (:y, :x)}(1.0, NamedTuple{(:y, :x), NTuple{2, Float64}}((1.0, 1.0)))
z = DiffCheck{Float64, (:z, )}(1.0, NamedTuple{(:z, ), NTuple{1, Float64}}((1.0, )))

#tests DiffCheck{T}(x::T) where T<:ReComp = DiffCheck{T, ()}(T(x), NamedTuple{(), NTuple{0,T}}(()))
@test DiffCheck(1.0) === constant
#tests DiffCheck(x::T, y::Symbol) where {T<:ReComp} = DiffCheck{T, (y, )}(x, NamedTuple{(y, ), NTuple{1, T}}((1,)))
@test DiffCheck(1.0, :x) === x
#tests DiffCheck(x::T, y::Array{Symbol}) where {T<:ReComp}
@test DiffCheck(1.0, [:y, :x]) === xy
# please be careful with the order of the symbols when comparing!!!
complex_z = DiffCheck(Complex{Float64}(1.0), :z)
complex_w = DiffCheck{Complex{Float64}}(0.0, :w)
complex_zw = DiffCheck(1.0+im*0, [:z, :w])
complex_zw2 = DiffCheck{Complex{Float64}}(1.0, [:z, :w])
###

@test complex_z === DiffCheck{Complex{Float64}, (:z, )}(1.0, NamedTuple{(:z, ), NTuple{1, Complex{Float64}}}((1.0, )))
@test complex_w === DiffCheck{Complex{Float64}, (:w, )}(0.0, NamedTuple{(:w, ), NTuple{1, Complex{Float64}}}((1.0, )))
@test complex_zw === DiffCheck{Complex{Float64}, (:w, :z)}(1.0, StaticTuple{Complex{Float64},(:w, :z),2}((1.0, 1.0)))
@test complex_zw2 === DiffCheck{Complex{Float64}, (:w, :z)}(1.0, StaticTuple{Complex{Float64},(:w, :z),2}((1.0, 1.0)))

end


text:
using ApproxFun, Test, StaticArrays, SpecialFunctions, LinearAlgebra
    import ApproxFun: testbandedbelowoperator, testbandedoperator, testspace, testtransforms, Vec,
                        maxspace, NoSpace, hasconversion, testfunctional,
                        jacobip, reverseorientation, ReverseOrientation

@testset "Jacobi" begin
    @testset "Basic" begin
        @test jacobip(0:5,2,0.5,0.1) ≈ [1.,0.975,-0.28031249999999996,-0.8636328125,-0.0022111816406250743,0.7397117980957031]

        testspace(Jacobi(.5,2.);haslineintegral=false)

        f=Fun(exp,Jacobi(.5,2.))
        @test f(.1) ≈ exp(.1)

        f=Fun(x->cos(100x),Jacobi(.5,2.124),500)
        @test f(.1) ≈ cos(100*.1)

        sp=Jacobi(.5,2.124)
        @time f=Fun(exp,sp)
        sp2=Jacobi(1.5,2.124)
        f2=Fun(exp,sp2)
        sp3=Jacobi(1.5,3.124)
        f3=Fun(exp,sp3)
        sp4=Jacobi(2.5,4.124)
        f4=Fun(exp,sp4)
        @test norm((Fun(f,sp2)-f2).coefficients)<100eps()
        @test norm((Fun(f,sp3)-f3).coefficients)<100eps()
        @test norm((Fun(f,sp4)-f4).coefficients)<200eps()

        m=20
        @time testtransforms(JacobiWeight(0.,m,Jacobi(0.,2m+1)))
        f=Fun(x->((1-x)/2).^m.*exp(x),JacobiWeight(0.,m,Jacobi(0.,2m+1)))
        @test abs(f(.1)-(x->((1-x)/2).^m.*exp(x))(.1))<10eps()

        m=10
        @time f=Fun(x->besselj(m,m*(1-x)),JacobiWeight(0.,m,Jacobi(0.,2m+1)))
        @test f(0.) ≈ besselj(m,m)
    end

    @testset "Conversion" begin
        testtransforms(Jacobi(-0.5,-0.5))
        @test norm(Fun(Fun(exp),Jacobi(-.5,-.5))-Fun(exp,Jacobi(-.5,-.5))) < 300eps()

        x=Fun(identity)
        ri=0.5/(1-x)
        @test ((1-x)/2 .* Fun(exp,JacobiWeight(0.,0.,Jacobi(0.,1.))))(.1) ≈ (1-.1)./2*exp(.1)

        @test ((1-x)/2 .* Fun(exp,JacobiWeight(0.,0.,Jacobi(0.,1.))))(.1) ≈ (1-.1)./2*exp(.1)

        @test (ri*Fun(exp,JacobiWeight(0.,0.,Jacobi(0.,1.))))(.1) ≈ .5/(1-.1)*exp(.1)
    end

    @testset "Derivative" begin
        D=Derivative(Jacobi(0.,1.,Segment(1.,0.)))
        @time testbandedoperator(D)

        S=JacobiWeight(0.,0.,Jacobi(0.,1.,Segment(1.,0.)))
        D=Derivative(S)
        testbandedoperator(D)

        f=Fun(exp,domainspace(D))
        @test (D*f-f).coefficients|>norm < eps(100000.)
        @test (f'-f).coefficients|>norm < eps(100000.)
        @test (D^2*f-f).coefficients|>norm < eps(100000000.)
        @test (D*(D*f)-f).coefficients|>norm < eps(100000000.)

        S=JacobiWeight(1,1,Ultraspherical(1))

        f=Fun(S,[1.,2.,3.])
        @test (Derivative(S,2)*f)(0.1) ≈ f''(0.1)
    end

    @testset "Jacobi multiplication" begin
        x=Fun(identity,Jacobi(0.,0.))
        f=Fun(exp,Jacobi(0.,0.))

        @test (x*f)(.1) ≈ .1exp(.1)

        x=Fun(identity,Jacobi(12.324,0.123))
        f=Fun(exp,Jacobi(0.,0.))

        @test (x*f)(.1) ≈ .1exp(.1)

        x=Fun(identity,Jacobi(12.324,0.123))
        f=Fun(exp,Jacobi(0.590,0.213))

        @test (x*f)(.1) ≈ .1exp(.1)

        g=Fun(cos,Jacobi(12.324,0.123))
        f=Fun(exp,Jacobi(0.590,0.213))

        @test (g*f)(.1) ≈ cos(.1)*exp(.1)
    end

    @testset "Jacobi integrate and sum" begin
        testtransforms(Legendre(0..2))
        @test sum(Fun(exp,Legendre(0..2))) ≈ sum(Fun(exp,0..2))

        a=Arc(0.,.1,0.,π/2)
        g=Fun(exp,Legendre(a))

        @test sum(g) ≈ sum(Fun(exp,a))
    end

    @testset "special derivative" begin
        x=Fun()
        f=exp(x)*sqrt(1-x^2)
        D=Derivative(WeightedJacobi(.5,.5))

        testtransforms(WeightedJacobi(.5,.5))
        testbandedoperator(D)

        @time g=(D*Fun(f,domainspace(D)))
        @test f'(0.1) ≈ g(0.1)
    end

    @testset "implementation of conversion between Chebyshev and Jacobi spaces using FastTransforms" begin
        f = Fun(x->cospi(1000x))
        g = Fun(f,Legendre())
        h = Fun(g,Chebyshev())
        @test norm(f.coefficients-h.coefficients,Inf) < 100eps()
        @time h = Fun(h,Legendre())
        @test norm(g.coefficients-h.coefficients,Inf) < 1000eps()
    end

    @testset "==" begin
        @test WeightedJacobi(0.1,0.2) == WeightedJacobi(0.1+eps(),0.2)
    end

    @testset "subspace bug" begin
        f=Fun(WeightedJacobi(0.1,0.2),rand(10))  # convert to Legendre expansion

        g=(f|(2:ApproxFun.∞))

        @test ApproxFun.coefficients(g.coefficients,space(g),ApproxFun.canonicalspace(g))[1] ==0.
        @test norm((Fun(g,space(f))|(2:ApproxFun.∞)-g).coefficients) < 10eps()
    end

    @testset "conversion for non-compatible paramters" begin
        S=Jacobi(1.2,0.1)
        x=Fun()

        p=(S,k)->Fun(S,[zeros(k);1.])
        n=1;
        @test norm(x*p(S,n-1)-(ApproxFun.recα(Float64,S,n)*p(S,n-1) + ApproxFun.recβ(Float64,S,n)*p(S,n))) < 10eps()
    end

    @testset "Log with squareroot singularities" begin
        a = 1.0; b = 2.0+im
        d = Segment(a,b)
        z = Fun(d)

        f = real(exp(z) / (sqrt(z-a)*sqrt(b-z)))
        S=space(f)
        x=4.0+2im;
        @test linesum(f*log(abs(x-z))) ≈ 13.740676344264614
    end

    @testset "Line sum for legendre" begin
        x = Fun(Legendre())
        @test sum(x+1) ≈ linesum(x+1)
        x=Fun(Legendre(Segment(2,1)))
        @test sum(x+1) ≈ -linesum(x+1)

        x=Fun(Segment(1,1+im))
        @test sum(x+1) ≈ im*linesum(x+1)

        x=Fun(Legendre(Segment(1,1+im)))
        @test sum(x+1) ≈ im*linesum(x+1)

        x=Fun(Legendre(Segment(im,1)))
        @test sum(x+1) ≈ (1-im)/sqrt(2)*linesum(x+1)
    end

    @testset "vector valued case" begin
        f=Fun((x,y)->real(exp(x+im*y)), Legendre(Segment(Vec(0.,0.),Vec(1.,1.))))
        @test f(0.1,0.1) ≈ real(exp(0.1+0.1im))
    end

    @testset "integer, float mixed" begin
        C=Conversion(Legendre(),Jacobi(1,0))
        testbandedoperator(C)
    end

    @testset "Addition of piecewise Legendre bug" begin
        f = Fun(exp,Legendre())
        f1 = Fun(exp,Legendre(-1..0))
        f2 = Fun(exp,Legendre(0..1))
        fp = f1+f2
        @test space(fp) isa PiecewiseSpace
        @test fp(0.1) ≈ exp(0.1)
        @test fp(0.) ≈ exp(0.)
        @test fp(-0.1) ≈ exp(-0.1)
    end

    @testset "JacobiWeight cumsum bug Issue #557" begin
        x = Fun(0.0..1.0)
        ν = 2
        @time f = x^(ν/2-1) * exp(-x/2) # 0.05s
        @test cumsum(f)' ≈ f
        @test cumsum(f)(1.0) ≈ 0.7869386805747332 # Mathematic

        x = Fun(Ray())
        ν = 2
        @time f = x^(ν/2-1) * exp(-x/2) # 0.05s
        @test cumsum(f)' ≈ f
        @test cumsum(f)(1.0) ≈ 0.7869386805747332
    end

    @testset "Jacobi–Chebyshev conversion" begin
        a,b = (Jacobi(-0.5,-0.5), Legendre())
        @test maxspace(a,b) == NoSpace()
        @test union(a,b) == a
        @test !hasconversion(a,b)

        a,b = (Chebyshev(), Legendre())
        @test maxspace(a,b) == NoSpace()
        @test union(a,b) == Jacobi(-0.5,-0.5)
        @test !hasconversion(a,b)
    end

    @testset "WeightedLaguerre cumsum" begin
        α = 2.7
        f = Fun(WeightedLaguerre(α), [1.0]);
        f = Fun(f, JacobiWeight(α,0,Chebyshev(Ray())));
        g = integrate(f)
        g(3.0) - cumsum(Fun(x -> f(x), 0..6))(3.0)
    end

    @testset "Definite integral" begin
        for S in (WeightedJacobi(0,0), JacobiWeight(0,0, Legendre(1.1..2.3)), Legendre())
            B = DefiniteIntegral(S)
            testfunctional(B)
            @test ApproxFun.rowstop(B,1) == 1
            B[1] == arclength(domain(S))
            f = Fun(exp, S)
            B*f == sum(Fun(exp,domain(S)))
        end
    end

    @testset "Reverse orientation" begin
        S = Jacobi(0.1,0.2)

        @test_throws ArgumentError Conversion(S, Jacobi(1.1,1.2,0..1))

        f = Fun(S, randn(10))
        @test f(0.1) ≈ (ReverseOrientation(S)*f)(0.1) ≈ reverseorientation(f)(0.1)
        @test rangespace(ReverseOrientation(S)) == space(reverseorientation(f)) ==
                    Jacobi(0.2,0.1,Segment(1,-1))

        R = Conversion(S, reverseorientation(S))
        @test (R*f)(0.1) ≈ f(0.1) ≈ reverseorientation(f)(0.1)

        S = Legendre()
        f = Fun(S, randn(10))
        @test f(0.1) ≈ (ReverseOrientation(S)*f)(0.1) ≈ reverseorientation(f)(0.1)
        @test rangespace(ReverseOrientation(S)) == space(reverseorientation(f)) ==
                    Legendre(Segment(1,-1))

        R = Conversion(S, reverseorientation(S))
        @test rangespace(R) == reverseorientation(S) ==
            space(reverseorientation(f))
        @test f(0.1) ≈ (R*f)(0.1) ≈ reverseorientation(f)(0.1)
    end
end


text:
using Colors
set_palette([c"#ca624b", c"#f7dba7", c"#daddd8", c"#61d095", c"#c7d59f"])
for i in 1:5
    plot!(randn(10), lw=0.05)
end


text:
using jInv.Mesh;
using ForwardHelmholtz
import FactoredEikonalFastMarching.getAnalytic2DeikonalSolution
#import EikonalInv.expandModelNearest
#import EikonalInv.addAbsorbingLayer
plotting = false;
if plotting
	using PyPlot;
	close("all");
end

include("My_sub2ind.jl");

n_tup = (512,512)
m = ones(n_tup);
n = collect(n_tup)
Minv = getRegularMesh([0.0,10.0,0.0,10.0],n-1);

pad = 35;
f = 0.5;
fmax = 2.5;
f = fmax;
w = 2*pi*f

src = div.(n,2);
q = zeros(Complex128,n_tup)
q[My_sub2ind(n,src)] = 1/(Minv.h[1]^2);

r = getAnalytic2DeikonalSolution(n,Minv.h,src)[1];
a = 1./(sqrt.(r+Minv.h[1]));
s_analytic = real(exp.(1im*pi/4).*a.*exp.(1im*w.*m.*r));

alpha = 0.5*2*pi;
maxOmega = getMaximalFrequency(m,Minv);
ABLamp = maxOmega;

gamma = getABL(Minv,false,[pad,pad],ABLamp);

H = GetHelmholtzOperator(Minv, m, w, gamma+alpha, false,false);

s1 = H\q[:];
s1 = real(reshape(s1,n_tup));

H = 0;
H = GetHelmholtzOperator(Minv, m, w-1im*alpha./2, gamma[:], false,false);

s2 = H\q[:];
s2 = real(reshape(s2,n_tup));

s_analytic = s_analytic./maximum(s_analytic);
s_analytic.*=exp.(-(alpha/2).*m.*r)
s_analytic.*=maximum(s1);

if plotting
	figure();
	imshow(s1');title("Helmholtz Solution");
	figure();
	imshow(s_analytic');title("Analytic Solution");
end

if plotting
	figure();
	plot(s1[src[1],:],"b");
	plot(s2[src[1],:],"r");
	plot(s_analytic[src[1],:],"g");
	legend(("s1","s2","Analytic"));

end


text:
using NCTiles, NCDatasets, MeshArrays, Dates

function writetestfile(fname,field,package)

    if isfile(fname)
        rm(fname)
    end

    README = ["This file is written as a test for NCTiles.jl."]
    if isa(field,Dict) # TileData Example
        savenames = joinpath.(ncfiltile2d*".".*lpad.(string.(1:field["data2d"].values.numtiles),4,"0").*".nc")
        dims = unique(vcat([field[v].dims for v in keys(field)]...))
        datasets = [createfile(savenames[tidx],field,README,ntile = length(savenames), itile = tidx) for tidx in 1:length(savenames)]

        ds = [x[1] for x in datasets]
        fldvars = [x[2] for x in datasets]
        #dims = [x[3] for x in datasets]

        for k in keys(field)
            if isa(field[k].values,TileData)
                addData(fldvars,field[k])
            else
                tmpfldvars = [fv[findfirst(isequal(k),name.(fv))] for fv in fldvars]
                addData.(tmpfldvars,Ref(field[k]))
            end
        end

        for dim in dims
            addDimData.(ds,Ref(dim))
        end
        close.(ds)

    else
        ds,fldvar,dimlist = createfile(fname,field,README)

        addData(fldvar,field)
        if package == NCDatasets
            addDimData.(Ref(ds),field.dims)
        end
        # Close the file
        close(ds)
    end
end

function testfile(fname,checkfld)
    pass = true
    ds = Dataset(fname)
    varcheck = Dict()
    varcheck["dims"] = Dict()
    for d in checkfld.dims
        fildim = ds[d.name]
        dimcheck = Dict()
        # check size
        dimcheck["size"] = length(fildim) == length(d.values)
        if isa(fildim[1],DateTime)
            dimvals = [Hour.(ds["tim"][t] .- DateTime(replace(d.units,"days since "=>""),"yyyy-mm-dd HH:MM:SS")).value for t in 1:length(ds["tim"])]/24
            dimcheck["values"] = d.values == dimvals
        else
            dimcheck["values"] = d.values == fildim
            if length(size(d.values)) > 1
                dimcheck["values"] = d.values[:,1] == fildim
            end
        end
        pass = pass && dimcheck["size"] && dimcheck["values"]
        # check attributes
        for k in keys(fildim.attrib)
            if k == "units"
                dimcheck["units"] = fildim.attrib[k] == d.units
                pass = pass && dimcheck["units"]
            else
                dimcheck[k] = fildim.attrib[k] == d.atts[k]
                pass = pass && dimcheck[k]
            end
        end

        varcheck["dims"][d.name] = dimcheck
    end

    dsvar = ds[checkfld.name]

    for k in keys(dsvar.attrib)
        if k == "units"
            varcheck["units"] = dsvar.attrib[k] == checkfld.units
            pass = pass && varcheck["units"]
        else
            varcheck[k] = dsvar.attrib[k] == checkfld.atts[k]
            pass = pass && varcheck[k]
        end
    end

    # Check values of test data
    varcheck["values"] = false
    if isa(checkfld.values,Array)
        varcheck["values"] = true
        if isa(checkfld.values[1],Array)
            if NCTiles.hastimedim(checkfld)
                for t in 1:length(checkfld.values)
                    if length(size(dsvar)) == 3
                        varcheck["values"] = varcheck["values"] && checkfld.values[t] == dsvar[:,:,t]
                    elseif length(size(dsvar)) == 4
                        varcheck["values"] = varcheck["values"] && checkfld.values[t] == dsvar[:,:,:,t]
                    elseif length(size(dsvar)) == 2
                        varcheck["values"] = varcheck["values"] && checkfld.values[t] == dsvar[:,t]
                    end
                end
            end
        else
            varcheck["values"] = varcheck["values"] && checkfld.values == dsvar
        end
    elseif isa(checkfld.values,BinData)
        varcheck["values"] = true
        for t in 1:length(checkfld.values.fnames)
            testdata = readbin(checkfld.values,t)
            if length(size(dsvar)) == 3
                varcheck["values"] = varcheck["values"] && testdata == dsvar[:,:,t]
            elseif length(size(dsvar)) == 4
                varcheck["values"] = varcheck["values"] && testdata == dsvar[:,:,:,t]
            elseif length(size(dsvar)) == 2
                varcheck["values"] = varcheck["values"] && testdata == dsvar[:,t]
            end
        end
    elseif isa(checkfld.values,NCData)
        varcheck["values"] = true
        testfld = Dataset(checkfld.values.fname)[checkfld.values.varname]
        for t in 1:size(testfld)[end]
            if length(size(dsvar)) == 3
                varcheck["values"] = varcheck["values"] && testfld[:,:,t] == dsvar[:,:,t]
            elseif length(size(dsvar)) == 4
                varcheck["values"] = varcheck["values"] && testfld[:,:,:,t] == dsvar[:,:,:,t]
            elseif length(size(dsvar)) == 2
                varcheck["values"] = varcheck["values"] && testfld[:,t] == dsvar[:,t]
            end
        end
    elseif isa(checkfld.values,TileData)
        itile = Int(ds.attrib["itile"])
        varcheck["values"] = true
        if isa(checkfld.values.vals,BinData)
            for t in 1:length(checkfld.values.vals.fnames)
                data = read(readbin(checkfld.values.vals.fnames[t],
                                    checkfld.values.vals.precision,checkfld.values.vals.iosize,
                                    checkfld.values.vals.fldidx),checkfld.values.tileinfo["XC"])
                testdata = NCTiles.gettile(data,checkfld.values.tileinfo,checkfld.values.tilesize,itile)
                if length(size(dsvar)) == 3
                    varcheck["values"] = varcheck["values"] && testdata == dsvar[:,:,t]
                elseif length(size(dsvar)) == 4
                    varcheck["values"] = varcheck["values"] && testdata == dsvar[:,:,:,t]
                elseif length(size(dsvar)) == 2
                    varcheck["values"] = varcheck["values"] && testdata == dsvar[:,t]
                end
            end
        else
            testdata = NCTiles.gettile(checkfld.values.vals,checkfld.values.tileinfo,checkfld.values.tilesize,itile)
            varcheck["values"] = varcheck["values"] && ((testdata .== dsvar) + isnan.(dsvar)) == ones(Int,size(dsvar)...)
        end
    end
    pass = pass && varcheck["values"]

    if ~pass
        missmatchstring = ""
        for k in keys(varcheck)
            if k == "dims"
                for d in keys(varcheck[k])
                    for a in keys(varcheck[k][d])
                        if ~varcheck[k][d][a]
                            missmatchstring = join([missmatchstring,join([k,d,a],'-')],", ")
                        end
                    end
                end
            else
                if ~varcheck[k]
                    missmatchstring = join([missmatchstring,k],", ")
                end
            end
        end
        println("The following file elements do not match the field elements:\n"*missmatchstring[3:end])
    end
    return pass
end

function maketestdata()
    lon=-180:20:180; lat=-90:20:90;
    depth = 5:10:100
    n1,n2,n3 = (length(lon),length(lat),length(depth))

    timeinterval = 1
    nsteps = 5
    time_steps = timeinterval/2:timeinterval:timeinterval*nsteps
    time_units = "days from 1992-01-01"

    data2d = [ones(Float32,n1,n2) for t in 1:nsteps]
    data3d = [ones(Float32,n1,n2,n3) for t in 1:nsteps]
    units = "test_units"
    longname2d = "Two Dimensional Test Data"
    longname3d = "Three Dimensional Test Data"

    dims = [NCvar("lon_c","degrees_east",n1,lon,Dict("long_name" => "longitude"),NCDatasets),
            NCvar("lat_c","degrees_north",n2,lat,Dict("long_name" => "latitude"),NCDatasets),
            NCvar("dep_c","m",n3,depth,Dict("long_name" => "depth","standard_name" => "depth","positive" => "down"),NCDatasets),
            NCvar("time","days from 1992-01-01",Inf,collect(time_steps),Dict(("long_name" => "time","standard_name" => "time")),NCDatasets)
            ]

    #fnames2d = ["testdata/data2d."*lpad(t,2,"0")*".data" for t in 1:length(data2d)]
    #fnames3d = ["testdata/data3d."*lpad(t,2,"0")*".data" for t in 1:length(data3d)]

    fnames2d = [mktemp()[1] for t in 1:length(data2d)]
    fnames3d = [mktemp()[1] for t in 1:length(data3d)]
    fnamestile2d = [mktemp()[1] for t in 1:length(data2d)]

    for t in 1:length(data2d)
        write(fnames2d[t],hton.(data2d[t]))
    end

    for t in 1:length(data3d)
        write(fnames3d[t],hton.(data3d[t]))
    end

    tilesize = (16,16)
    grid = GridSpec("CS32",getgrid()*"/")
    for f in fnamestile2d
        cp(joinpath(grid.path,"Depth.data"),f,force=true)
    end
    dep = -1 .* collect([-5 -15 -25 -35 -45 -55 -65 -75.0049972534180 -85.0250015258789 -95.0950012207031 -105.309997558594 -115.870002746582 -127.150001525879 -139.740005493164 -154.470001220703 -172.399993896484 -194.735000610352 -222.710006713867 -257.470001220703 -299.929992675781 -350.679992675781 -409.929992675781 -477.470001220703 -552.710021972656 -634.734985351563 -722.400024414063 -814.469970703125 -909.739990234375 -1007.15502929688 -1105.90502929688 -1205.53503417969 -1306.20495605469 -1409.15002441406  -1517.09497070313 -1634.17504882813 -1765.13500976563 -1914.15002441406 -2084.03491210938 -2276.22509765625 -2491.25000000000 -2729.25000000000 -2990.25000000000 -3274.25000000000 -3581.25000000000 -3911.25000000000 -4264.25000000000 -4640.25000000000 -5039.25000000000 -5461.25000000000 -5906.25000000000]')
    tile_ex = Dict(["tilesize" => tilesize,
                    "grid" => grid,
                    "dims" =>   [NCvar("i_c","1",tilesize[1],1:tilesize[1],Dict("long_name" => "Cartesian coordinate 1"),NCDatasets),
                                NCvar("j_c","1",tilesize[2],1:tilesize[2],Dict("long_name" => "Cartesian coordinate 2"),NCDatasets),
                                NCvar("dep_c","m",size(dep),dep,Dict("long_name" => "depth","standard_name" => "depth","positive" => "down"),NCDatasets),
                                NCvar("tim","days since 1992-1-1 0:0:0",Inf,time_steps,Dict(("long_name" => "time","standard_name" => "time")),NCDatasets)],
                    "fnamestile2d" => fnamestile2d
                    ])

    return Dict(["data2d" => data2d, "data3d" => data3d,
                 "units" => units, "dims" => dims,
                 "longname2d" => longname2d,"longname3d" => longname3d,
                 "fnames2d" => fnames2d, "fnames3d" => fnames3d, "tile_ex" => tile_ex])

end

function getgrid()

    testdir = abspath(joinpath(dirname(pathof(NCTiles)),"..","test"))
    griddir = joinpath(testdir,"GRID_CS32")
    if ~ispath(griddir)
        try
            file = joinpath(testdir,"GRID_CS32.tar.gz")
            download("https://github.com/gaelforget/GRID_CS32/archive/master.tar.gz",file)
            run(`tar -xzf $file -C $testdir`)
            olddir = joinpath(testdir,"GRID_CS32-master")
            run(`mv $olddir $griddir`)
            run(`rm $file`)
        catch e
            println("Could not download and extract grid. To fix, get the grid from https://github.com/gaelforget/GRID_CS32/ and put it in the "*testdir*" directory.")
        end
    end

    return testdir
end


text:
using DataDeps
using Statistics
using Oceananigans
using Oceananigans.Units
using ParameterEstimocean.Transformations: Transformation

fields_by_case = Dict(
   "free_convection" => (:b, :e),
   "weak_wind_strong_cooling" => (:b, :u, :v, :e),
   "strong_wind_weak_cooling" => (:b, :u, :v, :e),
   "strong_wind" => (:b, :u, :v, :e),
   "strong_wind_no_rotation" => (:b, :u, :e)
)

transformation = (b = Transformation(normalization=ZScore()),
                  u = Transformation(normalization=ZScore()),
                  v = Transformation(normalization=ZScore()),
                  e = Transformation(normalization=RescaledZScore(1e-1)))

function SyntheticObservationsBatch(path_fn, times, Nz; transformation=transformation, datadep = true, architecture = CPU(), field_names = (:b, :u, :v, :e), fields_by_case=fields_by_case)

   observations = Vector{SyntheticObservations}()

   for (case, forward_map_names) in zip(keys(fields_by_case), values(fields_by_case))

      data_path = datadep ? (@datadep_str path_fn(case)) : path_fn(case)
      observation = SyntheticObservations(data_path; transformation, times, field_names, forward_map_names, architecture, regrid=(1, 1, Nz))

      push!(observations, observation)
   end

   return observations
end

# Nz = 256
two_day_suite_path_1m(case) = "two_day_suite_1m/$(case)_instantaneous_statistics.jld2"
four_day_suite_path_1m(case) = "four_day_suite_1m/$(case)_instantaneous_statistics.jld2"
six_day_suite_path_1m(case) = "six_day_suite_1m/$(case)_instantaneous_statistics.jld2"

# Nz = 128
two_day_suite_path_2m(case) = "two_day_suite_2m/$(case)_instantaneous_statistics.jld2"
four_day_suite_path_2m(case) = "four_day_suite_2m/$(case)_instantaneous_statistics.jld2"
six_day_suite_path_2m(case) = "six_day_suite_2m/$(case)_instantaneous_statistics.jld2"

# Nz = 64
two_day_suite_path_4m(case) = "two_day_suite_4m/$(case)_instantaneous_statistics.jld2"
four_day_suite_path_4m(case) = "four_day_suite_4m/$(case)_instantaneous_statistics.jld2"
six_day_suite_path_4m(case) = "six_day_suite_4m/$(case)_instantaneous_statistics.jld2"

TwoDaySuite(; transformation=transformation, times=[2hours, 12hours, 1days, 36hours, 2days],
            Nz=64, architecture=CPU(), field_names = (:b, :u, :v, :e)) = SyntheticObservationsBatch(two_day_suite_path_2m, times, Nz; architecture, transformation, field_names)

FourDaySuite(; transformation=transformation, times=[2hours, 1days, 2days, 3days, 4days],
            Nz=64, architecture=CPU(), field_names = (:b, :u, :v, :e)) = SyntheticObservationsBatch(four_day_suite_path_2m, times, Nz; architecture, transformation, field_names)

SixDaySuite(; transformation=transformation, times=[2hours, 1.5days, 3days, 4.5days, 6days],
            Nz=64, architecture=CPU(), field_names = (:b, :u, :v, :e)) = SyntheticObservationsBatch(six_day_suite_path_2m, times, Nz; architecture, transformation, field_names)

function lesbrary_ensemble_simulation(observations;
                                             Nensemble = 30,
                                             architecture = CPU(),
                                             closure = ConvectiveAdjustmentVerticalDiffusivity(),
                                             Δt = 10.0
                                    )

    simulation = ensemble_column_model_simulation(observations;
                                                  Nensemble,
                                                  architecture,
                                                  tracers = (:b, :e),
                                                  closure)

    simulation.Δt = Δt

    Qᵘ = simulation.model.velocities.u.boundary_conditions.top.condition
    Qᵇ = simulation.model.tracers.b.boundary_conditions.top.condition
    N² = simulation.model.tracers.b.boundary_conditions.bottom.condition

    for (case, obs) in enumerate(observations)
      try
         f = obs.metadata.parameters.coriolis_parameter
         view(Qᵘ, :, case) .= obs.metadata.parameters.momentum_flux
         view(Qᵇ, :, case) .= obs.metadata.parameters.buoyancy_flux
         view(N², :, case) .= obs.metadata.parameters.N²_deep
         view(simulation.model.coriolis, :, case) .= Ref(FPlane(f=f))
      catch
         f = obs.metadata.coriolis.f
         view(Qᵘ, :, case) .= obs.metadata.parameters.Qᵘ
         view(Qᵇ, :, case) .= obs.metadata.parameters.Qᵇ
         view(N², :, case) .= obs.metadata.parameters.N²
         view(simulation.model.coriolis, :, case) .= Ref(FPlane(f=f))
      end
    end

    return simulation
end

function estimate_η_covariance(output_map, observations)

   @assert length(observations) > 2 "A two-sample covariance matrix has rank one and is therefore singular.
                                      Please increase the number of `observations` to at least 3."
   obs_maps = hcat([observation_map(output_map, obs) for obs in observations]...)
   return cov(transpose(obs_maps), corrected=false)
end


text:
using Coverage, Pkg
# process '*.cov' files
Pkg.test(; coverage=true)
coverage = process_folder() # defaults to src/; alternatively, supply the folder name as argument
LCOV.writefile("lcov.info", coverage)
clean_folder("src")


text:
using VUMPS
using VUMPS:qrpos,lqpos,leftorth,rightorth,leftenv,FLmap,rightenv,FRmap,ACenv,ACmap,Cenv,Cmap,LRtoC,ALCtoAC,ACCtoALAR,error,bigleftenv,BgFLmap,bigrightenv,BgFRmap, obs_FL, obs_FR, norm_FL, norm_FR, norm_FLmap, norm_FRmap
using CUDA
using LinearAlgebra
using Random
using Test
using OMEinsum
CUDA.allowscalar(false)

@testset "qr with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64]
    Random.seed!(100)
    A = atype(rand(dtype, 4,4))
    Q, R = qrpos(A)
    @test Array(Q*R) ≈ Array(A)
    @test all(real.(diag(R)) .> 0)
    @test all(imag.(diag(R)) .≈ 0)
end

@testset "lq with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64]
    Random.seed!(100)
    A = atype(rand(dtype, 4,4))
    L, Q = lqpos(A)
    @test Array(L*Q) ≈ Array(A)
    @test all(real.(diag(L)) .> 0)
    @test all(imag.(diag(L)) .≈ 0)
end

@testset "leftorth and rightorth with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64], Ni = [2], Nj = [2]
    Random.seed!(100)
    D, d = 3, 2
    A = Array{atype{dtype,3},2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        A[i,j] = atype(rand(dtype, D, d, D))
    end
    AL, L, λ = leftorth(A)
    R, AR, λ = rightorth(A)

    for j = 1:Nj,i = 1:Ni
        M = ein"cda,cdb -> ab"(AL[i,j],conj(AL[i,j]))
        @test (Array(M) ≈ I(D))

        LA = reshape(L[i,j] * reshape(A[i,j], D, d*D), d*D, D)
        ALL = reshape(AL[i,j], d*D, D) * L[i,j] * λ[i,j]
        @test (Array(ALL) ≈ Array(LA))

        M = ein"acd,bcd -> ab"(AR[i,j],conj(AR[i,j]))
        @test (Array(M) ≈ I(D))

        AxR = reshape(reshape(A[i,j], d*D, D)*R[i,j], D, d*D)
        RAR = R[i,j] * reshape(AR[i,j], D, d*D) * λ[i,j]
        @test (Array(RAR) ≈ Array(AxR))
    end
end

@testset "leftenv and rightenv with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64], Ni = [2], Nj = [2]
    Random.seed!(100)
    D, d = 3, 2
    A = Array{atype{dtype,3},2}(undef, Ni, Nj)
    M = Array{atype{dtype,4},2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        A[i,j] = atype(rand(dtype, D, d, D))
        M[i,j] = atype(rand(dtype, d, d, d, d))
    end

    AL, = leftorth(A)
    λL,FL = leftenv(AL, AL, M)
    _, AR, = rightorth(A)
    λR,FR = rightenv(AR, AR, M)

    for j = 1:Nj, i = 1:Ni
        ir = Ni + 1 - i
        @test λL[i,j] * FL[i,j] ≈ FLmap(AL[i,:], conj(AL[ir,:]), M[i,:], FL[i,j], j)
        @test λR[i,j] * FR[i,j] ≈ FRmap(AR[i,:], conj(AR[ir,:]), M[i,:], FR[i,j], j)
    end
end

@testset "observable leftenv and rightenv with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64], Ni = [2], Nj = [2]
    Random.seed!(100)
    D, d = 3, 2
    A = Array{atype{dtype,3},2}(undef, Ni, Nj)
    M = Array{atype{dtype,4},2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        A[i,j] = atype(rand(dtype, D, d, D))
        M[i,j] = atype(rand(dtype, d, d, d, d))
    end

    AL, = leftorth(A)
    λL,FL = obs_FL(AL, AL, M)
    _, AR, = rightorth(A)
    λR,FR = obs_FR(AR, AR, M)

    for j = 1:Nj, i = 1:Ni
        ir = Ni + 1 - i
        @test λL[i,j] * FL[i,j] ≈ FLmap(AL[i,:], AL[ir,:], M[i,:], FL[i,j], j)
        @test λR[i,j] * FR[i,j] ≈ FRmap(AR[i,:], AR[ir,:], M[i,:], FR[i,j], j)
    end
end

@testset "ACenv and Cenv with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64], Ni = [2], Nj = [2]
    Random.seed!(100)
    D, d = 3, 2
    A = Array{atype{dtype,3},2}(undef, Ni, Nj)
    M = Array{atype{dtype,4},2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        A[i,j] = atype(rand(dtype, D, d, D))
        M[i,j] = atype(rand(dtype, d, d, d, d))
    end

    AL, L = leftorth(A)
    λL,FL = leftenv(AL, AL, M)
    R, AR, = rightorth(A)
    λR,FR = rightenv(AR, AR, M)

    C = LRtoC(L, R)
    AC = ALCtoAC(AL, C)

    λAC, AC = ACenv(AC, FL, M, FR)
    λC, C = Cenv(C, FL, FR)
    for j = 1:Nj, i = 1:Ni
        jr = j + 1 - Nj * (j==Nj)
        @test λAC[i,j] * AC[i,j] ≈ ACmap(AC[i,j], FL[:,j], FR[:,j], M[:,j], i)
        @test λC[i,j] * C[i,j] ≈ Cmap(C[i,j], FL[:,jr], FR[:,j], i)
    end
end

@testset "bcvumps unit test with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64], Ni = [2], Nj = [2]
    Random.seed!(100)
    D, d = 3, 2
    A = Array{atype{dtype,3},2}(undef, Ni, Nj)
    M = Array{atype{dtype,4},2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        A[i,j] = atype(rand(dtype, D, d, D))
        M[i,j] = atype(rand(dtype, d, d, d, d))
    end

    AL, L = leftorth(A)
    λL,FL = leftenv(AL, AL, M)
    R, AR, = rightorth(A)
    λR,FR = rightenv(AR, AR, M)

    C = LRtoC(L,R)
    AC = ALCtoAC(AL,C)

    λAC, AC = ACenv(AC, FL, M, FR)
    λC, C = Cenv(C, FL, FR)
    AL, AR = ACCtoALAR(AC, C)
    err = error(AL,C,AR,FL,M,FR)
    @test err !== nothing
end

@testset "bigleftenv and bigrightenv with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64], Ni = [2], Nj = [2]
    Random.seed!(100)
    D, d = 3, 2
    A = Array{atype{dtype,3},2}(undef, Ni, Nj)
    M = Array{atype{dtype,4},2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        A[i,j] = atype(rand(dtype, D, d, D))
        M[i,j] = atype(rand(dtype, d, d, d, d))
    end

    ALu, = leftorth(A)
    ALd, = leftorth(A)
    λL,BgFL = bigleftenv(ALu, ALd, M)
    _, ARu, = rightorth(A)
    _, ARd, = rightorth(A)
    λR,BgFR = bigrightenv(ARu, ARd, M)

    for j = 1:Nj, i = 1:Ni
        ir = i + 1 - Ni * (i==Ni)
        irr = i + 2 - Ni * (i + 2 > Ni)
        @test λL[i,j] * BgFL[i,j] ≈ BgFLmap(ALu[i,:], ALd[irr,:], M[i,:], M[ir,:], BgFL[i,j], j)
        @test λR[i,j] * BgFR[i,j] ≈ BgFRmap(ARu[i,:], ARd[irr,:], M[i,:], M[ir,:], BgFR[i,j], j)
    end
end

@testset "norm leftenv and rightenv with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64], Ni = [2], Nj = [2]
    Random.seed!(100)
    D, d = 3, 2
    A = Array{atype{dtype,3},2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        A[i,j] = atype(rand(dtype, D, d, D))
    end

    ALu, = leftorth(A)
    ALd, = leftorth(A)
    λL, FL_norm = norm_FL(ALu, ALd)
    _, ARu, = rightorth(A)
    _, ARd, = rightorth(A)
    λR, FR_norm = norm_FR(ARu, ARd)

    for j = 1:Nj, i = 1:Ni
        @test λL[i,j] * FL_norm[i,j] ≈ norm_FLmap(ALu[i,:], ALu[i,:], FL_norm[i,j], j)
        @test λR[i,j] * FR_norm[i,j] ≈ norm_FRmap(ARu[i,:], ARd[i,:], FR_norm[i,j], j)
    end
end


text:
using quantumoptics

N = 500000
xmin = -62.5
xmax = 70.1

basis_position = quantumoptics.particle.PositionBasis(xmin, xmax, N)
basis_momentum = quantumoptics.particle.MomentumBasis(basis_position)

x0 = 5.1
p0 = -3.2
sigma = 1.
sigma_x = sigma/sqrt(2)
sigma_p = 1./(sigma*sqrt(2))

psix0 = quantumoptics.particle.gaussianstate(basis_position, x0, p0, sigma)
psip0 = quantumoptics.particle.gaussianstate(basis_momentum, x0, p0, sigma)

function transformation(b1::quantumoptics.particle.MomentumBasis, b2::quantumoptics.particle.PositionBasis, psi::Ket)
    Lp = (b1.pmax - b1.pmin)
    dx = quantumoptics.particle.spacing(b2)
    if b1.N != b2.N || abs(2*pi/dx - Lp)/Lp > 1e-12
        throw(IncompatibleBases())
    end
    N = b1.N
    psi_shifted = exp(1im*b2.xmin*(quantumoptics.particle.samplepoints(b1)-b1.pmin)).*psi.data
    psi_fft = exp(1im*b1.pmin*quantumoptics.particle.samplepoints(b2)).*ifft(psi_shifted)*sqrt(N)
    return Ket(b2, psi_fft)
end

function transformation(b1::quantumoptics.particle.PositionBasis, b2::quantumoptics.particle.MomentumBasis, psi::Ket)
    Lx = (b1.xmax - b1.xmin)
    dp = quantumoptics.particle.spacing(b2)
    if b1.N != b2.N || abs(2*pi/dp - Lx)/Lx > 1e-12
        throw(IncompatibleBases())
    end
    N = b1.N
    psi_shifted = exp(-1im*b2.pmin*(quantumoptics.particle.samplepoints(b1)-b1.xmin)).*psi.data
    psi_fft = exp(-1im*b1.xmin*quantumoptics.particle.samplepoints(b2)).*fft(psi_shifted)/sqrt(N)
    return Ket(b2, psi_fft)
end

@time psix0_fft = transformation(basis_position, basis_momentum, psix0)
@time psix0_fft = transformation(basis_position, basis_momentum, psix0)

function transformation!(b1::quantumoptics.particle.PositionBasis, b2::quantumoptics.particle.MomentumBasis, psi::Ket)
    N::Int = b1.N
    dp = quantumoptics.particle.spacing(b2)
    dx = quantumoptics.particle.spacing(b1)
    @inbounds for i=1:N
        psi.data[i] .*= exp(-1im*b2.pmin*(i-1)*dx)
    end
    fft!(psi.data)
    @inbounds for i=1:N
        psi.data[i] .*= exp(-1im*b1.xmin*(b2.pmin + (i-1)*dp))/sqrt(N)
    end
    return psi
end

psi_ = deepcopy(psix0)
@time transformation!(basis_position, basis_momentum, psi_)
psi_ = deepcopy(psix0)
@time transformation!(basis_position, basis_momentum, psi_)

function transformation!(b1::quantumoptics.particle.PositionBasis, b2::quantumoptics.particle.MomentumBasis, psi::Ket, fft_plan!, mul_before, mul_after)
    N::Int = b1.N
    @inbounds for i=1:N
        psi.data[i] *= mul_before[i]
    end
    fft_plan!(psi.data)
    @inbounds for i=1:N
        psi.data[i] *= mul_after[i]
    end
    return psi
end

mul_before = exp(-1im*basis_momentum.pmin*(quantumoptics.particle.samplepoints(basis_position)-basis_position.xmin))
mul_after = exp(-1im*basis_position.xmin*quantumoptics.particle.samplepoints(basis_momentum))/sqrt(N)

A = plan_fft!(psix0.data)
psi_ = deepcopy(psix0)
@time transformation!(basis_position, basis_momentum, psi_, A, mul_before, mul_after)
A = plan_fft!(psix0.data)
psi_ = deepcopy(psix0)
@time transformation!(basis_position, basis_momentum, psi_, A, mul_before, mul_after)

T = quantumoptics.particle.FFTOperator(basis_momentum, basis_position)
@time T*psix0
@time T*psix0

psi_ = deepcopy(psix0)
@time quantumoptics.operators.gemv!(Complex(1.), T, psix0, Complex(0.), psi_)
psi_ = deepcopy(psix0)
@time quantumoptics.operators.gemv!(Complex(1.), T, psix0, Complex(0.), psi_)


text:
"""
    match_num_point!(poly_1::Vector{Point}, poly_2::Vector{Point})

This is a helper function for [`morph`](@ref).
Given two polygons `poly_1` and `poly_2` points are added to the polygon with less points
until both polygons have the same number of points.
The polygon with less points gets mutated during this process.

# Arguments
- `poly_1::Vector{Point}`: The points which define the first polygon
- `poly_2::Vector{Point}`: The points which define the second polygon
"""
function match_num_point!(poly_1::Vector{Point}, poly_2::Vector{Point})
    l1 = length(poly_1)
    l2 = length(poly_2)
    # if both have the same number of points => we are done
    l1 == l2 && return

    # poly_1 should have less points than poly_2 so we flip if this is not the case
    flipped = false
    if l1 > l2
        poly_1, poly_2 = poly_2, poly_1
        l1, l2 = l2, l1
        flipped = true
    end

    # the difference of the length of points
    diff = l2 - l1

    points_per_edge = div(diff, l1)
    # how many extra points do we need
    points_per_edge_extra = rem(diff, l1)
    # => will add them to the first `points_per_edge_extra` edges

    # index is the index where the next point is added
    index = 2
    poly_1_orig = copy(poly_1)
    for i in 1:l1
        # p1 is the current point in the original polygon
        p1 = poly_1_orig[i]
        # p2 is the next point (which is the first of the polygon in the last iteration)
        if i + 1 > l1
            p2 = poly_1_orig[1]
        else
            p2 = poly_1_orig[i + 1]
        end
        # if we need 5 points and have only 4 edges we add 2 points for the first edge
        rem = 0
        if i <= points_per_edge_extra
            rem = 1
        end
        for j in 1:(points_per_edge + rem)
            # create the interpolated point between p1 and p2
            t = j / (points_per_edge + rem + 1)
            new_point = p1 + t * (p2 - p1)
            insert!(poly_1, index, new_point)
            index += 1
        end
        index += 1
    end

    if flipped
        poly_1, poly_2 = poly_2, poly_1
    end
    @assert length(poly_1) == length(poly_2)
end

"""
    morph(from_func::Function, to_func::Function; action=:stroke)

A closure for the [`_morph`](@ref) function.
This makes it easier to write the function inside an `Action`.

Currently morphing is quite simple and only works for basic shapes.
It especially does not work with functions which produce more than one polygon
or which produce filled polygons.
Blending between fills of polygons is definitely coming at a later stage.

**Important:** The functions itself should not draw the polygon
i.e. use `circle(Point(100,100), 50)` instead of `circle(Point(100,100), 50, :stroke)`

# Arguments
- `from_func::Function`: The function that creates the path for the first polygon.
- `to_func::Function`: Same as `from_func` but it defines the "result" polygon,
                       which will be displayed at the end of the Action

# Keywords
- `action::Symbol` defines whether the object has a fill or just a stroke. Defaults to stroke.

# Example

This creates a star that morphs into a circle and back.

```
using Javis

astar(args...) = star(O, 50)
acirc(args...) = circle(Point(100,100), 50)

video = Video(500, 500)
javis(video, [
    Action(1:100, ground),
    Action(1:50, morph(astar, acirc)),
    Action(51:100, morph(acirc, astar))
], creategif=true, tempdirectory="images",
    pathname="star2circle.gif", deletetemp=true)
```
"""
function morph(from_func::Function, to_func::Function; action = :stroke)
    return (video, scene_action, frame) ->
        _morph(video, scene_action, frame, from_func, to_func; draw_action = action)
end

"""
    save_morph_polygons!(action::Action, from_func::Function, to_func::Function)

Converts the paths created by the functions to polygons and calls [`match_num_point!`](@ref)
such that both polygons have the same number of points.
This is done once inside [`_morph`](@ref).
Saves the two polygons inside `action.opts[:from_poly]` and `action.opts[:to_poly]`.

**Assumption:** Both functions create only a single polygon each.
"""
function save_morph_polygons!(action::Action, from_func::Function, to_func::Function)
    newpath()
    from_func()
    closepath()
    from_poly = pathtopoly()[1]

    newpath()
    to_func()
    closepath()
    to_poly = pathtopoly()[1]

    match_num_point!(from_poly, to_poly)

    # find the smallest morphing distance to match the points in a more natural way
    # smallest_i holds the best starting point of from_path
    smallest_i = 1
    smallest_distance = typemax(Float64)

    for i in 1:length(from_poly)
        overall_distance = 0.0
        for j in 1:length(from_poly)
            p1 = from_poly[(j + i - 1) % length(from_poly) + 1]
            p2 = to_poly[j]
            overall_distance += distance(p1, p2)
        end
        if overall_distance < smallest_distance
            smallest_distance = overall_distance
            smallest_i = i
        end
    end
    new_from_poly = copy(from_poly)
    for i in 1:length(from_poly)
        new_from_poly[i] = from_poly[(i + smallest_i - 1) % length(from_poly) + 1]
    end

    action.opts[:from_poly] = new_from_poly
    action.opts[:to_poly] = to_poly
    action.opts[:points] = Vector{Point}(undef, length(new_from_poly))
end

"""
    _morph(video::Video, action::Action, frame, from_func::Function, to_func::Function; draw_action=:stroke)

Internal version of [`morph`](@ref) but described there.
"""
function _morph(
    video::Video,
    action::Action,
    frame,
    from_func::Function,
    to_func::Function;
    draw_action = :stroke,
)
    # computation of the polygons and the best way to morph in the first frame
    if frame == first(get_frames(action))
        save_morph_polygons!(action, from_func, to_func)
    end

    # obtain the computed polygons. These polygons have the same number of points.
    from_poly = action.opts[:from_poly]
    to_poly = action.opts[:to_poly]
    points = action.opts[:points]

    # compute the interpolation variable `t` for the current frame
    t = get_interpolation(action, frame)

    for (i, p1, p2) in zip(1:length(from_poly), from_poly, to_poly)
        new_point = p1 + t * (p2 - p1)
        points[i] = new_point
    end

    poly(points, draw_action; close = true)
end


text:
module BenchAlpha

using BenchmarkTools
using Ripserer

include(joinpath(@__DIR__, "utils.jl"))

benchmarks = [(file="torus10k.pts", dim_max=2), (file="klein500.pts", dim_max=3)]

suite = BenchmarkGroup()

for b in benchmarks
    name = b.file
    data = load_data(joinpath(@__DIR__, "data", b.file))
    suite[name] = @benchmarkable ripserer(Alpha($data); dim_max=$(b.dim_max)) samples = 2
end

end
BenchAlpha.suite


text:
import Splittings:UniformMesh
import Splittings:@Strang

"""

   exact(tf, mesh1, mesh2)

   Julia function to compute exact solution "

```math
\frac{d f}{dt} +  (y \frac{df}{delta1} - x \frac{df}{delta2}) = 0
```

"""
function exact(tf, mesh1::UniformMesh, mesh2::UniformMesh)

    f = zeros(Float64,(mesh1.length,mesh2.length))
    for (i, x) in enumerate(mesh1.points), (j, y) in enumerate(mesh2.points)
        xn = cos(tf) * x - sin(tf) * y
        yn = sin(tf) * x + cos(tf) * y
        f[i,j] = exp(-(xn-1)*(xn-1)/0.1)*exp(-(yn-1)*(yn-1)/0.1)
    end

    f

end

" Function to compute error "
function error1(f, f_exact)
    maximum(abs.(f .- f_exact))
end

import FFTW: FFTWPlan, plan_fft, fft!, ifft!

struct Fourier

    px :: FFTWPlan
    py :: FFTWPlan
    kx :: Vector{Float64}
    ky :: Vector{Float64}

    function Fourier( mesh1 :: UniformMesh,
		      mesh2 :: UniformMesh,
		      f  :: Array{Complex{Float64},2},
                      fᵗ :: Array{Complex{Float64},2})

        px = plan_fft(f,  1)
        py = plan_fft(fᵗ, 1)

        n1   = mesh1.length
        x1min = mesh1.start
        x1max = mesh1.stop
        kx = 2π/(x1max-x1min)*[0:n1÷2-1;n1÷2-n1:-1]

        n2   = mesh2.length
        ymin = mesh2.start
        ymax = mesh2.stop
        ky = 2π/(ymax-ymin)*[0:n2÷2-1;n2÷2-n2:-1]

	new( px, py, kx, ky)

    end
end

function advection_y!( f ::Array{Complex{Float64},2},
		       fᵗ::Array{Complex{Float64},2},
		       f̂ᵗ::Array{Complex{Float64},2},
		       mesh1::UniformMesh,
		       dt::Float64,
		       fourier::Fourier)

    exky = exp.( 1im * dt * fourier.ky .* transpose(mesh1.points))

    transpose!(fᵗ,f)
    mul!(f̂ᵗ,  fourier.py, fᵗ)
    f̂ᵗ .= f̂ᵗ .* exky
    ldiv!(fᵗ, fourier.py, f̂ᵗ)
    transpose!(f,fᵗ)

end

function advection_x!( f ::Array{Complex{Float64},2},
		       f̂ ::Array{Complex{Float64},2},
		       mesh2::UniformMesh,
		       dt::Float64,
		       fourier::Fourier)

    ekxy = exp.(-1im * dt * fourier.kx .* transpose(mesh2.points))

    mul!(f̂,  fourier.px, f)
    f̂ .= f̂ .* ekxy
    ldiv!(f, fourier.px, f̂)

end

function rotation_2d_fft(tf, nt, mesh1::UniformMesh, mesh2::UniformMesh)

    dt = tf/nt

    n1 = mesh1.length
    x1min, x1max = mesh1.start, mesh1.stop
    delta1 = mesh1.step

    n2 = mesh2.length
    ymin, ymax = mesh2.start, mesh2.stop
    delta2 = mesh2.step

    f  = zeros(Complex{Float64},(n1,n2))
    f̂  = similar(f)
    fᵗ = zeros(Complex{Float64},(n2,n1))
    f̂ᵗ = similar(fᵗ)

    fourier = Fourier( mesh1, mesh2, f, fᵗ)

    f .= exact(0.0, mesh1, mesh2)

    for n=1:nt

        @Strang( advection_y!( f, fᵗ, f̂ᵗ, mesh1, tan(dt), fourier),
		 advection_x!( f, f̂, mesh2, sin(dt), fourier))

    end
    real(f)
end

tf, nt = 200π, 1000

mesh1 = UniformMesh(-π, π, 128; endpoint=false)
mesh2 = UniformMesh(-π, π, 256; endpoint=false)

fc = rotation_2d_fft(tf, nt, mesh1, mesh2)
fe = exact(tf, mesh1, mesh2)

@testset "Rotation test with Fourier advections " begin

println(error1(fc, fe))
@test rotation_2d_fft(tf, nt, mesh1, mesh2) ≈ fe atol = 1e-10

end


text:
## Index convention:
## [i, j, k], where k = 1, 2 is the component:
## 1 - gas,
## 2 - liquid;
## [i, j] -> x, y

## -------------------------- BC -------------------------- ##
## (a f(x, y) + b ∂F / ∂n)|_{∂Ω} = g
##
## We keep a and b in an unmutable struct because they are
## set as a parameter in the beginning of the simulation,
## and they don't change.
## Therefore, they are kept with the other simulation parameters
##
## g, on the other hand, will get updated on every step,
## so it seems logical to keep together with all the other
## system parameters that change over time.
## We will denote it as f_RHS, since it is then easy to infer
## that we are refering the right hand side of the BC equation.
##
## The struct BoundaryCondtions holds two 3D arrays, a and b.
##
## Each of this holds for each component 4 arrays,
## one for each border: left, right, top, bottom.
##
## This is stored in a 3D array, where the last index
## k corresponds to the component, and the second index
## denotes the border.
##
## h[:][1][k] -> left
##
## h[:][2][k] -> right
##
## h[:][3][k] -> top
##
## h[:][4][k] -> bottom
##
## Pressure is the same for both components, but for
## consistency we will also use the same structure, but
## use only the first index k = 1.

struct BoundaryCondition{T}
    a::Array{T, 3}
    b::Array{T, 3}
end

## Here we keep the parameters a and b for each variable.
struct BoundaryConditions{T}
    u::BoundaryCondition{T}
    v::BoundaryCondition{T}
    ρ::BoundaryCondition{T}
    s::BoundaryCondition{T}
    P::BoundaryCondition{T}
end

## f_RHS == g
mutable struct BoundaryRHS{T}
    u::Array{T, 3}
    v::Array{T, 3}
    ρ::Array{T, 3}
    s::Array{T, 3}
    P::Array{T, 3}
end

Base.copy(bc_RHS::BoundaryRHS) = BoundaryRHS(copy.((bc_RHS.u,
                                                    bc_RHS.v,
                                                    bc_RHS.ρ,
                                                    bc_RHS.s,
                                                    bc_RHS.P,))... )

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

mutable struct System{T<:AbstractFloat}
    u::Array{T, 3}
    v::Array{T, 3}
    ρ::Array{T, 3}
    F::Array{T, 3}
    s::Array{T, 3}
    P::Array{T, 2}
    bc_RHS::BoundaryRHS{T}
end

Base.copy(sys::System) = System(copy.((sys.u, sys.v, sys.ρ,
                                       sys.F, sys.s, sys.P,
                                       sys.bc_RHS))... )

struct Parameters{T<:AbstractFloat}
    nsteps::Integer
    nx::Integer
    ny::Integer
    L::Integer
    Δx::T
    Δy::T
    Δt::T
    φ::T
    K::T
    Pin::T
    Pout::T
    P₀::T
    ψ::T
    ψ₀::T
    M::Vector{T}
    μ::Vector{T}
    bc::BoundaryConditions{T}
    wall::UnitRange{Int}
    inlet::UnitRange{Int}
    relax_steps::Integer
end

struct TaitEoS{T}
    B::T
    P₀::T
    C::T
    ρ₀::T
end

struct IdealGasEoS{T}
    T::T
    M::T
end

max_stars_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_issues_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_forks_repo_name: sofiabelen/Two-Phase-Filtration


text:
using CoordinatedSupplyChains
using Test

@testset "CoordinatedSupplyChains.jl" begin
    # Write your tests here
    dir = joinpath(@__DIR__,"TestV01")

    A,N,P,D,S,T,L,Sets,Pars = CoordinatedSupplyChains.LoadSSCaseData(dir)
    @test A.ID[1] == "A01"

    Output, Stats = OptimizeSSCase(A,N,P,D,S,T,L,Sets,Pars; CaseDataDirectory=dir)
    @test Output.z == 1.450100720657891e6
    @test Stats.NumVars == 3490

    SSRecordMaker(A,N,P,D,S,T,L,Output,Stats; CaseDataDirectory=dir)
    @test isdir(joinpath(dir,"records"))==true
    @test isfile(joinpath(dir,"records/model_stats.txt"))==true

    SSNetworkPlot(A,N,P,D,S,L,Output.f; CaseDataDirectory=dir, TestMode=true)
    @test isdir(joinpath(dir,"plots"))==true
    #@test isfile(joinpath(dir,"plots/NetworkPlot.png"))==true

    # Delete everything added by the code once the tests are run
    rm(joinpath(dir,"_Model.txt"))
    rm(joinpath(dir,"_SolutionData.txt"))
    rm(joinpath(dir,"records"), recursive=true)
    rm(joinpath(dir,"plots"), recursive=true)

    # RunSSCase calls the four previous functions again; the optional keyword options are reversed here to improve coverage
    RunSSCase(dir; PrintModel=true, WriteReport=false, PrintOutput=false, TestMode=true)
    @test isfile(joinpath(dir,"records/model_stats.txt"))==true
    #@test isfile(joinpath(dir,"plots/NetworkPlot.png"))==true

    # Delete everything added by the code once the tests are run; again
    rm(joinpath(dir,"_Model.txt"))
    rm(joinpath(dir,"records"), recursive=true)
    rm(joinpath(dir,"plots"), recursive=true)
end