text:
# This file contains code that was formerly a part of Julia. License is MIT: http://julialang.org/license

using Base.Math: @horner, libm
using Base.MPFR: ROUNDING_MODE

for f in (:erf, :erfc)
    @eval begin
        ($f)(x::Float64) = ccall(($(string(f)),libm), Float64, (Float64,), x)
        ($f)(x::Float32) = ccall(($(string(f,"f")),libm), Float32, (Float32,), x)
        ($f)(x::Real) = ($f)(float(x))
        ($f)(a::Float16) = Float16($f(Float32(a)))
        ($f)(a::Complex{Float16}) = Complex{Float16}($f(Complex{Float32}(a)))
        function ($f)(x::BigFloat)
            z = BigFloat()
            ccall(($(string(:mpfr_,f)), :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32), z, x, ROUNDING_MODE[])
            return z
        end
        ($f)(x::AbstractFloat) = error("not implemented for ", typeof(x))
    end
end

for f in (:erf, :erfc, :erfcx, :erfi, :Dawson)
    fname = (f === :Dawson) ? :dawson : f
    @eval begin
        ($fname)(z::Complex{Float64}) = Complex{Float64}(ccall(($(string("Faddeeva_",f)),openspecfun), Complex{Float64}, (Complex{Float64}, Float64), z, zero(Float64)))
        ($fname)(z::Complex{Float32}) = Complex{Float32}(ccall(($(string("Faddeeva_",f)),openspecfun), Complex{Float64}, (Complex{Float64}, Float64), Complex{Float64}(z), Float64(eps(Float32))))

        ($fname)(z::Complex) = ($fname)(float(z))
        ($fname)(z::Complex{<:AbstractFloat}) = throw(MethodError($fname,(z,)))
    end
end

for f in (:erfcx, :erfi, :Dawson)
    fname = (f === :Dawson) ? :dawson : f
    @eval begin
        ($fname)(x::Float64) = ccall(($(string("Faddeeva_",f,"_re")),openspecfun), Float64, (Float64,), x)
        ($fname)(x::Float32) = Float32(ccall(($(string("Faddeeva_",f,"_re")),openspecfun), Float64, (Float64,), Float64(x)))

        ($fname)(x::Real) = ($fname)(float(x))
        ($fname)(x::AbstractFloat) = throw(MethodError($fname,(x,)))
    end
end

"""
    erf(x)
Compute the error function of `x`, defined by ``\\frac{2}{\\sqrt{\\pi}} \\int_0^x e^{-t^2} dt``
for arbitrary complex `x`.
"""
erf

"""
    erfc(x)
Compute the complementary error function of `x`, defined by ``1 - \\operatorname{erf}(x)``.
"""
erfc

"""
    erfcx(x)

Compute the scaled complementary error function of `x`, defined by ``e^{x^2} \\operatorname{erfc}(x)``.
Note also that ``\\operatorname{erfcx}(-ix)`` computes the Faddeeva function ``w(x)``.
"""
erfcx

"""
    erfi(x)

Compute the imaginary error function of `x`, defined by ``-i \\operatorname{erf}(ix)``.
"""
erfi

"""
    dawson(x)

Compute the Dawson function (scaled imaginary error function) of `x`, defined by
``\\frac{\\sqrt{\\pi}}{2} e^{-x^2} \\operatorname{erfi}(x)``.
"""
dawson

# Compute the inverse of the error function: erf(erfinv(x)) == x,
# using the rational approximants tabulated in:
#     J. M. Blair, C. A. Edwards, and J. H. Johnson, "Rational Chebyshev
#     approximations for the inverse of the error function," Math. Comp. 30,
#     pp. 827--830 (1976).
#         http://dx.doi.org/10.1090/S0025-5718-1976-0421040-7
#         http://www.jstor.org/stable/2005402
"""
    erfinv(x)

Compute the inverse error function of a real `x`, defined by ``\\operatorname{erf}(\\operatorname{erfinv}(x)) = x``.
"""
function erfinv(x::Float64)
    a = abs(x)
    if a >= 1.0
        if x == 1.0
            return Inf
        elseif x == -1.0
            return -Inf
        end
        throw(DomainError(a, "`abs(x)` cannot be greater than 1."))
    elseif a <= 0.75 # Table 17 in Blair et al.
        t = x*x - 0.5625
        return x * @horner(t, 0.16030_49558_44066_229311e2,
                             -0.90784_95926_29603_26650e2,
                              0.18644_91486_16209_87391e3,
                             -0.16900_14273_46423_82420e3,
                              0.65454_66284_79448_7048e2,
                             -0.86421_30115_87247_794e1,
                              0.17605_87821_39059_0) /
                   @horner(t, 0.14780_64707_15138_316110e2,
                             -0.91374_16702_42603_13936e2,
                              0.21015_79048_62053_17714e3,
                             -0.22210_25412_18551_32366e3,
                              0.10760_45391_60551_23830e3,
                             -0.20601_07303_28265_443e2,
                              0.1e1)
    elseif a <= 0.9375 # Table 37 in Blair et al.
        t = x*x - 0.87890625
        return x * @horner(t, -0.15238_92634_40726_128e-1,
                               0.34445_56924_13612_5216,
                              -0.29344_39867_25424_78687e1,
                               0.11763_50570_52178_27302e2,
                              -0.22655_29282_31011_04193e2,
                               0.19121_33439_65803_30163e2,
                              -0.54789_27619_59831_8769e1,
                               0.23751_66890_24448) /
                   @horner(t, -0.10846_51696_02059_954e-1,
                               0.26106_28885_84307_8511,
                              -0.24068_31810_43937_57995e1,
                               0.10695_12997_33870_14469e2,
                              -0.23716_71552_15965_81025e2,
                               0.24640_15894_39172_84883e2,
                              -0.10014_37634_97830_70835e2,
                               0.1e1)
    else # Table 57 in Blair et al.
        t = 1.0 / sqrt(-log(1.0 - a))
        return @horner(t, 0.10501_31152_37334_38116e-3,
                          0.10532_61131_42333_38164_25e-1,
                          0.26987_80273_62432_83544_516,
                          0.23268_69578_89196_90806_414e1,
                          0.71678_54794_91079_96810_001e1,
                          0.85475_61182_21678_27825_185e1,
                          0.68738_08807_35438_39802_913e1,
                          0.36270_02483_09587_08930_02e1,
                          0.88606_27392_96515_46814_9) /
              (copysign(t, x) *
               @horner(t, 0.10501_26668_70303_37690e-3,
                          0.10532_86230_09333_27531_11e-1,
                          0.27019_86237_37515_54845_553,
                          0.23501_43639_79702_53259_123e1,
                          0.76078_02878_58012_77064_351e1,
                          0.11181_58610_40569_07827_3451e2,
                          0.11948_78791_84353_96667_8438e2,
                          0.81922_40974_72699_07893_913e1,
                          0.40993_87907_63680_15361_45e1,
                          0.1e1))
    end
end

function erfinv(x::Float32)
    a = abs(x)
    if a >= 1.0f0
        if x == 1.0f0
            return Inf32
        elseif x == -1.0f0
            return -Inf32
        end
        throw(DomainError(a, "`abs(x)` cannot be greater than 1."))
    elseif a <= 0.75f0 # Table 10 in Blair et al.
        t = x*x - 0.5625f0
        return x * @horner(t, -0.13095_99674_22f2,
                               0.26785_22576_0f2,
                              -0.92890_57365f1) /
                   @horner(t, -0.12074_94262_97f2,
                               0.30960_61452_9f2,
                              -0.17149_97799_1f2,
                               0.1f1)
    elseif a <= 0.9375f0 # Table 29 in Blair et al.
        t = x*x - 0.87890625f0
        return x * @horner(t, -0.12402_56522_1f0,
                               0.10688_05957_4f1,
                              -0.19594_55607_8f1,
                               0.42305_81357f0) /
                   @horner(t, -0.88276_97997f-1,
                               0.89007_43359f0,
                              -0.21757_03119_6f1,
                               0.1f1)
    else # Table 50 in Blair et al.
        t = 1.0f0 / sqrt(-log(1.0f0 - a))
        return @horner(t, 0.15504_70003_116f0,
                          0.13827_19649_631f1,
                          0.69096_93488_87f0,
                         -0.11280_81391_617f1,
                          0.68054_42468_25f0,
                         -0.16444_15679_1f0) /
              (copysign(t, x) *
               @horner(t, 0.15502_48498_22f0,
                          0.13852_28141_995f1,
                          0.1f1))
    end
end

erfinv(x::Integer) = erfinv(float(x))

# Inverse complementary error function: use Blair tables for y = 1-x,
# exploiting the greater accuracy of y (vs. x) when y is small.
"""
    erfcinv(x)

Compute the inverse error complementary function of a real `x`, defined by
``\\operatorname{erfc}(\\operatorname{erfcinv}(x)) = x``.
"""
function erfcinv(y::Float64)
    if y > 0.0625
        return erfinv(1.0 - y)
    elseif y <= 0.0
        if y == 0.0
            return Inf
        end
        throw(DomainError(y, "`y` must be nonnegative."))
    elseif y >= 1e-100 # Table 57 in Blair et al.
        t = 1.0 / sqrt(-log(y))
        return @horner(t, 0.10501_31152_37334_38116e-3,
                          0.10532_61131_42333_38164_25e-1,
                          0.26987_80273_62432_83544_516,
                          0.23268_69578_89196_90806_414e1,
                          0.71678_54794_91079_96810_001e1,
                          0.85475_61182_21678_27825_185e1,
                          0.68738_08807_35438_39802_913e1,
                          0.36270_02483_09587_08930_02e1,
                          0.88606_27392_96515_46814_9) /
              (t *
               @horner(t, 0.10501_26668_70303_37690e-3,
                          0.10532_86230_09333_27531_11e-1,
                          0.27019_86237_37515_54845_553,
                          0.23501_43639_79702_53259_123e1,
                          0.76078_02878_58012_77064_351e1,
                          0.11181_58610_40569_07827_3451e2,
                          0.11948_78791_84353_96667_8438e2,
                          0.81922_40974_72699_07893_913e1,
                          0.40993_87907_63680_15361_45e1,
                          0.1e1))
    else # Table 80 in Blair et al.
        t = 1.0 / sqrt(-log(y))
        return @horner(t, 0.34654_29858_80863_50177e-9,
                          0.25084_67920_24075_70520_55e-6,
                          0.47378_13196_37286_02986_534e-4,
                          0.31312_60375_97786_96408_3388e-2,
                          0.77948_76454_41435_36994_854e-1,
                          0.70045_68123_35816_43868_271e0,
                          0.18710_42034_21679_31668_683e1,
                          0.71452_54774_31351_45428_3e0) /
          (t * @horner(t, 0.34654_29567_31595_11156e-9,
                          0.25084_69079_75880_27114_87e-6,
                          0.47379_53129_59749_13536_339e-4,
                          0.31320_63536_46177_68848_0813e-2,
                          0.78073_48906_27648_97214_733e-1,
                          0.70715_04479_95337_58619_993e0,
                          0.19998_51543_49112_15105_214e1,
                          0.15072_90269_27316_80008_56e1,
                          0.1e1))
    end
end

function erfcinv(y::Float32)
    if y > 0.0625f0
        return erfinv(1.0f0 - y)
    elseif y <= 0.0f0
        if y == 0.0f0
            return Inf32
        end
        throw(DomainError(y, "`y` must be nonnegative."))
    else # Table 50 in Blair et al.
        t = 1.0f0 / sqrt(-log(y))
        return @horner(t, 0.15504_70003_116f0,
                          0.13827_19649_631f1,
                          0.69096_93488_87f0,
                         -0.11280_81391_617f1,
                          0.68054_42468_25f0,
                         -0.16444_15679_1f0) /
        (t * @horner(t, 0.15502_48498_22f0,
                        0.13852_28141_995f1,
                        0.1f1))
    end
end

erfcinv(x::Integer) = erfcinv(float(x))

# MPFR has an open TODO item for this function
# until then, we use [DLMF 7.12.1](https://dlmf.nist.gov/7.12.1) for the tail
function erfcx(x::BigFloat)
    if x <= (Clong == Int32 ? 0x1p15 : 0x1p30)
        # any larger gives internal overflow
        return exp(x^2)*erfc(x)
    elseif !isfinite(x)
        return 1/x
    else
        # asymptotic series
        # starts to diverge at iteration i = 2^30 or 2^60
        # final term will be < Γ(2*i+1)/(2^i * Γ(i+1)) / (2^(i+1))
        # so good to (lgamma(2*i+1) - lgamma(i+1))/log(2) - 2*i - 1
        #            ≈ 3.07e10 or 6.75e19 bits
        # which is larger than the memory of the respective machines
        ϵ = eps(BigFloat)/4
        v = 1/(2*x*x)
        k = 1
        s = w = -k*v
        while abs(w) > ϵ
            k += 2
            w *= -k*v
            s += w
        end
        return (1+s)/(x*sqrt(oftype(x,pi)))
    end
end


text:
function fib(n::Int)
  if n < 2
    return n
  end
  return fib(n - 2) + fib(n - 1)
end

println(fib(35))

ext: jl

max_stars_repo_name: peridot/peridot

max_issues_repo_name: peridot/peridot

max_forks_repo_name: peridot/peridot


text:
#!/cvmfs/soft.computecanada.ca/easybuild/software/2020/avx2/Core/julia/1.5.2/bin/julia
#SBATCH --mail-user=mkschleg@ualberta.ca
#SBATCH --mail-type=ALL
#SBATCH -o reproduce_test.out # Standard output
#SBATCH -e reproduce_test.err # Standard error
#SBATCH --mem-per-cpu=2000M # Memory request of 2 GB
#SBATCH --time=0:05:00 # Running time of 12 hours
#SBATCH --ntasks=4
#SBATCH --account=rrg-whitem

