text:
"""
    Reduced_lattice

Struct containing symmetry irreducible sites of a lattice graph.
* `name     :: String`                  : name of the original lattice
* `size     :: Int64`                   : bond truncation of the original lattice
* `model    :: String`                  : name of the initialized model
* `J        :: Vector{Vector{Float64}}` : coupling vector of the initialized model
* `sites    :: Vector{Site}`            : list of symmetry irreducible sites
* `overlap  :: Vector{Matrix{Int64}}`   : pairs of irreducible sites with their respective multiplicity in range of origin and another irreducible site
* `mult     :: Vector{Int64}`           : multiplicities of irreducible sites
* `exchange :: Vector{Int64}`           : images of the pair (origin, irreducible site) under site exchange
* `project  :: Matrix{Int64}`           : projections of pairs (site1, site2) of the original lattice to pair (origin, irreducible site)
"""
struct Reduced_lattice
    name     :: String
    size     :: Int64
    model    :: String
    J        :: Vector{Vector{Float64}}
    sites    :: Vector{Site}
    overlap  :: Vector{Matrix{Int64}}
    mult     :: Vector{Int64}
    exchange :: Vector{Int64}
    project  :: Matrix{Int64}
end

# check if matrix in list of matrices within numerical tolerance
function is_in(
    e    :: SMatrix{3, 3, Float64},
    list :: Vector{SMatrix{3, 3, Float64}}
    )    :: Bool

    in = false

    for item in list
        if maximum(abs.(item .- e)) <= 1e-8
            in = true
            break
        end
    end

    return in
end

# rotate vector onto a reference vector using Rodrigues formula
function get_rotation(
    vec :: SVector{3, Float64},
    ref :: SVector{3, Float64}
    )   :: SMatrix{3, 3, Float64}

    # buffer geometric information
    a = vec ./ norm(vec)
    b = ref ./ norm(ref)
    n = cross(a, b)
    s = norm(n)
    c = dot(a, b)

    # check if vectors are collinear
    if s < 1e-8
        # if vector are antiparallel do inversion
        if c < 0.0
            return SMatrix{3, 3, Float64}(-1.0, 0.0, 0.0, 0.0, -1.0, 0.0, 0.0, 0.0, -1.0)
        # if vectors are parallel return unity
        else
            return SMatrix{3, 3, Float64}(+1.0, 0.0, 0.0, 0.0, +1.0, 0.0, 0.0, 0.0, +1.0)
        end
    # if vectors are non-collinear use Rodrigues formula
    else
        n    = n ./ s
        temp = SMatrix{3, 3, Float64}(0.0, n[3], -n[2], -n[3], 0.0, n[1], n[2], -n[1], 0.0)
        mat  = SMatrix{3, 3, Float64}(+1.0, 0.0, 0.0, 0.0, +1.0, 0.0, 0.0, 0.0, +1.0)
        mat  = mat .+ s .* temp .+ (1.0 - c) .* temp * temp

        return mat
    end
end

