text:
function run_up_down_tests()

    TESTDIR = joinpath(tempdir(), "TEST_")
    rm(TESTDIR; force = true, recursive = true)

    try
        # create files
        _createfile(TESTDIR, "tesfile")
        _createfile(TESTDIR, "sub1", "tesfile")
        _createfile(TESTDIR, "sub1", "sub11", "tesfile")
        _createfile(TESTDIR, "sub1", "sub11", "sub111", "tesfile")
        _createfile(TESTDIR, "sub1", "sub12", "tesfile")
        _createfile(TESTDIR, "sub1", "sub12", "sub121", "tesfile")
        _createfile(TESTDIR, "sub2", "tesfile")
        _createfile(TESTDIR, "sub2", "sub21", "tesfile")
        _createfile(TESTDIR, "sub2", "sub21", "sub211", "tesfile")
        _createfile(TESTDIR, "sub2", "sub22", "tesfile")
        _createfile(TESTDIR, "sub2", "sub22", "sub221", "tesfile")

        @info("Testing walkdown", nthreads())

        # all
        for (thfrec, keepout, ref_filecount) in [
                (0.0, (dir) -> false, 11),
                (0.5, (dir) -> false, 11),
                (1.0, (dir) -> false, 11),
                (0.0, (dir) -> (basename(dir) == "sub2"), 6),
                (0.5, (dir) -> (basename(dir) == "sub2"), 6),
                (1.0, (dir) -> (basename(dir) == "sub2"), 6),
            ]

            @info("At", thfrec, ref_filecount)

            for _ in 1:10
                # all threaded
                filecounts = Dict{Int, Int}()
                walkdown(TESTDIR; keepout, thfrec) do path
                    thid = threadid()
                    get!(filecounts, thid, 0)

                    isfile(path) && (filecounts[thid] += 1)
                    thfrec > 0.0 && sleep(0.1)
                    false
                end
                filecount = sum(values(filecounts))
                @show filecounts
                @test filecount == ref_filecount
            end

            println()
        end

        @info("Testing filterdown", nthreads())

        for (thfrec, filter_, keepout, ref_filecount) in [
                (0.0, isdir, (dir) -> false, 10),
                (0.5, isdir, (dir) -> false, 10),
                (1.0, isdir, (dir) -> false, 10),
                (0.0, isdir, (dir) -> (basename(dir) == "sub2"), 6),
                (0.5, isdir, (dir) -> (basename(dir) == "sub2"), 6),
                (1.0, isdir, (dir) -> (basename(dir) == "sub2"), 6),

                (0.0, isfile, (dir) -> false, 11),
                (0.5, isfile, (dir) -> false, 11),
                (1.0, isfile, (dir) -> false, 11),
                (0.0, isfile, (dir) -> (basename(dir) == "sub2"), 6),
                (0.5, isfile, (dir) -> (basename(dir) == "sub2"), 6),
                (1.0, isfile, (dir) -> (basename(dir) == "sub2"), 6),
            ]

            @info("At", thfrec, nameof(filter_), ref_filecount)

            for _ in 1:10
                arr = filterdown(filter_, TESTDIR; thfrec, keepout)
                @show isfile.(arr)
                @test length(arr) == ref_filecount
            end

            println()
        end

        @info("Testing walkup and fiterup", nthreads())

        # all
        for (dir0, ref_filecount) in [
                (joinpath(TESTDIR, "sub2", "sub21", "sub211"), 4),
                (joinpath(TESTDIR, "sub2", "sub21"), 3),
                (joinpath(TESTDIR, "sub2"), 2),
                (joinpath(TESTDIR), 1),

            ]

            @info("At", dir0, ref_filecount)

            # all threaded
            filecount = 0
            walkup(dir0; root = TESTDIR) do path
                isfile(path) && (filecount += 1)
                false
            end
            @show filecount
            @test filecount == ref_filecount

            founds = filterup(isfile, dir0; root = TESTDIR)
            filecount = length(founds)
            @test filecount == ref_filecount

            println()
        end

    finally
        rm(TESTDIR; force = true, recursive = true)
    end

end
run_up_down_tests()


text:
using GR
using SpecialFunctions

x = 0:0.1:20
j0, j1, j2, j3 = [besselj.(order, x) for order in 0:3]

plot(x, j0, x, j1, x, j2, x, j3,
     labels=["J_0", "J_1", "J_2", "J_3"], location=11)

ext: jl


text:
#!/usr/bin/env julia

const STYLE="""

"""

function figure(src, caption)
"""

$caption

"""
end

function link(src)
"""
$src
"""
end

function makelink(filename)
    if isimagefile(filename)
        figure(basename(filename), basename(filename))
    else
        link(basename(filename))
    end
end

function isimagefile(f)
    f = lowercase(f)
    endswith(f, ".png") || endswith(f, ".jpg") || endswith(f,".jpeg")
end

function mybasename(x)
    y =basename(x)
    i = first(search(y,"."))
    y[1:i-1]
end

function relatedfiles(f)
    b = mybasename(f)
    filenames =  filter(x->startswith(x,b), readdir("."))
end

# Make webpages to allow navigation of results

