text:
 using RecursiveArrayTools, Parameters # for bifurcation point handling in ContRes
import Base: show, length		# simplified display method for ContRes

"""
	options = ContinuationPar(dsmin = 1e-4,...)

Returns a variable containing parameters to affect the `continuation` algorithm when solving `F(x,p) = 0`.

# Arguments
- `dsmin, dsmax` are the minimum, maximum arclength allowed value. It controls the density of points in the computed branch of solutions.
- `ds` is the initial arclength.
- `theta` is a parameter in the arclength constraint. It is very **important** to tune it. See the docs of [`continuation`](@ref).
- `pMin, pMax` allowed parameter range for `p`
- `maxSteps` maximum number of continuation steps
- `newtonOptions::NewtonPar`: options for the Newton algorithm
- `saveToFile = false`: save to file. A name is automatically generated.
- `saveSolEveryNsteps::Int64 = 0` at which continuation steps do we save the current solution`
- `plotEveryNsteps = 3`

Handling eigen elements
- `computeEigenValues = false` compute eigenvalues / eigenvectors
- `nev = 3` number of eigenvalues to be computed. It is automatically increased to have at least `nev` unstable eigenvalues. To be set for proper  bifurcation detection. See [Detection of bifurcation points](@ref) for more informations.
- `saveEigEveryNsteps = 1`	record eigen vectors every specified steps. **Important** for memory limited ressource, *e.g.* GPU.
- `saveEigenvectors	= true`	**Important** for memory limited ressource, *e.g.* GPU.

Handling bifurcation detection
- `precisionStability = 1e-10` lower bound on the real part of the eigenvalues to test for stability of equilibria and periodic orbits
- `detectFold = true` detect Fold bifurcations? It is a useful option although the detection of Fold is cheap. Indeed, it may happen that there is a lot of Fold points and this can saturate the memory in memory limited devices (e.g. on GPU)
- `detectBifurcation::Int` ∈ {0, 1, 2} detect bifurcations? If set to 0, bifurcation are not detected. If set to 1, bifurcation are detected along the continuation run, but not located precisely. If set to 2, a bisection algorithm is used to locate the bifurcations (slower). The possibility to switch off detection is a useful option. Indeed, it may happen that there are a lot of bifurcation points and this can saturate the memory in memory limited devices (e.g. on GPU)
- `dsminBisection` dsmin for the bisection algorithm when locating bifurcation points
- `nInversion` number of sign inversions in bisection algorithm
- `maxBisectionSteps` maximum number of bisection steps
- `tolBisectionEigenvalue` tolerance on real part of eigenvalue to detect bifurcation points in the bisection steps

Handling `ds` adaptation (see [`continuation`](@ref) for more information)
- `a  = 0.5` aggressiveness factor. It is used to adapt `ds` in order to have a number of newton iteration per continuation step roughly the same. The higher `a` is, the larger the step size `ds` is changed at each continuation step.
- `thetaMin = 1.0e-3` minimum value of `theta`
- `doArcLengthScaling` trigger further adaptation of `theta`

 Misc
- `finDiffEps::T  = 1e-9` ε used in finite differences computations

!!! tip "Mutating"
    For performance reasons, we decided to use an immutable structure to hold the parameters. One can use the package `Setfield.jl` to drastically simplify the mutation of different fields. See tutorials for more examples.
"""
@with_kw struct ContinuationPar{T, S <: AbstractLinearSolver, E <: AbstractEigenSolver}
	# parameters for arclength continuation
	dsmin::T	= 1e-3
	dsmax::T	= 2e-2;		@assert dsmax >= dsmin
	ds::T		= 1e-3;		@assert dsmax >= abs(ds);	@assert abs(ds) >= dsmin
	@assert dsmin > 0
	@assert dsmax > 0

	# parameters for scaling arclength step size
	theta::T					= 0.5 	# parameter in the dot product used for the extended system
	doArcLengthScaling::Bool  	= false
	gGoal::T					= 0.5
	gMax::T						= 0.8
	thetaMin::T					= 1.0e-3
	a::T						= 0.5  	# aggressiveness factor
	tangentFactorExponent::T 	= 1.5

	# parameters bound
	pMin::T	= -1.0
	pMax::T	=  1.0; 			@assert pMax >= pMin

	# maximum number of continuation steps
	maxSteps::Int64  = 100

	# Newton solver parameters
	finDiffEps::T  = 1e-9 				#constant for finite differences
	newtonOptions::NewtonPar{T, S, E} = NewtonPar()

	saveToFile::Bool = false 			# save to file?
	saveSolEveryNsteps::Int64 = 0		# what steps do we save the current solution

	# parameters for eigenvalues
 	computeEigenValues::Bool = false
	nev::Int64 = 3 						# number of eigenvalues
	saveEigEveryNsteps::Int64 = 1		# what steps do we keep the eigenvectors
	saveEigenvectors::Bool	= true		# useful options because if puts a high memory pressure

	plotEveryNsteps::Int64 = 10

	# handling bifucation points
	precisionStability::T = 1e-10		# lower bound for stability of equilibria and periodic orbits
	detectFold::Bool = true 			# detect fold points?
	detectBifurcation::Int64 = 0		# detect bifurcation points?
	dsminBisection::T = 1e-5			# dsmin for the bisection algorithm when locating bifurcation points
	nInversion::Int64 = 2				# number of sign inversions in bisection algorithm
	maxBisectionSteps::Int64 = 15		# maximum number of bisection steps
	tolBisectionEigenvalue::Float64 = 1e-5 # tolerance on real part of eigenvalue to detect bifurcation points in the bisection steps
	@assert iseven(nInversion) "The option `nInversion` number must be even"
	@assert detectBifurcation < 3 "The option `detectBifurcation` must belong to {0,1,2}"
end

# check the logic of the parameters
function check(contParams::ContinuationPar)
	if contParams.detectBifurcation > 0
		@set contParams.computeEigenValues = true
	else
		contParams
	end
end
####################################################################################################
# Structure to hold result
"""
$(TYPEDEF)

Structure which holds the results after a call to [`continuation`](@ref).

$(TYPEDFIELDS)

          A vector holding the set of bifurcation points detected during the computation of the branch. Each entry of the vector contains a tuple with fields:
-
    - `type` bifurcation type, `:hopf, :bp...`,
    - `idx` is the index in `eig` (see above) for which the bifurcation occurs.
    - `param` parameter value at the bifurcation point
    - `norm` norm of the equilibrium at the bifurcation point
    - `printsol = printSolution(x, param)`
    - `x` equilibrium at the bifurcation point
    - `tau` tangent along the branch at the bifurcation point
    - `ind_bif` is the eigenvalue index responsible for the bifurcation (if applicable)
    - `step` is the continuation step at which the bifurcation occurs,
    - `status ∈ {:converged, :guess}` indicates if the bisection algorithm was successful in detecting the bifurcation point
    - `δ = (δr, δi)` where δr indicates the change in the number of unstable eigenvalues and δi indicates the change in the number of unstable eigenvalues with nonzero imaginary part. `abs(δr)` is thus an estimate of the dimension of the kernel of the Jacobian at the bifurcation point.
"""
@with_kw_noshow struct ContResult{T, Teigvals, Teigvec, Biftype, Ts, Tfunc, Tpar, Tl <: Lens}
	"holds the low-dimensional information about the branch. More precisely, `branch[:,i]` contains the following information `(param, printSolution(u, param), Newton iterations, ds, i)` for each continuation step `i`."
	branch::VectorOfArray{T, 2, Array{Vector{T}, 1}}

	"A vector with eigen-elements at each continuation step."
	eig::Vector{NamedTuple{(:eigenvals, :eigenvec, :step), Tuple{Teigvals, Teigvec, Int64}}}

	"A vector holding the set of fold points detected during the computation of the branch."
	foldpoint::Vector{Biftype}

	"A `Vector{Bool}` holding the stability of the computed solution for each continuation step"
	stability::Vector{Bool}

	"A `Vector{Int64}` holding the number of eigenvalues with positive real part and non zero imaginary part for each continuation step (to detect Hopf bifurcation)"
	n_imag::Vector{Int64}

	"A `Vector{Int64}` holding the number of eigenvalues with positive real part for each continuation step (to detect stationary bifurcation)"
	n_unstable::Vector{Int64}

	"Vector of solutions sampled along the branch. This is set by the argument `saveSolEveryNsteps::Int64` (default 0) in [`ContinuationPar`](@ref)"
	sol::Ts

	"The parameters used for the call to `continuation` which produced this branch."
	contparams::ContinuationPar

	"Type of solutions computed in this branch."
	type::Symbol = :Equilibrium

	"Structure associated to the functional, useful for branch switching. For example, when computing periodic orbits, the functional `PeriodicOrbitTrapProblem`, `ShootingProblem`... will be saved here."
	functional::Tfunc = nothing

	"Parameters passed to continuation and used in the equation F(x, par) = 0"
	params::Tpar = nothing

	"Parameter axis used for computing the branch"
	param_lens::Tl

	# the explanation for this field is given after. Always put this field in last position
	""
	bifpoint::Vector{Biftype}

end

length(br::ContResult) = length(br.branch[1, :])

_show(io, bp, ii) = @printf(io, "- %3i, %7s point around p ≈ %4.8f, step = %3i, idx = %3i, ind_bif = %3i [%9s], δ = (%2i, %2i)\n", ii, bp.type, bp.param, bp.step, bp.idx, bp.ind_bif, bp.status, bp.δ[1], bp.δ[2])

function show(io::IO, br::ContResult)
	println(io, "Branch number of points: ", length(br.branch))
	println(io, "Branch of ", br.type)
	if length(br.bifpoint) > 0
		println(io, "Bifurcation points:")
		for ii in eachindex(br.bifpoint)
			_show(io, br.bifpoint[ii], ii)
		end
	end
	if length(br.foldpoint) > 0
		println(io, "Fold points:")
		for ii in eachindex(br.foldpoint)
			_show(io, br.foldpoint[ii], ii)
		end
	end
end

"""
This function is used to initialize the composite type `ContResult` according to the options contained in `contParams`
"""
function initContRes(br, x0, par, lens::Lens, evsol, contParams::ContinuationPar{T, S, E}) where {T, S, E}
	bif0 = (type = :none, idx = 1, param = T(0), norm  = T(0), printsol = T(0), x = x0, tau = x0, ind_bif = 0, step = 0, status = :guess, δ = (0,0))
	sol = contParams.saveSolEveryNsteps > 0 ? [(x = copy(x0), p = br[1,1], step = 0)] : nothing
	n_unstable = 0
	n_imag = 0
	stability = true

	if contParams.computeEigenValues
		evvectors = contParams.saveEigenvectors ? evsol[2] : nothing
		stability, n_unstable, n_imag = is_stable(contParams, evsol[1])
		_evvectors = (eigenvals = evsol[1], eigenvec = evvectors, step = 0)
	else
		_evvectors = (eigenvals = Complex{T}(1), eigenvec = nothing, step = 0)
	end

	return ContResult(
		branch = br,
		bifpoint = [bif0],
		foldpoint = [bif0],
		n_imag = [n_imag],
		n_unstable = [n_unstable],
		stability = [stability],
		eig = [_evvectors],
		sol = sol,
		contparams =  contParams,
		params = par,
		param_lens = lens)
end
####################################################################################################
# Iterator interface

@with_kw struct PALCIterable{TF, TJ, Tv, Tp, Tlens, T, S, E, Ttangent, Tlinear, Tplotsolution, Tprintsolution, TnormC, Tdot, Tfinalisesolution, Tcallback, Tfilename}
	F::TF
	J::TJ

	x0::Tv							# initial guess
	par::Tp							# reference to parameter, so no worry if this one is big
	param_lens::Tlens				# param axis to be considered specified by a ::Lens

	contParams::ContinuationPar{T, S, E}

	tangentAlgo::Ttangent
	linearAlgo::Tlinear

	plot::Bool = false
	plotSolution::Tplotsolution
	printSolution::Tprintsolution

	normC::TnormC
	dottheta::Tdot
	finaliseSolution::Tfinalisesolution
	callbackN::Tcallback

	verbosity::Int64 = 2

	filename::Tfilename
