text:
function mapreduce_index(f::Vector{<:Function}, op, itr, init)
    y = iterate(itr)
    y === nothing && return init
    v = op(init, (y[1]), f[y[2][2]])
    while true
        y = iterate(itr, y[2])
        y === nothing && break
        v = op(v, (y[1]), f[y[2][2]])
    end
    return v
end

# these functions are experimental and later must be moved to row_functions or hp_row_functions
function row_any_multi(ds::AbstractDataset, f::Vector{<:Function}, cols = :)
    colsidx = index(ds)[cols]
    @assert length(f) == length(colsidx) "number of provided functions must match the number of selected columns"
    _op_bool_add_multi(x::Bool,y::Bool) = x | y ? true : false
    op_for_any_multi!(x, y, f) = x .= _op_bool_add_multi.(x, f.(y))
    # mapreduce(identity, op_for_anymissing!, eachcol(ds)[colsidx[sel_colsidx]], init = zeros(Bool, size(ds,1)))
    mapreduce_index(f, op_for_any_multi!, view(_columns(ds),colsidx), zeros(Bool, size(ds,1)))
end

function row_all_multi(ds::AbstractDataset, f::Vector{<:Function}, cols = :)
    colsidx = index(ds)[cols]
    @assert length(f) == length(colsidx) "number of provided functions must match the number of selected columns"
    _op_bool_mult_multi(x::Bool,y::Bool) = x & y ? true : false
    op_for_all_multi!(x, y, f) = x .= _op_bool_mult_multi.(x, f.(y))
    # mapreduce(identity, op_for_anymissing!, eachcol(ds)[colsidx[sel_colsidx]], init = zeros(Bool, size(ds,1)))
    mapreduce_index(f, op_for_all_multi!, view(_columns(ds),colsidx), ones(Bool, size(ds,1)))
end

_op_bool_add_multi(x::Bool,y::Bool) = x || y ? true : false

function hp_bool_add_multi!(x, y, f)
    Threads.@threads for i in 1:length(x)
        @inbounds x[i] = _op_bool_add_multi(x[i], f(y[i]))
    end
    x
end

function hp_row_any_multi(ds::AbstractDataset, f::Vector{<:Function}, cols = :)
    colsidx = index(ds)[cols]
    @assert length(f) == length(colsidx) "number of provided functions must match the number of selected columns"
    _hp_op_for_any_multi!(x, y, f) = x .= hp_bool_add_multi!(x, y, f)
    # mapreduce(identity, op_for_anymissing!, eachcol(ds)[colsidx[sel_colsidx]], init = zeros(Bool, size(ds,1)))
    mapreduce_index(f, _hp_op_for_any_multi!, view(_columns(ds),colsidx), zeros(Bool, size(ds,1)))
end

_op_bool_mult_multi(x::Bool,y::Bool) = x && y ? true : false

function hp_bool_mult_multi!(x, y, f)
    Threads.@threads for i in 1:length(x)
        @inbounds x[i] = _op_bool_mult_multi(x[i], f(y[i]))
    end
    x
end

function hp_row_all_multi(ds::AbstractDataset, f::Vector{<:Function}, cols = :)
    colsidx = index(ds)[cols]
    @assert length(f) == length(colsidx) "number of provided functions must match the number of selected columns"
    _hp_op_for_all_multi!(x, y, f) = x .= hp_bool_mult_multi!(x, y, f)
    # mapreduce(identity, op_for_anymissing!, eachcol(ds)[colsidx[sel_colsidx]], init = zeros(Bool, size(ds,1)))
    mapreduce_index(f, _hp_op_for_all_multi!, view(_columns(ds),colsidx), ones(Bool, size(ds,1)))
end

__STRING(x) = string(x)
__STRING(x::AbstractString) = x
__STRING(::Missing) = ""
__STRING(x::Bool) = x ? "1" : "0"
__STRING(::Nothing) = ""

__CODEUNIT(x) = Base.CodeUnits(__STRING(x))

#
# function _op_for_row_join(x, y, f, delim, quotechar, idx, p)
#     idx[] += 1
#     if quotechar === nothing
#         if idx[] < p
#             x .*= STRING.(f[idx[]].(y))
#             x .*= delim
#         else
#             x .*= STRING.(f[idx[]].(y))
#             x .*= '\n'
#         end
#     else
#         if idx[] < p
#             x .*= quotechar
#             x .*= STRING.(f[idx[]].(y))
#             x .*= quotechar
#             x .*= delim
#         else
#             x .*= quotechar
#             x .*= STRING.(f[idx[]].(y))
#             x .*= quotechar
#             x .*= '\n'
#         end
#     end
#     x
# end
#
# function row_join(ds::AbstractDataset, f::Vector{<:Function}, cols  = :; delim = ',', quotechar = nothing)
#     colsidx = index(ds)[cols]
#     idx = Ref{Int}(0)
#     p = length(colsidx)
#     init0 = fill("", nrow(ds))
#     mapreduce(identity, (x,y)->_op_for_row_join(x,y,f, delim, quotechar, idx, p), view(_columns(ds), colsidx), init = init0)
# end

# some experimental ideas

function write_vals!(a, pos, x::Integer)
    n_positive, neg = Base.split_sign(x)
    #always we assume base = 10
    needed_space = neg + ndigits(x)
    available_space = length(a)-pos+1
    needed_space > available_space && throw(ArgumentError("not enough space in buffer to write value into it"))
    _base!(a, pos, 10, n_positive, 0, neg)
    pos+needed_space
end

#TODO this doesn't work - since in writeshortest() a must be a vector{UInt8} but we pass a view of a matrix in DLMReader
function write_vals!(a, pos, x::Union{Float16, Float32, Float64})
    needed_space = Base.Ryu.neededdigits(typeof(x))
    available_space = length(a)-pos+1
    # needed_space > available_space && throw(ArgumentError("not enough space in buffer to write value into it"))
    _writeshortest(a, pos, x)
end

function write_vals!(a, pos, x::Bool)
    needed_space = 1
    available_space = length(a)-pos+1
    needed_space > available_space && throw(ArgumentError("not enough space in buffer to write value into it"))
    x ? a[pos] = UInt8('1') : a[pos] = UInt8('0')
    pos+1
end

function write_vals!(a, pos, ::Missing)
    pos
end
function write_vals!(a, pos, ::Nothing)
    pos
end

function write_vals!(a, pos, x::AbstractString)
    needed_space = length(x)
    available_space = length(a)-pos+1
    needed_space > available_space && throw(ArgumentError("not enough space in buffer to write value into it"))
    z = Base.CodeUnits(x)
    for i in 1:length(x)
        a[i+pos-1] =z[i]
    end
    pos+length(x)
end

function write_vals!(a, pos, x)
    y = string(x)
    needed_space = length(y)
    available_space = length(a)-pos+1
    needed_space > available_space && throw(ArgumentError("not enough space in buffer to write value into it"))
    for i in 1:length(y)
        a[i+pos-1] = Base.CodeUnits(y)[i]
    end
    pos+length(y)
end

function write_delim!(a, pos, delim)
    needed_space = length(delim)
    available_space = length(a)-pos+1
    needed_space > available_space && throw(ArgumentError("not enough space in buffer to write value into it"))
    for i in 1:length(delim)
        a[i+pos-1] = delim[i]
    end
    pos+length(delim)
end

function write_eol!(a, pos)
    needed_space = 1
    available_space = length(a)-pos+1
    needed_space > available_space && throw(ArgumentError("not enough space in buffer to write value into it"))
    a[pos] = UInt8('\n')
    pos+1
end

function write_quotechar!(a, pos, quotechar)
    needed_space = 1
    available_space = length(a)-pos+1
    needed_space > available_space && throw(ArgumentError("not enough space in buffer to write value into it"))
    a[pos] = quotechar
    pos+1
end
function write_quotechar!(a, pos, ::Nothing)
    pos
end

function _base!(a, pos, base::Integer, x::Integer, pad::Int, neg::Bool)
    (x >= 0) | (base < 0) || throw(DomainError(x, "For negative `x`, `base` must be negative."))
    2 <= abs(base) <= 62 || throw(DomainError(base, "base must satisfy 2 ≤ abs(base) ≤ 62"))
    b = (base % Int)::Int
    digits = abs(b) <= 36 ? Base.base36digits : Base.base62digits
    # pad = 0 makes issue when x == 0 (n will be 0)
    n = neg + ndigits(x, base=b)
    i = n
    @inbounds while i > neg
        if b > 0
            a[i+pos-1] = digits[1 + (rem(x, b) % Int)::Int]
            x = div(x,b)
        else
            a[i+pos-1] = digits[1 + (mod(x, -b) % Int)::Int]
            x = cld(x,b)
        end
        i -= 1
    end
    if neg; @inbounds a[pos]=0x2d; end
    a
end

function _op_for_row_join!(buffer, currentpos, y, f, delim, quotechar, idx, p, lo, hi)
    idx[] += 1
    if quotechar === nothing
        if idx[] < p
            for i in lo:hi
                currentpos[i] = write_vals!(view(buffer, :, i), currentpos[i], f[idx[]](y[i]))
                currentpos[i] = write_delim!(view(buffer, :, i), currentpos[i], delim)
            end
        else
            for i in lo:hi
                currentpos[i] = write_vals!(view(buffer, :, i), currentpos[i], f[idx[]](y[i]))
                currentpos[i] = write_eol!(view(buffer, :, i), currentpos[i])
            end
        end
    else
        if nonmissingtype(eltype(y)) <: AbstractString
            quotecharval = UInt8(quotechar)
        else
            quotecharval = nothing
        end
        if idx[]<p
            for i in lo:hi
                currentpos[i] = write_quotechar!(view(buffer, :, i), currentpos[i], quotecharval)
                currentpos[i] = write_vals!(view(buffer, :, i), currentpos[i], f[idx[]](y[i]))
                currentpos[i] = write_quotechar!(view(buffer, :, i), currentpos[i], quotecharval)
                currentpos[i] = write_delim!(view(buffer, :, i), currentpos[i], delim)
            end
        else
            for i in lo:hi
                currentpos[i] = write_quotechar!(view(buffer, :, i), currentpos[i], quotecharval)
                currentpos[i] = write_vals!(view(buffer, :, i), currentpos[i], f[idx[]](y[i]))
                currentpos[i] = write_quotechar!(view(buffer, :, i), currentpos[i], quotecharval)
                currentpos[i] = write_eol!(view(buffer, :, i), currentpos[i])
            end
        end
    end
    currentpos
end

# buffer = Matrix{UInt8}(undef, lsize , nrow(ds))
# currentpos = ones(Int, nrow(ds))
function row_join!(buffer, currentpos, ds::AbstractDataset, f::Vector{<:Function}, cols = :; delim = ',', quotechar = nothing, threads = true)
    colsidx = multiple_getindex(index(ds), cols)
    p = length(colsidx)
    dlm = UInt8.(delim)
    if threads
        cz = div(nrow(ds), __NCORES)
        idx = [Ref{Int}(0) for _ in 1:__NCORES]
        Threads.@threads for i in 1:__NCORES
            lo = (i-1)*cz+1
            i == __NCORES ? hi = nrow(ds) : hi = i*cz
            mapreduce(identity, (x,y) -> _op_for_row_join!(buffer, x, y, f, dlm, quotechar, idx[i], p, lo, hi), view(_columns(ds),colsidx), init = currentpos)
        end
    else
        idx = Ref{Int64}(0)
        mapreduce(identity, (x,y) -> _op_for_row_join!(buffer, x, y, f, dlm, quotechar, idx, p, 1, nrow(ds)), view(_columns(ds),colsidx), init = currentpos)
    end
    currentpos, buffer
end

@inline function append_sign(x, plus, space, buf, pos)
    if signbit(x) && !isnan(x)  # suppress minus sign for signaling NaNs
        buf[pos] = UInt8('-')
        pos += 1
    elseif plus
        buf[pos] = UInt8('+')
        pos += 1
    elseif space
        buf[pos] = UInt8(' ')
        pos += 1
    end
    return pos
end

