text:
function getAdjointBasis(Lstar, λs, eigsolver; nev = 3, verbose = false)
	# we compute the eigen-elements of the adjoint of L
	λstar, evstar = eigsolver(Lstar, nev)
	verbose && Base.display(λstar)
	# vectors to hold eigen-elements for the adjoint of L
	λstars = Vector{eltype(λs)}()
	# TODO This is a horrible hack to get  the type of the left eigenvectors
	ζstars = Vector{typeof(geteigenvector(eigsolver, evstar, 1))}()

	for (idvp, λ) in enumerate(λs)
		I = argmin(abs.(λstar .- λ))
		abs(real(λstar[I])) > 1e-2 && @warn "Did not converge to the requested eigenvalues. We found $(real(λstar[I])) !≈ 0. This might not lead to precise normal form computation. You can perhaps increase the argument `nev`."
		verbose && println("--> VP[$idvp] paired with VPstar[$I]")
		ζstar = geteigenvector(eigsolver, evstar, I)
		push!(ζstars, copy(ζstar))
		push!(λstars, λstar[I])
		# we change λstar so that it is not used twice
		λstar[I] = 1e9
	end
	return ζstars, λstars
end

"""
$(SIGNATURES)

Return a left eigenvector for an eigenvalue closest to λ. `nev` indicates how many eigenvalues must be computed by the eigensolver. Indeed, for iterative solvers, it may be needed to compute more eigenvalues than necessary.
"""
function getAdjointBasis(Lstar, λ::Number, eigsolver; nev = 3, verbose = false)
	λstar, evstar = eigsolver(Lstar, nev)
	I = argmin(abs.(λstar .- λ))
	verbose && (println("--> left eigenvalues = "); display(λstar))
	verbose && println("--> right eigenvalue = ", λ, ", left eigenvalue = ", λstar[I])
	abs(real(λstar[I])) > 1e-2 && @warn "The bifurcating eigenvalue is not that close to Re = 0. We found $(real(λstar[I])) !≈ 0.  You can perhaps increase the argument `nev`."
	ζstar = geteigenvector(eigsolver ,evstar, I)
	return copy(ζstar), λstar[I]
end

"""
$(SIGNATURES)

Compute a normal form based on Golubitsky, Martin, David G Schaeffer, and Ian Stewart. Singularities and Groups in Bifurcation Theory. New York: Springer-Verlag, 1985, VI.1.d page 295.
"""
function computeNormalForm1d(F, dF, d2F, d3F, br::ContResult, ind_bif::Int; δ = 1e-8, nev = length(eigenvalsfrombif(br, ind_bif)), Jᵗ = nothing, verbose = false, lens = br.lens, issymmetric = false, Teigvec = vectortype(br), tolFold = 1e-3, scaleζ = norm)
	bifpt = br.specialpoint[ind_bif]
	@assert bifpt.type == :bp "The provided index does not refer to a Branch Point with 1d kernel. The type of the bifurcation is $(bifpt.type). The bifurcation point is $bifpt."
	@assert abs(bifpt.δ[1]) == 1 "We only provide normal form computation for simple bifurcation points e.g when the kernel of the jacobian is 1d. Here, the dimension of the kernel is $(abs(bifpt.δ[1]))."

	verbose && println("#"^53*"\n--> Normal form Computation for 1d kernel")
	verbose && println("--> analyse bifurcation at p = ", bifpt.param)

	# Newton parameters
	options = br.contparams.newtonOptions

	# bifurcation point
	# we need this conversion when running on GPU and loading the branch from the disk
	x0 = convert(Teigvec, bifpt.x)
	p = bifpt.param

	# parameter for vector field
	parbif = set(br.params, lens, p)

	# jacobian at bifurcation point
	L = dF(x0, parbif)

	# linear solver
	ls = options.linsolver

	# "zero" eigenvalue at bifurcation point, it must be real
	λ = real(br.eig[bifpt.idx].eigenvals[bifpt.ind_ev])
	verbose && println("--> smallest eigenvalue at bifurcation = ", λ)

	# corresponding eigenvector, it must be real
	if haseigenvector(br) == false
		# we recompute the eigen-elements if there were not saved during the computation of the branch
		@info "Eigen-elements not saved in the branch. Recomputing them..."
		_λ, _ev, _ = options.eigsolver(L, bifpt.ind_ev + 2)
		@assert _λ[bifpt.ind_ev] ≈ λ "We did not find the correct eigenvalue $λ. We found $(_λ)"
		ζ = real.(geteigenvector(options.eigsolver, _ev, bifpt.ind_ev))
	else
		ζ = real.(geteigenvector(options.eigsolver, br.eig[bifpt.idx].eigenvec, bifpt.ind_ev))
	end
	ζ ./= scaleζ(ζ)

	# extract eigen-elements for adjoint(L), needed to build spectral projector
	if issymmetric
		λstar = br.eig[bifpt.idx].eigenvals[bifpt.ind_ev]
		ζstar = copy(ζ)
	else
		_Jt = isnothing(Jᵗ) ? adjoint(L) : Jᵗ(x0, parbif)
		ζstar, λstar = getAdjointBasis(_Jt, conj(λ), options.eigsolver; nev = nev, verbose = verbose)
	end

	ζstar = real.(ζstar); λstar = real.(λstar)

	@assert abs(dot(ζ, ζstar)) > 1e-10 "We got ζ⋅ζstar = $((dot(ζ, ζstar))). This dot product should not be zero. Perhaps, you can increase `nev` which is currently $nev."
	ζstar ./= dot(ζ, ζstar)

	# differentials and projector on Range(L), there are real valued
	R2 = ((dx1, dx2)      -> d2F(x0, parbif, dx1, dx2))
	R3 = ((dx1, dx2, dx3) -> d3F(x0, parbif, dx1, dx2, dx3))
	E = x -> x .- dot(x, ζstar) .* ζ

	# we compute the reduced equation: a⋅(p - pbif) + x⋅(b1⋅(p - pbif) + b2⋅x/2 + b3⋅x^2/6)
	# coefficient of p
	R01 = (F(x0, set(parbif, lens, p + δ)) .- F(x0, set(parbif, lens, p - δ))) ./ (2δ)
	a = dot(R01, ζstar)
	verbose && println("--> a = ", a)

	# coefficient of x*p
	R11 = (apply(dF(x0, set(parbif, lens, p + δ)), ζ) - apply(dF(x0, set(parbif, lens, p - δ)), ζ)) ./ (2δ)
	Ψ01, _ = ls(L, E(R01))

	b1 = dot(R11 .- R2(ζ, Ψ01), ζstar)
	verbose && println("--> b1 = ", b1)

	# coefficient of x^2
	b2v = R2(ζ, ζ)
	b2 = dot(b2v, ζstar)
	verbose && println("--> b2/2 = ", b2/2)

	# coefficient of x^3, recall b2v = R2(ζ, ζ)
	wst, _ = ls(L, E(b2v))
	b3v = R3(ζ, ζ, ζ) .- 3 .* R2(ζ, wst)
	b3 = dot(b3v, ζstar)
	verbose && println("--> b3/6 = ", b3/6)

	bp = (x0, p, parbif, lens, ζ, ζstar, (a=a, b1=b1, b2=b2, b3=b3), :NA)
	if abs(a) < tolFold
		return 100abs(b2/2) < abs(b3/6) ? Pitchfork(bp[1:end-1]...) : Transcritical(bp...)
	else
		return Fold(bp...)
	end
	# we should never hit this
	return nothing
end

computeNormalForm1d(F, dF, d2F, d3F, br::Branch, ind_bif::Int; kwargs...) = computeNormalForm1d(F, dF, d2F, d3F, getContResult(br), ind_bif; kwargs...)

"""
$(SIGNATURES)

This function provides prediction for what the zeros of the Transcritical bifurcation point.

# Arguments
- `bp::Transcritical` the bifurcation point
- `ds` at with distance relative to the bifurcation point do you want the prediction. Can be negative. Basically the parameter is `p = bp.p + ds`

# Optional arguments
- `verbose`	display information
- `ampfactor = 1` factor multiplying prediction
"""
function predictor(bp::Transcritical, ds::T; verbose = false, ampfactor = T(1)) where T
	nf = bp.nf
	a, b1, b2, b3 = nf
	pnew = bp.p + ds
	# we solve b1 * ds + b2 * amp / 2 = 0
	amp = -2ds * b1 / b2 * ampfactor
	dsfactor = T(1)

	verbose && println("--> Prediction from Normal form, δp = $(pnew - bp.p), amp = $amp")
	return (x = bp.x0 .+ amp .* real.(bp.ζ), p = pnew, dsfactor = dsfactor)
end

"""
$(SIGNATURES)

This function provides prediction for what the zeros of the Pitchfork bifurcation point.

# Arguments
- `bp::Pitchfork` the bifurcation point
- `ds` at with distance relative to the bifurcation point do you want the prediction. Based on the criticality of the Picthfork, its sign is enforced no matter what you pass. Basically the parameter is `bp.p + abs(ds) * dsfactor` where `dsfactor = ±1` depending on the criticality.

# Optional arguments
- `verbose`	display information
- `ampfactor = 1` factor multiplying prediction
"""
function predictor(bp::Pitchfork, ds::T; verbose = false, ampfactor = T(1)) where T
	nf = bp.nf
	a, b1, b2, b3 = nf

	# we need to find the type, supercritical or subcritical
	dsfactor = b1 * b3 < 0 ? T(1) : T(-1)
	if 1==1
		# we solve b1 * ds + b3 * amp^2 / 6 = 0
		amp = ampfactor * sqrt(-6abs(ds) * dsfactor * b1 / b3)
		pnew = bp.p + abs(ds) * dsfactor
	# else
	# 	# we solve b1 * ds + b3 * amp^2 / 6 = 0
	# 	amp = ampfactor * abs(ds)
	# 	pnew = bp.p + dsfactor * ds^2 * abs(b3/b1/6)
	end
	verbose && println("--> Prediction from Normal form, δp = $(pnew - bp.p), amp = $amp")
	return (x = bp.x0 .+ amp .* real.(bp.ζ), p = pnew, dsfactor = dsfactor)
end

function predictor(bp::Fold, ds::T; verbose = false, ampfactor = T(1)) where T
	@info "It seems the point is a Saddle-Node bifurcation. The normal form is $(bp.nf)."
	return nothing
end
####################################################################################################
# type for bifurcation point Nd kernel for the jacobian
function factor3d(i,j,k)
	if i == j == k
		return 1/6
	else
		_power = length(unique((i,j,k)))
		if _power == 1
			factor = 1/6 /2
		elseif _power == 2
			factor = 1/2 / 3
		else
			factor = 1.0
		end
		return factor
	end
end

function (bp::NdBranchPoint)(::Val{:reducedForm}, x, p::T) where T
	# formula from https://fr.qwe.wiki/wiki/Taylor's_theorem
	# dimension of the kernel
	N = length(bp.ζ)
	# for the output
	out = zero(x)
	# normal form
	nf = bp.nf
	# coefficient p
	out .= p .* nf.a
	# factor to account for factorials
	factor = T(1)

	@inbounds for ii in 1:N
		factor = T(1)
		out[ii] = 0
		# coefficient x*p
		for jj in 1:N
			# coefficient x*p
			out[ii] += p * nf.b1[ii , jj] * x[jj]

			for kk in 1:N
				# coefficients of x^2
				factor = jj == kk ? 1/2 : 1
				out[ii] += nf.b2[ii, jj, kk] * x[jj] * x[kk] * factor / 2

				for ll in 1:N
					# coefficients of x^3
					factor = factor3d(ii, jj, kk)
					out[ii] += nf.b3[ii, jj, kk, ll] * x[jj] * x[kk]  * x[ll] * factor
				end
			end
		end
	end
	return out
end

function (bp::NdBranchPoint)(x, δp::Real)
	out = bp.x0 .+ x[1] .* bp.ζ[1]
	for ii in 2:length(x)
		out .+= x[ii] .* bp.ζ[ii]
	end
	return out
end

"""
$(SIGNATURES)

Print the normal form `bp` with a nice string.
"""
function nf(bp::NdBranchPoint; tol = 1e-6, digits = 4)
	superDigits = [c for c in "⁰ ²³⁴⁵⁶⁷⁸⁹"]

	nf = bp.nf
	N = length(nf.a)
	out = ["" for _ in 1:N]

	for ii = 1:N
		if abs(nf.a[ii]) > tol
			out[ii] *= "$(round(nf.a[ii],digits=digits)) ⋅ p"
		end
		for jj in 1:N
			coeff = round(nf.b1[ii,jj],digits=digits)
			if abs(coeff) > tol
				out[ii] *= " + $coeff * x$jj ⋅ p"
			end
		end

		for jj in 1:N
			for kk in jj:N
				coeff = round(nf.b2[ii,jj,kk] / 2,digits=digits)
				if abs(coeff) > tol
					if jj == kk
						out[ii] *= " + $coeff ⋅ x$(jj)²"
					else
						out[ii] *= " + $(round(2coeff,digits=digits)) ⋅ x$jj ⋅ x$kk"
					end
				end

				for ll in kk:N
					coeff = round(nf.b3[ii,jj,kk,ll] / 6,digits=digits)
					_pow = zeros(Int64,N)
					_pow[jj] += 1;_pow[kk] += 1;_pow[ll] += 1;

					if abs(coeff) > tol
						if jj == kk == ll
							out[ii] *= " + $coeff"
						else
							out[ii] *= " + $(round(3coeff,digits=digits))"
						end
						for mm in 1:N
							if _pow[mm] > 1
								out[ii] *= " ⋅ x$mm" * (superDigits[_pow[mm]+1])
							elseif _pow[mm] == 1
								out[ii] *= " ⋅ x$mm"
							end
						end
					end
				end
			end
		end
	end
	return out