# try to rotate vector onto a reference vector around a given axis (return matrix of zeros in case of failure)
function get_rotation_for_axis(
    vec  :: SVector{3, Float64},
    ref  :: SVector{3, Float64},
    axis :: SVector{3, Float64}
    )    :: SMatrix{3, 3, Float64}

    # normalize vectors and axis
    a  = vec ./ norm(vec)
    b  = ref ./ norm(ref)
    ax = axis ./ norm(axis)

    # determine othogonal projections onto axis
    ovec = vec .- dot(vec, ax) .* ax
    ovec = ovec ./ norm(ovec)
    oref = ref .- dot(ref, ax) .* ax
    oref = oref ./ norm(oref)

    # buffer geometric information
    s = norm(cross(ovec, oref))
    c = dot(ovec, oref)

    # allocate rotation matrix
    mat = SMatrix{3, 3, Float64}(ax[1]^2 + (1.0 - ax[1]^2) * c,
                                 ax[1] * ax[2] * (1.0 - c) + ax[3] * s,
                                 ax[1] * ax[3] * (1.0 - c) - ax[2] * s,
                                 ax[1] * ax[2] * (1.0 - c) - ax[3] * s,
                                 ax[2]^2 + (1.0 - ax[2]^2) * c,
                                 ax[2] * ax[3] * (1.0 - c) + ax[1] * s,
                                 ax[1] * ax[3] * (1.0 - c) + ax[2] * s,
                                 ax[2] * ax[3] * (1.0 - c) - ax[1] * s,
                                 ax[3]^2 + (1.0 - ax[3]^2) * c)

    # check if rotation works as expected
    if norm(mat * a .- b) > 1e-8
        # check if sense of rotation has to be inverted
        if norm(transpose(mat) * a .- b) < 1e-8
            return transpose(mat)
        # if algorithm fails return matrix of zeros
        else
            return SMatrix{3, 3, Float64}(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
        end
    # if sanity check is passed return rotation matrix
    else
        return mat
    end
end

# try to obtain a rotation of (vec1, vec2) onto (ref1, ref2) (return matrix of zeros in case of failure)
function rotate(
    vec1 :: SVector{3, Float64},
    vec2 :: SVector{3, Float64},
    ref1 :: SVector{3, Float64},
    ref2 :: SVector{3, Float64}
    )    :: SMatrix{3, 3, Float64}

    # rotate vec1 onto ref1
    mat = get_rotation(vec1, ref1)

    # try to rotate vec2 around ref1 onto ref2
    mat = get_rotation_for_axis(mat * vec2, ref2, ref1) * mat

    return mat
end

"""
    get_trafos_orig(
        l :: Lattice
        ) :: Vector{SMatrix{3, 3, Float64}}

Compute transformations which leave the origin of the lattice invariant (point group symmetries).
"""
function get_trafos_orig(
    l :: Lattice
    ) :: Vector{SMatrix{3, 3, Float64}}

    # allocate list for transformations, set reference site and its bonds
    trafos = SMatrix{3, 3, Float64}[]
    ref    = l.sites[1]
    con    = l.uc.bonds[1]

    # get a pair of non-collinear neighbors of origin
    for i1 in eachindex(con)
        for i2 in eachindex(con)
            # check that sites are unequal
            if i1 == i2
                continue
            end

            # get connected sites in Bravais representation
            int_i1 = ref.int .+ con[i1]
            int_i2 = ref.int .+ con[i2]

            # get connected sites in real space representation
            vec_i1 = get_vec(int_i1, l.uc)
            vec_i2 = get_vec(int_i2, l.uc)

            # check that sites are non-collinear
            if norm(cross(vec_i1, vec_i2)) <= 1e-8
                continue
            end

            # get site structs
            site_i1 = Site(int_i1, vec_i1)
            site_i2 = Site(int_i2, vec_i2)

            # get bonds
            bond_i1 = get_bond(ref, site_i1, l)
            bond_i2 = get_bond(ref, site_i2, l)

            # get another pair of non-collinear neighbors of origin
            for j1 in eachindex(con)
                for j2 in eachindex(con)
                    # check that sites are unequal
                    if j1 == j2
                        continue
                    end

                    # get connected sites in Bravais representation
                    int_j1 = ref.int .+ con[j1]
                    int_j2 = ref.int .+ con[j2]

                    # get connected sites in real space representation
                    vec_j1 = get_vec(int_j1, l.uc)
                    vec_j2 = get_vec(int_j2, l.uc)

                    # check that sites are non-collinear
                    if norm(cross(vec_j1, vec_j2)) <= 1e-8
                        continue
                    end

                    # get site structs
                    site_j1 = Site(int_j1, vec_j1)
                    site_j2 = Site(int_j2, vec_j2)

                    # get bonds
                    bond_j1 = get_bond(ref, site_j1, l)
                    bond_j2 = get_bond(ref, site_j2, l)

                    # try to obtain rotation
                    mat = rotate(site_i1.vec, site_i2.vec, site_j1.vec, site_j2.vec)

                    # if successful, verify rotation on test set before saving it
                    if maximum(abs.(mat)) > 1e-8
                        # check if rotation is already known
                        if is_in(mat, trafos) == false
                            # verify rotation
                            valid = true

                            for n in eachindex(l.test_sites)
                                # apply transformation to lattice site
                                mapped_vec = mat * l.test_sites[n].vec
                                orig_bond  = get_bond(ref, l.test_sites[n], l)
                                mapped_ind = get_site(mapped_vec, l)

                                # check if resulting site is in lattice and if the bonds match
                                if mapped_ind == 0
                                    valid = false
                                else
                                    valid = are_equal(orig_bond, get_bond(ref, l.sites[mapped_ind], l))
                                end

                                # break if test fails for one test site
                                if valid == false
                                    break
                                end
                            end

                            # save transformation
                            if valid
                                push!(trafos, mat)
                            end
                        end

                        # check if rotation combined with inversion is already known
                        if is_in(-mat, trafos) == false
                            # verify rotation
                            valid = true

                            for n in eachindex(l.test_sites)
                                # apply transformation to lattice site
                                mapped_vec = -mat * l.test_sites[n].vec
                                orig_bond  = get_bond(ref, l.test_sites[n], l)
                                mapped_ind = get_site(mapped_vec, l)

                                # check if resulting site is in lattice and if the bonds match
                                if mapped_ind == 0
                                    valid = false
                                else
                                    valid = are_equal(orig_bond, get_bond(ref, l.sites[mapped_ind], l))
                                end

                                # break if test fails for one test site
                                if valid == false
                                    break
                                end
                            end

                            # save transformation
                            if valid
                                push!(trafos, -mat)
                            end
                        end
                    end
                end
            end
        end
    end

    return trafos
end

# compute reduced representation of the lattice
function get_reduced(
    l :: Lattice
    ) :: Vector{Int64}

    # allocate a list of indices
    reduced = Int64[i for i in eachindex(l.sites)]

    # get transformations
    trafos = get_trafos_orig(l)

    # iterate over sites and try to find sites which are symmetry equivalent
    for i in 2 : length(reduced)
        # continue if index has already been replaced
        if reduced[i] < i
            continue
        end

        # apply all available transformations to current site and see where it is mapped
        for j in eachindex(trafos)
            mapped_vec = trafos[j] * l.sites[i].vec

            # check that site is not mapped to itself
            if norm(mapped_vec .- l.sites[i].vec) > 1e-8
                # determine image of site
                index = get_site(mapped_vec, l)

                # assert that the image is valid, otherwise this is not a valid transformation and our algorithm failed
                @assert index > 0 "Validity on test set could not be generalized."

                # replace site in list if bonds match
                if are_equal(l.bonds[1, i], get_bond(l.sites[1], l.sites[index], l))
                    reduced[index] = i
                end
            end
        end
    end

    return reduced
end

"""
    get_trafos_uc(
        l :: Lattice
        ) :: Vector{Tuple{SMatrix{3, 3, Float64}, Bool}}

Compute mappings of a lattice's basis sites to the origin.
The mappings consist of a transformation matrix and a boolean indicating if an inversion was used or not.
"""
function get_trafos_uc(
    l :: Lattice
    ) :: Vector{Tuple{SMatrix{3, 3, Float64}, Bool}}

    # allocate list for transformations, set reference site and its bonds
    trafos = Vector{Tuple{SMatrix{3, 3, Float64}, Bool}}(undef, length(l.uc.basis) - 1)
    ref    = l.sites[1]
    con    = l.uc.bonds[1]

    # iterate over basis sites and find a symmetry for each of them
    for b in 2 : length(l.uc.basis)
        # set basis site and connections
        int       = SVector{4, Int64}(0, 0, 0, b)
        basis     = Site(int, get_vec(int, l.uc))
        con_basis = l.uc.bonds[b]

        # get a pair of non-collinear neighbors of basis
        for b1 in eachindex(con_basis)
            for b2 in eachindex(con_basis)
                # check that sites are unequal
                if b1 == b2
                    continue
                end

                # get connected sites in Bravais representation
                int_b1 = basis.int .+ con_basis[b1]
                int_b2 = basis.int .+ con_basis[b2]

                # get connected sites in real space representation
                vec_b1 = get_vec(int_b1, l.uc)
                vec_b2 = get_vec(int_b2, l.uc)

                # check that sites are non-collinear
                if norm(cross(vec_b1 .- basis.vec, vec_b2 .- basis.vec)) <= 1e-8
                    continue
                end

                # get site structs
                site_b1 = Site(int_b1, vec_b1)
                site_b2 = Site(int_b2, vec_b2)

                # get bonds
                bond_b1 = get_bond(basis, site_b1, l)
                bond_b2 = get_bond(basis, site_b2, l)

                # get a pair of non-collinear neigbors of origin
                for ref1 in eachindex(con)
                    for ref2 in eachindex(con)
                        # check that sites are unequal
                        if ref1 == ref2
                            continue
                        end

                        # get connected sites in Bravais representation
                        int_ref1 = ref.int .+ con[ref1]
                        int_ref2 = ref.int .+ con[ref2]

                        # get connected sites in real space representation
                        vec_ref1 = get_vec(int_ref1, l.uc)
                        vec_ref2 = get_vec(int_ref2, l.uc)

                        # check that sites are non-collinear
                        if norm(cross(vec_ref1, vec_ref2)) <= 1e-8
                            continue
                        end

                        # get site structs
                        site_ref1 = Site(int_ref1, vec_ref1)
                        site_ref2 = Site(int_ref2, vec_ref2)

                        # get bonds
                        bond_ref1 = get_bond(ref, site_ref1, l)
                        bond_ref2 = get_bond(ref, site_ref2, l)

                        # check that the bonds match
                        if are_equal(bond_b1, bond_ref1) == false || are_equal(bond_b2, bond_ref2) == false
                            continue
                        end

                        # try shift -> rotation
                        mat = rotate(site_b1.vec .- basis.vec, site_b2.vec .- basis.vec, site_ref1.vec, site_ref2.vec)

                        # if successful, verify tranformation on test set before saving it
                        if maximum(abs.(mat)) > 1e-8
                            valid = true

                            for n in eachindex(l.test_sites)
                                # test only those sites which are in range of basis
                                if get_metric(l.test_sites[n], basis, l.uc) <= l.size
                                    # apply transformation to lattice site
                                    mapped_vec = mat * (l.test_sites[n].vec .- basis.vec)
                                    orig_bond  = get_bond(basis, l.test_sites[n], l)
                                    mapped_ind = get_site(mapped_vec, l)

                                    # check if resulting site is in lattice and if bonds match
                                    if mapped_ind == 0
                                        valid = false
                                    else
                                        valid = are_equal(orig_bond, get_bond(ref, l.sites[mapped_ind], l))
                                    end

                                    # break if test fails for one test site
                                    if valid == false
                                        break
                                    end
                                end
                            end

                            # save transformation
                            if valid
                                trafos[b - 1] = (mat, false)
                                break
                                break
                                break
                                break
                            end
                        end

                        # try inversion -> shift -> rotation
                        mat = rotate(-(site_b1.vec .- basis.vec), -(site_b2.vec .- basis.vec), site_ref1.vec, site_ref2.vec)

                        # if successful, verify tranformation on test set before saving it
                        if maximum(abs.(mat)) > 1e-8
                            valid = true

                            for n in eachindex(l.test_sites)
                                # test only those sites which are in range of basis
                                if get_metric(l.test_sites[n], basis, l.uc) <= l.size
                                    # apply transformation to lattice site
                                    mapped_vec = mat * (-l.test_sites[n].vec .+ basis.vec)
                                    orig_bond  = get_bond(basis, l.test_sites[n], l)
                                    mapped_ind = get_site(mapped_vec, l)

                                    # check if resulting site is in lattice and if bonds match
                                    if mapped_ind == 0
                                        valid = false
                                    else
                                        valid = are_equal(orig_bond, get_bond(ref, l.sites[mapped_ind], l))
                                    end

                                    # break if test fails for one test site
                                    if valid == false
                                        break
                                    end
                                end
                            end

                            # save transformation
                            if valid
                                trafos[b - 1] = (mat, true)
                                break
                                break
                                break
                                break
                            end
                        end
                    end
                end
            end
        end
    end

    return trafos
end

# apply transformation (i, j) -> (i0, j*), where i0 is the origin, to site j
function apply_trafo(
    s      :: Site,
    b      :: Int64,
    shift  :: SVector{3, Float64},
    trafos :: Vector{Tuple{SMatrix{3, 3, Float64}, Bool}},
    l      :: Lattice
    )      :: SVector{3, Float64}

    # check if equivalent to origin
    if b != 1
        # if inequivalent to origin use transformation inside unitcell
        if trafos[b - 1][2]
            return trafos[b - 1][1] * (shift .- s.vec .+ l.uc.basis[b])
        else
            return trafos[b - 1][1] * (s.vec .- shift .- l.uc.basis[b])
        end
    # if equivalent to origin, only shift along translation vectors needs to be performed
    else
        return s.vec .- shift
    end
end

# compute mappings onto reduced lattice
function get_mappings(
    l       :: Lattice,
    reduced :: Vector{Int64}
    )       :: Matrix{Int64}

    # allocate matrix
    num = length(l.sites)
    mat = Matrix{Int64}(I, num, num)

    # get transformations inside unitcell
    trafos = get_trafos_uc(l)

    # get distances to origin
    dists = Float64[norm(l.sites[i].vec) for i in eachindex(l.sites)]
    d_max = maximum(dists)

    # group sites in shells
    shell_dists = unique(trunc.(dists, digits = 8))
    shells      = Vector{Vector{Int64}}(undef, length(shell_dists))

    for i in eachindex(shells)
        shell = Int64[]

        for j in eachindex(dists)
            if abs(dists[j] - shell_dists[i]) < 1e-8
                push!(shell, j)
            end
        end

        shells[i] = shell
    end

    # compute entries of matrix
    Threads.@threads for i in eachindex(l.sites)
        # determine shift for second site
        si    = l.sites[i]
        b     = si.int[4]
        shift = get_vec(SVector{4, Int64}(si.int[1], si.int[2], si.int[3], 1), l.uc)

        for j in eachindex(l.sites)
            # compute only off-diagonal entries
            if i != j
                # apply symmetry transformation
                mapped_vec = apply_trafo(l.sites[j], b, shift, trafos, l)

                # locate transformed site in lattice
                index = 0
                d_map = norm(mapped_vec)

                # check if transformed site can be in lattice
                if d_max - d_map > -1e-8
                    # find respective shell
                    shell = shells[argmin(abs.(shell_dists .- d_map))]

                    # find matching site in shell
                    for k in shell
                        if norm(mapped_vec .- l.sites[k].vec) < 1e-8
                            index = k
                            break
                        end
                    end

                    if index != 0
                        mat[i, j] = reduced[index]
                    end
                end
            end
        end
    end

    return mat
end

# compute irreducible sites in overlap of two sites
function get_overlap(
    l           :: Lattice,
    reduced     :: Vector{Int64},
    irreducible :: Vector{Int64},
    mappings    :: Matrix{Int64}
    )           :: Vector{Matrix{Int64}}

    # allocate overlap
    overlap = Vector{Matrix{Int64}}(undef, length(irreducible))

    for i in eachindex(irreducible)
        # collect all sites in range of irreducible and origin
        temp = NTuple{2, Int64}[]

        for j in eachindex(l.sites)
            if mappings[irreducible[i], j] != 0
                push!(temp, (reduced[j], mappings[j, irreducible[i]]))
            end
        end

        # determine how often a certain pair occurs
        pairs = unique(temp)
        table = zeros(Int64, length(pairs), 3)

        for j in eachindex(pairs)
            # convert from original lattice index to new "irreducible" index
            table[j, 1] = findall(index -> index == pairs[j][1], irreducible)[1]
            table[j, 2] = findall(index -> index == pairs[j][2], irreducible)[1]

            # count multiplicity
            for k in eachindex(temp)
                if pairs[j] == temp[k]
                    table[j, 3] += 1
                end
            end
        end

        # save into overlap
        overlap[i] = table
    end

    return overlap
end

# compute multiplicity of irreducible sites
function get_mult(
    reduced     :: Vector{Int64},
    irreducible :: Vector{Int64}
    )           :: Vector{Int64}

    # allocate vector to store multiplicities
    mult = zeros(Float64, length(irreducible))

    for i in eachindex(mult)
        for j in eachindex(reduced)
            if reduced[j] == irreducible[i]
                mult[i] += 1
            end
        end
    end

    return mult
end

# compute mappings under site exchange
function get_exchange(
    irreducible :: Vector{Int64},
    mappings    :: Matrix{Int64}
    )           :: Vector{Int64}

    # allocate exchange list
    exchange = zeros(Int64, length(irreducible))

    # determine mappings under site exchange
    for i in eachindex(exchange)
        exchange[i] = findall(index -> index == mappings[irreducible[i], 1], irreducible)[1]
    end

    return exchange
end

# convert mapping table entries to irreducible site indices
function get_project(
    l           :: Lattice,
    irreducible :: Vector{Int64},
    mappings    :: Matrix{Int64}
    )           :: Matrix{Int64}

    # allocate projections
    project = zeros(Int64, length(l.sites), length(l.sites))

    for i in eachindex(l.sites)
        for j in eachindex(l.sites)
            if mappings[i, j] != 0
                project[i, j] = findall(index -> index == mappings[i, j], irreducible)[1]
            end
        end
    end

    return project
end

"""
    get_reduced_lattice(
        model   :: String,
        J       :: Vector{Vector{Float64}},
        l       :: Lattice
        ;
        verbose :: Bool = true
        )       :: Reduced_lattice

Compute symmetry reduced representation of a given lattice graph with spin interactions between sites.
The interactions are defined by passing a model's name and coupling vector.
"""
function get_reduced_lattice(
    model   :: String,
    J       :: Vector{Vector{Float64}},
    l       :: Lattice
    ;
    verbose :: Bool = true
    )       :: Reduced_lattice

    if verbose
        println("Performing symmetry reduction ...")
    end

    # initialize model by modifying bond matrix of lattice
    init_model!(model, J, l)

    # get reduced representation of lattice
    reduced     = get_reduced(l)
    irreducible = unique(reduced)
    sites       = Site[Site(l.sites[i].int, get_vec(l.sites[i].int, l.uc)) for i in irreducible]

    # get mapping table
    mappings = get_mappings(l, reduced)

    # get overlap
    overlap = get_overlap(l, reduced, irreducible, mappings)

    # get multiplicities
    mult = get_mult(reduced, irreducible)

    # get exchanges
    exchange = get_exchange(irreducible, mappings)

    # get projections
    project = get_project(l, irreducible, mappings)

    # build reduced lattice
    r = Reduced_lattice(l.name, l.size, model, J, sites, overlap, mult, exchange, project)

    if verbose
        println("Done. Reduced lattice has $(length(r.sites)) sites.")
    end

    return r
end


text:
function entropy(d::Truncated{Normal})
    s = std(d.untruncated)
    a = d.lower
    b = d.upper
    phi_a = pdf(d.untruncated, a) * s
    phi_b = pdf(d.untruncated, b) * s
    a_phi_a = a == -Inf ? 0.0 : a * phi_a
    b_phi_b = b == Inf ? 0.0 : b * phi_b
    z = d.nc
    return entropy(d.untruncated) + log(z) +
           0.5 * (a_phi_a - b_phi_b) / z - 0.5 * ((phi_a - phi_b) / z)^2
end

function mean(d::Truncated{Normal})
    delta = pdf(d.untruncated, d.lower) - pdf(d.untruncated, d.upper)
    return mean(d.untruncated) + delta * var(d.untruncated) / d.nc
end

function modes(d::Truncated{Normal})
    mu = mean(d.untruncated)
    if d.upper < mu
        return [d.upper]
    elseif d.lower > mu
        return [d.lower]
    else
        return [mu]
    end
end

function rand(d::Truncated{Normal})
    mu = mean(d.untruncated)
    sigma = std(d.untruncated)
    z = randnt((d.lower - mu) / sigma, (d.upper - mu) / sigma)
    return mu + sigma * z
end

function var(d::Truncated{Normal})
    s = std(d.untruncated)
    a = d.lower
    b = d.upper
    phi_a = pdf(d.untruncated, a) * s
    phi_b = pdf(d.untruncated, b) * s
    a_phi_a = a == -Inf ? 0.0 : a * phi_a
    b_phi_b = b == Inf ? 0.0 : b * phi_b
    z = d.nc
    return s^2 * (1 + (a_phi_a - b_phi_b) / z - ((phi_a - phi_b) / z)^2)
end

# Rejection sampler based on algorithm from Robert (1992)
#  - Available at http://arxiv.org/abs/0907.4010
function randnt(lower::Real, upper::Real)
    if (lower <= 0 && upper == Inf) ||
       (upper >= 0 && lower == Inf) ||
       (lower <= 0 && upper >= 0 && upper - lower > √2π)
        while true
            r = randn()
            if r > lower && r < upper
                return r
            end
        end
    elseif lower > 0 && upper - lower > 2.0 / (lower + sqrt(lower^2 + 4.0)) * exp((lower^2 - lower * sqrt(lower^2 + 4.0)) / 4.0)
        a = (lower + sqrt(lower^2 + 4.0))/2.0
        while true
            r = rand(Exponential(1.0 / a)) + lower
            u = rand()
            if u < exp(-0.5 * (r - a)^2) && r < upper
                return r
            end
        end
    elseif upper < 0 && upper - lower > 2.0 / (-upper + sqrt(upper^2 + 4.0)) * exp((upper^2 + upper * sqrt(upper^2 + 4.0)) / 4.0)
        a = (-upper + sqrt(upper^2 + 4.0)) / 2.0
        while true
            r = rand(Exponential(1.0 / a)) - upper
            u = rand()
            if u < exp(-0.5 * (r - a)^2) && r < -lower
                return -r
            end
        end
    else
        while true
            r = lower + rand() * (upper - lower)
            u = rand()
            if lower > 0
                rho = exp((lower^2 - r^2) * 0.5)
            elseif upper < 0
                rho = exp((upper^2 - r^2) * 0.5)
            else
                rho = exp(-r^2 * 0.5)
            end
            if u < rho
                return r
            end
        end
    end
end

ext: jl


text:
################################# JSON Parsers #################################
JSON2.@format Range keywordargs begin
    # lb => (default=0.0,)
    # ub => (default=MAXNUMBER,)
end

JSON2.@format Flexibility keywordargs begin
    # # flexibility_level => (default=0,)
    # fixed_price => (default=0.0,)
    # unit_price => (default=0.0,)
end

JSON2.@format FlexibleRange keywordargs begin
    # # soft_range => (default=Range(),)
    # hard_range => (default=Range(),)
    # lb_flex => (default=Flexibility(),)
    # ub_flex => (default=Flexibility(),)
end

JSON2.@format Location keywordargs begin
    # # id => (default="",)
    # # index => (default=-1,)
    # lat_y => (default=MAXNUMBER,)
    # long_x => (default=MAXNUMBER,)
    # opening_time_windows => (default=[Range()],)
    # energy_fixed_cost => (default=0.0,)
    # energy_unit_cost => (default=0.0,)
    # energy_recharging_speeds => (default=Float64[],)
end

# JSON2.@format LocationGroup keywordargs begin
#     # id => (default="",)
#     # location_ids => (default=String[],)
# end

JSON2.@format ProductCompatibilityClass keywordargs begin
    # # id => (default="",)
    # conflict_compatib_class_ids => (default=String[],)
    # prohibited_predecessor_compatib_class_ids => (default=String[],)
end

JSON2.@format ProductSharingClass keywordargs begin
    # # id => (default="",)
    # infinite_pickup_availabitilies => (default=true,)
    # infinite_delivery_availabitilies => (default=true,)
    # pickup_availabitilies_at_location_ids => (default=Dict{String,Float64}(),)
    # delivery_capacities_at_location_ids => (default=Dict{String,Float64}(),)
end

JSON2.@format ProductSpecificationClass keywordargs begin
    # # id => (default="",)
    # capacity_consumptions => (default=Dict{String,Tuple{Float64,Float64}}(),)
    # property_requirements => (default=Dict{String,Float64}(),)
end

JSON2.@format Request keywordargs begin
    # # id => (default="",)
    # product_compatibility_class_id => (default="default_id",)
    # product_sharing_class_id => (default="default_id",)
    # product_specification_class_id => (default="default_id",)
    # split_fulfillment => (default=false,)
    # request_flexibility => (default=Flexibility(),)
    # precedence_status => (default=0,)
    # product_quantity_range => (default=Range(lb = 1, ub = 1),)
    # pickup_location_group_id => (default="default_id",)
    # delivery_location_group_id => (default="default_id",)
    # pickup_service_time => (default=0.0,)
    # delivery_service_time => (default=0.0,)
    # max_duration => (default=MAXNUMBER,)
    # duration_unit_cost => (default=0.0,)
    # pickup_time_windows => (default=[FlexibleRange()],)
    # delivery_time_windows => (default=[FlexibleRange()],)
end

JSON2.@format VehicleCharacteristics keywordargs begin
    # of_vehicle => (default=Dict{String,Float64}(),)
    # of_compartments => (default=Dict{String,Dict{String,Float64}}(),)
end

JSON2.@format VehicleCategory keywordargs begin
    # # id => (default="",)
    # capacities => (default=VehicleCharacteristics(),)
    # properties => (default=VehicleCharacteristics(),)
    # loading_option => (default=0,)
    # energy_interval_lengths => (default=Float64[],)
end

JSON2.@format CostPeriod keywordargs begin
    # period => (default=Range(),)
    # travel_distance_unit_cost => (default=0.0,)
    # travel_time_unit_cost => (default=0.0,)
    # service_time_unit_cost => (default=0.0,)
    # waiting_time_unit_cost => (default=0.0,)
    # fixed_cost => (default=0.0,)
end

JSON2.@format HomogeneousVehicleSet keywordargs begin
    # # id => (default="",)
    # vehicle_category_id => (default="default_id",)
    # # departure_location_group_id => (default="",)
    # # arrival_location_group_id => (default="",)
    # cost_periods => (default=[CostPeriod()],)
    # work_periods => (default=[FlexibleRange()],)
    # initial_energy_charge => (default=MAXNUMBER,)
    # max_working_time => (default=MAXNUMBER,)
    # max_travel_distance => (default=MAXNUMBER,)
    # allow_shipment_over_multiple_work_periods => (default=true,)
    # nb_of_vehicles_range => (default=FlexibleRange(),)
end

JSON2.@format RvrpInstance keywordargs begin
    # # id => (default="",)
    # distance_mode => (default=0,)
    # coordinate_mode => (default=0,)
    # # travel_matrix_periods => (default=[Range()],)
    # # period_to_matrix_id => (default=Dict{Range,String}(Range() => "default_mat"),)
    # travel_time_matrices => (default=Dict{String,Array{Float64,2}}(),)
    # travel_distance_matrices => (default=Dict{String,Array{Float64,2}}(),)
    # energy_consumption_matrices => (default=Dict{String,Array{Float64,2}}(),)
    # locations => (default=Location[],)
    # location_groups => (default=LocationGroup[],)
    # product_compatibility_classes => (default=ProductCompatibilityClass[],)
    # product_sharing_classes => (default=ProductSharingClass[],)
    # product_specification_classes => (default=ProductSpecificationClass[],)
    # requests => (default=Request[],)
    # vehicle_categories => (default=VehicleCategory[],)
    # vehicle_sets => (default=HomogeneousVehicleSet[],)
end

function parse_to_json(data::RvrpInstance, file_path::String)
    json2_string = JSON2.write(data)
    io = open(file_path, "w")
    write(io, json2_string * "\n")
    close(io)
end

parse_from_json_string(s::String) = JSON2.read(s, RvrpInstance)

function parse_from_json(file_path::String)
    io = open(file_path, "r")
    s = read(io, String)
    close(io)
    return parse_from_json_string(s)
end


text:
using NumericalMethodsforEngineers

A = Float64[0 0 16; 0 4 5; 8 -4 22]
b = Float64[4, 2, 5]

println("\nCholesky factorization using banded storage\n")
println("\nBanded Coefficient Matrix:\n$A\n")
println("\nRight Hand side:\n$b\n")

cholin!(A)
println("\nL in Band Form:\n$A\n")

chobac!(A, b)
println("\nSolution Vector:\n$b")

A = Float64[0 0 16; 0 4 5; 8 -4 22]
b = Float64[4, 2, 5]

c = A \ b

@test round(c, 14) == b


text:
@doc """  if pos_gene have this position, return gene; else return ""
          Foundational function
""" ->
function query_gene(chr::ASCIIString,pos::Int64)
    info("query_gene for a $chr,$pos")
    if !isdefined(GeneMisc,:chr_sted_gene)
        warn("chr_sted_gene is not defined,load_index starts")
        load_index()
    end
    sted_gene = chr_sted_gene[chr]

    steds = collect(keys(sted_gene))
    # sorted by acclerator speed
    starts = sort(map(x->x[1], steds))
    ends   = sort(map(x->x[2], steds))
    debug = true
    function check_sted(starts,ends)
        l = length(starts)
        @assert length(ends) == l
        for i in 1:l
            @assert starts[i] < ends[i]
        end
    end
    if debug check_sted(starts,ends) end
    idx_st = searchsortedfirst(starts, pos)
    idx_ed = searchsortedfirst(ends, pos)
    function result(st_idx,ed_idx)
        info("$st_idx, $ed_idx")
        st = starts[st_idx]
        ed = ends[ed_idx]
        sted_gene[(st,ed)]
    end
    if idx_st == 1
        if starts[1] < pos
            return ""
        end
        @assert starts[1] >= pos
        return result(1,1)
    end
    if idx_st == idx_ed
        #if ends[idx_ed] > pos
        return ""
        #end
        #return result(idx_st)
    else
        @assert idx_st > idx_ed
        return result(idx_st-1,idx_ed)
    end

    #=
    for (st,ed) in steds
        if st<=pos<=ed
            return sted_gene[(st,ed)]
        end
    end
    =#

    ""
end

@doc """ if pos_gene have this position, return gene; else return ""
""" ->
function query_gene(chr::ASCIIString,pos::ASCIIString)
    pos_int = parse(Int64,pos)
    query_gene(chr,pos_int)
end

@doc """ query a lot of positions
""" ->
function query_gene(data::Array{ASCIIString,2})
    info("query gene location for you ...")
    @assert size(data,2) == 2
    n = size(data,1)
    genes = Array{ASCIIString,1}(n)
    for i = 1:n
        chr,pos = data[i,:]
        genes[i] = query_gene(chr,pos)
    end
    genes
end

@doc """ Given a gene name, find its location: chr, st, ed
""" ->
function query_geneloc(genename::ASCIIString)
    if !isdefined(GeneMisc, :gene_chrsted)
        load_index()
    end
    gene_chrsted[genename]
end


text:
# Example 5.4 of
# A. Ahmadi, R. Jungers, P. Parrilo and M. Roozbehani,
# "Joint spectral radius and path-complete graph Lyapunov functions."
# SIAM J. CONTROL OPTIM 52(1), 687-717, 2014.
# The JSR is 3.917384715148

@testset "[AJPR14] Example 5.4" begin
    s = discreteswitchedsystem([[-1 -1; -4 0], [3 3; -2 1]])
    @test sprint(show, s) == "Hybrid System with automaton OneStateAutomaton(2)"
    qub = 4.31959610746
    @test quicklb(s) ≈ 3
    @test quickub(s) ≈ qub
    qb = quickb(s)
    @test typeof(qb) == NTuple{2,Float64}
    @test qb[1] ≈ 3
    @test qb[2] ≈ qub
    # pnorm=Inf: [17.2409,7.42377,5.65853,4.92192,4.52597,4.28099,4.11275,3.99139,3.8991,3.827,3.76891,3.72121,3.6813,3.64744,3.61834,3.59307,3.57092,3.55135,3.53393,3.51832,3.50426,3.49152,3.47994]
    # pnorm=1:   [15.1245,7.35044,5.47402,4.84722,4.45847,4.23242,4.0709,3.95644,3.8687,3.80015,3.74485,3.69944,3.66141,3.62914,3.6014,3.5773,3.55616,3.53748,3.52085,3.50595,3.49253,3.48037,3.4693,3.45919]
    # pnorm=2:   [14.1535,6.81871,5.32651,4.7096,4.37067,4.15723,4.01117,3.90488,3.82407,3.7606,3.70945,3.66735,3.6321,3.60215,3.57639,3.554,3.53437,3.517,3.50154,3.48768,3.47519]

    @test pradius(s, 2) ≈ 3.234535151244 rtol=1e-9
    @test pradius(s, 2, BruteForce(), pnorm=2) ≈ 3.632097274822649

    pρlb = [3.23453515151, 3.42756015606, 3.54374928704]
    pρub = [4.57432347886, 4.07607892469, 3.97772408342]
    @testset "p-radius with $algo" for algo in (VeroneseLift(), KroneckerLift())
        for d in 1:3
            lb, ub = pradiusb(s, 2d, algo)
            @test isapprox(lb, pρlb[d])
            @test isapprox(ub, pρub[d])
        end
    end

    tmp = periodicswitching(s, [1, 2])
    tmp = periodicswitching(tmp.s, tmp.period, round(tmp.growthrate, digits=4)) # Avoid difference of floating point rounding
    @test sprint(show, tmp) == "PSW(3.9174, [1, 2])" #"Periodic switching of growth rate 3.9174 and modes: [1, 2] for the transitions [1, 2]"
    smp = periodicswitching(s, [1, 2])
    #@test smp != tmp
    @test smp != periodicswitching(s, [2, 1, 2, 1])
    @test smp != periodicswitching(s, [1, 2, 1, 2])
    @test smp != periodicswitching(s, [1, 1])
    #@test smp != periodicswitching(discreteswitchedsystem([[-1 -1; -4 0],[3 3; -2 1]]), [1, 2])
    @testset "JSR with $factory" for factory in sdp_factories
        sosdata(s).lb = 0
        tol = ismosek(factory) ? 1e-5 : 5e-4
        lb, ub = soslyapb(s, 1, factory=factory, tol=tol)
        @test log(lb) ≈ log(2.814640557) rtol=tol
        @test log(ub) ≈ log(3.980502849) rtol=tol
        lb, ub = soslyapb(s, 2, factory=factory, tol=tol)
        @test log(lb) ≈ log(3.299750624) rtol=tol
        @test log(ub) ≈ log(3.924086919) rtol=tol
        @test_throws ArgumentError sosbuildsequence(s, 1, v_0 = 2)
        for d in (1, 2)
            for p_0 in (:Primal, :Random)
                seq = sosbuildsequence(s, d, p_0=p_0)
                psw = findsmp(seq)
                @test psw !== nothing
                @test psw == smp
            end
        end
    end
    if isempty(sdp_factories)
        @test getub(s) ≈ pρub[end]
    else
        @test log(getub(s)) ≈ log(3.924086919) rtol=5e-4
        @test getsmp(s) == smp
        @test getlb(s) ≈ 3.917384715148 rtol=1e-12
    end
end


text:
"""
$(TYPEDEF)

Message prior on all dimensions of a variable node in the factor graph.

Notes
- Only temporary existance during CSM operations.
"""
struct MsgPrior{T <: SamplableBelief} <: AbstractPrior
  Z::T
  infoPerCoord::Vector{Float64}
  M
end

# MsgPrior{T}() where {T} = new{T}()
# MsgPrior{T}(z::T, infd::R) where {T <: SamplableBelief, R <: Real} = new{T}(z, infd)
# function MsgPrior(z::T, infd::R) where {T <: SamplableBelief, R <: Real}
#     MsgPrior{T}(z, infd)
# end
function getSample(cf::CalcFactor{<:MsgPrior})
  return rand(cf.factor.Z, 1)
end

#TODO check these for manifolds, may need updating to samplePoint
# MKD already returns a vector of points
function getSample(cf::CalcFactor{<:MsgPrior{<:ManifoldKernelDensity}})
  mkd = cf.factor.Z
  return samplePoint(mkd.manifold, mkd)
end

getManifold(mp::MsgPrior{<:ManifoldKernelDensity}) = mp.Z.manifold
getManifold(mp::MsgPrior) = mp.M

(cfo::CalcFactor{<:MsgPrior})(z, x1) = z .- x1

Base.@kwdef struct PackedMsgPrior <: AbstractPackedFactor
  Z::PackedSamplableBelief
  infoPerCoord::Vector{Float64}
end

function convert(::Type{PackedMsgPrior}, d::MsgPrior)
  PackedMsgPrior(convert(PackedSamplableBelief, d.Z), d.infoPerCoord)
end
function convert(::Type{<:MsgPrior}, d::PackedMsgPrior)
  MsgPrior(convert(SamplableBelief, d.Z), d.infoPerCoord)
end


text:
require("deCasteljau.jl")

function exercise_one ()
    controlPoints = [
                     [1.0 1 0];
                     [3.0 4 2];
                     [5.0 6 1];
                     [7.0 8 9];
                     [10.0 2 8];
                     [1.0 1 0];
                     ]

    bezierCurvePoints = drawCurve(controlPoints,
				  linspace(0,1,200))

    writeArrayForGnuplot(bezierCurvePoints, "bezier-exercise-one.coordinates")
end

function exercise_three()
    x = t -> 1 + t + t^2
    y = t -> t^3
    points = zeros(200, 2)
    params = linspace(0,1,200)
    for i=1:200
        t = params[i]
        points[i,:] =
            [x(t) y(t)]
    end
    writeArrayForGnuplot(points, "test-exercise-three.coordinates")
end

function exercise_two ()

    maximaControlPoints = [
                           [1.0 0];
                           [4/3 0];
                           [2.0 0];
                           [3.0 1];
                           ]

    maximaBezierCurvePoints = drawCurve(maximaControlPoints,
				  linspace(0,1,200))

    writeArrayForGnuplot(maximaControlPoints,
                         "control-poly-exercise-three-maxima.coordinates")

    writeArrayForGnuplot(maximaBezierCurvePoints,
                         "bezier-exercise-three-maxima.coordinates")

    controlPoints = parametricCoordinates([
			                   x -> 1 + x + x^2;
			                   y -> y^3
		                           ], 200)

    ## bezierCurvePoints = drawCurve(controlPoints,
    ##     			  linspace(0,1,200))

    ## writeArrayForGnuplot(bezierCurvePoints, "bezier-exercise-two.coordinates")
    writeArrayForGnuplot(controlPoints, "parametric-curve-exercise-three.coordinates")

    controlPoints = parametricCoordinates([
        		                   x -> 1 + x + x^2;
        		                   y -> y^3;
        		                   z -> z^2 - pi/2
        	                           ], 6)

    bezierCurvePoints = drawCurve(controlPoints,
        			  linspace(0,1,200))

    writeArrayForGnuplot(bezierCurvePoints, "bezier-three-axes-exercise-two.coordinates")
    writeArrayForGnuplot(controlPoints, "control-three-axes-poly-exercise-two.coordinates")
end

function exercise_four ()
    outerControlPoints = [
                          [0.0 0];
                          [1.0 3];
                          [2.5 4];
                          [5 4.5];
                          [4.5 2.5];
                          [6.0 1];
                          ]
    writeArrayForGnuplot(outerControlPoints, "control-poly-exercise-four.coordinates")
    t_star = .25
    q_value, Diagonal, Last = decasteljau(outerControlPoints, t_star)
    writeArrayForGnuplot(Diagonal, "cyan-control-poly-exercise-four.coordinates")
    writeArrayForGnuplot(Last, "magenta-control-poly-exercise-four.coordinates")

    #diagonalBezierCurvePoints = drawCurve(Diagonal, [t/t_star for t in linspace(0,1,100)])
    diagonalBezierCurvePoints = drawCurve(Diagonal, linspace(0,1,100))
    lastBezierCurvePoints = drawCurve(Last, linspace(0,1,100))
    writeArrayForGnuplot(diagonalBezierCurvePoints, "bezier-diagonal-exercise-four.coordinates")
    writeArrayForGnuplot(lastBezierCurvePoints, "bezier-last-exercise-four.coordinates")
end

function exercise_six()
    originalControlPoints = [
                             [1.0 1];
                             [3.0 4];
                             [5.0 6];
                             [7.0 1.5];
                             [10 2]
                             ]
    writeArrayForGnuplot(originalControlPoints, "control-poly-exercise-six.coordinates")
    originalBezierCurvePoints = drawCurve(originalControlPoints, linspace(0,1,200))
    writeArrayForGnuplot(originalBezierCurvePoints, "bezier-original-exercise-six.coordinates")

    oneMoreDegreeControlPoints = increaseDegree(originalControlPoints)
    oneMoreDegreeBezierCurvePoints = drawCurve(oneMoreDegreeControlPoints, linspace(0,1,200))
    writeArrayForGnuplot(oneMoreDegreeControlPoints, "control-poly-exercise-six-one-more-degree.coordinates")
    writeArrayForGnuplot(oneMoreDegreeBezierCurvePoints, "bezier-one-more-degree-exercise-six.coordinates")

    twoMoreDegreeControlPoints = increaseDegree(oneMoreDegreeControlPoints)
    twoMoreDegreeBezierCurvePoints = drawCurve(twoMoreDegreeControlPoints, linspace(0,1,200))
    writeArrayForGnuplot(twoMoreDegreeControlPoints, "control-poly-exercise-six-two-more-degree.coordinates")
    writeArrayForGnuplot(twoMoreDegreeBezierCurvePoints, "bezier-two-more-degree-exercise-six.coordinates")

    threeMoreDegreeControlPoints = increaseDegree(twoMoreDegreeControlPoints)
    threeMoreDegreeBezierCurvePoints = drawCurve(threeMoreDegreeControlPoints, linspace(0,1,200))
    writeArrayForGnuplot(threeMoreDegreeControlPoints, "control-poly-exercise-six-three-more-degree.coordinates")
    writeArrayForGnuplot(threeMoreDegreeBezierCurvePoints, "bezier-three-more-degree-exercise-six.coordinates")
end

function exercise_five()
    PointToRepeat = Float64[10 1]

    originalControlPoints = [
		             [2.0 4];
		             [6 12];
		             PointToRepeat;
		             [12 12]
	                     ]

    repeatedControlPoints = [
		             [2.0 4];
		             [6 12];
		             PointToRepeat;
		             PointToRepeat;
		             PointToRepeat;
		             PointToRepeat;
		             [12 12]
	                     ]

    originalBezierCurvePoints = drawCurve(originalControlPoints, linspace(0,1,200))
    repeatedBezierCurvePoints = drawCurve(repeatedControlPoints, linspace(0,1,200))
    writeArrayForGnuplot(originalBezierCurvePoints, "bezier-original-exercise-five.coordinates")
    writeArrayForGnuplot(repeatedBezierCurvePoints, "bezier-repeated-exercise-five.coordinates")
end

function exercise_seven()

    controlPoints = [
		     [0.0 0 0];
		     [1.0 2 0];
		     [3.0 2 0];
		     [6.0 -1 0];
	             ]

    continuityPoint, tangentPoint, obsculatorPoint, a_succ_i =
        junction_right(controlPoints, 1, 1)

    restControlPoints = [
                         [15.0 3 0];
		         [15.0 3 0];
		         [15.0 3 0];
		         [15.0 3 0];
		         [23.0 -4 0];
		         [24.0 0 0];
	                 ]

    continuityControlPoints = [
                               continuityPoint;
                               restControlPoints;
                               ]

    tangentControlPoints = [
                            continuityPoint;
                            tangentPoint;
                            restControlPoints;
                            ]

    obsculatingControlPoints = [
                                continuityPoint;
                                tangentPoint;
                                obsculatorPoint;
                                restControlPoints;
                                ]

    baseBezierCurvePoints = drawCurve(controlPoints,
                                      linspace(0,1,200))
    continuityBezierCurvePoints = drawCurve(continuityControlPoints,
                                            linspace(0,1,200))
    tangentBezierCurvePoints = drawCurve(tangentControlPoints,
                                         linspace(0,1,200))
    obsculatingBezierCurvePoints = drawCurve(obsculatingControlPoints,
                                             linspace(0,1,200))

    writeArrayForGnuplot(controlPoints,
                         "control-poly-base-exercise-seven-exercise-seven.coordinates")
    writeArrayForGnuplot(continuityControlPoints,
		         "control-poly-continuity-exercise-seven.coordinates")
    writeArrayForGnuplot(tangentControlPoints,
		         "control-poly-tangent-exercise-seven.coordinates")
    writeArrayForGnuplot(obsculatingControlPoints,
		         "control-poly-obsculating-exercise-seven.coordinates")

    writeArrayForGnuplot(baseBezierCurvePoints, "bezier-base-exercise-seven.coordinates")
    writeArrayForGnuplot(continuityBezierCurvePoints, "bezier-continuity-exercise-seven.coordinates")
    writeArrayForGnuplot(tangentBezierCurvePoints, "bezier-tangent-exercise-seven.coordinates")
    writeArrayForGnuplot(obsculatingBezierCurvePoints, "bezier-obsculating-exercise-seven.coordinates")

    writeArrayForGnuplot([
                          controlPoints;
                          obsculatingControlPoints;
                          ],
		         "control-poly-base-obsculating-exercise-seven.coordinates")

    writeArrayForGnuplot([
                          controlPoints[end-1,:];
                          a_succ_i;
                          obsculatingControlPoints[2,:];
                          ],
		         "control-poly-base-obsculating-with-a-succ-i-exercise-seven.coordinates")

end

function exercise_seven_left()

    controlPoints = [
		     [10.0 12 0];
		     [13.0 8.5 0];
		     [16.0 -2 0];
		     [26.0 0 0];
	             ]

    continuityPoint, tangentPoint, obsculatorPoint, a_i =
        junction_left(controlPoints, 1, 1)

    restControlPoints = [
		         [0.0 3 0];
		         [1.0 -4 0];
		         [5.0 0 0];
	                 ]

    continuityControlPoints = [
                               restControlPoints;
                               continuityPoint;
                               ]

    tangentControlPoints = [
                            restControlPoints;
                            tangentPoint;
                            continuityPoint;
                            ]

    obsculatingControlPoints = [
                                restControlPoints;
                                obsculatorPoint;
                                tangentPoint;
                                continuityPoint;
                                ]

    baseBezierCurvePoints = drawCurve(controlPoints,
                                      linspace(0,1,200))
    continuityBezierCurvePoints = drawCurve(continuityControlPoints,
                                            linspace(0,1,200))
    tangentBezierCurvePoints = drawCurve(tangentControlPoints,
                                         linspace(0,1,200))
    obsculatingBezierCurvePoints = drawCurve(obsculatingControlPoints,
                                             linspace(0,1,200))

    writeArrayForGnuplot(controlPoints,
                         "control-poly-base-left-exercise-seven-exercise-seven.coordinates")
    writeArrayForGnuplot(continuityControlPoints,
		         "control-poly-continuity-left-exercise-seven.coordinates")
    writeArrayForGnuplot(tangentControlPoints,
		         "control-poly-tangent-left-exercise-seven.coordinates")
    writeArrayForGnuplot(obsculatingControlPoints,
		         "control-poly-obsculating-left-exercise-seven.coordinates")

    writeArrayForGnuplot(baseBezierCurvePoints, "bezier-base-left-exercise-seven.coordinates")
    writeArrayForGnuplot(continuityBezierCurvePoints, "bezier-continuity-left-exercise-seven.coordinates")
    writeArrayForGnuplot(tangentBezierCurvePoints, "bezier-tangent-left-exercise-seven.coordinates")
    writeArrayForGnuplot(obsculatingBezierCurvePoints, "bezier-obsculating-left-exercise-seven.coordinates")

    writeArrayForGnuplot([
                          obsculatingControlPoints;
                          controlPoints;
                          ],
		         "control-poly-base-obsculating-left-exercise-seven.coordinates")

    writeArrayForGnuplot([
                          obsculatingControlPoints[end-1,:];
                          a_i;
                          controlPoints[2,:];
                          ],
		         "control-poly-base-obsculating-with-a-i-left-exercise-seven.coordinates")

end

function exercise_derivative()
    controlPoints = [
                     [1.0 1 0];
                     [3.0 4 2];
                     [5.0 6 1];
                     [7.0 8 9];
                     [10.0 2 8];
                     [1.0 1 0]
                     ]

    params = 200
    differencePoints = difference(controlPoints)
    derivativePoints = zeros(params, 3)

    j = 1
    for i=linspace(0,1,200)
        derivativePoints[j,:] = derivative(differencePoints,i)
        j += 1
    end

    writeArrayForGnuplot(drawCurve([
                                    [0.0 0 0];
                                    differencePoints;
                                    [0.0 0 0];
                                    ],
                                   linspace(0,1,200)),
		         "differences-of-exercise-one.coordinates")

    writeArrayForGnuplot(derivativePoints,
		         "derivative-of-exercise-one.coordinates")

end

function exercise_polar()
    controlPoints = [
                     [0.0 0.0];
                     [0.0 4.0];
                     [2.0 4.0];
                     [1.0 2.0];
                     ]

    n = length(controlPoints[:,1])
    params = 200
    differencePoints = difference(controlPoints)
    derivativePoints = zeros(params, 2)
    polarPoints = zeros(params, 2)
    bezierPoints = zeros(params, 2)

    j = 1
    for i=linspace(0,1,params)
        derivativePoints[j,:] = derivative(differencePoints,i)
        aPoint, Diagonal, Last = decasteljau(controlPoints, i)
        bezierPoints[j,:] = aPoint
        ## polarPoints[j,:] = aPoint + ((0.4 - i)/n)*derivativePoints[j,:]
        polarPoints[j,:] = polar(controlPoints, i, 0.4)
        j += 1
    end

    writeArrayForGnuplot(controlPoints,
                         "control-poly-polar.coordinates")

    ## bezierCurvePoints = drawCurve(controlPoints,
    ##     			  linspace(0,1,200))

    writeArrayForGnuplot(bezierPoints, "bezier-polar.coordinates")

    ## writeArrayForGnuplot(drawCurve([
    ##                                 [0.0 0 0];
    ##                                 differencePoints;
    ##                                 [0.0 0 0];
    ##                                 ],
    ##                                linspace(0,1,200)),
    ##     	         "differences-of-exercise-one.coordinates")

    writeArrayForGnuplot(polarPoints,
		         "derivative-polar.coordinates")

end

max_stars_repo_name: massimo-nocentini/cagd

max_issues_repo_name: massimo-nocentini/cagd

max_forks_repo_name: massimo-nocentini/cagd


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

include("api.jl")
include("common.jl")
include("generation.jl")
include("evaluation.jl")
include("spconstructs.jl")
# structures
include("deterministic_equivalent/generation.jl")
include("deterministic_equivalent/evaluation.jl")
include("deterministic_equivalent/optimization.jl")
include("deterministic_equivalent/spconstructs.jl")

include("stage_decomposition/generation.jl")
include("stage_decomposition/evaluation.jl")
include("stage_decomposition/optimization.jl")
include("stage_decomposition/spconstructs.jl")

include("scenario_decomposition/generation.jl")
include("scenario_decomposition/evaluation.jl")
include("scenario_decomposition/optimization.jl")
include("scenario_decomposition/spconstructs.jl")

include("util.jl")


text:
# Functions to perform stepper motor calibration
using LsqFit, JLD, StepperControl

"""
    create_xmodel(R, xorig)

## Description
Create a model for x displacements.
## Arguments
- R: bell crank length
- xorig: horizontal distance from the center of rotation to the origin
## Output
A function to convert steps, in x direction, into displacements with the following arguments:
- steps: number of steps (integer)
- p: parameters -> steps per revolution (to be find)
"""
function create_xmodel(R=80.0, xorig=-41.5)

    θ0 = asin(xorig/R)
    @. model(steps, p) = R*(sin(θ0 + 2π*steps/p[1]) - sin(θ0))

    return model
end

"""
    create_yxmodel(R, xorig)

## Description
Create a model for y displacements due to rotation of stepper motor responsible for x displacements.
## Arguments
- R: bell crank length
- xorig: horizontal distance from the center of rotation to the origin
## Output
A function to convert steps, in x direction, into displacements with the following arguments:
- steps: number of steps (integer)
- p: parameters -> steps per revolution (to be find)
"""
function create_yxmodel(R=80.0, xorig=-41.5)

    θ0 = asin(xorig/R)
    @. model(steps, p) = - R*(cos(θ0 + 2π*steps/p[1]) - cos(θ0))

    return model
end

"""
    create_ymodel(r)

## Description
Create a model for x displacements.
## Arguments
- r: gear radius
## Output
A function to convert steps, in y direction, into displacements with the following arguments:
- steps: number of steps (integer)
- p: parameters -> steps per revolution (to be find)
"""
function create_ymodel(r=12.5)

    @. model(steps, p) = 2π*steps*r/p[1]

    return model
end

"""
    iterative_calibr(dev::StepperSystem, steps::AbstractArray)

## Description
Iterative calibration process.
## Arguments
- dev: stepper motor device (see StepperControl)
- steps: array with the desired steps.
## Obs
The step values must be given relative to the origin (absolute displacements).
"""
function iterative_calibr(dev::StepperSystem, steps::AbstractArray)

    measurement = Float64[]
    println("Press enter to start calibration.\n")
    readline()

    for i in steps
        msg = stepper_move!(dev, i, relat=false)
        println(msg)
        println("Enter displacement for ", i, " steps: ")
        push!(measurement, parse(Float64, readline()))
    end

    println("\nEnd of calibration!")
    return measurement
end

"""
    calibr1DOF(sname::AbstractString="calibr.jld"; id, model, p0, steps::AbstractArray; measurement=missing)

## Description
Calibration function.
## Arguments
- sname: filename to be saved
- id: stepper ID (see StepperControl)
- model: model to be fit
- p0: initial paramenters (according to the model, elements must be of float type)
- steps: array with the desired steps (absolute displacement will be applied)
- measurement: measured data for each step. If is missing, then a iterative calibration is called
"""
function calibr1DOF(sname::AbstractString="calibr.jld"; id, model, p0, steps::AbstractArray, measurement=missing)

    n = length(steps)
    steps = Int.(round.(steps))

    if ismissing(measurement)
        r = stepper_open(1)
        stepper_config!(r, motorID=[id])
        println("Open the serial monitor and press enter to start.")
        readline()
        measurement = iterative_calibr(r, steps)
    end

    fit = curve_fit(model, steps, measurement, p0)
    save(sname, "fit", fit, "steps", steps, "measurement", measurement)

    return measurement, fit
end


text:
using LinearAlgebra

struct ZeroPolynomialInAlgebraicSetBridge{T, F <: MOI.AbstractVectorFunction,
                                          BT <: AbstractPolynomialBasis,
                                          DT <: AbstractSemialgebraicSet,
                                          MT <: AbstractMonomial,
                                          MVT <: AbstractVector{MT}} <: MOIB.AbstractBridge
    zero_constraint::MOI.ConstraintIndex{F, ZeroPolynomialSet{FullSpace, BT, MT, MVT}}
    domain::DT
    monomials::MVT
end

function ZeroPolynomialInAlgebraicSetBridge{T, F, BT, DT, MT, MVT}(model::MOI.ModelLike,
                                                                   f::MOI.AbstractVectorFunction,
                                                                   s::ZeroPolynomialSet{<:AbstractAlgebraicSet}) where {T, F <: MOI.AbstractVectorFunction,
                                                                                BT <: AbstractPolynomialBasis,
                                                                                DT <: AbstractSemialgebraicSet,
                                                                                MT <: AbstractMonomial,
                                                                                MVT <: AbstractVector{MT}}
    p = polynomial(collect(MOI.Utilities.eachscalar(f)), s.monomials)
    # As `*(::MOI.ScalarAffineFunction{T}, ::S)` is only defined if `S == T`, we
    # need to call `changecoefficienttype`. This is critical since `T` is
    # `Float64` when used with JuMP and the coefficient type is often `Int` with
    # `FixedVariablesSet`.
    # FIXME convert needed because the coefficient type of `r` is `Any` otherwise if `domain` is `AlgebraicSet`
    r = convert(typeof(p), rem(p, ideal(MultivariatePolynomials.changecoefficienttype(s.domain, T))))
    zero_constraint = MOI.add_constraint(model, MOIU.vectorize(coefficients(r)),
                                         ZeroPolynomialSet(FullSpace(), s.basis,
                                                           monomials(r)))
    return ZeroPolynomialInAlgebraicSetBridge{T, F, BT, DT, MT, MVT}(zero_constraint, s.domain, s.monomials)
end

function MOI.supports_constraint(::Type{ZeroPolynomialInAlgebraicSetBridge{T}},
                                 ::Type{<:MOI.AbstractVectorFunction},
                                 ::Type{<:ZeroPolynomialSet{<:AbstractAlgebraicSet}}) where T
    return true
end
function MOIB.added_constraint_types(::Type{<:ZeroPolynomialInAlgebraicSetBridge{T, F, BT, DT, MT, MVT}}) where {T, F, BT, DT, MT, MVT}
    return [(F, ZeroPolynomialSet{FullSpace, BT, MT, MVT})]
end
function MOIB.concrete_bridge_type(::Type{<:ZeroPolynomialInAlgebraicSetBridge{T}},
                                   F::Type{<:MOI.AbstractVectorFunction},
                                   ::Type{<:ZeroPolynomialSet{DT, BT, MT, MVT}}) where {T, BT, DT<:AbstractAlgebraicSet, MT, MVT}
    G = MOI.Utilities.promote_operation(-, T, F, F)
    return ZeroPolynomialInAlgebraicSetBridge{T, G, BT, DT, MT, MVT}
end

# Attributes, Bridge acting as an model
function MOI.get(::ZeroPolynomialInAlgebraicSetBridge{T, F, BT, DT, MT, MVT},
                 ::MOI.NumberOfConstraints{F, ZeroPolynomialSet{FullSpace, BT, MT, MVT}}) where {T, F, BT, DT, MT, MVT}
    return 1
end
function MOI.get(b::ZeroPolynomialInAlgebraicSetBridge{T, F, BT, DT, MT, MVT},
                 ::MOI.ListOfConstraintIndices{F, ZeroPolynomialSet{FullSpace, BT, MT, MVT}}) where {T, F, BT, DT, MT, MVT}
    return [b.zero_constraint]
end

# Indices
function MOI.delete(model::MOI.ModelLike, c::ZeroPolynomialInAlgebraicSetBridge)
    MOI.delete(model, c.zero_constraint)
end

# Attributes, Bridge acting as a constraint

# TODO ConstraintPrimal

# Let A be the linear map corresponding to A(p) = rem(p, ideal(set.domain))
# We want to compute A*(μ) (where A* is the conjugate of A); see
# slide 9 of https://www.youtube.com/watch?v=C8dHxJCUHYw.
# That is, for every monomial of the basis `mono`, we want to compute
# `⟨mono, A*(μ)⟩`. By definition of the conjugacy, we have `⟨mono, A*(μ)⟩ = ⟨A(mono), μ⟩`
# so we can compute compute it by `⟨rem(mono, ideal(set.domain)), μ⟩`
function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintDual,
                 bridge::ZeroPolynomialInAlgebraicSetBridge)
    dual = MOI.get(model, attr, bridge.zero_constraint)
    set = MOI.get(model, MOI.ConstraintSet(), bridge.zero_constraint)
    μ = measure(dual, set.monomials)
    I = ideal(bridge.domain)
    return [dot(rem(mono, I), μ) for mono in bridge.monomials]
