text:
export summarize_pdf

function bin_location(v::Float64, start::Int64)::Int64
    return round(abs(v) - abs(start)) + 1
end

# Current only use psd_mean to get pdf_mean
function summarize_pdf(psd_mean::Array{<:Real, 2}; min_db::Int64=-200, max_db::Int64=-50)
    pdf_bins = abs(max_db - min_db) + 1
    num_onehoursegments, freq_len = size(psd_mean)
    pdf_mean = Base.zeros(eltype(psd_mean), pdf_bins, freq_len)
    #pdf_max = zeros(pdf_bins, freq_len)
    #pdf_min = zeros(pdf_bins, freq_len)
    #pdf_medain = zeros(pdf_bins, freq_len)

    # Histogram accumulation
    for i in 1:freq_len
        for j in 1:num_onehoursegments
            pdf_mean[bin_location(psd_mean[j, i], max_db), i] += 1
            #pdf_max[] += 1
            #pdf_min[] += 1
            #pdf_median[] += 1
        end
    end

    # Calculate the PDF, i.e., probability
    pdf_mean ./= num_onehoursegments

    return pdf_mean
end


text:
@testset "RNN: vanilla" begin
    # gradient test for RNNCell
    m = RNNCell(10, 5); x = rand(10, 4); h = rand(5, 4)
    @test gradcheck((W_ih, W_hh, b_ih, b_hh, x, h) -> sum(rnn_forward(W_ih, W_hh, b_ih, b_hh, x, h)),
                    m.W_ih, m.W_hh, m.b_ih, m.b_hh, x, h)
    # smoke tests for RNN
    m = RNN(10 => 5); x_seq = rand(10, 4, 2); h = init_hidden(m, 4)
    grad((m, x_seq, h) -> begin h_all, h = m(x_seq, h); sum(h) end, m, x_seq, h)
    grad((m, x_seq, h) -> begin h_all, h = m(x_seq, h); sum(h_all) end, m, x_seq, h)
end

@testset "RNN: LSTM" begin
    # gradient test for LSTMCell
    m = LSTMCell(10, 5); x = rand(10, 4); h = rand(5, 4); c = rand(5, 4)
    @test gradcheck((W_ih, W_hh, b_ih, b_hh, x, h, c) -> begin
                    h, c = lstm_forward(W_ih, W_hh, b_ih, b_hh, x, h, c)
                    sum(h)
                    end,
                    m.W_ih, m.W_hh, m.b_ih, m.b_hh, x, h, c)
    # smokes for LSTM
    m = LSTM(10 => 5); x_seq = rand(10, 4, 2); h, c = init_hidden(m, 4)
    grad((m, x_seq, h, c) ->
         begin h_all, h, c = m(x_seq, h, c); sum(h) end, m, x_seq, h, c)
    grad((m, x_seq, h, c) ->
         begin h_all, h, c = m(x_seq, h, c); sum(h_all) end, m, x_seq, h, c)
end

@testset "RNN: GRU" begin
    # gradient test for GRUCell
    m = GRUCell(10, 5); x = rand(10, 4); h = rand(5, 4)
    @test gradcheck((W_ih, W_hh, b_ih, b_hh, x, h) -> sum(gru_forward(W_ih, W_hh, b_ih, b_hh, x, h)),
                    m.W_ih, m.W_hh, m.b_ih, m.b_hh, x, h)
    # smoke tests for GRU
    m = GRU(10 => 5); x_seq = rand(10, 4, 2); h = init_hidden(m, 4)
    grad((m, x_seq, h) -> begin h_all, h = m(x_seq, h); sum(h) end, m, x_seq, h)
    grad((m, x_seq, h) -> begin h_all, h = m(x_seq, h); sum(h_all) end, m, x_seq, h)
end

ext: jl


text:
module MatrixMarket

using Compat.SparseArrays
using Compat.LinearAlgebra

export mmread, mmwrite

struct ParseError
    error :: String
end

_parseint(x) = parse(Int, x)

"""
### mmread(filename, infoonly::Bool=false, retcoord::Bool=false)

Read the contents of the Matrix Market file 'filename' into a matrix,
which will be either sparse or dense, depending on the Matrix Market format
indicated by 'coordinate' (coordinate sparse storage), or 'array' (dense
array storage).

If infoonly is true (default: false), only information on the size and
structure is returned from reading the header. The actual data for the
matrix elements are not parsed.

If retcoord is true (default: false), the rows, column and value vectors
are returned, if it is a sparse matrix, along with the header information.
"""
function mmread(filename, infoonly::Bool=false, retcoord::Bool=false)
    open(filename,"r") do mmfile
        # Read first line
        firstline = chomp(readline(mmfile))
        tokens = split(firstline)
        if length(tokens) != 5
            throw(ParseError(string("Not enough words on first line: ", firstline)))
        end
        if tokens[1] != "%%MatrixMarket"
            throw(ParseError(string("Expected start of header `%%MatrixMarket`, got `$(tokens[1])`")))
        end
        (head1, rep, field, symm) = map(lowercase, tokens[2:5])
        if head1 != "matrix"
            throw(ParseError("Unknown MatrixMarket data type: $head1 (only \"matrix\" is supported)"))
        end

        eltype = field == "real" ? Float64 :
                 field == "complex" ? ComplexF64 :
                 field == "integer" ? Int64 :
                 field == "pattern" ? Bool :
                 throw(ParseError("Unsupported field $field (only real and complex are supported)"))

        symlabel = symm == "general" ? identity :
                   symm == "symmetric" ? symmetric! :
                   symm == "hermitian" ? hermitian! :
                   symm == "skew-symmetric" ? skewsymmetric! :
                   throw(ParseError("Unknown matrix symmetry: $symm (only general, symmetric, skew-symmetric and hermitian are supported)"))

        # Skip all comments and empty lines
        ll   = readline(mmfile)
        while length(chomp(ll))==0 || (length(ll) > 0 && ll[1] == '%')
            ll = readline(mmfile)
        end
        # Read matrix dimensions (and number of entries) from first non-comment line
        dd = map(_parseint, split(ll))
        if length(dd) < (rep == "coordinate" ? 3 : 2)
            throw(ParseError(string("Could not read in matrix dimensions from line: ", ll)))
        end
        rows = dd[1]
        cols = dd[2]
        entries = (rep == "coordinate") ? dd[3] : (rows * cols)
        infoonly && return (rows, cols, entries, rep, field, symm)

        rep == "coordinate" ||
            return symlabel(reshape([parse(Float64, readline(mmfile)) for i in 1:entries],
                                    (rows,cols)))

        rr = Vector{Int}(undef, entries)
        cc = Vector{Int}(undef, entries)
        xx = Vector{eltype}(undef, entries)
        for i in 1:entries
            line = readline(mmfile)
            splits = find_splits(line, eltype == ComplexF64 ? 3 : (eltype == Bool ? 1 : 2))
            rr[i] = _parseint(line[1:splits[1]])
            cc[i] = _parseint(eltype == Bool
                              ? line[splits[1]:end]
                              : line[splits[1]:splits[2]])
            if eltype == ComplexF64
                real = parse(Float64, line[splits[2]:splits[3]])
                imag = parse(Float64, line[splits[3]:length(line)])
                xx[i] = ComplexF64(real, imag)
            elseif eltype == Bool
                xx[i] = true
            else
                xx[i] = parse(eltype, line[splits[2]:length(line)])
            end
        end
        (retcoord
         ? (rr, cc, xx, rows, cols, entries, rep, field, symm)
         : symlabel(sparse(rr, cc, xx, rows, cols)))
    end
end

function find_splits(s::String, num)
    splits = Vector{Int}(undef, num)
    cur = 1
    in_space = s[1] == '\t' || s[1] == ' '
    @inbounds for i in 1:length(s)
        if s[i] == '\t' || s[i] == ' '
            if !in_space
                in_space = true
                splits[cur] = i
                cur += 1
                cur > num && break
            end
        else
            in_space = false
        end
    end

    splits
end

# Hack to represent skew-symmetric matrix as an ordinary matrix with duplicated elements
function skewsymmetric!(M::AbstractMatrix)
    m,n = size(M)
    m == n || throw(DimensionMismatch())
    return M .- transpose(tril(M, -1))
end

function symmetric!(M::AbstractMatrix)
    m,n = size(M)
    m == n || throw(DimensionMismatch())
    if eltype(M) == Bool
        return M .| transpose(tril(M, -1))
    else
        return M .+ transpose(tril(M, -1))
    end
end

function hermitian!(M::AbstractMatrix)
    m,n = size(M)
    m == n || throw(DimensionMismatch())
    if eltype(M) == Bool
        return M .| conj(transpose(tril(M, -1)))
    else
        return M .+ conj(transpose(tril(M, -1)))
    end
end

"""
### mmwrite(filename, matrix::SparseMatrixCSC)

Write a sparse matrix to file 'filename'.
"""
function mmwrite(filename, matrix :: SparseMatrixCSC)
  open(filename, "w") do file
    elem = eltype(matrix) <: Bool ? "pattern" :
           eltype(matrix) <: Integer ?  "integer" :
           eltype(matrix) <: AbstractFloat ? "real" :
           eltype(matrix) <: Complex ? "complex" :
           error("Invalid matrix type")
      sym = issymmetric(matrix) ? "symmetric" :
            ishermitian(matrix) ? "hermitian" :
            "general"

      # write mm header
      write(file, "%%MatrixMarket matrix coordinate $elem $sym\n")

      # only use lower triangular part of symmetric and Hermitian matrices
      if issymmetric(matrix) || ishermitian(matrix)
          matrix = tril(matrix)
      end

      # write matrix size and number of nonzeros
      write(file, "$(size(matrix, 1)) $(size(matrix, 2)) $(nnz(matrix))\n")

      rows = rowvals(matrix)
      vals = nonzeros(matrix)
      for i in 1:size(matrix, 2)
          for j in nzrange(matrix, i)
              write(file, "$(rows[j]) $i")
              if elem == "pattern" # omit values on pattern matrices
              elseif elem == "complex"
                  write(file, " $(real(vals[j])) $(imag(vals[j]))")
              else
                  write(file, " $(vals[j])")
              end
              write(file, "\n")
          end
      end
  end
end

end # module