end

"""
$(SIGNATURES)

Bi-orthogonalise the arguments.
"""
function biorthogonalise(ζs, ζstars, verbose)
	# change only the ζstars to have bi-orthogonal left/right eigenvectors
	# we could use projector P=A(A^{T}A)^{-1}A^{T}
	# we use Gram-Schmidt algorithm instead
	G = [ dot(ζ, ζstar) for ζ in ζs, ζstar in ζstars]
	@assert abs(det(G)) >1e-14 "The Gram matrix is not invertible! det(G) = $(det(G)), G = \n$G $(display(G))"

	# save those in case the first algo fails
	_ζs = deepcopy(ζs)
	_ζstars = deepcopy(ζstars)

	# first algo
	tmp = copy(ζstars[1])
	for ii in eachindex(ζstars)
		tmp .= ζstars[ii]
		for jj in eachindex(ζs)
			if ii != jj
				tmp .-= dot(tmp, ζs[jj]) .* ζs[jj] ./ dot(ζs[jj], ζs[jj])
			end
		end
		ζstars[ii] .= tmp ./ dot(tmp, ζs[ii])
	end

	G = [ dot(ζ, ζstar) for ζ in ζs, ζstar in ζstars]

	# we switch to another algo if the above fails
	if norm(G - LinearAlgebra.I, Inf) >= 1e-5
		@warn "Gram matrix not equal to idendity. Switching to LU algorithm."
		println("G (det = $(det(G))) = "); display(G)
		G = [ dot(ζ, ζstar) for ζ in _ζs, ζstar in _ζstars]
		_F = lu(G; check = true)
		display(inv(_F.L) * inv(_F.P) * G * inv(_F.U))
		ζs = inv(_F.L) * inv(_F.P) * _ζs
		ζstars = inv(_F.U)' * _ζstars
	end

	# test the bi-orthogonalization
	G = [ dot(ζ, ζstar) for ζ in ζs, ζstar in ζstars]
	verbose && (printstyled(color=:green, "--> Gram matrix = \n");Base.display(G))
	@assert norm(G - LinearAlgebra.I, Inf) < 1e-5 "Failure in bi-orthogonalisation of the right / left eigenvectors. The left eigenvectors do not form a basis. You may want to increase `nev`, G = \n $(display(G))"
	return ζs, ζstars
end

"""
$(SIGNATURES)

Compute the normal form of the bifurcation point located at `br.specialpoint[ind_bif]`.

# Arguments
- `F, dF, d2F, d3F` vector field `(x, p) -> F(x, p)` and its derivatives w.r.t. `x`.
- `br` result from a call to [`continuation`](@ref)
- `ind_bif` index of the bifurcation point in `br.specialpoint`

# Optional arguments
- `δ` used to compute ∂pF with finite differences
- `nev` number of eigenvalues used to compute the spectral projection. This number has to be adjusted when used with iterative methods.
- `Jᵗ = (x,p) -> ...` jacobian adjoint, it should be implemented in an efficient manner. For matrix-free methods, `transpose` is not readily available and the user must provide a dedicated method. In the case of sparse based jacobian, `Jᵗ` should not be passed as it is computed internally more efficiently, i.e. it avoids recomputing the jacobian as it would be if you pass `Jᵗ = (x, p) -> transpose(dF(x, p))`.
- `verbose` whether to display information
- `ζs` list of vectors spanning the kernel of `dF` at the bifurcation point. Useful to enforce the basis for the normal form.
- `lens::Lens` specify which parameter to take the partial derivative ∂pF
- `issymmetric` whether the Jacobian is Symmetric, avoid computing the left eigenvectors.
- `scaleζ` function to normalise the kernel basis. Indeed, when used with large vectors and `norm`, it results in ζs and the normal form coefficient being super small.

Based on Golubitsky, Martin, David G Schaeffer, and Ian Stewart. Singularities and Groups in Bifurcation Theory. New York: Springer-Verlag, 1985, VI.1.d page 295.
"""
function computeNormalForm(F, dF, d2F, d3F,
			br::ContResult, id_bif::Int ;
			δ = 1e-8,
			nev = length(eigenvalsfrombif(br, id_bif)),
			Jᵗ = nothing,
			verbose = false,
			ζs = nothing,
			lens = br.lens,
			issymmetric = false,
			Teigvec = getvectortype(br),
			scaleζ = norm)
	bifpt = br.specialpoint[id_bif]
	if abs(bifpt.δ[2]) > 0 # we try a Hopf point
		return hopfNormalForm(F, dF, d2F, d3F, br, id_bif; δ = δ, nev = nev, Jᵗ = Jᵗ, verbose = verbose, lens = lens, Teigvec = Teigvec, scaleζ = scaleζ)
	elseif abs(bifpt.δ[1]) == 1 # simple branch point
		return computeNormalForm1d(F, dF, d2F, d3F, br, id_bif ; δ = δ, nev = nev, Jᵗ = Jᵗ, verbose = verbose, lens = lens, issymmetric = issymmetric, Teigvec = Teigvec, scaleζ = scaleζ)
	end
	# kernel dimension:
	N = abs(bifpt.δ[1])
	# in case nev = 0 (number of unstable eigenvalues), we increase nev to avoid bug
	nev = max(2N, nev)
	verbose && println("#"^53*"\n--> Normal form Computation for a $N-d kernel")
	verbose && println("--> analyse bifurcation at p = ", bifpt.param)

	# Newton parameters
	options = br.contparams.newtonOptions

	# bifurcation point
	x0 = convert(Teigvec, bifpt.x)
	p = bifpt.param

	# parameter for vector field
	parbif = set(br.params, br.lens, p)

	# jacobian at bifurcation point
	L = dF(x0, parbif)

	# we invert L repeatdly, so we try to factorize it
	Linv = L isa AbstractMatrix ? factorize(L) : L

	# linear solver
	ls = options.linsolver

	# "zero" eigenvalues at bifurcation point
	rightEv = br.eig[bifpt.idx].eigenvals
	indev = br.specialpoint[id_bif].ind_ev
	λs = rightEv[indev-N+1:indev]
	verbose && println("--> smallest eigenvalues at bifurcation = ", real.(λs))

	# and corresponding eigenvectors
	if isnothing(ζs) # do we have a basis for the kernel?
		if haseigenvector(br) == false # are the eigenvector saved in the branch?
			@info "No eigenvector recorded, computing them on the fly"
			# we recompute the eigen-elements if there were not saved during the computation of the branch
			_λ, _ev, _ = options.eigsolver(L, length(rightEv))
			verbose && (println("--> (λs, λs (recomputed)) = "); display(hcat(rightEv, _λ[1:length(rightEv)])))
			if norm(_λ[1:length(rightEv)] - rightEv, Inf) > br.contparams.precisionStability
			@warn "We did not find the correct eigenvalues (see 1st col). We found the eigenvalues displayed in the second column:\n $(display(hcat(rightEv, _λ[1:length(rightEv)]))).\n Difference between the eigenvalues:" display(_λ[1:length(rightEv)] - rightEv)
			end
			ζs = [copy(geteigenvector(options.eigsolver, _ev, ii)) for ii in indev-N+1:indev]
		else
			ζs = [copy(geteigenvector(options.eigsolver, br.eig[bifpt.idx].eigenvec, ii)) for ii in indev-N+1:indev]
		end
	end

	# extract eigen-elements for transpose(L), needed to build spectral projector
	# it is OK to re-scale at this stage as the basis ζs is not touched anymore, we
	# only adjust ζstars
	for ζ in ζs; ζ ./= scaleζ(ζ); end
	if issymmetric
		λstars = copy(λs)
		ζstars = copy.(ζs)
	else
		_Jt = isnothing(Jᵗ) ? transpose(L) : Jᵗ(x0, parbif)
		ζstars, λstars = getAdjointBasis(_Jt, conj.(λs), options.eigsolver; nev = nev, verbose = verbose)
	end
	ζstars = real.(ζstars); λstars = real.(λstars)
	ζs = real.(ζs); λs = real.(λs)
	verbose && println("--> VP     = ", λs, "\n--> VPstar = ", λstars)

	ζs, ζstars = biorthogonalise(ζs, ζstars, verbose)

	# differentials should work as we are looking at reals
	R2 = ((dx1, dx2)      -> d2F(x0, parbif, dx1, dx2))
	R3 = ((dx1, dx2, dx3) -> d3F(x0, parbif, dx1, dx2, dx3))

	# projector on Range(L)
	function E(x)
		out = copy(x)
		for ii in 1:N
			out .= out .- dot(x, ζstars[ii]) .* ζs[ii]
		end
		out
	end

	# vector eltype
	Tvec = eltype(ζs[1])

	# coefficients of p
	dgidp = Vector{Tvec}(undef, N)
	R01 = (F(x0, set(parbif, lens, p + δ)) .- F(x0, set(parbif, lens, p - δ))) ./ (2δ)
	for ii in 1:N
		dgidp[ii] = dot(R01, ζstars[ii])
	end
	verbose && printstyled(color=:green,"--> a (∂/∂p) = ", dgidp, "\n")

	# coefficients of x*p
	d2gidxjdpk = zeros(Tvec, N, N)
	for ii in 1:N, jj in 1:N
		R11 = (apply(dF(x0, set(parbif, lens, p + δ)), ζs[jj]) .- apply(dF(x0, set(parbif, lens, p - δ)), ζs[jj])) ./ (2δ)
		Ψ01, flag = ls(Linv, E(R01))
		~flag && @warn "linear solver did not converge"
		d2gidxjdpk[ii,jj] = dot(R11 .- R2(ζs[jj], Ψ01), ζstars[ii])
	end
	verbose && (printstyled(color=:green, "\n--> b1 (∂²/∂x∂p)  = \n"); Base.display( d2gidxjdpk ))

	# coefficients of x^2
	d2gidxjdxk = zeros(Tvec, N, N, N)
	for ii in 1:N, jj in 1:N, kk in 1:N
		b2v = R2(ζs[jj], ζs[kk])
		d2gidxjdxk[ii, jj, kk] = dot(b2v, ζstars[ii])
	end

	if verbose
		printstyled(color=:green, "\n--> b2 (∂²/∂x²) = \n")
		for ii in 1:N
			printstyled(color=:blue, "--> component $ii\n")
			Base.display( d2gidxjdxk[ii,:,:] ./ 2)
		end
	end

	# coefficient of x^3
	d3gidxjdxkdxl = zeros(Tvec, N, N, N, N)
	for jj in 1:N, kk in 1:N, ll in 1:N
		b3v = R3(ζs[jj], ζs[kk], ζs[ll])
		# d3gidxjdxkdxl[ii,jj,kk,ll] = dot(b3v, ζstars[ii])

		wst, flag = ls(Linv, E(R2(ζs[ll], ζs[kk])))
		~flag && @warn "linear solver did not converge"
		b3v .-= R2(ζs[jj], wst)
		# d3gidxjdxkdxl[ii,jj,kk,ll] -= dot(R2(ζs[jj], wst), ζstars[ii])

		wst, flag = ls(Linv, E(R2(ζs[ll], ζs[jj])))
		~flag && @warn "linear solver did not converge"
		b3v .-= R2(ζs[kk], wst)
		# d3gidxjdxkdxl[ii,jj,kk,ll] -= dot(R2(ζs[kk], wst), ζstars[ii])

		wst, flag = ls(Linv, E(R2(ζs[kk], ζs[jj])))
		~flag && @warn "linear solver did not converge"
		b3v .-= R2(ζs[ll], wst)
		# d3gidxjdxkdxl[ii,jj,kk,ll] -= dot(R2(ζs[ll], wst), ζstars[ii])
		for ii in 1:N
			d3gidxjdxkdxl[ii, jj, kk, ll] = dot(b3v, ζstars[ii])
		end
	end
	if verbose
		printstyled(color=:green, "\n--> b3 (∂³/∂x³) = \n")
		for ii in 1:N
			printstyled(color=:blue, "--> component $ii\n")
			Base.display( d3gidxjdxkdxl[ii,:,:, :] ./ 6 )
		end
	end

	return NdBranchPoint(x0, p, parbif, lens, ζs, ζstars, (a=dgidp, b1=d2gidxjdpk, b2=d2gidxjdxk, b3=d3gidxjdxkdxl), Symbol("$N-d"))

end

computeNormalForm(F, dF, d2F, d3F, br::Branch, id_bif::Int; kwargs...) = computeNormalForm(F, dF, d2F, d3F, getContResult(br), id_bif; kwargs...)

