text:
"""
Various linear and nonlinear eigenvalue problems.

  * LEP: *L*inear *E*igenvalue *P*roblem, i.e. no bulk dispersion, fixed boundary conditions (can use with PMLs/cPMLs for resonance or RSM calculations)
  * CF: *C*onstant *F*lux states, i.e. at a given frequency, find eigen-susceptibilities and associated eigenfunctions. Can use either matched boundaries or PMLs. The CF problem is a particular kind of linear eigenvalue problem, though here it is treated separately.
  * NEP: *N*onlinear *E*igenvalue *P*roblem, i.e. with bulk dispersion, open boundary conditions, or both. Used to compute resonances/RSMs, no need for PMLs.

API: [`MaxwellLEP`](@ref), [`MaxwellCF`](@ref), [`MaxwellNEP`](@ref), for use with [`maxwelleigen`](@ref) (high-level),
or directly with [`ArnoldiMethod`](https://github.com/haampie/ArnoldiMethod.jl),
[`Arpack`](https://github.com/JuliaLinearAlgebra/Arpack.jl),
[`IterativeSolvers`](https://github.com/JuliaMath/IterativeSolvers.jl),
[`KrylovKit`](https://github.com/Jutho/KrylovKit.jl),
[`LinearAlgebra`](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/),
[`NonlinearEigenproblems`](https://github.com/nep-pack/NonlinearEigenproblems.jl) (low-level).
"""
module Spectral

export HelmholtzLEP
export HelmholtzCF
export HelmholtzNEP
export MaxwellLEP
export MaxwellCF
export MaxwellNEP
export helmholtzeigen
export maxwelleigen

using ..Common
using NonlinearEigenproblems

dimensions = (
		"1D/Spectral1D.jl",
        # "2D/Spectral2D.jl",
        # "3D/Spectral3D.jl"
        )

interfaces = (
	# "Interfaces/ArnoldiMethod.jl",
	"Interfaces/Arpack.jl",
	# "Interfaces/KrylovKit.jl",
	"Interfaces/LinearAlgebra.jl",
	"Interfaces/NonlinearEigenproblems.jl"
	)

# defin abstract containers

"""
	AbstractEigenproblem{N}
"""
abstract type AbstractEigenproblem{N} end

"""
	AbstractLinearEigenproblem{N} <: AbstractEigenproblem{N}
"""
abstract type AbstractLinearEigenproblem{N} <: AbstractEigenproblem{N} end

"""
	AbstractCFEigenproblem{N} <: AbstractEigenproblem{N}
"""
abstract type AbstractCFEigenproblem{N} <: AbstractEigenproblem{N} end

"""
	AbstractNonlinearEigenproblem{N} <: AbstractEigenproblem{N}
"""
abstract type AbstractNonlinearEigenproblem{N} <: AbstractEigenproblem{N} end

helmholtzeigen() = nothing
maxwelleigen() = nothing

################################################################################
# Define Eigenproblem Structures for Maxwell & Helmholtz

using SparseArrays

LEPs = (:Helmholtz,:Maxwell)
leps = (:helmholtz,:maxwell)
rep = (1,3) # num of times to repeat F for helmholtz v. maxwell
for (L,l,r) ∈ zip(LEPs,leps,rep)
	@eval begin
		"""
			$($(string(L,"LEP"))){N}

		`N`-dimensional $($(string(L))) Linear Eigenvalue Problem
		"""
		struct $(Symbol(L,"LEP")){N,T} <: AbstractLinearEigenproblem{N}
		    operator::T
			αεpFχ::SparseMatrixCSC{ComplexF64,Int}
			saturated::Base.RefValue{Bool}

		    $(Symbol(L,"LEP"))($l::$L{N}) where N = new{N,typeof($l)}($l,$l.αεpFχ,Ref(false))
		end

		"""
			$($(string(L,"CF"))){N}

		`N`-dimensional $($(string(L))) Constant Flux problem
		"""
		struct $(Symbol(L,"CF")){N,T} <: AbstractCFEigenproblem{N}
		    operator::T
			F::SparseMatrixCSC{Float64,Int}
			saturated::Base.RefValue{Bool}
		    $(Symbol(L,"CF"))($l::$L{N}) where N = new{N,typeof($l)}($l,spdiagm(0=>repeat($l.simulation.F,$r)),Ref(false))
		end

		"""
			$($(string(L,"NEP"))){N}

		`N`-dimensional $($(string(L))) Constant Flux problem
		"""
		struct $(Symbol(L,"NEP")){N,T,TNEP} <: AbstractNonlinearEigenproblem{N}
		    operator::T
		    nep::TNEP
			Fs::Vector{SparseMatrixCSC{ComplexF64,Int}}
			saturated::Base.RefValue{Bool}

		    $(Symbol(L,"NEP"))($l::$L{N},NEP,F_inds) where N = new{N,typeof($l),typeof(NEP)}($l,NEP,NEP.A[F_inds],Ref(false))
		end
	end
end

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

# include dimensions
foreach(include,dimensions)

################################################################################
# add convenience constructors for VectorFields

# ScalarFields for Helmholtz
for ep ∈ (HelmholtzLEP,HelmholtzCF,HelmholtzNEP)
	@eval begin
		Common.VectorField(hep::$(ep),n::Integer=1) = ScalarField(hep,n)
		Common.VectorField(arg,hep::$(ep)) = ScalarField(hep,arg)
		Common.VectorField(hep::$(ep),arg) = ScalarField(hep,arg)
		Common.ScalarField(hep::$(ep),n::Integer=1) = ScalarField(hep.simulation.x,n)
		Common.ScalarField(arg,hep::$(ep)) = ScalarField(hep.simulation.x,arg)
		Common.ScalarField(hep::$(ep),arg) = ScalarField(hep.simulation.x,arg)
	end
end

# ElectricFields for Maxwell
for ep ∈ (MaxwellLEP,MaxwellCF,MaxwellNEP)
# 	@eval begin
# 		# """
# 		# 	ElectricField(::$($ep),[m=1]) -> ElectricField
# 		# 	ElectricField(::$($ep),values) -> ElectricField
# 		#
# 		# arguments can be provided in any order
# 		# """
# 		Common.ElectricField(mep::$(ep)) = ElectricField(mep.simulation)
# 		Common.ElectricField(arg,mep::$(ep)) = ElectricField(mep.simulation,arg)
# 		Common.ElectricField(mep::$(ep),arg) = ElectricField(mep.simulation,arg)
# 	end
end

################################################################################
# LEP

"""
    HelmholtzLEP(sim) -> lep

`sim` is a `Simulation`, `lep` can be passed to `helmholtzeigen`
"""
HelmholtzLEP(args...; kwargs...) = HelmholtzLEP(Helmholtz(args...; kwargs...))

"""
    MaxwellLEP(sim; ky=0, kz=0) -> lep

`sim` is a `Simulation`, `lep` can be passed to `maxwelleigen`
"""
MaxwellLEP(args...; kwargs...) = MaxwellLEP(Maxwell(args...; kwargs...))

# generate Ax=λBx matrices, with or without saturation for Helmholtz
"""
    lep(ω) -> ∇², -ε-χ(ω)

	lep(ω,ωs,ψs) -> ∇², -ε-χ(ω,ωs,ψs)

`ωs` is an array of frequencies, `ψs` is a `ScalarField` containing one field for
each element of `ωs`
"""
function (lep::HelmholtzLEP)(ω, args...)
    lep.helmholtz(ω, args...)
	lep.saturated[] = isempty(args) ? false : true
    return lep.helmholtz.D², -lep.helmholtz.αεpFχ
end

# generate Ax=λBx matrices, with or without saturation for Maxwell
"""
    lep(ω) -> ∇×∇×, ε+χ(ω)

	lep(ω,ωs,ψs) -> ∇×∇×, ε+χ(ω,ωs,ψs)

`ωs` is an array of frequencies, `ψs` is an `ElectricField` containing one field for
each element of `ωs`
"""
function (lep::MaxwellLEP)(ω, args...)
    lep.maxwell(ω, args...)
	lep.saturated[] = isempty(args) ? false : true
    return lep.maxwell.D², lep.maxwell.αεpFχ
end

################################################################################
# CF

"""
    HelmholtzCF(::Simulation) -> lep
"""
HelmholtzCF(args...; kwargs...) = HelmholtzCF(Helmholtz(args...; kwargs...))
"""
    MaxwellCF(::Simulation; ky=0, kz=0) -> lep
"""
MaxwellCF(args...; kwargs...) = MaxwellCF(Maxwell(args...; kwargs...))

"""
	cf(ω) -> ∇² + εω² , -Fω²

	cf(ω,ωs,ψs) -> ∇² + εω² , -Fω²/(1 + Σᵢ|γ(ωsᵢ)ψs(i)|²)

`ωs` is an array of frequencies, `ψs` is a ScalarField containing one field for
each element of `ωs`
"""
function (cf::HelmholtzCF)(ω,args...)
	cf.helmholtz(ω,args...) # initialize
	_computeF!(cf.F,cf.simulation,ω,args...)
	cf.saturated[] = isempty(args) ? false : true
	ω² = ω^2
    rows = rowvals(cf.helmholtz.A)
    vals = nonzeros(cf.helmholtz.A)
    _, n = size(cf.helmholtz.A)
    for i ∈ 1:n
        col = i
        for j ∈ nzrange(cf.helmholtz.A, i)
            row = rows[j]
            vals[j] = cf.helmholtz.D²[row,col] + ω²*cf.helmholtz.αε[row,col]
        end
    end
	return cf.helmholtz.A, -ω²*cf.F
end

"""
	cf(ω) -> ∇×∇× - εω² , Fω²

	cf(ω,ωs,ψs) -> ∇×∇× - εω² , Fω²/(1 + Σᵢ|γ(ωsᵢ)ψs(i)|²)

`ωs` is an array of frequencies, `ψs` is an ElectricField containing one field for
each element of `ωs`
"""
function (cf::MaxwellCF)(ω,args...)
	cf.maxwell(ω,args...) # initialize
	_computeF!(cf.F,cf.simulation,ω,args...)
	cf.saturated[] = isempty(args) ? false : true
	ω² = ω^2
    rows = rowvals(cf.maxwell.A)
    vals = nonzeros(cf.maxwell.A)
    _, n = size(cf.maxwell.A)
    for i ∈ 1:n
        col = i
        for j ∈ nzrange(cf.maxwell.A, i)
            row = rows[j]
            vals[j] = cf.maxwell.D²[row,col] - ω²*cf.maxwell.αε[row,col]
        end
    end
	return cf.maxwell.A, ω²*cf.F
end

################################################################################
# NEP
# see dimensional files for dimension-specific constructors

"""
	nep(ω) -> A

	nep(ω,ωs,ψs) -> A

Full helmholtz operator `A` at frequency `ω` for nonlinear eigenproblem `nep`.
Problem saturated by `ωs::Vector` and `ψs::ScalarField`, if provided.
"""
function (mnep::HelmholtzNEP)(ω,args...)
	# saturate if extra args are provided
	if !isempty(args)
		mnep.helmholtz(ω,args...)
		mnep.saturated[] = true
	else
		mnep.saturated[] = false
	end
	# compute F, including saturation
	_computeF!(mnep.Fs, mnep.simulation, ω, args...)
	return compute_Mder(mnep.nep,ω)::SparseMatrixCSC{ComplexF64,Int}
end

"""
	nep(ω) -> A

	nep(ω,ωs,ψs) -> A

Full maxwell operator `A` at frequency `ω` for nonlinear eigenproblem `nep`.
Problem saturated by `ωs::Vector` and `ψs::ScalarField`, if provided.
"""
function (mnep::MaxwellNEP)(ω,args...)
	# saturate if extra args are provided
	if !isempty(args)
		mnep.maxwell(ω,args...) #
		mnep.saturated[] = true
	else
		mnep.saturated[] = false
	end
	# compute F, including saturation
	_computeF!(mnep.Fs, mnep.simulation, ω, args...)
	return compute_Mder(mnep.nep,ω)::SparseMatrixCSC{ComplexF64,Int}
end

################################################################################
# getproperty and propertynames for LEP, CF, NEP

for (L,l) ∈ zip(LEPs,leps)
	for ep ∈ (:LEP,:CF,:NEP)
		@eval begin
			function Base.getproperty(lep::$(Symbol(L,ep)), sym::Symbol)
				if Base.sym_in(sym, propertynames(getfield(lep,:operator)))
			        return getproperty(getfield(lep,:operator),sym)
				elseif Base.sym_in(sym, (:helmholtz, :maxwell))
					return getfield(lep,:operator)
				else
					return getfield(lep,sym)
			    end
			end

			Base.propertynames(lep::$(Symbol(L,ep)), private=false) = private ? fieldnames($(Symbol(L,ep))) : propertynames(lep.operator, false)
		end
	end
end

# Base.propertynames(cf::HelmholtzCF, private=false) = private ? fieldnames(HelmholtzCF) : propertynames(cf.helmholtz)

################################################################################
# saturate F

# CF: No saturation, reset to sim.F
function _computeF!(F::SparseMatrixCSC,sim::Simulation,ω)
	for i ∈ 1:size(F,1)
		j = mod1(i,length(sim))
		# either field is zero if no dispersive domain, or RESET to sim.F value
		if iszero(sim.dispersive_domain_indices[j])
			F.nzval[i] = 0
		else
			F.nzval[i] = sim.F[j]
		end
	end
	return nothing
end

# NEP: No saturation, reset to sim.Fs
function _computeF!(F::Vector{T},sim::Simulation,ω) where T<:AbstractSparseMatrix
	for d ∈ eachindex(F)
		for i ∈ 1:size(F[1],1)
			j = mod1(i,length(sim))
			if sim.dispersive_domain_indices[j] == d
				F[d].nzval[i] = sim.Fs[d][j]
			else
				F[d].nzval[i] = 0
			end
		end
	end
	return nothing
end

# CF: Saturate F (in cf.F, not in sim.F) according to background fields
function _computeF!(F::SparseMatrixCSC, sim::Simulation, ω, ωs::Array, ψs::VectorField)
	for i ∈ 1:size(F,1)
		j = mod1(i,length(sim))
		if typeof(sim.χ[j])<:TwoLevelSystem
			χ::TwoLevelSystem = sim.χ[j]
			if iszero(χ.D0)
				F.nzval[i] = 0
			else
				if iszero(sim.dispersive_domain_indices[j])
					F.nzval[i] = 0
				else
					F.nzval[i] = sim.F[j]*χ.hp1⁻¹[i]
				end
			end
		end
	end
	return nothing
