text:
"""
    AbstractBasis{F}

Abstract type that represents a basis on a manifold or a subset of it.

The type parameter `F` denotes the [`AbstractNumbers`](@ref) that will be used as scalars.
"""
abstract type AbstractBasis{F} end

"""
    number_system(::AbstractBasis)

The number system used as scalars in the given basis.
"""
number_system(::AbstractBasis{F}) where {F} = F

"""
    AbstractOrthonormalBasis{F}

Abstract type that represents an orthonormal basis on a manifold or a subset of it.

The type parameter `F` denotes the [`AbstractNumbers`](@ref) that will be used as scalars.
"""
abstract type AbstractOrthonormalBasis{F} <: AbstractBasis{F} end

"""
    AbstractPrecomputedOrthonormalBasis{F}

Abstract type that represents an orthonormal basis of the tangent space at a point
on a manifold. Tangent vectors can be obtained using function [`get_vectors`](@ref).

The vectors are not always fully precomputed because a partially precomputed
basis may be enough for implementing [`get_vector`](@ref) and [`get_coordinates`](@ref).

The type parameter `F` denotes the [`AbstractNumbers`](@ref) that will be used as scalars.
"""
abstract type AbstractPrecomputedOrthonormalBasis{F} <: AbstractOrthonormalBasis{F} end

"""
    ArbitraryOrthonormalBasis(F::AbstractNumbers = ℝ)

An arbitrary orthonormal basis on a manifold. This will usually
be the fastest orthonormal basis available for a manifold.

The type parameter `F` denotes the [`AbstractNumbers`](@ref) that will be used as scalars.
"""
struct ArbitraryOrthonormalBasis{F} <: AbstractOrthonormalBasis{F} end

ArbitraryOrthonormalBasis(F::AbstractNumbers = ℝ) = ArbitraryOrthonormalBasis{F}()

"""
    ProjectedOrthonormalBasis(method::Symbol, F::AbstractNumbers = ℝ)

An orthonormal basis that comes from orthonormalization of basis vectors
of the ambient space projected onto the subspace representing the tangent space
at a given point.

The type parameter `F` denotes the [`AbstractNumbers`](@ref) that will be used as scalars.

Available methods:
  - `:gram_schmidt` uses a modified Gram-Schmidt orthonormalization.
  - `:svd` uses SVD decomposition to orthogonalize projected vectors.
    The SVD-based method should be more numerically stable at the cost of
    an additional assumption (local metric tensor at a point where the
    basis is calculated has to be diagonal).
"""
struct ProjectedOrthonormalBasis{Method,F} <: AbstractOrthonormalBasis{F} end

function ProjectedOrthonormalBasis(method::Symbol, F::AbstractNumbers = ℝ)
    return ProjectedOrthonormalBasis{method,F}()
end

@doc doc"""
    DiagonalizingOrthonormalBasis(v, F::AbstractNumbers = ℝ)

An orthonormal basis `Ξ` as a vector of tangent vectors (of length determined by
[`manifold_dimension`](@ref)) in the tangent space that diagonalizes the curvature
tensor $R(u,v)w$ and where the direction `v` has curvature `0`.

The type parameter `F` denotes the [`AbstractNumbers`](@ref) that will be used as scalars.
"""
struct DiagonalizingOrthonormalBasis{TV,F} <: AbstractOrthonormalBasis{F}
    v::TV
end

function DiagonalizingOrthonormalBasis(v, F::AbstractNumbers = ℝ)
    return DiagonalizingOrthonormalBasis{typeof(v),F}(v)
end

const ArbitraryOrDiagonalizingBasis =
    Union{ArbitraryOrthonormalBasis,DiagonalizingOrthonormalBasis}

"""
    PrecomputedOrthonormalBasis(vectors::AbstractVector, F::AbstractNumbers = ℝ)

A precomputed orthonormal basis at a point on a manifold.

The type parameter `F` denotes the [`AbstractNumbers`](@ref) that will be used as scalars.
"""
struct PrecomputedOrthonormalBasis{TV<:AbstractVector,F} <:
       AbstractPrecomputedOrthonormalBasis{F}
    vectors::TV
end

function PrecomputedOrthonormalBasis(vectors::AbstractVector, F::AbstractNumbers = ℝ)
    return PrecomputedOrthonormalBasis{typeof(vectors),F}(vectors)
end

@doc doc"""
    DiagonalizingOrthonormalBasis(vectors, kappas, F::AbstractNumbers = ℝ)

A precomputed orthonormal basis `Ξ` as a vector of tangent vectors (of length determined
by [`manifold_dimension`](@ref)) in the tangent space that diagonalizes the curvature
tensor $R(u,v)w$ with eigenvalues `kappas` and where the direction `v` has curvature `0`.

The type parameter `F` denotes the [`AbstractNumbers`](@ref) that will be used as scalars.
"""
struct PrecomputedDiagonalizingOrthonormalBasis{TV<:AbstractVector,TK<:AbstractVector,F} <:
       AbstractPrecomputedOrthonormalBasis{F}
    vectors::TV
    kappas::TK
end

function PrecomputedDiagonalizingOrthonormalBasis(
    vectors::AbstractVector,
    kappas::AbstractVector,
    F::AbstractNumbers = ℝ,
)
    return PrecomputedDiagonalizingOrthonormalBasis{typeof(vectors),typeof(kappas),F}(
        vectors,
        kappas,
    )
end

"""
    get_coordinates(M::Manifold, x, v, B::AbstractBasis)

Compute a one-dimensional vector of coefficients of the tangent vector `v`
at point denoted by `x` on manifold `M` in basis `B`.

Depending on the basis, `x` may not directly represent a point on the manifold.
For example if a basis transported along a curve is used, `x` may be the coordinate
along the curve.

See also: [`get_vector`](@ref), [`get_basis`](@ref)
"""
function get_coordinates(M::Manifold, x, v, B::AbstractBasis)
    error("get_coordinates not implemented for manifold of type $(typeof(M)) a point of type $(typeof(x)), tangent vector of type $(typeof(v)) and basis of type $(typeof(B)).")
end
function get_coordinates(M::Manifold, x, v, B::AbstractPrecomputedOrthonormalBasis{ℝ})
    return map(vb -> real(inner(M, x, v, vb)), get_vectors(M, x, B))
end
function get_coordinates(M::Manifold, x, v, B::AbstractPrecomputedOrthonormalBasis)
    return map(vb -> inner(M, x, v, vb), get_vectors(M, x, B))
end

"""
    get_vector(M::Manifold, x, v, B::AbstractBasis)

Convert a one-dimensional vector of coefficients in a basis `B` of
the tangent space at `x` on manifold `M` to a tangent vector `v` at `x`.

Depending on the basis, `x` may not directly represent a point on the manifold.
For example if a basis transported along a curve is used, `x` may be the coordinate
along the curve.

See also: [`get_coordinates`](@ref), [`get_basis`](@ref)
"""
function get_vector(M::Manifold, x, v, B::AbstractBasis)
    error("get_vector not implemented for manifold of type $(typeof(M)) a point of type $(typeof(x)), tangent vector of type $(typeof(v)) and basis of type $(typeof(B)).")
end
function get_vector(M::Manifold, x, v, B::AbstractPrecomputedOrthonormalBasis)
    # quite convoluted but:
    #  1) preserves the correct `eltype`
    #  2) guarantees a reasonable array type `vout`
    #     (for example scalar * `SizedArray` is an `SArray`)
    bvectors = get_vectors(M, x, B)
    if isa(bvectors[1], ProductRepr)
        vt = v[1] * bvectors[1]
        vout = allocate(bvectors[1], eltype(vt))
        copyto!(vout, vt)
        for i = 2:length(v)
            vout += v[i] * bvectors[i]
        end
        return vout
    else
        vt = v[1] .* bvectors[1]
        vout = allocate(bvectors[1], eltype(vt))
        copyto!(vout, vt)
        for i = 2:length(v)
            vout .+= v[i] .* bvectors[i]
        end
        return vout
    end
end

"""
    get_basis(M::Manifold, x, B::AbstractBasis) -> AbstractBasis

Compute the basis vectors of the tangent space at a point on manifold `M`
represented by `x`.

Returned object derives from [`AbstractBasis`](@ref) and may have a field `.vectors`
that stores tangent vectors or it may store them implicitly, in which case
the function [`get_vectors`](@ref) needs to be used to retrieve the basis vectors.

See also: [`get_coordinates`](@ref), [`get_vector`](@ref)
"""
function get_basis(M::Manifold, x, B::AbstractBasis)
    error("get_basis not implemented for manifold of type $(typeof(M)) a point of type $(typeof(x)) and basis of type $(typeof(B)).")
end
"""
    get_basis(M::Manifold, x, B::ArbitraryOrthonormalBasis)

Compute the basis vectors of an [`ArbitraryOrthonormalBasis`](@ref).
"""
function get_basis(M::Manifold, x, B::ArbitraryOrthonormalBasis)
    dim = manifold_dimension(M)
    return PrecomputedOrthonormalBasis([
        get_vector(M, x, [ifelse(i == j, 1, 0) for j = 1:dim], B) for i = 1:dim
    ])
end
get_basis(M::Manifold, x, B::AbstractPrecomputedOrthonormalBasis) = B
function get_basis(M::ArrayManifold, x, B::AbstractPrecomputedOrthonormalBasis{ℝ})
    bvectors = get_vectors(M, x, B)
    N = length(bvectors)
    M_dim = manifold_dimension(M)
    if N != M_dim
        throw(ArgumentError("Incorrect number of basis vectors; expected: $M_dim, given: $N"))
    end
    for i = 1:N
        vi_norm = norm(M, x, bvectors[i])
        if !isapprox(vi_norm, 1)
            throw(ArgumentError("vector number $i is not normalized (norm = $vi_norm)"))
        end
        for j = i+1:N
            dot_val = real(inner(M, x, bvectors[i], bvectors[j]))
            if !isapprox(dot_val, 0; atol = eps(eltype(x)))
                throw(ArgumentError("vectors number $i and $j are not orthonormal (inner product = $dot_val)"))
            end
        end
    end
    return B
end

function get_coordinates(M::ArrayManifold, x, v, B::AbstractBasis; kwargs...)
    is_tangent_vector(M, x, v, true; kwargs...)
    return get_coordinates(M.manifold, x, v, B)
end

function get_vector(M::ArrayManifold, x, v, B::AbstractBasis; kwargs...)
    is_manifold_point(M, x, true; kwargs...)
    size(v) == (manifold_dimension(M),) || error("Incorrect size of vector v")
    return get_vector(M.manifold, x, v, B)
end

function _euclidean_basis_vector(x, i)
    y = zero(x)
    y[i] = 1
    return y
end