end
function MOI.get(model::MOI.ModelLike, attr::MomentsAttribute,
                 bridge::ZeroPolynomialInAlgebraicSetBridge)
    return MOI.get(model, attr, bridge.zero_constraint)
end


text:
# TODO: comment everything
module AsciiMutator
  using Helper

  export mutate
  export Code

  type Code
    str::String
    #
    # Start string should contain 3 symbols: 1 unchanged, 2 empty and 3 unchanged.
    # String should be a correct julia code.
    #
    function Code() new("; ;") end
  end
  #
  # TODO: optimize this method as deep aspossible
  # Produces one add/change/del mutation on code and returns it's modified version.
  # Probabilities array is an array of three items: add,change,del values. e.g.:
  # [1,2,1] means that add and del mutations will be twice less then change one.
  # Character for adding and changing is in visible ASCII range - 32:126. see ASCII
  # table for details. Change mutation is the most frequent, so it should be first.
  # @param code Code we have to modify
  # @param prob Probability array (add,change,del). e.g.: [X,X,X]
  # @return {AsciiMutator.Code}
  #
  function mutate(code::Code, prob::Array{Int})
    pIndex = Helper.getProbIndex(prob)
    len    = length(code.str) - 1
    i      = rand(2:len)

    # 1 - add, 2 - change, 3 - del
    if pIndex === 2 code.str     = string(code.str[1:i-1], Char(rand(32:126)), code.str[i+1:end])
    elseif pIndex === 1 code.str = string(code.str[1:i-1], Char(rand(32:126)), code.str[i:end])
    elseif len > 2 code.str      = string(code.str[1:i-1], code.str[i+1:end]) end
    code
  end
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
#------------------------------------------------------------------------------

