text:
import Base: iterate
abstract type AbstractContinuationIterable end
abstract type AbstractContinuationState end
####################################################################################################
# Iterator interface
"""
$(TYPEDEF)

# Useful functions
- `setParam(iter, p)` set parameter with lens `iter.lens` to `p`
- `isEventActive(iter)` whether the event detection is active
- `computeEigenElements(iter)` whether to compute eigen elements
- `saveEigenvectors(iter)` whether to save eigen vectors
- `getParams(iter)` get full list of params
- `length(iter)`
- `isInDomain(iter, p)` whether `p` in is domain [pMin, pMax]. (See [`ContinuationPar`](@ref))
- `isOnBoundary(iter, p)` whether `p` in is {pMin, pMax}
"""
@with_kw_noshow struct ContIterable{TF, TJ, Tv, Tp, Tlens, T, S, E, Ttangent, Tlinear, Tplotsolution, Tprintsolution, TnormC, Tdot, Tfinalisesolution, TcallbackN, Tevent, Tfilename} <: AbstractContinuationIterable
	F::TF
	J::TJ

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

	contParams::ContinuationPar{T, S, E}

	tangentAlgo::Ttangent
	linearAlgo::Tlinear

	plot::Bool = false
	plotSolution::Tplotsolution
	recordFromSolution::Tprintsolution
	event::Tevent = nothing

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

	verbosity::Int64 = 2

	filename::Tfilename
end

# default finalizer
finaliseDefault(z, tau, step, contResult; k...) = true

# constructor
function ContIterable(Fhandle, Jhandle,
					x0, par, lens::Lens,
					contParams::ContinuationPar{T, S, E},
					linearAlgo::AbstractBorderedLinearSolver = BorderingBLS(DefaultLS());
					filename = "branch-" * string(Dates.now()),
					tangentAlgo = SecantPred(),
					plot = false,
					plotSolution = (x, p; kwargs...) -> nothing,
					recordFromSolution = (x, p; kwargs...) -> norm(x),
					normC = norm,
					dotPALC = (x,y) -> dot(x,y) / length(x),
					finaliseSolution = finaliseDefault,
					callbackN = cbDefault,
					event = nothing,
					verbosity = 0, kwargs...
					) where {T <: Real, S, E}

	return ContIterable(F = Fhandle, J = Jhandle,
				x0 = x0, par = par, lens = lens,
				contParams = contParams,
				tangentAlgo = tangentAlgo,
				linearAlgo = linearAlgo,
				plot = plot,
				plotSolution = plotSolution,
				recordFromSolution = recordFromSolution,
				normC = normC,
				dottheta = DotTheta(dotPALC),
				finaliseSolution = finaliseSolution,
				callbackN = callbackN,
				event = event,
				verbosity = verbosity,
				filename = filename)
end

Base.eltype(it::ContIterable{TF, TJ, Tv, Tp, Tlens, T, S, E, Ttangent, Tlinear, Tplotsolution, Tprintsolution, TnormC, Tdot, Tfinalisesolution, TcallbackN, Tevent, Tfilename}) where {TF, TJ, Tv, Tp, Tlens, T, S, E, Ttangent, Tlinear, Tplotsolution, Tprintsolution, TnormC, Tdot, Tfinalisesolution, TcallbackN, Tevent, Tfilename} = T

setParam(it::ContIterable{TF, TJ, Tv, Tp, Tlens, T, S, E, Ttangent, Tlinear, Tplotsolution, Tprintsolution, TnormC, Tdot, Tfinalisesolution, TcallbackN, Tevent, Tfilename}, p0::T) where {TF, TJ, Tv, Tp, Tlens, T, S, E, Ttangent, Tlinear, Tplotsolution, Tprintsolution, TnormC, Tdot, Tfinalisesolution, TcallbackN, Tevent, Tfilename} = set(it.par, it.lens, p0)

@inline isEventActive(it::ContIterable) = !isnothing(it.event) && it.contParams.detectEvent > 0
@inline computeEigenElements(it::ContIterable) = computeEigenElements(it.contParams) || (isEventActive(it) && computeEigenElements(it.event))
@inline saveEigenvectors(it::ContIterable) = saveEigenvectors(it.contParams)

@inline getParams(it::ContIterable) = it.contParams
Base.length(it::ContIterable) = it.contParams.maxSteps
@inline isInDomain(it::ContIterable, p) = it.contParams.pMin < p < it.contParams.pMax
@inline isOnBoundary(it::ContIterable, p) = (it.contParams.pMin == p) || (p == it.contParams.pMax)
# clamp p value
clampPredp(p::Number, it::AbstractContinuationIterable) = clamp(p, it.contParams.pMin, it.contParams.pMax)
####################################################################################################
"""
	state = ContState(ds = 1e-4,...)

Returns a variable containing the state of the continuation procedure. The fields are meant to change during the continuation procedure.

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

# Useful functions
- `copy(state)` returns a copy of `state`
- `copyto!(dest, state)` returns a copy of `state`
- `getSolution(state)` returns the current solution (x, p)
- `getx(state)` returns the x component of the current solution
- `getp(state)` returns the p component of the current solution
- `getpreviousp(state)` returns the p component of the previous solution
- `isStable(state)` whether the current state is stable
"""
@with_kw_noshow mutable struct ContState{Tv, T, Teigvals, Teigvec, Tcb} <: AbstractContinuationState
	z_pred::Tv								# predictor solution
	tau::Tv									# tangent predictor
	z_old::Tv								# current solution

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

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

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

	# the following values encode the current, previous number of unstable (resp. imaginary) eigen values
	# it is initialized as -1 when unknown
	n_unstable::Tuple{Int64,Int64}  = (-1, -1)	# (current, previous)
	n_imag::Tuple{Int64,Int64} 		= (-1, -1)	# (current, previous)

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

	eventValue::Tcb = nothing
end

function Base.copy(state::ContState)
	return ContState(
		z_pred 	= _copy(state.z_pred),
		tau = _copy(state.tau),
		z_old 	= _copy(state.z_old),
		isconverged = state.isconverged,
		itnewton 	= state.itnewton,
		step 		= state.step,
		ds 			= state.ds,
		theta 		= state.theta,
		stopcontinuation = state.stopcontinuation,
		stepsizecontrol  = state.stepsizecontrol,
		n_unstable 		 = state.n_unstable,
		n_imag 			 = state.n_imag,
		eventValue		 = state.eventValue,
		eigvals			 = state.eigvals,
		eigvecs			 = state.eigvecs # can be removed? to save memory?
	)
end

function Base.copyto!(dest::ContState, src::ContState)
		copyto!(dest.z_pred , src.z_pred)
		copyto!(dest.tau , src.tau)
		copyto!(dest.z_old , src.z_old)
		src.isconverged = src.isconverged
		src.itnewton 	= src.itnewton
		src.step 		= src.step
		src.ds 			= src.ds
		src.theta 		= src.theta
		src.stopcontinuation = src.stopcontinuation
		src.stepsizecontrol  = src.stepsizecontrol
		src.n_unstable 		 = src.n_unstable
		src.n_imag 			 = src.n_imag
		src.eventValue		 = src.eventValue
		src.eigvals			 = src.eigvals
		src.eigvecs			 = src.eigvecs # can be removed? to save memory?
	return dest
end

# getters
getSolution(state::AbstractContinuationState) 		= state.z_old
getx(state::AbstractContinuationState) = state.z_old.u
@inline getp(state::AbstractContinuationState) = state.z_old.p
@inline getpreviousp(state::AbstractContinuationState) = state.z_pred.p
@inline isStable(state::AbstractContinuationState) = state.n_unstable[1] == 0
@inline stepsizecontrol(state::AbstractContinuationState) = state.stepsizecontrol
####################################################################################################
# condition for halting the continuation procedure (i.e. when returning false)
@inline done(it::ContIterable, state::ContState) =
			(state.step <= it.contParams.maxSteps) &&
			(isInDomain(it, state.z_old.p) || state.step == 0) &&
			(state.stopcontinuation == false)

function getStateSummary(it, state)
	x = getx(state); p = getp(state)
	pt = it.recordFromSolution(x, p)
	stable = computeEigenElements(it) ? isStable(state) : nothing
	return mergefromuser(pt, (param = p, itnewton = state.itnewton, itlinear = state.itlinear, ds = state.ds, theta = state.theta, n_unstable = state.n_unstable[1], n_imag = state.n_imag[1], stable = stable, step = state.step))
end

function updateStability!(state::ContState, n_unstable::Int, n_imag::Int)
	state.n_unstable = (n_unstable, state.n_unstable[1])
	state.n_imag = (n_imag, state.n_imag[1])
end

function save!(br::ContResult, it::AbstractContinuationIterable, state::AbstractContinuationState)
	# update branch field
	push!(br.branch, getStateSummary(it, state))
	# save solution
	if it.contParams.saveSolEveryStep > 0 && (modCounter(state.step, it.contParams.saveSolEveryStep) || ~done(it, state))
		push!(br.sol, (x = copy(getx(state)), p = getp(state), step = state.step))
	end
	# save eigen elements
	if computeEigenElements(it)
		if mod(state.step, it.contParams.saveEigEveryStep) == 0
			push!(br.eig, (eigenvals = state.eigvals, eigenvec = state.eigvecs, step = state.step))
		end
	end
end

function plotBranchCont(contres::ContResult, state::AbstractContinuationState, iter::ContIterable)
	if iter.plot && mod(state.step, getParams(iter).plotEveryStep) == 0
		return plotBranchCont(contres, getSolution(state), getParams(iter), iter.plotSolution)
	end
end

function ContResult(it::AbstractContinuationIterable, state::AbstractContinuationState)
	x0 = getx(state); p0 = getp(state)
	pt = it.recordFromSolution(x0, p0)
	eiginfo = computeEigenElements(it) ? (state.eigvals, state.eigvecs) : nothing
	return _ContResult(pt, getStateSummary(it, state), x0, setParam(it, p0), it.lens, eiginfo, getParams(it), computeEigenElements(it))
end

# function to update the state according to the event
function updateEvent!(it::ContIterable, state::ContState)
	# if the event is not active, we return false (not detected)
	if (isEventActive(it) == false) return false; end
	outcb = it.event(it, state)
	state.eventValue = (outcb, state.eventValue[1])
	# update number of positive values
	return isEventCrossed(it.event, it, state)
end
####################################################################################################
# Continuation Iterator
#
# function called at the beginning of the continuation
# used to determine first point on branch and tangent at this point
function Base.iterate(it::ContIterable; _verbosity = it.verbosity)
	# the keyword argument is to overwrite verbosity behaviour, like when locating bifurcations
	verbose = min(it.verbosity, _verbosity) > 0
	p0 = get(it.par, it.lens)
	T = eltype(it)

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

	# Get parameters
	@unpack pMin, pMax, maxSteps, newtonOptions, η, ds = it.contParams
	if !(pMin <= p0 <= pMax)
		@error "Initial parameter $p0 must be within bounds [$pMin, $pMax]"
		return nothing
	end

	# Converge initial guess
	verbose && printstyled("─────────────────  INITIAL GUESS ────────────────────", bold = true, color = :magenta)

	# we pass additional kwargs to newton so that it is sent to the newton callback
	u0, _, isconverged, _, _ = newton(it.F, it.J, it.x0, it.par, newtonOptions; normN = it.normC, callback = it.callbackN, iterationC = 0, p = p0)
	@assert isconverged "Newton failed to converge initial guess on the branch."
	verbose && (print("\n--> convergence of initial guess = ");printstyled("OK\n\n", color=:green))
	verbose && println("--> parameter = ", p0, ", initial step")
	verbose && printstyled("\n───────────────── INITIAL TANGENT ───────────────────", bold = true, color = :magenta)
	u_pred, _, isconverged, _, _ = newton(it.F, it.J,
			u0, setParam(it, p0 + ds / η), newtonOptions; normN = it.normC, callback = it.callbackN, iterationC = 0, p = p0 + ds / η)
	@assert isconverged "Newton failed to converge. Required for the computation of the initial tangent."
	verbose && (print("\n--> convergence of initial guess = ");printstyled("OK\n\n", color=:green))
	verbose && println("--> parameter = ", p0 + ds/η, ", initial step (bis)")
	return iterateFromTwoPoints(it, u0, p0, u_pred, p0 + ds / η; _verbosity = _verbosity)
end