function makehtml(filenames)

    n = length(filenames)

    for i in 1:n
        filename = filenames[i]
        info("filename $filename")
        name, ext = splitext(mybasename(filename))
        out =  "$name.html"

        info("Writing $out ...")

        open(out, "w") do f
            write(f, """

$STYLE

""")

            if i > 1
                x, ext = splitext(mybasename(filenames[i-1]))
                write(f, "previous\n")
            end
            if i < n
                x, ext = splitext(mybasename(filenames[i+1]))
                write(f, "next\n")
            end
            write(f, "
$name
\n")

            for rf in relatedfiles(filename)
                write(f, makelink(rf))
                write(f,"
")
            end

            write(f, "\n")
        end
    end
end

function makeindex(filenames)
    out = "index.html"
    info("Writing $out ...")

    open(out, "w") do f

        write(f, "$STYLE\n")

        for filename in filenames
            #name, ext = splitext(basename(filename))
            name = mybasename(filename)
            out2=  "$name.html"
            write(f, "
$name
\n")
        end
        write(f, "\n")
    end
end

function main(args)

    info("Writing HTML ...")

    filenames = args

    makehtml(filenames)
    makeindex(filenames)

    info("Done")

end

main(ARGS)


text:
###############################################################################
#
#   fmpz_mat.jl : Flint matrices over fmpz
#
###############################################################################

export fmpz_mat, FmpzMatSpace, getindex, getindex!, setindex!, rows, cols,
       charpoly, det, det_divisor, det_given_divisor, gram, hadamard,
       ishadamard, hnf, ishnf, hnf_with_transform, hnf_modular, lll, lll_ctx,
       lll_gram, lll_with_transform, lll_gram_with_transform, lll_with_removal,
       lll_with_removal_transform, nullspace, rank, rref, reduce_mod, similar,
       snf, snf_diagonal, issnf, solve, solve_rational, cansolve,
       cansolve_with_nullspace, solve_dixon, trace, transpose, content, hcat,
       vcat, addmul!, zero!, window, pseudo_inv, hnf_modular_eldiv,
       nullspace_right_rational

###############################################################################
#
#   Data type and parent object methods
#
###############################################################################

elem_type(::Type{FmpzMatSpace}) = fmpz_mat

parent_type(::Type{fmpz_mat}) = FmpzMatSpace

base_ring(a::FmpzMatSpace) = a.base_ring

parent(a::fmpz_mat, cached::Bool = true) =
    FmpzMatSpace(rows(a), cols(a), cached)

function check_parent(a::fmpz_mat, b::fmpz_mat)
   (rows(a) != rows(b) || cols(a) != cols(b)) && error("Incompatible matrices")
end

###############################################################################
#
#   Similar
#
###############################################################################

function similar(x::fmpz_mat)
   z = fmpz_mat(rows(x), cols(x))
   z.base_ring = x.base_ring
   return z
end

function similar(x::fmpz_mat, r::Int, c::Int)
   z = fmpz_mat(r, c)
   z.base_ring = x.base_ring
   return z
end

###############################################################################
#
#   Windows - handle with care!!!
#
###############################################################################

function Base.view(x::fmpz_mat, r1::Int, c1::Int, r2::Int, c2::Int)
  Generic._checkbounds(x, r1, c1)
  Generic._checkbounds(x, r2, c2)
  (r1 > r2 || c1 > c2) && error("Invalid parameters")
  b = fmpz_mat()
  b.base_ring = FlintZZ
  ccall((:fmpz_mat_window_init, :libflint), Void,
        (Ref{fmpz_mat}, Ref{fmpz_mat}, Int, Int, Int, Int),
            b, x, r1 - 1, c1 - 1, r2, c2)
  finalizer(b, _fmpz_mat_window_clear_fn)
  return b
end

function Base.view(x::fmpz_mat, r::UnitRange{Int}, c::UnitRange{Int})
  return Base.view(x, r.start, c.start, r.stop, c.stop)
end

function _fmpz_mat_window_clear_fn(a::fmpz_mat)
   ccall((:fmpz_mat_window_clear, :libflint), Void, (Ref{fmpz_mat},), a)
end

size(x::fmpz_mat) = tuple(rows(x), cols(x))

size(t::fmpz_mat, d) = d <= 2 ? size(t)[d] : 1

###############################################################################
#
#   Basic manipulation
#
###############################################################################

function getindex!(v::fmpz, a::fmpz_mat, r::Int, c::Int)
   z = ccall((:fmpz_mat_entry, :libflint), Ref{fmpz},
             (Ref{fmpz_mat}, Int, Int), a, r - 1, c - 1)
   ccall((:fmpz_set, :libflint), Void, (Ref{fmpz}, Ref{fmpz}), v, z)
end

@inline function getindex(a::fmpz_mat, r::Int, c::Int)
   @boundscheck Generic._checkbounds(a, r, c)
   v = fmpz()
   z = ccall((:fmpz_mat_entry, :libflint), Ref{fmpz},
             (Ref{fmpz_mat}, Int, Int), a, r - 1, c - 1)
   ccall((:fmpz_set, :libflint), Void, (Ref{fmpz}, Ref{fmpz}), v, z)
   return v
end

@inline function setindex!(a::fmpz_mat, d::fmpz, r::Int, c::Int)
   @boundscheck Generic._checkbounds(a, r, c)
   z = ccall((:fmpz_mat_entry, :libflint), Ref{fmpz},
             (Ref{fmpz_mat}, Int, Int), a, r - 1, c - 1)
   ccall((:fmpz_set, :libflint), Void, (Ref{fmpz}, Ref{fmpz}), z, d)
end

@inline setindex!(a::fmpz_mat, d::Integer, r::Int, c::Int) = setindex!(a, fmpz(d), r, c)

@inline function setindex!(a::fmpz_mat, d::Int, r::Int, c::Int)
   @boundscheck Generic._checkbounds(a, r, c)
   z = ccall((:fmpz_mat_entry, :libflint), Ref{fmpz},
             (Ref{fmpz_mat}, Int, Int), a, r - 1, c - 1)
   ccall((:fmpz_set_si, :libflint), Void, (Ref{fmpz}, Int), z, d)
end

@inline rows(a::fmpz_mat) = a.r

@inline cols(a::fmpz_mat) = a.c

zero(a::FmpzMatSpace) = a()

one(a::FmpzMatSpace) = a(1)

iszero(a::fmpz_mat) = ccall((:fmpz_mat_is_zero, :libflint), Bool,
                            (Ref{fmpz_mat},), a)

isone(a::fmpz_mat) = ccall((:fmpz_mat_is_one, :libflint), Bool,
                           (Ref{fmpz_mat},), a)

function deepcopy_internal(d::fmpz_mat, dict::ObjectIdDict)
   z = fmpz_mat(d)
   z.base_ring = d.base_ring
   return z
end

###############################################################################
#
#   Canonicalisation
#
###############################################################################

canonical_unit(a::fmpz_mat) = canonical_unit(a[1, 1])

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

function show(io::IO, a::FmpzMatSpace)
   print(io, "Matrix Space of ")
   print(io, a.rows, " rows and ", a.cols, " columns over ")
   print(io, "Integer Ring")
end

function show(io::IO, a::fmpz_mat)
   r = rows(a)
   c = cols(a)
   for i = 1:r
      print(io, "[")
      for j = 1:c
         print(io, a[i, j])
         if j != c
            print(io, " ")
         end
      end
      print(io, "]")
      if i != r
         println(io, "")
      end
   end
end

show_minus_one(::Type{fmpz_mat}) = show_minus_one(fmpz)

###############################################################################
#
#   Unary operations
#
###############################################################################

function -(x::fmpz_mat)
   z = similar(x)
   ccall((:fmpz_mat_neg, :libflint), Void,
         (Ref{fmpz_mat}, Ref{fmpz_mat}), z, x)
   return z
end

###############################################################################
#
#   transpose
#
###############################################################################

function transpose(x::fmpz_mat)
   z = similar(x, cols(x), rows(x))
   ccall((:fmpz_mat_transpose, :libflint), Void,
         (Ref{fmpz_mat}, Ref{fmpz_mat}), z, x)
   return z
end

###############################################################################
#
#   Binary operations
#
###############################################################################

function +(x::fmpz_mat, y::fmpz_mat)
   check_parent(x, y)
   z = similar(x)
   ccall((:fmpz_mat_add, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat},  Ref{fmpz_mat}),
               z, x, y)
   return z
end

function -(x::fmpz_mat, y::fmpz_mat)
   check_parent(x, y)
   z = similar(x)
   ccall((:fmpz_mat_sub, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat},  Ref{fmpz_mat}),
               z, x, y)
   return z
end

function *(x::fmpz_mat, y::fmpz_mat)
   cols(x) != rows(y) && error("Incompatible matrix dimensions")
   if rows(x) == cols(y) && rows(x) == cols(x)
      z = similar(x)
   else
      z = similar(x, rows(x), cols(y))
   end
   ccall((:fmpz_mat_mul, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat},  Ref{fmpz_mat}),
               z, x, y)
   return z
end

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

function *(x::Int, y::fmpz_mat)
   z = similar(y)
   ccall((:fmpz_mat_scalar_mul_si, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat}, Int), z, y, x)
   return z
end

function *(x::fmpz, y::fmpz_mat)
   z = similar(y)
   ccall((:fmpz_mat_scalar_mul_fmpz, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz}), z, y, x)
   return z
end

*(x::fmpz_mat, y::Int) = y*x

*(x::fmpz_mat, y::fmpz) = y*x

*(x::Integer, y::fmpz_mat) = fmpz(x)*y

*(x::fmpz_mat, y::Integer) = fmpz(y)*x

function +(x::fmpz_mat, y::Integer)
   z = deepcopy(x)
   for i = 1:min(rows(x), cols(x))
      z[i, i] += y
   end
   return z
end

function +(x::fmpz_mat, y::fmpz)
   z = deepcopy(x)
   for i = 1:min(rows(x), cols(x))
      z[i, i] = addeq!(z[i, i], y)
   end
   return z
end

+(x::Integer, y::fmpz_mat) = y + x

+(x::fmpz, y::fmpz_mat) = y + x

-(x::fmpz_mat, y::Integer) = x + (-y)

-(x::fmpz_mat, y::fmpz) = x + (-y)

function -(x::Integer, y::fmpz_mat)
   z = -y
   for i = 1:min(rows(y), cols(y))
      z[i, i] += x
   end
   return z
end

function -(x::fmpz, y::fmpz_mat)
   z = -y
   for i = 1:min(rows(y), cols(y))
      z[i, i] = addeq!(z[i, i], x)
   end
   return z
end

###############################################################################
#
#   Scaling
#
###############################################################################

doc"""
    <<(x::fmpz_mat, y::Int)
> Return $2^yx$.
"""
function <<(x::fmpz_mat, y::Int)
   y < 0 && throw(DomainError())
   z = similar(x)
   ccall((:fmpz_mat_scalar_mul_2exp, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat}, Int),
               z, x, y)
   return z
end

doc"""
    >>(x::fmpz_mat, y::Int)
> Return $x/2^y$ where rounding is towards zero.
"""
function >>(x::fmpz_mat, y::Int)
   y < 0 && throw(DomainError())
   z = similar(x)
   ccall((:fmpz_mat_scalar_tdiv_q_2exp, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat}, Int),
               z, x, y)
   return z
end

###############################################################################
#
#   Powering
#
###############################################################################

function ^(x::fmpz_mat, y::Int)
   y < 0 && throw(DomainError())
   rows(x) != cols(x) && error("Incompatible matrix dimensions")
   z = similar(x)
   ccall((:fmpz_mat_pow, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat}, Int),
               z, x, y)
   return z
end

###############################################################################
#
#   Comparisons
#
###############################################################################

function ==(x::fmpz_mat, y::fmpz_mat)
   check_parent(x, y)
   ccall((:fmpz_mat_equal, :libflint), Bool,
                                       (Ref{fmpz_mat}, Ref{fmpz_mat}), x, y)
end

###############################################################################
#
#   Ad hoc comparisons
#
###############################################################################

function ==(x::fmpz_mat, y::Integer)
   for i = 1:min(rows(x), cols(x))
      if x[i, i] != y
         return false
      end
   end
   for i = 1:rows(x)
      for j = 1:cols(x)
         if i != j && x[i, j] != 0
            return false
         end
      end
   end
   return true
end

==(x::Integer, y::fmpz_mat) = y == x

==(x::fmpz_mat, y::fmpz) = x == parent(x)(y)

==(x::fmpz, y::fmpz_mat) = parent(y)(x) == y

###############################################################################
#
#   Inversion
#
###############################################################################

function inv(x::fmpz_mat)
   z = similar(x)
   d = fmpz()
   ccall((:fmpz_mat_inv, :libflint), Void,
         (Ref{fmpz_mat}, Ref{fmpz}, Ref{fmpz_mat}), z, d, x)
   if d == 1
      return z
   end
   if d == -1
      return -z
   end
   error("Matrix not invertible")
end

###############################################################################
#
#   Pseudo inversion
#
###############################################################################

doc"""
    pseudo_inv(x::fmpz_mat)
> Return a tuple $(z, d)$ consisting of a matrix $z$ and denominator $d$ such
> that $z/d$ is the inverse of $x$.
"""
function pseudo_inv(x::fmpz_mat)
   z = similar(x)
   d = fmpz()
   ccall((:fmpz_mat_inv, :libflint), Void,
         (Ref{fmpz_mat}, Ref{fmpz}, Ref{fmpz_mat}), z, d, x)
   if !iszero(d)
      return (z, d)
   end
   error("Matrix is singular")