using Pkg
Pkg.activate(".")

using Reproduce
const config_file = "examples/arg_iter_config.toml"

function test_experiment()

    experiment = Experiment(config_file)

    pre_experiment(experiment; tldr="hello tldr")
    ret = job(experiment; num_workers=6)
    post_experiment(experiment, ret)

end

test_experiment()


text:
function parserule(rule)
    ss = split(rule, " bags contain ")
    rule = ss[1]
    edges = Vector{typeof(rule)}()
    rest = split(ss[2], ", ")
    rx = r"([0-9]+) ([a-z]+ [a-z]+)"
    for s in rest
        m = match(rx, s)
        if m == nothing
            continue
        end
        num = parse(Int, m.captures[1])
        append!(edges, repeat([m.captures[2]], num))
    end
    (rule, edges)
end

function parserules(rules)
    ss = split(rules, "\n", keepempty=false)
    map(parserule, ss)
end

#=
parserule("light red bags contain 1 bright white bag, 2 muted yellow bags.")
parserule("bright white bags contain 1 shiny gold bag.")
parserule("faded blue bags contain no other bags.")

rules = """
light red bags contain 1 bright white bag, 2 muted yellow bags.
dark orange bags contain 3 bright white bags, 4 muted yellow bags.
bright white bags contain 1 shiny gold bag.
muted yellow bags contain 2 shiny gold bags, 9 faded blue bags.
shiny gold bags contain 1 dark olive bag, 2 vibrant plum bags.
dark olive bags contain 3 faded blue bags, 4 dotted black bags.
vibrant plum bags contain 5 faded blue bags, 6 dotted black bags.
faded blue bags contain no other bags.
dotted black bags contain no other bags.
"""
parsed_rules = parserules(rules)
=#

function build_graph(parsed_rules::Vector{Tuple{SubString{String}, Vector{SubString{String}}}})
    m = length(parsed_rules)
    # Rows indicate what bags 'this' bag may contain. Columns indicate what bags may
    # contain 'this' bag.
    M = zeros(Int64, (m, m))
    idx_tracker = Dict{AbstractString, Int64}()
    # Just build index tracker first to make updating matrix easy.
    for (i, (bag, _)) in enumerate(parsed_rules)
        idx_tracker[bag] = i
    end
    # Build up matrix.
    for (i, (bag, contains_bags)) in enumerate(parsed_rules)
        for contains_bag in contains_bags
            j = idx_tracker[contains_bag]
            M[i, j] += 1
        end
    end
    (idx_tracker, M)
end

#=
(idxs, M) = build_graph(parsed_rules)
idxs
M
=#

function dfs_walk(f, M, col_idx)
    for (i, val) in enumerate(M[:, col_idx])
        if val == 0
            continue
        end
        f(i, col_idx, val)
        dfs_walk(f, M, i)
    end
end

function coalesce_bag_vals(M, row_idx)
    acc = sum(M[row_idx, :])
    if acc == 0
        return 0
    end
    for (j, val) in enumerate(M[row_idx, :])
        if val == 0
            continue
        end
        acc += val * coalesce_bag_vals(M, j)
    end
    acc
end

#=
dfs_walk((i, j, val)->println("walked $(i) $(j) $(val)"), M, 5)
coalesce_bag_vals(M, 5)
=#

function num_bags_contain(M, idxs, bagstr)
    if !haskey(idxs, bagstr)
        error("No index for bag")
    end
    idx = get(idxs, bagstr, 0)
    walked_bags = Vector{Int64}(undef, 0)
    dfs_walk((i, _, _) -> push!(walked_bags, i), M, idx)
    length(unique(walked_bags))
end

#=
num_bags_contain(M, idxs, "shiny gold")
bag_vals = num_bags_contain_val(M, idxs, "shiny gold")
=#

function day07part1(s)
    parsed_rules = parserules(s)
    (idxs, M) = build_graph(parsed_rules)
    num_bags_contain(M, idxs, "shiny gold")
end

function day07part2(s)
    parsed_rules = parserules(s)
    (idxs, M) = build_graph(parsed_rules)
    coalesce_bag_vals(M, get(idxs, "shiny gold", 0))
end

function day07(file)
    s = read(file, String)
    (day07part1(s), day07part2(s))
end

#=
day07("input.txt")
=#


text:
#
# implementation of Spline-based pair potentials
# and other functionality
#

import Dierckx
using Dierckx: Spline1D
using ForwardDiff

# https://github.com/kbarbary/Dierckx.jl

"""
`type SplinePairPotential`

Pair potential defined via splines. In the construction of these splines
atomic units (Angstrom) are assumed. If this is used in different units,
then something will probably break.

### Constructors:
```
SplinePairPotential(xdat, ydat; kwargs...)   # fit spline  from data points
SplinePairPotential(fname; kwargs...)        # load data points from file
```
Keyword arguments:
* `s = 0`: balances smoothness versus fit; by default this matches the fits obtained from ASE and LAMMPS
* `fixcutoff = true`: if true, then the fit will be artificially modified at the cut-off to ensure that the transition to zero is smooth
* `order = 3`: can use lower or higher order splines (0 <= order <= 5) but only 3 is tested
* `w = (1.0 + ydat).^(-2)`: this gives relative weights to datapoints to ensure a good in the important regions; the intuition behind the default choice is that is prevents overfitting at very high energies which are not physical anyhow, but ensure that sufficient data points are used in the low energy region.
"""
mutable struct SplinePairPotential <: SimplePairPotential
   spl::Spline1D          # The actual spline object
   rcut::Float64          # cutoff radius (??? could just use spl.t[end] ???)
   wrk::Vector{Float64}   # a work array for faster evaluation of derivatives
end

@pot SplinePairPotential

SplinePairPotential(spl::Spline1D) =
   SplinePairPotential(spl, maximum(spl.t), Vector{Float64}(undef, length(spl.t)))

cutoff(V::SplinePairPotential) = V.rcut

evaluate(V::SplinePairPotential, r::Number) = _evalspl(V.spl, r)
evaluate_d(V::SplinePairPotential, r::Number) = _evalspl_d(V.spl, r)
evaluate_dd(V::SplinePairPotential, r::Number) = _evalspl_dd(V.spl, r)

# a few interface routines so we can define AD in a nice way.
_evalspl(s::Spline1D, r) = s(r)
_evalspl_d(s::Spline1D, r) = Dierckx.derivative(s, r, 1)
_evalspl_dd(s::Spline1D, r) = Dierckx.derivative(s, r, 2)

function _evalspl(s::Spline1D, d::ForwardDiff.Dual{T}) where T
   x = ForwardDiff.value(d)
   ForwardDiff.Dual{T}( _evalspl(s, x),
                        _evalspl_d(s, x) * ForwardDiff.partials(d) )
end

function _evalspl_d(s::Spline1D, d::ForwardDiff.Dual{T}) where T
   x = ForwardDiff.value(d)
   ForwardDiff.Dual{T}( _evalspl_d(s, x),
                        _evalspl_dd(s, x) * ForwardDiff.partials(d) )
end

function SplinePairPotential(xdat, ydat; s = 0, fixcutoff=true, order=3,
                             w = (1.0 .+ abs.(ydat)).^(-2))
   # this creates a "fit" with s determining the balance between smoothness
   # and fitting the data (basically an error bound)
   spl = Spline1D(xdat, ydat; bc = "zero", s = s, k = order, w = w)
   # set the last few spline data-points to 0 to get a guaranteed
   # smooth transition to 0 at the cutoff.
   if fixcutoff
      spl.c[end-order+1:end] .= 0.0
   end
   return SplinePairPotential(spl)
end

SplinePairPotential(fname::AbstractString; kwargs...) =
   SplinePairPotential(load_ppfile(fname)...; kwargs...)

function load_ppfile(fname::AbstractString)
   if fname[end-3:end] == ".plt"
      return load_plt(fname)
   else
      error("unknown file type: $(fname[end-3:end])")
   end
   nothing
end

function load_plt(fname::AbstractString)
   data = readdlm(fname, Float64; comments=true)
   return data[:, 1], data[:, 2]
end


text:
"""

Find centers of the intervals defining locations of the noisy
observations recursively.

l       : left boundary
r       : right boundary
min_len : max length of the interval

"""
function _recmidpoints(l, r, min_len)
    if r - l <= min_len
        [r]
    else
        m = div(r-l, 2)
        vcat(_recmidpoints(l, l + m, min_len), _recmidpoints(l + m, r, min_len))
    end
end

function recmidpoints(l, r, min_len)
    mids = _recmidpoints(l, r, min_len)
    mids[1:end-1]
end

if 1 == 0
    ret = recmidpoints(1, 200, 50)
    println(ret)
end

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
mutable struct WPFBelief{S, A, O} <: AbstractParticleBelief{S}
    particles::Vector{S}
    weights::Vector{Float64}
    weight_sum::Float64
    belief::Int
    depth::Int
    tree::AdaOPSTree{S,A,O}
    _obs::O
    _probs::Union{Nothing, Dict{S,Float64}}
    _hist::Union{Nothing, Array{NamedTuple,1}}
end

WPFBelief(particles::Vector{S}, weights::Vector{Float64}, weight_sum::N, belief::Int, depth::Int, tree::AdaOPSTree{S,A,O}, obs::O) where {S,A,O,N<:Real} = WPFBelief{S,A,O}(particles, weights, convert(Float64, weight_sum), belief, depth, tree, obs, nothing, nothing)
WPFBelief(particles, weights, belief, depth, tree, obs) = WPFBelief(particles, weights, sum(weights), belief, depth, tree, obs, nothing, nothing)

function ParticleFilters.rand(rng::RNG, b::WPFBelief{S}) where {S,RNG<:AbstractRNG}
    t = rand(rng) * b.weight_sum
    i = 1
    cum_weight = b.weights[1]
    while cum_weight < t
        i += 1
        cum_weight += b.weights[i]
    end
    return b.particles[i]
end

ParticleFilters.particles(b::WPFBelief{S}) where S = b.particles
ParticleFilters.n_particles(b::WPFBelief) = length(b.particles)
ParticleFilters.weight(b::WPFBelief, i::Int) = b.weights[i]
ParticleFilters.particle(b::WPFBelief{S}, i::Int) where S = b.particles[i]
ParticleFilters.weight_sum(b::WPFBelief) = b.weight_sum
ParticleFilters.weights(b::WPFBelief) = b.weights
ParticleFilters.weighted_particles(b::WPFBelief{S}) where S = (b.particles[i]=>b.weights[i] for i in 1:length(b.weights))

function POMDPs.mean(b::WPFBelief{S}) where S
    mean_s = zero(eltype(particles(b)))
    for (w, s) in weighted_particles(b)
        mean_s += w * s
    end
    return mean_s / weight_sum(b)
end
POMDPs.currentobs(b::WPFBelief{S,A,O}) where {S,A,O} = b._obs

function POMDPs.history(belief::WPFBelief{S,A,O}) where {S,A,O}
    if belief._hist === nothing
        belief._hist = Vector{NamedTuple}(undef, belief.depth+1)
        D = belief.tree
        if belief.belief !== 1
            b = belief.belief
            depth = belief.depth + 1
            while depth !== 1
                ba = D.parent[b]
                belief._hist[depth] = (o=D.obs[b], a=D.ba_action[ba])
                b = D.ba_parent[ba]
                depth -= 1
            end
        end
        belief._hist[1] = (o=D.obs[1],)
    end
    return belief._hist
end

function resample!(resampled::WeightedParticleBelief, b::AbstractParticleBelief, rng::AbstractRNG)
    m = n_particles(resampled)
    step = weight_sum(b)/m
    U = rand(rng)*step
    c = weight(b,1) # accumulate sum of weights
    i = 1
    P = particles(b)
    P_resampled = particles(resampled)
    @inbounds for j in 1:m
        while U > c
            i += 1
            c += weight(b, i)
        end
        U += step
        P_resampled[j] = P[i]
    end
    return resampled
end

function resample!(resampled::WeightedParticleBelief, b, pomdp::POMDP, rng::AbstractRNG)
    m = n_particles(resampled)
    i = 0
    P_resampled = particles(resampled)
    @inbounds for i in 1:m
        s = rand(rng, b)
        while isterminal(pomdp, s)
            s = rand(rng, b)
        end
        P_resampled[i] = s
    end
    return resampled
end

function switch_to_sibling!(b::WPFBelief{S,A,O}, obs::O, weights::Array{Float64,1}) where {S,A,O}
    b.weights = weights
    b.weight_sum = sum(weights)
    b._obs = obs
    b._probs = nothing
    if b._hist !== nothing
        if length(b._hist) > 1
            b._hist[end] = (a=b._hist[end].a, o=obs)
        else
            b._hist[end] = (o=obs,)
        end
    end
end

function ParticleFilters.probdict(b::WPFBelief{S,A,O}) where {S,A,O}
    if b._probs === nothing
        # update the cache
        probs = Dict{S, Float64}()
        for (i,p) in enumerate(particles(b))
            if haskey(probs, p)
                probs[p] += weight(b, i)/weight_sum(b)
            else
                probs[p] = weight(b, i)/weight_sum(b)
            end
        end
        b._probs = probs
    end
    return b._probs
end

# Extract Specific type of belief from WPFBelief
POMDPs.initialize_belief(::NothingUpdater, b::WPFBelief) = nothing
POMDPs.initialize_belief(::PreviousObservationUpdater, b::WPFBelief{S,A,O}) where {S,A,O} = b._obs

function POMDPs.initialize_belief(up::KMarkovUpdater, b::WPFBelief{S,A,O}) where {S,A,O}
    hist = history(b)
    if length(hist) > up.k
        [tuple[:o] for tuple in hist[end-up.k+1:end]]
    else
        [tuple[:o] for tuple in hist]
    end