#------------------------------------------------------------------------------
####### Hashtable ######

# TODO: make immutable
#=
    stores hash of a monomial,
    corresponding divmask to speed up divisibility checks,
    index for position matrix (defaults to zero),
    and the todal degree
=#
mutable struct Hashvalue
    hash::UInt32

    #=

    =#
    divmask::UInt32

    idx::Int
    deg::Int
end

function copy_hashvalue(x::Hashvalue)
    Hashvalue(x.hash, x.divmask, x.idx, x.deg)
end

#=
    Hashtable implementing linear probing
    and designed to store and operate with monomials
=#
mutable struct MonomialHashtable
    exponents::Vector{ExponentVector}

    # maps exponent hash to its position in exponents array
    hashtable::Vector{Int}

    # stores hashes, division masks,
    # and other valuable info
    # for each hashtable enrty
    hashdata::Vector{Hashvalue}

    # values to hash exponents with, i.e
    # hash(e) = sum(hasher .* e)
    hasher::Vector{UInt32}

    #= Ring information =#
    # number of variables
    nvars::Int
    # raw length of exponent vector
    explen::Int
    # ring monomial ordering
    ord::Symbol

    #= Divisibility =#
    # divisor map to check divisibility faster
    divmap::Vector{UInt32}
    # variables selected for divmap
    divvars::Vector{Int}
    # count of divmap variables
    ndivvars::Int
    # bits per div variable
    ndivbits::Int

    size::Int
    # elements added
    load::Int
    #
    offset::Int