end

###############################################################################
#
#   Exact division
#
###############################################################################

function divexact(x::fmpz_mat, y::fmpz_mat)
   cols(x) != cols(y) && error("Incompatible matrix dimensions")
   x*inv(y)
end

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

function divexact(x::fmpz_mat, y::Int)
   z = similar(x)
   ccall((:fmpz_mat_scalar_divexact_si, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat}, Int), z, x, y)
   return z
end

function divexact(x::fmpz_mat, y::fmpz)
   z = similar(x)
   ccall((:fmpz_mat_scalar_divexact_fmpz, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz}), z, x, y)
   return z
end

divexact(x::fmpz_mat, y::Integer) = divexact(x, fmpz(y))

###############################################################################
#
#   Modular reduction
#
###############################################################################

doc"""
    reduce_mod(x::fmpz_mat, y::fmpz)
> Reduce the entries of $x$ modulo $y$ and return the result.
"""
function reduce_mod(x::fmpz_mat, y::fmpz)
   z = similar(x)
   ccall((:fmpz_mat_scalar_mod_fmpz, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz}), z, x, y)
   return z
end

doc"""
    reduce_mod(x::fmpz_mat, y::Integer)
> Reduce the entries of $x$ modulo $y$ and return the result.
"""
reduce_mod(x::fmpz_mat, y::Integer) = reduce_mod(x, fmpz(y))

###############################################################################
#
#   Characteristic polynomial
#
###############################################################################

function charpoly(R::FmpzPolyRing, x::fmpz_mat)
   rows(x) != cols(x) && error("Non-square")
   z = R()
   ccall((:fmpz_mat_charpoly, :libflint), Void,
                (Ref{fmpz_poly}, Ref{fmpz_mat}), z, x)
   return z
end

###############################################################################
#
#   Minimal polynomial
#
###############################################################################

function minpoly(R::FmpzPolyRing, x::fmpz_mat)
   rows(x) != cols(x) && error("Non-square")
   z = R()
   ccall((:fmpz_mat_minpoly, :libflint), Void,
                (Ref{fmpz_poly}, Ref{fmpz_mat}), z, x)
   return z
end

###############################################################################
#
#   Determinant
#
###############################################################################

function det(x::fmpz_mat)
   rows(x) != cols(x) && error("Non-square matrix")
   z = fmpz()
   ccall((:fmpz_mat_det, :libflint), Void,
                (Ref{fmpz}, Ref{fmpz_mat}), z, x)
   return z
end

doc"""
    det_divisor(x::fmpz_mat)
> Return some positive divisor of the determinant of $x$, if the determinant
> is nonzero, otherwise return zero.
"""
function det_divisor(x::fmpz_mat)
   rows(x) != cols(x) && error("Non-square matrix")
   z = fmpz()
   ccall((:fmpz_mat_det_divisor, :libflint), Void,
                (Ref{fmpz}, Ref{fmpz_mat}), z, x)
   return z
end

doc"""
    det_given_divisor(x::fmpz_mat, d::fmpz, proved=true)
> Return the determinant of $x$ given a positive divisor of its determinant. If
> `proved == true` (the default), the output is guaranteed to be correct,
> otherwise a heuristic algorithm is used.
"""
function det_given_divisor(x::fmpz_mat, d::fmpz, proved=true)
   rows(x) != cols(x) && error("Non-square")
   z = fmpz()
   ccall((:fmpz_mat_det_modular_given_divisor, :libflint), Void,
               (Ref{fmpz}, Ref{fmpz_mat}, Ref{fmpz}, Cint), z, x, d, proved)
   return z
end

doc"""
    det_given_divisor(x::fmpz_mat, d::Integer, proved=true)
> Return the determinant of $x$ given a positive divisor of its determinant. If
> `proved == true` (the default), the output is guaranteed to be correct,
> otherwise a heuristic algorithm is used.
"""
function det_given_divisor(x::fmpz_mat, d::Integer, proved=true)
   return det_given_divisor(x, fmpz(d), proved)
end

###############################################################################
#
#   Gram matrix
#
###############################################################################

function gram(x::fmpz_mat)
   z = similar(x, rows(x), rows(x))
   ccall((:fmpz_mat_gram, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat}), z, x)
   return z
end

###############################################################################
#
#   Hadamard matrix
#
###############################################################################

doc"""
    hadamard(R::FmpzMatSpace)
> Return the Hadamard matrix for the given matrix space. The number of rows and
> columns must be equal.
"""
function hadamard(R::FmpzMatSpace)
   R.rows != R.cols && error("Unable to create Hadamard matrix")
   z = R()
   success = ccall((:fmpz_mat_hadamard, :libflint), Bool,
                   (Ref{fmpz_mat},), z)
   !success && error("Unable to create Hadamard matrix")
   return z
end

doc"""
    ishadamard(x::fmpz_mat)
> Return `true` if the given matrix is Hadamard, otherwise return `false`.
"""
function ishadamard(x::fmpz_mat)
   return ccall((:fmpz_mat_is_hadamard, :libflint), Bool,
                   (Ref{fmpz_mat},), x)
end

###############################################################################
#
#   Hermite normal form
#
###############################################################################

doc"""
    hnf(x::fmpz_mat)
> Return the Hermite Normal Form of $x$.
"""
function hnf(x::fmpz_mat)
   z = similar(x)
   ccall((:fmpz_mat_hnf, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat}), z, x)
   return z
end

doc"""
    hnf_with_transform(x::fmpz_mat)
> Compute a tuple $(H, T)$ where $H$ is the Hermite normal form of $x$ and $T$
> is a transformation matrix so that $H = Tx$.
"""
function hnf_with_transform(x::fmpz_mat)
   z = similar(x)
   u = similar(x, rows(x), rows(x))
   ccall((:fmpz_mat_hnf_transform, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz_mat}), z, u, x)
   return z, u
end

doc"""
    hnf_modular(x::fmpz_mat, d::fmpz)
> Compute the Hermite normal form of $x$ given that $d$ is a multiple of the
> determinant of the nonzero rows of $x$.
"""
function hnf_modular(x::fmpz_mat, d::fmpz)
   z = similar(x)
   ccall((:fmpz_mat_hnf_modular, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz}), z, x, d)
   return z
end

doc"""
    hnf_modular_eldiv(x::fmpz_mat, d::fmpz)
> Compute the Hermite normal form of $x$ given that $d$ is a multiple of the
> largest elementary divisor of $x$. The matrix $x$ must have full rank.
"""
function hnf_modular_eldiv(x::fmpz_mat, d::fmpz)
   (rows(x) < cols(x)) &&
                error("Matrix must have at least as many rows as columns")
   z = deepcopy(x)
   ccall((:fmpz_mat_hnf_modular_eldiv, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz}), z, d)
   return z
end

doc"""
    ishnf(x::fmpz_mat)
> Return `true` if the given matrix is in Hermite Normal Form, otherwise return
> `false`.
"""
function ishnf(x::fmpz_mat)
   return ccall((:fmpz_mat_is_in_hnf, :libflint), Bool,
                   (Ref{fmpz_mat},), x)
end

###############################################################################
#
#   LLL
#
###############################################################################

mutable struct lll_ctx
   delta::Float64
   eta::Float64
   rep_type::Int
   gram_type::Int

   function lll_ctx(delta::Float64, eta::Float64,
                    rep::Symbol = :zbasis, gram::Symbol = :approx)
      rt = rep == :zbasis ? 1 : 0
      gt = gram == :approx ? 0 : 1
      return new(delta, eta, rt, gt)
   end
end

doc"""
> Compute a tuple $(L, T)$ where $L$ is the LLL reduction of $a$ and $T$ is a
> transformation matrix so that $L = Ta$. All the default parameters can be
> overridden by supplying an optional context object.
"""
function lll_with_transform(x::fmpz_mat, ctx::lll_ctx = lll_ctx(0.99, 0.51))
   z = deepcopy(x)
   u = similar(x, rows(x), rows(x))
   for i in 1:rows(u)
      u[i, i] = 1
   end
   ccall((:fmpz_lll, :libflint), Void,
         (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{lll_ctx}), z, u, ctx)
   return z, u
end

doc"""
    lll(x::fmpz_mat, ctx=lll_ctx(0.99, 0.51))
> Return the LLL reduction of the matrix $x$. By default the matrix $x$ is a
> $\mathbb{Z}$-basis and the Gram matrix is maintained throughout in
> approximate form. The LLL is performed with reduction parameters
> $\delta = 0.99$ and $\eta = 0.51$. All of these defaults can be overridden by
> specifying an optional context object.
"""
function lll(x::fmpz_mat, ctx::lll_ctx = lll_ctx(0.99, 0.51))
   z = deepcopy(x)
   if rows(z) == 0
     return z
   end
   u = similar(x, rows(x), rows(x))
   ccall((:fmpz_lll, :libflint), Void,
         (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{lll_ctx}), z, u, ctx)
   return z
end