end

# NEP: Saturate F (not in sim.F) according to background fields
function _computeF!(F::Vector{T}, sim::Simulation, ω, ωs::Array, ψs::VectorField) where T<:AbstractSparseMatrix
	for d ∈ eachindex(F)
		for i ∈ 1:size(F[1],1)
			j = mod1(i,length(sim))
			if typeof(sim.χ[j])<:TwoLevelSystem
				χ::TwoLevelSystem = sim.χ[j]
				if iszero(χ.D0)
					F[d].nzval[i] = 0
				else
					if sim.dispersive_domain_indices[j] == d
						F[d].nzval[i] = sim.F[j]*χ.hp1⁻¹[i]
					else
						F[d].nzval[i] = 0
					end
				end
			end
		end
	end
	return nothing
end

################################################################################
# ORTHOGONALIZE

using LinearAlgebra
import ..Common.ORTHOGONALIZE_OVERLAP_THRESHOLD

"""
	orthogonalize!(ψ,sim,η,B,ky,kz,[ind])

Check for nearly degenerate eigenpairs, if overlap integral exceeds $ORTHOGONALIZE_OVERLAP_THRESHOLD,
orthoganlize and renormalize.
"""
function orthogonalize!(u::ScalarField{N},sim::Simulation{N},η,B,args...) where N
	for i ∈ eachindex(η)
		for j ∈ 1:(i-1)
			if abs(2(η[i]-η[j])/(η[i]+η[j])) < 1e-3
				δ = sum(u[:,i].*u[:,j].*diag(B))*_measure(sim)
				if abs(δ) > ORTHOGONALIZE_OVERLAP_THRESHOLD
					β = 1/sqrt(1-δ^2)
					α = -δ*β
					u[:,j] = α*u[:,i] + β*u[:,j]
				end
			end
		end
	end
	return nothing
end

function orthogonalize!(u::ElectricField{N},sim::Simulation{N},η,B,ky,kz,args...) where N
	for i ∈ eachindex(η)
		for j ∈ 1:(i-1)
			if abs(2(η[i]-η[j])/(η[i]+η[j])) < 1e-3
				δ = sum(u[:,i].*u[:,j].*diag(B))*_measure(sim)
				if abs(δ) > ORTHOGONALIZE_OVERLAP_THRESHOLD
					if all(iszero,(ky,kz))
						u.z[:,j] .= 0
						u.y[:,i] .= 0
					else
						β = 1/sqrt(1-δ^2)
						α = -δ*β
						u[:,j] = α*u[:,i] + β*u[:,j]
					end
				end
			end
		end
	end
	if all(iszero,(ky,kz)) normalize!(sim,u,B,args...) end
	return nothing
end

################################################################################
# Pretty Printing

import ..Common.PRINTED_COLOR_DARK
import ..Common.PRINTED_COLOR_VARIABLE
import ..Common.PRINTED_COLOR_INSTRUCTION
import ..Common.PRINTED_COLOR_WARN

# LEP
function Base.show(io::IO,lep::AbstractLinearEigenproblem{N}) where N
    print(io,N,"D ")
	lep.saturated[] ? printstyled(io,"saturated ",color=PRINTED_COLOR_WARN) : nothing#print(io,"unsaturated ")
	if typeof(lep) <: MaxwellLEP
	    printstyled(io,"MaxwellLEP ",color=PRINTED_COLOR_DARK)
	elseif typeof(lep) <: HelmholtzLEP
	    printstyled(io,"HelmholtzLEP ",color=PRINTED_COLOR_DARK)
	end
    printstyled(io,"(pass to ",color=PRINTED_COLOR_INSTRUCTION)
	if typeof(lep) <: MaxwellLEP
	    printstyled(io,"maxwelleigen",color=PRINTED_COLOR_VARIABLE)
	elseif typeof(lep) <: HelmholtzLEP
	    printstyled(io,"helmholtzeigen",color=PRINTED_COLOR_VARIABLE)
	end
    printstyled(io,", or call w/ args ",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"ω",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,",[",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"ωs",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,",",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"ψs",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,"]->",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"A",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,",",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"B",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,", where ",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"Aψ=ω²Bψ",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,")",color=PRINTED_COLOR_INSTRUCTION)
end

# CF
function Base.show(io::IO,cf::AbstractCFEigenproblem{N}) where N
    print(io,N,"D ")
	cf.saturated[] ? printstyled(io,"saturated ",color=PRINTED_COLOR_WARN) : nothing#print(io,"unsaturated ")
	if typeof(cf) <: MaxwellCF
		printstyled(io,"MaxwellCF ",color=PRINTED_COLOR_DARK)
	elseif typeof(cf) <: HelmholtzCF
		printstyled(io,"HelmholtzCF ",color=PRINTED_COLOR_DARK)
	end
    printstyled(io,"(pass to ",color=PRINTED_COLOR_INSTRUCTION)
	if typeof(cf) <: MaxwellCF
		printstyled(io,"maxwelleigen",color=PRINTED_COLOR_VARIABLE)
	elseif typeof(cf) <: HelmholtzCF
		printstyled(io,"helmholtzeigen",color=PRINTED_COLOR_VARIABLE)
	end
    printstyled(io,", or call w/ args ",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"ω",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,",[",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"ωs",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,",",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"ψs",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,"]->",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"A",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,",",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"B",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,", where ",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"Au=ηBu",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,")",color=PRINTED_COLOR_INSTRUCTION)
end

# NEP
function Base.show(io::IO,nep::AbstractNonlinearEigenproblem{N}) where N
    print(io,N,"D ")
	nep.saturated[] ? printstyled(io,"saturated ",color=PRINTED_COLOR_WARN) : nothing#print(io,"unsaturated ")
	if typeof(nep) <: MaxwellNEP
	    printstyled(io,"MaxwellNEP ",color=PRINTED_COLOR_DARK)
	elseif typeof(nep) <: HelmholtzNEP
		printstyled(io,"HelmholtzNEP ",color=PRINTED_COLOR_DARK)
	end
    printstyled(io,"(pass to ", color=PRINTED_COLOR_INSTRUCTION)
	if typeof(nep) <: MaxwellNEP
	    printstyled(io,"maxwelleigen",color=PRINTED_COLOR_VARIABLE)
	elseif typeof(nep) <: HelmholtzNEP
		printstyled(io,"helmholtzeigen",color=PRINTED_COLOR_VARIABLE)
	end
	printstyled(io,", or NEP solver like ",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"iar_chebyshev",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,", or call w/ arg ",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"ω",color=PRINTED_COLOR_VARIABLE)
	printstyled(io,",[",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"ωs",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,",",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"ψs",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,"]->",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"A",color=PRINTED_COLOR_VARIABLE)
	printstyled(io,", where ",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"Aψ=0",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,")",color=PRINTED_COLOR_INSTRUCTION)
end

################################################################################
# DOCSTRINGS
# to be paired with default interface for `helmholtzeigen` and `maxwelleigen`

DOC_LEP_H = "
	helmholtzeigen(lep, Ω, [ωs, ψs; kwargs...]) -> ω, ψ

Linear eigenvalues `ω` closest to `Ω`, eigenfields `ψ::ScalarField`.

Optional `ωs::Vector` and `ψs::ScalarField` must be provided together, and define
nonlinearity-inducing background fields.

Keywords:"

DOC_LEP_M = "
	maxwelleigen(lep, Ω, [ωs, ψs; kwargs...]) -> ω, ψ

Linear eigenvalues `ω` closest to `Ω`, eigenfields `ψ::ElectricField`.

Optional `ωs::Vector` and `ψs::ElectricField` must be provided together, and define
nonlinearity-inducing background fields.

Keywords:"

DOC_CF_H = "
	helmholtzeigen(cf, Ω, [ωs, ψs; η=0, kwargs...]) -> H, u

Linear eigen-susceptibilities `H` closest to `η`, eigenfields `u::ScalarField`.
Optional `ωs::Vector` and `ψs::ScalarField` must be provided together, and define
nonlinearity-inducing background fields.

Keywords:
`η` Target CF eigenvalue (`0`);"

DOC_CF_M = "
	maxwelleigen(cf, Ω, [ωs, ψs; η=0, kwargs...]) -> H, u

Linear eigen-susceptibilities `H` closest to `η`, eigenfields `u::ElectricField`.
Optional `ωs::Vector` and `ψs::ElectricField` must be provided together, and define
nonlinearity-inducing background fields.

Keywords:
`η` Target CF eigenvalue (`0`);"

DOC_NEP_H = "
	helmholtzeigen(nep, Ω, [ωs, ψs; kwargs...]) -> ω, ψ

Nonlinear eigenvalues `ω` closest to `Ω`, eigenfields `ψ::ScalarField`.
Optional `ωs::Vector` and `ψs::ScalarField` must be provided together, and define
nonlinearity-inducing background fields.

Keywords:"

DOC_NEP_M = "
	maxwelleigen(nep, Ω, [ωs, ψs; kwargs...]) -> ω, ψ

Nonlinear eigenvalues `ω` closest to `Ω`, eigenfields `ψ::ElectricField`.
Optional `ωs::Vector` and `ψs::ElectricField` must be provided together, and define
nonlinearity-inducing background fields.

Keywords:"

################################################################################
# Interfaces

import ..Common.DEFAULT_LINEAR_EIGENSOLVER
import ..Common.DEFAULT_NONLINEAR_EIGENSOLVER
import ..Common.INDEX_OFFSET

foreach(include,interfaces)

end # module


text:
using Base: @_noinline_meta, @_inline_meta
using Core: is_top_bit_set
using Core.Intrinsics: bitcast, trunc_int, sext_int, zext_int, sle_int, eq_int, and_int

# onehot erorr
struct OneHotEncodeError <: Exception
    K
    val
    OneHotEncodeError(@nospecialize(K), @nospecialize(val)) = (@_noinline_meta; new(K, val))
end

function Base.showerror(io::IO, e::OneHotEncodeError)
    print(io, "OneHotEncodeError: cannot encode ")
    print(io, e.val)
    print(io, " with OneHot{")
    print(io, e.K)
    print(io, '}')
end

throw_onehotencode_error(K, val) = (@_noinline_meta; throw(OneHotEncodeError(K, val)))

# onehot encode
primitive type OneHot{K} <: AbstractOneHotArray{1} 32 end

OneHot(k) = OneHot{UInt32(k)}
OneHot{K}(x) where K = convert(OneHot(K), x)
OneHot(k, x) = OneHot{k}(x)

onehotsize(::OneHot{K}) where K = Int(K)

# array interface

Base.size(o::OneHot) = (onehotsize(o),)
Base.@propagate_inbounds function Base.getindex(o::OneHot, i::I) where {I<:Integer}
    @boundscheck checkbounds(o, i)
    return convert(I, o) == i
end

Base.@propagate_inbounds Base.getindex(o::OneHot, i::Colon) = o

Base.argmax(o::OneHot) = Int(o)
Base.argmin(o::OneHot) = isone(o) ? 2 : 1

# printing

function Base.showarg(io::IO, x::OneHot, toplevel)
    toplevel || print(io, "::")
    print(io, "OneHot{", onehotsize(x), '}')
end

# convert

Base.UInt32(o::OneHot) = bitcast(UInt32, o)
Base.UInt64(o::OneHot) = zext_int(UInt64, o)
Base.Int32(o::OneHot) = bitcast(Int32, o)
Base.Int64(o::OneHot) = zext_int(Int64, o)

Base.convert(::Type{Any}, o::OneHot) = o
Base.convert(::Type{OneHot{K}}, o::OneHot{K}) where {K} = o
Base.convert(::Type{OneHot{K}}, o::OneHot) where {K} = OneHot(K, UInt32(o))
Base.convert(::Type{UInt32},  o::OneHot) = UInt32(o)
Base.convert(::Type{To}, o::OneHot) where {To<:Number} = convert(To, UInt32(o))
Base.convert(::Type{T}, o::OneHot) where {T<:Array} = T(o)

Base.reinterpret(::Type{T}, o::OneHot) where {T} = reinterpret(T, UInt32(o))

# one

Base.one(o::O) where {O<:OneHot} = convert(O, 0x00000001)
Base.one(::Type{<:OneHot{K}}) where K = OneHot(K, 1)
Base.isone(o::OneHot) = isone(convert(UInt32, o))

# zero

Base.zero(o::O) where {O<:OneHot} = convert(O, 0x00000000)
Base.zero(::Type{<:OneHot{K}}) where K = OneHot(K, 0)
Base.iszero(o::OneHot) = iszero(convert(UInt32, o))

# number

Base.typemin(::Type{OneHot{K}}) where K = OneHot(K, 0)
Base.typemax(::Type{OneHot{K}}) where K = OneHot(K, K)

Base.isless(a::OneHot, b::OneHot) = isless(UInt32(a), UInt32(b))
Base.:(==)(a::OneHot, b::OneHot) = UInt32(a) == UInt32(b)

# bit-op

function check_onehot_top_bit(::Type{OneHot{K}}, x) where {K}
    @_inline_meta
    is_top_bit_set(x) && throw_onehotencode_error(K, x)
    x
end

function check_onehot_encode(ot::Type{OneHot{K}}, x) where {K}
    @_inline_meta
    sle_int(x, K) || throw_onehotencode_error(K, x)
    bitcast(ot, x)
end

function checked_onehot_trunc_sint(ot::Type{OneHot{K}}, x::From) where {K, From}
    @_inline_meta
    y = trunc_int(UInt32, x)
    back = sext_int(From, y)
    eq_int(x, back) || throw_onehotencode_error(K, x)
    check_onehot_encode(ot, y)
end

function checked_onehot_trunc_uint(ot::Type{OneHot{K}}, x::From) where {K, From}
    @_inline_meta
    y = trunc_int(UInt32, x)
    back = zext_int(From, y)
    eq_int(x, back) || throw_onehotencode_error(K, x)
    check_onehot_encode(ot, y)
end