function get_basis(M::Manifold, x, B::ProjectedOrthonormalBasis{:svd,ℝ})
    S = representation_size(M)
    PS = prod(S)
    dim = manifold_dimension(M)
    # projection
    # TODO: find a better way to obtain a basis of the ambient space
    vs = [
        convert(Vector, reshape(project_tangent(M, x, _euclidean_basis_vector(x, i)), PS))
        for i in eachindex(x)
    ]
    O = reduce(hcat, vs)
    # orthogonalization
    # TODO: try using rank-revealing QR here
    decomp = svd(O)
    rotated = Diagonal(decomp.S) * decomp.Vt
    vecs = [collect(reshape(rotated[i, :], S)) for i = 1:dim]
    # normalization
    for i = 1:dim
        i_norm = norm(M, x, vecs[i])
        vecs[i] /= i_norm
    end
    return PrecomputedOrthonormalBasis(vecs)
end

"""
    get_vectors(M::Manifold, x, B::AbstractBasis)

Get the basis vectors of basis `B` of the tangent space at point `x`.
"""
function get_vectors(M::Manifold, x, B::AbstractBasis)
    error("get_vectors not implemented for manifold of type $(typeof(M)) a point of type $(typeof(x)) and basis of type $(typeof(B)).")
end

get_vectors(::Manifold, x, B::PrecomputedOrthonormalBasis) = B.vectors
get_vectors(::Manifold, x, B::PrecomputedDiagonalizingOrthonormalBasis) = B.vectors

# related to DefaultManifold; to be moved to ManifoldsBase.jl in the future
function get_coordinates(M::DefaultManifold, x, v, ::ArbitraryOrthonormalBasis)
    return reshape(v, manifold_dimension(M))
end

function get_vector(M::DefaultManifold, x, v, ::ArbitraryOrthonormalBasis)
    return reshape(v, representation_size(M))
end

function get_basis(M::DefaultManifold, x, ::ArbitraryOrthonormalBasis)
    return PrecomputedOrthonormalBasis([
        _euclidean_basis_vector(x, i) for i in eachindex(x)
    ])
end

function get_basis(M::Manifold, x, B::ProjectedOrthonormalBasis{:gram_schmidt,ℝ}; kwargs...)
    E = [_euclidean_basis_vector(x, i) for i in eachindex(x)]
    N = length(E)
    Ξ = empty(E)
    dim = manifold_dimension(M)
    N < dim && @warn "Input only has $(N) vectors, but manifold dimension is $(dim)."
    K = 0
    @inbounds for n = 1:N
        Ξₙ = project_tangent(M, x, E[n])
        for k = 1:K
            Ξₙ .-= real(inner(M, x, Ξ[k], Ξₙ)) .* Ξ[k]
        end
        nrmΞₙ = norm(M, x, Ξₙ)
        if nrmΞₙ == 0
            @warn "Input vector $(n) has length 0."
            @goto skip
        end
        Ξₙ ./= nrmΞₙ
        for k = 1:K
            if !isapprox(real(inner(M, x, Ξ[k], Ξₙ)), 0; kwargs...)
                @warn "Input vector $(n) is not linearly independent of output basis vector $(k)."
                @goto skip
            end
        end
        push!(Ξ, Ξₙ)
        K += 1
        K * real_dimension(number_system(B)) == dim && return PrecomputedOrthonormalBasis(Ξ)
        @label skip
    end
    @warn "get_basis with bases $(typeof(B)) only found $(K) orthonormal basis vectors, but manifold dimension is $(dim)."
    return PrecomputedOrthonormalBasis(Ξ)
end


text:
export FaultedObservable, faulted

import Base: ==
import Base: show

"""
    FaultedObservable{D, H}(err, scheduler::H)

Observable that emits no items to the Actor and just sends an error notification on subscription.

# Constructor arguments
- `err`: error to emit on subscription
- `scheduler`: scheduler-like object

See also: [`faulted`](@ref)
"""
struct FaultedObservable{D, H} <: ScheduledSubscribable{D}
    err
    scheduler :: H
end

getscheduler(observable::FaultedObservable) = observable.scheduler

function on_subscribe!(observable::FaultedObservable, actor, scheduler)
    error!(actor, observable.err, scheduler)
    return voidTeardown
end

"""
    faulted(err; scheduler::H = AsapScheduler())            where { H <: AbstractScheduler }
    faulted(::Type{T}, err; scheduler::H = AsapScheduler()) where { T, H <: AbstractScheduler }

Creation operator for the `FaultedObservable` that emits no items to the Actor and immediately sends an error notification.

# Arguments
- `err`: the particular Error to pass to the error notification.
- `T`: type of output data source, optional, `Any` by default
- `scheduler`: optional, scheduler-like object

# Examples

```jldoctest
using Rocket

source = faulted("Error!")
subscribe!(source, logger())
;

# output

[LogActor] Error: Error!

```

See also: [`FaultedObservable`](@ref), [`subscribe!`](@ref)
"""
function faulted end

faulted(::Type{T}; scheduler::H = AsapScheduler())      where { T, H <: AbstractScheduler } = error("Missing error value in faulted() constructor.")
faulted(err; scheduler::H = AsapScheduler())            where { H <: AbstractScheduler }    = FaultedObservable{Any, H}(err, scheduler)
faulted(::Type{T}, err; scheduler::H = AsapScheduler()) where { T, H <: AbstractScheduler } = FaultedObservable{T, H}(err, scheduler)

Base.:(==)(e1::FaultedObservable{D, H},  e2::FaultedObservable{D, H}) where { D, H } = e1.err == e2.err
Base.:(==)(e1::FaultedObservable,     e2::FaultedObservable)                         = false

Base.show(io::IO, ::FaultedObservable{D, H}) where { D, H } = print(io, "FaultedObservable($D, $H)")

@deprecate throwError(T)    faulted(T)
@deprecate throwError(T, e) faulted(T, e)

ext: jl


text:
using plan_graph
using Test

@testset "plan_graph.jl" begin
    # Write your tests here.
end
#f = Compound(Symbol("on-table"), [Const(:a)])


text:
"""
Two different temporal plots, one with memory and one with some percolation rate.
While varying the percolation rate, we compare the ratio of the depths of the
graphs
"""

using QuNet
using LightGraphs, SimpleWeightedGraphs, GraphPlot, MetaGraphs, Plots, Colors, Statistics
using LaTeXStrings
using JLD
using Parameters

datafile = "data/bandwidth_wtih_memory"

# Params
num_trials = 10::Int64
perc_range = (0.05, 0.05, 0.95)::Tuple{Float64, Float64, Float64}
grid_size = 10::Int64
time_depth = 25::Int64
num_pairs = 50::Int64
asynchronus_weight = 100::Int64

generate_new_data = true
if generate_new_data == true

    # Generate ixi graph and extend it in time
    G = GridNetwork(grid_size, grid_size)

    # Extend graph in time without memory
    T = QuNet.TemporalGraph(G, time_depth, memory_prob=0.0)

    plot_data = []
    error_data = []
    for i in perc_range[1]:perc_range[2]:perc_range[3]
        println("Collecting for memory percolation rate : $i")
        raw_data = []
        for j in 1:num_trials

            # Make a network with some probability of memory
            T_mem = QuNet.TemporalGraph(G, time_depth, memory_prob=i, memory_costs = unit_costvector())
            # Make a copy of the network without memory
            T_copy = deepcopy(T)

            # Get i random userpairs with asynchronus src and dst nodes.
            mem_user_pairs = make_user_pairs(T_mem, num_pairs, src_layer=1, dst_layer=-1)
            user_pairs = make_user_pairs(T, num_pairs, src_layer=-1, dst_layer=-1)

            # Add async nodes
            QuNet.add_async_nodes!(T_mem, user_pairs, ϵ=asynchronus_weight)
            QuNet.add_async_nodes!(T_copy, user_pairs, ϵ=asynchronus_weight)

            # Get pathset data
            # NOTE: the line below used to have mem_user_pairs as an argument.
            pathset_mem, dum1, dum2 = QuNet.greedy_multi_path!(T_mem, QuNet.purify, mem_user_pairs, 1)
            pathset, dum1, dum2 = QuNet.greedy_multi_path!(T_copy, QuNet.purify, user_pairs, 1)
            # Pathset is an array of vectors containing edges describing paths between end-user pairs
            # Objective: find the largest timedepth used in the pathsets

            max_depth_mem = QuNet.max_timedepth(pathset_mem, T_mem)
            max_depth = QuNet.max_timedepth(pathset, T_copy)

            # Get the bandwidth of this quantity
            push!(raw_data, max_depth / max_depth_mem)
        end

        # Average the raw data, add it to plot data:
        push!(plot_data, mean(raw_data))
        # Get standard error
        push!(error_data, std(raw_data)/sqrt(num_trials - 1))
    end

    # Save data
    d = Dict{Symbol, Any}()
    @pack! d = num_trials, perc_range, grid_size, time_depth, num_pairs, asynchronus_weight, plot_data, error_data
    save("$datafile.jld", "data", d)

else
    # Load data
    if !isfile("$datafile.jld")
        error("$datafile.jld not found in working directory")
    end
    d = load("$datafile.jld")["data"]
    @unpack num_trials, perc_range, grid_size, time_depth, num_pairs, asynchronus_weight, plot_data, error_data = d
end

# Plot
x = collect(perc_range[1]:perc_range[2]:perc_range[3])
plot(x, plot_data, yerr = error_data, legend = false)
xaxis!(L"$\textrm{Proportion of Nodes with Quantum Memory}$")
yaxis!(L"$\textrm{R}$")

savefig("plots/bandwidth_with_memory_rate.png")
savefig("plots/bandwidth_with_memory_rate.pdf")

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
@testset "LoopPool" begin
    dim, N = 3, 4
    loopPool = ExprTree.LoopPool(:K, dim, N, Float64)
    basis1 = [1.0, 0.0, 0.0, 1.0]
    basis2 = [1.0, 1.0, 0.0, 0.0]
    basis3 = [1.0, 0.0, -1.0, 1.0]
    idx1 = ExprTree.append(loopPool, basis1)
    idx2 = ExprTree.append(loopPool, basis2)
    idx3 = ExprTree.append(loopPool, basis2)
    idx4 = ExprTree.append(loopPool, basis1)
    idx5 = ExprTree.append(loopPool, basis3)
    @test length(loopPool) == 3
    @test idx1 == idx4
    @test idx2 == idx3

    varK = rand(dim, N)
    ExprTree.update(loopPool, varK)
    @test ExprTree.current(loopPool, 1) ≈ varK * basis1
    @test ExprTree.current(loopPool, 2) ≈ varK * basis2
    @test ExprTree.current(loopPool, 3) ≈ varK * basis3

end