# same as the previous function but when two (initial guesses) points  are provided
function iterateFromTwoPoints(it::ContIterable, u0, p0::T, u1, p1::T; _verbosity = it.verbosity) where T
	theta = it.contParams.theta
	ds = it.contParams.ds
	# this is the last (first) point on the branch
	z_old   = BorderedArray(_copy(u0), p0)
	# this is a predictor for the next point on the branch, we could have used z_old as well
	z_pred	= BorderedArray(_copy(u1), p1)
	tau  = _copy(z_pred)

	# compute the tangent using Secant predictor
	getTangent!(tau, z_pred, z_old, it, ds, theta, SecantPred(), _verbosity)

	# compute eigenvalues to get the type. Necessary to give a ContResult
	if computeEigenElements(it)
		eigvals, eigvecs, = computeEigenvalues(it, u0, it.par, it.contParams.nev)
		if ~saveEigenvectors(it)
			eigvecs = nothing
		end
	else
		eigvals, eigvecs = nothing, nothing
	end

	# compute event value and store into state
	cbval = isEventActive(it) ? initialize(it.event, T) : nothing # event result
	state = ContState(z_pred = z_pred, tau = tau, z_old = z_old, isconverged = true, ds = it.contParams.ds, theta = it.contParams.theta, eigvals = eigvals, eigvecs = eigvecs, eventValue = (cbval, cbval))

	# update stability
	if computeEigenElements(it)
		_, n_unstable, n_imag = isStable(getParams(it), eigvals)
		updateStability!(state, n_unstable, n_imag)
	end

	# we update the event function result
	updateEvent!(it, state)
	return state, state
end

function Base.iterate(it::ContIterable, state::ContState; _verbosity = it.verbosity)
	if !done(it, state) return nothing end
	# next line is to overwrite verbosity behaviour, like when locating bifurcations
	verbosity = min(it.verbosity, _verbosity)
	verbose = verbosity > 0
	verbose1 = verbosity > 1

	@unpack step, ds, theta = state

	# Predictor: state.z_pred. The following method only mutates z_pred
	getPredictor!(state, it)

	if verbose
		printstyled("──"^35*"\nContinuation Step $step \n", bold= true);
		@printf("Step size = %2.4e\n", ds); print("Parameter ", getLensSymbol(it.lens))
		@printf(" = %2.4e ⟶  %2.4e [guess]\n", state.z_old.p, clampPredp(state.z_pred.p, it))
	end

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

	# Successful step
	if state.isconverged
		if verbose
			verbose1 && printstyled("--> Step Converged in $(state.itnewton) Nonlinear Iteration(s)\n", color=:green)
			print("Parameter ", getLensSymbol(it.lens))
			@printf(" = %2.4e ⟶  %2.4e \n", state.z_old.p, z_newton.p)
		end

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

		# record previous parameter (cheap) and update current solution
		state.z_pred.p = state.z_old.p
		copyto!(state.z_old, z_newton)

		# Eigen-elements computation, they are stored in state
		if computeEigenElements(it)
			# this compute eigen-elements, store them in state and update the stab indices in state
			iteigen = computeEigenvalues!(it, state)
			verbose1 && printstyled(color=:green,"--> Computed ", length(state.eigvals), " eigenvalues in ", iteigen, " iterations, #unstable = ", state.n_unstable[1],"\n")
		end
		state.step += 1
	else
		verbose && printstyled("Newton correction failed\n", color=:red)
		verbose && (println("--> Newton Residuals history = ");display(fval))
	end

	# Step size control
	if ~state.stopcontinuation && stepsizecontrol(state)
		# we update the PALC parameters ds and theta, they are in the state variable
		state.ds, state.theta, state.stopcontinuation = stepSizeControl(ds, theta, it.contParams, state.isconverged, state.itnewton, state.tau, it.tangentAlgo, verbosity)
	end

	return state, state
end

function continuation!(it::ContIterable, state::ContState, contRes::ContResult)
	contParams = getParams(it)
	verbose = it.verbosity > 0
	verbose1 = it.verbosity > 1

	next = (state, state)

	while ~isnothing(next)
		# we get the current state
		_, state = next
		########################################################################################
		# the new solution has been successfully computed
		# we perform saving, plotting, computation of eigenvalues...
		# the case state.step = 0 was just done above
		if state.isconverged && (state.step <= it.contParams.maxSteps) && (state.step > 0)
			# Detection of fold points based on parameter monotony, mutates contRes.specialpoint
			# if we detect bifurcations based on eigenvalues, we disable fold detection to avoid duplicates
			if contParams.detectFold && contParams.detectBifurcation < 2
				foldetected = locateFold!(contRes, it, state)
			end

			if contParams.detectBifurcation > 1 && detectBifucation(state)
				status::Symbol = :guess
				_T = eltype(it)
				interval::Tuple{_T, _T} = getinterval(state.z_pred.p, getp(state))
				# if the detected bifurcation point involves a parameter values with is on
				# the boundary of the parameter domain, we disable bisection because it would
				# lead to infinite looping. Indeed, clamping messes up the `ds`
				if contParams.detectBifurcation > 2 && ~isOnBoundary(it, getp(state))
					verbose1 && printstyled(color=:red, "--> Bifurcation detected before p = ", getp(state), "\n")
					# locate bifurcations with bisection, mutates state so that it stays very close to the bifurcation point. It also updates the eigenelements at the current state. The call returns :guess or :converged
					status, interval = locateBifurcation!(it, state, it.verbosity > 2)
				end
				# we double-ckeck that the previous line, which mutated `state`, did not remove the bifurcation point
				if detectBifucation(state)
					_, bifpt = getBifurcationType(contParams, state, it.normC, it.recordFromSolution, it.verbosity, status, interval)
					if bifpt.type != :none; push!(contRes.specialpoint, bifpt); end
					# detect loop in the branch
					# contParams.detectLoop && (state.stopcontinuation = detectLoop(contRes, bifpt))
				end
			end

			if isEventActive(it)
				# check if an event occured between the 2 continuation steps
				eveDetected = updateEvent!(it, state)
				verbose1 && printstyled(color = :blue, "--> Event values: ", state.eventValue[2], " --> ", state.eventValue[1],"\n")
				eveDetected && (verbose && printstyled(color=:red, "--> Event detected before p = ", getp(state), "\n"))
				# save the event if detected and / or use bisection to locate it precisely
				if eveDetected
					_T = eltype(it); status = :guess; intervalevent = (_T(0),_T(0))
					if contParams.detectEvent > 1
						# interval = getinterval(state.z_pred.p, getp(state))::Tuple{_T, _T}
						status, intervalevent = locateEvent!(it.event, it, state, it.verbosity > 2)
					end
					success, bifpt = getEventType(it.event, it, state, it.verbosity, status, intervalevent)
					state.stopcontinuation |= ~success
					if bifpt.type != :none; push!(contRes.specialpoint, bifpt); end
				end
			end

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

			# Call user saved finaliseSolution function. If returns false, stop continuation
			# we put a OR to stop continuation if the stop was required before
			state.stopcontinuation |= ~it.finaliseSolution(state.z_old, state.tau, state.step, contRes; state = state, iter = it)

			# Save current state in the branch
			save!(contRes, it, state)

			# Plotting
			plotBranchCont(contRes, state, it)
		end
		########################################################################################
		# body
		next = iterate(it, state)
	end

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

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

function continuation(it::ContIterable)
	## !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	# The return type of this method, e.g. ContResult
	# is not known at compile time so we
	# need a function barrier to resolve it
	## !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

	# we compute the cache for the continuation, i.e. state::ContState
	# In this call, we also compute the initial point on the branch (and its stability) and the initial tangent
	states = iterate(it)
	isnothing(states) && return nothing, nothing, nothing

	# variable to hold the result from continuation, i.e. a branch
	contRes = ContResult(it, states[1])

	# perform the continuation
	return continuation!(it, states[1], contRes)
end

function continuation(Fhandle, Jhandle, x0, par, lens::Lens, contParams::ContinuationPar, linearAlgo::AbstractBorderedLinearSolver; bothside = false, kwargs...)
	it = ContIterable(Fhandle, Jhandle, x0, par, lens, contParams, linearAlgo; kwargs...)
	if bothside
		res1 = continuation(it)
		@set! it.contParams.ds = -contParams.ds
		res2 = continuation(it)
		contresult = _merge(res1[1],res2[1])
		return contresult, res1[2], res1[3]

	else
		return continuation(it)
	end
end

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

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

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

# Arguments:
- `F` is a function with input arguments `(x, p)`, where `p` is the set of parameters passed to `F`, and returning a vector `r` which represents the functional. For type stability, the types of `x` and `r` should match. In particular, it is not **inplace**,
- `J` is the jacobian of `F` at `(x, p)`. It can assume three forms.
    1. Either `J` is a function and `J(x, p)` returns a `::AbstractMatrix`. In this case, you can use the default arguments of `contParams::ContinuationPar` for `continuation`.
    2. Or `J` is a function and `J(x, p)` returns a function taking one argument `dx` and returning `dr` of the same type as `dx`. In our notation, `dr = J * dx`. In this case, you cannot use the default parameters of `contParams::ContinuationPar` for `continuation` and you have to use a Matrix Free linear solver, for example `GMRESIterativeSolvers`,
    3. Or `J` is a function and `J(x, p)` returns a variable `j` which can assume any type. Then, you must implement a linear solver `ls` as a composite type, subtype of `AbstractLinearSolver` which is called like `ls(j, rhs)` and which returns the solution of the jacobian linear system. See for example `examples/SH2d-fronts-cuda.jl`. This linear solver is passed to `NewtonPar(linsolver = ls)` which itself is passed to `ContinuationPar`. Similarly, you have to implement an eigensolver `eig` as a composite type, subtype of `AbstractEigenSolver`.
- `x0` initial guess,
- `par` initial set of parameters,
- `lens::Lens` specifies which parameter axis among `par` is used for continuation. For example, if `par = (α = 1.0, β = 1)`, we can perform continuation w.r.t. `α` by using `lens = (@lens _.α)`. If you have an array `par = [ 1.0, 2.0]` and want to perform continuation w.r.t. the first variable, you can use `lens = (@lens _[1])`. For more information, we refer to `SetField.jl`.
- `contParams` parameters for continuation. See [`ContinuationPar`](@ref) for more information about the options

# Optional Arguments:
- `plot = false` whether to plot the solution while computing
- `recordFromSolution = (x, p) -> norm(x)` function used record a few indicators about the solution. It could be `norm` or `(x, p) -> x[1]`. This is also useful when saving several huge vectors is not possible for memory reasons (for example on GPU...). This function can return pretty much everything but you should keep it small. For example, you can do `(x, p) -> (x1 = x[1], x2 = x[2], nrm = norm(x))` or simply `(x, p) -> (sum(x), 1)`. This will be stored in `contres.branch` (see below). Finally, the first component is used to plot in the continuation curve.
- `plotSolution = (x, p; kwargs...) -> nothing` function implementing the plot of the solution. For example, you can pass something like `(x, p; kwargs...) -> plot(x; kwargs...)`.
- `finaliseSolution = (z, tau, step, contResult; kwargs...) -> true` Function called at the end of each continuation step. Can be used to alter the continuation procedure (stop it by returning `false`), saving personal data, plotting... The notations are ``z=(x, p)``, `tau` is the tangent at `z` (see below), `step` is the index of the current continuation step and `ContResult` is the current branch. For advanced use, the current `state::ContState` of the continuation is passed in `kwargs`. Note that you can have a better control over the continuation procedure by using an iterator, see [Iterator Interface](@ref).
- `callbackN` callback for newton iterations. See docs for [`newton`](@ref). Can be used to change preconditioners
- `tangentAlgo = SecantPred()` controls the algorithm used to predict the tangents along the curve of solutions or the corrector. Can be `NaturalPred`, `SecantPred` or `BorderedPred`. See below for more information.
- `linearAlgo = BorderingBLS()`. Used to control the way the extended linear system associated to the continuation problem is solved. Can be `MatrixBLS`, `BorderingBLS` or `MatrixFreeBLS`.
- `verbosity::Int = 0` controls the amount of information printed during the continuation process. Must belong to `{0,1,2,3}`. In case `contParams.newtonOptions.verbose = false`, the following is valid (Otherwise the newton iterations are shown). Each case prints more information then the previous one:
    - case 0: print nothing
    - case 1: print basic information about the continuation: used predictor, step size and parameter values
    - case 2: print newton iterations number, stability of solution, detected bifurcations / events
    - case 3: print information during bisection to detect bifurcation / events
- `normC = norm` norm used in the different Newton solves
- `dotPALC = (x, y) -> dot(x, y) / length(x)`, dot product used to define the weighted dot product (resp. norm) ``\\|(x, p)\\|^2_\\theta`` in the constraint ``N(x, p)`` (see below). This argument can be used to remove the factor `1/length(x)` for example in problems where the dimension of the state space changes (mesh adaptation, ...)
- `filename` name of a file to save the computed branch during continuation. The identifier .jld2 will be appended to this filename
- `bothside=true` compute the branches on the two sides of `p0`, merge them and return it.

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

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

!!! tip "Continuing the branch in the opposite direction"
    Just change the sign of `ds` in `ContinuationPar`.