end

import Base: eltype

eltype(it::PALCIterable{TF, TJ, Tv, Tp, Tlens, T, S, E, Ttangent, Tlinear, Tplotsolution, Tprintsolution, TnormC, Tdot, Tfinalisesolution, Tcallback, Tfilename}) where {TF, TJ, Tv, Tp, Tlens, T, S, E, Ttangent, Tlinear, Tplotsolution, Tprintsolution, TnormC, Tdot, Tfinalisesolution, Tcallback, Tfilename} = T

function PALCIterable(Fhandle, Jhandle,
					x0, par, lens::Lens,
					contParams::ContinuationPar{T, S, E},
					linearAlgo::AbstractBorderedLinearSolver = BorderingBLS()
					;
					filename = "branch-" * string(Dates.now()),
					tangentAlgo = SecantPred(),
					plot = false,
					plotSolution = (x, p; kwargs...) -> nothing,
					printSolution = (x, p) -> norm(x),
					normC = norm,
					dotPALC = (x,y) -> dot(x,y) / length(x),
					finaliseSolution = (z, tau, step, contResult) -> true,
					callbackN = (x, f, J, res, iteration, itlinear, optionsN; kwargs...) -> true,
					verbosity = 0
					) where {T <: Real, S, E}

	return PALCIterable(F = Fhandle, J = Jhandle, x0 = x0, par = par, param_lens = lens, contParams = check(contParams), tangentAlgo = tangentAlgo, linearAlgo = linearAlgo, plot = plot, plotSolution = plotSolution, printSolution = printSolution, normC = normC, dottheta = DotTheta(dotPALC), finaliseSolution = finaliseSolution, callbackN = callbackN, verbosity = verbosity, filename = filename)
end

"""
	state = PALCStateVariables(ds = 1e-4,...)

Returns a variable containing the state of the continuation procedure.

# Arguments
- `z_pred` current solution on the branch
- `tau` tangent predictor
- `z_old` previous solution
- `isconverged` Boolean for newton correction
- `itnewton` Number of newton iteration (in corrector)
- `step` current continuation step
- `ds` step size
- `theta` theta parameter for constraint equation in PALC
- `stopcontinuation` Boolean to stop continuation

# Useful functions
- `copy(state)` returns a copy of `state`
- `solution(state)` returns the current solution (x, p)
- `getx(state)` returns the x component of the current solution
- `getp(state)` returns the p component of the current solution
"""
@with_kw mutable struct PALCStateVariables{Tv, T, Teigvals, Teigvec}
	z_pred::Tv								# current solution
	tau::Tv									# tangent predictor

	z_old::Tv								# previous solution

	isconverged::Bool						# Boolean for newton correction
	itnewton::Int64							# Number of newton iteration (in corrector)

	step::Int64 = 0							# current continuation step
	ds::T									# step size
	theta::T								# theta parameter for constraint equation in PALC

	stopcontinuation::Bool = false			# Boolean to stop continuation
	stepsizecontrol::Bool = true			# Perform step size adaptation

	n_unstable::Tuple{Int64,Int64}  = (-1, -1)	# (current, previous)
	n_imag::Tuple{Int64,Int64} 		= (-1, -1)	# (current, previous)

	eigvals::Teigvals = nothing				# current eigenvalues
	eigvecs::Teigvec = nothing				# current eigenvectors
end

import Base: copy

function copy(state::PALCStateVariables)
	return PALCStateVariables(
		z_pred 	= _copy(state.z_pred),
		tau = _copy(state.tau),
		z_old 	= _copy(state.z_old),
		isconverged = state.isconverged,
		itnewton 	= state.itnewton,
		step 		= state.step,
		ds 			= state.ds,
		theta 		= state.theta,
		stopcontinuation = state.stopcontinuation,
		stepsizecontrol  = state.stepsizecontrol,
		n_unstable 		 = state.n_unstable,
		n_imag 			 = state.n_imag
	)
end

solution(state::PALCStateVariables) = state.z_old
getx(state::PALCStateVariables) = state.z_old.u
@inline getp(state::PALCStateVariables) = state.z_old.p
isstable(state::PALCStateVariables) = state.n_unstable[1] == 0

# condition for halting the continuation procedure
@inline done(it::PALCIterable, state::PALCStateVariables) =
			(state.step <= it.contParams.maxSteps) &&
			(it.contParams.pMin <= state.z_old.p <= it.contParams.pMax) &&
			(state.stopcontinuation == false)

function updatestability!(state::PALCStateVariables, n_unstable, n_imag)
	state.n_unstable = (n_unstable, state.n_unstable[1])
	state.n_imag = (n_imag, state.n_imag[1])
end

# we detect a bifurcation by a change in the number of unstable eigenvalues
function detectBifucation(state::PALCStateVariables)
	n1, n2 = state.n_unstable
	# deals with missing value encoded by n_unstable = -1
	if n1 == -1 || n2 == -1; return false; end
	return n1 !== n2
end

function save!(contres::ContResult, it::PALCIterable, state::PALCStateVariables)
	push!(contres.branch, getStateSummary(it, state))

	if state.n_unstable[1] >= 0 # whether to deal with n_unstable = -1
		push!(contres.n_unstable, state.n_unstable[1])
		push!(contres.stability, isstable(state))
	end

	# if condition to deal with n_imag = -1
	if state.n_imag[1] >= 0; push!(contres.n_imag, state.n_imag[1]); end

	# save solution
	if it.contParams.saveSolEveryNsteps > 0 &&
			mod(state.step, it.contParams.saveSolEveryNsteps) == 0
		push!(contres.sol, (x = copy(getx(state)), p = getp(state), step = state.step))
	end
	# save eigen elements
	if it.contParams.computeEigenValues
		if mod(state.step, it.contParams.saveEigEveryNsteps) == 0
			if it.contParams.saveEigenvectors
				push!(contres.eig, (eigenvals = state.eigvals, eigenvec = state.eigvecs, step = state.step))
			else
				push!(contres.eig, (eigenvals = state.eigvals, eigenvec = state.eigvecs, step = state.step))
			end
		end
	end
end

function initContRes(it::PALCIterable, state::PALCStateVariables)
	x0 = getx(state)
	p0 = getp(state)
	contParams = it.contParams

	if contParams.computeEigenValues
		eiginfo = contParams.newtonOptions.eigsolver(it.J(x0, set(it.par, it.param_lens, p0)), contParams.nev)
		_, n_unstable, n_imag = is_stable(contParams, eiginfo[1])
		updatestability!(state, n_unstable, n_imag)
		return initContRes(VectorOfArray([getStateSummary(it, state)]), x0, it.par, it.param_lens, eiginfo, contParams)
	else
		T = eltype(it)
		eiginfo = (Complex{T}(1), nothing, false, 0)
		return initContRes(VectorOfArray([getStateSummary(it, state)]), x0, it.par, it.param_lens, nothing, contParams)
	end
end

import Base: iterate

function iterate(it::PALCIterable; _verbosity = it.verbosity)
	# this is to overwrite verbosity behaviour, like when locating bifurcations
	verbosity = min(it.verbosity, _verbosity)
	pas0 = it.par
	p0 = get(it.par, it.param_lens)
	ds = it.contParams.ds
	theta = it.contParams.theta
	T = eltype(it)

	(verbosity > 0) && printstyled("#"^53*"\n********** Pseudo-Arclength Continuation ************\n\n", bold = true, color = :red)

	# Get parameters
	@unpack pMin, pMax, maxSteps, newtonOptions = it.contParams
	epsi = it.contParams.finDiffEps

	# Converge initial guess
	(verbosity > 0) && printstyled("*********** CONVERGE INITIAL GUESS *************", bold = true, color = :magenta)
	# we pass additional kwargs to newton so that it is sent to the newton callback
	u0, fval, isconverged, itnewton = newton(it.F, it.J, it.x0, it.par, newtonOptions; normN = it.normC, callback = it.callbackN, iterationC = 0, p = p0)
	@assert isconverged "Newton failed to converge initial guess"
	(verbosity > 0) && println("\n--> convergence of initial guess = OK")
	(verbosity > 0) && println("--> parameter = $(p0), initial step")
	(verbosity > 0) && printstyled("\n******* COMPUTING INITIAL TANGENT *************", bold = true, color = :magenta)
	η = T(150)
	u_pred, fval, isconverged, itnewton = newton(it.F, it.J,
			u0, set(it.par, it.param_lens, p0 + ds / η), newtonOptions; normN = it.normC, callback = it.callbackN, iterationC = 0, p = p0 + ds / η)
	@assert isconverged "Newton failed to converge for the computation of the initial tangent"
	(verbosity > 0) && (print("\n--> convergence of initial guess = ");printstyled("OK\n\n", color=:green))
	(verbosity > 0) && println("--> parameter = $(p0 + ds/η), initial step (bis)")
	z_old   = BorderedArray(copyto!(similar(u0), u0), p0)
	z_pred	= BorderedArray(copyto!(similar(u_pred), u_pred), p0 + ds / η)
	tau  = copy(z_pred)

	# compute the tangents
	getTangent!(tau, z_pred, z_old, it, ds, theta, SecantPred(), _verbosity)

	# compute eigenvalues to get the type
	if it.contParams.computeEigenValues
		eigvals, eigvecs, _, _ = it.contParams.newtonOptions.eigsolver(it.J(u0, it.par), it.contParams.nev)
		if it.contParams.saveEigenvectors == false
			eigvecs = nothing
		end
	else
		eigvals = nothing
		eigvecs = nothing
	end

	# return the state
	state = PALCStateVariables(z_pred = z_pred, tau = tau, z_old = z_old, isconverged = true, ds = it.contParams.ds, theta = it.contParams.theta, itnewton = 0, eigvals = eigvals, eigvecs = eigvecs)	# previous tangent
	return state, state
end

function iterate(it::PALCIterable, state::PALCStateVariables; _verbosity = it.verbosity)
	if !done(it, state) return nothing end
	# this is to overwrite verbosity behaviour, like when locating bifurcations
	verbosity = min(it.verbosity, _verbosity)

	@unpack step, ds, theta = state

	# Predictor: z_pred, following method only mutates z_pred
	getPredictor!(state.z_pred, state.z_old, state.tau, ds, it.tangentAlgo)
	(verbosity > 0) && println("#"^35)
	(verbosity > 0) && @printf("Start of Continuation Step %d:\nParameter = %2.4e ⟶  %2.4e [guess]\n", step, state.z_old.p, state.z_pred.p)
	(verbosity > 0) && @printf("Step size = %2.4e\n", ds)

	# Corrector, ie newton correction. This does not mutate the arguments
	z_newton, fval, state.isconverged, state.itnewton  = corrector(it,
			state.z_old, state.tau, state.z_pred,
			ds, theta,
			it.tangentAlgo, it.linearAlgo,
			normC = it.normC, callback = it.callbackN, iterationC = step, p = state.z_old.p)

	# Successful step
	if state.isconverged
		(verbosity > 0) && printstyled("--> Step Converged in $(state.itnewton) Nonlinear Iterations\n", color=:green)

		# Get predictor, it only mutates tau
		getTangent!(state.tau, z_newton, state.z_old, it,
					ds, theta, it.tangentAlgo, verbosity)

		# update current solution
		copyto!(state.z_old, z_newton)
	else
		(verbosity > 0) && printstyled("Newton correction failed\n", color=:red)
		(verbosity > 0) && println("--> Newton Residuals history = ", fval)
	end

	if state.stopcontinuation == false && state.stepsizecontrol == true
		# we update the PALC paramters ds and theta, they are in the state variable
		state.ds, state.theta, state.stopcontinuation = stepSizeControl(ds, theta, it.contParams, state.isconverged, state.itnewton, state.tau, verbosity)
	end

	state.step += 1
	return state, state