"""
$(SIGNATURES)

This function provides prediction for what the zeros of the reduced equation / normal form should be. The algorithm to find these zeros is based on deflated newton.
"""
function predictor(bp::NdBranchPoint, δp::T;
		verbose::Bool = false,
		ampfactor = T(1),
		nbfailures = 30,
		maxiter = 100,
		perturb = identity,
		J = nothing,
		normN = x -> norm(x, Inf),
		optn = NewtonPar(maxIter = maxiter, verbose = verbose)) where T
	# dimension of the kernel
	n = length(bp.ζ)

	# find zeros for the normal on each side of the bifurcation point
	function getRootsNf(_ds)
		deflationOp = DeflationOperator(2.0, dot, 1.0, [zeros(n)])
		failures = 0
		# we allow for 10 failures of nonlinear deflation
		outdef1 = rand(n)
		while failures < nbfailures
			if isnothing(J)
				jac = (x,p) -> ForwardDiff.jacobian( z->perturb(bp(Val(:reducedForm), z, p)), x)
				outdef1, hist, flag, _ = newton((x, p) -> perturb(bp(Val(:reducedForm), x, p)), jac, outdef1 .+ 0.1rand(n), _ds, optn, deflationOp; normN = normN)
			else
				outdef1, hist, flag, _ = newton((x, p) -> perturb(bp(Val(:reducedForm), x, p)), J, outdef1 .+ 0.1rand(n), _ds, optn, deflationOp; normN = normN)
			end
			flag && push!(deflationOp, ampfactor .* outdef1)
			~flag && (failures += 1)
		end
		return deflationOp.roots
	end
	rootsNFm =  getRootsNf(-abs(δp))
	rootsNFp =  getRootsNf(abs(δp))
	println("\n--> BS from Non simple branch point")
	printstyled(color=:green, "--> we find $(length(rootsNFm)) (resp. $(length(rootsNFp))) roots on the left (resp. right) of the bifurcation point counting the trivial solution (Reduced equation).\n")
	return (before = rootsNFm, after = rootsNFp)
end
####################################################################################################
"""
$(SIGNATURES)

Compute the Hopf normal form.

# Arguments
- `F, dF, d2F, d3F`: function `(x,p) -> F(x,p)` and its differentials `(x,p,dx) -> d1F(x,p,dx)`, `(x,p,dx1,dx2) -> d2F(x,p,dx1,dx2)`...
- `pt::Hopf` Hopf bifurcation point
- `ls` linear solver

# Optional arguments
- `δ = 1e-8` used for finite differences
- `verbose` bool to print information
"""
function hopfNormalForm(F, dF, d2F, d3F, pt::Hopf, ls; δ = 1e-8, verbose::Bool = false)
	x0 = pt.x0
	p = pt.p
	lens = pt.lens
	parbif = set(pt.params, lens, p)
	ω = pt.ω
	ζ = pt.ζ
	cζ = conj.(pt.ζ)
	ζstar = pt.ζstar

	# jacobian at the bifurcation point
	L = dF(x0, parbif)

	# we use BilinearMap to be able to call on complex valued arrays
	R2 = BilinearMap( (dx1, dx2)      -> d2F(x0, parbif, dx1, dx2) ./2)
	R3 = TrilinearMap((dx1, dx2, dx3) -> d3F(x0, parbif, dx1, dx2, dx3) ./6 )

	# −LΨ001 = R01
	R01 = (F(x0, set(parbif, lens, p + δ)) .- F(x0, set(parbif, lens, p - δ))) ./ (2δ)
	Ψ001, _ = ls(L, -R01)

	# (2iω−L)Ψ200 = R20(ζ,ζ)
	R20 = R2(ζ, ζ)
	Ψ200, _ = ls(L, R20; a₀ = Complex(0, 2ω), a₁ = -1)
	# @assert Ψ200 ≈ (Complex(0, 2ω)*I - L) \ R20

	# −LΨ110 = 2R20(ζ,cζ).
	R20 = 2 .* R2(ζ, cζ)
	Ψ110, _ = ls(L, -R20)

	# a = ⟨R11(ζ) + 2R20(ζ,Ψ001),ζ∗⟩
	av = (apply(dF(x0, set(parbif, lens, p + δ)), ζ) .- apply(dF(x0, set(parbif, lens, p - δ)), ζ)) ./ (2δ)
	av .+= 2 .* R2(ζ, Ψ001)
	a = dot(av, ζstar)

	# b = ⟨2R20(ζ,Ψ110) + 2R20(cζ,Ψ200) + 3R30(ζ,ζ,cζ), ζ∗⟩)
	bv = 2 .* R2(ζ, Ψ110) .+ 2 .* R2(cζ, Ψ200) .+ 3 .* R3(ζ, ζ, cζ)
	b = dot(bv, ζstar)

	# return coefficients of the normal form
	verbose && println((a = a, b = b))
	pt.nf = (a = a, b = b)
	if real(a) * real(b) < 0
		pt.type = :SuperCritical
	elseif real(a) * real(b) > 0
		pt.type = :SubCritical
	else
		pt.type = :Singular
	end
	verbose && printstyled(color = :red,"--> Hopf bifurcation point is: ", pt.type, "\n")
	return pt
end

"""
$(SIGNATURES)

Compute the Hopf normal form.

# Arguments
- `F, dF, d2F, d3F`: function `(x,p) -> F(x,p)` and its differentials `(x,p,dx) -> d1F(x,p,dx)`, `(x,p,dx1,dx2) -> d2F(x,p,dx1,dx2)`...
- `br` branch result from a call to [`continuation`](@ref)
- `ind_hopf` index of the bifurcation point in `br`
- `options` options for the Newton solver

# Optional arguments
- `Jᵗ` is the jacobian adjoint, used for computation of the eigen-elements of the jacobian adjoint, needed to compute the spectral projector
- `δ = 1e-8` used for finite differences
- `nev = 5` number of eigenvalues to compute to estimate the spectral projector
- `verbose` bool to print information
"""
function hopfNormalForm(F, dF, d2F, d3F, br::AbstractBranchResult, ind_hopf::Int; Jᵗ = nothing, δ = 1e-8, nev = length(eigenvalsfrombif(br, id_bif)), verbose::Bool = false, lens = br.lens, Teigvec = getvectortype(br), scaleζ = norm)
	@assert br.specialpoint[ind_hopf].type == :hopf "The provided index does not refer to a Hopf Point"
	verbose && println("#"^53*"\n--> Hopf Normal form computation")

	# Newton parameters
	options = br.contparams.newtonOptions

	# bifurcation point
	bifpt = br.specialpoint[ind_hopf]
	eigRes = br.eig

	# eigenvalue
	λ = eigRes[bifpt.idx].eigenvals[bifpt.ind_ev]
	ω = imag(λ)

	# parameter for vector field
	p = bifpt.param
	parbif = set(br.params, lens, p)

	# jacobian at bifurcation point
	L = dF(convert(Teigvec, bifpt.x), parbif)

	# right eigenvector
	if haseigenvector(br) == false
		# we recompute the eigen-elements if there were not saved during the computation of the branch
		_λ, _ev, _ = options.eigsolver(L, bifpt.ind_ev + 2)
		@assert _λ[bifpt.ind_ev] ≈ λ "We did not find the correct eigenvalue $λ. We found $(_λ)"
		ζ = geteigenvector(options.eigsolver, _ev, bifpt.ind_ev)
	else
		ζ = copy(geteigenvector(options.eigsolver ,br.eig[bifpt.idx].eigenvec, bifpt.ind_ev))
	end
	ζ ./= scaleζ(ζ)

	# left eigen-elements
	_Jt = isnothing(Jᵗ) ? adjoint(L) : Jᵗ(x, p)
	ζstar, λstar = getAdjointBasis(_Jt, conj(λ), options.eigsolver; nev = nev, verbose = verbose)

	# check that λstar ≈ conj(λ)
	abs(λ + λstar) > 1e-2 && @warn "We did not find the left eigenvalue for the Hopf point to be very close to the imaginary part, $λ ≈ $(λstar) and $(abs(λ + λstar)) ≈ 0?\n You can perhaps increase the number of computed eigenvalues, the number is nev = $nev"

	# normalise left eigenvector
	ζstar ./= dot(ζ, ζstar)
	@assert dot(ζ, ζstar) ≈ 1

	hopfpt = Hopf(
		bifpt.x,
		bifpt.param,
		ω,
		parbif,
		lens,
		ζ,
		ζstar,
		(a = 0. + 0im, b = 0. + 0im),
		:SuperCritical
	)
	return hopfNormalForm(F, dF, d2F, d3F, hopfpt, options.linsolver ; δ = δ, verbose = verbose)
end

function predictor(hp::Hopf, ds::T; verbose = false, ampfactor = T(1) ) where T
	# get the normal form
	nf = hp.nf
	a = nf.a
	b = nf.b

	# we need to find the type, supercritical or subcritical
	dsfactor = real(a) * real(b) < 0 ? T(1) : T(-1)
	pnew = hp.p + abs(ds) * dsfactor

	# we solve a * ds + b * amp^2 = 0
	amp = ampfactor * sqrt(-abs(ds) * dsfactor * real(a) / real(b))

	# o(1) correction to Hopf Frequency
	ω = hp.ω + (imag(a) - imag(b) * real(a) / real(b)) * ds

	return (orbit = t -> hp.x0 .+ 2amp .* real.(hp.ζ .* exp(complex(0, t))),
			amp = 2amp,
			ω = ω,
			p = pnew,
			dsfactor = dsfactor)
end


text:
using Gtk
params = HubertAmplifierParams()
rob = HubertAmplifier("hubertAmp", params)

if ask_dialog("You are trying to start the Hubert amplifier hardware test. Please ensure that the amp is turned on and a load is connected. Please continue only if the robot is safe to move!", "Cancel", "Start test")
    @error "Hubert amp test has not been implemented yet!"
    @test_broken 0

    @test ask_dialog("Did everything seem normal?")
else
    @test_broken 0
end

ext: jl


text:
module Tester

using Revise
using Logging
using Genie, Genie.Util, Genie.Configuration
using Reexport
@reexport using HTTP

"""
    setenv() :: Nothing

Switches Genie to the test env for the duration of the current execution.
"""
function setenv(environment = TEST) :: Nothing
  if Genie.config.tests_force_test_env
    Genie.config.app_env = environment
    @warn "Switched app to >> $(uppercase(Genie.config.app_env)) << env"

    Genie.Configuration.env() == environment || error("Could not switch env")
  end

  nothing
end

end


text:
using Test
using LinearAlgebra
using Distributed
using Logging
include("/usr/share/julia/test/testenv.jl")
addprocs_with_testenv(3)
@test nworkers() == 3

@everywhere using StochasticPrograms
using ProgressiveHedgingSolvers
using JuMP
using GLPKMathProgInterface
using OSQP

τ = 1e-2
reference_solver = GLPKSolverLP()
osqp = OSQP.OSQPMathProgBaseInterface.OSQPSolver(verbose=0)

executors = [Serial(),
             Synchronous()]

penalties = [Fixed(),
             Adaptive(θ = 1.01)]

problems = Vector{Tuple{<:StochasticProgram,String}}()
@info "Loading test problems..."
@info "Loading simple..."
include("simple.jl")
@info "Loading farmer..."
include("farmer.jl")
@info "Loading infeasible..."
include("infeasible.jl")
@info "Test problems loaded. Starting test sequence."

@testset "Distributed solver" begin
    @testset "$(solverstr(ph)): $name" for ph in [ProgressiveHedgingSolver(osqp,
                                                                           execution = executor,
                                                                           penalty = penalty,
                                                                           τ = 1e-3,
                                                                           log = false)
                                                  for executor in executors, penalty in penalties], (sp,name) in problems
        @testset "Distributed data" begin
            optimize!(sp, solver=reference_solver)
            x̄ = optimal_decision(sp)
            Q̄ = optimal_value(sp)
            with_logger(NullLogger()) do
                optimize!(sp, solver=ph)
            end
            @test abs(optimal_value(sp) - Q̄)/(1e-10+abs(Q̄)) <= τ
            @test norm(optimal_decision(sp) - x̄)/(1e-10+norm(x̄)) <= sqrt(τ)
        end
        @testset "Data on single remote node" begin
            sp_onenode = copy(sp)
            add_scenarios!(sp_onenode, scenarios(sp), workers()[1])
            optimize!(sp_onenode, solver=reference_solver)
            x̄ = optimal_decision(sp_onenode)
            Q̄ = optimal_value(sp_onenode)
            with_logger(NullLogger()) do
                optimize!(sp, solver=ph)
            end
            @test abs(optimal_value(sp_onenode) - Q̄)/(1e-10+abs(Q̄)) <= τ
            @test norm(optimal_decision(sp_onenode) - x̄)/(1e-10+norm(x̄)) <= sqrt(τ)
        end
        @testset "Local data" begin
            sp_nondist = copy(sp, procs = [1])
            add_scenarios!(sp_nondist, scenarios(sp))
            optimize!(sp_nondist, solver=reference_solver)
            x̄ = optimal_decision(sp_nondist)
            Q̄ = optimal_value(sp_nondist)
            with_logger(NullLogger()) do
                optimize!(sp, solver=ph)
            end
            @test abs(optimal_value(sp_nondist) - Q̄)/(1e-10+abs(Q̄)) <= τ
            @test norm(optimal_decision(sp_nondist) - x̄)/(1e-10+norm(x̄)) <= sqrt(τ)
        end
    end
end


text:
using Distributions
using DataStructures
using Statistics

"""
    SobolData

A struct which holds all information needed for the sampling and analysis of a
specific problem using Sobol Analysis:

`params::Union{OrderedDict{Symbol, <:Any}, Nothing} = nothing`: a dictionary mapping parameter names to their Distribution
`calc_second_order::Bool = true`: whether or not to calculate second order sensitivity indices
`N::Int = 1000`: the number of runs
"""
mutable struct SobolData
    params::Union{OrderedDict{Symbol, <:Any}, Nothing}
    calc_second_order::Bool
    N::Int

    function SobolData(;params= nothing, calc_second_order = true, N = 1000)
        return new(params, calc_second_order, N)
    end
end

"""
    DeltaData

A struct which holds all information needed for the sampling and analysis of a
specific problem using Delta Analysis:

`params::Union{OrderedDict{Symbol, <:Any}, Nothing} = nothing`: a dictionary mapping parameter names to their Distribution
`N::Int = 1000`: the number of runs
"""
mutable struct DeltaData
    params::Union{OrderedDict{Symbol, <:Any}, Nothing}
    N::Int

    function DeltaData(;params= nothing, N = 1000)
        return new(params, N)
    end