doc"""
    lll_gram_with_transform(x::fmpz_mat, ctx=lll_ctx(0.99, 0.51, :gram))
> Given the Gram matrix $x$ of a matrix $M$, compute a tuple $(L, T)$ where
> $L$ is the gram matrix of the LLL reduction of the matrix and $T$ is a
> transformation matrix so that $L = TM$.
"""
function lll_gram_with_transform(x::fmpz_mat, ctx::lll_ctx = lll_ctx(0.99, 0.51, :gram))
   z = deepcopy(x)
   u = similar(x, rows(x), rows(x))
   for i in 1:rows(u)
      u[i, i] = 1
   end
   ccall((:fmpz_lll, :libflint), Void,
         (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{lll_ctx}), z, u, ctx)
   return z, u
end

doc"""
    lll_gram(x::fmpz_mat, ctx=lll_ctx(0.99, 0.51, :gram))
> Given the Gram matrix $x$ of a matrix, compute the Gram matrix of its LLL
> reduction.
"""
function lll_gram(x::fmpz_mat, ctx::lll_ctx = lll_ctx(0.99, 0.51, :gram))
   z = deepcopy(x)
   u = similar(x, rows(x), rows(x))
   ccall((:fmpz_lll, :libflint), Void,
         (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{lll_ctx}), z, u, ctx)
   return z
end

doc"""
    lll_with_removal_transform(x::fmpz_mat, b::fmpz, ctx=lll_ctx(0.99, 0.51))
> Compute a tuple $(r, L, T)$ where the first $r$ rows of $L$ are those
> remaining from the LLL reduction after removal of vectors with norm exceeding
> the bound $b$ and $T$ is a transformation matrix so that $L = Tx$.
"""
function lll_with_removal_transform(x::fmpz_mat, b::fmpz, ctx::lll_ctx = lll_ctx(0.99, 0.51))
   z = deepcopy(x)
   u = similar(x, rows(x), rows(x))
   for i in 1:rows(u)
      u[i, i] = 1
   end
   d = Int(ccall((:fmpz_lll_with_removal, :libflint), Cint,
    (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz}, Ref{lll_ctx}), z, u, b, ctx))
   return d, z, u
end

doc"""
    lll_with_removal(x::fmpz_mat, b::fmpz, ctx=lll_ctx(0.99, 0.51))
> Compute the LLL reduction of $x$ and throw away rows whose norm exceeds
> the given bound $b$. Return a tuple $(r, L)$ where the first $r$ rows of $L$
> are the rows remaining after removal.
"""
function lll_with_removal(x::fmpz_mat, b::fmpz, ctx::lll_ctx = lll_ctx(0.99, 0.51))
   z = deepcopy(x)
   u = similar(x, rows(x), rows(x))
   d = Int(ccall((:fmpz_lll_with_removal, :libflint), Cint,
    (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz}, Ref{lll_ctx}), z, u, b, ctx))
   return d, z
end

###############################################################################
#
#   Nullspace
#
###############################################################################

function nullspace(x::fmpz_mat)
  H, T = hnf_with_transform(x')
  for i = rows(H):-1:1
    for j = 1:cols(H)
      if !iszero(H[i, j])
        N = similar(x, cols(x), rows(H) - i)
        for k = 1:rows(N)
          for l = 1:cols(N)
            N[k, l] = T[rows(T) - l + 1, k]
          end
        end
        return N, cols(N)
      end
    end
  end
  return similar(x, cols(x), 0), 0
end

doc"""
    nullspace_right_rational(x::fmpz_mat)
> Return the right rational nullspace of $x$, i.e. a set of vectors over
> $\mathbb{Z}$ giving a $\mathbb{Q}$-basis for the nullspace of $x$
> considered as a matrix over $\mathbb{Q}$.
"""
function nullspace_right_rational(x::fmpz_mat)
   z = similar(x)
   u = similar(x, cols(x), cols(x))
   rank = ccall((:fmpz_mat_nullspace, :libflint), Cint,
                (Ref{fmpz_mat}, Ref{fmpz_mat}), u, x)
   return u, rank
end

###############################################################################
#
#   Rank
#
###############################################################################

function rank(x::fmpz_mat)
   return ccall((:fmpz_mat_rank, :libflint), Int,
                (Ref{fmpz_mat},), x)
end

###############################################################################
#
#   Reduced row echelon form
#
###############################################################################

function rref(x::fmpz_mat)
   z = similar(x)
   d = fmpz()
   ccall((:fmpz_mat_rref, :libflint), Void,
         (Ref{fmpz_mat}, Ref{fmpz}, Ref{fmpz_mat}), z, d, x)
   return z, d
end

###############################################################################
#
#   Smith normal form
#
###############################################################################

doc"""
    snf(x::fmpz_mat)
> Compute the Smith normal form of $x$.
"""
function snf(x::fmpz_mat)
   z = similar(x)
   ccall((:fmpz_mat_snf, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat}), z, x)
   return z
end

doc"""
    snf_diagonal(x::fmpz_mat)
> Given a diagonal matrix $x$ compute the Smith normal form of $x$.
"""
function snf_diagonal(x::fmpz_mat)
   z = similar(x)
   ccall((:fmpz_mat_snf_diagonal, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat}), z, x)
   return z
end

doc"""
    issnf(x::fmpz_mat)
> Return `true` if $x$ is in Smith normal form, otherwise return `false`.
"""
function issnf(x::fmpz_mat)
   return ccall((:fmpz_mat_is_in_snf, :libflint), Bool,
                   (Ref{fmpz_mat},), x)
end

###############################################################################
#
#   Linear solving
#
###############################################################################

doc"""
    solve(a::fmpz_mat, b::fmpz_mat) -> fmpz_mat
> Return a matrix $x$ such that $ax = b$. An exception is raised
> if this is not possible.
"""
function solve(a::fmpz_mat, b::fmpz_mat)
   rows(b) != rows(a) && error("Incompatible dimensions in solve")
   fl, z = cansolve(a, b)
   if !fl
     error("system is inconsistent")
   end
   return z
end

doc"""
    cansolve(a::fmpz_mat, b::fmpz_mat) -> Bool, fmpz_mat
> Return true and a matrix $x$ such that $ax = b$, or false and some matrix
> in case $x$ does not exist.
"""
function cansolve(a::fmpz_mat, b::fmpz_mat)
   rows(b) != rows(a) && error("Incompatible dimensions in cansolve")
   H, T = hnf_with_transform(a')
   b = deepcopy(b)
   z = similar(a, cols(b), cols(a))
   l = min(rows(a), cols(a))
   for i = 1:cols(b)
     for j = 1:l
       k = 1
       while k <= cols(H) && iszero(H[j, k])
         k += 1
       end
       if k > cols(H)
         continue
       end
       q, r = divrem(b[k, i], H[j, k])
       if !iszero(r)
         return false, b
       end
       for h = k:cols(H)
         b[h, i] -= q*H[j, h]
       end
       z[i, j] = q
     end
   end
   if !iszero(b)
     return false, b
   end
   return true, (z*T)'
end

doc"""
    cansolve_with_nullspace(a::fmpz_mat, b::fmpz_mat) -> Bool, fmpz_mat, fmpz_mat
> Return true, a matrix $x$ and a matrix $k$ such that $ax = b$ and the columns
> of $k$ form a basis for the nullspace of $a$. In case $x$ does not exist, false
> and two arbitrary matrices are returned.
"""
function cansolve_with_nullspace(a::fmpz_mat, b::fmpz_mat)
   rows(b) != rows(a) && error("Incompatible dimensions in cansolve_with_nullspace")
   H, T = hnf_with_transform(a')
   z = similar(a, cols(b), cols(a))
   l = min(rows(a), cols(a))
   for i=1:cols(b)
     for j=1:l
       k = 1
       while k <= cols(H) && iszero(H[j, k])
         k += 1
       end
       if k > cols(H)
         continue
       end
       q, r = divrem(b[k, i], H[j, k])
       if !iszero(r)
         return false, b, b
       end
       for h=k:cols(H)
         b[h, i] -= q*H[j, h]
       end
       z[i, k] = q
     end
   end
   if !iszero(b)
     return false, b, b
   end

   for i = rows(H):-1:1
     for j = 1:cols(H)
       if !iszero(H[i,j])
         N = similar(a, cols(a), rows(H) - i)
         for k = 1:rows(N)
           for l = 1:cols(N)
             N[k,l] = T[rows(T) - l + 1, k]
           end
         end
         return true, (z*T)', N
       end
     end
   end
   N =  similar(a, cols(a), 0)

   return true, (z*T), N
end

doc"""
    solve_rational(a::fmpz_mat, b::fmpz_mat)
> If it exists, return a tuple $(x, d)$ consisting of a column vector $x$ such
> that $ax = db$. The element $b$ must be a column vector with the same number
> of rows as $a$ and $a$ must be a square matrix. If these conditions are not
> met or $(x, d)$ does not exist, an exception is raised.
"""
function solve_rational(a::fmpz_mat, b::fmpz_mat)
   rows(a) != cols(a) && error("Not a square matrix in solve_rational")
   rows(b) != rows(a) && error("Incompatible dimensions in solve_rational")
   z = similar(b)
   d = fmpz()
   nonsing = ccall((:fmpz_mat_solve, :libflint), Bool,
      (Ref{fmpz_mat}, Ref{fmpz}, Ref{fmpz_mat}, Ref{fmpz_mat}), z, d, a, b)
   !nonsing && error("Singular matrix in solve_rational")
   return z, d
end

function Generic.solve_with_det(a::fmpz_mat, b::fmpz_mat)
   return solve_rational(a, b)
end

doc"""
    solve_dixon(a::fmpz_mat, b::fmpz_mat)
> Return a tuple $(x, m)$ consisting of a column vector $x$ such that $ax = b
> \pmod{m}$. The element  $b$ must be a column vector with the same number > of
> rows as $a$ and $a$ must be a square matrix. If these conditions are not met
> or $(x, d)$ does not exist, an exception is raised.
"""
function solve_dixon(a::fmpz_mat, b::fmpz_mat)
   rows(a) != cols(a) && error("Not a square matrix in solve")
   rows(b) != rows(a) && error("Incompatible dimensions in solve")
   z = similar(b)
   d = fmpz()
   nonsing = ccall((:fmpz_mat_solve_dixon, :libflint), Bool,
      (Ref{fmpz_mat}, Ref{fmpz}, Ref{fmpz_mat}, Ref{fmpz_mat}), z, d, a, b)
   !nonsing && error("Singular matrix in solve")
   return z, d
end

###############################################################################
#
#   Trace
#
###############################################################################

function trace(x::fmpz_mat)
   rows(x) != cols(x) && error("Not a square matrix in trace")
   d = fmpz()
   ccall((:fmpz_mat_trace, :libflint), Int,
                (Ref{fmpz}, Ref{fmpz_mat}), d, x)
   return d
end

###############################################################################
#
#   Content
#
###############################################################################

function content(x::fmpz_mat)
  d = fmpz()
  ccall((:fmpz_mat_content, :libflint), Void,
        (Ref{fmpz}, Ref{fmpz_mat}), d, x)
  return d
end

###############################################################################
#
#   Concatenation
#
###############################################################################

function hcat(a::fmpz_mat, b::fmpz_mat)
  rows(a) != rows(b) && error("Incompatible number of rows in hcat")
  c = similar(a, rows(a), cols(a) + cols(b))
  ccall((:fmpz_mat_concat_horizontal, :libflint), Void,
        (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz_mat}), c, a, b)
  return c
end

function vcat(a::fmpz_mat, b::fmpz_mat)
  cols(a) != cols(b) && error("Incompatible number of columns in vcat")
  c = similar(a, rows(a) + rows(b), cols(a))
  ccall((:fmpz_mat_concat_vertical, :libflint), Void,
        (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz_mat}), c, a, b)
  return c