end

function getStateSummary(it, state)
	x = getx(state)
	p = getp(state)
	vcat(p, it.printSolution(x, p), state.itnewton, state.ds, state.theta, state.step)
end

function continuation!(it::PALCIterable, state::PALCStateVariables, contRes::ContResult)
	contParams = it.contParams
	verbosity = it.verbosity

	next = (state, state)

	while next !== nothing
		# we get the current state
		(i, state) = next
		########################################################################################
		# the new solution has been successfully computed
		# we perform saving, plotting, computation of eigenvalues...
		# the case state.step = 0 was just done above
		if state.isconverged && (state.step <= it.contParams.maxSteps) && (state.step > 0)

			# Eigen-elements computation, they are stored on state
			if contParams.computeEigenValues
				itnewton = computeEigenvalues!(it, state)

				(verbosity > 0) && printstyled(color=:green,"--> Computed ", length(state.eigvals), " eigenvalues in ", itnewton, " iterations, #unstable = ", state.n_unstable[1],"\n")
			end

			# Detection of fold points based on parameter monotony, mutates contRes.foldpoint
			if contParams.detectFold; locateFold!(contRes, it, state); end

			if contParams.detectBifurcation > 0 && detectBifucation(state)
				status::Symbol = :guess
				if contParams.detectBifurcation > 1
					(verbosity > 0) && printstyled(color=:red, "--> Bifurcation detected before p = ", getp(state), "\n")

					# locate bifurcations, mutates state so that it stays very close to the bifurcation point. It also updates the eigenelements at the current state
					status = locateBifurcation!(it, state, verbosity > 2)
				end
				if contParams.detectBifurcation>0 && detectBifucation(state)
					_, bifpt = getBifurcationType(contParams, state, it.normC, it.printSolution, verbosity, status)
					if bifpt.type != :none; push!(contRes.bifpoint, bifpt); end
				end
			end

			# Plotting
			(it.plot && mod(state.step, contParams.plotEveryNsteps) == 0 ) && plotBranchCont(contRes, state.z_old, contParams, it.plotSolution)

			# Saving Solution to File
			if contParams.saveToFile
				saveToFile(it.filename, getx(state), getp(state), state.step, contRes, contParams)
			end

			# Call user saved finaliseSolution function. If returns false, stop continuation
			if it.finaliseSolution(state.z_old, state.tau, state.step, contRes) == false
				state.stopcontinuation = true
			end

			# Save solution
			save!(contRes, it, state)
		end
		########################################################################################
		# body
		next = iterate(it, state)
	end

	# We remove the initial guesses which are meaningless
	popfirst!(contRes.bifpoint)
	popfirst!(contRes.foldpoint)

	# we remove the first element of branch, it was just to initialize it
	# popfirst!(contRes.branch.u)

	it.plot && plotBranchCont(contRes, state.z_old, contParams, it.plotSolution)

	# return current solution in case the corrector did not converge
	return contRes, state.z_old, state.tau
end

function continuation(it::PALCIterable)
	## !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	# contres is not known at compile time so we
	# need a function barrier to resolve its type
	#############################################

	# we compute the cache pour the continuation, i.e. state::PALCStateVariables
	state = iterate(it)[1]

	# variable to hold results from continuation
	contRes = initContRes(it, state)

	# perform the continuation
	return continuation!(it, state, contRes)
end

function continuation(Fhandle, Jhandle,
					x0, par, lens::Lens,
					contParams::ContinuationPar{T, S, E},
					linearAlgo::AbstractBorderedLinearSolver;
					tangentAlgo = SecantPred(),
					plot = false,
					plotSolution = (x, p; kwargs...) -> nothing,
					printSolution = (x, p) -> norm(x),
					normC = norm,
					dotPALC = (x,y) -> dot(x,y) / length(x),
					finaliseSolution = (z, tau, step, contResult) -> true,
					callbackN = (x, f, J, res, iteration, itlinear, optionsN; kwargs...) -> true,
					filename = "branch-" * string(Dates.now()),
					verbosity = 0) where {T <: Real, S, E}

	it = PALCIterable(Fhandle, Jhandle, x0, par, lens::Lens, contParams, linearAlgo;
						tangentAlgo = tangentAlgo, plot = plot, plotSolution = plotSolution, printSolution = printSolution, normC = normC, dotPALC = dotPALC, finaliseSolution = finaliseSolution, callbackN = callbackN, verbosity = verbosity, filename = filename)
	return continuation(it)
end

####################################################################################################

"""
	continuation(F, J, x0, par, lens::Lens, contParams::ContinuationPar; plot = false, normC = norm, dotPALC = (x,y) -> dot(x,y) / length(x), printSolution = norm, plotSolution = (x, p; kwargs...)->nothing, finaliseSolution = (z, tau, step, contResult) -> true, callbackN = (x, f, J, res, iteration, itlinear, options; kwargs...) -> true, linearAlgo = BorderingBLS(), tangentAlgo = SecantPred(), verbosity = 0)

Compute the continuation curve associated to the functional `F` and its jacobian `J`.

# Arguments:
- `F = (x, p) -> F(x, p)` where `p` is the set of parameters passed to `F`.
- `J = (x, p) -> d_xF(x, p)` its associated jacobian. It can be a matrix, a function or a callable struct.
- `x0` initial guess
- `par` initial set of parameters.
- `lens::Lens` specifies which parameter axis among `par` is used for continuation. For example, if `par = (α = 1.0, β = 1)`, we can perform continuation w.r.t. `α` by using `lens = (@lens _.α)`. If you have an array `par = [ 1.0, 2.0]` and want to perform continuation w.r.t. the first variable, you can use `lens = (@lens _.[1])`. For more information, we refer to `SetField.jl`.
- `contParams` parameters for continuation. See [`ContinuationPar`](@ref) for more information about the options
- `plot = false` whether to plot the solution while computing
- `printSolution = (x, p) -> norm(x)` function used to plot in the continuation curve. It is also used in the way results are saved. It could be `norm` or `x -> x[1]`. This is also useful when saving several huge vectors is not possible for memory reasons (for example on GPU...).
- `plotSolution = (x, p; kwargs...) -> nothing` function implementing the plot of the solution.
- `finaliseSolution = (z, tau, step, contResult) -> true` Function called at the end of each continuation step. Can be used to alter the continuation procedure (stop it by returning false), saving personal data, plotting... The notations are ``z=(x,p)``, `tau` is the tangent at `z` (see below), `step` is the index of the current continuation step and `ContResult` is the current branch. Note that you can have a better control over the continuation procedure by using an iterator, see [Iterator Interface](@ref).
- `callbackN` callback for newton iterations. see docs for `newton`. Can be used to change preconditioners
- `tangentAlgo = SecantPred()` controls the algorithm used to predict the tangents along the curve of solutions or the corrector. Can be `NaturalPred`, `SecantPred` or `BorderedPred`. See below for more information.
- `linearAlgo = BorderingBLS()`. Used to control the way the extended linear system associated to the continuation problem is solved. Can be `MatrixBLS`, `BorderingBLS` or `MatrixFreeBLS`.
- `verbosity::Int` controls the amount of information printed during the continuation process. Must belong to `{0,1,2,3}`
- `normC = norm` norm used in the different Newton solves
- `dotPALC = (x, y) -> dot(x, y) / length(x)`, dot product used to define the weighted dot product (resp. norm) ``\\|(x, p)\\|^2_\\theta`` in the constraint ``N(x, p)`` (see below). This arguement can be used to remove the factor `1/length(x)` for example in problems where the dimension of the state space changes (mesh adaptation, ...)
- `filename` name of a file to save the computed branch during continuation. The identifier .jld2 will be appended to this filename

# Outputs:
- `contres::ContResult` composite type which contains the computed branch. See [`ContResult`](@ref) for more information.
- `u::BorderedArray` the last solution computed on the branch

!!! tip "Controlling the argument `linearAlgo`"
    In this simplified interface to `continuation`, the argument `linearAlgo` is internally overwritten to provide a valid argument to the algorithm. If you do not want this to happen, call directly `continuation(F, J, x0, par, lens, contParams, linearAlgo; kwargs...)`.

# Simplified call:
You can also use the following call for which the jacobian is computed internally using Finite Differences

	continuation(Fhandle, x0, par, lens, contParams::ContinuationPar; kwargs...)

# Method

## Bordered system of equations

In what follows, we abuse of notations, `p` refers to the scalar value of the parameter we perform continuation with. Hence, it should be `p = get(par, lens)`.

The pseudo-arclength continuation method solves the equation ``F(x, p) = 0`` (of dimension N) together with the pseudo-arclength constraint ``N(x, p) = \\frac{\\theta}{length(x)} \\langle x - x_0, dx_0\\rangle + (1 - \\theta)\\cdot(p - p_0)\\cdot dp_0 - ds = 0`` and ``\\theta\\in[0,1]``. In practice, a curve ``\\gamma`` of solutions is sought and is parametrised by ``s``: ``\\gamma(s) = (x(s), p(s))`` is a curve of solutions to ``F(x, p)``. This formulation allows to pass turning points (where the implicit theorem fails). In the previous formula, ``(x_0, p_0)`` is a solution for a given ``s_0``, ``\\tau_0\\equiv(dx_0, dp_0)`` is the tangent to the curve ``\\gamma`` at ``s_0``. Hence, to compute the curve of solutions, we need to solve an equation of dimension N+1 which is called a Bordered system.

!!! warning "Parameter `theta`"
    The parameter `theta` in the struct `ContinuationPar`is very important. It should be tuned for the continuation to work properly especially in the case of large problems where the ``\\langle x - x_0, dx_0\\rangle`` component in the constraint might be favoured too much. Also, large `theta`s favour `p` as the corresponding term in ``N``

The parameter ds is adjusted internally depending on the number of Newton iterations and other factors. See the function `stepSizeControl` for more information. An important parameter to adjust the magnitude of this adaptation is the parameter `a` in the struct `ContinuationPar`.

## Algorithm

The algorithm works as follows:
0. Start from a known solution ``(x_0, p_0)`` with tangent to the curve of solutions: ``(dx_0 ,dp_0)``
1. **Predictor:** set ``(x_1, p_1) = (x_0, p_0) + ds\\cdot (dx_0, dp_0)``
2. **Corrector:** solve ``F(x, p)=0,\\ N(x, p)=0`` with a (Bordered) Newton Solver with initial guess ``(x_1, p_1)``.
    - if Newton in 3. did not converge, update ds/2 ⟶ ds in ``N`` and go to 1.
3. **New tangent:** Compute a new tangent (see below) ``(dx_1, dp_1)`` and update ``N`` with it. Set ``(x_0, p_0, dx_0, dp_0) = (x_1, p_1, dx_1, dp_1)`` and return to step 2

## Natural continuation

We speak of *natural* continuation when we do not consider the constraint ``N(x, p)=0``. Knowing ``(x_0, p_0)``, we use ``x_0`` as a guess for solving ``F(x, p_1)=0`` with ``p_1`` close to ``p_0``. Again, this fails at Turning points but it can be faster to compute than the constrained case. This is set by the option `tangentAlgo = NaturalPred()` in `continuation`.

## Tangent computation (step 4)
There are various ways to compute ``(dx_1, dp_1)``. The first one is called secant and is parametrised by the option `tangentAlgo = SecantPred()`. It is computed by ``(dx_1, dp_1) = (z_1, p_1) - (z_0, p_0)`` and normalised by the norm ``\\|(x, p)\\|^2_\\theta = \\frac{\\theta}{length(x)} \\langle x,x\\rangle + (1 - \\theta)\\cdot p^2``. Another method is to compute ``(dx_1, dp_1)`` by solving solving the bordered linear system ``\\begin{bmatrix} F_x & F_p	; \\ \\frac{\\theta}{length(x)}dx_0 & (1-\\theta)dp_0\\end{bmatrix}\\begin{bmatrix}dx_1 ;  dp_1\\end{bmatrix} =\\begin{bmatrix}0 ; 1\\end{bmatrix}`` ; it is set by the option `tangentAlgo = BorderedPred()`.

## Bordered linear solver

When solving the Bordered system ``F(x, p) = 0,\\ N(x, p)=0``, one faces the issue of solving the Bordered linear system ``\\begin{bmatrix} J & a	; b^T & c\\end{bmatrix}\\begin{bmatrix}X ;  y\\end{bmatrix} =\\begin{bmatrix}R ; n\\end{bmatrix}``. This can be solved in many ways via bordering (which requires two Jacobian inverses), by forming the bordered matrix (which works well for sparse matrices) or by using a full Matrix Free formulation. The choice of method is set by the argument `linearAlgo`. Have a look at the struct `linearBorderedSolver` for more information.

## Linear Algebra

Let us discuss here more about the norm and dot product. First, the option `normC` gives a norm that is used to evaluate the residual in the following way: ``max(normC(F(x,p)), \\|N(x,p)\\|)<tol``. It is thus used as a stopping criterion for a Newton algorithm. The dot product (resp. norm) used in ``N`` and in the (iterative) linear solvers is `LinearAlgebra.dot` (resp. `LinearAlgebra.norm`). It can be changed by importing these functions and redefining it. Not that by default, the ``L^2`` norm is used. These details are important because of the constraint ``N`` which incorporates the factor `length`. For some custom composite type implementing a Vector space, the dot product could already incorporates the `length` factor in which case you should either redefine the dot product or change ``\\theta``.

## Step size control

As explained above, each time the corrector phased failed, the step size ``ds`` is halved. This has the disavantage of having lost Newton iterations (which costs time) and impose small steps (which can be slow as well). To prevent this, the step size is controlled internally with the idea of having a constant number of Newton iterations per point. This is in part controlled by the aggressiveness factor `a` in `ContinuationPar`. Further tuning is performed by using `doArcLengthScaling=true` in `ContinuationPar`. This adjusts internally ``\\theta`` so that the relative contributions of ``x`` and ``p`` are balanced in the constraint ``N``.
"""
function continuation(Fhandle, Jhandle,
					x0, par, lens::Lens,
					contParams::ContinuationPar{T, S, E};
					tangentAlgo = SecantPred(),
					linearAlgo  = BorderingBLS(),
					plot = false,
					printSolution = (x, p) -> norm(x),
					normC = norm,
					dotPALC = (x,y) -> dot(x, y) / length(x),
					plotSolution = (x, p; kwargs...) -> nothing,
					finaliseSolution = (z, tau, step, contResult) -> true,
					callbackN = (x, f, J, res, iteration, itlinear, optionsN; kwargs...) -> true,
					filename = "branch-" * string(Dates.now()),
					verbosity = 0) where {T <: Real, S <: AbstractLinearSolver, E <: AbstractEigenSolver}

	# Create a bordered linear solver using newton linear solver
	_linearAlgo = @set linearAlgo.solver = contParams.newtonOptions.linsolver

	return continuation(Fhandle, Jhandle, x0, par, lens, contParams, _linearAlgo; tangentAlgo = tangentAlgo, plot = plot, printSolution = printSolution, normC = normC, dotPALC = dotPALC, plotSolution = plotSolution, finaliseSolution = finaliseSolution, callbackN = callbackN, filename = filename, verbosity = verbosity)

