text:
#Partitions

export
    integer_partitions,
    ncpartitions,
    partitions,
    prevprod
    #nextprod,

#integer partitions

struct IntegerPartitions
    n::Int
end

Base.length(p::IntegerPartitions) = npartitions(p.n)

function Base.iterate(p::IntegerPartitions, xs = Int[])
    length(xs) == p.n && return
    xs = nextpartition(p.n,xs)
    (xs, xs)
end

"""
    partitions(n)

Generate all integer arrays that sum to `n`. Because the number of partitions can be very
large, this function returns an iterator object. Use `collect(partitions(n))` to get an
array of all partitions. The number of partitions to generate can be efficiently computed
using `length(partitions(n))`.
"""
partitions(n::Integer) = IntegerPartitions(n)

function nextpartition(n, as)
    isempty(as) && return Int[n]

    xs = similar(as, 0)
    sizehint!(xs, length(as) + 1)

    for i = 1:length(as)-1
        if as[i+1] == 1
            x = as[i]-1
            push!(xs, x)
            n -= x
            while n > x
                push!(xs, x)
                n -= x
            end
            push!(xs, n)

            return xs
        end
        push!(xs, as[i])
        n -= as[i]
    end
    push!(xs, as[end]-1)
    push!(xs, 1)

    xs
end

let _npartitions = Dict{Int,Int}()
    global npartitions
    function npartitions(n::Int)
        if n < 0
            0
        elseif n < 2
            1
        elseif (np = get(_npartitions, n, 0)) > 0
            np
        else
            np = 0
            sgn = 1
            for k = 1:n
                np += sgn * (npartitions(n - (k*(3k-1)) >> 1) + npartitions(n - (k*(3k+1)) >> 1))
                sgn = -sgn
            end
            _npartitions[n] = np
        end
    end
end

# Algorithm H from TAoCP 7.2.1.4
# Partition n into m parts
# in colex order (lexicographic by reflected sequence)

struct FixedPartitions
    n::Int
    m::Int
end

Base.length(f::FixedPartitions) = npartitions(f.n,f.m)

"""
    partitions(n, m)

Generate all arrays of `m` integers that sum to `n`. Because the number of partitions can
be very large, this function returns an iterator object. Use `collect(partitions(n, m))` to
get an array of all partitions. The number of partitions to generate can be efficiently
computed using `length(partitions(n, m))`.
"""
partitions(n::Integer, m::Integer) =
    n >= 1 && m >= 1 ?
        FixedPartitions(n, m) :
        throw(DomainError((n, m), "n and m must be positive"))

function Base.iterate(f::FixedPartitions, s::Vector{Int} = Int[])
    f.m <= f.n || return
    if !isempty(s)
        (f.m == 1 || s[1]-1 <= s[end]) && return
    end

    xs = nextfixedpartition(f.n,f.m,s)
    (xs, xs)
end

function nextfixedpartition(n, m, bs)
    as = copy(bs)
    if isempty(as)
        # First iteration
        as = [n-m+1; ones(Int, m-1)]
    elseif as[2] < as[1]-1
        # Most common iteration
        as[1] -= 1
        as[2] += 1
    else
        # Iterate
        j = 0
        s = as[1]+as[2]-1
        for jj = 3:m # TODO: use `for outer j = ...` on 0.7
            j = jj
            as[jj] < as[1]-1 && break
            s += as[jj]
        end
        x = as[j] += 1
        for k = j-1:-1:2
            as[k] = x
            s -= x
        end
        as[1] = s
    end

    return as
end

let _nipartitions = Dict{Tuple{Int,Int},Int}()
    global npartitions
    function npartitions(n::Int, m::Int)
        if n < m || m == 0
            0
        elseif n == m
            1
        elseif (np = get(_nipartitions, (n,m), 0)) > 0
            np
        else
            _nipartitions[(n, m)] = npartitions(n-1, m-1) + npartitions(n-m, m)
        end
    end
end

# Algorithm H from TAoCP 7.2.1.5
# Set partitions

struct SetPartitions{T<:AbstractVector}
    s::T
end

Base.length(p::SetPartitions) = nsetpartitions(length(p.s))

"""
    partitions(s::AbstractVector)

Generate all set partitions of the elements of an array `s`, represented as arrays of
arrays. Because the number of partitions can be very large, this function returns an
iterator object. Use `collect(partitions(s))` to get an array of all partitions. The
number of partitions to generate can be efficiently computed using
`length(partitions(s))`.
"""
partitions(s::AbstractVector) = SetPartitions(s)

function Base.iterate(p::SetPartitions)
    n = length(p.s)
    iterate(p, (zeros(Int32, n), ones(Int32, n-1), n, 1))
end
function Base.iterate(p::SetPartitions, s)
    s[1][1] > 0 && return
    nextsetpartition(p.s, s...)
end

function nextsetpartition(s::AbstractVector, a, b, n, m)
    function makeparts(s, a, m)
        temp = [similar(s, 0) for k = 0:m]
        for i = 1:n
            push!(temp[a[i]+1], s[i])
        end
        filter!(!isempty, temp)
    end

    isempty(s) && return ([s], ([1], Int[], n, 1))

    part = makeparts(s,a,m)

    if a[end] != m
        a[end] += 1
    else
        j = 0
        for jj = n-1:-1:1
            j = jj
            a[jj] == b[jj] || break
        end
        a[j] += 1
        m = b[j] + (a[j] == b[j])
        for k = j+1:n-1
            a[k] = 0
            b[k] = m
        end
        a[end] = 0
    end

    return (part, (a,b,n,m))
end

let _nsetpartitions = Dict{Int,Int}()
    global nsetpartitions
    function nsetpartitions(n::Int)
        if n < 0
            0
        elseif n < 2
            1
        elseif (wn = get(_nsetpartitions, n, 0)) > 0
            wn
        else
            wn = 0
            for k = 0:n-1
                wn += binomial(n-1,k)*nsetpartitions(n-1-k)
            end
            _nsetpartitions[n] = wn
        end
    end
end

struct FixedSetPartitions{T<:AbstractVector}
    s::T
    m::Int
end

Base.length(p::FixedSetPartitions) = nfixedsetpartitions(length(p.s),p.m)

"""
    partitions(s::AbstractVector, m::Int)

Generate all set partitions of the elements of an array `s` into exactly `m` subsets,
represented as arrays of arrays. Because the number of partitions can be very large,
this function returns an iterator object. Use `collect(partitions(s, m))` to get
an array of all partitions. The number of partitions into `m` subsets is equal to the
Stirling number of the second kind, and can be efficiently computed using
`length(partitions(s, m))`.
"""
partitions(s::AbstractVector, m::Int) =
    length(s) >= 1 && m >= 1 ?
        FixedSetPartitions(s, m) :
        throw(DomainError((length(s), m), "length(s) and m must be positive"))

function Base.iterate(p::FixedSetPartitions)
    n = length(p.s)
    m = p.m
    state = m <= n ? (vcat(ones(Int, n-m),1:m), vcat(1,n-m+2:n), n) : (Int[], Int[], n)
    # state consists of:
    # vector a of length n describing to which partition every element of s belongs
    # vector b of length n describing the first index b[i] that belongs to partition i
    # integer n

    iterate(p, state)
end

function Base.iterate(p::FixedSetPartitions, s)
    (isempty(s[1]) || s[1][1] > 1) && return
    nextfixedsetpartition(p.s,p.m, s...)
end

function nextfixedsetpartition(s::AbstractVector, m, a, b, n)
    function makeparts(s, a)
        part = [similar(s, 0) for k = 1:m]
        for i = 1:n
            push!(part[a[i]], s[i])
        end
        return part
    end

    part = makeparts(s,a)

    if m == 1
        a[1] = 2
        return (part, (a, b, n))
    end

    if a[end] != m
        a[end] += 1
    else
        j = k = 0
        for jj = n-1:-1:1
            j = jj
            if a[j] < m && b[a[j]+1] < j
                break
            end
        end
        if j > 1
            a[j] += 1
            for p = j+1:n
                if b[a[p]] != p
                    a[p] = 1
                end
            end
        else
            for kk = m:-1:2
                k = kk
                if b[k-1] < b[k] - 1
                    break
                end
            end
            b[k] -= 1
            b[k+1:m] = n-m+k+1:n
            a[1:n] .= 1
            a[b] = 1:m
        end
    end

    return (part, (a,b,n))
end

function nfixedsetpartitions(n::Int, m::Int)
    numpart = 0
    for k = 0:m
        numpart += (-1)^(m-k) * binomial(m, k) * (k^n)
    end
    numpart = div(numpart, factorial(m))
    return numpart
end

# TODO: Base.DSP is no longer a thing in Julia 0.7
#This function is still defined in Base because it is being used by Base.DSP
#"""
#Next integer not less than `n` that can be written as $\prod k_i^{p_i}$ for integers $p_1$, $p_2$, etc.
#
#For a list of integers i1, i2, i3, find the smallest
#    i1^n1 * i2^n2 * i3^n3 >= x
#for integer n1, n2, n3
#"""
#function nextprod(a::Vector{Int}, x)
#    if x > typemax(Int)
#        throw(ArgumentError("unsafe for x > typemax(Int), got $x"))
#    end
#    k = length(a)
#    v = ones(Int, k)                  # current value of each counter
#    mx = [nextpow(ai,x) for ai in a]  # maximum value of each counter
#    v[1] = mx[1]                      # start at first case that is >= x
#    p::widen(Int) = mx[1]             # initial value of product in this case
#    best = p
#    icarry = 1
#
#    while v[end] < mx[end]
#        if p >= x
#            best = p < best ? p : best  # keep the best found yet
#            carrytest = true
#            while carrytest
#                p = div(p, v[icarry])
#                v[icarry] = 1
#                icarry += 1
#                p *= a[icarry]
#                v[icarry] *= a[icarry]
#                carrytest = v[icarry] > mx[icarry] && icarry < k
#            end
#            if p < x
#                icarry = 1
#            end
#        else
#            while p < x
#                p *= a[1]
#                v[1] *= a[1]
#            end
#        end
#    end
#    best = mx[end] < best ? mx[end] : best
#    return Int(best)  # could overflow, but best to have predictable return type
#end

"""
    prevprod(a::Vector{Int}, x)

Previous integer not greater than `x` that can be written as ``\\prod k_i^{p_i}`` for
integers ``p_1``, ``p_2``, etc.

For integers ``i_1``, ``i_2``, ``i_3``, this is equivalent to finding the largest ``x``
such that
```math
i_1^n_1 i_2^n_2 i_3^n_3 \\leq x
```
for integers ``n_1``, ``n_2``, ``n_3``.
"""
function prevprod(a::Vector{Int}, x)
    if x > typemax(Int)
        throw(ArgumentError("unsafe for x > typemax(Int), got $x"))
    end
    k = length(a)
    v = ones(Int, k)                  # current value of each counter
    mx = [nextpow(ai,x) for ai in a]  # allow each counter to exceed p (sentinel)
    first = Int(prevpow(a[1], x))     # start at best case in first factor
    v[1] = first
    p::widen(Int) = first
    best = p
    icarry = 1

    while v[end] < mx[end]
        while p <= x
            best = p > best ? p : best
            p *= a[1]
            v[1] *= a[1]
        end
        if p > x
            carrytest = true
            while carrytest
                p = div(p, v[icarry])
                v[icarry] = 1
                icarry += 1
                p *= a[icarry]
                v[icarry] *= a[icarry]
                carrytest = v[icarry] > mx[icarry] && icarry < k
            end
            if p <= x
                icarry = 1
            end
        end
    end
    best = x >= p > best ? p : best
    return Int(best)
end

"""
    integer_partitions(n)

List the partitions of the integer `n`.

!!! note
    The order of the resulting array is consistent with that produced by the computational
    discrete algebra software GAP.
"""
function integer_partitions(n::Integer)
    if n < 0
        throw(DomainError(n, "n must be nonnegative"))
    elseif n == 0
        return Vector{Int}[]
    elseif n == 1
        return Vector{Int}[[1]]
    end

    list = Vector{Int}[]

    for p in integer_partitions(n-1)
        push!(list, [p; 1])
        if length(p) == 1 || p[end] < p[end-1]
            push!(list, [p[1:end-1]; p[end]+1])
        end
    end

    list
end

#Noncrossing partitions

if VERSION >= v"0.7.0-DEV.3251"
    const _cmp = cmp
else
    const _cmp = lexcmp
end

#Produces noncrossing partitions of length n
function ncpartitions(n::Int)
    partitions = Vector{Vector{Int}}[]
    _ncpart!(1,n,n,Vector{Int}[], partitions)
    partitions
end

function _ncpart!(a::Int, b::Int, nn::Int, x::Vector, partitions::Vector)
    n = b - a + 1
    for k = 1:n, root in CoolLexCombinations(n, k)
        root .+= a - 1
        #Abort if construction is out of lex order
        !isempty(x) && _cmp(x[end], root) == 1 && return

        #Save if we've filled all the holes
        sofar = Vector{Int}[x..., root]
        ssofaru = sort!(union(sofar...))
        if length(ssofaru) == nn && ssofaru == collect(1:nn)
            push!(partitions, sofar)
            return
        end

        #otherwise patch all remaining holes
        blob = [ssofaru; nn + 1]
        for l = 1:length(blob)-1
            ap, bp = blob[l] + 1, blob[l+1] - 1
            ap <= bp && _ncpart!(ap, bp, nn, sofar, partitions)
        end
    end
end


text:
# Term Relations
# http://www.geneontology.org/page/ontology-relations

# The following relations are defined in this library:
# * is a
# * part of
# * has part
# * regulates
# * negatively regulates
# * positively regulates

