text:
#
# MCTS (Monte Carlo Tree Search) Choice Model
#
# Note that this implementation has the following features:
#
#   - the "reward" function is assumed to provide a value from 0 to 1, with 1 the best value (other limits will probably work,
#     though); note that higher being better is the opposite direction from the NMCS implementation's "fitness" function
#
#   - since in GödelTest somce choice points have infinite support, it is not possible to exhaustively enumerate all children
#	  therefore we sample (with replacement) according to the policy choice model (the number of samoples controlled by the
#     samplesize parameter); this is for consistency how this "finitization" problem is handled with the NMCS implementation;
#     see the julia3_spikes directory for alternative implementations that (in a prototypical manner) finitise by simply placing
#     an upper bound on the number of children and so can enumerate without sampling [note: these spikes Julia 3 code, and may
#     also have some other functional differences]
#
#	- while the standard UCB1 metric is used to control exploration vs exploitation, the godel number (the next "move") is returned
#     by taking the child node of the root that found the best reward (rather than best average reward); this is because GödelTes
#     is a single-player "game" and so we have the control to realise this reward
#

type MCTSNode
	trace::Dict
	totalreward::Real
	maxreward::Real
	visitcount::Int
	samplecount::Int
	childnodes::Dict{Real, MCTSNode}
	function MCTSNode(trace::Dict)
		new(trace, 0.0, -Inf, 0, 0, Dict{Real,MCTSNode}())
	end
end

#DEBUG function printnode(node::MCTSNode, action::AbstractString="root", indent::Int=0)
	#DEBUG println(" "^(indent*2) * "$(action): totalreward: $(node.totalreward) maxreward: $(node.maxreward) visitcount: $(node.visitcount) samplecount: $(node.samplecount)")
#DEBUG	for godelnumber in keys(node.childnodes)
#DEBUG		printnode(node.childnodes[godelnumber], "$(godelnumber)", indent+1)
#DEBUG 	end
#DEBUG end

type MCTSChoiceModel <: ChoiceModel
	policychoicemodel::ChoiceModel
	rewardfunction::Function						# higher is better
	visitbudget::Int
	samplesize::Int 								# the number of samples to take
	cp::Real										# the constant C_p in the calculation of UCB1
	rootgodelsequence::Vector{Real} 				# defines the state of the current "root" node
	rootnode::MCTSNode 								# current root node
	function MCTSChoiceModel(policychoicemodel::ChoiceModel, rewardfunction::Function, visitbudget::Int=1, samplesize::Int=1, cp::Real=1/sqrt(2.0))
		new(deepcopy(policychoicemodel), rewardfunction, visitbudget, samplesize, cp, (Real)[], MCTSNode(Dict()))
	end
end

function MCSTChoiceModel(g::Generator, rewardfunction::Function, visitbudget::Int=1, samplesize::Int=1, cp::Real=1/sqrt(2.0))
	MCSTChoiceModel(choicemodel(g), rewardfunction, visitbudget, samplesize, cp)
end

function setmctschoicemodel!(g::Generator, rewardfunction::Function, visitbudget::Int=1, samplesize::Int=1, cp::Real=1/sqrt(2.0))
	Base.depwarn("Use setchoicemodel!(g, MCTSChoiceModel(g, ...)) in place of setmtcschoicemodel!(g, ...)", :setmctschoicemodel!)
	setchoicemodel!(g, MCTSChoiceModel(g, rewardfunction, visitbudget, samplesize, cp))
end

# reset any state
function resetstate!(cm::MCTSChoiceModel)
	resetstate!(cm.policychoicemodel)
	cm.rootgodelsequence = (Real)[]
	cm.rootnode = MCTSNode(Dict())
end

#DEBUG function godelsequenceasstring(godelsequence::Vector{Real})
#DEBUG 	s = ""
#DEBUG 	for godelnumber in godelsequence
#DEBUG 		s = s * "$(godelnumber) "
#DEBUG 	end
#DEBUG 	s
#DEBUG end

#DEBUG function printmodel(cm::MCTSChoiceModel)
	#DEBUG println("- START model")
	#DEBUG println("rootgodelsequence: $(godelsequenceasstring(cm.rootgodelsequence))")
#DEBUG 	printnode(cm.rootnode)
	#DEBUG println("- END model")
#DEBUG 	readline()
#DEBUG end

# note: returns the godelnumber (i.e. action) of the best child node is returned, rather than the node itseld
function bestchildgodelnumber(node::MCTSNode, cp::Real)

	#DEBUG println("- START bestchildgodelnumber")
	#DEBUG println("using cp: $(cp)")

	@assert length(node.childnodes) >= 1

	bestucb1 = -Inf
	bestgodelnumbers = (Real)[]
	for childgodelnumber in keys(node.childnodes)
		childnode = node.childnodes[childgodelnumber]
		@assert childnode.visitcount >= 1
		ucb1 = (childnode.totalreward / childnode.visitcount) + cp * sqrt(2 * Base.log(node.visitcount) / childnode.visitcount)
		if ucb1 > bestucb1
			bestucb1 = ucb1
			bestgodelnumbers = [childgodelnumber]
		elseif ucb1 == bestucb1
			push!(bestgodelnumbers, childgodelnumber)
		end
	end

	# return one of the 'best' godelnumbers
	@assert length(bestgodelnumbers) >= 1
	best = bestgodelnumbers[rand(1:length(bestgodelnumbers))]

	#DEBUG println("choose: $(best)")
	#DEBUG println("- END bestchildgodelnumber")
	best

end

function godelnumber(cm::MCTSChoiceModel, cc::ChoiceContext)

	#DEBUG println("- START godelnumber")

	#DEBUG println("initial model:")
	#DEBUG printmodel(cm)

	@assert isempty(cm.rootgodelsequence) || (cm.rootgodelsequence == cc.derivationstate.godelsequence)
	cm.rootgodelsequence = deepcopy(cc.derivationstate.godelsequence)

	# keep on expanind and running simulations until the rootnode has been visited the required number of times

	while cm.rootnode.visitcount < cm.visitbudget

		currentnode = cm.rootnode
		currentgodelsequence = deepcopy(cm.rootgodelsequence)
		nodestoreward = [currentnode]
		reward = 0.0

		# when current node can't be expanded, move to best child, until we find an expandable (or terminal) node
		while (currentnode.samplecount == cm.samplesize) && !isempty(currentnode.childnodes)
			#DEBUG println("- START expandable node iteration")
			childgodelnumber = bestchildgodelnumber(currentnode, cm.cp)
			childnode = currentnode.childnodes[childgodelnumber]
			push!(nodestoreward, childnode)
			currentnode = childnode
			push!(currentgodelsequence, childgodelnumber)
			#DEBUG println("new godel sequence: $(godelsequenceasstring(currentgodelsequence))")
			#DEBUG println("- END expandable node iteration")
		end

		# now at a child node that can be expanded
		# perform a simulation that:
		#  (1) expand current node according to probability distribution of policy choice model
		#      - this is because we don't know set of possible actions from current state (and
		#      indeed this set may be uncountable), so we sample in a similar way to the NMCS
		#      implementation
		#  (2) apply the default policy to the end state

		#DEBUG println("- START simulation")

		policychoicemodel = deepcopy(cm.policychoicemodel)
		simulationgenerator = deepcopy(cc.derivationstate.generator)
		presetgodelsequence = deepcopy(currentgodelsequence)
		simulationcm = MCTSSimulationChoiceModel(policychoicemodel, presetgodelsequence)
		setchoicemodel!(simulationgeneration, simulationcm)
		result, state = nothing, nothing

		#DEBUG println("simulation using sequence: $(godelsequenceasstring(presetgodelsequence))")

		try
			result, state = generate(simulationgenerator; maxchoices=cc.derivationstate.maxchoices)
			reward = cm.rewardfunction(result)
		catch e
  			if !isa(e, GenerationTerminatedException)
				throw(e)
			end
			# TODO could punish further, e.g. negative reward?
		end

		#DEBUG println("simulation result: $(result)")
		#DEBUG println("simulation reward: $(reward)")
		#DEBUG println("- END simulation")

		if state != nothing

			#DEBUG println("- START applying simulation results")

			#DEBUG println("simulation returned sequence: $(godelsequenceasstring(state.godelsequence))")

			currentnode.samplecount += 1

			@assert length(state.godelsequence) >= length(currentgodelsequence)
			@assert state.godelsequence[1:length(currentgodelsequence)] == currentgodelsequence

			# check that a new godel number was returned (i.e. currentnode is not a terminal state)
			if length(state.godelsequence) > length(currentgodelsequence)

				# check which action (godelnumber) was taken next in the simulation
				childgodelnumber = state.godelsequence[length(currentgodelsequence)+1]

				#DEBUG println("first godel number in simulation: $(childgodelnumber)")

				# check to see if node for this action has already been created
				if haskey(currentnode.childnodes, childgodelnumber)
					#DEBUG println("first godel number is for existing child node")
					childnode = currentnode.childnodes[childgodelnumber]
				else
					#DEBUG println("first godel number is for new child node")
					# get choice point trace info for the child godel number: this will be the trace
					# info returned should this godel number be selected
					cptrace = state.cmtrace[length(currentgodelsequence)+1]
					childnode = MCTSNode(cptrace[2]) # cptrace is tuple of cp id and trace info for the returned gn, so need the second element of the tuple
					currentnode.childnodes[childgodelnumber] = childnode
				end

				push!(nodestoreward, childnode)

				currentnode = childnode # not necessary?
				push!(currentgodelsequence, childgodelnumber) # not necessary?

			end

			# now backpropogate reward
			for nodetoreward in nodestoreward
				nodetoreward.totalreward += reward
				nodetoreward.visitcount += 1
				if reward > nodetoreward.maxreward
					nodetoreward.maxreward = reward
				end
			end

			#DEBUG println("model after applying simulation results:")
			#DEBUG printmodel(cm)

			#DEBUG println("- END applying simulation results")

		end

		#DEBUG println("- END visit iteration")

	end

	# Taking desired godelnumber as that of the best child of the current root
	# is an ALTERNATIVE to the code below that finds max reward (and may be better for
	# the "single-player" game used here as we can ensure we makes moves to get this
	# max reward rather than the average reward)
	# gn = bestchildgodelnumber(cm.rootnode, 0.0) # NB with C_p of 0

	# desired godelnumber is that of the child with the largest maxreward
	bestmaxreward = -Inf
	bestgodelnumbers = (Real)[]
	for childgodelnumber in keys(cm.rootnode.childnodes)
		childnode = cm.rootnode.childnodes[childgodelnumber]
		if childnode.maxreward > bestmaxreward
			bestmaxreward = childnode.maxreward
			bestgodelnumbers = [childgodelnumber]
		elseif childnode.maxreward == bestmaxreward
			push!(bestgodelnumbers, childgodelnumber)
		end
	end

	# return one of the 'best' godelnumbers
	@assert length(bestgodelnumbers) >= 1
	gn = bestgodelnumbers[rand(1:length(bestgodelnumbers))]
	trace = cm.rootnode.childnodes[gn].trace

	# the child with this godelnumber now becomes the root
	cm.rootnode = cm.rootnode.childnodes[gn]
	push!(cm.rootgodelsequence, gn) # not strictly necessary since set from cc on entry (but needed now for assert)

	#DEBUG println("bounds: $(cc.lowerbound) $(cc.upperbound)")
	#DEBUG println("returning godelnumber: $(gn)")

	gn, trace

end

setparams!(cm::MCTSChoiceModel, params) = setparams!(cm.policychoicemodel, params)
getparams(cm::MCTSChoiceModel) = getparams(cm.policychoicemodel)
paramranges(cm::MCTSChoiceModel) = paramranges(cm.policychoicemodel)

type MCTSSimulationChoiceModel <: ChoiceModel
	policychoicemodel::ChoiceModel
	presetgodelsequence::Vector{Real}
end

function godelnumber(cm::MCTSSimulationChoiceModel, cc::ChoiceContext)
	if isempty(cm.presetgodelsequence)
		gn, trace = godelnumber(cm.policychoicemodel, cc)
	else
		gn, trace = shift!(cm.presetgodelsequence), Dict()
	end
	gn, trace
end

# reset any state
function resetstate!(cm::MCTSSimulationChoiceModel)
	resetstate!(cm.policychoicemodel)
end

show(io::IO, cm::MCTSChoiceModel) = print(io, "MCTS choice model (policy: ", cm.policychoicemodel, ")")


text:
# RwaJur4

export VERTICES_RWAJUR4,
		RwaJur4Exposure,
		RwaJur4Parameters,
		getLongExposureArray,
		getShortExposureArray,
		#getRwaJur4NetExposure,
		#getRwaJur4VerticalGap,
		#getRwaJur4Zones,
		#getRwaJur4HorizontalGapWithinZones,
		#getRwaJur4HorizontalGapBetweenZones,
		getPJur4,
		getRwaJur4

"""
Standardized vertices for mapping interest rate coupon rate exposures (RWAJUR4)
"""
const VERTICES_RWAJUR4 = VERTICES_RWAJUR2