end

continuation(Fhandle, x0, par, lens::Lens, contParams::ContinuationPar; kwargs...) = continuation(Fhandle, (x, p) -> finiteDifferences(u -> Fhandle(u, p), x), x0, par, lens, contParams; kwargs...)


text:
module Calculus
    using Compat
    export check_derivative,
           check_gradient,
           check_hessian,
           check_second_derivative,
           deparse,
           derivative,
           differentiate,
           hessian,
           jacobian,
           second_derivative

    # TODO: Debate type system more carefully
    # abstract type BundledFunction end
    # abstract type ScalarFunction end
    # abstract type VectorFunction end
    # abstract type ForwardDifference end
    # abstract type CentralDifference end
    # abstract type ComplexDifference end
    # abstract type GradientEstimator end
    # abstract type HessianEstimator end

    # const NonDifferentiableFunction = Function
    # type DifferentiableFunction
    #   f
    #   g
    # end
    # type TwiceDifferentiableFunction
    #   f
    #   g
    #   h
    # end
    # type NonDifferentiableBundledFunction <: BundledFunction
    #   f
    #   fstorage::Any
    # end
    # type DifferentiableBundledFunction <: BundledFunction
    #   f
    #   g
    #   fstorage::Any
    #   gstorage::Any
    # end
    # type TwiceDifferentiableBundledFunction <: BundledFunction
    #   f
    #   g
    #   h
    #   fstorage::Any
    #   gstorage::Any
    #   hstorage::Any
    # end

    include("finite_difference.jl")
    include("derivative.jl")
    include("check_derivative.jl")
    Base.@deprecate integrate(f,a,b) quadgk(f,a,b)[1]
    Base.@deprecate integrate(f,a,b,method) quadgk(f,a,b)[1]
    include("symbolic.jl")
    include("differentiate.jl")
    include("deparse.jl")
end

ext: jl


text:
"""
    get_culvert_coords(culvert_shapefile_path,
                       culvert_id_field),
                       culvert_id)

Get the X and Y coordinates for a culvert by ID

Parameters:
- `culvert_shapefile_path`: The file path to the shapefile containing culverts
- `culvert_id_field`: A Symbol for the field in the culvert shapefile to use for
uniquely identifying culverts. e.g. `Symbol("UID_FIELD")`
- `culvert_id`: The unique ID (of the field specified by `culvert_id_field`)
of the desired culvert

Returns a Tuple in the form of (X\\_COORDINATE, Y\\_COORDINATE)
"""
function get_culvert_coords(culvert_shapefile_path::String,
                            culvert_id_field::Symbol,
                            culvert_id::String)
    # Load the culverts shapefile and grab Table representation
    culverts_table = Shapefile.Table(culvert_shapefile_path)
    idxs = getproperty(culverts_table, culvert_id_field) .== culvert_id

    n_records = sum(idxs)
    if n_records == 0
        @error("No culvert exists with $(:culvert_id_field) = $(culvert_id)")
    elseif n_records > 1
        @warn("More than one culvert exists with $(:culvert_id_field) = $(culvert_id). Using the first matching record.")
    end

    culvert = Shapefile.shapes(culverts_table)[findfirst(idxs)]

    coord_x_column = culvert.x
    coord_y_row = culvert.y

    return coord_x_column, coord_y_row
end

"""
    extract_landcover(raster_path, side_length, x_coordinate, y_coordinate)

Extract a square from a landcover raster of side length `side_length`, centered on
the point at [y_coordinate, x_coordinate]. `side_length` will be automatically
adjusted (by no more that the width of single pixel) if needed to ensure that
the result is an odd number of pixels in width and height. This makes it play
nicer with `Omniscape.clip`.

Parameters:
- `raster_path`: The file path to the GeoTIFF raster containing landcover info.
- `side_length`: The side length of the square to extract in units of the projection of the raster.
- `x_coordinate`: X coordinate of the center of the square in units of the projection of the raster.
- `y_coordinate`: Y coordinate of the center of the square in units of the projection of the raster.

Returns a `GeoData.GeoArray`
"""
function extract_landcover(raster_path::String,
                           side_length::Number,
                           x_coordinate::Number,
                           y_coordinate::Number)
    # Load the GDALarray (GeoData.jl)
    landcover = GDALarray(raster_path)
    # Get dims to use for subsetting
    lc_dims = dims(landcover)
    lc_size = size(landcover)

    # Get dim values
    x_dim = lc_dims[1].val
    y_dim = lc_dims[2].val

    if iseven(Int(round(side_length / step(x_dim))))
        side_length = side_length + step(x_dim)
    end

    # Extract full extent of raster
    # Need to add step to get extent because entries in
    # x_dim and y_dim correspond to left and top edges of pixel, so need to
    # incorporate width of last row and column of pixels to get full extent
    west_bound = minimum(x_dim)
    east_bound = maximum(x_dim) + step(x_dim)
    south_bound = minimum(y_dim) + step(y_dim)
    north_bound = maximum(y_dim)

    # snap center point of square (that will be extracted from `landcover`)
    # to the nearest pixel center
    pixel_center_coord_x_column = maximum(x_dim[x_dim .< x_coordinate]) + step(x_dim)/2
    pixel_center_coord_y_row = maximum(y_dim[y_dim .< y_coordinate]) + step(y_dim)/2

    # Get extent boundary for subsetting the landcover raster
    # Use min and max to ensure extent isn't beyond the full raster's extent
    landcover_west_bound = max(pixel_center_coord_x_column - side_length / 2, west_bound)
    landcover_east_bound = min(pixel_center_coord_x_column +  side_length / 2, east_bound)
    landcover_south_bound = max(pixel_center_coord_y_row -  side_length / 2, south_bound)
    landcover_north_bound = min(pixel_center_coord_y_row +  side_length / 2, north_bound)

    # Add tiny padding to extent bounds and subset to get a GeoData.GeoArray
    # The padding ensures all pixels within the bounds are extracted
    landcover_ga = landcover[Lat(Between(landcover_south_bound - EXTENT_NUDGE, landcover_north_bound + EXTENT_NUDGE)),
                             Lon(Between(landcover_west_bound - EXTENT_NUDGE, landcover_east_bound + EXTENT_NUDGE))]

    # landcover_array = Array{Float64, 2}(landcover_ga.data[:,:,1])

    return permutedims(landcover_ga, (Lat, Lon, Band)), abs(step(y_dim)), pixel_center_coord_x_column, pixel_center_coord_y_row
end

"""
    prep_source_strength(source_strength_path, landcover_ga)

Parameters:
- `source_strength_path`: The file path to the layer representing source
strength.
- `landcover_ga`: The landcover GeoData.GeoArray to use as a snap raster for
resampling and masking the source strength raster.
"""
function prep_source_strength(source_strength_path::String,
                              landcover_ga::GeoArray)
    wkt = crs(landcover_ga).val
    lc_dims = dims(landcover_ga)

    # Get dim values
    # expects a north-side-up GeoArray, so 1st (rows) dim should
    # correspond to Latitude, y
    y_dim = lc_dims[1].val
    x_dim = lc_dims[2].val

    # Extract bounds
    # Need to add step to last entries to get extent because entries in
    # x_dim and y_dim correspond to left and top edge of pixel, so need to
    # incorporate width of last row and column of pixels to get full extent
    west_bound = minimum(x_dim)
    east_bound = maximum(x_dim) + step(x_dim)
    south_bound = minimum(y_dim) + step(y_dim)
    north_bound = maximum(y_dim)

    ArchGDAL.read(source_strength_path) do source
        ArchGDAL.gdalwarp([source], ["-t_srs", "$(wkt)",
                                     "-te", "$(west_bound + EXTENT_NUDGE)",
                                     "$(south_bound + EXTENT_NUDGE)",
                                     "$(east_bound + EXTENT_NUDGE)",
                                     "$(north_bound + EXTENT_NUDGE)",
                                     "-te_srs", "$(wkt)",
                                     "-tr", "$(step(x_dim))", "$(step(y_dim))",
                                     "-r", "cubic"]) do warped
            band = ArchGDAL.getband(warped, 1)

            geotransform = ArchGDAL.getgeotransform(warped)

            # Extract the array
            array_t = ArchGDAL.read(band)

            # This handles UInt tiff rasters that can still have negative NoData values
            # Need to convert the NoData value to Int64 in these cases
            if eltype(array_t) <: Integer
                ras_type = Int64
            else
                ras_type = eltype(array_t)
            end

            # Extract no data value, first converting it to the proper type (based on
            # the raster). Then, need to convert to Float64. Weird, yes,
            # but it's the only way I could get it to work for all raster types... -VL
            nodata_val = convert(Float64, convert(ras_type, ArchGDAL.getnodatavalue(band)))

            # Transpose the array -- ArchGDAL returns a x by y array, need y by x
            # Line to handle NaNs in datasets read from tifs
            array_t[isnan.(array_t)] .= nodata_val

            array = convert(Array{Union{Float64, Missing}, 2}, permutedims(array_t, [2, 1]))

            array[array .== nodata_val] .= missing

            array, geotransform
        end
    end