end

function POMDPs.initialize_belief(up::BasicParticleFilter, b::WPFBelief{S}) where S
    ParticleFilters.resample(up.resampler, b, up.rng)
end


text:
function factorial(n::Int64)
    if n == 0 || n == 1
        return 1
    else
        return n * factorial(n-1)
    end
end
n = parse(Int64, readline())
println(factorial(n))

max_stars_repo_path: Hackerrank/30 Days of Code/Julia/day 09.jl

max_stars_repo_name: Next-Gen-UI/Code-Dynamics

max_issues_repo_name: Next-Gen-UI/Code-Dynamics

max_forks_repo_name: Next-Gen-UI/Code-Dynamics


text:
#
# ropsargs
#

function reverse_pass!(∂ls::∂LoopSet)
    N = length(∂ls.opsparentsfirst)
    # @show ∂ls.opsparentsfirst
    # @show 5, ∂ls.ropsargs
    for i ∈ ∂ls.lsold.outer_reductions
        if ∂ls.tracked_ops[i]
            add_reverse_outer_reduct!(∂ls, i)
        end
    end
    # @show 11, ∂ls.ropsargs ∂ls.opsparentsfirst
    for n ∈ 0:N-1
        i = ∂ls.opsparentsfirst[N - n]
        ∂ls.tracked_ops[i] && add_reverse_operation!(∂ls, i)
    end
end

function add_reverse_outer_reduct!(∂ls::∂LoopSet, i::Int)
    @unpack rls, lsold, ropsargs = ∂ls
    oldop = operations(lsold)[i]
    constop = LoopVectorization.add_constant!(rls, diffsym(name(oldop)), 8)
    push!(ropsargs[i], constop)
    nothing
end

function add_reverse_operation!(∂ls::∂LoopSet, i::Int)
    @unpack fls, lsold = ∂ls
    ops = operations(lsold)
    op = ops[i]
    # @show op
    # @show 31, ∂ls.ropsargs, op
    if isload(op)
        add_load_reverse!(∂ls, i)
        # if istracked(∂ls, op)
        #     adjop = adj(name(op))
        #     ∂ref = ArrayReferenceMeta( ArrayReference( adjop, getindices(op) ), op.ref.loopedindex )
        #     ∂op = Operation(
        #         - 1, adjop, 8, :setindex!, memstore, loopdeps, NODEPENDENCY, Operation[], ∂ref, NODEPENDENCY
        #     )
        #     # backlogged store; wait until it has a parent to add
        #     ∂newops[i] = ∂op
        # end
    elseif iscompute(op) && !isreductcombineinstr(op)
        add_compute_reverse!(∂ls, i)
    elseif isstore(op)
        add_store_reverse!(∂ls, i)
        # if istracked(∂ls, op)
        #     adjop = adj(name(op))
        #     ∂ref = ArrayReferenceMeta( ArrayReference( adjop, getindices(op) ), op.ref.loopedindex )
        #     ∂op = Operation(
        #         - 1, adjop, 8, :getindex, memload, loopdeps, NODEPENDENCY, Operation[], ∂ref, NODEPENDENCY
        #     )
        #     # backlogged store; wait until it has a parent to add
        #     ∂newops[i] = ∂op
        # end
    # else#if isconstant(op)
        # add_constant_reverse!(∂ls, i)
    end
    nothing
end

function adjref(ref::ArrayReferenceMeta)
    LoopVectorization.ArrayReferenceMeta(LoopVectorization.ArrayReference(diffsym(ref.ref.array), ref.ref.indices), ref.loopedindex)
end
adjref(op::Operation) = adjref(op.ref)

isreductcombineinstr(op::Operation) = op.instruction.instr === :identity
function add_reverse_input!(ropsargs::Vector{Vector{Operation}}, opold::Operation, opnew::Operation)
    if isreductcombineinstr(opold) # changed how this is handled...
        add_reverse_input_parent!(ropsargs, opold, opnew, 1)
    else
        push!(ropsargs[identifier(opold)], opnew)
    end
    nothing
end
# Inserts the parent in ropsargs vector so that it may be accessed
function add_reverse_input_parent!(ropsargs::Vector{Vector{Operation}}, opold::Operation, opnew::Operation, i::Int)
    add_reverse_input!(ropsargs, parents(opold)[i], opnew)
end

# function add_constant_reverse!(∂ls::∂LoopSet, i::Int) end

# Store an adjoint
function add_load_reverse!(∂ls::∂LoopSet, i::Int)
    # istracked(∂ls, i) || return
    @unpack rls, lsold, ropsargs = ∂ls
    opold = operations(lsold)[i]
    rops = operations(rls)
    ref = adjref(opold)
    vparents = [get_adj_input_op!(∂ls, i)]
    opnew = Operation(
        length(rops), diffsym(name(opold)), 8, :setindex!, memstore, loopdependencies(opold), NODEPENDENCY, vparents, ref
    )
    add_store!(rls, opnew)
    nothing
end
# Reverse stores have no parents.
function add_store_reverse!(∂ls::∂LoopSet, i::Int)
    # istracked(∂ls, i) || return
    @unpack rls, lsold, ropsargs = ∂ls
    opold = operations(lsold)[i]
    rops = operations(rls)
    ref = adjref(opold)
    opnew = Operation(
        length(rops), diffsym(name(opold)), 8, :getindex, memload, loopdependencies(opold), NODEPENDENCY, NOPARENTS, ref, NODEPENDENCY
    )
    # push!(rops, opnew)
    add_load!(rls, opnew)
    add_reverse_input_parent!(ropsargs, opold, opnew, 1)
    nothing
end

function issubset(sub, sup)
    length(sub) < length(sup)
end

function make_updating_failed!(op₁::Operation, op₂::Operation)
    iscompute(op₁) || return true
    instr = instruction(op₁)
    if instr === :reduce_to_add
        op₁.instruction = :reduced_add
        parents(op₁)[2] = op₂
        return false
    end
    updateinstr = get(ReverseDiffExpressionsBase.MAKEUPDATING, instr, nothing)
    updateinstr === nothing && return true
    op₁.instruction = updateinstr
    pushfirst!(parents(op₁), op₂)
    false
end
function make_updating!(ls::LoopSet, op₁::Operation, op₂::Operation, loopdeps = loopdependencies(op₂))
    if make_updating_failed!(op₁, op₂)
        op₃ = Operation(
            length(operations(ls)), gensym(:combineadjoints), 8, :vadd, compute, loopdeps, Symbol[], [op₁, op₂], NOTAREFERENCE
        )
        push!(operations(ls), op₃)
        op₃
    else
        op₁
    end
end

function get_adj_input_op!(∂ls::∂LoopSet, ropsargsᵢ::Vector{Operation}, i::Int)
    @unpack rls, lsold = ∂ls
    oldop = operations(lsold)[i]
    rop₁ = first(ropsargsᵢ) # Throws boundserror if length 0
    for j ∈ 2:length(ropsargsᵢ)
        ropⱼ = ropsargsᵢ[j]
        if make_updating_failed!(rop₁, ropⱼ)
            rop₁ = make_updating!(rls, ropⱼ, rop₁, loopdependencies(oldop))
        end
    end
    rop₁
end

# function LoopVectorization.matches(v1, v2)
# function matches(v1, v2)
#     length(v1) == length(v2) || return false
#     if length(v1) == 0
#         return true
#     elseif length(v1) == 1
#         first(v1) == first(v2) && return true
#     elseif length(v1) == 2
#         v11 = v1[1]
#         v12 = v1[2]
#         v21 = v2[1]
#         v22 = v2[2]
#         if v11 == v21 && v12 == v22
#             return true
#         elseif v11 == v22 && v12 == v21
#             return true
#         else
#             return false
#         end
#     elseif v1 == v2
#         return true
#     else# trying to avoid sort
#         return sort(v1) == sort(v2)
#     end
# end

# function reductzero!(∂ls::∂LoopSet, op, reducedc)
#     instr = instruction(op)
#     # If there are different numbers of arguments...
#     vparents = parents(op)
#     length(vparents) == 2 || return :zero
#     parent1 = vparents[1]
#     parent2 = vparents[2]
#     mulops = :(:vmul, :*, :evmul)
#     if isconstant(parent1) && matches(reducedchildren(parent1), reducedc)

#     elseif isconstant(parent2) && matches(reducedchildren(parent2), reducedc)

#     end
#     :zero
# end

# This method means that we have a reduction.
# We must determine which variables are the reductions, and update parents accordingly
# If a parent performing the reduction is in the `c = a * b` family, we should swap `∂a` defintion of `b` with `c / a`.
# For reduction, ops must also be made updating.
#
# Strategy of this function is to go through ropsargs, and create reduced/reduce_to ops for all reductions,
# and then to call get_adj_input_op_noreduct!(∂ls::∂LoopSet, ropsargsᵢ, i::Int)
function get_adj_input_op_reduct!(∂ls::∂LoopSet, ropsargsᵢ::Vector{Operation}, i::Int)
    @unpack rls, diffops, lsold = ∂ls
    oldop = operations(lsold)[i]
    loopdeps = loopdependencies(oldop)
    # We need a simple reduce_to target, shareable by all reductions
    # mCt = gensym(:targetzero)
    # targetzero = LoopVectorization.add_constant!(rls, gensym(:targetzero), loopdeps, mCt, 8, :numericconstant)
    # push!(rls.preamble_zeros, (identifier(targetzero), LoopVectorization.IntOrFloat))
    isouterreduction = iszero(length(loopdeps))
    # And a reduction-specific zero init
    for j ∈ eachindex(ropsargsᵢ)
        ropⱼ = ropsargsᵢ[j]
        ropⱼdeps = loopdependencies(ropⱼ)
        length(ropⱼdeps) == length(loopdeps) && continue
        reducedc = [s for s ∈ ropⱼdeps if s ∉ loopdeps]
        # Now, we must
        # [x] 1. make ropⱼ a reduction, adding across the reducedc loops
        # [x] 2. Create a reduction zero initializer
        # [x] 3. Create a reduction finalizer (i.e., reduced_/reduce_to)
        append!(reduceddependencies(ropⱼ), reducedc)
        reductinit = LoopVectorization.add_constant!(rls, gensym(:reductzero), loopdeps, name(ropⱼ), 8, :numericconstant)
        push!(rls.preamble_zeros, (identifier(reductinit), LoopVectorization.IntOrFloat))
        reductinit.reduced_children = reducedc
        opcomb = make_updating!(rls, ropⱼ, reductinit)
        if isouterreduction
            push!(rls.outer_reductions, identifier(opcomb))
            ropsargsᵢ[j] = opcomb
        else
            # ropⱼr = Operation(length(operations(rls)), gensym(:reduct), 8, :reduce_to_add, compute, loopdeps, reducedc, [opcomb, targetzero])
            ropⱼr = Operation(length(operations(rls)), gensym(:reduct), 8, :identity, compute, loopdeps, reducedc, [opcomb])
            push!(operations(rls), ropⱼr)
            ropsargsᵢ[j] = ropⱼr
        end
    end
    # We've reduced all the reduction variables, so no we use the noreduct version.
    get_adj_input_op!(∂ls, ropsargsᵢ, i)
end

function get_adj_input_op_expand!(∂ls::∂LoopSet, ropsargsᵢ::Vector{Operation}, i::Int)
    @unpack rls, diffops, lsold = ∂ls
    oldop = operations(lsold)[i]
    loopdeps = loopdependencies(oldop)
    reduceddeps = reduceddependencies(oldop)
    for j ∈ eachindex(ropsargsᵢ)
        ropⱼ = ropsargsᵢ[j]
        instrⱼ = instruction(ropⱼ)
        if instrⱼ.instr ∉ (:*, :vmul, :evmul)
            continue
        end
        # I don't think I have to copy oldop? Why not just reuse it?
        push!(operations(rls), oldop)
        op = Operation(length(operations(rls)), gensym(:prodgrad), 8, :/, compute, loopdeps, reduceddeps, [ropⱼ, oldop])
        push!(operations(rls), op)
        ropsargsᵢ[j] = op
    end
    get_adj_input_op!(∂ls, ropsargsᵢ, i)
end

function update_store_tracker!(∂ls::∂LoopSet, i::Int, oldop::Operation)
    @unpack ropsargs, tracked_ops, rls = ∂ls
    adjarrayname = diffsym(name(oldop.ref))
    # @unpack model, ropsargs, tracked_ops, rls = ∂ls
    # if !isdefined(model, adjarrayname)
    #     define!(model, adjarrayname)
    #     return
    # end
    # adjarrayname is already defined, but now we load from it.
    loadadj = Operation(length(operations(rls)), gensym(:adjointload), 8, :getindex, memload, loopdependencies(oldop), reduceddependencies(oldop), NOPARENTS, adjref(oldop), reducedchildren(oldop))
    push!(ropsargs[i], loadadj)
    add_load!(rls, loadadj)
    nothing
end

# This function handles combining multiples
# and determining reductions
function get_adj_input_op!(∂ls::∂LoopSet, i::Int)
    @unpack rls, lsold, ropsargs, diffops = ∂ls
    oldop = operations(lsold)[i]
    if isload(oldop)
        update_store_tracker!(∂ls, i, oldop)
    end
    oldop_ld = loopdependencies(oldop)
    ndeps = length(oldop_ld)
    # This function must compute reduction
    ropsargsᵢ = ropsargs[i]
    # @show rls.operations
    # @show ropsargs oldop ropsargsᵢ, i
    if maximum(length ∘ loopdependencies, ropsargsᵢ) > ndeps
        get_adj_input_op_reduct!(∂ls, ropsargsᵢ, i)
    elseif minimum(length ∘ loopdependencies, ropsargsᵢ) < ndeps
        get_adj_input_op_expand!(∂ls, ropsargsᵢ, i)
    else
        get_adj_input_op!(∂ls, ropsargsᵢ, i)
    end