# Simplified call:
You can also use the following call for which the jacobian **matrix** (beware of large systems of equations!) is computed internally using Finite Differences

	continuation(Fhandle, x0, par, lens, contParams::ContinuationPar; kwargs...)
"""
function continuation(Fhandle, Jhandle, x0, par, lens::Lens, contParams::ContinuationPar;
					linearAlgo = nothing, kwargs...)
	# Create a bordered linear solver using the newton linear solver provided by the user
	if isnothing(linearAlgo)
		_linearAlgo = BorderingBLS(contParams.newtonOptions.linsolver)
	else
		# no linear solver has been specified
		if isnothing(linearAlgo.solver)
			_linearAlgo = @set linearAlgo.solver = contParams.newtonOptions.linsolver
		else
			_linearAlgo = linearAlgo
		end
	end

	return continuation(Fhandle, Jhandle, x0, par, lens, contParams, _linearAlgo; kwargs...)
end

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


text:
export DAQMeasurementProtocol, DAQMeasurementProtocolParams, sequenceName, sequence, mdf, prepareMDF

Base.@kwdef struct DAQMeasurementProtocolParams <: ProtocolParams
  sequenceName::AbstractString
end
DAQMeasurementProtocolParams(dict::Dict) = params_from_dict(DAQMeasurementProtocolParams, dict)

Base.@kwdef mutable struct DAQMeasurementProtocol <: Protocol
  name::AbstractString
  description::AbstractString
  scanner::MPIScanner
  params::DAQMeasurementProtocolParams

  sequence::Union{Sequence, Missing} = missing
  mdf::Union{MDFv2InMemory, Missing} = missing
  filename::AbstractString = ""
end

sequenceName(protocol::DAQMeasurementProtocol) = protocol.params.sequenceName
sequence(protocol::DAQMeasurementProtocol) = protocol.sequence
mdf(protocol::DAQMeasurementProtocol) = protocol.mdf

#TODO: This has currently no link to an MDF store. How should we integrate it?
function prepareMDF(protocol::DAQMeasurementProtocol, filename::AbstractString, study::MDFv2Study, experiment::MDFv2Experiment, operator::AbstractString="anonymous")
  protocol.mdf = MDFv2InMemory()
  protocol.mdf.root = defaultMDFv2Root()
  protocol.mdf.study = study
  protocol.mdf.experiment = experiment
  protocol.mdf.scanner = MDFv2Scanner(
    boreSize = ustrip(u"m", scannerBoreSize(protocol.scanner)),
    facility = scannerFacility(protocol.scanner),
    manufacturer = scannerManufacturer(protocol.scanner),
    name = scannerName(protocol.scanner),
    operator = operator,
    topology = scannerTopology(protocol.scanner)
  )

  protocol.filename = filename
end

function init(protocol::DAQMeasurementProtocol)
  scanner_ = scanner(protocol)
  configDir_ = configDir(scanner_)
  sequenceName_ = sequenceName(protocol)
  filename = joinpath(configDir_, "Sequences", "$sequenceName_.toml")
  protocol.sequence = Sequence(filename)
end

function execute(protocol::DAQMeasurementProtocol)
  scanner_ = scanner(protocol)
  seqCont = getSequenceController(scanner_)

  setupSequence(seqCont, sequence(protocol))
  startSequence(seqCont)

  # Give the acquisition thread some time, so we don't have mixed up messages when prompting for input
  waitTriggerReady(seqCont)

  numTriggers = length(acqNumFrames(seqCont.sequence))
  if numTriggers > 2
    @error "The DAQMeasurementProtocol is only designed to acquire with one trigger "*
           "(just foreground frames) or two triggers (background and foreground frames). "*
           "Please ckeck the associated sequence configuration."
  elseif numTriggers == 1
    @info "This configuration of the DAQMeasurementProtocol only reads foreground frames."
    decision = askConfirmation("Do you have everything prepared to start the measurement? "*
                               "If the answer is no, the protocol will be stopped.")
    if decision
      trigger(seqCont)
    else
      finish(seqCont)
      wait(seqCont)
      @info "Protocol stopped"
      return
    end
  elseif numTriggers == 2
    decision = askConfirmation("Do you have removed the sample in order to take the background measurement "*
                                "and do you want to start?\nIf the answer is no, the protocol will be stopped ")
    if decision
      trigger(seqCont, true)
    else
      finish(seqCont)
      wait(seqCont)
      @info "Protocol stopped"
      return
    end

    waitTriggerReady(seqCont)

    decision = askConfirmation("Do you have everything prepared to start the foreground measurement?\n"*
                                "If the answer is no, the protocol will be stopped.")
    if decision
      trigger(seqCont, false)
    else
      finish(seqCont)
      wait(seqCont)
      @info "Protocol stopped"
      return
    end
  end

  @info "All triggers applied. Now finishing sequence."
  finish(seqCont)
  wait(seqCont)
  if !ismissing(protocol.mdf)
    @info "Sequence finished. Now saving to MDF."
    fillMDF(seqCont, protocol.mdf)
    saveasMDF(protocol.filename, protocol.mdf)
  else
    @warn "No MDF defined and thus, no data is saved. If this is a mistake "*
          "please run `prepareMDF` prior to calling `runProtocol`."
  end
  @info "Protocol finished."
end

function cleanup(protocol::DAQMeasurementProtocol)

end

ext: jl


text:
# Autogenerated wrapper script for utf8proc_jll for x86_64-w64-mingw32
export libutf8proc

JLLWrappers.@generate_wrapper_header("utf8proc")
JLLWrappers.@declare_library_product(libutf8proc, "libutf8proc.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libutf8proc,
        "bin\\libutf8proc.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Base.Test

const ROLLER_PREFIX = tempname()
println("Path to ROLLER_PREFIX: $ROLLER_PREFIX")

configure(; modes = ["debug", "info", "warn", "error", "crazy"])
add_truck(Lumberjack.JsonTruck(Lumberjack.FileRoller(ROLLER_PREFIX)))

for i = 1:50000
    log("debug", "some-msg")
    log("info", "some-msg")
    log("warn", "some-msg")
    log("error", "some-msg")
    log("crazy", "some-msg")
end

@test success(`rm $ROLLER_PREFIX.0001`)
@test success(`rm $ROLLER_PREFIX.0002`)


text:
"""
    abbrev_symbol(s::Symbol, n::Int=4)

Abbreviate the symbol `s` to an `n`-character string.
"""
function abbrev_symbol(s::Symbol, n::Int=4)
    str = string(s)
    if length(str) ≤ n
        return str
    else
        return str[1:n]
    end
end

function sorted_list_insert!{T}(v::Vector{T}, x::T)
    insert_index = 1
    for val in v
        if x<val
            break
        else
            insert_index += 1
        end
    end
    insert!(v,insert_index,x)
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 7a8e5f68-beb1-41d4-af1f-4b89e5bbf388
begin
	using PlutoUI
	md"🦆 : $(@bind ducksize Slider(20:500, default=200))"
end

# ╔═╡ e703ea17-554a-4932-a20c-5e5975f4372a
begin
using HypertextLiteral
md"""Select your **settings** below:\
*Background color*: $(@bind bg_color html"<input type=color value='#00eeff'>")\
*Font*: $(@bind font_type Select(["cursive","serif","sans-serif","monospace","fantasy"]))\
*Font size*: $(@bind font_size NumberField(10:40, default=20))\
"""
end

# ╔═╡ b5094370-ee11-11eb-3e48-b3a77979f99a
md"""# An Introduction To Pluto
Welcome, dear reader, to this short introduction to Pluto. Knowing how to use Pluto is your first step towards working in an interactive environment with the `Julia` language!

In the rest of this notebook you'll find some handy **tips** and **tricks** to get you up and running in no time.

Good luck!
"""

# ╔═╡ 2dd30617-8deb-4762-9cc3-55f6ef341505
md"""## The Very Beginning
If you're here that means that you have already done this part!
But just for completion let us repeat what you need to do to start a `Pluto` notebook.
1. Open your `Julia` PERL (the application you downloaded from the Julia site)
2. Type `using Pluto` : This will load the Pluto module for you.
3. Type `Pluto.run()` : This starts your Pluto environment in your standard browser.
4. Select the notebook you want to open by typing the **path** to it or by inserting the **url** to the notebook.
"""

# ╔═╡ bc216d83-20ba-4efb-b201-bff73a83f5af
md"""## Interacting With Pluto
In this section we will quickly explain to you how to run notebooks and interact with them.

To begin with when you open a Pluto notebook, every cell is loaded and executed. So you may want to wait a couple of seconds before typing/working in a notebook.
"""

# ╔═╡ 8d9a48e4-382f-49a1-9f51-ab549a5b7039
md"""If you want to load a cell yourself you have two options. The first one is to just click the grey ▶ button in the bottom right corner. The second way is by pressing `Ctrl + Shift`. This will cause your current cell to be executed but it will also cause all the *related* cells to update. This way Pluto always keep your whole notebook up to date. Try to run the cell below. See how the other cell changes its output as well?"""

# ╔═╡ 288f2701-aa8e-4a81-b334-b58bcdc74ae7
duckcount = 2

# ╔═╡ 5bb6c26d-4d0f-47b4-8a47-e7d436be4892
repeat("🦆",duckcount)

# ╔═╡ 75f401fc-03eb-465a-a218-68b6ccc544d1
md"""Want to see the code behind a certain cell? Easy, just click on the 👁 button to the left of this cell. Go on! Try it with this cell."""

# Congrats! See how easy this is?

# ╔═╡ 22a6d095-e48f-4a56-9587-bca63beedf99
md"""If you're feeling ready to start typing yourself then you can create a new cell by clicking on the `+` icon on the upper or lower left corner of this cell.

Lastly we want to mention that if you want to remove a cell or just disable it temporary, you can so so by clicking on the `...` circle on the left side of your cell. Not seeing this button? That's normal! You'll only see this symbol when you are inside a code cell. (remember you can do this with the 👁 button).
"""

# ╔═╡ 1bffea5a-dd8d-4156-88d4-b6b614436881
md"""## Markdown and Pluto
You may have noticed that there is a bit of markup inside this notebook. The nice thing in Pluto is that it's really easy to do this. You can just use the `Markdown` language inside these notebooks. To start a `Markdown` part, just type `md"your text here"`. Want to type a big chunck of markdown text? Use `md\"""Your big chunk of text\"""` this way you can also use backspaces inside your markdown.