ext: jl


text:
module LinearMaps

export LinearMap

using LinearAlgebra
using SparseArrays

abstract type LinearMap{T} end

Base.eltype(::LinearMap{T}) where {T} = T

Base.isreal(A::LinearMap) = eltype(A) <: Real
LinearAlgebra.issymmetric(::LinearMap) = false # default assumptions
LinearAlgebra.ishermitian(A::LinearMap{<:Real}) = issymmetric(A)
LinearAlgebra.ishermitian(::LinearMap) = false # default assumptions
LinearAlgebra.isposdef(::LinearMap) = false # default assumptions

Base.ndims(::LinearMap) = 2
Base.size(A::LinearMap, n) = (n==1 || n==2 ? size(A)[n] : error("LinearMap objects have only 2 dimensions"))
Base.length(A::LinearMap) = size(A)[1] * size(A)[2]

Base.:(*)(A::LinearMap, x::AbstractVector) = mul!(similar(x, promote_type(eltype(A), eltype(x)), size(A, 1)), A, x)
function LinearAlgebra.mul!(y::AbstractVector, A::LinearMap{T}, x::AbstractVector, α::Number=one(T), β::Number=zero(T)) where {T}
    length(y) == size(A, 1) || throw(DimensionMismatch("mul!"))
    if α == 1
        β == 0 && (A_mul_B!(y, A, x); return y)
        β == 1 && (y .+= A * x; return y)
        # β != 0, 1
        rmul!(y, β)
        y .+= A * x
        return y
    elseif α == 0
        β == 0 && (fill!(y, zero(eltype(y))); return y)
        β == 1 && return y
        # β != 0, 1
        rmul!(y, β)
        return y
    else # α != 0, 1
        β == 0 && (A_mul_B!(y, A, x); rmul!(y, α); return y)
        β == 1 && (y .+= rmul!(A * x, α); return y)
        # β != 0, 1
        rmul!(y, β)
        y .+= rmul!(A * x, α)
        return y
    end
end
# the following is of interest in, e.g., subspace-iteration methods
function LinearAlgebra.mul!(Y::AbstractMatrix, A::LinearMap{T}, X::AbstractMatrix, α::Number=one(T), β::Number=zero(T)) where {T}
    (size(Y, 1) == size(A, 1) && size(X, 1) == size(A, 2) && size(Y, 2) == size(X, 2)) || throw(DimensionMismatch("mul!"))
    @inbounds @views for i = 1:size(X, 2)
        mul!(Y[:, i], A, X[:, i], α, β)
    end
    # starting from Julia v1.1, we could use the `eachcol` iterator
    # for (Xi, Yi) in zip(eachcol(X), eachcol(Y))
    #     mul!(Yi, A, Xi, α, β)
    # end
    return Y
end

A_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector)  = mul!(y, A, x)
At_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, transpose(A), x)
Ac_mul_B!(y::AbstractVector, A::AbstractMatrix, x::AbstractVector) = mul!(y, adjoint(A), x)

# Matrix: create matrix representation of LinearMap
function Base.Matrix(A::LinearMap)
    M, N = size(A)
    T = eltype(A)
    mat = Matrix{T}(undef, (M, N))
    v = fill(zero(T), N)
    @inbounds for i = 1:N
        v[i] = one(T)
        mul!(view(mat, :, i), A, v)
        v[i] = zero(T)
    end
    return mat
end

Base.Array(A::LinearMap) = Matrix(A)
Base.convert(::Type{Matrix}, A::LinearMap) = Matrix(A)
Base.convert(::Type{Array}, A::LinearMap) = Matrix(A)
Base.convert(::Type{SparseMatrixCSC}, A::LinearMap) = sparse(A)

# sparse: create sparse matrix representation of LinearMap
function SparseArrays.sparse(A::LinearMap{T}) where {T}
    M, N = size(A)
    rowind = Int[]
    nzval = T[]
    colptr = Vector{Int}(undef, N+1)
    v = fill(zero(T), N)
    Av = Vector{T}(undef, M)

    for i = 1:N
        v[i] = one(T)
        mul!(Av, A, v)
        js = findall(!iszero, Av)
        colptr[i] = length(nzval) + 1
        if length(js) > 0
            append!(rowind, js)
            append!(nzval, Av[js])
        end
        v[i] = zero(T)
    end
    colptr[N+1] = length(nzval) + 1

    return SparseMatrixCSC(M, N, colptr, rowind, nzval)
end

include("transpose.jl") # transposing linear maps
include("linearcombination.jl") # defining linear combinations of linear maps
include("composition.jl") # composition of linear maps
include("wrappedmap.jl") # wrap a matrix of linear map in a new type, thereby allowing to alter its properties
include("uniformscalingmap.jl") # the uniform scaling map, to be able to make linear combinations of LinearMap objects and multiples of I
include("functionmap.jl") # using a function as linear map
include("blockmap.jl") # block linear maps

"""
    LinearMap(A; kwargs...)
    LinearMap{T=Float64}(f, [fc,], M::Int, N::Int = M; kwargs...)

Construct a linear map object, either from an existing `LinearMap` or `AbstractMatrix` `A`,
with the purpose of redefining its properties via the keyword arguments `kwargs`, or
from a function or callable object `f`. In the latter case, one also needs to specify
the size of the equivalent matrix representation `(M, N)`, i.e. for functions `f` acting
on length `N` vectors and producing length `M` vectors (with default value `N=M`). Preferably,
also the `eltype` `T` of the corresponding matrix representation needs to be specified, i.e.
whether the action of `f` on a vector will be similar to e.g. multiplying by numbers of type `T`.
If not specified, the devault value `T=Float64` will be assumed. Optionally, a corresponding
function `fc` can be specified that implements the transpose/adjoint of `f`.

The keyword arguments and their default values for functions `f` are
*   issymmetric::Bool = false : whether `A` or `f` acts as a symmetric matrix
*   ishermitian::Bool = issymmetric & T<:Real : whether `A` or `f` acts as a Hermitian matrix
*   isposdef::Bool = false : whether `A` or `f` acts as a positive definite matrix.
For existing linear maps or matrices `A`, the default values will be taken by calling
`issymmetric`, `ishermitian` and `isposdef` on the existing object `A`.

For functions `f`, there is one more keyword arguments
*   ismutating::Bool : flags whether the function acts as a mutating matrix multiplication
    `f(y,x)` where the result vector `y` is the first argument (in case of `true`),
    or as a normal matrix multiplication that is called as `y=f(x)` (in case of `false`).
    The default value is guessed by looking at the number of arguments of the first occurence
    of `f` in the method table.
"""
LinearMap(A::Union{AbstractMatrix, LinearMap}; kwargs...) = WrappedMap(A; kwargs...)
LinearMap(f, M::Int; kwargs...) = LinearMap{Float64}(f, M; kwargs...)
LinearMap(f, M::Int, N::Int; kwargs...) = LinearMap{Float64}(f, M, N; kwargs...)
LinearMap(f, fc, M::Int; kwargs...) = LinearMap{Float64}(f, fc, M; kwargs...)
LinearMap(f, fc, M::Int, N::Int; kwargs...) = LinearMap{Float64}(f, fc, M, N; kwargs...)

LinearMap{T}(A::Union{AbstractMatrix, LinearMap}; kwargs...) where {T} = WrappedMap{T}(A; kwargs...)
LinearMap{T}(f, args...; kwargs...) where {T} = FunctionMap{T}(f, args...; kwargs...)

end # module


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

using Test, Printf

@testset "Printf" begin

@testset "%p" begin

    # pointers
    if Sys.WORD_SIZE == 64
        @test (Printf.@sprintf "%20p" 0) == "  0x0000000000000000"
        @test (Printf.@sprintf "%-20p" 0) == "0x0000000000000000  "
        @test (Printf.@sprintf "%20p" C_NULL) == "  0x0000000000000000"
        @test (@sprintf "%-20p" C_NULL) == "0x0000000000000000  "
    elseif Sys.WORD_SIZE == 32
        @test (Printf.@sprintf "%20p" 0) == "          0x00000000"
        @test (Printf.@sprintf "%-20p" 0) == "0x00000000          "
        @test (@sprintf "%20p" C_NULL) == "          0x00000000"
        @test (@sprintf "%-20p" C_NULL) == "0x00000000          "
    end

end

@testset "%a" begin

    # hex float
    @test (Printf.@sprintf "%a" 1.5) == "0x1.8p+0"
    @test (Printf.@sprintf "%a" 1.5f0) == "0x1.8p+0"
    @test (Printf.@sprintf "%a" big"1.5") == "0x1.8p+0"
    @test (Printf.@sprintf "%#.0a" 1.5) == "0x2.p+0"
    @test (Printf.@sprintf "%+30a" 1/3) == "         +0x1.5555555555555p-2"

    @test Printf.@sprintf("%a", 1.5) == "0x1.8p+0"
    @test Printf.@sprintf("%a", 3.14) == "0x1.91eb851eb851fp+1"
    @test Printf.@sprintf("%.0a", 3.14) == "0x2p+1"
    @test Printf.@sprintf("%.1a", 3.14) == "0x1.9p+1"
    @test Printf.@sprintf("%.2a", 3.14) == "0x1.92p+1"
    @test Printf.@sprintf("%#a", 3.14) == "0x1.91eb851eb851fp+1"
    @test Printf.@sprintf("%#.0a", 3.14) == "0x2.p+1"
    @test Printf.@sprintf("%#.1a", 3.14) == "0x1.9p+1"
    @test Printf.@sprintf("%#.2a", 3.14) == "0x1.92p+1"
    @test Printf.@sprintf("%.6a", 1.5) == "0x1.800000p+0"

end

@testset "%g" begin

    # %g
    for (val, res) in ((12345678., "1.23457e+07"),
                    (1234567.8, "1.23457e+06"),
                    (123456.78, "123457"),
                    (12345.678, "12345.7"),
                    (12340000.0, "1.234e+07"))
        @test (Printf.@sprintf("%.6g", val) == res)
    end
    for (val, res) in ((big"12345678.", "1.23457e+07"),
                    (big"1234567.8", "1.23457e+06"),
                    (big"123456.78", "123457"),
                    (big"12345.678", "12345.7"))
        @test (Printf.@sprintf("%.6g", val) == res)
    end
    for (fmt, val) in (("%10.5g", "     123.4"),
                    ("%+10.5g", "    +123.4"),
                    ("% 10.5g","     123.4"),
                    ("%#10.5g", "    123.40"),
                    ("%-10.5g", "123.4     "),
                    ("%-+10.5g", "+123.4    "),
                    ("%010.5g", "00000123.4")),
        num in (123.4, big"123.4")
        @test Printf.format(Printf.Format(fmt), num) == val
    end
    @test( Printf.@sprintf( "%10.5g", -123.4 ) == "    -123.4")
    @test( Printf.@sprintf( "%010.5g", -123.4 ) == "-0000123.4")
    @test( Printf.@sprintf( "%.6g", 12340000.0 ) == "1.234e+07")
    @test( Printf.@sprintf( "%#.6g", 12340000.0 ) == "1.23400e+07")
    @test( Printf.@sprintf( "%10.5g", big"-123.4" ) == "    -123.4")
    @test( Printf.@sprintf( "%010.5g", big"-123.4" ) == "-0000123.4")
    @test( Printf.@sprintf( "%.6g", big"12340000.0" ) == "1.234e+07")
    @test( Printf.@sprintf( "%#.6g", big"12340000.0") == "1.23400e+07")

    # %g regression gh #14331
    @test( Printf.@sprintf( "%.5g", 42) == "42")
    @test( Printf.@sprintf( "%#.2g", 42) == "42.")
    @test( Printf.@sprintf( "%#.5g", 42) == "42.000")

    @test Printf.@sprintf("%g", 0.00012) == "0.00012"
    @test Printf.@sprintf("%g", 0.000012) == "1.2e-05"
    @test Printf.@sprintf("%g", 123456.7) == "123457"
    @test Printf.@sprintf("%g", 1234567.8) == "1.23457e+06"

    # zeros
    @test Printf.@sprintf("%.15g", 0) == "0"
    @test Printf.@sprintf("%#.15g", 0) == "0.00000000000000"

end