end

# This function indexes the operations vector of the DiffRuleOperation to return
# the operation corresponding to dep.
# It checks if that operation is from the reverse-loop. If it is not (if it is from the forward loop)
# it add a load operation to the reverse loop and return that.
# If the operation must be loaded in the reverse loop, it also checks if it was already stored in the forward loop.
# If it wasn't, it adds a store operation to the forward loop.
function get_op_from_deriv_op_vec!(∂ls::∂LoopSet, dro::DiffRuleOperation, dep::Int)
    @unpack fls, rls, lsold, ropsargs, diffops = ∂ls
    sects = sections(dro)
    drops = operations(dro)
    fls_boundary = first(sects[2])
    dop = drops[dep]
    if !iszero(dep) && dep < fls_boundary
        # Then we need to check if the operation has been stored. If not, we must add a store to the previous loop.
        dopid = identifier(dop)
        storeid = ∂ls.stored_ops[dopid]
        # @show storeid dop
        ldref = loopdependencies(dop)
        if isload(dop)
            mref = dop.ref
        elseif isconstant(dop)
            constop = Operation(length(operations(rls)), name(dop), 8, instruction(dop), constant, loopdependencies(dop), reduceddependencies(dop), NOPARENTS, NOTAREFERENCE, reducedchildren(dop))
            dopid = identifier(dop)
            constid = identifier(constop)
            if (id = findfirst(ids -> first(ids) == dopid, lsold.preamble_symsym); !isnothing(id))
                push!(rls.preamble_symsym, (constid, last(lsold.preamble_symsym[id])))
            elseif (id = findfirst(ids -> first(ids) == dopid, lsold.preamble_symint); !isnothing(id))
                push!(rls.preamble_symint, (constid, last(lsold.preamble_symint[id])))
            elseif (id = findfirst(ids -> first(ids) == dopid, lsold.preamble_symfloat); !isnothing(id))
                push!(rls.preamble_symfloat, (constid, last(lsold.preamble_symfloat[id])))
            elseif (id = findfirst(ids -> first(ids) == dopid, lsold.preamble_zeros); !isnothing(id))
                push!(rls.preamble_zeros, (constid, last(lsold.preamble_zeros[id])))
            else
                @show dop
                throw("Constant value not found.")
            end
            return LoopVectorization.pushop!(rls, constop)
        elseif storeid == -1
            storeop = LoopVectorization.add_simple_store!(fls, dop, ArrayReference(gensym(:temporaryarray),ldref), 8)
            mref = storeop.ref
            push!(∂ls.temparrays, mref)
            LoopVectorization.maybeaddref!(fls, storeop)
            storeid = identifier(storeop)
            ∂ls.stored_ops[identifier(dop)] = storeid
        else
            mref = operations(fls)[storeid].ref
        end
        # Then we must cse-load it.
        dop = LoopVectorization.add_simple_load!(rls, name(dop), mref, ldref, 8)
        LoopVectorization.maybeaddref!(rls, dop)
    end
    dop
end
function mergedeps(ops::Vector{Operation})
    ldref = copy(loopdependencies(first(ops)))
    rdref = copy(reduceddependencies(first(ops)))
    for op ∈ @view(ops[2:end])
        LoopVectorization.mergesetv!(ldref, loopdependencies(op))
        LoopVectorization.mergesetv!(rdref, reduceddependencies(op))
    end
    ldref, rdref
end

function add_compute_reverse!(∂ls::∂LoopSet, i::Int)
    istracked(∂ls, i) || return
    @unpack rls, lsold, ropsargs, diffops = ∂ls
    oldop = operations(lsold)[i]
    rops = operations(rls)
    dro = diffops[i]
    # @show i, dro
    retinds = returned_inds(dro)
    # dro 0 can now be filled in
    sects = sections(dro)
    num_reverse_sections = length(sects) - 1
    nargs = length(sects) - 2
    drops = operations(dro)
    drops[0] = get_adj_input_op!(∂ls, i)
    # println("")
    # @show rops
    # @show sects
    for k ∈ 2:length(sects)
        # check if the final return of the section is tracked, if not, continue
        # k > 2 && @show drops[k - nargs - 3]
        (k > 2 && !istracked(∂ls, drops[k - nargs - 3])) && continue
        for j ∈ sects[k]
            instrⱼ, depsⱼ = dro[j]
            # @show j, instrⱼ, depsⱼ
            if instrⱼ === :adjoint
                drops[j] = drops[first(depsⱼ)]
                continue
            end
            parentsⱼ = Vector{Operation}(undef, length(depsⱼ))
            for (d,dep) ∈ enumerate(depsⱼ)
                parentsⱼ[d] = get_op_from_deriv_op_vec!(∂ls, dro, dep)
            end
            ldref, rdref = mergedeps(parentsⱼ)
            op = Operation(length(rops), gensym(:reverseop), 8, instrⱼ, compute, ldref, rdref, parentsⱼ)
            drops[j] = LoopVectorization.pushop!(rls, op)
        end
        if k > 2
            # We have to get_op_from_deriv_op_vec, because the return may have been computed in the first pass; e.g. exp
            add_reverse_input_parent!(ropsargs, oldop, get_op_from_deriv_op_vec!(∂ls, dro, retinds[k-1]), k-2)
        end
    end
end


text:
using Logging
using Crayons.Box
using LightGraphs, SimpleWeightedGraphs
using JuMP
using GLPK, SCIP, CPLEX
using Random
using ProgressMeter
using PrettyTables

include("common.jl")
include("heuristic.jl")
include("util.jl")
include("mtz.jl")
include("scf.jl")
include("mcf.jl")
include("cec.jl")
include("dcc.jl")

@enum HeaderName Filename Graphsize Treesize Opt Sol

function colot(name::HeaderName;
                        f::Union{Int, Nothing}=nothing,
                        ks::Union{Int, Nothing}=nothing,
                        k::Union{Int, Nothing}=nothing,
                        ms::Union{Int, Nothing}=nothing,
                        m::Union{Int, Nothing}=nothing,
                        sols::Union{Int, Nothing}=nothing,
                        sol::Union{Int, Nothing}=nothing,
                        line::Union{Int, Nothing}=nothing)
    return columnoftable(name, f=f, ks=ks, k=k, ms=ms, m=m, sols=sols, sol=sol, line=line)
end

function columnoftable(name::HeaderName;
                        f::Union{Int, Nothing}=nothing,
                        ks::Union{Int, Nothing}=nothing,
                        k::Union{Int, Nothing}=nothing,
                        ms::Union{Int, Nothing}=nothing,
                        m::Union{Int, Nothing}=nothing,
                        sols::Union{Int, Nothing}=nothing,
                        sol::Union{Int, Nothing}=nothing,
                        line::Union{Int, Nothing}=nothing)
    if name == Filename
        return 1
    elseif name == Graphsize
        return 2
    elseif name == Treesize
        return 3
    elseif name == Opt
        return 4
    elseif name == Sol
        @assert !isnothing(sols) "sols is not defined."
        @assert !isnothing(sol) "sol is not defined."
        @assert !isnothing(ms) "ms is not defined."
        @assert !isnothing(m) "m is not defined."
        return 4 + Int(sols != 1) + (Int(sols != 1) + ms) * (sol - 1) + m
    else
        @assert false "Defensive assert: Illegal case."
    end
end

function rowot(name::HeaderName;
                    f::Union{Int, Nothing}=nothing,
                    ks::Union{Int, Nothing}=nothing,
                    k::Union{Int, Nothing}=nothing,
                    ms::Union{Int, Nothing}=nothing,
                    m::Union{Int, Nothing}=nothing,
                    sols::Union{Int, Nothing}=nothing,
                    sol::Union{Int, Nothing}=nothing,
                    line::Union{Int, Nothing}=nothing)
    return rowoftable(name, f=f, ks=ks, k=k, ms=ms, m=m, sols=sols, sol=sol, line=line)
end

function rowoftable(name::HeaderName;
                    f::Union{Int, Nothing}=nothing,
                    ks::Union{Int, Nothing}=nothing,
                    k::Union{Int, Nothing}=nothing,
                    ms::Union{Int, Nothing}=nothing,
                    m::Union{Int, Nothing}=nothing,
                    sols::Union{Int, Nothing}=nothing,
                    sol::Union{Int, Nothing}=nothing,
                    line::Union{Int, Nothing}=nothing)
    lines = 4
    if name == Filename || name == Graphsize
        @assert !isnothing(f) "f is not defined."
        @assert !isnothing(ks) "ks is not defined."
        return 1 + (f - 1) * ks * lines
    elseif name == Treesize || name == Opt
        @assert !isnothing(f) "f is not defined."
        @assert !isnothing(ks) "ks is not defined."
        @assert !isnothing(k) "k is not defined."
        return 1 + (f - 1) * ks * lines + (k - 1) * lines
    elseif name == Sol
        @assert !isnothing(f) "f is not defined."
        @assert !isnothing(ks) "ks is not defined."
        @assert !isnothing(k) "k is not defined."
        @assert !isnothing(line) "line is not defined."
        return 1 + (f - 1) * ks * lines + (k - 1) * lines + (line - 1)
    else
        @assert false "Defensive assert: Illegal case."
    end
end

function solver_from_string(x::AbstractString)
    i::Int = -1
    xl = lowercase(x)
    for solver in instances(Solver)
        if string(solver) == xl
            i = Int(solver)
            break
        end
    end
    if i == -1
        error("No solver with name $x")
    end
    return Solver(i)
end

function solvingmode_from_string(x::AbstractString)
    i::Int = -1
    for mode in instances(SolvingMode)
        if string(mode) == x
            i = Int(mode)
            break
        end
    end
    if i == -1
        error("No mode with name $x")
    end
    return SolvingMode(i)
end

function main(;files::Vector{String}=[],
                directories::Vector{String}=[],
                modestrings::Vector{String}=["mtz"],
                kstrings::Vector{String}=["1.0"],
                solverstrings::Vector{String}=["cplex"],
                timeout_sec::Float64=Inf64,
                generate_timeout_sec::Float64=Inf64,
                preprocessinstances::Bool=true,
                preprocesstimeout::Float64=Inf64,
                printsolutiongraphs::Bool=true,
                printintermediatetable::Bool=true,
                debugmodels::Bool=true,
                intermediatetableinterval::Float64=60.0,
                tablestyle::String="text"
                )
    fl::Int = length(files)
    dl::Int = length(directories)
    modes::Vector{SolvingMode} = map(solvingmode_from_string, modestrings)
    ml::Int = length(modes)
    ks::Vector{Union{Float64,Int}} = map(s -> begin
        v::Float64 = parse(Float64, s)
        if v > 1 && v == trunc(v)
        return Int(trunc(Int, v))
    end
        if !(0 <= v <= 1)
        error("Specified illegal k=$s. Should be Float ∈ [0,1] or Int if > 1.")
    end
        return v
    end, kstrings)
    kl::Int = length(ks)
    solvers::Vector{Solver} = map(solver_from_string, solverstrings)
    sl::Int = length(solvers)
    @debug "Got $(fl) file$(fl == 1 ? "" : "s"), $(dl) $(dl == 1 ? "directory" : "directories"), $(ml) mode$(ml == 1 ? "" : "s"), $(kl) size$(kl == 1 ? "" : "s"), $(sl) solver$(sl == 1 ? "" : "s"), $(timeout == Inf64 ? "no" : format_seconds_readable(timeout)) timeout."
    fl = enumeratefiles!(files, directories, fl=fl)
    if fl == 0
        return
    end
    @assert fl > 0
    tablestyle = lowercase(tablestyle)
    if tablestyle == "text"
        tablebackend = :text
    elseif tablestyle == "html"
        tablebackend = :html
    elseif tablestyle == "latex"
        tablebackend = :latex
    else
        error("Unknown tablestyle: $tablestyle")
    end
    instances::Int = fl * kl * ml * sl
    if timeout != Inf64
        @debug "$instances instance$(instances == 1 ? "" : "s") to solve, $(format_seconds_readable(instances * timeout)) runtime upperbound."
    end
    (header, data) = generatetable(files=files, fl=fl, solverstrings=solverstrings, sl=sl, modestrings=modestrings, ml=ml, ks=ks, kl=kl)
    kMSTs!(files=files,
            fl=fl,
            sizes=ks,
            kl=kl,
            solvers=solvers,
            sl=sl,
            modes=modes,
            ml=ml,
            timeout_sec=timeout_sec,
            generate_timeout_sec=generate_timeout_sec,
            header=header,
            data=data,
            preprocessinstances=preprocessinstances,
            preprocesstimeout=preprocesstimeout,
            printintermediatetable=printintermediatetable,
            intermediatetableinterval=intermediatetableinterval,
            printsolutiongraphs=printsolutiongraphs,
            tablebackend=tablebackend)
    return
end

function enumeratefiles!(files::Vector{String}=[], directories::Vector{String}=[];fl::Int=length(files))::Int
    enumfiles = ProgressUnknown(desc="Enumerating files:", dt=0.5, enabled=isdebug())
    ProgressMeter.update!(enumfiles, fl)
    while !isempty(directories)
        directory = pop!(directories)
        for (root, dirs, fs) in walkdir(directory)
            for dir in dirs
                push!(directories, dir)
            end
            for file in fs
                if endswith(file, ".dat")
                    push!(files, joinpath(root, file))
                    ProgressMeter.next!(enumfiles)
                end
    end
        end
    end
    ProgressMeter.finish!(enumfiles)
    fl = length(files)
    @debug "Found $fl file$(fl == 1 ? "" : "s")."
    return fl