end

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function mul!(z::fmpz_mat, x::fmpz_mat, y::fmpz_mat)
   ccall((:fmpz_mat_mul, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz_mat}), z, x, y)
   return z
end

function mul!(y::fmpz_mat, x::Int)
   ccall((:fmpz_mat_scalar_mul_si, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat}, Int), y, y, x)
   return y
end

function mul!(y::fmpz_mat, x::fmpz)
   ccall((:fmpz_mat_scalar_mul_fmpz, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz}), y, y, x)
   return y
end

function addmul!(z::fmpz_mat, y::fmpz_mat, x::fmpz)
   ccall((:fmpz_mat_scalar_addmul_fmpz, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz}), z, y, x)
   return y
end

function addmul!(z::fmpz_mat, y::fmpz_mat, x::Int)
   ccall((:fmpz_mat_scalar_addmul_si, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat}, Int), z, y, x)
   return y
end

function addeq!(z::fmpz_mat, x::fmpz_mat)
   ccall((:fmpz_mat_add, :libflint), Void,
                (Ref{fmpz_mat}, Ref{fmpz_mat}, Ref{fmpz_mat}), z, z, x)
   return z
end

function zero!(z::fmpz_mat)
   ccall((:fmpz_mat_zero, :libflint), Void,
                (Ref{fmpz_mat},), z)
   return z
end

###############################################################################
#
#   Parent object call overloads
#
###############################################################################

function (a::FmpzMatSpace)()
   z = fmpz_mat(a.rows, a.cols)
   z.base_ring = FlintZZ
   return z
end

function (a::FmpzMatSpace)(arr::Array{fmpz, 2})
   _check_dim(a.rows, a.cols, arr)
   z = fmpz_mat(a.rows, a.cols, arr)
   z.base_ring = FlintZZ
   return z
end

function (a::FmpzMatSpace)(arr::Array{T, 2}) where {T <: Integer}
   _check_dim(a.rows, a.cols, arr)
   z = fmpz_mat(a.rows, a.cols, arr)
   z.base_ring = FlintZZ
   return z
end

function (a::FmpzMatSpace)(arr::Array{fmpz, 1})
   _check_dim(a.rows, a.cols, arr)
   z = fmpz_mat(a.rows, a.cols, arr)
   z.base_ring = FlintZZ
   return z
end

function (a::FmpzMatSpace)(arr::Array{T, 1}) where {T <: Integer}
   _check_dim(a.rows, a.cols, arr)
   z = fmpz_mat(a.rows, a.cols, arr)
   z.base_ring = FlintZZ
   return z
end

function (a::FmpzMatSpace)(d::fmpz)
   z = fmpz_mat(a.rows, a.cols, d)
   z.base_ring = FlintZZ
   return z
end

function (a::FmpzMatSpace)(d::Integer)
   z = fmpz_mat(a.rows, a.cols, fmpz(d))
   z.base_ring = FlintZZ
   return z
end

(a::FmpzMatSpace)(d::fmpz_mat) = d

###############################################################################
#
#   Conversions and promotions
#
###############################################################################

promote_rule(::Type{fmpz_mat}, ::Type{T}) where {T <: Integer} = fmpz_mat

promote_rule(::Type{fmpz_mat}, ::Type{fmpz}) = fmpz_mat

function convert(::Type{Matrix{Int}}, A::fmpz_mat)
    m, n = size(A)

    fittable = [fits(Int, A[i, j]) for i in 1:m, j in 1:n]
    if !all(fittable)
        error("When trying to convert a fmpz_mat to a Matrix{Int}, some elements were too large to fit into Int: try to convert to a matrix of BigInt.")
    end

    mat::Matrix{Int} = Int[A[i, j] for i in 1:m, j in 1:n]
    return mat
end

function convert(::Type{Matrix{BigInt}}, A::fmpz_mat)
    m, n = size(A)
    # No check: always ensured to fit a BigInt.
    mat::Matrix{BigInt} = BigInt[A[i, j] for i in 1:m, j in 1:n]
    return mat
end

###############################################################################
#
#   Matrix constructor
#
###############################################################################

function matrix(R::FlintIntegerRing, arr::Array{fmpz, 2})
   z = fmpz_mat(size(arr, 1), size(arr, 2), arr)
   z.base_ring = FlintZZ
   return z
end

function matrix(R::FlintIntegerRing, arr::Array{<: Integer, 2})
   z = fmpz_mat(size(arr, 1), size(arr, 2), arr)
   z.base_ring = FlintZZ
   return z
end

function matrix(R::FlintIntegerRing, r::Int, c::Int, arr::Array{fmpz, 1})
   _check_dim(r, c, arr)
   z = fmpz_mat(r, c, arr)
   z.base_ring = FlintZZ
   return z
end

function matrix(R::FlintIntegerRing, r::Int, c::Int, arr::Array{<: Integer, 1})
   _check_dim(r, c, arr)
   z = fmpz_mat(r, c, arr)
   z.base_ring = FlintZZ
   return z
end

###############################################################################
#
#  Zero matrix
#
###############################################################################

function zero_matrix(R::FlintIntegerRing, r::Int, c::Int)
   z = fmpz_mat(r, c)
   z.base_ring = FlintZZ
   return z
end

###############################################################################
#
#  Identity matrix
#
###############################################################################

function identity_matrix(R::FlintIntegerRing, n::Int)
   z = fmpz_mat(n, n)
   ccall((:fmpz_mat_one, :libflint), Void, (Ref{fmpz_mat}, ), z)
   z.base_ring = FlintZZ
   return z
end

###############################################################################
#
#   MatrixSpace constructor
#
###############################################################################

function MatrixSpace(R::FlintIntegerRing, r::Int, c::Int, cached::Bool = true)
   return FmpzMatSpace(r, c, cached)
end


text:
function json_headers!(headers::AbstractDict)
    headers["Accept"] = "application/json"
    headers["Content-Type"] = "application/json"
    return headers
end

function authentication_headers!(headers::AbstractDict, client::Client)::Nothing
    auth = get_auth(client)
    authentication_headers!(headers, auth)
    return nothing
end

function authentication_headers!(headers::AbstractDict, auth::AnonymousAuth)::Nothing
    return nothing
end

function authentication_headers!(headers::AbstractDict, oauth2::OAuth2)::Nothing
    headers["Authorization"] = "Bearer $(get_token(oauth2))"
    return nothing
end

function authentication_headers!(headers::AbstractDict, jwt_auth::JWTAuth)::Nothing
    headers["Authorization"] = "Bearer $(get_token(jwt_auth))"
    return nothing
end

function authentication_headers!(headers::AbstractDict, password_auth::UsernamePassAuth)::Nothing
    headers["Authorization"] = "Basic $(Base64.base64encode(string(get_username(password_auth), ':', get_password(password_auth))))"
    return nothing
end


text:
@symbol_func function K_VI(cur_reactor::AbstractReactor)
  cur_K = cur_reactor.l_i / 2.0

  cur_K -= 3/2

  cur_K += log( 8.0 / cur_reactor.epsilon )

  cur_K
end


text:
version = "3.15.08.02"
os_string = @osx? "darwin" : "linux"
extension = @osx? "zip" : "tar.gz"