end

"""
Prep inputs and run Omniscape for a single culvert
"""
function compute_connectivity(culvert_shapefile_path::String,
                              culvert_id_field::Symbol,
                              culvert_id::String,
                              landcover_path::String,
                              side_length_meters::Number,
                              source_strength_path::String,
                              mask_values::Array{Float64, 1},
                              os_args::Dict{String, String},
                              reclass_table::Array{Union{Float64, Missing}, 2}
                              )
    # Get culvert coordinates
    x_coordinate, y_coordinate = get_culvert_coords(
        culvert_shapefile_path,
        culvert_id_field,
        culvert_id
    )

    # Prep landcover
    lc_output = extract_landcover(landcover_path, side_length_meters,
                                  x_coordinate, y_coordinate)

    landcover_ga, resolution, pixel_x_column, pixel_y_row =  lc_output

    ### read in source and warp
    source_strength, geotransform = prep_source_strength(source_strength_path, landcover_ga)

    # Convert landcover to Array{Union{Float64, Missing}, 2} for Omniscape
    landcover = Array{Union{Missing, Float64}, 2}(landcover_ga.data[:,:,1])
    lc_missingval = Float64(landcover_ga.missingval)
    landcover[landcover .== lc_missingval] .= missing

    # Mask source strength with landcover, and give source strength a missing val where
    # landcover is missing
    source_strength[(!).(coalesce.(map(x -> in(x, mask_values), landcover), true))] .= missing
    source_strength[ismissing.(landcover)] .= missing
    # Grab wkt
    wkt = crs(landcover_ga).val

    # Objects to pass to Omniscape.jl, landcover (resistance), source_strength,
    # geotransform, wkt

    # Clip inputs
    center_pixel = Int(size(landcover)[1] / 2 + 0.5)
    landcover_clipped = Omniscape.clip(landcover, x = center_pixel, y = center_pixel, distance = center_pixel - 1)
    sources_clipped = Omniscape.clip(source_strength, x = center_pixel, y = center_pixel, distance = center_pixel - 1)

    if isdir("$(os_args["project_name"])")
        rm("$(os_args["project_name"])", force = true, recursive = true)
    end

    cum_current = run_omniscape(os_args, landcover_clipped,
                                reclass_table = reclass_table,
                                source_strength = sources_clipped,
                                wkt = wkt,
                                geotransform = geotransform,
                                write_outputs = true)

    cum_current, geotransform, pixel_x_column, pixel_y_row
end


text:
module Freenect

using libfreenect_jll
using DocStringExtensions

include("./export_utils.jl")
include("./freenect_sync.jl")
include("./pointcloud.jl")

end


text:
module ImageBinarization

using ImageContrastAdjustment
using ColorTypes
using ColorVectorSpace
using LinearAlgebra
using HistogramThresholding
using Polynomials

abstract type BinarizationAlgorithm end
struct Otsu <: BinarizationAlgorithm end
struct Polysegment <: BinarizationAlgorithm end
struct MinimumIntermodes <: BinarizationAlgorithm end
struct Intermodes <: BinarizationAlgorithm end

include("otsu.jl")
include("polysegment.jl")
include("minimum.jl")
include("intermodes.jl")

export
	# main functions
    binarize,
	binarize!,
	Otsu,
	Polysegment,
  	MinimumIntermodes,
  	Intermodes
end # module


text:
module CrystalAtoms
using Crystals.Structures: Crystal, position_for_crystal
export eachatom
import DataFrames

""" Wrapper around a row/atom in a Crystal for iteration """
immutable CrystalAtom{PARENT <: Crystal}
    """ Parent crystal """
    parent::PARENT
    """ Index in parent """
    index::Int64
end

Base.getindex(a::CrystalAtom, c::Any) = getindex(a.parent, a.index, c)
Base.getindex(a::CrystalAtom, c::Symbol, d::Integer) = getindex(a.parent, a.index, c, d)
Base.setindex!(a::CrystalAtom, v::Any, c::Any) = setindex!(a.parent, v, a.index, c)
function Base.setindex!(atom::CrystalAtom, v::Any, c::Symbol, d::Integer)
    setindex!(atom.parent, v, atom.index, c, d)
end
Base.names(a::CrystalAtom) = names(a.parent)

""" Wraps crystal for iteration purposes """
immutable AtomicIterator{PARENT <: Crystal}
    """ Parent crystal structure """
    parent::PARENT
end

""" Iterator over each atom in the crystal """
eachatom(crystal::Crystal) = AtomicIterator(crystal)
""" Iterator over each atom in the crystal """
DataFrames.eachrow(crystal::Crystal) = eachatom(crystal)

Base.start(itr::AtomicIterator) = 1
Base.done(itr::AtomicIterator, i::Integer) = i > size(itr.parent, 1)
Base.next(itr::AtomicIterator, i::Integer) = (CrystalAtom(itr.parent, i), i + 1)
Base.size(itr::AtomicIterator) = (length(itr.parent), )
Base.length(itr::AtomicIterator) = size(itr.parent, 1)
Base.getindex(itr::AtomicIterator, i::Any) = CrystalAtom(itr.parent, i)
Base.eltype{T <: Crystal}(::Type{AtomicIterator{T}}) = CrystalAtom{T}
end


text:
using Random
using JLD2

"""
    AgentsIO.to_serializable(t)

Return the serializable form of the passed value. This defaults to the value itself,
unless a more specific method is defined. Define a method for this function and for
[`AgentsIO.from_serializable`](@ref) if you need custom serialization for model
properties. This also enables passing keyword arguments to [`AgentsIO.load_checkpoint`](@ref)
and having access to them during deserialization of the properties. Some possible
scenarios where this may be required are:

- Your properties contain functions (or any type not supported by JLD2.jl). These may not
  be (de)serialized correctly. This could result in checkpoint files that cannot be loaded
  back in, or contain reconstructed types that do not retain their data/functionality.
- Your properties contain data that can be recalculated during deserialization. Omitting
  such properties can reduce the size of the checkpoint file, at the expense of some extra
  computation at deserialization.

If your model properties do not fall in the above scenarios, you do not need to use this
function.

This function is not called recursively on every type/value during serialization. The final
serialization functionality is enabled by JLD2.jl. To define custom serialization
for every occurence of a specific type (such as agent structs), refer to the
Custom Serialization section of JLD2.jl documentation.
"""
to_serializable(t) = t

"""
    AgentsIO.from_serializable(t; kwargs...)

Given a value in its serializable form, return the original version. This defaults
to the value itself, unless a more specific method is defined. Define a method for
this function and for [`AgentsIO.to_serializable`](@ref) if you need custom
serialization for model properties. This also enables passing keyword arguments
to [`AgentsIO.load_checkpoint`](@ref) and having access to them through `kwargs`.

Refer to [`AgentsIO.to_serializable`](@ref) to check when you need to define this function.

This function is not called recursively on every type/value during deserialization. The final
serialization functionality is enabled by JLD2.jl. To define custom serialization
for every occurence of a specific type (such as agent structs), refer to the
Custom Serialization section of JLD2.jl documentation.
"""
from_serializable(t; kwargs...) = t

struct SerializableABM{S,A<:AbstractAgent,P,R<:AbstractRNG}
    agents::Vector{A}
    space::S
    properties::P
    rng::R
    maxid::Int64
end

struct SerializableGridSpace{D,P}
    dims::NTuple{D,Int}
    metric::Symbol
end

struct SerializableContinuousSpace{D,P,T<:AbstractFloat}
    grid::SerializableGridSpace{D,P}
    dims::NTuple{D,Int}
    spacing::T
    extent::NTuple{D,T}
end

struct SerializableGraphSpace{G}
    graph::G
end

struct OSMAgentPositionData
    id::Int
    pos::Tuple{NTuple{2,Float64},NTuple{2,Float64},Float64}
    dest::Tuple{NTuple{2,Float64},NTuple{2,Float64},Float64}
    route::Vector{NTuple{2,Float64}}
end

struct SerializableOSMSpace
    agents::Vector{OSMAgentPositionData}
end

struct SerializableAStar{D,P,M,T,C}
    agent_paths::Vector{Tuple{Int,Vector{NTuple{D,T}}}}
    dims::NTuple{D,T}
    neighborhood::Vector{Dims{D}}
    admissibility::Float64
    walkmap::BitArray{D}
    cost_metric::C
end

JLD2.writeas(::Type{Pathfinding.AStar{D,P,M,T,C}}) where {D,P,M,T,C} = SerializableAStar{D,P,M,T,C}

function to_serializable(t::ABM{S}) where {S}
    sabm = SerializableABM(
        collect(allagents(t)),
        to_serializable(t.space),
        to_serializable(t.properties),
        t.rng,
        t.maxid.x,
    )
    if S <: OSM.OpenStreetMapSpace
        for i in 1:nagents(t)
            sabm.agents[i] = typeof(sabm.agents[i])(
                (
                    getproperty(sabm.agents[i], x) for x in fieldnames(typeof(sabm.agents[i]))
                )...,
            )
            sabm.agents[i].route = []
        end

        for a in allagents(t)
            push!(
                sabm.space.agents,
                OSMAgentPositionData(
                    a.id,
                    (OSM.latlon(a.pos[1], t), OSM.latlon(a.pos[2], t), a.pos[3]),
                    (
                        OSM.latlon(a.destination[1], t),
                        OSM.latlon(a.destination[2], t),
                        a.destination[3],
                    ),
                    [OSM.latlon(i, t) for i in a.route],
                ),
            )
        end
    end
    return sabm
end

to_serializable(t::GridSpace{D,P}) where {D,P} =
    SerializableGridSpace{D,P}(size(t.s), t.metric)

to_serializable(t::ContinuousSpace{D,P,T}) where {D,P,T} =
    SerializableContinuousSpace{D,P,T}(to_serializable(t.grid), t.dims, t.spacing, t.extent)

to_serializable(t::GraphSpace{G}) where {G} = SerializableGraphSpace{G}(t.graph)

to_serializable(t::OSM.OpenStreetMapSpace) = SerializableOSMSpace([])

JLD2.wconvert(::Type{SerializableAStar{D,P,M,T,C}}, t::Pathfinding.AStar{D,P,M,T,C}) where {D,P,M,T,C} =
    SerializableAStar{D,P,M,T,C}(
        [(k, collect(v)) for (k, v) in t.agent_paths],
        t.dims,
        map(Tuple, t.neighborhood),
        t.admissibility,
        t.walkmap,
        t.cost_metric,
    )

function from_serializable(t::SerializableABM{S,A}; kwargs...) where {S,A}
    abm = ABM(
        A,
        from_serializable(t.space; kwargs...);
        scheduler = get(kwargs, :scheduler, Schedulers.fastest),
        properties = from_serializable(t.properties; kwargs...),
        rng = t.rng,
        warn = get(kwargs, :warn, true),
    )
    abm.maxid[] = t.maxid

    if S <: SerializableOSMSpace
        agentdata = Dict(a.id => a for a in t.space.agents)
        for a in t.agents
            a.pos = (
                OSM.intersection(agentdata[a.id].pos[1], abm)[1],
                OSM.intersection(agentdata[a.id].pos[2], abm)[1],
                agentdata[a.id].pos[3],
            )
            a.destination = (
                OSM.intersection(agentdata[a.id].dest[1], abm)[1],
                OSM.intersection(agentdata[a.id].dest[2], abm)[1],
                agentdata[a.id].dest[3],
            )
            a.route = [OSM.intersection(i, abm)[1] for i in agentdata[a.id].route]
        end
    end

    for a in t.agents
        add_agent_pos!(a, abm)
    end
    return abm