### From Base.Ryu, because we need buf to be View of an array not vector (maybe we should change it in Ryu?)
function _writeshortest(buf, pos, x::T,
                       plus=false, space=false, hash=true,
                       precision=-1, expchar=UInt8('e'), padexp=false, decchar=UInt8('.'),
                       typed=false, compact=false) where {T}
    @assert 0 < pos <= length(buf)
    # special cases
    if x == 0
        if typed && x isa Float16
            buf[pos] = UInt8('F')
            buf[pos + 1] = UInt8('l')
            buf[pos + 2] = UInt8('o')
            buf[pos + 3] = UInt8('a')
            buf[pos + 4] = UInt8('t')
            buf[pos + 5] = UInt8('1')
            buf[pos + 6] = UInt8('6')
            buf[pos + 7] = UInt8('(')
            pos += 8
        end
        pos = append_sign(x, plus, space, buf, pos)
        buf[pos] = UInt8('0')
        pos += 1
        if hash
            buf[pos] = decchar
            pos += 1
        end
        if precision == -1
            buf[pos] = UInt8('0')
            pos += 1
            if typed && x isa Float32
                buf[pos] = UInt8('f')
                buf[pos + 1] = UInt8('0')
                pos += 2
            end
            if typed && x isa Float16
                buf[pos] = UInt8(')')
                pos += 1
            end
            return pos
        end
        while hash && precision > 1
            buf[pos] = UInt8('0')
            pos += 1
            precision -= 1
        end
        if typed && x isa Float32
            buf[pos] = UInt8('f')
            buf[pos + 1] = UInt8('0')
            pos += 2
        end
        if typed && x isa Float16
            buf[pos] = UInt8(')')
            pos += 1
        end
        return pos
    elseif isnan(x)
        pos = append_sign(x, plus, space, buf, pos)
        buf[pos] = UInt8('N')
        buf[pos + 1] = UInt8('a')
        buf[pos + 2] = UInt8('N')
        if typed
            if x isa Float32
                buf[pos + 3] = UInt8('3')
                buf[pos + 4] = UInt8('2')
            elseif x isa Float16
                buf[pos + 3] = UInt8('1')
                buf[pos + 4] = UInt8('6')
            end
        end
        return pos + 3 + (typed && x isa Union{Float32, Float16} ? 2 : 0)
    elseif !isfinite(x)
        pos = append_sign(x, plus, space, buf, pos)
        buf[pos] = UInt8('I')
        buf[pos + 1] = UInt8('n')
        buf[pos + 2] = UInt8('f')
        if typed
            if x isa Float32
                buf[pos + 3] = UInt8('3')
                buf[pos + 4] = UInt8('2')
            elseif x isa Float16
                buf[pos + 3] = UInt8('1')
                buf[pos + 4] = UInt8('6')
            end
        end
        return pos + 3 + (typed && x isa Union{Float32, Float16} ? 2 : 0)
    end

    output, nexp = Base.Ryu.reduce_shortest(x, compact ? 999_999 : nothing)

    if typed && x isa Float16
        buf[pos] = UInt8('F')
        buf[pos + 1] = UInt8('l')
        buf[pos + 2] = UInt8('o')
        buf[pos + 3] = UInt8('a')
        buf[pos + 4] = UInt8('t')
        buf[pos + 5] = UInt8('1')
        buf[pos + 6] = UInt8('6')
        buf[pos + 7] = UInt8('(')
        pos += 8
    end
    pos = append_sign(x, plus, space, buf, pos)

    olength = Base.Ryu.decimallength(output)
    exp_form = true
    pt = nexp + olength
    if -4 < pt <= (precision == -1 ? (T == Float16 ? 3 : 6) : precision) &&
        !(pt >= olength && abs(mod(x + 0.05, 10^(pt - olength)) - 0.05) > 0.05)
        exp_form = false
        if pt <= 0
            buf[pos] = UInt8('0')
            pos += 1
            buf[pos] = decchar
            pos += 1
            for _ = 1:abs(pt)
                buf[pos] = UInt8('0')
                pos += 1
            end
            # elseif pt >= olength
            # nothing to do at this point
            # else
            # nothing to do at this point
        end
    else
        pos += 1
    end
    i = 0
    ptr = pointer(buf)
    ptr2 = pointer(Base.Ryu.DIGIT_TABLE)
    if (output >> 32) != 0
        q = output ÷ 100000000
        output2 = (output % UInt32) - UInt32(100000000) * (q % UInt32)
        output = q

        c = output2 % UInt32(10000)
        output2 = div(output2, UInt32(10000))
        d = output2 % UInt32(10000)
        c0 = (c % 100) << 1
        c1 = (c ÷ 100) << 1
        d0 = (d % 100) << 1
        d1 = (d ÷ 100) << 1
        Base.Ryu.memcpy(ptr, pos + olength - 2, ptr2, c0 + 1, 2)
        Base.Ryu.memcpy(ptr, pos + olength - 4, ptr2, c1 + 1, 2)
        Base.Ryu.memcpy(ptr, pos + olength - 6, ptr2, d0 + 1, 2)
        Base.Ryu.memcpy(ptr, pos + olength - 8, ptr2, d1 + 1, 2)
        i += 8
    end
    output2 = output % UInt32
    while output2 >= 10000
        c = output2 % UInt32(10000)
        output2 = div(output2, UInt32(10000))
        c0 = (c % 100) << 1
        c1 = (c ÷ 100) << 1
        Base.Ryu.memcpy(ptr, pos + olength - i - 2, ptr2, c0 + 1, 2)
        Base.Ryu.memcpy(ptr, pos + olength - i - 4, ptr2, c1 + 1, 2)
        i += 4
    end
    if output2 >= 100
        c = (output2 % UInt32(100)) << 1
        output2 = div(output2, UInt32(100))
        Base.Ryu.memcpy(ptr, pos + olength - i - 2, ptr2, c + 1, 2)
        i += 2
    end
    if output2 >= 10
        c = output2 << 1
        buf[pos + 1] = Base.Ryu.DIGIT_TABLE[c + 2]
        buf[pos - exp_form] = Base.Ryu.DIGIT_TABLE[c + 1]
    else
        buf[pos - exp_form] = UInt8('0') + (output2 % UInt8)
    end

    if !exp_form
        if pt <= 0
            pos += olength
            precision -= olength
            while hash && precision > 0
                buf[pos] = UInt8('0')
                pos += 1
                precision -= 1
            end
        elseif pt >= olength
            pos += olength
            precision -= olength
            for _ = 1:nexp
                buf[pos] = UInt8('0')
                pos += 1
                precision -= 1
            end
            if hash
                buf[pos] = decchar
                pos += 1
                if precision < 0
                    buf[pos] = UInt8('0')
                    pos += 1
                end
                while precision > 0
                    buf[pos] = UInt8('0')
                    pos += 1
                    precision -= 1
                end
            end
        else
            pointoff = olength - abs(nexp)
            Base.Ryu.memmove(ptr, pos + pointoff + 1, ptr, pos + pointoff, olength - pointoff + 1)
            buf[pos + pointoff] = decchar
            pos += olength + 1
            precision -= olength
            while hash && precision > 0
                buf[pos] = UInt8('0')
                pos += 1
                precision -= 1
            end
        end
        if typed && x isa Float32
            buf[pos] = UInt8('f')
            buf[pos + 1] = UInt8('0')
            pos += 2
        end
    else
        if olength > 1 || hash
            buf[pos] = decchar
            pos += olength
            precision -= olength
        end
        if hash && olength == 1
            buf[pos] = UInt8('0')
            pos += 1
        end
        while hash && precision > 0
            buf[pos] = UInt8('0')
            pos += 1
            precision -= 1
        end

        buf[pos] = expchar
        pos += 1
        exp2 = nexp + olength - 1
        if exp2 < 0
            buf[pos] = UInt8('-')
            pos += 1
            exp2 = -exp2
        elseif padexp
            buf[pos] = UInt8('+')
            pos += 1
        end

        if exp2 >= 100
            c = exp2 % 10
            Base.Ryu.memcpy(ptr, pos, ptr2, 2 * div(exp2, 10) + 1, 2)
            buf[pos + 2] = UInt8('0') + (c % UInt8)
            pos += 3
        elseif exp2 >= 10
            Base.Ryu.memcpy(ptr, pos, ptr2, 2 * exp2 + 1, 2)
            pos += 2
        else
            if padexp
                buf[pos] = UInt8('0')
                pos += 1
            end
            buf[pos] = UInt8('0') + (exp2 % UInt8)
            pos += 1
        end
    end
    if typed && x isa Float16
        buf[pos] = UInt8(')')
        pos += 1
    end

    return pos
end


text:
using PyPlot
using PyCall
mpl = pyimport("tikzplotlib")
using JLD2

rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
font0 = Dict(
        "font.size" => 16,
        "axes.labelsize" => 16,
        "xtick.labelsize" => 16,
        "ytick.labelsize" => 16,
        "legend.fontsize" => 16,
)
merge!(rcParams, font0)

close("all")
for tid = 1:5
    @show tid
@load "Data/domain$(tid).jld2" domain

strain = hcat(domain.history["strain"]...)
stress = hcat(domain.history["stress"]...)
i=8
markevery=5
plot(strain[i,1:markevery:end], stress[i,1:markevery:end], "--o", fillstyle="none", label="Reference-tid$(tid)")

end

xlabel("Strain")
ylabel("Stress (MPa)")
legend()
PyPlot.tight_layout()
savefig("all_truss1d_strain_stress.pdf")
mpl.save("all_truss1d_strain_stress.tex")

ext: jl


text:
# This computes the scaling of the TTVs with the mass-ratios
# of the planets to the star, "epsilon", as well as the stepsize, h.
#

include("../../src/ttv.jl")
include("/Users/ericagol/Computer/Julia/regress.jl")

using PyPlot

# This routine takes derivative of transit times with respect
# to the initial orbital elements.
n = 3
n_body = n
t0 = 7257.93115525
h  = 0.07
tmax = 100.0

# Read in initial conditions:
elements = readdlm("elements.txt",',')

ntt = zeros(Int64,n)

# Make an array, tt,  to hold transit times:
# First, though, make sure it is large enough:
for i=2:n
  ntt[i] = ceil(Int64,tmax/elements[i,2])+3
end
dtdq0 = zeros(n,maximum(ntt),7,n)
tt  = zeros(n,maximum(ntt))
tt1 = zeros(n,maximum(ntt))

# Make an array to save the variables as a function of step
# size and mass:
nh = 5
nmass = 5
tt_save = zeros(BigFloat,nh,nmass,n,maximum(ntt))
# Save a counter for the actual number of transit times of each planet:
count = zeros(Int64,n)
count1 = zeros(Int64,n)
# Call the ttv function:
rstar = 1e12

mass_fac = big.([1.0,0.5,0.25,0.125,0.0625])
h_list = [h,h/2.,h/4.,h/8.,h/32.]
t0big = big(t0)
tmaxbig = big(tmax)
tt1big = big.(tt1)
rstarbig = big(rstar)
for i=1:nh, j=1:nmass
  tic()
# Create BigFloat versions of the input variables:
  elements_big = convert(Array{BigFloat,2},elements)
  elements_big[2:n,1] *= mass_fac[j]
  hbig = big(h_list[i])
  # Call the routine
  dqbig = ttv_elements!(n,t0big,hbig,tmaxbig,elements_big,tt1big,count,big(0.0),0,0,rstarbig)
  # Save this:
#  tt_save[i,j,:,:]=convert(Array{Float64,2},tt1big)
  tt_save[i,j,:,:]=tt1big
  println("finished h: ",convert(Float64,h_list[i])," mass_fac: ",convert(Float64,mass_fac[j])," elapsed: ",toq())
end

# Make a plot of transit time errors versus stepsize:
ntrans = sum(count)
clf()
tab = 0
hlabel = ["h-h/32","h/2-h/32","h/4-h/32","h/8-h/32"]
ch = ["black","red","green","blue","orange"]

# Make a plot of timing errors versus stepsize and epsilon:
ntrans = sum(count)
clf()
sigt = zeros(n-1,nh-1,nmass)
tab = 0
hlabel = ["h-h/32","h/2-h/32","h/4-h/32","h/8-h/32","h/16-big(h/32)"]
ch = ["black","red","green","blue","orange"]
for i=2:n
  for k=1:nmass
#    fn = zeros(Float64,2,count[i])
    fn = zeros(BigFloat,2,count[i])
#    sig = ones(count[i])
    sig = ones(BigFloat,count[i])
    tti16 = tt_save[nh,k,i,1:count[i]]
    for j=1:count[i]
#      fn[1,j] = 1.0
      fn[1,j] = big(1.0)
      fn[2,j] = round(Int64,(tti16[j]-elements[i,3])/elements[i,2])
    end
    for j=1:nh-1
      tti1 = tt_save[j,k,i,1:count[i]]
      coeff,cov = regress(fn,tti1-tti16,sig)
      diff = tti1-tti16-coeff[1]-coeff[2]*fn[2,:]
      sigt[i-1,j,k]=std(diff)
      if i == n
        plot(tti1,diff/h_list[j]^4/mass_fac[k],linestyle="dashed",c=ch[j])
      else
        plot(tti1,diff/h_list[j]^4/mass_fac[k],label=hlabel[j],c=ch[j])
      end
    end
  end
end
legend(loc="lower left")

read(STDIN,Char)
clf()

# Make a plot of some TTVs vs mass scaling:
# Try p=1:
p = 1
for j=1:nh-1
  loglog(mass_fac[1:nmass],sigt[1,j,:]*24.*3600./(h_list[j]/h[1]).^4,".",markersize=15,label="Inner planet")
  loglog(mass_fac[1:nmass],sigt[1,j,1]*24.*3600./(h_list[j]/h[1]).^4*mass_fac[1:nmass].^p,label=L"$\propto m^p$")
  loglog(mass_fac[1:nmass],sigt[2,j,:]*24.*3600./(h_list[j]/h[1]).^4,".",markersize=15,label="Outer planet")
  loglog(mass_fac[1:nmass],sigt[2,j,1]*24.*3600./(h_list[j]/h[1]).^4*mass_fac[1:nmass].^p,label=L"$\propto m^p$")
end
legend(loc = "upper left")
ylabel("RMS timing error [sec]")
xlabel("Mass ratio [day]")

PyPlot.savefig("timing_error_vs_mass_big.pdf",bbox_inches="tight")