end

#------------------------------------------------------------------------------

function hashnextindex(h::UInt32, j::UInt32, mod::UInt32)
     (h + j) & mod + UInt32(1)
end

#------------------------------------------------------------------------------

# initialize and set fields for basis hashtable
function initialize_basis_hash_table(
        ring::PolyRing,
        rng::Random.AbstractRNG;
        initial_size::Int=2^16)

    exponents = Vector{Vector{UInt16}}(undef, initial_size)
    hashdata  = Vector{Hashvalue}(undef, initial_size)
    hashtable = zeros(Int, initial_size)

    nvars = ring.nvars
    explen = ring.explen
    ord = ring.ord

    # initialize hashing vector
    hasher = zeros(UInt32, explen)
    for i in 1:explen
        # we don't want hash vector components to be zero
        while hasher[i] == 0
            hasher[i] = rand(rng, UInt32)
        end
    end

    # exponents[1:load] cover all stored exponents
    # , also exponents[1] is zeroed by default
    load = 1
    size = initial_size

    # exponents array starts from index offset,
    # We store buffer array at index 1
    offset = 2

    # initialize fast divisibility params
    charbit    = 8 # TODO ??
    int32bits  = charbit * sizeof(Int32)
    int32bits != 32 && error("Strange story with ints")
    ndivbits   = div(int32bits, nvars)
    # division mask stores at least 1 bit
    # per each of first charbit*sizeof(Int32) variables
    ndivbits == 0 && (ndivbits += 1)
    ndivvars = nvars < int32bits ? nvars : int32bits
    divvars  = Vector{Int}(undef, ndivvars)
    divmap   = Vector{UInt32}(undef, ndivvars * ndivbits)
    # count only first ndivvars variables for divisibility checks
    for i in 1:ndivvars
        divvars[i] = i
    end

    # first stored exponent used as buffer lately
    exponents[1] = zeros(UInt16, explen)

    MonomialHashtable(
        exponents, hashtable, hashdata, hasher,
        nvars, explen, ord,
        divmap, divvars, ndivvars, ndivbits,
        size, load, offset)