Not really familiar with the `Markdown` language? Don't worry it is really easy to learn, all the info you'll ever need can be found [here](https://www.markdownguide.org/cheat-sheet/).

Another nice thing to know is that you can use `Latex` and `unicode` symbols in Pluto. `Latex` gives you easy access to symbols like `α` (just type `\alpha<TAB>`). `Unicode` makes it so you can use all the smileys you love in your text and even as variable names! (big 🧠 time).
Again if you don't know any Latex or Unicode symbols you can just look them up on there respective pages. Just click [here](https://oeis.org/wiki/List_of_LaTeX_mathematical_symbols) for **Latex** and [here](https://unicode.org/emoji/charts/full-emoji-list.html) for **Unicode**.
"""

# ╔═╡ 4b0bca9a-426f-4157-935d-cf67cd03f201
md"""## Export Modes
If you want to export your current notebook there are a couple of ways and formats to do this with.

### Formats
There are 3 different formats as which you can export a notebook:
1. A notebook
2. A static HTML file (if you want to make static sites in Pluto)
3. A PDF file

### How to do it?
There are (at the moment) two ways to do this:
1. The most easy way to do this is by just scrolling to the top of your notebook, clicking the export button (looks like a triangle on top of a circle).

2. To automatically export your notebook to a pdf you first need to install the PlutoPDF package. This can be done with the following code `import Pkg; Pkg.add(url="https://github.com/JuliaPluto/PlutoPDF.jl")`. Next you just need to run the following code in Julia: `import PlutoPDF` followed by `PlutoPDF.pluto_to_pdf("notebook.jl")`. More info about this package can be found [here](https://github.com/JuliaPluto/PlutoPDF.jl).

"""

# ╔═╡ c6363449-3b51-4f25-ab41-c0cde7f0565c
md"""## Making It Interactive
To end the introduction, you will learn how to give your notebook that extra portion of interactivity.

Pluto allows you to use `Javascript` inside your notebook. Most people don't know `Javascript` so that's where `PlutoUI` comes to the rescue! Activate this by typing `using PlutoUI` at the top of your new notebook. By simply typing `@bind varname Slider(5:10)` you can make a nice slider. The value of the variable `varname` will now change depending on where you put the slider. Inside the `Slider` function you put `(minimal value : maximal value)`. Sliders are of course not the only thing you can use, there is a whole range interaction tools. You can learn about all of them inside the sample notebook *Sample Interactivity*. These sample notebooks can be found by clicking on the **Pluto.jl** logo at the top of your notebook and then selecting the *sample notebooks* option at the start menu.

To give you a little example, just play with the slider below and see what happens.
"""

# ╔═╡ 34971682-5962-4e58-aa18-bf155256c49f
Resource("https://upload.wikimedia.org/wikipedia/commons/thumb/2/24/Male_mallard_duck_2.jpg/1920px-Male_mallard_duck_2.jpg", :width => ducksize)

# ╔═╡ b3b85049-733d-4d05-b734-051949a61468
md"""# The Advanced Stuff
If you just want to use Pluto as you Julia environment then you know everything you need to know by now 🥳. However if you would like to know how to change the layout of your Pluto notebooks, add extra functionality, make a presentation in Pluto and how to publish it only as an interactive notebook. Then the next parts are for you!
"""

# ╔═╡ f5630f3c-b0eb-4b60-a428-da21b8a68416
md"""## Presentation Mode
Yes you are reading this correctly! Next to being very interactive and easy to learn you can also make presentations in Pluto.
This is not yet implemented in a 'nice' way but we can access it through a backdoor.
By typing the following code in an empty cell you'll create a button that allows you to enter presentation mode. `html"<button onclick='present()'>present</button>"`

When in presentation mode, all second order headings or higher will start a new slide. So to start a new slide just type `md"## cool title"` in an empty cell and you are good to go!

To exit presentation mode, just click again on the same button.
"""

# ╔═╡ ec19c690-5784-4cc4-a63b-5d2485eeafa9
html"<button onclick='present()'>present</button>"

# ╔═╡ 2a83d1e3-79c2-4a36-93b2-1bf2f89b27f3
md"## This starts a new slide"

# ╔═╡ 5a1f85e6-f9c6-44b7-a2d4-6b70d65aeaef
md"### This doesn't"

# ╔═╡ 6b20de32-f606-45d1-a800-b4cdadc33017
md"""## Javascript, HTML and CSS
As the title already tells you, it is indeed possible to tailor Pluto notebooks to your own preferences. Because you can actually use `Javascript, HTML and CSS` in Pluto notebooks. To do this you just install the module called `HypertextLiteral` by typing `using HypertextLiteral` at the beginning of your notebook. This allows you to use the `@htl()` macro in which you can start typing your Hypertext code encapsulated by \""".  If you would like some more info about this or a sneakpeak of its applications, then take a look at the sample notebook called *sample JavaScript* that you can find at the sample notebook section.

Below you will find an example of what you can do by combining these languages in Pluto.
"""

# ╔═╡ d1739e3b-6363-4a77-8d7c-61b6ba8978bd
@htl("""

		This is a little example of how you can use the automatic update function of Pluto to your advantage. By binding the right variables you can completely let the user choose how their notebook looks! Isn't that great?

""")

# ╔═╡ 3acc7d2e-7c54-423c-aac8-1f8f61963d77
md"""And this is how the above cell looks in `HTML` and `CSS` code:
```htmlmixed
<article class="learning">
	<p>
		This is a little example of how you can use the automatic update function of Pluto to your advantage. By binding the right variables you can completely let the user choose how their notebook looks! Isn't that great?
	</p>
</article>

<style>

	article.learning {
		background: $(bg_color);
		padding: 1em;
		border-radius: 5px;
		font-size: $(font_size)px;
		font-family: $(font_type);
	}
</style>
```

"""

# ╔═╡ 4c1ea155-df57-4b1b-94f7-c638e7f747c2
md"""## Publishing Your Notebook Online
Lastly we will take a small look at how you could post your interactive notebook online so that other people can also access it without using Julia.

At the moment this still isn't implemented in Pluto. It is however as we mentioned earlier possible to make a static version of your notebook.
But if you really want to have an up and running version of your notebook that everyone can access, then there are two main ways at the moment.

1. The first and easiest option is to use [Binder](https://github.com/fonsp/pluto-on-binder). This gives an easy and free way of running your notebook online. However it can take a long time to load.

2. The second way is a bit more complicated. In short you will run your notebook through `GitHub` and use a (payed) server to take care of the interactive part. If you don't have a server laying around this may not be for you. If you want to find out more about this you can read all about it [here](https://cotangent.dev/how-to-publish-pluto-jl-notebooks-online-interactively/)
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
HypertextLiteral = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
HypertextLiteral = "~0.9.0"
PlutoUI = "~0.7.9"
"""

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

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

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

[[HypertextLiteral]]
git-tree-sha1 = "72053798e1be56026b81d4e2682dbe58922e5ec9"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.0"

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

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "81690084b6198a2e1da36fcfda16eeca9f9f24e4"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.1"

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

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

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "c8abc88faa3f7a3950832ac5d6e690881590d6dc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "1.1.0"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

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

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

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

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

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

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

# ╔═╡ Cell order:
# ╟─b5094370-ee11-11eb-3e48-b3a77979f99a
# ╟─2dd30617-8deb-4762-9cc3-55f6ef341505
# ╟─bc216d83-20ba-4efb-b201-bff73a83f5af
# ╟─8d9a48e4-382f-49a1-9f51-ab549a5b7039
# ╠═288f2701-aa8e-4a81-b334-b58bcdc74ae7
# ╟─5bb6c26d-4d0f-47b4-8a47-e7d436be4892
# ╟─75f401fc-03eb-465a-a218-68b6ccc544d1
# ╟─22a6d095-e48f-4a56-9587-bca63beedf99
# ╟─1bffea5a-dd8d-4156-88d4-b6b614436881
# ╟─4b0bca9a-426f-4157-935d-cf67cd03f201
# ╟─c6363449-3b51-4f25-ab41-c0cde7f0565c
# ╟─7a8e5f68-beb1-41d4-af1f-4b89e5bbf388
# ╟─34971682-5962-4e58-aa18-bf155256c49f
# ╟─b3b85049-733d-4d05-b734-051949a61468
# ╟─f5630f3c-b0eb-4b60-a428-da21b8a68416
# ╟─ec19c690-5784-4cc4-a63b-5d2485eeafa9
# ╟─2a83d1e3-79c2-4a36-93b2-1bf2f89b27f3
# ╟─5a1f85e6-f9c6-44b7-a2d4-6b70d65aeaef
# ╟─6b20de32-f606-45d1-a800-b4cdadc33017
# ╟─e703ea17-554a-4932-a20c-5e5975f4372a
# ╟─d1739e3b-6363-4a77-8d7c-61b6ba8978bd
# ╟─3acc7d2e-7c54-423c-aac8-1f8f61963d77
# ╟─4c1ea155-df57-4b1b-94f7-c638e7f747c2
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
function distribute(f)
    function fs(dict :: Dict, xs...)
        for i in values(dict)
            fs(i, xs...)
        end
    end
    function fs(arr :: Array, xs...)
        for i in arr
            fs(i, xs...)
        end
    end
    function fs(xs...)
        f(xs...)
    end
    fs
end

function _inject(sym :: Symbol, dict :: Dict)
    dict[string(sym)] = eval(sym)
end

inject = distribute(_inject)

function combination{T}(xs :: Array{T, 1})
    if length(xs) < 2
        []
    elseif length(xs) == 2
        Array{Tuple{T, T}, 1}([(xs[1], xs[2])])
    else
        Array{Tuple{T, T}, 1}([[(xs[1], x) for x in xs[2:end]]..., combination(xs[2:end])...])
    end
end

function combination(xs)
    if length(xs) < 2
        []
    elseif length(xs) == 2
        [(xs[1], xs[2])]
    else
        [[(xs[1], x) for x in xs[2:end]]..., combination(xs[2:end])...]
    end
end

function pair_computation(func :: Function, xs)
    function compute(x)
        print(string("Result for pair ", x[1], " and ", x[2], " is ", func(x[1], x[2]), ".\n"))
    end
    foreach(compute, combination(xs))
end

max_stars_repo_name: Non-Contradiction/NextGeneration

max_issues_repo_name: Non-Contradiction/NextGeneration

max_forks_repo_name: Non-Contradiction/NextGeneration


text:
module Sources

using Adapt
using DocStringExtensions

export AbstractSourceLW,
    SourceLWNoScat,
    SourceLW2Str,
    SourceSW2Str,
    source_func_longwave,
    source_func_shortwave

"""
    AbstractSourceLW{FT,FTA1D,FTA2D}

Abstract longwave source for no-scattering and two stream longwave solvers.
"""
abstract type AbstractSourceLW{FT,FTA1D,FTA2D} end

"""
    SourceLWNoScat{FT,FTA1D,FTA2D} <: AbstractSourceLW{FT,FTA1D,FTA2D}

Longwave sources: computed at layer center, layer edges,
and at the surface for no scattering calculations

# Fields

$(DocStringExtensions.FIELDS)
"""
struct SourceLWNoScat{
    FT<:AbstractFloat,
    FTA1D<:AbstractArray{FT,1},
    FTA2D<:AbstractArray{FT,2},
} <: AbstractSourceLW{FT,FTA1D,FTA2D}
    "Planck source at layer average temperature `[W/m2]` `(nlay, ncol)`"
    lay_source::FTA2D
    "Planck source at layer edge in increasing ilay direction `[W/m2]` `(nlay, ncol)`, includes spectral weighting that accounts for state-dependent frequency to g-space mapping"
    lev_source_inc::FTA2D
    "Planck source at layer edge in decreasing ilay direction `[W/m2]` `(nlay, ncol)`, includes spectral weighting that accounts for state-dependent frequency to g-space mapping"
    lev_source_dec::FTA2D
    "Surface source `[W/m2]` `(ncol)`"
    sfc_source::FTA1D
    "temporary storage array used in no scattering GPU calculations `(nlay, ncol)`"
    src_up::FTA2D
    "temporary storage array used in no scattering GPU calculations `(nlay, ncol)`"
    src_dn::FTA2D
end
Adapt.@adapt_structure SourceLWNoScat

function SourceLWNoScat(
    ::Type{FT},
    ::Type{DA},
    nlay::I,
    ncol::I,
) where {FT<:AbstractFloat,I<:Int,DA}
    FTA1D = DA{FT,1}
    FTA2D = DA{FT,2}

    lay_source = FTA2D(undef, nlay, ncol) # lay_source
    lev_source_inc = FTA2D(undef, nlay, ncol) # lev_source_inc
    lev_source_dec = FTA2D(undef, nlay, ncol) # lev_source_dec
    sfc_source = FTA1D(undef, ncol)       # sfc_source
    src_up = FTA2D(undef, nlay, ncol) # src_up
    src_dn = FTA2D(undef, nlay, ncol) # src_dn

    return SourceLWNoScat{
        eltype(lay_source),
        typeof(sfc_source),
        typeof(lay_source),
    }(
        lay_source,
        lev_source_inc,
        lev_source_dec,
        sfc_source,
        src_up,
        src_dn,
    )
end

"""
    SourceLW2Str{FT,FTA1D,FTA2D} <: AbstractSourceLW{FT,FTA1D,FTA2D}

Longwave sources: computed at layer center, layer edges,
and at the surface for 2-stream calculations

# Fields

$(DocStringExtensions.FIELDS)
"""
struct SourceLW2Str{
    FT<:AbstractFloat,
    FTA1D<:AbstractArray{FT,1},
    FTA2D<:AbstractArray{FT,2},
} <: AbstractSourceLW{FT,FTA1D,FTA2D}
    "Planck source at layer average temperature `[W/m2]` `(nlay, ncol)`"
    lay_source::FTA2D
    "Planck source at layer edge in increasing ilay direction `[W/m2]` `(nlay, ncol)`, includes spectral weighting that accounts for state-dependent frequency to g-space mapping"
    lev_source_inc::FTA2D
    "Planck source at layer edge in decreasing ilay direction `[W/m2]` `(nlay, ncol)`, includes spectral weighting that accounts for state-dependent frequency to g-space mapping"
    lev_source_dec::FTA2D
    "Surface source `[W/m2]` `(ncol)`"
    sfc_source::FTA1D
    "level source `[W/m2]` `(nlay+1, ncol)`, used in 2 stream calculations"
    lev_source::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay, ncol)`"
    src_up::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay, ncol)`"
    src_dn::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay, ncol)`"
    Rdif::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay, ncol)`"
    Tdif::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay + 1, ncol)`"
    albedo::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay + 1, ncol)`"
    src::FTA2D
end
Adapt.@adapt_structure SourceLW2Str