end

function generatetable(;files::Vector{String}=[],
                        fl::Int=length(files),
                        solverstrings::Vector{String}=[],
                        sl::Int=length(solverstrings),
                        modestrings::Vector{String}=[],
                        ml::Int=length(modestrings),
                        ks::Vector{Union{Int, Float64}} = [],
                        kl::Int=length(ks))::Tuple{Vector{String},Matrix{Union{Int, Float64,String}}}
    header = ["Graph", "|V|", "k", "OPT"]
    if sl > 1
        for s in solverstrings
            push!(header, s)
            append!(header, modestrings)
        end
    else
        append!(header, modestrings)
    end
    lines = 4
    data::Matrix{Union{String, Int, Float64}} = fill("", fl * kl * lines, length(header))
    for (f, file) in enumerate(files)
        data[rowot(Filename, f=f, ks=kl), colot(Filename)] = file
        for (k, size) in enumerate(ks)
            data[rowot(Treesize, f=f, k=k, ks=kl), colot(Treesize)] = size
        end
    end
    return (header, data)
end

function kMSTs!(;
    files::Vector{String} = [],
    fl::Int=length(files),
    sizes::Vector{Union{Float64, Int}} = [],
    kl::Int=length(sizes),
    solvers::Vector{Solver} = [],
    sl::Int=length(solvers),
    modes::Vector{SolvingMode} = [],
    ml::Int=length(modes),
    timeout_sec::Float64 = Inf64,
    generate_timeout_sec::Float64 = Inf64,
    header::Vector{String}=[],
    data::Matrix{Union{String, Int, Float64}}=[],
    preprocessinstances::Bool=true,
    preprocesstimeout::Float64=Inf64,
    printintermediatetable::Bool = true,
    intermediatetableinterval::Float64 = 60.0,
    printsolutiongraphs::Bool=true,
    tablebackend::Symbol=:text
    )
    currenttime::Float64 = 0.0
    for (f, file) in enumerate(files)
        graph::SimpleWeightedGraph = read_file_as_simplegraph(file)
        n::Int = nv(graph) - 1
        data[rowot(Graphsize, f=f, ks=kl), colot(Graphsize, f=f, ks=kl)] = n
        for (k, size) in enumerate(sizes)
            treesize :: Int = 0
            if 0 <= size <= 1
                treesize = Int(floor(Int, (n * size)))
            else
                treesize = size
            end
            data[rowot(Treesize, f=f, k=k, ks=kl), colot(Treesize, f=f, k=k, ks=kl)] = treesize
            initialsolution::Union{Vector{Edge{Int}}, Nothing} = nothing
            lb::Int = 0
            ub::Int = typemax(Int)
            if preprocessinstances
                @debug "Preprocessing instance $file ($treesize nodes)."
                preprocesstime = @elapsed (initialsolution, lb, ub) = preprocess(graph, treesize, preprocesstimeout)
                @debug "Preprocessed instance $file ($treesize nodes) with lb: $(string(CYAN_FG(string(lb)))), ub: $(string(MAGENTA_FG(string(ub)))), rg: $(string(GREEN_FG(format_ratio_readable(1-(lb/ub))))), in $(format_seconds_readable(preprocesstime))."
            end
            opt::Union{Nothing,Int} = nothing
            for (s, solver) in enumerate(solvers)
                for (m, mode) in enumerate(modes)
                    line1::Union{String, Int, Float64} = ""
                    line2::Union{String, Int, Float64} = ""
                    line3::Union{String, Int, Float64} = ""
                    line4::Union{String, Int, Float64} = ""
                    @debug "Starting instance $file, with size $treesize, solver $solver and mode $mode."
                    kmstsolutionreport :: KMSTSolutionReport = kMST(graph, mode, treesize, solver, timeout_sec=timeout_sec, generate_timeout_sec=generate_timeout_sec, initialsolution=initialsolution, lowerbound=lb, upperbound=ub, printsolutiongraphs=printsolutiongraphs, debugmodels=debugmodels)
                    infostring = "Solved instance $file, with size $treesize, solver $solver and mode $mode."
                    if kmstsolutionreport.termination_status != MOI.OPTIMAL
                        infostring *=  " " * string(RED_FG("Did not find an optimal solution."))
                        if !isnothing(kmstsolutionreport.objective_value)
                            line1 = kmstsolutionreport.objective_value
                        end
                        if !isnothing(kmstsolutionreport.relative_gap)
                            line2 = format_ratio_readable(kmstsolutionreport.relative_gap)
                        end
                    else
                        line1 = format_seconds_readable(kmstsolutionreport.solve_time_sec)
                        @assert isnothing(opt) || opt == kmstsolutionreport.objective_value "Optimal value changed from $opt to $(kmstsolutionreport.objective_value) at $file with size $(treesize), solver $solver, mode $mode."
                        opt = kmstsolutionreport.objective_value
                    end
                    infostring *= " In $(format_seconds_readable(kmstsolutionreport.solve_time_sec)), with weight $(kmstsolutionreport.objective_value)"
                    if !isnothing(kmstsolutionreport.bnb_nodes)
                        line3 = kmstsolutionreport.bnb_nodes
                        infostring *= ", with $(kmstsolutionreport.bnb_nodes) BNB node$(kmstsolutionreport.bnb_nodes == 1 ? "" : "s")"
                    end
                    if !isnothing(kmstsolutionreport.violated_ineq)
                        line4 = kmstsolutionreport.violated_ineq
                        infostring *= ", with $(kmstsolutionreport.violated_ineq) violated inequalit$(kmstsolutionreport.bnb_nodes == 1 ? "y" : "ies")"
                    end
                    infostring *= "."
                    @info infostring
                    data[rowot(Sol, f=f, k=k, ks=kl, sol=s, sols=sl, m=m, ms=ml, line=1), colot(Sol, f=f, k=k, ks=kl, sol=s, sols=sl, m=m, ms=ml, line=1)] = line1
                    data[rowot(Sol, f=f, k=k, ks=kl, sol=s, sols=sl, m=m, ms=ml, line=2), colot(Sol, f=f, k=k, ks=kl, sol=s, sols=sl, m=m, ms=ml, line=2)] = line2
                    data[rowot(Sol, f=f, k=k, ks=kl, sol=s, sols=sl, m=m, ms=ml, line=3), colot(Sol, f=f, k=k, ks=kl, sol=s, sols=sl, m=m, ms=ml, line=3)] = line3
                    data[rowot(Sol, f=f, k=k, ks=kl, sol=s, sols=sl, m=m, ms=ml, line=4), colot(Sol, f=f, k=k, ks=kl, sol=s, sols=sl, m=m, ms=ml, line=4)] = line4
                    currenttime += kmstsolutionreport.solve_time_sec
                    if printintermediatetable && currenttime > intermediatetableinterval && (f != fl || k != kl || s != sl || m != ml)
                        pretty_table(data, header=header, backend=Val(tablebackend))
                        currenttime = 0.0
                    end
                    Base.GC.gc()
                end
            end
            o = "Unkn."
            if !isnothing(opt)
                o = opt
            end
            data[rowot(Opt, f=f, k=k, ks=kl, ms=ml, sols=sl), colot(Opt, f=f, k=k, ks=kl, ms=ml, sols=sl)] = o
        end
    end
    if tablebackend == :html
        pretty_table(data, header=header, backend=Val(tablebackend))
    else
        pretty_table(data, header=header, backend=Val(tablebackend), body_hlines=Int[4 * i for i in 1:(fl*kl)])
    end
end

struct KMSTSolutionReport
    termination_status :: MOI.TerminationStatusCode
    objective_value :: Union{Int, Nothing}
    relative_gap :: Union{Float64, Nothing}
    solve_time_sec :: Float64
    bnb_nodes :: Union{Int, Nothing}
    violated_ineq :: Union{Int, Nothing}
end

function kMST(graph :: SimpleWeightedGraph,
    mode::SolvingMode,
    k::Int,
    solver::Solver
    ;
    timeout_sec::Float64=Inf64,
    generate_timeout_sec::Float64=Inf64,
    initialsolution :: Union{Vector{Edge{Int}}, Nothing}=nothing,
    lowerbound::Int=0,
    upperbound::Int=typemax(Int),
    printsolutiongraphs::Bool=true,
    debugmodels::Bool=true) :: KMSTSolutionReport
    n::Int = nv(graph) - 1
    prelude_time::Float64 = 0
    prelude_time += @elapsed begin
        model = generate_model(graph, mode, k, solver, timeout_sec=timeout_sec, generate_timeout_sec=generate_timeout_sec, lowerbound=lowerbound, upperbound=upperbound,debugmodels=debugmodels)
    end
    if generate_timeout_sec != Inf64 && prelude_time > generate_timeout_sec
        return KMSTSolutionReport(MOI.TIME_LIMIT, upperbound, 1-(lowerbound/upperbound), prelude_time, 0, nothing)
    end
    nvars::Int = num_variables(model)
    if !isnothing(initialsolution) && solver != glpk
        prelude_time += @elapsed warmstart_model!(model, graph, mode, k, initialsolution)
    elseif isdebug()
        @debug "Model has $nvars variable$(nvars == 1 ? "" : "s")."
    end
    if isdebug()
        all_vars::Vector{VariableRef} = all_variables(model)
        bound::Int = count(v -> is_fixed(v) || (has_lower_bound(v) && has_upper_bound(v)), all_vars)
        if bound == nvars
            @debug "Model has no unbound variables."
        else
            free::Int = count(v -> !is_fixed(v) && !has_lower_bound(v) && !has_upper_bound(v), all_vars)
            onlyupperbound::Int = count(v -> !is_fixed(v) && !has_lower_bound(v) && has_upper_bound(v), all_vars)
            onlylowerbound::Int = count(v -> !is_fixed(v) && !has_upper_bound(v) && has_lower_bound(v), all_vars)
            @debug "Model has $bound ($(format_ratio_readable(bound,nvars))) bound variable$(bound == 1 ? "" : "s" ), $free ($(format_ratio_readable(free, nvars))) free variable$(free == 1 ? "" : "s"), $onlylowerbound ($(format_ratio_readable(onlylowerbound, nvars))) variable$(onlylowerbound == 1 ? "" : "s") with no upperbound, $onlyupperbound ($(format_ratio_readable(onlyupperbound, nvars))) variable$(onlyupperbound == 1 ? "" : "s") with no lowerbound."
        end

        constypes::Vector{Tuple{DataType, DataType}} = list_of_constraint_types(model)
        ncons::Int = 0
        debugstring::String = ""
        numcons::Dict{DataType, Int} = SortedDict()
        for (reftype, constype) in constypes
            c::Int = num_constraints(model, reftype, constype)
            if !(constype in keys(numcons))
                numcons[constype] = 0
            end
            numcons[constype] += c

            ncons += c
        end
        for (constype, c) in numcons
            if constype == MOI.LessThan{Float64}
                debugstring *= " $c ($(format_ratio_readable(c, ncons))) of type '≤'."
            elseif constype == MOI.GreaterThan{Float64}
                debugstring *= " $c ($(format_ratio_readable(c, ncons))) of type '≥'."
            elseif constype == MOI.EqualTo{Float64}
                debugstring *= " $c ($(format_ratio_readable(c, ncons))) of type '='."
            elseif constype == MOI.ZeroOne
                debugstring *= " $c ($(format_ratio_readable(c, ncons))) of type 'binary'."
            elseif constype == MOI.Integer
                debugstring *= " $c ($(format_ratio_readable(c, ncons))) of type 'integer'."
            elseif constype == MOI.Interval
                debugstring *= " $c ($(format_ratio_readable(c, ncons))) of type 'interval'."
            else
                debugstring *= " $c ($(format_ratio_readable(c, ncons))) of type '$constype'."
            end
        end
        @debug "Model has $ncons constraint$(ncons == 1 ? "" : "s")." * debugstring
    end

    @debug "Generated model in $(format_seconds_readable(prelude_time))."
    if generate_timeout_sec != Inf64 && prelude_time > generate_timeout_sec
        return KMSTSolutionReport(MOI.TIME_LIMIT, upperbound, 1-(lowerbound/upperbound), prelude_time, 0, nothing)
    end

    if mode == mcf && n >= 205
        return KMSTSolutionReport(MOI.MEMORY_LIMIT, upperbound, 1-(lowerbound/upperbound), prelude_time, 0, nothing)
    end

    kmstsolution = solve!(model, graph, k, mode)
    if printsolutiongraphs && !isnothing(kmstsolution.solution_graph)
        print_weighted_graph(kmstsolution.solution_graph)
    end

    objective_value :: Union{Int, Float64, Nothing} = kmstsolution.objective_value
    if objective_value == Inf64
        objective_value = nothing
    else
        objective_value = Int(round(Int, objective_value))
    end

    return KMSTSolutionReport(
        kmstsolution.termination_status,
        objective_value,
        kmstsolution.relative_gap,
        kmstsolution.solve_time_sec,
        kmstsolution.bnb_nodes,
        kmstsolution.violated_ineq)
end