end

function from_serializable(t::SerializableGridSpace{D,P}; kwargs...) where {D,P}
    s = Array{Vector{Int},D}(undef, t.dims)
    for i in eachindex(s)
        s[i] = Int[]
    end
    return GridSpace{D,P}(s, t.metric, Dict(), Dict())
end

function from_serializable(t::SerializableContinuousSpace{D,P,T}; kwargs...) where {D,P,T}
    update_vel! = get(kwargs, :update_vel!, Agents.defvel)
    ContinuousSpace(
        from_serializable(t.grid; kwargs...),
        update_vel!,
        t.dims,
        t.spacing,
        t.extent,
    )
end

from_serializable(t::SerializableGraphSpace; kwargs...) = GraphSpace(t.graph)

function from_serializable(t::SerializableOSMSpace; kwargs...)
    @assert haskey(kwargs, :map) "Path to OpenStreetMap not provided"

    OSM.OpenStreetMapSpace(
        get(kwargs, :map, OSM.TEST_MAP);   # Should never need default value
        use_cache = get(kwargs, :use_cache, false),
        trim_to_connected_graph = get(kwargs, :trim_to_connected_graph, true),
    )
end

JLD2.rconvert(::Type{Pathfinding.AStar{D,P,M,T,C}}, t::SerializableAStar{D,P,M,T,C}) where {D,P,M,T,C} =
    Pathfinding.AStar{D,P,M,T,C}(
        Dict{Int,Pathfinding.Path{D,T}}(
            k => Pathfinding.Path{D,T}(v...) for (k, v) in t.agent_paths
        ),
        t.dims,
        map(CartesianIndex, t.neighborhood),
        t.admissibility,
        t.walkmap,
        t.cost_metric,
    )

"""
    AgentsIO.save_checkpoint(filename, model::ABM)

Write the entire `model` to file specified by `filename`. The following points
should be considered before using this functionality:

- OpenStreetMap data is not saved. The path to the map should be specified when loading
  the model using the `map` keyword of [`AgentsIO.load_checkpoint`](@ref).
- Functions are not saved, including stepping functions, schedulers, and `update_vel!`.
  The last two can be provided to [`AgentsIO.load_checkpoint`](@ref) using the appropriate
  keyword arguments.
"""
function save_checkpoint(filename, model::ABM)
    model = to_serializable(model)
    @save filename model
end

"""
    AgentsIO.load_checkpoint(filename; kwargs...)

Load the model saved to the file specified by `filename`.

## Keywords
- `scheduler = Schedulers.fastest` specifies what scheduler should
  be used for the model.
- `warn = true` can be used to disable warnings from type checks on the
    agent type.
[`ContinuousSpace`](@ref) specific:
- `update_vel!` specifies a function that should be used to
  update each agent's velocity before it is moved. Refer to [`ContinuousSpace`](@ref) for
  details.
[`OpenStreetMapSpace`](@ref) specific:
- `map` is a path to the OpenStreetMap to be used for the space. This is a required parameter
  if the space is [`OpenStreetMapSpace`](@ref).
- `use_cache = false`, `trim_to_connected_graph = true` refer to [`OpenStreetMapSpace`](@ref)
"""
function load_checkpoint(filename; kwargs...)
    @load filename model
    return from_serializable(model; kwargs...)
end


text:
abstract type AbstractMatrixFreeSolver <: AbstractDisplacementSolver end

@params mutable struct StaticMatrixFreeDisplacementSolver{T,dim,TP<:AbstractPenalty{T}} <:
                       AbstractDisplacementSolver
    elementinfo::ElementFEAInfo{dim}
    problem::StiffnessTopOptProblem{dim}
    f::AbstractVector{T}
    meandiag::T
    u::AbstractVector{T}
    lhs::AbstractVector{T}
    rhs::AbstractVector{T}
    vars::AbstractVector{T}
    xes::Any
    fixed_dofs::Any
    free_dofs::Any
    penalty::TP
    prev_penalty::TP
    xmin::T
    cg_max_iter::Integer
    tol::T
    cg_statevars::CGStateVariables{T}
    preconditioner::Any
    preconditioner_initialized::Base.RefValue{Bool}
    conv::Any
end
function Base.show(
    ::IO, ::MIME{Symbol("text/plain")}, x::StaticMatrixFreeDisplacementSolver
)
    return println("TopOpt matrix free conjugate gradient iterative solver")
end
function StaticMatrixFreeDisplacementSolver(sp, args...; kwargs...)
    return StaticMatrixFreeDisplacementSolver(whichdevice(sp), sp, args...; kwargs...)
end

function StaticMatrixFreeDisplacementSolver(
    ::CPU,
    sp::StiffnessTopOptProblem{dim,T};
    conv=DefaultCriteria(),
    xmin=one(T) / 1000,
    cg_max_iter=700,
    tol=xmin,
    penalty=PowerPenalty{T}(1),
    prev_penalty=deepcopy(penalty),
    preconditioner=identity,
    quad_order=2,
) where {dim,T}
    elementinfo = ElementFEAInfo(sp, quad_order, Val{:Static})
    if eltype(elementinfo.Kes) <: Symmetric
        f = x -> sumdiag(rawmatrix(x).data)
    else
        f = x -> sumdiag(rawmatrix(x))
    end
    meandiag = mapreduce(f, +, elementinfo.Kes; init=zero(T))
    xes = deepcopy(elementinfo.fes)

    u = zeros(T, ndofs(sp.ch.dh))
    f = similar(u)
    lhs = similar(u)
    rhs = similar(u)
    vars = fill(one(T), getncells(sp.ch.dh.grid) - sum(sp.black) - sum(sp.white))
    cg_statevars = CGStateVariables{eltype(u),typeof(u)}(copy(u), similar(u), similar(u))

    fixed_dofs = sp.ch.prescribed_dofs
    free_dofs = setdiff(1:length(u), fixed_dofs)
    return StaticMatrixFreeDisplacementSolver(
        elementinfo,
        sp,
        f,
        meandiag,
        u,
        lhs,
        rhs,
        vars,
        xes,
        fixed_dofs,
        free_dofs,
        penalty,
        prev_penalty,
        xmin,
        cg_max_iter,
        tol,
        cg_statevars,
        preconditioner,
        Ref(false),
        conv,
    )
end

MatrixFreeOperator(solver::StaticMatrixFreeDisplacementSolver) = buildoperator(solver)
function buildoperator(solver::StaticMatrixFreeDisplacementSolver)
    penalty = getpenalty(solver)
    @unpack elementinfo, meandiag, vars, xmin, fixed_dofs, free_dofs, xes, conv = solver
    return MatrixFreeOperator(
        solver.f,
        elementinfo,
        meandiag,
        vars,
        xes,
        fixed_dofs,
        free_dofs,
        xmin,
        penalty,
        conv,
    )
end

function (s::StaticMatrixFreeDisplacementSolver)(;
    assemble_f=true, rhs=assemble_f ? s.f : s.rhs, lhs=assemble_f ? s.u : s.lhs, kwargs...
)
    if assemble_f
        assemble_f!(s.f, s.problem, s.elementinfo, s.vars, getpenalty(s), s.xmin)
    end
    matrix_free_apply2f!(
        rhs, s.elementinfo, s.meandiag, s.vars, s.problem, getpenalty(s), s.xmin
    )

    @unpack cg_max_iter, cg_statevars = s
    @unpack preconditioner_initialized, preconditioner, tol = s
    operator = buildoperator(s)

    if !(preconditioner === identity)
        if !preconditioner_initialized[]
            UpdatePreconditioner!(preconditioner, operator)
            preconditioner_initialized[] = true
        end
    end
    if preconditioner === identity
        return cg!(
            lhs,
            operator,
            rhs;
            tol=tol,
            maxiter=cg_max_iter,
            log=false,
            statevars=cg_statevars,
            initially_zero=false,
        )
    else
        return cg!(
            lhs,
            operator,
            rhs;
            tol=tol,
            maxiter=cg_max_iter,
            log=false,
            statevars=cg_statevars,
            initially_zero=false,
            Pl=preconditioner,
        )
    end
end


text:
using Test

@testset "All tests" begin
    include("models_1d.jl")
    include("inbuilt_models.jl")
    include("basic_properties.jl")
    include("io.jl")
    include("conversion.jl")
end


text:
type Setting
    img
    row::Int
    col::Int
end

function mod_fun(row_map, arr, col, scores)
    j = rand(1:size(arr)[1])        # We randomly choose a column.
    while arr[j, 1] == -1 || j == col
        j = rand(1:size(arr)[1])
    end
    elts = row_map[j]               # We get the elements on j-th column.
    nb_elt_on_col = size(elts)[1]   # We get the number of elements.
    start_index = arr[j, 3]         # We get the nearest index from the start point.
    # Randomly selecting a white pixel according to the normal distribution.
    rnd = randn()
    index = int(rnd * (nb_elt_on_col / 2) + start_index)
    index = max(index, 1)
    index = min(index, nb_elt_on_col)

    i = elts[index]

    return (i,j,index)
end

function cost_fun(orig, new, arr, pos, start)
    # Compute the local cost of the modification
    # Orig : the original image
    # new : the modified image
    # v : the candidate value for position `pos`
    # pos : position of the modification
    # start : hint about the start of the line given to the program

    rows, cols = size(new)
    row, col = start
    r, c = pos

    alpha = 1
    beta = 20

    cost = 0
    # We score the distance between the start point and the pixel.
    if c > 1 && arr[c - 1, 1] == -1
        cost += alpha * abs(r - row)
    end
    # We score the distance between the previous marked pixel and the current
    # one, weighted with the square of root of the distance with the start
    # pixel. We want the relation between neighbors to be important, but we
    # also want to be close to the line to match.
    if c > 1 && arr[c - 1, 1] != -1
        cost += beta * abs(r - arr[c - 1, 1]) * (abs(r - row) + 1) ^ (1/2)
    end

    cost
end

function extract(file_path, row, col, out_path)

    println("extracting line from: $(file_path)")
    println("starting point: (row = $(row), col = $(col))")

    # Read input image to sparse matrix
    img = convert(Array, float64(imread(file_path)))
    img[row, col] = 255                       # In order to be sure there is at least one pixel on this column.
    orig = sparse(img)

    # Get number of rows and columns of the image.
    rows, cols = size(orig)

    # Candidate
    I, J, V = findnz(orig)
    arr = zeros(cols, 3)                      # arr contains the tuples (row, index, start_index) marked, one by column.
    row_map = map((x) -> I[J .== x], 1:cols)  # row_map contains the indexes of the white pixels.

    # Initial random candidate
    for j = 1:cols
        elts = row_map[j]                     # We get elements on j-th column.
        nb_elt_on_col = size(elts)[1]         # We get the number of elements.

        start_index = 0
        # If there isn't any white pixel on the column, we don't mark it.
        if (nb_elt_on_col == 0)
            i = -1
            index = -1
        else
            # We initialize with the nearest (row speaking) pixel from the
            # start pixel.
            start_index = max(size(row_map[j][row_map[j] .<= row], 1), 1)
            # Then we mark a random one, according to the normal distribution.
            rnd = randn()
            index = int(rnd * (nb_elt_on_col / 2) + start_index)
            index = max(index, 1)
            index = min(index, nb_elt_on_col)

            i = elts[index]
        end
        arr[j,:] = [i, index, start_index]    # We mark this pixel.
    end
    # Set start point.
    arr[col, 1] = row

    obj = sparse(I, J, V)

    ###########################
    # Run simulated annealing #
    ###########################

    # Simulated annealing settings
    t = 2500.0
    t_step = 0.999999
    t_stop = 5
    max_epoc = 10000000000
    epoc = 0

    println("Compute score matrix")

    score = zeros(Float64, cols)
    for i = 1:cols
        score[i] = cost_fun(orig, obj, arr, (arr[i], i), (row, col))
    end

    tot_score = sum(score)
    println("Total score: $(tot_score)")

    println("Run Simulated Annealing")

    while epoc < max_epoc && t > t_stop

        # Choose neighbor
        r,c,index = mod_fun(row_map, arr, col, score)

        # Compute local score
        score_mod = cost_fun(orig, obj, arr, (r, c), (row, col))

        # Compute energy delta
        mod_e = score_mod - score[c]

        # Keep the candidate or not
        if mod_e < 0 || rand() <= exp(-mod_e / t)
            arr[c,:] = [r, index, arr[c,3]]
            tot_score = tot_score - score[c] + score_mod
            score[c] = score_mod
        end

        # Decrease temperature
        if epoc % 100 == 0
            t *= t_step
        end
        if epoc % 1000000 == 0
            println("Iteration number: $(epoc); cost: $(tot_score), temp: $(t)")
        end

        epoc += 1
    end

    println("temperature: $t")
    println("energy: $tot_score")
    println("epocs: $epoc")

    ###################
    # Display results #
    ###################

    img *= 0.0005
    for i = 1:cols
        if (arr[i] != -1)
            r = int(arr[i])
            img[r, i] = 255
        end
    end

    # Display result
    img2 = convert(Image, img)
    display(img2)
    #imwrite(img2, out_path)