"""
    RwaJur4Exposure

Holds exposures in interest rate coupon rates.

# Fields
* `date::Date`: date of evaluation
* `VXX::Float64`: exposure to vertex XX,
given by the sum of present value in BRL
of all cashflows that matures in XX business days (BD).

# Constructors
    RwaJur4Exposure(date::Date)

Create an instance with all exposures equal to zero.

	RwaJur4Exposure(date::Date, currency::AbstractString, longExposures::Array{Float64, 1}, shortExposures::Array{Float64, 1})

Create an instance with exposures given by two vectors of length 10 containing long and short exposures.
Exposures are assigned in sequence from vertex 1 BD to vertex 2520 BD.

"""
type RwaJur4Exposure
	# Date of evaluation
	date::Date
	# Currency
	currency::AbstractString
	# Exposures in standard vertices: VXX, where XX is a maturity calculated in business days
	# LONG EXPOSURES
	V1Long::Float64
	V21Long::Float64
	V42Long::Float64
	V63Long::Float64
	V126Long::Float64
	V252Long::Float64
	V504Long::Float64
	V756Long::Float64
	V1008Long::Float64
	V1260Long::Float64
	V2520Long::Float64
	# SHORT EXPOSURES
	V1Short::Float64
	V21Short::Float64
	V42Short::Float64
	V63Short::Float64
	V126Short::Float64
	V252Short::Float64
	V504Short::Float64
	V756Short::Float64
	V1008Short::Float64
	V1260Short::Float64
	V2520Short::Float64

	# Constructors
	RwaJur4Exposure(date::Date, currency::AbstractString) = new(date, currency, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)

	function RwaJur4Exposure(date::Date, currency::AbstractString, longExposures::Array{Float64, 1}, shortExposures::Array{Float64, 1})
		if size(longExposures, 1) != 11
			error("long exposures must have 11 values")
		end #if
		if any(longExposures .< 0)
			error("all long exposures must be positive")
		end #if

		if size(shortExposures, 1) != 11
			error("short exposures must have 11 values")
		end #if
		if any(shortExposures .< 0)
			error("all short exposures must be positive")
		end #if

		new(date,
			currency,
		    longExposures[1], longExposures[2], longExposures[3], longExposures[4],
			longExposures[5], longExposures[6], longExposures[7], longExposures[8],
			longExposures[9], longExposures[10], longExposures[11],
			shortExposures[1], shortExposures[2], shortExposures[3], shortExposures[4],
			shortExposures[5], shortExposures[6], shortExposures[7], shortExposures[8],
			shortExposures[9], shortExposures[10], shortExposures[11]
			)
	end #function

end #type

"""
	RWAJUR4Parameters

Holds BCB parameters to calculate interest rate coupon stantardized capital (RWA JUR4)

"""
type RwaJur4Parameters
	date::Date
	M4::Float64
end #type

"""
    getLongExposureArray(rExp::RwaJur4Exposure)

Return a vector with all long exposures in sequence,
from vertex 1 BD to vertex 2520 BD.

"""
getLongExposureArray(rExp::RwaJur4Exposure) = [rExp.V1Long, rExp.V21Long, rExp.V42Long, rExp.V63Long, rExp.V126Long, rExp.V252Long, rExp.V504Long, rExp.V756Long, rExp.V1008Long, rExp.V1260Long, rExp.V2520Long]

"""
    getShortExposureArray(rExp::RwaJur4Exposure)

Return a vector with all short exposures in sequence,
from vertex 1 BD to vertex 2520 BD.

"""
getShortExposureArray(rExp::RwaJur4Exposure) = [rExp.V1Short, rExp.V21Short, rExp.V42Short, rExp.V63Short, rExp.V126Short, rExp.V252Short, rExp.V504Short, rExp.V756Short, rExp.V1008Short, rExp.V1260Short, rExp.V2520Short]

# Casts
"""
    convert(::Type{RwaJur2Exposure}, x::RwaJur4Exposure)

Casts an RwaJur4Exposure into an RwaJur2Exposure.
"""
convert(::Type{RwaJur2Exposure}, x::RwaJur4Exposure) = RwaJur2Exposure(x.date, x.currency, getLongExposureArray(x), getShortExposureArray(x))

"""
    convert(::Type{RwaJur2Parameters}, x::RwaJur4Parameters)

Casts an RwaJur4Parameters into an RwaJur2Parameters.
"""
convert(::Type{RwaJur2Parameters}, x::RwaJur4Parameters) = RwaJur2Parameters(x.date, x.M4)

"""
    getRwaJur4NetExposure(rExp::RwaJur4Exposure)

Return a vector containing RWA JUR4 weighted net exposures.

"""
getRwaJur4NetExposure(rExp::RwaJur4Exposure) = getRwaJur2NetExposure(convert(RwaJur2Exposure, rExp))

"""
    getRwaJur4VerticalGap(rExp::RwaJur4Exposure)

Return a vector containing RWA JUR4 vertical gaps.

"""
getRwaJur4VerticalGap(rExp::RwaJur4Exposure) = getRwaJur2VerticalGap(convert(RwaJur2Exposure, rExp))

"""
    getRwaJur4Zones(rExp::RwaJur4Exposure)

Return a matrix containing RWA JUR4 zone exposures.
Column 1 contains long exposures, column2 contains short exposures.
Zones are represented in each of the 3 rows.

"""
getRwaJur4Zones(rExp::RwaJur4Exposure) = getRwaJur2Zones(convert(RwaJur2Exposure, rExp))

"""
    getRwaJur3HorizontalGapWithinZones(rExp::RwaJur4Exposure)

Return a vector containing RWA JUR4 horizontal gap within zones.

"""
getRwaJur4HorizontalGapWithinZones(rExp::RwaJur4Exposure) = getRwaJur2HorizontalGapWithinZones(convert(RwaJur2Exposure, rExp))

"""
    getRwaJur4HorizontalGapBetweenZones(rExp::RwaJur4Exposure)

Return a vector containing RWA JUR4 horizontal gap between zones.

"""
getRwaJur4HorizontalGapBetweenZones(rExp::RwaJur4Exposure) = getRwaJur2HorizontalGapBetweenZones(convert(RwaJur2Exposure, rExp))

"""
	getPJur3(rExp::RwaJur4Exposure, par::RwaJur4Parameters)

Calculate PJUR4 using maturity ladder methodology.

"""
getPJur4(rExp::RwaJur4Exposure, par::RwaJur4Parameters) = getPJur2(convert(RwaJur2Exposure, rExp), convert(RwaJur2Parameters, par))

"""
	getRwaJur4(rExp::RwaJur4Exposure, par::RwaJur4Parameters)

Calculate risk-weighted assets for interest rate coupon exposures (RWA JUR4),
using maturity ladder methodology.

"""
getRwaJur4(rExp::RwaJur4Exposure, par::RwaJur4Parameters) = getRwaJur2(convert(RwaJur2Exposure, rExp), convert(RwaJur2Parameters, par))

ext: jl


text:
module LaTeX

using SHA, Compat, Pkg, Dates, Requires
import Images

Sys.iswindows() && include("wincall.jl")

export Section, Table, Tabular, Figure, Image, ImageFileData, Code, TOC,
    Abstract, report, makepdf, writepdf, openpdf, document, DocumentClass, Title, Author

mutable struct TOC
end

mutable struct Abstract
    content
end

mutable struct Section
    title
    content
end

mutable struct Table
    caption
    content
end

mutable struct Tabular
    content::Array
end

mutable struct Figure
    caption
    content
end

mutable struct ImageFileData
    data::Array{UInt8}
    typ::Symbol
end

mutable struct Image
    height
    width
    data::ImageFileData
end