# @testset "CachedPool" begin
#     objType, weightType = Int, Float64
#     pool = ExprTree.CachedPool(:P, objType, weightType)
#     idx1 = ExprTree.append(pool, 1)
#     idx2 = ExprTree.append(pool, 2)
#     idx3 = ExprTree.append(pool, 2)
#     idx4 = ExprTree.append(pool, 1)
#     @test length(pool) == 2
#     @test idx1 == idx4
#     @test idx2 == idx3
# end

@testset "Generic Diagrams" begin

    """
        k1-k3                     k2+k3
        |                         |
    t1.L ↑     t1.L       t2.L     ↑ t2.L
        |-------------->----------|
        |       |  k3+k4   |      |
        |   v   |          |  v   |
        |       |    k4    |      |
        |--------------<----------|
    t1.L ↑    t1.L        t2.L     ↑ t2.L
        |                         |
        k1                        k2
    """
    # We only consider the direct part of the above diagram
    Gtype, Wtype = 1, 2
    spin = 2.0
    D = 3
    kF, β, mass2 = 1.919, 0.5, 1.0

    # varK = [rand(D) for i = 1:4] #k1, k2, k3, k4

    varK = rand(D, 4)
    varT = [rand() * β, rand() * β]

    K0 = [0.0, 0.0, 0.0]
    T0 = 0.0

    calcK(para, basis) = sum([para[i] .* basis[i] for i = 1:length(para)])
    calcT(para, basis) = para[basis[2]] - para[basis[1]]

    gorder, vorder = 0, 1

    weightType = Float64

    # function LoopPool(name::Symbol, dim::Int, N::Int, type::DataType)
    MomPool = ExprTree.LoopPool(:K, D, 4)

    diag = ExprTree.ExpressionTree(loopBasis = MomPool, nodePara = Int, weight = weightType)

    # #construct the propagator table
    gK = [[0.0, 0.0, 1.0, 1.0], [0.0, 0.0, 0.0, 1.0]]
    gT = [(1, 2), (2, 1)]
    g = [ExprTree.addpropagator!(diag, :G; site = gT[i], loop = gK[i], para = 1) for i = 1:2]

    vdK = [[0.0, 0.0, 1.0, 0.0], [0.0, 0.0, 1.0, 0.0]]
    vdT = [[1, 1], [2, 2]]
    vd = [ExprTree.addpropagator!(diag, :Vd; loop = vdK[i], para = 2) for i = 1:2]

    veK = [[1, 0, -1, -1], [0, 1, 0, -1]]
    veT = [[1, 1], [2, 2]]
    ve = [ExprTree.addpropagator!(diag, :Ve; loop = veK[i], para = 2) for i = 1:2]
    # ve = [ExprTree.addPropagator!(diag, Wtype, 1, veK[i], veT[i], Wsym)[1] for i = 1:2]
    # # W order is 1

    # # contruct the tree
    MUL, ADD = ExprTree.MUL, ExprTree.ADD
    ggn = ExprTree.addnode!(diag, MUL, :gxg, [g[1], g[2]], 1.0, para = 0)
    vdd = ExprTree.addnode!(diag, MUL, :dxd, [vd[1], vd[2]], spin, para = 0)
    vde = ExprTree.addnode!(diag, MUL, :dxe, [vd[1], ve[2]], -1.0, para = 0)
    ved = ExprTree.addnode!(diag, MUL, :exd, [ve[1], vd[2]], -1.0, para = 0)
    vsum = ExprTree.addnode!(diag, ADD, :sum, [vdd, vde, ved], 1.0, para = 0)
    root = ExprTree.addnode!(diag, MUL, :root, [ggn, vsum], 1.0, para = 0)
    push!(diag.root, root)
    ExprTree.initialize!(diag.node)

    # printBasisPool(diag)
    # printPropagator(diag)
    # ExprTree.printNodes(diag)
    # ExprTree.showTree(diag, diag.root[1])

    # #make sure the total number of diagrams are correct
    let
        DiagTree.eval(para, K, Tbasis, varT) = 1.0
        ExprTree.evalNaive!(diag, varK, varT)
        @test diag[1] ≈ -2 + 1 * spin
    end

    # #more sophisticated test of the weight evaluation
    let
        function evalG(K, τBasis, varT)
            ϵ = dot(K, K) / 2 - kF^2
            τ = varT[τBasis[2]] - varT[τBasis[1]]
            return Spectral.kernelFermiT(τ, ϵ, β)
        end

        evalV(K) = 8π / (dot(K, K) + mass2)

        function DiagTree.eval(para, K, Tbasis, varT)
            if para[1] == 1
                return evalG(K, Tbasis, varT)
            elseif para[1] == 2
                return evalV(K)
            else
                error("not implemented")
            end
        end

        # getK(basis, varK) = sum([basis[i] * K for (i, K) in enumerate(varK)])
        getK(basis, varK) = varK * basis

        gw = [evalG(getK(gK[i], varK), gT[i], varT) for i = 1:2]
        vdw = [evalV(getK(vdK[i], varK)) for i = 1:2]
        vew = [evalV(getK(veK[i], varK)) for i = 1:2]

        Vweight = spin * vdw[1] * vdw[2] - vdw[1] * vew[2] - vew[1] * vdw[2]
        Weight = gw[1] * gw[2] * Vweight

        # println(ExprTree.printPropagator(diag))
        ExprTree.evalNaive!(diag, varK, varT)
        @test diag[1] ≈ Weight
    end

end


text:
# Parent file to run all scripts which may generate
# some output that you want to display on the website.
# this can be used as a tester to check that all the code
# on your website runs properly.

dir = @__DIR__

"""
    genplain(s)

Small helper function to run some code and redirect the output (stdout) to a file.
"""
function genplain(s::String)
    open(joinpath(dir, "scripts/output", "$(splitext(s)[1]).out"), "w") do outf
        redirect_stdout(outf) do
            include(joinpath(dir, "scripts", s))
        end
    end
end

# output
using Glob
files = glob("./scripts/*.jl")
for f in files #[30:end]
    println("Generating output for "*f[11:end])
    genplain(f[11:end])
end

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
__precompile__(true)
module Unitful

import Base: ==, <, <=, +, -, *, /, .+, .-, .*, ./, .\, //, ^, .^
import Base: show, convert
import Base: abs, abs2, float, inv, sqrt
import Base: min, max, floor, ceil, log, log10, real, imag, conj

import Base: mod, rem, div, fld, cld, trunc, round, sign, signbit
import Base: isless, isapprox, isinteger, isreal, isinf, isfinite, isnan
import Base: copysign, flipsign
import Base: prevfloat, nextfloat, maxintfloat, rat, step #, linspace
import Base: promote_op, promote_array_type, promote_rule, unsafe_getindex
import Base: length, float, start, done, next, last, one, zero, colon#, range
import Base: getindex, eltype, step, last, first, frexp
import Base: Rational, typemin, typemax
import Base: steprange_last, unitrange_last, unsigned

export unit, dimension, uconvert, ustrip
export @dimension, @derived_dimension, @refunit, @unit, @u_str
export DimensionedQuantity, Quantity
export DimensionlessQuantity
export NoUnits, NoDims

include("Types.jl")
include("User.jl")
const NoUnits = Units{(), Dimensions{()}}()
const NoDims = Dimensions{()}()

(y::Units)(x::Number) = uconvert(y,x)

"""
```
type DimensionError <: Exception end
```

Thrown when dimensions don't match in an operation that demands they do.
"""
type DimensionError <: Exception end

"""
```
ustrip(x::Number)
```

Returns the number out in front of any units. This may be different from the value
in the case of dimensionless quantities. See [`uconvert`](@ref) and the example
below. Because the units are removed, information may be lost and this should
be used with some care.

This function is just calling `x/unit(x)`, which is as fast as directly
accessing the `val` field of `x::Quantity`, but also works for any other kind
of number.

This function is mainly intended for compatibility with packages that don't know
how to handle quantities. This function may be deprecated in the future.

```jldoctest
julia> ustrip(2u"μm/m") == 2
true

julia> uconvert(NoUnits, 2u"μm/m") == 2//1000000
true
```
"""
ustrip(x::Number) = x/unit(x)

"""
```
ustrip{T,D,U}(x::Array{Quantity{T,D,U}})
```

Strip units from an `Array` by reinterpreting to type `T`. The resulting
`Array` is a "unit free view" into array `x`. Because the units are
removed, information may be lost and this should be used with some care.

This function is provided primarily for compatibility purposes; you could pass
the result to PyPlot, for example. This function may be deprecated in the future.

```jldoctest
julia> a = [1u"m", 2u"m"]
2-element Array{Quantity{Int64, Dimensions:{𝐋}, Units:{m}},1}:
 1 m
 2 m

julia> b = ustrip(a)
2-element Array{Int64,1}:
 1
 2

julia> a[1] = 3u"m"; b
2-element Array{Int64,1}:
 3
 2
```
"""
ustrip{T,D,U}(x::Array{Quantity{T,D,U}}) = reinterpret(T, x)

"""
```
ustrip{T<:Number}(x::Array{T})
```

Fall-back that returns `x`.
"""
ustrip{T<:Number}(x::Array{T}) = x

"""
```
unit{T,D,U}(x::Quantity{T,D,U})
```

Returns the units associated with a quantity, `U()`.

Examples:

```jldoctest
julia> unit(1.0u"m") == u"m"
true

julia> typeof(u"m")
Unitful.Units{(Unitful.Unit{:Meter}(0,1//1),),Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}}
```
"""
unit{T,D,U}(x::Quantity{T,D,U}) = U()

"""
```
unit{T,D,U}(x::Type{Quantity{T,D,U}})
```

Returns the units associated with a quantity type, `U()`.

Examples:

```jldoctest
julia> unit(typeof(1.0u"m")) == u"m"
true
```
"""
unit{T,D,U}(::Type{Quantity{T,D,U}}) = U()

"""
```
unit(x::Number)
```

Returns a `Unitful.Units{(), Dimensions{()}}` object to indicate that ordinary
numbers have no units. This is a singleton, which we export as `NoUnits`.
The unit is displayed as an empty string.

Examples:

```jldoctest
julia> typeof(unit(1.0))
Unitful.Units{(),Unitful.Dimensions{()}}
julia> typeof(unit(Float64))
Unitful.Units{(),Unitful.Dimensions{()}}
julia> unit(1.0) == NoUnits
true
```
"""
unit(x::Number) = NoUnits
unit{T<:Number}(x::Type{T}) = NoUnits

"""
```
dimension(x::Number)
dimension{T<:Number}(x::Type{T})
```

Returns a `Unitful.Dimensions{()}` object to indicate that ordinary
numbers are dimensionless. This is a singleton, which we export as `NoDims`.
The dimension is displayed as an empty string.

Examples:

```jldoctest
julia> typeof(dimension(1.0))
Unitful.Dimensions{()}
julia> typeof(dimension(Float64))
Unitful.Dimensions{()}
julia> dimension(1.0) == NoDims
true
```
"""
dimension(x::Number) = NoDims
dimension{T<:Number}(x::Type{T}) = NoDims