max_forks_repo_name: langfzac/nbg-papers


text:
using KrylovKit
using LinearAlgebra
using DelimitedFiles
using Plots

# Load the file that contains the Polyharmonic interpolator function
include("polyHarmonicInterp.jl")
include("bisectionMethod.jl")

function calculateParents(offspring::Vector{T}, minChrom::Int,
	maxChrom::Int,stepChrom::Int) where T<:Real

	[ (x = copy(offspring); x[idx] = v; x)
		for idx in 1 : length(offspring)
			for v in max(offspring[idx]-1,minChrom):stepChrom:min(offspring[idx]+1,maxChrom)
				 if v != offspring[idx] ]

end

function q(parent::Vector{T},offspring::Vector{T},misRate::Float64,nChrom::Int) where T <: Real

	# parentCN = collect(linspaces[k][parent[k]] for k in 1:nChrom)
	# offspringCN = collect(linspaces[k][offspring[k]] for k in 1:nChrom)

	parentCN = parent
	offspringCN = offspring

	CNdist = norm(parentCN .- offspringCN,1)

	# No mutation occurred
	if (CNdist == 0.0)
		val = 1.0 - misRate
	else
		# FIXME: Should be a function dependent on parent ploidy
		val = misRate/nChrom
	end

	return val

end