file_name = "dReal-$version-$os_string-shared-libs.$extension"
file_url = "https://github.com/dreal/dreal3/releases/download/v$version/$file_name"
deps_dir = joinpath(joinpath(Pkg.dir("DReal"),"deps"))
prefix = joinpath(deps_dir,"usr")

println("Cleanup old build files")

try
  run(`rm  -r $(joinpath(deps_dir,"dReal*"))`)
catch end
try
  run(`rm -r $(joinpath(deps_dir,"usr"))`)
catch end

download(file_url,joinpath(deps_dir,file_name))
@osx? begin
  run(`unzip $file_name`)
end : begin
  run(`tar -xvf $file_name`)
end
run(`mv $(joinpath(deps_dir, "dReal-$version-$os_string-shared-libs")) usr`)


text:
using OrdinaryDiffEq
using Distributions
using DiffEqFlux, Flux
using Random
using Plots;

Random.seed!(123);

function sir_ode(u,p,t)
    (S,I,C) = u
    (β,γ) = p
    dS = -β*S*I
    dI = β*S*I - γ*I
    dC = β*S*I
    [dS,dI,dC]
end;

solver = RadauIIA3();

N = 1000.0
p = [0.5,0.25]
u0 = [0.99, 0.01, 0.0]
tspan = (0., 40.)
δt = 1;

sir_prob = ODEProblem(sir_ode, u0, tspan, p)
sir_sol = solve(sir_prob, solver, saveat = δt);

plot(sir_sol,
     xlabel = "Time",
     ylabel = "Proportion",
     labels = ["S" "I" "R"])

train_time = 30.0
tsdata = Array(sir_sol(0:δt:train_time))
cdata = diff(tsdata[3,:])
noisy_data = rand.(Poisson.(N .* cdata));

plot(1:δt:train_time, N .* cdata,
     xlabel = "Time",
     ylabel = "New cases per day",
     label = "True value")
scatter!(1:δt:train_time, noisy_data, label="Data")

foi1 = FastDense(1, 1, relu, bias=false)
p1_ = Float64.(initial_params(foi1))
length(p1_)

function sir_ude(u,p_,t,foi)
    S,I,C = u
    β,γ = p
    λ = foi([I],p_)[1]
    dS = -λ*S
    dI = λ*S - γ*I
    dC = λ*S
    [dS, dI, dC]
end;

tspan_train = (0,train_time)
sir_ude1 = (u,p_,t) -> sir_ude(u,p_,t,foi1)
prob_ude1 = ODEProblem(sir_ude1,
                      u0,
                      tspan_train,
                      p1_);

function predict(θ, prob)
    Array(solve(prob,
                solver;
                u0 = u0,
                p = θ,
                saveat = δt,
                sensealg = InterpolatingAdjoint(autojacvec=ReverseDiffVJP())))
end;

function loss(θ, prob)
    pred = predict(θ, prob)
    cpred = abs.(N*diff(pred[3,:]))
    Flux.poisson_loss(cpred, float.(noisy_data)), cpred
end;

loss(prob_ude1.p, prob_ude1);

const losses1 = []
callback1 = function (p, l, pred)
    push!(losses1, l)
    numloss = length(losses1)
    if numloss % 10 == 0
        display("Epoch: " * string(numloss) * " Loss: " * string(l))
    end
    return false
end;

res_ude1 = DiffEqFlux.sciml_train((θ)->loss(θ,prob_ude1),
                                  p1_,
                                  cb=callback1);

res_ude1.minimizer, losses1[end]

plot(losses1, xaxis = :log, xlabel = "Iterations", ylabel = "Loss", legend=false)

prob_ude1_fit = ODEProblem(sir_ude1, u0, tspan, res_ude1.minimizer)
sol_ude1_fit = solve(prob_ude1_fit, solver, saveat = δt)
scatter(sir_sol, label=["True Susceptible" "True Infected" "True Recovered"],title="Fitted true model")
plot!(sol_ude1_fit, label=["Estimated Susceptible" "Estimated Infected" "Estimated Recovered"])

Imax = maximum(tsdata[2,:])
Igrid = 0:0.01:0.5
λ = [foi1([I],res_ude1.minimizer)[1] for I in Igrid]
scatter(Igrid,λ,xlabel="Proportion of population infected, I",ylab="Force of infection, λ",label="Neural network prediction")
Plots.abline!(p[1],0,label="True value")
Plots.vline!([Imax],label="Upper bound of training data")

Random.seed!(1234)
nhidden = 4
foi2 = FastChain(FastDense(1, nhidden, relu),
                     FastDense(nhidden, nhidden, relu),
                     FastDense(nhidden, 1, relu))
p2_ = Float64.(initial_params(foi2))
length(p2_)

sir_ude2 = (u,p_,t) -> sir_ude(u,p_,t,foi2)
prob_ude2 = ODEProblem(sir_ude2,
                      u0,
                      tspan_train,
                      p2_);

const losses2 = []
callback2 = function (p, l, pred)
    push!(losses2, l)
    numloss = length(losses2)
    if numloss % 10 == 0
        display("Epoch: " * string(numloss) * " Loss: " * string(l))
    end
    return false
end;

res_ude2 = DiffEqFlux.sciml_train((θ)->loss(θ,prob_ude2),
                                  p2_,
                                  cb = callback2);

losses1[end],losses2[end]

prob_ude2_fit = ODEProblem(sir_ude2, u0, tspan, res_ude2.minimizer)
sol_ude2_fit = solve(prob_ude2_fit, solver, saveat = δt)
scatter(sir_sol, label=["True Susceptible" "True Infected" "True Recovered"],title="Fitted UDE model")
plot!(sol_ude2_fit, label=["Estimated Susceptible" "Estimated Infected" "Estimated Recovered"])

λ = [foi2([I],res_ude2.minimizer)[1] for I in Igrid]
scatter(Igrid, λ, xlabel="Proportion of population infected, i", ylab="Force of infection, λ", label="Neural network prediction")
Plots.abline!(p[1], 0,label="True value")
Plots.vline!([Imax], label="Upper bound of training data")


text:
import JSON

export deserialize

function deserialize(req::String, type)
    if type == Supply
        deserializeSupply(req)
    elseif type == Demand
        deserializeDemand(req)
    else
        throw("Not Supported")
    end
end

# TODO mapping using reflection aka fields in stead of explicit mapping

function deserializeSupply(req::String)::Supply
    reqDict = JSON.parse(req)
    Supply(
        reqDict["itemId"],
        reqDict["uom"],
        reqDict["nodeId"],
        reqDict["type"],
        reqDict["quantity"])
end

function deserializeDemand(req::String)::Demand
    reqDict = JSON.parse(req)
    Demand(
        reqDict["itemId"],
        reqDict["uom"],
        reqDict["nodeId"],
        reqDict["type"],
        reqDict["quantity"])
end

max_stars_repo_name: kawadesumit/availability

max_issues_repo_name: kawadesumit/availability

max_forks_repo_name: kawadesumit/availability


text:
# Create a version of the radon dataset with some counties removed.

using DrWatson
@quickactivate "EPExperiments"

using Random
using CSV
using DataFrames

const NEW_CSV_NAME = "small_radon_num_counties=20.csv"
const NEW_NUM_COUNTIES = 20 # Number of counties in the new dataset

function create_small_radon()
    Random.seed!(1234)

    radon_df = DataFrame(CSV.File(datadir("exp_raw", "radon.csv")))

    num_counties = length(unique(radon_df[:,:county]))
    ixs = randperm(num_counties)[1:NEW_NUM_COUNTIES]
    new_counties = unique(radon_df[:,:county])[ixs]

    small_radon_df = filter(row -> in(row.county, new_counties), radon_df)
    CSV.write(datadir("exp_raw", NEW_CSV_NAME), small_radon_df)
end

create_small_radon()


text:
module HonlyHelperWingsRefill

using ..Ahorn, Maple

@mapdef Entity "HonlyHelper/WingsRefill" WingsUpgrade(x::Integer, y::Integer, oneUse::Bool = false, feathers::Integer = 1)

const placements = Ahorn.PlacementDict(
    "Wings Refill (HonlyHelper)" => Ahorn.EntityPlacement(
        WingsUpgrade
    )
)

sprite = "objects/refill/idle00"

function Ahorn.selection(entity::WingsUpgrade)
    x, y = Ahorn.position(entity)
    return Ahorn.getSpriteRectangle(sprite, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::WingsUpgrade, room::Maple.Room)
    Ahorn.drawSprite(ctx, sprite, 0, 0)
end

end

max_stars_repo_name: HollyMagala/HonlyHelper

max_issues_repo_name: HollyMagala/HonlyHelper

max_forks_repo_name: HollyMagala/HonlyHelper


text:
import MPI
include("03-reduce-impl.jl")

function main()
    MPI.Init()

    do_reduce()

    MPI.Finalize()
end

main()


text:
"""
$(SIGNATURES)

Convenience function to add an id attribute to a html element
"""
attr(name::Symbol, val::AS) = ifelse(isempty(val), "", " $name=\"$val\"")

"""
$(SIGNATURES)

Convenience function for a list item
"""
html_li(in::AS) = "
$(in)
"

"""
$SIGNATURES

Convenience function for a sup item
"""
html_sup(id::String, in::AS) =  "$in"