end

"""
    scale_samples!(sequence::AbstractArray{<:Number, N1}, dists::AbstractArray{T, N2})

Rescale a Sobol `sequence` of parameters from the 0-to-1 range to their corresponding
univariate distributions `dists`.
"""
function scale_samples!(sequence::AbstractArray{<:Number, N1}, dists::AbstractArray{T, N2}) where T where N1 where N2
    D = length(dists) # number of parameters
    duplicate = size(sequence,2) == 2*D # sobol requires twice as many columns as D

    for param in 1:D
        dist = dists[param]
        if duplicate
            if typeof(dist) <: UnivariateDistribution
                sequence[:, [param, param + D]] = quantile.(dist, sequence[:, [param, param + D]])
            # temporary fix to problem calling quantile on a EmpiricalDistribution from Mimi
            else
                sequence[:, [param, param + D]] = quantile(dist, sequence[:, [param, param + D]])
            end
        else
            if typeof(dist) <: UnivariateDistribution
                sequence[:, param] = quantile.(dist, sequence[:, param])
            # temporary fix to problem calling quantile on a EmpiricalDistribution from Mimi
            else
                sequence[:, param] = quantile(dist, sequence[:, param])
            end
        end
    end
end

"""
    _check_conf_flag(num_resamples::Union{Nothing, Int}, conf_level::Union{Nothing, Number})
Check to see if confdience interval should be calculated based on the provided
`num_resamples` and `conf_level`.  Error if only one has a value and the other is Nothing.
"""
function _check_conf_flag(num_resamples::Union{Nothing, Int}, conf_level::Union{Nothing, Number})
    num_nothings = (num_resamples === nothing) + (conf_level === nothing)
    if num_nothings == 1
        error("Number of resamples is $num_resamples, while confidence level is $conf_level ... either none or both must be nothing")
    elseif num_nothings == 2
        conf_flag = false
    else
        conf_flag = true
    end
    return conf_flag
end

max_stars_repo_name: lrennels/SAJulia

max_issues_repo_name: lrennels/SAJulia


text:
const Pt{T<:Real} = Union{AbstractVector{T},NTuple{<:Any,T}}
const Poly = AbstractVector{<:Union{NTuple{<:Any,<:Real},AbstractVector{<:Real}}}

const DEFAULT_POINT_ORDER = (XDim, YDim)

# _fill_geometry!
# Fill a raster with `fill` where it interacts with a geometry.
# We reduce all geometries to vectors with `coordinates` so that a simple
# vector can also represent a polygon. In future we should do this lazily.
function _fill_geometry!(B::AbstractRaster, coll::GI.AbstractFeatureCollection; order, kw...)
    if bbox_might_overlap(B, coll, order)
        foreach(coll.features) do feature
            _fill_geometry!(B, feature; order, kw...)
        end
    end
    return B
end
function _fill_geometry!(B::AbstractRaster, feature::GI.AbstractFeature; order, kw...)
    if bbox_might_overlap(B, feature, order)
        _fill_geometry!(B, GI.geometry(feature); order, kw...)
    end
    return B
end
function _fill_geometry!(B::AbstractRaster, geoms::AbstractVector{<:Union{Missing,<:GI.AbstractGeometry,<:GI.AbstractFeature}}; order, kw...)
    any(map(g -> bbox_might_overlap(B, g, order), geoms)) || return B
    if all(g -> typeof(g) == typeof(first(geoms)), geoms)
        coords = GI.coordinates.(GI.geometry.(geoms))
        _fill_geometry!(B, coords; order, shape=_geom_shape(GI.geometry(first(geoms))), kw...)
    else
        _fill_mixed_geometries!(B, geoms; order, kw...)
    end
end
function _fill_geometry!(B::AbstractRaster, geom::GI.AbstractGeometry; order, kw...)
    # Dont do anything if the bbox doesn't overlap
    bbox_might_overlap(B, geom, order) || return B
    _fill_geometry!(B, GI.coordinates(geom); order, shape=_geom_shape(geom), kw...)
end
function _fill_geometry!(B::AbstractRaster, geom::AbstractVector; shape=:polygon, order, kw...)
    gbounds = _geom_bounds(geom, order)
    abounds = bounds(dims(B, order)) # Only mask if the gemoetry bounding box overlaps the array bounding box
    bounds_overlap(gbounds, abounds) || return B
    if shape === :polygon
        _fill_polygon!(B, geom; polybounds=gbounds, order, kw...)
    elseif shape === :line
        _fill_linestring!(B, geom; order, kw...)
    elseif shape === :point
        _fill_points!(B, geom; order, kw...)
    else
        _shape_error(shape)
    end
    return B
end

# Multiple geometry types. Rasterize by groups of types separately
@noinline function _fill_mixed_geometries!(B, geoms; order, kw...)
    geomgroups = Dict{Type,Vector}()
    for geom in geoms
        key = typeof(geom)
        if haskey(geomgroups, key)
            push!(geomgroups[key], geom)
        else
            geomgroups[key] = [geom]
        end
    end
    for group in values(geomgroups)
        _fill_geometry!(B, GI.coordinates.(group); order, shape=_geom_shape(first(group)), kw...)
    end
end

# _fill_polygon!
# Fill a raster with `fill` where pixels are inside a polygon
# `boundary` determines how edges are handled
function _fill_polygon!(B::AbstractRaster, poly; polybounds, order, fill=true, boundary=:center, kw...)
    # TODO take a view of B for the polygon
    # We need a tuple of all the dims in `order`
    # We also need the index locus to be the center so we are
    # only selecting cells more than half inside the polygon
    shifted_dims = map(dims(B, order)) do d
        d = DD.maybeshiftlocus(Center(), d)
        # DimPoints are faster if we materialise the dims
        modify(Array, d)
    end
    pts = DimPoints(shifted_dims)
    inpoly = if any(map(d -> DD.order(d) isa Unordered, shifted_dims))
        inpolygon(vec(pts), poly)
    else
        pointbounds = map(shifted_dims) do d
            f = first(d)
            l = last(d)
            f < l ? (f, l) : (l, f)
        end
        # Precalculate vmin, vmax and iyperm permutation vector
        # This is much faster than calling `sortperm` in PolygonInbounds.jl
        vmin = [first.(pointbounds)...]'
        vmax = [last.(pointbounds)...]'
        pmin = [first.(polybounds)...]'
        pmax = [last.(polybounds)...]'
        pts = DimPoints(shifted_dims)
        iyperm = _iyperm(shifted_dims)
        inpolygon(vec(pts), poly; vmin, vmax, pmin, pmax, iyperm)
    end
    inpolydims = dims(B, order)
    reshaped = Raster(reshape(inpoly, size(inpolydims)), inpolydims)
    return _inner_fill_polygon!(B, poly, inpoly, reshaped; order, fill, boundary)
end

# split to make a type stability function barrier
function _inner_fill_polygon!(B::AbstractRaster, poly, inpoly, reshaped; order, fill=true, boundary=:center, kw...)
    # Get the array as points
    # Use the first column of the output - the points in the polygon,
    # and reshape to match `A`
    for D in DimIndices(B)
        @inbounds if reshaped[D...]
            @inbounds B[D...] = fill
        end
    end
    if boundary === :touches
        # Add the line pixels
        _fill_linestring!(B, poly; order, fill)
    elseif boundary === :inside
        # Remove the line pixels
        _fill_linestring!(B, poly; order, fill=!fill)
    elseif boundary !== :center
        _boundary_error(boundary)
    end
    return B
end

function _iyperm(dims::Tuple{<:Dimension,<:Dimension})
    a1, a2 = map(dims) do d
        l = parent(d)
        LA.ordered_firstindex(l):_order_step(l):LA.ordered_lastindex(l)
    end
    iyperm = Array{Int}(undef, length(a1) * length(a2))
    lis = (LinearIndices(size(dims))[i, j] for j in a2 for i in a1)
    for (i, li) in enumerate(lis)
        iyperm[i] = li
    end
    return iyperm
end
function _iyperm(dims::Tuple{<:Dimension,<:Dimension,<:Dimension})
    a1, a2, a3 = map(dims) do d
        l = parent(d)
        LA.ordered_firstindex(l):_order_step(l):LA.ordered_lastindex(l)
    end
    iyperm = Array{Int}(undef, length(a1) * length(a2) * length(a3))
    lis = (LinearIndices(size(dims))[i, j, k] for k in a3 for j in a2 for i in a1)
    for (i, li) in enumerate(lis)
        Iyperm[i] = li
    end
    return iyperm
end

_order_step(x) = _order_step(order(x))
_order_step(::ReverseOrdered) = -1
_order_step(::ForwardOrdered) = 1

# _fill_point!
# Fill a raster with `fill` where points are inside raster pixels
function _fill_points!(B::AbstractRaster, points; kw...)
    # Just find which pixels contian the points, and set them to true
    _without_mapped_crs(B) do B1
        _fill_points1!(B, points; kw...)
    end
    return B
end

function _fill_points1!(B::AbstractRaster, points::AbstractVector; kw...)
    for point in points
        _fill_points1!(B, point; kw...)
    end
end
function _fill_points1!(B::AbstractRaster, point::Union{Tuple,<:AbstractVector{<:Real}};
    order, fill=true, atol=nothing, kw...
)
    selectors = map(dims(B, order), ntuple(i -> i, length(order))) do d, i
        _at_or_contains(d, point[i], atol)
    end
    if hasselection(B, selectors)
        B[selectors...] = fill
    end
end

# _fill_linestring!
# Fill a raster with `fill` where pixels touch lines in a linestring
# Separated for a type stability function barrier
function _fill_linestring!(B::AbstractRaster, linestring::AbstractVector{<:AbstractVector{<:AbstractVector}}, args...; kw...)
    for ls in linestring
        _fill_linestring!(B, ls, args...; kw...)
    end
end
function _fill_linestring!(B::AbstractRaster, linestring::Poly; order, fill=true, kw...)
    linestring = collect(linestring)
    # Flip the order with a view to keep our alg simple
    forward_ordered_B = reduce(dims(B); init=B) do A, d
        if DD.order(d) isa ReverseOrdered
            A = view(A, rebuild(d, lastindex(d):-1:firstindex(d)))
            set(A, d => reverse(d))
        else
            A
        end
    end
    # For each line segment, burn the line into B with `fill` values
    _fill_line_segments!(forward_ordered_B, linestring, fill, order)
    return B
end

function _fill_line_segments!(B::AbstractArray, linestring::Poly, fill, order::Tuple{<:Dimension,<:Dimension})
    isfirst = true
    local firstpoint, lastpoint
    for point in linestring
        if isfirst
            isfirst = false
            firstpoint = point
            lastpoint = point
            continue
        end
        if point == firstpoint
            isfirst = true
        end
        line = (
            start=(x=lastpoint[1], y=lastpoint[2]),
            stop=(x=point[1], y=point[2])
        )
        _fill_line!(B, line, fill, order)
        lastpoint = point
    end
end

# fill_line!
# Fill a raster with `fill` where pixels touch a line
# TODO: generalise to 3d and Irregular spacing?
function _fill_line!(A::AbstractRaster, line, fill, order::Tuple{<:Dimension,<:Dimension})
    regular = map(dims(A, order)) do d
        lookup(d) isa AbstractSampled && span(d) isa Regular
    end
    msg = """
        Can only fill lines where dimensions are regular.
        Consider using `boundary=center` reprojecting the crs,
        or make an issue in Rasters.jl on github if you need this to work.
        """
    all(regular) || throw(ArgumentError(msg))

    xd, yd = order
    x_scale = abs(step(span(A, X)))
    y_scale = abs(step(span(A, Y)))
    raw_x_offset = bounds(A, xd)[1]
    raw_y_offset = bounds(A, yd)[1]
    raw_start, raw_stop = line.start, line.stop # Float
    start = (; x=(raw_start.x - raw_x_offset)/x_scale, y=(raw_start.y - raw_y_offset)/y_scale)
    stop = (; x=(raw_stop.x - raw_x_offset)/x_scale, y=(raw_stop.y - raw_y_offset)/y_scale)
    x, y = floor(Int, start.x) + 1, floor(Int, start.y) + 1 # Int

    diff_x = stop.x - start.x
    diff_y = stop.y - start.y
    step_x = signbit(diff_x) * -2 + 1
    step_y = signbit(diff_y) * -2 + 1
    # Ray/Slope calculations
    # Straight distance to the first vertical/horizontal grid boundaries
    xoffset = stop.x > start.x ? (ceil(start.x) - start.x) : (start.x - floor(start.x))
    yoffset = stop.y > start.y ? (ceil(start.y) - start.y) : (start.y - floor(start.y))
    # Angle of ray/slope.
    angle = atan(-diff_y, diff_x)
    # max: How far to move along the ray to cross the first cell boundary.
    # delta: How far to move along the ray to move 1 grid cell.
    cs = cos(angle)
    si = sin(angle)
    max_x, delta_x = if isapprox(cs, zero(cs); atol=1e-10)
        -Inf, Inf
    else
        1.0 / cs, xoffset / cs
    end
    max_y, delta_y = if isapprox(si, zero(si); atol=1e-10)
        -Inf, Inf
    else
        1.0 / si, yoffset / si
    end
    # Travel one grid cell at a time.
    manhattan_distance = floor(Int, abs(floor(start.x) - floor(stop.x)) + abs(floor(start.y) - floor(stop.y)))
    # For arbitrary dimension indexing
    dimconstructors = map(DD.basetypeof, dims(A, order))
    for t in 0:manhattan_distance
        D = map((d, o) -> d(o), dimconstructors, (x, y))
        if checkbounds(Bool, A, D...)
            if fill isa Function
                @inbounds A[D...] = fill(A[D...])
            else
                @inbounds A[D...] = fill
            end
        end
        # Only move in either X or Y coordinates, not both.
        if abs(max_x) <= abs(max_y)
            max_x += delta_x
            x += step_x
        else
            max_y += delta_y
            y += step_y
        end
    end
    return A