function preprocess(graph::SimpleWeightedGraph, k::Int, timeout_sec::Float64=Inf64) :: Tuple{Vector{Edge{Int}}, Int, Int}
    nthreads :: Int = Base.Threads.nthreads()
    calc_time :: Float64 = 0

    lb::Int = 0
    ub::Int = typemax(Int)

    n::Int = nv(graph)

    lbfuture = Base.Threads.@spawn begin
    checked :: Set{Edge{Int}} = Set()
    tocheck_lock = Base.ReentrantLock()
    tocheck :: Vector{Edge{Int}} = []
    lowerbound_lock = Base.ReentrantLock()
    (lbtree, lb, s) = kruskal_heuristic(graph, k)
    append!(tocheck, setdiff(s, checked))
    for e in lbtree
        push!(checked, e)
    end

    while !isempty(tocheck)
        Threads.@threads for i = 1:max(1,div(nthreads-1,2))
            if !isempty(tocheck)
                unskippable = nothing
                lock(tocheck_lock)
                try
                    if !isempty(tocheck)
                        unskippable = pop!(tocheck)
                        while unskippable in checked && !isempty(tocheck)
                            unskippable = pop!(tocheck)
                        end
                        if unskippable in checked
                            unskippable = nothing
                        end
                    end
                finally
                    unlock(tocheck_lock)
                end
                if !isnothing(unskippable)
                    push!(checked, unskippable)
                    (lbtreecand, lbcand, s) = kruskal_heuristic(graph, k, unskippable=unskippable, bound=lb)
                    if !isnothing(lbtreecand)
                        for e in lbtreecand
                            push!(checked, e)
                        end
                    end
                    append!(tocheck, setdiff(s, checked))
                    if lbcand < lb
                        lock(lowerbound_lock)
                        try
                            if lbcand < lb
                                lbtree = lbtree
                                lb = lbcand
                            end
                        finally
                            unlock(lowerbound_lock)
                        end
                    end
                    if isempty(s)
                        lock(tocheck_lock)
                        try
                            empty!(tocheck)
                        finally
                            unlock(tocheck_lock)
                        end
                    end
                end
            end
        end
    end
    @assert !isnothing(lbtree) "No lower bound tree found."
    lb
    end

    ubtree::Vector{Edge{Int}} = []

    startnodes_lock = Base.ReentrantLock()
    upperbound_lock = Base.ReentrantLock()

    ubfuture = Base.Threads.@spawn begin
    startnodes::Vector{Int} = collect(Int,2:n)
    shuffle!(startnodes)
    checked :: Set{Int} = Set()

    while calc_time < timeout_sec && !isempty(startnodes)
        calc_time += @elapsed begin
            Threads.@threads for i = 1:max(1,div(nthreads-1,2)+((nthreads-1)%2))
                if !isempty(startnodes)
                    lock(startnodes_lock)
                    startnode = nothing
                    try
                        if !isempty(startnodes)
                            startnode = pop!(startnodes)
                            while startnode in checked && !isempty(startnodes)
                                startnode = pop!(startnodes)
                            end
                            if startnode in checked
                                startnode = nothing
                            end
                        end
                    finally
                        unlock(startnodes_lock)
                    end
                    if !isnothing(startnode)
                        heur::Union{Tuple{Vector{Edge{Int}}, Int}, Nothing} = prim_heuristic(graph, k, startnode=startnode, upperbound=ub)
                        if !isnothing(heur)
                            for e in heur[1]
                                for v in [src(e), dst(e)]
                                    push!(checked, v)
                                end
                            end
                            if heur[2] < ub
                                lock(upperbound_lock)
                                try
                                    if heur[2] < ub
                                        (ubtree, ub) = heur
                                    end
                                finally
                                    unlock(upperbound_lock)
                                end
                            end
                        end
                    end
                end
            end
        end
    end
    (ubtree,ub)
    end

    lb = fetch(lbfuture)
    (ubtree,ub) = fetch(ubfuture)

    @assert lb <= ub "Lowerbound has to be smaller than upperbound. !($lb < $ub)"

    return (ubtree, lb, ub)
end