"""
$(SIGNATURES)

Convenience function for a header
"""
html_hk(hk::String, header::AS; id::String="") = "<$hk$(attr(:id, id))>$header</$hk>"

"""
$(SIGNATURES)

Convenience function to introduce a hyper reference.
"""
function html_ahref(link::AS, name::Union{Int,AS};
                    title::AS="", class::AS="")
    a  = "$name"
    a
end

"""
$(SIGNATURES)

Convenience function to introduce a hyper reference relative to a key (local hyperref).
"""
function html_ahref_key(key::AS, name::Union{Int,AS})
    return html_ahref(url_curpage() * "#$key", name)
end

"""
$(SIGNATURES)

Convenience function to introduce a div block.
"""
html_div(name::AS, in::AS) = "
$in
"

"""
$(SIGNATURES)

Convenience function to introduce an image. The alt is escaped just in case the user adds quotation
marks in the alt string.
"""
html_img(src::AS, alt::AS="") = ""

"""
$(SIGNATURES)

Convenience function to introduce a code block.
"""
function html_code(c::AS, lang::AS="")
    isempty(c) && return ""
    isempty(lang) && return "
$c
"
    return "
$c
"
end

"""
$(SIGNATURES)

Convenience function to introduce inline code.
"""
html_code_inline(c::AS) = "$c"

"""
$(SIGNATURES)

Insertion of a visible red message in HTML to show there was a problem.
"""
html_err(mess::String="") = "
// $mess //
"

"""
$(SIGNATURES)

Helper function to get the relative url of the current page.
"""
function url_curpage()
    # go from /pages/.../something.md to /pub/.../something.html note that if
    # on windows then it would be \\ whence the PATH_SEP
    rp = replace(CUR_PATH[], Regex("^pages$(escape_string(PATH_SEP))")=>"pub$(PATH_SEP)")
    rp = unixify(rp)
    rp = splitext(rp)[1] * ".html"
    startswith(rp, "/") || (rp = "/" * rp)
    return rp
end


text:
using GeoArrayOps
using Test

@testset "pmf" begin
    # Write your own tests here.
    A = rand(25, 25)
    A[2, 2] = NaN
    B, flags = pmf(A)
    @test (A .<= B) == (flags .== 0.0)
    B, flags = pmf(A; circular = true)
    @test (A .<= B) == (flags .== 0.0)
end
@testset "smf" begin
    # Write your own tests here.
    A = rand(25, 25)
    A[2, 2] = NaN
    B = smf(A)
end
@testset "pssm" begin
    B = pssm(rand(25, 25))
end
@testset "psf" begin
    B = psf(rand(25, 25))
end
@testset "pitremoval" begin
    B = pitremoval(rand(25, 25))
end
@testset "spread" begin
    points = [0.0 0 0 0 2; 0 0 0 0 0; 0 0 0 0 0; 0 1 0 0 0; 0 0 0 0 0]
    initial = [8.0 8 8 8 4; 8 8 8 8 8; 8 8 8 8 8; 0 0 8 8 8; 0 0 8 8 8]
    friction = [1.0 200 1 1 1; 200 1 1 4 4; 1 1 4 4 4; 1 1 3 200 200; 1 Inf 3 200 4]
    @test spread(points, initial, friction) == spread2(points, initial, friction)
end
@testset "terrain" begin
    A = rand(25, 25)
    TRI(A)
    TPI(A)
    roughness(A)
end


text:
# Deprecated functionality

@deprecate CuDevice(ctx::CuContext) device(ctx)
@deprecate CuCurrentDevice() current_device()
@deprecate CuCurrentContext() current_context()
@deprecate CuContext(ptr::Union{Ptr,CuPtr}) context(ptr)
@deprecate CuDevice(ptr::Union{Ptr,CuPtr}) device(ptr)

@deprecate CuDefaultStream() default_stream()
@deprecate CuStreamLegacy() legacy_stream()
@deprecate CuStreamPerThread() per_thread_stream()
@deprecate query(s::CuStream) isdone(s)
@deprecate query(e::CuEvent) isdone(e)


text:
#===============================================================================

    Subroutine for the Newtonian Method

    Based on code on Ben Moll's webite:
            www.princton.edu/~moll/HACTproject.htm

===============================================================================#

function huggett_subroutine(r_t, data)

    v_st = data[1]
    gg_0 = data[2]
    Aswitch = data[3]

    SS, dS = [zeros(N,1) for i in 1:2]

    dVf, dVb = [zeros(H,2) for j in 1:2]

    V=v_st

    v_n = zeros(H,2,N)

        for n in N:-1:1
            v[:,:,n]=V
            #forward difference
            dVf[1:H-1,:] = (V[2:H,:]-V[1:H-1,:])/da
            dVf[H,:] = (z .+ r_t[n].*amax).^(-s) #will never be used, but impose state constraint a<=amax just in case
            #backward difference
            dVb[2:H,:] = (V[2:H,:]-V[1:H-1,:])/da
            dVb[1,:] = (z .+ r_t[n].*amin).^(-s) #state constraint boundary condition

            I_concave = dVb .> dVf #indicator whether value function is concave (problems arise if this is not the case)

            # consumption and savings with forward difference
            cf = dVf.^(-1/s)
            ssf = zz + r_t[n].*aa - cf
            # consumption and savings with backward difference
            cb = dVb.^(-1/s)
            ssb = zz + r_t[n].*aa - cb
            # consumption and derivative of value function at steady state
            c0 = zz + r_t[n].*aa
            dV0 = c0.^(-s)

            #= dV_upwind makes a choice of forward or backward differences based on
             the sign of the drift    =#
            If = ssf .> 0 #positive drift → forward difference
            Ib = ssb .< 0 #negative drift → backward difference
            I0 = (1.0 .-If-Ib) #at steady state

            dV_Upwind = dVf.*If + dVb.*Ib + dV0.*I0  #important to include third term
            c = dV_Upwind.^(-1/s)
            u = c.^(1-s)/(1-s)

            # CONSTRUCT MATRIX
            X = - min.(ssb,0)/da
            Y = - max.(ssf,0)/da + min.(ssb,0)/da
            Z = max.(ssf,0)/da

            A1=spdiagm(0 => Y[:,1], -1 => X[2:H,1], 1 =>Z[1:H-1,1])
            A2=spdiagm(0=> Y[:,2], -1 =>X[2:H,2], 1 => Z[1:H-1,2])

            global A = [A1 spzeros(H,H);spzeros(H,H) A2] + Aswitch

            push!(A_t,A)
            B = (1/dt + ρ)*sparse(I, 2*H, 2*H) - A

            u_stacked = [u[:,1];u[:,2]]
            V_stacked = [V[:,1];V[:,2]]

            b = u_stacked + V_stacked/dt
            V_stacked = B\b #SOLVE SYSTEM OF EQUATIONS

            V = [V_stacked[1:H] V_stacked[H+1:2*H]]

            global ss = zz+ r_t[n].*aa -c
        end

        gg[1,:,:] = gg0

        #loop for KFE, solve forward
        for q in 1:N
            AT=A_t[end-(q-1)]' #the indexing is backward from original code.
            gg[q+1,:,:] = (sparse(I,2*H,2*H)-AT*dt)\gg[q,:,:]
            SS[q] = (gg[q,:,:]'*aa[:].*da)[1]
            dS[q] = (gg[q+1,:,:]'*aa[:]*da - gg[q,:,:]'*aa[:]*da)[1]
        end

return dS
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct LabelsApi <: SwaggerApi
    client::Swagger.Client
end

function _swaggerinternal_apiV1FleetLabelsGet(_api::LabelsApi; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Any, "/api/v1/fleet/labels", ["bearerAuth"])
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? [] : [_mediaType])
    return _ctx
end