end
function _fill_line!(A::AbstractRaster, line, fill, order::Tuple{Vararg{<:Dimension}})
    msg = """"
        Converting a `:line` geometry to raster is currently only implemented for 2d lines.
        Make a Rasters.jl github issue if you need this for more dimensions.
        """
    throw(ArgumentError(msg))
end

# PolygonInbounds.jl setup

# _flat_nodes
# Convert a geometry/nested vectors to a flat iterator of point nodes for PolygonInbounds
_flat_nodes(A::GI.AbstractGeometry) = _flat_nodes(GI.coordinates(A))
_flat_nodes(A::GI.AbstractFeature) = _flat_nodes(GI.geometry(A))
_flat_nodes(A::AbstractVector{<:GI.AbstractGeometry}) = Iterators.flatten(map(_flat_nodes, A))
_flat_nodes(A::AbstractVector{<:GI.AbstractFeature}) = Iterators.flatten(map(_flat_nodes, A))
function _flat_nodes(A::AbstractVector{<:Union{Missing,<:GI.AbstractGeometry}})
    Iterators.flatten(map(_flat_nodes, A))
end
function _flat_nodes(A::AbstractVector{<:AbstractVector{<:AbstractVector}})
    Iterators.flatten(map(_flat_nodes, A))
end
_flat_nodes(A::AbstractVector{<:AbstractVector{<:Real}}) = A
_flat_nodes(A::AbstractVector{<:Real}) = A
_flat_nodes(A::AbstractVector{<:Tuple}) = A
_flat_nodes(iter::Base.Iterators.Flatten) = iter

_node_vecs(A::GI.AbstractGeometry) = GI.coordinates(A)
_node_vecs(A::GI.AbstractFeature) = _node_vecs(GI.geometry(A))
function _node_vecs(A::AbstractVector{<:Union{Missing,<:GI.AbstractGeometry,<:GI.AbstractFeature}})
    map(_node_vecs, A)
end
function _node_vecs(A::AbstractVector{<:AbstractVector{<:AbstractVector}})
    A
end
_node_vecs(A::AbstractVector{<:AbstractVector{<:AbstractFloat}}) = A
_node_vecs(A::AbstractVector{<:Tuple}) = A
_node_vecs(iter::Base.Iterators.Flatten) = collect(iter)

# _to_edges
# Convert a polygon to the `edges` needed by PolygonInbounds
function to_edges_and_nodes(poly)
    p1 = first(_flat_nodes(poly))
    edges = Vector{Tuple{Int,Int}}(undef, 0)
    nodes = Vector{typeof(p1)}(undef, 0)
    nodenum = 0
    _to_edges!(edges, nodes, nodenum, poly)
    edges = PermutedDimsArray(reinterpret(reshape, Int, edges), (2, 1))
    lastx = 0
    for (i, x) in enumerate(view(edges, :, 1))
        x == lastx + 1 || error("edge dup at $x $lastx $i")
        lastx = x
    end
    return edges, nodes
end

# _to_edges!(edges, edgenum, poly)
# fill edges vector with edges from polygon, numbered starting at `edgenum`
function _to_edges!(edges, nodes, pointnum, poly::AbstractVector{<:Union{Missing,<:GI.AbstractGeometry}})
    foldl(poly; init=pointnum) do n, p
        _to_edges!(edges, nodes, n, p)
    end
end
function _to_edges!(edges, nodes, pointnum, poly::GI.AbstractGeometry)
    _to_edges!(edges, nodes, pointnum, GI.coordinates(poly))
end
function _to_edges!(edges, nodes, pointnum, poly::AbstractVector{<:AbstractVector})
    foldl(poly; init=pointnum) do n, p
        _to_edges!(edges, nodes, n, p)
    end
end
# Analyse a single polygon
function _to_edges!(
    edges, nodes, pointnum,
    poly::AbstractVector{<:Union{<:NTuple{<:Any,T},<:AbstractVector{T}}}
) where T<:Real
    # Any polygon may actually be made up of multiple sub-polygons,
    # indicated by returning to the first point in the sub-polygon.
    # So we keep track of the first point, and insert closing edges
    # wherever we find it repeated.
    fresh_start = true
    startpoint = first(poly)
    start_pointnum = pointnum
    added_nodes = 0
    for (n, point) in enumerate(poly)
        if fresh_start
            # The first edge in the sub-polygon
            start_pointnum = pointnum + n
            startedge = (start_pointnum, start_pointnum + 1)
            push!(edges, startedge)
            startpoint = point
            fresh_start = false
        elseif point == startpoint || n == length(poly)
            # The closing edge of a sub-polygon
            closingedge = (pointnum + n, start_pointnum)
            push!(edges, closingedge)
            fresh_start = true
        else
            # A regular edge somewhere in a sub-polygon
            edge = (pointnum + n, pointnum + n + 1)
            push!(edges, edge)
        end
        # Track the total number of nodes we have added
        added_nodes = n
    end
    append!(nodes, poly)
    nextpoint = pointnum + added_nodes
    return nextpoint
end

# _geom_bounds
# Get the bounds of a geometry
function _geom_bounds(geom, order)
    nodes = _flat_nodes(geom)
    bounds = map(ntuple(identity, length(order))) do i
        extrema(p[i] for p in nodes)
    end
    return bounds
end

function _wrapped_geom_bounds(x, geom, order)
    wrapped_bounds = map(rebuild, dims(x, order), _geom_bounds(geom, order))
    return dims(wrapped_bounds, dims(x))
end

# bounds_overlap
# Check if two bounds tuples overlap
function bounds_overlap(a, b)
    axes_cross = map(a, b) do (a_min, a_max), (b_min, b_max)
        if a_max >= b_max
            a_min <= b_max
        else
            a_max >= b_min
        end
    end
    return all(axes_cross)
end

# bbox_might_overlap
# Check if there is a bbox for the geometry
function bbox_might_overlap(x, geom, order)
    x_bnds = bounds(x, order)
    p_bbox = GI.bbox(geom)
    # If there is no bbox available just act as if it
    # overlaps and let the checks later on sort it out
    p_bbox isa Nothing && return true
    bbox_dims = length(p_bbox) ÷ 2
    p_bnds = [(p_bbox[i], p_bbox[i+bbox_dims]) for i in 1:bbox_dims]

    return bounds_overlap(p_bnds, x_bnds)
end

# _geom_shape
# Get the shape category for a geometry
function _geom_shape(geom)
    typ = GI.geotype(geom)
    if typ in (:Point, :MultiPoint)
        return :point
    elseif typ in (:LineString, :MultiLineString)
        return :line
    elseif typ in (:Polygon, :MultiPolygon)
        return :polygon
    else
        throw(ArgumentError("Geometry type $typ not known"))
    end
    # TODO: What to do with :GeometryCollection
end

# Copied from PolygonInbounds, to add extra keyword arguments
# PR to include these when this has solidied
function inpoly2(vert, node, edge=zeros(Int);
    atol::T=0.0, rtol::T=NaN, iyperm=nothing,
    vmin=nothing, vmax=nothing, pmin=nothing, pmax=nothing
) where T<:AbstractFloat
    rtol = !isnan(rtol) ? rtol : iszero(atol) ? eps(T)^0.85 : zero(T)
    poly = PolygonInbounds.PolygonMesh(node, edge)
    points = PolygonInbounds.PointsInbound(vert)
    npoints = length(points)

    vmin = isnothing(vmin) ? minimum(points) : vmin
    vmax = isnothing(vmax) ? maximum(points) : vmax
    pmin = isnothing(pmin) ? minimum(poly) : pmin
    pmax = isnothing(pmax) ? maximum(poly) : pmax

    lbar = sum(pmax - pmin)
    tol = max(abs(rtol * lbar), abs(atol))

    ac = PolygonInbounds.areacount(poly)
    stat = ac > 1 ? falses(npoints, 2, ac) : falses(npoints, 2)
    # flip coordinates so expected effort is minimal

    dvert = vmax .- vmin
    if isnothing(iyperm)
        ix = dvert[1] < dvert[2] ? 1 : 2
        iyperm = sortperm(points, 3 - ix)
    else
        ix = 1
    end

    PolygonInbounds.inpoly2!(points, iyperm, poly, ix, tol, stat)
    return stat
end

_shape_error(shape) = throw(ArgumentError("`shape` must be :point, :line or :polygon, got $shape"))
_boundary_error(boundary) = throw(ArgumentError("`boundary` can be :touches, :inside, or :center, got $boundary"))


text:
module Adam
    using LoopVectorization

    function fit(;model::Any, input_data::Array{Float32}, output_data::Array{Float32}, loss_function::Any, monitor::Any, α::Float64=0.001, epochs::Int64=20, batch::Real=32, β₁::Float64=0.9, β₂::Float64=0.999, ϵ::Float64=1e-8)
        model.initialize(model, 1)
        input_shape = size(input_data)[1:end-1]
        output_shape = size(output_data)[1:end-1]
        batch_size = ceil(Int64, size(input_data, length(input_shape)+1)/batch)

        current_input_data = zeros(Float32, input_shape..., 1)
        current_output_data = zeros(Float32, output_shape..., 1)

        for e in 1:epochs
            print("Epoch ", e, "\n[")
            @time begin
                for t in 1:batch_size
                    index = rand(1:size(input_data)[end])
                    current_input_data .= reshape(selectdim(input_data, length(input_shape)+1, index), input_shape..., 1)
                    current_output_data .= reshape(selectdim(output_data, length(output_shape)+1, index), output_shape..., 1)

                    if t%ceil(batch_size/50)==0
                        print("=")
                    end
                    model.update(model, current_input_data, current_output_data, loss_function, monitor, "Adam", α, t, β₁, β₂, ϵ)
                end
                print("] with loss ", model.loss, "\nTime usage: ")
                model.loss = 0.0
            end
        end
    end
end


text:
@doc raw"""
    dominoes(imgname, colr, constrnt)

Approximate a grayscale image with complete sets of double-nine dominoes.

# Resources
Watch Bob Bosch of Oberlin College talk at Google about domino mosaics, and his interpretation of this art form as an integer program.
https://www.youtube.com/watch?v=g3CiVrN-BnY

IP modeling details can be found in this article.
http://www.optimization-online.org/DB_FILE/2003/09/722.pdf

# Examples
```julia-repl
julia> include("ModuleDominoes.jl")
julia> using .Dominoes
julia> createMosaic("Frankenstein", "white", "constrained")
```

# Author
Kalyani Nagaraj
March 2021
"""
module Dominoes

using JuMP
using GLPK
using DelimitedFiles
using MAT
using ImageView, Images, FileIO

include("dominoes.jl")
include("createMosaic.jl")

export createMosaic

end

max_stars_repo_name: kalyaninagaraj/Domino-Mosaics

max_issues_repo_name: kalyaninagaraj/Domino-Mosaics

max_forks_repo_name: kalyaninagaraj/Domino-Mosaics


text:
abstract type Shape end

struct Ellipse{T} <: Shape
    x::T
    y::T
    A::T
    B::T
    α::T
    ρ::T
end

function isinside(e::Ellipse, x, y)
    # https://stackoverflow.com/questions/7946187/point-and-ellipse-rotated-position-test-algorithm
    (cosd(e.α) * (x - e.x) + sind(e.α) * (y - e.y))^2 / e.A^2 +
    (sind(e.α) * (x - e.x) - cosd(e.α) * (y - e.y))^2 / e.B^2 ≤ one(x)
end

function projection(e::Ellipse, θᵢ, tᵢ)
    # p. 56
    s = √(e.x^2 + e.y^2)
    γ = atan(e.y, e.x)
    θ = θᵢ - deg2rad(e.α)
    t = tᵢ - s * cos(γ - θᵢ)
    a²θ = e.A^2 * cos(θ)^2 + e.B^2 * sin(θ)^2
    abs(t) ≤ √(a²θ) ? 2 * e.ρ * e.A * e.B * √(a²θ - t^2) / a²θ : zero(a²θ)
end

function projection(
    ellipses::AbstractVector{<:Ellipse},
    θ::AbstractVector,
    t::AbstractVector,
)
    [sum(projection(e, θᵢ, tᵢ) for e in ellipses) for tᵢ in t, θᵢ in θ]
end

# Table 3.1 summary of parameters for tomography simulations
SheppLoganPhantom = [
    Ellipse(0.0, 0.0, 0.92, 0.69, 90.0, 2.0),
    Ellipse(0.0, -0.0184, 0.874, 0.6624, 90.0, -0.98),
    Ellipse(0.22, 0.0, 0.31, 0.11, 72.0, -0.02),
    Ellipse(-0.22, 0.0, 0.41, 0.16, 108.0, -0.02),
    Ellipse(0.0, 0.35, 0.25, 0.21, 90.0, 0.01),
    Ellipse(0.0, 0.1, 0.046, 0.046, 0.0, 0.01),
    Ellipse(0.0, -0.1, 0.046, 0.046, 0.0, 0.01),
    Ellipse(-0.08, -0.605, 0.046, 0.023, 0.0, 0.01),
    Ellipse(0.0, -0.605, 0.023, 0.023, 0.0, 0.01),
    Ellipse(0.06, -0.605, 0.046, 0.023, 90.0, 0.01),
]

