text:
# Module that exports a few functions that when called return DataCell datasets
module DataGenerator
	import j4pr;
	export iris, quakes, boston, normal2d, fish, fishr, fnoise;

	"""
		iris()

	DataCell version of the `iris` dataset. The `RDatasets` package is required.
	"""
	# Iris dataset
	function iris()
		irisdf = j4pr.rdataset("iris");
		return j4pr.datacell(Array(irisdf[:,1:4])', Array(irisdf[:,5]), name="Iris Dataset")
	end

	"""
		quakes(thresh::Float64=5.0)

	DataCell version of the `quakes` dataset. The `thresh` parameter is used to threshold
	a continuous variable, in this case the quake magnitude to obtain the labels
	(e.g. 0 for < `thresh`, 1 >= `thresh`). The `RDatasets` package is required.
	"""
	# Quakes dataset
	function quakes(thresh=5.0)
		quakedf = j4pr.rdataset("quakes")
		labels = Array(quakedf[:,4].>=thresh)
		return j4pr.datacell(Array(quakedf[:,[1,2,3,5]])', labels, name="Quakes Dataset (threshold >=$thresh)")
	end

	"""
		boston()

	DataCell version of the `boston` dataset. The `RDatasets` package is required.
	"""
	# Boston dataset
	function boston()
		bostondf = j4pr.rdataset("Boston")
		return j4pr.datacell(Array{typeof(1.0)}(bostondf[:,1:13])', Array(bostondf[:,14]), name="Boston Dataset")
	end

	"""
		normal2d(N::Int=100)

	Dataset containing two classes and two variables  with `N` samples (default `100`), normally distributed
	so that the class distributions look perpendicular and overlapping; the mean vectors for each class
	and variable vary a little each time a new dataset is generated.
	"""
	# Two gaussians
	function normal2d(N::Int=100)
		N0 = round(Int, N/2)
		N1 = N - N0
		X = [	randn(N0)+round.(rand([1,2,3])) randn(N0)+round.(rand([1,2,3])); 		# class 0
       			5*randn(N1)+round.(rand([4,5,8])) 0.1*randn(N1)+round.(rand([1,2,3])); 		# class 1
		]'
		y = [zeros(N0);ones(N1)]
		order = shuffle(1:N)
		return j4pr.datacell(X[:,order], y[order], name="Two 2-D Gaussians")
	end

	"""
		fish(N::Int=100, noise::Float64=0.1)

	The `Fish dataset` is a 2-D dataset portraying a fish which represents one class, with the background as other class; noise
	can be added. `N` is the size of the grid (default `100`, e.g. grid of 100×100 points) and `noise` (default `0.1`) is the
	fraction of samples with their labels flipped (e.g. noisy).

	Dataset based on the "Fish data" from `Kuncheva L. "Combining Pattern Classifiers", 2-nd Ed. ISBN 978-1-118-31523-1`
	"""
	# Fish
	function fish(N::Int=100, noise::Float64=0.01)

		data = zeros(2,N^2)

		# generate the 2-D grid
		lab1 = BitVector(N*N);
		lab2 = BitVector(N*N);
		for (i,y) in enumerate(1:N), (j,x) in enumerate(1:N)
			xn = x/N; yn = y/N;
			lab1[(i-1)*N+j] = xn^3 -2xn*yn + 1.6yn^2 < 0.4
			lab2[(i-1)*N+j] = -xn^3 + 2*sin(xn)*yn + yn < 0.7
			data[1,(i-1)*N+j] = x	# x dimension 1'st variable
			data[2,(i-1)*N+j] = y	# y dimension 2'nd variable
		end
		l = xor.(lab1, lab2)
		nnoisy = round(Int, N^2*noise)
		if nnoisy > 0.0
			nidx = randperm(N^2)[1:nnoisy] # get indexes of noisy samples
			l[nidx]=.!l[nidx]
		end

		labels = map(x->x ? "fish" : "background", l)
		return j4pr.datacell(data, labels[:], name="Fish Dataset, classification")
	end

	"""
		fishr(N::Int=100, noise::Float64=0.1)

	Version of the `fish` dataset for regression problems.

	Dataset based on the "Fish data" from `Kuncheva L. "Combining Pattern Classifiers", 2-nd Ed. ISBN 978-1-118-31523-1`
	"""
	function fishr(N::Int=100, noise::Float64=0.01)

		data = zeros(2,N^2)

		# generate the 2-D grid
		lab1 = falses(N*N);
		lab2 = falses(N*N);
		for (i,y) in enumerate(1:N), (j,x) in enumerate(1:N)
			xn = x/N; yn = y/N;
			lab1[(i-1)*N+j] = xn^3 -2xn*yn + 1.6yn^2 < 0.4
			lab2[(i-1)*N+j] = -xn^3 + 2*sin(xn)*yn + yn < 0.7
			data[1,(i-1)*N+j] = x	# x dimension 1'st variable
			data[2,(i-1)*N+j] = y	# y dimension 2'nd variable
		end
		l = xor.(lab1, lab2)
		nnoisy = round(Int, N^2*noise)
		labels = map(x->x ? 1.0 : 0.0, l)
		if nnoisy > 0.0
			nidx = randperm(N^2)[1:nnoisy] # get indexes of noisy samples
			labels[nidx]+=rand(nnoisy)
		end

		return j4pr.datacell(data, labels[:], name="Fish Dataset, regression")
	end

	"""
		fnoise(points, f, fnoise)

	Generates a dataset containing the `points` as data and the function `f`
	applied to these points as values. Noise is added through `fnoise`. The
	full formula for data generation is `y = f.(points) + fnoise.(f.(points))`

	# Arguments
	  * `points::AbstractVector` is the vector of points
	  * `f` is a function applied to each point (default `sin`)
	  * `fnoise` is a function applied to each value of the function (default `x->x+rand()`)

	# Examples
	```
	julia> data = DataGenerator.fnoise(0:0.01:pi, x->sin(x), x->x*rand())
	Function #71, with noise, 315 obs, 1 vars, 1 target(s)/obs, 315 distinct values

	julia> data |> x->UnicodePlots.lineplot(collect(x.x),x.y)
	     ┌────────────────────────────────────────┐
	   2 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⢠⡆⢠⠀⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡇⣾⡇⣸⢰⣧⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣰⡆⣇⣿⣧⣿⣾⣿⣴⣇⡄⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	     │⠀⠀⠀⠀⠀⠀⠀⠀⢠⣰⣿⣷⣿⣿⣿⣿⣿⣿⣿⣿⣿⡇⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	     │⠀⠀⠀⠀⠀⠀⠀⢠⢸⣿⣿⣿⣿⡟⢿⣿⣿⣿⣿⣿⣿⣧⣿⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	     │⠀⠀⠀⠀⠀⠀⠀⣿⣿⣿⣿⢻⣿⡇⢸⣿⣿⣿⣿⣿⣿⣿⣿⡆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	     │⠀⠀⠀⠀⠀⠀⣾⣿⣿⣿⣿⢸⢸⠇⠈⠉⡇⢿⣿⣿⣿⣿⣿⣿⣤⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	     │⠀⠀⠀⠀⠀⢠⣿⣿⣿⣿⢿⠀⢸⠀⠀⠀⠁⠈⠙⠻⠙⢹⣿⣿⣿⣷⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	     │⠀⠀⠀⠀⣠⣿⣿⣿⡿⡿⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⣿⣿⣿⣿⣿⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	     │⠀⠀⠀⢀⣿⣿⠙⠋⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠻⢻⣿⣿⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	     │⠀⠀⠀⣿⣿⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⠋⠏⣿⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	     │⠀⠀⣼⣿⡏⠃⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢿⡇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	     │⠀⢠⣿⠈⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢷⣷⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	     │⢀⡿⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠹⣄⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	   0 │⡼⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠹⡄⠀⠀⠀⠀⠀⠀⠀⠀│
	     └────────────────────────────────────────┘
	     0                                        4

	julia> data = DataGenerator.fnoise(0:0.01:2pi, x->cos(x), x->0)
	Function #59, with noise, 629 obs, 1 vars, 1 target(s)/obs, 629 distinct values

	julia> data |> x->UnicodePlots.lineplot(collect(x.x),x.y)
	      ┌────────────────────────────────────────┐
	    1 │⠉⠓⢦⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡴⠋⠉⠀⠀⠀⠀│
	      │⠀⠀⠀⠳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢠⠞⠀⠀⠀⠀⠀⠀⠀│
	      │⠀⠀⠀⠀⠹⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢠⠏⠀⠀⠀⠀⠀⠀⠀⠀│
	      │⠀⠀⠀⠀⠀⠱⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢠⠏⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	      │⠀⠀⠀⠀⠀⠀⢳⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⡎⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	      │⠀⠀⠀⠀⠀⠀⠀⢧⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡜⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	      │⠀⠀⠀⠀⠀⠀⠀⠈⣇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡸⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	      │⠤⠤⠤⠤⠤⠤⠤⠤⠼⡤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⢴⠧⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠤⠄│
	      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⢹⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠏⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢳⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡞⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢇⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡼⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⣆⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡼⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⣆⠀⠀⠀⠀⠀⠀⠀⠀⡰⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⢦⠀⠀⠀⠀⠀⢀⡼⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	   -1 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠈⠳⢄⣀⣀⡤⠞⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
	      └────────────────────────────────────────┘
	      0                                        7
	```
	"""
	function fnoise(points::AbstractVector, f::Function=x->x, fnoise::Function=x->x+rand())
		X = points
		y = f.(points)
		y = y + fnoise.(y)
		return j4pr.datacell(X,y, name="Function with noise")
	end

end


text:
print("##############################################\n")
print("JaynesCummings.jl Benchmarks - Single Threaded\n")
print("Date: ", Libc.strftime(time()), "\n")
print("Revision:\n")
run(`git rev-parse HEAD`)
print("##############################################\n\n")

versioninfo()
print("\n\n")

# Including the module
using JaynesCummings

# This first run is just used to compile the functions with a basic case
cutoffN = 4
initket = ["|g,1>"]
qubit_freq = 2*pi*6.57E9
resonator_freq = qubit_freq
coupling_freq = 0.5*2*pi*19E6
finalTime = 500E-9
samples = 300
wigner_samples = 10

print("## Execution times for basic case, including build time.\nParameters:\ncutoffN = $cutoffN\nwigner_samples = $wigner_samples\ninitialstate = $initket\n\n")
print("gen_initialstate: ")
initialstate = @time gen_initialstate(cutoffN,initket)
print("gen_hamiltonian: ")
hamiltonian = @time gen_hamiltonian(qubit_freq,resonator_freq,coupling_freq,cutoffN)
print("gen_timeevoarray: ")
time_vec, time_evo_array = @time gen_timeevoarray(hamiltonian,finalTime,samples)
print("calc_qubittimeevo: ")
excited_prob = @time calc_qubittimeevo(initialstate,time_evo_array)
print("calc_photonnumbers: ")
photons = @time calc_photonnumbers(time_vec,real(excited_prob),cutoffN,coupling_freq)
print("calc_densitymatrix_resonator: ")
densitymatrix = @time calc_densitymatrix_resonator(cutoffN,coupling_freq,initialstate,time_vec,time_evo_array)
print("calc_wignerfunction_resonator: ")
wignerfunction = @time calc_wignerfunction_resonator(cutoffN,wigner_samples,coupling_freq,initialstate,time_vec,time_evo_array)

# This second run is much more computationally intensive
cutoffN = 32
wigner_samples = 50
initket = ["|g,0>","0.707106781+0.907106781im|g,3>","1.6025403-0.5im|g,4>","-1.6|g,6>"]
print("\n\n## Execution times for test case, does NOT include build time.\nParameters:\ncutoffN = $cutoffN\nwigner_samples = $wigner_samples\ninitialstate = $initket\n\n")
print("gen_initialstate: ")
initialstate = @time gen_initialstate(cutoffN,initket)
print("gen_hamiltonian: ")
hamiltonian = @time gen_hamiltonian(qubit_freq,resonator_freq,coupling_freq,cutoffN)
print("gen_timeevoarray: ")
time_vec, time_evo_array = @time gen_timeevoarray(hamiltonian,finalTime,samples)
print("calc_qubittimeevo: ")
excited_prob = @time calc_qubittimeevo(initialstate,time_evo_array)
print("calc_photonnumbers: ")
photons = @time calc_photonnumbers(time_vec,real(excited_prob),cutoffN,coupling_freq)
print("calc_densitymatrix_resonator: ")
densitymatrix = @time calc_densitymatrix_resonator(cutoffN,coupling_freq,initialstate,time_vec,time_evo_array)
print("calc_wignerfunction_resonator: ")
wignerfunction = @time calc_wignerfunction_resonator(cutoffN,wigner_samples,coupling_freq,initialstate,time_vec,time_evo_array)

ext: jl


text:
function get_commandline_sh_script(exe_dir, executable_path)
"""
#!/bin/sh

osascript -e "
set RootPath to POSIX path of (path to me)
tell application id \\\"com.apple.terminal\\\"
    do script (\\\"exec '\\\" & RootPath & \\\"$exe_dir/$executable_path'\\\")
    activate
end tell
"
"""
end

# Build a wrapper app that opens a terminal and runs the provided binary.
# Returns a new applet name only if binary_name is already "applet" (to prevent collision).
function build_commandline_app_bundle(builddir, binary_name, appname, verbose)
    println("~~~~~~ Creating commandline-app wrapper applet. ~~~~~~~")

    mkpath(builddir)  # Create builddir if it doesn't already exist.

    app_path = "$builddir/$appname.app"
    exe_dir = "Contents/MacOS"  # Put the binaries next to the applet in MacOS.
    applet_name = "applet"
    if binary_name == applet_name  # Prevent collisions.
        applet_name = "applet_wrapper"
    end
    applet_path = "$app_path/$exe_dir/$applet_name"
    mkpath(dirname(applet_path))

    verbose && println("    Creating wrapper script: $applet_path")
    write(applet_path, get_commandline_sh_script(exe_dir, binary_name))
    run(`chmod +x $applet_path`)

    return applet_name
end


text:
struct MarginalBuffer
    ϕ::Vector{Float64}
    θ::Vector{Float64}
    β::Vector{Float64}
    ind::Vector{Int}
    w::Vector{Float64}
    v::Vector{Float64}
    V::Matrix{Float64}
    Vβ::Vector{Float64}
    Vfβ::Vector{Float64}
    δ::Vector{Float64}
    nbuff1::Vector{Float64}
    ∇β::Vector{Float64}
    Λ::Vector{Float64}
    Q::Vector{Float64}
    α::Matrix{Float64}
    ∇α::Vector{Float64}
    init::Vector{Float64}
    n::Int
end

function MarginalBuffer(n::Int)
    ϕ = Vector{Float64}(9)
    θ = Vector{Float64}(7)#Representations from nested cubic.
    β = Vector{Float64}(10)#Representation multiplied by V.
    ind = Vector{Int}(n)
    w = Vector{Float64}(n)
    v = Vector{Float64}(n)
    V =  ones(10,n)
    Vβ = Vector{Float64}(n)
    Vfβ = Vector{Float64}(n)
#    VFβ = Vector{Float64}(n)
    δ = Vector{Float64}(n)
    nbuff1 = Vector{Float64}(n)
    ∇β = Vector{Float64}(10)
    Λ = Vector{Float64}(2)
    Q = Vector{Float64}(3)
    α = zeros(7,10)
    ∇α = Vector{Float64}(7)
    init = [0.0, 0.0, 0.0, 0.0, -2.0, 0.5, -0.2, -5.0, 3.5]

    MarginalBuffer(ϕ, θ, β, ind, w, v, V, Vβ, Vfβ, δ, nbuff1, ∇β, Λ, Q, α, ∇α, init, n)
end


text:
using NWNTools, StaticArrays, BenchmarkTools, LinearAlgebra, SparseArrays, LightGraphs, SimpleWeightedGraphs
dims = SA[20.0,20.0];
ρₙ = 0.4;
l = 7.0;
ρₑ = 22.2/1000; #put in micro Ohm metres
D = 0.02;
V₊ = 5.0;
Rⱼ = 11.0;
Rₑ = 7.0;
lines = rand_lines(dims, round(Int, prod(dims)*ρₙ), l);
props = repeat([WireProp(ρₑ,D)],length(lines));
grnds = [Line(SA[dims[1],0.0],dims)];
elecs = [Line(SA[0.0,0.0],SA[0.0,dims[2]])];

nwn = NWN_circuit{:MNR}(lines, props, elecs, [V₊], grnds, dims, Rⱼ, Rₑ);


text:
## Internal API
# get names
"""
    names(env::Fym)

Get the symbols of sub-environments.

```jldoctest
# Examples
julia> using LazyFym
julia> struct Env1 <: Fym
       end

julia> struct Env2 <: Fym
       end

julia> struct Env3 <: Fym
       end

julia> struct Env12 <: Fym
           env1::Env1
           env2::Env2
       end

julia> struct Env12_3 <: Fym
           env12::Env12
           env3::Env3
       end

julia> env12_3 = Env12_3(Env12(Env1(), Env2()), Env3())
Env12_3(Env12(Env1(), Env2()), Env3())

julia> LazyFym.names(env12_3)
2-element Array{Symbol,1}:
 :env12
 :env3

julia> LazyFym.names(env12_3.env12)
2-element Array{Symbol,1}:
 :env1
 :env2
```
"""
function Base.names(env::Fym)
    return [name for name in fieldnames(typeof(env)) if typeof(getfield(env, name)) <: Fym]
end

function preprocess(env::Fym, x0)
    env_size_nt = size(env, x0)
    env_flatten_length = flatten_length(env, x0)
    env_index_nt = index(env, x0, 1:env_flatten_length)
    return env_index_nt, env_size_nt
end

# raw view
function raw(env, x)
    env_names = names(env)
    if env_names == []
        return x
    else
        _x = env_names |> MapCat(name -> raw(getfield(env, name), x[name])) |> collect
        return _x
    end
end

# processed view
function process(_x, env_index_nt, env_size_nt)
    if typeof(env_index_nt) <: AbstractRange
        if env_size_nt == ()
            return _x[env_index_nt][1]  # scalar
        else
            return reshape(_x[env_index_nt], env_size_nt...)
        end
    else
        index_names = keys(env_index_nt)
        processed_values = index_names |> Map(name -> process(_x, env_index_nt[name], env_size_nt[name]))
        return (; zip(index_names, processed_values)...)
    end
end

# size
function Base.size(env::Fym, x0)
    env_names = names(env)
    if env_names == []
        return size(x0)
    else
        env_sizes = env_names |> Map(name -> size(getfield(env, name), x0[name]))
        return (; zip(env_names, env_sizes)...)  # NamedTuple
    end
end

# flatten length
function flatten_length(env::Fym, x0)
    env_names = names(env)
    if env_names == []
        return prod(size(x0))
    else
        return env_names |> Map(name -> flatten_length(getfield(env, name), x0[name])) |> sum
    end
end

# index
function index(env::Fym, x0, _range)
    env_names = names(env)
    if env_names == []
        return _range
    else
        env_accumulated_flatten_lengths = env_names |> Map(name -> flatten_length(getfield(env, name), x0[name])) |> Scan(+) |> collect
        range_first = first(_range)
        env_ranges_tmp = (range_first-1) .+ [0, env_accumulated_flatten_lengths...] |> Consecutive(length(env_accumulated_flatten_lengths); step=1)
        env_ranges = zip(env_ranges_tmp...) |> MapSplat((x, y) -> x+1:y)
        env_indices = zip(env_names, env_ranges) |> MapSplat((name, range) -> index(getfield(env, name), x0[name], range))
        return (; zip(env_names, env_indices)...)  # NamedTuple
    end
end


text:
using ADMPS
using ADMPS: num_grad
using ADMPS: leftenv,rightenv,norm_FL,norm_FR
using ChainRulesCore
using CUDA
using KrylovKit
using LinearAlgebra
using OMEinsum
using Random
using Test
using Zygote
CUDA.allowscalar(false)

@testset "Zygote with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64]
    a = atype(randn(2,2))
    @test Zygote.gradient(norm, a)[1] ≈ num_grad(norm, a)

    foo1 = x -> sum(atype(Float64[x 2x; 3x 4x]))
    @test Zygote.gradient(foo1, 1)[1] ≈ num_grad(foo1, 1)
end

@testset "Zygote.@ignore" begin
    function foo2(x)
        return x^2
    end
    function foo3(x)
        return x^2 + Zygote.@ignore x^3
    end
    @test foo2(1) != foo3(1)
    @test Zygote.gradient(foo2,1)[1] ≈ Zygote.gradient(foo3,1)[1]
end

@testset "linsolve with $atype{$dtype}" for atype in [Array], dtype in [Float64]
    Random.seed!(100)
    D,d = 2^2,2
    A = atype(rand(D,d,D))
    工 = ein"asc,bsd -> abcd"(A,conj(A))
    λLs, Ls, info = eigsolve(L -> ein"ab,abcd -> cd"(L,工), atype(rand(D,D)), 1, :LM)
    λL, L = λLs[1], Ls[1]
    λRs, Rs, info = eigsolve(R -> ein"abcd,cd -> ab"(工,R), atype(rand(D,D)), 1, :LM)
    λR, R = λRs[1], Rs[1]

    dL = atype(rand(D,D))
    dL -= ein"ab,ab -> "(L,dL)[] * L
    @test ein"ab,ab ->  "(L,dL)[] ≈ 0 atol = 1e-9
    ξL, info = linsolve(R -> ein"abcd,cd -> ab"(工,R), dL, -λL, 1)
    @test ein"ab,ab -> "(ξL,L)[] ≈ 0 atol = 1e-9

    dR = atype(rand(D,D))
    dR -= ein"ab,ab -> "(R,dR)[] * R
    @test ein"ab,ab -> "(R,dR)[] ≈ 0 atol = 1e-9
    ξR, info = linsolve(L -> ein"ab,abcd -> cd"(L,工), dR, -λR, 1)
    @test ein"ab,ab -> "(ξR,R)[] ≈ 0 atol = 1e-9
end

@testset "loop_einsum mistake with $atype" for atype in [Array, CuArray]
    Random.seed!(100)
    D = 10
    A = atype(rand(D,D,D))
    B = atype(rand(D,D))
    function foo(x)
        C = A * x
        D = B * x
        E = ein"abc,abc -> "(C,C)
        F = ein"ab,ab -> "(D,D)
        return Array(E)[]/Array(F)[]
    end
    @test Zygote.gradient(foo, 1)[1] ≈ num_grad(foo, 1) atol = 1e-8
end

@testset "leftenv and rightenv with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64]
    Random.seed!(100)
    d = 2
    D = 3
    Au = atype(rand(dtype,D,d,D))
    Ad = atype(rand(dtype,D,d,D))

    S = atype(rand(ComplexF64,D,d,D,D,d,D))
    M = atype(rand(dtype,d,d,d,d))
    function foo1(M)
        _,FL = leftenv(Au, Ad, M)
        A = ein"(abc,abcdef),def -> "(FL,S,FL)
        B = ein"abc,abc -> "(FL,FL)
        return norm(Array(A)[]/Array(B)[])
    end
    @test Zygote.gradient(foo1, M)[1] ≈ num_grad(foo1, M) atol = 1e-8

    function foo2(M)
        _,FR = rightenv(Au, Ad, M)
        A = ein"(abc,abcdef),def -> "(FR,S,FR)
        B = ein"abc,abc -> "(FR,FR)
        return norm(Array(A)[]/Array(B)[])
    end
    @test Zygote.gradient(foo2, M)[1] ≈ num_grad(foo2, M) atol = 1e-8
end

@testset "norm_FL and norm_FR with $atype{$dtype}" for atype in [Array, CuArray], dtype in [Float64, ComplexF64]
    Random.seed!(100)
    d = 2
    D = 3
    Au = atype(rand(dtype,D,d,D))
    Ad = atype(rand(dtype,D,d,D))

    S = atype(rand(dtype,D,D,D,D))
    function foo1(Au)
        _,FL = norm_FL(Au, Ad)
        A = ein"(ab,abcd),cd -> "(FL,S,FL)
        B = ein"ab,ab -> "(FL,FL)
        return norm(Array(A)[]/Array(B)[])
    end
    @test Zygote.gradient(foo1, Au)[1] ≈ num_grad(foo1, Au) atol = 1e-8

    function foo2(Ad)
        _,FR = norm_FR(Au, Ad)
        A = ein"(ab,abcd),cd -> "(FR,S,FR)
        B = ein"ab,ab -> "(FR,FR)
        return norm(Array(A)[]/Array(B)[])
    end
    @test Zygote.gradient(foo2, Ad)[1] ≈ num_grad(foo2, Ad) atol = 1e-8
end


text:
streams_spec =
[
	(:_type, AbstractString, nothing),
	(:expr, AbstractString, nothing),
	(:scale, AbstractString, nothing),
	#(:predicate, AbstractString, nothing)
]

primitivefactory(:VegaStream, streams_spec)

max_stars_repo_name: jpthor/vega


text:
using SparseArrays
"""
```
steadystate!(m::TwoAssetHANK)
```

Calculates the model's steady-state values. `steadystate!(m)` must be called whenever
the parameters of `m` are updated.
"""
function steadystate!(m::TwoAssetHANK)
    # Read in parameters
    aalpha = m[:aalpha].value
    ddelta = m[:ddelta].value
    ddeath = m[:ddeath].value
    rrho   = m[:rrho].value
    chi0   = m[:chi0].value
    chi1   = m[:chi1].value
    chi2   = m[:chi2].value
    a_lb   = m[:a_lb].value
    pam    = m[:pam].value
    xxi    = m[:xxi].value
    ggamma = m[:ggamma].value
    tau_I  = m[:tau_I].value
    trans  = m[:trans].value

    # Set liquid rates
    r_b      = m[:r_b_SS].value
    r_b_borr = m[:r_b_borr_SS].value
    lambda   = get_setting(m, :lambda)
    K_liquid = get_setting(m, :K_liquid)

    # Output flags
    aggregate_variables        = get_setting(m, :aggregate_variables)
    distributional_variables   = get_setting(m, :distributional_variables)
    distributional_variables_1 = get_setting(m, :distributional_variables_1)

    # Approximation parameters
    maxit_HJB = get_setting(m, :maxit_HJB)
    crit_HJB  = get_setting(m, :crit_HJB)
    Delta     = get_setting(m, :Delta)

    maxit_HIS = get_setting(m, :maxit_HIS)
    crit_HIS  = get_setting(m, :crit_HIS)
    start_HIS = get_setting(m, :start_HIS)

    maxit_KFE = get_setting(m, :maxit_KFE)
    crit_KFE  = get_setting(m, :crit_KFE)
    Δ_KFE     = get_setting(m, :Delta_KFE)

    maxit_KL  = get_setting(m, :maxit_KL)
    crit_KL   = get_setting(m, :crit_KL)
    relax_KL  = get_setting(m, :relax_KL)

    I      = get_setting(m, :I)
    J      = get_setting(m, :J)
    I_g    = get_setting(m, :I_g)
    J_g    = get_setting(m, :J_g)
    N      = get_setting(m, :N)

    # Read in grids
    a      = get_setting(m, :a)
    b      = get_setting(m, :b)
    a_g    = get_setting(m, :a_g)
    b_g    = get_setting(m, :b_g)
    y      = real(vec(get_setting(m, :y)))
    y_dist = get_setting(m, :y_dist)
    y_mean = get_setting(m, :y_mean)
    KL     = get_setting(m, :KL_0)
    interp_decision = get_setting(m, :interp_decision)

    # Compute prices associated with initial guess of KL
    w	= (1 - aalpha) * (KL ^ aalpha)
    r_a	= aalpha * (KL ^ (aalpha - 1)) - ddelta

    # Construct problem functions
    util, deposit, cost = construct_problem_functions(ggamma, chi0, chi1, chi2, a_lb)

    # Initialize difference grids
    a_g_grid, b_g_grid, dab_g_tilde_grid, dab_g_tilde = set_grids(a, b, a_g, b_g, N)

    dab_g_tilde_vec = repeat(dab_g_tilde, N)
    r_b_vec         = r_b .* (b .>= 0) + r_b_borr .* (b .< 0)

    b_vec = repeat(repeat(b, outer=J), outer=N)
    a_vec = repeat(repeat(a, inner=I), outer=N)
    y_vec = repeat(y, inner=I*J)

    y_grid   = reshape(repeat(y, inner=I*J), I, J, N)
    y_g_grid = reshape(repeat(y, inner=I_g*J_g), I_g, J_g, N)

    # Initialization for consumption and value function
    c_0 = (1-xxi) * w * y_vec .+ (r_a + ddeath*pam) .* a_vec +
        (r_b_borr + ddeath * pam) .* b_vec .+ trans - tau_I * w * y_vec

    V_0	= reshape(1/(rrho + ddeath) * util.(c_0), I, J, N)

    # Initial distribution
    gg0 = zeros(I_g, J_g, N)
    gg0[b .== 0, 1, :] = vec(y_dist)
    gg0 = vec(gg0)
    gg0 = (gg0 ./ sum(gg0)) ./ dab_g_tilde_vec # ensure integration to 1
    gg  = gg0

    AT  = Array{Float64}(undef, 0, 0)
    Vn1 = Array{Float64}(undef, I, J, N)
    Vn  = Array{Float64}(undef, I, J, N)

    gg_tilde = Array{Float64}(undef, 0, 0)

    g   = Array{Float64}(undef, 0, 0)
    c   = Array{Float64,3}(undef, I, J, N)
    s   = Array{Float64,3}(undef, I, J, N)
    u   = Array{Float64,3}(undef, I, J, N)
    d   = Array{Float64,3}(undef, I, J, N)
    c_g = Array{Float64,3}(undef, I_g, J_g, N)

    K_supply, L_supply, KL_supply = 0.0, 0.0, 0.0

    #----------------------------------------------------------------
    # Iterate on KL to find steady state
    #----------------------------------------------------------------
    for ii = 1 : maxit_KL

        # Derive aggregates, given KL
	    w   = (1 - aalpha) * (KL ^ aalpha)
 	    r_a = aalpha * (KL ^ (aalpha - 1)) - ddelta

	    # Store current value function
	    Vn = V_0

	    #----------------------------------------------------------------
	    # Solve HJB
	    #----------------------------------------------------------------
        @time for nn = 1 : maxit_HJB
            c, s, d = solve_hjb(Vn, a_lb, ggamma, ddeath, pam, trans, xxi,
                                tau_I, 0.0, w, r_b_vec, y, a, b, cost, util, deposit)
            u = util.(c)

            # Interpolate
            d_g = reshape(interp_decision * vec(d), I_g, J_g, N)
            s_g = reshape(interp_decision * vec(s), I_g, J_g, N)
            c_g = reshape(interp_decision * vec(c), I_g, J_g, N)

            A, AT = transition(ddeath, pam, xxi, w, a_lb, 0.0, d, d_g, s, s_g,
                               r_a, a, a_g, b, b_g, y, cost)

            #------------------------------------------------------------
            # Update value function, Howard improvement steps
            #------------------------------------------------------------
            Vn1 = update_value_fn(A, Vn, lambda, u, Delta, rrho, ddeath)

            if nn >= start_HIS
                for jj = 1:maxit_HIS
                    Vn2    = update_value_fn(A, Vn1, lambda, u, Delta, rrho, ddeath)
                    VHIS_Δ = Vn2 - Vn1
                    Vn1    = Vn2
                    dist   = maximum(abs.(VHIS_Δ))
                    if dist < crit_HIS
                        break
                    end
                end
            end

            # Check for convergence
            V_Δ  = Vn1 - Vn
            Vn   = Vn1
            dist = maximum(abs.(V_Δ))
            @show dist

            if dist < crit_HJB
                 println("Value Function Converged, Iteration = ", nn)
                break
            end
        end

        # Store final value function
        V_0 = Vn

        #----------------------------------------------------------------
        # Compute stationary distribution and update aggregate KL
        #----------------------------------------------------------------

        # Find new stationary distribution associated with decision rules
        A   = AT
        λ0  = lambda - diagm(0 => diag(lambda)) # transition matrix with diagonal killed
        λ0p = λ0'

        gg_tilde = dab_g_tilde_vec .* gg
        gg1      = Array{Float64,2}(undef, I_g * J_g, N)
        g        = Array{Float64,3}(undef, I_g,  J_g, N)

        K_supply, L_supply, KL_supply = 0.0, 0.0, 0.0

        #----------------------------------------------------------------
        # Iteration to solve Kolmogorov Forward equation
        #----------------------------------------------------------------
        @time for nn = 1:maxit_KFE

            gg_tilde  = dab_g_tilde_vec .* gg
            gg_tilde2 = reshape(gg_tilde, I_g*J_g, N)
            gg1       = Array{Float64,2}(undef, I_g * J_g, N)

            for kk = 1:N
                Ak = A[1+(kk-1)*(I_g*J_g):kk*(I_g*J_g), 1+(kk-1)*(I_g*J_g):kk*(I_g*J_g)]

                death_inflow = zeros(I_g, J_g)
                death_inflow[b .== 0, 1] .= sum(gg_tilde2[:,kk])

                gk_sum = vec(sum(broadcast(*, λ0p[kk,:], gg_tilde2'), dims=1))

                gg1[:,kk] = add_diag(add_diag(-Δ_KFE*Ak, 1.0),
                                     -Δ_KFE * (lambda[kk,kk] - ddeath))' \
                                         (gg_tilde2[:,kk] + Δ_KFE * gk_sum +
                                          Δ_KFE * ddeath * vec(death_inflow))
            end

            gg11 = (vec(gg1) ./ sum(gg1)) ./ dab_g_tilde_vec
            dist = maximum(abs.(gg11-gg))

            if mod(nn,100) == 0
                println("Iteration ", nn, ", distance is ", dist)
            end
            if dist < crit_KFE
                gg = gg11
                g  = reshape(gg, I_g, J_g, N)

                println("Distribution Converged, Iteration = ", nn)
                break
            end
            gg = gg11
            g  = reshape(gg, I_g, J_g, N)
        end

        #------------------------------------------------------------
        # Update guess of Capital, Labor Supply
        #------------------------------------------------------------
        ggdab_g = g .* dab_g_tilde_grid

        K_supply  = sum((a_g_grid .+ (K_liquid ? b_g_grid : 0.0)) .* ggdab_g)
        L_supply  = sum(y_g_grid .* ggdab_g)
        KL_supply = K_supply / L_supply

        # Check for convergence and update
        KL_Δ = KL - KL_supply
        println("The current gap is ", KL_Δ)

        if abs(KL_Δ) > crit_KL
            KL = relax_KL * KL_supply + (1 - relax_KL) * KL
        else
            println("I have found the steady state, Iteration = ", ii)
            break
        end
    end

    ccu                  = kron(lambda, my_speye(I_g * J_g))
    A                    = AT + ccu
    dab_small            = reshape(dab_g_tilde_grid[:,:,1], I_g*J_g, 1)
    loc                  = findall(!iszero, b .== 0)
    dab_small            = dab_small ./ dab_small[loc] * ddeath
    dab_small[loc]      .= 0.0
    death_process        = -ddeath * my_speye(I_g*J_g)
    death_process[loc,:] = vec(dab_small)
    death_process        = kron(my_speye(N), death_process)

    g_new   = A' * gg_tilde
    g_new   = (1.0 ./ dab_g_tilde_vec) .* g_new + death_process * gg
    g_new   = reshape(g_new ,I_g, J_g, N)
    g_new_a = sum(sum(g_new, dims=1), dims=3)
    #compute_savings()

    # Rename variables in steady state
    m[:KL_SS]           = real(KL)
    m[:r_a_SS]          = aalpha * (m[:KL_SS].value ^ (aalpha - 1)) - ddelta
    m[:w_SS]            = (1 - aalpha) * (m[:KL_SS].value ^ aalpha)
    m[:K_SS]            = real(m[:KL_SS].value * L_supply)
    m[:u_SS]            = u
    m[:c_SS]            = c
    m[:d_SS]            = d
    m[:V_SS]            = Vn
    m[:g_SS]            = g
    m[:C_SS]            = sum(g .* c_g .* dab_g_tilde_grid)
    m[:C_Var_SS]        = sum(g .* log.(c_g).^2 .* dab_g_tilde_grid) -
                                 sum(g .* log.(c_g) .* dab_g_tilde_grid)^2
    m[:I_SS]            = dot(g_new_a, a_g)
    m[:B_SS]            = sum(g .* b_g_grid .* dab_g_tilde_grid)
    m[:Y_SS]            = real((K_supply ^ aalpha) * (L_supply ^ (1 - aalpha)))
    m[:n_SS]            = real(L_supply)

    r_b_g_grid          = r_b .* (b_g_grid .>= 0) + r_b_borr .* (b_g_grid .< 0)
    m[:earn_SS]         = log.((1-tau_I) * w * y_g_grid + b_g_grid .*
                               (r_b_g_grid .+ ddeath*pam) .+ trans .+ a_g_grid .*
                               (r_a .+ ddeath*pam))

    m[:earn_Var_SS]     = sum(g .* m[:earn_SS].value .^ 2 .* dab_g_tilde_grid) -
                                sum(g .* m[:earn_SS].value .* dab_g_tilde_grid)^2

    a_g_0pos = get_setting(m, :a_g_0pos)
    b_g_0pos = get_setting(m, :b_g_0pos)

    ###
    # Consumption by hand-to-mouth status
    ###
    WHTM_indicator      = zeros(I_g, J_g, N)
    WHTM_indicator[b_g_0pos:b_g_0pos+1, a_g_0pos+2:end,:] .= 1.0

    m[:WHTM_indicator] = WHTM_indicator
    m[:WHTM_SS]        = sum(g[:] .* WHTM_indicator[:] .* dab_g_tilde_grid[:])
    m[:C_WHTM_SS]      = sum(WHTM_indicator[:] .* c_g[:] .* g[:] .* dab_g_tilde_grid[:])

    PHTM_indicator     = zeros(I_g, J_g, N)
    PHTM_indicator[b_g_0pos:b_g_0pos+1,a_g_0pos:a_g_0pos+2:end,:] .= 1.0
    m[:PHTM_indicator] = PHTM_indicator
    m[:PHTM_SS]        = sum(g[:] .* PHTM_indicator[:] .* dab_g_tilde_grid[:])
    m[:C_PHTM_SS]      = sum(c_g[:] .* g[:] .* PHTM_indicator[:] .* dab_g_tilde_grid[:])

    NHTM_indicator      = zeros(I_g,J_g,N)
    NHTM_indicator[(b_g_0pos+3):end,2:end,:] .= 1.0

    m[:NHTM_indicator] = NHTM_indicator
    m[:NHTM_SS]        = sum(g[:] .* NHTM_indicator[:] .* dab_g_tilde_grid[:])
    m[:C_NHTM_SS]      = sum(c_g[:] .* g[:] .* NHTM_indicator[:] .* dab_g_tilde_grid[:]) /
        m[:NHTM_SS].value

    n_v = get_setting(m, :n_v)
    n_g = get_setting(m, :n_g)
    n_p = get_setting(m, :n_p)
    n_Z = get_setting(m, :n_Z)

    # Collect variables into vector
    vars_SS                  = zeros(n_v + n_g + n_p + n_Z)
    vars_SS[1:n_v]           = vec(m[:V_SS].value)
    gg_SS                    = vec(m[:g_SS].value)
    m[:gg_SS]                = gg_SS

    vars_SS[n_v+1:n_v+n_g] = gg_SS[1:I_g*J_g*N-1]
    vars_SS[n_v+n_g+1]     = m[:K_SS].value
    vars_SS[n_v+n_g+2]     = m[:r_b_SS].value

    if aggregate_variables == 1
        vars_SS[n_v+n_g+3] = m[:Y_SS].value
        vars_SS[n_v+n_g+4] = m[:C_SS].value
    elseif distributional_variables == 1
        vars_SS[n_v+n_g+3] = m[:C_Var_SS].value
        vars_SS[n_v+n_g+4] = m[:earn_Var_SS].value
    elseif distributional_variables_1 == 1
        vars_SS[n_v+n_g+3] = m[:C_WHTM_SS].value
        vars_SS[n_v+n_g+4] = m[:C_PHTM_SS].value
    end
    vars_SS[n_v+n_g+n_p+1,1] = 0.0

    # Compute illiquid wedge
    m[:illiquid_wedge] = m[:r_a_SS].value - m[:r_b_SS].value

    # Save SS variables
    m[:vars_SS] = vars_SS

    return m
end


text:
module CSDP

# Try to load the binary dependency
if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("CSDP not properly installed. Please run Pkg.build(\"CSDP\")")
end

export Blockmatrix

include("blockmat.h.jl")
include("blockmat.jl")
include("declarations.h.jl")
include("declarations.jl")
include("debug-mat.jl")
include("CSDPSolverInterface.jl")

end # module


text:
function CursorSchema(stmt::QueryStmt)

    columns_count = ncol(stmt)
    column_query_info = columns_info(stmt)

    column_names_index = Dict{String, Int}()

    for (column_index, query_info) in enumerate(column_query_info)
        column_names_index[column_name(query_info)] = column_index
    end

    return CursorSchema(column_query_info, column_names_index)
end

@inline ncol(schema::CursorSchema) = length(schema.column_names_index)
@inline ncol(row::ResultSetRow) = ncol(row.cursor)
@inline ncol(cursor::Cursor) = ncol(cursor.stmt)
@inline ncol(rs::ResultSet) = ncol(rs.schema)

@inline nrow(rs::ResultSet) = length(rs.rows)

@inline Base.isempty(rs::ResultSet) = isempty(rs.rows)
@inline Base.size(rs::ResultSet) = ( nrow(rs), ncol(rs) )

"""
    names(schema::CursorSchema) :: Vector{String}
    names(row::ResultSetRow) :: Vector{String}
    names(cursor::Cursor) :: Vector{String}
    names(rs::ResultSet) :: Vector{String}

Return a freshly allocated `Vector{String}` of names of columns contained in schema, row, cursor, or result set.
"""
Base.names(schema::CursorSchema) = [ column_name(orainfo) for orainfo in schema.column_query_info ]
Base.names(cursor::Cursor) = names(cursor.schema)
Base.names(row::ResultSetRow) = names(row.schema)
Base.names(rs::ResultSet) = names(rs.schema)

@inline stmt(cursor::Cursor) = cursor.stmt

function Cursor(stmt::QueryStmt)
    schema = CursorSchema(stmt)
    return Cursor(stmt, schema)
end

function Base.iterate(cursor::Cursor, nil::Nothing=nothing)
    fetch_result = fetch(cursor.stmt)
    if fetch_result.found
        current_row = ResultSetRow(cursor)
        return (current_row, nothing)
    else
        return nothing
    end
end

function query(f::Function, conn::Connection, sql::AbstractString;
               scrollable::Bool=false,
               tag::AbstractString="",
               fetch_array_size::Integer=ORA_DEFAULT_FETCH_ARRAY_SIZE,
               exec_mode::OraExecMode=ORA_MODE_EXEC_DEFAULT)

    stmt(conn, sql, scrollable=scrollable, tag=tag, fetch_array_size=fetch_array_size) do stmt
        execute(stmt, exec_mode=exec_mode)
        cursor = Cursor(stmt)
        f(cursor)
    end

    nothing
end

function execute_and_fetch_all!(stmt::Stmt; exec_mode::OraExecMode=ORA_MODE_EXEC_DEFAULT) :: ResultSet
    local schema::CursorSchema
    rows = Vector{ResultSetRow}()

    query(stmt, exec_mode=exec_mode) do cursor
        schema = cursor.schema

        for row in cursor
            push!(rows, row)
        end
    end

    return ResultSet(schema, rows)
end

function query(f::Function, stmt::Stmt; exec_mode::OraExecMode=ORA_MODE_EXEC_DEFAULT)
    execute(stmt, exec_mode=exec_mode)
    f(Cursor(stmt))
    nothing
end

function query(stmt::Stmt; exec_mode::OraExecMode=ORA_MODE_EXEC_DEFAULT) :: ResultSet
    return execute_and_fetch_all!(stmt, exec_mode=exec_mode)
end

function query(conn::Connection, sql::AbstractString;
               scrollable::Bool=false,
               tag::AbstractString="",
               fetch_array_size::Integer=ORA_DEFAULT_FETCH_ARRAY_SIZE,
               exec_mode::OraExecMode=ORA_MODE_EXEC_DEFAULT) :: ResultSet

    local rs::ResultSet

    stmt(conn, sql, scrollable=scrollable, tag=tag, fetch_array_size=fetch_array_size) do stmt
        rs = execute_and_fetch_all!(stmt, exec_mode=exec_mode)
    end

    return rs
end

function ResultSetRow(cursor::Cursor)
    data = Vector{Any}()

    # parse data from Cursor
    for column_index in 1:ncol(cursor)
        oracle_value = query_oracle_value(stmt(cursor), column_index)
        push!(data, parse_oracle_value(oracle_value))
    end

    return ResultSetRow(cursor.schema, data)
end

Base.length(row::ResultSetRow) = length(row.data)
Base.isempty(row::ResultSetRow) = isempty(row.data)

@inline function Base.getindex(row::ResultSetRow, column_index::Integer)
    check_inbounds(row, column_index)
    @inbounds return row.data[column_index]
end

@inline function Base.getindex(row::ResultSetRow, column_name::AbstractString)
    check_inbounds(row, column_name)
    column_index = row.schema.column_names_index[column_name]
    @inbounds return row.data[column_index]
end

@inline function check_inbounds(schema::CursorSchema, column::Integer)
    @assert 0 < column <= ncol(schema) "Column $column not found."
end

@inline function check_inbounds(schema::CursorSchema, column::AbstractString)
    @assert haskey(schema.column_names_index, column) "Column $column not found."
end

@inline check_inbounds(row::ResultSetRow, column::Union{AbstractString,Integer}) = check_inbounds(row.schema, column)

@inline function check_inbounds(rs::ResultSet, row::Integer, column::Union{AbstractString, Integer})
    check_inbounds(rs.schema, column)
    @assert 0 < row <= nrow(rs) "Row $row not found."
end

has_possibly_more_rows(r::FetchRowsResult) = Bool(r.more_rows)

Base.IteratorSize(::Cursor) = Base.SizeUnknown()
Base.eltype(::Cursor) = ResultSetRow

function fetchrow(stmt::QueryStmt) :: Union{Nothing, ResultSetRow}
    fetch_result = fetch(stmt)

    if fetch_result.found
        return ResultSetRow(Cursor(stmt))
    else
        return nothing
    end
end

function Base.getindex(rs::ResultSet, row::Integer, column::Union{AbstractString,Integer})
    check_inbounds(rs, row, column)
    return rs.rows[row][column]
end

function Base.getindex(rs::ResultSet, ::Colon, column::Union{AbstractString,Integer})
    column_data = Vector()

    !isempty(rs) && check_inbounds(rs.rows[1], column)

    for row in rs.rows
        push!(column_data, row[column])
    end

    return column_data
end

function Base.lastindex(rs::ResultSet, d::Integer)
    @assert d == 1 || d == 2 "Invalid dimension for ResultSet: $d."

    # rows
    if d == 1
        return nrow(rs)
    else
        # columns
        return ncol(rs)
    end
end

function Base.getindex(rs::ResultSet, row_range::UnitRange{T}, column::Union{AbstractString,Integer}) where {T<:Integer}
    column_data = Vector()

    for r in row_range
        check_inbounds(rs, r, column)
        push!(column_data, rs[r, column])
    end

    return column_data
end


text:
@cache struct EulerCache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
end

@cache struct SplitEulerCache{uType,rateType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  k::rateType
  fsalfirst::rateType
end

function alg_cache(alg::SplitEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  SplitEulerCache(u,uprev,similar(u),zero(rate_prototype),zero(rate_prototype))
end

struct SplitEulerConstantCache <: OrdinaryDiffEqConstantCache end

alg_cache(alg::SplitEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = SplitEulerConstantCache()

function alg_cache(alg::Euler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  EulerCache(u,uprev,similar(u),zero(rate_prototype),zero(rate_prototype))
end

struct EulerConstantCache <: OrdinaryDiffEqConstantCache end

alg_cache(alg::Euler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = EulerConstantCache()

@cache struct HeunCache{uType,rateType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  atmp::uNoUnitsType
  k::rateType
  fsalfirst::rateType
end

@cache struct RalstonCache{uType,rateType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  tmp::uType
  atmp::uNoUnitsType
  k::rateType
  fsalfirst::rateType
end

function alg_cache(alg::Heun,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  HeunCache(u,uprev,similar(u),similar(u, uEltypeNoUnits),zero(rate_prototype),
            zero(rate_prototype))
end

function alg_cache(alg::Ralston,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  RalstonCache(u,uprev,similar(u),similar(u, uEltypeNoUnits),zero(rate_prototype),
               zero(rate_prototype))
end

struct HeunConstantCache <: OrdinaryDiffEqConstantCache end

alg_cache(alg::Heun,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = HeunConstantCache()

struct RalstonConstantCache <: OrdinaryDiffEqConstantCache end

alg_cache(alg::Ralston,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = RalstonConstantCache()

@cache struct MidpointCache{uType,rateType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  tmp::uType
  atmp::uNoUnitsType
  fsalfirst::rateType
end

struct MidpointConstantCache <: OrdinaryDiffEqConstantCache end

function alg_cache(alg::Midpoint,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u); atmp = similar(u, uEltypeNoUnits)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  MidpointCache(u,uprev,k,tmp,atmp,fsalfirst)
end

alg_cache(alg::Midpoint,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = MidpointConstantCache()

@cache struct RK4Cache{uType,rateType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k₂::rateType
  k₃::rateType
  k₄::rateType
  k::rateType
  tmp::uType
  atmp::uNoUnitsType
end

struct RK4ConstantCache <: OrdinaryDiffEqConstantCache end

function alg_cache(alg::RK4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  k₁ = zero(rate_prototype)
  k₂ = zero(rate_prototype)
  k₃ = zero(rate_prototype)
  k₄ = zero(rate_prototype)
  k  = zero(rate_prototype)
  tmp = similar(u); atmp = similar(u, uEltypeNoUnits)
  RK4Cache(u,uprev,k₁,k₂,k₃,k₄,k,tmp,atmp)
end

alg_cache(alg::RK4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = RK4ConstantCache()

@cache struct BS3Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  fsalfirst::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::BS3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = BS3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  BS3Cache(u,uprev,k1,k2,k3,k4,utilde,tmp,atmp,tab)
end

alg_cache(alg::BS3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = BS3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct OwrenZen3Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::OwrenZen3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = OwrenZen3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  OwrenZen3Cache(u,uprev,k1,k2,k3,k4,utilde,tmp,atmp,tab)
end

alg_cache(alg::OwrenZen3,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = OwrenZen3ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct OwrenZen4Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::OwrenZen4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = OwrenZen4ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  OwrenZen4Cache(u,uprev,k1,k2,k3,k4,k5,k6,utilde,tmp,atmp,tab)
end

alg_cache(alg::OwrenZen4,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = OwrenZen4ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct OwrenZen5Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::OwrenZen5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = OwrenZen5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  k7 = zero(rate_prototype)
  k8 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  OwrenZen5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,utilde,tmp,atmp,tab)
end

alg_cache(alg::OwrenZen5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = OwrenZen5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct BS5Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  k8::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::BS5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = BS5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  k7 = zero(rate_prototype)
  k8 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  BS5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,k8,utilde,tmp,atmp,tab)
end

alg_cache(alg::BS5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = BS5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct Tsit5Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

@cache struct RK46NLCache{uType,rateType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  tmp::uType
  fsalfirst::rateType
  tab::TabType
end

struct RK46NLConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  α2::T
  α3::T
  α4::T
  α5::T
  α6::T
  β1::T
  β2::T
  β3::T
  β4::T
  β5::T
  β6::T
  c2::T2
  c3::T2
  c4::T2
  c5::T2
  c6::T2

  function RK46NLConstantCache(::Type{T}, ::Type{T2}) where {T,T2}
    α2 = T(-0.737101392796)
    α3 = T(-1.634740794343)
    α4 = T(-0.744739003780)
    α5 = T(-1.469897351522)
    α6 = T(-2.813971388035)
    β1 = T(0.032918605146)
    β2 = T(0.823256998200)
    β3 = T(0.381530948900)
    β4 = T(0.200092213184)
    β5 = T(1.718581042715)
    β6 = T(0.27)
    c2 = T2(0.032918605146)
    c3 = T2(0.249351723343)
    c4 = T2(0.466911705055)
    c5 = T2(0.582030414044)
    c6 = T2(0.847252983783)
    new{T,T2}(α2, α3, α4, α5, α6, β1, β2, β3, β4, β5, β6, c2, c3, c4, c5, c6)
  end
end

function alg_cache(alg::RK46NL,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = RK46NLConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  RK46NLCache(u,uprev,k,tmp,fsalfirst,tab)
end

function alg_cache(alg::RK46NL,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  RK46NLConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function alg_cache(alg::Tsit5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = Tsit5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  k7 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  Tsit5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,utilde,tmp,atmp,tab)
end

alg_cache(alg::Tsit5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = Tsit5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct DP5Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  dense_tmp3::rateType
  dense_tmp4::rateType
  update::rateType
  bspl::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::DP5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = k2
  k7 = zero(rate_prototype) # This is FSAL'd to k1
  dense_tmp3 = k2
  dense_tmp4 = k5
  bspl = k3

  tmp = similar(u) # has to be separate for FSAL
  utilde = tmp

  if eltype(u) != uEltypeNoUnits || calck
    update = zero(rate_prototype)
    atmp = similar(u,uEltypeNoUnits)
  else
    update = k7
    atmp = k3
  end

  tab = DP5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  cache = DP5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,dense_tmp3,dense_tmp4,update,bspl,utilde,tmp,atmp,tab)
  cache
end

alg_cache(alg::DP5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = DP5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct DP5ThreadedCache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  dense_tmp3::rateType
  dense_tmp4::rateType
  update::rateType
  bspl::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::DP5Threaded,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  k7 = zero(rate_prototype)
  dense_tmp3 = zero(rate_prototype)
  dense_tmp4 = zero(rate_prototype)
  update = zero(rate_prototype)
  bspl = zero(rate_prototype)
  utilde = similar(u)
  tmp = similar(u); atmp = similar(u,uEltypeNoUnits)
  tab = DP5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  DP5ThreadedCache(u,uprev,k1,k2,k3,k4,k5,k6,k7,dense_tmp3,dense_tmp4,update,bspl,utilde,tmp,atmp,tab)
end

@cache struct Anas5Cache{uType,rateType,uNoUnitsType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k1::rateType
  k2::rateType
  k3::rateType
  k4::rateType
  k5::rateType
  k6::rateType
  k7::rateType
  utilde::uType
  tmp::uType
  atmp::uNoUnitsType
  tab::TabType
end

function alg_cache(alg::Anas5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tab = Anas5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))
  k1 = zero(rate_prototype)
  k2 = zero(rate_prototype)
  k3 = zero(rate_prototype)
  k4 = zero(rate_prototype)
  k5 = zero(rate_prototype)
  k6 = zero(rate_prototype)
  k7 = zero(rate_prototype)
  utilde = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)
  Anas5Cache(u,uprev,k1,k2,k3,k4,k5,k6,k7,utilde,tmp,atmp,tab)
end

alg_cache(alg::Anas5,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = Anas5ConstantCache(real(uBottomEltypeNoUnits),real(tTypeNoUnits))

@cache struct KYK2014DGSSPRK_3S2_Cache{uType,rateType,TabType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  fsalfirst::rateType
  tab::TabType
  #temporary values for Shu-Osher
  u_1::uType
  u_2::uType
  kk_1::rateType
  kk_2::rateType
end

struct KYK2014DGSSPRK_3S2_ConstantCache{T,T2} <: OrdinaryDiffEqConstantCache
  #These are not α and β for RK but for Shu-Osher
  #see top of page 317 in
  #Optimal Strong-Stability-Preserving Runge–Kutta Time Discretizations for
  #Discontinuous Garlekin Methods, Kubatko, Yaeger, Ketcheson 2014
  α_10::T
  α_20::T
  α_21::T
  α_30::T
  α_32::T
  β_10::T
  β_21::T
  β_30::T
  β_32::T
  #Shu-Osher is normally stated for autonomous systems, the times
  #are calculated by hand for this scheme
  c_1::T
  c_2::T

  function KYK2014DGSSPRK_3S2_ConstantCache(::Type{T}, ::Type{T2}) where {T,T2}
    α_10 = T(1.0)
    α_20 = T(0.087353119859156)
    α_21 = T(0.912646880140844)
    α_30 = T(0.344956917166841)
    α_32 = T(0.655043082833159)
    β_10 = T(0.528005024856522)
    β_21 = T(0.481882138633993)
    β_30 = T(0.022826837460491)
    β_32 = T(0.345866039233415)
    c_1 = β_10
    c_2 = α_21 * β_10 + β_21 # ==0.96376427726
    new{T,T2}(α_10, α_20, α_21, α_30, α_32, β_10, β_21, β_30, β_32, c_1, c_2)
  end
end

function alg_cache(alg::KYK2014DGSSPRK_3S2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  u_1 = similar(u)
  u_2 = similar(u)
  kk_1 = zero(rate_prototype)
  kk_2 = zero(rate_prototype)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  tab = KYK2014DGSSPRK_3S2_ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  KYK2014DGSSPRK_3S2_Cache(u,uprev,k,fsalfirst,tab, u_1, u_2, kk_1, kk_2)
end

function alg_cache(alg::KYK2014DGSSPRK_3S2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  KYK2014DGSSPRK_3S2_ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end


text:
using HttpServer
server = Server() do req, res
    "Goodbye, World!"
end
run(server, 8080)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#!/usr/bin/env julia
using Test, DynamicBoundspODEsIneq, DynamicBoundsBase

const DBB = DynamicBoundsBase

@testset "DifferentialInequality Relaxation" begin

    x0(p) = [34.0; 20.0; 0.0; 0.0; 16.0; 0.0]
    function f!(du, u, p, t)
        du[1] = -p[1]*u[1]*u[2] + p[2]*u[3] + p[6]*u[6]
        du[2] = -p[1]*u[1]*u[2] + p[2]*u[3] + p[3]*u[3]
        du[3] =  p[1]*u[1]*u[2] - p[2]*u[3] - p[3]*u[3]
        du[4] =  p[3]*u[3] - p[4]*u[4]*u[5] + p[5]*u[6]
        du[5] = -p[4]*u[4]*u[5] + p[5]*u[6] + p[6]*u[6]
        du[6] =  p[4]*u[4]*u[5] - p[5]*u[6] - p[6]*u[6]
        return
    end

    tspan = (0.0, 0.1)
    pL = [0.1; 0.033; 16.0; 5.0; 0.5; 0.3]
    pU = 10.0*pL

    prob = DynamicBoundsBase.ODERelaxProb(f!, tspan, x0, pL, pU)

    A = [0.0 -1.0 -1.0  0.0  0.0  0.0;
         0.0  0.0  0.0  0.0 -1.0 -1.0;
         1.0 -1.0 0.0 1.0 -1.0 0.0]
    b = [-20.0; -16.0; -2.0]
    set!(prob, PolyhedralConstraint(A, b))

    xL = zeros(6)
    xU = [34.0; 20.0; 20.0; 34.0; 16.0; 16.0]
    set!(prob, ConstantStateBounds(xL,xU))

    integrator = DifferentialInequality(prob, calculate_relax = true,
                                        calculate_subgradient = true)

    relax!(integrator)

    @test isapprox(integrator.relax_lo[6,77], 0.021504949016377915, atol = 1E-5)
    @test isapprox(integrator.relax_hi[6,77], 10.049461417090038, atol = 1E-3)
    @test isapprox(integrator.relax_cv[6,77], 0.02150494901637792, atol = 1E-5)
    @test isapprox(integrator.relax_cc[6,77], 5.173156379620435, atol = 1E-3)
    @test isapprox(integrator.relax_cv_grad[6,77][6], -0.00767, atol = 1E-5)
    @test isapprox(integrator.relax_cc_grad[6,77][5], 0.004196170550054147, atol = 1E-5)
end

@testset "Access functions" begin
    x0(p) = [34.0; 20.0; 0.0; 0.0; 16.0; 0.0]
    function f!(du, u, p, t)
        du[1] = -p[1]*u[1]*u[2] + p[2]*u[3] + p[6]*u[6]
        du[2] = -p[1]*u[1]*u[2] + p[2]*u[3] + p[3]*u[3]
        du[3] =  p[1]*u[1]*u[2] - p[2]*u[3] - p[3]*u[3]
        du[4] =  p[3]*u[3] - p[4]*u[4]*u[5] + p[5]*u[6]
        du[5] = -p[4]*u[4]*u[5] + p[5]*u[6] + p[6]*u[6]
        du[6] =  p[4]*u[4]*u[5] - p[5]*u[6] - p[6]*u[6]
        return
    end

    tspan = (0.0, 0.1)
    pL = [0.1; 0.033; 16.0; 5.0; 0.5; 0.3]
    pU = 10.0*pL

    prob = DynamicBoundsBase.ODERelaxProb(f!, tspan, x0, pL, pU)

    set!(prob, DBB.SupportSet([i for i in 0.0:0.005:0.1]))

    A = [0.0 -1.0 -1.0  0.0  0.0  0.0;
         0.0  0.0  0.0  0.0 -1.0 -1.0;
         1.0 -1.0 0.0 1.0 -1.0 0.0]
    b = [-20.0; -16.0; -2.0]
    set!(prob, PolyhedralConstraint(A, b))

    xL = zeros(6)
    xU = [34.0; 20.0; 20.0; 34.0; 16.0; 16.0]
    set!(prob, ConstantStateBounds(xL,xU))

    integrator = DifferentialInequality(prob, calculate_relax = true,
                                        calculate_subgradient = true,
                                        calculate_local_sensitivity = true)

    @test DBB.supports(integrator, DBB.IntegratorName())
    @test DBB.supports(integrator, DBB.Gradient())
    @test DBB.supports(integrator, DBB.Subgradient())
    @test DBB.supports(integrator, DBB.Bound{Lower}())
    @test DBB.supports(integrator, DBB.Bound{Upper}())
    @test DBB.supports(integrator, DBB.Relaxation{Lower}())
    @test DBB.supports(integrator, DBB.Relaxation{Upper}())
    @test DBB.supports(integrator, DBB.IsNumeric())
    @test DBB.supports(integrator, DBB.IsSolutionSet())
    @test DBB.supports(integrator, DBB.TerminationStatus())
    @test DBB.supports(integrator, DBB.Value())
    @test DBB.supports(integrator, DBB.ParameterValue())

    @test DBB.get(integrator, DBB.IntegratorName()) === "DifferentialInequality Integrator"
    @test !DBB.get(integrator, DBB.IsNumeric())
    @test DBB.get(integrator, DBB.IsSolutionSet())
    @test DBB.get(integrator, DBB.TerminationStatus()) === RELAXATION_NOT_CALLED

    DBB.relax!(integrator)
    DBB.integrate!(integrator)

    support_set = DBB.get(integrator, DBB.SupportSet())
    @test isapprox(support_set.s[2], 0.005)

    vout = zeros(6, size(integrator.local_problem_storage.x, 2))
    DBB.getall!(vout, integrator, Value())
    @test isapprox(vout[5,end], 1.1438878551545235, atol = 1E-5)
    @test isapprox(vout[6,end], 14.856112144845488, atol = 1E-5)

    out2 = Matrix{Float64}[]
    for i = 1:6
        push!(out2, zeros(6, size(integrator.local_problem_storage.dxdp[1], 2)))
    end
    DBB.getall!(out2, integrator, DBB.Gradient{Nominal}())
    @test isapprox(out2[2][3, 20], -0.007927203440720683, atol = 1E-5)

    out3 = Matrix{Float64}[]
    for i = 1:6
        push!(out3, zeros(6, size(integrator.local_problem_storage.dxdp[1], 2)))
    end
    @test_throws ErrorException DBB.getall!(out3, integrator, DBB.Gradient{Lower}())

    out4 = Matrix{Float64}[]
    for i = 1:6
        push!(out4, zeros(6, size(integrator.local_problem_storage.dxdp[1], 2)))
    end
    @test_throws ErrorException DBB.getall!(out4, integrator, DBB.Gradient{Upper}())

    out5 = Matrix{Float64}[]
    for i = 1:6
        push!(out5, zeros(6, size(integrator.local_problem_storage.dxdp[1], 2)))
    end
    DBB.getall!(out5, integrator, DBB.Subgradient{Lower}())
    @test isapprox(out5[2][3, 20], 0.0, atol = 1E-8)

    out6 = Matrix{Float64}[]
    for i = 1:6
        push!(out6, zeros(6, size(integrator.local_problem_storage.dxdp[1], 2)))
    end
    DBB.getall!(out6, integrator, DBB.Subgradient{Upper}())
    @test isapprox(out6[2][3, 20], 0.0, atol = 1E-8)

    out7 = zeros(6, size(integrator.relax_lo,2))
    DBB.getall!(out7, integrator, DBB.Bound{Lower}())
    @test isapprox(out7[6, 21], 0.09008717985737755, atol = 1E-5)

    out8 = zeros(6, size(integrator.relax_lo,2))
    DBB.getall!(out8, integrator, DBB.Bound{Upper}())
    @test isapprox(out8[6, 21], 1066.5815332523712, atol = 1E-5)

    out11 = zeros(6, size(integrator.relax_lo,2))
    DBB.getall!(out11, integrator, DBB.Relaxation{Lower}())
    @test isapprox(out11[6, 21], 0.09008717985737755, atol = 1E-5)

    out12 = zeros(6, size(integrator.relax_lo,2))
    DBB.getall!(out12, integrator, DBB.Relaxation{Upper}())
    @test isapprox(out12[6, 21], 976.3394543605548, atol = 1E-5)

    val1 = zeros(6) .- 0.1
    DBB.setall!(integrator, DBB.ParameterBound{Lower}(), val1)
    @test integrator.pL[1] == -0.1

    val2 = zeros(6) .+ 0.2
    DBB.setall!(integrator, DBB.ParameterBound{Upper}(), val2)
    @test integrator.pU[1] == 0.2

    val3 = zeros(6) .+ 0.1
    DBB.setall!(integrator, DBB.ParameterValue(), val3)
    @test integrator.p[1] == 0.1

    val01 = -0.3
    DBB.set!(integrator, DBB.ParameterBound{Lower}(2), val01)
    @test integrator.pL[2] == -0.3

    val02 = 0.5
    DBB.set!(integrator, DBB.ParameterBound{Upper}(2), val02)
    @test integrator.pU[2] == 0.5

    val03 = 0.2
    DBB.set!(integrator, DBB.ParameterValue(2), val03)
    @test integrator.p[2] == 0.2

    out15 = zeros(size(integrator.relax_lo,2))
    DBB.getall!(out15, integrator, DBB.Bound{Lower}())
    @test isapprox(out15[12], 11.452655899207926, atol = 1E-5)

    out16 = zeros(size(integrator.relax_lo,2))
    DBB.getall!(out16, integrator, DBB.Bound{Upper}())
    @test isapprox(out16[12], 32.8336928873473, atol = 1E-5)
end


text:
import Pkg
Pkg.activate("examples")
using Revise, GpABC, DifferentialEquations, Distances, Distributions

# seed!(4)

#
# ABC settings
#

n_particles = 1000
threshold_schedule = [3.0, 2.0, 1.0, 0.5, 0.2]
progress_every = 1000

#
# Emulation settings
#
n_design_points = 100

#
# True parameters
#
true_params =  [2.0, 2.0, 15.0, 1.0, 1.0, 1.0, 100.0, 1.0, 1.0, 1.0]
priors = [Uniform(0.2, 5.), Uniform(0.2, 5.), Uniform(10., 20.),
            Uniform(0.2, 2.), Uniform(0.2, 2.), Uniform(0.2, 2.),
            Uniform(75., 125.),
            Uniform(0.2, 2.), Uniform(0.2, 2.), Uniform(0.2, 2.)]
param_indices = [1,2,3]
n_var_params = length(param_indices)

#
# ODE solver settings
#
Tspan = (0.0, 10.0)
x0 = [0.0, 0.0, 0.0]
solver = RK4()
saveat = 0.1

function ODE_3GeneReg(dx, x, pars, t)
  dx[1] = pars[1]/(1+pars[7]*x[3]) - pars[4]*x[1]
  dx[2] = pars[2]*pars[8]*x[1]/(1+pars[8]*x[1]) - pars[5]*x[2]
  dx[3] = pars[3]*pars[9]*x[1]*pars[10]*x[2]/((1+pars[9]*x[1])*(1+pars[10]*x[2])) - pars[6]*x[3]
end

#
# Returns the solution to the toy model as solved by DifferentialEquations
#
function GeneReg(params::AbstractArray{Float64,1},
    Tspan::Tuple{Float64,Float64}, x0::AbstractArray{Float64,1},
    solver::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm, saveat::Float64)

  prob = ODEProblem(ODE_3GeneReg, x0 ,Tspan, params)
  Obs = solve(prob, solver, saveat=saveat)

  return Array{Float64, 2}(Obs)
end

reference_data = GeneReg(true_params, Tspan, x0, solver, saveat)
# reference_data += randn(size(reference_data)) * 0.1

function simulator_function(var_param_idx::Int, var_param::Real)
    params = copy(true_params)
    params[var_param_idx] = var_param
    GeneReg(params, Tspan, x0, solver, saveat)
end

function simulator_function(var_param_indices::AbstractArray{Int, 1}, var_params::AbstractArray{T, 1}) where {T<:Real}
    params = copy(true_params)
    params[var_param_indices] .= var_params
    GeneReg(params, Tspan, x0, solver, saveat)
end

function simulator_function(var_params::AbstractArray{T, 1}) where {T<:Real}
    simulator_function(param_indices, var_params)
    params = copy(true_params)
    params[param_indices] .= var_params
    GeneReg(params, Tspan, x0, solver, saveat)
end

function simulate_distance(var_params)
    distance_metric(summary_stats(reference_data),
        summary_stats(simulator_function(var_params)))
end

function simulate_distance(var_param_idx::Int, var_param::Real)
    distance_metric(summary_stats(reference_data),
        summary_stats(simulator_function(var_param_idx, var_param)))
end

function simulate_distance(var_param_indices::AbstractArray{Int, 1}, var_params::AbstractArray{T, 1}) where {T<:Real}
    distance_metric(summary_stats(reference_data),
        summary_stats(simulator_function(var_param_indices, var_params)))
end

println("EMULATION")
emu_out = EmulatedABCSMC(reference_data,
    simulator_function,
    priors[param_indices],
    threshold_schedule,
    n_particles,
    n_design_points;
    batch_size=1000,
    emulator_retraining = PreviousPopulationThresholdRetraining(n_design_points, min(100, n_design_points/2), 10),
    emulated_particle_selection = MeanVarEmulatedParticleSelection()
    )

println("SIMULATION")
sim_out = SimulatedABCSMC(reference_data,
    simulator_function,
    priors[param_indices],
    threshold_schedule,
    n_particles)


text:
using InplaceDHMC
using Test

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


text:
function createmediadict()
  Dict(
   :mud=>11.35,
   :water=>8.339,
   :air=>0.0,
   :heavymud=>12.0,
   :lightmud=>10.0
  )
end


text:
using ExecutableSpecifications
using ExecutableSpecifications: findstepdefinition, NonUniqueStepDefinition, StepDefinitionLocation, NoMatchingStepDefinition
using ExecutableSpecifications: FromMacroStepDefinitionMatcher, CompositeStepDefinitionMatcher
using ExecutableSpecifications.Gherkin
using ExecutableSpecifications.Gherkin: Given, When, Then

@testset "Step definitions" begin
    @testset "Find a step definition" begin
        @testset "Find a step definition; A matching given step; A step is found" begin
            given = ExecutableSpecifications.Gherkin.Given("some definition")
            stepdef_matcher = ExecutableSpecifications.FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications.@given

                @given "some definition" begin
                    x = 1
                end
            """)

            stepdefinition = ExecutableSpecifications.findstepdefinition(stepdef_matcher, given)
            @test isa(stepdefinition.definition, Function)
        end

        @testset "Find a step definition; A non-matching given step; No step definition found" begin
            given = ExecutableSpecifications.Gherkin.Given("some other definition")
            stepdef_matcher = ExecutableSpecifications.FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications.@given

                @given "some definition" begin
                    x = 1
                end
            """)

            @test_throws ExecutableSpecifications.NoMatchingStepDefinition ExecutableSpecifications.findstepdefinition(stepdef_matcher, given)
        end

        @testset "Find a step definition; A matching given step with another description; A step is found" begin
            given = ExecutableSpecifications.Gherkin.Given("some other definition")
            stepdef_matcher = ExecutableSpecifications.FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications.@given

                @given "some other definition" begin
                    x = 1
                end
            """)

            stepdefinition = ExecutableSpecifications.findstepdefinition(stepdef_matcher, given)
            @test isa(stepdefinition.definition, Function)
        end

        @testset "Find a step definition in another matcher; The other matcher has no matching step; No step is found" begin
            # This test ensures that step definitions are local to a single matcher, so that they aren't
            # kept globally.
            given = ExecutableSpecifications.Gherkin.Given("some definition")
            stepdef_matcher = ExecutableSpecifications.FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications.@given

                @given "some definition" begin
                    x = 1
                end
            """)

            # There is no step definitions here, so it should not find any matching definitions.
            empty_matcher = ExecutableSpecifications.FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications.@given
            """)

            @test_throws ExecutableSpecifications.NoMatchingStepDefinition ExecutableSpecifications.findstepdefinition(empty_matcher, given)
        end
    end

    @testset "Execute a step definition" begin
        @testset "Execute a step definition; Store an int in context; Context stores the value" begin
            given = ExecutableSpecifications.Gherkin.Given("some definition")
            stepdef_matcher = ExecutableSpecifications.FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications.@given

                @given "some definition" begin
                    context[:x] = 1
                end
            """)

            context = ExecutableSpecifications.StepDefinitionContext()
            stepdefinition = ExecutableSpecifications.findstepdefinition(stepdef_matcher, given)
            stepdefinition.definition(context)

            @test context[:x] == 1
        end

        @testset "Execute a step definition; Store a string in context; Context stores the value" begin
            given = ExecutableSpecifications.Gherkin.Given("some definition")
            stepdef_matcher = ExecutableSpecifications.FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications.@given

                @given "some definition" begin
                    context[:x] = "Some string"
                end
            """)

            context = ExecutableSpecifications.StepDefinitionContext()
            stepdefinition = ExecutableSpecifications.findstepdefinition(stepdef_matcher, given)
            stepdefinition.definition(context)

            @test context[:x] == "Some string"
        end

        @testset "Execute a step definition; Retrieve a value from the context; Context value is present" begin
            given = ExecutableSpecifications.Gherkin.Then(":x has value 1")
            stepdef_matcher = ExecutableSpecifications.FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications: @then, @expect

                @then ":x has value 1" begin
                    @expect context[:x] == 1
                end
            """)

            context = ExecutableSpecifications.StepDefinitionContext()
            context[:x] = 1
            stepdefinition = ExecutableSpecifications.findstepdefinition(stepdef_matcher, given)

            @test stepdefinition.definition(context) == ExecutableSpecifications.SuccessfulStepExecution()
        end

        @testset "Execute a step definition; An empty step definition; Success is returned" begin
            given = ExecutableSpecifications.Gherkin.Given("some definition")
            stepdef_matcher = ExecutableSpecifications.FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications.@given

                @given "some definition" begin

                end
            """)

            context = ExecutableSpecifications.StepDefinitionContext()
            stepdefinition = ExecutableSpecifications.findstepdefinition(stepdef_matcher, given)
            @test stepdefinition.definition(context) == ExecutableSpecifications.SuccessfulStepExecution()
        end

        @testset "Execute a step definition; An assert fails; StepFailed is returned" begin
            given = ExecutableSpecifications.Gherkin.Given("some definition")
            stepdef_matcher = ExecutableSpecifications.FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications: @given, @expect

                @given "some definition" begin
                    @expect 1 == 2
                end
            """)

            context = ExecutableSpecifications.StepDefinitionContext()
            stepdefinition = ExecutableSpecifications.findstepdefinition(stepdef_matcher, given)
            @test stepdefinition.definition(context) isa ExecutableSpecifications.StepFailed
        end

        @testset "Execute a step definition; An empty When step; Success is returned" begin
            when = When("some action")
            stepdef_matcher = ExecutableSpecifications.FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications: @when

                @when "some action" begin

                end
            """)

            context = ExecutableSpecifications.StepDefinitionContext()
            stepdefinition = ExecutableSpecifications.findstepdefinition(stepdef_matcher, when)
            @test stepdefinition.definition(context) == ExecutableSpecifications.SuccessfulStepExecution()
        end

        @testset "Execute a step definition; An empty Then step; Success is returned" begin
            then = Then("some postcondition")
            stepdef_matcher = ExecutableSpecifications.FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications: @then

                @then "some postcondition" begin

                end
            """)

            context = ExecutableSpecifications.StepDefinitionContext()
            stepdefinition = ExecutableSpecifications.findstepdefinition(stepdef_matcher, then)
            @test stepdefinition.definition(context) == ExecutableSpecifications.SuccessfulStepExecution()
        end

        @testset "Execute a step definition; Step throws an exception; The error is not caught" begin
            given = Given("some precondition")
            stepdef_matcher = ExecutableSpecifications.FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications: @given

                @given "some precondition" begin
                    throw(ErrorException("Some error"))
                end
            """)

            context = ExecutableSpecifications.StepDefinitionContext()
            stepdefinition = findstepdefinition(stepdef_matcher, given)
            @test_throws ErrorException stepdefinition.definition(context)
        end

        @testset "Execute a step definition; Call a method defined in the steps file; Method is in scope" begin
            when = ExecutableSpecifications.Gherkin.When("calling empty function foo")
            stepdef_matcher = ExecutableSpecifications.FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications.@when

                foo() = nothing

                @when "calling empty function foo" begin
                    foo()
                end
            """)

            context = ExecutableSpecifications.StepDefinitionContext()
            stepdefinition = ExecutableSpecifications.findstepdefinition(stepdef_matcher, when)
            @test stepdefinition.definition(context) == ExecutableSpecifications.SuccessfulStepExecution()
        end
    end

    @testset "Non-unique step definitions" begin
        @testset "Find a step definition; Two step definitions have the same description; NonUniqueStepDefinition is thrown" begin
            given = ExecutableSpecifications.Gherkin.Given("some definition")
            stepdef_matcher = ExecutableSpecifications.FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications.@given

                @given "some definition" begin
                end

                @given "some definition" begin
                end
            """)

            context = ExecutableSpecifications.StepDefinitionContext()
            @test_throws NonUniqueStepDefinition ExecutableSpecifications.findstepdefinition(stepdef_matcher, given)
        end

        @testset "Find a step definition; Two step definitions have the same description; File is reported for both" begin
            given = ExecutableSpecifications.Gherkin.Given("some definition")
            stepdef_matcher = ExecutableSpecifications.FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications.@given

                @given "some definition" begin
                end

                @given "some definition" begin
                end
            """; filename="steps.jl")

            context = ExecutableSpecifications.StepDefinitionContext()
            exception_thrown = false
            try
                ExecutableSpecifications.findstepdefinition(stepdef_matcher, given)
            catch ex
                if ex isa NonUniqueStepDefinition
                    exception_thrown = true
                    @test ex.locations[1].filename == "steps.jl"
                    @test ex.locations[2].filename == "steps.jl"
                else
                    rethrow()
                end
            end

            @assert exception_thrown "No NonUniqueStepDefinition exception was thrown!"
        end

        @testset "Find a step definition; Two step definitions have the same description; Another file is reported for both" begin
            given = ExecutableSpecifications.Gherkin.Given("some definition")
            stepdef_matcher = ExecutableSpecifications.FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications.@given

                @given "some definition" begin
                end

                @given "some definition" begin
                end
            """; filename="othersteps.jl")

            context = ExecutableSpecifications.StepDefinitionContext()
            exception_thrown = false
            try
                ExecutableSpecifications.findstepdefinition(stepdef_matcher, given)
            catch ex
                if ex isa NonUniqueStepDefinition
                    exception_thrown = true
                    @test ex.locations[1].filename == "othersteps.jl"
                    @test ex.locations[2].filename == "othersteps.jl"
                else
                    rethrow()
                end
            end

            @assert exception_thrown "No NonUniqueStepDefinition exception was thrown!"
        end
    end

    @testset "Composite matcher" begin
        @testset "Find a step definition from a composite; First matcher has the definition; Definition is found" begin
            given = Given("some precondition")
            matcher1 = FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications: @given

                @given "some precondition" begin

                end
            """)

            compositematcher = CompositeStepDefinitionMatcher(matcher1)

            stepdefinition = findstepdefinition(compositematcher, given)
            @test stepdefinition.definition isa Function
            @test stepdefinition.description == "some precondition"
        end

        @testset "Find a step definition from a composite; Second matcher has the definition; Definition is found" begin
            given = Given("some other precondition")
            matcher1 = FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications: @given
            """)
            matcher2 = FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications: @given

                @given "some other precondition" begin

                end
            """)

            compositematcher = CompositeStepDefinitionMatcher(matcher1, matcher2)

            stepdefinition = findstepdefinition(compositematcher, given)
            @test stepdefinition.definition isa Function
            @test stepdefinition.description == "some other precondition"
        end

        @testset "Find two step definitions from a composite; Both exist in a matcher; Definitions are found" begin
            given = Given("some other precondition")
            when = When("some action")
            matcher1 = FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications: @given

            """)
            matcher2 = FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications: @given

                @given "some other precondition" begin

                end
            """)
            matcher3 = FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications: @when

                @when "some action" begin

                end
            """)

            compositematcher = CompositeStepDefinitionMatcher(matcher1, matcher2, matcher3)

            stepdefinition = findstepdefinition(compositematcher, given)
            @test stepdefinition.definition isa Function
            @test stepdefinition.description == "some other precondition"

            stepdefinition = findstepdefinition(compositematcher, when)
            @test stepdefinition.description == "some action"
        end

        @testset "Find a step definition from a composite; Matching two definitions; Non unique step exception thrown" begin
            given = Given("some precondition")
            matcher1 = FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications: @given

                @given "some precondition" begin

                end
            """)
            matcher2 = FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications: @given

                @given "some precondition" begin

                end
            """)

            compositematcher = CompositeStepDefinitionMatcher(matcher1, matcher2)

            @test_throws NonUniqueStepDefinition findstepdefinition(compositematcher, given)
        end

        @testset "Find a step definition from a composite; Matching two definitions in one matcher; Non unique step exception thrown" begin
            given = Given("some precondition")
            matcher1 = FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications: @given

                @given "some precondition" begin

                end

                @given "some precondition" begin

                end
            """)
            matcher2 = FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications: @given
            """)

            compositematcher = CompositeStepDefinitionMatcher(matcher1, matcher2)

            @test_throws NonUniqueStepDefinition findstepdefinition(compositematcher, given)
        end

        @testset "Find a step definition from a composite; Matches in both matchers; Non unique locations indicates both matchers" begin
            given = Given("some precondition")
            matcher1 = FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications: @given

                @given "some precondition" begin

                end
            """; filename="matcher1.jl")
            matcher2 = FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications: @given

                @given "some precondition" begin

                end
            """; filename="matcher2.jl")

            compositematcher = CompositeStepDefinitionMatcher(matcher1, matcher2)

            exception_is_thrown = false
            try
                findstepdefinition(compositematcher, given)
            catch ex
                if ex isa NonUniqueStepDefinition
                    exception_is_thrown = true
                    location_filenames = [location.filename for location in ex.locations]
                    @test "matcher1.jl" in location_filenames
                    @test "matcher2.jl" in location_filenames
                else
                    rethrow()
                end
            end

            @test exception_is_thrown
        end

        @testset "Find a step definition from a composite; Two matchings in one and one matching in second matcher; Non unique locations indicates both matchers" begin
            given = Given("some precondition")
            matcher1 = FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications: @given

                @given "some precondition" begin

                end

                @given "some precondition" begin

                end
            """; filename="matcher1.jl")
            matcher2 = FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications: @given

                @given "some precondition" begin

                end
            """; filename="matcher2.jl")

            compositematcher = CompositeStepDefinitionMatcher(matcher1, matcher2)

            exception_is_thrown = false
            try
                findstepdefinition(compositematcher, given)
            catch ex
                if ex isa NonUniqueStepDefinition
                    exception_is_thrown = true
                    location_filenames = [location.filename for location in ex.locations]
                    @test "matcher1.jl" in location_filenames
                    @test "matcher2.jl" in location_filenames
                else
                    rethrow()
                end
            end

            @test exception_is_thrown
        end

        @testset "Find a step definition from a composite; No matches found; NoMatchingStepDefinition thrown" begin
            given = Given("some precondition")
            matcher1 = FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications: @given
            """)
            matcher2 = FromMacroStepDefinitionMatcher("""
                using ExecutableSpecifications: @given
            """)

            compositematcher = CompositeStepDefinitionMatcher(matcher1, matcher2)

            @test_throws NoMatchingStepDefinition findstepdefinition(compositematcher, given)
        end
    end
end


text:
# =========================================================================================
"""
    readhmo{T <: AbstractFloat}(
        stream::IOStream,
              ::Type{T}=Float64
    )

Reads a complete surface model from the given HMO file.

# Return type
`Model{T, Triangle{T}}`

# Alias

    readhmo{T}(fname::String, ::Type{T}=Float64)

Reads the model using a file name rather than a `IOStream` object.
"""
function readhmo(
        stream::IOStream,
              ::Type{T}=Float64
    ) where T <: AbstractFloat
    nodes = readhmo_nodes(stream, T)
    Model(nodes, readhmo_elements(stream, nodes), readhmo_charges(stream, T))
end

function readhmo(
        fname::String,
             ::Type{T}=Float64
    ) where T
    open(fh -> readhmo(fh, T), fname)
end

# =========================================================================================
"""
    readhmo_nodes{T <: AbstractFloat}(
        stream::IOStream,
              ::Type{T}=Float64
    )

Reads all nodes from the given HMO file.

# Return type
`Vector{Vector{T}}`
"""
function readhmo_nodes(
        stream::IOStream,
              ::Type{T}=Float64
    ) where T <: AbstractFloat
    nodes = Vector{T}[]
    seek(stream, "BEG_NODL_DATA")
    readline(stream) # skip first line
    for line in eachline(stream)
        line == "END_NODL_DATA" && break
        push!(nodes, [parse(T, a) for a in split(line)[2:end]])
    end
    nodes
end

# =========================================================================================
"""
    readhmo_elements{T <: AbstractFloat}(
        stream::IOStream,
        nodes ::Vector{Vector{T}}
    )

Reads all elements from the given HMO file.

# Return type
`Vector{Triangle{T}}`
"""
function readhmo_elements(
        stream::IOStream,
        nodes ::Vector{Vector{T}}
    ) where T <: AbstractFloat
    elements = Triangle{T}[]
    seek(stream, "BEG_ELEM_DATA")
    readline(stream) # skip first line
    for line in eachline(stream)
        line == "END_ELEM_DATA" && break
        push!(elements, Triangle([nodes[parse(Int, a)] for a in split(line)[4:end]]...))
    end
    elements
end

# =========================================================================================
"""
    readhmo_charges{T <: AbstractFloat}(
        stream::IOStream,
              ::Type{T}=Float64
    )

Reads all charges from the given HMO file.

# Return type
`Vector{Charge{T}}`
"""
function readhmo_charges(
        stream::IOStream,
              ::Type{T}=Float64
    ) where T <: AbstractFloat
    charges = Charge{T}[]
    seek(stream, "BEG_CHARGE_DATA")
    readline(stream) # skip first line
    for line in eachline(stream)
        line == "END_CHARGE_DATA" && break
        push!(charges, Charge([parse(T, a) for a in split(line)[2:end]]...))
    end
    charges
end

# =========================================================================================
"""
    writehmo{T}(
        stream::IOStream,
        model ::Model{T, Triangle{T}}
    )

Creates a HMO file from the given surface model.

# Return type
`Void`

# Alias

    writehmo{T}(
        fname::String,
        model::Model{T, Triangle{T}}
    )

Creates the HMO file by name rather than `IOStream` object.
"""
function writehmo(
        stream::IOStream,
        model ::Model{T, Triangle{T}}
    ) where T
    write(stream, "BEG_NODL_DATA\n")
    write(stream, "$(length(model.nodes))\n")
    for (i, node) in enumerate(model.nodes)
        join(stream, Any[i, node...], "\t")
        write(stream, "\n")
    end
    write(stream, "END_NODL_DATA\n\n")

    getidx = v -> findfirst(isequal(v), model.nodes)

    write(stream, "BEG_ELEM_DATA\n")
    join(stream, [length(model.elements), zeros(Int, 12)...], "\t")
    write(stream, "\n")
    for (i, elem) in enumerate(model.elements)
        join(stream, Any[i, 1, 103, getidx(elem.v1), getidx(elem.v2), getidx(elem.v3)], "\t")
        write(stream, "\n")
    end
    write(stream, "END_ELEM_DATA\n\n")

    write(stream, "BEG_CHARGE_DATA\n")
    write(stream, "$(length(model.charges))\n")
    for (i, charge) in enumerate(model.charges)
        join(stream, Any[i, charge.pos..., charge.val], "\t")
        write(stream, "\n")
    end
    write(stream, "END_CHARGE_DATA\n")
end

function writehmo(fname::String, model::Model{T, Triangle{T}}) where T
    open(fh -> writehmo(fh, model), fname, "w")
end