text:
#========================================================================================#
#	Laboratory 400
#
# Welcome to course 400: Introduction to Quantum Computation!
#
# Author: Niall Palfreyman, 24/04/2022
#========================================================================================#
[
	Activity(
		"""
		Welcome to course 400: Introduction to Quantum Computation!

		In this course we will explore the wonderful world of quantum physics, quantum
		computing and quantum information. In this introductory lab, we will first gather the
		various mathematical tools we need for dealing with quantum systems - in particular,
		we will need to extend our knowledge of Julia to include complex numbers :)

		You might be asking yourself: Why do we need complex numbers? The answer is that quantum
		theory teaches us that particles possess certain probabilities that can interfere with
		each other, and the only way we can calculate these interferences is by using complex
		numbers. So let's start ...

		Recall that a complex number z has the general form x+iy, where x and y are real numbers,
		and i^2==-1. Check this now: In Julia, the number i is written "im", so tell me the
		value of

			im^2
		""",
		"Notice the way in which Julia writes this number - does it make sense to you?",
		x -> x == -1
	),
	Activity(
		"""
		Any complex number has a REAL part x and an IMAGINARY part y. Form the following two
		complex numbers and tell me their sum:

			z1 = 4+5im
			z2 = 3-5im
			z1 + z2
		""",
		"",
		x -> x == 7
	),
	Activity(
		"""
		What is the value of the difference z1 - z2?
		""",
		"",
		x -> x == 1+10im
	),
	Activity(
		"""
		We calculate the CONJUGATE of a complex number by simply reversing the sign of its
		imaginary part. Calculate the complex conjugate of 5+2im:

			conj(5+2im)
		""",
		"",
		x -> x == 5-2im
	),
	Activity(
		"""
		We can also use the operator ' to take the conjugate of a complex number. Calculate
		the value of the complex conjugate of z1:

			z1'

		""",
		"",
		x -> x == 4-5im
	),
	Activity(
		"""
		The conjugate of a product of complex numbers is equal to the product of their
		individual conjugates. Check that this is true:

			(z1*z2)' == a1' * z2'
		""",
		"",
		x -> x == true
	),
	Activity(
		"""
		If we think of a complex number as a vector in the Gaussian plane, its real part is
		its coordinate along the horizontal real axis, and its imaginary part is its coordinate
		along the vertical imaginary axis. In this case, we can calculate the NORM of a complex
		number by using Pythagoras: We square its real and imaginary parts, add them together,
		then take the square-root of the result. Check that this is true by using the Julia
		function abs():

			abs(z1)^2 == real(z1)^2 + imag(z1)^2
		""",
		"",
		x -> x == true
	),
	Activity(
		"""
		We can also compute the norm of a complex number by multiplying the number by its
		conjugate and taking its square-root. Check that this is correct:

			a = abs(z2)
			b = sqrt(z2 * conj(z2))
			a == b
		""",
		"",
		x -> x == true
	),
	Activity(
		"""
		We calculate complex exponents using Euler's famous formula:

			θ = 3
			a = exp(θ*im)
			b = cos(θ) + sin(θ)*im
			a == b
		""",
		"",
		x -> x == true
	),
	Activity(
		"""
		In general, we can express any complex number as the product of its norm and a complex
		exponentiation:

			abs(z1) * exp(angle(z1)*im)

		If you calculate this value, you will see that rounding errors make it not quite equal to
		z1, but it is very close. Approximately how big is the absolute magnitude of the difference
		between these two complex numbers?
		""",
		"Notice we are asking whether z1 == abs(z1) * exp(angle(z1)*im)",
		x -> x < 1e-15
	),
	Activity(
		"""
		In classical computing, the state of a bit in computer memory can only have ONE value:
		EITHER 0 OR 1 (but not both!). The very crucial difference in quantum computing is that
		the state of a bit in a quantum computer memory can SIMULTANEOUSLY have the two different
		values 0 and 1 with certain probabilities. We call this probabilistic state of a quantum
		bit a QUBIT.

		We think of a qubit as an N-dimensional column vector of complex numbers, where N is
		usually a power of 2. We call such a column vector a KET, written mathematically as "|x>".
		In Julia, we will use small letters to identify kets, for example:

			x = [-1, 2im,  1]

		Use Julia to calculate the norm of the ket x by squaring the absolute value of each
		complex component, adding these together, then taking the square root:
		""",
		"sqrt(sum(abs.(x).^2))",
		x -> 2.44 < x < 2.45
	),
	Activity(
		"""
		We know how to use ' to calculate the conjugate of a complex number; now we will see that
		this operator has a more general purpose. Use Julia to check that the size of x is (3,1),
		then find the size of x' :
		""",
		"size(x')",
		x -> x == (1,3)
	),
	Activity(
		"""
		x' is the ADJOINT of x; we calculate it by transposing x (i.e., swapping rows and columns),
		then taking the complex conjugate of each element of x. We can calculate the adjoint of
		any matrix. What is the adjoint of the following matrix?

			M = [
				1+4im 2-3im -2+im
				2-5im 2+3im    5
				  2     3    3-5im
			]
		""",
		"",
		x -> x == [1-4im  2+5im  2+0im; 2+3im  2-3im  3+0im; -2-1im  5+0im  3+5im]
	),
	Activity(
		"""
		We call a row vector of complex numbers a BRA, written mathematically as "<x|". In Julia,
		we will use small letters followed by a bang (!) to identify bras, for example:

			y! = [ 1,   0, im]'

		Use Julia to calculate the norm of the bra y! by squaring the absolute value of each
		complex component, adding these together, then taking the square root:
		""",
		"sqrt(sum(abs.(y!).^2))",
		x -> 1.41 < x < 1.42
	),
	Activity(
		"""
		We can use the adjoint operation to turn a ket into a bra or to turn a bra into a ket.
		Look at the values of x and x', and check that they are respectively a ket and a bra:

			x
			x'
		""",
		"",
		x -> true
	),
	Activity(
		"""
		Verify that taking the adjoint twice returns us to the original ket:

			x'' == x
		""",
		"",
		x -> true
	),
	Activity(
		"""
		Verify that taking the adjoint twice returns us to the original bra:

			y!'' == y!
		""",
		"",
		x -> true
	),
	Activity(
		"""
		We calculate the INNER product (also called the dot product) by matrix-multiplying a
		bra and a ket: <y|.|x> = <y||x> = <y|x> . Paul Dirac suggested the names bra and ket
		because together they form the BRAcKETs <y|x> that we use to calculate probabilities
		in quantum theory. Calculate the following inner product:

			y! * x
		""",
		"",
		x -> x == -1-im
	),
	Activity(
		"""
		Verify that <y|x> is not the same as <x|y>:

			y!*x != x' * y!'
		""",
		"",
		x -> true
	),
	Activity(
		"""
		Do you think the following rule is generally true: <x|y>' == <y|x> ?
		""",
		"This rule is indeed correct - can you verify it using x and y (= y!')?",
		x -> occursin('y',lowercase(x))
	),
	Activity(
		"""
		Two non-zero vectors are ORTHOGONAL if their inner product is zero. For example,
		define the following vector:

			z = [1,im,-1]

		Now verify that z is orthogonal to x:

			<z|x> == <x|z> == 0
		""",
		"z'*x == x'*z == 0",
		x -> true
	),
	Activity(
		"""
		A bra or ket whose inner product with itself is equal to 1.0 is called NORMALISED. In
		quantum computing, qubits are state vectors that represent probability distributions,
		so the sum of the individual probabilities must equal 1. Because of this, we always use
		normalised vectors to represent qubits.

		Verify that this ket is normalised:

			n = [0.6,0.8im]
		""",
		"<n|n> == n'*n == 1.0",
		x -> true
	),
	Activity(
		"""
		We form the OUTER product of two vectors by matrix-multiplying a ket and a bra:

			|x><y| = x*y!

		NOTICE: We calculate an INNER product as (bra*ket), but we calculate the OUTER product
		as (ket*bra). The size of our ket x is (3,1), and the size of our bra y! is (1,3). What
		will be the size of the outer product x*y! ?
		""",
		"Remember the rules for multiplying matrices!",
		x -> x == (3,3)
	),
	Activity(
		"""
		What is the value of the matrix |x><y| ?
		""",
		"x*y!",
		x -> x == [-1  0  im; 2im  0  2; 1  0  -im]
	),
	Activity(
		"""
		A square matrix A is HERMITIAN if A is equal to its own adjoint:

			A == A'

		Verify that the following matrix is Hermitian:

			A = [
				1		3+2im
				3-2im	0
			]
		""",
		"",
		x -> x == -1
	),
	Activity(
		"""
		The special point about Hermitian matrices is that their eigenvalues are always real
		numbers. In quantum theory, eigenvalues are the values that we can measure in physical
		experiments, so it makes sense that we want these to be real numbers. Consequently,
		Hermitian matrices represent physical measurements like: "Is this qubit equal to 1?"

		Is this matrix Hermitian?

			B = [
				-5			99-5im
				99+5im		2
			]
		""",
		"Remember: Transpose B, then take the complex conjugate",
		x -> occursin('y',lowercase(x))
	),
	Activity(
		"""
		A matrix U is UNITARY if the adjoint of U is equal to its own inverse:

			U * U' == U' * U == I

		(where I is the identity matrix). Is the following matrix U unitary?

			U = [0 im;-im 0]
		""",
		"Test whether U'*U==I, and whether U*U'==I",
		x -> occursin('y',lowercase(x))
	),
	Activity(
		"""
		Did you notice that U is not only unitary, but also Hermitian? Verify this now:
		""",
		"",
		x -> true
	),
	Activity(
		"""
		The special point about unitary matrices is that they do not change the norm of a vector.
		They are like rotations in a complex vector space that change the orientation of the
		vector, but not its norm. We will see that this is important for describing how a qubit
		changes over time. Is the following matrix S unitary?

			S = [1 0;0 exp(im)]
		""",
		"",
		x -> occursin('y',lowercase(x))
	),
	Activity(
		"""
		Is the matrix X Hermitian?
		""",
		"",
		x -> occursin('n',lowercase(x))
	),
	Activity(
		"""
		To summarise our findings so far, adjoining turns kets into bras, turns bras into kets,
		and reverses the order of multiplications:

			x! = x'
			y  = y!'
			<y|x>' == <x|y>
			(|x><y|)' == |y>'<x|'
			(A * x)' == x' * A'
			(A * S)' == S' * A'

		Verify each of these rules now ...
		""",
		"Remember that <x|y> simply means x'*y",
		x -> true
	),
	Activity(
		"""
		In general, multiplying a matrix M with a ket x will change the direction of x; however,
		there is often a special set of kets |ψ> (the EIGENVECTORS of M) and complex numbers λ (the
		EIGENVALUES), for which the following equation applies:

			M |ψ> == λ |ψ>

		That is, M only changes the magnitude of |ψ> (multiplies it by λ), but doesn't change the
		direction of |ψ>. For example, here is one of the famous Pauli spin matrices:

			S_y = [0 -im;im 0]

		Verify that the vector

			ψ_yp = [1,im]

		is an eigenvector of S_y, and tell me the corresponding eigenvalue:
		""",
		"Remember, we require: S_y ψ_yp == λ ψ_yp",
		x -> x == 1
	),
	Activity(
		"""
		Now verify that ψ_ym = [1,-im] is also an eigenvector of S_y, and tell me the corresponding
		eigenvalue:
		""",
		"",
		x -> x == -1
	),
	Activity(
		"""
		To end this mathematical introduction, let me remind you that the TRACE of a matrix A is
		sum of the elements along its leading diagonal. What is the trace of the following matrix?

			M = [1 2;3 4]
		""",
		"",
		x -> x == 5
	),
	Activity(
		"""
		The special point about the trace of a matrix is that the trace of an outer product is
		equal to the corresponding inner product. Let's verify this using the two Pauli spin
		eigenkets ψ_yp and ψ_ym. Suppose x is one of these two kets, and y is again either of
		the two (y might even be the same as x), then we can combine these two kets either as an
		inner or an outer product:

			<x|y> == (x' * y) , or
			|x><y| == (x * y')

		The first of these expressions is a number - the INNER product of x and y - while the
		second expression is a matrix - the OUTER product of x and y. Now take the trace of the
		outer product - is it equal to the inner product?
		""",
		"",
		x -> x == occursin('y',lowercase(x))
	),
	Activity(
		"""
		OK, now we have all the mathematical apparatus we need to get started with quantum
		computation! The first thing we will do in the following lab is to build Our Very
		Own Quantum Computer simulator, so that we can get started writing quantum algorithms
		to run on it. What fun! :)
		""",
		"",
		x -> true
	),
]


text:
using Distributed, Dates
t1 = now()
addprocs()