function construct(ellipses::AbstractVector{<:Ellipse}, N::Integer)
    f = zeros(N, N)
    for (i, xᵢ) in enumerate(range(-1, 1, length = N)),
        (j, yⱼ) in enumerate(range(-1, 1, length = N))

        @inbounds f[j, i] = sum([isinside(e, xᵢ, yⱼ) ? e.ρ : zero(e.ρ) for e in ellipses])
    end
    f
end


text:
##                                                                            ##
##              Problem set #5 for Structural modelling course                ##
##                               Nicolás Urdaneta                             ##
##                                                                            ##

Pkg.add("DataFramesMeta")

using Random
using LinearAlgebra
using Statistics
using Optim
using DataFrames
using DataFramesMeta
using CSV
using HTTP
using GLM

# read in function to create state transitions for dynamic model
include("create_grids.jl")

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 1: reshaping the data
#:::::::::::::::::::::::::::::::::::::::::::::::::::
# load in the data

url = "https://raw.githubusercontent.com/OU-PhD-Econometrics/fall-2020/master/ProblemSets/PS5-ddc/busdataBeta0.csv"
df=CSV.read(HTTP.get(url).body)

# create bus id variable
df = @transform(df, bus_id = 1:size(df,1))

#---------------------------------------------------
# reshape from wide to long (must do this twice be-
# cause DataFrames.stack() requires doing it one
# variable at a time)
#---------------------------------------------------
# first reshape the decision variable
dfy = @select(df, :bus_id,:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20,:RouteUsage,:Branded)
dfy_long = DataFrames.stack(dfy, Not([:bus_id,:RouteUsage,:Branded]))
rename!(dfy_long, :value => :Y)
dfy_long = @transform(dfy_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfy_long, Not(:variable))

# next reshape the odometer variable
dfx = @select(df, :bus_id,:Odo1,:Odo2,:Odo3,:Odo4,:Odo5,:Odo6,:Odo7,:Odo8,:Odo9,:Odo10,:Odo11,:Odo12,:Odo13,:Odo14,:Odo15,:Odo16,:Odo17,:Odo18,:Odo19,:Odo20)
dfx_long = DataFrames.stack(dfx, Not([:bus_id]))
rename!(dfx_long, :value => :Odometer)
dfx_long = @transform(dfx_long, time = kron(collect([1:20]...),ones(size(df,1))))
select!(dfx_long, Not(:variable))

# join reshaped df's back together
df_long = leftjoin(dfy_long, dfx_long, on = [:bus_id,:time])
sort!(df_long,[:bus_id,:time])

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 2: estimate a static version of the model
#:::::::::::::::::::::::::::::::::::::::::::::::::::

# [Your turn]

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 3a: read in data for dynamic model
#:::::::::::::::::::::::::::::::::::::::::::::::::::

# [ Load in the data ]

Y = Matrix(df[:,[:Y1,:Y2,:Y3,:Y4,:Y5,:Y6,:Y7,:Y8,:Y9,:Y10,:Y11,:Y12,:Y13,:Y14,:Y15,:Y16,:Y17,:Y18,:Y19,:Y20]])

# [ convert other data frame columns to matrices]

#:::::::::::::::::::::::::::::::::::::::::::::::::::
# Question 3b: generate state transition matrices
#:::::::::::::::::::::::::::::::::::::::::::::::::::
zval,zbin,xval,xbin,xtran = create_grids()


text:
# Outer approximation algorithm

mutable struct Optimizer <: MOI.AbstractOptimizer
    verbose::Bool
    tol_feas::Float64
    oa_solver::Union{Nothing, MOI.OptimizerWithAttributes}
    use_iterative_method::Bool

    oa_opt::Union{Nothing, MOI.ModelLike}
    con_vars::Vector{Vector{VI}}
    con_sets::Vector{MOI.AbstractVectorSet}
    num_cuts::Int

    # for iterative method
    num_iters::Int
    status::Symbol

    # for OA-solver-driven method
    has_discrete_constraint::Bool
    dummy_discrete_constraint::Union{Nothing, Tuple{VI, CI}}
    num_callbacks::Int

    function Optimizer()
        opt = new()
        opt.verbose = true
        opt.tol_feas = 1e-7
        opt.oa_solver = nothing
        opt.use_iterative_method = false
        return _empty(opt)
    end
end

function _empty(opt::Optimizer)
    opt.oa_opt = nothing
    opt.con_vars = Vector{VI}[]
    opt.con_sets = MOI.AbstractVectorSet[]
    opt.num_cuts = 0
    opt.num_iters = 0
    opt.status = :Empty
    opt.has_discrete_constraint = false
    opt.dummy_discrete_constraint = nothing
    opt.num_callbacks = 0
    return opt
end

function MOI.optimize!(opt::Optimizer)
    # add initial cuts
    for k in eachindex(opt.con_sets)
        add_init_cuts(k, opt)
    end
    if opt.verbose
        println("added $(opt.num_cuts) initial cuts; starting OA solver")
    end

    if opt.use_iterative_method
        iterative_method(opt)
    else
        oa_solver_driven_method(opt)
    end

    if opt.verbose
        oa_status = MOI.get(opt.oa_opt, MOI.TerminationStatus())
        println("OA solver finished with status $oa_status, after $(opt.num_cuts) cuts")
        if opt.use_iterative_method
            println("iterative method used $(opt.num_iters) iterations\n")
        else
            println("OA solver driven method used $(opt.num_callbacks) callbacks\n")
        end
    end
    return
end

# TODO handle statuses properly
function iterative_method(opt::Optimizer)
    while true
        MOI.optimize!(opt.oa_opt)

        oa_opt_status = MOI.get(opt.oa_opt, MOI.TerminationStatus())
        if oa_opt_status == MOI.OPTIMAL
            obj_bound = MOI.get(opt.oa_opt, MOI.ObjectiveBound())
            if opt.verbose
                @printf("%5d %8d %12.4e\n", opt.num_iters, opt.num_cuts, obj_bound)
            end
        elseif oa_opt_status == MOI.DUAL_INFEASIBLE
            if opt.verbose
                println("OA solver status is $oa_opt_status; relaxation likely unbounded")
            end
            opt.status = :DualInfeasibleRelaxation
            break
        elseif oa_opt_status == MOI.INFEASIBLE
            if opt.verbose
                println("infeasibility detected; terminating")
            end
            opt.status = :PrimalInfeasible
            break
        elseif oa_opt_status == MOI.INFEASIBLE_OR_UNBOUNDED
            if opt.verbose
                println("OA solver status is $oa_opt_status; assuming infeasibility")
            end
            opt.status = :PrimalInfeasible
            break
        else
            error("OA solver status $oa_opt_status not handled")
        end

        # TODO
        # if opt.num_iters == opt.max_iters
        #     opt.verbose && println("iteration limit reached; terminating")
        #     opt.status = :IterationLimit
        #     break
        # end

        # if time() - start_time >= opt.time_limit
        #     opt.verbose && println("time limit reached; terminating")
        #     opt.status = :TimeLimit
        #     break
        # end

        is_cut_off = any(add_sep_cuts(k, opt) for k in eachindex(opt.con_sets))
        if !is_cut_off
            if opt.verbose
                println("no cuts were added; terminating")
            end
            opt.status = :Optimal
            break
        end

        opt.num_iters += 1
        flush(stdout)
    end
    return
end

function oa_solver_driven_method(opt::Optimizer)
    if !opt.has_discrete_constraint
        # no discrete constraints, so add one to force OA solver to use lazy callbacks
        if opt.verbose
            println("model has no discrete constraints; adding a dummy discrete constraint")
        end
        opt.dummy_discrete_constraint = MOI.add_constrained_variable(opt, MOI.Integer())
    end

    function lazy_callback(cb)
        opt.num_callbacks += 1
        is_cut_off = any(add_sep_cuts(k, opt, cb) for k in eachindex(opt.con_sets))
        if !is_cut_off && opt.verbose
            println("no cuts were added during callback")
        end
    end
    MOI.set(opt.oa_opt, MOI.LazyConstraintCallback(), lazy_callback)

    MOI.optimize!(opt.oa_opt)

    if !isnothing(opt.dummy_discrete_constraint)
        # delete dummy discrete constraint
        (vi, ci) = opt.dummy_discrete_constraint
        MOI.delete(opt.oa_opt, ci)
        MOI.delete(opt.oa_opt, vi)
        opt.dummy_discrete_constraint = nothing
    end

    return
end

function add_init_cuts(k::Int, opt::Optimizer)
    con_vars = opt.con_vars[k]
    con_set = opt.con_sets[k]

    cuts = Cuts.get_init_cuts(con_set)
    for cut in cuts
        cut_expr = SAF(MOI.ScalarAffineTerm.(cut, con_vars), 0.0)
        MOI.add_constraint(opt.oa_opt, cut_expr, MOI.GreaterThan(0.0))
    end

    opt.num_cuts += length(cuts)
    return nothing
end

function add_sep_cuts(k::Int, opt::Optimizer, cb = nothing)
    vars = opt.con_vars[k]
    point = _get_var_primal(vars, opt.oa_opt, cb)
    @assert !any(isnan, point)

    cuts = Cuts.get_sep_cuts(point, opt.con_sets[k], opt.tol_feas)
    for cut in cuts
        cut_expr = SAF(MOI.ScalarAffineTerm.(cut, vars), 0.0)
        _add_constraint(cut_expr, opt.oa_opt, cb)
    end

    opt.num_cuts += length(cuts)
    return !isempty(cuts)
end

function _add_constraint(cut_expr::SAF, oa_opt::MOI.ModelLike, ::Nothing)
    return MOI.add_constraint(oa_opt, cut_expr, MOI.GreaterThan(0.0))
end

function _add_constraint(cut_expr::SAF, oa_opt::MOI.ModelLike, cb)
    return MOI.submit(oa_opt, MOI.LazyConstraint(cb), cut_expr, MOI.GreaterThan(0.0))
end

function _get_var_primal(vars::Vector{VI}, oa_opt::MOI.ModelLike, ::Nothing)
    return MOI.get(oa_opt, MOI.VariablePrimal(), vars)
end

function _get_var_primal(vars::Vector{VI}, oa_opt::MOI.ModelLike, cb)
    return MOI.get(oa_opt, MOI.CallbackVariablePrimal(cb), vars)
end


text:
module Streams

export Stream, closebody, isaborted,
       header, hasheader,
       setstatus, setheader

import ..HTTP
using ..IOExtras
using ..Parsers
using ..Messages
import ..Messages: header, hasheader, setheader,
                   writeheaders, writestartline
import ..ConnectionPool: getrawstream, byteview
import ..@require, ..precondition_error
import ..@ensure, ..postcondition_error
import ..@debug, ..DEBUG_LEVEL

mutable struct Stream{M <: Message, S <: IO} <: IO
    message::M
    parser::Parser
    stream::S
    writechunked::Bool
    readchunked::Bool
    ntoread::Int
end

"""
    Stream(::IO, ::Request, ::Parser)

Creates a `HTTP.Stream` that wraps an existing `IO` stream.

 - `startwrite(::Stream)` sends the `Request` headers to the `IO` stream.
 - `write(::Stream, body)` sends the `body` (or a chunk of the body).
 - `closewrite(::Stream)` sends the final `0` chunk (if needed) and calls
   `closewrite` on the `IO` stream. When the `IO` stream is a
   [`HTTP.ConnectionPool.Transaction`](@ref), calling `closewrite` releases
   the [`HTTP.ConnectionPool.Connection`](@ref) back into the pool for use by the
   next pipelined request.

 - `startread(::Stream)` calls `startread` on the `IO` stream then
    reads and parses the `Response` headers.  When the `IO` stream is a
   [`HTTP.ConnectionPool.Transaction`](@ref), calling `startread` waits for other
   pipelined responses to be read from the [`HTTP.ConnectionPool.Connection`](@ref).
 - `eof(::Stream)` and `readavailable(::Stream)` parse the body from the `IO`
    stream.
 - `closeread(::Stream)` reads the trailers and calls `closeread` on the `IO`
    stream.  When the `IO` stream is a [`HTTP.ConnectionPool.Transaction`](@ref),
    calling `closeread` releases the readlock and allows the next pipelined
    response to be read by another `Stream` that is waiting in `startread`.
    If the `Parser` has not recieved a complete response, `closeread` throws
    an `EOFError`.
"""

Stream(r::M, p::Parser, io::S) where {M, S} =
    Stream{M,S}(r, p, io, false, false, 0)

header(http::Stream, a...) = header(http.message, a...)
setstatus(http::Stream, status) = (http.message.response.status = status)
setheader(http::Stream, a...) = setheader(http.message.response, a...)
getrawstream(http::Stream) = getrawstream(http.stream)

IOExtras.isopen(http::Stream) = isopen(http.stream)

# Writing HTTP Messages

messagetowrite(http::Stream{Response}) = http.message.request
messagetowrite(http::Stream{Request}) = http.message.response

IOExtras.iswritable(http::Stream) = iswritable(http.stream)

function IOExtras.startwrite(http::Stream)
    if !iswritable(http.stream)
        startwrite(http.stream)
    end
    m = messagetowrite(http)
    if !hasheader(m, "Content-Length") &&
       !hasheader(m, "Transfer-Encoding") &&
       !hasheader(m, "Upgrade")
        http.writechunked = true
        setheader(m, "Transfer-Encoding" => "chunked")
    else
        http.writechunked = ischunked(m)
    end
    writeheaders(http.stream, m)
end