Base.convert(ot::Type{OneHot{K}}, x::Int8) where {K} = check_onehot_encode(ot, sext_int(UInt32, check_onehot_top_bit(ot, x)))
Base.convert(ot::Type{OneHot{K}}, x::Int16) where {K} = check_onehot_encode(ot, sext_int(UInt32, check_onehot_top_bit(ot, x)))
Base.convert(ot::Type{OneHot{K}}, x::Int32) where {K} = check_onehot_encode(ot, bitcast(UInt32, check_onehot_top_bit(ot, x)))
Base.convert(ot::Type{OneHot{K}}, x::Int64) where {K} = checked_onehot_trunc_sint(ot, check_onehot_top_bit(ot, x))
Base.convert(ot::Type{OneHot{K}}, x::Int128) where {K} = checked_onehot_trunc_sint(ot, check_onehot_top_bit(ot, x))
Base.convert(ot::Type{OneHot{K}}, x::UInt8) where {K} = check_onehot_encode(ot, zext_int(UInt32, x))
Base.convert(ot::Type{OneHot{K}}, x::UInt16) where {K} = check_onehot_encode(ot, zext_int(UInt32, x))
Base.convert(ot::Type{OneHot{K}}, x::UInt32) where {K} = check_onehot_encode(ot, x)
Base.convert(ot::Type{OneHot{K}}, x::UInt64) where {K} = checked_onehot_trunc_uint(ot, x)
Base.convert(ot::Type{OneHot{K}}, x::UInt128) where {K} = checked_onehot_trunc_uint(ot, x)
Base.convert(ot::Type{OneHot{K}}, x::Bool) where {K} = and_int(zext_int(ot, x), Base.convert(ot, 0x1))

ext: jl


text:
# * Derivatives
# ** Three-point stencils

# α = super-/subdiagonal of second derivative
# α̃ = subdiagonal of second derivative, in case not same as superdiagonal
# β = -1/2 diagonal of second derivative
# γ = diagonal of first derivative
# δ = super-/subdiagonal of first derivative

# *** Cartesian finite-differences

α(::FiniteDifferences{T}, ::Integer) where T = one(T)
β(::FiniteDifferences{T}, ::Integer) where T = one(T)
γ(::FiniteDifferences{T}, ::Integer) where T = zero(T)
δ(::FiniteDifferences{T}, ::Integer) where T = one(T)

# *** Staggered finite-differences

# Eq. (20), Schafer (2000)
α(::Uniform,  ::StaggeredFiniteDifferences{T}, j::Integer) where T = j^2/(j^2 - one(T)/4)
β(::Uniform, B::StaggeredFiniteDifferences{T}, j::Integer) where T = (j^2 - j + one(T)/2)/(j^2 - j + one(T)/4)

α(d::Uniform, B::StaggeredFiniteDifferences) = α.(Ref(d), Ref(B), 1:length(B.r)-1)
β(d::Uniform, B::StaggeredFiniteDifferences) = β.(Ref(d), Ref(B), 1:length(B.r))
δ(d::Uniform, B::StaggeredFiniteDifferences) = α(d, B)

function staggered_stencil_common(fun, B::StaggeredFiniteDifferences{T}, δN=0) where T
    r = B.r
    N = length(r)
    v = zeros(T, N+δN)

    r̃ = vcat(r, 2r[end]-r[end-1])
    a = 2r[1]-r[2]
    b,c,d = r[1],r[2],r[3]

    for j = 1:N+δN
        if j < N
            d = r̃[j+2]
        end
        v[j] = fun(a,b,c,d)

        a,b,c = b,c,d
    end
    v
end

function α(::NonUniform, B::StaggeredFiniteDifferences)
    staggered_stencil_common(B, -1) do a,b,c,d
        # Eq. (A13) Krause 1999
        2/((c-b)*√((d-b)*(c-a)))*((b+c)/2)^2/(b*c)
    end
end

function α(::NonUniform, B::StaggeredFiniteDifferences, f::Function)
    staggered_stencil_common(B, -1) do a,b,c,d
        rⱼ₊₁₂ = (b+c)/2
        # Eq. (A13) Krause 1999
        2f(rⱼ₊₁₂)/((c-b)*√((d-b)*(c-a)))*rⱼ₊₁₂^2/(b*c)
    end
end

function β(::NonUniform, B::StaggeredFiniteDifferences)
    staggered_stencil_common(B) do a,b,c,d
        # Eq. (A14) Krause 1999
        fp = ((c+b)/2b)^2
        fm = ((b+a)/2b)^2
        1/(c-a)*(1/(c-b)*fp + 1/(b-a)*fm)
    end
end

function β(::NonUniform, B::StaggeredFiniteDifferences, f::Function)
    staggered_stencil_common(B) do a,b,c,d
        # Eq. (A14) Krause 1999
        b⁻² = inv(b^2)
        rⱼ₊₁₂ = (b+c)/2
        rⱼ₋₁₂ = (a+b)/2
        fp = rⱼ₊₁₂^2 * b⁻²
        fm = rⱼ₋₁₂^2 * b⁻²
        1/(c-a)*(f(rⱼ₊₁₂)/(c-b)*fp + f(rⱼ₋₁₂)/(b-a)*fm)
    end
end

function δ(::NonUniform, B::StaggeredFiniteDifferences)
    staggered_stencil_common(B, -1) do a,b,c,d
        2/(√((d-b)*(c-a)))*((b+c)/2)^2/(b*c)
    end
end

α(B::StaggeredFiniteDifferences, args...)   = α(distribution(B), B, args...)
β(B::StaggeredFiniteDifferences, args...)   = β(distribution(B), B, args...)
γ(B::StaggeredFiniteDifferences{T}) where T = zeros(T, length(B.r))
δ(B::StaggeredFiniteDifferences)            = δ(distribution(B), B)

# *** Implicit finite-differences

α( ::Uniform, B::ImplicitFiniteDifferences{T}) where T = ones(T, length(B.r)-1)
α̃(d::Uniform, B::ImplicitFiniteDifferences) = α(d, B)
β( ::Uniform, B::ImplicitFiniteDifferences{T}) where T = ones(T, length(B.r))
γ( ::Uniform, B::ImplicitFiniteDifferences{T}) where T = zeros(T, length(B.r))
δ( ::Uniform, B::ImplicitFiniteDifferences{T}) where T = ones(T, length(B.r)-1)
δ̃(d::Uniform, B::ImplicitFiniteDifferences) = δ(d, B)

# Non-uniform, implicit finite differences, as described by
# Patchkovskii (2016)

function implicit_stencil_common(fun, B::ImplicitFiniteDifferences{T}, s=1, δN=0) where T
    r = B.r
    N = length(r)
    v = zeros(T, N+δN)

    r̃ = vcat(r[s:end], 2r[end]-r[end-1])
    a = r̃[1]-(s > 1 ? r[s-1] : 0)
    b = r̃[2]-r̃[1]
    c = zero(T)
    for j = 1:N+δN
        if j < N - (s-1)
            c = r̃[j+2]-r̃[j+1]
        end

        v[j] = fun(a,b,c)

        a,b = b,c
    end
    v
end

function α̃(::NonUniform, B::ImplicitFiniteDifferences)
    implicit_stencil_common(B, 2, -1) do a,b,c
        2/(a*(a+b)) # Eq. (33)
    end
end

function α(::NonUniform, B::ImplicitFiniteDifferences)
    implicit_stencil_common(B, 1, -1) do a,b,c
        2/(b*(a+b)) # Eq. (33)
    end
end

function β(::NonUniform, B::ImplicitFiniteDifferences)
    implicit_stencil_common(B) do a,b,c
        1/(a*b) # Eq. (33), remember β is multiplied by -2
    end
end

α(B::ImplicitFiniteDifferences) = α(distribution(B), B)
α̃(B::ImplicitFiniteDifferences) = α̃(distribution(B), B)
β(B::ImplicitFiniteDifferences) = β(distribution(B), B)
γ(B::ImplicitFiniteDifferences) = γ(distribution(B), B)
δ(B::ImplicitFiniteDifferences) = δ(distribution(B), B)
δ̃(B::ImplicitFiniteDifferences) = δ̃(distribution(B), B)

# *** Common

α(B::AbstractFiniteDifferences) = α.(Ref(B), 1:size(B,2)-1)
α̃(B::AbstractFiniteDifferences) = α(B)
β(B::AbstractFiniteDifferences) = β.(Ref(B), 1:size(B,2))
γ(B::AbstractFiniteDifferences) = γ.(Ref(B), 1:size(B,2))
δ(B::AbstractFiniteDifferences) = δ.(Ref(B), 1:size(B,2)-1)
δ̃(B::AbstractFiniteDifferences) = δ(B)

# *** Restrictions

for f in [:α, :α̃, :δ, :δ̃]
    @eval begin
        function $f(B̃::RestrictedFiniteDifferences)
            j = indices(B̃,2)
            $f(parent(B̃))[j[1]:(j[end]-1)]
        end
    end
end

for f in [:β, :γ]
    @eval begin
        function $f(B̃::RestrictedFiniteDifferences)
            $f(parent(B̃))[indices(B̃,2)]
        end
    end
end

for f in [:α, :α̃, :β, :γ, :δ, :δ̃]
    @eval begin
        function $f(Ã::BasisOrRestricted{<:AbstractFiniteDifferences},
                    B̃::BasisOrRestricted{<:AbstractFiniteDifferences},
                    band=0)
            i,i′ = extrema(indices(Ã,2))
            j,j′ = extrema(indices(B̃,2))
            k = max(i+min(0,band),j-max(0,band)):min(i′+min(0,band),j′-max(0,band))
            $f(parent(B̃))[k]
        end
    end
end

# ** Materialization

non_symmetric_laplacian(::Uniform, B::ImplicitFiniteDifferences) = false
non_symmetric_laplacian(::NonUniform, B::ImplicitFiniteDifferences) = true
non_symmetric_laplacian(B::ImplicitFiniteDifferences) = non_symmetric_laplacian(distribution(B), B)
non_symmetric_laplacian(::AbstractFiniteDifferences) = false

function derivative_matrix(::Type{T}, Ac, B, diff_order) where T
    A = parent(Ac)
    parent(A) == parent(B) ||
        throw(ArgumentError("Cannot multiply functions on different grids"))
    Ai,Bi = last(indices(A)), last(indices(B))
    if Ai == Bi
        n = length(Ai)
        if diff_order == 1 || non_symmetric_laplacian(parent(B))
            Tridiagonal(Vector{T}(undef, n-1),
                        Vector{T}(undef, n),
                        Vector{T}(undef, n-1))
        elseif diff_order == 2
            SymTridiagonal(Vector{T}(undef, n),
                           Vector{T}(undef, n-1))
        end
    else
        m,n = length(Ai),length(Bi)
        offset = Ai[1]-Bi[1]
        BandedMatrix{T}(undef, (m,n), (1-offset,1+offset))
    end
end

@materialize function *(Ac::AdjointBasisOrRestricted{<:AbstractFiniteDifferences},
                        D::Derivative,
                        B::BasisOrRestricted{<:AbstractFiniteDifferences})
    T -> begin
        derivative_matrix(T, Ac, B, 1)
    end
    dest::Tridiagonal{T} -> begin
        A = parent(Ac)
        parent(A) == parent(B) ||
            throw(ArgumentError("Cannot multiply functions on different grids"))

        # Central difference approximation
        μ = 2step(B)
        dest.dl .= -δ̃(B)/μ
        dest.d .= γ(B)/μ
        dest.du .= δ(B)/μ
    end
    dest::BandedMatrix{T} -> begin
        A = parent(Ac)
        parent(A) == parent(B) ||
            throw(ArgumentError("Cannot multiply functions on different grids"))

        dl,d,du = bandrange(dest)
        # Central difference approximation
        μ = 2step(B)
        dest[Band(dl)] .= -δ̃(A,B,-1)/μ
        dest[Band(d)] .= γ(A,B)/μ
        dest[Band(du)] .= δ(A,B,1)/μ
    end
end

@materialize function *(Ac::AdjointBasisOrRestricted{<:AbstractFiniteDifferences},
                        Dc::QuasiAdjoint{<:Any,<:Derivative},
                        D::Derivative,
                        B::BasisOrRestricted{<:AbstractFiniteDifferences})
    T -> begin
        derivative_matrix(T, Ac, B, 2)
    end
    dest::SymTridiagonal{T} -> begin
        A = parent(Ac)
        parent(A) == parent(B) ||
            throw(ArgumentError("Cannot multiply functions on different grids"))

        δ² = step(B)^2
        dest.dv .= -2β(B)/δ²
        dest.ev .= α(B)/δ²
    end
    dest::Tridiagonal{T} -> begin
        A = parent(Ac)
        parent(A) == parent(B) ||
            throw(ArgumentError("Cannot multiply functions on different grids"))

        δ² = step(B)^2
        dest.dl .= α̃(B)/δ²
        dest.d .= -2β(B)/δ²
        dest.du .= α(B)/δ²
    end
    dest::BandedMatrix{T} -> begin
        A = parent(Ac)
        parent(A) == parent(B) ||
            throw(ArgumentError("Cannot multiply functions on different grids"))

        dl,d,du = bandrange(dest)
        δ² = step(B)^2
        dest[Band(dl)] .= α̃(A,B,-1)/δ²
        dest[Band(d)] .= -2β(A,B)/δ²
        dest[Band(du)] .= α(A,B,1)/δ²
    end
end

@materialize function *(Ac::AdjointBasisOrRestricted{<:StaggeredFiniteDifferences},
                        Dc::QuasiAdjoint{<:Any,<:Derivative},
                        O::QuasiDiagonal,
                        D::Derivative,
                        B::BasisOrRestricted{<:StaggeredFiniteDifferences})
    T -> begin
        derivative_matrix(T, Ac, B, 2)
    end
    dest::SymTridiagonal{T} -> begin
        A = parent(Ac)
        parent(A) == parent(B) ||
            throw(ArgumentError("Cannot multiply functions on different grids"))

        δ² = step(B)^2
        f = r -> O.diag[max(zero(T),r)]
        dest.dv .= -2β(B,f)/δ²
        dest.ev .= α(B,f)/δ²
    end
    dest::Tridiagonal{T} -> begin
        A = parent(Ac)
        parent(A) == parent(B) ||
            throw(ArgumentError("Cannot multiply functions on different grids"))

        δ² = step(B)^2
        dest.dl .= α̃(B)/δ²
        dest.d .= -2β(B, f)/δ²
        dest.du .= α(B)/δ²
    end
    # dest::BandedMatrix{T} -> begin
    #     A = parent(Ac)
    #     parent(A) == parent(B) ||
    #         throw(ArgumentError("Cannot multiply functions on different grids"))

    #     dl,d,du = bandrange(dest)
    #     δ² = step(B)^2
    #     dest[Band(dl)] .= α̃(A,B,-1)/δ²
    #     # Need to include f here
    #     dest[Band(d)] .= -2β(A,B)/δ²
    #     dest[Band(du)] .= α(A,B,1)/δ²
    # end