end

function copy_hashtable(ht::MonomialHashtable)
    exps = Vector{Vector{UInt16}}(undef, ht.size)
    table = Vector{Int}(undef, ht.size)
    data = Vector{Hashvalue}(undef, ht.size)
    exps[1] = zeros(UInt16, ht.explen)

    @inbounds for i in 2:ht.load
        # TODO: ELIDE COPY
        exps[i] = copy(ht.exponents[i])
        table[i] = ht.hashtable[i]
        data[i] = copy_hashvalue(ht.hashdata[i])
    end

    MonomialHashtable(
        ht.exponents, ht.hashtable, ht.hashdata, ht.hasher,
        ht.nvars, ht.explen, ht.ord,
        ht.divmap, ht.divvars, ht.ndivvars, ht.ndivbits,
        ht.size, ht.load, ht.offset)
end

#------------------------------------------------------------------------------

# initialize hashtable either for `symbolic_preprocessing` or for `update` functions
# These are of the same purpose and structure as basis hashtable,
# but are more local oriented
function initialize_secondary_hash_table(basis_ht::MonomialHashtable)

    # 2^6 seems to be the best out of 2^5, 2^6, 2^7
    initial_size = 2^6

    exponents = Vector{Vector{UInt16}}(undef, initial_size)
    hashdata  = Vector{Hashvalue}(undef, initial_size)
    hashtable = zeros(Int, initial_size)

    # preserve ring info
    explen = basis_ht.explen
    nvars  = basis_ht.nvars
    ord    = basis_ht.ord

    # preserve division info
    divmap   = basis_ht.divmap
    divvars  = basis_ht.divvars
    ndivvars = basis_ht.ndivvars
    ndivbits = basis_ht.ndivbits

    # preserve hasher
    hasher = basis_ht.hasher

    load = 1
    size = initial_size
    offset = 2

    exponents[1] = zeros(UInt16, explen)

    MonomialHashtable(
        exponents, hashtable, hashdata, hasher,
        nvars, explen, ord,
        divmap, divvars, ndivvars, ndivbits,
        size, load, offset)
end

#------------------------------------------------------------------------------

# resizes (if needed) ht so that it can store `size` elements,
# and clears all previoud data
function reinitialize_hash_table!(ht::MonomialHashtable, size::Int)
    if size > ht.size
        while size > ht.size
            ht.size *= 2
        end
        # TODO
        resize!(ht.hashdata, ht.size)
        resize!(ht.exponents, ht.size)
    end
    ht.hashtable = zeros(Int, ht.size)
    ht.load = 1
end

# doubles the size of storage in `ht`,
# and rehashes all elements
function enlarge_hash_table!(ht::MonomialHashtable)
    ht.size *= 2
    resize!(ht.hashdata,  ht.size)
    resize!(ht.exponents, ht.size)

    ht.hashtable = zeros(Int, ht.size)

    mod = UInt32(ht.size - 1)
    for i in ht.offset:ht.load
        # hash for this elem is already computed
        he = ht.hashdata[i].hash
        hidx = he
        @inbounds for j in UInt32(1):UInt32(ht.size)
            hidx = hashnextindex(he, j, mod)
            ht.hashtable[hidx] != 0 && continue
            ht.hashtable[hidx] = i
            break
        end
    end
end

#------------------------------------------------------------------------------
####### Basis #######

# s-pair, a pair of polynomials
struct SPair
    # first generator as index from the basis array
    poly1::Int
    # second generator -//-
    poly2::Int
    # position of lcm(poly1, poly2) in hashtable
    lcm::Int
    # total degree of lcm
    deg::UInt
end

mutable struct Pairset
    pairs::Vector{SPair}
    # number of filled pairs,
    # Initially zero
    load::Int
end

function initialize_pairset(; initial_size=2^6) # TODO: why 64?
    pairs = Vector{SPair}(undef, initial_size)
    return Pairset(pairs, 0)
end

function Base.isempty(ps::Pairset)
    return ps.load == 0
end

function check_enlarge_pairset!(ps::Pairset, added::Int)
    sz = length(ps.pairs)
    # TODO: shrink pairset ?
    # @warn ps.load + added sz
    if ps.load + added >= sz
        # @error "resuze"  ps.load + added sz ps.load + added >= sz
        newsz = max(2 * sz, ps.load + added)
        resize!(ps.pairs, newsz)
    end
end

#------------------------------------------------------------------------------

mutable struct Basis{T}
    # vector of polynomials, each polynomial is a vector of monomials,
    # each monomial is represented with it's position in hashtable
    gens::Vector{Vector{Int}}
    # polynomial coefficients
    coeffs::Vector{Vector{T}}

    #= Keeping track of sizes   =#
    #=  ndone <= ntotal <= size =#
    # total allocated size,
    # size == length(gens) is always true
    size::Int
    # number of processed polynomials in `gens`
    # Iitially a zero
    ndone::Int
    # total number of polys filled in `gens`
    # (these will be handled during next update! call)
    # Iitially a zero
    ntotal::Int

    #= Keeping track of redundancy =#
    #= invariant =#
    #= length(lead) == length(nonred) == count(isred) == nlead =#
    # if element of the basis
    # is redundant
    isred::Vector{Int8}
    # positions of non-redundant elements in the basis
    nonred::Vector{Int}
    # division masks of leading monomials of
    # non redundant basis elements
    lead::Vector{UInt32}
    # number of filled elements in lead
    nlead::Int

    # characteristic of ground field
    ch::UInt64
end

function initialize_basis(ring::PolyRing, ngens::Int, ::Type{T}) where {T<:Coeff}
    #=
        always true
        length(gens) == length(coeffs) == length(isred) == size
    =#

    sz     = ngens * 2 # hmmm
    ndone  = 0
    ntotal = 0
    nlead  = 0

    gens   = Vector{Vector{Int}}(undef, sz)
    coeffs = Vector{Vector{T}}(undef, sz)
    isred  = zeros(Int8, sz)
    nonred = Vector{Int}(undef, sz)
    lead   = Vector{UInt32}(undef, sz)

    ch = ring.ch

    Basis(gens, coeffs, sz, ndone, ntotal, isred, nonred, lead, nlead, ch)
end

function initialize_basis(ring::PolyRing, hashedexps, coeffs::Vector{Vector{T}}) where {T<:Coeff}
    sz     = length(hashedexps) # hmmm
    ndone  = 0
    ntotal = 0
    nlead  = 0

    isred  = zeros(Int8, sz)
    nonred = Vector{Int}(undef, sz)
    lead   = Vector{UInt32}(undef, sz)

    ch = ring.ch

    Basis(hashedexps, coeffs, sz, ndone, ntotal, isred, nonred, lead, nlead, ch)
end

#------------------------------------------------------------------------------