"""
```
dimension{U,D}(u::Units{U,D})
```

Returns a [`Unitful.Dimensions`](@ref) object corresponding to the dimensions
of the units, `D()`. For a dimensionless combination of units, a
`Unitful.Dimensions{()}` object is returned.

Examples:

```jldoctest
julia> dimension(u"m")
𝐋

julia> typeof(dimension(u"m"))
Unitful.Dimensions{(Unitful.Dimension{:Length}(1//1),)}

julia> typeof(dimension(u"m/km"))
Unitful.Dimensions{()}
```
"""
dimension{U,D}(u::Units{U,D}) = D()

"""
```
dimension{D}(x::DimensionedQuantity{D})
```

Returns a [`Unitful.Dimensions`](@ref) object `D()` corresponding to the
dimensions of quantity `x`. For a dimensionless [`Unitful.Quantity`](@ref), a
`Unitful.Dimensions{()}` object is returned.

Examples:

```jldoctest
julia> dimension(1.0u"m")
𝐋

julia> typeof(dimension(1.0u"m/μm"))
Unitful.Dimensions{()}
```
"""
dimension{D}(x::DimensionedQuantity{D}) = D()
dimension{D}(x::Type{DimensionedQuantity{D}}) = D()
dimension{T,D,U}(::Type{Quantity{T,D,Units{U,D}}}) = D()

"""
```
dimension{T<:Number}(x::AbstractArray{T})
```

Just calls `map(dimension, x)`.
"""
dimension{T<:Number}(x::AbstractArray{T}) = map(dimension, x)

"""
```
dimension{T<:Units}(x::AbstractArray{T})
```

Just calls `map(dimension, x)`.
"""
dimension{T<:Units}(x::AbstractArray{T}) = map(dimension, x)