@doc raw"""Labels

Params:

Return: Any
"""
function apiV1FleetLabelsGet(_api::LabelsApi; _mediaType=nothing)
    _ctx = _swaggerinternal_apiV1FleetLabelsGet(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function apiV1FleetLabelsGet(_api::LabelsApi, response_stream::Channel; _mediaType=nothing)
    _ctx = _swaggerinternal_apiV1FleetLabelsGet(_api; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_apiV1FleetLabelsLabelIdGet(_api::LabelsApi, label_id::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Any, "/api/v1/fleet/labels/{label_id}", ["bearerAuth"])
    Swagger.set_param(_ctx.path, "label_id", label_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? [] : [_mediaType])
    return _ctx
end

@doc raw"""Label by ID

Params:
- label_id::String (required)

Return: Any
"""
function apiV1FleetLabelsLabelIdGet(_api::LabelsApi, label_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_apiV1FleetLabelsLabelIdGet(_api, label_id; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function apiV1FleetLabelsLabelIdGet(_api::LabelsApi, response_stream::Channel, label_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_apiV1FleetLabelsLabelIdGet(_api, label_id; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export apiV1FleetLabelsGet, apiV1FleetLabelsLabelIdGet


text:
mark_spec =
[
	(:name, AbstractString, nothing),
	(:key, AbstractString, nothing),
	(:_type, AbstractString, "symbol"),
	(:from, VegaMarkFrom, nothing),
	(:properties, VegaMarkProperties, nothing),
	(:delay, VegaValueRef, nothing),
	(:ease, AbstractString, nothing),
	(:marks, Vector, nothing),
	(:scales, Vector{VegaScale}, nothing),
	(:axes, Vector{VegaAxis}, nothing),
	(:interactive, Bool, nothing)
]

primitivefactory(:VegaMark, mark_spec)

max_stars_repo_name: jpthor/vega


text:
function normalize!{T}(s::Array{T, 2})
# Normalize s so that the squared L2 vector norm of each row is at most 1

  dictSize = size(s, 1);
  for k=1:dictSize;
    norm_s = norm(s[k, :]);
    if norm_s*norm_s>1;
      s[k, :] /= norm_s;
    end
  end
end

max_stars_repo_name: bkong/FSC

max_issues_repo_name: bkong/FSC

max_forks_repo_name: bkong/FSC


text:
# This file is part of Kpax3. License is MIT.

function initializepartition(settings::KSettings;
                             kset::UnitRange{Int}=1:0)
  if settings.verbose
    Printf.@printf("Computing pairwise distances... ")
  end

  tmp = zeros(UInt8, length(settings.miss))
  idx = 0
  for c in 1:length(settings.miss)
    if settings.miss[c] != UInt8('-')
      idx += 1
      tmp[idx] = settings.miss[c]
    end
  end

  miss = if idx > 0
           copyto!(zeros(UInt8, idx), 1, tmp, 1, idx)
         else
           zeros(UInt8, 1)
         end

  (data, id, ref) = readfasta(settings.ifile, settings.protein, miss,
                              settings.l, false, 0)

  n = size(data, 2)

  d = if settings.protein
        distaamtn84(data, ref)
      else
        distntmtn93(data, ref)
      end

  D = zeros(Float64, n, n)
  idx = 0
  for j in 1:(n - 1), i in (j + 1):n
    idx += 1
    D[i, j] = D[j, i] = d[idx]
  end

  if settings.verbose
    Printf.@printf("done\n")
  end

  # expected number of cluster approximately between cbrt(n) and sqrt(n)
  g = ceil(Int, n^(2 / 5))

  kset = if length(kset) == 0
           max(1, g - 20):min(n, g + 20)
         elseif kset[1] > 0
           if kset[end] > n
             if kset[1] == 1
               2:n
             else
               kset[1]:n
             end
           else
             if kset[1] == 1
               2:kset[end]
             else
               kset[1]:kset[end]
             end
           end
         else
           throw(KDomainError("First element of 'kset' is less than one."))
         end

  if settings.verbose
    Printf.@printf("Initializing partition...\n")
  end

  missval = if (length(settings.miss) == 1) && (settings.miss[1] == 0x00)
              0x00
            else
              UInt8('?')
            end

  (data1, id1, ref1) = readfasta(settings.ifile, settings.protein,
                                 settings.miss, settings.l, false, 0)
  (bindata, val, key) = categorical2binary(data1, UInt8(127), missval)

  priorR = EwensPitman(settings.α, settings.θ)
  priorC = AminoAcidPriorCol(bindata, settings.γ, settings.r)

  s = initializestate(bindata, D, kset, priorR, priorC, settings)

  if settings.verbose
    Printf.@printf("Initialization done.\n")
  end

  normalizepartition(s.R, n)
end

function initializepartition(data::Matrix{UInt8},
                             D::Matrix{Float64},
                             settings::KSettings;
                             kset::UnitRange{Int}=1:0)
  n = size(data, 2)

  # expected number of cluster approximately between cbrt(n) and sqrt(n)
  g = ceil(Int, n^(2 / 5))

  kset = if length(kset) == 0
           max(1, g - 20):min(n, g + 20)
         elseif kset[1] > 0
           if kset[end] > n
             if kset[1] == 1
               2:n
             else
               kset[1]:n
             end
           else
             if kset[1] == 1
               2:kset[end]
             else
               kset[1]:kset[end]
             end
           end
         else
           throw(KDomainError("First element of 'kset' is less than one."))
         end

  if settings.verbose
    Printf.@printf("Initializing partition...\n")
  end

  priorR = EwensPitman(settings.α, settings.θ)
  priorC = AminoAcidPriorCol(data, settings.γ, settings.r)

  s = initializestate(data, D, kset, priorR, priorC, settings)

  if settings.verbose
    Printf.@printf("Initialization done.\n")
  end

  normalizepartition(s.R, n)
end

"""
remove "gaps" and non-positive values from the partition
Example: [1, 1, 0, 1, -2, 4, 0] -> [3, 3, 2, 3, 1, 4, 2]
"""
function normalizepartition(partition::Vector{Int},
                            n::Int)
  if length(partition) != n
    throw(KInputError(string("Length of argument 'partition' is not equal to ",
                             "the sample size: ", length(partition),
                             " instead of ", n, ".")))
  end

  map(Int, indexin(partition, sort(unique(partition))))
end

function normalizepartition(ifile::String,
                            n::Int)
  d = DelimitedFiles.readdlm(ifile, ',', Int)

  if size(d, 2) != 1
    throw(KInputError(string("Too many columns found in file ", ifile, ".")))
  end

  if length(d) != n
    throw(KInputError(string("Partition length is not equal to the sample ",
                             "size: ", length(d), " instead of ", n, ".")))
  end

  map(Int, indexin(d[:, 1], sort(unique(d[:, 1]))))
end

function normalizepartition(ifile::String,
                            id::Vector{String})
  d = DelimitedFiles.readdlm(ifile, ',', String)

  if size(d, 1) != length(id)
    throw(KInputError(string("Partition length is not equal to the sample ",
                             "size: ", size(d, 1), " instead of ", length(id),
                             ".")))
  end

  if size(d, 2) == 1
    partition = [parse(Int, x) for x in d[:, 1]]
    map(Int, indexin(partition, sort(unique(partition))))
  elseif size(d, 2) == 2
    idx = map(Int, indexin(id, [x::String for x in d[:, 1]]))
    partition = [parse(Int, x) for x in d[:, 2]]

    for i in 1:length(idx)
      if idx[i] == 0
        throw(KInputError(string("Missing ids in file ", ifile, ".")))
      end
    end

    map(Int, indexin(partition, sort(unique(partition)))[idx])
  else
    throw(KInputError(string("Too many columns found in file ", ifile, ".")))
  end
end

function decodepartition!(R::Vector{Int})
  for a in 1:length(R)
    R[a] = 1 + div(R[a] - a, length(R))
  end

  nothing
end

function encodepartition!(R::Vector{Int},
                          S::Vector{Int})
  n = length(S)

  lidx = zeros(Int, n)

  g = 0
  for a in 1:n
    g = S[a]

    if lidx[g] == 0
      lidx[g] = a
    end

    # linear index of unit a in the n-by-n adjacency matrix
    R[a] = a + n * (lidx[g] - 1)
  end

  nothing
end

function modifypartition!(R::Vector{Int},
                          k::Int)
  n = length(R)
  q = k

  if (k > 0) && (k < n + 1)
    q = rand(max(1, k - 10):min(n, k + 10))

    if q < k
      modifymerge!(R, k, q)
    elseif q > k
      modifysplit!(R, k, q)
    else
      modifyscramble!(R, k)
    end
  end

  q
end

function modifymerge!(R::Vector{Int},
                      k::Int,
                      q::Int)
  n = length(R)

  cset = zeros(Int, k)
  empty = trues(n)
  l = 1
  for a in 1:n
    if empty[R[a]]
      empty[R[a]] = false
      cset[l] = R[a]
      l += 1
    end
  end

  c = zeros(Int, 2)

  while k != q
    StatsBase.sample!(cset, c, replace=false, ordered=false)

    for a in 1:n
      if R[a] == c[2]
        R[a] = c[1]
      end
    end

    l = 1
    while cset[l] != c[2]
      l += 1
    end

    deleteat!(cset, l)
    k -= 1
  end

  nothing
end

function modifysplit!(R::Vector{Int},
                      k::Int,
                      q::Int)
  n = length(R)

  t = zeros(Int, n)
  for a in 1:n
    t[R[a]] += 1
  end

  g = 0

  while k != q
    k += 1

    w = StatsBase.ProbabilityWeights(Float64[t[a] > 0 ? t[a] - 1 : 0
                                             for a in 1:n])
    g = StatsBase.sample(w)

    for a in 1:n
      if (R[a] == g) && ((t[k] == 0) || (rand() <= 0.25))
        R[a] = k
        t[g] -= 1
        t[k] += 1
      end
    end
  end

  nothing
end

function modifyscramble!(R::Vector{Int},
                         k::Int)
  n = length(R)

  t = zeros(Int, n)
  for a in 1:n
    t[R[a]] += 1
  end

  v = zeros(Int, n)
  moved = false

  l = 1
  g = 1
  h = 1
  while g < k
    if t[l] > 0
      t[l] = 0
      v[l] = 0

      for a in (l + 1):n
        v[a] = t[a] > 0 ? t[a] - 1 : 0
      end

      w = StatsBase.ProbabilityWeights(v)
      h = StatsBase.sample(w)

      keepgoing = true
      moved = false
      a = 1
      while keepgoing
        if R[a] == h
          if moved
            if rand() <= 0.05
              R[a] = g
              t[h] -= 1

              if t[h] == 1
                keepgoing = false
              end
            end
          else
            moved = true
            R[a] = g
            t[h] -= 1

            if t[h] == 1
              keepgoing = false
            end
          end
        end

        a += 1
        if a > n
          keepgoing = false
        end
      end

      g += 1
    end

    l += 1
  end

  nothing
end