text:
@muladd function DiffEqBase.addsteps!(k,t,uprev,u,dt,f,p,cache::Vern6Cache,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if length(k) < 9 || always_calc_begin
    @unpack c1,c2,c3,c4,c5,c6,a21,a31,a32,a41,a43,a51,a53,a54,a61,a63,a64,a65,a71,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,a91,a94,a95,a96,a97,a98 = cache.tab
    @unpack k1,k2,k3,k4,k5,k6,k7,k8,k9,tmp = cache
    @.. tmp = uprev+dt*(a21*k1)
    f(k2,tmp,p,t+c1*dt)
    @.. tmp = uprev+dt*(a31*k1+a32*k2)
    f(k3,tmp,p,t+c2*dt)
    @.. tmp = uprev+dt*(a41*k1+a43*k3)
    f(k4,tmp,p,t+c3*dt)
    @.. tmp = uprev+dt*(a51*k1+a53*k3+a54*k4)
    f(k5,tmp,p,t+c4*dt)
    @.. tmp = uprev+dt*(a61*k1+a63*k3+a64*k4+a65*k5)
    f(k6,tmp,p,t+c5*dt)
    @.. tmp = uprev+dt*(a71*k1+a73*k3+a74*k4+a75*k5+a76*k6)
    f(k7,tmp,p,t+c6*dt)
    @.. tmp = uprev+dt*(a81*k1+a83*k3+a84*k4+a85*k5+a86*k6+a87*k7)
    f(k8,tmp,p,t+dt)
    @.. tmp = uprev+dt*(a91*k1+a94*k4+a95*k5+a96*k6+a97*k7+a98*k8)
    f(k9,tmp,p,t+dt)
    copyat_or_push!(k,1,k1)
    copyat_or_push!(k,2,k2)
    copyat_or_push!(k,3,k3)
    copyat_or_push!(k,4,k4)
    copyat_or_push!(k,5,k5)
    copyat_or_push!(k,6,k6)
    copyat_or_push!(k,7,k7)
    copyat_or_push!(k,8,k8)
    copyat_or_push!(k,9,k9)
  end
  if (allow_calc_end && length(k)< 12) || force_calc_end # Have not added the extra stages yet
    @unpack c10,a1001,a1004,a1005,a1006,a1007,a1008,a1009,c11,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,c12,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211 = cache.tab.extra
    @unpack tmp = cache
    rtmp = similar(cache.k1)
    uidx = eachindex(uprev)
    @.. tmp = uprev+dt*(a1001*k[1]+a1004*k[4]+a1005*k[5]+a1006*k[6]+a1007*k[7]+a1008*k[8]+a1009*k[9])
    f(rtmp,tmp,p,t+c10*dt); copyat_or_push!(k,10,rtmp)
    @.. tmp = uprev+dt*(a1101*k[1]+a1104*k[4]+a1105*k[5]+a1106*k[6]+a1107*k[7]+a1108*k[8]+a1109*k[9]+a1110*k[10])
    f(rtmp,tmp,p,t+c11*dt); copyat_or_push!(k,11,rtmp)
    @.. tmp = uprev+dt*(a1201*k[1]+a1204*k[4]+a1205*k[5]+a1206*k[6]+a1207*k[7]+a1208*k[8]+a1209*k[9]+a1210*k[10]+a1211*k[11])
    f(rtmp,tmp,p,t+c12*dt); copyat_or_push!(k,12,rtmp)
  end
  nothing
end

@muladd function DiffEqBase.addsteps!(k,t,uprev,u,dt,f,p,cache::Vern7Cache,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if length(k) < 10 || always_calc_begin
    @unpack c2,c3,c4,c5,c6,c7,c8,a021,a031,a032,a041,a043,a051,a053,a054,a061,a063,a064,a065,a071,a073,a074,a075,a076,a081,a083,a084,a085,a086,a087,a091,a093,a094,a095,a096,a097,a098,a101,a103,a104,a105,a106,a107 = cache.tab
    @unpack k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,tmp = cache
    f(k1,uprev,p,t)
    @.. tmp = uprev+dt*(a021*k1)
    f(k2,tmp,p,t+c2*dt)
    @.. tmp = uprev+dt*(a031*k1+a032*k2)
    f(k3,tmp,p,t+c3*dt)
    @.. tmp = uprev+dt*(a041*k1+a043*k3)
    f(k4,tmp,p,t+c4*dt)
    @.. tmp = uprev+dt*(a051*k1+a053*k3+a054*k4)
    f(k5,tmp,p,t+c5*dt)
    @.. tmp = uprev+dt*(a061*k1+a063*k3+a064*k4+a065*k5)
    f(k6,tmp,p,t+c6*dt)
    @.. tmp = uprev+dt*(a071*k1+a073*k3+a074*k4+a075*k5+a076*k6)
    f(k7,tmp,p,t+c7*dt)
    @.. tmp = uprev+dt*(a081*k1+a083*k3+a084*k4+a085*k5+a086*k6+a087*k7)
    f(k8,tmp,p,t+c8*dt)
    @.. tmp = uprev+dt*(a091*k1+a093*k3+a094*k4+a095*k5+a096*k6+a097*k7+a098*k8)
    f(k9,tmp,p,t+dt)
    @.. tmp = uprev+dt*(a101*k1+a103*k3+a104*k4+a105*k5+a106*k6+a107*k7)
    f(k10,tmp,p,t+dt)
    copyat_or_push!(k,1,k1)
    copyat_or_push!(k,2,k2)
    copyat_or_push!(k,3,k3)
    copyat_or_push!(k,4,k4)
    copyat_or_push!(k,5,k5)
    copyat_or_push!(k,6,k6)
    copyat_or_push!(k,7,k7)
    copyat_or_push!(k,8,k8)
    copyat_or_push!(k,9,k9)
    copyat_or_push!(k,10,k10)
  end
  if (allow_calc_end && length(k)< 16) || force_calc_end # Have not added the extra stages yet
    @unpack tmp = cache
    rtmp = similar(cache.k1)
    @unpack c11,a1101,a1104,a1105,a1106,a1107,a1108,a1109,c12,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1211,c13,a1301,a1304,a1305,a1306,a1307,a1308,a1309,a1311,a1312,c14,a1401,a1404,a1405,a1406,a1407,a1408,a1409,a1411,a1412,a1413,c15,a1501,a1504,a1505,a1506,a1507,a1508,a1509,a1511,a1512,a1513,c16,a1601,a1604,a1605,a1606,a1607,a1608,a1609,a1611,a1612,a1613 = cache.tab.extra
    @.. tmp = uprev+dt*(a1101*k[1]+a1104*k[4]+a1105*k[5]+a1106*k[6]+a1107*k[7]+a1108*k[8]+a1109*k[9])
    f(rtmp,tmp,p,t+c11*dt); copyat_or_push!(k,11,rtmp)
    @.. tmp = uprev+dt*(a1201*k[1]+a1204*k[4]+a1205*k[5]+a1206*k[6]+a1207*k[7]+a1208*k[8]+a1209*k[9]+a1211*k[11])
    f(rtmp,tmp,p,t+c12*dt); copyat_or_push!(k,12,rtmp)
    @.. tmp = uprev+dt*(a1301*k[1]+a1304*k[4]+a1305*k[5]+a1306*k[6]+a1307*k[7]+a1308*k[8]+a1309*k[9]+a1311*k[11]+a1312*k[12])
    f(rtmp,tmp,p,t+c13*dt); copyat_or_push!(k,13,rtmp)
    @.. tmp = uprev+dt*(a1401*k[1]+a1404*k[4]+a1405*k[5]+a1406*k[6]+a1407*k[7]+a1408*k[8]+a1409*k[9]+a1411*k[11]+a1412*k[12]+a1413*k[13])
    f(rtmp,tmp,p,t+c14*dt); copyat_or_push!(k,14,rtmp)
    @.. tmp=  uprev+dt*(a1501*k[1]+a1504*k[4]+a1505*k[5]+a1506*k[6]+a1507*k[7]+a1508*k[8]+a1509*k[9]+a1511*k[11]+a1512*k[12]+a1513*k[13])
    f(rtmp,tmp,p,t+c15*dt); copyat_or_push!(k,15,rtmp)
    @.. tmp = uprev+dt*(a1601*k[1]+a1604*k[4]+a1605*k[5]+a1606*k[6]+a1607*k[7]+a1608*k[8]+a1609*k[9]+a1611*k[11]+a1612*k[12]+a1613*k[13])
    f(rtmp,tmp,p,t+c16*dt); copyat_or_push!(k,16,rtmp)
  end
  nothing
end

@muladd function DiffEqBase.addsteps!(k,t,uprev,u,dt,f,p,cache::Vern8Cache,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if length(k) <13 || always_calc_begin
    @unpack c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0804,a0805,a0806,a0807,a0901,a0904,a0905,a0906,a0907,a0908,a1001,a1004,a1005,a1006,a1007,a1008,a1009,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211,a1301,a1304,a1305,a1306,a1307,a1308,a1309,a1310 = cache.tab
    @unpack k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,tmp = cache
    f(k1,uprev,p,t)
    @.. tmp = uprev+dt*(a0201*k1)
    f(k2,tmp,p,t+c2*dt)
    @.. tmp = uprev+dt*(a0301*k1+a0302*k2)
    f(k3,tmp,p,t+c3*dt)
    @.. tmp = uprev+dt*(a0401*k1+a0403*k3)
    f(k4,tmp,p,t+c4*dt)
    @.. tmp = uprev+dt*(a0501*k1+a0503*k3+a0504*k4)
    f(k5,tmp,p,t+c5*dt)
    @.. tmp = uprev+dt*(a0601*k1+a0604*k4+a0605*k5)
    f(k6,tmp,p,t+c6*dt)
    @.. tmp = uprev+dt*(a0701*k1+a0704*k4+a0705*k5+a0706*k6)
    f(k7,tmp,p,t+c7*dt)
    @.. tmp = uprev+dt*(a0801*k1+a0804*k4+a0805*k5+a0806*k6+a0807*k7)
    f(k8,tmp,p,t+c8*dt)
    @.. tmp = uprev+dt*(a0901*k1+a0904*k4+a0905*k5+a0906*k6+a0907*k7+a0908*k8)
    f(k9,tmp,p,t+c9*dt)
    @.. tmp = uprev+dt*(a1001*k1+a1004*k4+a1005*k5+a1006*k6+a1007*k7+a1008*k8+a1009*k9)
    f(k10,tmp,p,t+c10*dt)
    @.. tmp = uprev+dt*(a1101*k1+a1104*k4+a1105*k5+a1106*k6+a1107*k7+a1108*k8+a1109*k9+a1110*k10)
    f(k11,tmp,p,t+c11*dt)
    @.. tmp = uprev+dt*(a1201*k1+a1204*k4+a1205*k5+a1206*k6+a1207*k7+a1208*k8+a1209*k9+a1210*k10+a1211*k11)
    f(k12,tmp,p,t+dt)
    @.. tmp = uprev+dt*(a1301*k1+a1304*k4+a1305*k5+a1306*k6+a1307*k7+a1308*k8+a1309*k9+a1310*k10)
    f(k13,tmp,p,t+dt)
    copyat_or_push!(k,1,k1)
    copyat_or_push!(k,2,k2)
    copyat_or_push!(k,3,k3)
    copyat_or_push!(k,4,k4)
    copyat_or_push!(k,5,k5)
    copyat_or_push!(k,6,k6)
    copyat_or_push!(k,7,k7)
    copyat_or_push!(k,8,k8)
    copyat_or_push!(k,9,k9)
    copyat_or_push!(k,10,k10)
    copyat_or_push!(k,11,k11)
    copyat_or_push!(k,12,k12)
    copyat_or_push!(k,13,k13)
  end
  if (allow_calc_end && length(k)< 21) || force_calc_end # Have not added the extra stages yet
    rtmp = similar(cache.k1)
    @unpack c14,a1401,a1406,a1407,a1408,a1409,a1410,a1411,a1412,c15,a1501,a1506,a1507,a1508,a1509,a1510,a1511,a1512,a1514,c16,a1601,a1606,a1607,a1608,a1609,a1610,a1611,a1612,a1614,a1615,c17,a1701,a1706,a1707,a1708,a1709,a1710,a1711,a1712,a1714,a1715,a1716,c18,a1801,a1806,a1807,a1808,a1809,a1810,a1811,a1812,a1814,a1815,a1816,a1817,c19,a1901,a1906,a1907,a1908,a1909,a1910,a1911,a1912,a1914,a1915,a1916,a1917,c20,a2001,a2006,a2007,a2008,a2009,a2010,a2011,a2012,a2014,a2015,a2016,a2017,c21,a2101,a2106,a2107,a2108,a2109,a2110,a2111,a2112,a2114,a2115,a2116,a2117 = cache.tab.extra
    @unpack tmp = cache
    @.. tmp = uprev+dt*(a1401*k[1]+a1406*k[6]+a1407*k[7]+a1408*k[8]+a1409*k[9]+a1410*k[10]+a1411*k[11]+a1412*k[12])
    f(rtmp,tmp,p,t+c14*dt); copyat_or_push!(k,14,rtmp)
    @.. tmp = uprev+dt*(a1501*k[1]+a1506*k[6]+a1507*k[7]+a1508*k[8]+a1509*k[9]+a1510*k[10]+a1511*k[11]+a1512*k[12]+a1514*k[14])
    f(rtmp,tmp,p,t+c15*dt); copyat_or_push!(k,15,rtmp)
    @.. tmp = uprev+dt*(a1601*k[1]+a1606*k[6]+a1607*k[7]+a1608*k[8]+a1609*k[9]+a1610*k[10]+a1611*k[11]+a1612*k[12]+a1614*k[14]+a1615*k[15])
    f(rtmp,tmp,p,t+c16*dt); copyat_or_push!(k,16,rtmp)
    @.. tmp = uprev+dt*(a1701*k[1]+a1706*k[6]+a1707*k[7]+a1708*k[8]+a1709*k[9]+a1710*k[10]+a1711*k[11]+a1712*k[12]+a1714*k[14]+a1715*k[15]+a1716*k[16])
    f(rtmp,tmp,p,t+c17*dt); copyat_or_push!(k,17,rtmp)
    @.. tmp = uprev+dt*(a1801*k[1]+a1806*k[6]+a1807*k[7]+a1808*k[8]+a1809*k[9]+a1810*k[10]+a1811*k[11]+a1812*k[12]+a1814*k[14]+a1815*k[15]+a1816*k[16]+a1817*k[17])
    f(rtmp,tmp,p,t+c18*dt); copyat_or_push!(k,18,rtmp)
    @.. tmp = uprev+dt*(a1901*k[1]+a1906*k[6]+a1907*k[7]+a1908*k[8]+a1909*k[9]+a1910*k[10]+a1911*k[11]+a1912*k[12]+a1914*k[14]+a1915*k[15]+a1916*k[16]+a1917*k[17])
    f(rtmp,tmp,p,t+c19*dt); copyat_or_push!(k,19,rtmp)
    @.. tmp = uprev+dt*(a2001*k[1]+a2006*k[6]+a2007*k[7]+a2008*k[8]+a2009*k[9]+a2010*k[10]+a2011*k[11]+a2012*k[12]+a2014*k[14]+a2015*k[15]+a2016*k[16]+a2017*k[17])
    f(rtmp,tmp,p,t+c20*dt); copyat_or_push!(k,20,rtmp)
    @.. tmp = uprev+dt*(a2101*k[1]+a2106*k[6]+a2107*k[7]+a2108*k[8]+a2109*k[9]+a2110*k[10]+a2111*k[11]+a2112*k[12]+a2114*k[14]+a2115*k[15]+a2116*k[16]+a2117*k[17])
    f(rtmp,tmp,p,t+c21*dt); copyat_or_push!(k,21,rtmp)
  end
  nothing
end

@muladd function DiffEqBase.addsteps!(k,t,uprev,u,dt,f,p,cache::Vern9Cache,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if length(k) < 10 || always_calc_begin
    @unpack c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0806,a0807,a0901,a0906,a0907,a0908,a1001,a1006,a1007,a1008,a1009,a1101,a1106,a1107,a1108,a1109,a1110,a1201,a1206,a1207,a1208,a1209,a1210,a1211,a1301,a1306,a1307,a1308,a1309,a1310,a1311,a1312,a1401,a1406,a1407,a1408,a1409,a1410,a1411,a1412,a1413,a1501,a1506,a1507,a1508,a1509,a1510,a1511,a1512,a1513,a1514,a1601,a1606,a1607,a1608,a1609,a1610,a1611,a1612,a1613 = cache.tab
    @unpack k1,k2,k3,k4,k5,k6,k7,k8,k9,k10,k11,k12,k13,k14,k15,k16,tmp = cache
    uidx = eachindex(uprev)
    f(k1,uprev,p,t)
    @.. tmp = uprev+dt*(a0201*k1)
    f(k2,tmp,p,t+c1*dt)
    @.. tmp = uprev+dt*(a0301*k1+a0302*k2)
    f(k3,tmp,p,t+c2*dt)
    @.. tmp = uprev+dt*(a0401*k1+a0403*k3)
    f(k4,tmp,p,t+c3*dt)
    @.. tmp = uprev+dt*(a0501*k1+a0503*k3+a0504*k4)
    f(k5,tmp,p,t+c4*dt)
    @.. tmp = uprev+dt*(a0601*k1+a0604*k4+a0605*k5)
    f(k6,tmp,p,t+c5*dt)
    @.. tmp = uprev+dt*(a0701*k1+a0704*k4+a0705*k5+a0706*k6)
    f(k7,tmp,p,t+c6*dt)
    @.. tmp = uprev+dt*(a0801*k1+a0806*k6+a0807*k7)
    f(k8,tmp,p,t+c7*dt)
    @.. tmp = uprev+dt*(a0901*k1+a0906*k6+a0907*k7+a0908*k8)
    f(k9,tmp,p,t+c8*dt)
    @.. tmp = uprev+dt*(a1001*k1+a1006*k6+a1007*k7+a1008*k8+a1009*k9)
    f(k10,tmp,p,t+c9*dt)
    @.. tmp = uprev+dt*(a1101*k1+a1106*k6+a1107*k7+a1108*k8+a1109*k9+a1110*k10)
    f(k11,tmp,p,t+c10*dt)
    @.. tmp = uprev+dt*(a1201*k1+a1206*k6+a1207*k7+a1208*k8+a1209*k9+a1210*k10+a1211*k11)
    f(k12,tmp,p,t+c11*dt)
    @.. tmp = uprev+dt*(a1301*k1+a1306*k6+a1307*k7+a1308*k8+a1309*k9+a1310*k10+a1311*k11+a1312*k12)
    f(k13,tmp,p,t+c12*dt)
    @.. tmp = uprev+dt*(a1401*k1+a1406*k6+a1407*k7+a1408*k8+a1409*k9+a1410*k10+a1411*k11+a1412*k12+a1413*k13)
    f(k14,tmp,p,t+c13*dt)
    @.. tmp = uprev+dt*(a1501*k1+a1506*k6+a1507*k7+a1508*k8+a1509*k9+a1510*k10+a1511*k11+a1512*k12+a1513*k13+a1514*k14)
    f(k15,tmp,p,t+dt)
    @.. tmp = uprev+dt*(a1601*k1+a1606*k6+a1607*k7+a1608*k8+a1609*k9+a1610*k10+a1611*k11+a1612*k12+a1613*k13)
    f(k16,tmp,p,t+dt)
    copyat_or_push!(k,1,k1)
    copyat_or_push!(k,2,k8)
    copyat_or_push!(k,3,k9)
    copyat_or_push!(k,4,k10)
    copyat_or_push!(k,5,k11)
    copyat_or_push!(k,6,k12)
    copyat_or_push!(k,7,k13)
    copyat_or_push!(k,8,k14)
    copyat_or_push!(k,9,k15)
    copyat_or_push!(k,10,k16)
  end
  if (allow_calc_end && length(k)< 20) || force_calc_end # Have not added the extra stages yet
    rtmp = similar(cache.k1)
    uidx = eachindex(uprev)
    @unpack tmp = cache
    @unpack c17,a1701,a1708,a1709,a1710,a1711,a1712,a1713,a1714,a1715,c18,a1801,a1808,a1809,a1810,a1811,a1812,a1813,a1814,a1815,a1817,c19,a1901,a1908,a1909,a1910,a1911,a1912,a1913,a1914,a1915,a1917,a1918,c20,a2001,a2008,a2009,a2010,a2011,a2012,a2013,a2014,a2015,a2017,a2018,a2019,c21,a2101,a2108,a2109,a2110,a2111,a2112,a2113,a2114,a2115,a2117,a2118,a2119,a2120,c22,a2201,a2208,a2209,a2210,a2211,a2212,a2213,a2214,a2215,a2217,a2218,a2219,a2220,a2221,c23,a2301,a2308,a2309,a2310,a2311,a2312,a2313,a2314,a2315,a2317,a2318,a2319,a2320,a2321,c24,a2401,a2408,a2409,a2410,a2411,a2412,a2413,a2414,a2415,a2417,a2418,a2419,a2420,a2421,c25,a2501,a2508,a2509,a2510,a2511,a2512,a2513,a2514,a2515,a2517,a2518,a2519,a2520,a2521,c26,a2601,a2608,a2609,a2610,a2611,a2612,a2613,a2614,a2615,a2617,a2618,a2619,a2620,a2621 = cache.tab.extra
    @.. tmp = uprev+dt*(a1701*k[1]+a1708*k[2]+a1709*k[3]+a1710*k[4]+a1711*k[5]+a1712*k[6]+a1713*k[7]+a1714*k[8]+a1715*k[9])
    f(rtmp,tmp,p,t+c17*dt); copyat_or_push!(k,11,rtmp)
    @.. tmp = uprev+dt*(a1801*k[1]+a1808*k[2]+a1809*k[3]+a1810*k[4]+a1811*k[5]+a1812*k[6]+a1813*k[7]+a1814*k[8]+a1815*k[9]+a1817*k[11])
    f(rtmp,tmp,p,t+c18*dt); copyat_or_push!(k,12,rtmp)
    @.. tmp = uprev+dt*(a1901*k[1]+a1908*k[2]+a1909*k[3]+a1910*k[4]+a1911*k[5]+a1912*k[6]+a1913*k[7]+a1914*k[8]+a1915*k[9]+a1917*k[11]+a1918*k[12])
    f(rtmp,tmp,p,t+c19*dt); copyat_or_push!(k,13,rtmp)
    @.. tmp = uprev+dt*(a2001*k[1]+a2008*k[2]+a2009*k[3]+a2010*k[4]+a2011*k[5]+a2012*k[6]+a2013*k[7]+a2014*k[8]+a2015*k[9]+a2017*k[11]+a2018*k[12]+a2019*k[13])
    f(rtmp,tmp,p,t+c20*dt); copyat_or_push!(k,14,rtmp)
    @.. tmp = uprev+dt*(a2101*k[1]+a2108*k[2]+a2109*k[3]+a2110*k[4]+a2111*k[5]+a2112*k[6]+a2113*k[7]+a2114*k[8]+a2115*k[9]+a2117*k[11]+a2118*k[12]+a2119*k[13]+a2120*k[14])
    f(rtmp,tmp,p,t+c21*dt); copyat_or_push!(k,15,rtmp)
    @.. tmp = uprev+dt*(a2201*k[1]+a2208*k[2]+a2209*k[3]+a2210*k[4]+a2211*k[5]+a2212*k[6]+a2213*k[7]+a2214*k[8]+a2215*k[9]+a2217*k[11]+a2218*k[12]+a2219*k[13]+a2220*k[14]+a2221*k[15])
    f(rtmp,tmp,p,t+c22*dt); copyat_or_push!(k,16,rtmp)
    @.. tmp = uprev+dt*(a2301*k[1]+a2308*k[2]+a2309*k[3]+a2310*k[4]+a2311*k[5]+a2312*k[6]+a2313*k[7]+a2314*k[8]+a2315*k[9]+a2317*k[11]+a2318*k[12]+a2319*k[13]+a2320*k[14]+a2321*k[15])
    f(rtmp,tmp,p,t+c23*dt); copyat_or_push!(k,17,rtmp)
    @.. tmp  = uprev+dt*(a2401*k[1]+a2408*k[2]+a2409*k[3]+a2410*k[4]+a2411*k[5]+a2412*k[6]+a2413*k[7]+a2414*k[8]+a2415*k[9]+a2417*k[11]+a2418*k[12]+a2419*k[13]+a2420*k[14]+a2421*k[15])
    f(rtmp,tmp,p,t+c24*dt); copyat_or_push!(k,18,rtmp)
    @.. tmp = uprev+dt*(a2501*k[1]+a2508*k[2]+a2509*k[3]+a2510*k[4]+a2511*k[5]+a2512*k[6]+a2513*k[7]+a2514*k[8]+a2515*k[9]+a2517*k[11]+a2518*k[12]+a2519*k[13]+a2520*k[14]+a2521*k[15])
    f(rtmp,tmp,p,t+c25*dt); copyat_or_push!(k,19,rtmp)
    @.. tmp = uprev+dt*(a2601*k[1]+a2608*k[2]+a2609*k[3]+a2610*k[4]+a2611*k[5]+a2612*k[6]+a2613*k[7]+a2614*k[8]+a2615*k[9]+a2617*k[11]+a2618*k[12]+a2619*k[13]+a2620*k[14]+a2621*k[15])
    f(rtmp,tmp,p,t+c26*dt); copyat_or_push!(k,20,rtmp)
  end
  nothing
end

@muladd function DiffEqBase.addsteps!(k,t,uprev,u,dt,f,p,cache::Vern6ConstantCache,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if length(k) < 9 || always_calc_begin
    @unpack c1,c2,c3,c4,c5,c6,a21,a31,a32,a41,a43,a51,a53,a54,a61,a63,a64,a65,a71,a73,a74,a75,a76,a81,a83,a84,a85,a86,a87,a91,a94,a95,a96,a97,a98 = cache.tab
    copyat_or_push!(k,1,f(uprev,p,t))
    copyat_or_push!(k,2,f(uprev+dt*(a21*k[1]),p,t+c1*dt))
    copyat_or_push!(k,3,f(uprev+dt*(a31*k[1]+a32*k[2]),p,t+c2*dt))
    copyat_or_push!(k,4,f(uprev+dt*(a41*k[1]       +a43*k[3]),p,t+c3*dt))
    copyat_or_push!(k,5,f(uprev+dt*(a51*k[1]       +a53*k[3]+a54*k[4]),p,t+c4*dt))
    copyat_or_push!(k,6,f(uprev+dt*(a61*k[1]       +a63*k[3]+a64*k[4]+a65*k[5]),p,t+c5*dt))
    copyat_or_push!(k,7,f(uprev+dt*(a71*k[1]       +a73*k[3]+a74*k[4]+a75*k[5]+a76*k[6]),p,t+c6*dt))
    copyat_or_push!(k,8,f(uprev+dt*(a81*k[1]       +a83*k[3]+a84*k[4]+a85*k[5]+a86*k[6]+a87*k[7]),p,t+dt))
    copyat_or_push!(k,9,f(uprev+dt*(a91*k[1]+a94*k[4]+a95*k[5]+a96*k[6]+a97*k[7]+a98*k[8]),p,t+dt))
  end
  if (allow_calc_end && length(k)< 12) || force_calc_end # Have not added the extra stages yet
    @unpack c10,a1001,a1004,a1005,a1006,a1007,a1008,a1009,c11,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,c12,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211 = cache.tab.extra
    copyat_or_push!(k,10,f(uprev+dt*(a1001*k[1]+a1004*k[4]+a1005*k[5]+a1006*k[6]+a1007*k[7]+a1008*k[8]+a1009*k[9]),p,t+c10*dt))
    copyat_or_push!(k,11,f(uprev+dt*(a1101*k[1]+a1104*k[4]+a1105*k[5]+a1106*k[6]+a1107*k[7]+a1108*k[8]+a1109*k[9]+a1110*k[10]),p,t+c11*dt))
    copyat_or_push!(k,12,f(uprev+dt*(a1201*k[1]+a1204*k[4]+a1205*k[5]+a1206*k[6]+a1207*k[7]+a1208*k[8]+a1209*k[9]+a1210*k[10]+a1211*k[11]),p,t+c12*dt))
  end
  nothing
end

@muladd function DiffEqBase.addsteps!(k,t,uprev,u,dt,f,p,cache::Vern7ConstantCache,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if length(k) < 10 || always_calc_begin
    @unpack c2,c3,c4,c5,c6,c7,c8,a021,a031,a032,a041,a043,a051,a053,a054,a061,a063,a064,a065,a071,a073,a074,a075,a076,a081,a083,a084,a085,a086,a087,a091,a093,a094,a095,a096,a097,a098,a101,a103,a104,a105,a106,a107 = cache.tab
    copyat_or_push!(k,1,f(uprev,p,t))
    copyat_or_push!(k,2,f(uprev+dt*(a021*k[1]),p,t+c2*dt))
    copyat_or_push!(k,3,f(uprev+dt*(a031*k[1]+a032*k[2]),p,t+c3*dt))
    copyat_or_push!(k,4,f(uprev+dt*(a041*k[1]       +a043*k[3]),p,t+c4*dt))
    copyat_or_push!(k,5,f(uprev+dt*(a051*k[1]       +a053*k[3]+a054*k[4]),p,t+c5*dt))
    copyat_or_push!(k,6,f(uprev+dt*(a061*k[1]       +a063*k[3]+a064*k[4]+a065*k[5]),p,t+c6*dt))
    copyat_or_push!(k,7,f(uprev+dt*(a071*k[1]       +a073*k[3]+a074*k[4]+a075*k[5]+a076*k[6]),p,t+c7*dt))
    copyat_or_push!(k,8,f(uprev+dt*(a081*k[1]       +a083*k[3]+a084*k[4]+a085*k[5]+a086*k[6]+a087*k[7]),p,t+c8*dt))
    copyat_or_push!(k,9,f(uprev+dt*(a091*k[1]          +a093*k[3]+a094*k[4]+a095*k[5]+a096*k[6]+a097*k[7]+a098*k[8]),p,t+dt))
    copyat_or_push!(k,10,f(uprev+dt*(a101*k[1]          +a103*k[3]+a104*k[4]+a105*k[5]+a106*k[6]+a107*k[7]),p,t+dt))
  end
  if (allow_calc_end && length(k)< 16) || force_calc_end # Have not added the extra stages yet
    @unpack c11,a1101,a1104,a1105,a1106,a1107,a1108,a1109,c12,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1211,c13,a1301,a1304,a1305,a1306,a1307,a1308,a1309,a1311,a1312,c14,a1401,a1404,a1405,a1406,a1407,a1408,a1409,a1411,a1412,a1413,c15,a1501,a1504,a1505,a1506,a1507,a1508,a1509,a1511,a1512,a1513,c16,a1601,a1604,a1605,a1606,a1607,a1608,a1609,a1611,a1612,a1613 = cache.tab.extra
    copyat_or_push!(k,11,f(uprev+dt*(a1101*k[1]+a1104*k[4]+a1105*k[5]+a1106*k[6]+a1107*k[7]+a1108*k[8]+a1109*k[9]),p,t+c11*dt))
    copyat_or_push!(k,12,f(uprev+dt*(a1201*k[1]+a1204*k[4]+a1205*k[5]+a1206*k[6]+a1207*k[7]+a1208*k[8]+a1209*k[9]+a1211*k[11]),p,t+c12*dt))
    copyat_or_push!(k,13,f(uprev+dt*(a1301*k[1]+a1304*k[4]+a1305*k[5]+a1306*k[6]+a1307*k[7]+a1308*k[8]+a1309*k[9]+a1311*k[11]+a1312*k[12]),p,t+c13*dt))
    copyat_or_push!(k,14,f(uprev+dt*(a1401*k[1]+a1404*k[4]+a1405*k[5]+a1406*k[6]+a1407*k[7]+a1408*k[8]+a1409*k[9]+a1411*k[11]+a1412*k[12]+a1413*k[13]),p,t+c14*dt))
    copyat_or_push!(k,15,f(uprev+dt*(a1501*k[1]+a1504*k[4]+a1505*k[5]+a1506*k[6]+a1507*k[7]+a1508*k[8]+a1509*k[9]+a1511*k[11]+a1512*k[12]+a1513*k[13]),p,t+c15*dt))
    copyat_or_push!(k,16,f(uprev+dt*(a1601*k[1]+a1604*k[4]+a1605*k[5]+a1606*k[6]+a1607*k[7]+a1608*k[8]+a1609*k[9]+a1611*k[11]+a1612*k[12]+a1613*k[13]),p,t+c16*dt))
  end
  nothing
end

@muladd function DiffEqBase.addsteps!(k,t,uprev,u,dt,f,p,cache::Vern8ConstantCache,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if length(k) <13 || always_calc_begin
    @unpack c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0804,a0805,a0806,a0807,a0901,a0904,a0905,a0906,a0907,a0908,a1001,a1004,a1005,a1006,a1007,a1008,a1009,a1101,a1104,a1105,a1106,a1107,a1108,a1109,a1110,a1201,a1204,a1205,a1206,a1207,a1208,a1209,a1210,a1211,a1301,a1304,a1305,a1306,a1307,a1308,a1309,a1310 = cache.tab
    copyat_or_push!(k,1,f(uprev,p,t))
    copyat_or_push!(k,2,f(uprev+dt*(a0201*k[1]),p,t+c2*dt ))
    copyat_or_push!(k,3,f(uprev+dt*(a0301*k[1]+a0302*k[2]),p,t+c3*dt ))
    copyat_or_push!(k,4,f(uprev+dt*(a0401*k[1]       +a0403*k[3]),p,t+c4*dt ))
    copyat_or_push!(k,5,f(uprev+dt*(a0501*k[1]       +a0503*k[3]+a0504*k[4]),p,t+c5*dt ))
    copyat_or_push!(k,6,f(uprev+dt*(a0601*k[1]                +a0604*k[4]+a0605*k[5]),p,t+c6*dt ))
    copyat_or_push!(k,7,f(uprev+dt*(a0701*k[1]                +a0704*k[4]+a0705*k[5]+a0706*k[6]),p,t+c7*dt ))
    copyat_or_push!(k,8,f(uprev+dt*(a0801*k[1]                +a0804*k[4]+a0805*k[5]+a0806*k[6]+a0807*k[7]),p,t+c8*dt ))
    copyat_or_push!(k,9,f(uprev+dt*(a0901*k[1]                +a0904*k[4]+a0905*k[5]+a0906*k[6]+a0907*k[7]+a0908*k[8]),p,t+c9*dt ))
    copyat_or_push!(k,10,f(uprev+dt*(a1001*k[1]                +a1004*k[4]+a1005*k[5]+a1006*k[6]+a1007*k[7]+a1008*k[8]+a1009*k[9]),p,t+c10*dt))
    copyat_or_push!(k,11,f(uprev+dt*(a1101*k[1]                +a1104*k[4]+a1105*k[5]+a1106*k[6]+a1107*k[7]+a1108*k[8]+a1109*k[9]+a1110*k[10]),p,t+c11*dt))
    copyat_or_push!(k,12,f(uprev+dt*(a1201*k[1]                +a1204*k[4]+a1205*k[5]+a1206*k[6]+a1207*k[7]+a1208*k[8]+a1209*k[9]+a1210*k[10]+a1211*k[11]),p,t+    dt))
    copyat_or_push!(k,13,f(uprev+dt*(a1301*k[1]                +a1304*k[4]+a1305*k[5]+a1306*k[6]+a1307*k[7]+a1308*k[8]+a1309*k[9]+a1310*k[10]),p,t+    dt))
  end
  if (allow_calc_end && length(k)< 21) || force_calc_end # Have not added the extra stages yet
    @unpack c14,a1401,a1406,a1407,a1408,a1409,a1410,a1411,a1412,c15,a1501,a1506,a1507,a1508,a1509,a1510,a1511,a1512,a1514,c16,a1601,a1606,a1607,a1608,a1609,a1610,a1611,a1612,a1614,a1615,c17,a1701,a1706,a1707,a1708,a1709,a1710,a1711,a1712,a1714,a1715,a1716,c18,a1801,a1806,a1807,a1808,a1809,a1810,a1811,a1812,a1814,a1815,a1816,a1817,c19,a1901,a1906,a1907,a1908,a1909,a1910,a1911,a1912,a1914,a1915,a1916,a1917,c20,a2001,a2006,a2007,a2008,a2009,a2010,a2011,a2012,a2014,a2015,a2016,a2017,c21,a2101,a2106,a2107,a2108,a2109,a2110,a2111,a2112,a2114,a2115,a2116,a2117 = cache.tab.extra
    copyat_or_push!(k,14,f(uprev+dt*(a1401*k[1]+a1406*k[6]+a1407*k[7]+a1408*k[8]+a1409*k[9]+a1410*k[10]+a1411*k[11]+a1412*k[12]),p,t+c14*dt))
    copyat_or_push!(k,15,f(uprev+dt*(a1501*k[1]+a1506*k[6]+a1507*k[7]+a1508*k[8]+a1509*k[9]+a1510*k[10]+a1511*k[11]+a1512*k[12]+a1514*k[14]),p,t+c15*dt))
    copyat_or_push!(k,16,f(uprev+dt*(a1601*k[1]+a1606*k[6]+a1607*k[7]+a1608*k[8]+a1609*k[9]+a1610*k[10]+a1611*k[11]+a1612*k[12]+a1614*k[14]+a1615*k[15]),p,t+c16*dt))
    copyat_or_push!(k,17,f(uprev+dt*(a1701*k[1]+a1706*k[6]+a1707*k[7]+a1708*k[8]+a1709*k[9]+a1710*k[10]+a1711*k[11]+a1712*k[12]+a1714*k[14]+a1715*k[15]+a1716*k[16]),p,t+c17*dt))
    copyat_or_push!(k,18,f(uprev+dt*(a1801*k[1]+a1806*k[6]+a1807*k[7]+a1808*k[8]+a1809*k[9]+a1810*k[10]+a1811*k[11]+a1812*k[12]+a1814*k[14]+a1815*k[15]+a1816*k[16]+a1817*k[17]),p,t+c18*dt))
    copyat_or_push!(k,19,f(uprev+dt*(a1901*k[1]+a1906*k[6]+a1907*k[7]+a1908*k[8]+a1909*k[9]+a1910*k[10]+a1911*k[11]+a1912*k[12]+a1914*k[14]+a1915*k[15]+a1916*k[16]+a1917*k[17]),p,t+c19*dt))
    copyat_or_push!(k,20,f(uprev+dt*(a2001*k[1]+a2006*k[6]+a2007*k[7]+a2008*k[8]+a2009*k[9]+a2010*k[10]+a2011*k[11]+a2012*k[12]+a2014*k[14]+a2015*k[15]+a2016*k[16]+a2017*k[17]),p,t+c20*dt))
    copyat_or_push!(k,21,f(uprev+dt*(a2101*k[1]+a2106*k[6]+a2107*k[7]+a2108*k[8]+a2109*k[9]+a2110*k[10]+a2111*k[11]+a2112*k[12]+a2114*k[14]+a2115*k[15]+a2116*k[16]+a2117*k[17]),p,t+c21*dt))
  end
  nothing
end

@muladd function DiffEqBase.addsteps!(k,t,uprev,u,dt,f,p,cache::Vern9ConstantCache,always_calc_begin = false,allow_calc_end = true,force_calc_end = false)
  if length(k) < 10 || always_calc_begin
    @unpack c1,c2,c3,c4,c5,c6,c7,c8,c9,c10,c11,c12,c13,a0201,a0301,a0302,a0401,a0403,a0501,a0503,a0504,a0601,a0604,a0605,a0701,a0704,a0705,a0706,a0801,a0806,a0807,a0901,a0906,a0907,a0908,a1001,a1006,a1007,a1008,a1009,a1101,a1106,a1107,a1108,a1109,a1110,a1201,a1206,a1207,a1208,a1209,a1210,a1211,a1301,a1306,a1307,a1308,a1309,a1310,a1311,a1312,a1401,a1406,a1407,a1408,a1409,a1410,a1411,a1412,a1413,a1501,a1506,a1507,a1508,a1509,a1510,a1511,a1512,a1513,a1514,a1601,a1606,a1607,a1608,a1609,a1610,a1611,a1612,a1613 = cache.tab
    copyat_or_push!(k,1,f(uprev,p,t))
    copyat_or_push!(k,2,f(uprev+dt*(a0201*k[1]),p,t+c1*dt))
    copyat_or_push!(k,3,f(uprev+dt*(a0301*k[1]+a0302*k[2]),p,t+c2*dt))
    copyat_or_push!(k,4,f(uprev+dt*(a0401*k[1]       +a0403*k[3]),p,t+c3*dt))
    copyat_or_push!(k,5,f(uprev+dt*(a0501*k[1]       +a0503*k[3]+a0504*k[4]),p,t+c4*dt))
    copyat_or_push!(k,6,f(uprev+dt*(a0601*k[1]                +a0604*k[4]+a0605*k[5]),p,t+c5*dt))
    copyat_or_push!(k,7,f(uprev+dt*(a0701*k[1]                +a0704*k[4]+a0705*k[5]+a0706*k[6]),p,t+c6*dt))
    copyat_or_push!(k,2,f(uprev+dt*(a0801*k[1]                                  +a0806*k[6]+a0807*k[7]),p,t+c7*dt))
    copyat_or_push!(k,3,f(uprev+dt*(a0901*k[1]                                  +a0906*k[6]+a0907*k[7]+a0908*k[2]),p,t+c8*dt))
    copyat_or_push!(k,4,f(uprev+dt*(a1001*k[1]                                  +a1006*k[6]+a1007*k[7]+a1008*k[2]+a1009*k[3]),p,t+c9*dt))
    copyat_or_push!(k,5,f(uprev+dt*(a1101*k[1]                                  +a1106*k[6]+a1107*k[7]+a1108*k[2]+a1109*k[3]+a1110*k[4]),p,t+c10*dt))
    temp6 = recursivecopy(k[6])
    temp7 = recursivecopy(k[7])
    copyat_or_push!(k,6,f(uprev+dt*(a1201*k[1]                                  +a1206*temp6+a1207*temp7+a1208*k[2]+a1209*k[3]+a1210*k[4]+a1211*k[5]),p,t+c11*dt))
    copyat_or_push!(k,7,f(uprev+dt*(a1301*k[1]                                  +a1306*temp6+a1307*temp7+a1308*k[2]+a1309*k[3]+a1310*k[4]+a1311*k[5]+a1312*k[6]),p,t+c12*dt))
    copyat_or_push!(k,8,f(uprev+dt*(a1401*k[1]                                  +a1406*temp6+a1407*temp7+a1408*k[2]+a1409*k[3]+a1410*k[4]+a1411*k[5]+a1412*k[6]+a1413*k[7]),p,t+c13*dt))
    copyat_or_push!(k,9,f(uprev+dt*(a1501*k[1]                                  +a1506*temp6+a1507*temp7+a1508*k[2]+a1509*k[3]+a1510*k[4]+a1511*k[5]+a1512*k[6]+a1513*k[7]+a1514*k[8]),p,t+dt))
    copyat_or_push!(k,10,f(uprev+dt*(a1601*k[1]                                  +a1606*temp6+a1607*temp7+a1608*k[2]+a1609*k[3]+a1610*k[4]+a1611*k[5]+a1612*k[6]+a1613*k[7]),p,t+dt))
  end
  if (allow_calc_end && length(k)< 20) || force_calc_end # Have not added the extra stages yet
    @unpack c17,a1701,a1708,a1709,a1710,a1711,a1712,a1713,a1714,a1715,c18,a1801,a1808,a1809,a1810,a1811,a1812,a1813,a1814,a1815,a1817,c19,a1901,a1908,a1909,a1910,a1911,a1912,a1913,a1914,a1915,a1917,a1918,c20,a2001,a2008,a2009,a2010,a2011,a2012,a2013,a2014,a2015,a2017,a2018,a2019,c21,a2101,a2108,a2109,a2110,a2111,a2112,a2113,a2114,a2115,a2117,a2118,a2119,a2120,c22,a2201,a2208,a2209,a2210,a2211,a2212,a2213,a2214,a2215,a2217,a2218,a2219,a2220,a2221,c23,a2301,a2308,a2309,a2310,a2311,a2312,a2313,a2314,a2315,a2317,a2318,a2319,a2320,a2321,c24,a2401,a2408,a2409,a2410,a2411,a2412,a2413,a2414,a2415,a2417,a2418,a2419,a2420,a2421,c25,a2501,a2508,a2509,a2510,a2511,a2512,a2513,a2514,a2515,a2517,a2518,a2519,a2520,a2521,c26,a2601,a2608,a2609,a2610,a2611,a2612,a2613,a2614,a2615,a2617,a2618,a2619,a2620,a2621 = cache.tab.extra
    copyat_or_push!(k,11,f(uprev+dt*(a1701*k[1]+a1708*k[2]+a1709*k[3]+a1710*k[4]+a1711*k[5]+a1712*k[6]+a1713*k[7]+a1714*k[8]+a1715*k[9]),p,t+c17*dt))
    copyat_or_push!(k,12,f(uprev+dt*(a1801*k[1]+a1808*k[2]+a1809*k[3]+a1810*k[4]+a1811*k[5]+a1812*k[6]+a1813*k[7]+a1814*k[8]+a1815*k[9]+a1817*k[11]),p,t+c18*dt))
    copyat_or_push!(k,13,f(uprev+dt*(a1901*k[1]+a1908*k[2]+a1909*k[3]+a1910*k[4]+a1911*k[5]+a1912*k[6]+a1913*k[7]+a1914*k[8]+a1915*k[9]+a1917*k[11]+a1918*k[12]),p,t+c19*dt))
    copyat_or_push!(k,14,f(uprev+dt*(a2001*k[1]+a2008*k[2]+a2009*k[3]+a2010*k[4]+a2011*k[5]+a2012*k[6]+a2013*k[7]+a2014*k[8]+a2015*k[9]+a2017*k[11]+a2018*k[12]+a2019*k[13]),p,t+c20*dt))
    copyat_or_push!(k,15,f(uprev+dt*(a2101*k[1]+a2108*k[2]+a2109*k[3]+a2110*k[4]+a2111*k[5]+a2112*k[6]+a2113*k[7]+a2114*k[8]+a2115*k[9]+a2117*k[11]+a2118*k[12]+a2119*k[13]+a2120*k[14]),p,t+c21*dt))
    copyat_or_push!(k,16,f(uprev+dt*(a2201*k[1]+a2208*k[2]+a2209*k[3]+a2210*k[4]+a2211*k[5]+a2212*k[6]+a2213*k[7]+a2214*k[8]+a2215*k[9]+a2217*k[11]+a2218*k[12]+a2219*k[13]+a2220*k[14]+a2221*k[15]),p,t+c22*dt))
    copyat_or_push!(k,17,f(uprev+dt*(a2301*k[1]+a2308*k[2]+a2309*k[3]+a2310*k[4]+a2311*k[5]+a2312*k[6]+a2313*k[7]+a2314*k[8]+a2315*k[9]+a2317*k[11]+a2318*k[12]+a2319*k[13]+a2320*k[14]+a2321*k[15]),p,t+c23*dt))
    copyat_or_push!(k,18,f(uprev+dt*(a2401*k[1]+a2408*k[2]+a2409*k[3]+a2410*k[4]+a2411*k[5]+a2412*k[6]+a2413*k[7]+a2414*k[8]+a2415*k[9]+a2417*k[11]+a2418*k[12]+a2419*k[13]+a2420*k[14]+a2421*k[15]),p,t+c24*dt))
    copyat_or_push!(k,19,f(uprev+dt*(a2501*k[1]+a2508*k[2]+a2509*k[3]+a2510*k[4]+a2511*k[5]+a2512*k[6]+a2513*k[7]+a2514*k[8]+a2515*k[9]+a2517*k[11]+a2518*k[12]+a2519*k[13]+a2520*k[14]+a2521*k[15]),p,t+c25*dt))
    copyat_or_push!(k,20,f(uprev+dt*(a2601*k[1]+a2608*k[2]+a2609*k[3]+a2610*k[4]+a2611*k[5]+a2612*k[6]+a2613*k[7]+a2614*k[8]+a2615*k[9]+a2617*k[11]+a2618*k[12]+a2619*k[13]+a2620*k[14]+a2621*k[15]),p,t+c26*dt))
  end
  nothing
end


text:
mutable struct CholeskyFactorization{Tv, Ti} <: AbstractLUFactorization{Tv,Ti}
    A::Union{ExtendableSparseMatrix{Tv,Ti},Nothing}
    fact::Union{SuiteSparse.CHOLMOD.Factor{Tv},Nothing}
    phash::UInt64
    A64
end

"""
CholeskyFactorization(;valuetype=Float64, indextype=Int64)
CholeskyFactorization(matrix)

Default Cholesky factorization via cholmod.
"""
CholeskyFactorization(;valuetype::Type=Float64, indextype::Type=Int64)=CholeskyFactorization{valuetype,indextype}(nothing,nothing,0,nothing)

function update!(cholfact::CholeskyFactorization)
    A=cholfact.A
    flush!(A)
    if A.phash!=cholfact.phash
        cholfact.A64=Symmetric(SparseMatrixCSC{Float64,Int64}(A.cscmatrix))
        cholfact.fact=cholesky(cholfact.A64)
        cholfact.phash=A.phash
    else
        cholfact.A64.data.nzval.=A.cscmatrix.nzval
        cholfact.fact=cholesky!(cholfact.fact,cholfact.A64)
    end
    cholfact
end

LinearAlgebra.ldiv!(fact::CholeskyFactorization, v)=fact.fact\v
LinearAlgebra.ldiv!(u,fact::CholeskyFactorization, v)=u.=fact.fact\v


text:
#========================================================================================#
#	Sample startup.jl file
#
# Insert this file into your directory ~/.julia/config to have it start automatically
# whenever you start Julia.
#
# Author: Niall Palfreyman, 06/02/2022
#========================================================================================#

# This function moves to the indicated directory and includes Ingolstadt, allowing you to
# start working on your current laboratory. Remember to replace ??? by your own
# directory specification:
function ingo()
	cd("C:/Users/???/Ingolstadt")
	include("src/Ingolstadt.jl")
end


text:
module Grids

export Center, Face
export AbstractTopology, Periodic, Bounded, Flat, FullyConnected, LeftConnected, RightConnected, topology

export AbstractGrid, AbstractUnderlyingGrid, halo_size, total_size
export AbstractRectilinearGrid, RectilinearGrid
export XRegRectilinearGrid, YRegRectilinearGrid, ZRegRectilinearGrid, HRegRectilinearGrid, RegRectilinearGrid
export AbstractCurvilinearGrid, AbstractHorizontallyCurvilinearGrid
export LatitudeLongitudeGrid, XRegLatLonGrid, YRegLatLonGrid, ZRegLatLonGrid
export ConformalCubedSphereFaceGrid, ConformalCubedSphereGrid
export node, xnode, ynode, znode, xnodes, ynodes, znodes, nodes
export offset_data, new_data
export on_architecture

using CUDA
using Adapt
using OffsetArrays

using Oceananigans
using Oceananigans.Architectures

import Base: size, length, eltype, show
import Oceananigans.Architectures: architecture

#####
##### Abstract types
#####

"""
    Center

A type describing the location at the center of a grid cell.
"""
struct Center end

"""
	Face

A type describing the location at the face of a grid cell.
"""
struct Face end

"""
    AbstractTopology

Abstract supertype for grid topologies.
"""
abstract type AbstractTopology end

"""
    Periodic

Grid topology for periodic dimensions.
"""
struct Periodic <: AbstractTopology end

"""
    Bounded

Grid topology for bounded dimensions, e.g., wall-bounded dimensions.
"""
struct Bounded <: AbstractTopology end

"""
    Flat

Grid topology for flat dimensions, generally with one grid point, along which the solution
is uniform and does not vary.
"""
struct Flat <: AbstractTopology end

"""
    FullyConnected

Grid topology for dimensions that are connected to other models or domains.
"""
struct FullyConnected <: AbstractTopology end

"""
    LeftConnected

Grid topology for dimensions that are connected to other models or domains only on the left (the other direction is bounded)
"""
struct LeftConnected <: AbstractTopology end

"""
    RightConnected

Grid topology for dimensions that are connected to other models or domains only on the right (the other direction is bounded)
"""
struct RightConnected <: AbstractTopology end

"""
    AbstractGrid{FT, TX, TY, TZ}

Abstract supertype for grids with elements of type `FT` and topology `{TX, TY, TZ}`.
"""
abstract type AbstractGrid{FT, TX, TY, TZ, Arch} end

"""
    AbstractUnderlyingGrid{FT, TX, TY, TZ}

Abstract supertype for "primary" grids (as opposed to grids with immersed boundaries)
with elements of type `FT` and topology `{TX, TY, TZ}`.
"""
abstract type AbstractUnderlyingGrid{FT, TX, TY, TZ, Arch} <: AbstractGrid{FT, TX, TY, TZ, Arch} end

"""
    AbstractRectilinearGrid{FT, TX, TY, TZ}

Abstract supertype for rectilinear grids with elements of type `FT` and topology `{TX, TY, TZ}`.
"""
abstract type AbstractRectilinearGrid{FT, TX, TY, TZ, Arch} <: AbstractUnderlyingGrid{FT, TX, TY, TZ, Arch} end

"""
    AbstractCurvilinearGrid{FT, TX, TY, TZ}

Abstract supertype for curvilinear grids with elements of type `FT` and topology `{TX, TY, TZ}`.
"""
abstract type AbstractCurvilinearGrid{FT, TX, TY, TZ, Arch} <: AbstractUnderlyingGrid{FT, TX, TY, TZ, Arch} end

"""
    AbstractHorizontallyCurvilinearGrid{FT, TX, TY, TZ}

Abstract supertype for horizontally-curvilinear grids with elements of type `FT` and topology `{TX, TY, TZ}`.
"""
abstract type AbstractHorizontallyCurvilinearGrid{FT, TX, TY, TZ, Arch} <: AbstractCurvilinearGrid{FT, TX, TY, TZ, Arch} end

isrectilinear(grid) = false

include("grid_utils.jl")
include("zeros.jl")
include("new_data.jl")
include("inactive_node.jl")
include("automatic_halo_sizing.jl")
include("input_validation.jl")
include("grid_generation.jl")
include("rectilinear_grid.jl")
include("conformal_cubed_sphere_face_grid.jl")
include("latitude_longitude_grid.jl")

end # module


text:
using Test, YaoArrayRegister, YaoBase

@testset "select" begin
    reg = product_state(4, 6; nbatch = 2)
    # println(focus!(reg, [1,3]))
    r1 = select!(focus!(copy(reg), [2, 3]), 0b11) |> relax!(to_nactive = 2)
    r2 = select(focus!(copy(reg), [2, 3]), 0b11) |> relax!(to_nactive = 2)
    r3 = copy(reg) |> focus!(2, 3) |> select!(0b11) |> relax!(to_nactive = 2)

    @test r1' * r1 ≈ ones(2)
    @test r1 ≈ r2
    @test r3 ≈ r2
end

@testset "measure and resetto/remove" begin
    reg = rand_state(4)
    res = measure!(YaoBase.ResetTo(0), reg, (4,))
    @test isnormalized(reg)
    result = measure(reg; nshots = 10)
    @test all(result .< 8)
    @test ndims(res) == 0

    reg = rand_state(6) |> focus!(1, 4, 3)
    reg0 = copy(reg)
    res = measure!(YaoBase.RemoveMeasured(), reg)
    @test nqubits(reg) == 3
    select(reg0, res)
    @test select(reg0, res) |> normalize! ≈ reg
    @test ndims(res) == 0

    r = rand_state(10)
    r1 = copy(r) |> focus!(1, 4, 3)
    res = measure!(YaoBase.RemoveMeasured(), r, (1, 4, 3))
    r2 = select(r1, res)
    r2 = relax!(r2, (); to_nactive = nqubits(r2))
    @test normalize!(r2) ≈ r
    @test ndims(res) == 0

    reg = rand_state(6, nbatch = 5) |> focus!((1:5)...)
    res = measure!(YaoBase.ResetTo(0), reg, 1)
    @test nactive(reg) == 5
    @test ndims(res) == 1
end

@testset "fix measure kwargs error" begin
    r = rand_state(10)
    @test length(measure(r; nshots = 10)) == 10
    @test_throws MethodError measure!(r; nshots = 10)
    @test_throws MethodError measure!(YaoBase.RemoveMeasured(), r; nshots = 10)
end

@testset "fix measure output type error" begin
    res = measure(rand_state(1; nbatch = 10))
    @test res isa Matrix{BitStr64{1}}
end

@testset "fix measure kwargs error" begin
    r = rand_state(10)
    @test length(measure(r; nshots = 10)) == 10
    @test_throws MethodError measure!(r; nshots = 10)
    @test_throws MethodError measure!(YaoBase.RemoveMeasured(), r; nshots = 10)
end


text:
import YAML
using Twitter
using Printf

SECRET_KEYS_FILENAME = "secret_keys.yaml"

if !ispath(SECRET_KEYS_FILENAME)
    println(@sprintf "%s file not found! Please make sure to create the appropriate file containing your TWitter Developer keys." SECRET_KEYS_FILENAME)

else
    secrets = YAML.load_file(SECRET_KEYS_FILENAME)

    twitterauth(secrets["consumer_key"],
        secrets["consumer_token"],
        secrets["oauth_token"],
        secrets["oauth_secret"])

    user_favorites_list = get_favorites_list(screen_name = secrets["twitter_handle"], count=200)
    number_likes_for_deletion = length(user_favorites_list)

    println(@sprintf "%d favorites/likes found from account %s... Would you like to delete these favorites? (Y/N)" number_likes_for_deletion secrets["twitter_handle"])
    response = readline()

    if (response == "Y")
        println(@sprintf """'%s' response detected. Likes/Favorites will now be deleted...""" response)
        for favorite in user_favorites_list
            favorite_id = favorite.id
            post_favorites_destroy(id=favorite_id)
        end
        println("The tweets have now been deleted. Exiting...")
    else
        println(@sprintf """ '%s' response detected. No tweets will be deleted. Exiting...""" response)
    end
end


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

module Printf
using Base: Grisu, GMP
export @printf, @sprintf

### printf formatter generation ###
const SmallFloatingPoint = Union{Float64,Float32,Float16}
const SmallNumber = Union{SmallFloatingPoint,Base.BitInteger}

function gen(s::AbstractString)
    args = []
    blk = Expr(:block, :(local neg, pt, len, exp, do_out, args))
    for x in parse(s)
        if isa(x,AbstractString)
            push!(blk.args, :(write(out, $(length(x)==1 ? x[1] : x))))
        else
            c = lowercase(x[end])
            f = c=='f' ? gen_f :
                c=='e' ? gen_e :
                c=='a' ? gen_a :
                c=='g' ? gen_g :
                c=='c' ? gen_c :
                c=='s' ? gen_s :
                c=='p' ? gen_p :
                         gen_d
            arg, ex = f(x...)
            push!(args, arg)
            push!(blk.args, ex)
        end
    end
    push!(blk.args, :nothing)
    return args, blk
end

### printf format string parsing ###

function parse(s::AbstractString)
    # parse format string in to stings and format tuples
    list = []
    i = j = start(s)
    while !done(s,j)
        c, k = next(s,j)
        if c == '%'
            isempty(s[i:j-1]) || push!(list, s[i:j-1])
            flags, width, precision, conversion, k = parse1(s,k)
            '\'' in flags && error("printf format flag ' not yet supported")
            conversion == 'n'    && error("printf feature %n not supported")
            push!(list, conversion == '%' ? "%" : (flags,width,precision,conversion))
            i = j = k
        else
            j = k
        end
    end
    isempty(s[i:end]) || push!(list, s[i:end])
    # coalesce adjacent strings
    i = 1
    while i < length(list)
        if isa(list[i],AbstractString)
            for outer j = i+1:length(list)
                if !isa(list[j],AbstractString)
                    j -= 1
                    break
                end
                list[i] *= list[j]
            end
            deleteat!(list,i+1:j)
        end
        i += 1
    end
    return list
end

## parse a single printf specifier ##

# printf specifiers:
#   %                       # start
#   (\d+\$)?                # arg (not supported)
#   [\-\+#0' ]*             # flags
#   (\d+)?                  # width
#   (\.\d*)?                # precision
#   (h|hh|l|ll|L|j|t|z|q)?  # modifier (ignored)
#   [diouxXeEfFgGaAcCsSp%]  # conversion

next_or_die(s::AbstractString, k) = !done(s,k) ? next(s,k) :
    throw(ArgumentError("invalid printf format string: $(repr(s))"))

function parse1(s::AbstractString, k::Integer)
    j = k
    width = 0
    precision = -1
    c, k = next_or_die(s,k)
    # handle %%
    if c == '%'
        return "", width, precision, c, k
    end
    # parse flags
    while c in "#0- + '"
        c, k = next_or_die(s,k)
    end
    flags = String(s[j:k-2])
    # parse width
    while '0' <= c <= '9'
        width = 10*width + c-'0'
        c, k = next_or_die(s,k)
    end
    # parse precision
    if c == '.'
        c, k = next_or_die(s,k)
        if '0' <= c <= '9'
            precision = 0
            while '0' <= c <= '9'
                precision = 10*precision + c-'0'
                c, k = next_or_die(s,k)
            end
        end
    end
    # parse length modifer (ignored)
    if c == 'h' || c == 'l'
        prev = c
        c, k = next_or_die(s,k)
        if c == prev
            c, k = next_or_die(s,k)
        end
    elseif c in "Ljqtz"
        c, k = next_or_die(s,k)
    end
    # validate conversion
    if !(c in "diouxXDOUeEfFgGaAcCsSpn")
        throw(ArgumentError("invalid printf format string: $(repr(s))"))
    end
    # TODO: warn about silly flag/conversion combinations
    flags, width, precision, c, k
end

### printf formatter generation ###

function special_handler(flags::String, width::Int)
    @gensym x
    blk = Expr(:block)
    pad = '-' in flags ? rpad : lpad
    pos = '+' in flags ? "+" :
          ' ' in flags ? " " : ""
    abn = quote
        isnan($x) ? $(pad("NaN", width)) :
         $x < 0   ? $(pad("-Inf", width)) :
                    $(pad("$(pos)Inf", width))
    end
    ex = :(isfinite($x) ? $blk : write(out, $abn))
    x, ex, blk
end

function pad(m::Int, n, c::Char)
    if m <= 1
        :($n > 0 && write(out,$c))
    else
        @gensym i
        quote
            $i = $n
            while $i > 0
                write(out,$c)
                $i -= 1
            end
        end
    end
end

function dynamic_pad(m, val, c::Char)
    @gensym i
    quote
        if $m <= 1
            $val > 0 && write(out,$c)
        else
            $i = $val
            while $i > 0
                write(out,$c)
                $i -= 1
            end
        end
    end
end

function print_fixed(out, precision, pt, ndigits, trailingzeros=true)
    pdigits = pointer(DIGITS)
    if pt <= 0
        # 0.0dddd0
        write(out, '0')
        write(out, '.')
        precision += pt
        while pt < 0
            write(out, '0')
            pt += 1
        end
        unsafe_write(out, pdigits, ndigits)
        precision -= ndigits
    elseif ndigits <= pt
        # dddd000.000000
        unsafe_write(out, pdigits, ndigits)
        while ndigits < pt
            write(out, '0')
            ndigits += 1
        end
        if trailingzeros
            write(out, '.')
        end
    else # 0 < pt < ndigits
        # dd.dd0000
        ndigits -= pt
        unsafe_write(out, pdigits, pt)
        write(out, '.')
        unsafe_write(out, pdigits+pt, ndigits)
        precision -= ndigits
    end
    if trailingzeros
        while precision > 0
            write(out, '0')
            precision -= 1
        end
    end
end

function print_exp_e(out, exp::Integer)
    write(out, exp < 0 ? '-' : '+')
    exp = abs(exp)
    d = div(exp,100)
    if d > 0
        if d >= 10
            print(out, exp)
            return
        end
        write(out, Char('0'+d))
    end
    exp = rem(exp,100)
    write(out, Char('0'+div(exp,10)))
    write(out, Char('0'+rem(exp,10)))
end

function print_exp_a(out, exp::Integer)
    write(out, exp < 0 ? '-' : '+')
    exp = abs(exp)
    print(out, exp)
end

function gen_d(flags::String, width::Int, precision::Int, c::Char)
    # print integer:
    #  [dDiu]: print decimal digits
    #  [o]:    print octal digits
    #  [x]:    print hex digits, lowercase
    #  [X]:    print hex digits, uppercase
    #
    # flags:
    #  (#): prefix hex with 0x/0X; octal leads with 0
    #  (0): pad left with zeros
    #  (-): left justify
    #  ( ): precede non-negative values with " "
    #  (+): precede non-negative values with "+"
    #
    x, ex, blk = special_handler(flags,width)
    # interpret the number
    prefix = ""
    if lowercase(c)=='o'
        fn = '#' in flags ? :decode_0ct : :decode_oct
    elseif c=='x'
        '#' in flags && (prefix = "0x")
        fn = :decode_hex
    elseif c=='X'
        '#' in flags && (prefix = "0X")
        fn = :decode_HEX
    else
        fn = :decode_dec
    end
    push!(blk.args, :((do_out, args) = $fn(out, $x, $flags, $width, $precision, $c)))
    ifblk = Expr(:if, :do_out, Expr(:block))
    push!(blk.args, ifblk)
    blk = ifblk.args[2]
    push!(blk.args, :((len, pt, neg) = args))
    # calculate padding
    width -= length(prefix)
    space_pad = width > max(1,precision) && '-' in flags ||
                precision < 0 && width > 1 && !('0' in flags) ||
                precision >= 0 && width > precision
    padding = nothing
    if precision < 1; precision = 1; end
    if space_pad
        if '+' in flags || ' ' in flags
            width -= 1
            if width > precision
                padding = :($width-(pt > $precision ? pt : $precision))
            end
        else
            if width > precision
                padding = :($width-neg-(pt > $precision ? pt : $precision))
            end
        end
    end
    # print space padding
    if padding !== nothing && !('-' in flags)
        push!(blk.args, pad(width-precision, padding, ' '))
    end
    # print sign
    '+' in flags ? push!(blk.args, :(write(out, neg ? '-' : '+'))) :
    ' ' in flags ? push!(blk.args, :(write(out, neg ? '-' : ' '))) :
                   push!(blk.args, :(neg && write(out, '-')))
    # print prefix
    for ch in prefix
        push!(blk.args, :(write(out, $ch)))
    end
    # print zero padding & leading zeros
    if space_pad && precision > 1
        push!(blk.args, pad(precision-1, :($precision-pt), '0'))
    elseif !space_pad && width > 1
        zeros = '+' in flags || ' ' in flags ? :($(width-1)-pt) : :($width-neg-pt)
        push!(blk.args, pad(width-1, zeros, '0'))
    end
    # print integer
    push!(blk.args, :(unsafe_write(out, pointer(DIGITS), pt)))
    # print padding
    if padding !== nothing && '-' in flags
        push!(blk.args, pad(width-precision, padding, ' '))
    end
    # return arg, expr
    :(($x)::Real), ex
end

function gen_f(flags::String, width::Int, precision::Int, c::Char)
    # print to fixed trailing precision
    #  [fF]: the only choice
    #
    # flags
    #  (#): always print a decimal point
    #  (0): pad left with zeros
    #  (-): left justify
    #  ( ): precede non-negative values with " "
    #  (+): precede non-negative values with "+"
    #
    x, ex, blk = special_handler(flags,width)
    # interpret the number
    if precision < 0; precision = 6; end
    push!(blk.args, :((do_out, args) = fix_dec(out, $x, $flags, $width, $precision, $c)))
    ifblk = Expr(:if, :do_out, Expr(:block))
    push!(blk.args, ifblk)
    blk = ifblk.args[2]
    push!(blk.args, :((len, pt, neg) = args))
    # calculate padding
    padding = nothing
    if precision > 0 || '#' in flags
        width -= precision+1
    end
    if '+' in flags || ' ' in flags
        width -= 1
        if width > 1
            padding = :($width-(pt > 0 ? pt : 1))
        end
    else
        if width > 1
            padding = :($width-(pt > 0 ? pt : 1)-neg)
        end
    end
    # print space padding
    if padding !== nothing && !('-' in flags) && !('0' in flags)
        push!(blk.args, pad(width-1, padding, ' '))
    end
    # print sign
    '+' in flags ? push!(blk.args, :(write(out, neg ? '-' : '+'))) :
    ' ' in flags ? push!(blk.args, :(write(out, neg ? '-' : ' '))) :
                   push!(blk.args, :(neg && write(out, '-')))
    # print zero padding
    if padding !== nothing && !('-' in flags) && '0' in flags
        push!(blk.args, pad(width-1, padding, '0'))
    end
    # print digits
    if precision > 0
        push!(blk.args, :(print_fixed(out,$precision,pt,len)))
    else
        push!(blk.args, :(unsafe_write(out, pointer(DIGITS), len)))
        push!(blk.args, :(while pt >= (len+=1) write(out,'0') end))
        '#' in flags && push!(blk.args, :(write(out, '.')))
    end
    # print space padding
    if padding !== nothing && '-' in flags
        push!(blk.args, pad(width-1, padding, ' '))
    end
    # return arg, expr
    :(($x)::Real), ex
end

function gen_e(flags::String, width::Int, precision::Int, c::Char, inside_g::Bool=false)
    # print float in scientific form:
    #  [e]: use 'e' to introduce exponent
    #  [E]: use 'E' to introduce exponent
    #
    # flags:
    #  (#): always print a decimal point
    #  (0): pad left with zeros
    #  (-): left justify
    #  ( ): precede non-negative values with " "
    #  (+): precede non-negative values with "+"
    #
    x, ex, blk = if inside_g
        @gensym x
        blk = Expr(:block)
        x, blk, blk
    else
        special_handler(flags,width)
    end
    # interpret the number
    if precision < 0; precision = 6; end
    ndigits = min(precision+1,length(DIGITS)-1)
    push!(blk.args, :((do_out, args) = ini_dec(out,$x,$ndigits, $flags, $width, $precision, $c)))
    ifblk = Expr(:if, :do_out, Expr(:block))
    push!(blk.args, ifblk)
    blk = ifblk.args[2]
    push!(blk.args, :((len, pt, neg) = args))
    push!(blk.args, :(exp = pt-1))
    expmark = isupper(c) ? "E" : "e"
    if precision==0 && '#' in flags
        expmark = string(".",expmark)
    end
    # calculate padding
    padding = nothing
    width -= precision+length(expmark)+(precision>0)+4
    # 4 = leading + expsign + 2 exp digits
    if '+' in flags || ' ' in flags
        width -= 1 # for the sign indicator
        if width > 0
            padding = quote
                padn=$width
                if (exp<=-100)|(100<=exp)
                    if isa($x,SmallNumber)
                        padn -= 1
                    else
                        padn -= Base.ndigits0z(exp) - 2
                    end
                end
                padn
            end
        end
    else
        if width > 0
            padding = quote
                padn=$width-neg
                if (exp<=-100)|(100<=exp)
                    if isa($x,SmallNumber)
                        padn -= 1
                    else
                        padn -= Base.ndigits0z(exp) - 2
                    end
                end
                padn
            end
        end
    end
    # print space padding
    if padding !== nothing && !('-' in flags) && !('0' in flags)
        push!(blk.args, pad(width, padding, ' '))
    end
    # print sign
    '+' in flags ? push!(blk.args, :(write(out, neg ? '-' : '+'))) :
    ' ' in flags ? push!(blk.args, :(write(out, neg ? '-' : ' '))) :
                   push!(blk.args, :(neg && write(out, '-')))
    # print zero padding
    if padding !== nothing && !('-' in flags) && '0' in flags
        push!(blk.args, pad(width, padding, '0'))
    end
    # print digits
    push!(blk.args, :(write(out, DIGITS[1])))
    if precision > 0
        if inside_g && !('#' in flags)
            push!(blk.args, :(endidx = $ndigits;
                              while endidx > 1 && DIGITS[endidx] == UInt8('0')
                                  endidx -= 1
                              end;
                              if endidx > 1
                                  write(out, '.')
                                  unsafe_write(out, pointer(DIGITS)+1, endidx-1)
                              end
                              ))
        else
            push!(blk.args, :(write(out, '.')))
            push!(blk.args, :(unsafe_write(out, pointer(DIGITS)+1, $(ndigits-1))))
            if ndigits < precision+1
                n = precision+1-ndigits
                push!(blk.args, pad(n, n, '0'))
            end
        end
    end
    for ch in expmark
        push!(blk.args, :(write(out, $ch)))
    end
    push!(blk.args, :(print_exp_e(out, exp)))
    # print space padding
    if padding !== nothing && '-' in flags
        push!(blk.args, pad(width, padding, ' '))
    end
    # return arg, expr
    :(($x)::Real), ex
end

function gen_a(flags::String, width::Int, precision::Int, c::Char)
    # print float in hexadecimal format
    #  [a]: lowercase hex float, e.g. -0x1.cfp-2
    #  [A]: uppercase hex float, e.g. -0X1.CFP-2
    #
    # flags:
    #  (#): always print a decimal point
    #  (0): pad left with zeros
    #  (-): left justify
    #  ( ): precede non-negative values with " "
    #  (+): precede non-negative values with "+"
    #
    x, ex, blk = special_handler(flags,width)
    if c == 'A'
        hexmark, expmark = "0X", "P"
        fn = :ini_HEX
    else
        hexmark, expmark = "0x", "p"
        fn = :ini_hex
    end
    # if no precision, print max non-zero
    if precision < 0
        push!(blk.args, :((do_out, args) = $fn(out,$x, $flags, $width, $precision, $c)))
    else
        ndigits = min(precision+1,length(DIGITS)-1)
        push!(blk.args, :((do_out, args) = $fn(out,$x,$ndigits, $flags, $width, $precision, $c)))
    end
    ifblk = Expr(:if, :do_out, Expr(:block))
    push!(blk.args, ifblk)
    blk = ifblk.args[2]
    push!(blk.args, :((len, exp, neg) = args))
    if precision==0 && '#' in flags
        expmark = string(".",expmark)
    end
    # calculate padding
    padding = nothing
    if precision > 0
        width -= precision+length(hexmark)+length(expmark)+4
        # 4 = leading + expsign + 1 exp digit + decimal
    else
        width -= length(hexmark)+length(expmark)+3+(precision<0 && '#' in flags)
        # 3 = leading + expsign + 1 exp digit
    end
    if '+' in flags || ' ' in flags
        width -= 1 # for the sign indicator
        if width > 0
            padding = :($(width+1) - Base.ndigits(exp))
        end
    else
        if width > 0
            padding = :($(width+1) - neg - Base.ndigits(exp))
        end
    end
    if precision < 0 && width > 0
        if '#' in flags
            padding = :($padding - (len-1))
        else
            padding = :($padding - (len>1 ? len : 0))
        end
    end
    # print space padding
    if padding !== nothing && !('-' in flags) && !('0' in flags)
        push!(blk.args, pad(width, padding, ' '))
    end
    # print sign
    '+' in flags ? push!(blk.args, :(write(out, neg ? '-' : '+'))) :
    ' ' in flags ? push!(blk.args, :(write(out, neg ? '-' : ' '))) :
                    push!(blk.args, :(neg && write(out, '-')))
    # hex prefix
    for ch in hexmark
        push!(blk.args, :(write(out, $ch)))
    end
    # print zero padding
    if padding !== nothing && !('-' in flags) && '0' in flags
        push!(blk.args, pad(width, padding, '0'))
    end
    # print digits
    push!(blk.args, :(write(out, DIGITS[1])))
    if precision > 0
        push!(blk.args, :(write(out, '.')))
        push!(blk.args, :(unsafe_write(out, pointer(DIGITS)+1, $(ndigits-1))))
        if ndigits < precision+1
            n = precision+1-ndigits
            push!(blk.args, pad(n, n, '0'))
        end
    elseif precision < 0
        ifvpblk = Expr(:if, :(len > 1), Expr(:block))
        vpblk = ifvpblk.args[2]
        if '#' in flags
            push!(blk.args, :(write(out, '.')))
        else
            push!(vpblk.args, :(write(out, '.')))
        end
        push!(vpblk.args, :(unsafe_write(out, pointer(DIGITS)+1, len-1)))
        push!(blk.args, ifvpblk)
    end
    for ch in expmark
        push!(blk.args, :(write(out, $ch)))
    end
    push!(blk.args, :(print_exp_a(out, exp)))
    # print space padding
    if padding !== nothing && '-' in flags
        push!(blk.args, pad(width, padding, ' '))
    end
    # return arg, expr
    :(($x)::Real), ex
end

function gen_c(flags::String, width::Int, precision::Int, c::Char)
    # print a character:
    #  [cC]: both the same for us (Unicode)
    #
    # flags:
    #  (0): pad left with zeros
    #  (-): left justify
    #
    @gensym x
    blk = Expr(:block, :($x = Char($x)))
    if width > 1 && !('-' in flags)
        p = '0' in flags ? '0' : ' '
        push!(blk.args, pad(width-1, :($width-charwidth($x)), p))
    end
    push!(blk.args, :(write(out, $x)))
    if width > 1 && '-' in flags
        push!(blk.args, pad(width-1, :($width-charwidth($x)), ' '))
    end
    :(($x)::Integer), blk
end

function _limit(s, prec)
    prec >= sizeof(s) && return s
    p = prevind(s, prec+1)
    n = nextind(s, p)-1
    s[1:(prec>=n ? n : prevind(s,p))]
end

function gen_s(flags::String, width::Int, precision::Int, c::Char)
    # print a string:
    #  [sS]: both the same for us (Unicode)
    #
    # flags:
    #  (-): left justify
    #
    @gensym x
    blk = Expr(:block)
    if width > 0
        if !('#' in flags)
            push!(blk.args, :($x = string($x)))
        else
            push!(blk.args, :($x = repr($x)))
        end
        if precision!=-1
            push!(blk.args, :($x = _limit($x, $precision)))
        end
        if !('-' in flags)
            push!(blk.args, pad(width, :($width-strwidth($x)), ' '))
        end
        push!(blk.args, :(write(out, $x)))
        if '-' in flags
            push!(blk.args, pad(width, :($width-strwidth($x)), ' '))
        end
    else
        if precision!=-1
            push!(blk.args, :(io = IOBuffer()))
        else
            push!(blk.args, :(io = out))
        end
        if !('#' in flags)
            push!(blk.args, :(print(io, $x)))
        else
            push!(blk.args, :(show(io, $x)))
        end
        if precision!=-1
            push!(blk.args, :(write(out, _limit(String(take!(io)), $precision))))
        end
    end
    :(($x)::Any), blk
end

# TODO: faster pointer printing.

function gen_p(flags::String, width::Int, precision::Int, c::Char)
    # print pointer:
    #  [p]: the only option
    #
    @gensym x
    blk = Expr(:block)
    ptrwidth = Sys.WORD_SIZE>>2
    width -= ptrwidth+2
    if width > 0 && !('-' in flags)
        push!(blk.args, pad(width, width, ' '))
    end
    push!(blk.args, :(write(out, '0')))
    push!(blk.args, :(write(out, 'x')))
    push!(blk.args, :(write(out, String(hex(unsigned($x), $ptrwidth)))))
    if width > 0 && '-' in flags
        push!(blk.args, pad(width, width, ' '))
    end
    :(($x)::Ptr), blk
end

function gen_g(flags::String, width::Int, precision::Int, c::Char)
    # print to fixed trailing precision
    #  [g]: lower case e on scientific
    #  [G]: Upper case e on scientific
    #
    # flags
    #  (#): always print a decimal point
    #  (0): pad left with zeros
    #  (-): left justify
    #  ( ): precede non-negative values with " "
    #  (+): precede non-negative values with "+"
    #
    x, ex, blk = special_handler(flags,width)
    if precision < 0; precision = 6; end
    ndigits = min(precision+1,length(DIGITS)-1)
    # See if anyone else wants to handle it
    push!(blk.args, :((do_out, args) = ini_dec(out,$x,$ndigits, $flags, $width, $precision, $c)))
    ifblk = Expr(:if, :do_out, Expr(:block))
    push!(blk.args, ifblk)
    blk = ifblk.args[2]
    push!(blk.args, :((len, pt, neg) = args))
    push!(blk.args, :(exp = pt-1))
    push!(blk.args, :(do_f = $precision > exp >= -4)) # Should we interpret like %f or %e?
    feblk = Expr(:if, :do_f, Expr(:block), Expr(:block))
    push!(blk.args, feblk)
    fblk = feblk.args[2]
    eblk = feblk.args[3]

    ### %f branch
    # Follow the same logic as gen_f() but more work has to be deferred until runtime
    # because precision is unknown until then.
    push!(fblk.args, :(fprec = $precision - (exp+1)))
    push!(fblk.args, :((do_out, args) = fix_dec(out, $x, $flags, $width, fprec, $c - 1)))
    fifblk = Expr(:if, :do_out, Expr(:block))
    push!(fblk.args, fifblk)
    blk = fifblk.args[2]
    push!(blk.args, :((len, pt, neg) = args))
    push!(blk.args, :(padding = nothing))
    push!(blk.args, :(width = $width))
    # need to compute value before left-padding since trailing zeros are elided
    push!(blk.args, :(tmpout = IOBuffer()))
    push!(blk.args, :(print_fixed(tmpout,fprec,pt,len,$('#' in flags))))
    push!(blk.args, :(tmpstr = String(take!(tmpout))))
    push!(blk.args, :(width -= length(tmpstr)))
    if '+' in flags || ' ' in flags
        push!(blk.args, :(width -= 1))
    else
        push!(blk.args, :(if neg width -= 1; end))
    end
    push!(blk.args, :(if width >= 1 padding = width; end))
    # print space padding
    if !('-' in flags) && !('0' in flags)
        padexpr = dynamic_pad(:width, :padding, ' ')
        push!(blk.args, :(if padding !== nothing
                          $padexpr; end))
    end
    # print sign
    '+' in flags ? push!(blk.args, :(write(out, neg ? '-' : '+'))) :
    ' ' in flags ? push!(blk.args, :(write(out, neg ? '-' : ' '))) :
                   push!(blk.args, :(neg && write(out, '-')))
    # print zero padding
    if !('-' in flags) && '0' in flags
        padexpr = dynamic_pad(:width, :padding, '0')
        push!(blk.args, :(if padding !== nothing
                          $padexpr; end))
    end
    # finally print value
    push!(blk.args, :(write(out,tmpstr)))
    # print space padding
    if '-' in flags
        padexpr = dynamic_pad(:width, :padding, ' ')
        push!(blk.args, :(if padding !== nothing
                          $padexpr; end))
    end

    ### %e branch
    # Here we can do all the work at macro expansion time
    var, eex = gen_e(flags, width, precision-1, c, true)
    push!(eblk.args, :($(var.args[1]) = $x))
    push!(eblk.args, eex)

    :(($x)::Real), ex
end

### core unsigned integer decoding functions ###

macro handle_zero(ex)
    quote
        if $(esc(ex)) == 0
            DIGITS[1] = '0'
            return Int32(1), Int32(1), $(esc(:neg))
        end
    end
end

decode_oct(out, d, flags::String, width::Int, precision::Int, c::Char) = (true, decode_oct(d))
decode_0ct(out, d, flags::String, width::Int, precision::Int, c::Char) = (true, decode_0ct(d))
decode_dec(out, d, flags::String, width::Int, precision::Int, c::Char) = (true, decode_dec(d))
decode_hex(out, d, flags::String, width::Int, precision::Int, c::Char) = (true, decode_hex(d))
decode_HEX(out, d, flags::String, width::Int, precision::Int, c::Char) = (true, decode_HEX(d))
fix_dec(out, d, flags::String, width::Int, precision::Int, c::Char) = (true, fix_dec(d, precision))
ini_dec(out, d, ndigits::Int, flags::String, width::Int, precision::Int, c::Char) = (true, ini_dec(d, ndigits))
ini_hex(out, d, ndigits::Int, flags::String, width::Int, precision::Int, c::Char) = (true, ini_hex(d, ndigits))
ini_HEX(out, d, ndigits::Int, flags::String, width::Int, precision::Int, c::Char) = (true, ini_HEX(d, ndigits))
ini_hex(out, d, flags::String, width::Int, precision::Int, c::Char) = (true, ini_hex(d))
ini_HEX(out, d, flags::String, width::Int, precision::Int, c::Char) = (true, ini_HEX(d))

# fallbacks for Real types without explicit decode_* implementation
decode_oct(d::Real) = decode_oct(Integer(d))
decode_0ct(d::Real) = decode_0ct(Integer(d))
decode_dec(d::Real) = decode_dec(Integer(d))
decode_hex(d::Real) = decode_hex(Integer(d))
decode_HEX(d::Real) = decode_HEX(Integer(d))

handlenegative(d::Unsigned) = (false, d)
function handlenegative(d::Integer)
    if d < 0
        return true, unsigned(oftype(d,-d))
    else
        return false, unsigned(d)
    end
end

function decode_oct(d::Integer)
    neg, x = handlenegative(d)
    @handle_zero x
    pt = i = div((sizeof(x)<<3)-leading_zeros(x)+2,3)
    while i > 0
        DIGITS[i] = '0'+(x&0x7)
        x >>= 3
        i -= 1
    end
    return Int32(pt), Int32(pt), neg
end

function decode_0ct(d::Integer)
    neg, x = handlenegative(d)
    # doesn't need special handling for zero
    pt = i = div((sizeof(x)<<3)-leading_zeros(x)+5,3)
    while i > 0
        DIGITS[i] = '0'+(x&0x7)
        x >>= 3
        i -= 1
    end
    return Int32(pt), Int32(pt), neg
end

function decode_dec(d::Integer)
    neg, x = handlenegative(d)
    @handle_zero x
    pt = i = Base.ndigits0z(x)
    while i > 0
        DIGITS[i] = '0'+rem(x,10)
        x = div(x,10)
        i -= 1
    end
    return Int32(pt), Int32(pt), neg
end

function decode_hex(d::Integer, symbols::Array{UInt8,1})
    neg, x = handlenegative(d)
    @handle_zero x
    pt = i = (sizeof(x)<<1)-(leading_zeros(x)>>2)
    while i > 0
        DIGITS[i] = symbols[(x&0xf)+1]
        x >>= 4
        i -= 1
    end
    return Int32(pt), Int32(pt), neg
end

const hex_symbols = b"0123456789abcdef"
const HEX_symbols = b"0123456789ABCDEF"

decode_hex(x::Integer) = decode_hex(x,hex_symbols)
decode_HEX(x::Integer) = decode_hex(x,HEX_symbols)

function decode(b::Int, x::BigInt)
    neg = x.size < 0
    pt = Base.ndigits(x, abs(b))
    length(DIGITS) < pt+1 && resize!(DIGITS, pt+1)
    neg && (x.size = -x.size)
    GMP.MPZ.get_str!(DIGITS, b, x)
    neg && (x.size = -x.size)
    return Int32(pt), Int32(pt), neg
end
decode_oct(x::BigInt) = decode(8, x)
decode_dec(x::BigInt) = decode(10, x)
decode_hex(x::BigInt) = decode(16, x)
decode_HEX(x::BigInt) = decode(-16, x)

function decode_0ct(x::BigInt)
    neg = x.size < 0
    DIGITS[1] = '0'
    if x.size == 0
        return Int32(1), Int32(1), neg
    end
    pt = Base.ndigits0z(x, 8) + 1
    length(DIGITS) < pt+1 && resize!(DIGITS, pt+1)
    neg && (x.size = -x.size)
    p = convert(Ptr{UInt8}, DIGITS) + 1
    GMP.MPZ.get_str!(p, 8, x)
    neg && (x.size = -x.size)
    return neg, Int32(pt), Int32(pt)
end

### decoding functions directly used by printf generated code ###

# decode_*(x)=> fixed precision, to 0th place, filled out
# fix_*(x,n) => fixed precision, to nth place, not filled out
# ini_*(x,n) => n initial digits, filled out

# alternate versions:
#   *_0ct(x,n) => ensure that the first octal digits is zero
#   *_HEX(x,n) => use uppercase digits for hexadecimal

# - returns (len, point, neg)
# - implies len = point
#

function decode_dec(x::SmallFloatingPoint)
    if x == 0.0
        DIGITS[1] = '0'
        return (Int32(1), Int32(1), false)
    end
    len,pt,neg = grisu(x,Grisu.FIXED,0)
    if len == 0
        DIGITS[1] = '0'
        return (Int32(1), Int32(1), false)
    else
        for i = len+1:pt
            DIGITS[i] = '0'
        end
    end
    return Int32(len), Int32(pt), neg
end
# TODO: implement decode_oct, decode_0ct, decode_hex, decode_HEX for SmallFloatingPoint

## fix decoding functions ##
#
# - returns (neg, point, len)
# - if len less than point, trailing zeros implied
#

# fallback for Real types without explicit fix_dec implementation
fix_dec(x::Real, n::Int) = fix_dec(float(x),n)

fix_dec(x::Integer, n::Int) = decode_dec(x)

function fix_dec(x::SmallFloatingPoint, n::Int)
    if n > length(DIGITS)-1; n = length(DIGITS)-1; end
    len,pt,neg = grisu(x,Grisu.FIXED,n)
    if len == 0
        DIGITS[1] = '0'
        return (Int32(1), Int32(1), neg)
    end
    return Int32(len), Int32(pt), neg
end

## ini decoding functions ##
#
# - returns (neg, point, len)
# - implies len = n (requested digits)
#

# fallback for Real types without explicit fix_dec implementation
ini_dec(x::Real, n::Int) = ini_dec(float(x),n)

function ini_dec(d::Integer, n::Int)
    neg, x = handlenegative(d)
    k = ndigits(x)
    if k <= n
        pt = k
        for i = k:-1:1
            DIGITS[i] = '0'+rem(x,10)
            x = div(x,10)
        end
        for i = k+1:n
            DIGITS[i] = '0'
        end
    else
        p = Base.powers_of_ten[k-n+1]
        r = rem(x,p)
        if r >= (p>>1)
            x += p
            if x >= Base.powers_of_ten[k+1]
                p *= 10
                k += 1
            end
        end
        pt = k
        x = div(x,p)
        for i = n:-1:1
            DIGITS[i] = '0'+rem(x,10)
            x = div(x,10)
        end
    end
    return n, pt, neg
end

function ini_dec(x::SmallFloatingPoint, n::Int)
    if x == 0.0
        ccall(:memset, Ptr{Void}, (Ptr{Void}, Cint, Csize_t), DIGITS, '0', n)
        return Int32(1), Int32(1), signbit(x)
    else
        len,pt,neg = grisu(x,Grisu.PRECISION,n)
    end
    return Int32(len), Int32(pt), neg
end

function ini_dec(x::BigInt, n::Int)
    if x.size == 0
        ccall(:memset, Ptr{Void}, (Ptr{Void}, Cint, Csize_t), DIGITS, '0', n)
        return Int32(1), Int32(1), false
    end
    d = Base.ndigits0z(x)
    if d <= n
        info = decode_dec(x)
        d == n && return info
        p = convert(Ptr{Void}, DIGITS) + info[2]
        ccall(:memset, Ptr{Void}, (Ptr{Void}, Cint, Csize_t), p, '0', n - info[2])
        return info
    end
    return (n, d, decode_dec(round(BigInt,x/big(10)^(d-n)))[3])
end

ini_hex(x::Real, n::Int) = ini_hex(x,n,hex_symbols)
ini_HEX(x::Real, n::Int) = ini_hex(x,n,HEX_symbols)

ini_hex(x::Real) = ini_hex(x,hex_symbols)
ini_HEX(x::Real) = ini_hex(x,HEX_symbols)

ini_hex(x::Real, n::Int, symbols::Array{UInt8,1}) = ini_hex(float(x), n, symbols)
ini_hex(x::Real, symbols::Array{UInt8,1}) = ini_hex(float(x), symbols)

function ini_hex(x::SmallFloatingPoint, n::Int, symbols::Array{UInt8,1})
    x = Float64(x)
    if x == 0.0
        ccall(:memset, Ptr{Void}, (Ptr{Void}, Cint, Csize_t), DIGITS, '0', n)
        return Int32(1), Int32(0), signbit(x)
    else
        s, p = frexp(x)
        sigbits = 4*min(n-1,13)
        s = 0.25*round(ldexp(s,1+sigbits))
        # ensure last 2 exponent bits either 01 or 10
        u = (reinterpret(UInt64,s) & 0x003f_ffff_ffff_ffff) >> (52-sigbits)
        if n > 14
            ccall(:memset, Ptr{Void}, (Ptr{Void}, Cint, Csize_t), DIGITS, '0', n)
        end
        i = (sizeof(u)<<1)-(leading_zeros(u)>>2)
        while i > 0
            DIGITS[i] = symbols[(u&0xf)+1]
            u >>= 4
            i -= 1
        end
        # pt is the binary exponent
        return Int32(n), Int32(p-1), x < 0.0
    end
end

function ini_hex(x::SmallFloatingPoint, symbols::Array{UInt8,1})
    x = Float64(x)
    if x == 0.0
        ccall(:memset, Ptr{Void}, (Ptr{Void}, Cint, Csize_t), DIGITS, '0', 1)
        return Int32(1), Int32(0), signbit(x)
    else
        s, p = frexp(x)
        s *= 2.0
        u = (reinterpret(UInt64,s) & 0x001f_ffff_ffff_ffff)
        t = (trailing_zeros(u) >> 2)
        u >>= (t<<2)
        n = 14-t
        for i = n:-1:1
            DIGITS[i] = symbols[(u&0xf)+1]
            u >>= 4
        end
        # pt is the binary exponent
        return Int32(n), Int32(p-1), x < 0.0
    end
end

function ini_hex(x::Integer)
    len,pt,neg = decode_hex(x)
    pt = (len-1)<<2
    len,pt,neg
end
function ini_HEX(x::Integer)
    len,pt,neg = decode_HEX(x)
    pt = (len-1)<<2
    len,pt,neg
end

# not implemented
ini_hex(x::Integer,ndigits::Int) = throw(MethodError(ini_hex,(x,ndigits)))

#BigFloat
fix_dec(out, d::BigFloat, flags::String, width::Int, precision::Int, c::Char) = bigfloat_printf(out, d, flags, width, precision, c)
ini_dec(out, d::BigFloat, ndigits::Int, flags::String, width::Int, precision::Int, c::Char) = bigfloat_printf(out, d, flags, width, precision, c)
ini_hex(out, d::BigFloat, ndigits::Int, flags::String, width::Int, precision::Int, c::Char) = bigfloat_printf(out, d, flags, width, precision, c)
ini_HEX(out, d::BigFloat, ndigits::Int, flags::String, width::Int, precision::Int, c::Char) = bigfloat_printf(out, d, flags, width, precision, c)
ini_hex(out, d::BigFloat, flags::String, width::Int, precision::Int, c::Char) = bigfloat_printf(out, d, flags, width, precision, c)
ini_HEX(out, d::BigFloat, flags::String, width::Int, precision::Int, c::Char) = bigfloat_printf(out, d, flags, width, precision, c)
function bigfloat_printf(out, d::BigFloat, flags::String, width::Int, precision::Int, c::Char)
    fmt_len = sizeof(flags)+4
    if width > 0
        fmt_len += ndigits(width)
    end
    if precision >= 0
        fmt_len += ndigits(precision)+1
    end
    fmt = IOBuffer(fmt_len)
    write(fmt, '%')
    write(fmt, flags)
    if width > 0
        print(fmt, width)
    end
    if precision == 0
        write(fmt, '.')
        write(fmt, '0')
    elseif precision > 0
        write(fmt, '.')
        print(fmt, precision)
    end
    write(fmt, 'R')
    write(fmt, c)
    write(fmt, UInt8(0))
    printf_fmt = take!(fmt)
    @assert length(printf_fmt) == fmt_len
    bufsiz = length(DIGITS)
    lng = ccall((:mpfr_snprintf,:libmpfr), Int32,
                (Ptr{UInt8}, Culong, Ptr{UInt8}, Ref{BigFloat}...),
                DIGITS, bufsiz, printf_fmt, d)
    lng > 0 || error("invalid printf formatting for BigFloat")
    unsafe_write(out, pointer(DIGITS), min(lng, bufsiz-1))
    return (false, ())
end

### external printf interface ###

is_str_expr(ex) =
    isa(ex,Expr) && (ex.head == :string || (ex.head == :macrocall && isa(ex.args[1],Symbol) &&
    endswith(string(ex.args[1]),"str")))

function _printf(macroname, io, fmt, args)
    isa(fmt, AbstractString) || throw(ArgumentError("$macroname: format must be a plain static string (no interpolation or prefix)"))
    sym_args, blk = gen(fmt)

    has_splatting = false
    for arg in args
       if isa(arg, Expr) && arg.head == :...
          has_splatting = true
          break
       end
    end

    #
    #  Immediately check for corresponding arguments if there is no splatting
    #
    if !has_splatting && length(sym_args) != length(args)
       throw(ArgumentError("$macroname: wrong number of arguments ($(length(args))) should be ($(length(sym_args)))"))
    end

    for i = length(sym_args):-1:1
        var = sym_args[i].args[1]
        if has_splatting
           unshift!(blk.args, :($var = G[$i]))
        else
           unshift!(blk.args, :($var = $(esc(args[i]))))
        end
    end

    #
    #  Delay generation of argument list and check until evaluation time instead of macro
    #  expansion time if there is splatting.
    #
    if has_splatting
       x = Expr(:call,:tuple,args...)
       unshift!(blk.args,
          quote
             G = $(esc(x))
             if length(G) != $(length(sym_args))
                throw(ArgumentError($macroname,": wrong number of arguments (",length(G),") should be (",$(length(sym_args)),")"))
             end
          end
       )
    end

    unshift!(blk.args, :(out = $io))
    Expr(:let, Expr(:block), blk)
end

"""
    @printf([io::IOStream], "%Fmt", args...)

Print `args` using C `printf` style format specification string, with some caveats:
`Inf` and `NaN` are printed consistently as `Inf` and `NaN` for flags `%a`, `%A`,
`%e`, `%E`, `%f`, `%F`, `%g`, and `%G`. Furthermore, if a floating point number is
equally close to the numeric values of two possible output strings, the output
string further away from zero is chosen.

Optionally, an `IOStream`
may be passed as the first argument to redirect output.

# Examples
```jldoctest
julia> @printf("%f %F %f %F\\n", Inf, Inf, NaN, NaN)
Inf Inf NaN NaN\n

julia> @printf "%.0f %.1f %f\\n" 0.5 0.025 -0.0078125
1 0.0 -0.007813
```
"""
macro printf(args...)
    isempty(args) && throw(ArgumentError("@printf: called with no arguments"))
    if isa(args[1], AbstractString) || is_str_expr(args[1])
        _printf("@printf", :STDOUT, args[1], args[2:end])
    else
        (length(args) >= 2 && (isa(args[2], AbstractString) || is_str_expr(args[2]))) ||
            throw(ArgumentError("@printf: first or second argument must be a format string"))
        _printf("@printf", esc(args[1]), args[2], args[3:end])
    end
end

"""
    @sprintf("%Fmt", args...)

Return `@printf` formatted output as string.

# Examples
```jldoctest
julia> s = @sprintf "this is a %s %15.1f" "test" 34.567;

julia> println(s)
this is a test            34.6
```
"""
macro sprintf(args...)
    isempty(args) && throw(ArgumentError("@sprintf: called with zero arguments"))
    isa(args[1], AbstractString) || is_str_expr(args[1]) ||
        throw(ArgumentError("@sprintf: first argument must be a format string"))
    letexpr = _printf("@sprintf", :(IOBuffer()), args[1], args[2:end])
    push!(letexpr.args[2].args, :(String(take!(out))))
    letexpr
end

end # module

max_stars_repo_name: visr/julia

max_issues_repo_name: visr/julia

max_forks_repo_name: visr/julia


text:
#=
  - Julia Version: 1.4.0
  - Author       : Alexander Bills and Abhinav Misalkar
  - Date         : 12th May 2020
  - Title        : Taxonomy for Cost Inputs
=#

mutable struct struct_electrode_costs
    AM::Array{Float64,1}
    cond::Array{Float64,1}
    bind::Array{Float64,1}
    solvent::Array{Float64,1}
    baseline_AM_required::Array{Float64,1}
    baseline_solvent_evap::Array{Float64,1}
    bind_mult::Array{Float64,1}
    CC::Array{Float64,1}
    terminal_assembly::Array{Float64,1}
    material_yield::Array{Float64,1}
    cc_yield::Array{Float64,1}
end

mutable struct struct_cell_costs
    container::Array{Float64,1}
    seperator_cost
    seperator_yield
    electrolyte_cost
    electrolyte_yield
    cell_yield::Float64
    cathode::struct_electrode_costs
    anode::struct_electrode_costs

end

mutable struct struct_manufacturing
    receiving::Array{Float64,2}
    pos_mat_prep::Array{Float64, 2}
    neg_mat_prep::Array{Float64, 2}
    pos_electrode_coat::Array{Float64, 2}
    neg_electrode_coat::Array{Float64, 2}
    bind_slvnt_recov::Array{Float64, 2}
    pos_calendr::Array{Float64, 2}
    neg_calendr::Array{Float64, 2}
    intr_process_mat_handling::Array{Float64, 2}
    electrode_slitting::Array{Float64, 2}
    vacuum_drying::Array{Float64, 2}
    control_lab::Array{Float64, 2}
    cell_stacking::Array{Float64, 2}
    CC_welding::Array{Float64, 2}
    cell_insertion_container::Array{Float64, 2}
    electrolyte_filling_and_cell_sealing::Array{Float64, 2}
    dry_room_control::Array{Float64, 2}
    formation_cycling::Array{Float64, 2}
    cell_sealing::Array{Float64, 2}
    charge_retention_testing::Array{Float64, 2}
    rejected_cell_scrap_recycle_assembly::Array{Float64, 2}
    shipping::Array{Float64, 2}
end

mutable struct struct_general_costs
        pack_int::Float64                                # Pack Integration Cost                                  # USD/Pack
        building::Float64                                # Building cost                                          # USD/m2
        dir_labor::Float64                               # Labor cost                                             # $/hr
        no_units_mfg::Float64
        pos_bindr_slvnt_recovery::Float64                # Positive Binder Solvent Recovery percentage
        manufacturing::struct_manufacturing
        formation_cycling_mult::Float64
        cell_stacking_mult::Float64
        pos_electrode_coating_solvent_SR::Float64
        neg_electrode_coating_solvent_SR::Float64
end

mutable struct struct_baseline
    energy::Float64                          # Baseline annual energy                                   # kWh/year
    cells_adjusted_yield::Float64            # Baseline annual cells adjusted for yield                 # cells/year
    electrode_area::Float64                  # Baseline annual electrode area                           # m2/year
    dry_room_area::Float64                   # Baseline annual dry room area required                   # m2
end

mutable struct struct_pack_level
    no_cells_per_pack::Float64

end

mutable struct struct_ovhd_rate
    launch_dir_annual_materials::Float64             # Launch cost percent direct annual materials
    launch_dir_labor_var_ovhd::Float64               # Launch cost percent direct labor plus variable overhead
    depr_capital_equip::Float64                      # Depreciation rate capital equipment
    depr_building_inv::Float64                       # Depreciation rate building investment
    GSA_dir_labor_var_ovhd::Float64                  # General, Sales and Administration rate direct labor plus variable overhead
    GSA_depr::Float64                                # General, Sales and Administration rate depreciation
    RnD_depr::Float64                                # R & D rate depreciation
    var_ovhd_dir_labor::Float64                      # Variable overhead rate direct labor
    var_ovhd_depr::Float64                           # Variable overhead rate depreciation
    profit_after_taxes::Float64                      # Profit rate after taxes
    working_capital_rate::Float64                    # Working capital rate
    warranty::Float64                                # Warranty cost added to the price rate
end

mutable struct struct_cost
    cell_costs::struct_cell_costs
    manufacturing_costs::struct_manufacturing
    general_costs::struct_general_costs
    baseline::struct_baseline
    ovhd_rate::struct_ovhd_rate

end

max_issues_repo_name: battmodels/BatteryCost

max_forks_repo_name: battmodels/BatteryCost


text:
# Code that deals with objects and child objects.

PVectorCollection() = PVectorCollection(OrderedDict{ObjectId, ParamVector}());

function PVectorCollection(v :: Vector{ParamVector})
    pvecV = PVectorCollection();
    for pvec in v
        push!(pvecV, pvec);
    end
    return pvecV
end

Base.show(io :: IO, pvecV :: PVectorCollection) =
    print(io, "PVectorCollection of length $(length(pvecV))");

has_pvector(pvecV :: PVectorCollection, pv :: ParamVector) =
    has_pvector(pvecV, get_object_id(pv));

has_pvector(pvecV :: PVectorCollection, objId :: ObjectId) =
    haskey(pvecV.d, objId);

function Base.push!(pvecV :: PVectorCollection, pv :: ParamVector)
    @assert !has_pvector(pvecV, pv)  "$pv already exists.";
    pvecV.d[get_object_id(pv)] = pv;
end

Lazy.@forward PVectorCollection.d (
    Base.iterate, Base.isempty, Base.length
    );

function validate_pvec(pvecV :: PVectorCollection; silent = true)
    isValid = true;
    for (_, pvec) in pvecV
        isValid = isValid  &&  validate_pvec(pvec; silent);
    end
    return isValid
end

"""
	$(SIGNATURES)

Collect all `ParamVector`s in an object and its child objects.
Returns an empty `ParamVector` if nothing found.
"""
function collect_pvectors(o :: ModelObject)
    pvecV = PVectorCollection();
    objV = collect_model_objects(o);
    if !isempty(objV)
        for i1 = 1 : length(objV)
            if has_pvector(objV[i1])
                push!(pvecV, get_pvector(objV[i1]));
            end
        end
    end
    return pvecV :: PVectorCollection
end

"""
    $(SIGNATURES)

Find a `ParamVector` for a given ObjectId.
Returns the `ParamVector`. Or `nothing` if not found.
"""
find_pvector(o :: ModelObject, objId) = find_pvector(collect_pvectors(o), objId);

function find_pvector(pvv :: PVectorCollection, objId :: ObjectId)
    if has_pvector(pvv, objId)
        return pvv.d[objId];
    else
        return nothing;
    end
end

"""
	$(SIGNATURES)

Find all occurrences of a given parameter name.
Returns a `Dict` that maps `ObjectId => Param`. Or is empty.

# Example
```julia
d = find_param(o, :alpha);
for (objId, p) in d
    # Since a reference to the `Param` is returned, this changes the calibration status in `o`.
    calibrate!(p);
    pvec = find_pvector(o, objId);
end
```
"""
function find_param(o :: ModelObject, pName)
    pvv = collect_pvectors(o);
    return find_param(pvv, pName);
end

function find_param(pvv :: PVectorCollection, pName)
    d = Dict{ObjectId, Param}();
    for (objId, pvec) in pvv
        if param_exists(pvec, pName)
            p = retrieve(pvec, pName);
            d[objId] = p;
        end
    end
    return d
end

## --------------  Vectors and Dicts

"""
	$(SIGNATURES)

Make a `PVectorCollection` into a `Dict{String, Dict{Symbol, Any}}`.
Each entry is one `ParamVector`.
They key is the `ObjectId` of each `ParamVector` made into a `String`. Such as "parent > child > grandchild[2, 1]".
The value is the `ParamVector` made into a Dict.
Objects without calibrated params have no entries.

This is a format that can be saved without using user defined types. There is hope this can be serialized.

All values are stored for parameters where some values are calibrated and others are fixed.
"""
function make_dict(pvv :: PVectorCollection; isCalibrated :: Bool = true)
    n = length(pvv);
    if n < 1
        return nothing
    end
    isCalibrated  ?  (vType = :value)  :  (vType = :defaultValue);

    d = nothing;
    for (objId, pv) in pvv.d
        # pv = pvv[j];
        key = make_string(objId);
        pd = make_dict(pv; isCalibrated, valueType = vType);
        if isnothing(d)
            d = Dict([key => pd]);
        else
            push!(d, key => pd);
        end
    end

    return d
end

## -------------  Setting values

"""
	$(SIGNATURES)

Set model values from `Dict{String, Dict}` generated by `make_dict`.
Includes child objects.
"""
function set_values_from_dicts!(
    x :: ModelObject,
    pvDict :: D1;
    isCalibrated :: Bool = true,
    skipInvalidSize :: Bool = false
    ) where D1 <: AbstractDict

    # Collect all model objects
    mObjV = collect_model_objects(x);

    # Loop over `ParamVector`s, represented as Dicts
    # `nameStr` is the `ObjectId` converted into a `String`
    # `pd` is the `ParamVector` converted into a `Dict`
    for (nameStr,  pd) in pvDict
        # Make string into ObjectId. Allow for `nameStr` to be something that can be made into `String`
        oId = make_object_id(string(nameStr));
        # Find the matching model object
        obj = find_object(x, oId);
        if !isnothing(obj)  &&  has_pvector(obj)
            set_values_from_dict!(get_pvector(obj), pd; skipInvalidSize);
            set_own_values_from_pvec!(obj, isCalibrated);
        end
    end
    return nothing
end

"""
	$(SIGNATURES)

Copy all values from a vector of `ParamVector` into an object, including child objects.
Only changes values that are `isCalibrated` in object and `v`.

# Arguments
- skipInvalidSize: skip parameters that have the wrong size. Otherwise throw error.
"""
function set_values_from_pvectors!(
    x :: ModelObject,
    v :: PVectorCollection,
    isCalibrated :: Bool;
    skipInvalidSize :: Bool = false
    )
    # Collect all model objects
    mObjV = collect_model_objects(x);

    # Loop over `ParamVector`s
    for (objId, pvec) in v.d
        # Find the matching model object
        obj = find_object(x, objId);
        if !isnothing(obj)  &&  has_pvector(obj)
            set_own_values_from_pvec!(get_pvector(obj), pvec, isCalibrated;
                skipInvalidSize);
            set_own_values_from_pvec!(obj, isCalibrated);
        end
    end
    return nothing
end

## --------  Testing

function make_test_pvector_collection(n)
    pvecV = PVectorCollection();
    for j = 1 : n
        oId = Symbol("id$j");
        push!(pvecV, make_test_pvector(n; objId = oId));
    end
    return pvecV
end

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