end


text:
using FundStabABM.Params

@testset "Parameters" begin
    @test Params.default().bign >= Params.default().bigk >= 1
    @test Params.default().bigt >= 1
    @test -1< Params.default().drift <1
    @test Params.default().marketvol >= 0
    @test 0 < Params.default().perfwindow[1] <=
        Params.default().perfwindow[end] < Params.default().bigt
    @test Params.default().marketstartval > 0
    @test Params.default().marketvol >= 0
    @test Params.default().betamean > 0
    @test Params.default().betastd >= 0
    @test 0 <= Params.default().stockvolrange[1] <=
        Params.default().stockvolrange[end]
    @test 0 <=Params.default().invcaprange[1] <= Params.default().invcaprange[2]
    @test Params.default().stockstartval >= 0
    @test -1 < Params.default().thresholdmean < 1
    @test Params.default().thresholdstd >= 0
    @test 1 <=Params.default().portfsizerange[1] <=
        Params.default().portfsizerange[end] <=Params.default().bigm
    @test 0 <= Params.default().impactrange[1] <=Params.default().impactrange[2]

end


text:
using DataStructures
using Compat
using Compat: @warn
using IterableTables: getiterator
using DataFrames

import Phylo.API: _nodetype, _branchtype, _getnodes, _getbranches
import Phylo.API: _getnodenames, _getbranchnames, _getleafnames
import Phylo.API: _getleafinfo, _setleafinfo!, _leafinfotype
import Phylo.API: _resetleaves!, _getnoderecord, _setnoderecord!
import Phylo.API: _addnode!, _deletenode!, _addbranch!, _deletebranch!, _validate
import Phylo.API: _hasrootheight, _getrootheight, _setrootheight!, _clearrootheight!
import Phylo.API: _getnode, _getbranch, _setnode!, _setbranch!, _nleaves

_branchtype(::AbstractBranchTree{NL, BL}) where {NL, BL} = Branch{NL}

"""
    BinaryTree

Binary phylogenetic tree object with known leaves and per node data
"""
mutable struct BinaryTree{LI, ND} <: AbstractBranchTree{String, Int}
    nodes::OrderedDict{String, BinaryNode{Int}}
    branches::Dict{Int, Branch{String}}
    leafinfos::LI
    noderecords::OrderedDict{String, ND}
    rootheight::Float64
end

_leafinfotype(::BinaryTree{LI, ND}) where {LI, ND} = LI
_nleaves(tree::BinaryTree{LI, ND}) where {LI, ND} =
    length(nodefilter(_isleaf, tree))

function BinaryTree(lt::BinaryTree{LI, ND};
                    copyinfo=false, empty=true) where {LI, ND}
    validate(lt) || error("Tree to copy is not valid")
    leafnames = getleafnames(lt)
    # Leaf records may be conserved across trees, as could be invariant?
    leafinfos = copyinfo ? deepcopy(lt.leafinfos) : lt.leafinfos
    if empty # Empty out everything else
        nodes = OrderedDict(map(leaf -> leaf => BinaryNode{Int}(), leafnames))
        branches = Dict{Int, Branch{String}}()
        noderecords = OrderedDict(map(leaf -> leaf => ND(), leafnames))
    else # Make copies of everything
        nodes = deepcopy(getnodes(lt))
        noderecords = deepcopy(lt.noderecords)
        branches = deepcopy(getbranches(lt))
    end
    return BinaryTree{LI, ND}(nodes, branches,
                              leafinfos, noderecords, lt.rootheight)
end

function BinaryTree{LI, ND}(leaves::Vector{String},
                            treetype::Type{BinaryTree{LI, ND}} =
                            BinaryTree{LI, ND};
                            rootheight::Float64 = NaN) where {LI, ND}
    nodes = OrderedDict(map(leaf -> leaf => BinaryNode{Int}(), leaves))
    leafinfos = LI()
    noderecords = OrderedDict(map(leaf -> leaf => ND(), leaves))
    return BinaryTree{LI, ND}(nodes, OrderedDict{Int, Branch{String}}(),
                              leafinfos, noderecords, rootheight)
end

function BinaryTree{LI, ND}(numleaves::Int = 0,
                            treetype::Type{BinaryTree{LI, ND}} =
                            BinaryTree{LI, ND};
                            rootheight::Float64 = NaN) where {LI, ND}
    leaves = map(num -> "Leaf $num", 1:numleaves)
    nodes = OrderedDict(map(leaf -> leaf => BinaryNode{Int}(), leaves))
    leafinfos = LI()
    noderecords = OrderedDict(map(leaf -> leaf => ND(), leaves))
    return BinaryTree{LI, ND}(nodes, OrderedDict{Int, Branch{String}}(),
                              leafinfos, noderecords, rootheight)
end

function BinaryTree{LI, ND}(leafinfos::LI; rootheight::Float64 = NaN) where {LI, ND}
    leafnames = unique(collect(map(info -> info[1], getiterator(leafinfos))))
    nodes = OrderedDict(map(leaf -> leaf => BinaryNode{Int}(), leafnames))
    branches = Dict{Int, Branch{String}}()
    noderecords = OrderedDict(map(leaf -> leaf => ND(),
                                  leafnames))
    return BinaryTree{LI,
                      Dict{String, Any}}(nodes, branches, leafinfos,
                                         noderecords, rootheight)
end

BinaryTree(leafinfos::LI; rootheight::Float64 = NaN) where LI =
    BinaryTree{LI, Dict{String, Any}}(leafinfos; rootheight = rootheight)

_nodetype(::BinaryTree) = BinaryNode{Int}

function _getnodes(bt::BinaryTree)
    return bt.nodes
end

function _getbranches(bt::BinaryTree)
    return bt.branches
end

function _getleafinfo(bt::BinaryTree)
    return bt.leafinfos
end

function _getleafinfo(bt::BinaryTree, leafname)
    return Iterators.filter(info -> info[1] == leafname,
                            getiterator(bt.leafinfos))
end

function _setleafinfo!(bt::BinaryTree, info)
    bt.leafinfos = info
end

function _resetleaves!(bt::BinaryTree)
    bt.leafinfos = empty!(bt.leafinfos)
    return bt
end

function _getnoderecord(bt::BinaryTree, nodename)
    return bt.noderecords[nodename]
end

function _setnoderecord!(bt::BinaryTree, nodename, value)
    bt.noderecords[nodename] = value
end

function _addnode!(tree::BinaryTree{LI, NR}, nodename) where {LI, NR}
    !_hasnode(tree, nodename) ||
        error("Node $nodename already present in tree")
    _setnode!(tree, nodename, BinaryNode{Int}())
    setnoderecord!(tree, nodename, NR())
    return nodename
end

function _deletenode!(tree::BinaryTree, nodename)
    node = getnode(tree, nodename)
    if _hasinbound(node)
        deletebranch!(tree, _getinbound(node))
    end
    for b in _getoutbounds(node)
        deletebranch!(tree, b)
    end
    delete!(_getnodes(tree), nodename)
    delete!(tree.noderecords, nodename)
    return nodename
end

function _validate(tree::BinaryTree)
    if !isempty(tree.leafinfos) && length(getiterator(tree.leafinfos)) > 0
        if Set(map(info -> info[1], getiterator(tree.leafinfos))) !=
            Set(_getleafnames(tree))
            @warn "LeafInfo names do not match actual leaves of tree"
            return false
        end
    end

    if Set(nodenamefilter(_hasoutboundspace, tree)) !=
        Set(nodenamefilter(_isleaf, tree))
        @warn "Nodes must have two or zero outbound connections."
        return false
    end

    if Set(keys(tree.noderecords)) != Set(keys(getnodes(tree)))
        @warn "Leaf records do not match node records of tree"
        return false
    end

    return true
end

function _hasrootheight(tree::BinaryTree)
    return !isnan(tree.rootheight)
end

function _getrootheight(tree::BinaryTree)
    return tree.rootheight
end

function _setrootheight!(tree::BinaryTree, height::Float64)
    tree.rootheight = height
    return height
end

function _clearrootheight!(tree::BinaryTree)
    tree.rootheight = NaN
end

"""
    NamedBinaryTree

Binary phylogenetic tree object with known leaves
"""
const NamedTree = NamedBinaryTree = BinaryTree{DataFrame, Dict{String, Any}}

"""
    PolytomousTree

Phylogenetic tree object with polytomous branching, and known leaves and per node data
"""
mutable struct PolytomousTree{LI, ND} <: AbstractBranchTree{String, Int}
    nodes::OrderedDict{String, Node{Int}}
    branches::Dict{Int, Branch{String}}
    leafinfos::LI
    noderecords::OrderedDict{String, ND}
    rootheight::Float64
end

_leafinfotype(::PolytomousTree{LI, ND}) where {LI, ND} = LI
_nleaves(tree::PolytomousTree{LI, ND}) where {LI, ND} =
    length(nodefilter(_isleaf, tree))

function PolytomousTree(lt::PolytomousTree{LI, ND};
                        copyinfo=false, empty=true) where {LI, ND}
    validate(lt) || error("Tree to copy is not valid")
    leafnames = getleafnames(lt)
    # Leaf records may be conserved across trees, as could be invariant?
    leafinfos = copyinfo ? deepcopy(lt.leafinfos) : lt.leafinfos
    if empty # Empty out everything else
        nodes = OrderedDict(map(leaf -> leaf => Node{Int}(), leafnames))
        branches = Dict{Int, Branch{String}}()
        noderecords = OrderedDict(map(leaf -> leaf => ND(), leafnames))
    else # Make copies of everything
        nodes = deepcopy(getnodes(lt))
        noderecords = deepcopy(lt.noderecords)
        branches = deepcopy(getbranches(lt))
    end
    return PolytomousTree{LI, ND}(nodes, branches, leafinfos, noderecords,
                            lt.rootheight)
end

function PolytomousTree{LI, ND}(leafinfos::LI,
                                treetype::Type{PolytomousTree{LI, ND}} =
                                PolytomousTree{LI, Dict{String, Any}};
                                rootheight::Float64 = NaN) where {LI, ND}
    leafnames = unique(collect(map(info -> info[1], getiterator(leafinfos))))
    nodes = OrderedDict(map(leaf -> leaf => Node{Int}(), leafnames))
    branches = Dict{Int, Branch{String}}()
    noderecords = OrderedDict(map(leaf -> leaf => ND(), leafnames))
    return PolytomousTree{LI, ND}(nodes, branches,
                                  leafinfos, noderecords, rootheight)