@everywhere begin

    using Dates, Glob, SeisIO, SeisNoise, Serialization

    # directories
    BASEDIR = "/media/four/"
    DATADIR = joinpath(BASEDIR,"LASSO/SAC/")
    FFTDIR = joinpath(BASEDIR,"LASSO/FFT/")
    XMLDIR = joinpath(BASEDIR,"LASSO/XML")
    files = glob("*.SAC",DATADIR)

    if !isdir(FFTDIR)
        mkpath(FFTDIR)
    end

    # parameters
    fs = 250.
    cc_len = 900.
    cc_step = 450.
    freqmin = 0.1
    freqmax = 20.

    # fast saving of FFT
    function serialize_fft(F::FFTData,FFTDIR::String)
        filename = joinpath(FFTDIR,F.name)
        if isa(F.fft,SeisNoise.AbstractGPUArray)
            F = F |> cpu
        end
        serialize(filename,F)
        return nothing
    end

    function seis2fft(
        file::String,
        XMLDIR::String,
        FFTDIR::String,
        fs::Real,
        cc_len::Real,
        cc_step::Real,
        freqmin::Float64,
        freqmax::Float64,
    )
        println("Reading file $(basename(file))")
        S = read_data("sac",file)
        S.fs .= round.(S.fs)
        detrend!(S)
        taper!(S)
        ungap!(S)
        filtfilt!(S,rt="Lowpass",fh=fs/2,np=3)
        resample!(S,fs=fs)

        # get xmlfile
        net,sta,loc,chan = split(S[1].id,'.')
        xmlfile = joinpath(XMLDIR,net*'.'*sta*".xml")
        Resp = read_meta("sxml",xmlfile)
        S.loc[1] = Resp.loc[1]
        translate_resp!(S,Resp[1].resp)
        remove_resp!(S)
        filtfilt!(S,fl=freqmin,fh=freqmax,np=2)
        R = RawData(S[1],cc_len,cc_step)
        demean!(R)
        taper!(R)
        F = rfft(R)
        serialize_fft(F,FFTDIR)
        return nothing
    end
end

t2 = now()
pmap((file -> seis2fft(file,XMLDIR,FFTDIR,fs,cc_len,cc_step,freqmin,freqmax)),files)
t3 = now()

# time for stuff
compute = Dates.canonicalize(Dates.CompoundPeriod(t3-t2))
total = Dates.canonicalize(Dates.CompoundPeriod(t3-t1))
println("Total computation took $total")
println("Computation took $compute")

ext: jl


text:
function n_call!(ys::YASStyle, fst::FST, s::State)
    style = getstyle(ys)

    f = n -> n.typ === PLACEHOLDER || n.typ === NEWLINE

    for (i, n) in enumerate(fst.nodes)
        if i == 3
            # The indent is set here to handle the edge
            # case where the first argument of Call is
            # nestable.
            # ref https://github.com/domluna/JuliaFormatter.jl/issues/387
            fst.indent = s.line_offset
        end

        if n.typ === NEWLINE
            s.line_offset = fst.indent
        elseif n.typ === PLACEHOLDER
            si = findnext(f, fst.nodes, i + 1)
            nest_if_over_margin!(style, fst, s, i; stop_idx = si)
        elseif n.typ === TRAILINGSEMICOLON
            n.val = ""
            n.len = 0
            nest!(style, n, s)
        elseif is_gen(n)
            n.indent = fst.indent
            n.extra_margin = 1
            nest!(style, n, s)
        else
            diff = fst.indent - fst[i].indent
            add_indent!(n, s, diff)
            n.extra_margin = 1
            nest!(style, n, s)
        end
    end
end
@inline n_curly!(ys::YASStyle, fst::FST, s::State) = n_call!(ys, fst, s)
@inline n_ref!(ys::YASStyle, fst::FST, s::State) = n_call!(ys, fst, s)
@inline n_macrocall!(ys::YASStyle, fst::FST, s::State) = n_call!(ys, fst, s)
@inline n_typedcomprehension!(ys::YASStyle, fst::FST, s::State) = n_call!(ys, fst, s)
@inline n_typedvcat!(ys::YASStyle, fst::FST, s::State) = n_call!(ys, fst, s)

function n_tuple!(ys::YASStyle, fst::FST, s::State)
    style = getstyle(ys)
    fst.indent = s.line_offset
    length(fst.nodes) > 0 && is_opener(fst[1]) && (fst.indent += 1)

    f = n -> n.typ === PLACEHOLDER || n.typ === NEWLINE

    for (i, n) in enumerate(fst.nodes)
        if n.typ === NEWLINE
            s.line_offset = fst.indent
        elseif n.typ === PLACEHOLDER
            si = findnext(f, fst.nodes, i + 1)
            nest_if_over_margin!(style, fst, s, i; stop_idx = si)
        elseif n.typ === TRAILINGSEMICOLON
            n.val = ""
            n.len = 0
            nest!(style, n, s)
        elseif is_gen(n)
            n.indent = fst.indent
            n.extra_margin = 1
            nest!(style, n, s)
        else
            diff = fst.indent - fst[i].indent
            add_indent!(n, s, diff)
            n.extra_margin = 1
            nest!(style, n, s)
        end
    end
end
@inline n_braces!(ys::YASStyle, fst::FST, s::State) = n_tuple!(ys, fst, s)
@inline n_vect!(ys::YASStyle, fst::FST, s::State) = n_tuple!(ys, fst, s)
@inline n_parameters!(ys::YASStyle, fst::FST, s::State) = n_tuple!(ys, fst, s)
@inline n_invisbrackets!(ys::YASStyle, fst::FST, s::State) = n_tuple!(ys, fst, s)
@inline n_comprehension!(ys::YASStyle, fst::FST, s::State) = n_tuple!(ys, fst, s)
@inline n_vcat!(ys::YASStyle, fst::FST, s::State) = n_tuple!(ys, fst, s)
@inline n_bracescat!(ys::YASStyle, fst::FST, s::State) = n_tuple!(ys, fst, s)

function n_generator!(ys::YASStyle, fst::FST, s::State)
    style = getstyle(ys)
    diff = s.line_offset - fst[1].indent

    # if the first argument is not a leaf
    # aligns it to be inside the generator
    # expression
    add_indent!(fst[1], s, diff)

    for (i, n) in enumerate(fst.nodes)
        if n.typ === NEWLINE
            s.line_offset = fst.indent
        elseif n.typ === PLACEHOLDER
            si = findnext(n -> n.typ === PLACEHOLDER, fst.nodes, i + 1)
            nest_if_over_margin!(style, fst, s, i; stop_idx = si)
        elseif is_gen(n)
            n.indent = fst.indent
            n.extra_margin = 1
            nest!(style, n, s)
        else
            n.extra_margin = 1
            nest!(style, n, s)
        end
    end
end
@inline n_filter!(ys::YASStyle, fst::FST, s::State) = n_generator!(ys, fst, s)
@inline n_flatten!(ys::YASStyle, fst::FST, s::State) = n_generator!(ys, fst, s)

function n_whereopcall!(ys::YASStyle, fst::FST, s::State)
    style = getstyle(ys)
    fst.indent = s.line_offset
    # after "A where "
    Blen = sum(length.(fst[2:end]))
    fst[1].extra_margin = Blen + fst.extra_margin

    for (i, n) in enumerate(fst.nodes)
        if n.typ === NEWLINE
            s.line_offset = fst.indent
        elseif n.typ === PLACEHOLDER
            si = findnext(n -> n.typ === PLACEHOLDER, fst.nodes, i + 1)
            nest_if_over_margin!(style, fst, s, i; stop_idx = si)
        elseif is_opener(n) && n.val == "{"
            fst.indent = s.line_offset + 1
            nest!(style, n, s)
        elseif i == 1 || i == length(fst.nodes)
            nest!(style, n, s)
        else
            n.extra_margin = 1
            nest!(style, n, s)
        end
    end
end

function n_using!(ys::YASStyle, fst::FST, s::State)
    style = getstyle(ys)
    idx = findfirst(n -> n.val == ":", fst.nodes)
    fst.indent = s.line_offset
    if idx === nothing
        fst.indent += sum(length.(fst[1:2]))
    else
        fst.indent += sum(length.(fst[1:idx+1]))
    end
    for (i, n) in enumerate(fst.nodes)
        if n.typ === PLACEHOLDER
            si = findnext(n -> n.typ === PLACEHOLDER, fst.nodes, i + 1)
            nest_if_over_margin!(style, fst, s, i; stop_idx = si)
        elseif n.typ === NEWLINE
            s.line_offset = fst.indent
        else
            nest!(style, n, s)
        end
    end
end
@inline n_export!(ys::YASStyle, fst::FST, s::State) = n_using!(ys, fst, s)
@inline n_import!(ys::YASStyle, fst::FST, s::State) = n_using!(ys, fst, s)

function n_chainopcall!(ys::YASStyle, fst::FST, s::State)
    style = getstyle(ys)
    n_block!(DefaultStyle(style), fst, s, indent = s.line_offset)
end

function n_comparison!(ys::YASStyle, fst::FST, s::State)
    style = getstyle(ys)
    n_block!(DefaultStyle(style), fst, s, indent = s.line_offset)
end

function n_binaryopcall!(ys::YASStyle, fst::FST, s::State; indent::Int = -1)
    style = getstyle(ys)
    if findfirst(n -> n.typ === PLACEHOLDER, fst.nodes) !== nothing
        n_binaryopcall!(DefaultStyle(style), fst, s; indent = indent)
        return
    end

    start_line_offset = s.line_offset
    walk(increment_line_offset!, fst.nodes[1:end-1], s, fst.indent)
    nest!(style, fst[end], s)
end

function n_for!(ys::YASStyle, fst::FST, s::State)
    style = getstyle(ys)
    block_idx = findfirst(n -> !is_leaf(n), fst.nodes)
    if block_idx === nothing || length(fst[block_idx]) == 0
        nest!(style, fst.nodes, s, fst.indent, extra_margin = fst.extra_margin)
        return
    end

    ph_idx = findfirst(n -> n.typ === PLACEHOLDER, fst[block_idx].nodes)
    for (i, n) in enumerate(fst.nodes)
        if n.typ === NEWLINE && fst.nodes[i+1].typ === Block
            s.line_offset = fst.nodes[i+1].indent
        elseif n.typ === NOTCODE && fst.nodes[i+1].typ === Block
            s.line_offset = fst.nodes[i+1].indent
        elseif n.typ === NEWLINE
            s.line_offset = fst.indent
        else
            n.extra_margin = fst.extra_margin
            if i == 3 && n.typ === Block
                n_block!(style, n, s, indent = s.line_offset)
            else
                nest!(style, n, s)
            end
        end
    end

    # return if the argument block was nested
    ph_idx !== nothing && fst[3][ph_idx].typ === NEWLINE && return

    idx = 5
    n = fst[idx]
    if n.typ === NOTCODE && n.startline == n.endline
        res = get(s.doc.comments, n.startline, (0, ""))
        res == (0, "") && (fst[idx-1] = Whitespace(0))
    end
end