function Base.unsafe_write(http::Stream, p::Ptr{UInt8}, n::UInt)
    if n == 0
        return 0
    end
    if !iswritable(http) && isopen(http.stream)
        startwrite(http)
    end
    if !http.writechunked
        return unsafe_write(http.stream, p, n)
    end
    return write(http.stream, hex(n), "\r\n") +
           unsafe_write(http.stream, p, n) +
           write(http.stream, "\r\n")
end

"""
    closebody(::Stream)

Write the final `0` chunk if needed.
"""

function closebody(http::Stream)
    if http.writechunked
        http.writechunked = false
        write(http.stream, "0\r\n\r\n")
    end
end

function IOExtras.closewrite(http::Stream{Response})
    if !iswritable(http)
        return
    end
    closebody(http)
    closewrite(http.stream)
end

function IOExtras.closewrite(http::Stream{Request})
    @require iswritable(http)

    closebody(http)
    closewrite(http.stream)

    if hasheader(http.message, "Connection", "close") ||
       http.message.version < v"1.1" &&
      !hasheader(http.message, "Connection", "keep-alive")

        @debug 1 "✋  \"Connection: close\": $(http.stream)"
        close(http.stream)
    end
end

# Reading HTTP Messages

IOExtras.isreadable(http::Stream) = isreadable(http.stream)

function IOExtras.startread(http::Stream)

    startread(http.stream)

    reset!(http.parser)
    readheaders(http.stream, http.parser, http.message)
    handle_continue(http)

    http.readchunked = ischunked(http.message)
    http.ntoread = bodylength(http.message)

    return http.message
end

"""
100 Continue
https://tools.ietf.org/html/rfc7230#section-5.6
https://tools.ietf.org/html/rfc7231#section-6.2.1
"""

function handle_continue(http::Stream{Response})
    if http.message.status == 100
        @debug 1 "✅  Continue:   $(http.stream)"
        reset!(http.parser)
        readheaders(http.stream, http.parser, http.message)
    end

end

function handle_continue(http::Stream{Request})
    if hasheader(http.message, "Expect", "100-continue")
        if !iswritable(http.stream)
            startwrite(http.stream)
        end
        @debug 1 "✅  Continue:   $(http.stream)"
        writeheaders(http.stream, Response(100))
    end
end

function Base.eof(http::Stream)
    if !headerscomplete(http.message)
        startread(http)
    end
    if http.ntoread == 0
        return true
    end
    if eof(http.stream)
        return true
    end
    return false
end

function Base.readavailable(http::Stream)::ByteView
    @require headerscomplete(http.message)

    if http.ntoread == 0
        return nobytes
    end
    if nb_available(http.stream) > http.ntoread
        bytes = byteview(read(http.stream, http.ntoread))
    else
        bytes = readavailable(http.stream)
    end
    l = length(bytes)
    if l == 0
        return nobytes
    end
    if http.readchunked
        bytes, excess = parsebody(http.parser, bytes)
        unread!(http, excess)
        if bodycomplete(http.parser)
            http.ntoread = 0
        end
    end
    if http.ntoread != unknown_length
        http.ntoread -= length(bytes)
    end
    @ensure http.ntoread >= 0
    return bytes
end

IOExtras.unread!(http::Stream, excess) = unread!(http.stream, excess)

function Base.read(http::Stream)
    buf = IOBuffer()
    write(buf, http)
    return take!(buf)
end

"""
    isaborted(::Stream{Response})

Has the server signaled that it does not wish to receive the message body?

"If [the response] indicates the server does not wish to receive the
 message body and is closing the connection, the client SHOULD
 immediately cease transmitting the body and close the connection."
[RFC7230, 6.5](https://tools.ietf.org/html/rfc7230#section-6.5)
"""

function isaborted(http::Stream{Response})

    if iswritable(http.stream) &&
       iserror(http.message) &&
       hasheader(http.message, "Connection", "close")
        @debug 1 "✋  Abort on $(sprint(writestartline, http.message)): " *
                 "$(http.stream)"
        @debug 2 "✋  $(http.message)"
        return true
    end
    return false
end

function IOExtras.closeread(http::Stream{Response})

    # Discard body bytes that were not read...
    while !eof(http)
        readavailable(http)
    end

    # Read trailers...
    if bodycomplete(http.parser) && !messagecomplete(http.parser)
        readtrailers(http.stream, http.parser, http.message)
    end

    if http.ntoread != unknown_length && http.ntoread > 0
        # Error if Message is not complete...
        close(http.stream)
        throw(EOFError())
    elseif hasheader(http.message, "Connection", "close")
        # Close conncetion if server sent "Connection: close"...
        @debug 1 "✋  \"Connection: close\": $(http.stream)"
        close(http.stream)
    elseif isreadable(http.stream)
        closeread(http.stream)
    end

    return http.message
end

function IOExtras.closeread(http::Stream{Request})
    if http.ntoread != unknown_length && http.ntoread > 0
        # Error if Message is not complete...
        close(http.stream)
        throw(EOFError())
    end
    if isreadable(http)
        closeread(http.stream)
    end
end

end #module Streams


text:
using Primes, Printf
function strfactor(n::Integer)
    n > -2 || return "-1 × " * strfactor(-n)
    isprime(n) || n < 2 && return dec(n)
    f = factor(Vector{typeof(n)}, n)
    return join(f, " × ")
end

lo, hi = -4, 40
println("Factor print $lo to $hi:")
for n in lo:hi
    @printf("%5d = %s\n", n, strfactor(n))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module BSR

using MATLAB

export
    gpb,
    gpb_owt_ucm,
    gpb_owt_ucm_strength

include("matutil.jl")
include("bsrfun.jl")

end


text:
function strip_maybe(T)
    T === Missing && return Union{}
    !isa(T, Union) && return T
    return Union{strip_maybe(T.a), strip_maybe(T.b)}
end

optional_eltype(::Type{Union{T, Void}}) where {T} = T
optional_eltype(_) = Union{}

macro from_json(T, renamings)
    syms = Symbol[]
    body = Expr[]
    for (N,fT) in zip(fieldnames(T), T.types)
        sym = gensym()
        push!(syms, sym)
        name = string(N)
        if haskey(renamings, N)
            val = renamings[N]
            if isa(val, String)
                name = val
            else
                @assert isa(val, Expr)
                push!(body, :($sym = $val))
                continue
            end
        end
        sfT = strip_maybe(fT)
        oeT = optional_eltype(sfT)
        if oeT !== Union{}
            sfT = oeT
        end
        if sfT <: QboRef
            name *= "Ref"
            op = :($sym = from_json($sfT, data[$(name)]))
        elseif sfT <: Vector
            op = :($sym = map(x->from_json($(eltype(sfT)), x),data[$(name)]))
        else
            op = :($sym = from_json($sfT, data[$(name)]))
        end
        if oeT !== Union{}
            push!(body, quote
                if haskey(data, $(name))
                    $op
                else
                    $sym = nothing
                end
            end)
        else
            push!(body, op)
        end
    end
    esc(quote
        $(body...)
        $T($(syms...))
    end)
end

macro to_json(T, renamings)
    body = Expr[:(ret = Dict{String, Any}())]
    for (N,fT) in zip(fieldnames(T), T.types)
        name = string(N)
        if haskey(renamings, N)
            val = renamings[N]
            if isa(val, String)
                name = val
            else
                @assert isa(val, Expr)
                push!(body, val)
                continue
            end
        end
        sfT = strip_maybe(fT)
        oeT = optional_eltype(sfT)
        if oeT !== Union{}
            sfT = oeT
        end
        if sfT <: QboRef
            name *= "Ref"
            op = :(ret[$name] = to_json($sfT, getfield(data, $(quot(N)))))
        elseif sfT <: Vector
            op = :(ret[$name] = map(to_json,getfield(data, $(quot(N)))))
        else
            op = :(ret[$name] = to_json(getfield(data, $(quot(N)))))
        end
        push!(body, quote
            if !ismissing(getfield(data, $(quot(N)))) && getfield(data, $(quot(N))) !== nothing
                $op
            end
        end)
    end
    a = esc(quote
        $(body...)
        ret
    end)
end

from_json(::Type{Int}, s::String) = parse(Int, s)
from_json(::Type{String}, s::String) = s
from_json(::Type{Date}, s::String) = Date(s)
from_json(::Type{Decimal}, d::Decimal) = d
from_json(::Type{Bool}, b::Bool) = b

to_json(s::String) = s
to_json(d::Date) = repr(d)
to_json(d::Decimal) = d


text:
mutable struct MultiGPCovars{KEY}
    D::Array{Float64,2}
    y::Vector{Float64}
    groupKeys::Vector{KEY}
    mgp::Dict{KEY,GPE}
    groupIndices::Dict{KEY,UnitRange{Int64}}
    p::Int
    dim::Int
    nobs::Int
    logNoise::Float64
    mean::Mean
    kernel::Kernel
    βkern::Kernel
    βdata::KernelData
    # Auxiliary data
    cK::PDMats.PDMat        # (k + obsNoise)
    alpha::Vector{Float64}  # (k + obsNoise)⁻¹y
    mll::Float64            # Marginal log-likelihood
    dmll::Vector{Float64}   # Gradient marginal log-likelihood
    function MultiGPCovars{KEY}(
            D::Array{Float64,2},
            y::Vector{Float64},
            groupKeys::Vector{KEY},
            mgp::Dict{KEY,GPE},
            groupIndices::Dict{KEY,UnitRange{Int64}},
            p::Int,
            dim::Int,
            nobs::Int,
            logNoise::Float64,
            mean::Mean,
            kernel::Kernel,
            βkern::Kernel
            ) where {KEY}
        size(D, 2) == nobs || throw("incompatible dimensions of covariates matrix and gaussian processes")
        size(D, 1) == p || throw("first dimension of D is not p")
        length(groupKeys) == length(groupIndices) || throw("groupKeys and groupIndices should have same length")
        length(groupKeys) == length(mgp) || throw("groupKeys and mgp should have same length")
        βdata=KernelData(βkern, D, D)
        mgpcv = new(D, y, groupKeys, mgp, groupIndices, p, dim, nobs, logNoise, mean, kernel, βkern, βdata)
        propagate_params!(mgpcv)
        initialise_mll!(mgpcv)
        return mgpcv
    end
end
function MultiGPCovars(Dlist::Vector{M}, gpList::Vector{GPE}, groupKeys::Vector{KEY}, βkern::Kernel) where {M<:AbstractMatrix{Float64}, KEY}
    total_nobs = sum([gp.nobs for gp in gpList])
    D = hcat(Dlist...)
    size(D, 2) == total_nobs || throw(ArgumentError("incompatible dimensions of covariates matrix and gaussian processes"))
    p = size(D,1)

    # get kernel and parameters from first GP in list
    first_gp = gpList[1]
    dim = first_gp.dim
    logNoise = first_gp.logNoise
    kern = first_gp.kernel
    mean = first_gp.mean
    ngroups = length(groupKeys)
    @assert ngroups == length(gpList)
    gpDict = Dict{KEY, GPE}()
    groupIndices = Dict{KEY, UnitRange{Int64}}()
    istart = 1
    for j in 1:ngroups
        gp = gpList[j]
        key = groupKeys[j]
        nobs = gp.nobs
        gp.nobs > 0 || throw("empty group")
        groupIndices[key] = istart:istart+nobs-1
        gpDict[key] = gp
        istart += nobs
    end
    y = vcat([gp.y for gp in gpList]...)
    mgpcv = MultiGPCovars{KEY}(
                D, y, groupKeys, gpDict, groupIndices,
                p, dim, total_nobs, logNoise, mean, kern, βkern)
    return mgpcv
end
function propagate_params!(mgpcv::MultiGPCovars)
    for gp in values(mgpcv.mgp)
        # harmonize parameters
        gp.kernel = mgpcv.kernel
        gp.mean = mgpcv.mean
        gp.logNoise = mgpcv.logNoise
    end
end

function update_mll!(mgpcv::MultiGPCovars)
    propagate_params!(mgpcv)
    cK = mgpcv.cK.mat
    cov!(cK, mgpcv.βkern, mgpcv.D, mgpcv.βdata)
    μ = Array{Float64}(undef, mgpcv.nobs)
    for key in mgpcv.groupKeys
        gp = mgpcv.mgp[key]
        slice = mgpcv.groupIndices[key]
        μ[slice] = mean(mgpcv.mean,gp.x)
        addcov!(view(cK, slice, slice), mgpcv.kernel, gp.x, gp.data)
    end
    add_diag!(cK, max(exp(2*mgpcv.logNoise),1e-8))
    mgpcv.cK = update_chol!(mgpcv.cK)
    mgpcv.alpha = mgpcv.cK \ (mgpcv.y - μ)
    mgpcv.mll = -dot((mgpcv.y - μ),mgpcv.alpha)/2.0 - logdet(mgpcv.cK)/2.0 - mgpcv.nobs*log(2π)/2.0 # Marginal log-likelihood
end

function initialise_mll!(mgpcv::MultiGPCovars)
    cK = Array{Float64}(undef, mgpcv.nobs, mgpcv.nobs)
    propagate_params!(mgpcv)
    cov!(cK, mgpcv.βkern, mgpcv.D, mgpcv.βdata)
    μ = Array{Float64}(undef, mgpcv.nobs)
    for key in mgpcv.groupKeys
        gp = mgpcv.mgp[key]
        slice = mgpcv.groupIndices[key]
        μ[slice] = mean(mgpcv.mean, gp.x)
        addcov!(view(cK, slice, slice), mgpcv.kernel, gp.x, gp.data)
    end
    add_diag!(cK, max(exp(2*mgpcv.logNoise),1e-8))
    mgpcv.cK = PDMats.PDMat(cK)
    mgpcv.alpha = mgpcv.cK \ (mgpcv.y .- μ)
    mgpcv.mll = -dot((mgpcv.y-μ),mgpcv.alpha)/2.0 - logdet(mgpcv.cK)/2.0 - mgpcv.nobs*log(2π)/2.0