Image(height, width, data::Array) = Image(height, width, size(data,3) == 3 ? Images.colorim(data) : Images.grayim(data'))
function Image(height, width, image)
    filename = tempname()*".png"
    Images.save(filename, image)
    r = read(filename)
    rm(filename)
    Image(height, width, ImageFileData(r, :png))
end

mutable struct Code
    code
end

# declarations: non-displayed items controlling document metadata
abstract type AbstractDecl end

mutable struct DocumentClass <: AbstractDecl
    class::AbstractString
    settings::Vector{AbstractString}
end

mutable struct Author <: AbstractDecl
    text::AbstractString
end

mutable struct Title <: AbstractDecl
    text::AbstractString
end

mutable struct Style <: AbstractDecl
    section::AbstractString
end

isdecl(::AbstractDecl) = true
isdecl(::Date) = true
isdecl(_) = false

function processdecl(d::DocumentClass)
    if isempty(d.settings)
        "\\documentclass{$(d.class)}"
    else
        "\\documentclass[$(join(d.settings, ','))]{$(d.class)}"
    end
end
processdecl(a::Author) = "\\author{$(a.text)}"
processdecl(t::Title) = "\\title{$(t.text)}"
processdecl(d::Date) = "\\date{$d}"
processdecl(s::Style) = "\\allsectionsfont{$(s.section)}"

"""
    makepdf(latex)

Build pdf document in temporary directory using pdflatex and returns
its path.
"""
function makepdf(latex)
    dirname = "$(tempname()).d"
    mkdir(dirname)
    texname = joinpath(dirname, "document.tex")
    pdfname = joinpath(dirname, "document.pdf")
    open(texname, "w") do file
        write(file, latex)
    end
    cd(dirname) do
        for i in 1:2
            output = read(`pdflatex -shell-escape -halt-on-error $texname`,String)
            occursin("Error:", output) && println(output)
        end
    end
    pdfname
end

"""
    writepdf(latex, filename)

Build pdf document and copy it to `filename`.
"""
function writepdf(latex, filename)
    pdfname = makepdf(latex)
    cp(pdfname,filename,remove_destination=false)
    nothing
end

"""
    openpdf(latex)

Build pdf document and open it.
"""
function openpdf(latex)
    pdfname = makepdf(latex)
    if Sys.iswindows()
        command = `cmd /K start \"\" $pdfname`
        run(command)
    else
        run(`open $pdfname`)
    end
    pdfname
end

flatten(a) = flatten(Any[],a)
function flatten(r, a)
    for x in a
        isa(x,Array) ? flatten(r,x) : push!(r,x)
    end
    r
end

processitem(p, item::T, indent) where {T<:AbstractString} = [item]
processitem(p, item::T, indent) where {T<:Number} = ["$item"]
processitem(p, decl, indent) = []

function processitem(p, items::Array, indent)
    isempty(items) && return [""]
    map(x -> processitem(p, x, indent), items)
end

function processitem(p, item::TOC, indent)
    ["\\tableofcontents"]
end

function processitem(p, item::Abstract, indent)
    r = Any["\\begin{abstract}"]
    append!(r, processitem(p, item.content, indent))
    push!(r, "\\end{abstract}")
end

function processitem(p, item::Section, indent)
    if indent > p[:maxdepth]
        cmd = last(p[:sectioncommands])
    else
        cmd = p[:sectioncommands][indent]
    end

    r = Any["\\$cmd{$(item.title)}\\nopagebreak"]
    append!(r, processitem(p, item.content, indent+1))
end

function processitem(p, item::Figure, indent)
    r = Any["\\begin{figure}[!ht]"]
    append!(r, processitem(p, item.content, indent))
    append!(r, ["\\caption{$(item.caption)}", "\\end{figure}"])
end

processitem(p, item::Code, indent) = [
    "\\usestyle{default}",
    "\\begin{pygmented}{jl}",
    split(item.code, '\n')...,
    "\\end{pygmented}"]

function processitem(p, item::Image, indent)
    filename = joinpath(tempdir(), bytes2hex(sha256(item.data.data))*".$(item.data.typ)")
    open(filename, "w") do file
        write(file, item.data.data)
    end

    r = Any["\\includegraphics["]
    if !isempty(item.width)
        if item.width <= 1
            push!(r, "width=$(item.width)\\textwidth,")
        else
            push!(r, "width=$(item.width)cm,")
        end
    end
    if !isempty(item.height)
        if item.height <= 1
            push!(r, "height=$(item.height)\\textheight,")
        else
            push!(r, "height=$(item.height)cm")
        end
    end
	escaped_filename = replace(filename,"\\" => "/")
    push!(r, "]{$escaped_filename}")
    flatten(r)
end

function processitem(p, item::Table, indent)
    r = Any["\\begin{table}[!ht]"]
    push!(r, processitem(p, item.content, indent))
    push!(r, ["\\caption{$(item.caption)}", "\\end{table}"])
end

function processitem(p, item::Tabular, indent)
    if ndims(item.content) == 1
        item.content = reshape(item.content, (1, length(item.content)))
    end
    sm, sn = size(item.content)
    r = Any["\\begin{tabular}[!ht]{$(repeat("c", sn))}"]
    for m = 1:sm
        for n = 1:sn
            push!(r, processitem(p, item.content[m,n], indent))
            if n < sn
                push!(r, " & ")
            end
        end
        push!(r, " \\\\")
    end
    push!(r, "\\end{tabular}")
end

# get required packages and settings for this element and all children
getrequirements(item) = Dict()
function getrequirements(items::Array)
    base = Dict()
    for i in items
        requirements = getrequirements(i)
        for (requirement, settings) in requirements
            if haskey(base, requirement)
                union!(base[requirement], settings)
            else
                base[requirement] = settings
            end
        end
    end
    base
end
getrequirements(::Code) = Dict("texments" => Set([]))
getrequirements(::Image) = Dict("graphicx" => Set([]))
getrequirements(a::Abstract) = getrequirements(a.content)
getrequirements(s::Section) = getrequirements(Any[s.title, s.content])
getrequirements(t::Table) = getrequirements(Any[t.caption, t.content])
getrequirements(t::Tabular) = getrequirements(t.content)
getrequirements(f::Figure) = getrequirements(Any[f.caption, f.content])

getrequirements(d::Style) = Dict("sectsty" => Set([]))

# add information about document class in dictionary.
function inform!(p, d::DocumentClass)
    if d.class == "report"
        p[:sectioncommands] = ["chapter", "section", "subsection",
            "subsubsection","paragraph"]
    else  # article & others
        p[:sectioncommands] = ["section", "subsection", "subsubsection",
            "paragraph"]
    end
end

document(items) = document(Dict(), items)
function document(p, items)
    # make required path
    p = merge((Dict(:maxdepth => 3, :tmppath => tempdir())), p)

    mkpath(p[:tmppath])

    # document properties
    dclass = nothing
    decls = Dict()

    # parse items tree for declarations (currently to depth 1)
    for tr in items
        if isa(tr, DocumentClass)
            dclass = tr
        elseif isdecl(tr)
            decls[typeof(tr)] = tr
        end
    end

    # modify preamble dict to reflect docclass, e.g. section headers
    inform!(p, dclass)

    # global required packages (keep small)
    require = Dict(
        "inputenc" => Set(["latin1"]),
        "fullpage" => Set(["cm"]),
        "morefloats" => Set([]),
        "placeins" => Set(["section"]))

    # add required packages (dynamic based off items)
    docrequire = getrequirements(items)

    preamble = AbstractString[]
    docbody = processitem(p, items, 1)

    # build up the document starting from the beginning
    push!(preamble, processdecl(dclass))
    for (package, settings) in merge(require, docrequire)
        if isempty(settings)
            push!(preamble, "\\usepackage{$package}")
        else
            push!(preamble, "\\usepackage[$(join(settings, ','))]{$package}")
        end
    end

    # some other default stuff
    append!(preamble, [
        [processdecl(dl) for dl in values(decls)]...,
        "\\begin{document}",  # technically not preamble anymore
        "\\maketitle"])

    # append the document
    r = vcat(preamble, docbody)
    push!(r, "\\end{document}")

    r = join(flatten(r), "\n")
    r
end

report(items; kargs...) = report(Dict(), items; kargs...)
function report(p, items; author = "", title = "Report", date = "", toc = false,
    theabstract = "")

    doctree = Any[DocumentClass("report", ["11pt", "a4paper"])]
    # make date, author, etc. if available
    isempty(date) || push!(doctree, Date(date))
    isempty(author) || push!(doctree, Author(author))
    push!(doctree, Title(title))
    isempty(theabstract) || push!(doctree, Abstract(theabstract))
    toc && push!(doctree, TOC())

    # default styling
    push!(doctree, Style("\\normalfont\\sffamily\\bfseries"))

    if isa(items, Vector)
        append!(doctree, items)
    else
        push!(doctree, items)
    end
    document(p, doctree)
end

function __init__()
    @require Winston = "bd07be1c-e76f-5ff0-9c0b-f51ef45303c6" include("winston.jl")
    @require Gadfly = "c91e804a-d5a3-530f-b6f0-dfbca275c004" include("gadfly.jl")
    @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("plots.jl")
end

#openpdf(report(Dict(),[]))
#openpdf(report(Dict(),["Test"]))
#openpdf(report(Section("Section header", "content"), toc = false))

#sec1 = Section("sec1", "content1")
#sec2 = Section("sec", "content1")
#ch1 = Section("ch1", [sec1, sec2])

#openpdf(report(ch1, toc=true))

#openpdf(report(Figure("caption", "content")))
#openpdf(report(Table("caption", "content")))
#openpdf(report(Tabular(['1' '2'])))

end # module


text:
using Primes
using Printf

"""
Stores global and per-rank information about a grid imposed over a simulation
domain.
"""
abstract type AbstractGrid end

"""
    UniformGrid(lower_dims, upper_dims, num_cells, num_guard_cells, periodic)

Stores global and per-rank information about a uniform `N`-dimensional spatial
grid.
"""
struct UniformGrid{N, FT, IT} <: AbstractGrid
    global_num_cells::NTuple{N, IT}
    num_cells::NTuple{N, IT}
    lower_cell_bound::NTuple{N, IT}
    upper_cell_bound::NTuple{N, IT}
    num_guard_cells::IT

    global_lower_dimensions::NTuple{N, FT}
    lower_dimensions::NTuple{N, FT}
    global_upper_dimensions::NTuple{N, FT}
    upper_dimensions::NTuple{N, FT}
    global_simulation_size::NTuple{N, FT}
    simulation_size::NTuple{N, FT}
    cell_size::NTuple{N, FT}

    periodic::NTuple{N, Bool}
    proc_dims::NTuple{N, IT}
end

function UniformGrid(global_lower_dimensions::NTuple{N, FT},
              global_upper_dimensions::NTuple{N, FT},
              global_num_cells::NTuple{N, IT},
              num_guard_cells::IT=zero(IT),
              periodic::NTuple{N, Bool}=ntuple(_ -> true, N),
             ) where {N, FT, IT}

    proc_dims = Tuple(compute_processor_dimensions(global_num_cells,
                                                MPI.Comm_size(MPI.COMM_WORLD)))

    lower_cell_bound, upper_cell_bound =
        decompose_domain(global_num_cells,
                         MPI.Comm_size(MPI.COMM_WORLD),
                         MPI.Comm_rank(MPI.COMM_WORLD))
    num_cells = upper_cell_bound .- lower_cell_bound .+ 1

    global_simulation_size  = global_upper_dimensions .- global_lower_dimensions
    cell_size = global_simulation_size ./ global_num_cells

    lower_dimensions = (lower_cell_bound .- 1) .* cell_size
    upper_dimensions = upper_cell_bound .* cell_size
    simulation_size  = upper_dimensions .- lower_dimensions

    return UniformGrid{N, FT, IT}(global_num_cells,
                           num_cells,
                           lower_cell_bound,
                           upper_cell_bound,
                           num_guard_cells,

                           global_lower_dimensions,
                           lower_dimensions,
                           global_upper_dimensions,
                           upper_dimensions,
                           global_simulation_size,
                           simulation_size,
                           cell_size,

                           periodic,
                           proc_dims
                          )
end

function Base.show(io::IO, g::UniformGrid{N, FT, IT}) where {N, FT, IT}
    print(io, "UniformGrid{$N, $FT, $IT}(global=(")
    for i in 1:N
        str = @sprintf("%.4g:%.4g:%.4g", g.global_lower_dimensions[i], g.cell_size[i],
                       g.global_upper_dimensions[i])
        print(io, str)
        i != N && print(io, ", ")
    end
    print(io, "), local=(")
    for i in 1:N
        str = @sprintf("%.4g:%.4g:%.4g", g.lower_dimensions[i], g.cell_size[i],
                       g.upper_dimensions[i])
        print(io, str)
        i != N && print(io, ", ")
    end
    print(io, ")), $(g.lower_cell_bound), $(g.upper_cell_bound)")
end

"""
    compute_processor_dimensions(num_cells::NTuple{N}, num_procs)

Given the desired number of cells in each dimension, and a number of processors,
returns the number of subdivisions along each dimension that minimize the amount
of inter-rank communication.
"""
function compute_processor_dimensions(num_cells::NTuple{N}, num_procs) where N
    factors = factor(Vector, num_procs)

    proc_dims = fill(1, N)

    # For each prime factor of num_procs, we cut the simulation domain in the
    # dimension with the most cells remaining.
    for factor in Iterators.reverse(factors)
        # Find the current largest dimension
        largest_value = num_cells[1] / proc_dims[1]
        largest_index = 1
        for j in 2:N
            if num_cells[j] / proc_dims[j] > largest_value
                largest_value = num_cells[j] / proc_dims[j]
                largest_index = j
            end
        end

        # Cut along that dimension to make it smaller
        proc_dims[largest_index] *= factor
    end

    # We should have make exactly num_procs subdomains
    @assert reduce(*, proc_dims) == num_procs

    return proc_dims
end

"""
    compute_proc_indices(proc_dims, proc)

Computes the unique indices in the processor grid for `proc`.
"""
function compute_proc_indices(proc_dims, proc)
    indices = Vector{Int}(undef, length(proc_dims))
    factor = 1
    sum = 0
    for i in 1:length(proc_dims)
        indices[i] = div(mod(proc - sum, factor * proc_dims[i]), factor)
        factor *= proc_dims[i]
        sum += indices[i]
    end
    return Tuple(indices)
end

"""
    compute_proc_rank(proc_dims, indices)

Computes the rank associated with the processor at `indices`.
"""
function compute_proc_rank(proc_dims, indices)
    rank = 0
    prod = 1
    for i in 1:length(indices)
        rank += prod * mod(indices[i], proc_dims[i])
        prod *= proc_dims[i]
    end
    return rank
end

"""
    compute_proc_bounds(num_cells, proc_dims, proc)

Computes the range of cells that `proc` is responsible for.
"""
function compute_proc_bounds(num_cells, proc_dims, proc)
    indices = compute_proc_indices(proc_dims, proc)

    num_subcells = Tuple(div.(num_cells, proc_dims))
    extra_cells  = Tuple(mod.(num_cells, proc_dims))

    lower_bound = num_subcells .* indices .+ 1 .+
        ifelse.(indices .< extra_cells, indices, extra_cells)
    upper_bound = lower_bound .+ num_subcells .- 1 .+
        ifelse.(indices .< extra_cells, 1, 0)

    return lower_bound, upper_bound
end

"""
    decompose_domain(num_cells::NTuple{N}, num_procs[, proc])

Given a target number of pieces `num_procs` and an `N`-cube simulation domain
with `num_cells` along each dimension, computes a decomposition of the domain.
If `proc` is specified, returns a pair of Tuples corresponding to the lower and
upper bounds of the cells assigned to `proc`. Otherwise, returns a pair of
vectors of the lower and upper bounds for each processor.

NOTE: `proc` is zero indexed.
"""
function decompose_domain(num_cells::NTuple{N}, num_procs) where N
    proc_dims = compute_processor_dimensions(num_cells, num_procs)

    lower_bounds = Vector{typeof(num_cells)}(undef, num_procs)
    upper_bounds = Vector{typeof(num_cells)}(undef, num_procs)

    for i in 1:num_procs
        lower_bounds[i], upper_bounds[i] = compute_proc_bounds(num_cells, proc_dims, i - 1)
    end

    return lower_bounds, upper_bounds
end

function decompose_domain(num_cells::NTuple{N}, num_procs, proc) where N
    proc_dims = compute_processor_dimensions(num_cells, num_procs)

    return compute_proc_bounds(num_cells, proc_dims, proc)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Testing

@testset "test initialize scalar field problem" begin
    el = Element(Seg2, [1, 2])
    pr = Problem(Heat, 1)
    push!(pr, el)
    initialize!(pr)
    @test haskey(el, "temperature")
    # one timestep in field "temperature"
    @test length(el["temperature"]) == 1
    # this way we access to field at default time t=0.0, it's different than ^!
    @test length(el("temperature")) == 2
    # length of single increment
    @test length(el("temperature", 0.0)) == 2
    @test length(last(el, "temperature").data) == 2
end

@testset "test initialize vector field problem" begin
    el = Element(Seg2, [1, 2])
    pr = Problem(Elasticity, 2)
    push!(pr, el)
    initialize!(pr)
    @test haskey(el, "displacement")
    @test length(el["displacement"]) == 1
    # this way we access to field at default time t=0.0, it's different than ^!
    @test length(el("displacement")) == 2
    # length of single increment
    @test length(el("displacement", 0.0)) == 2
    @test length(last(el, "displacement").data) == 2
end

@testset "test initialize boundary problem" begin
    el = Element(Seg2, [1, 2])
    pr = Problem(Dirichlet, "bc", 1, "temperature")
    push!(pr, el)
    initialize!(pr)
    @test haskey(el, "reaction force")
    @test haskey(el, "temperature")
end

#=
@testset "dict field depending from problems" begin
    p1 = Problem(Elasticity, "Body 1", 2)
    p2 = Problem(Elasticity, "Body 2", 2)
    X = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0],
        3 => [1.0, 1.0],
        4 => [0.0, 1.0])
    update!([p1, p2], "geometry", 0.0 => X)
    @test isapprox(p1("geometry", 0.0)[1], [0.0, 0.0])
    @test isapprox(p2("geometry", 0.0)[1], [0.0, 0.0])
    p1("geometry", 0.0)[1] = [1.0, 2.0]
    @test isapprox(p2("geometry", 0.0)[1], [1.0, 2.0])
end
=#

@testset "dict field depending from problems" begin
    X = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0],
        3 => [1.0, 1.0],
        4 => [0.0, 1.0],
        5 => [0.0, 2.0],
        6 => [1.0, 2.0],
        7 => [1.0, 3.0],
        8 => [0.0, 3.0])
    p1 = Problem(Elasticity, "Body 1", 2)
    p2 = Problem(Elasticity, "Body 2", 2)
    e1 = Element(Quad4, 1, [1, 2, 3, 4])
    e2 = Element(Quad4, 2, [5, 6, 7, 8])
    push!(p1, e1)
    push!(p2, e2)
    update!(p1, "geometry", 0.0 => X)
    update!(p2, "geometry", 0.0 => X)
    @test isapprox(p1("geometry", 0.0)[1], [0.0, 0.0])
    @test isapprox(p2("geometry", 0.0)[1], [0.0, 0.0])
    p1("geometry", 0.0)[1] = [1.0, 2.0]
    @test isapprox(p2("geometry", 0.0)[1], [1.0, 2.0])
    @test isapprox(e1("geometry", 0.0)[1], [1.0, 2.0])
end


text:
function resolve(path)
    if islink(path)
        dir = dirname(path)
        resolve(joinpath(dir, readlink(path)))
    else
        path
    end
end

# return a list of valid directories, resolving symlinks and pruning duplicates
function valid_dirs(dirs)
    map!(resolve, dirs, dirs)
    filter(isdir, unique(dirs))
end

## generic discovery routines