end

# *** Implicit derivatives

#=
This is an implementation of finite difference scheme described in

- Muller, H. G. (1999). An Efficient Propagation Scheme for the
  Time-Dependent Schrödinger equation in the Velocity Gauge. Laser
  Physics, 9(1), 138–148.

where the first derivative is approximated as

\[\partial f =
\left(1+\frac{h^2}{6}\Delta_2\right)^{-1}
\Delta_1 f
\equiv
M_1^{-1}\tilde{\Delta}_1 f,\]
where
\[M_1 \equiv
\frac{1}{6}
\begin{bmatrix}
4+\lambda' & 1 &\\
1 & 4 & 1\\
& 1 & 4 & \ddots\\
&&\ddots&\ddots\\
\end{bmatrix},\]
and
\[\tilde{\Delta}_1 \equiv
\frac{1}{2h}
\begin{bmatrix}
\lambda & 1 &\\
-1 &  & 1\\
& -1 &  & \ddots\\
&&\ddots&\ddots\\
\end{bmatrix},\]

where \(\lambda=\lambda'=\sqrt{3}-2\) for problems with a singularity
at the boundary \(r=0\) and zero otherwise; and the second derivative
as

\[\partial^2 f =
\left(1+\frac{h^2}{12}\Delta_2\right)^{-1}
\Delta_2 f
\equiv
-2M_2^{-1}\Delta_2 f,\]
where
\[M_2 \equiv
-\frac{1}{6}
\begin{bmatrix}
10-2\delta\beta_1 & 1 &\\
1 & 10 & 1\\
& 1 & 10 & \ddots\\
&&\ddots&\ddots\\
\end{bmatrix},\]
and
\[\Delta_2 \equiv
\frac{1}{h^2}
\begin{bmatrix}
-2(1+\delta\beta_1) & 1 &\\
1 & -2 & 1\\
& 1 & -2 & \ddots\\
&&\ddots&\ddots\\
\end{bmatrix},\]

where, again, \(\delta\beta_1 = -Zh[12-10Zh]^{-1}\) is a correction
introduced for problems singular at the origin.

=#

function implicit_lhs(::Uniform, B::BasisOrRestricted{<:ImplicitFiniteDifferences{T}}, difforder) where T
    m = size(B,2)
    M = SymTridiagonal(Vector{T}(undef, m), Vector{T}(undef, m-1))

    j₁ = first(indices(B, 2))

    if difforder == 1
        M.dv .= 4
        M.ev .= 1

        # f'  =   M₁⁻¹*Δ₁*f   Eq. (19)
        M /= 6
    elseif difforder == 2
        M.dv .= 10
        M.ev .= 1

        # f'' = -2M₂⁻¹*Δ₂*f   Eq. (13)
        M /= (-6)*(-2)
    end
end

function implicit_lhs(::NonUniform, B::BasisOrRestricted{<:ImplicitFiniteDifferences{T}}, difforder) where T
    dl, d, du = if difforder == 1
    elseif difforder == 2
        # Eq. (33), Patchkovskii (2016)
        dl = implicit_stencil_common(B, 2, -1) do a,b,c
            (a^2 + a*b - b^2)/(6*a*(a+b))
        end
        d = implicit_stencil_common(B) do a,b,c
            (a^2 + 3a*b + b^2)/(6*a*b)
        end
        du = implicit_stencil_common(B, 1, -1) do a,b,c
            (-a^2 + a*b + b^2)/(6*b*(a+b))
        end
        dl, d, du
    end

    Tridiagonal(dl, d, du)
end

implicit_lhs(B::BasisOrRestricted{<:ImplicitFiniteDifferences}, args...) =
    implicit_lhs(distribution(B), B, args...)

function implicit_derivative(::Type{T}, M, difforder) where T
    A = parent(first(M.args))
    B = last(M.args)
    A == B || throw(DimensionMismatch("Implicit finite-differences does not support different restrictions"))
    ImplicitDerivative(copyto!(similar(M, T), M),
                       implicit_lhs(last(M.args), difforder))
end

LazyArrays.simplifiable(::typeof(*), ::AdjointBasisOrRestricted{B}, ::Derivative, ::BasisOrRestricted{B}) where {T,B<:ImplicitFiniteDifferences{T}} = Val(true)
LazyArrays._simplify(::typeof(*), A::AdjointBasisOrRestricted{B}, D::Derivative, C::BasisOrRestricted{B})  where {T,B<:ImplicitFiniteDifferences{T}} =
        implicit_derivative(T, ApplyQuasiArray(*,A,D,C), 1)

LazyArrays.simplifiable(::typeof(*), ::AdjointBasisOrRestricted{B}, ::QuasiAdjoint{T,<:Derivative}, ::Derivative, ::BasisOrRestricted{B}) where {T,B<:ImplicitFiniteDifferences{T}} = Val(true)
LazyArrays._simplify(::typeof(*), Ac::AdjointBasisOrRestricted{B}, Dc::QuasiAdjoint{T,<:Derivative}, D::Derivative, A::BasisOrRestricted{B}) where {T,B<:ImplicitFiniteDifferences{T}} =
        implicit_derivative(T, ApplyQuasiArray(*,Ac,Dc,D,A), 2)

# *** Coulomb derivatives

# Below, we implement some corrections to the derivative matrices in
# case of an attractive Coulomb problem [-∂²/2 - Z/r + ℓ(ℓ+1)/2r^2].
#
# For StaggeredFiniteDifferences, the correction is just a heuristic
# quoted by
#
# - Schafer, K. J. (2009). Numerical Methods in Strong Field Physics. In
#   T. Brabec (Eds.), (pp. 111–145). Springer.
#
# whereas for ImplicitFiniteDifferences, it is derived via Taylor
# expansions, for the uniform case, see
#
# - Muller, H. G. (1999). An Efficient Propagation Scheme for the
#   Time-Dependent Schrödinger equation in the Velocity Gauge. Laser
#     Physics, 9(1), 138–148.
#
# and for the non-uniform case, see
#
# - Patchkovskii, S., & Muller, H. (2016). Simple, Accurate, and
#   Efficient Implementation of 1-Electron Atomic Time-Dependent
#   Schrödinger Equation in Spherical Coordinates. Computer Physics
#   Communications, 199,
#   153–169. http://dx.doi.org/10.1016/j.cpc.2015.10.014

gradient_correction!(∂, B::AbstractFiniteDifferences, Z, ℓ) = ∂

function gradient_correction!(∂::ImplicitDerivative, ::Uniform, B::ImplicitFiniteDifferences, Z, ℓ)
    # This is a correction for the gradient on the half-line to
    # approximately restore anti-Hermiticity, cf. Eq. (20ff), Muller
    # (1999).

    ρ = step(B)

    ∂.Δ[1,1] += (√3 - 2)/2ρ
    ∂.M[1,1] += (√3 - 2)/6

    ∂.M⁻¹ = factorize(∂.M)
    ∂
end

function gradient_correction!(∂::ImplicitDerivative, ::NonUniform, B::ImplicitFiniteDifferences, Z, ℓ)
    Z == 0 && return ∂
    a = B.r[1]
    b = B.r[2] - a
    if ℓ == 0
        # Eq. (31), Patchkovskii (2016)
        num = (a^2 + b*a + b^2)*(2*(a+b)*(a+2b)*a*Z - 6a^2 - 15b*a - 8b^2)
        ∂.Δ[1,1] = -(a+b)^2*(a+2b)*(2a^3*Z + b*a*(3 - 2b*Z) - 6a^2 + b^2)/(a*b*num)
        ∂.M[1,1] = (a+b)^2*(a+2b)*((a+b)*a*Z - 3a - b)/num
        ∂.M[1,2] = a^2*(a+2b)*((a+b)*a*Z - 3a - 2b)/num
    end
    ∂.M⁻¹ = factorize(∂.M)
    ∂
end

gradient_correction!(∂, B::ImplicitFiniteDifferences, Z, ℓ) =
    gradient_correction!(∂, distribution(B), B::ImplicitFiniteDifferences, Z, ℓ)

@materialize function *(Ac::AdjointBasisOrRestricted{<:AbstractFiniteDifferences},
                        CD::CoulombDerivative,
                        B::BasisOrRestricted{<:AbstractFiniteDifferences})
    T -> begin
        derivative_matrix(T, Ac, B, 1)
    end
    dest::AbstractMatrix{T} -> begin
        A = parent(Ac)
        parent(A) == parent(B) ||
            throw(ArgumentError("Cannot multiply functions on different grids"))

        D = CD.D

        copyto!(dest, ApplyQuasiArray(*,Ac,D,A))
        if first(indices(parent(Ac),2)) == first(indices(B,2)) == 1
            gradient_correction!(dest, parent(B), CD.Z, CD.ℓ)
        end
        dest
    end
end

function laplacian_correction!(∂²::AbstractMatrix, B::StaggeredFiniteDifferences, Z, ℓ)
    if ℓ == 0 && Z ≠ 0
        ρ = local_step(B,1)
        # Eq. (22) Schafer (2009)
        δβ₁ = 1/ρ^2*Z*ρ/8*(1 + Z*ρ)
        ∂²[1,1] -= 2δβ₁
    end

    ∂²
end

function laplacian_correction!(∂²::ImplicitDerivative, ::Uniform, B::ImplicitFiniteDifferences, Z, ℓ)
    if ℓ == 0 && Z ≠ 0
        ρ = step(B)
        # Eq. (17), Muller (1999)
        δβ₁ = -Z*ρ/(12 - 10Z*ρ)
        ∂².Δ[1,1] -= 2δβ₁/ρ^2
        ∂².M[1,1] -= 2δβ₁/((-6)*(-2))

        ∂².M⁻¹ = factorize(∂².M)
    end
    ∂²
end

function laplacian_correction!(∂²::ImplicitDerivative, ::NonUniform, B::ImplicitFiniteDifferences, Z, ℓ)
    Z == 0 && return ∂²
    a = B.r[1]
    b = B.r[2] - a
    if ℓ == 0
        # Eq. (34), Patchkovskii (2016)
        num = b*(a*(3a+4b)*Z - 6*(a+b))
        ∂².Δ[1,1] = 2*(a+b)*(6a - (3a - b)*(a+b)*Z)/(a^2*num)
        ∂².M[1,1] = (a+b)*(a*(a+b)*(a+3b)*Z - 3*(a^2 + 3b*a + b^2))/(3a*num)
        ∂².M[1,2] = (a^3*(-Z) + a^2*(b*Z+3) + b*a*(2*b*Z-3) - 3*b^2)/(3*num)
    else
        a = B.r[1]
        b = B.r[2] - a

        # Eq. (35)
        num = b*(3a+4b)
        ∂².Δ[1,1] = (3a - b)*(a+b)^2/(a^3*numb)
        ∂².M[1,1] = (a+b)^2*(a+3b)/(3a*num)
        ∂².M[1,2] = -(a-2b)*(a+b)/(3num)
    end
    ∂².M⁻¹ = factorize(∂².M)
    ∂²
end

laplacian_correction!(∂²::ImplicitDerivative, B::ImplicitFiniteDifferences, Z, ℓ) =
    laplacian_correction!(∂², distribution(B), B, Z, ℓ)

laplacian_correction!(∂², B::AbstractFiniteDifferences, Z, ℓ) = ∂²

@materialize function *(Ac::AdjointBasisOrRestricted{<:AbstractFiniteDifferences},
                        CDc::QuasiAdjoint{<:Any,<:CoulombDerivative},
                        CD::CoulombDerivative,
                        B::BasisOrRestricted{<:AbstractFiniteDifferences})
    T -> begin
        derivative_matrix(T, Ac, B, 2)
    end
    dest::AbstractMatrix{T} -> begin
        A = parent(Ac)
        parent(A) == parent(B) ||
            throw(ArgumentError("Cannot multiply functions on different grids"))

        @assert parent(CDc).Z == CD.Z
        Dc = adjoint(parent(CDc).D)
        D = CD.D

        copyto!(dest, ApplyQuasiArray(*,Ac,Dc,D,A))
        if first(indices(parent(Ac),2)) == first(indices(B,2)) == 1
            laplacian_correction!(dest, parent(B), CD.Z, CD.ℓ)
        end
        dest
    end
end

@materialize function *(Ac::AdjointBasisOrRestricted{<:AbstractFiniteDifferences},
                        CDc::QuasiAdjoint{<:Any,<:CoulombDerivative},
                        O::QuasiDiagonal,
                        CD::CoulombDerivative,
                        B::BasisOrRestricted{<:AbstractFiniteDifferences})
    T -> begin
        derivative_matrix(T, Ac, B, 2)
    end
    dest::AbstractMatrix{T} -> begin
        A = parent(Ac)
        parent(A) == parent(B) ||
            throw(ArgumentError("Cannot multiply functions on different grids"))

        @assert parent(CDc).Z == CD.Z
        Dc = adjoint(parent(CDc).D)
        D = CD.D

        copyto!(dest, ApplyQuasiArray(*,Ac,Dc,O,D,A))
        if first(indices(parent(Ac),2)) == first(indices(B,2)) == 1
            laplacian_correction!(dest, parent(B), CD.Z, CD.ℓ)
        end
        dest
    end
end

LazyArrays.simplifiable(::typeof(*), ::AdjointBasisOrRestricted{R}, ::CoulombDerivative, ::BasisOrRestricted{R}) where {T,R<:ImplicitFiniteDifferences{T}} = Val(true)
function LazyArrays._simplify(::typeof(*), Ac::AdjointBasisOrRestricted{R}, CD::CoulombDerivative, B::BasisOrRestricted{R}) where {T,R<:ImplicitFiniteDifferences{T}}
    D = CD.D
    ∂ = implicit_derivative(T, ApplyQuasiArray(*,Ac,D,B), 1)
    if first(indices(parent(Ac),2)) == first(indices(B,2)) == 1
        gradient_correction!(∂, parent(B), CD.Z, CD.ℓ)
    end
    ∂
end

LazyArrays.simplifiable(::typeof(*), ::AdjointBasisOrRestricted{R}, ::QuasiAdjoint{T,<:CoulombDerivative}, ::CoulombDerivative, ::BasisOrRestricted{R}) where {T,R<:ImplicitFiniteDifferences{T}} = Val(true)
function LazyArrays._simplify(::typeof(*), Ac::AdjointBasisOrRestricted{R}, CDc::QuasiAdjoint{T,<:CoulombDerivative}, CD::CoulombDerivative, B::BasisOrRestricted{R}) where {T,R<:ImplicitFiniteDifferences{T}}
    @assert parent(CDc).Z == CD.Z
    Dc = adjoint(parent(CDc).D)
    D = CD.D

    ∂² = implicit_derivative(T, ApplyQuasiArray(*,Ac,Dc,D,B), 2)
    if first(indices(parent(Ac),2)) == first(indices(B,2)) == 1
        laplacian_correction!(∂², parent(B), CD.Z, CD.ℓ)
    end
    ∂²
end


text:
test_delivery_unassigned_env = {
    # By hour starting at 8AM
    "coverage"  => Array[
        [3.0, 3.0, 3.0, 3.0, 3.0, 4.0, 4.0, 6.0, 6.0, 7.0, 7.0, 5.0, 2.0, 1.0],
        [2.0, 2.0, 4.0, 5.0, 3.0, 3.0, 4.0, 4.0, 7.0, 6.0, 5.0, 4.0, 2.0, 0.0],
        [3.0, 3.0, 5.0, 4.0, 5.0, 4.0, 5.0, 5.0, 7.0, 6.0, 6.0, 5.0, 3.0, 1.0],
        [2.0, 3.0, 3.0, 5.0, 5.0, 2.0, 5.0, 6.0, 7.0, 6.0, 6.0, 5.0, 3.0, 1.0],
        [3.0, 4.0, 4.0, 3.0, 5.0, 4.0, 6.0, 5.0, 6.0, 8.0, 5.0, 4.0, 4.0, 3.0],
        [3.0, 4.0, 5.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 6.0, 6.0, 5.0, 2.0, 2.0],
        [3.0, 4.0, 5.0, 5.0, 5.0, 4.0, 6.0, 6.0, 11.0, 8.0, 6.0, 4.0, 2.0, 0.0]
    ],
    "cycle_length" => 14,
    "shift_time_min" => 4,
    "shift_time_max" => 8,
    "time_between_coverage" => (10),
}

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
# This file was generated, do not modify it.

using MLJ
using PyPlot
ioff() # hide
import Statistics

Xs = source()
ys = source()

atom = @load DecisionTreeRegressor
atom.n_subfeatures = 4 # to ensure diversity among trained atomic models

machines = (machine(atom, Xs, ys) for i in 1:100)

Statistics.mean(v...) = mean(v)
Statistics.mean(v::AbstractVector{<:AbstractNode}) = node(mean, v...)

yhat = mean([predict(m, Xs) for  m in machines]);

surrogate = Deterministic()
mach = machine(surrogate, Xs, ys; predict=yhat)

@from_network mach begin
    mutable struct OneHundredModels
        atom=atom
    end
end

one_hundred_models = OneHundredModels()

X, y = @load_boston;

r = range(atom,
          :min_samples_split,
          lower=2,
          upper=100, scale=:log)

mach = machine(atom, X, y)

curve = learning_curve!(mach,
                        range=r,
                        measure=mav,
                        resampling=CV(nfolds=9),
                        verbosity=0)

plot(curve.parameter_values, curve.measurements)
xlabel(curve.parameter_name)

savefig(joinpath(@OUTPUT, "e1.svg")) # hide

r = range(one_hundred_models,
          :(atom.min_samples_split),
          lower=2,
          upper=100, scale=:log)

mach = machine(one_hundred_models, X, y)

curve = learning_curve!(mach,
                        range=r,
                        measure=mav,
                        resampling=CV(nfolds=9),
                        verbosity=0)

plot(curve.parameter_values, curve.measurements)
xlabel(curve.parameter_name)

savefig(joinpath(@OUTPUT, "e2.svg")) # hide

PyPlot.close_figs() # hide


text:
using DataFrames, CSV, Statistics, StatsBase
data = CSV.read("../../../externals/Core.Math.Samples/data/Pejcic_318.csv", copycols = true)

# describe(data[:ATT])

for col in eachcol(data,true);
    if (eltype(col[2]) <: Real) println(describe(col[2]));
    end;
end;


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

using Markdown
using InteractiveUtils

# ╔═╡ ed7c391e-d937-11ea-31e0-d5884018ca29
md"# Orthogonal Procrustes"

# ╔═╡ fd1fbc70-d936-11ea-2712-ed19c93c0c9f
begin
	using DelimitedFiles
	using Gadfly, Cairo, Colors
	using LinearAlgebra
	using Statistics

	set_default_plot_size(15cm,13cm)
	color_palette = Scale.color_discrete().f(7)[[1, 4, 6]]
	themes = [Theme(default_color=color,
			        panel_fill="white",
			        grid_color="lightgray",
                    background_color=parse(Colorant, "white")
			        ) for color in color_palette];
end

# ╔═╡ 3b0cfa3e-d937-11ea-2d53-255142465b91
B = readdlm("point_cloud.txt")

# ╔═╡ 3d667abc-d937-11ea-338d-e194f0379319
p1 = plot(x=B[1, :], y=B[2, :], Geom.point, themes[1],
         Guide.title("a point set, B"), themes[1],
         Coord.cartesian(fixed=true),
        )

# ╔═╡ 53c23896-d937-11ea-2e06-6b6331e0fe4c
draw(PNG("ptset_B.png", 5inch, 4inch, dpi=300), p1)

# ╔═╡ 7ca103a0-d937-11ea-0d5a-af0e3a9141a6
rotation_matrix2d(θ::Float64) = [cos(θ) -sin(θ); sin(θ) cos(θ)]

# ╔═╡ 83746456-d937-11ea-04a1-f15bf81a17d4
begin
	θ = π * 0.8
	R_known = rotation_matrix2d(θ)
end

# ╔═╡ 8b9acb16-d937-11ea-0e02-7f6bcec57932
begin
	ϵ = 0.01 # noise
	A = R_known * B .+ ϵ * randn(size(B)...)

	p2 = plot(x=A[1, :], y=A[2, :], Geom.point, themes[3],
	         Guide.title("point set, A"), themes[3],
	         Coord.cartesian(fixed=true),
	    )
end

# ╔═╡ a85a6086-d937-11ea-3fed-21f714cad837
draw(PNG("ptset_A.png", 5inch, 4inch, dpi=300), p2)

# ╔═╡ b04a9cfa-d937-11ea-38cd-e1a2a9d588e3
p3 = plot(layer(x=A[1, :], y=A[2, :], Geom.point, themes[3]),
         layer(x=B[1, :], y=B[2, :], Geom.point, themes[1]),
         Guide.title("orthogonal Procrustes"), themes[1],
         Coord.cartesian(fixed=true),
         Guide.manual_color_key("point set", ["A", "B"], [color_palette[3], color_palette[1]])
    )

# ╔═╡ b7f97a90-d937-11ea-1d20-a3752908fe8c
draw(PNG("before_alignment.png", 5inch, 4inch, dpi=300), p3)

# ╔═╡ bcb90346-d937-11ea-1bbb-2f481b298a95
p4 = plot(layer(x=A[1, :], y=A[2, :], Geom.point, themes[3]),
         layer(x=B[1, :], y=B[2, :], Geom.point, themes[1]),
         [layer(x=[A[1, i], B[1, i]], y=[A[2, i], B[2, i]],
            Geom.line, Theme(default_color=colorant"gray")) for i = 1:size(A)[2]]...,
         Guide.title("orthogonal Procrustes"), themes[1],
         Coord.cartesian(fixed=true),
         Guide.manual_color_key("point set", ["A", "B"], [color_palette[3], color_palette[1]])
    )

# ╔═╡ c17420c0-d937-11ea-3fad-0b62648de42f
draw(PNG("correspondence.png", 5inch, 4inch, dpi=300), p4)

# ╔═╡ c61c4986-d937-11ea-2423-51c2154ee671
F = svd(A * B')

# ╔═╡ c9fae814-d937-11ea-3e2a-8d010deb8f54
R = F.V * F.U'

# ╔═╡ ce49e546-d937-11ea-3c90-2150e15a7ce8
isapprox(R, rotation_matrix2d(-θ), rtol=0.01)

# ╔═╡ d4242594-d937-11ea-1a8e-c1b01e9563f8
A_transformed = R * A

# ╔═╡ d86f9a84-d937-11ea-2807-6de9156292c8
p5 = plot(layer(x=A_transformed[1, :], y=A_transformed[2, :], Geom.point, themes[3]),
         layer(x=B[1, :], y=B[2, :], Geom.point, themes[1]),
         Guide.title("orthogonal Procrustes"), themes[1],
         Coord.cartesian(fixed=true),
         Guide.manual_color_key("point set", ["RA", "B"], [color_palette[3], color_palette[1]])
    )

# ╔═╡ e461e8cc-d937-11ea-07f3-a1a8abd49594
draw(PNG("after_alignment.png", 5inch, 4inch, dpi=300), p5)

# ╔═╡ Cell order:
# ╟─ed7c391e-d937-11ea-31e0-d5884018ca29
# ╠═fd1fbc70-d936-11ea-2712-ed19c93c0c9f
# ╠═3b0cfa3e-d937-11ea-2d53-255142465b91
# ╠═3d667abc-d937-11ea-338d-e194f0379319
# ╠═53c23896-d937-11ea-2e06-6b6331e0fe4c
# ╠═7ca103a0-d937-11ea-0d5a-af0e3a9141a6
# ╠═83746456-d937-11ea-04a1-f15bf81a17d4
# ╠═8b9acb16-d937-11ea-0e02-7f6bcec57932
# ╠═a85a6086-d937-11ea-3fed-21f714cad837
# ╠═b04a9cfa-d937-11ea-38cd-e1a2a9d588e3
# ╠═b7f97a90-d937-11ea-1d20-a3752908fe8c
# ╠═bcb90346-d937-11ea-1bbb-2f481b298a95
# ╠═c17420c0-d937-11ea-3fad-0b62648de42f
# ╠═c61c4986-d937-11ea-2423-51c2154ee671
# ╠═c9fae814-d937-11ea-3e2a-8d010deb8f54
# ╠═ce49e546-d937-11ea-3c90-2150e15a7ce8
# ╠═d4242594-d937-11ea-1a8e-c1b01e9563f8
# ╠═d86f9a84-d937-11ea-2807-6de9156292c8
# ╠═e461e8cc-d937-11ea-07f3-a1a8abd49594

max_stars_repo_name: SimonEnsemble/catbook

max_issues_repo_name: SimonEnsemble/catbook

max_forks_repo_name: SimonEnsemble/catbook


text:
using FuzzySets, Test

@test FuzzySets.name() == "FuzzySets.jl"

step_size = 0.005
levels = collect(0:step_size:1)

A = FuzzyNumber(levels, number=0.0, width=0.5)
@test A.levels == levels
@test A.grades[1] == Interval(-0.5, 0.5)
@test A[1] == Interval(-0.5, 0.5)
@test A(-0.2) == 0.6
@test A(0) == 1.0
@test A(0.2) == 0.6
@test A(0.5) == 0.0
@test A(1) == 0.0
@test length(A) == 201

@test support(A) == Interval(-0.5, 0.5)
@test core(A) == Interval(0) == Interval(0.0)
@test height(A) == 1.0
@test cut(A, 0) == Interval(-0.5, 0.5)
@test cut(A, 0.5) == Interval(-0.25, 0.25)
@test cut(A, 1) == Interval(0)


text:
abstract type AbstractAggregation end
abstract type AbstractAggregator end
"""
    RawAggregationParameter

An optimizer attribute used for raw parameters of the aggregator. Defers to `RawParameter`.
"""
struct RawAggregationParameter <: AggregationParameter
    name::Any
end

function MOI.get(aggregator::AbstractAggregator, param::RawAggregationParameter)
    name = Symbol(param.name)
    if !(name in fieldnames(typeof(aggregator)))
        error("Unrecognized parameter name: $(name) for aggregator $(typeof(aggregator)).")
    end
    return getfield(aggregator, name)
end

function MOI.set(aggregator::AbstractAggregator, param::RawAggregationParameter, value)
    name = Symbol(param.name)
    if !(name in fieldnames(typeof(aggregator)))
        error("Unrecognized parameter name: $(name) for aggregator $(typeof(aggregator)).")
    end
    setfield!(aggregator, name, value)
    return nothing
end

# ------------------------------------------------------------
include("cut_collection.jl")
include("distance_measures.jl")
include("selection_rules.jl")
include("cluster_rules.jl")
include("lock_rules.jl")
include("no_aggregation.jl")
include("partial_aggregation.jl")
include("dynamic_aggregation.jl")
include("cluster_aggregation.jl")
include("hybrid_aggregation.jl")


text:
println("Hello World!")

println('ggg')
println('hhhee')

println('nnn')

println('en plus')

max_stars_repo_name: nicolas-hbt/testGit

max_issues_repo_name: nicolas-hbt/testGit

max_forks_repo_name: nicolas-hbt/testGit


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ObjectReplicationPolicyRule <: SwaggerModel
    ruleId::Any # spec type: Union{ Nothing, String } # spec name: ruleId
    sourceContainer::Any # spec type: Union{ Nothing, String } # spec name: sourceContainer
    destinationContainer::Any # spec type: Union{ Nothing, String } # spec name: destinationContainer
    filters::Any # spec type: Union{ Nothing, ObjectReplicationPolicyFilter } # spec name: filters

    function ObjectReplicationPolicyRule(;ruleId=nothing, sourceContainer=nothing, destinationContainer=nothing, filters=nothing)
        o = new()
        validate_property(ObjectReplicationPolicyRule, Symbol("ruleId"), ruleId)
        setfield!(o, Symbol("ruleId"), ruleId)
        validate_property(ObjectReplicationPolicyRule, Symbol("sourceContainer"), sourceContainer)
        setfield!(o, Symbol("sourceContainer"), sourceContainer)
        validate_property(ObjectReplicationPolicyRule, Symbol("destinationContainer"), destinationContainer)
        setfield!(o, Symbol("destinationContainer"), destinationContainer)
        validate_property(ObjectReplicationPolicyRule, Symbol("filters"), filters)
        setfield!(o, Symbol("filters"), filters)
        o
    end
end # type ObjectReplicationPolicyRule

const _property_map_ObjectReplicationPolicyRule = Dict{Symbol,Symbol}(Symbol("ruleId")=>Symbol("ruleId"), Symbol("sourceContainer")=>Symbol("sourceContainer"), Symbol("destinationContainer")=>Symbol("destinationContainer"), Symbol("filters")=>Symbol("filters"))
const _property_types_ObjectReplicationPolicyRule = Dict{Symbol,String}(Symbol("ruleId")=>"String", Symbol("sourceContainer")=>"String", Symbol("destinationContainer")=>"String", Symbol("filters")=>"ObjectReplicationPolicyFilter")
Base.propertynames(::Type{ ObjectReplicationPolicyRule }) = collect(keys(_property_map_ObjectReplicationPolicyRule))
Swagger.property_type(::Type{ ObjectReplicationPolicyRule }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ObjectReplicationPolicyRule[name]))}
Swagger.field_name(::Type{ ObjectReplicationPolicyRule }, property_name::Symbol) =  _property_map_ObjectReplicationPolicyRule[property_name]