text:
function example_domain_H(parent = RealField(precision(BigFloat)))
    angles =
        fmpq.([
            3 // 2,
            1 // 2,
            1 // 2,
            1 // 2,
            1 // 2,
            3 // 2,
            3 // 2,
            1 // 2,
            1 // 2,
            1 // 2,
            1 // 2,
            3 // 2,
        ])
    vertices = [
        parent.(v) for v in [
            (0, 0),
            (0, 1),
            (-1, 1),
            (-1, -2),
            (0, -2),
            (0, -1),
            (1, -1),
            (1, -2),
            (2, -2),
            (2, 1),
            (1, 1),
            (1, 0),
        ]
    ]

    domain = Polygon(angles, vertices, parent)

    u1 = StandaloneVertexEigenfunction(vertex(domain, 1), fmpq(1 // 2), fmpq(3 // 2))
    u2 = StandaloneVertexEigenfunction(vertex(domain, 6), fmpq(0), fmpq(3 // 2))
    u3 = StandaloneVertexEigenfunction(vertex(domain, 7), fmpq(-1 // 2), fmpq(3 // 2))
    u4 = StandaloneVertexEigenfunction(vertex(domain, 12), fmpq(1), fmpq(3 // 2))

    u = CombinedEigenfunction(
        domain,
        [u1, u2, u3, u4],
        us_to_boundary = [
            setdiff(boundaries(domain), (12, 1)),
            setdiff(boundaries(domain), (5, 6)),
            setdiff(boundaries(domain), (6, 7)),
            setdiff(boundaries(domain), (11, 12)),
        ],
    )

    return domain, u
end

function example_domain_ngon(
    n,
    parent = RealField(precision(BigFloat));
    lightning = false,
    linked = false,
)
    θ = fmpq((n - 2) // n)
    angles = fill(θ, n)
    vertices = [(cos(θ), sin(θ)) for θ in (2parent(π) / n) .* (0:n-1)]
    domain = Polygon(angles, vertices, parent)

    if !lightning
        if linked
            us = [
                LinkedEigenfunction(
                    [
                        StandaloneVertexEigenfunction(
                            vertex(domain, i),
                            i * (1 - θ) + θ * 1 // 2,
                            θ,
                        ) for i in boundaries(domain)
                    ],
                    excluded_boundaries = [
                        BitSet([mod1(i - 1, n), i]) for i in boundaries(domain)
                    ],
                ),
            ]
            us_to_boundary = fill(BitSet([1]), length(us))
        else
            us = [
                StandaloneVertexEigenfunction(
                    vertex(domain, i),
                    i * (1 - θ) + θ * 1 // 2,
                    θ,
                ) for i in boundaries(domain)
            ]
            us_to_boundary = [
                setdiff(boundaries(domain), (i, mod1(i - 1, length(boundaries(domain))))) for i in eachindex(us)
            ]
        end
    else
        if linked
            us = [
                LinkedEigenfunction(
                    [
                        StandaloneLightningEigenfunction(
                            vertex(domain, i),
                            i * (1 - θ) + θ * 1 // 2,
                            θ,
                        ) for i in boundaries(domain)
                    ],
                    ones(n),
                ),
                StandaloneInteriorEigenfunction(domain, stride = n),
            ]

            us_to_boundary = fill(BitSet([1]), length(us))
        else
            us = [
                [
                    StandaloneLightningEigenfunction(
                        vertex(domain, i),
                        i * (1 - θ) + θ * 1 // 2,
                        θ,
                    ) for i in boundaries(domain)
                ]
                StandaloneInteriorEigenfunction(domain)
            ]

            us_to_boundary = fill(boundaries(domain), length(us))
        end
    end

    u = CombinedEigenfunction(domain, us, us_to_boundary = us_to_boundary)

    return domain, u
end

function example_domain_triangle_in_triangle(parent = RealField(precision(BigFloat)))
    domain1 = Triangle(fmpq(1 // 3), fmpq(1 // 3), parent)
    domain2 = TransformedDomain(
        Triangle(fmpq(1 // 3), fmpq(1 // 3), parent),
        fmpq(0),
        parent(0.25),
        SVector(cospi(fmpq(1 // 6), parent), sinpi(fmpq(1 // 6), parent)) / 4,
    )

    domain = IntersectedDomain(domain1, domain2)

    u1 = StandaloneLightningEigenfunction(domain1, 1)

    u2 = LinkedEigenfunction([StandaloneLightningEigenfunction(domain1, i) for i = 2:3])

    u3 = StandaloneLightningEigenfunction(domain2, 1, l = parent(0.15), outside = true)

    u4 = LinkedEigenfunction([
        StandaloneLightningEigenfunction(domain2, i, l = parent(0.15), outside = true)
        for i = 2:3
    ])

    u5 = StandaloneInteriorEigenfunction(domain)

    us = [u1, u2, u3, u4, u5]
    #us = vcat([u1], u2.us, [u3], u4.us, [u5])
    orders = [1, 1, 1, 1, 1]
    us_to_boundary = fill(1:6, length(us))

    u = CombinedEigenfunction(domain, us, orders, us_to_boundary = us_to_boundary)

    return domain, u
end

function example_domain_ngon_in_ngon(
    n1,
    n2,
    parent = RealField(precision(BigFloat));
    linked = true,
    even = true,
    T = arb,
)
    θ1 = fmpq((n1 - 2) // n1)
    angles1 = fill(θ1, n1)
    vertices1 = [(cos(θ), sin(θ)) for θ in (2parent(π) / n1) .* (0:n1-1)]
    domain1 = Polygon(angles1, vertices1, parent)

    θ2 = fmpq((n2 - 2) // n2)
    angles2 = fill(θ2, n2)
    vertices2 = [(cos(θ), sin(θ)) for θ in (2parent(π) / n2) .* (0:n2-1)]
    domain2 = TransformedDomain(
        Polygon(angles2, vertices2, parent),
        fmpq(0),
        parent(0.5),
        SVector(parent(0), parent(0)),
    )

    domain = IntersectedDomain(domain1, domain2)

    if linked
        us = [
            LinkedEigenfunction([
                StandaloneLightningEigenfunction{T,arb}(domain1, i; even) for
                i in boundaries(domain1)
            ]),
            LinkedEigenfunction([
                StandaloneLightningEigenfunction{T,arb}(
                    domain2,
                    i,
                    outside = true,
                    l = parent(0.4);
                    even,
                ) for i in boundaries(domain2)
            ]),
            StandaloneInteriorEigenfunction(domain, stride = 4; even),
        ]

        orders = ifelse(even, [2, 2, 1], [3, 3, 2])
        even_boundaries = ifelse(even, [1, 5], [])

        u = CombinedEigenfunction(
            domain,
            us,
            orders,
            us_to_boundary = fill(BitSet([1, 5]), length(us));
            even_boundaries,
        )
    else
        us = vcat(
            [
                StandaloneLightningEigenfunction{T,arb}(domain1, i) for
                i in boundaries(domain1)
            ],
            [
                StandaloneLightningEigenfunction{T,arb}(
                    domain2,
                    i,
                    outside = true,
                    l = parent(0.4),
                ) for i in boundaries(domain2)
            ],
            [StandaloneInteriorEigenfunction(domain)],
        )

        u = CombinedEigenfunction(domain, us, [3, 3, 3, 3, 3, 3, 3, 3, 2])
    end

    return domain, u
end

"""
    example_domain_goal_v1(N, d, h, parent; ...)

Return the domain given by a hexagon with 6 triangles cut out as well
as a corresponding eigenfunction.

The placement of the holes is determined by `(N, d, h)`.

The keyword `symmetry_class` determines which symmetries are assumed
by the eigenfunction. `symmetry_class = 1` assumes it's six-fold
symmetric and even, this is satisfied by the first and second
eigenvalue. `symmetry_class = 2` assumes it's even w.r.t. the x-axis
and odd w.r.t. the y-axis, `symmetry_class = 3` does the opposite,
this is used by the third and fourth eigenvalue.

Computing the first few eigenvalues in Matlab using a finite element
method gives (with `(N, d, h) = (22, 9, 6)`)
```
   31.5943
   75.8995
   75.9633
   75.9633
```
though it's likely not accurate to more than 1 or 2 decimals. If `(N,
d, h) = (27, 11, 6)` then they are
```
   31.0432
   63.2104
   63.7259
   63.7259
   68.2629
   68.2629
```
For the second eigenvalue the value `63.2083359862688427940` might be
slightly more accurate.
"""
function example_domain_goal_v1(
    N::Integer = 27,
    d::Integer = 11,
    h::Integer = 6,
    parent = RealField(precision(BigFloat));
    T = arb,
    lightning = false,
    inner_expansion = true,
    outer_expansion = false,
    symmetry_class = 1,
    even = true,
    reversed = true,
)
    ### Compute the domain ###
    # The main domain is a hexagon
    n = 6
    θ = fmpq((n - 2) // n)
    angles = fill(θ, n)
    vertices = [(cospi(θ, parent), sinpi(θ, parent)) for θ in fmpq(2 // n) .* (0:n-1)]
    exterior = Polygon(angles, vertices, parent)

    # The interior domains are triangles
    @assert mod(h, 3) == 0
    points = [
        SVector(parent(d // N), -parent(h // N) / sqrt(parent(3))),
        SVector(parent((d + h) // N), parent(0)),
        SVector(parent(d // N), parent(h // N) / sqrt(parent(3))),
    ]

    interior_angles = [1 // 3, 1 // 3, 1 // 3]

    interiors = [
        TransformedDomain(Polygon(interior_angles, points, parent), i // 3, 1, [0, 0])
        for i = 0:5
    ]

    domain = IntersectedDomain(exterior, interiors)

    ### Compute the eigenfunction ###
    # Expansions from the vertices of the hexagon
    if symmetry_class == 1
        if lightning
            outer_eigenfunction =
                i -> StandaloneLightningEigenfunction(
                    vertex(exterior, i),
                    fmpq(mod(Rational(1 - θ // 2 + (i - 1) * (1 - θ)), 2)),
                    θ;
                    even,
                )
            outer_excluded_boundary = i -> BitSet()
            order1 = ifelse(even, 2, 3)

        else
            outer_eigenfunction =
                i -> StandaloneVertexEigenfunction(
                    vertex(exterior, i),
                    i * (1 - θ) + θ * 1 // 2,
                    θ,
                    stride = ifelse(even, 2, 1),
                )
            outer_excluded_boundary =
                i -> BitSet([mod1(i - 1, length(boundaries(exterior))), i])
            order1 = ifelse(even, 1, 2)
        end

        u1 = LinkedEigenfunction(
            [outer_eigenfunction(i) for i in boundaries(exterior)],
            excluded_boundaries = [
                outer_excluded_boundary(i) for i in boundaries(exterior)
            ],
        )

        # Expansions from the outer tip of the triangles
        u2 = LinkedEigenfunction([
            StandaloneLightningEigenfunction{T,fmpq}(
                vertex(d, 2),
                7 // 6 + d.rotation,
                2 - d.original.angles[2],
                l = parent(h // 2N),
                even = even,
            ) for d in interiors
        ])
        if even && reversed
            order2 = 3 * 2
        elseif even
            order2 = 2
        else
            order2 = 3
        end

        # Expansions from the inner tips of the triangles
        u3 = LinkedEigenfunction(
            [
                StandaloneLightningEigenfunction{T,fmpq}(
                    vertex(d, i),
                    ifelse(i == 1, 1 // 2, -1 // 6) + d.rotation,
                    2 - d.original.angles[i],
                    l = parent(h // 2N),
                    even = even && !reversed,
                    reversed = reversed && i == 3,
                ) for d in interiors, i in [1, 3]
            ][:],
        )
        if even && reversed
            order3 = 3 * 3
        elseif even
            order3 = 2
        else
            order3 = 3
        end

        # Expansion from the center
        u4 = StandaloneInteriorEigenfunction(SVector(domain.parent(0), domain.parent(0)), stride = 6; even)
        order4 = ifelse(even, 1, 2)

        # Interior expansions closer to outer boundaries
        u5 = LinkedEigenfunction([
            StandaloneInteriorEigenfunction(
                parent(3 // 4) .*
                SVector(cospi(fmpq(i // 3), parent), sinpi(fmpq(i // 3), parent)),
                fmpq(i // 3),
                stride = 6,
                even = true,
            ) for i = 0:5
        ])
        order5 = ifelse(even, 1, 2)

        us = AbstractPlanarEigenfunction[u1, u2, u3]
        orders = Int[order1, order2, order3]

        if inner_expansion
            push!(us, u4)
            push!(orders, order4)
        end
        if outer_expansion
            push!(us, u5)
            push!(orders, order5)
        end

        us_to_boundary = fill(BitSet([1, 8, 9]), length(us))
        even_boundaries = ifelse(even, Int[1, 9], Int[])
    elseif symmetry_class == 2
        ####################################################################
        is_even = i -> i == 1 || i == 4
        is_reversed = i -> i == 3 || i == 6

        if lightning
            outer_eigenfunction =
                i -> StandaloneLightningEigenfunction(
                    vertex(exterior, i),
                    fmpq(mod(Rational(1 - θ // 2 + (i - 1) * (1 - θ)), 2)),
                    θ;
                    even = is_even(i),
                    reversed = is_reversed(i),
                )
            outer_excluded_boundary = i -> BitSet()
            order11 = 2
            order12 = 3
        else
            outer_eigenfunction =
                i -> StandaloneVertexEigenfunction(
                    vertex(exterior, i),
                    i * (1 - θ) + θ * 1 // 2,
                    θ,
                    stride = ifelse(is_even(i), 2, 1),
                    reversed = is_reversed(i),
                )
            outer_excluded_boundary =
                i -> BitSet([mod1(i - 1, length(boundaries(exterior))), i])
            order11 = 1
            order12 = 1
        end

        u11 = LinkedEigenfunction(
            [outer_eigenfunction(i) for i in (1, 4)],
            [parent(1), parent(-1)],
            excluded_boundaries = [outer_excluded_boundary(i) for i in (1, 4)],
        )

        u12 = LinkedEigenfunction(
            [outer_eigenfunction(i) for i in (2, 3, 5, 6)],
            [parent(1), parent(-1), parent(-1), parent(1)],
            excluded_boundaries = [outer_excluded_boundary(i) for i in (2, 3, 5, 6)],
        )

        # Expansions from the outer tip of the triangles
        u21 = LinkedEigenfunction(
            [
                StandaloneLightningEigenfunction{T,fmpq}(
                    vertex(d, 2),
                    7 // 6 + d.rotation,
                    2 - d.original.angles[2],
                    l = parent(h // 2N),
                    even = is_even(i),
                ) for (i, d) in ((1, interiors[1]), (4, interiors[4]))
            ],
            [parent(1), parent(-1)],
        )
        order21 = 2

        u22 = LinkedEigenfunction(
            [
                StandaloneLightningEigenfunction{T,fmpq}(
                    vertex(d, 2),
                    7 // 6 + d.rotation,
                    2 - d.original.angles[2],
                    l = parent(h // 2N),
                    even = is_even(i),
                    reversed = is_reversed(i),
                ) for (i, d) in (
                    (2, interiors[2]),
                    (3, interiors[3]),
                    (5, interiors[5]),
                    (6, interiors[6]),
                )
            ],
            [parent(1), parent(-1), parent(-1), parent(1)],
        )
        order22 = 3

        # Expansions from the inner tips of the triangles
        u31 = LinkedEigenfunction(
            [
                StandaloneLightningEigenfunction{T,fmpq}(
                    vertex(d, i),
                    ifelse(i == 1, 1 // 2, -1 // 6) + d.rotation,
                    2 - d.original.angles[i],
                    l = parent(h // 2N),
                    reversed = i == 3,
                ) for (j, d) in ((1, interiors[1]), (4, interiors[4])), i in [1, 3]
            ][:],
            [parent(1), parent(-1), parent(1), parent(-1)],
        )
        order31 = 3

        u32 = LinkedEigenfunction(
            [
                StandaloneLightningEigenfunction{T,fmpq}(
                    vertex(d, i),
                    ifelse(i == 1, 1 // 2, -1 // 6) + d.rotation,
                    2 - d.original.angles[i],
                    l = parent(h // 2N),
                    reversed = i == 3,
                ) for (j, d, i) in (
                    (2, interiors[2], 1),
                    (3, interiors[3], 3),
                    (5, interiors[5], 1),
                    (6, interiors[6], 3),
                )
            ],
            [parent(1), parent(-1), parent(-1), parent(1)],
        )
        order32 = 3

        u33 = LinkedEigenfunction(
            [
                StandaloneLightningEigenfunction{T,fmpq}(
                    vertex(d, i),
                    ifelse(i == 1, 1 // 2, -1 // 6) + d.rotation,
                    2 - d.original.angles[i],
                    l = parent(h // 2N),
                    reversed = i == 3,
                ) for (j, d, i) in (
                    (2, interiors[2], 3),
                    (3, interiors[3], 1),
                    (5, interiors[5], 3),
                    (6, interiors[6], 1),
                )
            ],
            [parent(1), parent(-1), parent(-1), parent(1)],
        )
        order33 = 3

        # Expansion from the center
        u4 = StandaloneInteriorEigenfunction(domain, stride = 2, offset = 1; even)
        order4 = 1

        us = AbstractPlanarEigenfunction[u11, u12, u21, u22, u31, u32, u33, u4]
        orders =
            Int[order11, order12, 3order21, 3order22, 3order31, 3order32, 3order33, order4]

        us_to_boundary = fill(BitSet([1, 2, 8, 9, 10, 11, 12]), length(us))
        even_boundaries = Int[2, 9]
    elseif symmetry_class == 3
        ####################################################################
        is_odd = i -> i == 1 || i == 4
        is_reversed = i -> i == 3 || i == 4 || i == 6

        if lightning
            outer_eigenfunction =
                i -> StandaloneLightningEigenfunction(
                    vertex(exterior, i),
                    fmpq(mod(Rational(1 - θ // 2 + (i - 1) * (1 - θ)), 2)),
                    θ;
                    odd = is_odd(i),
                    reversed = is_reversed(i),
                )
            outer_excluded_boundary = i -> BitSet()
            order11 = 1
            order12 = 3
        else
            outer_eigenfunction =
                i -> StandaloneVertexEigenfunction(
                    vertex(exterior, i),
                    i * (1 - θ) + θ * 1 // 2,
                    θ,
                    stride = ifelse(is_odd(i), 2, 1),
                    offset = ifelse(is_odd(i), 1, 0),
                    reversed = is_reversed(i),
                )
            outer_excluded_boundary =
                i -> BitSet([mod1(i - 1, length(boundaries(exterior))), i])
            order11 = 1
            order12 = 1
        end

        u11 = LinkedEigenfunction(
            [outer_eigenfunction(i) for i in (1, 4)],
            [parent(1), parent(1)],
            excluded_boundaries = [outer_excluded_boundary(i) for i in (1, 4)],
        )

        u12 = LinkedEigenfunction(
            [outer_eigenfunction(i) for i in (2, 3, 5, 6)],
            [parent(1), parent(1), parent(-1), parent(-1)],
            excluded_boundaries = [outer_excluded_boundary(i) for i in (2, 3, 5, 6)],
        )

        # Expansions from the outer tip of the triangles
        u21 = LinkedEigenfunction(
            [
                StandaloneLightningEigenfunction{T,fmpq}(
                    vertex(d, 2),
                    7 // 6 + d.rotation,
                    2 - d.original.angles[2],
                    l = parent(h // 2N),
                    odd = is_odd(i),
                    reversed = i == 4,
                ) for (i, d) in ((1, interiors[1]), (4, interiors[4]))
            ],
            [parent(1), parent(1)],
        )
        order21 = 1

        u22 = LinkedEigenfunction(
            [
                StandaloneLightningEigenfunction{T,fmpq}(
                    vertex(d, 2),
                    7 // 6 + d.rotation,
                    2 - d.original.angles[2],
                    l = parent(h // 2N),
                    odd = is_odd(i),
                    reversed = is_reversed(i),
                ) for (i, d) in (
                    (2, interiors[2]),
                    (3, interiors[3]),
                    (5, interiors[5]),
                    (6, interiors[6]),
                )
            ],
            [parent(1), parent(1), parent(-1), parent(-1)],
        )
        order22 = 3

        # Expansions from the inner tips of the triangles
        u31 = LinkedEigenfunction(
            [
                StandaloneLightningEigenfunction{T,fmpq}(
                    vertex(d, i),
                    ifelse(i == 1, 1 // 2, -1 // 6) + d.rotation,
                    2 - d.original.angles[i],
                    l = parent(h // 2N),
                    reversed = i == 3,
                ) for (j, d) in ((1, interiors[1]), (4, interiors[4])), i in [1, 3]
            ][:],
            [parent(1), parent(-1), parent(-1), parent(1)],
        )
        order31 = 3

        u32 = LinkedEigenfunction(
            [
                StandaloneLightningEigenfunction{T,fmpq}(
                    vertex(d, i),
                    ifelse(i == 1, 1 // 2, -1 // 6) + d.rotation,
                    2 - d.original.angles[i],
                    l = parent(h // 2N),
                    reversed = i == 3,
                ) for (j, d, i) in (
                    (2, interiors[2], 1),
                    (3, interiors[3], 3),
                    (5, interiors[5], 1),
                    (6, interiors[6], 3),
                )
            ],
            [parent(1), parent(1), parent(-1), parent(-1)],
        )
        order32 = 3

        u33 = LinkedEigenfunction(
            [
                StandaloneLightningEigenfunction{T,fmpq}(
                    vertex(d, i),
                    ifelse(i == 1, 1 // 2, -1 // 6) + d.rotation,
                    2 - d.original.angles[i],
                    l = parent(h // 2N),
                    reversed = i == 3,
                ) for (j, d, i) in (
                    (2, interiors[2], 3),
                    (3, interiors[3], 1),
                    (5, interiors[5], 3),
                    (6, interiors[6], 1),
                )
            ],
            [parent(1), parent(1), parent(-1), parent(-1)],
        )
        order33 = 3

        # Expansion from the center
        u4 = StandaloneInteriorEigenfunction(domain, stride = 2, odd = true)
        order4 = 1

        us = AbstractPlanarEigenfunction[u11, u12, u21, u22, u31, u32, u33, u4]

        orders =
            Int[order11, order12, 3order21, 3order22, 3order31, 3order32, 3order33, order4]

        us_to_boundary = fill(BitSet([1, 2, 8, 9, 10, 11, 12]), length(us))
        even_boundaries = Int[2, 9]
    else
        throw(ArgumentError("symmetry_class should be 1 or 2, got $symmetry_class"))
    end

    u = CombinedEigenfunction(domain, us, orders; us_to_boundary, even_boundaries)

    return domain, u
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    CoDaArray(table)

Construct an array of compositional data from `table`.
"""
struct CoDaArray{D,PARTS} <: AbstractVector{Composition{D,PARTS}}
  data::Matrix{Union{Float64,Missing}}
end

function CoDaArray(table)
  PARTS = Tables.columnnames(table)
  data  = Tables.matrix(table, transpose=true)
  CoDaArray{length(PARTS),Tuple(PARTS)}(data)
end

Base.getindex(array::CoDaArray{D,PARTS}, ind) where {D,PARTS} =
  Composition(PARTS, getfield(array, :data)[:,ind])

Base.size(array::CoDaArray) = (size(getfield(array, :data), 2),)

Base.IndexStyle(::Type{<:CoDaArray}) = IndexLinear()

function Base.getproperty(array::CoDaArray{D,PARTS}, part::Symbol) where {D,PARTS}
  ind = findfirst(isequal(part), PARTS)
  if isnothing(ind)
    throw(KeyError("$part"))
  else
    vec(getfield(array, :data)[ind,:])
  end
end

"""
    parts(array)

Parts in compositional `array`.
"""
parts(::CoDaArray{D,PARTS}) where {D,PARTS} = PARTS

"""
    compose(table, cols; keepcols=true, as=:coda)

Convert columns `cols` of `table` into parts of a
composition and save the result in a [`CoDaArray`](@ref).
If `keepcols` is set to `true`, then save the result `as`
a column in a new table with all other columns preserved.
"""
function compose(table, cols=Tables.columnnames(table);
                 keepcols=true, as=:coda)
  # construct compositional array from selected columns
  coda = table |> Select(cols) |> CoDaArray

  # different types of return
  if keepcols
    other = setdiff(Tables.columnnames(table), cols)
    osel  = table |> Select(other)
    ocol  = [o => Tables.getcolumn(osel, o) for o in other]
    # preserve input table type
    𝒯 = Tables.materializer(table)
    𝒯((; ocol..., as => coda))
  else
    coda
  end
end

# -----------------
# TABLES INTERFACE
# -----------------

# implement table interface for CoDaArray
Tables.istable(::Type{<:CoDaArray}) = true
Tables.rowaccess(::Type{<:CoDaArray}) = true
Tables.rows(array::CoDaArray) = array

# implement row interface for Composition
Tables.getcolumn(c::Composition, i::Int)    = getfield(c, :data)[i]
Tables.getcolumn(c::Composition, n::Symbol) = getfield(c, :data)[n]
Tables.columnnames(c::Composition{D,PARTS}) where {D,PARTS} = PARTS


text:
"""
convert tensor to matrix or matrix to tensor
storing the index info
"""
module TensorMatrices_lemon

	export AbstractTenmat, Tenmat, tensorSize, size, getindex, setindex!, similar, ndims, tensor2tenmat, tenmat2tensor, getMatrix

	include("tenmat.jl")

end


text:
# A lightly-modified version of the same function in Base
# Highlights argument types with color specified by highlighter(typ)
function show_tuple_as_call(@nospecialize(highlighter), io::IO, name::Symbol, sig::Type, demangle=false, kwargs=nothing)
    if sig === Tuple
        print(io, demangle ? Base.demangle_function_name(name) : name, "(...)")
        return
    end
    tv = Any[]
    env_io = io
    while isa(sig, UnionAll)
        push!(tv, sig.var)
        env_io = IOContext(env_io, :unionall_env => sig.var)
        sig = sig.body
    end
    sig = (sig::DataType).parameters

    ft = sig[1]
    uw = Base.unwrap_unionall(ft)
    if ft <: Function && isa(uw,DataType) && isempty(uw.parameters) &&
            isdefined(uw.name.module, uw.name.mt.name) &&
            ft == typeof(getfield(uw.name.module, uw.name.mt.name))
        print(env_io, (demangle ? demangle_function_name : identity)(uw.name.mt.name))
    elseif isa(ft, DataType) && ft.name === Type.body.name && !Core.Compiler.has_free_typevars(ft)
        f = ft.parameters[1]
        print(env_io, f)
    else
        print(env_io, "(::", ft, ")")
    end

    first = true
    print(io, "(")
    for i = 2:length(sig)  # fixme (iter): `eachindex` with offset?
        first || print(io, ", ")
        first = false
        printstyled(env_io, "::", sig[i], color=highlighter(sig[i]))
    end
    if kwargs !== nothing
        print(io, "; ")
        first = true
        for (k, t) in kwargs
            first || print(io, ", ")
            first = false
            print(io, k, "::")
            show(io, t)
        end
    end
    print(io, ")")
    Base.show_method_params(io, tv)
    nothing
end

function stripType(@nospecialize(typ))
    if isa(typ, UnionAll)
        typ = Base.unwrap_unionall(typ)
    elseif isa(typ, TypeVar) || isa(typ, Union)
        return typ
    end
    return typ <: Type && length(typ.parameters) == 1 ? typ.parameters[1] : typ
end
nonconcrete_red(@nospecialize(typ)) = isconcretetype(stripType(typ)) ? :nothing : :red

const _emptybackedges = MethodInstance[]
# Extension API
backedges(mi::MethodInstance) = isdefined(mi, :backedges) ? mi.backedges : _emptybackedges
method(mi::MethodInstance) = mi.def
specTypes(mi::MethodInstance) = mi.specTypes
instance(mi::MethodInstance) = mi

function callstring(io, mi)
    show_tuple_as_call(nonconcrete_red, IOContext(io, :color=>true), method(mi).name, specTypes(mi))
    return String(take!(io))
end

if has_treemenu
    struct Data{T}
        callstr::String
        nd::T
    end

    function treelist(mi)
        io = IOBuffer()
        str = callstring(io, mi)
        treelist!(Node(Data(str, mi)), io, mi, "", Base.IdSet{typeof(mi)}())
    end
    function treelist!(parent::Node, io, mi, indent::AbstractString, visited::Base.IdSet)
        mi ∈ visited && return parent
        push!(visited, mi)
        indent *= " "
        for edge in backedges(mi)
            str = indent * callstring(io, edge)
            child = Node(Data(str, edge), parent)
            treelist!(child, io, edge, indent, visited)
        end
        return parent
    end
else
    # TreeMenu can't be implemented, fallback to non-folding menu
    treelist(mi) = treelist!(String[], typeof(mi)[], IOBuffer(), mi, "", Base.IdSet{typeof(mi)}())

    function treelist!(strs, mis, io::IO, mi, indent::AbstractString, visited::Base.IdSet)
        mi ∈ visited && return strs, mis
        push!(visited, mi)
        str = indent * callstring(io, mi)
        push!(strs, str)
        push!(mis, mi)
        indent *= " "
        for edge in backedges(mi)
            treelist!(strs, mis, io, edge, indent, visited)
        end
        return strs, mis
    end
end


text:
"""
    AbstractKKTVector{T, VT}

Supertype for KKT's right-hand-side vectors ``(x, s, y, z, ν, w)``.
"""
abstract type AbstractKKTVector{T, VT} end

"""
    full(X::AbstractKKTVector)

Return the all the values stored inside the KKT vector `X`.
"""
function full end

Base.length(rhs::AbstractKKTVector) = length(full(rhs))

"""
    number_primal(X::AbstractKKTVector)

Get total number of primal values ``(x, s)`` in KKT vector `X`.
"""
number_primal(rhs::AbstractKKTVector) = length(primal(rhs))

"""
    number_dual(X::AbstractKKTVector)

Get total number of dual values ``(y, z)`` in KKT vector `X`.
"""
number_dual(rhs::AbstractKKTVector) = length(dual(rhs))

"""
    primal(X::AbstractKKTVector)

Return the primal values ``(x, s)`` stored in the KKT vector `X`.
"""
function primal end

"""
    dual(X::AbstractKKTVector)

Return the dual values ``(y, z)`` stored in the KKT vector `X`.
"""
function dual end

"""
    primal_dual(X::AbstractKKTVector)

Return both the primal and the dual values ``(x, s, y, z)`` stored in the KKT vector `X`.
"""
function primal_dual end

"""
    dual_lb(X::AbstractKKTVector)

Return the dual values ``ν`` associated to the lower-bound
stored in the KKT vector `X`.
"""
function dual_lb end

"""
    dual_ub(X::AbstractKKTVector)

Return the dual values ``w`` associated to the upper-bound
stored in the KKT vector `X`.
"""
function dual_ub end

function Base.fill!(rhs::AbstractKKTVector{T}, val::T) where T
    fill!(full(rhs), val)
end

# Overload basic BLAS operations.
norm(X::AbstractKKTVector, p::Real) = norm(full(X), p)
dot(X::AbstractKKTVector, Y::AbstractKKTVector) = dot(full(X), full(Y))
function axpy!(a::Number, X::AbstractKKTVector, Y::AbstractKKTVector)
    axpy!(a, full(X), full(Y))
end

#=
    ReducedKKTVector
=#

"""
    ReducedKKTVector{T, VT<:AbstractVector{T}} <: AbstractKKTVector{T, VT}

KKT vector ``(x, s, y, z)``, associated to a [`AbstractReducedKKTSystem`](@ref).

Compared to [`UnreducedKKTVector`](@ref), it does not store
the dual values associated to the primal's lower and upper bounds.
"""
struct ReducedKKTVector{T, VT<:AbstractVector{T}} <: AbstractKKTVector{T, VT}
    values::VT
    xp::VT # unsafe view
    xl::VT # unsafe view
end

ReducedKKTVector(n::Int, m::Int, nlb::Int, nub::Int) = ReducedKKTVector(n, m)
function ReducedKKTVector(n::Int, m::Int)
    x = Vector{Float64}(undef, n + m)
    fill!(x, 0.0)
    # Wrap directly array x to avoid dealing with views
    pp = pointer(x)
    xp = unsafe_wrap(Vector{Float64}, pp, n)
    pd = pointer(x, n + 1)
    xl = unsafe_wrap(Vector{Float64}, pd, m)
    return ReducedKKTVector{Float64, Vector{Float64}}(x, xp, xl)
end
function ReducedKKTVector(rhs::AbstractKKTVector)
    return ReducedKKTVector(number_primal(rhs), number_dual(rhs))
end

full(rhs::ReducedKKTVector) = rhs.values
primal(rhs::ReducedKKTVector) = rhs.xp
dual(rhs::ReducedKKTVector) = rhs.xl
primal_dual(rhs::ReducedKKTVector) = rhs.values

#=
    UnreducedKKTVector
=#

"""
    UnreducedKKTVector{T, VT<:AbstractVector{T}} <: AbstractKKTVector{T, VT}

Full KKT vector ``(x, s, y, z, ν, w)``, associated to a [`AbstractUnreducedKKTSystem`](@ref).

"""
struct UnreducedKKTVector{T, VT<:AbstractVector{T}} <: AbstractKKTVector{T, VT}
    values::VT
    x::VT  # unsafe view
    xp::VT # unsafe view
    xl::VT # unsafe view
    xzl::VT # unsafe view
    xzu::VT # unsafe view
end

function UnreducedKKTVector(n::Int, m::Int, nlb::Int, nub::Int)
    values = Vector{Float64}(undef, n + m + nlb + nub)
    fill!(values, 0.0)
    # Wrap directly array x to avoid dealing with views
    pp = pointer(values)
    x = unsafe_wrap(Vector{Float64}, pp, n + m) # Primal-Dual
    xp = unsafe_wrap(Vector{Float64}, pp, n) # Primal
    pd = pointer(values, n + 1)
    xl = unsafe_wrap(Vector{Float64}, pd, m) # Dual
    pzl = pointer(values, n + m + 1)
    xzl = unsafe_wrap(Vector{Float64}, pzl, nlb) # Lower bound
    pzu = pointer(values, n + m + nlb + 1)
    xzu = unsafe_wrap(Vector{Float64}, pzu, nub) # Upper bound
    return UnreducedKKTVector{Float64, Vector{Float64}}(values, x, xp, xl, xzl, xzu)
end

full(rhs::UnreducedKKTVector) = rhs.values
primal(rhs::UnreducedKKTVector) = rhs.xp
dual(rhs::UnreducedKKTVector) = rhs.xl
primal_dual(rhs::UnreducedKKTVector) = rhs.x
dual_lb(rhs::UnreducedKKTVector) = rhs.xzl
dual_ub(rhs::UnreducedKKTVector) = rhs.xzu


text:
mutable struct Tree <: AbstractAgent
    id::Int
    pos::Dims{2}
    status::Symbol  #:green, :burning, :burnt
end

"""
``` julia
forest_fire(;
    density = 0.8,
    griddims = (100, 100)
)
```
Same as in [Forest fire model](@ref).
"""
function forest_fire(; density = 0.7, griddims = (100, 100))
    space = GridSpace(griddims; periodic = false, metric = :euclidean)
    forest = AgentBasedModel(Tree, space)
    for position in positions(forest)
        if rand(forest.rng) < density
            state = position[1] == 1 ? :burning : :green
            add_agent!(position, forest, state)
        end
    end
    return forest, forest_agent_step!, dummystep
end

function forest_agent_step!(tree, forest)
    if tree.status == :burning
        for neighbor in nearby_agents(tree, forest)
            if neighbor.status == :green
                neighbor.status = :burning
            end
        end
        tree.status = :burnt
    end
end


text:
# # 2D decaying turbulence
#
#md # This example can be viewed as a Jupyter notebook via [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/literated/twodnavierstokes_decaying.ipynb).
#
# A simulation of decaying two-dimensional turbulence.
#
# ## Install dependencies
#
# First let's make sure we have all required packages installed.

# ```julia
# using Pkg
# pkg"add GeophysicalFlows, Printf, Random, Plots"
# ```

# ## Let's begin
# Let's load `GeophysicalFlows.jl` and some other needed packages.
#
using GeophysicalFlows, Printf, Random, Plots

using Random: seed!
using GeophysicalFlows: peakedisotropicspectrum

# ## Choosing a device: CPU or GPU

dev = CPU()     # Device (CPU/GPU)
nothing # hide

# ## Numerical, domain, and simulation parameters
#
# First, we pick some numerical and physical parameters for our model.

n, L  = 128, 2π             # grid resolution and domain length
nothing # hide

# Then we pick the time-stepper parameters
    dt = 1e-2  # timestep
nsteps = 4000  # total number of steps
 nsubs = 20    # number of steps between each plot
nothing # hide

# ## Problem setup
# We initialize a `Problem` by providing a set of keyword arguments. The
# `stepper` keyword defines the time-stepper to be used.
prob = TwoDNavierStokes.Problem(dev; nx=n, Lx=L, ny=n, Ly=L, dt=dt, stepper="FilteredRK4")
nothing # hide

# Next we define some shortcuts for convenience.
sol, clock, vars, grid = prob.sol, prob.clock, prob.vars, prob.grid
x, y = grid.x, grid.y
nothing # hide

# ## Setting initial conditions

# Our initial condition tries to reproduce the initial condition used by McWilliams (_JFM_, 1984).
seed!(1234)
k₀, E₀ = 6, 0.5
ζ₀ = peakedisotropicspectrum(grid, k₀, E₀, mask=prob.timestepper.filter)
TwoDNavierStokes.set_ζ!(prob, ζ₀)
nothing # hide

# Let's plot the initial vorticity field. Note that when plotting, we decorate the variable
# to be plotted with `Array()` to make sure it is brought back on the CPU when `vars` live on
# the GPU.
heatmap(x, y, Array(vars.ζ'),
    aspectratio = 1,
              c = :balance,
           clim = (-40, 40),
          xlims = (-L/2, L/2),
          ylims = (-L/2, L/2),
         xticks = -3:3,
         yticks = -3:3,
         xlabel = "x",
         ylabel = "y",
          title = "initial vorticity",
     framestyle = :box)

# ## Diagnostics

# Create Diagnostics -- `energy` and `enstrophy` functions are imported at the top.
E = Diagnostic(TwoDNavierStokes.energy, prob; nsteps=nsteps)
Z = Diagnostic(TwoDNavierStokes.enstrophy, prob; nsteps=nsteps)
diags = [E, Z] # A list of Diagnostics types passed to "stepforward!" will  be updated every timestep.
nothing # hide

# ## Output

# We choose folder for outputing `.jld2` files and snapshots (`.png` files).
filepath = "."
plotpath = "./plots_decayingTwoDNavierStokes"
plotname = "snapshots"
filename = joinpath(filepath, "decayingTwoDNavierStokes.jld2")
nothing # hide

# Do some basic file management
if isfile(filename); rm(filename); end
if !isdir(plotpath); mkdir(plotpath); end
nothing # hide

# And then create Output
get_sol(prob) = prob.sol # extracts the Fourier-transformed solution
get_u(prob) = irfft(im * grid.l .* grid.invKrsq .* sol, grid.nx)
out = Output(prob, filename, (:sol, get_sol), (:u, get_u))
saveproblem(out)
nothing # hide

# ## Visualizing the simulation

# We initialize a plot with the vorticity field and the time-series of
# energy and enstrophy diagnostics.

p1 = heatmap(x, y, Array(vars.ζ'),
         aspectratio = 1,
                   c = :balance,
                clim = (-40, 40),
               xlims = (-L/2, L/2),
               ylims = (-L/2, L/2),
              xticks = -3:3,
              yticks = -3:3,
              xlabel = "x",
              ylabel = "y",
               title = "vorticity, t=" * @sprintf("%.2f", clock.t),
          framestyle = :box)

p2 = plot(2, # this means "a plot with two series"
               label = ["energy E(t)/E(0)" "enstrophy Z(t)/Z(0)"],
              legend = :right,
           linewidth = 2,
               alpha = 0.7,
              xlabel = "t",
               xlims = (0, 41),
               ylims = (0, 1.1))

l = @layout Plots.grid(1, 2)
p = plot(p1, p2, layout = l, size = (800, 360))

# ## Time-stepping the `Problem` forward

# We time-step the `Problem` forward in time.

startwalltime = time()

anim = @animate for j = 0:Int(nsteps/nsubs)
  if j % (1000 / nsubs) == 0
    cfl = clock.dt * maximum([maximum(vars.u) / grid.dx, maximum(vars.v) / grid.dy])

    log = @sprintf("step: %04d, t: %d, cfl: %.2f, ΔE: %.4f, ΔZ: %.4f, walltime: %.2f min",
        clock.step, clock.t, cfl, E.data[E.i]/E.data[1], Z.data[Z.i]/Z.data[1], (time()-startwalltime)/60)

    println(log)
  end

  p[1][1][:z] = Array(vars.ζ)
  p[1][:title] = "vorticity, t=" * @sprintf("%.2f", clock.t)
  push!(p[2][1], E.t[E.i], E.data[E.i]/E.data[1])
  push!(p[2][2], Z.t[Z.i], Z.data[Z.i]/Z.data[1])

  stepforward!(prob, diags, nsubs)
  TwoDNavierStokes.updatevars!(prob)
end

mp4(anim, "twodturb.mp4", fps=18)

# Last we can save the output by calling
# ```julia
# saveoutput(out)`
# ```

# ## Radial energy spectrum

# After the simulation is done we plot the instantaneous radial energy spectrum to illustrate
# how `FourierFlows.radialspectrum` can be used,

E  = @. 0.5 * (vars.u^2 + vars.v^2) # energy density
Eh = rfft(E)                  # Fourier transform of energy density
kr, Ehr = FourierFlows.radialspectrum(Eh, grid, refinement=1) # compute radial specturm of `Eh`
nothing # hide

# and we plot it.
plot(kr, abs.(Ehr),
    linewidth = 2,
        alpha = 0.7,
       xlabel = "kᵣ", ylabel = "∫ |Ê| kᵣ dk_θ",
        xlims = (5e-1, grid.nx),
       xscale = :log10, yscale = :log10,
        title = "Radial energy spectrum",
       legend = false)


text:
using Documenter, ManualScraper

makedocs()


text:
function get_ae(obj::Normalizer)
    return jcall(obj, "getAe", jdouble, ())
end

function get_max_degree(obj::Normalizer)
    return jcall(obj, "getMaxDegree", jint, ())
end

function get_max_order(obj::Normalizer)
    return jcall(obj, "getMaxOrder", jint, ())
end

function get_mu(obj::Normalizer)
    return jcall(obj, "getMu", jdouble, ())
end

function get_offset(obj::Normalizer, arg0::AbsoluteDate)
    return jcall(obj, "getOffset", jdouble, (AbsoluteDate,), arg0)
end

function get_reference_date(obj::Normalizer)
    return jcall(obj, "getReferenceDate", AbsoluteDate, ())
end

function get_tide_system(obj::Normalizer)
    return jcall(obj, "getTideSystem", TideSystem, ())
end

function on_date(obj::Normalizer, arg0::AbsoluteDate)
    return jcall(obj, "onDate", NormalizedSphericalHarmonicsProvider_NormalizedSphericalHarmonics, (AbsoluteDate,), arg0)
end


text:
using Test
import BenchmarkTools

@testset "solve: zero allocations" begin

    fs = (sin, cos, x -> -sin(x))
    x0 = (3,4)
    Ms = (Order1(), Order2(), Order5(), Order8(), Order16(),
          Roots.Order1B(), Roots.Order2B(),
          Roots.Bisection(), Roots.A42(), Roots.AlefeldPotraShi(), Roots.Brent(),
          Roots.Newton(), Roots.Halley(), Roots.Schroder()) # not Order0(), FalsePosition()

    for M ∈ Ms
        @test BenchmarkTools.@ballocated(solve(ZeroProblem($fs, $x0), $M)) == 0
    end

end

@testset "simple: zero allocations" begin
    @test BenchmarkTools.@ballocated( Roots.bisection(sin, 3,4) ) == 0
    @test BenchmarkTools.@ballocated( Roots.secant_method(sin, 3) ) == 0
    @test BenchmarkTools.@ballocated( Roots.muller(sin, 2.9,3.0,3.1) ) == 0
    @test BenchmarkTools.@ballocated( Roots.newton((sin,cos), 3) ) == 0
    @test BenchmarkTools.@ballocated( Roots.dfree(sin, 3) ) == 0
end


text:
"""
    copy_eigenvalues!(λs, A) → λs

Puts the eigenvalues of a quasi-upper triangular matrix A in the λs vector.
"""
function copy_eigenvalues!(λs, A::AbstractMatrix{T}, range = OneTo(size(A, 2)), tol = eps(real(T))) where {T}
    i = first(range)

    @inbounds while i < last(range)
        if is_offdiagonal_small(A, i, tol)
            λs[i] = A[i, i]
            i += 1
        else
            # Conjugate pair
            d = A[i,i] * A[i+1,i+1] - A[i,i+1] * A[i+1,i]
            x = (A[i,i] + A[i+1,i+1]) / 2
            y = sqrt(complex(x*x - d))
            λs[i + 0] = x + y
            λs[i + 1] = x - y
            i += 2
        end
    end

    @inbounds if i == last(range)
        λs[i] = A[i, i]
    end

    return λs
end

"""
    eigenvalue(R, i) → λ

Get the `i`th eigenvalue of R. NOTE: assumes `i` points to the start of
a block.
"""
function eigenvalue(R, i)
    n = minimum(size(R))

    @inbounds begin
        if i == n || iszero(R[i+1,i])
            return complex(R[i,i])
        else
            d = R[i,i] * R[i+1,i+1] - R[i,i+1] * R[i+1,i]
            x = (R[i,i] + R[i+1,i+1]) / 2
            y = sqrt(complex(x*x - d))
            return x + y
        end
    end
end

"""
    eigenvalues(A::AbstractMatrix{T}) → Vector{complex(T)}

Computes the eigenvalues of the matrix A. Assumes that A is quasi-upper triangular.
The eigenvalues are returned in complex arithmetic, even if their imaginary
part is 0.
"""
eigenvalues(A::AbstractMatrix{T}, tol = eps(real(T))) where {T} =
    copy_eigenvalues!(Vector{complex(T)}(undef, size(A, 2)), A, OneTo(size(A, 2)), tol)

"""
    partialeigen(P::PartialSchur) → (Vector{<:Union{Real,Complex}}, Matrix{<:Union{Real,Complex}})

Transforms a partial Schur decomposition into an eigendecomposition.

!!! note

    For real-symmetric and Hermitian matrices the Schur vectors coincide with
    the eigenvectors, and hence it is not necessary to call this function in
    that case.

The method still relies on LAPACK to compute the eigenvectors of the (quasi)
upper triangular matrix `R` from the partial Schur decomposition.

!!! note

    This method is currently type unstable for real matrices, since we have not
    yet decided how to deal with complex conjugate pairs of eigenvalues. E.g.
    if almost all eigenvalues are real, but there are just a few conjugate
    pairs, should all eigenvectors be complex-valued?
"""
function partialeigen(P::PartialSchur)
    vals, vecs = eigen(P.R)
    return vals, P.Q*vecs
end


text:
"""
    BSpline{T}
bspline is a structures which includes the knots and the
vertices of a BSpline.
"""
struct BSpline{T}
    knots::Vector{T}
    vertices::Vector{T}
end

"""
    knots(bs::BSpline)
Returns the field knots for a certain `BSpline`.
"""
knots(bs::BSpline) = bs.knots

"""
    vertices(bs::BSpline)
Returns the field vertices for a certain `BSpline`.
"""
vertices(bs::BSpline) = bs.vertices

"""
    BSpline2dim{T}
bspline2dim is a structures which includes the knots in x and y direction
and the vertices (in vector form) of a 2 dimensional BSpline.
"""
struct BSpline2dim{T}
    knotsx::Vector{T}
    knotsy::Vector{T}
    vertices::Vector{T}
end
## Document
knotsx(bsdim2::BSpline2dim) = bsdim2.knotsx
knotsy(bsdim2::BSpline2dim) = bsdim2.knotsy
vertices(bsdim2::BSpline2dim) = bsdim2.vertices

"""
    buildKnots(xin, xend, nvertices)
Returns the knots of a `BSpline` with beginning xin, end xend
with nvertices components. The knots returned by this function
are equispaced between xin and xend.

#Examples
```julia julia-repl
julia> buildKnots(1.,5.,10)
10-element Vector{Float64}:
 0.4285714285714286
 1.0
 1.5714285714285714
 2.142857142857143
 2.7142857142857144
 3.2857142857142856
 3.857142857142857
 4.428571428571429
 5.0
 5.571428571428571
```
"""
function buildKnots(xin, xend, nvertices)
    knotsIn = LinRange(xin, xend, nvertices-2)
    knots = vcat(knotsIn[1]-(knotsIn[2]-knotsIn[1]), knotsIn,  knotsIn[end]+(knotsIn[2]-knotsIn[1]))
    return knots
end


text:
using DynamicalSystemsBase
using DynamicalSystemsBase: CDS, DS, columns
using Test, StaticArrays, OrdinaryDiffEq

println("\nTesting continuous system evolution...")

@testset "Lorenz63 oop" begin

    lo11 = Systems.lorenz() #with Jac
    # lo33 = Systems.lorenz(big.([0.0, 10.0, 0.0]))
    # trajectory pure:
    ts1 = trajectory(lo11, 1.0, Δt = 0.01)
    # ts3 = trajectory(lo33, 1.0)
    # @test eltype(ts3[1]) == BigFloat
    @test size(ts1) == (101, 3)
    # @test ts1[1, :] ≈ ts3[end,:]
    @test ts1[1, :] == SVector{3}(lo11.u0)
    @test !any(x -> abs(x) > 1e3, ts1[end])
    # trajectory with diff_eq_kwargs and Δt:
    ts1 = trajectory(lo11, 1.0; Δt=0.1,
    diff_eq_kwargs=Dict(:abstol=>1e-8, :reltol=>1e-8))
    # ts3 = trajectory(lo33, 1.0; Δt=0.1,
    # diff_eq_kwargs=Dict(:abstol=>1e-8, :reltol=>1e-8))
    @test size(ts1) == (11, 3)
    @test ts1[1, :] == SVector{3}(lo11.u0)
    # @test ts1[end, :] ≈ ts3[end,:]

    data2 = trajectory(lo11, 100; save_idxs = 1:2)
    @test size(data2)[2] == 2
end

@testset "Lorenz96 iin" begin
    u = ones(5)
    lo11 = Systems.lorenz96(5, u)
    # lo33 = Systems.lorenz96(5, big.(ones(5)))
    # trajectory pure:
    ts1 = trajectory(lo11, 2.0, Δt = 0.01)
    # ts3 = trajectory(lo33, 2.0)
    # @test eltype(ts3[1]) == BigFloat
    @test size(ts1) == (201, 5)
    # @test ts1[end, :] ≈ ts3[end,:]
    @test ts1[1, :] == SVector{5}(u)
    # trajectory with diff_eq_kwargs and Δt:
    ts1 = trajectory(lo11, 2.0; Δt=0.1,
    diff_eq_kwargs=Dict(:abstol=>1e-9, :reltol=>1e-9))
    # ts3 = trajectory(lo33, 2.0; Δt=0.1,
    # diff_eq_kwargs=Dict(:abstol=>1e-9, :reltol=>1e-9))
    @test size(ts1) == (21, 5)
    # @test ts1[end, :] ≈ ts3[end,:]
    @test ts1[1, :] == SVector{5}(u)
    @test !any(x -> abs(x) > 1e3, ts1[end])
end


text:
using Pkg, Sockets, Serialization, Logging

import Base.Threads.@spawn

# assumes the julia executable starts in the default environment
const default_env = dirname(Pkg.project().path)

# both this script and the symbol server is going to use this shared environment
Pkg.activate("linter-julia", shared=true)

### variables ###

# packages needed and minimum versions
needed_packages = Dict(
    "StaticLint" => v"7.0.0",
    "SymbolServer" => v"6.0.1",
    "CSTParser" => v"3.1.0",
    "JSON" => v"0.21.1"
)

@inline function withtrace(e, msg=nothing)
    buffer = IOBuffer();
    if !isnothing(msg)
        println(buffer, msg)
    end
    println(buffer, e)
    st = stacktrace(catch_backtrace());
    for (idx, l) in enumerate(st)
        println(buffer, " [$idx] $l")
    end
    return(String(take!(buffer)))
end

const store_location = dirname(Pkg.project().path)
const port = ARGS[1]
const atom_pid = parse(Int32, split(port,"_")[end])

if ispath(port)
    # if another is already set up, just finish quietly
    exit()
end

# it is needed for the log file at the very first install
mkpath(store_location)

# remove previous logs
logs = readdir(store_location, join = true)
for f in logs
    if basename(splitext(f)[1]) == "symbolserver.log"
        rm(f)
    end
end

logger = open(joinpath(store_location, "symbolserver.log." * string(atom_pid)), "w")

global_logger( SimpleLogger(logger) )

### environment set-up ###

# check if the necessary packages are correctly installed

installed_packages = Dict(pkg.name=>pkg.version for pkg in values(Pkg.dependencies()) if pkg.is_direct_dep == true )

for (pkg, minver) in needed_packages
    try
        if !(pkg in keys(installed_packages))
            @info "installing $pkg"
            Pkg.add(pkg)
        elseif installed_packages[pkg] < minver
            @info "updating $pkg"
            Pkg.update(pkg)
        end
    catch err
        @error withtrace(err, "something went wrong while installing $pkg")
        rethrow()
    end
end

# load the packages here

using StaticLint, SymbolServer, CSTParser, JSON

# list of hints that are going to be classified as "info"

const infohints = [ StaticLint.ConstIfCondition, StaticLint.PointlessOR, StaticLint.PointlessAND,
    StaticLint.UnusedBinding, StaticLint.UnusedTypeParameter, StaticLint.TypePiracy,
    StaticLint.UnusedFunctionArgument, StaticLint.NotEqDef, StaticLint.InappropriateUseOfLiteral  ]

### types ###

mutable struct LintMsg
    filename::String
    message::String
    code::String
    severity::String
    startpos::Int64
    endpos::Int64
end

mutable struct SSData
    server::Union{StaticLint.FileServer, Nothing}
    project_mtime::Union{Float64,Nothing}
    manifest_mtime::Union{Float64,Nothing}
    # lock mechanism
    available::Channel{Bool}
    # generation message sent
    msgsent::Bool
end

### variables ###

# this maps environments to symbol servers
symbolservers = Dict{AbstractString, SSData}()
server_queue = Channel(Inf)

### functions ###

"""
    guess_environment(fname::AbstractString)::AbstractString

Guesses the most likely environment path for a .jl file.
If no Project.toml is found in parents, the default environment is returned.
"""
function guess_environment(fname::AbstractString)::AbstractString

    if isempty(fname)
        return default_env
    end

    try
        fname = realpath(fname)
        actdir = dirname(fname)

        while length(splitpath(actdir)) != 1
            act_candidate = joinpath(actdir, "Project.toml")
            if isfile(act_candidate)
                return dirname(act_candidate)
            end
            actdir = dirname(actdir)
        end
    catch e
        @warn "error while getting environment of file: $fname"
    end

    # return the default environment
    return default_env

end

function exit_if_atom_dies()

    PROCESS_QUERY_INFORMATION = 0x0400
    ERROR_INVALID_PARAMETER = 0x57
    STILL_ACTIVE = 259

    function CloseHandle(handle)
        Base.windowserror(:CloseHandle, 0 == ccall(:CloseHandle, stdcall, Cint, (Ptr{Cvoid},), handle))
        nothing
    end

    function OpenProcess(id::Integer, rights = PROCESS_QUERY_INFORMATION)
        proc = ccall((:OpenProcess, "kernel32"), stdcall, Ptr{Cvoid}, (UInt32, Cint, UInt32), rights, false, id)
        Base.windowserror(:OpenProcess, proc == C_NULL)
        proc
    end

    while true

        if !Sys.iswindows()

            if ccall(:kill, Int32, (Int32, Int32), atom_pid, 0) != 0
                exit()
            end

        else
            # an Absolutely Ridiculous Difference ...

            hProcess = try OpenProcess(atom_pid)
            catch err
                if err isa SystemError && err.extrainfo.errnum == ERROR_INVALID_PARAMETER
                    exit()
                end
            end

            exitCode = Ref{UInt32}()

            if ccall(:GetExitCodeProcess, stdcall, Cint, (Ptr{Cvoid}, Ref{UInt32}), hProcess, exitCode) != 0
                CloseHandle(hProcess)
                if exitCode[] != STILL_ACTIVE
                    exit()
                end
            else
                CloseHandle(hProcess)
            end

        end
        sleep(5)
    end
end

"""
    guess_package_root(environment::AbstractString)::Union{AbstractString,nothing}

Guesses the most likely root for a package.
Can return nothing if a likely candidate is not found.
"""
function guess_package_root(env::AbstractString)::Union{AbstractString,Nothing}

    try
        # guess a project name first
        have_name = false
        if isfile( joinpath(env, "Project.toml") )
            projectdef = Pkg.TOML.parsefile( joinpath(env, "Project.toml") )
            if "name" in keys(projectdef)
                projectname = projectdef["name"]
                have_name = true
            end
        end
        if !have_name
            projectname = splitext(basename(env))[1]
        end

        candidates = [
            joinpath(env, "src", projectname * ".jl"),
            joinpath(env, projectname * ".jl")
        ]

        for c in candidates
            if isfile(c)
                return c
            end
        end

    catch e
        @warn "error while getting package root of env: $env"
    end

    return nothing

end

function generate_messages( fname::AbstractString, code::AbstractString, env::AbstractString, ss )::Vector{LintMsg}

    if isnothing(ss.server)
        if ss.msgsent == false
            ss.msgsent = true
            @debug "waiting message generated for $fname"
            return [ LintMsg(fname, "$env", "I000", "info", 1, 1) ]
        else
            return []
        end
    end

    rootfile = guess_package_root(env)
    if isnothing(rootfile)
        rootfile = fname
    end

    # lock server
    take!(ss.available)

    # this is a hack from StaticLint.lint_file() and StaticLint.sslint_string()

    empty!(ss.server.files)
    root = !isempty(rootfile) ? StaticLint.loadfile(ss.server, rootfile) : nothing

    buffer = StaticLint.File(fname, code, CSTParser.parse(code, true), root, ss.server)
    if isnothing(root) || fname == rootfile
        StaticLint.setroot(buffer, buffer)
        root = buffer
    end
    StaticLint.setfile(ss.server, fname, buffer)

    StaticLint.semantic_pass(buffer)
    StaticLint.semantic_pass(root)

    for (p,f) in ss.server.files
        StaticLint.check_all(f.cst, StaticLint.LintOptions(), ss.server)
    end

    msgs = LintMsg[]

    for (path, file) in ss.server.files
        # skip hints unrelated to this file
        if path != fname
            continue
        end
        act_hints = StaticLint.collect_hints(file.cst, ss.server)
        for (offset, x) in act_hints
            if StaticLint.haserror(x)
                msg = StaticLint.LintCodeDescriptions[x.meta.error]
                if x.meta.error in infohints
                    code = "I"*string(Int32(x.meta.error), pad=3)
                    severity = "info"
                else
                    code = "W"*string(Int32(x.meta.error), pad=3)
                    severity = "warning"
                end
            elseif CSTParser.headof(x) === :errortoken
                msg = "Parsing error"
                code = "E000"
                severity = "error"
            else
                msg = "Missing reference: '$(CSTParser.valof(x))'"
                code = "W000"
                severity = "warning"
            end
            push!( msgs, LintMsg(path, msg, code, severity, offset, offset + x.fullspan) )
        end
    end

    # enable the server for other threads
    put!(ss.available, true)

    @debug "hints generated for $fname: $msgs"

    return msgs

end

function generate_symbolserver(ch, env::AbstractString)

    server = StaticLint.setup_server( env )
    put!(ch, (env, server))

end

function find_line_column(io::IOBuffer, position::Int64)

    line = 1
    column = 1

    try
        before = read(io, position)
        line = count(x -> x == UInt8('\n'), before) + 1
        prevreturn = findprev( x -> x == UInt8('\n'), before, length(before) )
        column = isnothing(prevreturn) ? position : position - prevreturn
    catch err
        @error withtrace(err, "error while determining line and column for position $position")
    end

    return line, column
end

function find_line_column_file(fname::AbstractString, position::Int64)
    open(fname, "r") do io
        return find_line_column(io, position)
    end
end

function find_line_column_string(str::AbstractString, position::Int64)
    return find_line_column(IOBuffer(str), position)
end

function get_last_charpos(str::AbstractString, position::Int64)::Int64
    prev = findprev( x -> !isspace(x), str, position )
    if isnothing(prev)
        prev = 1
    end
    return prev
end

function convertmsgtojson(msgs, code)
    output = Any[]
    for msg in msgs

        # determine line and column from the file
        @assert msg.startpos <= msg.endpos
        if msg.endpos > lastindex(code)
            @warn "endpos $(msg.endpos) is higher than code length $(lastindex(code))"
            msg.endpos = lastindex(code)
        end
        endpos = get_last_charpos(code, msg.endpos)
        startline, startcolumn = find_line_column_string(code, msg.startpos)
        endline, endcolumn = find_line_column_string(code, endpos)

        if startline > endline
            startline = endline
            startcolumn = endcolumn
        end
        if startline == endline && startcolumn > endcolumn
            startcolumn = endcolumn
        end

        # this looks like a bug in linter: as we type at the end of the buffer, the editor seems to propagate
        # the end position, despite linter returning the correct line/column. Fix this here by marking one less
        # character if we are the end of the buffer. A bit ugly when typing, but less ugly than dragging the error
        if endpos == length(code)
            endcolumn -= 1
        end

        # Atom index starts from zero thus minus one
        errorrange = Array[[startline-1, startcolumn], [endline-1, endcolumn]]

        # this seems to be by design, https://github.com/steelbrain/linter/issues/1235
        # weird is that it seems to call the linter, but does not display the returned data
        if isempty(msg.filename)
            loc = Dict("position" => errorrange)
        else
            loc = Dict("file" => msg.filename, "position" => errorrange)
        end

        push!(output, Dict( "severity" => msg.severity,
                            "location" => loc,
                            "excerpt" => msg.message,
                            "description" => msg.code) )

    end
    return output
end

function filtermsgs(msgs, data)
    if haskey(data,"ignore_warnings")
        if data["ignore_warnings"]
            msgs = filter(i -> i.severity != "warning", msgs)
        end
    end
    if haskey(data,"ignore_info")
        if data["ignore_info"]
            msgs = filter(i -> i.severity != "info", msgs)
        end
    end
    if haskey(data,"ignore_codes")
        msgs = filter(i -> !(i.code in data["ignore_codes"]), msgs)
    end

    return msgs
end

function handle_connection(conn)

    global symbolservers, server_queue

    try

        # insert waiting servers
        while isready(server_queue)
            # there might be a race condition lurking around right here,
            # but this is light without a yield, and we do not use threading here
            senv, server = take!(server_queue)
            symbolservers[senv].server = server
            @info "added server for $senv"
        end

        data = JSON.parse(conn)

        haskey(data, "file") ? fname = data["file"] : fname = ""

        @debug "request for file : $fname"

        env = guess_environment(fname)

        # check mtime of Project.toml and Manifest.toml
        project_mtime = isfile(joinpath(env, "Project.toml")) ? stat( joinpath(env, "Project.toml") ).mtime : nothing
        manifest_mtime = isfile(joinpath(env, "Manifest.toml")) ? stat( joinpath(env, "Manifest.toml") ).mtime : nothing

        # build symbol server if needed
        if !(env in keys(symbolservers)) || symbolservers[env].project_mtime != project_mtime || symbolservers[env].manifest_mtime != manifest_mtime
            @debug "launching server generation for $env"
            symbolservers[env] = SSData(nothing, project_mtime, manifest_mtime, Channel{Bool}(1), false)
            put!(symbolservers[env].available, true)
            @spawn generate_symbolserver(server_queue, env)
        end

        msgs = generate_messages( fname, data["code_str"], env, symbolservers[env] )

        msgs = filtermsgs(msgs, data)
        out = convertmsgtojson(msgs, data["code_str"])
        JSON.print(conn, out)

        @debug "symbolservers: " [ (k, !isnothing(v.server)) for (k,v) in symbolservers ]

        flush(logger)

    catch err
        @error withtrace(err, "connection ended with error")
    finally
        close(conn)
    end

end

### main script starts here ###

@spawn exit_if_atom_dies()

try
    server = listen(port)

    # this is looping until Atom dies
    while true
        conn = accept(server)
        # the pipe is fifo anyway and cannot deal with multiple connections
        @async handle_connection(conn)
    end

finally
    close(server)
end

max_stars_repo_name: takbal/linter-julia

max_issues_repo_name: takbal/linter-julia

max_forks_repo_name: takbal/linter-julia


text:
struct mv_gaussian <: distibution_sample
    μ::AbstractArray{Float64,1}
    Σ::AbstractArray{Float64,2}
    invΣ::AbstractArray{Float64,2}
    logdetΣ::Float64
    # mvn::MvNormal
end

function dinvquad!(r,a,x)
    dcolwise_dot!(r,x, a \ x)
end

function log_likelihood!(r::AbstractArray,x::AbstractArray, distibution_sample::mv_gaussian , group::Int64 = -1)
    z = x .- distibution_sample.μ
    dcolwise_dot!(r,z, distibution_sample.invΣ * z)
    r .= -((length(distibution_sample.Σ) * Float64(log(2π)) + logdet(distibution_sample.Σ))/2) .-r
end


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

# Breadth-first search
export BFS
"""
  `BFS(fst,initial)`

Returns an iterator that performs a breadth-first search (BFS) of `fst` starting from the state `initial`.

At each iteartion the tuple `(i,p)` is returned, where `i` is the current state and `p` a [`Path`](@ref).
"""
struct BFS{W,A,D,F <: WFST{W,A,D}}
  fst::F
  initial::D
  function BFS(fst::F,initial::D) where {W,A,D,F <: WFST{W,A,D}}
    if initial==0
      throw(ErrorException("Invalid initial state"))
    end
    new{W,A,D,F}(fst,initial)
  end
end

function init_bfs(fst::WFST{W,A,D}, initial::D) where {W,D,A<:Arc{W,D}}
  q = Queue{Tuple{D,Path{W,D}}}()
  isym = get_isym(fst)
  osym = get_osym(fst)
  w = get_initial_weight(fst,initial)
  p = Path(Vector{D}(), Vector{D}(), w, isym, osym)
  enqueue!(q, (initial, p))
  return q
end

function Base.iterate(it::BFS,
                      q = init_bfs(it.fst,it.initial))
  if isempty(q)
    return nothing
  else
    i, p = dequeue!(q)
    s = it.fst[i]
    for arc in s
      p_new = update_path(p, get_ilabel(arc), get_olabel(arc), get_weight(arc))
      enqueue!(q, (get_nextstate(arc), p_new))
    end
    return (i,p), q
  end
end


text:
export iauZp
"""
**  Zero a p-vector.
**
**  This function is part of the International Astronomical Union's
**  SOFA (Standards Of Fundamental Astronomy) software collection.
**
**  Status:  vector/matrix support function.
**
**  Returned:
**     p        double[3]      p-vector
**
**  This revision:  2013 June 18
**
**  SOFA release 2018-01-30
**
**  Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauZp(double p[3])

function iauZp()
   p = zeros(Float64, 3)

   ccall((:iauZp, libsofa_c), Cvoid,
           (Ptr{Cdouble},),
           p)

   # Transpose since C call return row-major operation
   return SVector{3}(p)
end


text:
"""
Validation PSSE/TGOV1:
This case study defines a three bus system with an infinite bus, GENROU+AC1A+TGOV1 and a load.
The fault drop the line connecting the infinite bus and GENROU.
"""

##################################################
############### SOLVE PROBLEM ####################
##################################################

raw_file = joinpath(dirname(@__FILE__), "benchmarks/psse/TGOV1/ThreeBusMulti.raw")
dyr_file = joinpath(dirname(@__FILE__), "benchmarks/psse/TGOV1/ThreeBus_TGOV1.dyr")
csv_file = joinpath(dirname(@__FILE__), "benchmarks/psse/TGOV1/TEST_TGOV1.csv")

#Construct system
sys = System(raw_file, dyr_file);

path = (joinpath(pwd(), "test-psse-tgov1"))
!isdir(path) && mkdir(path)
try
    sim = Simulation!(
        ImplicitModel,
        sys, #system
        path,
        (0.0, 20.0), #time span
        BranchTrip(1.0, "BUS 1-BUS 2-i_1"), #Type of Fault
    )

    #Solve problem in equilibrium
    execute!(sim, IDA(), dtmax = 0.005, saveat = 0.005)

    #Obtain small signal results for initial conditions
    #NOT WORKING DUE TO TYPES ON EXECUTE
    #small_sig = small_signal_analysis(sim)

    series = get_state_series(sim, ("generator-102-1", :δ))
    t = series[1]
    δ = series[2]

    #Obtain PSSE results
    t_psse, δ_psse = get_csv_delta(csv_file)

    diff = [0.0]
    res = get_init_values_for_comparison(sim)
    for (k, v) in test_psse_tgov1_init
        diff[1] += LinearAlgebra.norm(res[k] - v)
    end
    #Test Initial Condition
    @test (diff[1] < 1e-3)
    #Test Solution DiffEq
    @test sim.solution.retcode == :Success
    #Test Small Signal
    #@test small_sig.stable
    #Test Transient Simulation Results
    # PSSE results are in Degrees
    @test LinearAlgebra.norm(δ - (δ_psse .* pi / 180), Inf) <= 1e-2
    @test LinearAlgebra.norm(t - round.(t_psse, digits = 3)) == 0.0
finally
    @info("removing test files")
    rm(path, force = true, recursive = true)
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## UMFPACK

## Type of solve
const UMFPACK_A     =  0     # Ax=b
const UMFPACK_At    =  1     # adjoint(A)x=b
const UMFPACK_Aat   =  2     # transpose(A)x=b
const UMFPACK_Pt_L  =  3     # adjoint(P)Lx=b
const UMFPACK_L     =  4     # Lx=b
const UMFPACK_Lt_P  =  5     # adjoint(L)Px=b
const UMFPACK_Lat_P =  6     # transpose(L)Px=b
const UMFPACK_Lt    =  7     # adjoint(L)x=b
const UMFPACK_Lat   =  8     # transpose(L)x=b
const UMFPACK_U_Qt  =  9     # U*adjoint(Q)x=b
const UMFPACK_U     =  10    # Ux=b
const UMFPACK_Q_Ut  =  11    # Q*adjoint(U)x=b
const UMFPACK_Q_Uat =  12    # Q*transpose(U)x=b
const UMFPACK_Ut    =  13    # adjoint(U)x=b
const UMFPACK_Uat   =  14    # transpose(U)x=b

## Sizes of Control and Info arrays for returning information from solver
const UMFPACK_INFO = 90
const UMFPACK_CONTROL = 20
const UMFPACK_PRL = 1

## Status codes
const UMFPACK_OK = 0
const UMFPACK_WARNING_singular_matrix       = 1
const UMFPACK_WARNING_determinant_underflow = 2
const UMFPACK_WARNING_determinant_overflow  = 3
const UMFPACK_ERROR_out_of_memory           = -1
const UMFPACK_ERROR_invalid_Numeric_object  = -3
const UMFPACK_ERROR_invalid_Symbolic_object = -4
const UMFPACK_ERROR_argument_missing        = -5
const UMFPACK_ERROR_n_nonpositive           = -6
const UMFPACK_ERROR_invalid_matrix          = -8
const UMFPACK_ERROR_different_pattern       = -11
const UMFPACK_ERROR_invalid_system          = -13
const UMFPACK_ERROR_invalid_permutation     = -15
const UMFPACK_ERROR_internal_error          = -911
const UMFPACK_ERROR_file_IO                 = -17
const UMFPACK_ERROR_ordering_failed         = -18

max_stars_repo_name: greimel/julia

max_issues_repo_name: greimel/julia

max_forks_repo_name: greimel/julia


text:
push!(LOAD_PATH,"src");
push!(LOAD_PATH,"modules");
push!(LOAD_PATH,"submodules");


text:
# TODO: seed, copy

"""
A box in R^n.
I.e., each coordinate is bounded.

Example usage:
box_obj.action_space = Box(-10, 10, (1,))

Two kinds of valid input:
    Box(low=-1.0, high=1.0, shape=(3,4)) # low and high are scalars, and shape is provided
    Box(low=[-1.0,-2.0], high=[2.0,4.0]) # low and high are arrays of the same shape
"""
mutable struct Box <: AbstractSpace
    low::Array
    high::Array
    shape::Tuple
    dtype::DataType
    #seed::Int
end

function Box(low::Number, high::Number, shape::Union{Tuple, Array{Int64, 1}}, dtype::Union{DataType, Nothing}=nothing)
    if isnothing(dtype)
        dtype = high == 255 ? UInt8 : Float32
        @warn "dtype was autodetected as $(dtype). Please provide explicit data type."
    end
    if dtype <: Integer
        low = floor(dtype, low)
        high = floor(dtype, high)
    end
    Low = dtype(low) .+ zeros(dtype, shape)
    High = dtype(high) .+ zeros(dtype, shape)
    return Box(Low, High, shape, dtype)
end

function Box(low::Array, high::Array, dtype::Union{DataType, Nothing}=nothing)
    @assert size(low) == size(high)
    shape = size(low)
    if isnothing(dtype)
        dtype = all(high .== 255) ? UInt8 : Float32
        @warn "dtype was autodetected as $(dtype). Please provide explicit data type."
    end
    if dtype <: Integer
        low = floor.(dtype, low)
        high = floor.(dtype, high)
    else
        low = dtype.(low)
        high = dtype.(high)
    end
    return Box(low, high, shape, dtype)
end
#=
function seed!(box_obj::Box, seed::Int)
    box_obj.seed = seed
end
=#
function sample(box_obj::Box)
    box_obj.dtype <: AbstractFloat ?
        rand(box_obj.dtype, box_obj.shape) .* (box_obj.high .- box_obj.low) .+ box_obj.low :
        rand.(UnitRange.(box_obj.low, box_obj.high))
end

function contains(box_obj::Box, x)
    isa(x, Number) && box_obj.shape == (1,) && (x = [x])
    size(x) == box_obj.shape && all(box_obj.low .<= x .<= box_obj.high)
end

Base.:(==)(box_obj::Box, other::Box) = isapprox(box_obj.low, other.low) && isapprox(box_obj.high, other.high)


text:
function jackknifeEstimate(x, func, θ::Float64)
    n = length(x)

    je = zeros(n)
    for i in 0:n-1
        je[i+1] = func(x[1:end-i])
    end

    return mean(je) .- je
end

function jackknifeEstimate(x, y, func, θ::Float64)
    n = minimum([length(x), length(y)])

    je = zeros(n)
    for i in 0:n-1
		je[i+1] = func(x[1:end-i], y[1:end-i])
    end

	return mean(je) .- je
end


text:
function rouwenhorst(n, mu, sigma, rho)

    # grid
    width = sqrt((n-1) * sigma^2 / ( 1 - rho^2))
    grid = vcat(range( mu-width, mu + width, length=n))

    # transition matrix
    p0 = (1 + rho) / 2
    trans = [p0 1-p0; 1-p0 p0]

    if n > 2
        for i = 1:n-2
            cstr_temp = zeros(length(trans[:,1]), 1)
            trans = p0 .* [trans cstr_temp; cstr_temp' 0] + (1 - p0 ) .* [cstr_temp trans; cstr_temp' 0]  + (1 - p0 ) .*  [ cstr_temp' 0; trans cstr_temp] + p0 .* [ cstr_temp' 0; cstr_temp trans]
        end
        trans = trans ./ sum(trans,dims=2)
    end

    # ergodic distribution
    dist = ones(1,n)./n
    for i = 1:100
        dist = dist*(trans^i)
    end
    dist = vcat(dist...)

    return grid, trans, dist
end


text:
"""
    estimateKernelSize{T <: Real}(X::Array{T, 2}; sampleSize = 100)

    Estimate kernel size by taking median distance between points.
"""
function estimateKernelSize(X::Array{T, 2}; sampleSize = 100) where T <: Real

	M = size(X)[1]

	# set kernel size to median distance between points
	if M > sampleSize
		Xmed = X[shuffle(collect(1:M))[1:sampleSize],:]
		S = sampleSize
	else
		Xmed = X
		S = M
	end

	dists = pairwise(SqEuclidean(), Xmed, Xmed, dims=1)
    sig = sqrt(0.5 * median(dists))

    return sig
end

"""
    rbfDotProduct(X::Array{T, 2}, X::Array{T, 2}, kernelSize::Float64)

"""
function rbfDotProduct(X::Array{T, 2}, Y::Array{T, 2}, kernelSize) where T <: Real

    H = zeros(size(X,1), size(Y,1))
    @assert size(X,2) == size(Y,2)

    for i in 1:size(X,1)
        for j in 1:size(Y,1)
            @inbounds begin
                H[i,j] = first( (X[i,:].-Y[j,:]) * (X[i,:].-Y[j,:])' )
            end
        end
    end

    return exp.(-H/2/kernelSize^2)
end

eye(n) = Matrix(I, n, n)


text:
module Advanced

# Only exports advanced stuff which one typically can't be arsed to use from the REPl
using Reexport: @reexport
# Traits
@reexport using ..NaiveNASlib: NamedTrait, AfterΔSizeTrait, named, logged, validated

# Size strategies
@reexport using ..NaiveNASlib:  DefaultJuMPΔSizeStrategy, ThrowΔSizeFailError, ΔSizeFailNoOp, LogΔSizeExec, ΔNout, ΔNoutExact, ΔNoutRelaxed, ΔNin,
                                ΔNinExact, ΔNinRelaxed, AlignNinToNout, TruncateInIndsToValid, WithUtilityFun, WithKwargs, TimeLimitΔSizeStrategy,
                                TimeOutAction, UpperInsertBound, AfterΔSizeCallback, logafterΔsize, validateafterΔsize

# Align size strategies, e.g what to do with sizes of vertices connected to a removed vertex
@reexport using ..NaiveNASlib:  IncreaseSmaller, DecreaseBigger, AlignSizeBoth, ChangeNinOfOutputs, FailAlignSizeNoOp, FailAlignSizeError,
                                FailAlignSizeWarn, NoSizeChange, CheckAligned, CheckNoSizeCycle, CheckCreateEdgeNoSizeCycle, PostAlign

# Connect strategies
@reexport using ..NaiveNASlib: ConnectAll, ConnectNone

# For elementwise operations
@reexport using ..NaiveNASlib: VertexConf, traitconf, outwrapconf

@reexport using ..NaiveNASlib: RemoveStrategy

@reexport using ..NaiveNASlib: findterminating, all_in_graph, output!, ancestors, descendants, defaultutility

end