function library_versioned_names(name::String, version::Union{Nothing,VersionNumber,String}=nothing)
    names = String[]
    if Sys.iswindows()
        # Windows encodes the version in the filename
        if version isa VersionNumber
            append!(names, ["$(name)$(Sys.WORD_SIZE)_$(version.major)$(version.minor).$(Libdl.dlext)",
                            "$(name)$(Sys.WORD_SIZE)_$(version.major).$(Libdl.dlext)"])
        elseif version isa String
            push!(names, "$(name)$(Sys.WORD_SIZE)_$(version).$(Libdl.dlext)")
        elseif version === nothing
            push!(names, "$(name)$(Sys.WORD_SIZE).$(Libdl.dlext)")
        end

        # some libraries (e.g. CUTENSOR) are shipped without the word size-prefix
        if version isa VersionNumber
            append!(names, ["$(name)_$(version.major)$(version.minor).$(Libdl.dlext)",
                            "$(name)_$(version.major).$(Libdl.dlext)"])
        elseif version isa String
            push!(names, "$(name)_$(version).$(Libdl.dlext)")
        elseif version === nothing
            push!(names, "$(name).$(Libdl.dlext)")
        end
    elseif Sys.isapple()
        # macOS puts the version number before the dylib extension
        if version isa VersionNumber
            append!(names, ["lib$(name).$(version.major).$(version.minor).$(Libdl.dlext)",
                            "lib$(name).$(version.major).$(Libdl.dlext)"])
        elseif version isa String
            push!(names, "lib$(name).$(version).$(Libdl.dlext)")
        elseif version === nothing
            push!(names, "lib$(name).$(Libdl.dlext)")
        end
    elseif Sys.isunix()
        # most UNIX distributions ship versioned libraries (also see JuliaLang/julia#22828)
        if version isa VersionNumber
            append!(names, ["lib$(name).$(Libdl.dlext).$(version.major).$(version.minor).$(version.patch)",
                            "lib$(name).$(Libdl.dlext).$(version.major).$(version.minor)",
                            "lib$(name).$(Libdl.dlext).$(version.major)"])
        elseif version isa String
            push!(names, "lib$(name).$(Libdl.dlext).$(version)")
        elseif version === nothing
            push!(names, "lib$(name).$(Libdl.dlext)")
        end
    elseif version === nothing
        push!(names, "lib$name.$(Libdl.dlext)")
    end
    return names
end

"""
    find_library(name, version; locations=String[])

Wrapper for Libdl.find_library, performing a more exhaustive search:

- variants of the library name (including version numbers, platform-specific tags, etc);
- various subdirectories of the `locations` list, and finally system library directories.

Returns the full path to the library.
"""
function find_library(name::String, version::Union{Nothing,VersionNumber,String}=nothing;
                      locations::Vector{String}=String[])
    @debug "Request to look for library $name $version" locations

    # figure out names
    all_names = library_versioned_names(name, version)

    # figure out locations
    all_locations = String[]
    for location in locations
        push!(all_locations, location)
        push!(all_locations, joinpath(location, "lib"))
        if Sys.WORD_SIZE == 64
            push!(all_locations, joinpath(location, "lib64"))
            push!(all_locations, joinpath(location, "libx64"))
        end
        if Sys.iswindows()
            push!(all_locations, joinpath(location, "bin"))
            push!(all_locations, joinpath(location, "bin", Sys.WORD_SIZE==64 ? "x64" : "Win32"))
        end
    end

    @debug "Looking for library $(join(all_names, ", "))" locations=all_locations
    name_found = Libdl.find_library(all_names, all_locations)
    if isempty(name_found)
        return nothing
    end

    # find the full path of the library (which Libdl.find_library doesn't guarantee to return)
    path = Libdl.dlpath(name_found)
    @debug "Found library $(basename(path)) at $(dirname(path))"
    return path
end

"""
    find_binary(name; locations=String[])

Similar to `find_library`, performs an exhaustive search for a binary in various
subdirectories of `locations`, and finally PATH by using `Sys.which`.
"""
function find_binary(name::String; locations::Vector{String}=String[])
    @debug "Request to look for binary $name" locations

    # figure out locations
    all_locations = String[]
    for location in locations
        push!(all_locations, location)
        push!(all_locations, joinpath(location, "bin"))
    end
    # we look in PATH too by using `Sys.which` with unadorned names

    @debug "Looking for binary $name" locations=all_locations
    all_paths = [name; [joinpath(location, name) for location in all_locations]]
    for path in all_paths
        try
            program_path = Sys.which(path)
            if program_path !== nothing
                @debug "Found binary $path at $program_path"
                return program_path
            end
        catch
            # some system disallow `stat` on certain paths
        end
    end

    return nothing
end

## CUDA-specific discovery routines

const cuda_releases = [v"1.0", v"1.1",
                       v"2.0", v"2.1", v"2.2",
                       v"3.0", v"3.1", v"3.2",
                       v"4.0", v"4.1", v"4.2",
                       v"5.0", v"5.5",
                       v"6.0", v"6.5",
                       v"7.0", v"7.5",
                       v"8.0",
                       v"9.0", v"9.1", v"9.2",
                       v"10.0", v"10.1", v"10.2",
                       v"11.0", v"11.1", v"11.2", v"11.3", v"11.4"]

const cuda_library_versions = Dict(
    v"11.0.1" => Dict(
        # NOTE: encountered this version in a Docker container; not sure where it came from.
        "cupti"     => "2020.1.0", # wtf
        "nvtx"      => v"11.0.167",
        "cublas"    => v"11.0.0", #.191
        "cufft"     => v"10.1.3", #.191
        "curand"    => v"10.2.0", #.191
        "cusolver"  => v"10.4.0", #.191
        "cusparse"  => v"11.0.0", #.191
    ),
    v"11.0.2" => Dict(
        "cupti"     => "2020.1.0", # wtf
        "nvtx"      => v"11.0.167",
        "cublas"    => v"11.0.0", #.191
        "cufft"     => v"10.1.3", #.191
        "curand"    => v"10.2.0", #.191
        "cusolver"  => v"10.4.0", #.191
        "cusparse"  => v"11.0.0", #.191
    ),
    v"11.0.3" => Dict(
        "cupti"     => "2020.1.1", # docs mention 11.0.221
        "nvtx"      => v"11.0.167",
        "cublas"    => v"11.2.0", #.252
        "cufft"     => v"10.2.1", #.245
        "curand"    => v"10.2.1", #.245
        "cusolver"  => v"10.6.0", #.245
        "cusparse"  => v"11.1.1", #.245
    ),
    v"11.1.0" => Dict(
        "cupti"     => "2020.2.0", # docs mention 11.1.69
        "nvtx"      => v"11.1.74",
        "cublas"    => v"11.2.1", #.74
        "cufft"     => v"10.3.0", #.74
        "curand"    => v"10.2.2", #.74
        "cusolver"  => v"11.0.0", #.74
        "cusparse"  => v"11.2.0", #.275
    ),
    v"11.1.1" => Dict(
        "cupti"     => "2020.2.1", # docs mention 11.1.105
        "nvtx"      => v"11.1.74",
        "cublas"    => v"11.3.0", #.106
        "cufft"     => v"10.3.0", #.105
        "curand"    => v"10.2.2", #.105
        "cusolver"  => v"11.0.1", #.105
        "cusparse"  => v"11.3.0", #.10
    ),
    v"11.2.0" => Dict(
        "cupti"     => "2020.3.0", # docs mention 11.2.67
        "nvtx"      => v"11.2.67",
        "cublas"    => v"11.3.1", #.68
        "cufft"     => v"10.4.0", #.72
        "curand"    => v"10.2.3", #.68
        "cusolver"  => v"11.0.2", #.68
        "cusparse"  => v"11.3.1", #.68
    ),
    v"11.2.1" => Dict(
        "cupti"     => "2020.3.1", # docs mention 11.2.135
        "nvtx"      => v"11.2.67",
        "cublas"    => v"11.4.1", #.1026
        "cufft"     => v"10.4.0", #.135
        "curand"    => v"10.2.3", #.135
        "cusolver"  => v"11.1.0", #.135
        "cusparse"  => v"11.4.0", #.135
    ),
    v"11.2.2" => Dict(
        "cupti"     => "2020.3.1", # docs mention 11.2.152
        "nvtx"      => v"11.2.152",
        "cublas"    => v"11.4.1", #.1043
        "cufft"     => v"10.4.1", #.152
        "curand"    => v"10.2.3", #.152
        "cusolver"  => v"11.1.0", #.152
        "cusparse"  => v"11.4.1", #.1152
    ),
    v"11.3.0" => Dict(
        "cupti"     => "2021.1.0", # docs mention 11.3.58
        "nvtx"      => v"11.3.58",
        "cublas"    => v"11.4.2", #.10064
        "cufft"     => v"10.4.2", #.58
        "curand"    => v"10.2.4", #.58
        "cusolver"  => v"11.1.1", #.58
        "cusparse"  => v"11.5.0", #.58
    ),
    v"11.3.1" => Dict(
        "cupti"     => "2021.1.1", # docs mention 11.3.111
        "nvtx"      => v"11.3.109",
        "cublas"    => v"11.5.1", #.109
        "cufft"     => v"10.4.2", #.109
        "curand"    => v"10.2.4", #.109
        "cusolver"  => v"11.1.2", #.109
        "cusparse"  => v"11.6.0", #.109
    ),
    v"11.4.0" => Dict(
        "cupti"     => "2021.2.0", # docs mention 11.4.65
        "nvtx"      => v"11.4.43",
        "cublas"    => v"11.5.2", #.43
        "cufft"     => v"10.5.0", #.43
        "curand"    => v"10.2.5", #.43
        "cusolver"  => v"11.2.0", #.43
        "cusparse"  => v"11.6.0", #.43
    ),
)

function cuda_library_version(library, toolkit_version)
    if library == "nvtx"
        v"1"
    elseif toolkit_version >= v"11"
        # starting with CUDA 11, libraries are versioned independently
        if !haskey(cuda_library_versions, toolkit_version)
            error("CUDA.jl does not yet support CUDA $toolkit_version; please file an issue.")
        end
        # HACK: generalize this?
        if library == "cusolverMg"
            library = "cusolver"
        end
        if library == "cublasLt"
            library = "cublas"
        end
        cuda_library_versions[toolkit_version][library]
    else
        toolkit_version
    end
end

const cuda_library_names = Dict(
    "nvtx"      => "nvToolsExt"
)

# only for nvdisasm, to discover the CUDA toolkit version
const cuda_binary_versions = Dict(
    v"11.0.1" => Dict(
        # NOTE: encountered this version in a Docker container; not sure where it came from.
        "nvdisasm"  => v"11.0.167"
    ),
    v"11.0.2" => Dict(
        "nvdisasm"  => v"11.0.194"
    ),
    v"11.0.3" => Dict(
        "nvdisasm"  => v"11.0.221"
    ),
    v"11.1.0" => Dict(
        "nvdisasm"  => v"11.1.74"
    ),
    v"11.1.1" => Dict(
        "nvdisasm"  => v"11.1.74"   # ambiguous!
    ),
    v"11.2.0" => Dict(
        "nvdisasm"  => v"11.2.67"
    ),
    v"11.2.1" => Dict(
        "nvdisasm"  => v"11.2.135"
    ),
    v"11.2.2" => Dict(
        "nvdisasm"  => v"11.2.152"
    ),
    v"11.3.0" => Dict(
        "nvdisasm"  => v"11.3.58",
        "ptxas"     => v"11.3.58"
    ),
    v"11.3.1" => Dict(
        "nvdisasm"  => v"11.3.58",  # ambiguous!
        "ptxas"     => v"11.3.109"
    ),
    v"11.3.1" => Dict(
        "nvdisasm"  => v"11.4.43"
    ),
)

# simplified find_library/find_binary entry-points,
# looking up name aliases and known version numbers
# and passing the (optional) toolkit dirs as locations.
function find_cuda_library(library::String, toolkit_dirs::Vector{String},
                           toolkit_version::VersionNumber)
    toolkit_release = VersionNumber(toolkit_version.major, toolkit_version.minor)

    # figure out the location
    locations = toolkit_dirs
    ## CUPTI is in the "extras" directory of the toolkit
    if library == "cupti"
        toolkit_extras_dirs = filter(dir->isdir(joinpath(dir, "extras")), toolkit_dirs)
        cupti_dirs = map(dir->joinpath(dir, "extras", "CUPTI"), toolkit_extras_dirs)
        append!(locations, cupti_dirs)
    end
    ## NVTX is located in an entirely different location on Windows
    if library == "nvtx" && Sys.iswindows()
        if haskey(ENV, "NVTOOLSEXT_PATH")
            dir = ENV["NVTOOLSEXT_PATH"]
            @debug "Looking for NVTX library via environment variable" dir
        else
            program_files = ENV[Sys.WORD_SIZE == 64 ? "ProgramFiles" : "ProgramFiles(x86)"]
            dir = joinpath(program_files, "NVIDIA Corporation", "NvToolsExt")
            @debug "Looking for NVTX library in the default directory" dir
        end
        isdir(dir) && push!(locations, dir)
    end

    version = cuda_library_version(library, toolkit_version)
    name = get(cuda_library_names, library, library)
    find_library(name, version; locations=locations)
end
find_cuda_binary(name::String, toolkit_dirs::Vector{String}=String[]) =
    find_binary(name; locations=toolkit_dirs)