function check_required(o::ObjectReplicationPolicyRule)
    (getproperty(o, Symbol("sourceContainer")) === nothing) && (return false)
    (getproperty(o, Symbol("destinationContainer")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ObjectReplicationPolicyRule }, name::Symbol, val)
end


text:
#
# @delegate macro and support
#
using MacroTools

function delegated_args(args::Vector)
    newargs = []
    for a in args
        if a isa Symbol
            push!(newargs, a)

        elseif (@capture(a, var_::T_ = val_) || @capture(a, var_ = val_))
            push!(newargs, :($var = $var))

        elseif @capture(a, var_::T_)
            push!(newargs, var)
        else
            error("Unrecognized argument format: $a")
        end
    end
    return newargs
end

"""
Macro to define a method that simply delegate to a method with the same signature
but using the specified field name of the original first argument as the first arg
in the delegated call. That is,

    `@delegate compid(ci::CompositeComponentInstance, i::Int, f::Float64) => leaf`

expands to:

    `compid(ci::CompositeComponentInstance, i::Int, f::Float64) = compid(ci.leaf, i, f)`

If a second expression is given, it is spliced in, mainly to support the deprecated
decache(m)". We might delete this feature, but why bother?
"""
macro delegate(ex, other=nothing)
    result = nothing
    other = (other === nothing ? [] : [other])  # make vector so $(other...) disappears if empty

    if (@capture(ex, fname_(varname_::T_, args__; kwargs__) => rhs_) ||
        @capture(ex, fname_(varname_::T_, args__) => rhs_))
        # @info "args: $args"
        new_args = delegated_args(args)

        if kwargs === nothing
            kwargs = new_kwargs = []
        else
            new_kwargs = delegated_args(kwargs)
        end

        result = quote
            function $fname($varname::$T, $(args...); $(kwargs...))
                retval = $fname($varname.$rhs, $(new_args...); $(new_kwargs...))
                $(other...)
                return retval
            end
        end
    end

    if result === nothing
        error("Calls to @delegate must be of the form 'func(obj, args...) => X', where X is a field of obj to delegate to'. Expression was: $ex")
    end
    return esc(result)