@testset "%f" begin

    # Inf / NaN handling
    @test (Printf.@sprintf "%f" Inf) == "Inf"
    @test (Printf.@sprintf "%+f" Inf) == "+Inf"
    @test (Printf.@sprintf "% f" Inf) == " Inf"
    @test (Printf.@sprintf "% #f" Inf) == " Inf"
    @test (Printf.@sprintf "%f" -Inf) == "-Inf"
    @test (Printf.@sprintf "%+f" -Inf) == "-Inf"
    @test (Printf.@sprintf "%f" NaN) == "NaN"
    @test (Printf.@sprintf "%+f" NaN) == "NaN"
    @test (Printf.@sprintf "% f" NaN) == "NaN"
    @test (Printf.@sprintf "% #f" NaN) == "NaN"
    @test (Printf.@sprintf "%e" big"Inf") == "Inf"
    @test (Printf.@sprintf "%e" big"NaN") == "NaN"

    @test (Printf.@sprintf "%.0f" 3e142) == "29999999999999997463140672961703247153805615792184250659629251954072073858354858644285983761764971823910371920726635399393477049701891710124032"

    @test Printf.@sprintf("%f", 1.234) == "1.234000"
    @test Printf.@sprintf("%F", 1.234) == "1.234000"
    @test Printf.@sprintf("%+f", 1.234) == "+1.234000"
    @test Printf.@sprintf("% f", 1.234) == " 1.234000"
    @test Printf.@sprintf("%f", -1.234) == "-1.234000"
    @test Printf.@sprintf("%+f", -1.234) == "-1.234000"
    @test Printf.@sprintf("% f", -1.234) == "-1.234000"
    @test Printf.@sprintf("%#f", 1.234) == "1.234000"
    @test Printf.@sprintf("%.2f", 1.234) == "1.23"
    @test Printf.@sprintf("%.2f", 1.235) == "1.24"
    @test Printf.@sprintf("%.2f", 0.235) == "0.23"
    @test Printf.@sprintf("%4.1f", 1.234) == " 1.2"
    @test Printf.@sprintf("%8.1f", 1.234) == "     1.2"
    @test Printf.@sprintf("%+8.1f", 1.234) == "    +1.2"
    @test Printf.@sprintf("% 8.1f", 1.234) == "     1.2"
    @test Printf.@sprintf("% 7.1f", 1.234) == "    1.2"
    @test Printf.@sprintf("% 08.1f", 1.234) == " 00001.2"
    @test Printf.@sprintf("%08.1f", 1.234) == "000001.2"
    @test Printf.@sprintf("%-08.1f", 1.234) == "1.2     "
    @test Printf.@sprintf("%-8.1f", 1.234) == "1.2     "
    @test Printf.@sprintf("%08.1f", -1.234) == "-00001.2"
    @test Printf.@sprintf("%09.1f", -1.234) == "-000001.2"
    @test Printf.@sprintf("%09.1f", 1.234) == "0000001.2"
    @test Printf.@sprintf("%+09.1f", 1.234) == "+000001.2"
    @test Printf.@sprintf("% 09.1f", 1.234) == " 000001.2"
    @test Printf.@sprintf("%+ 09.1f", 1.234) == "+000001.2"
    @test Printf.@sprintf("%+ 09.1f", 1.234) == "+000001.2"
    @test Printf.@sprintf("%+ 09.0f", 1.234) == "+00000001"
    @test Printf.@sprintf("%+ #09.0f", 1.234) == "+0000001."
end

@testset "%e" begin

    # Inf / NaN handling
    @test (Printf.@sprintf "%e" Inf) == "Inf"
    @test (Printf.@sprintf "%+e" Inf) == "+Inf"
    @test (Printf.@sprintf "% e" Inf) == " Inf"
    @test (Printf.@sprintf "% #e" Inf) == " Inf"
    @test (Printf.@sprintf "%e" -Inf) == "-Inf"
    @test (Printf.@sprintf "%+e" -Inf) == "-Inf"
    @test (Printf.@sprintf "%e" NaN) == "NaN"
    @test (Printf.@sprintf "%+e" NaN) == "NaN"
    @test (Printf.@sprintf "% e" NaN) == "NaN"
    @test (Printf.@sprintf "% #e" NaN) == "NaN"
    @test (Printf.@sprintf "%e" big"Inf") == "Inf"
    @test (Printf.@sprintf "%e" big"NaN") == "NaN"

    # scientific notation
    @test (Printf.@sprintf "%.0e" 3e142) == "3e+142"
    @test (Printf.@sprintf "%#.0e" 3e142) == "3.e+142"
    @test (Printf.@sprintf "%.0e" big"3e142") == "3e+142"
    @test (Printf.@sprintf "%#.0e" big"3e142") == "3.e+142"

    @test (Printf.@sprintf "%.0e" big"3e1042") == "3e+1042"

    @test (Printf.@sprintf "%e" 3e42) == "3.000000e+42"
    @test (Printf.@sprintf "%E" 3e42) == "3.000000E+42"
    @test (Printf.@sprintf "%e" 3e-42) == "3.000000e-42"
    @test (Printf.@sprintf "%E" 3e-42) == "3.000000E-42"

    @test Printf.@sprintf("%e", 1.234) == "1.234000e+00"
    @test Printf.@sprintf("%E", 1.234) == "1.234000E+00"
    @test Printf.@sprintf("%+e", 1.234) == "+1.234000e+00"
    @test Printf.@sprintf("% e", 1.234) == " 1.234000e+00"
    @test Printf.@sprintf("%e", -1.234) == "-1.234000e+00"
    @test Printf.@sprintf("%+e", -1.234) == "-1.234000e+00"
    @test Printf.@sprintf("% e", -1.234) == "-1.234000e+00"
    @test Printf.@sprintf("%#e", 1.234) == "1.234000e+00"
    @test Printf.@sprintf("%.2e", 1.234) == "1.23e+00"
    @test Printf.@sprintf("%.2e", 1.235) == "1.24e+00"
    @test Printf.@sprintf("%.2e", 0.235) == "2.35e-01"
    @test Printf.@sprintf("%4.1e", 1.234) == "1.2e+00"
    @test Printf.@sprintf("%8.1e", 1.234) == " 1.2e+00"
    @test Printf.@sprintf("%+8.1e", 1.234) == "+1.2e+00"
    @test Printf.@sprintf("% 8.1e", 1.234) == " 1.2e+00"
    @test Printf.@sprintf("% 7.1e", 1.234) == " 1.2e+00"
    @test Printf.@sprintf("% 08.1e", 1.234) == " 1.2e+00"
    @test Printf.@sprintf("%08.1e", 1.234) == "01.2e+00"
    @test Printf.@sprintf("%-08.1e", 1.234) == "1.2e+00 "
    @test Printf.@sprintf("%-8.1e", 1.234) == "1.2e+00 "
    @test Printf.@sprintf("%-8.1e", 1.234) == "1.2e+00 "
    @test Printf.@sprintf("%08.1e", -1.234) == "-1.2e+00"
    @test Printf.@sprintf("%09.1e", -1.234) == "-01.2e+00"
    @test Printf.@sprintf("%09.1e", 1.234) == "001.2e+00"
    @test Printf.@sprintf("%+09.1e", 1.234) == "+01.2e+00"
    @test Printf.@sprintf("% 09.1e", 1.234) == " 01.2e+00"
    @test Printf.@sprintf("%+ 09.1e", 1.234) == "+01.2e+00"
    @test Printf.@sprintf("%+ 09.1e", 1.234) == "+01.2e+00"
    @test Printf.@sprintf("%+ 09.0e", 1.234) == "+0001e+00"
    @test Printf.@sprintf("%+ #09.0e", 1.234) == "+001.e+00"
end

@testset "strings" begin

    @test Printf.@sprintf("Hallo heimur") == "Hallo heimur"
    @test Printf.@sprintf("+%s+", "hello") == "+hello+"
    @test Printf.@sprintf("%.1s", "foo") == "f"
    @test Printf.@sprintf("%s", "%%%%") == "%%%%"
    @test Printf.@sprintf("%s", "Hallo heimur") == "Hallo heimur"
    @test Printf.@sprintf("%+s", "Hallo heimur") == "Hallo heimur"
    @test Printf.@sprintf("% s", "Hallo heimur") == "Hallo heimur"
    @test Printf.@sprintf("%+ s", "Hallo heimur") == "Hallo heimur"
    @test Printf.@sprintf("%1s", "Hallo heimur") == "Hallo heimur"
    @test Printf.@sprintf("%20s", "Hallo") == "               Hallo"
    @test Printf.@sprintf("%-20s", "Hallo") == "Hallo               "
    @test Printf.@sprintf("%0-20s", "Hallo") == "Hallo               "
    @test Printf.@sprintf("%.20s", "Hallo heimur") == "Hallo heimur"
    @test Printf.@sprintf("%20.5s", "Hallo heimur") == "               Hallo"
    @test Printf.@sprintf("%.0s", "Hallo heimur") == ""
    @test Printf.@sprintf("%20.0s", "Hallo heimur") == "                    "
    @test Printf.@sprintf("%.s", "Hallo heimur") == ""
    @test Printf.@sprintf("%20.s", "Hallo heimur") == "                    "
    @test (Printf.@sprintf "%s" "test") == "test"
    @test (Printf.@sprintf "%s" "tést") == "tést"
    @test Printf.@sprintf("ø%sø", "hey") == "øheyø"
    @test Printf.@sprintf("%4sø", "ø") == "   øø"
    @test Printf.@sprintf("%-4sø", "ø") == "ø   ø"

    @test (Printf.@sprintf "%8s" "test") == "    test"
    @test (Printf.@sprintf "%-8s" "test") == "test    "

    @test (Printf.@sprintf "%s" :test) == "test"
    @test (Printf.@sprintf "%#s" :test) == ":test"
    @test (Printf.@sprintf "%#8s" :test) == "   :test"
    @test (Printf.@sprintf "%#-8s" :test) == ":test   "

    @test (Printf.@sprintf "%8.3s" "test") == "     tes"
    @test (Printf.@sprintf "%#8.3s" "test") == "     \"te"
    @test (Printf.@sprintf "%-8.3s" "test") == "tes     "
    @test (Printf.@sprintf "%#-8.3s" "test") == "\"te     "
    @test (Printf.@sprintf "%.3s" "test") == "tes"
    @test (Printf.@sprintf "%#.3s" "test") == "\"te"
    @test (Printf.@sprintf "%-.3s" "test") == "tes"
    @test (Printf.@sprintf "%#-.3s" "test") == "\"te"

end

@testset "chars" begin

    @test Printf.@sprintf("%c", 'a') == "a"
    @test Printf.@sprintf("%c",  32) == " "
    @test Printf.@sprintf("%c",  36) == "\$"
    @test Printf.@sprintf("%3c", 'a') == "  a"
    @test Printf.@sprintf( "%c", 'x') == "x"
    @test Printf.@sprintf("%+c", 'x') == "x"
    @test Printf.@sprintf("% c", 'x') == "x"
    @test Printf.@sprintf("%+ c", 'x') == "x"
    @test Printf.@sprintf("%1c", 'x') == "x"
    @test Printf.@sprintf("%20c"  , 'x') == "                   x"
    @test Printf.@sprintf("%-20c" , 'x') == "x                   "
    @test Printf.@sprintf("%-020c", 'x') == "x                   "
    @test Printf.@sprintf("%c", 65) == "A"
    @test Printf.@sprintf("%c", 'A') == "A"
    @test Printf.@sprintf("%3c", 'A') == "  A"
    @test Printf.@sprintf("%-3c", 'A') == "A  "
    @test Printf.@sprintf("%c", 248) == "ø"
    @test Printf.@sprintf("%c", 'ø') == "ø"
    @test Printf.@sprintf("%c", "ø") == "ø"
    @test Printf.@sprintf("%c", '𐀀') == "𐀀"

end

function _test_flags(val, vflag::AbstractString, fmt::AbstractString, res::AbstractString, prefix::AbstractString)
    vflag = string("%", vflag)
    space_fmt = string(length(res) + length(prefix) + 3, fmt)
    fsign = string((val < 0 ? "-" : "+"), prefix)
    nsign = string((val < 0 ? "-" : " "), prefix)
    osign = val < 0 ? string("-", prefix) : string(prefix, "0")
    esign = string(val < 0 ? "-" : "", prefix)
    esignend = val < 0 ? "" : " "

    for (flag::AbstractString, ans::AbstractString) in (
            ("", string("  ", nsign, res)),
            ("+", string("  ", fsign, res)),
            (" ", string("  ", nsign, res)),
            ("0", string(osign, "00", res)),
            ("-", string(esign, res, "  ", esignend)),
            ("0+", string(fsign, "00", res)),
            ("0 ", string(nsign, "00", res)),
            ("-+", string(fsign, res, "  ")),
            ("- ", string(nsign, res, "  ")),
        )
        fmt_string = string(vflag, flag, space_fmt)
        fmtd = Printf.format(Printf.Format(fmt_string), val)
        @test fmtd == ans
    end
end