function copy_basis_thorough(basis::Basis{T}) where {T}
    #  That cost a day of debugging ////
    gens   = Vector{Vector{Int}}(undef, basis.size)
    coeffs = Vector{Vector{T}}(undef, basis.size)
    @inbounds for i in 1:basis.ntotal
        gens[i] = Vector{Int}(undef, length(basis.gens[i]))
        coeffs[i] = Vector{T}(undef, length(basis.coeffs[i]))

        @inbounds for j in 1:length(basis.gens[i])
            gens[i][j] = basis.gens[i][j]
            coeffs[i][j] = basis.coeffs[i][j]
        end
    end
    isred  = copy(basis.isred)
    nonred = copy(basis.nonred)
    lead = copy(basis.lead)
    Basis(gens, coeffs, basis.size, basis.ndone,
            basis.ntotal, isred, nonred, lead,
            basis.nlead, basis.ch)
end

#------------------------------------------------------------------------------

function check_enlarge_basis!(basis::Basis{T}, added::Int) where {T}
    if basis.ndone + added >= basis.size
        basis.size = max(basis.size * 2, basis.ndone + added)
        resize!(basis.gens, basis.size)
        resize!(basis.coeffs, basis.size)
        resize!(basis.isred, basis.size)
        basis.isred[basis.ndone+1:end] .= 0
        resize!(basis.nonred, basis.size)
        resize!(basis.lead, basis.size)
    end
end

#------------------------------------------------------------------------------

# Normalize each element of the input basis
# by dividing it by leading coefficient
function normalize_basis!(basis::Basis{CoeffFF})
    cfs = basis.coeffs
    @inbounds for i in 1:basis.ntotal
        # mul = inv(cfs[i][1])
        # hack for now, TODODO
        if !isassigned(cfs, i)
            continue
        end
        ch = basis.ch
        mul = invmod(cfs[i][1], ch) % ch
        @inbounds for j in 2:length(cfs[i])
            # cfs[i][j] *= mul
            # TODO: faster division
            cfs[i][j] = cfs[i][j]*mul % ch
        end
        cfs[i][1] = one(cfs[i][1])
    end
    basis
end

# Normalize each element of the input basis
# by dividing it by leading coefficient
function normalize_basis!(basis::Basis{CoeffQQ})
    cfs = basis.coeffs
    @inbounds for i in 1:basis.ntotal
        # mul = inv(cfs[i][1])
        # hack for now, TODODO
        if !isassigned(cfs, i)
            continue
        end
        mul = inv(cfs[i][1])
        @inbounds for j in 2:length(cfs[i])
            cfs[i][j] *= mul
        end
        cfs[i][1] = one(cfs[i][1])
    end
    basis
end

#------------------------------------------------------------------------------
####### Matrix #######

mutable struct MacaulayMatrix{T<:Coeff}
    #=
        Matrix of the following structure

        | A  B |
        | C  D |

        A contains known pivots of reducing rows,
        and CD are rows to be reduced by AB

    =#

    # rows from upper, AB part of the matrix,
    # stored as vectors of corresponding exponents (already hashed)
    uprows::Vector{Vector{Int}}
    # rows from lower, CD part of the matrix,
    # stored as vectors of corresponding exponents (already hashed)
    lowrows::Vector{Vector{Int}}

    # maps column idx {1 ... ncols} to monomial hash {2 ... ht.load}
    # in some (?) hashtable
    col2hash::Vector{Int}

    # row coefficients
    # (some of the rows are stored in the basis,
    #  and some are stored here)
    coeffs::Vector{Vector{T}}

    #= sizes info =#
    # total number of allocated rows
    size::Int
    # number of pivots,
    # ie new basis elements discovered after matrix reduction
    npivots::Int
    # number of filled rows, nrows <= size
    nrows::Int
    # number of columns
    ncols::Int

    # number of upper rows (in AB section)
    nup::Int
    # number of lower rows (in CD section)
    nlow::Int
    # number of left cols  (in AC section)
    nleft::Int
    # number of right cols (in BD section)
    nright::Int

    # maps column idx {1 ... ncols} to index of coefficient array
    # First nleft indices point to coefficients from AB part of the matrix
    # (these coefficients are owned by the basis struct)
    # Last nright indices point to coefficients from CD part of the matrix
    # (these are owned by this object and stored in coeffs array)
    # Essentially each coefficient array from the first part represents a reducer row,
    # and each array from the second part stands for a reduced *nonzero* row
    # should be row to coef
    up2coef::Vector{Int}
    low2coef::Vector{Int}
end

function refresh_matrix!(matrix)
    matrix.npivots = 0
    matrix.nrows   = 0
    matrix.ncols   = 0
    matrix.nup     = 0
    matrix.nlow    = 0
    matrix.nleft   = 0
    matrix.nright  = 0
end

function initialize_matrix(ring::PolyRing, ::Type{T}) where {T<:Coeff}
    uprows   = Vector{Vector{Int}}(undef, 0)
    lowrows  = Vector{Vector{Int}}(undef, 0)
    col2hash = Vector{Int}(undef, 0)
    coeffs   = Vector{Vector{T}}(undef, 0)

    size    = 0
    npivots = 0
    nrows   = 0
    ncols   = 0

    nup    = 0
    nlow   = 0
    nleft  = 0
    nright = 0

    up2coef   = Vector{Int}(undef, 0)
    low2coef   = Vector{Int}(undef, 0)

    MacaulayMatrix(uprows, lowrows, col2hash, coeffs,
            size, npivots, nrows, ncols,
            nup, nlow, nleft, nright,
            up2coef, low2coef)
end

# TODO: change semantic
function reinitialize_matrix!(matrix::MacaulayMatrix{T}, npairs::Int) where {T}
    resize!(matrix.uprows, npairs*2)
    resize!(matrix.lowrows, npairs*2)
    resize!(matrix.up2coef, npairs*2)
    resize!(matrix.low2coef, npairs*2)
    matrix.size = 2*npairs
    matrix.ncols = 0
    matrix.nleft = 0
    matrix.nright = 0
    matrix.nup   = 0
    matrix.nlow  = 0
    matrix
end


text:
int main()
{
   int x = 10;
   int y, z, u, v, w;
   y = 13;
   w = x + y;
   z = y + 5;
   u = z - w;
   w = w % (-u);
   w++;
   y = z * w;
   w++;
   v = 6 * y;
   v = u + v;
   z = 6 * y + 5 * u - 4 * (w + u + v / 2 + z / 7);
   x = z * 16 + u * 10;
   printInt(x);
   printInt(y);
   printInt(z);
   printInt(u);
   printInt(v);
   printInt(w);
   return 0;
}

max_stars_repo_name: lukaszcz/javalette

max_issues_repo_name: lukaszcz/javalette

max_forks_repo_name: lukaszcz/javalette


text:
const NeuroMetaArray{T,N,A<:AbstractArray{T,N},D} = ImageMeta{T,N,A,<:NeuroMetadata{D}}

NeuroMetaArray(x::AbstractArray) = ImageMeta(x, NeuroMetadata())
NeuroMetaArray(x::AbstractArray, p::NeuroMetadata) = NeuroMetaArray(x, p)
NeuroMetaArray(x::AbstractArray, p::AbstractDict) = NeuroMetaArray(x, NeuroMetadata(p))
NeuroMetaArray(x; kwargs...) = NeuroMetaArray(x, NeuroMetadata(; kwargs...))

# TODO: add to ImageMetadata
ImageCore.HasProperties(::Type{T}) where {T<:ImageMeta} = HasProperties{true}()

@inline Base.getproperty(x::NeuroMetaArray, s::Symbol) = neuroproperty(x, s)
@inline Base.setproperty!(x::NeuroMetaArray, s::Symbol, val) = neuroproperty!(x, s, val)


text:
module ConcatOperator

VERSION < v"0.4.0-dev+6294" && error("""
    Cannot load ConcatOperator module.
    This version of Julia does not support parsing "++"
    """)

export ++

"""
Concatenation operator for strings, characters, and vectors.
"x ++ y ++ z ++ ..." is transformed by the parser to "++(x,y,z,...)".
"""
++(ss::Union{AbstractString,Char}...) = string(ss...)
++(x1::Vector, xs::Vector...) = vcat(x1, xs...)

end


text:
# overlay spots on image

function setpixel!(img, x, y, col)
    setpixel!(img, round(Int, x), round(Int, y), col)
    setpixel!(img, round(Int, x)+1, round(Int, y), col)
    setpixel!(img, round(Int, x), round(Int, y)+1, col)
    setpixel!(img, round(Int, x)-1, round(Int, y), col)
    setpixel!(img, round(Int, x), round(Int, y)-1, col)
end
function setpixel!(img, x::Int, y, col)
    if checkbounds(Bool, img, x, y)
        img[x, y] = col
    end
end
function overlay_spots(A, df, col; kwargs...)
    A1 = copy(A)
    overlay_spots!(A1, df, col; kwargs...)
    A1
end
function overlay_spots!(A::AbstractArray{<:Any, 3}, df::CSVFiles.CSVFile, col; kwargs...)
    overlay_spots!(A, collect(df), col; kwargs...)
end
function extract_frame(df::AbstractArray{T}, frame) where {NAMES, T <: NamedTuple{NAMES}}
    if :Centroid_X in NAMES
        [(x.Centroid_Y, x.Centroid_X, x.AreaInPixels) for x in df if x.Frame==frame && x.Centroid_X != 0]
    elseif :area in NAMES
        [(x.y, x.x, x.area) for x in df if x.frame==frame && x.x != 0]
    end
end
function overlay_spots!(A::AbstractArray{<:Any, 3}, df::AbstractArray{<:NamedTuple}, col; kwargs...)
    #(:x, :y, :area, :frame)
    for frame in 1:size(A, 3)
        d = extract_frame(df, frame)
        overlay_spots!((@view A[:,:,frame]), d, col; kwargs...)
    end
end
function overlay_spots!(A::AbstractArray{T, 2}, xyarea, col::T; d=-1.0) where T
    #d = get(kwargs, :d, Inf)
    for (i, (x, y, area)) in enumerate(xyarea)
        r = sqrt(area/pi) + 0.5
        circlen = pi*r
        if tooclose(xyarea, i, d)
            col1 = T(0.5, 0.5, 1.0)
        else
            col1 = col
        end
        for angle in range(0, pi, length=ceil(Int, circlen))
            dx = r*cos(angle)
            dy = r*sin(angle)
            setpixel!(A, x+dx, y+dy, col1)
            setpixel!(A, x-dx, y-dy, col1)
        end
    end
end

function tooclose(p, i::Int, R)
    x0, y0 = p[i]
    R2 = Float64(R^2)
    for i1 in eachindex(p)
        if i != i1
            x1, y1 = p[i1]
            d2 = (x1-x0)^2 + (y1-y0)^2
            if d2 < R2
                return true
            end
        end
    end
    false
end


text:
# # Testing different approximations for the rate-of-change law

# Ultimately we want to evolve, in time, an equation of the form $u_t = F(u)$. For the 2D NSE, the most computationally demanding term is the bilinear one: $B(\omega) = \left(\partial_x^2 - \partial_y^2\right)(uv) + \partial_{xy}\left(v^2 - u^2\right)$, where $\mathbf{u} = (u(x, y), v(x, y))$ is the velocity field, which, in the two-dimensional periodic case, can easily be obtained from the vorticity $\omega = \boldsymbol{\nabla} \times \mathbf{u}$.

# Here, we test different ways of approximating $B(u)$

# Here are the packages we are gonna need.

using FFTW
using Plots
using LinearAlgebra: mul!
using Test
using BenchmarkTools
using Random

@info "Threads: $(FFTW.nthreads())"

# ## Operators

# We first define a method to build the various operators acting in spectral space:

function get_operators(N, κ₀)
    ## Differentiation in spectral space
    Dx_hat = im * κ₀ * [
        ifelse(k1 ≤ div(N, 2) + 1, k1 - 1, k1 - 1 - N) for k2 in 1:div(N, 2)+1, k1 in 1:N
    ]
    Dy_hat = im * κ₀ * [k2 - 1 for k2 in 1:div(N, 2)+1, k1 in 1:N]

    Delta_hat = - κ₀^2 * [
        ifelse(k1 ≤ div(N, 2) + 1, (k1 - 1)^2 + (k2 - 1)^2, (k1 - 1 - N)^2 + (k2 - 1)^2)
        for k2 in 1:div(N, 2)+1, k1 in 1:N
    ]

    ## Recovering of the velocity field from the vorticity
    Hu_hat = - Dy_hat ./ Delta_hat
    Hu_hat[1, 1] = 0.0
    Hv_hat = Dx_hat ./ Delta_hat
    Hv_hat[1, 1] = 0.0

    ## For the Basdevant formulation
    DxsqDysq_hat = Dx_hat.^2 .- Dy_hat.^2
    Dxy_hat = Dx_hat .* Dy_hat

    return Dx_hat, Dy_hat, Delta_hat, Hu_hat, Hv_hat, DxsqDysq_hat, Dxy_hat