end


text:
"""
    @test_forall argument_data, prop

Check statement of proposition over all tuples of values of variables, specified in `argument_data`. Generate Pass/Fail/Error output for Base.Test.

# Arguments
* `argument_data`: comma separated expression with arbitrarily many statements, containing information about ranges of variables. Statements must have form `x in iter`, where `x` is symbol and `iter` is any expression valid in `"for x in iter"` syntax used in construction of for-loops. In particular, `iter` can be `a:b`, where `a` & `b` are integer-valued expressions; or array/tuple-valued expression.
* `prop` is boolean-valued expression for proposition to be checked for arguments in specified ranges (see Examples below).

# Details
Code constructs a sequence of for-loops and checks that proposition evaluates to `true` for all values of variables specified. Equivalently, test is performed over cartesian product of ranges of variables.

Loops are generated from `argument_data` inductively: from outer-loop to inner-loop. If range of a variable depends on other variables, then their ranges must be defined in `argument_data` earlier.

#Examples
```julia
julia> @test_forall x in -1:1, x*(x-1)*(x+1) == 0
Test Passed
  Expression: (x in -1:1,x * (x - 1) * (x + 1) == 0)

julia> @test_forall x in -1:1, x*(x-2)*(x+1) == 0
Test Failed
  Expression: (x in -1:1,x * (x - 2) * (x + 1) == 0)

julia> @test_forall x in [0,1,2], y in x:4, (y+4>2*x)==true
Test Passed
  Expression: (x in [0,1,2],y in x:4,(y + 4 > 2x) == true)

julia> @test_forall x in (1-1,0+1,sqrt(4)), y in x:4, (y+4>2*x)==true
Test Passed
  Expression: (x in (1 - 1,0 + 1,sqrt(4)),y in x:4,(y + 4 > 2x) == true)

julia> @test_forall x in ["a","b"], y in ["z","w"], x*y in Set(["az","aw","bz","bw"])
Test Passed
  Expression: (x in ["a","b"],y in ["z","w"],x * y in Set(["az","aw","bz","bw"]))

# Bad design: throws UndefVarError: x not defined, because outer loop for y variable refers to the value of x coming after.
julia> @test_forall y in x:4,x in 0:2, (y+4>2*x)==true
Error During Test
  Test threw an exception of type UndefVarError
  Expression: (y in x:4,x in 0:2,(y + 4 > 2x) == true)
  UndefVarError: x not defined
```
"""
macro test_forall(factex::Expr)
  ex = begin
    try
    if factex.head != :tuple
        error("Expression of unsupported format for @test_forall: $factex")
    end
    # Defining ex
	iterate = :()
	for p in factex.args[1:end-1]
		iterate = Expr(:block,iterate.args...,:($(esc(p.args[2]))=$(esc(p.args[3]))))
	end
	prop = quote
		res = res && $(esc(factex.args[end]))
		if !res
			fail_data = $(esc(factex.args[end]))
			break
		end
	end
	ex = Expr(:for,iterate,prop)
  catch err
    ex =
        quote
		throw($err)
        end
  end # try
  ex
  end # begin
    #println(ex)
    quote
	res = true
	fail_data = false
	result = try
          $ex
          res ? Pass(:test,$(Expr(:quote, factex)), nothing, nothing) :
                      Fail(:test,$(Expr(:quote, factex)), fail_data, nothing)
        catch err
          Error(:test_error,$(Expr(:quote, factex)), err, catch_backtrace())
        end
        record(get_testset(), result)
    end
end


text:
module RobustFactorizations

using Statistics, LinearAlgebra
using ProximalOperators, TSVD

export rpca, rpca_fista, rpca_admm, lowrank, sparserandn

function lowrank(r,c,n)
    A = randn(r,c)
    U,S,V = svd(A)
    U[:,1:n]*Diagonal(S[1:n])*V[:,1:n]'
end

function sparserandn(r,c,γ=0.01)
    randn(r,c) .* (rand(r,c) .< γ)
end

struct RPCA{T<:AbstractMatrix}
    L::T
    S::T
    D::T
end

"""
    RPCA = rpca(A::Matrix; λ=1.0 / √(maximum(size(A))), iters=1000, tol=1.0e-7, ρ=1.5, verbose=false, nonnegL=false, nonnegS=false, nukeA=true)

minimize_{L,D,S} ||L||_* + λ||S||₁ + γ||D||₂² s.t. A = L+D+S

Ref: "The Augmented Lagrange Multiplier Method for Exact Recovery of Corrupted Low-Rank Matrices", Zhouchen Lin, Minming Chen, Leqin Wu, Yi Ma, https://people.eecs.berkeley.edu/~yima/psfile/Lin09-MP.pdf

#Arguments:
- `A`: Input matrix
- `λ`: Sparsity regularization
- `iters`: Maximum number of iterations
- `tol`: Tolerance
- `ρ`: Algorithm tuning param
- `verbose`: Print status
- `nonnegL`: Hard thresholding on A
- `nonnegS`: Hard thresholding on E
- `proxL`: NuclearNorm(1/2)
- `proxD`: nothing
- `proxS`: NormL1(λ))

To speed up convergence you may either increase the tolerance or increase `ρ`. Increasing `tol` is often the best solution.
"""
function rpca(A::AbstractMatrix{T};
                          λ              = real(T)(1.0/sqrt(maximum(size(A)))),
                          iters::Int     = 1000,
                          tol            = sqrt(eps(real(T))),
                          ρ              = real(T)(1.5),
                          verbose::Bool  = false,
                          nonnegL::Bool  = false,
                          nonnegS::Bool  = false,
                          proxL          = NuclearNorm(real(T)(nonnegL ? 1 : 1/2)),
                          proxD          = nothing,
                          proxS          = NormL1(λ)) where T
    RT        = real(T)
    M, N      = size(A)
    d         = min(M,N)
    L, S      = zeros(T, M, N), zeros(T, M, N)
    D         = zeros(T, M, N)
    Z         = similar(A)
    Y         = copy(A)
    norm²     = opnorm(Y)::RT # can be tuned
    norm∞     = norm(Y, Inf) / λ
    dual_norm = max(norm², norm∞)
    d_norm    = norm²
    Y       ./= dual_norm
    μ         = RT(1.25) / norm²
    μ̄         = μ  * RT(1.0e+7)
    sv        = 10
    for k = 1:iters
        if proxD !== nothing
            Z .= A .- L .- S .+ (1/μ) .* Y
            prox!(D, proxD, Z, 1/μ)
        end
        Z .= A .- L .- D .+ (1/μ) .* Y
        prox!(S, proxS, Z, 1/μ)
        if nonnegS
            S .= max.(S, 0)
        end
        Z .= A .- S .- D .+ (1/μ) .* Y
        prox!(L, proxL, Z, 1/μ)
        if nonnegL
            L .= max.(L, 0)
        end

        @. Z = A - L - S - D # Z are the reconstruction errors
        @. Y = Y + μ * Z # Y can not be moved below as it depends on μ which is changed below
        μ = min(μ*ρ, μ̄)

        cost = opnorm(Z) / d_norm
        verbose && println("$(k) cost: $(round(cost, sigdigits=4))")

        if cost < tol
            verbose && println("converged")
            break
        end
        k == iters && @warn "Maximum number of iterations reached, cost: $cost, tol: $tol"
    end

    RPCA(L, S, D)
end

function rpca_fista(A::AbstractMatrix{T};
                          λ              = real(T)(1.0/sqrt(maximum(size(A)))),
                          γ              = 0.5,
                          iters::Int     = 10000,
                          tol            = sqrt(eps(real(T))),
                          verbose::Bool  = false,
                          printerval     = 50,
                          nonnegL::Bool  = false,
                          nonnegS::Bool  = false,
                          proxL          = NuclearNorm(real(T)(nonnegL ? 1 : 1/2)),
                          proxD          = nothing,
                          proxS          = NormL1(λ)) where T

    RT             = real(T)
    M, N           = size(A)
    d              = min(M,N)
    L, S           = zeros(T, M, N), zeros(T, M, N)
    YL, YS         = zeros(T, M, N), zeros(T, M, N)
    L_extr, S_extr = zeros(T, M, N), zeros(T, M, N)
    D              = zeros(T, M, N)
    Z              = similar(A)
    S_prev         = copy(S)
    L_prev         = copy(L)
    prox           = SeparableSum(proxS, proxL)
    t              = RT(1.0)
    for k = 1:iters
        t1 = (1 + sqrt(1+4t^2))/2
        acc = (t-1)/t1
        t = t1
        S_extr .= S .+ acc.*(S .- S_prev)
        L_extr .= L .+ acc.*(L .- L_prev)
        @. D = A - S - L
        @. YS = S_extr + γ*D
        @. YL = L_extr + γ*D
        S_prev .= S
        L_prev .= L
        prox!((S, L), prox, (YS, YL), γ)
        if nonnegS
            S .= max.(S, 0)
        end
        if nonnegL
            L .= max.(L, 0)
        end
        ϵ = max(norm(S_extr-S, Inf), norm(L_extr-L, Inf))/γ
        cost = ϵ/(1+max(norm(S,Inf), norm(L,Inf)))
        verbose && k % printerval == 0 &&  println("$(k) cost: $(round(cost, sigdigits=4))")
        if cost <= tol
            verbose && println("converged, cost: $(round(cost, sigdigits=4))")
            break
        end
        k == iters && @warn "Maximum number of iterations reached, cost: $cost, tol: $tol"
    end
    return RPCA(L,S,D)
end

function rpca_admm(A::AbstractMatrix{T};
                          λ              = real(T)(1.0/sqrt(maximum(size(A)))),
                          ρ              = 0.5,
                          iters::Int     = 10000,
                          printerval     = 100,
                          tol            = sqrt(eps(real(T))),
                          verbose::Bool  = false,
                          nonnegL::Bool  = false,
                          nonnegS::Bool  = false,
                          proxL          = NuclearNorm(real(T)(nonnegL ? 1 : 1/2)),
                          proxD          = nothing,
                          proxS          = NormL1(λ)) where T

  d_norm         = opnorm(A)
  RT             = real(T)
  RT             = real(T)
  M, N           = size(A)
  d              = min(M,N)
  L, S           = zeros(T, M, N), zeros(T, M, N)
  Y              = zeros(T, M, N)
  D              = zeros(T, M, N)
  for k = 0:iters-1
      @. Y = A - S + D
      prox!(L, proxL, Y, 1/ρ)
      if nonnegL
          L .= max.(L, 0)
      end
      @. Y = A - L + D
      prox!(S, proxS, Y, 1/ρ)
      if nonnegS
          S .= max.(S, 0)
      end
      @. D = D + A - L - S
      cost = opnorm(D) / d_norm
      verbose && k % printerval == 0 && println("$(k) cost: $(round(cost, sigdigits=4))")
      if cost <= tol
          verbose && println("converged")
          break
      end
      k == iters && @warn "Maximum number of iterations reached, cost: $cost, tol: $tol"
  end
  return RPCA(L,S,D)
end

end # module


text:
import RowEchelon
import LinearAlgebra

#=
Homework 7
Section 3.3
Problems 19
=#


text:
module Day11

using ..ReTest, ..DocStringExtensions

import ..solve, ..PuzzleInfo, ..TEST_DATA_DIRECTORY

abstract type Day11PuzzleInfo <: PuzzleInfo end