function basefactorhelper(inex, ex, t, p)
    if isinteger(p)
        p = Integer(p)
    end

    can_exact = (ex < typemax(Int))
    can_exact &= (1/ex < typemax(Int))

    ex2 = 10.0^t*float(ex)^p
    can_exact &= (ex2 < typemax(Int))
    can_exact &= (1/ex2 < typemax(Int))
    can_exact &= isinteger(p)

    if can_exact
        (inex^p, (ex//1*10^t)^p)
    else
        ((inex * ex * 10.0^t)^p, 1)
    end
end

@generated function basefactor(x::Units)
    tunits = x.parameters[1]
    fact1 = map(basefactor, tunits)
    inex1 = mapreduce(x->getfield(x,1), *, 1.0, fact1)
    ex1   = mapreduce(x->getfield(x,2), *, 1, fact1)
    :(($inex1,$ex1))
end

function tensfactor(x::Unit)
    p = power(x)
    if isinteger(p)
        p = Integer(p)
    end
    tens(x)*p
end

@generated function tensfactor(x::Units)
    tunits = x.parameters[1]
    a = mapreduce(tensfactor, +, 0, tunits)
    :($a)
end

# Addition / subtraction
for op in [:+, :-]
    @eval ($op){S,T,D,U}(x::Quantity{S,D,U}, y::Quantity{T,D,U}) =
        Quantity(($op)(x.val,y.val), U())

    # If not generated, there are run-time allocations
    @eval @generated function ($op){S,T,D,SU,TU}(x::Quantity{S,D,SU},
            y::Quantity{T,D,TU})
        result_units = promote_type(SU,TU)()
        :($($op)(uconvert($result_units, x), uconvert($result_units, y)))
    end

    @eval ($op)(::Quantity, ::Quantity) = throw(DimensionError())
    @eval function ($op)(x::Quantity, y::Number)
        if isa(x, DimensionlessQuantity)
            ($op)(promote(x,y)...)
        else
            throw(DimensionError())
        end
    end
    @eval function ($op)(x::Number, y::Quantity)
        if isa(y, DimensionlessQuantity)
            ($op)(promote(x,y)...)
        else
            throw(DimensionError())
        end
    end

    @eval ($op)(x::Quantity) = Quantity(($op)(x.val),unit(x))
end

*(x::Number, y::Units, z::Units...) = Quantity(x,*(y,z...))

# Kind of weird, but okay, no need to make things noncommutative.
*(x::Units, y::Number) = *(y,x)

# These six are defined for use in `*(a0::Unitlike, a::Unitlike...)`
unit{S}(x::Unit{S}) = S
unit{S}(x::Dimension{S}) = S
tens(x::Unit) = x.tens
tens(x::Dimension) = 0
power(x::Unit) = x.power
power(x::Dimension) = x.power

"""
```
*(a0::Unitlike, a::Unitlike...)
```

Given however many unit-like objects, multiply them together. Both
[`Unitful.Dimensions`](@ref) and [`Unitful.Units`](@ref) objects are considered
to be `Unitlike` in the sense that you can multiply them, divide them, and
collect powers. This function will fail if there is an attempt to multiply a
unit by a dimension or vice versa.

Collect [`Unitful.Unit`](@ref) objects from the type parameter of the
[`Unitful.Units`](@ref) objects. For identical units including SI prefixes
(i.e. cm ≠ m), collect powers and sort uniquely by the name of the `Unit`.
The unique sorting permits easy unit comparisons.

Examples:

```jldoctest
julia> u"kg*m/s^2"
kg m s^-2

julia> u"m/s*kg/s"
kg m s^-2

julia> typeof(u"m/s*kg/s") == typeof(u"kg*m/s^2")
true
```
"""
@generated function *(a0::Unitlike, a::Unitlike...)

    # Sort the units uniquely. This is a generated function so that we
    # don't have to figure out the units each time.

    D = a0 <: Units ? Unit : Dimension
    b = Array{D,1}()
    a0p = a0.parameters[1]
    length(a0p) > 0 && append!(b, a0p)
    for x in a
        xp = x.parameters[1]
        length(xp) > 0 && append!(b, xp)
    end

    # b is an Array containing all of the Unit or Dimension objects that were
    # found in the type parameters of the Units or Dimensions object (a0, a...)

    sort!(b, by=x->power(x))
    D == Unit && sort!(b, by=x->tens(x))
    sort!(b, by=x->unit(x))

    # Units[m,m,cm,cm^2,cm^3,nm,m^4,µs,µs^2,s]
    # reordered as:
    # Units[nm,cm,cm^2,cm^3,m,m,m^4,µs,µs^2,s]

    # Collect powers of a given unit
    c = Array{D,1}()
    if !isempty(b)
        i = start(b)
        oldstate = b[i]
        p=0//1
        while !done(b, i)
            (state, i) = next(b, i)
            if tens(state) == tens(oldstate) && unit(state) == unit(oldstate)
                p += power(state)
            else
                if p != 0
                    push!(c, D{unit(oldstate)}(tens(oldstate),p))
                end
                p = power(state)
            end
            oldstate = state
        end
        if p != 0
            push!(c, D{unit(oldstate)}(tens(oldstate),p))
        end
    end
    # results in:
    # Units[nm,cm^6,m^6,µs^3,s]

    if a0 <: Units
        T = Units
        d = (c...)
        f = typeof(mapreduce(dimension,*,NoDims,c))
        :(($T){$d,$f}())
    else
        T = Dimensions
        d = (c...)
        :(($T){$d}())
    end
end

function *{T,D,U}(x::Quantity{T,D,U}, y::Units, z::Units...)
    result_units = *(U(),y,z...)
    Quantity(x.val,result_units)
end

function *(x::Quantity, y::Quantity)
    xunits = unit(x)
    yunits = unit(y)
    result_units = xunits*yunits
    z = x.val*y.val
    Quantity(z,result_units)
end

# Next two lines resolves some method ambiguity:
*{T<:Quantity}(x::Bool, y::T) =
    ifelse(x, y, ifelse(signbit(y), -zero(y), zero(y)))
*(x::Quantity, y::Bool) = Quantity(x.val*y, unit(x))

*(y::Number, x::Quantity) = *(x,y)
*(x::Quantity, y::Number) = Quantity(x.val*y, unit(x))

# See operators.jl
# Element-wise operations with units
for (f,F) in [(:./, :/), (:.*, :*), (:.+, :+), (:.-, :-)]
    @eval ($f)(x::Units, y::Units) = ($F)(x,y)
    @eval ($f)(x::Number, y::Units)   = ($F)(x,y)
    @eval ($f)(x::Units, y::Number)   = ($F)(x,y)
end
.\(x::Unitlike, y::Unitlike) = y./x
.\(x::Number, y::Units) = y./x
.\(x::Units, y::Number) = y./x

# See arraymath.jl
./(x::Units, Y::AbstractArray) =
    reshape([ x ./ y for y in Y ], size(Y))
./(X::AbstractArray, y::Units) =
    reshape([ x ./ y for x in X ], size(X))
.\(x::Units, Y::AbstractArray) =
    reshape([ x .\ y for y in Y ], size(Y))
.\(X::AbstractArray, y::Units) =
    reshape([ x .\ y for x in X ], size(X))

for f in (:.*,) # looked in arraymath.jl for similar code
    @eval begin
        function ($f){T}(A::Units, B::AbstractArray{T})
            F = similar(B, promote_op($f,typeof(A),T))
            for (iF, iB) in zip(eachindex(F), eachindex(B))
                @inbounds F[iF] = ($f)(A, B[iB])
            end
            return F
        end
        function ($f){T}(A::AbstractArray{T}, B::Units)
            F = similar(A, promote_op($f,T,typeof(B)))
            for (iF, iA) in zip(eachindex(F), eachindex(A))
                @inbounds F[iF] = ($f)(A[iA], B)
            end
            return F
        end
    end
end

# Division (units)

/(x::Unitlike, y::Unitlike) = *(x,inv(y))
/(x::Quantity, y::Units) = Quantity(x.val, unit(x) / y)
/(x::Units, y::Quantity) = Quantity(1/y.val, x / unit(y))
/(x::Number, y::Units) = Quantity(x,inv(y))
/(x::Units, y::Number) = (1/y) * x

//(x::Unitlike, y::Unitlike)  = x/y
//(x::Quantity, y::Units) = Quantity(x.val, unit(x) / y)
//(x::Units, y::Quantity) = Quantity(1//y.val, x / unit(y))
//(x::Number, y::Units) = Rational(x)/y
//(x::Units, y::Number) = (1//y) * x

# Division (quantities)

for op in (:/, ://)
    @eval begin
        ($op)(x::Quantity, y::Quantity) = Quantity(($op)(x.val, y.val), unit(x) / unit(y))
        ($op)(x::Quantity, y::Number) = Quantity(($op)(x.val, y), unit(x))
        ($op)(x::Number, y::Quantity) = Quantity(($op)(x, y.val), inv(unit(y)))
    end
end

# ambiguity resolution
//(x::Quantity, y::Complex) = Quantity(//(x.val, y), unit(x))

# Division (other functions)

for f in (:div, :fld, :cld)
    @eval function ($f)(x::Quantity, y::Quantity)
        z = uconvert(unit(y), x)
        ($f)(z.val,y.val)
    end
end

for f in (:mod, :rem)
    @eval function ($f)(x::Quantity, y::Quantity)
        z = uconvert(unit(y), x)
        Quantity(($f)(z.val,y.val), unit(y))
    end
end

# Exponentiation is not type stable.
# For now we define a special `inv` method to at least
# enable division to be fast.
@generated function inv(x::Dimensions)
    tup = x.parameters[1]
    length(tup) == 0 && return :(x)
    tup2 = map(x->x^-1,tup)
    y = *(Dimensions{tup2}())
    :($y)
end

@generated function inv(x::Units)
    tup = x.parameters[1]
    length(tup) == 0 && return :(x)
    tup2 = map(x->x^-1,tup)
    D = typeof(mapreduce(dimension, *, NoDims, tup2))
    y = *(Units{tup2, D}())
    :($y)
end

^{T}(x::Unit{T}, y::Integer) = Unit{T}(tens(x),power(x)*y)
^{T}(x::Unit{T}, y) = Unit{T}(tens(x),power(x)*y)

^{T}(x::Dimension{T}, y::Integer) = Dimension{T}(power(x)*y)
^{T}(x::Dimension{T}, y) = Dimension{T}(power(x)*y)

function ^{T}(x::Dimensions{T}, y::Integer)   # needed for ambiguity resolution
    *(Dimensions{map(a->a^y, T)}())
end

function ^{T}(x::Dimensions{T}, y)
    *(Dimensions{map(a->a^y, T)}())
end

function ^{U,D}(x::Units{U,D}, y::Integer)
    utup = map(a->a^y, U)
    *(Units{utup, ()}()) # dimensions get reconstructed anyway
    # would use mapreduce(dimension, *, NoDims, utup)
end

function ^{U,D}(x::Units{U,D}, y)
    utup = map(a->a^y, U)
    *(Units{utup, ()}())
end

# All of these are needed for ambiguity resolution
^{T,D,U}(x::Quantity{T,D,U}, y::Integer) = Quantity((x.val)^y, U()^y)
^{T,D,U}(x::Quantity{T,D,U}, y::Rational) = Quantity((x.val)^y, U()^y)
^{T,D,U}(x::Quantity{T,D,U}, y::Real) = Quantity((x.val)^y, U()^y)

# Other mathematical functions
"""
```
sqrt(x::Quantity)
```
"""
sqrt(x::Quantity) = Quantity(sqrt(x.val), sqrt(unit(x)))

# This is a generated function to ensure type stability and keep `sqrt` fast.
@generated function sqrt(x::Dimensions)
    tup = x.parameters[1]
    tup2 = map(x->x^(1//2),tup)
    y = *(Dimensions{tup2}())    # sort appropriately
    :($y)
end

# This is a generated function to ensure type stability and keep `sqrt` fast.
@generated function sqrt(x::Units)
    tup = x.parameters[1]
    tup2 = map(x->x^(1//2),tup)
    y = *(Units{tup2,()}())    # sort appropriately
    :($y)
end

for (f, F) in [(:min, :<), (:max, :>)]
    @eval @generated function ($f)(x::Quantity, y::Quantity)
        xdim = x.parameters[2]()
        ydim = y.parameters[2]()
        if xdim != ydim
            return :(throw(DimensionError()))
        end

        xunits = x.parameters[3].parameters[1]
        yunits = y.parameters[3].parameters[1]

        factx = mapreduce(.*, xunits) do x
            vcat(basefactor(x)...)
        end
        facty = mapreduce(.*, yunits) do x
            vcat(basefactor(x)...)
        end

        tensx = mapreduce(tensfactor, +, xunits)
        tensy = mapreduce(tensfactor, +, yunits)

        convx = *(factx..., (10.0)^tensx)
        convy = *(facty..., (10.0)^tensy)

        :($($F)(x.val*$convx, y.val*$convy) ? x : y)
    end

    @eval ($f)(x::Units, y::Units) =        # TODO remove
        unit(($f)(Quantity(1.0, x), Quantity(1.0, y)))
end

@vectorize_2arg Quantity max
@vectorize_2arg Quantity min

abs(x::Quantity) = Quantity(abs(x.val),  unit(x))
abs2(x::Quantity) = Quantity(abs2(x.val), unit(x)*unit(x))

trunc(x::Quantity) = Quantity(trunc(x.val), unit(x))
round(x::Quantity) = Quantity(round(x.val), unit(x))

copysign(x::Quantity, y::Number) = Quantity(copysign(x.val,y/unit(y)), unit(x))
flipsign(x::Quantity, y::Number) = Quantity(flipsign(x.val,y/unit(y)), unit(x))

isless{T,D,U}(x::Quantity{T,D,U}, y::Quantity{T,D,U}) = isless(x.val, y.val)
isless(x::Quantity, y::Quantity) = isless(uconvert(unit(y), x).val, y.val)
<{T,D,U}(x::Quantity{T,D,U}, y::Quantity{T,D,U}) = (x.val < y.val)
<(x::Quantity, y::Quantity) = <(uconvert(unit(y), x).val,y.val)

isapprox{T,D,U}(x::Quantity{T,D,U}, y::Quantity{T,D,U}) = isapprox(x.val, y.val)
isapprox(x::Quantity, y::Quantity) = isapprox(uconvert(unit(y), x).val, y.val)
isapprox(x::Quantity, y::Number) = isapprox(uconvert(NoUnits, x), y)
isapprox(x::Number, y::Quantity) = isapprox(y,x)

function isapprox{T1,D,U1,T2,U2}(x::AbstractArray{Quantity{T1,D,U1}},
        y::AbstractArray{Quantity{T2,D,U2}}; rtol::Real=Base.rtoldefault(T1,T2),
        atol=zero(Quantity{T1,D,U1}), norm::Function=vecnorm)

    d = norm(x - y)
    if isfinite(d)
        return d <= atol + rtol*max(norm(x), norm(y))
    else
        # Fall back to a component-wise approximate comparison
        return all(ab -> isapprox(ab[1], ab[2]; rtol=rtol, atol=atol), zip(x, y))
    end
end
isapprox{S<:Quantity,T<:Quantity}(x::AbstractArray{S}, y::AbstractArray{T};
    kwargs...) = false
function isapprox{S<:Quantity,N<:Number}(x::AbstractArray{S}, y::AbstractArray{N};
    kwargs...)
    if dimension(N) == dimension(S)
        isapprox(map(x->uconvert(NoUnits,x),x),y; kwargs...)
    else
        false
    end
end
isapprox{S<:Quantity,N<:Number}(y::AbstractArray{N}, x::AbstractArray{S};
    kwargs...) = isapprox(x,y; kwargs...)

=={S,T,D,U}(x::Quantity{S,D,U}, y::Quantity{T,D,U}) = (x.val == y.val)
function ==(x::Quantity, y::Quantity)
    dimension(x) != dimension(y) && return false
    uconvert(unit(y), x).val == y.val
end

function ==(x::Quantity, y::Number)
    if dimension(x) == NoDims
        uconvert(NoUnits, x) == y
    else
        false
    end
end
==(x::Number, y::Quantity) = ==(y,x)
<=(x::Quantity, y::Quantity) = <(x,y) || x==y

for f in (:zero, :floor, :ceil)
    @eval ($f)(x::Quantity) = Quantity(($f)(x.val), unit(x))
end
zero{T,D,U}(x::Type{Quantity{T,D,U}}) = zero(T)*U()

"""
```
one(x::Quantity)
```

Returns the multiplicative identity for `x`.
"""
one(x::Quantity) = one(x.val)

"""
```
one{T,D,U}(x::Type{Quantity{T,U}})
```

Returns the multiplicative identity for this type (it's `one(T)`).
"""
one{T,D,U}(x::Type{Quantity{T,D,U}}) = one(T)

isinteger(x::Quantity) = isinteger(x.val)
isreal(x::Quantity) = isreal(x.val)
isfinite(x::Quantity) = isfinite(x.val)
isinf(x::Quantity) = isinf(x.val)
isnan(x::Quantity) = isnan(x.val)

unsigned(x::Quantity) = Quantity(unsigned(x.val), unit(x))

log(x::DimensionlessQuantity) = log(uconvert(NoUnits, x))
log10(x::DimensionlessQuantity) = log10(uconvert(NoUnits, x))

real(x::Quantity) = Quantity(real(x.val), unit(x))
imag(x::Quantity) = Quantity(imag(x.val), unit(x))
conj(x::Quantity) = Quantity(conj(x.val), unit(x))

"""
```
sign(x::Quantity)
```

Returns the sign of `x`.
"""
sign(x::Quantity) = sign(x.val)

"""
```
signbit(x::Quantity)
```

Returns the sign bit of the underlying numeric value of `x`.
"""
signbit(x::Quantity) = signbit(x.val)

prevfloat{T<:AbstractFloat,D,U}(x::Quantity{T,D,U}) =
    Quantity(prevfloat(x.val), unit(x))
nextfloat{T<:AbstractFloat,D,U}(x::Quantity{T,D,U}) =
    Quantity(nextfloat(x.val), unit(x))

function frexp{T<:AbstractFloat,D,U}(x::Quantity{T,D,U})
    a,b = frexp(x.val)
    a *= unit(x)
    a,b
end

"""
```
float(x::Quantity)
```

Convert the numeric backing type of `x` to a floating-point representation.
Returns a `Quantity` with the same units.
"""
float(x::Quantity) = Quantity(float(x.val), unit(x))

"""
```
Integer(x::Quantity)
```

Convert the numeric backing type of `x` to an integer representation.
Returns a `Quantity` with the same units.
"""
Integer(x::Quantity) = Quantity(Integer(x.val), unit(x))

"""
```
Rational(x::Quantity)
```

Convert the numeric backing type of `x` to a rational number representation.
Returns a `Quantity` with the same units.
"""
Rational(x::Quantity) = Quantity(Rational(x.val), unit(x))

colon(start::Quantity, step::Quantity, stop::Quantity) =
    StepRange(promote(start, step, stop)...)

function Base.steprange_last{T<:Number,D,U}(start::Quantity{T,D,U}, step, stop)
    z = zero(step)
    step == z && throw(ArgumentError("step cannot be zero"))
    if stop == start
        last = stop
    else
        if (step > z) != (stop > start)
            last = start - step
        else
            diff = stop - start
            if T<:Signed && (diff > zero(diff)) != (stop > start)
                # handle overflowed subtraction with unsigned rem
                if diff > zero(diff)
                    remain = -convert(typeof(start), unsigned(-diff) % step)
                else
                    remain = convert(typeof(start), unsigned(diff) % step)
                end
            else
                remain = Base.steprem(start,stop,step)
            end
            last = stop - remain
        end
    end
    last
end

"""
```
offsettemp(::Unit)
```

For temperature units, this function is used to set the scale offset.
"""
offsettemp(::Unit) = 0

include("Display.jl")
include("Promotion.jl")
include("Conversion.jl")

defaults()
end


text:
mime_types = Dict(
    "css" => "text/css",
    "gif" => "image/gif",
    "ico" => "image/x-icon",
    "jpe" => "image/jpeg",
    "jpeg" => "image/jpeg",
    "jpg" => "image/jpeg",
    "js" => "text/javascript",
    "pdf" => "application/pdf",
    "svg" => "image/svg+xml",
    "txt" => "text/plain"
)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "ASLBuilder"
version = v"3.1.0"

# Collection of sources required to build ASLBuilder
sources = [
    "https://github.com/ampl/mp/archive/3.1.0.tar.gz" =>
    "587c1a88f4c8f57bef95b58a8586956145417c8039f59b1758365ccc5a309ae9",
    "./bundled"
]

# Bash recipe for building across all platforms
script = raw"""
# Use staticfloat's cross-compile trick for ASL https://github.com/ampl/mp/issues/115
cd $WORKSPACE/srcdir/mp-3.1.0
rm -rf thirdparty/benchmark
patch -p1 < $WORKSPACE/srcdir/mp-extra/no_benchmark.patch

# Build ASL

mkdir build
cd build

## STATIC BUILD START
if [ $target = "x86_64-apple-darwin14" ]; then
  export AR=/opt/x86_64-apple-darwin14/bin/x86_64-apple-darwin14-ar
fi

if [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then
   cmake -DCMAKE_C_FLAGS='-fPIC -DPIC' -DBUILD_SHARED_LIBS=OFF -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DRUN_HAVE_STD_REGEX=0 -DRUN_HAVE_STEADY_CLOCK=0 -DHAVE_ACCESS_DRIVER_EXITCODE=0 -DHAVE_EXCEL_DRIVER_EXITCODE=0 -DHAVE_ODBC_TEXT_DRIVER_EXITCODE=0    ../
else
   cmake -DCMAKE_C_FLAGS='-fPIC -DPIC' -DBUILD_SHARED_LIBS=OFF -DCMAKE_INSTALL_PREFIX=$prefix  -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain       -DRUN_HAVE_STD_REGEX=0       -DRUN_HAVE_STEADY_CLOCK=0       ../
fi
## STATIC BUILD END

## DYNAMIC BUILD START
#if [ $target = "x86_64-w64-mingw32" ] || [ $target = "i686-w64-mingw32" ]; then
#   cmake -DCMAKE_C_FLAGS='-fPIC -DPIC' -DBUILD_SHARED_LIBS=ON -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DRUN_HAVE_STD_REGEX=0 -DRUN_HAVE_STEADY_CLOCK=0 -DHAVE_ACCESS_DRIVER_EXITCODE=0 -DHAVE_EXCEL_DRIVER_EXITCODE=0 -DHAVE_ODBC_TEXT_DRIVER_EXITCODE=0    ../
#else
#   cmake -DCMAKE_C_FLAGS='-fPIC -DPIC' -DBUILD_SHARED_LIBS=ON -DCMAKE_INSTALL_PREFIX=$prefix  -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain       -DRUN_HAVE_STD_REGEX=0       -DRUN_HAVE_STEADY_CLOCK=0       ../
#fi
## DYNAMIC BUILD END

# Copy over pregenerated files after building arithchk, so as to fake out cmake,
# because cmake will delete our arith.h

## If this fails it is ok, we just want to prevend cmake deleting arith.h
set +e
make arith-h VERBOSE=1
set -e

mkdir -p src/asl
cp -v $WORKSPACE/srcdir/mp-extra/expr-info.cc ../src/expr-info.cc
cp -v $WORKSPACE/srcdir/mp-extra/arith.h.${target} src/asl/arith.h

# Build and install ASL

make -j${nproc} VERBOSE=1
make install VERBOSE=1
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libasl", :libasl)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using POMDPs
using DESPOT
using POMDPToolbox
using POMDPModels

import DESPOT: bounds

struct BabyBounds end
bounds{S}(::BabyBounds, p::BabyPOMDP, ::Vector{DESPOTParticle{S}}, ::DESPOTConfig) = (p.r_feed+p.r_hungry)/(1.0-discount(p)), 0

solver = DESPOTSolver{Bool, Bool, Bool, BabyBounds, MersenneStreamArray}(bounds = BabyBounds(),
                                                    random_streams=MersenneStreamArray(MersenneTwister(1)),
                                                    next_state=false,
                                                    curr_obs=false,
                                                    rng=Base.GLOBAL_RNG,
                                                    max_trials=100
                                                    )
problem = BabyPOMDP()

# test_solver is from POMDPToolbox
test_solver(solver, problem, updater=updater(problem))

blink(solver)


text:
function transform(qs, ks, vs, bags)
	d = size(qs,1)
	os = map(bags) do b
		isempty(b) && return(vs[:,0:-1])
		q, k, v = qs[:,b], ks[:,b], vs[:,b]
		v * softmax(q' * k ./ sqrt(d))
	end
	reduce(catobs, os)
end

struct SingleHeadTransform{Q,K,V}
	q::Q
	k::K
	v::V
	n::Flux.LayerNorm
end

SingleHeadTransform(d) = SingleHeadTransform(Dense(d,d),Dense(d,d),Dense(d,d), LayerNorm(d))

Flux.@functor(SingleHeadTransform)

(m::SingleHeadTransform)(::Missing, bags) = missing

function (m::SingleHeadTransform)(x, bags)
	o = transform(m.q(x), m.k(x), m.v(x), bags)
	m.n(x + o)
end


text:
"""
Table 1 from Efron-Tibshirani 1996
"""
const pain_scores = [
    3
    7
    6
    1
    2
    3
    3
    1
    7
    5
    4
    4
    1
    3
    3
    5
    0
    1
    0
    0
    2
    2
    0
    0
    0
    1
    0
    0
    1
    1
    1
    0
    0
    0
    0
    0
    0
    0
    0
    0
]


text:
"""
# Hybrid-λ (Daw et al. 2011)
# --------------------------

Hybrid: Forward Learner + SARSA-λ

See learner/forwardlearner.jl, learner/sarsaLlearner.jl
and learner/hybridconvexlearner.jl for the core implementation.

"""
function getLL_HybridSarsaLambda(data::Array{Int,2}, alpha::Float64, gamma::Float64,
                                lambda::Float64, eta::Float64, decay_offset::Float64,
                                decay_slope::Float64, temp::Float64)
    data_stats, nr_states, nr_actions, nr_samples = get_data_stats_and_dimensions(data)
    if data_stats.nr_subjects != 1
        error("Invalid argument: getLL_HybridSarsaLambda_biased supports only one subject")
    end

    LL_unbiased = 0.
    action_selection_probs = zeros(Float64, nr_samples)

    # signals:
    RPE_record = zeros(Float64, nr_samples)
    SPE_record = zeros(Float64, nr_samples)
    Qsa_hyb_record = zeros(Float64, nr_samples)
    learnerMB = []
    learnerMF = []
    learnerHybrid = []
    current_puzzle_id = -1
    # Puzzle = independent block of the experiment.
    # The code was made for general purposes, in this experiment we
    # have only ONE puzzle.
    isendofexperiment = false
    modulatingfactor = 0.
    for i = 1:nr_samples
        if current_puzzle_id != data[i,Int(c_puzzle)]
            # if new puzzle -> reset learning
            current_puzzle_id = data[i,Int(c_puzzle)]
            learnerMB = ForwardLearner(nr_states = nr_states,
                                nr_actions = nr_actions,
                                η = eta, γ = gamma,
                                init_R = 0.
                                )
            learnerMF = SarsaL(nr_states = nr_states,
                                nr_actions = nr_actions,
                                α = alpha,
                                γ = gamma,
                                λ = lambda,
                                initvalue = 0.
                                )
            learnerHybrid = HybridConvex(nr_states = nr_states,
                                nr_actions = nr_actions,
                                decay_offset = decay_offset,
                                decay_slope = decay_slope
                                )
        end
        state = data[i, Int(c_state)]
        action = data[i, Int(c_action)]
        reward = data[i, Int(c_reward)]
        next_state = data[i, Int(c_nextState)]

        # Hybrid learner
        update!(learnerHybrid, state, action, reward, learnerMB.QSA, learnerMF.QSA) #Q_hybrid =  w_MB*learnerMB.QSA + w_MF*learnerMF.QSA
        Qsa_hyb_record[i] = copy(learnerHybrid.QSA[state, action])

        p_actions = getSoftmaxProbs(learnerHybrid.QSA[state,:], temp, [0., 0.])
        action_selection_probs[i] = copy(p_actions[action])
        LL_unbiased += log(p_actions[action])

        # ForwardLearner
        update!(learnerMB, state, action, next_state, reward)
        SPE_record[i] = copy(learnerMB.spe[1])

        # Sort out if end of episode or end of experiment
        next_action = 0
        if reward==0 #non goal state
            if nr_samples > i
                next_action = data[i+1,Int(c_action)]
            else
                isendofexperiment = true
            end
        end

        # Sarsa-Lambda with replacing traces
        update!(learnerMF, state, action, next_state, next_action, reward,
                isendofexperiment, modulatingfactor)
        RPE_record[i] = copy(learnerMF.rpe[1])
    end

    LL = LL_unbiased
    trial_by_trial_p_sa = (action_selection_probs)
    algo = "Hybrid FWD+SarsaL"
    results_dict = Dict{String,Any}("Qsa_hyb_record" => Qsa_hyb_record,
                                    "finalQSA_Hybrid_record" => learnerHybrid.QSA,
                                    "finalQSA_MF_record" => learnerMF.QSA,
                                    "finalQSA_MB_record" => learnerMB.QSA,
                                    "final_Tsas" => learnerMB.TSAS)
    param_v = [alpha, lambda, gamma, eta, decay_offset, decay_slope, temp]
    parma_str = "alpha=$(round(alpha,digits=4)), lambda=$(round(lambda,digits=4)),
                gamma=$(round(gamma,digits=4)), eta=$(round(eta,digits=4)),
                decay_offset=$(round(decay_offset,digits=4)),
                decay_slope=$(round(decay_slope,digits=4)), temp=$(round(temp,digits=4))"

    signals = Step_by_step_signals(data_stats.all_subject_ids[1],
        RPE_record,
        SPE_record,
        trial_by_trial_p_sa,
        LL,
        data_stats,
        algo,
        param_v,
        parma_str,
        99, #biased_action_idx, dummy value
        results_dict)
        return LL, signals
end


text:
include("help.jl")
function mark_kross!(r::Robot)     #Основная функция
    for side in (Nord,West,Sud,Ost) #Идем во все стороны
        num_steps=get_num_steps_putmarkers!(r,side) #запомнимаем число шагов
        movements!(r,invers(side),num_steps) #Возвращаемся
    end
    putmarker!(r) #Ставим маркеры
end

invers(side::HorizonSide) = HorizonSide(mod(Int(side) + 2,4)) #Меняем направление на противоположное

function get_num_steps_putmarkers!(r::Robot,side::HorizonSide) #Идем в side и запоминаем шаги
    num_steps=0
    while isborder(r, side) == false #пока нет стенки
        move!(r,side)
        putmarker!(r)
        num_steps += 1
    end
    return num_steps #Возвращаем число шагов
end

max_stars_repo_name: MCPavlo/mirea-progs

max_issues_repo_name: MCPavlo/mirea-progs

max_forks_repo_name: MCPavlo/mirea-progs


text:
# module_rewrite.jl - rewriting the function as a verilog module generator

################################################################################
## MODULE VERSION.
## complete rewrite of the function expression
##

function strip_output_definition!(f::Expr)
  #if the function already has an expression type
  if f.args[1].head == :(::)
    f.args[1] = f.args[1].args[1]# = output_type
  end
end

function strip_wire_inputs!(f::Expr)
  #issue a list of places where we need extra wire macros.
  inputs_list = Vector{Tuple{Symbol, Any}}()

  #make a copy of the new arguments list.
  new_arguments = deepcopy(f.args[1].args[1:1])
  param_list = Symbol[]

  for argument in f.args[1].args[2:end]
    #a case where we have an untyped argument
    if isa(argument, Symbol)
      push!(param_list, argument)
      push!(new_arguments, argument)
      continue
    end

    #a bare wire statement... Presume a @input macro statement already exists.
    if argument.args[2] == :Wire
      #push a blank input symbol to the inputs_list stack.
      push!(inputs_list, (argument.args[1], nothing))
    elseif argument.args[2] == :SingleWire
      push!(inputs_list, (argument.args[1], :(0:0v)))
    elseif isa(argument.args[2], Expr) &&
      (argument.args[2].head == :curly) &&
      (argument.args[2].args[1] == :Wire) &&
      (argument.args[2].args[2].head == :(:))
      #save the needed symbol and wire range parameter.
      push!(inputs_list, (argument.args[1], argument.args[2].args[2]))
    else
      push!(param_list, argument.args[1])
      push!(new_arguments, argument)
    end
  end

  #take the old arguments list and make them the new arguments
  f.args[1].args = new_arguments

  inputs_list, param_list
end

function inject_inputs!(f::Expr, inputlist)
  #reverse it, because inject! prepends
  for input in reverse(inputlist)
    inject!(f, :(Verilog.@input $(input[1]) $(input[2])))
  end
end

function inject_verilog_generator!(f::Expr, params)
  fn_symbol = f.args[1].args[1]
  inject!(f, :(Verilog.@verigen($fn_symbol, $(params...))))
end

#this should be passed the function block.
function linebyline_adaptor!(block::Expr, input_list = nothing)
  newargs = []
  for argument in block.args
    if isa(argument, Symbol)
      identifier = argument
      push!(newargs, :(Verilog.@final $identifier))
    elseif (argument.head == :tuple)
      #check if we're trying to return a tuple of wire symbols.
      identifiers = argument.args
      push!(newargs, :(Verilog.@final $(identifiers...)))
    elseif (input_list != nothing) &&
        (argument.head == :macrocall)  &&
        (argument.args[1] == Symbol("@input"))

      #search through the input list for a corresponding input statement and set
      #the parameter value.
      for idx = 1:length(input_list)
        if (input_list[idx][1] == argument.args[2])
          input_list[idx] = (argument.args[2], argument.args[3])
        end
      end
    elseif (argument.head == :(=))
      #all assignments should be surrounded by check to update the list.
      identifier = argument.args[1]
      assignment = argument.args[2]
      push!(newargs, :(Verilog.@assign $identifier $assignment))
    elseif (argument.head == :for)
      #println("forloop block:")
      linebyline_adaptor!(argument.args[2])
      push!(newargs, argument)
    elseif (argument.head == :if)
      linebyline_adaptor!(argument.args[2])
      if length(argument.args) == 3
        linebyline_adaptor!(argument.args[3])
      end
      push!(newargs, argument)
    elseif (argument.head == :return)
      if isa(argument.args[1], Symbol)
        #check to see if we're trying to return a symbol.
        identifier = argument.args[1]
        push!(newargs, :(Verilog.@final $identifier))
      elseif isa(argument.args[1], Expr) && (argument.args[1].head == :tuple)
        #check if we're trying to return a tuple of wire symbols.
        identifiers = argument.args[1].args
        push!(newargs, :(Verilog.@final $(identifiers...)))
      else
        throw(ArgumentError("@verilog functions may only return wires or tuples of wires."))
      end
      push!(newargs, :(@goto fin))
    else
      push!(newargs, argument)
    end
  end
  #replace the old arguments with new arguments.
  block.args = newargs
  nothing
end

function inject_verilog_finisher!(f::Expr)
  push!(f.args[2].args, :(@label fin))
  push!(f.args[2].args, :(Verilog.@verifin))
end

function module_transform!(f_module::Expr, input_list, params_list)
  #do a line-by-line analysis and readapt the contents of f_module
  linebyline_adaptor!(f_module.args[2], input_list)

  #add the last line, which will finish the verilog analysis.
  inject_verilog_finisher!(f_module)

  #inject input macros.
  inject_inputs!(f_module, input_list)

  #next, inject the script generator

  inject_verilog_generator!(f_module, params_list)
end

function module_translate(f::Expr)
  #for now, do nothing.
  f_module = copy(f)
  #check if there's an existing output definition and send it to nothing.
  strip_output_definition!(f_module)
  #strip any inputs that happen to be wires
  input_list, params_list = strip_wire_inputs!(f_module)

  (length(input_list) > 0) || throw(ArgumentError("selected function does not appear to be a verilog module."))

  module_transform!(f_module, input_list, params_list)

  set_output_type!(f_module, :String)

  inject!(f_module, :(Verilog.@verimode :verilog))

  return f_module
end

################################################################################
## TEXT rewrite.
## what happens when you submit a "text rewrite" to the

function substitute_wire_inputs_as_wiretext!(f::Expr)
  #issue a list of places where we need extra wire macros.
  inputs_list = Vector{Tuple{Symbol, Any}}()
  param_list = Symbol[]

  for argument in f.args[1].args[2:end]
    if isa(argument, Symbol)
      push!(param_list, argument)
      continue
    end

    #a bare wire statement.
    if argument.args[2] == :Wire
      #reset it to Unsigned.
      argument.args[2] = :(Verilog.WireObject)
      #save the needed symbol and wire range parameter.
      push!(inputs_list, (argument.args[1], nothing))
    elseif argument.args[2] == :SingleWire
      argument.args[2] = :(Verilog.WireObject{0:0v})

      #save the needed symbol and wire range parameter.
      push!(inputs_list, (argument.args[1], :(0:0v)))
    elseif isa(argument.args[2], Expr) &&
      (argument.args[2].head == :curly) &&
      (argument.args[2].args[1] == :Wire) &&
      (argument.args[2].args[2].head == :(:))

      colondef = argument.args[2].args[2]

      #then reset it to WireObject
      argument.args[2] = :(Verilog.WireObject{$colondef})

      #save the needed symbol and wire range parameter.
      push!(inputs_list, (argument.args[1], argument.args[2].args[2]))
    else
      push!(param_list, argument.args[1])
    end
  end

  inputs_list, param_list
end

function wiretext_translate(f::Expr)
  f_wiretext = copy(f)

  input_list,param_list = substitute_wire_inputs_as_wiretext!(f_wiretext)

  module_transform!(f_wiretext, input_list, param_list)

  set_output_type!(f_wiretext, :(Verilog.ModuleObject))

  input_identifiers = [t[1] for t in input_list]

  inject!(f_wiretext, :(Verilog.@verimode :modulecall $(input_identifiers...)))

  f_wiretext
end


text:
# Boostingthecomputationofthematrixexponential

export graph_sid_exp;

"""
    (graph,cref)=graph_sid_exp(k;T=Float64)

Computes a polynomial evaluation approximating the exponential using `k` matrix
multiplications following the procedure in the reference. The coefficients are
directly copied from the paper.

The evaluation is embedded in the degopt-format, see [`graph_degopt`](@ref),
using the function [`graph_sastre_yks_degopt`](@ref). Moreover, for `k<=3` it
uses the Paterson–Stockmeyer method.

**Reference**

1. J. Sastre, J. Ibáñez, E. Defez. "Boosting the computation of the matrix
   exponential". Applied Mathematics of Computation, 340:206-220, 2019.
   DOI: [10.1016/j.amc.2018.08.017](https://doi.org/10.1016/j.amc.2018.08.017)
"""
function graph_sid_exp(k; T = Float64)
    if (k <= 3)
        if k == 3
            deg = 6
        elseif k == 2
            deg = 4
        elseif k == 1
            deg = 2
        else
            deg = 1
        end
        return graph_ps_degopt(convert.(T, 1 ./ factorial.(0:deg)))

    elseif (k == 4)
        c1  = 4.018761610201036 * 1e-4
        c2  = 2.945531440279683 * 1e-3
        c3  = -8.709066576837676 * 1e-3
        c4  = 4.017568440673568 * 1e-1
        c5  = 3.230762888122312 * 1e-2
        c6  = 5.768988513026145 * 1e0
        c7  = 2.338576034271299 * 1e-2
        c8  = 2.381070373870987 * 1e-1
        c9  = 2.224209172496374 * 1e0
        c10 = -5.792361707073261 * 1e0
        c11 = -4.130276365929783 * 1e-2
        c12 = 1.040801735231354 * 1e1
        c13 = -6.331712455883370 * 1e1
        c14 = 3.484665863364574 * 1e-1
        c15 = 1
        c16 = 1

        c = [
            [[c2, c1], [0.0, 0, 0]], # (10)
            [[1.0], [0, c4, c3], [1.0], [0, 0, c5], [c6], [0.0, 0, c7]], # (11)
            [
                [0.0, 1],
                [0, c9, c8],
                [c10, 1],
                [0, c11, 0],
                [c13, c12],
                [c16, c15, c14],
            ], # (12)
        ]

        c = convert.(Vector{Vector{T}}, c)

        (graph, cref) = graph_sastre_yks_degopt(2, 2, c)

    elseif (k == 5)
        c1  = 1.161658834444880e-06
        c2  = 4.500852739573010e-06
        c3  = 5.374708803114821e-05
        c4  = 2.005403977292901e-03
        c5  = 6.974348269544424e-02
        c6  = 9.418613214806352e-01
        c7  = 2.852960512714315e-03
        c8  = -7.544837153586671e-03
        c9  = 1.829773504500424e+00
        c10 = 3.151382711608315e-02
        c11 = 1.392249143769798e-01
        c12 = -2.269101241269351e-03
        c13 = -5.394098846866402e-02
        c14 = 3.112216227982407e-01
        c15 = 9.343851261938047e+00
        c16 = 6.865706355662834e-01
        c17 = 3.233370163085380e+00
        c18 = -5.726379787260966e+00
        c19 = -1.413550099309667e-02
        c20 = -1.638413114712016e-01

        c = [
            [[c3, c2, c1], [0.0, 0, 0, 0]], # (17)
            [
                [1.0],
                [0, c6, c5, c4],
                [1.0],
                [0, 0, c8, c7],
                [c9],
                [0, 0, c11, c10],
            ], # (18)
            [
                [0.0, 1],
                [0, c14, c13, c12],
                [c15, 1],
                [0, c16, 0, 0],
                [c18, c17],
                [1, 1, c20, c19],
            ], # (19)
        ]

        c = convert.(Vector{Vector{T}}, c)

        (graph, cref) = graph_sastre_yks_degopt(2, 3, c)

    elseif (k == 6)
        # Coefficients from http://personales.upv.es/~jorsasma/software/expmpol.m
        c = [
            1.172460202011541e-08
            9.379681616092325e-08
            1.406952242413849e-06
            2.294895435403922e-05
            2.024281516007681e-03
            1.430688980356062e-02
            1.952545843107103e-01
            2.865001388641538e+00
            -1.204349003694297e-03
            2.547056607231984e-03
            2.721930992200371e-02
            2.498969092549990e+02
            2.018492049443954e-02
            1.965098904519709e-01
            1.739158441630994e+00
            8.290085751394409e+00
            2.919349464582001e-04
            1.758035313846159e-04
            1.606091400855144e-02
            3.655234395347475e-02
            2.243394407902074e-03
            -3.005000525808178e-02
            1.969779342112314e-01
            1
            1
        ]

        # Evaluation from  http://personales.upv.es/~jorsasma/software/expmpol.m
        # y0s=Ap{4}*(c(1)*Ap{4}+c(2)*Ap{3}+c(3)*Ap{2}+c(4)*A);
        # y1s=(y0s+c(5)*Ap{4}+c(6)*Ap{3}+c(7)*Ap{2}+c(8)*A)*(y0s+c(9)*Ap{4}+c(10)*Ap{3}+c(11)*Ap{2})+c(12)*y0s+c(13)*Ap{4}+c(14)*Ap{3}+c(15)*Ap{2}+c(16)*A;
        # sol=y1s*(y0s+c(17)*Ap{4}+c(18)*Ap{3}+c(19)*Ap{2}+c(20)*A)+c(21)*Ap{4}+c(22)*Ap{3}+c(23)*Ap{2}+A+eye(n);
        CC = [
            [[c[4], c[3], c[2], c[1]], [0.0, 0, 0, 0, 0]],
            [
                [1.0],
                [0, c[8], c[7], c[6], c[5]],
                [1.0],
                [0, 0, c[11], c[10], c[9]],
                [c[12]],
                [0, c[16], c[15], c[14], c[13]],
            ],
            [
                [0.0, 1],
                [0.0, 0, 0, 0, 0],
                [1.0, 0],
                [0, c[20], c[19], c[18], c[17]],
                [0.0, 0],
                [1, 1, c[23], c[22], c[21]],
            ], # (19)
        ]

        CC = convert.(Vector{Vector{T}}, CC)

        (graph, cref) = graph_sastre_yks_degopt(2, 4, CC)

    elseif (k == 7)
        # Coefficients from http://personales.upv.es/~jorsasma/software/expmpol.m
        c = [
            1.556371639324141e-11
            1.556371639324141e-10
            2.957106114715868e-09
            6.204734935438909e-08
            1.313681421698863e-06
            3.501669195497238e-05
            1.283057135586989e-03
            2.479095151834799e-02
            4.155284057336423e-01
            5.951585263506065e+00
            3.753710741641900e-05
            2.100333647757715e-04
            2.630043177655382e-03
            3.306559506631931e-02
            6.175954247606858e+01
            2.742336655922557e-03
            3.005135891320298e-02
            2.857950268422422e-01
            2.991654767354374e+00
            1.110689398085882e+01
            8.572383602707347e-06
            9.027588625491207e-05
            1.121744731945438e-03
            8.139086096860678e-03
            -2.638236222337760e-04
            6.263526066651383e-05
            4.985549176118462e-03
            7.705596948494946e-02
            5.029302610017967e-01
            1
            1
        ]

        # Evaluation from  http://personales.upv.es/~jorsasma/software/expmpol.m
        # y0s=Ap{5}*(c(1)*Ap{5}+c(2)*Ap{4}+c(3)*Ap{3}+c(4)*Ap{2}+c(5)*A);
        # y1s=(y0s+c(6)*Ap{5}+c(7)*Ap{4}+c(8)*Ap{3}+c(9)*Ap{2}+c(10)*A)*(y0s+c(11)*Ap{5}+c(12)*Ap{4}+c(13)*Ap{3}+c(14)*Ap{2})+c(15)*y0s+c(16)*Ap{5}+c(17)*Ap{4}+c(18)*Ap{3}+c(19)*Ap{2}+c(20)*A;
        # sol=y1s*(y0s+c(21)*Ap{5}+c(22)*Ap{4}+c(23)*Ap{3}+c(24)*Ap{2}+c(25)*A)+c(26)*Ap{5}+c(27)*Ap{4}+c(28)*Ap{3}+c(29)*Ap{2}+A+eye(n);
        CC = [
            [[c[5], c[4], c[3], c[2], c[1]], [0.0, 0, 0, 0, 0, 0]],
            [
                [1.0],
                [0, c[10], c[9], c[8], c[7], c[6]],
                [1.0],
                [0, 0, c[14], c[13], c[12], c[11]],
                [c[15]],
                [0, c[20], c[19], c[18], c[17], c[16]],
            ],
            [
                [0.0, 1],
                [0.0, 0, 0, 0, 0, 0],
                [1.0, 0],
                [0, c[25], c[24], c[23], c[22], c[21]],
                [0.0, 0],
                [1, 1, c[29], c[28], c[27], c[26]],
            ],
        ]

        CC = convert.(Vector{Vector{T}}, CC)

        (graph, cref) = graph_sastre_yks_degopt(2, 5, CC)

    elseif (k > 7)
        (graph, cref) = graph_sid_exp(7; T = T)
        # Square it 7-k times
        for j = 8:k
            degopt = Degopt(graph)
            square!(degopt)
            scale!(degopt, convert(T, 1 / 2))
            (graph, cref) = graph_degopt(degopt)
        end
    else
        error("Not implemented for k=$k")
    end

    return (graph, cref)
end


text:
"""
    gmtinfo(cmd0::String="", arg1=nothing; kwargs...)

Reads files and finds the extreme values in each of the columns.

Full option list at [`gmtinfo`]($(GMTdoc)gmtinfo.html)

Parameters
----------

- **A** | **ranges** :: [Type => Str]

    Specify how the range should be reported.
    ($(GMTdoc)gmtinfo.html#a)
- **C** | **per_column** :: [Type => Bool]

    Report the min/max values per column in separate columns [Default uses <min/max> format].
    ($(GMTdoc)gmtinfo.html#c)
- **D** | **center** :: [Type => Bool]

    Modifies results obtained by -I by shifting the region to better align with the center of the data.
    ($(GMTdoc)gmtinfo.html#d)
- **E** | **get_record** :: [Type => Str | []]

    Returns the record whose column col contains the minimum (l) or maximum (h) value.
    ($(GMTdoc)gmtinfo.html#e)
- **F** | **counts** :: [Type => Str | []]

    Returns the counts of various records depending on the appended mode.
    ($(GMTdoc)gmtinfo.html#f)
- **I** | **inc** :: [Type => Str | Number]

    Report the min/max of the first n columns to the nearest multiple of the provided increments
    and output results in the form -Rw/e/s/n
    ($(GMTdoc)gmtinfo.html#i)
- **L** | **common_limits** :: [Type => Bool]

    Determines common limits across tables or segments.
    ($(GMTdoc)gmtinfo.html#l)
- **S** | **for_error_bars** :: [Type => Str | []]

    Add extra space for error bars. Useful together with I option and when later plotting with `plot E`.
    ($(GMTdoc)gmtinfo.html#s)
- **T** | **nearest_multiple** :: [Type => Str | Number]    ``Arg = dz[+ccol]``

    Report the min/max of the first (0’th) column to the nearest multiple of dz and output this as
    the string -Tzmin/zmax/dz.
    ($(GMTdoc)gmtinfo.html#t)
- $(GMT.opt_V)
- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_o)
- $(GMT.opt_r)
- $(GMT.opt_swap_xy)
"""
function gmtinfo(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("gmtinfo", cmd0, arg1)

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode

	cmd, = parse_common_opts(d, "", [:V_params :e :f :o :r :yx])
	(endswith(cmd, "-:")) && (cmd *= "i")    # Need to be -:i not -: to not swap output too
	cmd = parse_these_opts(cmd, d, [[:A :ranges], [:C :per_column], [:D :center], [:E :get_record], [:F :counts],
	                                [:L :common_limits], [:S :for_error_bars]])
	cmd = add_opt(d, cmd, 'I', [:I :inc],
	              (exact=("e", nothing, 1), polyg=("b", nothing, 1), surface=("s", nothing, 1), fft=("d", nothing, 1), inc=("", arg2str, 2)), false, true)
	cmd = add_opt(d, cmd, 'T', [:T :nearest_multiple], (dz="", col="+c", column="+c"))

	# If file name sent in, read it.
	cmd, arg1, = read_data(d, cmd0, cmd, arg1, " ")
	if (dbg_print_cmd(d, cmd) !== nothing)  return cmd  end
	isa(arg1, Tuple) ? gmt("gmtinfo " * cmd, arg1...) : gmt("gmtinfo " * cmd, arg1)
end

# ---------------------------------------------------------------------------------------------------
gmtinfo(arg1, cmd0::String=""; kw...) = gmtinfo(cmd0, arg1; kw...)