function SourceLW2Str(
    ::Type{FT},
    ::Type{DA},
    nlay::I,
    ncol::I,
) where {FT<:AbstractFloat,I<:Int,DA}
    FTA1D = DA{FT,1}
    FTA2D = DA{FT,2}

    lay_source = FTA2D(undef, nlay, ncol) # lay_source
    lev_source_inc = FTA2D(undef, nlay, ncol) # lev_source_inc
    lev_source_dec = FTA2D(undef, nlay, ncol) # lev_source_dec
    sfc_source = FTA1D(undef, ncol) # sfc_source
    lev_source = FTA2D(undef, nlay + 1, ncol) # lev_source
    src_up = FTA2D(undef, nlay, ncol) # src_up
    src_dn = FTA2D(undef, nlay, ncol) # src_dn
    Rdif = FTA2D(undef, nlay, ncol) # Rdif
    Tdif = FTA2D(undef, nlay, ncol) # Tdif
    albedo = FTA2D(undef, nlay + 1, ncol) # albedo
    src = FTA2D(undef, nlay + 1, ncol) # src

    return SourceLW2Str{
        eltype(lay_source),
        typeof(sfc_source),
        typeof(lay_source),
    }(
        lay_source,
        lev_source_inc,
        lev_source_dec,
        sfc_source,
        lev_source,
        src_up,
        src_dn,
        Rdif,
        Tdif,
        albedo,
        src,
    )
end

"""
    source_func_longwave(
        ::Type{FT},
        ncol::I,
        nlay::I,
        OPC::Symbol,
        ::Type{DA},
    ) where {FT<:AbstractFloat,I<:Int,DA}

Initializes the longwave source for one scalar and two stream simulations.
"""
function source_func_longwave(
    ::Type{FT},
    ncol::I,
    nlay::I,
    OPC::Symbol,
    ::Type{DA},
) where {FT<:AbstractFloat,I<:Int,DA}
    lay_source = DA{FT,2}(undef, nlay, ncol)
    lev_source_inc = DA{FT,2}(undef, nlay, ncol)
    lev_source_dec = DA{FT,2}(undef, nlay, ncol)
    sfc_source = DA{FT,1}(undef, ncol)
    src_up = DA{FT,2}(undef, nlay, ncol)
    src_dn = DA{FT,2}(undef, nlay, ncol)

    if OPC === :OneScalar
        return SourceLWNoScat{
            eltype(lay_source),
            typeof(sfc_source),
            typeof(src_up),
        }(
            lay_source,
            lev_source_inc,
            lev_source_dec,
            sfc_source,
            src_up,
            src_dn,
        )
    else
        lev_source = DA{FT,2}(undef, nlay + 1, ncol)
        Rdif = DA{FT,2}(undef, nlay, ncol)
        Tdif = DA{FT,2}(undef, nlay, ncol)
        albedo = DA{FT,2}(undef, nlay + 1, ncol)
        src = DA{FT,2}(undef, nlay + 1, ncol)
        return SourceLW2Str{FT,typeof(sfc_source),typeof(src_up)}(
            lay_source,
            lev_source_inc,
            lev_source_dec,
            sfc_source,
            lev_source,
            src_up,
            src_dn,
            Rdif,
            Tdif,
            albedo,
            src,
        )
    end
end