"""
    find_toolkit()::Vector{String}

Look for directories where (parts of) the CUDA toolkit might be installed. This returns a
(possibly empty) list of paths that can be used as an argument to other discovery functions.

The behavior of this function can be overridden by defining the `CUDA_PATH`, `CUDA_HOME` or
`CUDA_ROOT` environment variables, which should point to the root of the CUDA toolkit.
"""
function find_toolkit()
    dirs = String[]

    # look for environment variables to override discovery
    envvars = ["CUDA_PATH", "CUDA_HOME", "CUDA_ROOT"]
    filter!(var -> haskey(ENV, var) && ispath(ENV[var]), envvars)
    if !isempty(envvars)
        paths = unique(map(var->ENV[var], envvars))
        if length(paths) > 1
            @warn "Multiple CUDA environment variables set to different values: $(join(paths, ", "))"
        end

        @debug "Looking for CUDA toolkit via environment variables $(join(envvars, ", "))"
        append!(dirs, paths)
        return dirs
    end

    # look for the compiler binary (in the case PATH points to the installation)
    ptxas_path = find_binary("ptxas")
    if ptxas_path !== nothing
        ptxas_dir = dirname(ptxas_path)
        if occursin(r"^bin(32|64)?$", basename(ptxas_dir))
            ptxas_dir = dirname(ptxas_dir)
        end

        @debug "Looking for CUDA toolkit via ptxas binary" path=ptxas_path dir=ptxas_dir
        push!(dirs, ptxas_dir)
    end

    # look for the runtime library (in the case LD_LIBRARY_PATH points to the installation)
    libcudart_path = find_library("cudart")
    if libcudart_path !== nothing
        libcudart_dir = dirname(libcudart_path)
        if occursin(r"^(lib|bin)(32|64)?$", basename(libcudart_dir))
            libcudart_dir = dirname(libcudart_dir)
        end

        @debug "Looking for CUDA toolkit via CUDA runtime library" path=libcudart_path dir=libcudart_dir
        push!(dirs, libcudart_dir)
    end

    # look in default installation directories
    default_dirs = String[]
    if Sys.iswindows()
        # CUDA versions are installed in separate directories under a single base dir
        program_files = ENV[Sys.WORD_SIZE == 64 ? "ProgramFiles" : "ProgramFiles(x86)"]
        basedir = joinpath(program_files, "NVIDIA GPU Computing Toolkit", "CUDA")
        if isdir(basedir)
            entries = map(dir -> joinpath(basedir, dir), readdir(basedir))
            append!(default_dirs, entries)
        end
    else
        # CUDA versions are installed in unversioned dirs, or suffixed with the version
        basedirs = ["/usr/local/cuda", "/opt/cuda"]
        for ver in cuda_releases, dir in basedirs
            push!(default_dirs, "$dir-$(ver.major).$(ver.minor)")
        end
        append!(default_dirs, basedirs)
        push!(default_dirs, "/usr/lib/nvidia-cuda-toolkit")
        push!(default_dirs, "/usr/share/cuda")
    end
    reverse!(default_dirs) # we want to search starting from the newest CUDA version
    default_dirs = valid_dirs(default_dirs)
    if !isempty(default_dirs)
        @debug "Looking for CUDA toolkit via default installation directories" dirs=default_dirs
        append!(dirs, default_dirs)
    end

    # filter
    dirs = valid_dirs(dirs)
    @debug "Found CUDA toolkit at $(join(dirs, ", "))"
    return dirs
end

# figure out the CUDA toolkit version (by looking at the output of a tool like `nvdisasm`)
function parse_toolkit_version(tool, tool_path::String)
    # parse the version string
    verstr = withenv("LANG"=>"C") do
        read(`$tool_path --version`, String)
    end
    m = match(r"\bV(?<major>\d+).(?<minor>\d+).(?<patch>\d+)\b", verstr)
    if m === nothing
        @error "Could not parse CUDA version info (\"$verstr\"); please file an issue."
        return nothing
    end

    version = VersionNumber(parse(Int, m[:major]),
                            parse(Int, m[:minor]),
                            parse(Int, m[:patch]))

    if version >= v"11"
        # starting with CUDA 11, binaries are versioned independently
        # NOTE: we can't always tell, e.g. nvdisasm is the same in CUDA 11.1.0 and 11.1.1.
        #       return the lowest version to ensure compatibility.
        for toolkit_version in sort(collect(keys(cuda_binary_versions)))
            if haskey(cuda_binary_versions[toolkit_version], tool) &&
               cuda_binary_versions[toolkit_version][tool] == version
                @debug "CUDA toolkit identified as $toolkit_version (providing $tool $version)"
                return toolkit_version
            end
        end
        @error "CUDA.jl does not yet support CUDA with $tool $version; please file an issue."
        return nothing
    else
        @debug "CUDA toolkit identified as $version"
        return version
    end
end

"""
    find_libdevice(toolkit_dirs::Vector{String})

Look for the CUDA device library supporting `targets` in any of the CUDA toolkit directories
`toolkit_dirs`. On CUDA >= 9.0, a single library unified library is discovered and returned
as a string. On older toolkits, individual libraries for each of the targets are returned as
a vector of strings.
"""
function find_libdevice(toolkit_dirs)
    @debug "Request to look for libdevice" locations=toolkit_dirs

    # figure out locations
    dirs = String[]
    for toolkit_dir in toolkit_dirs
        push!(dirs, toolkit_dir)
        push!(dirs, joinpath(toolkit_dir, "libdevice"))
        push!(dirs, joinpath(toolkit_dir, "nvvm", "libdevice"))
    end

    # filter
    dirs = valid_dirs(dirs)
    @debug "Look for libdevice" locations=dirs

    for dir in dirs
        path = joinpath(dir, "libdevice.10.bc")
        if isfile(path)
            @debug "Found unified device library at $path"
            return path
        end
    end

    return nothing
end

"""
    find_libcudadevrt(toolkit_dirs::Vector{String})

Look for the CUDA device runtime library in any of the CUDA toolkit directories
`toolkit_dirs`.
"""
function find_libcudadevrt(toolkit_dirs)
    locations = toolkit_dirs
    @debug "Request to look for libcudadevrt " locations

    name = nothing
    if Sys.isunix()
        name = "libcudadevrt.a"
    elseif Sys.iswindows()
        name = "cudadevrt.lib"
    else
        error("No support for discovering the CUDA device runtime library on your platform, please file an issue.")
    end

    # figure out locations
    all_locations = String[]
    for location in locations
        push!(all_locations, location)
        if Sys.iswindows()
            if Sys.WORD_SIZE == 64
                push!(all_locations, joinpath(location, "lib", "x64"))
            elseif Sys.WORD_SIZE == 32
                push!(all_locations, joinpath(location, "lib", "Win32"))
            end
        else
            push!(all_locations, joinpath(location, "lib"))
            if Sys.WORD_SIZE == 64
                push!(all_locations, joinpath(location, "lib64"))
            end
        end
    end

    @debug "Looking for CUDA device runtime library $name" locations=all_locations
    paths = filter(isfile, map(location->joinpath(location, name), all_locations))

    if isempty(paths)
        return nothing
    else
        path = first(paths)
        @debug "Found CUDA device runtime library $(basename(path)) at $(dirname(path))"
        return path
    end
end


text:
@userplot MtxPowersPlot

"""
    mtxpowersplot(A, nmax)

Plot the matrix norm of `A^k` for `k` up to `nmax`

# Examples
```jldoctest
julia> opts = Dict{Symbol,Any}(:npts => 40,:ax => [-110,110,-110,110])
julia> A = diagm(1 => 100.0*ones(4)) + diagm(-2 => 1e-4*ones(3))
5×5 Matrix{Float64}:
 0.0     100.0       0.0       0.0    0.0
 0.0       0.0     100.0       0.0    0.0
 0.0001    0.0       0.0     100.0    0.0
 0.0       0.0001    0.0       0.0  100.0
 0.0       0.0       0.0001    0.0    0.0
 julia> ps_data = new_matrix(A,opts)
 [...]
 julia> driver!(ps_data, opts)
 julia> mtxpowersplot(ps_data, nmax=20)
```
"""
function mtxpowersplot end

RecipesBase.@recipe function f(p::MtxPowersPlot;
    nmax::Integer = 50,
    lbmethod = :none,
    lbdk = 0.25
    )
    if length(p.args) != 1 || !(typeof(p.args[1]) <: PSAStruct)
        error("mtxpowersplot must be given a PSAStruct. Got $(typeof(p.args))")
    end
    ps_data::PSAStruct = p.args[1]
    A = ps_data.input_matrix

    powers, transient = get_mtxpowersnorm(A, nmax)

    alp = maximum(abs.(ps_data.ps_dict[:ews]))

    # NB. We need Plots to use @layout
    #layout := @layout [a; b]
    #link := :x # link the x-axis of both subplots

    # normal plot
    @series begin
        seriestype := :path
        subplot := 1
        label := "Matrix norm"
        powers, transient
    end
    @series begin
        seriestype := :path
        subplot := 1
        label := "Spectral Radius"
        powers, alp.^powers
    end
    # log plot
    #=
    @series begin
        seriestype := :path
        subplot := 2
        label := ""
        yscale := :log10
        powers, transient
    end
    @series begin
        seriestype := :path
        subplot := 2
        label := ""
        yscale := :log10
        powers, alp.^powers
    end
    =#
end


text:
"""
    generate_evaluate(E, ::Type{T})

Generate the statements for the evaluation of the polynomial with exponents `E`.
This assumes that E is in reverse lexicographic order.
"""
function generate_evaluate(E, ::Type{T}, access_input) where T
    exprs = []
    coeffs = map(i -> :(c[$i]), 1:size(E, 2))
    last_expr = generate_evaluate!(exprs, E, T, size(E, 1), 1, coeffs, access_input)
    Expr(:block, exprs..., last_expr)
end

function generate_evaluate!(exprs, E, ::Type{T}, nvar, nterm, coeffs, access_input) where T
    m, n = size(E)

    if n == 1
        return first(monomial_product(T, E[:,1], coeffs[nterm], access_input=access_input))
    end

    if m == 1
        return evalpoly(T, E[1,:], coeffs[nterm:nterm+n-1], access_input(nvar))
    end

    # Recursive
    rec_coeffs = []
    sub_nterm = nterm
    degrees, submatrices = degrees_submatrices(E)
    for E_d in submatrices
        coeff_subexpr = generate_evaluate!(exprs, E_d, T, nvar - 1, sub_nterm, coeffs, access_input)
        # If the returned value is just a symbol propagate it
        if isa(coeff_subexpr, Symbol)
            push!(rec_coeffs, coeff_subexpr)
        else
            @gensym c
            push!(exprs, :($c = $coeff_subexpr))
            push!(rec_coeffs, c)
        end
        sub_nterm += size(E_d, 2)
    end

    return evalpoly(T, degrees, rec_coeffs, access_input(nvar))
end


text:
@testset "1094.car-pooling.jl" begin
    @test car_pooling([[2,1,5],[3,3,7]], 4) == false
    @test car_pooling([[2,1,5],[3,3,7]], 5) == true
    @test car_pooling([[2,1,5],[3,5,7]], 3) == true
    @test car_pooling([[3,2,7],[3,7,9],[8,3,9]], 11) == true
end


text:
#-----------------------------------------------------------------------# RegERM
mutable struct RegERM{L <: Loss, P <: Penalty} <: AbstractMLModel
    β0::Float64
    β::Vector{Float64}
    λ::Vector{Float64}
    loss::L
    penalty::P
end
RegERM(loss, penalty, λ=0.0) = RegERM(0.0, zeros(0), fill(λ, 1), loss, penalty)

# struct StatsModel{L, P, D, T}
#     β0::T
#     β::Vector{T}
#     λ::Vector{T}
#     loss::L
#     penalty::P
#     data::D
# end
# function StatsModel(data::D, loss::L, penalty::P; kw...) where {D<:StatsModelData,L,P}
#     T = eltype(data.x)
#     p = nvars(data)
#     StatsModel{L,P,D,T}(zero(T), zeros(T, p), zeros(T, p), loss, penalty, data)
# end
# function Base.show(io::IO, o::StatsModel)
#     println(io, "StatsModel")
#     println(io, "  > β0      : ", o.β0)
#     println(io, "  > β       : ", o.β)
#     println(io, "  > λ       : ", o.λ)
#     println(io, "  > loss    : ", o.loss)
#     println(io, "  > penalty : ", o.penalty)
#     print(io,   "  > data    : ", o.data)
# end

# function StatsModel(x, y, loss::Loss = L2DistLoss(), penalty::Penalty = NoPenalty())
#     StatsModel(StatsModelData(x, y), loss, penalty)
# end


text:
using Configurations
import TOML
import Pluto: Pluto, Token, Notebook
export NotebookSession,
    SliderServerSettings, ExportSettings, PlutoDeploySettings, get_configuration
using TerminalLoggers: TerminalLogger

###
# SESSION DEFINITION

abstract type RunResult end

Base.@kwdef struct RunningNotebook <: RunResult
    path::String
    notebook::Pluto.Notebook
    original_state::Any
    token::Token = Token()
    bond_connections::Dict{Symbol,Vector{Symbol}}
end
Base.@kwdef struct FinishedNotebook <: RunResult
    path::String
    original_state::Any
end

Base.@kwdef struct NotebookSession{
    C<:Union{Nothing,String},
    D<:Union{Nothing,String},
    R<:Union{Nothing,RunResult},
}
    path::String
    current_hash::C
    desired_hash::D
    run::R
end