end

function update_mll_and_dmll!(
        mgpcv::MultiGPCovars, ααinvcKI::AbstractMatrix
        ;
        noise::Bool=true,  # include gradient component for the logNoise term
        domean::Bool=true, # include gradient components for the mean parameters
        kern::Bool=true,   # include gradient components for the spatial kernel parameters
        beta::Bool=true,   # include gradient components for the linear regression prior terms
    )
    update_mll!(mgpcv)
    n_mean_params = num_params(mgpcv.mean)
    n_kern_params = num_params(mgpcv.kernel)
    n_beta_params = num_params(mgpcv.βkern)
    dmll = Array{Float64}(undef, noise + domean*n_mean_params + kern*n_kern_params + beta*n_beta_params)
    logNoise = mgpcv.logNoise
    get_ααinvcKI!(ααinvcKI, mgpcv.cK, mgpcv.alpha)
    i=1
    if noise
        dmll[i] = exp(2*logNoise)*trace(ααinvcKI)
        i+=1
    end
    if domean
        Mgrads = vcat([grad_stack(gp.mean, gp.x) for gp in values(mgpcv.mgp)]...)
        for j in 1:n_mean_params
            dmll[i] = dot(Mgrads[:,j],mgpcv.alpha)
            i+=1
        end
    end
    if kern
        dmll_k = @view(dmll[i:i+n_kern_params-1])
        fill!(dmll_k, 0.0)
        for key in mgpcv.groupKeys
            gp = mgpcv.mgp[key]
            slice = mgpcv.groupIndices[key]
            ααview = view(ααinvcKI, slice, slice)
            dmll_k_gp = Vector{Float64}(undef, n_kern_params)
            dmll_kern!(dmll_k_gp, gp.kernel, gp.x, gp.data, ααview)
            dmll_k .+= dmll_k_gp
        end
        i += n_kern_params
    end
    if beta
        dmll_β = @view(dmll[i:end])
        dmll_kern!(dmll_β, mgpcv.βkern, mgpcv.D, mgpcv.βdata, ααinvcKI)
        i+=1
    end
    mgpcv.dmll = dmll
end
function get_params(mgpcv::MultiGPCovars; noise::Bool=true, domean::Bool=true, kern::Bool=true, beta::Bool=true)
    params = Float64[]
    if noise; push!(params, mgpcv.logNoise); end
    if domean;  append!(params, get_params(mgpcv.mean)); end
    if kern; append!(params,  get_params(mgpcv.kernel)); end
    if beta; append!(params,  get_params(mgpcv.βkern)); end
    return params
end
function set_params!(mgpcv::MultiGPCovars, hyp::Vector{Float64};
                    noise::Bool=true, domean::Bool=true, kern::Bool=true, beta::Bool=true)
    i=1
    if noise
        mgpcv.logNoise = hyp[i]
        i+=1
    end
    if domean
        set_params!(mgpcv.mean, hyp[i:i+num_params(mgpcv.mean)-1])
        i+=num_params(mgpcv.mean)
    end
    if kern
        set_params!(mgpcv.kernel, hyp[i:i+num_params(mgpcv.kernel)-1])
        i+=num_params(mgpcv.kernel)
    end
    if beta
        set_params!(mgpcv.βkern, hyp[i:i+num_params(mgpcv.βkern)-1])
        i+=num_params(mgpcv.βkern)
    end
    propagate_params!(mgpcv)
end

function optimize!(mgpcv::MultiGPCovars; noise::Bool=true, domean::Bool=true, kern::Bool=true, beta::Bool=true,
                    method=ConjugateGradient(), options=Optim.Options())
    cK_buffer = Array{Float64}(undef, mgpcv.nobs, mgpcv.nobs)
    function mll(hyp::Vector{Float64})
        try
            set_params!(mgpcv, hyp; noise=noise, domean=domean, kern=kern, beta=beta)
            update_mll!(mgpcv)
            return -mgpcv.mll
        catch err
             if !all(isfinite.(hyp))
                println(err)
                return Inf
            elseif isa(err, ArgumentError)
                println(err)
                return Inf
            elseif isa(err, LinearAlgebra.PosDefException)
                println(err)
                return Inf
            else
                throw(err)
            end
        end
    end

    function mll_and_dmll!(grad::Vector{Float64}, hyp::Vector{Float64})
        try
            set_params!(mgpcv, hyp; noise=noise, domean=domean, kern=kern, beta=beta)
            update_mll_and_dmll!(mgpcv, cK_buffer; noise=noise, domean=domean, kern=kern, beta=beta)
            grad[:] = -mgpcv.dmll
            return -mgpcv.mll
        catch err
             if !all(isfinite.(hyp))
                println(err)
                return Inf
            elseif isa(err, ArgumentError)
                println(err)
                return Inf
            elseif isa(err, LinearAlgebra.PosDefException)
                println(err)
                return Inf
            else
                throw(err)
            end
        end
    end
    function dmll!(grad::Vector{Float64}, hyp::Vector{Float64})
        mll_and_dmll!(grad, hyp)
    end

    init = get_params(mgpcv;  noise=noise, domean=domean, kern=kern, beta=beta)  # Initial hyperparameter values
    func = Optim.OnceDifferentiable(mll, dmll!, mll_and_dmll!, init)
    results = Optim.optimize(func,init, method, options)
    set_params!(mgpcv, Optim.minimizer(results); noise=noise, domean=domean, kern=kern, beta=beta)
    mll(Optim.minimizer(results))
    return results
end

"""
    Obtain the full (block diagonal) spatial covariance matrix
"""
function spatial_cov(mgpcv::MultiGPCovars)
    cK = zeros(mgpcv.nobs, mgpcv.nobs)
    for key in mgpcv.groupKeys
        gp = mgpcv.mgp[key]
        slice = mgpcv.groupIndices[key]
        addcov!(view(cK, slice, slice), mgpcv.kernel, gp.x, gp.data)
    end
    return cK
end

"""
    Obtain the covariance matrix of the outcome conditional
    on the β coefficients.
    cov(Y ∣ β)
"""
function get_ΣYβ(mgpcv::MultiGPCovars)
    cK = spatial_cov(mgpcv)
    for i in 1:mgpcv.nobs
        cK[i,i] += max(exp(2*mgpcv.logNoise),1e-8)
    end
    return PDMats.PDMat(cK)
end

"""
    Evaluate the mean function.
"""
function mean(mgpcv::MultiGPCovars)
    μ = Array{Float64}(undef, mgpcv.nobs)
    for key in mgpcv.groupKeys
        gp = mgpcv.mgp[key]
        slice = mgpcv.groupIndices[key]
        μ[slice] = mean(mgpcv.mean,gp.x)
    end
    return μ
end

"""
    Get the posterior mean coefficient value E(β ∣ Y)
    Notation:
        ΣYβ = cov(Y ∣ β)
        Σβ = cov(β) = σβ² I # prior variance for β

    precision = Dᵀ ΣYβ⁻¹ D + Σβ⁻¹
    βhat = precision-weighted average
         = (precision⁻¹ D) (ΣYβ⁻¹ (Y-μ))

    See:
        Notes on “Analytical covariates”
        BDA3 equation 14.4
"""
function postmean_β(mgpcv::MultiGPCovars)
    ΣY_β = get_ΣYβ(mgpcv)
    precision = PDMats.X_invA_Xt(ΣY_β, mgpcv.D)
    for i in 1:mgpcv.p
        precision[i,i] += mgpcv.βkern.ℓ2
    end
    m = mean(mgpcv)
    βhat = (precision \ mgpcv.D) * (ΣY_β \ (mgpcv.y.-m))
    return βhat
end


text:
push!(LOAD_PATH,"../src/")

using Documenter, LatticeUtilities

DocMeta.setdocmeta!(LatticeUtilities, :DocTestSetup, :(using LatticeUtilities); recursive=true)

makedocs(
    sitename = "LatticeUtilities.jl",
    modules  = [LatticeUtilities],
    pages    = [
        "Home"            => "index.md",
        "Getting Started" => "getting_started.md",
        "Examples"        => "examples.md",
        "API"             => "api.md"
    ],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    )
)

deploydocs(
    repo = "github.com/cohensbw/LatticeUtilities.jl.git"
)


text:
module SessionActions

import ..Pluto: ServerSession, ClientSession, Notebook, emptynotebook, tamepath, new_notebooks_directory, without_dotjl, numbered_until_new, readwrite, move_notebook!, update_save_run!, putnotebookupdates!, putplutoupdates!, load_notebook, clientupdate_notebook_list, WorkspaceManager, @asynclog

struct NotebookIsRunningException <: Exception
    notebook::Notebook
end

abstract type AbstractUserError <: Exception end
struct UserError <: AbstractUserError
    msg::String
end
function Base.showerror(io::IO, e::UserError)
    print(io, e.msg)
end

function open_url(session::ServerSession, url::AbstractString; kwargs...)
    path = download(url, emptynotebook().path)
    open(session, path; kwargs...)
end

function open(session::ServerSession, path::AbstractString; run_async=true, compiler_options=nothing, as_sample=false)
    if as_sample
        new_filename = "sample " * without_dotjl(basename(path))
        new_path = numbered_until_new(joinpath(new_notebooks_directory(), new_filename); suffix=".jl")

        readwrite(path, new_path)
        path = new_path
    end

    for nb in values(session.notebooks)
        if realpath(nb.path) == realpath(tamepath(path))
            throw(NotebookIsRunningException(nb))
        end
    end

    nb = load_notebook(tamepath(path), session.options.evaluation.run_notebook_on_load)

    # overwrites the notebook environment if specified
    if compiler_options !== nothing
        nb.compiler_options = compiler_options
    end

    session.notebooks[nb.notebook_id] = nb
    if session.options.evaluation.run_notebook_on_load
        update_save_run!(session, nb, nb.cells; run_async=run_async, prerender_text=true)
        # TODO: send message when initial run completed
    end

    if run_async
        @asynclog putplutoupdates!(session, clientupdate_notebook_list(session.notebooks))
    else
        putplutoupdates!(session, clientupdate_notebook_list(session.notebooks))
    end

    nb
end

function new(session::ServerSession; run_async=true)
    nb = emptynotebook()
    update_save_run!(session, nb, nb.cells; run_async=run_async, prerender_text=true)
    session.notebooks[nb.notebook_id] = nb

    if run_async
        @asynclog putplutoupdates!(session, clientupdate_notebook_list(session.notebooks))
    else
        putplutoupdates!(session, clientupdate_notebook_list(session.notebooks))
    end

    nb
end

function shutdown(session::ServerSession, notebook::Notebook; keep_in_session=false, async=false)
    if !keep_in_session
        listeners = putnotebookupdates!(session, notebook) # TODO: shutdown message
        delete!(session.notebooks, notebook.notebook_id)
        putplutoupdates!(session, clientupdate_notebook_list(session.notebooks))
        for client in listeners
            @async close(client.stream)
        end
    end
    WorkspaceManager.unmake_workspace((session, notebook); async=async)
end

function attach_client(session::ServerSession, clientid::Symbol, stream::IO)
  client = get(session.connected_clients, clientid, nothing)
  if client === nothing
    client = ClientSession(clientid, stream)
    session.connected_clients[clientid] = client
    @info "Client connected $(clientid) (total $(length(session.connected_clients)))"
  else
    client.stream = stream # it might change when the same client reconnects
  end
  client
end

function detach_client(session::ServerSession, client::ClientSession)
  @info "Client disconnected $(client.id)"
  delete!(session.connected_clients, client.id)
  if client.connected_notebook !== nothing
    notebook = client.connected_notebook

    # cleanup client-owned cells
    to_remove = []
    for (cellid, cell) in notebook.cells_dict
      if cell.owner == client.id
        push!(to_remove, cellid)
      end
    end
    if !isempty(to_remove)
      for cellid in to_remove
        delete!(notebook.cells_dict, cellid)
      end
      notebook.cell_order = filter(
        cellid -> haskey(notebook.cells_dict, cellid),
        notebook.cell_order
      )
      # TODO(andreypopp): do we want to do that here?
      # update_save_run!(session, notebook, notebook.cells; run_async=true, save=false)
    end
  end
end

end


text:
mutable struct ReplayMemory{S<:Real, A<:Real, R<:Real}
    n::Int64
    maxsize::Int64
    batchsize::Int64
    s::Matrix{S}
    a::Matrix{A}
    r::Vector{R}
    d::Vector{Bool}
    function ReplayMemory{S, A, R}(obssize, actiondim, maxsize, batchsize) where {S<:Real, A<:Real, R<:Real}
        new(
            0,
            maxsize,
            batchsize,
            zeros(S, obssize, maxsize),
            zeros(A, actiondim, maxsize),
            zeros(R, maxsize),
            zeros(Bool, maxsize)
        )
    end
end

function append!(memory::ReplayMemory, s, a, r, d)
    ind = memory.n % memory.maxsize + 1
    memory.n += 1
    memory.s[:, ind] .= s
    memory.a[:, ind] .= a
    memory.r[ind] = r
    memory.d[ind] = d
end

function sample(memory::ReplayMemory)
    inds = sample(1:(min(memory.n, memory.maxsize)-1), min(memory.n-1, memory.batchsize), replace=false)
    return (
        memory.s[:, inds],
        memory.a[:, inds],
        memory.r[inds],
        memory.s[:, inds .+ 1],
        memory.d[inds]
    )
end