"""
    SourceSW2Str{FT,FTA1D,FTA2D}

Shortwave sources: computed at layer center, layer edges,
and at the surface for 2-stream calculations

# Fields

$(DocStringExtensions.FIELDS)
"""
struct SourceSW2Str{
    FT<:AbstractFloat,
    FTA1D<:AbstractArray{FT,1},
    FTA2D<:AbstractArray{FT,2},
}
    "surface source `(ncol)`"
    sfc_source::FTA1D
    "albedo `(nlay + 1, ncol)`"
    albedo::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay, ncol)`"
    src_up::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay, ncol)`"
    src_dn::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay, ncol)`"
    Rdif::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay, ncol)`"
    Tdif::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay, ncol)`"
    Rdir::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay, ncol)`"
    Tdir::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay, ncol)`"
    Tnoscat::FTA2D
    "temporary storage array, used in 2 stream calculations `(nlay + 1, ncol)`"
    src::FTA2D
end
Adapt.@adapt_structure SourceSW2Str

function SourceSW2Str(
    ::Type{FT},
    ::Type{DA},
    nlay::I,
    ncol::I,
) where {FT<:AbstractFloat,I<:Int,DA}
    FTA1D = DA{FT,1}
    FTA2D = DA{FT,2}

    sfc_source = FTA1D(undef, ncol)
    albedo = FTA2D(undef, nlay + 1, ncol)
    src_up = FTA2D(undef, nlay, ncol)
    src_dn = FTA2D(undef, nlay, ncol)
    Rdif = FTA2D(undef, nlay, ncol)
    Tdif = FTA2D(undef, nlay, ncol)
    Rdir = FTA2D(undef, nlay, ncol)
    Tdir = FTA2D(undef, nlay, ncol)
    Tnoscat = FTA2D(undef, nlay, ncol)
    src = FTA2D(undef, nlay + 1, ncol)

    return SourceSW2Str{eltype(sfc_source),typeof(sfc_source),typeof(albedo)}(
        sfc_source,
        albedo,
        src_up,
        src_dn,
        Rdif,
        Tdif,
        Rdir,
        Tdir,
        Tnoscat,
        src,
    )
end

"""
    source_func_shortwave(
        ::Type{FT},
        ncol::I,
        nlay::I,
        opc::Symbol,
        ::Type{DA},
    ) where {FT<:AbstractFloat,I<:Int,DA}

Initializes the shortwave source for one scalar and two stream simulations.
"""
function source_func_shortwave(
    ::Type{FT},
    ncol::I,
    nlay::I,
    opc::Symbol,
    ::Type{DA},
) where {FT<:AbstractFloat,I<:Int,DA}
    if opc == :OneScalar
        return nothing
    else
        sfc_source = DA{FT}(undef, ncol)
        albedo = DA{FT}(undef, nlay + 1, ncol)
        src_up = DA{FT}(undef, nlay, ncol)
        src_dn = DA{FT}(undef, nlay, ncol)
        Rdif = DA{FT}(undef, nlay, ncol)
        Tdif = DA{FT}(undef, nlay, ncol)
        Rdir = DA{FT}(undef, nlay, ncol)
        Tdir = DA{FT}(undef, nlay, ncol)
        Tnoscat = DA{FT}(undef, nlay, ncol)
        src = DA{FT}(undef, nlay + 1, ncol)
        return SourceSW2Str{FT,typeof(sfc_source),typeof(src_up)}(
            sfc_source,
            albedo,
            src_up,
            src_dn,
            Rdif,
            Tdif,
            Rdir,
            Tdir,
            Tnoscat,
            src,
        )
    end
end

end


text:
using Replay

instructions = raw"""
using OhMyREPL
using Printf
println("HelloWorld")
N = 10
@printf "%.2f %.2f %.2f\n" sin(pi/4) cos(pi/4) tan(pi/4)
@show N
"""

replay(instructions, use_ghostwriter=false, julia_project=@__DIR__)

replay(instructions, use_ghostwriter=true, julia_project=@__DIR__)


text:
# License is MIT: https://github.com/JuliaString/Unicode_Entities/LICENSE.md

__precompile__()

"""
# Public API (nothing is exported)

* lookupname(str)
* matchchar(char)
* matches(str)
* longestmatches(str)
* completions(str)
"""
module Unicode_Entities

using StrTables

VER = UInt32(1)

struct PackedNames{T,S,O} <: AbstractPackedTable{T}
    offsetvec::Vector{O}
    namtab::Vector{S}
    wrd1::StrTable{T}           # This has sorted words for 1-byte
    wrd2::StrTable{T}           # This has sorted words for 2-byte
end

Base.getindex(tab::PackedNames, ind::Integer) =
    _unpackword(tab.namtab[tab.offsetvec[ind] + 1 : tab.offsetvec[ind+1]], tab.wrd1, tab.wrd2)

struct Unicode_Table{T} <: AbstractEntityTable
    ver::UInt32
    tim::String
    inf::String
    base32::UInt32
    nam::PackedNames{T}
    ind::Vector{UInt16}
    val16::Vector{UInt16}
    ind16::Vector{UInt16}
    val32::Vector{UInt16}
    ind32::Vector{UInt16}
end

function __init__()
    (ver, tim, inf, base32, nam, ind, wrd1, wrd2, val16, ind16, val32, ind32) =
        StrTables.load(joinpath(@__DIR__, "../data/unicode.dat"))
    global default = Unicode_Table(ver, tim, inf, base32,
                                   PackedNames(nam.offsetvec, nam.namtab, wrd1, wrd2),
                                   ind, val16, ind16, val32, ind32)
end

"""
Internal function to unpack the packed Unicode entity names

Unicode has a very limited character set for the entity names, 0-9, A-Z, -, (, and ).
( and ) are only used for 4 names, so those are special cased to store in the word table.

0 represents a space, when needed between characters
1 represents a hyphen
2-11 represents the digits 0-9
12-37 represents the letters A-Z
38-53 represents up to 16*256 words stored in the wrd2 table
54-255 represent 202 words stored in the wrd1 table
"""
function _unpackword(v::Vector{UInt8}, w1, w2)
    io = IOBuffer()
    pos = 0
    len = length(v)
    prevw = false
    while (pos += 1) <= len
        ch = v[pos]
        if ch < 0x26 # single character (space, hyphen, 0-9, A-Z)
            if ch == 0x00
                write(io, ' ')
            elseif ch == 0x01
                write(io, '-')
            else
                pos != 1 && prevw && write(io, ' ')
                write(io, ch + (ch < 0x0c ? 0x2e : 0x35))
            end
            prevw = false
        else
            pos != 1 && (prevw || v[pos-1]>0x1) && write(io, ' ')
            if ch < 0x36
                write(io, w2[((((ch - 0x25)%UInt16)<<8) | v[pos+=1])-255])
            else
                write(io, w1[ch - 0x35])
            end
            prevw = true
        end
    end
    String(take!(io))
end

## Override methods

function StrTables._get_str(tab::Unicode_Table, ind)
    string(Char(ind <= tab.base32 ? tab.val16[ind] : tab.val32[ind - tab.base32] + 0x10000))
end

function StrTables.lookupname(tab::Unicode_Table, str::AbstractString)
    rng = searchsorted(tab.nam, uppercase(str))
    isempty(rng) ? StrTables._empty_str : StrTables._get_str(tab, tab.ind[rng.start])
end

StrTables.matches(tab::Unicode_Table, vec::Vector{T}) where {T} =
    length(vec) == 1 ? matchchar(tab, vec[1]) : StrTables._empty_str_vec

StrTables.longestmatches(tab::Unicode_Table, vec::Vector{T}) where {T} =
    isempty(vec) ? StrTables._empty_str_vec : matchchar(tab, uppercase(vec[1]))

function StrTables.completions(tab::Unicode_Table, str::AbstractString)
    up = uppercase(str)
    [nam for nam in tab.nam if startswith(nam, up)]
end

end # module Unicode_Entities


text:
#!/bin/julia
using Elf

function main()
    isempty(ARGS) && error("No file specifed.")

    bin = read(ARGS[1])

    ehdr = Elf64_Ehdr(bin)

    !iself(ehdr) && error("Not an ELF binary.")

    !(elfclass(ehdr) == ELFCLASS64) && error("Unknown class.")

    !(endian(ehdr) == ELFDATA2MSB) && error("Unknown endian.")

    # Print ELF file type
    println("ELF Type: $(elftype(ehdr))")

    # Print Section Names
    println("Sections:")
    for (i, s) in enumerate(collect(sections(bin)))
        println("\t[$i]:\t\"$(s.first)\"$(" " ^ (25 - length(s.first)))type: $(SHType(s.second.sh_type))")
    end
    println()

    # Print Segments Names
    println("Segments:")
    for (i, s) in enumerate(collect(segments(bin)))
        println("\t[$i]:\t$(s.first)")
    end

    # Print Symbols Names
    println("Symbols:")
    syms, rels = symbols(bin)

    for (i, s) in enumerate(collect(syms))
        println("\t[$i]:\t\"$(s.first)\"$(" " ^ (50 - length(s.first)))type: $(elf64_st_type(s.second.st_info))")
    end

    println()

    if isempty(rels)
        println("No relocations.")
    else
        println("Relocations:")
        for (i, s) in enumerate(collect(rels))
            println("\t[$i]:\t$(s.first)")
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
import ACE.OrthPolys: TransformedPolys, OrthPolyBasis
import ACE.Transforms: GetNorm

@doc raw"""
`Rn1pBasis`

One-particle basis of the form $R_n(r_{ij})$, i.e.,
no dependence on species or on $l$.

This does two things: (1) translates the `TransformedPolys` into a valid
one-particle basis; and (2) treat it as having vectorial input, i.e. value is
scalar but gradient is vectorial.

The default symbols are `:rr` for the state and `:n` for the index of the
basis function.
"""
function Rn1pBasis(  R::OrthPolyBasis, trans = nothing;
                     varsym = :rr, nsym = :n, label="R$nsym", )
   spec = [ NamedTuple{(nsym,)}((i,)) for i = 1:length(R) ]
   degrees = collect(0:(length(R)-1))
   R1 = chain(λ("rr -> norm(rr)"), trans, R)
   Rn = B1pComponent(R1, GetVal{varsym}(), spec, degrees, label)
   rl = ACE.inv_transform(r -> Base.invokelatest(trans, r), R.tl)
   rr = ACE.inv_transform(r -> Base.invokelatest(trans, r), R.tr)
   rin, rcut = extrema([rl, rr])
   Rn.meta["rin"] = rin
   Rn.meta["rcut"] = rcut
   return Rn
end


text:
"""
    jumping_rate([type_number, ]ctqw)

Return default value of the jumping rate for the adjacency matrix. This function
is not implemented for `:laplacian` or `:normalized_laplacian`. For `:adjacency` case the greatest eigenvalue
is returned. By default `type_number` is set to `ComplexF64`.
"""
function jumping_rate(::Type{T}, ctqw::AbstractCTQW) where T<:Number
   @assert ctqw.matrix == :adjacency "Default jumping rate known for adjacency matrix only"

   if ne(ctqw.graph) == 0 # if graph is empty
      return T(1.) #does not matter
   else
      return T(1. /eigs(adjacency_matrix(ctqw.graph), nev=1)[1][1])
   end
end

function jumping_rate(ctqw::AbstractCTQW)
   jumping_rate(ComplexF64, ctqw)
end

function normalized_laplacian(g::Graph, ::Type{T}) where T<:Number
  l = laplacian_matrix(g, T)
  d = Diagonal(sqrt.(1 ./(diag(l))))
  d*l*d
end

"""
    graph_hamiltonian(type, ctqw)

Returns default evolution matrix of `type`, adjacency, laplacian matrix or normalized laplacian depending on
`ctqw` parametrization.
"""
function graph_hamiltonian(::Type{T}, ctqw::AbstractCTQW) where T<:Number
   if ctqw.matrix == :adjacency
      return -adjacency_matrix(ctqw.graph, T)
   elseif ctqw.matrix == :laplacian
      return laplacian_matrix(ctqw.graph, T)
   elseif ctqw.matrix == :normalized_laplacian
      return normalized_laplacian(ctqw.graph, T)
   else
      throw(ErrorException("Model $ctqw poorly parametrized"))
   end
end
"""
    hamiltonian_evolution(hamiltonian, initstate, runtime)

Evolve `initstate` according to `hamiltonian` for time `runtime` according to
Schrödinger equation. For dense matrices matrix exponantation is calculated. For
sparse matrices `expmv` from `Expokit` package is used.
"""
function hamiltonian_evolution(hamiltonian::AbstractMatrix{<:Number},
                               initstate::Vector{<:Number},
                               runtime::Real)
   exp(-1im*hamiltonian*runtime)*initstate
end

function hamiltonian_evolution(hamiltonian::SparseMatrixCSC{<:Number},
                               initstate::Vector{<:Number},
                               runtime::Real)
   expmv(runtime, -1im*hamiltonian, initstate)
end


text:
Base.@kwdef mutable struct ConstraintParams{T}
	ϕ::T = 1.0 	      # penalty scaling parameter
	μ0::T = 10.0 	  # initial penalty parameter
	μ_max::T = 1e8    # max penalty parameter
	λ_max::T = 1e8    # max Lagrange multiplier
    use_default::MVector{4,Bool} = @MVector ones(Bool,4)
end

function ConstraintParams(ϕ::Real, μ0::Real, μ_max::Real, λ_max::Real)
    params = promote(ϕ, μ0, μ_max, λ_max)
	ConstraintParams{eltype(params)}(params...)
end

struct ALConVal{C,P,W,V,M} <: TO.AbstractConstraintValues{C}
    con::C
    inds::Vector{Int}
    vals::Vector{V}
    jac::Matrix{M}
    viol::Vector{V}
    ∇viol::Matrix{M}
    λ::Vector{V}
    μ::Vector{V}
    λbar::Vector{V}
    active::Vector{BitVector}
    c_max::Vector{Float64}
    is_const::BitArray{2}
    iserr::Bool
    params::ConstraintParams{Float64}

    tmp::M
    ∇proj::Vector{SizedMatrix{P,P,Float64,2,Matrix{Float64}}}   # (p,p) projection Jacobian
    ∇²proj::Vector{SizedMatrix{P,P,Float64,2,Matrix{Float64}}}  # (p,p) projection "Hessian," or Jacobian of ∇Π'Π
    grad::Vector{SizedVector{W,Float64,Vector{Float64}}}    # gradient of augmented Lagrangian wrt inputs
    hess::Vector{SizedMatrix{W,W,Float64,2,Matrix{Float64}}}    # hessian of augmented Lagrangian wrt inputs
    const_hess::BitVector

	function ALConVal(n::Int, m::Int, con::TO.AbstractConstraint, inds::AbstractVector{Int},
			jac, vals, iserr::Bool=false)
		if !iserr && size(TO.gen_jacobian(con)) != size(jac[1])
			throw(DimensionMismatch("size of jac[i] $(size(jac[1])) does not match the expected size of $(size(gen_jacobian(con)))"))
		end
        params = ConstraintParams()
        p = length(con)
        P = length(vals)

        viol = deepcopy(vals)
        ∇viol = deepcopy(jac)
        λ = deepcopy(vals)
        μ = [zero(v) .+ 1.0 for v in vals]
        λbar = deepcopy(vals)
        active = [ones(Bool,p) for i = 1:P]

        ix = 1:n
        iu = n .+ (1:m)
        c_max = zeros(P)
        is_const = BitArray(undef, size(jac))

        # jac = [jac; jac[end:end,:]]  # append extra for temporary array
        tmp = zero(jac[1])

        ni = size(jac[1],2)  # size of inputs to the constraint
        ∇proj  = [SizedMatrix{p,p}(zeros(p,p)) for i = 1:P]
        ∇²proj = [SizedMatrix{p,p}(zeros(p,p)) for i = 1:P]
        grad = [SizedVector{ni}(zeros(ni)) for i = 1:P]
        hess = [SizedMatrix{ni,ni}(zeros(ni,ni)) for i = 1:P]
        const_hess = BitArray(undef, P)

        new{typeof(con), p, ni, eltype(vals), eltype(jac)}(con,
            collect(inds), vals, jac, viol, ∇viol, λ, μ, λbar, active, c_max, is_const, iserr, params,
            tmp, ∇proj, ∇²proj, grad, hess, const_hess)
    end
end

function ALConVal(n::Int, m::Int, cval::ALConVal)
	# create a ConVal for the "raw" Jacobians, if needed
	# 	otherwise return the same ConVal
	if cval.iserr
		p = length(cval.con)
		ws = TO.widths(cval.con, n, m)
		jac = [SizedMatrix{p,w}(zeros(p,w)) for k in cval.inds, w in ws]
		ALConVal(n, m, cval.con, cval.inds, jac, cval.vals, false)
	else
		return cval
	end
end

function ALConVal(n::Int, m::Int, con::TO.AbstractConstraint, inds::UnitRange{Int}, iserr::Bool=false)
	C,c = TO.gen_convals(n,m,con,inds)
	ALConVal(n, m, con, inds, C, c)
end

TO.sense(cval::ALConVal) = TO.sense(cval.con)

"""
    violation!(::ALConVal)

Calculate the constraint violations. For equality constraints, this is simply the
constraint value. For conic constraints, this is the difference between the projected
and current constraint value. Assumes the constraint values have already been computed.
"""
function violation!(cval::ALConVal)
    cone = TO.sense(cval)
    for i = 1:length(cval.inds)
        cval.viol[i] .= TO.violation(cone, SVector(cval.vals[i]))
    end
end

"""
    ∇violation!(::ALConVal)

Calculate the Jacobian of the constraint violation. Assumes the constraint value and Jacobian
have already been computed.
"""
function ∇violation!(cval::ALConVal)
    cone = TO.sense(cval)
    for i = 1:length(cval.inds)
        TO.∇violation!(cone, cval.∇viol[i], cval.jac[i], SVector(cval.vals[i]), cval.∇proj[i])
    end
end

function TO.max_violation!(cval::ALConVal)
	s = TO.sense(cval.con)
    for i in eachindex(cval.inds)
        cval.c_max[i] = TO.max_violation(s, SVector(cval.vals[i]))
    end
end

function reset_duals!(con::ALConVal)
    for i in eachindex(con.λ)
        con.λ[i] .*= 0
    end
end

function reset_penalties!(con::ALConVal)
    for i in eachindex(con.μ)
        con.μ[i] .= con.params.μ0
    end
end

function set_params!(cval::ALConVal, opts)
    p = cval.params
    p.use_default[1] && (p.ϕ = opts.penalty_scaling)
    p.use_default[2] && (p.μ0 = opts.penalty_initial)
    p.use_default[3] && (p.μ_max = opts.penalty_max)
    p.use_default[4] && (p.λ_max = opts.dual_max)
end

function set_params!(cval::ALConVal;
        penalty_initial=NaN,
        penalty_scaling=NaN,
        penalty_max=NaN,
        dual_max=NaN
    )
    p = cval.params
    if !isnan(penalty_scaling)
        p.ϕ = penalty_scaling
        p.use_default[1] = false
    end
    if !isnan(penalty_initial)
        p.μ0 = penalty_initial
        p.use_default[2] = false
    end
    if !isnan(penalty_max)
        p.μ_max = penalty_max
        p.use_default[3] = false
    end
    if !isnan(dual_max)
        p.λ_max = dual_max
        p.use_default[4] = false
    end
end

function shift_fill!(cval::ALConVal, n=1)
    shift_fill!(cval.λ, n)
    shift_fill!(cval.μ, n)
end


text:
#
# Tensor the matrix V of order n up to order N,
# where V acts on qubits: j, j + 1,... k; j < k
#
# and |ψ₀> = |0>₁ ⊗ |0>₂ ⊗ ... ⊗ |0>(N)
#
# V(k) -> I₁ ⊗ ... ⊗ I(k - 1) ⊗ V(k) ⊗ I(k + n - 1) ⊗ ... ⊗ I(N)
#

module Tensor

using LinearAlgebra

export tensor, ⊗

const C32 = Complex{Float32}

function tensor(V::Matrix{C32}, indx::Tuple{Vararg{Int}}, N::Int)
    if length(indx) == 2
        j, k = indx
        if j > k
            # Qₖ ⊗ Qⱼ -> Qⱼ ⊗ Q(j - 1)
            Ṽ = tensor(V, j, N)
            R, L = σ(j, k, N)
        else
            # Qⱼ ⊗ Qₖ -> Q(k - 1) ⊗ Qₖ
            Ṽ = tensor(V, k, N)
            R, L = σ(k - 1, j, N)
        end
        return L * Ṽ * R
    else
        i, j, k = indx
        if i > j && j > k
            Ṽ = tensor(V, i, N)
            R₁, L₁ = σ(i, k, N)
            R₂, L₂ = σ(i - 1, j - 1, N)
        elseif i < j && j > k
            if i < k
                Ṽ = tensor(V, j, N)
                R₁, L₁ = σ(j, k, N)
                R₂, L₂ = σ(j - 2, i, N)
            else
                Ṽ = tensor(V, j, N)
                R₁, L₁ = σ(j, k, N)
                R₂, L₂ = σ(j - 2, i - 1, N)
            end
        elseif i > j && j < k
            if i < k
                Ṽ = tensor(V, k, N)
                R₁, L₁ = σ(k - 1, j, N)
                R₂, L₂ = σ(j - 2, i - 1, N)
            else
                Ṽ = tensor(V, i, N)
                R₁, L₁ = σ(i, k, N)
                R₂, L₂ = σ(i - 1, j, N)
            end
        else
            Ṽ = tensor(V, k, N)
            R₁, L₁ = σ(k - 1, j, N)
            R₂, L₂ = σ(k - 2, i, N)
        end
        return L₁ * L₂ * Ṽ * R₂ * R₁
    end
end

function tensor(V::Matrix{C32}, k::Int, N::Int)
    n = Int(log(2, size(V, 1)))
    L = diagm(ones(C32, 2^(k - n)))
    R = diagm(ones(C32, 2^(N - k)))
    L ⊗ V ⊗ R
end

⊗(x, y) = kron(x, y)

const SWAP = [1 0 0 0;
              0 0 1 0;
              0 1 0 0;
              0 0 0 1]

# Q(i) -> Q(i + 1)
τ(i, N) = tensor(C32.(SWAP), i + 1, N)

# k < j => Q(k) -> Q(j)
σ(j, k, N) = begin
    if k < j
        τs = [τ(k + j - i - 1, N) for i = k:j-1];
        return (*(τs...), *(reverse(τs)...))
    else
        return (C32(1), C32(1))
    end
end

end

max_stars_repo_name: aarontrowbridge/QuIPS

max_issues_repo_name: aarontrowbridge/QuIPS

max_forks_repo_name: aarontrowbridge/QuIPS


text:
using Zygote
include("chainrules_patch.jl")

import YaoExtensions, Random

c = YaoExtensions.variational_circuit(5)
h = YaoExtensions.heisenberg(5)

function loss(h, c, θ) where N
    # the assign is nessesary!
    c = dispatch(c, fill(θ, nparameters(c)))
    reg = apply(zero_state(nqubits(c)), c)
    real(expect(h, reg))
end

reg0 = zero_state(5)
zygote_grad = Zygote.gradient(θ->loss(h, c, θ), 0.5)[1]

# check gradients
using Test
dispatch!(c, fill(0.5, nparameters(c)))
greg, gparams = expect'(h, zero_state(5)=>c)
true_grad = sum(gparams)

@test zygote_grad ≈ true_grad

# the batched version
function loss2(h, c, θ) where N
    # the assign is nessesary!
    c = dispatch(c, fill(θ, nparameters(c)))
    reg = zero_state(nqubits(c),nbatch=2)
    reg = apply(reg, c)
    sum(real(expect(h, reg)))
end

reg0 = zero_state(5)
zygote_grad2 = Zygote.gradient(θ->loss2(h, c, θ), 0.5)[1]
true_grad2 = (loss2(h, c, 0.5+1e-5) - loss2(h, c, 0.5-1e-5)) / 2e-5
@test true_grad2 ≈ zygote_grad2


text:
module FillZY
export fillz, filly

using Statistics: mean
using StaticArrays: SMatrix, MVector
using OffsetArrays
using Unitful # for ustrip and u"m"
using LinearAlgebra: norm, ⋅
using ..Constants: μ₀, ϵ₀, c₀, twopi, fourpi, tol, tdigits
using ..Layers: Layer
using ..Sheets: RWGSheet
using ..RWG: RWGData
using ..PGF: c3_calc, d3_calc
using ..Zint: zint, filljk!, vtxcrd
using ..PGF: electric_modal_sum_funcs, magnetic_modal_sum_funcs
using ..Log: @logfile

const next = (2,3,1)
const prev = (3,1,2)
const third = SMatrix{3,3}([0 3 2; 3 0 1; 2 1 0])

"""
    fillz(k0,u,layers::AbstractVector{Layer},s,ψ₁,ψ₂,metal::RWGSheet,rwgdat::RWGData) -> zmat

Fill the generalized impedance matrix for an FSS of electric current type.

## Arguments:

- `k0`: Free-space wavenumber (rad/meter).
- `u`:  Green's function smoothing factor (1/meter).
- `layers`:  An array characterizing the dielectric layers surrounding the FSS sheet.
- `s`:  An integer indexing the interface within layers at which the FSS
          sheet is located. `s=1` implies the sheet is between `layers[1]` and
          `layers[2]`, etc.
- `ψ₁`,`ψ₂`:  Variables containing the unit cell incremental phase shifts (radians).
- `metal`:  A variable which characterizes the metalization region of the FSS/PSS.
- `rwgdat`:  A variable which defines the basis functions.

## Return value

- `zmat`: Complex array of size `(Nbf,Nbf)`, where `Nbf` is the
          number of basis functions. On exit, this array will have been
          filled with the generalized impedance matrix of the moment
          method formulation.

"""
function fillz(k0,u,layers::AbstractVector{Layer},s,ψ₁,ψ₂,metal::RWGSheet,rwgdat::RWGData)

    closed = true              # Always use singularity extraction.
    spinlock = Threads.SpinLock()
    nbf = size(rwgdat.bfe, 2)
    zmat = rwgdat.zorymat
    zmat .= zero(eltype(zmat))

    # Initialize the Green's functions expansion coefficients:
    c3 = c3_calc(k0, u, layers[s].μᵣ, layers[s].ϵᵣ, layers[s+1].μᵣ, layers[s+1].ϵᵣ)
    d3 = d3_calc(k0, u, layers[s].μᵣ, layers[s].ϵᵣ, layers[s+1].μᵣ, layers[s+1].ϵᵣ)

    # Calculate beta00, the fundamental transverse wave vector (1/meter)
    units_per_meter = ustrip(Float64, metal.units, 1u"m")
    β₁, β₂ = metal.β₁ * units_per_meter, metal.β₂ * units_per_meter
    β₀₀ = (ψ₁ * β₁ + ψ₂ * β₂) / twopi

    # Initialize functions for modal series:
    (Σm1_func, Σm2_func) = electric_modal_sum_funcs(k0, u, ψ₁, ψ₂, layers, s, β₁, β₂, β₀₀, 1e-7)

    # floquet_factor is indexed into using values in rwgdat.eci:
    floquet_factor = OffsetArray(zeros(ComplexF64, 5), 0:4)
    floquet_factor[0] = 1.0       # Edges not on unit cell boundary
    floquet_factor[1] = 1.0       # Edges at ξ=0 boundary.
    floquet_factor[2] = cis(-ψ₁)  # Edges at ξ=1 boundary.
    floquet_factor[3] = 1.0       # Edges at η=0 boundary.
    floquet_factor[4] = cis(-ψ₂)  # Edges at η=1 boundary.

    # Set up aliases into data structures for more convenient reference:
    bfe = rwgdat.bfe
    bff = rwgdat.bff
    ebf = rwgdat.ebf
    eci = rwgdat.eci

    nface = size(metal.fv,2)
    i2s = CartesianIndices((nface,nface))
    nbf = size(rwgdat.bfe,2)
    ϵᵣ₁ = layers[s].ϵᵣ
    μᵣ₁ = layers[s].μᵣ
    ϵᵣ₂ = layers[s+1].ϵᵣ
    μᵣ₂ = layers[s+1].μᵣ
    μ̃ = 2 * μᵣ₁ * μᵣ₂ / (μᵣ₁ + μᵣ₂) # Equation (4-11) (normalized to μ₀)
    ϵ̄ = (ϵᵣ₁ + ϵᵣ₂)/2              # Equation (4-23) (normalized to ϵ₀)
    ω = k0 * c₀                  # Radian frequency (Radians/second)

    # Check whether or not the frequency-independent face/face integrals are up to date:
    if ψ₁ ≠ metal.ψ₁ ||  ψ₂ ≠ metal.ψ₂ || metal.u == 0 ||
        abs((u / units_per_meter - metal.u) / metal.u) > tol
        # Set up the values of the face/face integrals' params stored in metal:
        metal.ψ₁ = ψ₁
        metal.ψ₂ = ψ₂
        metal.u = u / units_per_meter
        # Fill the frequency-independent face/face integrals:
        t_spatial = time()
        filljk!(metal, rwgdat, closed)
        t_spatial = time() - t_spatial
        @logfile "      $(round(t_spatial,digits=tdigits)) seconds for spatial face integrals"
    end

    t1 = time_ns()
    Threads.@threads for iufp in 1:rwgdat.nufp  # Loop over each unique face pair
        zcontrib = MVector{9,ComplexF64}(0im, 0im, 0im, 0im, 0im, 0im, 0im, 0im, 0im)
        mbfsave = MVector{9,Int}(0, 0, 0, 0, 0, 0, 0, 0, 0)
        sbfsave = MVector{9,Int}(0, 0, 0, 0, 0, 0, 0, 0, 0)
        ifmifs = rwgdat.ufp2fp[iufp][1]  # Obtain index into face/face matrix
        rowcol = i2s[ifmifs]
        ifm, ifs = rowcol[1], rowcol[2] # indices of match and source triangles
        is = @view metal.fe[:,ifs]     # Obtain the three edges of the source triangle.
        # Obtain the coordinates (in meters) of the source triangle's vertices:
        rs = vtxcrd(ifs, metal) ./ units_per_meter
        # Calculate the signed area of the source triangle:
        rs32 = rs[3] - rs[2]
        rs12 = rs[1] - rs[2]
        area = 0.5 * (rs32[1] * rs12[2] - rs32[2] * rs12[1])
        area48 = 48 * abs(area)  # Needed for surface loading

        self_tri = (ifm == ifs) # Source and match tri are the same?

        if self_tri
            # Obtain the length of each side of the source=match
            # triangle. To be used later in surface loading calculation.
            ls = (norm(rs32), norm(rs[1] - rs[3]), norm(rs12))
        end

        clsflg = self_tri || closed  # Extract if self or if always.

        i_m = @view metal.fe[:,ifm]   # Obtain the three edges of the match triangle.
        rm = vtxcrd(ifm, metal) ./ units_per_meter # Coordinates (m) of the match tri. vertices.
        rmc = mean(rm)        # Match face centroid (meters).

        # Perform frequency-dependent integrals over source triangle
        (I1, I1_ξ, I1_η, I2) = zint(Σm1_func, Σm2_func, rs, rmc)
        # Recall the frequency-independent spatial integrals:
        J = metal.J[iufp]
        J_ξ = metal.J_ξ[iufp]
        J_η = metal.J_η[iufp]
        K = metal.K[iufp]
        K_ξ = metal.K_ξ[iufp]
        K_η = metal.K_η[iufp]
        rinv = metal.rinv[iufp]
        ρ_r = metal.ρ_r[iufp]
        I1_ζ = I1 - I1_ξ - I1_η
        J_ζ = J - J_ξ - J_η
        K_ζ = K - K_ξ - K_η

        # Compute vector from each vertex to centroid of match triangle
        # (divided by 2) as in Eq. (7-15):
        ρc2 =  0.5 * [rmc - rm[i] for i in 1:3]

        # Loop over the face pairs in this equivalence class:
        for (i0, ifmifs2) in enumerate(rwgdat.ufp2fp[iufp])
            zcontrib .= zero(eltype(zcontrib))
            mbfsave .= one(eltype(mbfsave))
            sbfsave .= one(eltype(mbfsave))
            savecounter = 0

            rowcol = i2s[ifmifs2]
            ifm, ifs = rowcol[1], rowcol[2] # indices of match and source triangles
            is = @view metal.fe[:,ifs]      # Obtain the three edges of source tri.
            i_m = @view metal.fe[:,ifm]      # Obtain the three edges of the match tri.
            # Loop over each edge of the source triangle.
            for isl in 1:3   # isl is local edge index, is[isl] is global edge index.
                sbf = ebf[is[isl]]  # Source basis function for this source edge.
                sbf == 0 && continue  # no basis func. is defined
                if bff[1,sbf] == ifs # "Plus" source triangle:
                    source_flag = floquet_factor[eci[is[isl]]]
                elseif bff[2,sbf] == ifs  # "Minus" source triangle:
                    source_flag = -floquet_factor[eci[is[isl]]]
                else
                    error("Impossible situation in source edge loop")
                end
                A_source_flag = μ₀ / fourpi * source_flag * μ̃
                Φ_source_flag = im / ϵ̄ * (source_flag / (twopi*ω*ϵ₀))

                # Compute singular contribution for this edge (the middle term in
                # square brackets in Equation (7-21) using (B.2):
                Asing = ρ_r  +  u * rinv * (rmc - rs[isl])

                # Compute Eq. (7-26) (but correct sign is carried in source flags):
                I1_i = rs[1] * I1_ξ + rs[2] * I1_η + rs[3] * I1_ζ - rs[isl] * I1
                J_i = rs[1] * J_ξ + rs[2] * J_η + rs[3] * J_ζ - rs[isl] * J
                K_i = rs[1] * K_ξ + rs[2] * K_η + rs[3] * K_ζ - rs[isl] * K

                # Compute Equation (7-21):
                A_i = A_source_flag * (fourpi * I1_i + Asing + u * J_i + c3 / u * K_i)
                # Compute Equation (7-31):
                Φ_i = Φ_source_flag * (fourpi * I2 + u * (rinv + J) + d3 / u * K)

                # Now loop over each edge of the match triangle:
                for iml in 1:3 # iml is local edge index, i_m[iml] is global index.
                    savecounter += 1
                    self_edge = i_m[iml] == is[isl] # Match edge = Source edge?
                    mbf = ebf[i_m[iml]]  # Match basis function for this match edge.
                    mbf == 0 && continue  # no basis function defined for this match tri edge
                    if bff[1,mbf] == ifm   # "Plus" match triangle.
                        match_flag = conj(floquet_factor[eci[i_m[iml]]])
                    elseif bff[2,mbf] == ifm  # "Minus" match tri.
                        match_flag = -conj(floquet_factor[eci[i_m[iml]]])
                    else
                        error("Impossible situation in match edge loop")
                    end

                    # Compute one of the dot products in Eq (7-15) apart from sign:
                    dotprod = ρc2[iml] ⋅ A_i
                    # Add contribution to the impedance matrix
                    mbfsave[savecounter] = mbf
                    sbfsave[savecounter] = sbf
                    zcontrib[savecounter] += match_flag * (im*ω*dotprod - Φ_i)
                    #zmat[mbf,sbf] += match_flag * (im*ω*dotprod - Φ_i)
                    # Add surface loading, if applicable:
                    if self_tri && metal.fr[ifm] ≠ 0
                        if self_edge
                            Zload = metal.fr[ifs] / area48 *
                                (3*(ls[next[isl]]^2 + ls[prev[isl]]^2) - ls[isl]^2) * one(ComplexF64)  # Eq. (7-34)
                        else
                            Zload = metal.fr[ifs] / area48 * source_flag * match_flag *
                                (ls[isl]^2  +  ls[iml]^2  - 3 * ls[third[isl,iml]]^2) # Eq. (7-35)
                        end
                        zcontrib[savecounter] += Zload
                        #zmat[mbf,sbf] += Zload
                    end
                end
            end # loop over isl, source triangle edges
            Threads.lock(spinlock)
                @inbounds for kkk in 1:savecounter
                    zmat[mbfsave[kkk], sbfsave[kkk]] += zcontrib[kkk]
                end
            Threads.unlock(spinlock)
        end # face pairs in this equivalence class
    end # facepair loop
    t2 = time_ns()
    tsec = round((t2-t1)/1e9; digits=tdigits)
    @logfile "      $tsec seconds to fill $(size(zmat,1)) × $(size(zmat,2)) matrix entries"
    return zmat
end

"""
    filly(k0,u,layers::AbstractVector{Layer},s,ψ₁,ψ₂,apert::RWGSheet,rwgdat::RWGData) -> ymat

Fill the generalized impedance matrix for an FSS of electric current type.

## Arguments:

- `k0`: Free-space wavenumber (rad/meter).
- `u`:  Green's function smoothing factor (1/meter).
- `layers`:  An array characterizing the dielectric layers surrounding the FSS sheet.
- `s`:  An integer indexing the interface within layers at which the FSS
          sheet is located. `s=1` implies the sheet is between `layers[1]` and
          `layers[2]`, etc.
- `ψ₁`,`ψ₂`:  Variables containing the unit cell incremental phase shifts (radians).
- `apert`:  A variable which characterizes the aperture region of the FSS/PSS.
- `rwgdat`:  A variable which defines the basis functions.

## Return value

- `ymat`: Complex array of size `(Nbf,Nbf)`, where `Nbf` is the
          number of basis functions. On exit, this array will have been
          filled with the generalized admittance matrix of the moment
          method formulation.

"""
function filly(k0, u, layers::AbstractVector{Layer}, s, ψ₁, ψ₂, apert, rwgdat)

    closed = true              # Always use singularity extraction.
    spinlock = Threads.SpinLock()
    nbf = size(rwgdat.bfe, 2)
    ymat = rwgdat.zorymat
    ymat .= zero(eltype(ymat))

    # Initialize the Green's functions expansion coefficients:
    c3s = c3_calc(k0, u, layers[s].μᵣ, layers[s].ϵᵣ, layers[s].μᵣ, layers[s].ϵᵣ)
    c3sp1 = c3_calc(k0, u, layers[s+1].μᵣ, layers[s+1].ϵᵣ, layers[s+1].μᵣ, layers[s+1].ϵᵣ)
    d3s = d3_calc(k0, u, layers[s].μᵣ, layers[s].ϵᵣ, layers[s].μᵣ, layers[s].ϵᵣ)
    d3sp1 = d3_calc(k0, u, layers[s+1].μᵣ, layers[s+1].ϵᵣ, layers[s+1].μᵣ, layers[s+1].ϵᵣ)

    # Calculate beta00, the fundamental transverse wave vector (1/meter)
    units_per_meter = ustrip(Float64, apert.units, 1u"m")
    β₁, β₂ = apert.β₁ * units_per_meter, apert.β₂ * units_per_meter
    β₀₀ = (ψ₁ * β₁ + ψ₂ * β₂) / twopi

    # Initialize functions for modal series:
    (Σm1_func, Σm2_func) = magnetic_modal_sum_funcs(k0, u, ψ₁, ψ₂, layers, s, β₁, β₂, β₀₀,1e-7)

    # floquet_factor is indexed into using values in rwgdat.eci:
    floquet_factor = OffsetArray(zeros(ComplexF64, 5), 0:4)
    floquet_factor[0] = 1.0       # Edges not on unit cell boundary
    floquet_factor[1] = 1.0       # Edges at ξ=0 boundary.
    floquet_factor[2] = cis(-ψ₁)  # Edges at ξ=1 boundary.
    floquet_factor[3] = 1.0       # Edges at η=0 boundary.
    floquet_factor[4] = cis(-ψ₂)  # Edges at η=1 boundary.

    # Set up aliases into data structures for more convenient reference:
    bfe = rwgdat.bfe
    bff = rwgdat.bff
    ebf = rwgdat.ebf
    eci = rwgdat.eci

    nface = size(apert.fv,2)
    i2s = CartesianIndices((nface,nface))
    nbf = size(rwgdat.bfe,2)
    ϵᵣ₁ = layers[s].ϵᵣ
    μᵣ₁ = layers[s].μᵣ
    ϵᵣ₂ = layers[s+1].ϵᵣ
    μᵣ₂ = layers[s+1].μᵣ
    μ̃ = 2.0 * μᵣ₁ * μᵣ₂ / (μᵣ₁ + μᵣ₂) # Eq (4-11) (normalized to μ₀)
    ϵ̄ = 0.5 * (ϵᵣ₁ + ϵᵣ₂) # Equation (1-23) (normalized to ϵ₀)
    ω = k0 * c₀              # Radian frequency (Radians/second)
    I1fact = π / ϵ̄
    KFfact = (c3s*ϵᵣ₁ + c3sp1*ϵᵣ₂) / (2 * ϵ̄ * u)
    I2fact = π * μ̃
    KPfact = μ̃ / (2*u) * (d3s / μᵣ₁ + d3sp1 / μᵣ₂)

    # Check whether or not the apert face/face integrals are up to date:
    if ψ₁ ≠ apert.ψ₁ || ψ₂ ≠ apert.ψ₂ || apert.u == 0 ||
        abs((u / units_per_meter - apert.u) / apert.u) > tol
        # Set up the values of the face/face integrals' params stored in apert:
        apert.ψ₁ = ψ₁
        apert.ψ₂ = ψ₂
        apert.u = u / units_per_meter # Units are 1/(local length units)
        # Fill the frequency-independent face/face integrals:
        t_spatial = time()
        filljk!(apert, rwgdat, closed)
        t_spatial = time() - t_spatial
        @logfile "      $(round(t_spatial,digits=tdigits)) seconds for spatial face integrals"
    end

    t1 = time_ns()
    Threads.@threads for iufp in 1:rwgdat.nufp  # Loop over each unique face pair
        ycontrib = MVector{9,ComplexF64}(0im, 0im, 0im, 0im, 0im, 0im, 0im, 0im, 0im)
        mbfsave = MVector{9,Int}(0, 0, 0, 0, 0, 0, 0, 0, 0)
        sbfsave = MVector{9,Int}(0, 0, 0, 0, 0, 0, 0, 0, 0)
        ifmifs = rwgdat.ufp2fp[iufp][1]  # Obtain index into face/face matrix.
        rowcol = i2s[ifmifs]
        ifm, ifs = rowcol[1], rowcol[2] # indices of match and source triangles
        is = @view apert.fe[:,ifs]     # Obtain the three edges of the source triangle.
        # Obtain the coordinates (in meters) of the source triangle's vertices:
        rs = vtxcrd(ifs, apert) ./ units_per_meter
        # Calculate the signed area of the source triangle:
        rs32 = rs[3] - rs[2]
        rs12 = rs[1] - rs[2]
        area = 0.5 * (rs32[1] * rs12[2] - rs32[2] * rs12[1])
        area48 = 48 * abs(area)  # Needed for surface loading

        self_tri = (ifm == ifs) # Source and match tri are the same?
        clsflg = self_tri || closed  # Extract if self or if always.

        i_m = @view apert.fe[:,ifm]   # Obtain the three edges of the match triangle.
        rm = vtxcrd(ifm, apert) ./ units_per_meter  # Coords (m) of the match tri. vertices.
        rmc = mean(rm)        # Match face centroid (meters).

        # Perform frequency-dependent integrals over source triangle
        (I1, I1_ξ, I1_η, I2) = zint(Σm1_func, Σm2_func, rs, rmc)
        # Recall the frequency-independent spatial integrals:
        J = apert.J[iufp]
        J_ξ = apert.J_ξ[iufp]
        J_η = apert.J_η[iufp]
        K = apert.K[iufp]
        K_ξ = apert.K_ξ[iufp]
        K_η = apert.K_η[iufp]
        rinv = apert.rinv[iufp]
        ρ_r = apert.ρ_r[iufp]
        I1_ζ = I1 - I1_ξ - I1_η
        J_ζ = J - J_ξ - J_η
        K_ζ = K - K_ξ - K_η

        # Compute vector from each vertex to centroid of match triangle
        # (divided by 2) as in Eq. (7.53):
        ρc2 =  0.5 * [rmc - rm[i] for i in 1:3]

        # Loop over the face pairs this equivalence class:
        for (i0, ifmifs) in enumerate(rwgdat.ufp2fp[iufp])
            ycontrib .= zero(eltype(ycontrib))
            mbfsave .= one(eltype(mbfsave))
            sbfsave .= one(eltype(sbfsave))
            savecounter = 0

            rowcol = i2s[ifmifs]
            ifm, ifs = rowcol[1], rowcol[2] # indices of match and source triangles
            is = @view apert.fe[:,ifs]      # Obtain the three edges of source tri.
            i_m = @view apert.fe[:,ifm]      # Obtain the three edges of the match tri.
            # Loop over each edge of the source triangle.
            for isl in 1:3 # isl is local edge index, is[isl] is global edge index.
                sbf = ebf[is[isl]]  # Source basis function for this edge.
                sbf == 0 && continue  # no basis function is defined
                if bff[1,sbf] == ifs # "Plus" source triangle:
                    source_flag = floquet_factor[eci[is[isl]]]
                elseif bff[2,sbf] == ifs # "Minus" source triangle:
                    source_flag = -floquet_factor[eci[is[isl]]]
                else
                    error("Impossible situation in source edge loop")
                end
                F_source_flag = -ϵ₀ / π * ϵ̄ * source_flag
                Ψ_source_flag = 2im / (π*ω*μ₀*μ̃) * source_flag

                # Compute singular contribution for this edge (the middle term in
                # square brackets in Equation (7-57)) using (B-2):
                Fsing = ρ_r  +  u * rinv * (rmc - rs[isl])

                # Compute Eq. (7-26) (but correct sign is carried in source flags):
                I1_i = rs[1] * I1_ξ + rs[2] * I1_η + rs[3] * I1_ζ - rs[isl] * I1
                J_i = rs[1] * J_ξ + rs[2] * J_η + rs[3] * J_ζ - rs[isl] * J
                K_i = rs[1] * K_ξ + rs[2] * K_η + rs[3] * K_ζ - rs[isl] * K

                # Compute Equation (7-57):
                F_i = F_source_flag * (I1fact * I1_i + Fsing + u * J_i + KFfact * K_i)

                # Compute Equation (7-60):
                Ψ_i = Ψ_source_flag * (I2fact * I2 + u * (rinv  + J) + KPfact * K)

                # Now loop over each edge of the match triangle:
                for iml in 1:3 # iml is local edge index, im[iml] is global index.
                    savecounter += 1
                    self_edge = i_m[iml] == is[isl] # Match edge = Source edge?
                    mbf = ebf[i_m[iml]]  # Match basis function for this edge.
                    mbf == 0 && continue # no basis function defined for this edge
                    if bff[1,mbf] == ifm  # "Plus" match triangle.
                        match_flag = conj(floquet_factor[eci[i_m[iml]]])
                    elseif bff[2,mbf] == ifm # "Minus" match tri.
                        match_flag = -conj(floquet_factor[eci[i_m[iml]]])
                    else
                        error("Impossible situation in match edge loop")
                    end

                    # Compute one of the dot products in Eq (7-53) apart from sign:
                    dotprod = ρc2[iml] ⋅ F_i

                    # Add contribution to the admittance matrix
                    mbfsave[savecounter] = mbf
                    sbfsave[savecounter] = sbf
                    ycontrib[savecounter] += match_flag * (-im*ω*dotprod - Ψ_i)
                    #ymat[mbf,sbf] += match_flag * (-im*ω*dotprod - Ψ_i)
                end
            end # loop over source edges
            Threads.lock(spinlock)
                @inbounds for kkk in 1:savecounter
                    ymat[mbfsave[kkk], sbfsave[kkk]] += ycontrib[kkk]
                end
            Threads.unlock(spinlock)
        end
    end # loop over face pairs
    t2 = time_ns()
    tsec = round((t2-t1)/1e9; digits=tdigits)
    @logfile "      $tsec seconds to fill $(size(ymat,1)) × $(size(ymat,2)) matrix entries"

    return ymat
end

end # module