@testset "basics" begin

    @test Printf.@sprintf("%%") == "%"
    @test Printf.@sprintf("hey there") == "hey there"
    @test_throws ArgumentError Printf.Format("")
    @test_throws ArgumentError Printf.Format("%+")
    @test_throws ArgumentError Printf.Format("%.")
    @test_throws ArgumentError Printf.Format("%.0")
    @test isempty(Printf.Format("%%").formats)
    @test Printf.@sprintf("%d%d", 1, 2) == "12"
    @test (Printf.@sprintf "%d%d" [1 2]...) == "12"
    @test (Printf.@sprintf("X%d", 2)) == "X2"
    @test (Printf.@sprintf("\u00d0%d", 2)) == "\u00d02"
    @test (Printf.@sprintf("\u0f00%d", 2)) == "\u0f002"
    @test (Printf.@sprintf("\U0001ffff%d", 2)) == "\U0001ffff2"
    @test (Printf.@sprintf("%dX%d", 1, 2)) == "1X2"
    @test (Printf.@sprintf("%d\u00d0%d", 1, 2)) == "1\u00d02"
    @test (Printf.@sprintf("%d\u0f00%d", 1, 2)) == "1\u0f002"
    @test (Printf.@sprintf("%d\U0001ffff%d", 1, 2)) == "1\U0001ffff2"
    @test (Printf.@sprintf("%d\u2203%d\u0203", 1, 2)) == "1\u22032\u0203"
    @test_throws ArgumentError Printf.Format("%y%d")
    @test_throws ArgumentError Printf.Format("%\u00d0%d")
    @test_throws ArgumentError Printf.Format("%\u0f00%d")
    @test_throws ArgumentError Printf.Format("%\U0001ffff%d")
    @test Printf.@sprintf("%10.5d", 4) == "     00004"
    @test (Printf.@sprintf "%d" typemax(Int64)) == "9223372036854775807"

    for (fmt, val) in (("%7.2f", "   1.23"),
                   ("%-7.2f", "1.23   "),
                   ("%07.2f", "0001.23"),
                   ("%.0f", "1"),
                   ("%#.0f", "1."),
                   ("%.4e", "1.2345e+00"),
                   ("%.4E", "1.2345E+00"),
                   ("%.2a", "0x1.3cp+0"),
                   ("%.2A", "0X1.3CP+0")),
        num in (1.2345, big"1.2345")
        @test Printf.format(Printf.Format(fmt), num) == val
    end

    for (fmt, val) in (("%i", "42"),
                   ("%u", "42"),
                   ("Test: %i", "Test: 42"),
                   ("%#x", "0x2a"),
                   ("%x", "2a"),
                   ("%X", "2A"),
                   ("% i", " 42"),
                   ("%+i", "+42"),
                   ("%4i", "  42"),
                   ("%-4i", "42  "),
                   ("%f", "42.000000"),
                   ("%g", "42"),
                   ("%e", "4.200000e+01")),
        num in (UInt16(42), UInt32(42), UInt64(42), UInt128(42),
                Int16(42), Int32(42), Int64(42), Int128(42), big"42")
        @test Printf.format(Printf.Format(fmt), num) == val
    end

    for i in (
            (42, "", "i", "42", ""),
            (42, "", "d", "42", ""),

            (42, "", "u", "42", ""),
            (42, "", "x", "2a", ""),
            (42, "", "X", "2A", ""),
            (42, "", "o", "52", ""),

            (42, "#", "x", "2a", "0x"),
            (42, "#", "X", "2A", "0X"),
            (42, "#", "o", "052", ""),

            (1.2345, "", ".2f", "1.23", ""),
            (1.2345, "", ".2e", "1.23e+00", ""),
            (1.2345, "", ".2E", "1.23E+00", ""),

            (1.2345, "#", ".0f", "1.", ""),
            (1.2345, "#", ".0e", "1.e+00", ""),
            (1.2345, "#", ".0E", "1.E+00", ""),

            (1.2345, "", ".2a", "1.3cp+0", "0x"),
            (1.2345, "", ".2A", "1.3CP+0", "0X"),
        )
        _test_flags(i...)
        _test_flags(-i[1], i[2:5]...)
    end

    # reasonably complex
    @test (Printf.@sprintf "Test: %s%c%C%c%#-.0f." "t" 65 66 67 -42) == "Test: tABC-42.."

    # combo
    @test (Printf.@sprintf "%f %d %d %f" 1.0 [3 4]... 5) == "1.000000 3 4 5.000000"

    # multi
    @test (Printf.@sprintf "%s %f %9.5f %d %d %d %d%d%d%d" [1:6;]... [7,8,9,10]...) == "1 2.000000   3.00000 4 5 6 78910"

    # comprehension
    @test (Printf.@sprintf "%s %s %s %d %d %d %f %f %f" Any[10^x+y for x=1:3,y=1:3 ]...) == "11 101 1001 12 102 1002 13.000000 103.000000 1003.000000"

    # more than 16 formats/args
    @test (Printf.@sprintf "%s %s %s %d %d %d %f %f %f %s %s %s %d %d %d %f %f %f" Any[10*x+(x+1) for x=1:18 ]...) ==
        "12 23 34 45 56 67 78.000000 89.000000 100.000000 111 122 133 144 155 166 177.000000 188.000000 199.000000"

    # Check bug with trailing nul printing BigFloat
    @test (Printf.@sprintf("%.330f", BigFloat(1)))[end] != '\0'

    # issue #29662
    @test (Printf.@sprintf "%12.3e" pi*1e100) == "  3.142e+100"

    @test string(Printf.Format("%a").formats[1]) == "%a"
    @test string(Printf.Format("%a").formats[1]; modifier="R") == "%Ra"

    @test Printf.@sprintf("%d", 3.14) == "3"
    @test Printf.@sprintf("%2d", 3.14) == " 3"
    @test Printf.@sprintf("%2d", big(3.14)) == " 3"
    @test Printf.@sprintf("%s", 1) == "1"
    @test Printf.@sprintf("%f", 1) == "1.000000"
    @test Printf.@sprintf("%e", 1) == "1.000000e+00"
    @test Printf.@sprintf("%g", 1) == "1"

    # escaped '%'
    @test_throws ArgumentError @sprintf("%s%%%s", "a")
    @test @sprintf("%s%%%s", "a", "b") == "a%%b"

    # print float as %d uses round(x)
    @test @sprintf("%d", 25.5) == "26"

    # 37539
    @test @sprintf(" %.1e\n", 0.999) == " 1.0e+00\n"
    @test @sprintf("   %.1f", 9.999) == "   10.0"

    # 37552
    @test @sprintf("%d", 1.0e100) == "10000000000000000159028911097599180468360808563945281389781327557747838772170381060813469985856815104"
    @test @sprintf("%d", 3//1) == "3"
end

@testset "integers" begin

    @test Printf.@sprintf( "% d",  42) == " 42"
    @test Printf.@sprintf( "% d", -42) == "-42"
    @test Printf.@sprintf( "% 5d",  42) == "   42"
    @test Printf.@sprintf( "% 5d", -42) == "  -42"
    @test Printf.@sprintf( "% 15d",  42) == "             42"
    @test Printf.@sprintf( "% 15d", -42) == "            -42"
    @test Printf.@sprintf("%+d",  42) == "+42"
    @test Printf.@sprintf("%+d", -42) == "-42"
    @test Printf.@sprintf("%+5d",  42) == "  +42"
    @test Printf.@sprintf("%+5d", -42) == "  -42"
    @test Printf.@sprintf("%+15d",  42) == "            +42"
    @test Printf.@sprintf("%+15d", -42) == "            -42"
    @test Printf.@sprintf( "%0d",  42) == "42"
    @test Printf.@sprintf( "%0d", -42) == "-42"
    @test Printf.@sprintf( "%05d",  42) == "00042"
    @test Printf.@sprintf( "%05d", -42) == "-0042"
    @test Printf.@sprintf( "%015d",  42) == "000000000000042"
    @test Printf.@sprintf( "%015d", -42) == "-00000000000042"
    @test Printf.@sprintf("%-d",  42) == "42"
    @test Printf.@sprintf("%-d", -42) == "-42"
    @test Printf.@sprintf("%-5d",  42) == "42   "
    @test Printf.@sprintf("%-5d", -42) == "-42  "
    @test Printf.@sprintf("%-15d",  42) == "42             "
    @test Printf.@sprintf("%-15d", -42) == "-42            "
    @test Printf.@sprintf("%-0d",  42) == "42"
    @test Printf.@sprintf("%-0d", -42) == "-42"
    @test Printf.@sprintf("%-05d",  42) == "42   "
    @test Printf.@sprintf("%-05d", -42) == "-42  "
    @test Printf.@sprintf("%-015d",  42) == "42             "
    @test Printf.@sprintf("%-015d", -42) == "-42            "
    @test Printf.@sprintf( "%0-d",  42) == "42"
    @test Printf.@sprintf( "%0-d", -42) == "-42"
    @test Printf.@sprintf( "%0-5d",  42) == "42   "
    @test Printf.@sprintf( "%0-5d", -42) == "-42  "
    @test Printf.@sprintf( "%0-15d",  42) == "42             "
    @test Printf.@sprintf( "%0-15d", -42) == "-42            "
    @test_throws ArgumentError Printf.Format("%d %")

    @test Printf.@sprintf("%lld", 18446744065119617025) == "18446744065119617025"
    @test Printf.@sprintf("%+8lld", 100) == "    +100"
    @test Printf.@sprintf("%+.8lld", 100) == "+00000100"
    @test Printf.@sprintf("%+10.8lld", 100) == " +00000100"
    @test_throws ArgumentError Printf.Format("%_1lld")
    @test Printf.@sprintf("%-1.5lld", -100) == "-00100"
    @test Printf.@sprintf("%5lld", 100) == "  100"
    @test Printf.@sprintf("%5lld", -100) == " -100"
    @test Printf.@sprintf("%-5lld", 100) == "100  "
    @test Printf.@sprintf("%-5lld", -100) == "-100 "
    @test Printf.@sprintf("%-.5lld", 100) == "00100"
    @test Printf.@sprintf("%-.5lld", -100) == "-00100"
    @test Printf.@sprintf("%-8.5lld", 100) == "00100   "
    @test Printf.@sprintf("%-8.5lld", -100) == "-00100  "
    @test Printf.@sprintf("%05lld", 100) == "00100"
    @test Printf.@sprintf("%05lld", -100) == "-0100"
    @test Printf.@sprintf("% lld", 100) == " 100"
    @test Printf.@sprintf("% lld", -100) == "-100"
    @test Printf.@sprintf("% 5lld", 100) == "  100"
    @test Printf.@sprintf("% 5lld", -100) == " -100"
    @test Printf.@sprintf("% .5lld", 100) == " 00100"
    @test Printf.@sprintf("% .5lld", -100) == "-00100"
    @test Printf.@sprintf("% 8.5lld", 100) == "   00100"
    @test Printf.@sprintf("% 8.5lld", -100) == "  -00100"
    @test Printf.@sprintf("%.0lld", 0) == "0"
    @test Printf.@sprintf("%#+21.18llx", -100) == "-0x000000000000000064"
    @test Printf.@sprintf("%#.25llo", -100) == "-00000000000000000000000144"
    @test Printf.@sprintf("%#+24.20llo", -100) == "  -000000000000000000144"
    @test Printf.@sprintf("%#+18.21llX", -100) == "-0X000000000000000000064"
    @test Printf.@sprintf("%#+20.24llo", -100) == "-0000000000000000000000144"
    @test Printf.@sprintf("%#+25.22llu", -1) == "  -0000000000000000000001"
    @test Printf.@sprintf("%#+25.22llu", -1) == "  -0000000000000000000001"
    @test Printf.@sprintf("%#+30.25llu", -1) == "    -0000000000000000000000001"
    @test Printf.@sprintf("%+#25.22lld", -1) == "  -0000000000000000000001"
    @test Printf.@sprintf("%#-8.5llo", 100) == "000144  "
    @test Printf.@sprintf("%#-+ 08.5lld", 100) == "+00100  "
    @test Printf.@sprintf("%#-+ 08.5lld", 100) == "+00100  "
    @test Printf.@sprintf("%.40lld",  1) == "0000000000000000000000000000000000000001"
    @test Printf.@sprintf("% .40lld",  1) == " 0000000000000000000000000000000000000001"
    @test Printf.@sprintf("% .40d",  1) == " 0000000000000000000000000000000000000001"
    @test Printf.@sprintf("%lld",  18446744065119617025) == "18446744065119617025"

    @test Printf.@sprintf("+%d+",  10) == "+10+"
    @test Printf.@sprintf("%#012x",  1) == "0x0000000001"
    @test Printf.@sprintf("%#04.8x",  1) == "0x00000001"

    @test Printf.@sprintf("%#-08.2x",  1) == "0x01    "
    @test Printf.@sprintf("%#08o",  1) == "00000001"
    @test Printf.@sprintf("%d",  1024) == "1024"
    @test Printf.@sprintf("%d", -1024) == "-1024"
    @test Printf.@sprintf("%i",  1024) == "1024"
    @test Printf.@sprintf("%i", -1024) == "-1024"
    @test Printf.@sprintf("%u",  1024) == "1024"
    @test Printf.@sprintf("%u",  UInt(4294966272)) == "4294966272"
    @test Printf.@sprintf("%o",  511) == "777"
    @test Printf.@sprintf("%o",  UInt(4294966785)) == "37777777001"
    @test Printf.@sprintf("%x",  305441741) == "1234abcd"
    @test Printf.@sprintf("%x",  UInt(3989525555)) == "edcb5433"
    @test Printf.@sprintf("%X",  305441741) == "1234ABCD"
    @test Printf.@sprintf("%X",  UInt(3989525555)) == "EDCB5433"
    @test Printf.@sprintf("%+d",  1024) == "+1024"
    @test Printf.@sprintf("%+d", -1024) == "-1024"
    @test Printf.@sprintf("%+i",  1024) == "+1024"
    @test Printf.@sprintf("%+i", -1024) == "-1024"
    @test Printf.@sprintf("%+u",  1024) == "+1024"
    @test Printf.@sprintf("%+u",  UInt(4294966272)) == "+4294966272"
    @test Printf.@sprintf("%+o",  511) == "+777"
    @test Printf.@sprintf("%+o",  UInt(4294966785)) == "+37777777001"
    @test Printf.@sprintf("%+x",  305441741) == "+1234abcd"
    @test Printf.@sprintf("%+x",  UInt(3989525555)) == "+edcb5433"
    @test Printf.@sprintf("%+X",  305441741) == "+1234ABCD"
    @test Printf.@sprintf("%+X",  UInt(3989525555)) == "+EDCB5433"
    @test Printf.@sprintf("% d",  1024) == " 1024"
    @test Printf.@sprintf("% d", -1024) == "-1024"
    @test Printf.@sprintf("% i",  1024) == " 1024"
    @test Printf.@sprintf("% i", -1024) == "-1024"
    @test Printf.@sprintf("% u",  1024) == " 1024"
    @test Printf.@sprintf("% u",  UInt(4294966272)) == " 4294966272"
    @test Printf.@sprintf("% o",  511) == " 777"
    @test Printf.@sprintf("% o",  UInt(4294966785)) == " 37777777001"
    @test Printf.@sprintf("% x",  305441741) == " 1234abcd"
    @test Printf.@sprintf("% x",  UInt(3989525555)) == " edcb5433"
    @test Printf.@sprintf("% X",  305441741) == " 1234ABCD"
    @test Printf.@sprintf("% X",  UInt(3989525555)) == " EDCB5433"
    @test Printf.@sprintf("%+ d",  1024) == "+1024"
    @test Printf.@sprintf("%+ d", -1024) == "-1024"
    @test Printf.@sprintf("%+ i",  1024) == "+1024"
    @test Printf.@sprintf("%+ i", -1024) == "-1024"
    @test Printf.@sprintf("%+ u",  1024) == "+1024"
    @test Printf.@sprintf("%+ u",  UInt(4294966272)) == "+4294966272"
    @test Printf.@sprintf("%+ o",  511) == "+777"
    @test Printf.@sprintf("%+ o",  UInt(4294966785)) == "+37777777001"
    @test Printf.@sprintf("%+ x",  305441741) == "+1234abcd"
    @test Printf.@sprintf("%+ x",  UInt(3989525555)) == "+edcb5433"
    @test Printf.@sprintf("%+ X",  305441741) == "+1234ABCD"
    @test Printf.@sprintf("%+ X",  UInt(3989525555)) == "+EDCB5433"
    @test Printf.@sprintf("%#o",  511) == "0777"
    @test Printf.@sprintf("%#o",  UInt(4294966785)) == "037777777001"
    @test Printf.@sprintf("%#x",  305441741) == "0x1234abcd"
    @test Printf.@sprintf("%#x",  UInt(3989525555)) == "0xedcb5433"
    @test Printf.@sprintf("%#X",  305441741) == "0X1234ABCD"
    @test Printf.@sprintf("%#X",  UInt(3989525555)) == "0XEDCB5433"
    @test Printf.@sprintf("%#o",  UInt(0)) == "00"
    @test Printf.@sprintf("%#x",  UInt(0)) == "0x0"
    @test Printf.@sprintf("%#X",  UInt(0)) == "0X0"
    @test Printf.@sprintf("%1d",  1024) == "1024"
    @test Printf.@sprintf("%1d", -1024) == "-1024"
    @test Printf.@sprintf("%1i",  1024) == "1024"
    @test Printf.@sprintf("%1i", -1024) == "-1024"
    @test Printf.@sprintf("%1u",  1024) == "1024"
    @test Printf.@sprintf("%1u",  UInt(4294966272)) == "4294966272"
    @test Printf.@sprintf("%1o",  511) == "777"
    @test Printf.@sprintf("%1o",  UInt(4294966785)) == "37777777001"
    @test Printf.@sprintf("%1x",  305441741) == "1234abcd"
    @test Printf.@sprintf("%1x",  UInt(3989525555)) == "edcb5433"
    @test Printf.@sprintf("%1X",  305441741) == "1234ABCD"
    @test Printf.@sprintf("%1X",  UInt(3989525555)) == "EDCB5433"
    @test Printf.@sprintf("%20d",  1024) == "                1024"
    @test Printf.@sprintf("%20d", -1024) == "               -1024"
    @test Printf.@sprintf("%20i",  1024) == "                1024"
    @test Printf.@sprintf("%20i", -1024) == "               -1024"
    @test Printf.@sprintf("%20u",  1024) == "                1024"
    @test Printf.@sprintf("%20u",  UInt(4294966272)) == "          4294966272"
    @test Printf.@sprintf("%20o",  511) == "                 777"
    @test Printf.@sprintf("%20o",  UInt(4294966785)) == "         37777777001"
    @test Printf.@sprintf("%20x",  305441741) == "            1234abcd"
    @test Printf.@sprintf("%20x",  UInt(3989525555)) == "            edcb5433"
    @test Printf.@sprintf("%20X",  305441741) == "            1234ABCD"
    @test Printf.@sprintf("%20X",  UInt(3989525555)) == "            EDCB5433"
    @test Printf.@sprintf("%-20d",  1024) == "1024                "
    @test Printf.@sprintf("%-20d", -1024) == "-1024               "
    @test Printf.@sprintf("%-20i",  1024) == "1024                "
    @test Printf.@sprintf("%-20i", -1024) == "-1024               "
    @test Printf.@sprintf("%-20u",  1024) == "1024                "
    @test Printf.@sprintf("%-20u",  UInt(4294966272)) == "4294966272          "
    @test Printf.@sprintf("%-20o",  511) == "777                 "
    @test Printf.@sprintf("%-20o",  UInt(4294966785)) == "37777777001         "
    @test Printf.@sprintf("%-20x",  305441741) == "1234abcd            "
    @test Printf.@sprintf("%-20x",  UInt(3989525555)) == "edcb5433            "
    @test Printf.@sprintf("%-20X",  305441741) == "1234ABCD            "
    @test Printf.@sprintf("%-20X",  UInt(3989525555)) == "EDCB5433            "
    @test Printf.@sprintf("%020d",  1024) == "00000000000000001024"
    @test Printf.@sprintf("%020d", -1024) == "-0000000000000001024"
    @test Printf.@sprintf("%020i",  1024) == "00000000000000001024"
    @test Printf.@sprintf("%020i", -1024) == "-0000000000000001024"
    @test Printf.@sprintf("%020u",  1024) == "00000000000000001024"
    @test Printf.@sprintf("%020u",  UInt(4294966272)) == "00000000004294966272"
    @test Printf.@sprintf("%020o",  511) == "00000000000000000777"
    @test Printf.@sprintf("%020o",  UInt(4294966785)) == "00000000037777777001"
    @test Printf.@sprintf("%020x",  305441741) == "0000000000001234abcd"
    @test Printf.@sprintf("%020x",  UInt(3989525555)) == "000000000000edcb5433"
    @test Printf.@sprintf("%020X",  305441741) == "0000000000001234ABCD"
    @test Printf.@sprintf("%020X",  UInt(3989525555)) == "000000000000EDCB5433"
    @test Printf.@sprintf("%#20o",  511) == "                0777"
    @test Printf.@sprintf("%#20o",  UInt(4294966785)) == "        037777777001"
    @test Printf.@sprintf("%#20x",  305441741) == "          0x1234abcd"
    @test Printf.@sprintf("%#20x",  UInt(3989525555)) == "          0xedcb5433"
    @test Printf.@sprintf("%#20X",  305441741) == "          0X1234ABCD"
    @test Printf.@sprintf("%#20X",  UInt(3989525555)) == "          0XEDCB5433"
    @test Printf.@sprintf("%#020o",  511) == "00000000000000000777"
    @test Printf.@sprintf("%#020o",  UInt(4294966785)) == "00000000037777777001"
    @test Printf.@sprintf("%#020x",  305441741) == "0x00000000001234abcd"
    @test Printf.@sprintf("%#020x",  UInt(3989525555)) == "0x0000000000edcb5433"
    @test Printf.@sprintf("%#020X",  305441741) == "0X00000000001234ABCD"
    @test Printf.@sprintf("%#020X",  UInt(3989525555)) == "0X0000000000EDCB5433"
    @test Printf.@sprintf("%0-20d",  1024) == "1024                "
    @test Printf.@sprintf("%0-20d", -1024) == "-1024               "
    @test Printf.@sprintf("%0-20i",  1024) == "1024                "
    @test Printf.@sprintf("%0-20i", -1024) == "-1024               "
    @test Printf.@sprintf("%0-20u",  1024) == "1024                "
    @test Printf.@sprintf("%0-20u",  UInt(4294966272)) == "4294966272          "
    @test Printf.@sprintf("%-020o",  511) == "777                 "
    @test Printf.@sprintf("%-020o",  UInt(4294966785)) == "37777777001         "
    @test Printf.@sprintf("%-020x",  305441741) == "1234abcd            "
    @test Printf.@sprintf("%-020x",  UInt(3989525555)) == "edcb5433            "
    @test Printf.@sprintf("%-020X",  305441741) == "1234ABCD            "
    @test Printf.@sprintf("%-020X",  UInt(3989525555)) == "EDCB5433            "
    @test Printf.@sprintf("%.20d",  1024) == "00000000000000001024"
    @test Printf.@sprintf("%.20d", -1024) == "-00000000000000001024"
    @test Printf.@sprintf("%.20i",  1024) == "00000000000000001024"
    @test Printf.@sprintf("%.20i", -1024) == "-00000000000000001024"
    @test Printf.@sprintf("%.20u",  1024) == "00000000000000001024"
    @test Printf.@sprintf("%.20u",  UInt(4294966272)) == "00000000004294966272"
    @test Printf.@sprintf("%.20o",  511) == "00000000000000000777"
    @test Printf.@sprintf("%.20o",  UInt(4294966785)) == "00000000037777777001"
    @test Printf.@sprintf("%.20x",  305441741) == "0000000000001234abcd"
    @test Printf.@sprintf("%.20x",  UInt(3989525555)) == "000000000000edcb5433"
    @test Printf.@sprintf("%.20X",  305441741) == "0000000000001234ABCD"
    @test Printf.@sprintf("%.20X",  UInt(3989525555)) == "000000000000EDCB5433"
    @test Printf.@sprintf("%20.5d",  1024) == "               01024"
    @test Printf.@sprintf("%20.5d", -1024) == "              -01024"
    @test Printf.@sprintf("%20.5i",  1024) == "               01024"
    @test Printf.@sprintf("%20.5i", -1024) == "              -01024"
    @test Printf.@sprintf("%20.5u",  1024) == "               01024"
    @test Printf.@sprintf("%20.5u",  UInt(4294966272)) == "          4294966272"
    @test Printf.@sprintf("%20.5o",  511) == "               00777"
    @test Printf.@sprintf("%20.5o",  UInt(4294966785)) == "         37777777001"
    @test Printf.@sprintf("%20.5x",  305441741) == "            1234abcd"
    @test Printf.@sprintf("%20.10x",  UInt(3989525555)) == "          00edcb5433"
    @test Printf.@sprintf("%20.5X",  305441741) == "            1234ABCD"
    @test Printf.@sprintf("%20.10X",  UInt(3989525555)) == "          00EDCB5433"
    @test Printf.@sprintf("%020.5d",  1024) == "               01024"
    @test Printf.@sprintf("%020.5d", -1024) == "              -01024"
    @test Printf.@sprintf("%020.5i",  1024) == "               01024"
    @test Printf.@sprintf("%020.5i", -1024) == "              -01024"
    @test Printf.@sprintf("%020.5u",  1024) == "               01024"
    @test Printf.@sprintf("%020.5u",  UInt(4294966272)) == "          4294966272"
    @test Printf.@sprintf("%020.5o",  511) == "               00777"
    @test Printf.@sprintf("%020.5o",  UInt(4294966785)) == "         37777777001"
    @test Printf.@sprintf("%020.5x",  305441741) == "            1234abcd"
    @test Printf.@sprintf("%020.10x",  UInt(3989525555)) == "          00edcb5433"
    @test Printf.@sprintf("%020.5X",  305441741) == "            1234ABCD"
    @test Printf.@sprintf("%020.10X",  UInt(3989525555)) == "          00EDCB5433"
    @test Printf.@sprintf("%20.0d",  1024) == "                1024"
    @test Printf.@sprintf("%20.d", -1024) == "               -1024"
    @test Printf.@sprintf("%20.d",  0) == "                   0"
    @test Printf.@sprintf("%20.0i",  1024) == "                1024"
    @test Printf.@sprintf("%20.i", -1024) == "               -1024"
    @test Printf.@sprintf("%20.i",  0) == "                   0"
    @test Printf.@sprintf("%20.u",  1024) == "                1024"
    @test Printf.@sprintf("%20.0u",  UInt(4294966272)) == "          4294966272"
    @test Printf.@sprintf("%20.u",  UInt(0)) == "                   0"
    @test Printf.@sprintf("%20.o",  511) == "                 777"
    @test Printf.@sprintf("%20.0o",  UInt(4294966785)) == "         37777777001"
    @test Printf.@sprintf("%20.o",  UInt(0)) == "                   0"
    @test Printf.@sprintf("%20.x",  305441741) == "            1234abcd"
    @test Printf.@sprintf("%20.0x",  UInt(3989525555)) == "            edcb5433"
    @test Printf.@sprintf("%20.x",  UInt(0)) == "                   0"
    @test Printf.@sprintf("%20.X",  305441741) == "            1234ABCD"
    @test Printf.@sprintf("%20.0X",  UInt(3989525555)) == "            EDCB5433"
    @test Printf.@sprintf("%20.X",  UInt(0)) == "                   0"

end

end # @testset "Printf"

max_stars_repo_name: mlxd/julia

max_issues_repo_name: mlxd/julia

max_forks_repo_name: mlxd/julia


text:
function residual_violation(mechanism::Mechanism)
    violation = 0.0
    for joint in mechanism.joints
        res = constraint(mechanism, joint)
        shift = 0
        for element in [joint.translational, joint.rotational]
            Nλ = joint_length(element)
            Nb = limits_length(element)
            subres = res[shift + 2Nb .+ (1:Nλ)]
            violation = max(violation, norm(subres, Inf))
            shift += impulses_length(element)
        end
    end
    for body in mechanism.bodies
        res = constraint(mechanism, body)
        violation = max(violation, norm(res, Inf))
    end
    for contact in mechanism.contacts
        res = constraint(mechanism, contact)
        violation = max(violation, norm(res, Inf))
    end
    return violation
end

function bilinear_violation(mechanism::Mechanism)
    violation = 0.0
    for contact in mechanism.contacts
        comp = complementarity(mechanism, contact)
        violation = max(violation, norm(comp, Inf))
    end
    for joint in mechanism.joints
        comp = complementarity(mechanism, joint)
        violation = max(violation, norm(comp, Inf))
    end
    return violation
end


text:
function MCMCRun(Problem::MCMCSetup)

#Assume additive Guassian noise for each measurement
likelihood = (μ,σ) -> MvNormal(μ,σ*ones(length(μ)))
#Choose an MCMC sampler
sampler = Turing.NUTS(0.65) #0.65 is the desired acceptance ratio
#Define all of the parameters
getParNames = copy(Problem.MI.Parameters)
parTuple = Tuple{(push!(getParNames,:σ))...} #Tack on the std
# Define the default value when no data is passed
allData = [experiment[2] for experiment in Problem.Data]
defaults = (x = allData,)

modelFramework(varInfo,sampler,model) = begin
  #Set the accumulated log Posterior to zero
  varInfo.logp = 0.0
  #Set the default data for when no input is given
  x = model.defaults.x

  #Loop through the priors and assign them to the model
  θ = Vector(undef,length(Problem.Priors))

  for (prior,name,i) in zip(Problem.Priors,parNames,1:length(Problem.MI.Parameters))
    θ[i], logp = Turing.assume(sampler,prior,Turing.VarName([:mf,name],""),varInfo)
    varInfo.logp += logp
  end

  #Convert θ to array
  θ = convert(Array{typeof(first(θ))},θ) #Is this needed?

  #Assume normal random error in data
  likelihoodDistPriors = InverseGamma(2,3) #prior for data? Why this in particular?
  σ, logp = Turing.assume(sampler,likelihoodDistPriors,Turing.VarName([:mf,parNames[end]],""),varInfo)
  varInfo.logp += logp

  #Remake a new problem with the different parameter set
  newProb = remake(Problem.ODE,u0 = convert.(eltype(θ),(Problem.ODE.u0)), p=θ)
  sol = solve(newProb,Problem.Alg)

  #Need to loop through experiments to compare to model
  for (times,data) in Problem.Data

    if sol.retcode != :Success
      #println(newp) #what were the bad pars?
      varInfo.logp -= Inf
    else
      #Loop through time
      for (idx,t) in enumerate(times)
        #What states are being measured?
          measured = @. !ismissing(data[idx,:])
        #What are their concentrations?
          obs = ObserveTransform(sol,t,measured)
        #Update the likelihood with each measurement
        logp = Turing.observe(
          sampler,
          likelihood(obs,σ),
          data[idx,measured],
          varInfo)

        varInfo.logp += logp
      end #Time points
    end #If solve succeeds
  end #Loop through data

  return varInfo
end #modelFramework

# Instantiate a Model object.
model = Turing.Model{parTuple, Tuple{:x}}(modelFramework, allData, defaults)
chain = sample(model,sampler,Problem.Samples)
end #mcmcODE function

max_stars_repo_name: RobertGregg/TuringODE

max_issues_repo_name: RobertGregg/TuringODE

max_forks_repo_name: RobertGregg/TuringODE


text:
using ZChop
"""
    QFI(ρ, dρ [, abstol])

Numerically evaluate the quantum Fisher information for the matrix ρ given its derivative dρ wrt the parameter

This function is the implementation of Eq. (13) in Paris, Int. J. Quantum Inform. 7, 125 (2009).

# Arguments
    * `ρ`:  Density matrix
    * `dhro`: Derivative wrt the parameter to be estimated
    * `abstol = 1e-5`: tolerance in the denominator of the formula
"""
function QFI(ρ, dρ; abstol = 1e-5)
    # Get the eigenvalues and eigenvectors of the density matrix
    # We enforce its Hermiticity so that the algorithm is more efficient and returns real values
    eigval, eigvec = eigen(Matrix(zchop(ρ,1e-10)))
    eigval = real(eigval)
    dim = length(eigval)
    return real(2*sum( [( (eigval[n] + eigval[m] > abstol) ? (1. / (eigval[n] + eigval[m])) * abs(eigvec[:,n]' * dρ * eigvec[:,m])^2 : 0.) for n=1:dim, m=1:dim]))
end


text:
@testset "SimResource FIFO/LIFO" begin
    function rfifo(fr)
        tmp = []
        r = SimResource{Int}(quantity=5, lo=0, hi=5,
                             max_requests=3, fifo_requests=fr)
        s = Scheduler()
        @test_throws ErrorException SimResource{Int}(quantity=0, lo=1, hi=5)
        @test_throws ErrorException SimResource{Int}(quantity=3, lo=3, hi=2)

        @test request!(s, r, 2, x -> push!(tmp, 1))[1]
        s.now = 1.0
        @test request!(s, r, 2, x -> push!(tmp, 2))[1]
        @test request!(s, r, 2, x -> push!(tmp, 3))[1]
        @test request!(s, r, 2, x -> push!(tmp, 4))[1]
        @test request!(s, r, 2, x -> push!(tmp, 5))[1]
        @test !request!(s, r, 2, x -> push!(tmp, 6))[1]
        @test length(r.requests) == 3
        s.now = 2.0
        @test provide!(s, r, 1) == 1
        @test request!(s, r, 2, x -> push!(tmp, 7))[1]
        s.now = 3.0
        @test provide!(s, r, 10) == 5
        go!(s, 6.0)
        @test length(r.requests) == 1
        r.requests[1].request(s)
    end

    @test rfifo(true) == [1,2,3,4,5,7]
    @test rfifo(false) == [1,2,5,7,4,3]
end

@testset "SimResource waive! test" begin
    s = Scheduler()
    r = SimResource{Int}()
    @test !waive!(r, ResourceRequest(1, identity))

    rr1 = request!(s, r, 2, x -> error(2))[2]
    rr2 = request!(s, r, 2, x -> error(3))[2]
    rr3 = request!(s, r, 2, x -> error(4))[2]

    @test length(r.requests) == 3
    @test waive!(r, rr2)
    @test length(r.requests) == 2
    @test !waive!(r, rr2)
    @test waive!(r, rr1)
    @test waive!(r, rr3)
    @test isempty(r.requests)
end


text:
# Example file showing the classic 1D robot and four door example
# Best if run piece by piece in a Julia REPL or VSCode IDE

using IncrementalInference

# [OPTIONAL] plotting libraries are loaded separately
using Cairo, RoMEPlotting
Gadfly.set_default_plot_size(35cm,20cm)

## example parameters
# Number of kernels representing each marginal belief
N=100

# prior knowledge of four possible door locations
cv = 3.0
doorPrior = Mixture(Prior,
                    [Normal(-100,cv);Normal(0,cv);Normal(100,cv);Normal(300,cv)],
                    [1/4;1/4;1/4;1/4] )

## Build the factor graph object
fg = initfg()
 # WIP, will become default option in the future
getSolverParams(fg).useMsgLikelihoods = true

# first pose location
v1 = addVariable!(fg,:x1,ContinuousScalar,N=N)
# see a door for the first time
addFactor!(fg,[:x1], doorPrior)

# first solution with only one variable and factor (may take a few moments on first JIT compiling)
solveTree!(fg)
plotKDE(fg, :x1)

## drive to second pose location
addVariable!(fg,:x2, ContinuousScalar, N=N)
addFactor!(fg,[:x1;:x2],LinearRelative(Normal(50.0,2.0)))

# drive to third pose location
v3=addVariable!(fg,:x3,ContinuousScalar, N=N)
addFactor!(fg,[:x2;:x3], LinearRelative( Normal(50.0,4.0)))

# see a door for the second time
addFactor!(fg,[:x3], doorPrior)

# second solution should be much quicker
solveTree!(fg)
plotKDE(fg, [:x1; :x2; :x3])

# drive to forth and final pose location
addVariable!(fg,:x4,ContinuousScalar, N=N)
addFactor!(fg,[:x3;:x4], LinearRelative( Normal(200.0,4.0)))

# lets see the prediction of where pose :x4 might be
ensureAllInitialized!(fg)
plotKDE(fg, :x4)

## make a third door sighting
addFactor!(fg,[:x4], doorPrior)

# solve over all data
tree = solveTree!(fg)

# list variables and factors in fg
@show ls(fg) # |> sortDFG
@show lsf(fg)

## draw all beliefs
pl = plotKDE(fg, [:x1;:x2;:x3;:x4])

#save plot to file
pl |> PNG("4doors.png") # can also do SVG, PDF

## If interested, here is the Bayes/Junction tree too
# drawTree(tree, show=true) # using Graphviz and Linux evince for pdf

#


text:
include("cluster.jl")

Nx = Ny = 16
@show steps = 20000

Jc = 0.4406868                              # critical temp reduced coupling const
J = 1.0                                     # coupling constant in Kelvin
Tc = J/Jc                                   # critical temperature in Kelvin
println("Critical temperature Tc = $Tc")

# Try different lattice sizes
Ns = [16, 32, 128]
# 17 distinct temperatures
Ts = [0.5, 1.0, 1.25, 1.5:0.1:2.2, 2.22, 2.24, 2.26, 2.265, 2.27, 2.28, 2.3:0.1:3.0]
# magnetic fields
Bs = [0.00001, 0.0001, 0.001, 0.01, 0.1]

Tvals = [0.001:0.001:2.269]
Mvals = M(Tvals, J)

# maximum correlation distance
mcd = 12

print("Importing Gadfly... ")
import Gadfly
Gadfly.set_default_plot_size(20Gadfly.cm, 12Gadfly.cm)
println("done.")
# for plotting results
function MeansvsT(means, b, desc::String)

    layers = Array(Gadfly.Layer,0);
    for n in 1:length(Ns)
        newlayer = Gadfly.layer(
            x=Ts, y=means[:,n], Gadfly.Geom.point, Gadfly.Geom.line,
            color=["Simulated for N=$(Ns[n])"]);
        layers = cat(1, layers, newlayer)
    end

    return Gadfly.plot(
        layers,
        Gadfly.layer(
            x=Tvals, y=Mvals, Gadfly.Geom.line,
            color=["Predicted magnetization"]),
        Gadfly.layer(
            x=[Tc,Tc], y=[0,1], Gadfly.Geom.line,
            color=["Critical Temperature 2.269K"]),
        Gadfly.Guide.xticks(ticks=[0:0.5:2, 2.269, 2.5, 3]),
        Gadfly.Guide.yticks(ticks=[0:0.1:1]),
        Gadfly.Guide.xlabel("Temperature"),
        Gadfly.Guide.ylabel("Mean normalized absolute magnetization"),
        Gadfly.Guide.title("$desc vs. Temperature for B=$(Bs[b]) at Different Lattice Sizes")
    )
end

# for plotting results
μvsT(b) = Gadfly.plot(
    Gadfly.layer(
        x=Ts, y=μs[:,b], Gadfly.Geom.point, Gadfly.Geom.line,
        color=["Measured for B=$(Bs[b])"]),
    Gadfly.layer(
        x=Tvals, y=Mvals, Gadfly.Geom.line,
        color=["Predicted magnetization"]),
    Gadfly.Guide.xlabel("Temperature"),
    Gadfly.Guide.ylabel("Mean normalized absolute magnetization"),
    Gadfly.Guide.title("Magnetization vs. Temperature for B=$(Bs[b]) with $(Nx*Ny) Particles")
)

# see magnetization distribution
μhist(mags) = Gadfly.plot(
    Gadfly.Guide.title("Magnetization Histogram for B = 0.1"),
    Gadfly.layer(
        x=abs(mags), color=["Absolute Magnetization"],
        Gadfly.Geom.histogram(bincount=100, density=true), order=1),
    Gadfly.layer(
        x=mags, color=["Magnetization"],
        Gadfly.Geom.histogram(bincount=100, density=true), order=2)
)

# Make plots of the data
absmeans = Array(Float64, length(Ts), length(Ns), length(Bs))
absmeans[:,:,1] = readdlm("absmeans-1.0e-5.dat")
absmeans[:,:,2] = readdlm("absmeans-0.0001.dat")
absmeans[:,:,3] = readdlm("absmeans-0.001.dat")
absmeans[:,:,4] = readdlm("absmeans-0.01.dat")
absmeans[:,:,5] = readdlm("absmeans-0.1.dat")

means = Array(Float64, length(Ts), length(Ns), length(Bs))
means[:,:,1] = readdlm("means-1.0e-5.dat")
means[:,:,2] = readdlm("means-0.0001.dat")
means[:,:,3] = readdlm("means-0.001.dat")
means[:,:,4] = readdlm("means-0.01.dat")
means[:,:,5] = readdlm("means-0.1.dat")

plots = Array(Gadfly.Plot, length(Bs))
for i in 1:length(plots)
    plots[i] = MeansvsT(means[:,:,i], i, "Magnetization")
end

absplots = Array(Gadfly.Plot, length(Bs))
for i in 1:length(plots)
    absplots[i] = MeansvsT(absmeans[:,:,i], i, "Absolute Magnetization")
end


text:
"""
Module with graph theoretical methods for maximum matching, Pantelides index reduction and Tarjan's strongly connected components.

* Author: Hilding Elmqvist, Mogram AB
* Date: July-August 2016 (rewritten). Array version: Spring 2017
* License: MIT

A bipartite graph (E, V) is defined by an array of arrays. Each E-entry has an integer array of indices to V-entries.

Example bipartite graph:

    G = [
      [3, 5],
      [4, 6],
      [1, 7, 9],
      [2, 8, 9],
      [1, 2]
    ]
"""
module BLTandPantelides

export matching, pantelides!, BLT, checkAssign
export matchingArray, pantelidesArray!, BLTArray, compactify
export expandArray, expandArrayOfArray, expandLengths!, scalarToArray, collectDerivatives, invertInvertDer
export expandStringArray
using ..BLTandPantelidesUtilities

"Controls logging"
const log = false

"""
    function augmentPath!(G, i, assign, vColour, eColour, vActive)
Construction of augmenting path

Reference:
Pantelides, C.: The consistent initialization of differential-algebraic systems. SIAM Journal
of Scientific and Statistical Computing, 9(2), pp. 213–231 (1988).
"""
function augmentPath!(G, i, assign, vColour, eColour, vActive)
    # returns pathFound
    # assign: assign[j] contains the E-node to which V-node j is assigned or 0 if V-node j not assigned
    # i: E-node
    # vActive: set to false has the same effect as deleting V-node and corresponding edges
    # j: V-node

    if log
        println("augmentPath: equation $i")
    end

    pathFound = false
    eColour[i] = true

    # If a V-node j exists such that edge (i-j) exists and assign[j] == 0
    for j in G[i]
        if vActive[j] && assign[j] == 0
            pathFound = true
            assign[j] = i
            return pathFound
        end
    end

    # For every j such that edge (i-j) exists and j is uncoloured
    for j in G[i]
        if vActive[j] && !vColour[j]
            vColour[j] = true
            k = assign[j]
            pathFound = augmentPath!(G, k, assign, vColour, eColour, vActive)

            if pathFound
                assign[j] = i
                return pathFound
            end
        end
    end
    return pathFound
end

# function augmentPathArray!(G, i, assign, vColour, eColour, vActive, vLengths, eLengths)
# returns pathFound
# assign: assign[j] contains the E-node to which V-node j is assigned or 0 if V-node j not assigned
# i: E-node
# vActive: set to false has the same effect as deleting V-node and corresponding edges
# j: V-node
function augmentPathArray!(G::Array{Array{Int,1},1}, i::Tuple{Int,Int}, assign::Array{Array{Tuple{Int,Int},1},1}, vColour::Array{Array{Bool,1},1}, eColour::Array{Array{Bool,1},1}, vActive::Array{Bool,1}, vLengths::Array{Int,1}, eLengths::Array{Int,1})::Bool
    if log
        println("\naugmentPath: equation $i")
    end

    pathFound = false
    eColour[i[1]][i[2]] = true

    # If a V-node j exists such that edge (i-j) exists and assign[j] == 0
    for j in G[i[1]]
        # for j in reverse(G[i[1]])
        if vActive[j]
            # Check for all elements of variable j
            for jj in 1:vLengths[j]
                if assign[j][jj] == (0, 0)
                    pathFound = true
                    assign[j][jj] = i
                    if log
                        println("Simply assigned: assign[$j][$jj] = $i")
                    end
                    return pathFound
                end
            end
        end
    end

    # For every j such that edge (i-j) exists and j is uncoloured
    for j in G[i[1]]
        # for j in reverse(G[i[1]])
        if vActive[j]
            # Check for all elements of variable j
            for jj in 1:vLengths[j]
                if !vColour[j][jj]
                    vColour[j][jj] = true
                    k::Tuple{Int,Int} = assign[j][jj]
                    pathFound = augmentPathArray!(G, k, assign, vColour, eColour, vActive, vLengths, eLengths)
                    if pathFound
                        assign[j][jj] = i
                        if log
                            println("Assigned: assign[$j][$jj] = $i")
                        end
                        return pathFound
                    end
                end
            end
        end
    end
    return pathFound
end

function checkAssign(assign, VSizes, VTypes, ESizes, ETypes, equationsInfix, variableNames, A, vActive=[a == 0 for a in A])
    println("Checking assignment")
    assignmentOK = true
    for j in 1:length(assign)
        if vActive[j]
            i = assign[j]
            if i > 0 && VSizes[j] != ESizes[i]
                assignmentOK = false
                print("Error: Variable ")
                printList(variableNames, [j], A, newLine=false)
                println(" (($j)) with size=$(VSizes[j]) is assigned in equation (($i)) with size $(ESizes[i])")
            end
        end
    end

    if assignmentOK
        println("Assignment is OK")
    else
        # error("Assignment not OK")
  end
end

"""
    function matching(G, M)
Find maximum matching in bipartite graph

* `G`: bipartite graph
* `M`: number of V-nodes
* `vActive`: set to false has the same effect as deleting V-node and corresponding edges
* `return assign`: assign[j] contains the E-node to which V-node j is assigned or 0 if V-node j not assigned

Reference:
Pantelides, C.: The consistent initialization of differential-algebraic systems. SIAM Journal
of Scientific and Statistical Computing, 9(2), pp. 213–231 (1988).
"""
function matching(G, M, vActive=fill(true, M))
    assign::Array{Int,1} = fill(0, M)
    for i in 1:length(G)
        eColour::Array{Bool,1} = fill(false, length(G))
        vColour::Array{Bool,1} = fill(false, M)
        pathFound = augmentPath!(G, i, assign, vColour, eColour, vActive)
    end
    return assign
end

function matchingArray(G::Array{Array{Int,1},1}, vActive::Array{Bool,1}, vLengths::Array{Int,1}, eLengths::Array{Int,1})
    local assign::Array{Array{Tuple{Int,Int},1},1}

    assign = fillArrayOfArray((0, 0), vLengths)
    for i in 1:length(G)
        eColour::Array{Array{Bool,1},1} = fillArrayOfArray(false, eLengths)
        vColour::Array{Array{Bool,1},1} = fillArrayOfArray(false, vLengths)
        for ii in 1:eLengths[i]
            pathFound = augmentPathArray!(G, (i, ii), assign, vColour, eColour, vActive, vLengths, eLengths)
        end
    end
    return assign
end

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

"""
    function pantelides!(G, M, A)
Perform index reduction with Pantelides algorithm.

* `G`: bipartite graph (updated)
* `M`: number of V-nodes
* `A`: A[j] = if V[k] = der(V[j]) then k else 0
* `return assign`: assign[j] contains the E-node to which V-node j is assigned or 0 if V-node j not assigned
* `return A`: A[j] = if V[k] = der(V[j]) then k else 0
* `return B`: B[i] = if E[l] = der(E[i]) then l else 0

Reference:
Pantelides, C.: The consistent initialization of differential-algebraic systems. SIAM Journal
of Scientific and Statistical Computing, 9(2), pp. 213–231 (1988).
"""
function pantelides!(G, M, A)
    assign::Array{Int,1} = fill(0, M)
    B::Array{Int,1} = fill(0, length(G))
    N = length(G)
    N2 = N
    for k in 1:N2
        pathFound = false
        i = k
        while !pathFound
            # Delete all V-nodes with A[.] != 0 and all their incidence edges from the graph
            vActive::Array{Bool,1} = [a == 0 for a in A]
            # Designate all nodes as "uncoloured"
            eColour::Array{Bool,1} = fill(false, length(G))
            vColour::Array{Bool,1} = fill(false, M)
            pathFound = augmentPath!(G, i, assign, vColour, eColour, vActive)
            if !pathFound
                if log
                    println("\nDifferentiate:")
                end

                # For every coloured V-node j do
                for j in 1:length(vColour)
                    if vColour[j]
                        M += 1
                        if log
                            println("New variable derivative: var($M) = der($j)")
                        end
                        push!(A, 0)
                        A[j] = M
                        push!(assign, 0)
                    end
                end

                # For every coloured E-node l do
                for l in 1:N
                    if eColour[l]
                        N += 1
                        if log
                            println("New equation derivative:  equ($N) = DER($l)")
                        end

                        # Create new E-node N
                        push!(G, copy(G[l]))

                        # Create edges from E-node N to all V-nodes j and A[j] such that edge (l-j) exists
                        for m in 1:length(G[l])
                            j = G[l][m]
                            if !(A[j] in G[N])
                                push!(G[N], A[j])
                            end
                        end
                        push!(B, 0)

                        # Set B[l] = N
                        B[l] = N
                    end
                end

                # For every coloured V-node j
                for j in 1:length(vColour)
                    if vColour[j]
                        if log
                            println("Assigning derivative of variable $(A[j]) to derivative of equation: $(B[assign[j]])")
                        end
                        assign[A[j]] = B[assign[j]]
                    end
                end
                i = B[i]
            end
        end
    end

    return assign, A, B
end

function fillArrayOfArray(val, Sizes)
    a = Array{Any,1}[]
    for s in Sizes
        push!(a, fill(val, s))
    end
    return a
end

function expandArray(a, lengths)
    offsets = cumsum(lengths) # Offsets in expanded array for each unexpanded variable
    aa = []
    scalarToArray = [] #fill((), offsets[end])  # Not needed?
    for v in 1:length(a)
        d = a[v]
        for i in 1:lengths[v]
            if d != 0
                push!(aa, offsets[d] - lengths[d] + i)
            else
                push!(aa, 0)
            end
            push!(scalarToArray, v)
        end
    end

    # @show aa
    # @show scalarToArray
    return aa #, scalarToArray
end

function expandStringArray(a, lengths, par=false)
    aa = []
    for v in 1:length(a)
        for i in 1:lengths[v]
            if par
                push!(aa, "( " * a[v] * " ) [" * string(i) * "]")
            else
                push!(aa, a[v] * "[" * string(i) * "]")
            end
        end
    end
    return aa
end

function expandArrayOfArray(G, vLengths, eLengths)
    vOffsets = cumsum(vLengths)
    GG = []
    scalarToArrayEquation = []
    for e in 1:length(G)
        g = G[e]
        gg = []

        for v in g
            for i in 1:vLengths[v]
                push!(gg, vOffsets[v] - vLengths[v] + i)
            end
        end

        for j in 1:eLengths[e]
            push!(GG, gg)
            push!(scalarToArrayEquation, (e, j))
        end
    end

    return GG # , scalarToArrayEquation
end

function expandLengths!(lengths, AB)
    (orgIndex, derOrder) = invertDer(AB)
    # @show lengths
    # @show AB
    # @show orgIndex derOrder

    for i in length(lengths) + 1:length(AB)
        # @show i
        org = orgIndex[i]
        # @show org
        push!(lengths, lengths[scalarToArray(lengths)[org]])
    end
end

function scalarToArray(lengths)
    s2A = []
    for i in 1:length(lengths)
        for j in 1:lengths[i]
            push!(s2A, i)
        end
    end
    return s2A
end

function collectDerivatives(AABB, lengths)
    (orgIndex, derOrder) = invertDer(AABB)
    @show AABB lengths
    @show orgIndex derOrder
    n = count(derOrder .== 0)
    orgInd = orgIndex[n + 1:end]
    derOrd = derOrder[n + 1:end]
    p = sortperm(orgInd)
    orgInd2 = orgInd[p]
    derOrd2 = derOrd[p]
    p2 = sortperm(derOrd2)
    derOrd3 = derOrd2[p2]
    orgInd3 = orgInd2[p2]
    orgIndex2 = [1:n; orgInd3]
    derOrder2 = [fill(0, n); derOrd3]

    @show orgIndex2 derOrder2
    return orgIndex2, derOrder2
end

function collectDerivativesOld(AABB, lengths)
    # Collect derivative variables or equations together
    # AABB updated
    AB = [] # New A or B array for array variables or equations
    scalarToArray = fill(0, length(AABB))
    k = 0 # Position in AB array

    # Sample last array element from the AABB list
    # since if one element is differentiated, all other elements are also
    for i in cumsum(lengths)
        index = i
        aabb = AABB[i]
        if aabb == 0
            k += 1
            push!(AB, 0)
            scalarToArray[index] = k
        elseif aabb > 0
            # Place variable and all its derivatives after each other
            while aabb > 0
                k += 1
                push!(AB, k + 1)
                scalarToArray[index] = k
                # (AABB[aabb], aabb) = (-1, AABB[aabb])  # Mark variable as already collected
                index = aabb
                aabb = AABB[aabb]
                AABB[index] = -1
            end

            k += 1
            push!(AB, 0)
            scalarToArray[index] = k
        end
    end
    return AB, scalarToArray
end

"""
    function pantelidesArray!(G, M, A, vLengths, eLengths)
Perform index reduction for array equations and variables with Pantelides algorithm.

* `G`: bipartite graph (updated)
* `M`: number of V-nodes
* `A`: A[j] = if V[k] = der(V[j]) then k else 0 (updated)
* `vLengths`: number of elements of the (array) variables, i.e. product of the dimensions (updated)
* `eLengths`: number of elements of the (array) equations, i.e. product of the dimensions (updated)
* `return assign`: assign[j] contains the E-node to which V-node j is assigned or 0 if V-node j not assigned
* `return A`: A[j] = if V[k] = der(V[j]) then k else 0
* `return B`: B[i] = if E[l] = der(E[i]) then l else 0
"""
function pantelidesArray!(G::Array{Array{Int,1},1}, M::Int, A::Array{Int,1}, vLengths::Array{Int,1}, eLengths::Array{Int,1})
    local assign::Array{Array{Tuple{Int,Int},1},1}
    local B::Array{Int,1}

    assign = fillArrayOfArray((0, 0), vLengths)
    B = fill(0, length(G))
    N = length(G)
    N2 = N
    for k in 1:N2  ### Gives two split variables
        # for k in N2:-1:1
        for ii in 1:eLengths[k]
            # for ii in eLengths[k]:-1:1
            pathFound = false
            i = k
            while !pathFound
                # Delete all V-nodes with A[.] != 0 and all their incidence edges from the graph
                vActive::Array{Bool,1} = [a == 0 for a in A]

                # Designate all nodes as "uncoloured"
                eColour::Array{Array{Bool,1},1} = fillArrayOfArray(false, eLengths)
                vColour::Array{Array{Bool,1},1} = fillArrayOfArray(false, vLengths)
                pathFound = augmentPathArray!(G, (i, ii), assign, vColour, eColour, vActive, vLengths, eLengths)
                if !pathFound
                    if log
                        println("\nDifferentiate:")
                    end
                    # For every coloured V-node j do
                    for j in 1:length(vColour)
                        for jj in 1:vLengths[j]
                            if A[j] == 0 && vColour[j][jj]
                                # if A[j] == 0 && any(vColour[j])
                                M += 1
                                if log
                                    println("New variable derivative: var($M) = der($j)")
                                end
                                push!(A, 0)
                                A[j] = M
                                push!(vLengths, vLengths[j]) # The new differentiated variable has the same length as the original.
                                push!(assign, fill((0, 0), vLengths[j]))
                            end
                        end
                    end

                    # For every coloured E-node l do
                    for l in 1:N
                        for ll in 1:eLengths[l]
                            if B[l] == 0 && eColour[l][ll]
                                # if B[l] == 0 && any(eColour[l])
                                N += 1
                                if log
                                    println("New equation derivative:  equ($N) = DER($l)")
                                end

                                # Create new E-node N
                                push!(G, copy(G[l]))
                                # push!(G, [])  # Problem with wrong sorting was due to this temporary test :-(

                                # Create edges from E-node N to all V-nodes j and A[j] such that edge (l-j) exists
                                for m in 1:length(G[l])
                                    j = G[l][m]

                                    # if ! (A[j] in G[N])
                                    if !(A[j] in G[l])  ### Strange difference to function pantelides!
                                        push!(G[N], A[j])
                                    end
                                end
                                push!(B, 0)

                                # Set B[l] = N
                                B[l] = N
                                push!(eLengths, eLengths[l]) # The new differentiated equation has the same length as the original.
                            end
                        end
                    end

                    # For every coloured V-node j
                    # Assign each variable derivative to the differential of the equation to which the variable was originally assigned
                    for j in 1:length(vColour)
                        for jj in 1:vLengths[j]
                            if vColour[j][jj]
                                # if any(vColour[j])
                                # assign[A[j]] = B[assign[j]]   # Scalar version
                                assign[A[j]][jj] = (B[assign[j][jj][1]], assign[j][jj][2]) # Same element of the derivative array as in the original variable is assigned.
                                if log
                                    println("Assigning derivative of variable ($(A[j]), $jj) to derivative of equation: $((B[assign[j][jj][1]], assign[j][jj][2]))")
                                end
                            end
                        end
                    end
                    i = B[i]
                end
            end
        end
    end
    return assign, A, B
end

const notOnStack = 1000000000

"""
Find minimal systems of equations that have to be solved simultaneously.
Reference:
Tarjan, R. E. (1972), "Depth-first search and linear graph algorithms", SIAM Journal on Computing 1 (2): 146–160, doi:10.1137/0201010
"""
function strongConnect(G, assign, v, nextnode, stack, components, lowlink, number)
    # println("strongConnect: ", v)

    if v == 0
        return nextnode
    end

    nextnode += 1
    lowlink[v] = number[v] = nextnode
    push!(stack, v)

    for w in [assign[j] for j in G[v]] # for w in the adjacency list of v
        if w > 0   # Is assigned
            if number[w] == 0 # if not yet numbered
                nextnode = strongConnect(G, assign, w, nextnode, stack, components, lowlink, number)
                lowlink[v] = min(lowlink[v], lowlink[w])
            else
                if number[w] < number[v]
                    # (v, w) is a frond or cross-link
                    # if w is on the stack of points. Always valid since otherwise number[w]=notOnStack (a big number)
                    lowlink[v] = min(lowlink[v], number[w])
                end
            end
        end
    end

    if lowlink[v] == number[v]
        # v is the root of a component
        # start a new strongly connected component
        # println("start a new strongly connected component")
        # @show v lowlink[v] number[v]
        # @show stack
        # @show number[v]
        comp = []
        repeat = true
        while repeat
            # delete w from point stack and put w in the current component
            # println("delete w from point stack and put w in the current component")
            w = pop!(stack)
            # @show w number[w]
            number[w] = notOnStack
            push!(comp, w)
            # @show comp
            repeat = w != v
        end
        push!(components, comp)
    end
    return nextnode
end

const notOnStack = 1000000000

"""
    function strongConnectArray(G, assign, v, nextnode, stack, components, lowlink, number, vLengths, eLengths)

Find minimal systems of equations that have to be solved simultaneously.
Reference:
Tarjan, R. E. (1972), "Depth-first search and linear graph algorithms", SIAM Journal on Computing 1 (2): 146–160, doi:10.1137/0201010
"""
function strongConnectArray(G::Array{Array{Int,1},1}, assign::Array{Array{Tuple{Int,Int},1},1}, v::Tuple{Int,Int}, nextnode::Int, stack::Array{Tuple{Int,Int},1}, components::Array{Array{Tuple{Int,Int},1},1}, lowlink::Array{Array{Int,1},1}, number::Array{Array{Int,1},1}, vLengths::Array{Int,1}, eLengths::Array{Int,1})

    #  println("\nstrongConnect: ", v)

    (i, ii) = v
    # @show v, nextnode
    if v == (0, 0) # 0   ### When can this happen?
        return nextnode
    end

    nextnode += 1
    lowlink[i][ii] = number[i][ii] = nextnode
    push!(stack, v)

    for w in [assign[j][jj] for j in G[i] for jj in 1:vLengths[j]] # for w in the adjacency list of v
        if w != (0, 0)   # Is assigned
            k = w[1]
            kk = w[2]
            if number[k][kk] == 0 # if not yet numbered
                nextnode = strongConnectArray(G, assign, w, nextnode, stack, components, lowlink, number, vLengths, eLengths)
                lowlink[i][ii] = min(lowlink[i][ii], lowlink[k][kk])
            else
                if number[k][kk] < number[i][ii]
                    # (v, w) is a frond or cross-link
                    # if w is on the stack of points. Always valid since otherwise number[w]=notOnStack (a big number)
                    lowlink[i][ii] = min(lowlink[i][ii], number[k][kk])
                end
            end
        end
    end

    if lowlink[i][ii] == number[i][ii]
        # v is the root of a component
        # start a new strongly connected component
        comp = Tuple{Int,Int}[]
        repeat = true
        while repeat
            # delete w from point stack and put w in the current component
            w = pop!(stack)
            (k, kk) = w
            number[k][kk] = notOnStack
            push!(comp, w)
            if log
                println("Added equation $w to $comp")
            end
            repeat = w != v
        end
        push!(components, comp)
        if log
            println("Added component $comp to $components")
        end
    end
    return nextnode
end

"""
    function BLT(G, assign)

Find Block Lower Triangular structure for a bipartite graph `G` with assignment `assign`

* `G`: bipartite graph
* `assign`: assign[j] contains the E-node to which V-node j is assigned or 0 if V-node j not assigned
* `return components`: cell array of components. Each component is a list of indices to E-nodes
"""
function BLT(G, assign)
    nextnode::Int = 0
    stack = []
    components = []
    lowlink = fill(0, length(G))
    number = fill(0, length(G))

    for v in 1:length(G)
        if number[v] == 0
            nextnode = strongConnect(G, assign, v, nextnode, stack, components, lowlink, number)
        end
    end
    return components
end

"""
    function BLT(G, assign)

Find Block Lower Triangular structure for a bipartite graph `G` with assignment `assign`

* `G`: bipartite graph
* `assign`: assign[j] contains the E-node to which V-node j is assigned or 0 if V-node j not assigned
* `return components`: cell array of components. Each component is a list of indices to E-nodes
"""
function BLTArray(G::Array{Array{Int,1},1}, assign::Array{Array{Tuple{Int,Int},1},1}, vLengths::Array{Int,1}, eLengths::Array{Int,1})
    nextnode = 0
    stack = Tuple{Int,Int}[]
    components = Array{Tuple{Int,Int},1}[]
    lowlink::Array{Array{Int,1},1} = fillArrayOfArray(0, eLengths)
    number::Array{Array{Int,1},1} = fillArrayOfArray(0, eLengths)

    for i in 1:length(G)
        for ii in 1:eLengths[i]
            if number[i][ii] == 0
                nextnode = strongConnectArray(G, assign, (i, ii), nextnode, stack, components, lowlink, number, vLengths, eLengths)
            end
        end
    end
    return components
end

function compactify(components, assign, vLengths, eLengths, vNames, equationsInfix)
    # @show vLengths eLengths
    componentsArray = [unique([e[1] for e in c]) for c in components]
    # @show componentsArray

    # println("\nVerify that all array equation elements are within the same strongly connected component")
    for c in components
        es = [ee[1] for ee in c] # Collect all array equations involved in c
        for e in c
            if count(e[1] .== es) != eLengths[e[1]] # For each array equation in c, make sure all elements are involved in c.
                println("\nThe equation $e:")
                println(equationsInfix[e[1]])
                println("is split into several strongly connected components:")
                println(c)
            end
        end
    end
    compArr = collect(Iterators.flatten(componentsArray))
    # @show size(compArr)
    # @show size(unique(compArr))
    @assert size(compArr) == size(unique(compArr)) "Error: array split into different strongly connected components."

    # Collect all variables assigned in each component
    invAssign = zeros(Int, length(eLengths))
    assignArray = zeros(Int, length(vLengths))
    for c in components
        if length(c) > 9 ###
            println("Block")
        end

        # @show c
        # @show length(c)
        vars = []
        for e in c
            # @show e
            # println(equationsInfix[e[1]])

            # Search for all variables assigned in all equations of the component
            for i in 1:length(assign)
                for j in 1:length(assign[i])
                    # @show e[1] assign[i][j]
                    if e == assign[i][j]   # All equation elements are within the same component.
                        push!(vars, i)
                        # print("push ", vNames[i])
                        if length(assign[i]) > 1
                            # print("[$j]")
                        end
                        # println()
                    end
                end
            end
        end

        if length(c) > 9 ###
            println("end block")
        end
        # @show length(vars)
        vars = unique(vars) # Remove duplicates
        # @show length(vars)
        # @show vars
        cc = unique([e[1] for e in c])
        # @show cc
        # @show length(cc) length(vars)
        # @show sum(eLengths[e] for e in cc)
        # @show sum(vLengths[v] for v in vars)
        # assert(length(cc) == length(vars))
        scalarV = 0
        scalarE = 0
        for i in 1:length(vars)
            # @show cc[i] vars[i]
            invAssign[cc[i]] = vars[i]
            assignArray[vars[i]] = cc[i]
            scalarV += vLengths[vars[i]]
            scalarE += eLengths[cc[i]]
        end
        # @show scalarV scalarE
        if scalarV != scalarE
            for i in 1:length(vars)
                println(vLengths[vars[i]], " ", vNames[vars[i]])
            end
            for i in 1:length(vars)
                println(eLengths[cc[i]], " ", equationsInfix[cc[i]])
            end
        # assert(false)
        end
    end
    # @show invAssign

    for e in 1:length(invAssign)
        if invAssign[e] > 0
            # println("equation ", e, "  variable: ", invAssign[e], " ", vNames[invAssign[e]])
        else
            println("Error: equation ", e, " not assigned")
            # assert(false)
        end
    end
    return invAssign, assignArray, componentsArray
end

# Could not use using BLTandPantelidesUtilities!

function invertDer(A)
    orgIndex = [i for i in 1:length(A)]  # Index of original variable or equation
    derOrder = fill(0, length(A)) # Derivative order
    for i in 1:length(A)
        a = A[i]
        if a > 0
            derOrder[a] = derOrder[i] + 1
            orgIndex[a] = orgIndex[i]
        end
    end
    return orgIndex, derOrder
end

function invertInvertDer(orgIndex, derOrder)
    A = fill(0, length(orgIndex))
    for i in 1:length(orgIndex)
        if derOrder[i] > 0
            e = orgIndex[i]
            for d in 1:derOrder[i] - 1
                e = A[e]
            end
            A[e] = i
        end
    end
    return A
end

function invertAssign(assign, n=length(assign))
    invAssign = fill(0, n)
    unAssigned::Vector{Int} = []
    for j in 1:length(assign)
        i = assign[j]
        if i > 0
            invAssign[i] = j
        else
            push!(unAssigned, j)
        end
    end
    return invAssign, unAssigned
end

function printAssignedEquations(equations, variables, indices, assign, A, B)
    (orgIndexVar, derOrderVar) = invertDer(A)
    (orgIndexEqu, derOrderEqu) = invertDer(B)
    (assignedVar,) = invertAssign(assign)
    for i in indices
        j = assignedVar[i]
        if j > 0
            if derOrderVar[j] == 1
                print("der(")
            elseif derOrderVar[j] > 1
                print("der", derOrderVar[j], "(")
            end
            print(variables[orgIndexVar[j]])
            if derOrderVar[j] > 0
                print(")")
            end
            print("   (($j))   :=   ")
        end

        if derOrderEqu[i] == 1
            print("DER()   ")
        elseif derOrderEqu[i] > 1
            print("DER", derOrderEqu[i], "()   ")
        end
        println(equations[orgIndexEqu[i]], "   (($i))")
    end
end

function printList(infixes, indices, A, vertical=false; newLine=true)
    (orgIndex, derOrder) = invertDer(A)
    for ind in 1:length(indices)
        j = indices[ind]
        if j > 0
            if ind > 1
                if vertical
                    println()
                else
                    print(", ")
                end
            end

            if derOrder[j] > 0
                if vertical
                    print("DER")
                else
                    print("der")
                end

                if derOrder[j] > 1
                    print(derOrder[j])
                end
                print("(")
            end

            print(infixes[orgIndex[j]])

            if derOrder[j] > 0
                print(")")
            end
        end
    end

    if newLine
        println()
    end
end

end