text:
facts("Flatten") do

    a = Signal(0)
    b = Signal(1)

    c = Signal(a)

    d = flatten(c)
    cnt = foldp((x, y) -> x+1, 0, d)

    context("Signal{Signal} -> flat Signal") do
        # Flatten implies:
        @fact value(c) --> a
        @fact value(d) --> value(a)
    end

    context("Initial update count") do
        @fact value(cnt) --> 0
    end

    context("Current signal updates") do
        push!(a, 2)
        step()

        @fact value(cnt) --> 1
        @fact value(d) --> value(a)
    end

    context("Signal swap") do
        push!(c, b)
        step()
        @fact value(cnt) --> 2
        @fact value(d) --> value(b)

        push!(a, 3)
        step()
        @fact value(cnt) --> 2
        @fact value(d) --> value(b)

        push!(b, 3)
        step()

        @fact value(cnt) --> 3
        @fact value(d) --> value(b)
    end

    context("Subtle sig swap issue") do
        # When a node, a map (e) in this case, has a flatten as a parent, but also
        # a signal that is the flatten parent sigsig's (c's) current value ("a" here)
        # then when the sigsig gets pushed another value, you want the map to
        # still update on changes to a, even after the map is "rewired" when a
        # new value is pushed to c.

        a = Signal(1)
        b = Signal(2)
        c = Signal(a)
        d = flatten(c)
        e = map(*, a, d) #e is dependent on "a" directly and through d

        @fact value(e) --> 1
        push!(a, 3)
        step()
        @fact value(a) --> 3
        @fact value(d) --> 3
        @fact value(e) --> 9

        push!(c, b)
        @fact value(e) --> 9 # no change until step
        step()
        @fact value(d) --> 2 # d now takes b's value
        @fact value(e) --> 6 # e == d * a == 2 * 3 == 6

        # the push!(c, b) should have triggered a "rewiring" of the graph
        # so that updates to b affect d and e
        push!(b, 9)
        step()
        @fact value(d) --> 9  # d now takes b's value
        @fact value(e) --> 27 # e == d * a == 9 * 3 == 27

        # changes to a should still affect e (but not d)
        push!(a, 4)
        @fact value(e) --> 27 # no change until step
        step()
        @fact value(a) --> 4
        @fact value(d) --> 9 # no change to d
        @fact value(e) --> 36 # a*d == 4 * 9

        @fact queue_size() --> 0
    end

    context("Sigsig's value created after SigSig") do
        # This is why we need bind! in flatten implementation rather than just
        # setting the flatten's parents to (input, current_node) every time
        # input updates. That won't work if the current_node was created after
        # the flatten (e.g. in the below after pushing a new value to `a`),
        # because updates to the current_node happen further down the `nodes`
        # list than the flatten, so the flatten doesn't get updated.
        a = Signal(number())
        local b
        c = foreach(a) do av
                b = Signal(av)
                foreach(identity, b)
            end
        d = flatten(c)
        b_orig = b

        @fact d.value --> a.value

        push!(b, number())
        step()
        @fact d.value --> b.value

        push!(a, number())
        step()
        @fact d.value --> a.value
        @fact b_orig --> not(b)

        push!(b, number())
        step()
        @fact d.value + 1 --> b.value + 1
    end
end


text:
# ------------ Switch trace ------------ #

struct SwitchTrace{T} <: Trace
    gen_fn::GenerativeFunction{T}
    index::Int
    branch::Trace
    retval::T
    args::Tuple
    score::Float64
    noise::Float64
end

@inline get_choices(tr::SwitchTrace) = get_choices(tr.branch)
@inline get_retval(tr::SwitchTrace) = tr.retval
@inline get_args(tr::SwitchTrace) = tr.args
@inline get_score(tr::SwitchTrace) = tr.score
@inline get_gen_fn(tr::SwitchTrace) = tr.gen_fn
@inline Base.getindex(tr::SwitchTrace, addr) = Base.getindex(tr.branch, addr)
@inline project(tr::SwitchTrace, selection::Selection) = project(tr.branch, selection)
@inline project(tr::SwitchTrace, ::EmptySelection) = tr.noise

max_stars_repo_name: georgematheos/Gen

max_issues_repo_name: georgematheos/Gen

max_forks_repo_name: georgematheos/Gen


text:
export iar_chebyshev
export compute_y0_cheb

using IterativeSolvers
using LinearAlgebra
using Random

# Types specifying which way to compute y0 in chebyshev iar
abstract type ComputeY0Cheb end;
abstract type ComputeY0ChebDEP <: ComputeY0Cheb end;
abstract type ComputeY0ChebPEP <: ComputeY0Cheb end;
abstract type ComputeY0ChebSPMF_NEP <: ComputeY0Cheb end;
abstract type ComputeY0ChebAuto <: ComputeY0Cheb end;

# Data collected in a precomputation phase.
# These are made mutable (could be made immutable by appropriate modification in precompute_data)
abstract type AbstractPrecomputeData end
mutable struct PrecomputeDataDEP <: AbstractPrecomputeData
    Tc; Ttau;
end
mutable struct PrecomputeDataPEP <: AbstractPrecomputeData
    Tc; D;
end
mutable struct PrecomputeDataSPMF <: AbstractPrecomputeData
    Tc; DDf;
end
mutable struct PrecomputeDataNEP <: AbstractPrecomputeData
    P; P_inv; α; γ; σ;
end

"""
    iar_chebyshev(nep,[maxit=30,][σ=0,][γ=1,][linsolvecreator=DefaultLinSolverCreator(),][tolerance=eps()*10000,][neigs=6,][errmeasure,][v=rand(size(nep,1),1),][logger=0,][check_error_every=1,][orthmethod=DGKS][a=-1,][b=1,][compute_y0_method=ComputeY0ChebAuto])

Run the infinite Arnoldi method (Chebyshev version) on the nonlinear eigenvalue problem stored in `nep`.

The target `σ` is the center around which eiganvalues are computed. A Ritz pair `λ` and `v` is flagged a as converged (to an eigenpair) if `errmeasure` is less than `tol`. The vector
`v` is the starting vector for constructing the Krylov space. The orthogonalization method, used in contructing the orthogonal basis of the Krylov space, is specified by `orthmethod`, see the package `IterativeSolvers.jl`. The iteration
is continued until `neigs` Ritz pairs converge. This function throws a `NoConvergenceException` if the wanted eigenpairs are not computed after `maxit` iterations.
However, if `neigs` is set to `Inf` the iteration is continued until `maxit` iterations without an error being thrown.
The kwarg `compute_y0_method` specifying how the next vector of the Krylov space (in Chebyshev format) can be computed. See [`compute_y0_cheb`](@ref) in the module NEPSolver with the command `?NEPSolver.compute_y0_cheb`.

See [`augnewton`](@ref) for other parameters.

# Example
```julia-repl
julia> using NonlinearEigenproblems, LinearAlgebra
julia> nep=nep_gallery("dep0",100);
julia> v0=ones(size(nep,1));
julia> λ,v=iar_chebyshev(nep;v=v0,tol=1e-5,neigs=3);
julia> norm(compute_Mlincomb!(nep,λ[1],v[:,1])) # Is it an eigenvalue?
julia> λ    # print the computed eigenvalues
3-element Array{Complex{Float64},1}:
julia> norm(compute_Mlincomb(nep,λ[1],v[:,1]))
 0.050462487848960284 - 1.4289626573515395e-18im
 -0.07708779190301127 + 7.703053374113074e-18im
   0.1503856540695659 - 1.662582577182149e-17im
```

# References
* Algorithm 2 in Jarlebring, Michiels Meerbergen, A linear eigenvalue algorithm for the nonlinear eigenvalue problem, Numer. Math, 2012
"""
iar_chebyshev(nep::NEP;params...)=iar_chebyshev(ComplexF64,nep;params...)
function iar_chebyshev(
    ::Type{T},
    nep::NEP;
    orthmethod::Type{T_orth}=DGKS,
    maxit=30,
    linsolvercreator=DefaultLinSolverCreator(),
    tol=eps(real(T))*10000,
    neigs=6,
    errmeasure::ErrmeasureType = DefaultErrmeasure(nep),
    σ=zero(T),
    γ=one(T),
    v=randn(real(T),size(nep,1)),
    logger=0,
    check_error_every=1,
    compute_y0_method::Type{T_y0}=ComputeY0ChebAuto,
    a = isa(nep,DEP) ? -maximum(nep.tauv) : -1.0,
    b = isa(nep,DEP) ? 0.0 : 1.0
    )where{T,T_orth<:IterativeSolvers.OrthogonalizationMethod,T_y0<:ComputeY0Cheb}

    @parse_logger_param!(logger)

    if (compute_y0_method == ComputeY0ChebAuto)
        if (isa(nep,DEP))
            compute_y0_method=ComputeY0ChebDEP;
        elseif (isa(nep,PEP))
            compute_y0_method=ComputeY0ChebPEP;
        elseif  (isa(nep,SPMF_NEP))
            compute_y0_method=ComputeY0ChebSPMF_NEP;
        else
            compute_y0_method=ComputeY0Cheb;
        end
    end

    if ( σ!=zero(T) || γ!=one(T) ) && compute_y0_method<:Union{ComputeY0ChebDEP,ComputeY0ChebPEP}
            @warn "The problem will be explicitly shifted and scaled. The shift and scaling feature is not supported in the general version of iar_chebyshev."
            # TODO: use the original errmeasure and not compute_resnorm. I don't know why doesn't work
            #errmeasure=(μ,v)-> errmeasure(σ+γ*μ,v) #this is what we want but it does not work
            errmeasure=function (μ,v) compute_resnorm(nep,σ+γ*μ,v) end
            nep=shift_and_scale(nep,shift=σ,scale=γ);
            σ_orig=σ; γ_orig=γ
            σ=zero(T); γ=one(T)
    end
    push_info!(logger,"IAR Chebyshev with interval [$a,$b]");

    cc=(a+b)/(a-b);   kk=2/(b-a); # scale and shift parameters for the Chebyshev basis

    n = size(nep,1);
    m = maxit;

    # initialization
    V = zeros(T,n*(m+1),m+1);
    H = zeros(T,m+1,m);
    y = zeros(T,n,m+1);
    α=Vector{T}(γ.^(0:m)); α[1]=zero(T);
    local M0inv::LinSolver=create_linsolver(linsolvercreator,nep,σ)
    err = ones(m,m);
    λ=zeros(T,m+1); Q=zeros(T,n,m+1);

    vv=view(V,1:1:n,1); # next vector V[:,k+1]
    vv[:]=v; vv[:]=vv[:]/norm(vv);
    k=1; conv_eig=0;

    # hardcoded matrix L
    L=diagm(0 => vcat(2, 1 ./ (2:m))) + diagm(-2 => -vcat(1 ./ (1:(m-2))))
    L=L*(b-a)/4

    # precomputation for exploiting the structure DEP, PEP, GENERAL
    precomp=precompute_data(T,nep,compute_y0_method,a,b,maxit,γ,σ)

    while (k <= m) && (conv_eig<neigs)

        VV=view(V,1:1:n*(k+1),1:k); # extact subarrays, memory-CPU efficient
        vv=view(V,1:1:n*(k+1),k+1); # next vector V[:,k+1]

        # compute y (only steps different then standard IAR)
        y=zeros(T,n,k+1);
        if (compute_y0_method != ComputeY0Cheb)
            y[:,2:k+1]  = reshape(VV[1:1:n*k,k],n,k)*L[1:k,1:k];
        end
        y[:,1] = compute_y0_cheb(T,nep,compute_y0_method,reshape(VV[1:1:n*k,k],n,k),
                                 y,M0inv,precomp);

        vv[:]=reshape(y[:,1:k+1],(k+1)*n,1);
        # orthogonalization
        H[k+1,k] = orthogonalize_and_normalize!(VV, vv, view(H,1:k,k), orthmethod)

        # compute Ritz pairs (every check_error_every iterations)
        if ((rem(k,check_error_every)==0)||(k==m))&&(k>2)
            D,Z = eigen(H[1:k,1:k])
            VV=view(V,1:1:n,1:k);
            Q=VV*Z
            λ=σ .+ γ ./ D
            conv_eig=0;
            # compute the errors
            err[k,1:size(λ,1)]=
              map(s-> estimate_error(errmeasure,λ[s],Q[:,s]), 1:size(λ,1))
            # Log them and compute the converged
            push_iteration_info!(logger,2, k,err=err[k,1:size(λ,1)],
                                 continues=true);
            for s=1:size(λ,1)
                if err[k,s]<tol;
                    conv_eig=conv_eig+1;
                    push_info!(logger,"+", continues=true);
                elseif err[k,s]<tol*10
                    push_info!(logger,"=", continues=true);
                else
                    push_info!(logger,"-", continues=true);
                end
            end
            push_info!(logger,"");
            # Sort the errors
            idx=sortperm(err[k,1:k]); # sort the error
            err[k,1:k]=err[k,idx];

            # extract the converged Ritzpairs
            if (k==m)||(conv_eig>=neigs)
                nrof_eigs = Int(min(length(λ),neigs))
                λ=λ[idx[1:nrof_eigs]]
                Q=Q[:,idx[1:nrof_eigs]]
            end
        end
        k=k+1;
    end
    k=k-1
    # NoConvergenceException
    if conv_eig<neigs && neigs != Inf
        err=err[end,1:neigs];
        idx=sortperm(err); # sort the error
        λ=λ[idx];  Q=Q[:,idx]; err=err[idx];
        msg="Number of iterations exceeded. maxit=$(maxit)."
        if conv_eig<3
            msg=string(msg, " Check that σ is not an eigenvalue.")
        end
        throw(NoConvergenceException(λ,Q,err,msg))
    end

    # eventually shift and rescale the eigenvalues if the problem was shifted and rescaled
    # TODO: we shouldn't use the exception system for this, perhaps a Bool flag is a better choice
    try
        λ = σ_orig .+ γ_orig * λ
    catch
        # ignore
    end

    # extract the converged Ritzpairs
    λ=λ[1:min(length(λ),conv_eig)];
    Q=Q[:,1:min(size(Q,2),conv_eig)];
    return λ,Q,err[1:k,:],V[:,1:k],H[1:k,1:k]
end

# Contructors for the precomputed data
function PrecomputeDataInit(::Type{ComputeY0ChebDEP})
    return PrecomputeDataDEP(0,0);
end
function PrecomputeDataInit(::Type{ComputeY0ChebPEP})
    return PrecomputeDataPEP(0,0);