"""

	jacobianMap()

	This function defines our Jacobian map which is used by Krylov to calculate
	the maximum eigenvalues of the Jacobian matrix. Of interest is the point at
	which the maximum eigenvalue is 0 (the critical curve) which separates
	extinction from growth.

	INPUTS
		x: 			A vector which is acted on the map of J (J*x -> y)
		mutRate:	The mutation rate of the ploidy types
		birthRate:	The birth rate of the ploidy types
		deathRate:	The death rate of the ploidy types
		matSize:	The number of ploidy types

	OUTPUTS
		y: 			The result of J*x -> y

"""
function jacobianMap(x::Vector{Float64}, misRate::Float64, pars)

	# Initialize the output vector
	y = zeros(length(x))

	(interp,nChrom,chromArray,deathRate,maxY) = pars

	# Linear indexing
	lidx = LinearIndices(zeros(5,5,5,5,5))

	# Iterate over each compartment
	for (i,focal) in enumerate(
		Iterators.product((1:length(j) for j in chromArray)...))

		focalCN = collect(chromArray[k][focal[k]] for k in 1:nChrom)

		parentCNList = calculateParents(focalCN, 1,5,1)

		inflow = 0.0

		# Inflow from other parents
		for parentCN in parentCNList

			# Get birth rate
			birthRate = min(max(PolyharmonicInterpolation.polyharmonicSpline(interp,parentCN')[1],0.0),maxY)

			# Get flow rate from parentCN -> focalCN
			flowRate = q(parentCN,focalCN,misRate,nChrom)

			# Will need to have this be the index in the future
			inflow += birthRate*flowRate*x[lidx[parentCN...]]

		end

		# Grab the focal cells birth rate
		birthRate = min(max(PolyharmonicInterpolation.polyharmonicSpline(interp,focalCN')[1],0.0),maxY)

		inflow += birthRate*(1.0 - 2*misRate)*x[i]

		# Death term
		outflow = deathRate*x[i]

		y[i] = inflow - outflow

	end

	return y

end

function getParameterInfo()

	# Read the birth rate file to be used for the interpolation
	X = readdlm("CNVs brain cancer CLs.txt", '\t')
	Y = readdlm("GrowthRate brain cancer CLs.txt", '\t')

	# Prune data to keep only the numeric part
	X = float.(X[2:end,2:end]) .+ 2.0
	Y = float.(Y[:,2])

	# Polyharmonic interpolator
	interp = PolyharmonicInterpolation.PolyharmonicInterpolator(X,Y)

	nChrom = interp.dim

	minX,maxX = minimum(X,dims=1),maximum(X,dims=1)

	stepsize = 1

	chromArray = [1:stepsize:5 for i in 1:nChrom]

	# Total number of compartments
	matSize = prod(length(arr) for arr in chromArray)

	return interp,nChrom,chromArray,matSize,maximum(Y)

end

function getEigenvalues(misRate::Float64,matSize::Int,constParams::Tuple)

	# Find maximum eigenvalues
	evals,evecs,info = (eigsolve(matSize, 1,:LR; verbosity=0,issymmetric=false,
	krylovdim=50,tol=1e-4)
	do u
		# This computes y = J*u (the action of J on the vector u)
		jacobianMap(u,misRate,constParams)
	end)

	# Check if eigenvalue converged
	try
		return real(evals[1])
	catch
		error("Eigenvalue did not converge")
	end

end

# Parameters
deathRateVec = (0.0:0.05:0.7)
interp,nChrom,chromArray,matSize,maxY = getParameterInfo()

# Initialize solution vector
beta = Vector{Float64}()

for deathRate in reverse(deathRateVec)

	# Params for the linear map
	constParams = (interp,nChrom,chromArray,deathRate,maxY)

	# Define the function
	myfun = misRate -> getEigenvalues(misRate,matSize,constParams)

	# Execute bisection search method
	results = Bisection.bisectionSearch(0.0,1.0,myfun)

	# sgn(f(a)) = sgn(f(b)) - bisection does not work. We check to see the sign
	# individually to determine whether any misRate is acceptable or none is.
	if typeof(results) <: Tuple
		a,fa,b,fb = results

		# If f(b) < 0 - Any missegregation rate is acceptable
		if fb < 0
			push!(beta,0.0)
		else
			push!(beta,NaN)
		end
	# The signs are unequal and the solver went to work...
	else
		push!(beta,results.best)
	end

	@show deathRate,beta[end]

end

# # Initialize solution vector
# beta = Vector{Float64}()

# for deathRate in reverse(deathRateVec)

# 	result = (optimize(mRate -> lossFunction(mRate,
# 	matSize,(interp,nChrom,chromArray,deathRate)),[0.05], LBFGS()))

# 	@show deathRate, result.minimizer[1]

# 	push!(beta,result.minimizer[1])

# 	initGuess = result.minimizer[1]

# end

# # Beta is a probability and so we remove unphysical solutions by the solver
# # (Alternatively, we can do box-constrained optimization if necessary)
# beta[beta .> 1.0] .= 1.0
# beta[beta .< 0.0] .= 0.0

max_stars_repo_name: MathOnco/ploidyEvolution

max_issues_repo_name: MathOnco/ploidyEvolution

max_forks_repo_name: MathOnco/ploidyEvolution


text:
root = joinpath(@__DIR__, "..")
using Pkg; Pkg.activate(root)
using Literate

src = joinpath(root, "notebook_script")
out = joinpath(root, "notebook")

function preprocess(s)
    s = "using Pkg; Pkg.activate(\".\"); Pkg.instantiate()\n#-\n" * s
end

for f in ["Project.toml", "Manifest.toml"]
    cp(joinpath(root, f), joinpath(out, f), force=true)
end

for x in readdir(src)
    Literate.notebook(joinpath(src, x), out; preprocess=preprocess, execute=true, documenter=true)
end


text:
using LsqFit, Distributions, DSP, DelimitedFiles

struct input_data
    t::Array{Float64,1}
	DTOF::Array{Float64,1}
	IRF::Array{Float64,1}
	DTOFerrors::Array{Float64,1}
	ρ::Float64
	nmed::Float64
	ndet::Float64
	lambda::Float64
	normpeaks::Bool
end

struct model_params
	#yerrors::Array{Float64,1}          #uncertainties in dependent variable
    model::Function                    #model to be fitted against
	initparams::Array{Float64,1}
	lb::Array{Float64,1}
	ub::Array{Float64,1}
	risefactor::Float64  			#percent of peak to fit on rising edge
	tailfactor::Float64				#perfecnt of peak to fit on failing tail
	dof::Int64
end

struct fitresult
    #Fit inputs
    xdata::Array{Float64,1}     #independent variable
    ydata::Array{Float64,1}     #dependent variable
    yerrors::Array{Float64,1}   #uncertainties in dependent variable
	model::Function             #model to be fitted against

    #Fit results
    residuals::Array{Float64,1}  	 #weighted residuals
    params::Array{Float64,1}    		#best fit parameters for model
	marginerror::Array{Float64,1}        #product of standard error and the critical value of each parameter at a certain significance level (default is 5%) from t-distribution.
    perrors::Array{Float64,1}            #uncertainties in best fit parameters
    dof::Int64                  		#number of degrees of freedom in fit
    chisq::Float64             		 #chi^2 value of fit
    probability::Float64      		  #probability of chi^2 value
end

function DT_model(t, β::Array{Float64,1}, ρ::Float64, ndet::Float64, nmed::Float64)
	#optimzed semi-infinite model

    n::Float64 = nmed/ndet
    μa::Float64 = β[1]
    μsp::Float64 = β[2]
    D::Float64 = 1/3μsp
	ν::Float64 = 29.9792345/nmed

    Rt1 = Array{Float64}(undef, length(t))
    Rt2 = Array{Float64}(undef, length(t))

	if n > 1.0
		A = 504.332889 - 2641.00214n + 5923.699064n^2 - 7376.355814n^3 +
		 5507.53041n^4 - 2463.357945n^5 + 610.956547n^6 - 64.8047n^7
	elseif n < 1.0
		A = 3.084635 - 6.531194n + 8.357854n^2 - 5.082751n^3
	else
		A = 1.0
	end

	zs::Float64 = 1/μsp
	ze::Float64 = 2A*D

	z3m::Float64 = - zs
    z4m::Float64 = 2ze +zs

	Threads.@threads for n in eachindex(t)

   		Rt1[n] = -exp(-(ρ^2/(4D*ν*t[n])) - μa*ν*t[n])
   		Rt1[n] = Rt1[n]/(2*(4π*D*ν)^(3/2)*t[n]^(5/2))

   		Rt2[n] = z3m*exp(-(z3m^2/(4D*ν*t[n]))) - z4m*exp(-(z4m^2/(4D*ν*t[n])))

	end

   Rt = replace!(Rt1.*Rt2, NaN => 0)
   return Rt./maximum(Rt)

end

function conv_DT(t, β::Array{Float64,1}, data::input_data)

	RtDT = Array{Float64}(undef, length(data.t))
	convDT = Array{Float64}(undef, 2*length(data.t) - 1)

	RtDT = DT_model(data.t, β, data.ρ, data.nmed, data.ndet)

	convDT = conv(data.IRF, RtDT)
	convDT = convDT./maximum(convDT)

	tidx = findfirst(x -> x == t[1], data.t)

	if data.normpeaks == true
		peakindDTOF, peakindconvDT = findmax(data.DTOF)[2], findmax(convDT)[2]
		if peakindDTOF > peakindconvDT
			convDT = [zeros(peakindDTOF-peakindconvDT); convDT]
		elseif peakindDTOF < peakindconvDT
			convDT = convDT[peakindconvDT-peakindDTOF:end]
		end
	end

	convDT = convDT[tidx:tidx+length(t)-1]
	return log.(convDT)
end

function getfit(input_data, model_params; alpha=0.1)

	#get fitting window
	counts = input_data.DTOF
	maxvalue, maxindex = findmax(counts)

	ind1 = findfirst(x -> x > maxvalue*model_params.risefactor, counts)
	ind2 =  findlast(x -> x >maxvalue*model_params.tailfactor, counts)

	weights = 1 ./ input_data.DTOFerrors

	fit = curve_fit((t, β) -> conv_DT(t, β, input_data), input_data.t[ind1:ind2], log.(input_data.DTOF[ind1:ind2]), model_params.initparams, lower=model_params.lb, upper=model_params.ub)

	perrors = estimate_errors(fit,alpha)
	marginerror = margin_error(fit,alpha)

	chisq = sum(fit.resid.^2)
	chisqdist = Distributions.Chisq(2)
	probability = Distributions.ccdf(chisqdist,chisq)
	return fitresult(input_data.t[ind1:ind2],log.(input_data.DTOF[ind1:ind2]),input_data.DTOFerrors, DT_model, fit.resid,fit.param,marginerror,perrors,2,chisq,probability)
end

function plotfit(fit::fitresult, data::input_data)
	scatter(fit.xdata, fit.ydata, color="black", label="DTOF", markersize=3, alpha=0.8)

	#scatter(data.t[600:900], log.(data.DTOF[600:900]), color="black", label="DTOF", markersize=3, alpha=0.8)
	#plot!(data.t[600:900], conv_DT(data.t[600:900], fit.params, data), color="red", lw=3, label="DT-fit", alpha=0.8)

	plot!(fit.xdata, conv_DT(fit.xdata, fit.params, data), color="red", lw=3, label="DT-fit", alpha=0.8)
	xlabel!("time (ns)")
	ylabel!("Counts")

end

function plotfit_n(fit::fitresult1, data::input_data)

	xind1 = findmax(fit.ydata)[2]
	xind2 = findmax(conv_DT(fit.xdata, fit.params, data))[2]
	scatter(data.t[600:900].- fit.xdata[xind1], log.(data.DTOF[600:900]), color="black", label="DTOF", markersize=3, alpha=0.8)
	plot!(data.t[600:900].- fit.xdata[xind2], conv_DT(data.t[600:900], fit.params, data), color="red", lw=3, label="DT-fit", alpha=0.8)

	#scatter(fit.xdata .- fit.xdata[xind1], fit.ydata, color="black", label="DTOF", markersize=3, alpha=0.8)
	#plot!(fit.xdata .- fit.xdata[xind2], conv_DT(fit.xdata, fit.params, data), color="red", lw=3, label="DT-fit", alpha=0.8)
	xlabel!("time (ns)")
	ylabel!("Counts")
end

function loadData(filename)
    #loads data from asc
    data_array = readdlm(filename,skipstart=10)
    t = data_array[1:end-1,1]
    Rt = data_array[1:end-1,2]
    (time = convert(Array{Float64},t), counts = convert(Array{Float64},Rt))
end

max_stars_repo_name: sammzer/LightPropagation

max_issues_repo_name: sammzer/LightPropagation

max_forks_repo_name: sammzer/LightPropagation


text:
function Base.show(io::IO, s::Session)
    print(io, "Session services available: $(s.opened_services)")
end

function Base.show(io::IO, s::Service)
    print(io, "Service(\"$(s.name)\")")
end

# Operation is always top level
function Base.show(io::IO, op::Operation)
    println(io, "Operation: $(op.name)")
    println(io)
    println(io, "Request Definition")
    sep(io)
    show(io, op.request_definition)
    println(io)
    if isempty(op.response_definitions)
        println(io, "Response Definition is empty")
    else
        num_response_definitions = length(op.response_definitions)
        for (i, resp_def) in enumerate(op.response_definitions)
            println(io, "Response Definition $i of $num_response_definitions:")
            sep(io)
            show(io, resp_def)
            println(io)
        end
    end
end

Base.show(io::IO, schema::SchemaElementDefinition) = print_schema_tree(io, schema)
Base.show(io::IO, schema_type::AbstractSchemaTypeDefinition) = print_schema_tree(io, schema_type)

function print_schema_tree(io::IO, schema::SchemaElementDefinition, prefix="")
    print(io, "$(prefix)Schema Element `$(schema.name)` - $(schema.status) - ")
    print(io, "Values $(schema.min_values)..")
    if schema.max_values == typemax(schema.max_values)
        print(io, "Unbounded")
    else
        print(io, "$(schema.max_values)")
    end
    println(io)

    if !isempty(schema.alternate_names)
        println(io, "Alternate Names: $(join(schema.alternate_names, ", "))")
    end

    print_schema_tree(io, schema.schema_type, prefix)
end

function print_schema_tree(io::IO, schema::SimpleSchemaTypeDefinition, prefix="")
    println(io, "$(prefix)Simple Schema Type `$(schema.name)` - $(schema.description) - $(schema.status) <: $(schema.datatype)")
end

function print_schema_tree(io::IO, schema::ComplexSchemaTypeDefinition, prefix="")
    println(io, "$(prefix)Complex Schema Type `$(schema.name)` - $(schema.description) - $(schema.status) <: $(schema.datatype)")

    if isempty(schema.elements)
        println(io, "$(prefix)Empty.")
    else
        num_elements = length(schema.elements)
        for (i, element_def) in enumerate(schema.elements)
            println(io, "$(prefix)$(schema.name) element $(i) of $num_elements:")
            print_schema_tree(io::IO, element_def, "    $prefix")
        end
    end
end

function print_schema_tree(io::IO, schema::EnumerationSchemaTypeDefinition, prefix="")
    println(io, "$(prefix)Enumeration Schema Type `$(schema.name)` - $(schema.description) - $(schema.status) <: $(schema.datatype)")
    print_schema_tree(io, schema.enumeration, "    $prefix")
end

function sep(io::IO)
    println(io, "--------------------------")
    println(io)
end

Base.show(io::IO, blpname::BLPName) = print(io, "BLPName(\"$(blpname.symbol)\")")

function print_schema_tree(io::IO, list::BLPConstantList, prefix="")
    println(io, "$(prefix)Enumeration $(list.name) - $(list.description) - $(list.status) - $(list.datatype)")
    if isempty(list.list)
        println(io, "$(prefix)Enumeration is empty")
    else
        num_elements = length(list.list)

        for (i, blp_constant) in enumerate(list.list)
            print(io, "$(prefix)$(list.name) value $(i) of $num_elements: ")
            show(io, blp_constant)
            println(io)
        end
    end
end

function Base.show(io::IO, blp_const::BLPConstant)
    print(io, "$(blp_const.name) - $(blp_const.description) - $(blp_const.status) - $(blp_const.datatype). Julia Value: `$(blp_const.value)`.")
end

Base.show(io::IO, el::Element) = print_element_tree(io, el)

function print_element_tree(io::IO, root::Element{A,D}, prefix="") where {A,D}
    new_prefix = "    " * prefix

    if is_choice_datatype(D)
        println(io, "$prefix$(root.name) <: Element{$A, $D}")
        print_element_tree(io, get_choice(root), new_prefix)
    elseif D == BLPAPI_DATATYPE_SEQUENCE

        println(io, "$prefix$(root.name) <: Element{$A, $D}")

        if A # vector of sequences
            vector = get_element_value(root)
            for element in vector
                print_element_tree(io, element, new_prefix)
            end
        else
            # single sequence
            for child in each_child_element(root)
                print_element_tree(io, child, new_prefix)
            end
        end
    else
        @assert is_simple_datatype(D)
        println(io, "$prefix$(root.name) <: Element{$A, $D} = $(get_element_value(root))")
    end
end

function Base.show(io::IO, event::Event)
    println(io, "Event $(event.event_type)")

    for (i, message) in enumerate(each_message(event))
        println(io, "Message $i")
        println(io, message.correlation_ids)
        element = BLPData.Element(message)
        print_element_tree(io, element, "    ")
    end
end

max_stars_repo_licenses: OpenSSL MIT

max_issues_repo_licenses: OpenSSL MIT

max_forks_repo_licenses: OpenSSL MIT


text:
include("theta.jl")
include("delta.jl")

max_forks_repo_name: ohinder/Everest


text:
using FastGMM
using Distributions
using StatsBase

d = 2
n = 1000000

# start very simple
m = MvNormal([10, -10], [5, 0.2])

X = Float32.(rand(m, n))
ω = Float32.(repeat([1], n))

gmm = fitGMM(X, ω, 1)

K = 4

# build a model to sample from
m = MixtureModel(MvNormal, [
    ([0, 0], [1, 1]),
    ([5, 5], [0.2, 0.4]),
    ([10, 0], [0.8, 0.1]),
    ([1, 24], [0.2, 2])],
    [0.4, 0.2, 0.1, 0.3])

X = Float32.(rand(m, n))
ω = Float32.(repeat([1], n))

gmm = fitGMM(X, ω, K)


text:
export Sites, sitetype, data
export Bit, Spin, Half
export up, down

################################ Site Label ####################################
"""
    SiteLabel
Site Labels, type of each element, like Spin, Bit, etc.
"""
abstract type SiteLabel end

"""
    BitSiteLabel <: SiteLabel
Binary Site, like Spin, Bit, etc.
"""
abstract type BitSiteLabel <: SiteLabel end

"""
    nstate(label)

Returns number of different states of this label.
"""
nstate(::Type{<:BitSiteLabel}) = 2

#################################################################
# Interfaces

"""
    up(BitSiteLabel) -> tag
up tag for this label. e.g. `1` for `Bit`, `0.5` for `Half`.
"""
function up end

"""
    down(BitSiteLabel) -> tag
down tag for this label. e.g. `0` for `Bit`, `-0.5` for `Half`.
"""
function down end

"""
    Spin <: BitSiteLabel
Binary State, Spin. Has two state: -1, 1
"""
abstract type Spin <: BitSiteLabel end

"""
    Bit <: BitSiteLabel
Binary State, Bit. Has two state: 0, 1
"""
abstract type Bit <: BitSiteLabel end

"""
    Half <: BitSiteLabel
Binary State, Half. Has two state: -0.5, 0.5
"""
abstract type Half <: BitSiteLabel end

Base.eltype(::Type{Spin}) = Float64
up(::Type{Spin}) = convert(eltype(Spin), 1.0f0)
down(::Type{Spin}) = convert(eltype(Spin), -1.0f0)

Base.eltype(::Type{Bit}) = Float64
up(::Type{Bit}) = convert(eltype(Bit), 1.0f0)
down(::Type{Bit}) = convert(eltype(Bit), 0.0f0)

Base.eltype(::Type{Half}) = Float64
up(::Type{Half}) = convert(eltype(Half), 0.5f0)
down(::Type{Half}) = convert(eltype(Half), -0.5f0)

################################## Sites #######################################

"""
    AbstractSites{Label, T, N} <: AbstractArray{T, N}
Sites are arrays with certain labels
"""
abstract type AbstractSites{Label, T, N} <: AbstractArray{T, N} end

"""
    sitetype(sites)

get sitetype from subtypes of `AbstractSites`
"""
function sitetype end

"""
    data(sites)
get data of this `sites`
"""
function data end

# TODO: type promote to normal arrays when mixed with normal arrays.

"""
    Sites{L <: SiteLabel, T, N} <: AbstractSites{L, T, N}

Lattice Sites are Julia Arrays with certain Label, it is able to use array
interface.

# Example

    Sites(Bit, 2, 2)
    Sites(Bit, (2, 2))
    rand(Bit, 2, 2)
"""
mutable struct Sites{L <: SiteLabel, T, N} <: AbstractSites{L, T, N}
    data::Array{T, N}
end

Sites{L}(x::Array{T, N}) where {L, T, N} = Sites{L, T, N}(x)

# Enable type inference
Sites(::Type{L}, data::Array{T, N}) where {L, T, N} = Sites{L, T, N}(data)

Sites(::Type{L}, shape::Int...) where L <: SiteLabel = Sites(L, shape)
Sites(::Type{L}, shape::Tuple) where L <: SiteLabel =
    Sites{L, eltype(L), length(shape)}(fill(down(L), shape))

data(s::Sites) = s.data

# use array interface
Base.eltype(x::Sites{L, T}) where {L, T} = T
Base.length(x::Sites) = length(x.data)
Base.ndims(x::Sites) = ndims(x.data)
Base.size(x::Sites) = size(x.data)
Base.size(x::Sites, n::Integer) = size(x.data, n)
Base.axes(x::Sites) = axes(x.data)
Base.axes(x::Sites, d::Integer) = axes(x.data, d)
Base.eachindex(x::Sites) = eachindex(x.data)
Base.stride(x::Sites, k::Integer) = stride(x.data, k)
Base.strides(x::Sites) = strides(x.data)
Base.getindex(x::Sites, index::Integer...) = getindex(x.data, index...)
Base.getindex(x::Sites, index::NTuple{N, T}) where {N, T <: Integer} = getindex(x.data, index...)
Base.setindex!(x::Sites, val, index::Integer...) = setindex!(x.data, val, index...)
Base.setindex!(x::Sites, val, index::NTuple{N, T}) where {N, T <: Integer} = setindex!(x.data, val, index...)

# TODO: value check
Base.getindex(::Type{Bit}, xs...) = Sites{Bit}([xs...])
Base.getindex(::Type{Spin}, xs...) = Sites{Spin}([xs...])
Base.getindex(::Type{Half}, xs...) = Sites{Spin}([xs...])

@inline function Base.copy(b::Sites{L, T, N}) where {L, T, N}
    return Sites{L, T, N}(copy(b.data))
end

sitetype(::Type{Sites{L}}) where L = L
sitetype(x::Sites{L}) where L = L

export reset!, set!

@inline function reset!(b::Sites{L, T, N}) where {L, T, N}
    fill!(b.data, down(L))
    return b
end

@inline function set!(s::Sites{Bit, T, N}, n::Int) where {T, N}
    for i in eachindex(s)
        @inbounds s[i] = ((n>>(i-1)) & 0x1)
    end
    return s
end

@inline function set!(s::Sites{Spin, T, N}, n::Int) where {T, N}
    for i in eachindex(s)
        @inbounds s[i] = 2 * ((n>>(i-1)) & 0x1) - 1
    end
    return s
end

@inline function set!(s::Sites{Half, T, N}, n::Int) where {T, N}
    for i in eachindex(s)
        @inbounds s[i] = ((n>>(i-1)) & 0x1) - 0.5
    end
    return s
end

#######################
# randomize interface
#######################

using Random

function Random.rand!(rng::AbstractRNG, b::Sites{L, T, N}) where {L, T, N}
    rand!(rng, b.data, [up(L), down(L)])
    return b
end

Random.rand!(b::Sites{L, T, N}) where {L <: SiteLabel, T, N} = rand!(GLOBAL_RNG, b)
Random.rand(rng::AbstractRNG, ::Type{L}, shape::Dims) where {L <: SiteLabel} =
    rand!(rng, Sites(L, shape))

####################

@inline function flip!(b::Sites{L, T, N}, index::Integer...) where {L, T, N}
    if b[index...] == up(L)
        b[index...] = down(L)
    else
        b[index...] = up(L)
    end
    return b
end

@inline function randflip!(b::Sites{L, T, N}) where {L, T, N}
    offset = rand(1:length(b))
    return flip!(b, offset)
end

# carry bit
function carrybit!(a::Sites{L}) where L
    @inbounds for i in eachindex(a)
        if a[i] == up(L)
            a[i] = down(L)
        else
            a[i] = up(L)
            break
        end
    end
    a
end

"""
    ndowns(sites)

count the number of down tags, e.g number of `0`s
"""
function ndowns(s::Sites{L}) where {L <:BitSiteLabel}
    count = 0
    @inbounds for i in eachindex(s)
        if s[i] == down(L)
            count += 1
        end
    end
    count
end

function nups(s::Sites{L}) where {L <:BitSiteLabel}
    count = 0
    @inbounds for i in eachindex(s)
        if s[i] == up(L)
            count += 1
        end
    end
    count
end

function n_downs_ups(s::Sites{L}) where {L <:BitSiteLabel}
    count_down = 0
    count_up   = 0
    @inbounds for i in eachindex(s)
        s[i] == up(L) ? (count_up += 1) : (count_down += 1)
    end
    count_down, count_up
end

# TODO: whether this should have side effect?
# will side effect affects performance (extra instance)?
function Base.:(<<)(a::Sites, b::Int)
    for i = 1:b
        carrybit!(a)
    end
    return a
end

##############
# Conversions
##############

import Base: convert

function convert(::Type{Integer}, x::Sites{L, T, N}) where {L, T, N} end

for IntType in (:Int8, Int16, Int32, Int64, Int128, BigInt)
@eval begin
        function convert(::Type{$IntType}, x::Sites{Bit, T, N}) where {T, N}
            if sizeof($IntType) * 8 < length(x)
                throw(Compat.InexactError(:convert, $IntType, x))
            end

            sum($IntType(each) << (i-1) for (i, each) in enumerate(x))
        end

        function convert(::Type{$IntType}, x::Sites{Spin, T, N}) where {T, N}
            if sizeof($IntType) * 8 < length(x)
                throw(Compat.InexactError(:convert, $IntType, x))
            end

            sum($IntType(div(each+1, 2)) << (i-1) for (i, each) in enumerate(x))
        end

        function convert(::Type{$IntType}, x::Sites{Half, T, N}) where {T, N}
            if sizeof($IntType) * 8 < length(x)
                throw(Compat.InexactError(:convert, $IntType, x))
            end

            sum($IntType(each+0.5) << (i-1) for (i, each) in enumerate(x))
        end
    end
end

max_stars_repo_name: Roger-luo/QuHamiltonian

max_issues_repo_name: Roger-luo/QuHamiltonian

max_forks_repo_name: Roger-luo/QuHamiltonian


text:
using MatrixDepot
using Base.Test

include("test_cauchy.jl")
include("test_circul.jl")
include("test_hadamard.jl")
include("test_hilb.jl")
include("test_dingdong.jl")
include("test_frank.jl")
include("test_invhilb.jl")
include("test_forsythe.jl")
include("test_magic.jl")
include("test_grcar.jl")
include("test_triw.jl")
include("test_moler.jl")
include("test_pascal.jl")
include("test_kahan.jl")
include("test_pei.jl")
include("test_vand.jl")
include("test_invol.jl")
include("test_chebspec.jl")
include("test_lotkin.jl")
include("test_clement.jl")
include("test_fiedler.jl")
include("test_minij.jl")
include("test_binomial.jl")
include("test_tridiag.jl")
include("test_lehmer.jl")
include("test_parter.jl")
include("test_chow.jl")
include("test_randcorr.jl")
include("test_poisson.jl")
include("test_neumann.jl")
include("test_rosser.jl")
include("test_sampling.jl")
include("test_wilkinson.jl")
include("test_rando.jl")
include("test_randsvd.jl")
include("test_rohess.jl")
include("test_kms.jl")
include("test_wathen.jl")
include("test_oscillate.jl")
include("test_toeplitz.jl")
include("test_hankel.jl")
include("test_golub.jl")
include("test_companion.jl")
include("test_erdrey.jl")
include("test_gilbert.jl")
include("test_smallworld.jl")

tests = [
         "clean",
         "common",
         "include_generator",
         "download",
         "number",
         "property",
         "regu",
         ]

for t in tests
    tp = joinpath(Pkg.dir("MatrixDepot"), "test", "$(t).jl")
    println("running $(tp) ...")
    include(tp)
end

println("Success in all tests.")


text:
abstract type AbstractXAIMethod end
# All analyzers are implemented such that they return an array of explanations:
#   (method::AbstractXAIMethod)(input, ns::AbstractNeuronSelector)::Vector{Explanation}

const BATCHDIM_MISSING = ArgumentError(
    """The input is a 1D vector and therefore missing the required batch dimension.
    Call analyze with the keyword argument add_batch_dim=false."""
)

"""
    analyze(input, method)
    analyze(input, method, neuron_selection)

Return raw classifier output and explanation.
If `neuron_selection` is specified, the explanation will be calculated for that neuron.
Otherwise, the output neuron with the highest activation is automatically chosen.

## Keyword arguments
- `add_batch_dim`: add batch dimension to the input without allocating. Default is `false`.
"""
function analyze(
    input::AbstractArray{<:Real},
    method::AbstractXAIMethod,
    neuron_selection::Union{Integer,Tuple{<:Integer}};
    kwargs...,
)
    return _analyze(input, method, IndexSelector(neuron_selection); kwargs...)
end

function analyze(input::AbstractArray{<:Real}, method::AbstractXAIMethod; kwargs...)
    return _analyze(input, method, MaxActivationSelector(); kwargs...)
end

function (method::AbstractXAIMethod)(
    input::AbstractArray{<:Real},
    neuron_selection::Union{Integer,Tuple{<:Integer}};
    kwargs...,
)
    return _analyze(input, method, IndexSelector(neuron_selection); kwargs...)
end
function (method::AbstractXAIMethod)(input::AbstractArray{<:Real}; kwargs...)
    return _analyze(input, method, MaxActivationSelector(); kwargs...)
end

# lower-level call to method
function _analyze(
    input::AbstractArray{T,N},
    method::AbstractXAIMethod,
    sel::AbstractNeuronSelector;
    add_batch_dim::Bool=false,
    kwargs...,
) where {T<:Real,N}
    if add_batch_dim
        return method(batch_dim_view(input), sel; kwargs...)
    end
    N < 2 && throw(BATCHDIM_MISSING)
    return method(input, sel; kwargs...)
end

# for convenience, the anaylyzer can be called directly

# Explanations and outputs are returned in a wrapper.
# Metadata such as the analyzer allows dispatching on functions like `heatmap`.
struct Explanation{A,O,I,L}
    attribution::A
    output::O
    neuron_selection::I
    analyzer::Symbol
    layerwise_relevances::L
end


text:
using DataArrays, DataFrames, YAML, Interpolations

@enum(INDEXFN,
  Sellmeier =1,
  Sellmeier2 =2,
  Polynomial =3 ,
  RFInfo =4,
  Cauchy =5,
  Gases =6,
  Herzberger =7,
  Retro =8,
  Exotic =9
)

const λ_i = (0.36501)	#	Hg	UV
const λ_h = (0.40466)	#	Hg	violet
const λ_g = (0.43584)	#	Hg	blue
const λ_Fp = (0.47999)	#	Cd	blue
const λ_F = (0.48613)	#	H	blue
const λ_e = (0.54607)	#	Hg	green
const λ_d = (0.58756)	#	He	yellow
const λ_D = (0.5893)	#	Na	yellow
const λ_Cp = (0.64385)	#	Cd	red
const λ_C = (0.65627)	#	H	red
const λ_r = (0.70652)	#	He	red
const λ_Ap = (0.7682)	#	K	IR
const λ_s = (0.85211)	#	Cs	IR
const λ_t = (0.101398)	#	Hg	IR

_fraunhofer_lines = Dict(
"i" => ("Hg", 0.36501, "blue"),
"h" => ("Hg", 0.40466, "blue"),
"g" => ("Hg", 0.43585, "blue"),
"F'" => ("Cd", 0.47999, "green"),
"F" => ("Hβ", 0.486134, "infrared"),
"e" => ("Hg", 0.546073, "infrared"),
"d" => ("He", 0.58756, "infrared"),
"D" => ("Na", 0.5893, "red"),
"C'" => ("Cd", 0.64385, "red"),
"C" => ("Hα", 0.656281, "red"),
"r" => ("He", 0.70652, "ultraviolet"),
"A'" => ("K", 0.7682, "violet"),
"s" => ("Cs", 0.85211, "yellow"),
"t" => ("Hg", 1.01398, "yellow"),
"t2" => ("Ni", 0.299444, ""),
"T" => ("Fe", 0.302108, ""),
"P" => ("Ti+", 0.336112, ""),
"N" => ("Fe", 0.358121, ""),
"L" => ("Fe", 0.382044, ""),
"K" => ("Ca+", 0.3933666, ""),
"H" => ("Ca+", 0.396847, ""),
"h2" => ("Hδ", 0.410175, ""),
"G2" => ("Ca", 0.430774, ""),
"G" => ("Fe", 0.43079, ""),
"G'" => ("Hγ", 0.434047, ""),
"e2" => ("Fe", 0.438355, ""),
"d2" => ("Fe", 0.466814, ""),
"c" => ("Fe", 0.495761, ""),
"b4" => ("Mg", 0.516733, ""),
"b3" => ("Fe", 0.516891, ""),
"b2" => ("Mg", 0.51727, ""),
"b1" => ("Mg", 0.518362, ""),
"E2" => ("Fe", 0.527039, ""),
"D3" => ("He", 0.5875618, ""),
"D2" => ("Na", 0.588995, ""),
"D1" => ("Na", 0.589592, ""),
"a" => ("O2", 0.627661, ""),
"B" => ("O2", 0.686719, ""),
"A" => ("O2", 0.75937, ""),
"Z" => ("O2", 0.822696, ""),
"y" => ("O2", 0.898765, "")
)

fraunhofer_λum(line::AbstractString) = get(_fraunhofer_lines,line,NaN)[2]
fraunhofer_source(line::AbstractString) = get(_fraunhofer_lines,line,"")[1]
fraunhofer_color(line::AbstractString) = get(_fraunhofer_lines,line,"")[3]
fraunhofer_standards() = keys(filter((k,v)->!isempty(v[3]), _fraunhofer_lines))

function refindex{T<:Real}(::Type{Val{Sellmeier}}, C::Array{T,1}, λum::T)
  local N = length(C)
  @assert isodd(N) && N>=3 "index equation requires an odd number of coefficients"
  local λsq::T = λum^2
  local nsq::T = one(T) + C[1]
  for i = 2:2:N-1
    nsq += C[i] * λsq / (λsq - C[i+1]^2)
  end
  sqrt(nsq)
end

function refindex{T<:Real}(::Type{Val{Sellmeier2}}, C::Array{T,1}, λum::T)
  local N = length(C)
  @assert isodd(N) && N>=3 "index equation requires an odd number of coefficients"
  local λsq::T = λum^2
  local nsq::T = one(T) + C[1]
  for i = 2:2:N-1
    nsq += C[i] * λsq / (λsq - C[i+1])
  end
  sqrt(nsq)
end

function refindex{T<:Real}(::Type{Val{Polynomial}}, C::Array{T,1}, λum::T)
  local N = length(C)
  @assert isodd(N) && N>=3 "index equation requires an odd number of coefficients"
  local nsq::T = C[1]
  for i = 2:2:N-1
    nsq += C[i] * λum^C[i+1]
  end
  sqrt(nsq)
end

function refindex{T<:Real}(::Type{Val{RFInfo}}, C::Array{T,1}, λum::T)
  local N = length(C)
  @assert isodd(N) && N>=11 "index equation requires an odd number of coefficients ≥ 11"
  local λsq::T = λum^2
  local nsq::T = C[1]
  for i = 2:4:9
    nsq += C[i] * λum^C[i+2] / (λsq - C[i+3]^C[i+4])
  end
  for i = 10:2:N-1
    nsq += C[i]*λum^C[i+1]
  end
  sqrt(nsq)
end

function refindex{T<:Real}(::Type{Val{Cauchy}}, C::Array{T,1}, λum::T)
  local N = length(C)
  @assert isodd(N) && N>=3 "index equation requires an odd number of coefficients"
  local n::T = C[1]
  for i = 2:2:N-1
    n += C[i]*λum^C[i+1]
  end
  n
end

function refindex{T<:Real}(::Type{Val{Gases}}, C::Array{T,1}, λum::T)
  local N = length(C)
  @assert isodd(N) && N>=3 "index equation requires an odd number of coefficients"
  local n::T = one(T) + C[1]
  local invsqrtλ::T = λum^-2
  for i = 2:2:N-1
    n += C[i] / (C[i+1] - invsqrtλ)
  end
  n
end

function refindex{T<:Real}(::Type{Val{Herzberger}}, C::Array{T,1}, λum::T)
  local N = length(C)
  @assert N>=6 "index equation requires 6 coefficients"
  local n::T = C[1]
  local λsq::T = λum^2
  n += C[2] / (λsq - 0.028)
  n += C[3]*(1/(λsq - 0.028))^2
  for i = 4::N
    n += C[i] * λum^(2*(i-3))
  end
  n
end

function refindex{T<:Real}(::Type{Val{Retro}}, C::Array{T,1}, λum::T)
  local N = length(C)
  @assert N>=4 "index equation requires 4 coefficients"
  local λsq::T = λum^2
  local r::T = C[1] + C[2]*λsq/(λsq-C[3]) + C[4]*λsq
  sqrt(-2r - 1) / sqrt(r - 1)
end

function refindex{T<:Real}(::Type{Val{Exotic}}, C::Array{T,1}, λum::T)
  local N = length(C)
  @assert N>=6 "index equation requires 6 coefficients"
  local λsq::T = λum^2
  local nsq::T = C[1] + C[2]/(λsq-C[3]) + C[4]*(λum - C[5])/((λum - C[5])^2 + C[6])
  sqrt(nsq)
end

fusedSiCstr = "0 0.6961663 0.0684043 0.4079426 0.1162414 0.8974794 9.896161"
fusedSiC = readdlm(IOBuffer(fusedSiCstr),' ', Float64)

type Material
  name::AbstractString
  nD::Float64
  VD::Float64
  fnum::Int
  range::Tuple{Float64, Float64}
  C::Array{Float64,1}
  Ntable::Array{Float64,2}
  Ktable::Array{Float64,2}
  Ninterp::Nullable{Interpolations.GriddedInterpolation}
  Kinterp::Nullable{Interpolations.GriddedInterpolation}
  info::Dict{AbstractString,Any}
  comments::AbstractString
  references::AbstractString
  Material(name) = new(name)
end

function Material(name::AbstractString, path::AbstractString, file::AbstractString)
  matd = YAML.load_file(joinpath(path,file))
  m = Material(name)
  m.name = name
  m.nD = NaN
  m.VD = NaN
  m.fnum = -1
  m.range = (NaN, NaN,)
  m.C = Array(Float64,0)
  m.Ntable = Array(Float64,0,0)
  m.Ktable = Array(Float64,0,0)
  m.Ninterp = Nullable()
  m.Kinterp = Nullable()
  if haskey(matd, "DATA")
    for data in matd["DATA"]
      typestr = lowercase(get(data, "type", ""))
      if startswith(typestr,"formula") && haskey(data, "coefficients")
        m.fnum = parse(Int,split(typestr)[end])
        m.C = vec(readdlm(IOBuffer(data["coefficients"]), ' ', Float64)[1,:])
        range = haskey(data, "range") ? readdlm(IOBuffer(data["range"]),' ', Float64) : [NaN, NaN]
        m.range = tuple(range...)
      elseif startswith(typestr,"tabulated") && haskey(data,"data")
        tabnk = lowercase(split(typestr)[end])
        hasn = search(tabnk, 'n') > 0
        hask = search(tabnk, 'k') > 0
        tab = readdlm(IOBuffer(data["data"]), ' ', Float64)
        if any(isnan, m.range)
          m.range = extrema(tab[:,1])
        end
        if hasn && hask
          m.Ntable = tab[:,[1,2]]
          m.Ktable = tab[:,[1,3]]
        elseif hasn
          m.Ntable = tab
        else
          m.Ktable = tab
        end
      end
    end
  end
  if length(m.Ntable) > 0
    knots = (vec(m.Ntable[:,1]),)  # tuple of x coords
    vals = vec(m.Ntable[:,2])      # f(x)
    m.Ninterp = Nullable(interpolate(knots, vals, Gridded(Linear())))
  end
  if length(m.Ktable) > 0
    knots = (vec(m.Ktable[:,1]),)  # tuple of x coords
    vals = vec(m.Ktable[:,2])      # f(x)
    m.Kinterp = Nullable(interpolate(knots, vals, Gridded(Linear())))
  end

  m.info = get(matd, "INFO", Dict())
  m.comments = get(matd, "COMMENTS", "")
  m.references = get(matd, "REFERENCES","")
  m.nD = get(m.info, "n_d", NaN)
  m.VD = get(m.info, "V_d", NaN)
  #println(">> range: $(m.range[1]) to $(m.range[2])")
  if (isnan(m.nD) || isnan(m.VD)) && m.range[1] <= λ_F && λ_C <= m.range[2]
    nD = refindex(m, λ_D)
    nC = refindex(m, λ_C)
    nF = refindex(m, λ_F)
    #nF = refindex(Val{INDEXFN(m.fnum)}, m.C, λ_F)
    #println(">>Calculated nD=$nD, nC=$nC, nF=$nF")
    m.nD = nD
    m.VD = (m.nD - 1) / (nF - nC)
  end
  m
end

function refindex{T<:Real}(m::Material, λum::T)
  if 1 <= m.fnum <= 9 && length(m.C) > 0
    refindex(Val{INDEXFN(m.fnum)}, m.C, λum)
  elseif !isnull(m.Ninterp) && m.range[1] <= λum <= m.range[2]
    get(m.Ninterp)[λum]
  else
    NaN
  end
end

dbfolder = Pkg.dir("GlassCatalog","assets","refractiveindex.info-database","database")

println("=== L-LAH86 ===")
lah86 = Material("L-LAH86", dbfolder, "glass/ohara/L-LAH86.yml")
dump(lah86)
println()

println("=== BSC7 ===")
bsc7 = Material("BSC7", dbfolder, "glass/hoya/BSC7.yml")
dump(bsc7)
println()

println("=== GaInP ===")
gainp = Material("GaInP", dbfolder, "other/mixed crystals/GaInP/Schubert.yml")
dump(gainp)
println()


text: # import BuildExecutable


text:
struct SimpleEM <: DiffEqBase.AbstractSDEAlgorithm end
export SimpleEM

@muladd function DiffEqBase.solve(prob::SDEProblem{uType,tType,false},alg::SimpleEM,args...;
                          dt = error("dt required for SimpleEM")) where {uType,tType}

  f = prob.f
  g = prob.g
  u0 = prob.u0
  tspan = prob.tspan
  p = prob.p

  @inbounds begin
      n = Int((tspan[2] - tspan[1])/dt) + 1
      u = [u0 for i in 1:n]
      t = [tspan[1] + i*dt for i in 0:n-1]
      sqdt = sqrt(dt)
  end

  @inbounds for i in 2:n
      uprev = u[i-1]
      tprev = t[i-1]
      if u0 isa Number
          u[i] = uprev + f(uprev,p,tprev)*dt + sqdt*g(uprev,p,tprev)*randn(typeof(u0))
      else
          u[i] = uprev + f(uprev,p,tprev)*dt + sqdt*g(uprev,p,tprev).*randn(typeof(u0))
      end
  end

  sol = DiffEqBase.build_solution(prob,alg,t,u,
                       calculate_error = false)
end

@muladd function DiffEqBase.solve(prob::SDEProblem{uType,tType,true},alg::SimpleEM,args...;
                          dt = error("dt required for SimpleEM")) where {uType,tType}

  f = prob.f
  g = prob.g
  u0 = prob.u0
  tspan = prob.tspan
  p = prob.p
  ftmp = zero(u0)
  gtmp = DiffEqBase.is_diagonal_noise(prob) ? zero(u0) : zero(prob.noise_rate_prototype)
  gtmp2 = DiffEqBase.is_diagonal_noise(prob) ? nothing : zero(u0)
  dW = DiffEqBase.is_diagonal_noise(prob) ? zero(u0) : false .* prob.noise_rate_prototype[1,:]

  @inbounds begin
      n = Int((tspan[2] - tspan[1])/dt) + 1
      u = [copy(u0) for i in 1:n]
      t = [tspan[1] + i*dt for i in 0:n-1]
      sqdt = sqrt(dt)
  end

  @inbounds for i in 2:n
      uprev = u[i-1]
      tprev = t[i-1]
      f(ftmp,uprev,p,tprev)
      g(gtmp,uprev,p,tprev)
      @. dW = randn(eltype(dW))

      if DiffEqBase.is_diagonal_noise(prob)
          DiffEqBase.@.. u[i] = uprev + ftmp*dt + sqdt*gtmp*dW
      else
          mul!(gtmp2,gtmp,dW)
          DiffEqBase.@.. u[i] = uprev + ftmp*dt + sqdt*gtmp2
      end
  end

  sol = DiffEqBase.build_solution(prob,alg,t,u,
                       calculate_error = false)
end


text:
module HDModule

using HierarchicalStateMachineModule: Event,State
using FieldModule: Field,DirectionIndicators,maximum,set!,fill!,+,invert,repeat,get
using CoordsModule: Coords, DirectionIndicator,LatLonSectionCoords,
      is_ocean, is_outflow, is_truesink, is_lake
using GridModule: Grid, for_all,find_downstream_coords,for_section_with_line_breaks
using UserExceptionModule: UserError
import HierarchicalStateMachineModule: handle_event
using InteractiveUtils
using Printf: @printf

struct RunHD <: Event end

abstract type PrognosticFields <: State end

struct RiverParameters
  flow_directions::DirectionIndicators
  river_reservoir_nums::Field{Int64}
  overland_reservoir_nums::Field{Int64}
  base_reservoir_nums::Field{Int64}
  river_retention_coefficients::Field{Float64}
  overland_retention_coefficients::Field{Float64}
  base_retention_coefficients::Field{Float64}
  landsea_mask::Field{Bool}
  cascade_flag::Field{Bool}
  grid::Grid
  step_length::Float64
  function RiverParameters(flow_directions::DirectionIndicators,
                           river_reservoir_nums::Field{Int64},
                           overland_reservoir_nums::Field{Int64},
                           base_reservoir_nums::Field{Int64},
                           river_retention_coefficients::Field{Float64},
                           overland_retention_coefficients::Field{Float64},
                           base_retention_coefficients::Field{Float64},
                           landsea_mask::Field{Bool},
                           grid::Grid,
                           day_length::Float64,
                           step_length::Float64)
    cascade_flag::Field{Bool} = invert(landsea_mask)
    river_retention_coefficients = river_retention_coefficients *
                                   (day_length/step_length)
    overland_retention_coefficients = overland_retention_coefficients *
                                      (day_length/step_length)
    base_retention_coefficients = base_retention_coefficients *
                                  (day_length/step_length)
    return new(flow_directions,river_reservoir_nums,overland_reservoir_nums,
               base_reservoir_nums,river_retention_coefficients,
               overland_retention_coefficients,base_retention_coefficients,
               landsea_mask,cascade_flag,grid,step_length)
  end
end

RiverParameters(flow_directions::DirectionIndicators,
                river_reservoir_nums::Field{Int64},
                overland_reservoir_nums::Field{Int64},
                base_reservoir_nums::Field{Int64},
                river_retention_coefficients::Field{Float64},
                overland_retention_coefficients::Field{Float64},
                base_retention_coefficients::Field{Float64},
                landsea_mask::Field{Bool},
                grid::Grid) = RiverParameters(flow_directions::DirectionIndicators,
                                              river_reservoir_nums::Field{Int64},
                                              overland_reservoir_nums::Field{Int64},
                                              base_reservoir_nums::Field{Int64},
                                              river_retention_coefficients::Field{Float64},
                                              overland_retention_coefficients::Field{Float64},
                                              base_retention_coefficients::Field{Float64},
                                              landsea_mask::Field{Bool},
                                              grid::Grid,86400.0,86400.0)

mutable struct RiverPrognosticFields
  runoff::Field{Float64}
  drainage::Field{Float64}
  lake_evaporation::Field{Float64}
  river_inflow::Field{Float64}
  base_flow_reservoirs::Array{Field{Float64},1}
  overland_flow_reservoirs::Array{Field{Float64},1}
  river_flow_reservoirs::Array{Field{Float64},1}
  water_to_ocean::Field{Float64}
  function RiverPrognosticFields(river_parameters::RiverParameters)
    runoff = Field{Float64}(river_parameters.grid,0.0)
    drainage = Field{Float64}(river_parameters.grid,0.0)
    lake_evaporation = Field{Float64}(river_parameters.grid,0.0)
    river_inflow = Field{Float64}(river_parameters.grid,0.0)
    base_flow_reservoirs =      repeat(Field{Float64}(river_parameters.grid,0.0),
                                       maximum(river_parameters.base_reservoir_nums))
    overland_flow_reservoirs =  repeat(Field{Float64}(river_parameters.grid,0.0),
                                       maximum(river_parameters.overland_reservoir_nums))
    river_flow_reservoirs =     repeat(Field{Float64}(river_parameters.grid,0.0),
                                       maximum(river_parameters.river_reservoir_nums))
    water_to_ocean = Field{Float64}(river_parameters.grid,0.0)
    return new(runoff,drainage,lake_evaporation,river_inflow,base_flow_reservoirs,
               overland_flow_reservoirs,river_flow_reservoirs,water_to_ocean)
  end
end

mutable struct RiverDiagnosticFields
  runoff_to_rivers::Field{Float64}
  drainage_to_rivers::Field{Float64}
  river_outflow::Field{Float64}
  function RiverDiagnosticFields(river_parameters::RiverParameters)
    runoff_to_rivers::Field{Float64} =
      Field{Float64}(river_parameters.grid,0.0)
    drainage_to_rivers::Field{Float64} =
      Field{Float64}(river_parameters.grid,0.0)
    river_outflow::Field{Float64} =
      Field{Float64}(river_parameters.grid,0.0)
    new(runoff_to_rivers,drainage_to_rivers,river_outflow)
  end
end

mutable struct RiverDiagnosticOutputFields
  cumulative_river_flow::Field{Float64}
  function RiverDiagnosticOutputFields(river_parameters::RiverParameters)
    cumulative_river_flow::Field{Float64} =
      Field{Float64}(river_parameters.grid,0.0)
    new(cumulative_river_flow)
  end
end

get_river_parameters(obj::T) where {T <: PrognosticFields} =
  obj.river_parameters::RiverParameters

get_river_fields(obj::T) where {T <: PrognosticFields} =
  obj.river_fields::RiverPrognosticFields

get_river_diagnostic_fields(obj::T) where {T <: PrognosticFields} =
  obj.river_diagnostic_fields::RiverDiagnosticFields

get_river_diagnostic_output_fields(obj::T) where {T <: PrognosticFields} =
  obj.river_diagnostic_output_fields::RiverDiagnosticOutputFields

get_using_lakes(obj::T) where {T <: PrognosticFields} =
  obj.using_lakes::Bool

struct RiverPrognosticFieldsOnly <: PrognosticFields
  river_parameters::RiverParameters
  river_fields::RiverPrognosticFields
  river_diagnostic_fields::RiverDiagnosticFields
  river_diagnostic_output_fields::RiverDiagnosticOutputFields
  using_lakes::Bool
  RiverPrognosticFieldsOnly(river_parameters::RiverParameters,
                            river_fields::RiverPrognosticFields) =
    new(river_parameters,river_fields,RiverDiagnosticFields(river_parameters),
        RiverDiagnosticOutputFields(river_parameters),false)
end

function handle_event(prognostic_fields::PrognosticFields,
                      run_hd::RunHD)
  river_parameters::RiverParameters = get_river_parameters(prognostic_fields)
  river_fields::RiverPrognosticFields = get_river_fields(prognostic_fields)
  river_diagnostic_fields::RiverDiagnosticFields =
    get_river_diagnostic_fields(prognostic_fields)
  using_lakes::Bool =
              get_using_lakes(prognostic_fields)
  local lake_water_from_ocean::Field{Float64}
  if using_lakes
    lake_water_input::Field{Float64},lake_water_from_ocean =
      water_from_lakes(prognostic_fields,river_parameters.step_length)
    for_all(river_parameters.grid) do coords::Coords
      cell_lake_water_input::Float64 = lake_water_input(coords)
      if river_parameters.cascade_flag(coords)
        set!(river_fields.river_inflow,coords,
             get(river_fields.river_inflow,coords)+cell_lake_water_input)
      else
        set!(river_diagnostic_fields.river_outflow,
             coords,cell_lake_water_input)
      end
    end
  end
  cascade(river_fields.overland_flow_reservoirs,
          river_fields.runoff,
          river_diagnostic_fields.runoff_to_rivers,
          river_parameters.overland_retention_coefficients,
          river_parameters.base_reservoir_nums,
          river_parameters.cascade_flag,river_parameters.grid,
          river_parameters.step_length)
  cascade(river_fields.base_flow_reservoirs,
          river_fields.drainage,
          river_diagnostic_fields.drainage_to_rivers,
          river_parameters.base_retention_coefficients,
          river_parameters.base_reservoir_nums,
          river_parameters.cascade_flag,river_parameters.grid,
          river_parameters.step_length)
  cascade(river_fields.river_flow_reservoirs,
          river_fields.river_inflow,
          river_diagnostic_fields.river_outflow,
          river_parameters.river_retention_coefficients,
          river_parameters.river_reservoir_nums,
          river_parameters.cascade_flag,river_parameters.grid,
          river_parameters.step_length)
  fill!(river_fields.river_inflow,0.0)
  route(river_parameters.flow_directions,
        river_diagnostic_fields.river_outflow+
        river_diagnostic_fields.runoff_to_rivers+
        river_diagnostic_fields.drainage_to_rivers,
        river_fields.river_inflow,
        river_parameters.grid)
  fill!(river_diagnostic_fields.river_outflow,0.0)
  fill!(river_diagnostic_fields.runoff_to_rivers,0.0)
  fill!(river_diagnostic_fields.drainage_to_rivers,0.0)
  for_all(river_parameters.grid) do coords::Coords
            flow_direction::DirectionIndicator =
              get(river_parameters.flow_directions,coords)
            if is_ocean(flow_direction) || is_outflow(flow_direction) ||
               is_truesink(flow_direction)
                set!(river_fields.water_to_ocean,coords,
                      get(river_fields.river_inflow,coords) +
                      get(river_fields.runoff,coords) +
                      get(river_fields.drainage,coords))
                set!(river_fields.river_inflow,coords,0.0)
            elseif using_lakes && is_lake(flow_direction)
                water_to_lakes(prognostic_fields,coords,
                               get(river_fields.river_inflow,coords) +
                               get(river_fields.runoff,coords) +
                               get(river_fields.drainage,coords) -
                               get(river_fields.lake_evaporation,coords),
                               river_parameters.step_length)
                set!(river_fields.water_to_ocean,coords,
                     -1.0*get(lake_water_from_ocean,coords))
                set!(river_fields.river_inflow,coords,0.0)
            end
          end
  fill!(river_fields.runoff,0.0)
  fill!(river_fields.drainage,0.0)
  return prognostic_fields
end

function water_to_lakes(prognostic_fields::PrognosticFields,coords::Coords,inflow::Float64,
                        step_length::Float64)
  throw(UserError())
end

function water_from_lakes(prognostic_fields::PrognosticFields,step_length::Float64)
  throw(UserError())
end

function cascade(reservoirs::Array{Field{Float64},1},
                 inflow::Field{Float64},
                 outflow::Field{Float64},
                 retention_coefficients::Field{Float64},
                 reservoir_nums::Field{Int64},
                 cascade_flag::Field{Bool},grid::Grid,
                 step_length::Float64)
  for_all(grid) do coords::Coords
    if get(cascade_flag,coords)
      cascade_kernel(coords,
                    reservoirs,
                    inflow,
                    outflow,
                    retention_coefficients,
                    reservoir_nums,
                    step_length)
    end
    return
  end
end

function cascade(reservoirs::Array{Array{Field{Float64},1},1},
                 inflow::Array{Field{Float64},1},
                 outflow::Array{Field{Float64},1},
                 retention_coefficients::Array{Field{Float64},1},
                 reservoir_nums::Array{Field{Int64},1},
                 cascade_flag::Field{Bool},grid::Grid,
                 cascade_num::Int64,
                 step_length::Float64)
  for_all(grid) do coords::Coords
    if get(cascade_flag,coords)
      for i = 1:cascade_num
        reservoirs_i::Array{Field{Float64},1} = reservoirs[i]
        inflow_i::Field{Float64} = inflow[i]
        outflow_i::Field{Float64} = outflow[i]
        retention_coefficients_i::Field{Float64} = retention_coefficients[i]
        reservoir_nums_i::Field{Int64} = reservoir_nums[i]
        cascade_kernel(coords,
                       reservoirs_i,
                       inflow_i,
                       outflow_i,
                       retention_coefficients_i,
                       reservoir_nums_i,
                       step_length)
      end
    end
    return
  end
end

function cascade_kernel(coords::Coords,
                        reservoirs::Array{Field{Float64},1},
                        inflow::Field{Float64},
                        outflow::Field{Float64},
                        retention_coefficients::Field{Float64},
                        reservoir_nums::Field{Int64},
                        step_length::Float64)
  flow::Float64 = get(inflow,coords)*step_length
  for i = 1:(get(reservoir_nums,coords))
    reservoir::Field{Float64} = reservoirs[i]
    new_reservoir_value::Float64 = get(reservoir,coords) + flow
    flow = new_reservoir_value/(get(retention_coefficients,coords)+1.0)
    set!(reservoir,coords,new_reservoir_value - flow)
  end
  flow /= step_length
  set!(outflow,coords,flow)
  return
end

# macro fuse_cascades(exprs...)
#   cascade_num = 0
#   for expr in exprs
#     if expr[1] == :call && expr[1] == :cascade
#       cascade_num += 1
#       # figure out what to put in arrays
#     end
#   end
#   return quote
#             reservoirs_macro = Array{Field{Float64},1}[]
#             inflow_macro = Field{Float64}[]
#             outflow_macro = Field{Float64}[]
#             retention_coefficients_macro = Field{Float64}[]
#             reservoir_nums_macro = Field{Int64}[]
#             cascade_num_macro = $cascade_num
#             cascade(reservoirs_macro,
#                     inflow_macro,
#                     outflow_macro,
#                     retention_coefficients_macro,
#                     reservoir_nums_macro,
#                     cascade_flag,grid
#                     cascade_num_macro)
#          end
# end

function route(flow_directions::DirectionIndicators,
               flow_in::Field{Float64},
               flow_out::Field{Float64},
               grid::Grid)
  for_all(grid) do coords::Coords
    flow_in_local::Float64 = flow_in(coords)
    if flow_in_local != 0.0
      flow_direction::DirectionIndicator =
        get(flow_directions,coords)
      if is_truesink(flow_direction) || is_lake(flow_direction) ||
         is_ocean(flow_direction) || is_outflow(flow_direction)
        flow_in_local += get(flow_out,coords)
        set!(flow_out,coords,flow_in_local)
      else
        new_coords::Coords =
          find_downstream_coords(grid,
                                flow_direction,
                                coords)
        flow_in_local += get(flow_out,new_coords)
        set!(flow_out,new_coords,flow_in_local)
      end
    end
    return
  end
end

struct SetDrainage <: Event
  new_drainage::Field{Float64}
end

struct SetRunoff <: Event
  new_runoff::Field{Float64}
end

struct SetLakeEvaporation <: Event
  new_lake_evaporation::Field{Float64}
end

function handle_event(prognostic_fields::PrognosticFields,
                      set_drainage::SetDrainage)
  river_fields::RiverPrognosticFields = get_river_fields(prognostic_fields)
  river_fields.drainage = set_drainage.new_drainage
  return prognostic_fields
end

function handle_event(prognostic_fields::PrognosticFields,
                      set_runoff::SetRunoff)
  river_fields::RiverPrognosticFields = get_river_fields(prognostic_fields)
  river_fields.runoff = set_runoff.new_runoff
  return prognostic_fields
end

function handle_event(prognostic_fields::PrognosticFields,
                      set_lake_evaporation::SetLakeEvaporation)
  river_fields::RiverPrognosticFields = get_river_fields(prognostic_fields)
  river_fields.lake_evaporation = set_lake_evaporation.new_lake_evaporation
  return prognostic_fields
end

struct PrintResults <: Event
  timestep::Int64
end

struct PrintSection <: Event end

function handle_event(prognostic_fields::PrognosticFields,
                      print_results::PrintResults)
  println("Timestep: $(print_results.timestep)")
  print_river_results(prognostic_fields)
  return prognostic_fields
end

function print_river_results(prognostic_fields::PrognosticFields)
  river_fields::RiverPrognosticFields = get_river_fields(prognostic_fields)
  println("")
  println("River Flow")
  println(river_fields.river_inflow)
  println("")
  println("Water to Ocean")
  println(river_fields.water_to_ocean)
end

function handle_event(prognostic_fields::PrognosticFields,
                      print_results::PrintSection)
   print_river_results_section(prognostic_fields)
  return prognostic_fields
end

function print_river_results_section(prognostic_fields::PrognosticFields)
  section_coords::LatLonSectionCoords = LatLonSectionCoords(65,75,125,165)
  river_fields::RiverPrognosticFields = get_river_fields(prognostic_fields)
  river_parameters::RiverParameters = get_river_parameters(prognostic_fields)
  for_section_with_line_breaks(river_parameters.grid,section_coords) do coords::Coords
    @printf("%.2f ",river_fields.river_inflow(coords))
    flush(stdout)
  end
  println("")
end

struct WriteRiverInitialValues <: Event end

function handle_event(prognostic_fields::PrognosticFields,
                      print_results::WriteRiverInitialValues)
  river_fields::RiverPrognosticFields = get_river_fields(prognostic_fields)
  river_parameters::RiverParameters = get_river_parameters(prognostic_fields)
  hd_start_filepath::AbstractString = "/Users/thomasriddick/Documents/data/temp/river_model_out.nc"
  write_river_initial_values(hd_start_filepath,river_parameters,river_fields)
  return prognostic_fields
end

function write_river_initial_values(hd_start_filepath::AbstractString,
                                    river_parameters::RiverParameters,
                                    river_prognostic_fields::RiverPrognosticFields)
  throw(UserError())
end

struct WriteRiverFlow <: Event
  timestep::Int64
end

function handle_event(prognostic_fields::PrognosticFields,
                      write_river_flow::WriteRiverFlow)
  river_fields::RiverPrognosticFields = get_river_fields(prognostic_fields)
  river_parameters::RiverParameters = get_river_parameters(prognostic_fields)
  write_river_flow_field(river_parameters,river_fields.river_inflow,
                         timestep=write_river_flow.timestep)
  return prognostic_fields
end

function write_river_flow_field(river_parameters::RiverParameters,
                                river_flow_field::Field{Float64};
                                timestep::Int64=-1)
  throw(UserError())
end

struct AccumulateRiverFlow <: Event end

function handle_event(prognostic_fields::PrognosticFields,
                      accumulate_river_flow::AccumulateRiverFlow)
  river_fields::RiverPrognosticFields = get_river_fields(prognostic_fields)
  river_diagnostic_output_fields::RiverDiagnosticOutputFields =
    get_river_diagnostic_output_fields(prognostic_fields)
  river_diagnostic_output_fields.cumulative_river_flow += river_fields.river_inflow
  return prognostic_fields
end

struct ResetCumulativeRiverFlow <: Event end

function handle_event(prognostic_fields::PrognosticFields,
                      reset_cumulative_river_flow::ResetCumulativeRiverFlow)
  river_diagnostic_output_fields::RiverDiagnosticOutputFields =
    get_river_diagnostic_output_fields(prognostic_fields)
  fill!(river_diagnostic_output_fields.cumulative_river_flow,0.0)
  return prognostic_fields
end

struct WriteMeanRiverFlow <: Event
  timestep::Int64
  number_of_timesteps::Int64
end

function handle_event(prognostic_fields::PrognosticFields,
                      write_mean_river_flow::WriteMeanRiverFlow)
  river_diagnostic_output_fields::RiverDiagnosticOutputFields =
    get_river_diagnostic_output_fields(prognostic_fields)
  river_parameters::RiverParameters = get_river_parameters(prognostic_fields)
  mean_river_flow::Field{Float64} =
      Field{Float64}(river_parameters.grid,0.0)
  for_all(river_parameters.grid) do coords::Coords
    set!(mean_river_flow,coords,
         river_diagnostic_output_fields.cumulative_river_flow(coords)/
         convert(Float64,write_mean_river_flow.number_of_timesteps))
  end
  write_river_flow_field(river_parameters,mean_river_flow,
                         timestep=write_mean_river_flow.timestep)
  return prognostic_fields
end

struct PrintGlobalValues <: Event
end

function handle_event(prognostic_fields::PrognosticFields,
                      print_global_values::PrintGlobalValues)
  river_fields::RiverPrognosticFields = get_river_fields(prognostic_fields)
  river_parameters::RiverParameters = get_river_parameters(prognostic_fields)
  global_sum_base_flow_res::Float64 = 0.0
  global_sum_overland_flow_res::Float64 = 0.0
  global_sum_river_flow_res::Float64 = 0.0
  println("Global sums:")
  for_all(river_parameters.grid) do coords::Coords
    global_sum_base_flow_res += river_fields.base_flow_reservoirs[1](coords)
    global_sum_overland_flow_res += river_fields.overland_flow_reservoirs[1](coords)
    for i = 1:5
      global_sum_river_flow_res += river_fields.river_flow_reservoirs[i](coords)
    end
  end
  println("Base Flow Res Content: $(global_sum_base_flow_res)")
  println("Overland Flow Res Content: $(global_sum_overland_flow_res)")
  println("River Flow Res Content: $(global_sum_river_flow_res)")
  return prognostic_fields
end

end

max_stars_repo_name: ThomasRiddick/DynamicHD

max_issues_repo_name: ThomasRiddick/DynamicHD

max_forks_repo_name: ThomasRiddick/DynamicHD


text:
module RE2

include("regex.jl")

end # module


text:
# -- PRIVATE METHODS ----------------------------------------------------------- #
function _process_gene_location(gene_location::String)::String

    # does the gene_location contain a .x?
    dot_index = findfirst(x -> x == '.', gene_location)
    if (isnothing(dot_index))
        return gene_location
    end

    # ok, we have a dot -
    bare_gene_location = gene_location[1:dot_index - 1]

    # return -
    return bare_gene_location
end
# ------------------------------------------------------------------------------ #

# -- PUBLIC METHODS ------------------------------------------------------------ #
function get_sequence_for_gene(gene_array::Array{String,1})::Some

    try

        # initialize -
        seq_obj_array = Array{KEGGSequence,1}()
        for gene in gene_array
            seq_obj = get_sequence_for_gene(gene) |> check
            if (isnothing(seq_obj) == false)
                push!(seq_obj_array, seq_obj)
            end
        end

        # return -
        return Some(seq_obj_array)
    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

function get_sequence_for_gene(gene_location::String)::Some

    # TODO: check gene string -

    try

        # gene seq - cutoff the trailing .x -
        kegg_gene_location = _process_gene_location(gene_location)
        if (isnothing(kegg_gene_location) == true)
            return Some(nothing)
        end

        # setup the url string -
        url_string = "$(_KEGG_GET_URL)/$(kegg_gene_location)/ntseq"

        # get the sequence -
        # returns: {header}\n{sequence_body}\n{?}
        ntseq = http_get_call_with_url(url_string) |> check
        if (length(ntseq) == 0)
            return Some(nothing)
        end

        # split on the newline -
        fragment_array = split(ntseq, '\n')

        # build new sequence type -
        sequence = KEGGSequence()
        sequence.gene_location = kegg_gene_location
        sequence.type = :gene
        sequence.header = fragment_array[1]

        # build seq string -
        seq_buffer = ""
        [seq_buffer *= string(x) for x in fragment_array[2:end - 1]]
        sequence.sequence = LongDNASeq(seq_buffer)

        # return the sequence -
        return Some(sequence)

    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

function get_sequence_for_protein(gene_array::Array{String,1})::Some

    try

        # initialize -
        seq_obj_array = Array{KEGGSequence,1}()
        for gene in gene_array
            seq_obj = get_sequence_for_protein(gene) |> check
            if (isnothing(seq_obj) == false)
                push!(seq_obj_array, seq_obj)
            end
        end

        # return -
        return Some(seq_obj_array)
    catch error

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)
    end
end

function get_sequence_for_protein(gene_location::String)::Some

    # TODO: check gene string -

    try

        # gene seq - cutoff the trailing .x -
        kegg_gene_location = _process_gene_location(gene_location)
        if (isnothing(kegg_gene_location) == true)
            return Some(nothing)
        end

        # setup the url string -
        url_string = "$(_KEGG_GET_URL)/$(kegg_gene_location)/aaseq"

        # get the sequence -
        # returns: {header}\n{sequence_body}\n{?}
        aaseq = http_get_call_with_url(url_string) |> check
        if (length(aaseq) == 0)
            return Some(nothing)
        end

        # split on the newline -
        fragment_array = split(aaseq, '\n')

        # build new sequence type -
        sequence = KEGGSequence()
        sequence.gene_location = kegg_gene_location
        sequence.type = :protein
        sequence.header = fragment_array[1]

        # build seq string -
        seq_buffer = ""
        [seq_buffer *= string(x) for x in fragment_array[2:end - 1]]
        sequence.sequence = LongSequence{AminoAcidAlphabet}(seq_buffer)

        # return the sequence -
        return Some(sequence)

    catch eror

        # get the original error message -
        error_message = sprint(showerror, error, catch_backtrace())
        vl_error_obj = ErrorException(error_message)

        # Package the error -
        return Some(vl_error_obj)

    end
end
# ------------------------------------------------------------------------------ #


text:
using STBImage
using Test

@testset "STBImage.jl" begin
    # Write your own tests here.
end


text:
using SubspacePreservingSparsification
using Test

@testset "SubspacePreservingSparsification.jl" begin
    @testset "Bin pattern" begin
        include("test_bin_pattern.jl")
    end
    @testset "diffenrent types without nullspace" begin
        include("test_compare_with_matlab.jl")
    end
    @testset "with nullspace" begin
        include("test_compare_with_matlab_nullspace.jl")
    end
    @testset "Integer Matrix" begin
        include("test_int_matrix.jl")
    end
end


text:
module Vinyl

using MacroTools

include("hooks.jl")
include("interpret.jl")
include("continuations.jl")

end # module


text:
module Higher

using ..Omega: RandVar, ciid, randtuple, cond, ==ᵣ, ==ₛ
using DocStringExtensions

using ..Space: transfertags

"A `RandVar` whose input inherits the tags the tags of ω"
inherittags(rv::RandVar, ω) = ciid(ω_ -> rv(transfertags(ω_, ω)))

include("rcd.jl")          # Random Conditional Distribution
include("rid.jl")          # Random Interventional Distribution

export rid, rcd, ∥, ∥ₛ

end


text:
using Test
using AugmentedGaussianProcesses
using LinearAlgebra
using Statistics
using Distributions
using KernelFunctions
const AGP = AugmentedGaussianProcesses
include("testingtools.jl")

nData = 100; nDim = 2
k = SqExponentialKernel(10.0)
ν = 5.0
r = 10
K = 4
X = rand(nData,nDim)
f = ones(nData)
while !(maximum(f) > 0 && minimum(f) < 0)
    global f = rand(MvNormal(zeros(nData),kernelmatrix(k,X,obsdim=1)+1e-3I))
end
width = maximum(f)-minimum(f)
normf = (f.-minimum(f))/width*K

y = Dict("Regression"=>f,"Classification"=>sign.(f),"MultiClass"=>floor.(Int64,normf),"Poisson"=>rand.(Poisson.(2.0*AGP.logistic.(f))),"NegBinomial"=>rand.(NegativeBinomial.(r,AGP.logistic.(f))))
n_class = length(unique(y["MultiClass"]))
reg_likelihood = ["GaussianLikelihood","StudentTLikelihood","LaplaceLikelihood","HeteroscedasticLikelihood"]
class_likelihood = ["BayesianSVM","LogisticLikelihood"]
multiclass_likelihood = ["LogisticSoftMaxLikelihood","SoftMaxLikelihood"]
poisson_likelihood = ["PoissonLikelihood"]
negbin_likelihood = ["NegBinomialLikelihood"]
likelihood_types = [reg_likelihood,class_likelihood,multiclass_likelihood,poisson_likelihood,negbin_likelihood]
# likelihood_types = [negbin_likelihood]
likelihood_names = ["Regression","Classification","MultiClass","Poisson","NegBinomial"]
# likelihood_names = ["NegBinomial"]
# inferences = ["GibbsSampling"]#,"NumericalInference"]
inferences = ["AnalyticVI"]#,"GibbsSampling","QuadratureVI"]
floattypes = [Float64]
@testset "VGP" begin
    for (likelihoods,l_names) in zip(likelihood_types,likelihood_names)
        @testset "$l_names" begin
            for l in likelihoods
                @testset "$(string(l))" begin
                    for inference in inferences
                        @testset "$(string(inference))" begin
                            if in(inference,methods_implemented_VGP[l])
                                for floattype in floattypes
                                    @test typeof(VGP(X,y[l_names],k,eval(Meta.parse(l*"("*addlargument(l)*")")),eval(Meta.parse(inference*"("*addiargument(false,inference)*")")))) <: VGP{floattype,eval(Meta.parse(l*"{"*string(floattype)*"}")),eval(Meta.parse(inference*"{"*string(floattype)*","*nlatent(l)*"}")),AGP._VGP{floattype},eval(Meta.parse(nlatent(l)))}
                                    global model = VGP(X,y[l_names],k,eval(Meta.parse(l*"("*addlargument(l)*")")),eval(Meta.parse(inference*"("*addiargument(false,inference)*")")),verbose=2)
                                    @test train!(model,50)
                                    @test testconv(model,l_names,X,y[l_names])
                                end
                            else
                                @test_throws AssertionError VGP(X,y[l_names],k,eval(Meta.parse(l*"("*addlargument(l)*")")),eval(Meta.parse(inference*"()")))
                            end
                        end
                    end
                end
            end
        end
    end
end


text:
const MAXMEM = 100

@testset "Iteration test" begin
  Δt = 1e-2
  prob1, _, _ = ConstrainedSystems.basic_constrained_problem()
  integrator = ConstrainedSystems.init(prob1, LiskaIFHERK(),dt=Δt)
  @test needs_iteration(integrator.f,integrator.u,integrator.p,integrator.u) == false

  prob2, _, _ = ConstrainedSystems.cartesian_pendulum_problem()
  integrator = ConstrainedSystems.init(prob2, LiskaIFHERK(),dt=Δt)
  @test needs_iteration(integrator.f,integrator.u,integrator.p,integrator.u) == true

end

@testset "Recursive copy" begin

  struct MyStruct{T}
    a :: T
  end

  n = 1000
  a = rand(n,n)
  z = similar(a)

  s1 = MyStruct(a)
  s2 = MyStruct(z)

  recursivecopy!(s2,s1)
  @test s2.a == s1.a
  @test !(s2.a === s1.a)
  @test @allocated(recursivecopy!(s2,s1)) < MAXMEM

  ss1 = MyStruct(s1)
  ss2 = MyStruct(MyStruct(z))
  recursivecopy!(ss2,ss1)
  @test ss2.a.a == ss1.a.a
  @test !(ss2.a.a === ss1.a.a)
  @test @allocated(recursivecopy!(s2,s1)) < MAXMEM

  t1 = [s1,s1]
  t2 = [MyStruct(z),MyStruct(z)]
  recursivecopy!(t2,t1)
  @test t2[1].a == t1[1].a
  @test t2[2].a == t1[2].a
  @test !(t2[1].a === t1[1].a)
  @test @allocated(recursivecopy!(s2,s1)) < MAXMEM

end