end

function PolytomousTree{LI, ND}(leaves::Vector{String},
                                treetype::Type{PolytomousTree{LI, ND}} =
                                PolytomousTree{LI, ND};
                                rootheight::Float64 = NaN) where {LI, ND}
    nodes = OrderedDict(map(leaf -> leaf => Node{Int}(), leaves))
    leafinfos = LI()
    noderecords = OrderedDict(map(leaf -> leaf => ND(), leaves))
    return PolytomousTree{LI, ND}(nodes, OrderedDict{Int, Branch{String}}(),
                              leafinfos, noderecords, rootheight)
end

function PolytomousTree{LI, ND}(numleaves::Int = 0,
                                treetype::Type{PolytomousTree{LI, ND}} =
                                PolytomousTree{LI, Dict{String, Any}};
                                rootheight::Float64 = NaN) where {LI, ND}
    leaves = map(num -> "Leaf $num", 1:numleaves)
    nodes = OrderedDict(map(leaf -> leaf => Node{Int}(), leaves))
    leafinfos = LI()
    noderecords = OrderedDict(map(leaf -> leaf => ND(), leaves))
    return PolytomousTree{LI, ND}(nodes, OrderedDict{Int, Branch{String}}(),
                                  leafinfos, noderecords, rootheight)
end

PolytomousTree(leafinfos::LI; rootheight::Float64 = NaN) where LI =
    PolytomousTree{LI, Dict{String, Any}}(leafinfos; rootheight = rootheight)

_nodetype(::PolytomousTree) = Node{Int}

function _getnodes(pt::PolytomousTree)
    return pt.nodes
end

function _getbranches(pt::PolytomousTree)
    return pt.branches
end

function _getleafinfo(pt::PolytomousTree)
    return pt.leafinfos
end

function _getleafinfo(pt::PolytomousTree, leaf)
    return Iterators.filter(info -> info[1] == leafname,
        getiterator(pt.leafinfos))
end

function _setleafinfo!(pt::PolytomousTree, info)
    pt.leafinfos = info
end

function _resetleaves!(pt::PolytomousTree)
    pt.leafinfos = empty!(pt.leafinfos)
    return pt
end

function _getnoderecord(pt::PolytomousTree, nodename)
    return pt.noderecords[nodename]
end

function _setnoderecord!(pt::PolytomousTree, nodename, value)
    pt.noderecords[nodename] = value
end

function _addnode!(tree::PolytomousTree{LI, NR}, nodename) where {LI, NR}
    !_hasnode(tree, nodename) ||
        error("Node $nodename already present in tree")
    _setnode!(tree, nodename, Node{Int}())
    setnoderecord!(tree, nodename, NR())
    return nodename
end

function _deletenode!(tree::PolytomousTree, nodename)
    node = getnode(tree, nodename)
    if _hasinbound(node)
        deletebranch!(tree, _getinbound(node))
    end
    for b in _getoutbounds(node)
        deletebranch!(tree, b)
    end
    delete!(_getnodes(tree), nodename)
    delete!(tree.noderecords, nodename)
    return nodename
end

function _validate(tree::PolytomousTree)
    if length(getiterator(tree.leafinfos)) > 0
        if Set(map(info -> info[1], getiterator(tree.leafinfos))) !=
            Set(_getleafnames(tree))
            @warn "LeafInfo names do not match actual leaves of tree"
            return false
        end
    end
    if Set(keys(tree.noderecords)) != Set(keys(getnodes(tree)))
        @warn "Leaf records do not match node records of tree"
        return false
    end

    return true
end

function _hasrootheight(tree::PolytomousTree)
    return !isnan(tree.rootheight)
end

function _getrootheight(tree::PolytomousTree)
    return tree.rootheight
end

function _setrootheight!(tree::PolytomousTree, height::Float64)
    tree.rootheight = height
    return height
end

function _clearrootheight!(tree::PolytomousTree)
    tree.rootheight = NaN
end

"""
    NamedPolytomousTree

Polytomous phylogenetic tree object with known leaves
"""
const NamedPolytomousTree = PolytomousTree{DataFrame, Dict{String, Any}}

_getnodenames(tree::AbstractTree) = collect(keys(_getnodes(tree)))
_getbranchnames(tree::AbstractTree) = collect(keys(_getbranches(tree)))
#  - _hasnode()
_hasnode(tree::AbstractTree, label) = haskey(_getnodes(tree), label)
#  - _hasbranch()
_hasbranch(tree::AbstractTree, label) = haskey(_getbranches(tree), label)
#  - _addbranch!()
function _addbranch!(tree::AbstractBranchTree, source, destination, length::Float64, label)
    # Add the new branch
    _setbranch!(tree, label, Branch(source, destination, length))

    # Update the associated source and destination nodes
    _addoutbound!(getnode(tree, source), label)
    _setinbound!(getnode(tree, destination), label)

    # Return updated tree
    return label
end
#  - _deletebranch!()
function _deletebranch!(tree::AbstractBranchTree, label)
    # Find the branch
    branch = _getbranch(tree, label)
    # Remove branch reference from its source node
    _deleteoutbound!(_getnode(tree, _src(branch)), label)
    # Remove branch reference from its destination node
    _deleteinbound!(_getnode(tree, _dst(branch)), label)
    # Remove branch itself
    delete!(_getbranches(tree), label)
    # Return the branch label
    return label
end

"""
    clearrootheight(::AbstractTree)

Clears the tree's root height record.
"""
function clearrootheight!(tree::AbstractTree)
    _clearrootheight!(tree)
end

function _getnode(tree::AbstractTree, label)
    return _getnodes(tree)[label]
end

function _getbranch(tree::AbstractTree, label)
    return _getbranches(tree)[label]
end

function _setnode!(tree::AbstractTree, label, node)
    return _getnodes(tree)[label] = node
end

function _setbranch!(tree::AbstractTree, label, branch)
    _hasbranch(tree, label) &&
        error("Branch $label already exists")
    return _getbranches(tree)[label] = branch
end


text:
# using DataFrames
# using CSV
using JSON
using SHA

include("entities/AbstractDataEntries.jl")
# include("entities/BlobStores.jl")

include("services/AbstractDataEntries.jl")
include("services/DataEntryBlob.jl")
include("services/BlobStores.jl")

include("services/FileDataEntryBlob.jl")
include("services/InMemoryDataEntryBlob.jl")

export AbstractDataEntry, GeneralDataEntry, MongodbDataEntry, FileDataEntry

export getData, addData!, updateData!, deleteData!, listStoreEntries
export getDataBlob, addDataBlob!, updateDataBlob!, deleteDataBlob!, listDataBlobs
export copyStore

export getId, getHash, getCreatedTimestamp


text:
import Test
import PredictMD

a = PredictMD.package_directory()
Test.@test( isdir(a) )

b = PredictMD.package_directory(
    "assets",
    )
Test.@test( isdir(b) )
Test.@test( dirname(b) == a )

c = PredictMD.package_directory(
    "assets",
    "icd",
    )
Test.@test( isdir(c) )
Test.@test( dirname(c) == b )

d = PredictMD.package_directory(
    "assets",
    "icd",
    "icd9",
    )
Test.@test( isdir(d) )
Test.@test( dirname(d) == c )

e = PredictMD.package_directory(
    "assets",
    "icd",
    "icd9",
    "ccs"
    )
Test.@test( isdir(e) )
Test.@test( dirname(e) == d )

f = PredictMD.package_directory(
    "assets",
    "icd",
    "icd9",
    "ccs",
    "AppendixASingleDX.txt"
    )
Test.@test( dirname(f) == e )
Test.@test( isfile(f) )


text:
# demonstrates the use of just an y argument

using SixelGraphics

x = linspace(0,6*pi)
y = sin(x)

s = sixelplot(y, typ='b')


text:
# This file was generated, do not modify it. # hide
import RDatasets: dataset
import DataFrames: describe, select, Not, rename!
boston = dataset("MASS", "Boston")
first(boston, 3)


text:
using PlotlyJS, JSON

function ternary_markers()
    function make_ax(title, tickangle)
        attr(title=title, titlefont_size=20, tickangle=tickangle,
            tickfont_size=15, tickcolor="rgba(0, 0, 0, 0)", ticklen=5,
            showline=true, showgrid=true)
    end

    raw_data = [
        Dict(:journalist=>75, :developer=>:25, :designer=>0, :label=>"point 1"),
        Dict(:journalist=>70, :developer=>:10, :designer=>20, :label=>"point 2"),
        Dict(:journalist=>75, :developer=>:20, :designer=>5, :label=>"point 3"),
        Dict(:journalist=>5, :developer=>:60, :designer=>35, :label=>"point 4"),
        Dict(:journalist=>10, :developer=>:80, :designer=>10, :label=>"point 5"),
        Dict(:journalist=>10, :developer=>:90, :designer=>0, :label=>"point 6"),
        Dict(:journalist=>20, :developer=>:70, :designer=>10, :label=>"point 7"),
        Dict(:journalist=>10, :developer=>:20, :designer=>70, :label=>"point 8"),
        Dict(:journalist=>15, :developer=>:5, :designer=>80, :label=>"point 9"),
        Dict(:journalist=>10, :developer=>:10, :designer=>80, :label=>"point 10"),
        Dict(:journalist=>20, :developer=>:10, :designer=>70, :label=>"point 11")
    ]

    t = scatterternary(
        mode="markers",
        a=[_x[:journalist] for _x in raw_data],
        b=[_x[:developer] for _x in raw_data],
        c=[_x[:designer] for _x in raw_data],
        text=[_x[:label] for _x in raw_data],
        marker=attr(symbol=100, color="#DB7365", size=14, line_width=2)
    )
    layout = Layout(
        ternary=attr(
            sum=100,
            aaxis=make_ax("Journalist", 0),
            baxis=make_ax("Developer", 45),
            caxis=make_ax("Designer", -45),
            bgcolor="#fff1e0",
        ), annotations=attr(
            showarrow=false,
            text="Replica of Tom Pearson's block",
            x=1.0, y=1.3, font_size=15
        ),
        paper_bgcolor="#fff1e0"
    )
    plot(t, layout)
end

function filled_ternary()
    function make_ax(title)
        attr(
            title=title,
            ticksuffix="%",
            min=0.01,
            linewidth=2,
            ticks="outside",
            ticklen=8,
            showgrid=true
        )
    end

    fn = tempname()
    download("https://gist.githubusercontent.com/davenquinn/988167471993bc2ece29/raw/f38d9cb3dd86e315e237fde5d65e185c39c931c2/data.json", fn)
    raw_data = JSON.parsefile(fn)
    rm(fn)

    colors = [
        "#8dd3c7",
        "#ffffb3",
        "#bebada",
        "#fb8072",
        "#80b1d3",
        "#fdb462",
        "#b3de69",
        "#fccde5",
        "#d9d9d9",
        "#bc80bd",
        "#ccebc5",
        "#ffed6f"
    ]

    traces = Array{GenericTrace,1}(undef, length(raw_data))
    for (i, (k, v)) in enumerate(raw_data)
        traces[i] = scatterternary(mode="lines", name=k,
            a=[_x["clay"] for _x in v],
            b=[_x["sand"] for _x in v],
            c=[_x["silt"] for _x in v],
            line_color="#444",
            fill="toself",
            fillcolor=colors[i],
            hoveron="fills+points"
        )
    end
    layout = Layout(
        ternary=attr(
            sum=100,
            aaxis=make_ax("Clay"),
            baxis=make_ax("Sand"),
            caxis=make_ax("Slit")),
        showlegend=false,
        width=700,
        annotations=[attr(
            showarrow=false, x=0.15, y=1.1, text="Soil types fill plot"
        )]
    )
    plot(traces, layout)
end