end
function PrecomputeDataInit(::Type{ComputeY0ChebSPMF_NEP})
    return PrecomputeDataSPMF(0,0);
end
function PrecomputeDataInit(::Type{ComputeY0Cheb})
    return PrecomputeDataNEP(0,0,0,0,0);
end

# Precompute data (depending on NEP-type and y0 computation method)
function precompute_data(T,nep::NEPTypes.DEP,::Type{ComputeY0ChebDEP},a,b,m,γ,σ)
    if ( σ!=zero(T) || γ!=one(T) )
        error("This function does not support shift and scale parameters");
    end
    cc=(a+b)/(a-b);   kk=2/(b-a); # scale and shift parameters for the Chebyshev basis
    precomp=PrecomputeDataInit(ComputeY0ChebDEP);
    # creating the vector containing T_i(c)
    precomp.Tc=cos.((0:m)'.*acos(cc));  # vector containing T_i(c)
    # creating the matri  containing T_i(-kk*nep.tauv[j]+cc)
    precomp.Ttau=zeros(T,length(nep.tauv),m+2);    II=(0:m+1)';
    for j=1:length(nep.tauv)
        tauv_SS=-kk*nep.tauv[j]+cc;
        if abs(tauv_SS)<=1
            precomp.Ttau[j,:]=cos.(II.*acos(tauv_SS));
        elseif tauv_SS>=1
            precomp.Ttau[j,:]=cosh.(II.*acosh(tauv_SS));
        else
            precomp.Ttau[j,:]=((-1).^II).*cosh.(II.*acosh(-tauv_SS));
        end
    end
    return precomp;
end
function precompute_data(T,nep::NEPTypes.PEP,::Type{ComputeY0ChebPEP},a,b,m,γ,σ)
    # The matrix D is the derivation map, namely D_N [T_0(s),T_1(s), ..., T_{N-1}(s)]=[T_0'(s),T_1'(s), ..., T_{N-1}'(s)]
    # NOTE: The matrix D_N is the submatrix D_M[1:N+1,1:N+1] for M>=N (larger derivarive matrix)
    if ( σ!=zero(T) || γ!=one(T) )
        error("This function does not support shift and scale parameters");
    end
    cc=(a+b)/(a-b);   kk=2/(b-a); # scale and shift parameters for the Chebyshev basis
    precomp=PrecomputeDataInit(ComputeY0ChebPEP);
    precomp.Tc=cos.((0:m)'.*acos(cc));  # vector containing T_i(c)
    L=diagm(0 => vcat(2, 1 ./ (2:m))) + diagm(-2 => -vcat(1 ./ (1:(m-2))))
    L=L*(b-a)/4
    L=inv(L[1:m,1:m]); precomp.D=vcat(zeros(1,m),L[1:m-1,:]);
    return precomp;
end
function precompute_data(T,nep::NEPTypes.AbstractSPMF,::Type{ComputeY0ChebSPMF_NEP},a,b,m,γ,σ)
    # Precomputation for NEPS in SPMF format
    cc=(a+b)/(a-b);   kk=2/(b-a); # scale and shift parameters for the Chebyshev basis
    precomp=PrecomputeDataInit(ComputeY0ChebSPMF_NEP);
    precomp.Tc=cos.((0:m)'.*acos(cc));  # vector containing T_i(c)
    L=diagm(0 => vcat(2, 1 ./ (2:m))) + diagm(-2 => -vcat(1 ./ (1:(m-2))))
    L=L*(b-a)/4; L=inv(L[1:m,1:m]); # integration matrix
    D=vcat(zeros(1,m),L[1:m-1,:]); # detivation matrix

    fv, Av = get_fv(nep), get_Av(nep)
    DDf = Array{Array{T,2}}(undef, length(fv))
    for i = 1:length(fv)
        DDs = Matrix{T}(σ*I, size(D,1), size(D,1)) + γ*D
        DDf[i] = γ * DD0_mat_fun(T, fv[i], DDs, σ)
    end
    precomp.DDf = DDf
    return precomp
end
function precompute_data(T,nep::NEPTypes.NEP,::Type{ComputeY0Cheb},a,b,m,γ,σ)
    # Precomputation for a NEP that is not a DEP, PEP or SPMF
    @warn "The nep does not belong to the class of DEP, PEP or SPMF and the function compute_y0 is not provided. Check if the nep belongs to such classes and define it accordingly or provide the function compute_y0. If none of these options are possible, the method will be based on the convertsion between Chebyshev and monomial base and may be numerically unstable if many iterations are performed."
    cc=(a+b)/(a-b);   kk=2/(b-a); # scale and shift parameters for the Chebyshev basis
    precomp=PrecomputeDataInit(ComputeY0Cheb);
    precomp.P = mapslices(x->cheb2mon(T,kk,cc,x), Matrix{T}(I,m+1,m+1), dims = 1)'        # P maps chebyshev to monomials as matrix vector action
    precomp.P_inv = mapslices(x->mon2cheb(T,kk,cc,x), Matrix{T}(I,m+1,m+1), dims = 1)'    # P_inv maps monomials to chebyshev as matrix vector action
    precomp.σ=σ;
    precomp.γ=γ;
    precomp.α=γ.^(0:m);
    return precomp;
end
"""
    y0 = compute_y0_cheb([eltype],nep::NEPTypes.DEP,::Type{ComputeY0ChebPEP},X,Y,M0inv,precomp::AbstractPrecomputeData)

Computes the vector y0 used in [`iar_chebyshev`](@ref) given by
```math
 y_0 = \\sum_{i=1}^N T_{i-1}(γ) x_i - \\sum_{j=1}^m A_j \\left( \\sum_{i=1}^{N+1} T_{i-1}(-ρ \\tau_j+γ) y_i \\right )
```
where T(c) is the vector containing \$T_i(c)\$ as coefficients, where \$T_i\$ is the i-th Chebyshev polynomial of the first kind.
"""
function compute_y0_cheb(T,nep::NEPTypes.DEP,::Type{ComputeY0ChebDEP},X,Y,M0inv,precomp::AbstractPrecomputeData)
    Tc=precomp.Tc;
    Ttau=precomp.Ttau;
    Av=get_Av(nep);

    n,N=size(X)
    y0=sum(broadcast(*,X,view(Tc,1:1,1:N)), dims = 2) # \sum_{i=1}^N T_{i-1}(γ) x_i
    for j=1:length(nep.tauv) # - \sum_{j=1}^m A_j \left( \sum_{i=1}^{N+1} T_{i-1}(-ρ \tau_j+γ) y_i\right )
        y0-=Av[j+1]*sum(broadcast(*,Y,view(Ttau,j:j,1:N+1)), dims = 2)
    end
    y0=lin_solve(M0inv,y0)
    return y0
end
"""
    y0 = compute_y0_cheb([eltype],nep::NEPTypes.PEP,::Type{ComputeY0ChebPEP},X,Y,M0inv,precomp::AbstractPrecomputeData)

Computes the vector y0 used in [`iar_chebyshev`](@ref) given by
```math
 y_0 = \\sum_{j=0}^{d-1} A_{j+1} x D^j T(c) - y T(c)
```
where T(c) is the vector containing \$T_i(c)\$ as coefficients, where \$T_i\$ is the i-th Chebyshev polynomial of the first kind and \$D\$ is the derivation matrix in Chebyshev basis.
"""
function compute_y0_cheb(T,nep::NEPTypes.PEP,::Type{ComputeY0ChebPEP},X,Y,M0inv,precomp::AbstractPrecomputeData)
    Tc=precomp.Tc;
    n,N=size(X)
    d=length(nep.A)-1;  # degree of the PEP
    # sum for every coefficiet
    v=view(Tc,1:N);
    y0=zeros(T,n,1);
    for j=0:d-1
        y0+=nep.A[j+2]*(X*v);
        v=view(precomp.D,1:N,1:N)*v;
    end
    y0=-lin_solve(M0inv,y0)
    y0-=Y*(view(Tc,1:N+1));
    return y0
end
"""
    y0 = compute_y0_cheb([eltype],nep::NEPTypes.SPMF_NEP,::Type{ComputeY0ChebPEP},X,Y,M0inv,precomp::AbstractPrecomputeData)

Computes the vector y0 used in [`iar_chebyshev`](@ref) given by
```math
 y_0= \\sum_{j=0}^{m} M^{(j)}(\\mu) X b_j \\left( D_N \\right) T_N(c) - Y T_N(c)
```
where T(c) is the vector containing \$T_i(c)\$ as coefficients, where \$T_i\$ is the i-th Chebyshev polynomial of the first kind and \$b_j(\\lambda)=(f_j(0)-f_j(\\lambda))/\\lambda=f[\\lambda,0]\$ are divided differences.
"""
function compute_y0_cheb(T,nep::NEPTypes.AbstractSPMF,::Type{ComputeY0ChebSPMF_NEP},X,Y,M0inv,precomp::AbstractPrecomputeData)
    Tc=precomp.Tc;
    n,N=size(X);
    fv,Av=get_fv(nep),get_Av(nep)
    y0=zero(X)
    for i=1:length(fv)
        y0+=Av[i]*X*view(precomp.DDf[i],1:N,1:N)
    end
    y0=y0*(vec(view(Tc,1:1,1:N)))
    y0=-lin_solve(M0inv,y0)
    y0-=Y*(view(Tc,1:N+1));
    return y0
end
"""
    y0 = compute_y0_cheb([eltype],nep::NEPTypes.NEP,::Type{ComputeY0ChebNEP},X,Y,M0inv,precomp::AbstractPrecomputeData)

Computes the vector y0 used in [`iar_chebyshev`](@ref) defined as
```math
 y_0 =\\left( \\sum_{i=0}^{N-1} B \\left( \\frac{d}{d \\theta} \\right) \\hat T_i(\\theta) x_i \\right)(0) - \\sum_{i=0}^{N} T_i(c) y_i
```
where \$T_i\$ is the i-th Chebyshev polynomial of the first kind, \$ \\ hat T_i\$ is the i-th Chebyshev polynomial of the first kind for the interval [a,b]. For a generic `nep`, this quantity is computed by converting polynomials in monomial basis. This procedure may be numerical unstable if many iterations are required. If for the specific `nep` a closed formula is available, we suggest to overload this function.
"""
function compute_y0_cheb(T,nep::NEPTypes.NEP,::Type{ComputeY0Cheb},X,Y,M0inv,precomp::AbstractPrecomputeData)
    # compute_y0_cheb computes y0 for the NEP
    # This function convert the coefficients x=vec(X) from the Chebyshev basis to the monomial basis, then compute the vector y=vec(Y) as in iar (Taylor version) and then convert y=vec(Y) in Chebyshev basis
    k=size(X,2);
    α=precomp.α; σ=precomp.σ;
    Y[:,2:k+1] = X*view(precomp.P,1:k,1:k);
    broadcast!(/,view(Y,:,2:k+1),view(Y,:,2:k+1),(1:k)')
    Y[:,1] = compute_Mlincomb(nep,σ,view(Y,:,1:k+1),α[1:k+1]);
    Y[:,1] = -lin_solve(M0inv,Y[:,1]);
    Y[:,:]=Y*view(precomp.P_inv,1:k+1,1:k+1);
    return Y[:,1];
end
function mon2cheb(T,ρ,γ,a)
#cheb2mon: Monomial basis to shifted-and-scaled Chebyshev basis conversion.
#    c = cheb2mon(rho,gamma,a) converts a polynomial written in
#    monomial basis
#    p(x)=\sum_{j=0}^n a(j+1) x^j
#    to a polynomial written in shifted-and-scaled Chebyshev basis
#    p(x)=\sum_{j=0}^n c(j+1) T_j(rho*x+gamma)
#    where T_j(x) is the j-th Chebyshev polynomial of the first kind.
#
#    Example:
#    Suppose we have a polynomial in the monomial basis:
#    2-4*x-x^2
#    This polynomial is expressed in shifted-and-scaled Chebyshev basis as:
#    c(1)*T_0(2*x+3)+c(2)*T_1(2*x+3)+a(3)*T_2(2*x+3)
#    where
#    c=cheb2mon(2,3,[2 -4 -1])
#
#    This function is based on
#    MATH77 and mathc90, Release 6.0, Libraries of Mathematical Subprograms in Fortran 77 and C
#    Fred T. Krogh, Charles L. Lawson, W. Van Snyder
#    in chapter 11 section 3
#    available in the following link
#    http://www.netlib.org/math/docpdf/ch11-03.pdf
    n=length(a)-1;
    α=1/(2*ρ);    β=-γ/ρ;
    b=zeros(T,n+3,1);
    bb=zeros(T,n+3,1);

    for j=n:-1:0
        bb[1]=α*b[2]+β*b[1]+a[j+1];
        bb[2]=β*b[2]+α*b[3]+2*α*b[1];
        for k=3:n-j-1
            bb[k]=α*b[k-1]+β*b[k]+α*b[k+1];
        end
        if n-j>2
            bb[n-j]=α*b[n-j-1]+β*b[n-j];
        end
        if n-j+1>2
            bb[n-j+1]=α*b[n-j];
        end
        b=bb;
        bb=0*bb;
    end
    c=b[1:n+1];
    c=c[:,1];
end
function cheb2mon(T,ρ,γ,c)
#cheb2mon: shifted-and-scaled Chebyshev basis conversion to Monomial basis.
#    a = cheb2mon(T,rho,gamma,c) converts a polynomial written in
#    shifted-and-scaled Chebyshev basis
#    p(x)=\sum_{j=0}^n c(j+1) T_j(rho*x+gamma)
#    to a polynomial written in monomial basis
#    p(x)=\sum_{j=0}^n a(j+1) x^j
#    where T_j(x) is the j-th Chebyshev polynomial of the first kind.
#    T is the type, e.g., Float64
#
#    Example:
#    Suppose we have a polynomial in the shifted-and-scaled Chebyshev basis:
#    2*T_0(2x+3)-4*T_1(2x+3)-T_2(2x+3)
#    This polynomial is expressed in monomial basis as
#    a(1)+a(2)*x+a(3)*x^2
#    where
#    a=cheb2mon(Float64,2,3,[2 -4 -1])
#
#    This function is based on
#    MATH77 and mathc90, Release 6.0, Libraries of Mathematical Subprograms in Fortran 77 and C
#    Fred T. Krogh, Charles L. Lawson, W. Van Snyder
#    in chapter 11 section 3
#    available in the following link
#    http://www.netlib.org/math/docpdf/ch11-03.pdf
    n=length(c)-1;
    α=1/(2*ρ);    β=-γ/ρ;
    a=zeros(T,n+3,1);     b=zeros(T,n+3,1);
    bb=zeros(T,n+3,1);    bb[1:n+1]=c;
    for j=1:1:n+1
        for k=n-j+1:-1:2
            b[k]=(bb[k+1]-β*b[k+1]-α*b[k+2])/α;
        end
        b[1]=(bb[2]-β*b[2]-α*b[3])/(2*α);
        a[j]=bb[1]-α*b[2]-β*b[1];

        bb=b; b=0*b;
    end
    a=a[1:n+1,1];
end
function DD0_mat_fun(T,f,S,σ)
    # evaluete the divided differences matrix function
    # f[S,0] by using the equality
    #
    # f(S I) = (f(S) 		f[S,σI]  )
    #  (0 0)   (0			f(0)	)
    #
    # or, written in a more compact way, it is
    #
    # f([S I; 0 σI])=[f(S) f[S,σI]; 0 f(σ)I]
    #
    # Notice that f[S,0] is defined also for S singular.
    # If S is nonsingular it holds f[S,0]=S^(-1)-(f(S)-f(0))
    # Example:
    # n=10; S=rand(n,n); T=ComplexF64; f=x->exp(x)+x^2
    # Y1=DD0_mat_fun(T,f,S); Y2=inv(S)*(f(S)-f(zeros(S)));
    # opnorm(Y1-Y2)

    n=size(S,1);
    A=zeros(T,2*n,2*n);
    A[1:n, 1:n] = S
    A[1:n, n+1:2*n] = Matrix{T}(I, n, n)
    A[n.+(1:n), n.+(1:n)] = Matrix{T}(σ*I, n, n)
    return f(A)[1:n,n+1:end];
end


text:
using Quon
using Test

function apply_first_match!(r, q)
    ms = match(r, q)
    isempty(ms) && (return q)
    m = first(ms)
    Quon.rewrite!(q, m)
    return q
end

@testset "Genus in Yang-Baxter Eq" begin
    q = Tait(2)
    push_gate!(q, Val(:CZ), 1, 2)
    push_gate!(q, Val(:CZ), 1, 2)

    apply_first_match!(Quon.Rule(:yang_baxter_triangle), q)
    apply_first_match!(Quon.Rule(:yang_baxter_star), q)
    apply_first_match!(Quon.Rule(:perm_rz), q)
    apply_first_match!(Quon.Rule(:yang_baxter_triangle), q)
    apply_first_match!(Quon.Rule(:yang_baxter_star), q)
    apply_first_match!(Quon.Rule(:perm_rz), q)
    apply_first_match!(Quon.Rule(:perm_rz), q)
    apply_first_match!(Quon.Rule(:perm_rz), q)
    apply_first_match!(Quon.Rule(:yang_baxter_triangle), q)
    apply_first_match!(Quon.Rule(:yang_baxter_star), q)
    apply_first_match!(Quon.Rule(:perm_rz), q)
    apply_first_match!(Quon.Rule(:perm_rz), q)
    apply_first_match!(Quon.Rule(:perm_rz), q)
    @test length(match(Quon.Rule(:yang_baxter_star), q)) == 0
end

@testset "add_edge" begin
    q = Tait(2)
    push_gate!(q, Val(:CZ), 1, 2)
    push_gate!(q, Val(:CZ), 1, 2)

    apply_first_match!(Quon.Rule(:yang_baxter_triangle), q)
    apply_first_match!(Quon.Rule(:yang_baxter_triangle), q)
    apply_first_match!(Quon.Rule(:yang_baxter_star), q)
    apply_first_match!(Quon.Rule(:yang_baxter_star), q)
    apply_first_match!(Quon.Rule(:perm_rz), q)
    apply_first_match!(Quon.Rule(:yang_baxter_triangle), q)
    apply_first_match!(Quon.Rule(:perm_rz), q)
    apply_first_match!(Quon.Rule(:yang_baxter_star), q)
    apply_first_match!(Quon.Rule(:yang_baxter_triangle), q)
    apply_first_match!(Quon.Rule(:yang_baxter_star), q)
    apply_first_match!(Quon.Rule(:yang_baxter_triangle), q)
    apply_first_match!(Quon.Rule(:yang_baxter_star), q)
    apply_first_match!(Quon.Rule(:yang_baxter_triangle), q)
    apply_first_match!(Quon.Rule(:perm_rz), q)
    apply_first_match!(Quon.Rule(:yang_baxter_star), q)
    apply_first_match!(Quon.Rule(:yang_baxter_triangle), q)
    apply_first_match!(Quon.Rule(:yang_baxter_star), q)
    apply_first_match!(Quon.Rule(:yang_baxter_triangle), q)
    apply_first_match!(Quon.Rule(:yang_baxter_star), q)
    apply_first_match!(Quon.Rule(:yang_baxter_triangle), q)
    apply_first_match!(Quon.Rule(:yang_baxter_star), q)
    apply_first_match!(Quon.Rule(:yang_baxter_triangle), q)
    apply_first_match!(Quon.Rule(:yang_baxter_star), q)
    apply_first_match!(Quon.Rule(:yang_baxter_triangle), q)
    apply_first_match!(Quon.Rule(:yang_baxter_star), q)
    apply_first_match!(Quon.Rule(:yang_baxter_triangle), q)
    apply_first_match!(Quon.Rule(:yang_baxter_star), q)
    apply_first_match!(Quon.Rule(:yang_baxter_triangle), q)
    apply_first_match!(Quon.Rule(:yang_baxter_star), q)
    apply_first_match!(Quon.Rule(:yang_baxter_triangle), q)
    apply_first_match!(Quon.Rule(:yang_baxter_star), q)
    apply_first_match!(Quon.Rule(:yang_baxter_triangle), q)
    apply_first_match!(Quon.Rule(:yang_baxter_triangle), q)
    apply_first_match!(Quon.Rule(:yang_baxter_star), q)
    apply_first_match!(Quon.Rule(:yang_baxter_triangle), q)
    apply_first_match!(Quon.Rule(:x_fusion), q)
    apply_first_match!(Quon.Rule(:yang_baxter_star), q)
    apply_first_match!(Quon.Rule(:yang_baxter_triangle), q)
    apply_first_match!(Quon.Rule(:identity), q)
    apply_first_match!(Quon.Rule(:perm_rz), q)
    apply_first_match!(Quon.Rule(:perm_rz), q)
    apply_first_match!(Quon.Rule(:perm_rz), q)
    apply_first_match!(Quon.Rule(:perm_rz), q)
    apply_first_match!(Quon.Rule(:perm_rz), q)
    @test true
end


text:
function petk04()
#
# M. Petkovic testing polynomials, p139
#
    z = [3.0; -1*[1;1;1]; 2*im*[1;1;1]; (-2+im)*[1;1]; (-2-im)*[1;1];
         (2+im)*[1;1]; (2-im)*[1;1]];
    p = poly(z);
    z = [3 1; -1 3;  2*im 3; -2+im 2; -2-im 2; 2+im 2; 2-im 2];

    p, PolyZeros(z)
end


text:
using Test
using AutoPreallocation
using AutoPreallocation: AllocationRecord, record_alloc!

@testset "isequal and hash" begin
    record1 = AllocationRecord()
    record_alloc!(record1, ones(3,7))
    record_alloc!(record1, zeros(40))

    record2 = AllocationRecord()
    record_alloc!(record2, zeros(3,7))
    record_alloc!(record2, ones(40))

    @test record1 == record2
    @test isequal(record1, record2)
    @test hash(record1) === hash(record2)

    @test record1 != AllocationRecord()
    @test !isequal(record1, AllocationRecord())
end

@testset "show" begin
    record = AllocationRecord()
    record_alloc!(record, ones(3,7))
    record_alloc!(record, zeros(40))

    @test repr(record) == "AutoPreallocation.AllocationRecord(\n" *
        "    [Array{Float64,2}(undef, (3, 7)), Array{Float64,1}(undef, (40,))],\n" *
        "    [(3, 7), (40,)]\n" *
        ")"

    # Check canon requirement of a good repr:
    @test eval(Meta.parse(repr(record))) == record
end

@testset "reinitialize" begin
    @testset "resizing up: maxsize = $maxsize" for maxsize in (0, 64, 64^3)
        record = AllocationRecord()
        record_alloc!(record, ones(maxsize))

        # Access preallocated item direct to mimic resizing it
        @assert length(record.allocations[1]) == maxsize
        empty!(record.allocations[1]) # shrink it down
        @assert length(record.allocations[1]) == 0
        @test (@ballocated reinitialize!($record)) == 0
        @test length(record.allocations[1]) == maxsize
    end

    @testset "resizing down: maxsize = $maxsize" for maxsize in (0, 64, 64^3)
        record = AllocationRecord()
        record_alloc!(record, Int[])

        # Access preallocated item direct to mimic resizing it
        @assert length(record.allocations[1]) == 0
        append!(record.allocations[1], ones(maxsize)) # grow it
        @assert length(record.allocations[1]) == maxsize
        @test (@ballocated reinitialize!($record)) == 0
        @test length(record.allocations[1]) == 0
    end
end


text:
# call add_procs(n) to add extra processors if desired
# load packages
using KOSHELPER
@everywhere using KOS
using KOS_FIELD_SOLVE
#
# other Julia packages
@everywhere using Interpolations
using PyPlot
#
ion()
# distance function
@everywhere function d_ball_0_1( x::Array{Float64,1} )
  # out=d_ball_0_1(x)
  #
  # Distance of points to the domain boundary
  # Domain is a ball of radius 1 centred at (0,0)
  # x: point
  R=1.
  return R-sqrt(sum(x.^2))
end
#
@everywhere function vec_d_ball_0_1( xs )
  # out=d_ball_0_1(x)
  # same as above
  # points are rows of xs
  R=1.
  return (R-sqrt.(sum(xs.^2,length(size(xs)))))[:]
end
# Zero functions
@everywhere function fn_zero( x::Array{Float64,2})
  #
  return 0.
end
#
@everywhere function fn_zero( x::Array{Float64,1} )
  # Right-hand side designed to expose slow coupling rates.  #
  # x: point(s) as columns for evaluation
  return zeros(size(x)[1])
end
##########################################
function mwe2(alpha)
  # Main routine for computing eigenvalues on unit ball
  # Inputs: alpha, parameter for fractional Laplacian in (0,2)
  # lam_tol, tolerance for the eigenvalue
  lam_tol=5e-2
  # no_its, number of Arnoldi iterations
  no_its=5
  # relax_flag, use variable accuracy if True
  relax_flag=true
  #
  # Ouputs:lam, computed eignvalue
  # el_time, elapsed time in seconds
  # eig_res, eigenvalue residual
  println("\n==================\nGet smallest eigenvalue.\n")
  #
  # domain is set here (fn_zero functions don't do much)
  p=KOS.Params(fn_zero,fn_zero,vec_d_ball_0_1,alpha)
  #
  # set the parameters for the method
  dict=Dict("max_samples"=>Int32(1e5), "report"=>true,
    #
    "no_levels"=>7, "l0"=>5,  "solve_method"=>"ml",
    "lam_tol"=>lam_tol, # tolerance for eigenvalue solve
    "no_arnoldi_its"=>no_its,
    #
    "relax"=>relax_flag, # for inexact Arnolid
    "B"=>3,#see paper for B (large B smaller WOS tolerance)
   )
  println("Variable accuracy=", dict["relax"],
          ", no_levels=",dict["no_levels"], ", lam_tol=",dict["lam_tol"])
  # set up the mesh
  G=KOSHELPER.get_tri_mesh([-1.,1.,-1.,1.],
                            dict["no_levels"],p.d)
  # first Arnoldi run
  # initial vector for eigenvector
  Y=ones(G.n_int);
  for i=1:G.n_int
    Y[i]=1-(1-p.d(G.get_pt(i))[1])^2
  end
  arnoldi=Construct_Arnoldi(Y)
  #########################################
  lam,evec,delta_lam,arnoldi,el_time,eig_res=wos_ev(p,dict,G,arnoldi,1.)
  dyda_compare(alpha,  lam)
  println("-----------------------------------\n\n")
  #
  println("Variable accuracy=",dict["relax"],
          ", no_levels=",dict["no_levels"], ", lam_tol=",dict["lam_tol"])
  out=G.reshape(evec,KOS.vec_fn_zero)
  Z=G.move_to_TriMesh(out)
  #
  do_plot(G,Z,"evec", "evec.png")
  return lam, el_time, eig_res
end
#######################################
################ end eg_ev.jl
####################################


text:
using CentroidalVoronoiController
using GeometryBasics
using Test

@testset "CentroidalVoronoiController.jl" begin
    # Write your tests here.
    @testset "CentroidalVoronoiController.ControllerTests" begin
        include("test_controller.jl")
    end
end