"""
Day 11 part 1 puzzle information.

$(FIELDS)
"""
struct Part1PuzzleInfo <: Day11PuzzleInfo
    """Path to input data file."""
    datapath::String
end

"""
Day 11 part 2 puzzle information.

$(FIELDS)
"""
struct Part2PuzzleInfo <: Day11PuzzleInfo
    """Path to input data file."""
    datapath::String
end

"""
    $(FUNCTIONNAME)(input)

Read two-dimensional integer (1:9) energy level data from input stream `input`.
"""
function readenergylevels(input::IO)
    reduce(vcat, [parse(Int, c) for c in line]' for line in eachline(input))
end

"""
    $(FUNCTIONNAME)(index, shape)

Get iterator over valid offset indices at index `index` for 2D array with size `shape`.
"""
function offsetsat(index::CartesianIndex{2}, shape::Tuple{Integer, Integer})
    [
        CartesianIndex(i, j)
        for i in (index[1] > 1 ? -1 : 0):(index[1] < shape[1] ? 1 : 0)
        for j in (index[2] > 1 ? -1 : 0):(index[2] < shape[2] ? 1 : 0)
    ]
end

"""
    $(FUNCTIONNAME)(energylevels)

Step dumbo octopus flashing model forward given current energy levels `energylevels`.
Returns the set of indices corresponding to octopi which flashed during step.
"""
function stepmodel!(energylevels::Array{T, 2}) where T <: Integer
    energylevels .+= 1
    flashindices = findall(energylevels .> 9)
    energylevels[flashindices] .= 0
    flashed = Set(flashindices)
    while length(flashindices) > 0
        for flashindex in flashindices
            for offset in offsetsat(flashindex, size(energylevels))
                adjacentindex = flashindex + offset
                if adjacentindex ∉ flashed
                    energylevels[flashindex + offset] += 1
                end
            end
        end
        flashindices = findall(energylevels .> 9)
        energylevels[flashindices] .= 0
        union!(flashed, flashindices)
    end
    flashed
end

function solve(info::Part1PuzzleInfo)
    energylevels = open(readenergylevels, info.datapath)
    numflashes = 0
    for step in 1:100
        flashed = stepmodel!(energylevels)
        numflashes += length(flashed)
    end
    numflashes
end

function solve(info::Part2PuzzleInfo)
    energylevels = open(readenergylevels, info.datapath)
    step = 0
    while !all(energylevels .== 0)
        stepmodel!(energylevels)
        step += 1
    end
    step
end

@testset begin
    @test solve(Part1PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day11.txt"))) == 1656
    @test solve(Part2PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day11.txt"))) == 195
end

end # module


text:
# syntax: proto3
using Compat
using ProtoBuf
import ProtoBuf.meta
import Base: hash, isequal, ==

type _Any
    type_url::AbstractString
    value::Array{UInt8,1}
    _Any(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #type _Any
hash(v::_Any) = ProtoBuf.protohash(v)
isequal(v1::_Any, v2::_Any) = ProtoBuf.protoisequal(v1, v2)
==(v1::_Any, v2::_Any) = ProtoBuf.protoeq(v1, v2)

export _Any


text:
# ic.jl
# Needed to initialize a problem.

function ICConstant(mesh::AbstractMesh, sbp::AbstractSBP, eqn::SimpleODEData, opts,
                    u0::AbstractArray)

  for i=1:length(u0)
    u0[i] = 2
  end

  return nothing
end

@doc """
### SimpleODEMod.ICx2plust4

Computes the initial conditions for the state variable
    u = x^2 + t^4

**Inputs**

*  `mesh` : AbstractMesh type
*  `sbp`  : Summation-by-parts operator
*  `eqn`  : Simple ODE equation object
*  `opts` : Options dictionary
*  `u0`   : Array that stores inital state variables

**Outputs**

*  None

"""->
function ICx2plust4(mesh::AbstractMesh{Tmsh},
sbp::AbstractSBP{Tsbp}, eqn::SimpleODEData{Tsol, Tres, 2},
opts, u0::AbstractArray{Tsol}) where {Tmsh, Tsbp, Tsol, Tres}

  t = 0.0

  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      dofnums_j = sview(mesh.dofs, :, j, i)
#       x = mesh.coords[1, j, i]
#       y = mesh.coords[2, j, i]
#       u0[dofnums_j] = x^2 + t^2
      u0[dofnums_j] = calc_x2_4t3(mesh.coords[:, j, i], eqn.params, t)
    end
  end

  return nothing
end # end function ICx2plust4

@doc """
### SimpleODEMod.ICx2plust3

Computes the initial conditions for the state variable
    u = x^2 + t^3

**Inputs**

*  `mesh` : AbstractMesh type
*  `sbp`  : Summation-by-parts operator
*  `eqn`  : Simple ODE equation object
*  `opts` : Options dictionary
*  `u0`   : Array that stores inital state variables

**Outputs**

*  None

"""->
function ICx2plust3(mesh::AbstractMesh{Tmsh},
sbp::AbstractSBP{Tsbp}, eqn::SimpleODEData{Tsol, Tres, 2},
opts, u0::AbstractArray{Tsol}) where {Tmsh, Tsbp, Tsol, Tres}

  t = 0.0

  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      dofnums_j = sview(mesh.dofs, :, j, i)
#       x = mesh.coords[1, j, i]
#       y = mesh.coords[2, j, i]
#       u0[dofnums_j] = x^2 + t^2
      u0[dofnums_j] = calc_x2_3t2(mesh.coords[:, j, i], eqn.params, t)
    end
  end

  return nothing
end # end function ICx2plust3

@doc """
### SimpleODEMod.ICx2plust2

Computes the initial conditions for the state variable
    u = x^2 + t^2

**Inputs**

*  `mesh` : AbstractMesh type
*  `sbp`  : Summation-by-parts operator
*  `eqn`  : Simple ODE equation object
*  `opts` : Options dictionary
*  `u0`   : Array that stores inital state variables

**Outputs**

*  None

"""->

function ICx2plust2(mesh::AbstractMesh{Tmsh},
sbp::AbstractSBP{Tsbp}, eqn::SimpleODEData{Tsol, Tres, 2},
opts, u0::AbstractArray{Tsol}) where {Tmsh, Tsbp, Tsol, Tres}

  t = 0.0

  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      dofnums_j = sview(mesh.dofs, :, j, i)
#       x = mesh.coords[1, j, i]
#       y = mesh.coords[2, j, i]
#       u0[dofnums_j] = x^2 + t^2
      u0[dofnums_j] = calc_x2_t2(mesh.coords[:, j, i], eqn.params, t)
    end
  end

  return nothing
end # end function ICx2plust2

@doc """
### SimpleODEMod.ICx2plust2

Computes the initial conditions for the state variable
    u = x^2 + t^2

**Inputs**

*  `mesh` : AbstractMesh type
*  `sbp`  : Summation-by-parts operator
*  `eqn`  : Simple ODE equation object
*  `opts` : Options dictionary
*  `u0`   : Array that stores inital state variables

**Outputs**

*  None

"""->

function ICallzero(mesh::AbstractMesh{Tmsh},
sbp::AbstractSBP{Tsbp}, eqn::SimpleODEData{Tsol, Tres, 2},
opts, u0::AbstractArray{Tsol}) where {Tmsh, Tsbp, Tsol, Tres}
  t = 0.0

  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      dofnums_j = sview(mesh.dofs, :, j, i)
#       x = mesh.coords[1, j, i]
#       y = mesh.coords[2, j, i]
#       u0[dofnums_j] = x^2 + t^2
      u0[dofnums_j] = 0.0
    end
  end

  return nothing
end # end function ICallzero

function ICallOnes(mesh::AbstractMesh{Tmsh},
sbp::AbstractSBP{Tsbp}, eqn::SimpleODEData{Tsol, Tres, 2},
opts, u0::AbstractArray{Tsol}) where {Tmsh, Tsbp, Tsol, Tres}
  t = 0.0

  for i = 1:mesh.numEl
    for j = 1:mesh.numNodesPerElement
      dofnums_j = sview(mesh.dofs, :, j, i)
#       x = mesh.coords[1, j, i]
#       y = mesh.coords[2, j, i]
#       u0[dofnums_j] = x^2 + t^2
      u0[dofnums_j] = 1.0
    end
  end

  return nothing
end # end function ICallOnes

# TODO: rename ICallzero to ICallZeros

global const ICDict = Dict{Any, Function}(
"ICConstant" => ICConstant,
"ICx2plust4" => ICx2plust4,
"ICx2plust3" => ICx2plust3,
"ICx2plust2" => ICx2plust2,
"ICallzero" => ICallzero,
"ICallOnes" => ICallOnes,
)


text:
using Documenter, DimensionalData

makedocs(
    modules = [DimensionalData],
    sitename = "DimensionalData.jl",
)

deploydocs(
    repo = "github.com/rafaqz/DimensionalData.jl.git",
)


text:
function RTAU(DP)
      if (DP<=25.0)
          return rtauv=1.0
      elseif (DP>=50.0)
          return rtauv=0.1
      else
          return rtauv=1.9-.036*DP
      end
end


text:
module Bifurcation

using PyPlot

include("model/model.jl")
using .Model
include("./differential_equation.jl")
include("bifurcation/continuation.jl")

function analysis()
    global fp, br

    # BP: name(index) of bifurcation parameter
    global BP = C.mitogen

    global p = param_values()

    new_curve!(p)

    fp = readdlm("./Data/fp.dat",'\t',Float64,'\n')
    ev = readdlm("./Data/ev.dat",'\t',Float64,'\n')
    br = get_bistable_regime(ev)
end

function diagram()
    rc("figure",figsize = (8,6))
    rc("font",family = "Arial")
    rc("font",size = 20)
    rc("axes",linewidth = 1)
    rc("lines",linewidth = 2)
    rc("lines",markersize = 10)
    gca().spines["right"].set_visible(false)
    gca().spines["top"].set_visible(false)
    gca().yaxis.set_ticks_position("left")
    gca().xaxis.set_ticks_position("bottom")

    plot(fp[1:br[1]-1,VN+1],fp[1:br[1]-1,V.E2f+1],"k-")
    plot(fp[br,VN+1],fp[br,V.E2f+1],"k:")
    plot(fp[br[end]+1:end,VN+1],fp[br[end]+1:end,V.E2f+1],"k-")

    axvspan(fp[br[1],VN+1],fp[br[end],VN+1],facecolor="m",alpha=0.05)

    xlabel("Mitogen (a.u.)")
    ylabel("E2F activity (a.u.)")

    xlim(0,1)
    xticks([0,0.5,1])
    ylim(-0.05,1.05)
    yticks([0,0.5,1])

    show()
end

end  # module


text:
@testset "buffer" begin

data = rand(UInt8, 8)

let
    membuf = MemoryBuffer(data)
    dispose(membuf)
end

@test_throws LLVMException MemoryBufferFile("nonexisting")

MemoryBuffer(data) do membuf
    @test pointer(data) != pointer(membuf)
    @test length(membuf) == length(data)
    @test convert(Vector{UInt8}, membuf) == data
end

MemoryBuffer(data, "SomeBuffer", false) do membuf
    @test pointer(data) == pointer(membuf)
end

mktemp() do path, _
    let
        membuf = MemoryBufferFile(path)
        dispose(membuf)
    end

    MemoryBufferFile(path) do membuf
    end
end

end


text:
module Dante
# Reimplementation of finite volume MHD in Julia.

export solve, setInitRiemann, EulerExact

using Printf, TimerOutputs

# Set main timer for global timing of functions
const main_timer = TimerOutput()

# Always call timer() to hide implementation details
timer() = main_timer

include("utility.jl")
include("input.jl")
include("advance.jl")
include("divergence.jl")
include("facevalue.jl")
include("flux.jl")
include("state.jl")
include("boundary.jl")
include("source.jl")
include("time.jl")
include("output.jl")

"""
	solve(paramFile)

Run the model given input parameter file `paramFile`.
"""
function solve(paramFile="examples/PARAM_sods.toml")

	reset_timer!(timer())

	@timeit timer() "input" param = setParameters(paramFile)

	@timeit timer() "init" state_GV = set_init(param)

	@timeit timer() "advance" advance!(param, state_GV)

	# Compare with analytical solution, if possible
	if param.DoPlot && param.IC == "Riemann"
		@timeit timer() "analytic sol" plot_Riemann_exact(param)
	end

	show(timer())
	println()

	# Check if running test
	if occursin("test", paramFile)
		return param, state_GV
	end
end

end


text:
module Calc

import Base
import REPL
import REPL.LineEdit
import REPL.LineEdit: terminal
import Statistics

mutable struct CalcStack <: AbstractArray{Any, 1}
    x::Array{Any}
end

CalcStack() = CalcStack(Any[])
Base.getindex(s::CalcStack, i) = s.x[i]
Base.setindex!(s::CalcStack, x, i) = setindex!(s.x, x, i)
Base.size(s::CalcStack, o...) = size(s.x, o...)
Base.similar(s::CalcStack) = CalcStack()
Base.push!(s::CalcStack, x) = push!(s.x, x)
Base.splice!(s::CalcStack, i) = splice!(s.x, i)
Base.copy(s::CalcStack) = CalcStack(copy(s.x))
function Base.show(io::IO, s::CalcStack)
    println(io)
    println(io, "Stack [$(state.usedegrees ? "deg" : "rad")|$(state.usepolar ? "polr" : "rect")]")
    n = length(s)
    for i in 1:n
        print(io, n-i+1, ": ")
        printelement(io, s.x[i])
        println(io)
    end
    if n == 0
        println(io, ".")
    end
end

cs(x) = sprint(showcompact, x)
printelement(io::IO, x) = show(IOContext(io, :compact => true), x)
printelement(io::IO, x::Complex) = state.usepolar ? print(io, "$(cs(abs(x)))∠$(cs(rad2deg(angle(x))))°") : show(IOContext(io, :compact => true), x)

mutable struct CalcState
    history::Array{CalcStack}
    panel
    position::Int
    usedegrees::Bool
    usepolar::Bool
end

const state = CalcState(CalcStack[CalcStack()], 0, 1, true, false)

activestack() = state.history[state.position]

function advance(stack)
    if stack != activestack()
        state.position += 1
        if state.position > length(state.history)
             push!(state.history, stack)
        else
             state.history[state.position] = stack
             if state.position < length(state.history)
                 state.history = state.history[1:state.position]
             end
        end
    end
end

"""
    calcfun(fun, n = 0, splatoutput = false)

where

* `fun` : defines the operation using the stack.
* `n`   : number of stack elements to pass to `fun`. If negative, pass the whole
          stack.
* `splatoutput` : if `true`, splat the results back on to the stack. This is
          useful for vector operations and entering multiple values on the
          stack.

As part of this operation, `n` elements are removed from the stack.

With one argument (`n == 1`), the top of the stack `x` is passed to `fun`.
With `n > 1`, the top `n` stack elements are passed to `fun` in reverse order.
By convention, `x` is the top of the stack, and `y` is the second element, so
`function(y, x)` is the appropriate order for two arguments.

Here are examples taken from the core key definitions using `setkeys`:

```julia
Calc.setkeys(Dict(
    "P"  => Calc.calcfun(() -> pi, 0),
    "Q"  => Calc.calcfun(sqrt, 1),
    "-"  => Calc.calcfun(-, 2),
    "n"  => Calc.calcfun(-, 1),   # negate
    # tab - swap x & y on the stack
    "\t" => Calc.calcfun((y, x) -> Any[x, y], 2, true),
    "fs" => Calc.calcfun(sign, 1),
    # pack the stack into a vector
    "Vp" => Calc.calcfun(x -> Any[x], -1)
))
```
Returns a function for use in a keymap.
"""
function calcfun(fun, n = 0, splatoutput = false)
    (s, args...) -> begin
        println(terminal(s))
        stack = copy(activestack())
        b = LineEdit.buffer(s)
        newval = Base.eval(Main, Base.parse_input_line(String(take!(b))))
        if newval != nothing
            push!(stack, newval)
            advance(stack)
            stack = copy(stack)
        end
        if n ≥ 0
            ns = length(stack)
            args = splice!(stack, ns-n+1:ns)
            val = fun(args...)
            if val != nothing
                if splatoutput
                    push!(stack, val...)
                else
                    push!(stack, val)
                end
            end
        else       # Negative: pass and return the whole stack
            stack.x = fun(stack.x)
        end
        advance(stack)
        show(terminal(s), activestack())
        :done
    end
end

## This code was adapted from the following:
##    https://github.com/Keno/Gallium.jl/blob/1ef7111880495f3c5a7989d88a5fbc60e7eeb285/src/lldbrepl.jl
## Copyright (c) 2015: Keno Fischer. Licensed under the MIT "Expat" License:
##    https://github.com/Keno/Gallium.jl/blob/b4bc668a4cbd0f2d4f63fbdb0597a1264afd7b4d/LICENSE.md

function initiate_calc_repl(repl)
    mirepl = isdefined(repl,:mi) ? repl.mi : repl

    main_mode = mirepl.interface.modes[1]

    calc_launch_keymap = Dict{Any,Any}(
        '=' => function (s,args...)
            if isempty(s)
                if !haskey(s.mode_state, panel)
                    s.mode_state[panel] = LineEdit.init_state(repl.t, panel)
                end
                println(terminal(s))
                println(terminal(s))
                show(terminal(s), activestack())
                println(terminal(s))
                LineEdit.transition(s,panel)
            else
                LineEdit.edit_insert(s,'=')
            end
        end
    )

    # Setup the repl panel
    panel = LineEdit.Prompt("calc> ";
        # Copy colors from the prompt object
        prompt_prefix = Base.text_colors[:blue],
        prompt_suffix = main_mode.prompt_suffix,
        on_enter = REPL.return_callback)
    panel.on_done = REPL.respond(repl, panel; pass_empty = false) do line
        :(  )
    end
    # Setup the alternate repl panel for algebraic entry
    # Setup the panel for input
    inputpanel = LineEdit.Prompt(">>>>> ";
        prompt_prefix = Base.text_colors[:cyan],
        prompt_suffix = main_mode.prompt_suffix,
        on_enter = s -> true)
    inputpanel.on_done = REPL.respond(repl, panel; pass_empty = false) do line
        isempty(line) ? :() : line
    end

    hp = main_mode.hist
    hp.mode_mapping[:calc] = panel
    panel.hist = hp

    push!(mirepl.interface.modes, panel)

    search_prompt, skeymap = LineEdit.setup_search_keymap(hp)
    mk = REPL.mode_keymap(main_mode)

    b = Dict{Any,Any}[skeymap, mk, LineEdit.history_keymap, LineEdit.default_keymap, LineEdit.escape_defaults]

    function input(fun::Function, s, prompt::AbstractString)
        inputpanel.prompt = prompt
        inputpanel.on_done = REPL.respond(repl, panel; pass_empty = false) do line
            :( $(try fun(line) catch e warn(e) end) )
        end
        if !haskey(s.mode_state, inputpanel)
            s.mode_state[inputpanel] = LineEdit.init_state(repl.t, inputpanel)
        end
        LineEdit.transition(s, inputpanel)
    end

    calc_keymap = Dict{Any,Any}(
    # arithmetic
        "+" => calcfun(+, 2),
        "-" => calcfun(-, 2),
        "n" => calcfun(-, 1),   # negate
        "*" => calcfun(*, 2),
        "/" => calcfun(/, 2),
        "&" => calcfun(x -> 1/x, 1),
        "%" => calcfun(%, 2),
        "A" => calcfun(abs, 1),
        "fs" => calcfun(sign, 1),
        "fn" => calcfun(min, 2),
        "fx" => calcfun(max, 2),
        "f[" => calcfun(x -> x - 1, 1),
        "f]" => calcfun(x -> x + 1, 1),
    # algebraic
        "Q" => calcfun(sqrt, 1),
        "IQ" => calcfun(x -> x*x, 1),
        "L" => calcfun(log, 1),
        "E" => calcfun(exp, 1),
        "IL" => calcfun(exp, 1),
        "HL" => calcfun(log10, 1),
        "IHL" => calcfun(x -> 10^x, 1),
        "B" => calcfun((y, x) -> log(x, y), 2),
        "^" => calcfun(^, 2),
        "I^" => calcfun((y, x) -> y ^ (1/x), 2),
        "fh" => calcfun((y, x) -> sqrt(y^2 + x^2), 2),
    # trig
        "S" => calcfun(x -> state.usedegrees ? sind(x) : sin(x), 1),
        "C" => calcfun(x -> state.usedegrees ? cosd(x) : cos(x), 1),
        "T" => calcfun(x -> state.usedegrees ? tand(x) : tan(x), 1),
        "IS" => calcfun(x -> state.usedegrees ? asind(x) : asin(x), 1),
        "IC" => calcfun(x -> state.usedegrees ? acosd(x) : acos(x), 1),
        "IT" => calcfun(x -> state.usedegrees ? atand(x) : atan(x), 1),
        "P" => calcfun(() -> pi, 0),
    # settings
        "mr" =>  (s, o...) -> (println(terminal(s), "\nUsing radians..."); state.usedegrees = false; :done),
        "md" =>  (s, o...) -> (println(terminal(s), "\nUsing degrees..."); state.usedegrees = true; :done),
        "mp" =>  (s, o...) -> (state.usepolar = !state.usepolar; println(terminal(s), "\nUsing $(state.usepolar ? "polar" : "rectangular") coordinates..."); :done),
    # complex numbers
        "X" => calcfun(complex, 2),
        "IX" => calcfun(x -> [real(x), imag(x)], 1, true),
        # polar entry with y in degrees
        "Z" => calcfun((y, x) -> y * exp(1.0im * x * π / 180), 2),
        "IZ" => calcfun(x -> [abs(x), rad2deg(angle(x))], 1, true),
        "J" => calcfun(conj, 1),
        "G" => calcfun(x -> state.usedegrees ? rad2deg(angle(x)) : angle(x), 1),
        "fr" => calcfun(real, 1),
        "fi" => calcfun(imag, 1),
    # percentages
        "\e%" => calcfun(x -> x/100, 1),
        "c%" => calcfun(x -> 100x, 1),
        "b%" => calcfun((y, x) -> 100(x-y)/y, 2),
    # vectors
        "|" => calcfun(vcat, 2),
        "Vu" => calcfun(x -> x, 1, true),
        "Vp" => calcfun(x -> Any[x], -1),
    # statistics
        "u#" => calcfun(length, 1),
        "u+" => calcfun(sum, 1),
        "u*" => calcfun(prod, 1),
        "uX" => calcfun(maximum, 1),
        "uN" => calcfun(minimum, 1),
        "uM" => calcfun(Statistics.mean, 1),
        "HuM" => calcfun(Statistics.median, 1),
        "uS" => calcfun(Statistics.std, 1),
        "HuS" => calcfun(Statistics.var, 1),
    # storing/recalling
        # store x in the prompted variable
        "ss" => (s, o...) -> input(s, "Variable name> ") do x
                    Base.eval(Main, Expr(:(=), Symbol(x), activestack()[end]))
                end,
        # store copy of the whole stack in the prompted variable as Vector{Any}
        "sS" => (s, o...) -> input(s, "Variable name> ") do x
                    Base.eval(Main, Expr(:(=), Symbol(x), copy(activestack().x)))
                end,
    # general
        # Meta-k - Copy `x` to the clipboard
        "\ek" => calcfun(x -> (clipboard(x); x), 1),
        # Ctrl-k - Pop `x` to the clipboard -- not sure why "^K" doesn't work
        "\u0b"  => calcfun(x -> (clipboard(x); return nothing), 1),
        # delete
        "\e[3~" => (s,o...)-> eof(LineEdit.buffer(s)) ? calcfun((stack) -> stack[1:end-1], -1)(s,o...) : LineEdit.edit_delete(s),
        # tab - swap x & y on the stack
        "\t" => calcfun((y, x) -> Any[x, y], 2, true),
        # space / Enter for stack entry
        " " => calcfun(x -> x, -1),
        "\r" => LineEdit.KeyAlias(" "),
        "\n" => LineEdit.KeyAlias(" "),
        # undo
        "U" => (s, o...) -> begin
                    if state.position > 1
                        state.position -= 1
                    end
                    show(terminal(s), activestack())
                    :done
                end,
        # redo
        "D" => (s, o...) -> begin
                    if state.position < length(state.history)
                        state.position += 1
                    end
                    show(terminal(s), activestack())
                    :done
                end,
        # trigger algebraic entry
        "=" => (s, o...) -> input(s, "calc= ") do line
                    stack = copy(Calc.activestack())
                    push!(stack, Base.eval(Main, fixrefs(Base.parse_input_line(line))))
                    Calc.advance(stack)
                    show(terminal(s), activestack())
                end
    )

    main_mode.keymap_dict = LineEdit.keymap_merge(main_mode.keymap_dict, calc_launch_keymap)
    panel.keymap_dict = LineEdit.keymap_merge(LineEdit.keymap(b), calc_keymap)
    state.panel = panel

    # Finish the input repl panel
    hp.mode_mapping[:input] = inputpanel
    inputpanel.hist = hp

   # Convert _1, _2, ... to stack references
    fixrefs(x) = x
    fixrefs(e::Expr) = Expr(e.head, Any[fixrefs(a) for a in e.args]...)
    function fixrefs(s::Symbol)
        st = string(s)
        if first(st) != '_'
            return s
        end
        n = try
            parse(st[2:end])
        catch
            ""
        end
        if isa(n, Integer) && n > 0
            return :( $activestack()[end - $n + 1] )
        else
            return s
        end
    end

    push!(mirepl.interface.modes, inputpanel)
    inputpanel.keymap_dict = LineEdit.keymap(b)

    nothing
end

"""
    setkeys(keymap::Dict)

Merge `keymap` with the Calc keymap. Use to define or replace keys in the
calculator. `keymap` is a mapping from a key sequence to a function that
defines the operation. Keys and key sequences are defined using Unix-style
key definitions (I haven't found a good link). Use `^` for a CTRL key modifier
and `\e` for a Meta or ALT modifier.

See `Calc.calcfun` for a helper function for defining operations on the stack.
Here is an example:

```julia
Calc.setkeys(Dict("fp" => Calc.calcfun((y, x) -> 1 / (1/y + 1/x), 2)))
```
Returns the new keymap.
"""
function setkeys(keymap)
    state.panel.keymap_dict = LineEdit.keymap_merge(state.panel.keymap_dict, keymap)
end

function __init__()
    isdefined(Base, :active_repl) && initiate_calc_repl(Base.active_repl)
end

end # module


text:
using CatIndices
using Test

### BidirectionalVector
@testset "CatIndices" begin
    function checkvals(v)
        for i in axes(v,1)
            @test v[i] == i
        end
        nothing
    end
    @testset "BidirectionalVector" begin
        v = BidirectionalVector(1:3)
        @test axes(v,1) == 1:3
        checkvals(v)
        @test isa(push!(v, 4), BidirectionalVector)
        @test axes(v,1) == 1:4
        checkvals(v)
        @test isa(append!(v, 5:11), BidirectionalVector)
        @test axes(v,1) == 1:11
        checkvals(v)
        @test pop!(v) == 11
        @test axes(v,1) == 1:10
        checkvals(v)
        @test isa(deletetail!(v, 2), BidirectionalVector)
        @test axes(v,1) == 1:8
        checkvals(v)
        @test popfirst!(v) == 1
        @test axes(v,1) == 2:8
        checkvals(v)
        @test isa(pushfirst!(v, 0, 1), BidirectionalVector)
        @test isa(pushfirst!(v, -1), BidirectionalVector)
        @test axes(v,1) == -1:8
        checkvals(v)
        @test isa(prepend!(v, -5:-2), BidirectionalVector)
        @test axes(v,1) == -5:8
        checkvals(v)
        @test isa(deletehead!(v,2), BidirectionalVector)
        @test axes(v,1) == -3:8
        checkvals(v)

        v[0] = 200
        @test v[0] == 200

        @test axes(similar(v)) === axes(v)
        @test axes(similar(v, Float64)) === axes(v)
        @test axes(similar(v, Float64, 3)) === (Base.OneTo(3),)
        a = similar(Array{Float64}, axes(v))
        @test isa(a, BidirectionalVector)
        @test axes(a) === axes(v)
    end

    @testset "vcat" begin
        @test !CatIndices.is_pinned(1:3)
        @test CatIndices.is_pinned(PinIndices(1:3))
        v = vcat(1:3, PinIndices(4:5), 6:10)
        @test axes(v,1) == -2:7
        @test_throws ArgumentError vcat(1:3, PinIndices(4:5), PinIndices(6:10))
    end
end