end

# Methods to generate a scalar field from a certain list of excitable wavenumbers.

function field_from_spectrum(L, N, modes::Matrix{<:Integer}, amps::Matrix{<:Real})
    κ₀ = 2π/L
    x = y = (L/N):(L/N):L
    field = sum(
        [
            2κ₀^2 * (k[1]^2 + k[2]^2) * (
                a[1] * cos.(κ₀ * (k[1] * one.(y) * x' + k[2] * y * one.(x)'))
                - a[2] * sin.(κ₀ * (k[1] * one.(y) * x' + k[2] * y * one.(x)'))
            )
            for (k, a) in zip(eachrow(modes), eachrow(amps))
        ]
    )
    return field
end

function field_from_spectrum(rng::AbstractRNG, L, N, num_modes::Int)

    modes = rand(rng, 1:div(N,10), num_modes, 2)
    amps = rand(rng, num_modes, 2)

    field = field_from_spectrum(L, N, modes, amps)

    return field
end

field_from_spectrum(L, N, num_modes::Int) = field_from_spectrum(Xoshiro(), L, N, num_modes)

# Evolution steps

function step_naive!(vort_hat, dt, params)
    operators, vars = params
    Dx_hat, Dy_hat, Delta_hat, Hu_hat, Hv_hat,
        DxsqDysq_hat, Dxy_hat, Exp_nu_dt_Delta_hat = operators
    g_hat, N, Nsub = vars
    u_hat = Hu_hat .* vort_hat
    v_hat = Hv_hat .* vort_hat
    u = irfft(u_hat, N)
    v = irfft(v_hat, N)
    vort = irfft(vort_hat, N)
    wu_hat = rfft(vort .* u)
    wv_hat = rfft(vort .* v)
    vort_hat .= Exp_nu_dt_Delta_hat .* (
        vort_hat .+ dt * (
            g_hat .- Dx_hat .* wu_hat .- Dy_hat .* wv_hat
        )
    )
    ## dealiasing
    vort_hat[div(Nsub,2) + 1:end, :] .= 0.0im
    vort_hat[:, div(Nsub,2) + 1:div(N,2) + div(Nsub,2)] .= 0.0im
    return vort_hat
end

function step_Basdevant!(vort_hat, dt, params)
    operators, vars = params
    Dx_hat, Dy_hat, Delta_hat, Hu_hat, Hv_hat,
        DxsqDysq_hat, Dxy_hat, Exp_nu_dt_Delta_hat = operators
    g_hat, N, Nsub = vars
    u_hat = Hu_hat .* vort_hat
    v_hat = Hv_hat .* vort_hat
    u = irfft(u_hat, N)
    v = irfft(v_hat, N)
    uv_hat = rfft(u .* v)
    v2u2_hat = rfft(v.^2 - u.^2)
    vort_hat .= Exp_nu_dt_Delta_hat .* (
        vort_hat .+ dt * (
            g_hat .- DxsqDysq_hat .* uv_hat .- Dxy_hat .* v2u2_hat
        )
    )
    ## dealiasing
    vort_hat[div(Nsub,2) + 1:end, :] .= 0.0im
    vort_hat[:, div(Nsub,2) + 1:div(N,2) + div(Nsub,2)] .= 0.0im
    return vort_hat
end

function step_Basdevant_plan!(vort_hat, dt, params)
    operators, vars, auxs, plans = params

    Dx_hat, Dy_hat, Delta_hat, Hu_hat, Hv_hat,
        DxsqDysq_hat, Dxy_hat, Exp_nu_dt_Delta_hat = operators
    g_hat, N, Nsub = vars
    u_hat, v_hat, u, v, uv, v2u2, uv_hat, v2u2_hat = auxs
    plan, plan_inv = plans

    u_hat .= Hu_hat .* vort_hat
    v_hat .= Hv_hat .* vort_hat

    mul!(u, plan_inv, u_hat)
    mul!(v, plan_inv, v_hat)

    uv .= u .* v
    v2u2 .= v.^2 .- u.^2
    mul!(uv_hat, plan, uv)
    mul!(v2u2_hat, plan, v2u2)

    @. vort_hat = Exp_nu_dt_Delta_hat * (
        vort_hat + dt * (
            g_hat - DxsqDysq_hat * uv_hat - Dxy_hat * v2u2_hat
        )
    )

#=     vort_hat .= Exp_nu_dt_Delta_hat .* (
        vort_hat .+ dt .* (
            g_hat .- DxsqDysq_hat .* uv_hat .- Dxy_hat .* v2u2_hat
        )
    ) =#
    ## dealiasing
    vort_hat[div(Nsub,2) + 1:end, :] .= 0.0im
    vort_hat[:, div(Nsub,2) + 1:div(N,2) + div(Nsub,2)] .= 0.0im
    return vort_hat
end

function step_Basdevant_plan_loop!(vort_hat, dt, params)
    operators, vars, auxs, plans = params

    Dx_hat, Dy_hat, Delta_hat, Hu_hat, Hv_hat,
        DxsqDysq_hat, Dxy_hat, Exp_nu_dt_Delta_hat = operators
    g_hat, N, Nsub = vars
    u_hat, v_hat, u, v, uv, v2u2, uv_hat, v2u2_hat = auxs
    plan, plan_inv = plans

    @inbounds for i in eachindex(u_hat)
        u_hat[i] = Hu_hat[i] * vort_hat[i]
        v_hat[i] = Hv_hat[i] .* vort_hat[i]
    end

    mul!(u, plan_inv, u_hat)
    mul!(v, plan_inv, v_hat)

    @inbounds for i in eachindex(u)
        uv[i] = u[i] * v[i]
        v2u2[i] = v[i]^2 - u[i]^2
    end

    mul!(uv_hat, plan, uv)
    mul!(v2u2_hat, plan, v2u2)

    for i in eachindex(vort_hat)
        vort_hat[i] = Exp_nu_dt_Delta_hat[i] * (
            vort_hat[i] + dt * (
                g_hat[i] - DxsqDysq_hat[i] * uv_hat[i] - Dxy_hat[i] * v2u2_hat[i]
            )
        )
    end
    ## dealiasing
    @inbounds for j in 1:N, i in div(Nsub,2) + 1:div(N, 2) + 1
        vort_hat[i, j] = 0.0im
    end
    @inbounds for j in div(Nsub,2) + 1:div(N,2) + div(Nsub,2), i in 1:div(N, 2) + 1
        vort_hat[i, j] = 0.0im
    end
    return vort_hat
end

function step_Basdevant_plan_doubleloop!(vort_hat, dt, params)
    operators, vars, auxs, plans = params

    Dx_hat, Dy_hat, Delta_hat, Hu_hat, Hv_hat,
        DxsqDysq_hat, Dxy_hat, Exp_nu_dt_Delta_hat = operators
    g_hat, N, Nsub = vars
    u_hat, v_hat, u, v, uv, v2u2, uv_hat, v2u2_hat = auxs
    plan, plan_inv = plans

    M, N = size(vort_hat)

    @inbounds for j in 1:N, i in 1:M
        u_hat[i, j] = Hu_hat[i, j] * vort_hat[i, j]
        v_hat[i, j] = Hv_hat[i, j] .* vort_hat[i, j]
    end

    mul!(u, plan_inv, u_hat)
    mul!(v, plan_inv, v_hat)

    @inbounds for j in 1:N, i in 1:N
        uv[i, j] = u[i, j] * v[i, j]
        v2u2[i, j] = v[i, j]^2 - u[i, j]^2
    end

    mul!(uv_hat, plan, uv)
    mul!(v2u2_hat, plan, v2u2)

    for j in 1:N, i in 1:M
        vort_hat[i, j] = Exp_nu_dt_Delta_hat[i, j] * (
            vort_hat[i, j] + dt * (
                g_hat[i, j] - DxsqDysq_hat[i, j] * uv_hat[i, j] - Dxy_hat[i, j] * v2u2_hat[i, j]
            )
        )
    end
    ## dealiasing
    @inbounds for j in 1:N, i in div(Nsub,2) + 1:M
        vort_hat[i, j] = 0.0im
    end
    @inbounds for j in div(Nsub,2) + 1:div(N,2) + div(Nsub,2), i in 1:M
        vort_hat[i, j] = 0.0im
    end
    return vort_hat
end

# ## The spatial domain

L = 2π
κ₀ = 2π/L
N = 128
Nsub = 84
x = y = (L/N):(L/N):L

ν = 1.0

dt = 1.0e-5 # 1.0e-3 a 1.0e-5
t_final = 1.0

num_steps = Int(round(t_final / dt))

Dx_hat, Dy_hat, Delta_hat, Hu_hat, Hv_hat, DxsqDysq_hat, Dxy_hat = get_operators(N, κ₀)
Exp_nu_dt_Delta_hat = exp.(ν * dt * Delta_hat)

vort_init = field_from_spectrum(L, N, 4)
vort_init_hat = rfft(vort_init)

g_steady = field_from_spectrum(L, N, 1)
g_steady_hat = rfft(g_steady)

vort_steady_hat = - g_steady_hat ./ Delta_hat
vort_steady_hat[1, 1] = 0.0im
vort_steady = irfft(vort_steady_hat, N)

vort_hat = copy(vort_init_hat)
vort = irfft(vort_hat, N)

u_hat = zero(vort_hat)
v_hat = zero(vort_hat)
u = zeros(N, N)
v = zeros(N, N)
uv = zeros(N, N)
v2u2 = zeros(N, N)
uv_hat = zero(vort_hat)
v2u2_hat = zero(vort_hat)
plan = plan_rfft(vort, flags=FFTW.MEASURE)
plan_inv = plan_irfft(vort_hat, N, flags=FFTW.MEASURE)

operators =
    Dx_hat, Dy_hat, Delta_hat, Hu_hat, Hv_hat, DxsqDysq_hat, Dxy_hat, Exp_nu_dt_Delta_hat
vars = g_steady_hat, N, Nsub
auxs = u_hat, v_hat, u, v, uv, v2u2, uv_hat, v2u2_hat
plans = plan, plan_inv

params = (
    operators,
    vars,
    auxs,
    plans
)

vort_hat = copy(vort_init_hat)
vort = irfft(vort_hat, N)

println("Enstrophy convergence:")
println(sum(abs2, vort_hat - vort_steady_hat) * (1/N)^4)
for n in 1:num_steps
    ## step_naive!(vort_hat, dt, (operators, vars))
    ## step_Basdevant!(vort_hat, dt, (operators, vars))
    ## step_Basdevant_plan!(vort_hat, dt, (operators, vars, auxs, plans))
    step_Basdevant_plan_doubleloop!(vort_hat, dt, (operators, vars, auxs, plans))
    if rem(n, 1000) == 0
        println(sum(abs2, vort_hat - vort_steady_hat) * (1/N)^4)
    end
end

vort = irfft(vort_hat, N)
display(heatmap(x, y, vort, xlabel="x", ylabel="y", title="vorticity", titlefont=12))

display(surface(x, y, vort - vort_steady, xlabel="x", ylabel="y", zlabel="error", title="difference `vort_final .- vort_steady`", titlefont=12))

@info "step_naive!"
@btime step_naive!(vh, τ, p) setup = (vh = copy(vort_hat); τ = $dt; p = $(operators, vars));
@info "step_Basdevant!"
@btime step_Basdevant!(vh, τ, p) setup = (vh = copy(vort_hat); τ = $dt; p = $(operators, vars));
@info "step_Basdevant_plan!"
@btime step_Basdevant_plan!(vh, τ, p) setup = (vh = copy(vort_hat); τ = $dt; p = $(operators, vars, auxs, plans));
@info "step_Basdevant_plan_loop!"
@btime step_Basdevant_plan_loop!(vh, τ, p) setup = (vh = copy(vort_hat); τ = $dt; p = $(operators, vars, auxs, plans));
@info "step_Basdevant_plan_doubleloop!"
@btime step_Basdevant_plan_doubleloop!(vh, τ, p) setup = (vh = copy(vort_hat); τ = $dt; p = $(operators, vars, auxs, plans));

#=
[ Info: step_naive!
  261.500 μs (211 allocations: 2.16 MiB)
[ Info: step_Basdevant!
  222.791 μs (174 allocations: 2.03 MiB)
[ Info: step_Basdevant_plan!
  136.834 μs (0 allocations: 0 bytes)
[ Info: step_Basdevant_plan_loop!
  131.042 μs (0 allocations: 0 bytes)
[ Info: step_Basdevant_plan_doubleloop!
  137.000 μs (0 allocations: 0 bytes)
=#