function generate_model(graph::SimpleWeightedGraph,
        mode::SolvingMode,
        k::Int,
        solver::Solver;
        timeout_sec::Float64 = Inf64,
        generate_timeout_sec::Float64 = Inf64,
        lowerbound::Int = 0,
        upperbound::Int = typemax(Int),
        debugmodels::Bool=true)
    model = Model()
    if solver == glpk
        set_optimizer(model, GLPK.Optimizer)
    elseif solver == cplex
        set_optimizer(model, CPLEX.Optimizer)
    elseif solver == scip
        set_optimizer(model, SCIP.Optimizer)
    else
        @assert false "Defensive assert: Unhandled solver $solver."
    end
    if !isdebug()
        set_silent(model)
    else
        unset_silent(model)
        if solver == glpk
            set_optimizer_attribute(model, "msg_lev", GLPK.GLP_MSG_ON)
        end
    end
    if timeout_sec != Inf64
        if timeout < 0
            error("Timeout should be positive")
        end
        set_time_limit_sec(model, timeout)
    end
    if solver == cplex
        available_memory::Int = round(Int, Base.Sys.total_memory() // 2^20)
        work_memory::Int = max(1024, min(1024 * 8, div(available_memory, 3)))
        tree_limit::Int = max(1024 * 2, max(work_memory, div(9 * available_memory, 10)))
        @debug "There $(available_memory == 1 ? "is" : "are") $available_memory MiB available. Limiting work memory to $work_memory MiB, and setting tree limit to $tree_limit MiB."
        #set_optimizer_attribute(model, "CPX_PARAM_NODEFILEIND", 3)
        set_optimizer_attribute(model, "CPX_PARAM_TRELIM", tree_limit)
        set_optimizer_attribute(model, "CPX_PARAM_WORKMEM", work_memory)
        set_optimizer_attribute(model, "CPX_PARAM_THREADS", 1)
    end
    generate_time::Float64 = 0
    @debug "Generating common variables."
    generate_time += @elapsed basic_kmst!(model, graph, k, lowerbound=lowerbound, upperbound=upperbound, isdebug=isdebug())
    if generate_time > generate_timeout_sec + 1.0
        return model
    end
    generate_time += @elapsed begin
    @debug "Generating $mode variables."
    if mode == mtz
        miller_tucker_zemlin(model, graph, k)
    elseif mode == scf
        single_commodity_flow!(model, graph, k)
    elseif mode == mcf
        multi_commodity_flow!(model, graph, k, generate_timeout_sec=generate_timeout_sec-generate_time, isdebug=isdebug())
    elseif mode == cec
        cycle_elimination_constraints!(model, graph, k)
    elseif mode == dcc
        directed_cutset_constraints!(model, graph, k)
    end
    end
    if debugmodels && generate_time < generate_timeout_sec + 1.0
        @debug model
    end
    return model
end

function warmstart_model!(model, graph::SimpleWeightedGraph, mode::SolvingMode, k::Int, solution::Vector{Edge{Int}})
    basic_kmst_warmstart!(model, graph, k, solution)
    if mode == mtz
        miller_tucker_zemlin_warmstart(model, graph, k, solution)
    end

    filled = check_kmst_warmstart!(model, graph, k)

    if mode == mtz
        filled += check_miller_tucker_zemlin_warmstart(model, graph, k)
    end

    @debug "Found $filled variable$(filled == 1 ? "" : "s") that could be filled."

    if isdebug()
        num_vars :: Int = JuMP.num_variables(model)
        vars :: Vector{VariableRef} = all_variables(model)
        start_values = start_value.(vars)
        with_startvalue :: Int = num_vars - count(isnothing, start_values)
        are_fixed::Int = count(is_fixed, vars)
        known::Int = count(v -> is_fixed(v) || !isnothing(start_value(v)), vars)
        @debug "Model has $num_vars variables, $known ($(format_ratio_readable(known, num_vars))) $(known == 1 ? "is" : "are") known, $with_startvalue ($(format_ratio_readable(with_startvalue, num_vars))) $(with_startvalue == 1 ? "has" : "have") a start value, $are_fixed ($(format_ratio_readable(are_fixed, num_vars))) $(are_fixed == 1 ? "is" : "are") fixed."
    end
end

function read_file_as_simplegraph(path::String)
    graph::Union{Nothing,SimpleWeightedGraph} = nothing
    vertexsize::Int = 0
    edgesize::Union{AbstractFloat,Int} = Inf
    edgeprog = nothing
    totaltime, totallines = open(path) do f
        linecounter::Int = 0
        timetaken = @elapsed for line in eachline(f)
            if linecounter <= 1
                if linecounter == 0
                    vertexsize = parse(Int, line)
                    graph = SimpleWeightedGraph(vertexsize)
                else
                    edgesize = parse(Int, line)
                    edgeProg = Progress(edgesize, dt=0.5,
                    barglyphs=BarGlyphs('[', '#', ['-','~','+','*','=','>','#'], ' ', ']'),
                    desc="Parsing edges:",
                    enabled=isdebug())
                end
            else
                elemsstrings = split(line, " ", limit=4)
                elems = [parse(Int, v) for v in elemsstrings]
                index = elems[1]
                node1 = elems[2]
                node2 = elems[3]
                weight = elems[4]
                if weight > 0
                    add_edge!(graph, node1 + 1, node2 + 1, weight)
                else
                    add_edge!(graph, node1 + 1, node2 + 1, min(1.0 / edgesize, 0.5 - (eps() * 2.0)))
                end
            end
            linecounter += 1
            if linecounter >= (edgesize + 2)
                break
            end
        end
        (timetaken, linecounter)
    end

    if totallines != 1
        linesplural = "s"
    else
        linesplural = ""
    end

    @debug ("Read $(string(YELLOW_FG(string(totallines)))) line$(linesplural) in $(string(BLUE_FG(format_seconds_readable(totaltime, 1)))) at $(string(CYAN_FG(format_seconds_readable(totaltime / totallines, 2)))) per line")
    es = ne(graph)
    if es != edgesize
        @warn "Input file claims $edgesize edge$(edgesize == 1 ? "" : "s") but got $es:"
        print_weighted_graph(graph, nothing, nv(graph), es)
    end
    return graph
end

max_stars_repo_name: Entze/jkMST

max_issues_repo_name: Entze/jkMST

max_forks_repo_name: Entze/jkMST


text:
function johnson_shortest_paths(g::AbstractGraph{U},
distmx::AbstractMatrix{T}=weights(g)) where T <: Real where U <: Integer

    nvg = nv(g)
    type_distmx = typeof(distmx)
#Change when parallel implementation of Bellman Ford available
    wt_transform = bellman_ford_shortest_paths(g, vertices(g), distmx).dists

    if !type_distmx.mutable && type_distmx !=  LightGraphs.DefaultDistance
        distmx = sparse(distmx) #Change reference, not value
    end

#Weight transform not needed if all weights are positive.
    if type_distmx !=  LightGraphs.DefaultDistance
        for e in edges(g)
            distmx[src(e), dst(e)] += wt_transform[src(e)] - wt_transform[dst(e)]
        end
    end

    dijk_state = Parallel.dijkstra_shortest_paths(g, vertices(g), distmx)
    dists = dijk_state.dists
    parents = dijk_state.parents

    broadcast!(-, dists, dists, wt_transform)
    for v in vertices(g)
        dists[:, v] .+= wt_transform[v] #Vertical traversal prefered
    end

    if type_distmx.mutable
        for e in edges(g)
            distmx[src(e), dst(e)] += wt_transform[dst(e)] - wt_transform[src(e)]
        end
    end

    return JohnsonState(dists, parents)
end

function enumerate_paths(s::JohnsonState{T,U}, v::Integer) where T <: Real where U <: Integer
    pathinfo = s.parents[v, :]
    paths = Vector{Vector{U}}()
    for i in 1:length(pathinfo)
        if (i == v) || (s.dists[v, i] == typemax(T))
            push!(paths, Vector{U}())
        else
            path = Vector{U}()
            currpathindex = i
            while currpathindex != 0
                push!(path, currpathindex)
                currpathindex = pathinfo[currpathindex]
            end
            push!(paths, reverse(path))
        end
    end
    return paths
end

enumerate_paths(s::JohnsonState) = [enumerate_paths(s, v) for v in 1:size(s.parents, 1)]
enumerate_paths(st::JohnsonState, s::Integer, d::Integer) = enumerate_paths(st, s)[d]


text:
@test ip"127.0.0.1" == IPv4(127,0,0,1)
@test ip"192.0" == IPv4(192,0,0,0)
@test ip"192.0xFFF" == IPv4(192,0,15,255)
@test ip"192.0xFFFF" == IPv4(192,0,255,255)
@test ip"192.0xFFFFF" == IPv4(192,15,255,255)
@test ip"192.0xFFFFFF" == IPv4(192,255,255,255)
@test_throws ErrorException Base.parseipv4("192.0xFFFFFFF")
@test ip"022.0.0.1" == IPv4(18,0,0,1)

@test_throws ErrorException Base.parseipv4("192.0xFFFFFFFFF")
@test_throws ErrorException Base.parseipv4("192.")

@test ip"::1" == IPv6(1)
@test ip"2605:2700:0:3::4713:93e3" == IPv6(parseint(Uint128,"260527000000000300000000471393e3",16))

@test ip"2001:db8:0:0:0:0:2:1" == ip"2001:db8::2:1" == ip"2001:db8::0:2:1"

@test ip"0:0:0:0:0:ffff:127.0.0.1" == IPv6(0xffff7f000001)

# RFC 5952 Compliance

@test repr(ip"2001:db8:0:0:0:0:2:1") == "ip\"2001:db8::2:1\""
@test repr(ip"2001:0db8::0001") == "ip\"2001:db8::1\""
@test repr(ip"2001:db8::1:1:1:1:1") == "ip\"2001:db8:0:1:1:1:1:1\""
@test repr(ip"2001:db8:0:0:1:0:0:1") == "ip\"2001:db8::1:0:0:1\""
@test repr(ip"2001:0:0:1:0:0:0:1") == "ip\"2001:0:0:1::1\""

c = Base.Condition()
@async begin
	s = listen(2134)
	Base.notify(c)
	sock = accept(s)
	write(sock,"Hello World\n")
	close(s)
	close(sock)
end
wait(c)
@test readall(connect(2134)) == "Hello World\n"

socketname = @windows ? "\\\\.\\pipe\\uv-test" : "testsocket"
@unix_only isfile(socketname) && Base.FS.unlink(socketname)
@async begin
	s = listen(socketname)
	Base.notify(c)
	sock = accept(s)
	write(sock,"Hello World\n")
	close(s)
	close(sock)
end
wait(c)
@test readall(connect(socketname)) == "Hello World\n"

@test_throws Base.UVError getaddrinfo(".invalid")
@test_throws Base.UVError connect("localhost", 21452)

server = listen(2134)
@async @test_throws ErrorException accept(server)
sleep(0.1)
close(server)

server = listen(2134)
@async connect("localhost",2134)
s1 = accept(server)
@test_throws ErrorException accept(server,s1)
close(server)

@test_throws Base.UVError connect(".invalid",80)

a = UdpSocket()
b = UdpSocket()
bind(a,ip"127.0.0.1",2134)
bind(b,ip"127.0.0.1",2135)

c = Condition()
@async begin
    @test bytestring(recv(a)) == "Hello World"
    # Issue 6505
    @async begin
        @test bytestring(recv(a)) == "Hello World"
        notify(c)
    end
    send(b,ip"127.0.0.1",2134,"Hello World")
end
send(b,ip"127.0.0.1",2134,"Hello World")
wait(c)

@test_throws MethodError bind(UdpSocket(),2134)

close(a)
close(b)

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text: pseudos(sample) = transpose(hcat([StatsBase.ordinalrank(sample[i,:])./(size(sample,2)+1) for i in 1:size(sample,1)]...))


text:
using BinDeps
using Compat

@BinDeps.setup

libnlopt = library_dependency("libnlopt", aliases=["libnlopt_cxx", "libnlopt$(Sys.WORD_SIZE)"])

provides(AptGet, "libnlopt0", libnlopt)

provides(Sources,URI("http://ab-initio.mit.edu/nlopt/nlopt-2.4.tar.gz"), libnlopt)

provides(BuildProcess,Autotools(configure_options =
    ["--enable-shared", "--without-guile", "--without-python",
    "--without-octave", "--without-matlab","--with-cxx"],
    libtarget="libnlopt_cxx.la"),libnlopt, os = :Unix)

if is_apple()
    using Homebrew
    provides( Homebrew.HB, "homebrew/science/nlopt", libnlopt, os = :Darwin )
end

nloptname = "nlopt-2.4.1"

libdir = BinDeps.libdir(libnlopt)
srcdir = BinDeps.srcdir(libnlopt)
downloadsdir = BinDeps.downloadsdir(libnlopt)
extractdir(w) = joinpath(srcdir,"w$w")
destw(w) = joinpath(libdir,"libnlopt$(w).dll")

type FileCopyRule <: BinDeps.BuildStep
    src::AbstractString
    dest::AbstractString
end
Base.run(fc::FileCopyRule) = isfile(fc.dest) || cp(fc.src, fc.dest)

provides(BuildProcess,
	(@build_steps begin
		FileDownloader("http://ab-initio.mit.edu/nlopt/$(nloptname)-dll32.zip", joinpath(downloadsdir, "$(nloptname)-dll32.zip"))
		FileDownloader("http://ab-initio.mit.edu/nlopt/$(nloptname)-dll64.zip", joinpath(downloadsdir, "$(nloptname)-dll64.zip"))
		CreateDirectory(srcdir, true)
		CreateDirectory(joinpath(srcdir,"w32"), true)
		CreateDirectory(joinpath(srcdir,"w64"), true)
		FileUnpacker(joinpath(downloadsdir,"$(nloptname)-dll32.zip"), extractdir(32), joinpath(extractdir(32),"matlab"))
		FileUnpacker(joinpath(downloadsdir,"$(nloptname)-dll64.zip"), extractdir(64), joinpath(extractdir(64),"matlab"))
		CreateDirectory(libdir, true)
        FileCopyRule(joinpath(extractdir(32),"libnlopt-0.dll"), destw(32))
        FileCopyRule(joinpath(extractdir(64),"libnlopt-0.dll"), destw(64))
	end), libnlopt, os = :Windows)

if is_windows()
    push!(BinDeps.defaults, BuildProcess)
end

@BinDeps.install Dict(:libnlopt => :libnlopt)

if is_windows()
    pop!(BinDeps.defaults)
end


text:
using HDF5, DataFrames, CSV, JLD
using StatsBase, Polynomials, CurveFit
using PyPlot

kep_filename = joinpath(Pkg.dir("ExoplanetsSysSim"), "data", "q1_q17_dr25_stellar.csv")
gaia_filename = joinpath(Pkg.dir("ExoplanetsSysSim"), "data", "gaiadr2_keplerdr25_crossref.csv")
mast_filename = joinpath(Pkg.dir("ExoplanetsSysSim"), "data", "KeplerMAST_TargetProperties.csv")
tmass_filename = joinpath(Pkg.dir("ExoplanetsSysSim"), "data", "Gaia_Kepler_2MASS_mags.csv")
stellar_catalog_file_out = joinpath(Pkg.dir("ExoplanetsSysSim"), "data", "q1q17_dr25_gaia_m_2mass.jld")

kep_df = CSV.read(kep_filename)
gaia_df = CSV.read(gaia_filename, rows_for_type_detect=30000)
mast_df = CSV.read(mast_filename)
tmass_df = CSV.read(tmass_filename, rows_for_type_detect=30000)

dup_gaiaid = find(nonunique(DataFrame(x = gaia_df[:source_id])))
deleterows!(gaia_df, dup_gaiaid)

println("Total crossref target stars = ", length(gaia_df[:kepid]))

mag_diff = gaia_df[:phot_g_mean_mag].-gaia_df[:kepmag]
quant_arr = quantile(mag_diff, [0.067,0.933])   # 1.5-sigma cut
mag_match = find(x->quant_arr[1]<=x<=quant_arr[2], mag_diff)
gaia_df = gaia_df[mag_match,:]

gaia_col = [:kepid,:source_id,:parallax,:parallax_error,:astrometric_gof_al,:astrometric_excess_noise_sig,:phot_g_mean_mag,:bp_rp,:a_g_val,:priam_flags,:teff_val,:teff_percentile_lower,:teff_percentile_upper,:radius_val,:radius_percentile_lower,:radius_percentile_upper,:lum_val,:lum_percentile_lower,:lum_percentile_upper]
df = join(kep_df, gaia_df[:,gaia_col], on=:kepid)
df = join(df, mast_df, on=:kepid)
df = join(df, tmass_df, on=:source_id, makeunique=true)
kep_df = nothing
gaia_df = nothing

println("Total target stars (KOIs) matching magnitude = ", length(df[:kepid]), " (", sum(df[:nkoi]),")")

df[:teff] = df[:teff_val]
df[:teff_err1] = df[:teff_percentile_upper].-df[:teff_val]
df[:teff_err2] = df[:teff_percentile_lower].-df[:teff_val]
delete!(df, :teff_val)
delete!(df, :teff_percentile_upper)
delete!(df, :teff_percentile_lower)
for x in 1:length(df[:kepid])
    if !isnan(df[x,:radius_val])
        df[x,:radius_err1] = df[x,:radius_percentile_upper]-df[x,:radius_val]
        df[x,:radius_err2] = df[x,:radius_percentile_lower]-df[x,:radius_val]
        df[x,:radius] = df[x,:radius_val]
    end
end
delete!(df, :radius_val)
delete!(df, :radius_percentile_upper)
delete!(df, :radius_percentile_lower)

not_binary_suspect = (df[:astrometric_gof_al] .<= 20) .& (df[:astrometric_excess_noise_sig] .<= 5)
astrometry_good = []
for x in 1:length(df[:kepid])
    if !(ismissing(df[x,:priam_flags]))
        pflag = string(df[x,:priam_flags])
         if (pflag[2] == '0') & (pflag[3] == '0') # WARNING: Assumes flag had first '0' removed by crossref script
             push!(astrometry_good, true)
         else
             push!(astrometry_good, false)
         end
     else
         push!(astrometry_good, false)
     end
end
astrometry_good = astrometry_good .& (df[:parallax_error] .< 0.3*df[:parallax])
planet_search = df[:kepmag] .<= 16.

has_mass = .! (ismissing.(df[:mass]) .| ismissing.(df[:mass_err1]) .| ismissing.(df[:mass_err2]))
has_radius = .! (ismissing.(df[:radius]) .| ismissing.(df[:radius_err1]) .| ismissing.(df[:radius_err2]))# .| isnan.(df[:radius]) .| isnan.(df[:radius_err1]) .| isnan.(df[:radius_err2]))
#has_dens = .! (ismissing.(df[:dens]) .| ismissing.(df[:dens_err1]) .| ismissing.(df[:dens_err2]))
has_cdpp = .! (ismissing.(df[:rrmscdpp01p5]) .| ismissing.(df[:rrmscdpp02p0]) .| ismissing.(df[:rrmscdpp02p5]) .| ismissing.(df[:rrmscdpp03p0]) .| ismissing.(df[:rrmscdpp03p5]) .| ismissing.(df[:rrmscdpp04p5]) .| ismissing.(df[:rrmscdpp05p0]) .| ismissing.(df[:rrmscdpp06p0]) .| ismissing.(df[:rrmscdpp07p5]) .| ismissing.(df[:rrmscdpp09p0]) .| ismissing.(df[:rrmscdpp10p5]) .| ismissing.(df[:rrmscdpp12p0]) .| ismissing.(df[:rrmscdpp12p5]) .| ismissing.(df[:rrmscdpp15p0]))
has_rest = .! (ismissing.(df[:dataspan]) .| ismissing.(df[:dutycycle]))
has_limbd = .! (ismissing.(df[:limbdark_coeff1]) .| ismissing.(df[:limbdark_coeff2]) .| ismissing.(df[:limbdark_coeff3]) .| ismissing.(df[:limbdark_coeff4]))
has_tmass = .! (ismissing.(df[:ks_msigcom]) .| ismissing.(df[:ks_m]) .| ismissing.(df[:j_m]))

mast_cut =.&(df[:numLCEXqtrs].>0,df[:numLCqtrs].>4)

is_usable = .&(has_rest, has_cdpp, mast_cut, astrometry_good, has_tmass)#, not_binary_suspect)

df = df[find(is_usable),:]
println("Total stars (KOIs) with valid parameters = ", length(df[:kepid]), " (", sum(df[:nkoi]),")")

m_absg = df[:phot_g_mean_mag] + 5 + 5*log10.(df[:parallax]/1000.)
m_absk = df[:ks_m] + 5 + 5*log10.(df[:parallax]/1000.)

m_absk_err = sqrt.(df[:ks_msigcom].^2 .+ ((5./(df[:parallax]*log(10))).^2 .* df[:parallax_error].^2))
df[:ks_absm] = m_absk
df[:ks_absm_err] = m_absk_err
m_gcolor = (1.7 .<= df[:bp_rp])# .<= 5.0)
m_tcolor = (0.617 + 0.162) .<= (df[:j_m] .- df[:ks_m])

df[:radius] = 1.9515 - (0.352*m_absk) .+ (0.0168*m_absk.^2)
m_rad_err = abs.((-0.352 + 2*0.0168*m_absk) .* m_absk_err)
df[:radius_err1] = m_rad_err
df[:radius_err2] = -m_rad_err

df[:mass] = 0.5858 + (0.3872*m_absk) .- (0.1217*m_absk.^2) .+ (0.0106*m_absk.^3) .- (2.7262e-4*m_absk.^4)
m_mass_err = abs.((0.3872 - (2*0.1217*m_absk) .+ (3*0.0106*m_absk.^2) .- (4*2.7262e-4*m_absk.^3)).*m_absk_err)
df[:mass_err1] = m_mass_err
df[:mass_err2] = -m_mass_err

# for i in 1:length(df[:kepid])
#     println(i, " ", m_absk[i], " ", df[i,:radius], " ", df[i,:mass])
# end

m_teff = df[:teff] .<= 4000
m_logg = df[:logg] .> 3
#m_teff = (2320 .<= df[:teff] .<= 3870)
m_gmag = (7.55 .<= m_absg)# .<= 16.33)
m_kmag = (4.8 .<= m_absk)

M_samp = find(m_tcolor)
println("Total cool stars (KOIs) with valid parameters = ", length(M_samp), " (", sum(df[M_samp, :nkoi]),")")
M_samp = find(.&(m_tcolor,m_kmag))
println("Total M stars (KOIs) with valid parameters = ", length(M_samp), " (", sum(df[M_samp, :nkoi]),")")
# M_samp = find(.&(m_tcolor,m_kmag, df[:radius] .< 2))#, df[:mass] .> 0.5))
# println("Total MS M stars (KOIs) with valid parameters = ", length(M_samp), " (", sum(df[M_samp, :nkoi]),")")
#println("Maximum distance (pc) = ", maximum(1.0./(df[M_samp, :parallax]/1000.0)))
#println("Minimum Ks uncertainty = ", minimum(df[M_samp, :ks_msigcom]))

plot_samp = rand(1:length(df[:kepid]), 5000)
plt[:scatter](df[M_samp, :bp_rp], m_absg[M_samp], s=3, label = "M", color="red")
plt[:scatter](df[plot_samp, :bp_rp], m_absg[plot_samp], s=3, label = "All", color="black")
# plt[:scatter](df[M_samp, :j_m] .- df[M_samp, :ks_m], m_absk[M_samp], s=3, label = "M", color="red")
# plt[:scatter](df[plot_samp, :j_m] .- df[plot_samp, :ks_m], m_absk[plot_samp], s=3, label = "All", color="black")
# plt[:scatter](df[M_samp, :radius], df[M_samp, :radius_val], s=3, label = "M", color="black")

plt[:ylabel](L"$M_{G}$")
plt[:xlabel](L"$B_p - R_p$")
# plt[:ylabel](L"$M_{K}$")
# plt[:xlabel](L"$J - K$")
# plt[:xlabel]("Kepler R_*")
# plt[:ylabel]("Mann R_*")
# plt[:xlim]((0.3, 0.8))
# plt[:ylim]((0.3, 0.8))
#plt[:xlim]((0, 1))
plt[:ylim](reverse(plt[:ylim]()))#bottom=-4)))
plt[:legend]()
# plt[:savefig]("m-dwarf_samp_radiicomp.png")
plt[:savefig]("m-dwarf_samp_gaiamag.png")
plt[:savefig]("m-dwarf_samp_gaiamag.eps")
# plt[:savefig]("m-dwarf_samp_2massmag.png")
# plt[:savefig]("m-dwarf_samp_2massmag.eps")

num_less = 0
num_mid = 0
tmp_df = df[M_samp,:]
for i in 1:length(tmp_df[:kepid])
    maxmes = max(tmp_df[i,:mesthres01p5],tmp_df[i,:mesthres02p0],tmp_df[i,:mesthres02p5],tmp_df[i,:mesthres03p0],tmp_df[i,:mesthres03p5],tmp_df[i,:mesthres04p5],tmp_df[i,:mesthres05p0],tmp_df[i,:mesthres06p0],tmp_df[i,:mesthres07p5],tmp_df[i,:mesthres09p0],tmp_df[i,:mesthres10p5],tmp_df[i,:mesthres12p0],tmp_df[i,:mesthres12p5],tmp_df[i,:mesthres15p0])
    if maxmes > 8
        num_less += 1
    elseif maxmes > 7.1
        num_mid += 1
    end
end
println("Total # of M targets with MES threshold > 8 = ", num_less)
println("Total # of M targets with MES threshold between 7.1 and 8 = ", num_mid)

# # See options at: http://exoplanetarchive.ipac.caltech.edu/docs/API_keplerstellar_columns.html
# # TODO SCI DETAIL or IMPORTANT?: Read in all CDPP's, so can interpolate?
# symbols_to_keep = [ :kepid, :source_id, :mass, :mass_err1, :mass_err2, :radius, :radius_err1, :radius_err2, :dens, :dens_err1, :dens_err2, :teff, :phot_g_mean_mag, :bp_rp, :j_m, :ks_m, :lum_val, :rrmscdpp01p5, :rrmscdpp02p0, :rrmscdpp02p5, :rrmscdpp03p0, :rrmscdpp03p5, :rrmscdpp04p5, :rrmscdpp05p0, :rrmscdpp06p0, :rrmscdpp07p5, :rrmscdpp09p0, :rrmscdpp10p5, :rrmscdpp12p0, :rrmscdpp12p5, :rrmscdpp15p0, :dataspan, :dutycycle, :limbdark_coeff1, :limbdark_coeff2, :limbdark_coeff3, :limbdark_coeff4, :contam]
# # delete!(df, [~(x in symbols_to_keep) for x in names(df)])    # delete columns that we won't be using anyway
# df = df[M_samp, symbols_to_keep]
# tmp_df = DataFrame()
# for col in names(df)
#     tmp_df[col] = collect(skipmissing(df[col]))
# end
# df = tmp_df

# save(stellar_catalog_file_out,"stellar_catalog", df)


text:
function count_set(df, cols, setcols)
    count(eachrow(df)) do row
        any(ismissing, row[cols[setcols]])        && return false # all setcols must be non-missing
        !all(row[cols[setcols]])                  && return false # all set colls must be true
        any(skipmissing(row[cols[Not(setcols)]])) && return false # no non-set cols can be true
        return true
    end
end

function upset_dots!(ax, colsets, nsets=maximum(Iterators.flatten(colsets)))
    xlims!(ax, (0, length(colsets) + 1))
    hidexdecorations!(ax)
    hideydecorations!(ax, ticklabels=false)
    hidespines!(ax)

    for (x, p) in enumerate(colsets)
        scatter!(ax, fill(x, nsets - length(p)), collect(1:nsets)[Not(p)], markersize=20, color=:lightgray)
        scatter!(ax, fill(x, length(p)), p, markersize=20, color=:black)
        length(p) > 1 && lines!(ax, [x,x], [extrema(p)...], color=:black, linewidth=3)
    end
end

# https://github.com/JuliaStats/MultivariateStats.jl/pull/162
function loadings(M::MultivariateStats.MDS)
    ev = eigvals(M)
    return ev' .* projection(M)[:, 1:length(ev)]
end

# https://github.com/JuliaStats/MultivariateStats.jl/pull/162
function loadings(M::MultivariateStats.MDS, dim)
    l = loadings(M)
    return l[:, dim]
end

varexplained(M::MultivariateStats.MDS) = eigvals(M) ./ sum(eigvals(M))

mdsaxis(M::MultivariateStats.MDS, dim::Int) = "MDS$dim ($(round(varexplained(M)[dim] * 100, digits=2))%)"

varexpl(p::PERMANOVA.PSummary) = p.results[1, 3] * 100
pvalue(p::PERMANOVA.PSummary) = p.results[1, 5]

function permanovas(comm, metadatums)
    permdf = DataFrame()
    for md in metadatums
        com_md = get(comm, md)
        hasmd = findall(!ismissing, com_md)
        df = DataFrame(test = com_md[hasmd])
        disallowmissing!(df)
        size(df, 1) < 20 && @warn "Very small number of not missing $md"

        p = permanova(df, abundances(comm)[:, hasmd]', BrayCurtis, @formula(1 ~ test))
        push!(permdf, (; metadatum = md, varexpl=varexpl(p), pvalue=pvalue(p)))
    end
    return permdf

end

function plot_permanovas(comms, metadatums; commlabels=[], mdlabels=[], colormap=:blues, colorrange=(0,10))
    fig = Figure()
    ax = Axis(fig[1,1])
    hm = plot_permanovas!(ax, comms, metadatums; commlabels, mdlabels, colormap, colorrange)

    return fig, ax, hm
end

function plot_permanovas!(ax, comms, metadatums; commlabels=[], mdlabels=[], colormap=:blues, colorrange=(0,10))
    ps = [permanovas(comm, metadatums) for comm in comms]
    vmat = mapreduce(df-> df.varexpl, hcat, ps)
    pmat = mapreduce(df-> df.pvalue, hcat, ps)

    hm = heatmap!(ax, vmat'; colormap, colorrange)

    for ci in CartesianIndices(pmat)
        c = vmat[ci] < (colorrange[2] - colorrange[1]) / 2 ? :black : :lightgray
        text!(string(round(vmat[ci], digits=2), "%"); position=(ci[2],ci[1]), align=(:center, :center), color=c)
        p = pmat[ci]
        stars = p < 0.001 ? "***" : p < 0.01 ? "**" : p < 0.05 ? "*" : ""
        text!(stars; position=(ci[2],ci[1]), align=(:center, :bottom), color=c)
    end

    ax.xticks = (1:length(comms), isempty(commlabels) ? ["comm$i" for i in 1:length(comms)] : commlabels)
    ax.yticks = (1:length(metadatums), isempty(commlabels) ? string.(metadatums) : mdlabels)

    return hm
end

max_stars_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_issues_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis

max_forks_repo_name: Klepac-Ceraj-Lab/ResonanceAnalysis


text:
using LinearAlgebra

"""
    regularized_Picard(B, K, dpp_samples_ids, Fredholm_sample_ids, lambda, it_max, tol, use_inverse)

Runs regularized Picard iteration starting from an initial positive definite matrix B.

# Arguments
- `B::Array{Float64,2}`: an initial positive definite matrix B.
- `K::Array{Float64,2}`: full kernel matrix.
- `dpp_samples_ids:Array{Array{Int64,1},1}`: array of arrays containing indices of s DPP samples.
- `Fredholm_sample_ids:Array{Int64,1}`: indices of samples for approximating Fredholm determinant.
- `lambda:Float`: regularization parameter.
- `it_max::Integer`: largest number of iterations.
- `tol:Float`: relative objective variation desired for regularized Picard iteration.
- `use_inverse:Bool`: if true, inverse matrices are computed, otherwise only linear systems.

# Output
- `B::Array{Float64,2}`: approximate solution.
- `R::Array{Float64,2}`: Cholesky s.t. K = R'R.
- `obj`: array containing the objective values at each iteration.
- `i_stop`: number of iterations runned before stopping.

"""
function regularized_Picard(B::Array{Float64,2},K::Array{Float64,2}, dpp_samples_ids::Array{Array{Int64,1},1}, Fredholm_sample_ids::Array{Int64,1}, lambda::Float64, it_max::Int64 ,tol::Float64,use_inverse::Bool)

# number of dpp_samples_ids
nb_dpp_samples_ids = length(dpp_samples_ids);

# total number of points
m = size(K,1);

# Chol decomposition
R = cholesky(K).U;

# sampling matrix for uniformSample
nb_unif = length(Fredholm_sample_ids)

# initialization
obj = zeros(it_max,1);
i_stop = it_max;

# regularization
epsilon = 1e-10;

# iterations
for i in 1:it_max
    # construct  Delta
    BRDeltaRtB = zeros(m,m);
    RtBR = R'*B*R
    BR = B*R
    for l = 1:nb_dpp_samples_ids
        id = dpp_samples_ids[l];
        M = (RtBR + UniformScaling(epsilon))[id,id]
        BRid = BR[:,id]
        if use_inverse
            BRDeltaRtB +=  BRid * inv(M) * BRid';
        else
            BRDeltaRtB +=  BRid * (M\(BRid'));
        end
    end

    N = (RtBR + UniformScaling(nb_unif))[Fredholm_sample_ids,Fredholm_sample_ids]
    BRunif = BR[:,Fredholm_sample_ids]
    if use_inverse
        BRDeltaRtB = BRDeltaRtB/nb_dpp_samples_ids - BRunif*inv(N)*BRunif';
    else
        BRDeltaRtB = BRDeltaRtB/nb_dpp_samples_ids - BRunif*(N\(BRunif'));
    end

    BRDeltaRtB = 0.5*(BRDeltaRtB + BRDeltaRtB');

    # Picard iteration
    pB = B + BRDeltaRtB;

    # final expression
    B = (0.5/lambda)*(real(sqrt(I+4*lambda*pB))-I);
    B = 0.5*(B+B');

    # track the objective values
    obj_det,ob_reg = Picard_objective(B, dpp_samples_ids, Fredholm_sample_ids, R,lambda);
    obj[i] = obj_det + ob_reg;
    # printing in repl
    if i%100 == 0
        rel_variation = abs(obj[i]-obj[i-1])/abs(obj[i])
        print("---------------------------------------------------------------\n")
        print("$(i) / $it_max\n")
        print("relative objective variation $(rel_variation)\n")
        print("objective = $(obj[i]) \n")

    end
    # stopping criterion
    if i>1 && abs(obj[i]-obj[i-1])/abs(obj[i])< tol
        i_stop = i;
        print("---------------------------------------------------------------\n")
        print("Relative tolerance $(tol) attained after $(i) iterations.\n")
        print("Final objective= $(obj[i])\n")
        print("---------------------------------------------------------------\n")
        break
    end
    if i==it_max
        print("iteration has not yet converged.\n")
    end
end

return B, R, obj, i_stop

end

"""
    Picard_objective(B, dpp_samples_ids, Fredholm_samples_ids, R,lambda)

Computes Picard objective function.

# Arguments
- `B::Array{Float64,2}`: positive definite matrix B.
- `dpp_samples_ids:Array{Array{Int64,1},1}`: array of arrays containing indices of s DPP samples.
- `Fredholm_sample_ids:Array{Int64,1}`: indices of samples for approximating Fredholm determinant.
- `R`: Cholesky decomposition such that K = R'R.
- `lambda:Float`: regularization parameter.

# Outputs
- `ob_det:Float`: term of the objective including the determinants.
- `ob_reg:Float`: term of the objective including the regularization.

For the full objective do ob_det + ob_reg
"""
function Picard_objective(B, dpp_samples_ids, Fredholm_samples_ids, R,lambda)

    # number of dpp samples
    nb_dpp_samples_ids = length(dpp_samples_ids);

    # samples Fredholm
    nb_unif = length(Fredholm_samples_ids)

    PhiBPhi = R'*B*R;
    ob_det = 0
    for l = 1:nb_dpp_samples_ids
        id = dpp_samples_ids[l];
        ob_det -= logdet(PhiBPhi[id,id]);
        if ob_det==Inf
            error("singular determinant in objective\n")
        end
    end
    ob_det = ob_det/nb_dpp_samples_ids;
    ob_det += logdet(I + (1/nb_unif)*PhiBPhi[Fredholm_samples_ids,Fredholm_samples_ids]);
    ob_reg = lambda*tr(B);

    # full objective = ob_det + ob_reg
    return ob_det, ob_reg;
end


text:
#=
estep.jl

    Expectation step of the Expectation-Maximization (EM) algorithm to estimate
    the hyper parameters of a linear Gaussian State Space model.

@author: Quint Wiersma <q.wiersma@vu.nl>

@date: 2022/02/09
=#

"""
    estep!(smoother, filter, sys)

Perform expectation step of the joint log-likelihood of the linear Gaussian
State Space model as part of the EM algorithm.

#### Arguments
  - `sys::StateSpaceSystem` : state space system matrices

#### Returns
  - `filter::KalmanFilter`  : Kalman filter output
  - `smoother::Smoother`    : Kalman smoother output
"""
function estep!(smoother::Smoother, filter::KalmanFilter, sys::StateSpaceSystem)
	# Run filter
    kalman_filter!(filter, sys)
	# Run smoother
    kalman_smoother_cov!(smoother, filter, sys)

    return nothing
end