text:
using Convex,SCS

mutable struct SSPOR
    basis::Basis
    selected_sensors::Vector
    n_sensors::Int
    n_basis_modes::Int
    fit_::Function
    predict_::Function
    optimizer
end

struct Reconstructor
    fit::Function
    predict::Function
end

# function fit(model::SSPOR;quiet=true)
#     if model.fit_ === fit_lsq
#         model.fit_()
#     end
# end

# function predict(model::SSPOR,x::Vector;quiet=true)
#     if model.predict_ === predict_lsq
#         model.predict_(model.Ψ,x,model.n_basis_modes)
#     end

# end

function fit_lsq(Ψ,x,y,selected_sensors)
    y_sensed = y[selected_sensors];
    x_sensed = x[selected_sensors];
    Ψ_sensed = Ψ[:,selected_sensors]
    return Ψ_sensed,x_sensed,y_sensed
end

function predict_lsq(Ψ,x,y,n_basis_modes)
    X = Convex.Variable(size(x))
    prob = minimize(sumsquares(y-Ψ*X))
    solve!(prob,()-> SCS.Optimizer(verbose=false))
    prob.status
    return prob.optval
end


text:
include(joinpath(@__DIR__, "..", "src/dynamics", "hopper_2D", "visuals.jl"))
T = Float64
vis = Visualizer()
open(vis)

# simulation model
include(joinpath(module_dir(), "src/dynamics/hopper_2D/model.jl"))

model_sim = Hopper2D(Dimensions(nq, nu, nw, nc, nb),
			   1.2 * mb, ml, 1.2 * Jb, Jl,
			   μ_world, μ_joint, g,
			   BaseMethods(), DynamicsMethods(), ContactMethods(),
			   ResidualMethods(), ResidualMethods(),
			   SparseStructure(spzeros(0, 0), spzeros(0, 0)),
			   SVector{4}(zeros(4)),
			   environment_2D_flat())

expr_base = generate_base_expressions(model, M_analytical=true)
instantiate_base!(model_sim.base, expr_base)
expr_dyn = generate_dynamics_expressions(model_sim)
instantiate_dynamics!(model_sim.dyn, expr_dyn)
expr_res, rz_sp, rθ_sp = generate_residual_expressions(model_sim)
instantiate_residual!(model_sim.res, expr_res)
model_sim.spa.rz_sp = rz_sp
model_sim.spa.rθ_sp = rθ_sp

# get hopper model
model = get_model("hopper_2D")
nq = model.dim.q
nu = model.dim.u
nc = model.dim.c
nb = model.dim.b
nd = nq + nc + nb
nr = nq + nu + nc + nb + nd

# get trajectory
# ref_traj = get_trajectory("hopper_2D", "gait_in_place", load_type=:joint_traj)
ref_traj = get_trajectory("hopper_2D", "gait_forward", load_type=:joint_traj)
ref_traj_copy = deepcopy(ref_traj)

# time
H = ref_traj.H
h = ref_traj.h
N_sample = 5
H_mpc = 10
h_sim = h / N_sample
H_sim = 5000

# barrier parameter
κ_mpc = 1.0e-4

obj = TrackingObjective(H_mpc, model.dim,
    q = [Diagonal(1.0e-1 * [1,3,1,3])   for t = 1:H_mpc],
    u = [Diagonal(1.0e-0 * [1e-2, 1e0]) for t = 1:H_mpc],
    γ = [Diagonal(1.0e-100 * ones(model.dim.c)) for t = 1:H_mpc],
    b = [Diagonal(1.0e-100 * ones(model.dim.b)) for t = 1:H_mpc])

p = linearized_mpc_policy(ref_traj, model, obj,
    H_mpc = H_mpc,
    N_sample = N_sample,
    κ_mpc = κ_mpc,
    n_opts = NewtonOptions(
        r_tol = 3e-4,
        max_iter = 5),
    mpc_opts = LinearizedMPCOptions(
        altitude_update = true,
        altitude_impact_threshold = 0.5,
        altitude_verbose = true,
        )
    )

q1_ref = copy(ref_traj.q[2])
q0_ref = copy(ref_traj.q[1])
q1_sim = SVector{model.dim.q}(q1_ref)
q0_sim = SVector{model.dim.q}(copy(q1_sim - (q1_ref - q0_ref) / N_sample))
@assert norm((q1_sim - q0_sim) / h_sim - (q1_ref - q0_ref) / h) < 1.0e-8

sim = ContactImplicitMPC.simulator(model_sim, q0_sim, q1_sim, h_sim, H_sim,
    p = p,
	uL = [-0.25, -0.75], uU = [0.25, 0.75],
    ip_opts = ContactImplicitMPC.InteriorPointOptions(
        r_tol = 1.0e-8,
        κ_init = 1.0e-8,
        κ_tol = 2.0e-8),
    sim_opts = ContactImplicitMPC.SimulatorOptions(warmstart = true))

@time status = ContactImplicitMPC.simulate!(sim)

(sim.traj.q[end][1] - sim.traj.q[1][1]) / (h_sim * H_sim)
anim = visualize_robot!(vis, sim.model, sim.traj, sample=20)

plot(hcat(sim.traj.u...)[1:2, 1:10:end]')

ext: jl


text:
const rank  = split("A23456789TJQK", "")
const suit  = split("♣♦♥♠", "")
const deck  = Vector{String}()

const mslcg = [0]
rng() = (mslcg[1] = ((mslcg[1] * 214013 + 2531011) & 0x7fffffff)) >> 16

initdeck() = for r in rank, s in suit push!(deck, "$r$s") end

function deal(num = rand(UInt,1)[1] % 32000 + 1)
    initdeck()
    mslcg[1] = num
    println("\nGame # ", num)
    while length(deck) > 0
        choice = rng() % length(deck) + 1
        deck[choice], deck[end] = deck[end], deck[choice]
        print(" ", pop!(deck), length(deck) % 8 == 4 ? "\n" : "")
    end
end

deal(1)
deal(617)
deal()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
   setbits(x::Unsigned, y::Unsigned, mask::Unsigned)

If you have AVX512, setbits of vector-arguments will select bits according to mask `m`, selecting from `y` if 0 and from `x` if `1`.
For scalar arguments, or vector arguments without AVX512, `setbits` requires the additional restrictions on `y` that all bits for
which `m` is 1, `y` must be 0.
That is for scalar arguments or vector arguments without AVX512, it requires the restriction that
((y ⊻ m) & m) == m
"""
@inline setbits(x, y, m) = (x & m) | y
@inline function setbits(x::Vec{W,U}, y, m, ::True) where {W,U<:Union{UInt32,UInt64}}
    VectorizationBase.vpternlog(vbroadcast(Val{W}(), y), x, vbroadcast(Val{W}(), m), Val{216}())
end
@inline function setbits(x::VecUnroll{N,W,U}, y, m, ::True) where {N,W,U<:Union{UInt32,UInt64}}
    VectorizationBase.VecUnroll(VectorizationBase.fmap(setbits, data(x), y, m))
end
@inline setbits(x, y, m, ::False) = ((x & m) | y)
@inline setbits(x::VectorizationBase.AbstractSIMD, y, m) = setbits(x, y, m, VectorizationBase.has_feature(Val(:x86_64_avx512f)))

"""
The masks mask off bits so that a set of bits will be within the range 1.0 and prevfloat(2.0).
They do this via first using the bitwise-and to set all the exponent bits and the sign-bit to 0,
and then using the bitwise-or to set all but the first exponent bit to 1.
Neither operation touches the fraction-bits.

For example, you can see the binary representation of a double precision number on wikipedia:
https://en.wikipedia.org/wiki/Double-precision_floating-point_format
and then
julia> bitstring(0x000fffffffffffff)
"0000000000001111111111111111111111111111111111111111111111111111"

julia> bitstring(0x3ff0000000000000)
"0011111111110000000000000000000000000000000000000000000000000000"

so that `(x & 0x000fffffffffffff) | 0x3ff0000000000000` has the desired effect of setting the sign
bit to 0, and the exponent bit to be within the proper range.
Note that the second string has 10 "1"s, so that it is 2^10.
For double precision, the exponential contribution is 2^(e - 1023). With e = 2^10,
we have 2^(2^10 - 1023) = 2^(1024 - 1023) = 2^1 = 2.
Thus, the fraction specifies where in the range between 1.0 and prevfloat(2.0).
If the fraction bits are completely random, the distribution will be uniform between 1.0 and prevfloat(2.0).

The reason we target that range (1.0 to prevfloat(2.0)) is to get a nice, uniform distribution.
We can get uniform(0, 1] via 2 - u, or a uniform[0,1) via u - 1.
Choosing other fractional intervals, eg [0.5, 1.0) or [2.0, 4.0) would take more work to translate.
Alternatively, if we tried to span multiple fractional intervals, suddenly trying to get a fairly
uniform distribution would get complicated.
"""
@inline floatbitmask(x, ::Type{Float64}) = reinterpret(Float64, setbits(reinterpret(UInt64, x), 0x3ff0000000000000, 0x000fffffffffffff))
@inline floatbitmask(x, ::Type{Float32}) = reinterpret(Float32, setbits(reinterpret(UInt32, x), 0x3f800000, 0x007fffff))


text:
using Test
import LinearAlgebra: dot
using RationalSimplex

@testset "RationalSimplex" begin

@testset "Test standard form" begin
      # Test case 1
      # min -10x - 12y - 12z
      # st     x +  2y +  2z <= 20
      #       2x +   y +  2z <= 20
      #       2x +  2y +   z <= 20
      #       x, y, z >= 0
      @testset "Test case 1 (optimal)" begin
            c = [-10//1, -12//1, -12//1, 0//1, 0//1, 0//1]
            A = [  1//1    2//1    2//1  1//1  0//1  0//1;
                   2//1    1//1    2//1  0//1  1//1  0//1;
                   2//1    2//1    1//1  0//1  0//1  1//1]
            b = [20//1, 20//1, 20//1]
            status, x = simplex(c, A, b)
            @test status == :Optimal
            @test x[1] == 4//1
            @test x[2] == 4//1
            @test x[3] == 4//1
      end

      # Test case 2
      # min -90x - y
      # st    2x + y <= 40
      #      20x + y <= 100
      #       2x     <= 3
      #        x,  y >= 0
      @testset "Test case 2 (optimal)" begin
            c = [-90//1, -1//1, 0//1, 0//1, 0//1]
            A = [  2//1   1//1  1//1  0//1  0//1;
                  20//1   1//1  0//1  1//1  0//1;
                   2//1   0//1  0//1  0//1  1//1]
            b = [ 40//1, 100//1, 3//1]
            status, x = simplex(c, A, b)
            @test status == :Optimal
            @test x[1] == 3//2
            @test x[2] == 37//1
            @test x[3] == 0//1
            @test x[4] == 33//1
            @test x[5] == 0//1
      end

      # Test case 3
      # min x
      # st  x >= 3
      #     x <= 2
      @testset "Test case 3 (infeasible)" begin
            c = [1//1, 0//1, 0//1]
            A = [1//1 -1//1  0//1;
                 1//1  0//1  1//1]
            b = [3//1, 2//1]
            status, x = simplex(c, A, b)
            @test status == :Infeasible
      end

      # Test case 4
      # min -x
      # st  x >= 2
      @testset "Test Case 4 (unbounded)" begin
            c = [-1//1, 0//1]
            A = [ 1//1 -1//1;]
            b = [ 2//1]
            status, x = simplex(c, A, b)
            @test status == :Unbounded
      end

      # Check simplex detects nonneg rhs
      @testset "Check for nonnegative RHS" begin
            c = [-1//1, 0//1]
            A = [ 1//1 -1//1;]
            b = [-2//1]
            @test_throws AssertionError simplex(c, A, b)
      end
end

@testset "Test general form" begin
      # Test case 1
      # max  10x + 12y + 12z
      # st     x +  2y +  2z <= 20
      #       2x +   y +  2z <= 20
      #       2x +  2y +   z <= 20
      #       x, y, z >= 0
      @testset "Test case 1 (optimal)" begin
            c = [ 10//1,  12//1,  12//1]
            A = [  1//1    2//1    2//1;
                   2//1    1//1    2//1;
                   2//1    2//1    1//1]
            b = [20//1, 20//1, 20//1]
            status, x = simplex(c, :Max, A, b, ['<','<','<'])
            @test status == :Optimal
            @test x[1] == 4//1
            @test x[2] == 4//1
            @test x[3] == 4//1
      end

      # Test case 2
      # max  90x + y
      # st    2x + y <= 40
      #      20x + y <= 100
      #       2x     <= 3
      #        x,  y >= 0
      @testset "Test case 2 (optimal)" begin
            c = [ 90//1,  1//1]
            A = [  2//1   1//1;
                  20//1   1//1;
                   2//1   0//1]
            b = [ 40//1, 100//1, 3//1]
            status, x = simplex(c, :Max, A, b, ['<','<','<'])
            @test status == :Optimal
            @test x[1] == 3//2
            @test x[2] == 37//1
      end

      # Test case 3
      # min x
      # st  x >= 3
      #     x <= 2
      @testset "Test case 3 (infeasible)" begin
            c = [1//1, 0//1]
            A = [1//1 0//1;
                 1//1 0//1]
            b = [3//1, 2//1]
            status, x = simplex(c, :Min, A, b, ['>','<'])
            @test status == :Infeasible
      end

      # Test case 4
      # min x
      # st   x >=  3
      #     -x >= -2
      @testset "Test case 4 (infeasible)" begin
            c = [1//1, 0//1]
            A = [ 1//1 0//1;
                 -1//1 0//1]
            b = [3//1, -2//1]
            status, x = simplex(c, :Min, A, b, ['>','>'])
            @test status == :Infeasible
      end
end

@testset "Kantorovich distances" begin
      # Problem description from
      # http://stla.github.io/stlapblog/posts/KantorovichWithJulia.html
      c = [0//1, 1//1, 1//1, 1//1, 0//1, 1//1, 1//1, 1//1, 0//1]
      b = [1//7, 2//7, 4//7, 1//4, 1//4, 1//2]
      M = [1//1 1//1 1//1 0//1 0//1 0//1 0//1 0//1 0//1;
           0//1 0//1 0//1 1//1 1//1 1//1 0//1 0//1 0//1;
           0//1 0//1 0//1 0//1 0//1 0//1 1//1 1//1 1//1;
           1//1 0//1 0//1 1//1 0//1 0//1 1//1 0//1 0//1;
           0//1 1//1 0//1 0//1 1//1 0//1 0//1 1//1 0//1;
           0//1 0//1 1//1 0//1 0//1 1//1 0//1 0//1 1//1]

      status, x = simplex(c, :Min, M, b, ['=','=','=','=','=','='])
      @test dot(c,x) == 3//28
end

end  # RationalSimplex testset.


text:
# This file is a part of MGVI.jl, licensed under the MIT License (MIT).

struct PDLinMapWithChol{
    T,
    PMap<:LinearMaps.LinearMap{T},
    CMap<:LinearMaps.LinearMap{T}
} <: LinearMaps.LinearMap{T}
    _linmap::PMap
    _chol_L::CMap
end

function PDLinMapWithChol(A::AbstractMatrix, chol_A::AbstractMatrix)
    chol_L = LinearMap(
        chol_A,
        issymmetric = false, ishermitian = false, isposdef = false
    )

    wrapped_A = LinearMap(
        A,
        issymmetric = true, ishermitian = true, isposdef = true
    )

    PDLinMapWithChol(wrapped_A, chol_L)
end

PDLinMapWithChol(A::AbstractMatrix) = PDLinMapWithChol(A, cholesky_L(A))

without_chol(A::PDLinMapWithChol) = A._linmap
cholesky_L(A::PDLinMapWithChol) = A._chol_L

Base.Matrix(A::PDLinMapWithChol) = Matrix(without_chol(A))

Base.convert(::Type{AbstractMatrix}, A::PDLinMapWithChol) =
    convert(AbstractMatrix, without_chol(A))
Base.convert(::Type{Matrix}, A::PDLinMapWithChol) =
    convert(Matrix, without_chol(A))
Base.convert(::Type{SparseMatrixCSC}, A::PDLinMapWithChol) =
    convert(SparseMatrixCSC, without_chol(A))
SparseArrays.sparse(A::PDLinMapWithChol) =
    sparse(without_chol(A))

Base.size(A::PDLinMapWithChol) = size(without_chol(A))

LinearAlgebra.issymmetric(A::PDLinMapWithChol) = true
LinearAlgebra.ishermitian(A::PDLinMapWithChol) = true
LinearAlgebra.isposdef(A::PDLinMapWithChol) = true

LinearAlgebra.adjoint(A::PDLinMapWithChol) = A

LinearAlgebra.transpose(A::PDLinMapWithChol{<:Real}) = A
LinearAlgebra.transpose(A::PDLinMapWithChol) = PDLinMapWithChol(transpose(without_chol(A)), cholesky_L(A))

LinearMaps.MulStyle(A::PDLinMapWithChol) = LinearMaps.MulStyle(without_chol(A))

Base.@propagate_inbounds LinearAlgebra.mul!(
    y::AbstractVector, A::PDLinMapWithChol, x::AbstractVector, α::Number, β::Number
) = mul!(y, without_chol(A), x, α, β)

Base.@propagate_inbounds LinearAlgebra.mul!(
    y::AbstractMatrix, A::PDLinMapWithChol, x::AbstractMatrix, α::Number, β::Number
) = mul!(y, without_chol(A), x, α, β)

Base.:(*)(A::LinearMap, B::PDLinMapWithChol) = A * without_chol(B)
Base.:(*)(A::PDLinMapWithChol, B::LinearMap) = without_chol(A) * B
Base.:(*)(A::PDLinMapWithChol, B::PDLinMapWithChol) = without_chol(A) * without_chol(B)

function _blockdiag(As::NTuple{N,PDLinMapWithChol}) where N
    PDLinMapWithChol(blockdiag(map(without_chol, As)...), blockdiag(map(cholesky_L, As)...))
end

const DiagLinearMap{T} = LinearMaps.WrappedMap{T,<:Diagonal}
const DiagPDLinMapWithChol = PDLinMapWithChol{T,<:DiagLinearMap,<:DiagLinearMap} where T

function _blockdiag(As::NTuple{N, DiagPDLinMapWithChol}) where N
    PDLinMapWithChol(
        LinearMap(Diagonal(vcat(map(A -> get_diagonal(without_chol(A).lmap), As)...)), issymmetric = true, ishermitian = true, isposdef = true),
        LinearMap(Diagonal(vcat(map(A -> get_diagonal(cholesky_L(A).lmap), As)...)), issymmetric = true, ishermitian = true, isposdef = true)
    )
end

function _blockdiag(A::AbstractVector{<:DiagPDLinMapWithChol})
    d = _flatten_vec_of_vec(map(x -> get_diagonal(without_chol(x).lmap), A))
    PDLinMapWithChol(Diagonal(d))
end

blockdiag(As::PDLinMapWithChol...) = _blockdiag(As)


text:
# eq 4.9
# Gilbert, J.C., Lemaréchal, C. Some numerical experiments with variable-storage quasi-Newton algorithms.
# Mathematical Programming 45, 407–435 (1989). https://doi.org/10.1007/BF01589113
function gilbert_init(α, s, y)
    a = dot(y, Diagonal(α), y)
    b = dot(y, s)
    c = dot(s, inv(Diagonal(α)), s)
    return @. b / (a / α + y^2 - (a / c) * (s / α)^2)
end

"""
    lbfgs_inverse_hessians(θs, ∇logpθs; Hinit=gilbert_init, history_length=5, ϵ=1e-12) -> Vector{WoodburyPDMat}

From an L-BFGS trajectory and gradients, compute the inverse Hessian approximations at each point.

Given positions `θs` with gradients `∇logpθs`, construct LBFGS inverse Hessian
approximations with the provided `history_length`.
"""
function lbfgs_inverse_hessians(θs, ∇logpθs; Hinit=gilbert_init, history_length=5, ϵ=1e-12)
    L = length(θs) - 1
    θ = θs[1]
    ∇logpθ = ∇logpθs[1]
    n = length(θ)

    # allocate caches/containers
    history_ind = 0 # index of last set history entry
    history_length_effective = 0 # length of history so far
    s = similar(θ) # cache for BFGS update, i.e. sₗ = θₗ₊₁ - θₗ = -λ Hₗ ∇logpθₗ
    y = similar(∇logpθ) # cache for yₗ = ∇logpθₗ₊₁ - ∇logpθₗ = Hₗ₊₁ \ s₁ (secant equation)
    S = similar(s, n, min(history_length, L)) # history of s
    Y = similar(y, n, min(history_length, L)) # history of y
    α = fill!(similar(θ), true) # diag(H₀)
    H = lbfgs_inverse_hessian(Diagonal(α), S, Y, history_ind, history_length_effective) # H₀ = I
    Hs = [H] # trace of H

    for l in 1:L
        θlp1, ∇logpθlp1 = θs[l + 1], ∇logpθs[l + 1]
        (any(isnan, θlp1) || any(isnan, ∇logpθlp1)) && break
        s .= θlp1 .- θ
        y .= ∇logpθ .- ∇logpθlp1
        if dot(y, s) > ϵ * sum(abs2, y)  # curvature is positive, safe to update inverse Hessian
            # add s and y to history
            history_ind = mod1(history_ind + 1, history_length)
            history_length_effective = max(history_ind, history_length_effective)
            S[1:n, history_ind] .= s
            Y[1:n, history_ind] .= y

            # initial diagonal estimate of H
            α = Hinit(α, s, y)
        else
            @debug "Skipping inverse Hessian update from iteration $l to avoid negative curvature."
        end

        θ, ∇logpθ = θlp1, ∇logpθlp1
        H = lbfgs_inverse_hessian(Diagonal(α), S, Y, history_ind, history_length_effective)
        push!(Hs, H)
    end

    return Hs
end

"""
    lbfgs_inverse_hessian(H₀, S₀, Y₀, history_ind, history_length) -> WoodburyPDMat

Compute approximate inverse Hessian initialized from `H₀` from history stored in `S₀` and `Y₀`.

`history_ind` indicates the column in `S₀` and `Y₀` that was most recently added to the
history, while `history_length` indicates the number of first columns in `S₀` and `Y₀`
currently being used for storing history.
`S = S₀[:, history_ind+1:history_length; 1:history_ind]` reorders the columns of `₀` so that the
oldest is first and newest is last.

From Theorem 2.2 of [^Byrd1994], the expression for the inverse Hessian ``H`` is

```math
\\begin{align}
B &= \\begin{pmatrix}H_0 Y & S\\end{pmatrix}\\\\
R &= \\operatorname{triu}(S^\\mathrm{T} Y)\\\\
E &= I \\circ R\\\\
D &= \\begin{pmatrix}
    0 & -R^{-1}\\\\
    -R^{-\\mathrm{T}} & R^\\mathrm{-T} (E + Y^\\mathrm{T} H_0 Y ) R^\\mathrm{-1}\\\\
H &= H_0 + B D B^\\mathrm{T}
\\end{pmatrix}
\\end{align}
```

[^Byrd1994]: Byrd, R.H., Nocedal, J. & Schnabel, R.B.
             Representations of quasi-Newton matrices and their use in limited memory methods.
             Mathematical Programming 63, 129–156 (1994).
             doi: [10.1007/BF01582063](https://doi.org/10.1007/BF01582063)
"""
function lbfgs_inverse_hessian(H₀::Diagonal, S0, Y0, history_ind, history_length)
    J = history_length
    α = H₀.diag
    B = similar(α, size(α, 1), 2J)
    D = fill!(similar(α, 2J, 2J), false)
    iszero(J) && return WoodburyPDMat(H₀, B, D)

    hist_inds = [(history_ind + 1):history_length; 1:history_ind]
    @views begin
        S = S0[:, hist_inds]
        Y = Y0[:, hist_inds]
        B₁ = B[:, 1:J]
        B₂ = B[:, (J + 1):(2J)]
        D₁₁ = D[1:J, 1:J]
        D₁₂ = D[1:J, (J + 1):(2J)]
        D₂₁ = D[(J + 1):(2J), 1:J]
        D₂₂ = D[(J + 1):(2J), (J + 1):(2J)]
    end

    mul!(B₁, Diagonal(α), Y)
    copyto!(B₂, S)
    mul!(D₂₂, S', Y)
    triu!(D₂₂)
    R = UpperTriangular(D₂₂)
    nRinv = UpperTriangular(D₁₂)
    copyto!(nRinv, -I)
    ldiv!(R, nRinv)
    nRinv′ = LowerTriangular(copyto!(D₂₁, nRinv'))
    tril!(D₂₂) # eliminate all but diagonal
    mul!(D₂₂, Y', B₁, true, true)
    LinearAlgebra.copytri!(D₂₂, 'U', false, false)
    rmul!(D₂₂, nRinv)
    lmul!(nRinv′, D₂₂)

    return WoodburyPDMat(H₀, B, D)
end


text:
#!/usr/bin/env julia

using WaterMeanMeanForce

using ArgParse

s = ArgParseSettings()
s.autofix_names = true

@add_arg_table! s begin
    "-T"
        metavar = "T"
        help = "comma-separated list of temperatures (K)"
        required = true
end

c = parse_args(ARGS, s, as_symbols=true)

for T_str in split(c[:T], ",")
    T = parse(Float64, T_str)
    beta, slope = harmonic_exact_slope(2, T)

    println("$(T) $(beta) $(slope)")
end


text:
# Enumeration.

dist = @pp enum() do
    local
    a = flip(0.5),
    b = flip(0.5),
    c = flip(0.5)
    factor(a||b ? 0 : -1)
    a + b + c
end

@test_approx_eq(exp(score(dist,0)), 0.054615886286517)
@test_approx_eq(exp(score(dist,1)), 0.351538628762172)
@test_approx_eq(exp(score(dist,2)), 0.445384113713482)
@test_approx_eq(exp(score(dist,3)), 0.148461371237827)

@test_throws ErrorException @pp enum() do
    factor(-Inf)
    true
end

# Ensure the context is restored when an exception occurs during
# enumeration.
try @pp enum(()->foo())
catch e
    if isa(e, UndefVarError)
        @test Stochy.ctx == Stochy.Prior()
    else
        rethrow(e)
    end
end

# PMCMC.

dist = @pp pmcmc(5,5) do
    local x = flip()
    factor(x ? 0 : -1)
    x
end

@test length(support(dist)) == 2
@test true in support(dist)
@test false in support(dist)

try @pp pmcmc(()->foo(),1,1)
catch e
    if isa(e, UndefVarError)
        @test Stochy.ctx == Stochy.Prior()
    else
        rethrow(e)
    end
end

# MH.

dist = @pp mh(5) do
    local x = flip(0.5)
    factor(x ? 0 : -1)
    x
end

@test issubset(support(dist), [true, false])

try @pp mh(()->foo(),1)
catch e
    if isa(e, UndefVarError)
        @test Stochy.ctx == Stochy.Prior()
    else
        rethrow(e)
    end
end


text:
"""
    SubDataFrame{<:AbstractVector{Int}} <: AbstractDataFrame

A view of row subsets of an AbstractDataFrame

A `SubDataFrame` is meant to be constructed with `view`.  A
SubDataFrame is used frequently in split/apply sorts of operations.

```julia
view(d::AbstractDataFrame, rows)
```

### Arguments

* `d` : an AbstractDataFrame
* `rows` : any indexing type for rows, typically an Int,
  AbstractVector{Int}, AbstractVector{Bool}, or a Range

### Notes

A `SubDataFrame` is an AbstractDataFrame, so expect that most
DataFrame functions should work. Such methods include `describe`,
`dump`, `nrow`, `size`, `by`, `stack`, and `join`. Indexing is just
like a DataFrame; copies are returned.

To subset along columns, use standard column indexing as that creates
a view to the columns by default. To subset along rows and columns,
use column-based indexing with `view`.

### Examples

```julia
df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
               b = repeat([2, 1], outer=[4]),
               c = randn(8))
sdf1 = view(df, 1:6)
sdf2 = view(df, df[:a] .> 1)
sdf3 = view(df[[1,3]], df[:a] .> 1)  # row and column subsetting
sdf4 = groupby(df, :a)[1]  # indexing a GroupedDataFrame returns a SubDataFrame
sdf5 = view(sdf1, 1:3)
sdf1[:,[:a,:b]]
```

"""
struct SubDataFrame{T<:AbstractVector{Int}} <: AbstractDataFrame
    parent::DataFrame
    rows::T # maps from subdf row indexes to parent row indexes

    function SubDataFrame{T}(parent::DataFrame, rows::T) where {T <: AbstractVector{Int}}
        if length(rows) > 0
            rmin, rmax = extrema(rows)
            if rmin < 1 || rmax > size(parent, 1)
                throw(BoundsError())
            end
        end
        new(parent, rows)
    end
end

function SubDataFrame(parent::DataFrame, rows::T) where {T <: AbstractVector{Int}}
    return SubDataFrame{T}(parent, rows)
end

function SubDataFrame(parent::DataFrame, rows::Colon)
    return SubDataFrame(parent, 1:nrow(parent))
end

function SubDataFrame(parent::DataFrame, row::Integer)
    throw(ArgumentError("invalid row index: $row of type $(typeof(row))"))
end

function SubDataFrame(parent::DataFrame, rows::AbstractVector{<:Integer})
    if any(x -> x isa Bool, rows)
        throw(ArgumentError("invalid row index of type `Bool`"))
    end
    return SubDataFrame(parent, convert(Vector{Int}, rows))
end

function SubDataFrame(parent::DataFrame, rows::AbstractVector{Bool})
    if length(rows) != nrow(parent)
        throw(ArgumentError("invalid length of `AbstractVector{Bool}` row index" *
                            " (got $(length(rows)), expected $(nrow(parent)))"))
    end
    return SubDataFrame(parent, findall(rows))
end

function SubDataFrame(parent::DataFrame, rows::AbstractVector)
    if !all(x -> (x isa Integer) && !(x isa Bool), rows)
        throw(ArgumentError("only `Integer` indices are accepted in `rows`"))
    end
    return SubDataFrame(parent, convert(Vector{Int}, rows))
end

function SubDataFrame(sdf::SubDataFrame, rowinds)
    return SubDataFrame(parent(sdf), rows(sdf)[rowinds])
end

SubDataFrame(sdf::SubDataFrame, rowinds::Colon) = sdf

Base.parent(sdf::SubDataFrame) = getfield(sdf, :parent)
Base.parentindices(sdf::SubDataFrame) = (rows(sdf), axes(parent(sdf), 2))
rows(sdf::SubDataFrame) = getfield(sdf, :rows)

Base.view(adf::AbstractDataFrame, colinds) = view(adf, :, colinds)
Base.view(adf::AbstractDataFrame, rowinds, colind::ColumnIndex) =
    view(adf[colind], rowinds)
Base.view(adf::AbstractDataFrame, rowinds, colind::Bool) =
    throw(ArgumentError("invalid column index $colind of type `Bool`"))
Base.view(adf::AbstractDataFrame, rowinds, colinds) =
    SubDataFrame(adf[colinds], rowinds)
Base.view(adf::AbstractDataFrame, rowinds, ::Colon) = SubDataFrame(adf, rowinds)

##############################################################################
##
## AbstractDataFrame interface
##
##############################################################################

index(sdf::SubDataFrame) = index(parent(sdf))

# TODO: Remove these
nrow(sdf::SubDataFrame) = ncol(sdf) > 0 ? length(rows(sdf))::Int : 0
ncol(sdf::SubDataFrame) = length(index(sdf))

Base.getindex(sdf::SubDataFrame, colind::ColumnIndex) =
    view(parent(sdf)[colind], rows(sdf))
Base.getindex(sdf::SubDataFrame, colinds::AbstractVector) =
    SubDataFrame(parent(sdf)[colinds], rows(sdf))
Base.getindex(sdf::SubDataFrame, ::Colon) = sdf
Base.getindex(sdf::SubDataFrame, rowind::Integer, colind::ColumnIndex) =
    parent(sdf)[rows(sdf)[rowind], colind]
Base.getindex(sdf::SubDataFrame, rowinds::AbstractVector, colind::ColumnIndex) =
    parent(sdf)[rows(sdf)[rowinds], colind]
Base.getindex(sdf::SubDataFrame, ::Colon, colind::ColumnIndex) =
    parent(sdf)[rows(sdf), colind]
Base.getindex(sdf::SubDataFrame, ::Colon, colinds::AbstractVector) =
    parent(sdf)[rows(sdf), colinds]
Base.getindex(sdf::SubDataFrame, rowinds::AbstractVector, colinds::AbstractVector) =
    parent(sdf)[rows(sdf)[rowinds], colinds]
Base.getindex(sdf::SubDataFrame, rowinds::AbstractVector, ::Colon) =
    parent(sdf)[rows(sdf)[rowinds], :]
Base.getindex(sdf::SubDataFrame, ::Colon, ::Colon) =
    parent(sdf)[rows(sdf), :]

function Base.setindex!(sdf::SubDataFrame, val::Any, colinds::Any)
    parent(sdf)[rows(sdf), colinds] = val
    return sdf
end

function Base.setindex!(sdf::SubDataFrame, val::Any, rowinds::Any, colinds::Any)
    parent(sdf)[rows(sdf)[rowinds], colinds] = val
    return sdf
end

##############################################################################
##
## Miscellaneous
##
##############################################################################

Base.copy(sdf::SubDataFrame) = parent(sdf)[rows(sdf), :]

without(sdf::SubDataFrame, c) = view(without(parent(sdf), c), rows(sdf), :)
# Resolve a method ambiguity
without(sdf::SubDataFrame, c::Vector{<:Integer}) = view(without(parent(sdf), c), rows(sdf), :)

deleterows!(df::SubDataFrame, ind) =
    throw(ArgumentError("SubDataFrame does not support deleting rows"))


text:
abstract AbstractStateSpaceModel{T <: Real}

immutable StateSpaceModel{T} <: AbstractStateSpaceModel{T}

    # Process transition matrix, control matrix, and noise covariance
    A::Function
    B::Function
    V::Function

    # Observation matrix, feed-forward matrix, and noise covariance
    C::Function
    D::Function
    W::Function

    # Inital state and error covariance conditions
    x1::Vector{T}
    P1::Matrix{T}

    # Model dimensions
    nx::Int
    ny::Int
    nu::Int

    function StateSpaceModel{T}(A::Function, B::Function, V::Function,
                                C::Function, D::Function, W::Function,
                                x1::Vector{T}, P1::Matrix{T})

        ispossemidef(x::Matrix) = issym(x) && (eigmin(x) >= 0)
        @assert ispossemidef(V(1))
        @assert ispossemidef(W(1))
        @assert ispossemidef(P1)

        nx, ny, nu = confirm_matrix_sizes(A(1), B(1), V(1), C(1), D(1), W(1), x1, P1)
        new(A, B, V, C, D, W, x1, P1, nx, ny, nu)
    end
end

# Time-dependent constructor
StateSpaceModel{T}(A::Function, V::Function, C::Function, W::Function,
                          x1::Vector{T}, P1::Matrix{T};
                          B::Function=_->zeros(size(V(1), 1), 1),
                          D::Function=_->zeros(size(W(1), 1), 1)) =
	  StateSpaceModel{T}(A, B, V, C, D, W, x1, P1)

# Time-independent constructor
StateSpaceModel{T}(A::Matrix{T}, V::Matrix{T}, C::Matrix{T}, W::Matrix{T},
                          x1::Vector{T}, P1::Matrix{T};
                          B::Matrix{T}=zeros(size(A, 1), 1),
                          D::Matrix{T}=zeros(size(C, 1), 1)) =
	  StateSpaceModel{T}(_->A, _->B, _->V, _->C, _->D, _->W, x1, P1)

function show{T}(io::IO, mod::StateSpaceModel{T})
    dx, dy = mod.nx, mod.ny
    println("StateSpaceModel{$T}, $dx-D process x $dy-D observations")
    println("Process evolution matrix A:")
    show(mod.A(1))
    println("\n\nControl input matrix B:")
    show(mod.B(1))
    println("\n\nProcess error covariance V:")
    show(mod.V(1))
    println("\n\nObservation matrix C:")
    show(mod.C(1))
    println("\n\nFeed-forward matrix D:")
    show(mod.D(1))
    println("\n\nObseration error covariance W:")
    show(mod.W(1))
end

# Time-independent parametrized matrix type
# TODO: A more general ParametrizedArray (allowing vectors) could be preferable
immutable ParametrizedMatrix{T<:Real}

    f::Vector{T}
    D::Matrix{T}
    np::Int
    dims::Tuple{Int, Int}

    function ParametrizedMatrix(f::Vector{T}, D::Matrix{T}, dims::Tuple{Int, Int})
        @assert length(f) == size(D, 1)
        @assert length(f) == dims[1] * dims[2]
        new(f, D, size(D, 2), dims)
    end #ParametrizedMatrix

end #ParametrizedMatrix

ParametrizedMatrix{T}(f::Vector{T}, D::Matrix{T}, dims::Tuple{Int, Int}) = ParametrizedMatrix{T}(f, D, dims)

function show{T}(io::IO, cpm::ParametrizedMatrix{T})

    function combinestrelems(a, b)
        if (a != "") & (b != "")
          a * " + " * b
        elseif (a == "") & (b == "")
          "0"
        else
          a * b
        end #if
    end #combinestrelems

    nx, ny = cpm.dims
    println(nx, "x", ny, " ParametrizedMatrix{$T}")
    conststring = map(x -> x == 0 ? "" : string(x), cpm.f)
    paramstring = fill("", length(cpm.f))
    paramelems = fill("", cpm.np)

    for i = 1:size(cpm.D,1)
        for j = 1:cpm.np
            if  cpm.D[i,j] == 0
                paramelems[j] = ""
            elseif cpm.D[i,j] == 1
                paramelems[j] = string(Char(96+j))
            else
                paramelems[j] = string(cpm.D[i,j],Char(96+j))
            end #if
        end #for
        paramstring[i] = reduce(combinestrelems, paramelems)
    end #for

    finalstrings = map(combinestrelems, conststring, paramstring)

    showcompact(reshape(finalstrings, cpm.dims))
end

Base.length(cpm::ParametrizedMatrix) = cpm.dims[1] * cpm.dims[2]
Base.size(cpm::ParametrizedMatrix) = cpm.dims
Base.size(cpm::ParametrizedMatrix, dim::Int) = cpm.dims[dim]
Base.call(cpm::ParametrizedMatrix, params::Vector) = all(cpm.D .== 0) ?
        reshape(cpm.f, cpm.dims) : reshape(cpm.f + cpm.D * params, cpm.dims)

parametrize_full(X::Union{Vector, Matrix}) =
        ParametrizedMatrix(zeros(length(X)), eye(length(X)), size(X'')), collect(X)

function parametrize_diag(x::Vector)
    n = length(x)
    D = zeros(n^3)
    D[[1 + (i-1)*(n^2 + n + 1) for i in 1:n]] = 1
    D = reshape(D, n^2, n)
    f = zeros(n^2)
    return ParametrizedMatrix(f, D, (n,n)), x
end #parametrize_diag

parametrize_none{T}(X::Union{Vector{T}, Matrix{T}}) =
        ParametrizedMatrix(collect(X), zeros(length(X))'', size(X'')), T[]

# Time-independent parametrized state space model
immutable ParametrizedSSM{T} <: AbstractStateSpaceModel{T}

    # Transition equation and noise covariance

    A1::Function
    A2::ParametrizedMatrix{T}
    A3::Function

    B1::Function
    B2::ParametrizedMatrix{T}

    G::Function
    Q::ParametrizedMatrix{T}

    # Observation equation and noise covariance

    C1::Function
    C2::ParametrizedMatrix{T}
    C3::Function

    D1::Function
    D2::ParametrizedMatrix{T}

    H::Function
    R::ParametrizedMatrix{T}

    # Initial state and error covariance
    x1::ParametrizedMatrix{T}
    P1::ParametrizedMatrix{T} #TODO ... Eliminate?

    nx::Int
    ny::Int
    nu::Int
    nq::Int
    nr::Int

    function ParametrizedSSM{T}(A1::Function, A2::ParametrizedMatrix{T}, A3::Function,
                                B1::Function, B2::ParametrizedMatrix{T},
                                G::Function, Q::ParametrizedMatrix{T},
                                C1::Function, C2::ParametrizedMatrix{T}, C3::Function,
                                D1::Function, D2::ParametrizedMatrix{T},
                                H::Function, R::ParametrizedMatrix{T},
                                x1::ParametrizedMatrix{T}, P1::ParametrizedMatrix{T})
        nx, ny, nu, nq, nr = confirm_matrix_sizes(A1(1), A2, A3(1), B1(1), B2, G(1), Q,
                                                  C1(1), C2, C3(1), D1(1), D2, H(1), R, x1, P1)
        new(A1, A2, A3, B1, B2, G, Q, C1, C2, C3, D1, D2, H, R, x1, P1, nx, ny, nu, nq, nr)
    end

end #ParametrizedSSM

ParametrizedSSM{T}(A2::ParametrizedMatrix{T}, Q::ParametrizedMatrix{T},
                          C2::ParametrizedMatrix{T}, R::ParametrizedMatrix{T},
                          x1::ParametrizedMatrix{T}, P1::ParametrizedMatrix{T};
                          A1::Function=_->eye(size(A2,1)), A3::Function=_->eye(size(A2,2)),
                          B1::Function=_->eye(size(A1(1),1)),
                          B2::ParametrizedMatrix{T}=parametrize_none(zeros(size(A1(1),1), 1))[1],
                          G::Function=_->eye(size(A1(1),1)),
                          C1::Function=_->eye(size(C2, 1)), C3::Function=_->eye(size(C2,2)),
                          D1::Function=_->eye(size(C1(1),1)),
                          D2::ParametrizedMatrix{T}=parametrize_none(zeros(size(C1(1),1), 1))[1],
                          H::Function=_->eye(size(C1(1),1))) =
          ParametrizedSSM{T}(A1, A2, A3, B1, B2, G, Q, C1, C2, C3, D1, D2, H, R, x1, P1)

# State space model parameters
immutable SSMParameters{T}

    # Process transition and noise covariance
    A::Vector{T}
    B::Vector{T}
    Q::Vector{T}

    # Observation and noise covariance
    C::Vector{T}
    D::Vector{T}
    R::Vector{T}

    # Initial state and error covariance
    x1::Vector{T}
    P1::Vector{T}

end #SSMParameters

SSMParameters{T}(_::T; A::Vector{T}=T[], B::Vector{T}=T[], Q::Vector{T}=T[],
                  C::Vector{T}=T[], D::Vector{T}=T[], R::Vector{T}=T[],
                  x1::Vector{T}=T[], P1::Vector{T}=T[]) =
              SSMParameters{T}(A, B, Q, C, D, R, x1, P1)

function Base.call(m::ParametrizedSSM, p::SSMParameters)
    A(t) = m.A1(t) * m.A2(p.A) * m.A3(t)
    B(t) = m.B1(t) * m.B2(p.B)
    V(t) = m.G(t) * m.Q(p.Q) * m.G(t)'
    C(t) = m.C1(t) * m.C2(p.C) * m.C3(t)
    D(t) = m.D1(t) * m.D2(p.D)
    W(t) = m.H(t) * m.R(p.R) * m.H(t)'
    x1 = vec(m.x1(p.x1))
    P1 = m.P1(p.P1)
    return StateSpaceModel(A, V, C, W, x1, P1, B=B, D=D)
end #call

function confirm_matrix_sizes(F, B, V, G, D, W, x1, P1)

    nx = size(F, 1)
    nu = size(B, 2)
    ny = size(G, 1)

    @assert size(F) == (nx, nx)
    @assert size(B) == (nx, nu)
    @assert size(V) == (nx, nx)

    @assert size(G) == (ny, nx)
    @assert size(D) == (ny, nu)
    @assert size(W) == (ny, ny)

    @assert length(x1) == nx
    @assert size(P1) == (nx, nx)

    return nx, ny, nu

end #confirm_matrix_sizes

function confirm_matrix_sizes(A1, A2, A3, B1, B2, G, Q, C1, C2, C3, D1, D2, H, R, x1, P1)

    @assert size(B2, 2) == size(D2, 2)

    nx = size(A1, 1)
    ny = size(C1, 1)
    nu = size(B2, 2)

    na1, na2  = size(A1, 2), size(A2, 2)
    nb        = size(B1, 2)
    nc1, nc2  = size(C1, 2), size(C2, 2)
    nd        = size(D1, 2)

    nq = size(Q, 1)
    nr = size(R, 1)

    @assert size(A1) == (nx, na1)
    @assert size(A2) == (na1, na2)
    @assert size(A3) == (na2, nx)

    @assert size(B1) == (nx, nb)
    @assert size(B2) == (nb, nu)

    @assert size(G)  == (nx, nq)
    @assert size(Q)  == (nq, nq)

    @assert size(C1) == (ny, nc1)
    @assert size(C2) == (nc1, nc2)
    @assert size(C3) == (nc2, nx)

    @assert size(D1) == (ny, nd)
    @assert size(D2) == (nd, nu)

    @assert size(H)  == (ny, nr)
    @assert size(R)  == (nr, nr)

    @assert length(x1)  == nx
    @assert size(P1)    == (nx, nx)

    return nx, ny, nu, nq, nr

end #confirm_matrix_sizes

"
Generates a realization of a state space model.

Arguments:
  model : StateSpaceModel
  Model defining the process

  n : Int
  Number of steps to simulate.

  u : Array (optional)
  n x nu array of exogenous inputs
"
function simulate(model::StateSpaceModel, n::Int; u::Array=zeros(n, model.nu))
    @assert size(u, 1) == n
    @assert size(u, 2) == model.nu

    # create empty arrays to hold the state and observed series
    x = zeros(model.nx, n)
    y = zeros(model.ny, n)
    u = u'

    # Cholesky decompositions of the covariance matrices, for generating
    # random noise
    V_chol(t) = all(model.V(t) .== 0) ? model.V(t) : chol(model.V(t), Val{:L})
    W_chol(t) = all(model.W(t) .== 0) ? model.W(t) : chol(model.W(t), Val{:L})

    # Generate the series
    x[:, 1] = model.x1
    y[:, 1] = model.C(1) * model.x1 + model.D(1) * u[:, 1] + W_chol(1) * randn(model.ny)
    for i=2:n
        x[:, i] = model.A(i-1) * x[:, i-1] + model.B(i-1) * u[:, i-1] + V_chol(i-1) * randn(model.nx)
        y[:, i] = model.C(i)   * x[:, i]   + model.D(i)   * u[:, i]   + W_chol(i) * randn(model.ny)
    end

    return x', y'
end


text:
# A DBserver struct contains the connect information for a database.
struct DBserver
    instanceName::String
    host::String
    user::String
    pass::String
    dbType::String
    Graphulo
end

using JavaCall

# dbsetup calls the DBserver constructor given configuration information.
# When provided only the instance name, dbsetup assumes the default configuration directory.
#   The default location is that on the MIT Supercloud.
#   Assumes existance of two files: accumulo_user_password.txt and dnsname containing password and hostname, respectively.
#   Assumes username is AccumuloUser.
# When provided a directory, will assume configuration files are in that directory.
#   Assumes existance of two files: accumulo_user_password.txt and dnsname containing password and hostname, respectively.
#   Assumes username is AccumuloUser.
# When provided a configuration file, will read in and use the configuration in the file.
#   See example configuration file for formatting.

function dbinit()
    if ~JavaCall.isloaded()
        JavaCall.addClassPath(joinpath(dirname(pathof(D4M)),"..","libext","*"))
        JavaCall.addClassPath(joinpath(dirname(pathof(D4M)),"..","lib","graphulo-3.0.0.jar"))
        JavaCall.init()
    else
        println("JVM already initialized")
    end

    # Check if correct packages are on classpath
    System = @jimport java.lang.System
    cpath = jcall(System, "getProperty", JString, (JString,), "java.class.path")

    if ~occursin("graphulo", cpath)
        println("Required Graphulo libraries for database operations missing from Java classpath.")
        println("To fix, add the required libraries (see Database Use in the readme).")
        println("Then restart Julia, and intialize jvm using dbinit().")
    end
end

function dbsetup(instance::AbstractString, hostname::AbstractString, username::AbstractString, pword::AbstractString)
    # Setup with arguments; insecure but convenient.

    if ~JavaCall.isloaded()
        println("Starting up JVM for DB operations")
        dbinit()
    end

    g = @jimport "edu.mit.ll.graphulo.MatlabGraphulo"
    Graphulo = g((JString, JString, JString, JString,), instance, hostname, username, pword)
    DB = DBserver(instance,hostname,username,pword,"BigTableLike",Graphulo)
    return DB
end

function dbsetup(instance, config="/home/gridsan/tools/groups/")
    # Note default value for config location is the MIT Supercloud default location
    if isdir(config) # Config dir
        dbdir = config*"/databases/"*instance
        f = open(dbdir*"/accumulo_user_password.txt","r")
        pword = read(f, String)
        username = "AccumuloUser"
        f = open(dbdir*"/dnsname","r")
        hostname = replace(read(f, String),"\n" => "")*":2181"
    else # Config file
        f = open(config)
        conf = readlines(config)
        conf = Dict(l[1] => l[2] for l in split.(conf,"="))
        instance = conf["instance"]
        hostname = conf["hostname"]
        username = conf["username"]
        pword = conf["password"]
    end

    return dbsetup(instance, hostname, username, pword)
end

# ls returns a list of tables that exist in the DBserver DB.
function ls(DB::DBserver)
    dbInfo = @jimport "edu.mit.ll.d4m.db.cloud.D4mDbInfo"
    dbInfoObj = dbInfo((JString, JString, JString, JString,), DB.instanceName, DB.host, DB.user, DB.pass)

    tables = jcall(dbInfoObj, "getTableList", JString, (),)

    return split(tables[1:end-1],tables[end])
end

# getindex returns a binding to a table. If the table does it exist, it creates the table.
# A database table struct is returned.
function getindex(DB::DBserver,tableName::String)

    ops = @jimport "edu.mit.ll.d4m.db.cloud.D4mDbTableOperations"
    opsObj = ops((JString, JString, JString, JString,), DB.instanceName, DB.host, DB.user, DB.pass)

    if ~any(ls(DB) .== tableName) # Create new table if it doesn't exist
        println("Creating "*tableName*" in "*DB.instanceName);
        jcall(opsObj,"createTable", Nothing, (JString,), tableName)
    end

    d4mQuery = @jimport "edu.mit.ll.d4m.db.cloud.D4mDataSearch"
    queryObj = d4mQuery((JString, JString, JString, JString, JString,), DB.instanceName, DB.host, tableName, DB.user, DB.pass)

    return DBtable(DB, tableName, "", 0, 0, "", 5e5, queryObj,opsObj)

end

# getindex returns a binding to a table. If the table does it exist, it creates the table.
# When a second table name is provided, a database table pair is returned.
function getindex(DB::DBserver,tableName1::String,tableName2::String)

    ops = @jimport "edu.mit.ll.d4m.db.cloud.D4mDbTableOperations"
    opsObj = ops((JString, JString, JString, JString,), DB.instanceName, DB.host, DB.user, DB.pass)

    # Create new tables if they don't exist
    if ~any(ls(DB) .== tableName1) || ~any(ls(DB) .== tableName2)
        if ~any(ls(DB) .== tableName1)
            println("Creating "*tableName1*" in "*DB.instanceName)
            jcall(opsObj,"createTable", Nothing, (JString,), tableName1)
        end
        if ~any(ls(DB) .== tableName2)
            println("Creating "*tableName2*" in "*DB.instanceName)
            jcall(opsObj,"createTable", Nothing, (JString,), tableName2)
        end
    end

    d4mQuery = @jimport "edu.mit.ll.d4m.db.cloud.D4mDataSearch"
    queryObj = d4mQuery((JString, JString, JString, JString, JString,), DB.instanceName, DB.host, tableName1, DB.user, DB.pass)

    return DBtablePair(DB, tableName1, tableName2, "", 0, 0, "", 5e5, queryObj,opsObj)

end

# Returns whether a table with the given name exists in Accumulo
function ispresent(DB::DBserver, tableName1::String)
    return any(ls(DB) .== tableName1)
end

# Deletes a table with the given name.
function deletename(DB::DBserver, name::AbstractString)
    defaultnames = ["accumulo.metadata", "accumulo.replication", "accumulo.root", "trace"]
    if any(defaultnames .== name)
        println("Table with name " * name * " is a default table")
    elseif ispresent(DB, name)
        delete(DB[String(name)])
    else
        println("Table with name " * name * " does not exist in " * DB.instanceName)
    end
end

# Deletes all tables, except for the 4 default ones in the Accumulo database.
function deleteall(DB::DBserver)
    defaultnames = ["accumulo.metadata", "accumulo.replication", "accumulo.root", "trace"]
    for name in ls(DB)
        if ~(any(defaultnames .== name))
            delete(DB[String(name)])
        end
    end
end

# Deletes all tables, withh a given prefix, except for the 4 default ones in the Accumulo database.
function deleteprefix(DB::DBserver, prefix::AbstractString)
    defaultnames = ["accumulo.metadata", "accumulo.replication", "accumulo.root", "trace"]
    for name in ls(DB)
        if ~(any(defaultnames .== name)) && startswith(name, prefix)
                delete(DB[String(name)])
        end
    end
end


text:
using LinearAlgebra
function powermethod(A;maxiter=100000, tol=1e-6)
  x = normalize!(randn(size(A,1)))
  y = A*x
  λ = x'*y
  for i=1:maxiter
    if norm(y - λ*x)/(abs(λ)+1) <= tol
      break
    end
    x = normalize!(y)
    y = A*x
    λ = x'*y
  end
  if norm(y - λ*x)/(abs(λ)+1) >= tol
    @warn("powermethod did not converge, relative residual is ",
      norm(y - λ*x)/(abs(lam)+1))
  end
  return x, λ
end
##
using Random
Random.seed!(2)
A = randn(3,3)
v, lam = powermethod(A)
##
Lams, V = eigen(A)
##
@show lam
@show Lams[1]
##
function powermethod_explain(A,V;maxiter=9)
  x = normalize!(randn(size(A,1)))
  y = A*x
  for i=1:9
    x = normalize!(y)
    y = A*x
    println(i, ", ", V'*x)
  end
end
using Random
Random.seed!(2)
A = randn(3,3)
A = A+A' # make it symmetric
Lams, V = eigen(A)
powermethod_explain(A,V)
## example on population growth matrix
G = [1.0215 0.0127; 0.0627 1.0667 ]
v, lam = powermethod(G)
@show  g = v / sum(v), lam


text:
"""
    Graph(; kws...)

A type representing a graph (a set of nodes and edges) that can also store additional node and edge data.

Graph datasets in MLDatasets.jl contain one or more `Graph` objects.

# Fields

- num_nodes
- num_edges
- edge_index
- node_data
- edge_data
"""
Base.@kwdef struct Graph
    num_nodes::Int = 0
    num_edges::Int = 0
    edge_index::Tuple{Vector{Int}, Vector{Int}} = (Int[], Int[])
    node_data = nothing
    edge_data = nothing
end

function Base.show(io::IO, d::Graph)
    print(io, "Graph($(d.num_nodes), $(d.num_edges))")
end

function Base.show(io::IO, ::MIME"text/plain", d::Graph)
    recur_io = IOContext(io, :compact => false)
    print(io, "Graph:")
    for f in fieldnames(Graph)
        if !startswith(string(f), "_")
            fstring = leftalign(string(f), 10)
            print(recur_io, "\n  $fstring  =>    ")
            print(recur_io, "$(_summary(getfield(d, f)))")
        end
    end
end

function adjlist2edgeindex(adj; inneigs=true)
    s, t = Int[], Int[]
    for i in 1:length(adj)
        for j in adj[i]
            push!(s, i)
            push!(t, j)
        end
    end

    if inneigs
        return t, s
    else
        return s, t
    end
end


text:
using Bijectors: Exp
using WarpedRegressors: posterior

@testset "warped_regressor" begin
    @testset "Stheno" begin
        rng = MersenneTwister(123456)

        # Construct a GP.
        N = 3
        s = 0.1
        x = randn(rng, N)
        f = GP(SEKernel())
        fx = f(x, s)

        # Construct a WarpedGP.
        ϕ = Exp()
        ϕf = warp(f, ϕ)
        ϕfx = ϕf(x, s)

        # Ensure that sampling acts heirarchically. Not a phenominal test.
        z = rand(MersenneTwister(1), fx)
        y = ϕ.(z)
        @test y == rand(MersenneTwister(1), ϕfx)

        # Roughly check that the logpdf is correct. This really isn't a great test.
        manual_logpdf = sum(logabsdetjac.(Ref(inv(ϕ)), y)) + logpdf(fx, inv(ϕ).(y))
        @test manual_logpdf ≈ logpdf(ϕfx, y)

        # Check that gradients of the lml w.r.t. the parameters can be computed.
        function warped_regressor_logpdf(x, y, s)
            ϕf = warp(GP(SEKernel()), Exp())
            return logpdf(ϕf(x, s), y)
        end
        Zygote.gradient(warped_regressor_logpdf, x, y, s)

        # Check that posterior does the correct thing. This really also isn't a great test.
        manual_posterior = warp(posterior(fx, z), ϕ)
        post = posterior(ϕfx, y)
        @test logpdf(manual_posterior(x), y) ≈ logpdf(post(x), y)

        # Check that an affine transformation of a GP given the same as scaling the GP.
        let
            α = randn(rng)
            αf = GP(α^2 * SEKernel())
            ϕf = warp(f, Bijectors.Scale(α))
            y = rand(rng, αf(x))
            @test logpdf(αf(x), y) ≈ logpdf(ϕf(x), y)
        end
    end
    @testset "BayesianLinearRegressors" begin
        rng = MersenneTwister(123456)

        # Construct a regressor.
        N = 11
        D = 3
        s = 0.1
        x = randn(rng, D, N)
        f = BayesianLinearRegressor(ones(D), Matrix{Float64}(I, D, D))
        fx = f(x, s)

        # Construct a WarpedGP.
        ϕ = Exp()
        ϕf = warp(f, ϕ)
        ϕfx = ϕf(x, s)

        # Ensure that sampling acts heirarchically. Not a phenominal test.
        z = rand(MersenneTwister(1), fx)
        y = ϕ.(z)
        @test y == rand(MersenneTwister(1), ϕfx)

        # Roughly check that the logpdf is correct. This really isn't a great test.
        manual_logpdf = sum(logabsdetjac.(Ref(inv(ϕ)), y)) + logpdf(fx, inv(ϕ).(y))
        @test manual_logpdf ≈ logpdf(ϕfx, y)

        # Check that gradients of the lml w.r.t. the parameters can be computed.
        function warped_regressor_logpdf(x, y, s)
            f = BayesianLinearRegressor(ones(D), Matrix{Float64}(I, D, D))
            ϕf = warp(f, ϕ)
            return logpdf(ϕf(x, s), y)
        end
        Zygote.gradient(warped_regressor_logpdf, x, y, s)

        # Check that posterior does the correct thing. This really also isn't a great test.
        manual_posterior = warp(BayesianLinearRegressors.posterior(fx, z), ϕ)
        post = posterior(ϕfx, y)
        @test logpdf(manual_posterior(x, 1e-6), y) ≈ logpdf(post(x, 1e-6), y)
    end
end


text:
# greedy seeds heuristic by A. Catalan and M. Fisher (2012) doi:10. 2139/ssrn.2166687
function GREEDYSEEDS(trans::SparseMatrixCSC{Float64, Int64},
                     Q::Matrix{Int64},
                     capacity::Array{Int64,1})
    if CHECKCAPACITY(Q::Array{Int64,2}, capacity::Array{Int64,1}) == 1
        ## Calculate the coapperance matrix
        capacity_left = copy(capacity) .- 1
        ## Sort the SKUs by decreasing sales
        sales = SORTSALES(trans)
        ## Assign the top (best selling) SKU to the largest DC
        X = Array{Int64,2}(undef,size(trans,2),size(capacity,1)) .= 0
        X[sales[1,1],1] = 1
        ## Mark the SKU that is now allocated and the capacity used
        sales[1,4] = 1
        ## Start the assignment of the other SKUs
        GREEDYSEEDSTART!(sales::Array{Int64,2}, X::Array{Int64,2}, Q::Array{Int64,2}, capacity_left::Array{Int64,1})
        ## For each SKUs sorted by decreasing sales that has not been allocated yet
        GREEDYSEEDMAIN!(sales::Array{Int64,2}, X::Array{Int64,2}, Q::Array{Int64,2}, capacity_left::Array{Int64,1})
        ## While there is unallocated space in the DCs do
        X = FILLUP(X::Array{Int64,2}, Q::Array{Int64,2}, capacity_left::Array{Int64,1})
        ## Return the solution
        return X::Array{Int64,2}
    end
end

max_stars_repo_name: beyondatlas/splitdeliveries

max_issues_repo_name: beyondatlas/splitdeliveries

max_forks_repo_name: beyondatlas/splitdeliveries


text:
using LightXML

struct AnimatedTile
    name::String
    path::String

    origX::Number
    origY::Number

    posX::Number
    posY::Number

    delay::Number
end

struct Mask
    tiles::UInt8
    ignores::UInt8
end

# Relevant characters of the tile mask string
# We never need to check the centre tile, otherwise we wouldn't be in the mask check
const lookupmask = [1, 2, 3, 5, 7, 9, 10, 11]

# Convert mask string into a bitmask.
# ignoremask is 1 for every tile that should NOT be ignored.
function getMaskFromString(s::String)
    data = s[lookupmask]
    tilemask = UInt8(0)
    ignoremask = UInt8(0)
    for (i, c) in enumerate(data)
        tilemask |= UInt8(c == '1') << (8 - i)
        ignoremask |= UInt8(c != 'x') << (8 - i)
    end
    return Mask(tilemask, ignoremask)
end

function sortScore(mask::Mask)
    # Counts the number of 0-bits in the ignore string
    # by counting the number of 1-bits.
    return 8 - sum([(mask.ignores >> i) & 0b1 for i = 0:7])
end

# Get quads for a given mask
function convertTileString(s::String)
    parts = split(s, ";")
    res = Coord[]

    for part in parts
        x, y = parse.(Int, split(part, ","))

        push!(res, Coord(x, y))
    end

    return res
end

function loadAnimatedTilesXML(fn::String)
    xdoc = parse_file(fn)
    xroot = root(xdoc)

    animatedTiles = AnimatedTile[]

    for elem in child_elements(xroot)
        push!(animatedTiles, AnimatedTile(
            attribute(elem, "name"),
            attribute(elem, "path"),

            parseNumber(attribute(elem, "origX")),
            parseNumber(attribute(elem, "origY")),

            parseNumber(attribute(elem, "posX")),
            parseNumber(attribute(elem, "posY")),

            parseNumber(attribute(elem, "delay"))
        ))
    end

    return animatedTiles
end

struct Coord
    x::Int
    y::Int
end

struct MaskData
    mask::Mask
    coords::Array{Coord, 1}
    sprites::String
end

struct Ignore
    wildcard::Bool
    ignores::String

    Ignore(s::String) = new('*' in s, filter(c -> c != '*', s))
end

function loadTilesetXML(fn::String)
    xdoc = parse_file(fn)
    xroot = root(xdoc)

    paths = Dict{Char, String}()
    masks = Dict{Char, Array{MaskData, 1}}()
    padding = Dict{Char, Array{Coord, 1}}()
    center = Dict{Char, Array{Coord, 1}}()
    ignores = Dict{Char, Ignore}()

    orderedElements = sort(collect(child_elements(xroot)), by=set -> attribute(set, "copy") !== nothing)

    for set in orderedElements
        id = attribute(set, "id")[1] # Char doesn't work on 1 length string?
        path = attribute(set, "path")
        copyTileset = attribute(set, "copy")
        ignore = attribute(set, "ignores")

        paths[id] = "tilesets/$path"

        if ignore !== nothing
            ignores[id] = Ignore(ignore)
        end

        if copyTileset !== nothing
            tilesetId = copyTileset[1]
            padding[id] = deepcopy(padding[tilesetId])
            center[id] = deepcopy(center[tilesetId])
            masks[id] = deepcopy(masks[tilesetId])
        end

        currMasks = MaskData[]

        for c in child_elements(set)
            mask = attribute(c, "mask")
            tilesString = attribute(c, "tiles")
            quads = convertTileString(tilesString)

            if mask == "padding"
                padding[id] = quads

            elseif mask == "center"
                center[id] = quads

            else
                sprites = attribute(c, "sprites")
                push!(currMasks, MaskData(
                    getMaskFromString(mask),
                    quads,
                    isa(sprites, String) ? sprites : ""
                ))
            end
        end

        if !isempty(currMasks)
            masks[id] = currMasks
            sort!(masks[id], by=m -> sortScore(m.mask), rev=false)
        end
    end

    return paths, masks, padding, center, ignores
end

struct TilerMeta
    paths::Dict{Char, String}
    masks::Dict{Char, Array{MaskData, 1}}
    paddings::Dict{Char, Array{Coord, 1}}
    centers::Dict{Char, Array{Coord, 1}}
    ignores::Dict{Char, Ignore}
end

TilerMeta(s::String) = TilerMeta(loadTilesetXML(s)...)

fgTilerMeta = nothing
bgTilerMeta = nothing
animatedTilesMeta = nothing

function loadFgTilerMeta(side::Union{Side, Nothing}, fn::String, loadDefault::Bool=false)
    path = joinpath(storageDirectory, "XML", "ForegroundTiles.xml")

    if side !== nothing && !loadDefault
        meta = get(Dict{String, Any}, side.data, "meta")
        metaPath = get(meta, "ForegroundTiles", "")

        hasRoot, modRoot = getModRoot(fn)
        xmlPath = joinpath(modRoot, metaPath)

        if !isempty(metaPath) && hasRoot
            path = xmlPath
        end
    end

    try
        global fgTilerMeta = TilerMeta(path)

    catch e
        if !loadDefault
            println(Base.stderr, "Failed to load custom ForegroundTiles XML")
            println(Base.stderr, e)

            for (exc, bt) in Base.catch_stack()
                showerror(Base.stderr, exc, bt)
                println()
            end

            loadFgTilerMeta(side, fn, true)
        end
    end
end

function loadBgTilerMeta(side::Union{Side, Nothing}, fn::String, loadDefault::Bool=false)
    path = joinpath(storageDirectory, "XML", "BackgroundTiles.xml")

    if side !== nothing && !loadDefault
        meta = get(Dict{String, Any}, side.data, "meta")
        metaPath = get(meta, "BackgroundTiles", "")

        hasRoot, modRoot = getModRoot(fn)
        xmlPath = joinpath(modRoot, metaPath)

        if !isempty(metaPath) && hasRoot
            path = xmlPath
        end
    end

    try
        global bgTilerMeta = TilerMeta(path)

    catch e
        if !loadDefault
            println(Base.stderr, "Failed to load custom BackgroundTiles XML")
            println(Base.stderr, e)

            for (exc, bt) in Base.catch_stack()
                showerror(Base.stderr, exc, bt)
                println()
            end

            loadBgTilerMeta(side, fn, true)
        end
    end
end

function loadAnimatedTilesMeta(side::Union{Side, Nothing}, fn::String, loadDefault::Bool=false)
    path = joinpath(storageDirectory, "XML", "AnimatedTiles.xml")

    if side !== nothing && !loadDefault
        meta = get(Dict{String, Any}, side.data, "meta")
        metaPath = get(meta, "AnimatedTiles", "")

        hasRoot, modRoot = getModRoot(fn)
        xmlPath = joinpath(modRoot, metaPath)

        if !isempty(metaPath) && hasRoot
            path = xmlPath
        end
    end

    try
        global animatedTilesMeta = loadAnimatedTilesXML(path)

    catch e
        if !loadDefault
            println(Base.stderr, "Failed to load custom AnimatedTiles XML")
            println(Base.stderr, e)

            for (exc, bt) in Base.catch_stack()
                showerror(Base.stderr, exc, bt)
                println()
            end

            loadAnimatedTilesMeta(side, fn, true)
        end
    end
end

function loadXMLMeta()
    side = loadedState.side
    filename = loadedState.filename

    loadFgTilerMeta(side, filename)
    loadBgTilerMeta(side, filename)
    loadAnimatedTilesMeta(side, filename)
end

function getTile(tiles::Tiles, x::Int, y::Int)
    return get(tiles.data, (y, x), ' ')
end

function checkPadding(data::Array{Char, 2}, x::Int, y::Int, width::Int, height::Int)
    # Checks for '0' even though getIfInboundsKnownSize might return ' '
    # This is due to quirks/special flags in the actuall autotiler
    return getIfInboundsKnownSize(data, y - 2, x, height, width, ' ') == '0' ||
        getIfInboundsKnownSize(data, y + 2, x, height, width, ' ') == '0' ||
        getIfInboundsKnownSize(data, y, x - 2, height, width, ' ') == '0' ||
        getIfInboundsKnownSize(data, y, x + 2, height, width, ' ') == '0'
end

const ignoresDefault = Ignore("")

# Return the "mask" value for a tile
# Returns as UInt8 for performance
# Manual string iteratation for performance
function checkTile(value::Char, target::Char, ignore::Ignore=ignoresDefault)
    if target == '0' || (ignore.wildcard && value != target)
        return 0b0

    else
        for c in ignore.ignores
            if target == c
                return 0b0
            end
        end
    end

    return 0b1
end

function checkMask(mask::Mask, adj::UInt8)
    return ((adj ⊻ mask.tiles) & mask.ignores) == 0
end

function getIfInboundsKnownSize(data::Array{Char, 2}, y::Int, x::Int, height::Int, width::Int, default::Char='0')::Char
    if 1 <= y <= height && 1 <= x <= width
        return @inbounds data[y, x]

    else
        return default
    end
end

# Unrolled for performance
function getAdjacencyMask(data::Array{Char, 2}, y::Int, x::Int, height::Int, width::Int, default::Char, ignores::Ignore)::UInt8
    return checkTile(default, getIfInboundsKnownSize(data, y - 1, x - 1, height, width, default), ignores) << 7 |
        checkTile(default, getIfInboundsKnownSize(data, y - 1, x, height, width, default), ignores) << 6 |
        checkTile(default, getIfInboundsKnownSize(data, y - 1, x + 1, height, width, default), ignores) << 5 |

        checkTile(default, getIfInboundsKnownSize(data, y, x - 1, height, width, default), ignores) << 4 |
        checkTile(default, getIfInboundsKnownSize(data, y, x + 1, height, width, default), ignores) << 3 |

        checkTile(default, getIfInboundsKnownSize(data, y + 1, x - 1, height, width, default), ignores) << 2 |
        checkTile(default, getIfInboundsKnownSize(data, y + 1, x, height, width, default), ignores) << 1 |
        checkTile(default, getIfInboundsKnownSize(data, y + 1, x + 1, height, width, default), ignores)
end

const defaultPaddingCoords = Coord[Coord(0, 0)]
const maskDataDefault = MaskData[]

function getMaskQuads(x::Int, y::Int, tiles::Tiles, meta::TilerMeta)
    data = tiles.data
    height, width = size(data)
    value = @inbounds tiles.data[y, x]

    masks = get(meta.masks, value, maskDataDefault)
    ignores = get(meta.ignores, value, ignoresDefault)

    adjacent = getAdjacencyMask(data, y, x, height, width, value, ignores)

    for data in masks
        if checkMask(data.mask, adjacent)
            return data.coords, data.sprites
        end
    end

    if checkPadding(tiles.data, x, y, width, height)
        return get(meta.paddings, value, defaultPaddingCoords), ""

    else
        return get(meta.centers, value, defaultPaddingCoords), ""
    end
end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
using Javis
using FFTW
using FFTViews
using FileIO
using Images
using TravelingSalesmanHeuristics

function ground(args...)
    background("black")
    sethue("white")
end

function circ(; r = 10, vec = O, action = :stroke, color = "white")
    sethue(color)
    circle(O, r, action)
    my_arrow(O, vec)
    return vec
end

function my_arrow(start_pos, end_pos)
    arrow(
        start_pos,
        end_pos;
        linewidth = distance(start_pos, end_pos) / 100,
        arrowheadlength = 7,
    )
    return end_pos
end

function draw_line(
    p1 = O,
    p2 = O;
    color = "white",
    action = :stroke,
    edge = "solid",
    linewidth = 3,
)
    sethue(color)
    setdash(edge)
    setline(linewidth)
    line(p1, p2, action)
end

function draw_path!(path, pos, color)
    sethue(color)

    push!(path, pos)
    return draw_line.(path[2:end], path[1:(end - 1)]; color = color)
end

function get_points(img)
    findall(x -> x == 1, img) .|> x -> Point(x.I)
end

c2p(c::Complex) = Point(real(c), imag(c))

remap_idx(i::Int) = (-1)^i * floor(Int, i / 2)
remap_inv(n::Int) = 2n * sign(n) - 1 * (n > 0)

function animate_fourier(options)
    npoints = options.npoints
    nplay_frames = options.nplay_frames
    nruns = options.nruns
    nframes = nplay_frames + options.nend_frames

    # obtain points from julialogo
    points = get_points(load(File(format"PNG", "thin_img.png")))
    npoints = length(points)
    println("#points: $npoints")
    # solve tsp to reduce length of extra edges
    distmat = [distance(points[i], points[j]) for i = 1:npoints, j = 1:npoints]

    path, cost = solve_tsp(distmat; quality_factor = options.tsp_quality_factor)
    println("TSP cost: $cost")
    points = points[path] # tsp saves the last point again

    # optain the fft result and scale
    y = [p.x - options.width for p in points] ./ 3
    x = [p.y - options.height for p in points] ./ 3

    fs = FFTView(fft(complex.(x, y)))
    # normalize the points as fs isn't normalized
    fs ./= npoints
    npoints = length(fs)

    video = Video(options.width, options.height)
    Background(1:nframes, ground)

    circles = Object[]

    for i = 1:npoints
        ridx = remap_idx(i)

        push!(circles, Object((args...) -> circ(; r = abs(fs[ridx]), vec = c2p(fs[ridx]))))

        if i > 1
            # translate to the tip of the vector of the previous circle
            act!(circles[i], Action(1:1, anim_translate(circles[i - 1])))
        end
        ridx = remap_idx(i)
        act!(circles[i], Action(1:nplay_frames, anim_rotate(0.0, ridx * 2π * nruns)))
    end

    trace_points = Point[]
    Object(1:nframes, (args...) -> draw_path!(trace_points, pos(circles[end]), "pink"))

    return render(video; pathname = joinpath(@__DIR__, options.filename))
    # return render(video; liveview = true)
end

function main()
    # hd_options = (
    # npoints = 3001, # rough number of points for the shape => number of circles
    # nplay_frames = 1200, # number of frames for the animation of fourier
    # nruns = 2, # how often it's drawn
    # nend_frames = 200,  # number of frames in the end
    # width = 1920,
    # height = 1080,
    # shape_scale = 2.5, # scale factor for the logo
    # tsp_quality_factor = 50,
    # filename = "julia_hd.mp4",
    # )

    gif_options = (
        npoints = 1001, # rough number of points for the shape => number of circles
        nplay_frames = 600, # number of frames for the animation of fourier
        nruns = 1, # how often it's drawn
        nend_frames = 200,  # number of frames in the end
        width = 455,
        height = 490,
        tsp_quality_factor = 20,
        filename = "julia_fast.gif",
    )

    # gif_options = (
    # npoints = 651, # rough number of points for the shape => number of circles
    # nplay_frames = 600, # number of frames for the animation of fourier
    # nruns = 2, # how often it's drawn
    # nend_frames = 0,  # number of frames in the end
    # width = 350,
    # height = 219,
    # shape_scale = 0.8, # scale factor for the logo
    # tsp_quality_factor = 80,
    # filename = "julia_logo_dft.gif",
    # )
    return animate_fourier(gif_options)
end

max_stars_repo_name: Ved-Mahajan/JavisIsAwsome

max_issues_repo_name: Ved-Mahajan/JavisIsAwsome

max_forks_repo_name: Ved-Mahajan/JavisIsAwsome


text:
# algorithm type
struct ButuceaComte end

"""
    ButuceaComteOptions(; bandwidth = :auto)

The Butucea Comte estimator of linear functional in the deconvolution model. Given a `LinearEBayesTarget`
with characteristic function ``\\psi^*`` and samples convolved with Gaussian noise (with ``\\phi`` Standard Gaussian pdf),
it is defined as follows

```math
\\hat{L}_{\\text{BC},h_n} = \\frac{1}{2 \\pi n}\\sum_{i=1}^n \\int_{-1/h_n}^{1/h_n} \\exp(it Z_k) \\frac{\\psi^*(-t)}{\\varphi^*(t)}dt
```
``h`` is the bandwidth, the option `:auto` will pick it automatically.

## Reference:
  >Butucea, C. and Comte, F., 2009.
  >Adaptive estimation of linear functionals in the convolution model and applications.
  >Bernoulli, 15(1), pp.69-98.
"""
Base.@kwdef struct ButuceaComteOptions{S}
    bandwidth::S = :auto
end

fit(bcopt::ButuceaComteOptions, Zs) = bcopt #do nothing

struct ButuceaComteEstimator{EBT<:LinearEBayesTarget}
    target::EBT
    h::Float64
    f_real #ApproxFun based
    f_imag #ApproxFun based
end

# Make theoretical results of BC practical.
Base.@kwdef struct ButuceaComteNoiseDensityParams{T <: Real}
    γ::T
    α::T
    ρ::T
    κ0upper::T
    κ0lower::T
end

function ButuceaComteNoiseDensityParams(Z::StandardNormalSample)
    ButuceaComteNoiseDensityParams(γ = 0.0, ρ = 2.0, α = 0.5, κ0upper=1.0, κ0lower=1.0)
end

Base.@kwdef struct ButuceaComteSobolevParams{T<:Real}
   b::T
   L::T        # constant
   a::T = 0.0
   r::T = 0.0
end

Base.@kwdef struct ButuceaComteTargetParams{T<:Real}
   B::T        #polynomial decay (1+x^2)^B
   A::T
   R::T
   Cψ::T       #constant
end

function ButuceaComteTargetParams(target::PriorDensityTarget)
    ButuceaComteTargetParams(A=0.0, B=0.0, R=0.0, Cψ=1.0)
end

# only accounting for b,B for now
function butucea_comte_squared_bias_bound(noise::ButuceaComteNoiseDensityParams,
                                  sobolev::ButuceaComteSobolevParams,
                                  target::ButuceaComteTargetParams,
                                  h::Real)
    R = target.R
    r = sobolev.r

    a = target
    B = target.B
    b = sobolev.b

    L = sobolev.L
    Cψ = target.Cψ

    L*Cψ/π/(2b + 2B -1)*h^(1-2b-2B)
end

# only do case PriorDensityTarget for now
function butucea_comte_unit_var_proxy(h::Real)
    1/π^2*exp(h^2/2)/h^2
end

function default_bandwidth(::Type{ButuceaComteEstimator},
                           target::MarginalDensityTarget{<:StandardNormalSample},
                           n)
    sqrt(log(n))
end

default_bandwidth(::Type{ButuceaComteEstimator}, target::LFSRNumerator, n) = sqrt(log(n)*2/3)

default_bandwidth(::Type{ButuceaComteEstimator}, target::PosteriorMeanNumerator, n) = sqrt(log(n))

function default_bandwidth(::Type{ButuceaComteEstimator},
                           target::PriorDensityTarget,
                           n)
    # should make this also depend on noise distribution etc ...
    # π*m_n  = (log(n)/(2α + 1))^{1/ρ}
    #gaussian case α=1//2, ρ=2
    sqrt(log(n)/2)
end

function ButuceaComteEstimator(target::EBayesTarget; n::Integer)
    h = default_bandwidth(ButuceaComteEstimator, target, n)
    ButuceaComteEstimator(target, h)
end

function ButuceaComteEstimator(target::EBayesTarget, h::Float64)
    #todo: remove hardcoding of Normal
    _f(t) = cf(target,-t)/cf(Normal(),t)
    _set = Interval(-h,h)
    _f_real = Fun(real ∘ _f, _set)
    _f_imag = Fun(imag ∘ _f, _set)
    ButuceaComteEstimator(target, h, _f_real, _f_imag)
end

function (cb::ButuceaComteEstimator)(z)
    real_part = fourier(cb.f_real, z)
    imag_part = fourier(cb.f_imag, z)

    real(real_part + im*imag_part)/2/π
    #h = cb.h
    ##target = cb.target
    #f(t) = real(exp(im*t*z)*cf(target,-t)/cf(Normal(),t))
    #res, _ = quadgk(f, -h, +h) res/2π
end

function DiscretizedAffineEstimator(mhist, cb::ButuceaComteEstimator)
    DiscretizedAffineEstimator(mhist, z -> cb(z))
end

"""
    estimate(target::LinearEBayesTarget, bcopt::ButuceaComteOptions, Zs)

Estimate `target` based on samples `Zs` with the Butucea-Comte estimator with
[`ButuceaComteOptions`](@ref) bcopt.
"""
function estimate(target::LinearEBayesTarget, bcopt::ButuceaComteOptions,
                  Zs::AbstractVector{<:StandardNormalSample})
    (bcopt.bandwidth == :auto) || error("only :auto bandwidth for now")
    bcest = ButuceaComteEstimator(target;  n=nobs(Zs))
    mean(bcest.(response.(Zs)))
end

#function estimate(Xs,cb::ComteButucea)
#    mean(cb.Q.(Xs))
#end

#struct DiscretizedAffineEstimator{MH<:MCEBHistogram}
#    mhist::MH
#    Q::Vector{Float64}
#    Qo::Float64 #offset
#end

#default_bandwidth(::Type{ComteButucea}, n) = sqrt(log(n))
#default_bandwidth(::Type{ComteButucea}, target, n) = default_bandwidth(ComteButucea, n)


text:
cradixsort!(x) = ccall(
    (:radix_sort,"libuntitled3"),
     Void,
     (Ptr{Ptr{UInt8}},UInt),
    pointer(x), length(x)
)

# test on variable length 8
const M=100_000_000; const K=100
using FastGroupBy
srand(1);
svec1 = rand([string(rand(Char.(32:126), rand(1:8))...) for k in 1:M÷K], M);
@time cradixsort!(svec1) # 50 seconds
issorted(svec1)

srand(1);
svec1 = rand([string(rand(Char.(32:126), rand(1:8))...) for k in 1:M÷K], M);
@time radixsort!(svec1) # 18 seconds
issorted(svec1)

# test on fixed length 10
const M = 100_000_000; const K = 100;
srand(1);
svec = rand(["id"*dec(k,10) for k in 1:M÷K], M);

@time cradixsort(svec) # 2.5 minutes
issorted(svec)

using FastGroupBy
srand(1);
svec = rand(["id"*dec(k,10) for k in 1:M÷K], M);
@time radixsort!(svec) # 18 seconds
issorted(svec)


text:
# internal method essentially the same as Base.replace!(y, (z .=> r)...)
# but more efficient.
# Similar to the behaviour of `Base.replace!` if `z` contain repetions of values in
# `y` then only the transformation corresponding to the first occurence is performed
# i.e `_replace!([1,5,3], [1,4], 4:5)` would return `[4,5,3]` rather than `[5,5,3]`
# (which replaces `1=>4` and then `4=>5`)
function _replace!(y::AbstractVector, z::AbstractVector, r::AbstractVector)
    length(r) == length(z) ||
     throw(DimensionMismatch("`z` and `r` has to be of the same length"))
    @inbounds for i in eachindex(y)
        for j in eachindex(z)
            isequal(z[j], y[i]) && (y[i] = r[j]; break)
        end
    end
    return y
end

"""
 softmax(X::AbstractMatrix{<:Real})
Return the softmax computed in a numerically stable way:
``σ(x) = exp.(x) ./ sum(exp.(x))``
Implementation taken from NNlib.jl.
"""
function softmax(X::AbstractMatrix{<:Real})
    S = copyto!(similar(X, eltype(X)), X)
    return softmax!(S)
end

"""
 softmax!(X::AbstractMatrix{<:Real})
Return the softmax computed in a numerically stable way:
``σ(x) = exp.(x) ./ sum(exp.(x))`` and store the result in the
input matrix
Implementation taken from NNlib.jl.
"""
function softmax!(X::AbstractMatrix{<:Real})
    max_ = maximum(X, dims=2)
    X .= exp.(X .- max_)
    X ./= sum(X, dims=2)
    return X
end


text:
n = 0

while n < 10
    global n = n + 1
    println(n)

    for i = 1:5
        println("oi")
    end
end

a = zeros(n,n)

for i = 1:n
  for j = 1:n
    if i==j

      a[i,j] = -4

      elseif i == (j-3) || i == (j+3)

        a[i,j] = 1

      elseif (mod(i,3) != 0 && j == (i+1))

        a[i,j] = 1

      elseif (mod(i,3) != 1 && j == (i-1))

        a[i,j] = 1

      end
  end
end

println("The A matrix is: ")

println(a)

for i = 1:n

  println(a[i,:])

end

max_stars_repo_path: Linear System/Julia/useful_info.jl

max_stars_repo_name: arturofburgos/Comparing-Languages

max_issues_repo_path: Linear System/Julia/useful_info.jl

max_forks_repo_path: Linear System/Julia/useful_info.jl


text:
export Atoms

"""
    Atoms{S,T<:AbstractFloat}

Basic atomic parameters: element symbols, numbers and masses

Masses are converted to atomic units.
Constructed using either element symbols or masses.

```jldoctest
julia> Atoms(:H)
Atoms{1, Float64}([:H], UInt8[0x01], [1837.4715941070515])

julia> Atoms([:H, :H, :H, :C])
Atoms{4, Float64}([:H, :H, :H, :C], UInt8[0x01, 0x01, 0x01, 0x06], [1837.4715941070515, 1837.4715941070515, 1837.4715941070515, 21894.713607956142])

julia> Atoms([100, 200])
Atoms{2, Float64}([:X, :X], UInt8[0x00, 0x00], [100.0, 200.0])
```
"""
struct Atoms{S,T<:AbstractFloat}
    types::SVector{S,Symbol}
    numbers::SVector{S,UInt8}
    masses::SVector{S,T}
end

function Atoms{T}(atom_types::AbstractVector{Symbol}) where {T}
    S = length(atom_types)
    types = SVector{S,Symbol}(atom_types)
    numbers = SVector{S,UInt8}([element.number for element in elements[atom_types]])
    masses = SVector{S,T}([austrip(element.atomic_mass) for element in elements[atom_types]])
    Atoms{S,T}(types, numbers, masses)
end

function Atoms{T}(masses::AbstractVector) where {T}
    S = length(masses)
    types = SVector{S,Symbol}(fill(:X, S))
    numbers = SVector{S,UInt8}(zeros(S))
    masses = SVector{S,T}(austrip.(masses))
    Atoms{S,T}(types, numbers, masses)
end

Atoms{T}(atom_type) where {T} = Atoms{T}([atom_type])
Atoms(atom_types) = Atoms{Float64}(atom_types)

Base.length(::Atoms{S}) where {S} = S
Base.range(::Atoms{S}) where {S} = range(1; length=S)

Base.IndexStyle(::Atoms) = IndexLinear()
Base.getindex(A::Atoms{S,T}, i::Int) where {S,T} = Atoms{T}(A.types[i])
Base.getindex(A::Atoms{S,T}, i::AbstractRange) where {S,T} = Atoms{T}(A.types[i])


text:
@testset "Alligators" begin

@info "$(now()) Starting test set: Alligators"

#=
  Stomach contents of alligators in different lakes.
  https://bookdown.org/mpfoley1973/statistics/generalized-linear-models-glm.html#multinomial-logistic-regression
=#

################################################################################
# Supporting functions

function construct_alligator_dataframe(data, xnames, ylevels)
    result = fill("", 219, 5)
    rownum = 0
    for profile = 1:16
        pdata = view(data, :, profile)
        for outcome = 1:5
            nobs = pdata[4 + outcome]
            y    = ylevels[outcome]
            for i = 1:nobs
                rownum += 1
                result[rownum, 1] = y
                for (j, xname) in enumerate(xnames)
                    result[rownum, 1 + j] = pdata[j]
                end
            end
        end
    end
    colnames = vcat("StomachContents", xnames)
    DataFrame(result, colnames)
end

################################################################################
# Data and target parameters (see link above)

ylevels = ["Fish", "Invertebrate", "Reptile", "Bird", "Other"]  # Stomach contents
xnames  = ["profile", "Gender", "Size", "Lake"]                 # Profile of alligators
rawdata = ["1", "f", "<2.3", "george", 3, 9, 1, 0, 1,
"2", "m", "<2.3", "george", 13, 10, 0, 2, 2,
"3", "f", ">2.3", "george", 8, 1, 0, 0, 1,
"4", "m", ">2.3", "george", 9, 0, 0, 1, 2,
"5", "f", "<2.3", "hancock", 16, 3, 2, 2, 3,
"6", "m", "<2.3", "hancock", 7, 1, 0, 0, 5,
"7", "f", ">2.3", "hancock", 3, 0, 1, 2, 3,
"8", "m", ">2.3", "hancock", 4, 0, 0, 1, 2,
"9", "f", "<2.3", "oklawaha", 3, 9, 1, 0, 2,
"10", "m", "<2.3", "oklawaha", 2, 2, 0, 0, 1,
"11", "f", ">2.3", "oklawaha", 0, 1, 0, 1, 0,
"12", "m", ">2.3", "oklawaha", 13, 7, 6, 0, 0,
"13", "f", "<2.3", "trafford", 2, 4, 1, 1, 4,
"14", "m", "<2.3", "trafford", 3, 7, 1, 0, 1,
"15", "f", ">2.3", "trafford", 0, 1, 0, 0, 0,
"16", "m", ">2.3", "trafford", 8, 6, 6, 3, 5]

# Each row is [parameter name, parameter value, stderror]
target = ["(Intercept):1"   -1.8568     0.5813;
"(Intercept):2"   -1.6115     0.5508;
"(Intercept):3"   -2.2866     0.6566;
"(Intercept):4"   -0.6642     0.3802;
"Lakegeorge:1"    -0.5753     0.7952;
"Lakegeorge:2"     1.7805     0.6232;
"Lakegeorge:3"    -1.1295     1.1928;
"Lakegeorge:4"    -0.7666     0.5686;
"Lakeoklawaha:1"  -1.1256     1.1923;
"Lakeoklawaha:2"   2.6937     0.6693;
"Lakeoklawaha:3"   1.4008     0.8105;
"Lakeoklawaha:4"  -0.7405     0.7421;
"Laketrafford:1"   0.6617     0.8461;
"Laketrafford:2"   2.9363     0.6874;
"Laketrafford:3"   1.9316     0.8253;
"Laketrafford:4"   0.7912     0.5879;
"Size>2.3:1"       0.7302     0.6523;
"Size>2.3:2"      -1.3363     0.4112;
"Size>2.3:3"       0.5570     0.6466;
"Size>2.3:4"      -0.2906     0.4599;
"Genderm:1"       -0.6064     0.6888;
"Genderm:2"       -0.4630     0.3955;
"Genderm:3"       -0.6276     0.6853;
"Genderm:4"       -0.2526     0.4663]

################################################################################
# Script

# Dummy-encode categorical data
rawdata        = reshape(rawdata, 9, 16)  # Each column is a vector of data with rownames = vcat(xnames, ylevels)
data           = construct_alligator_dataframe(rawdata, xnames, ylevels)
data.intercept = fill(1.0, nrow(data))
data.male      = [x == "m"        ? 1.0 : 0.0 for x in data.Gender]  # levels = ["f", "m"]
data.large     = [x == ">2.3"     ? 1.0 : 0.0 for x in data.Size]    # levels = ["<2.3", ">2.3"]
data.george    = [x == "george"   ? 1.0 : 0.0 for x in data.Lake]    # levels = ["hancock", "george", "oklawaha", "trafford"]
data.oklawaha  = [x == "oklawaha" ? 1.0 : 0.0 for x in data.Lake]
data.trafford  = [x == "trafford" ? 1.0 : 0.0 for x in data.Lake]

# Construct training data
ylevels_new = Dict("Fish" => 1, "Bird" => 2, "Invertebrate" => 3, "Reptile" => 4, "Other" => 5)
y = [ylevels_new[x] for x in data.StomachContents]
X = Matrix(data[:, ["intercept", "george", "oklawaha", "trafford", "large", "male"]])

# Train
fitted = fit(y, X)

# Assess fitted parameters and standard errors
target_params   = transpose(reshape(target[:, 2], 4, 6))
fitted_params   = fitted.params
target_stderror = transpose(reshape(target[:, 3], 4, 6))
fitted_stderror = fitted.stderror
@test maximum(abs.(fitted_params .- target_params)) <= 0.0001
@test maximum(abs.(fitted_stderror .- target_stderror)) <= 0.0001

end


text:
using BSON: @save, @load

using CuArrays
using Flux:
    Chain, Dense, relu, σ, softmax, logitcrossentropy, onehotbatch, flatten,
    params, gradient, gpu, cpu, loadparams!, trainmode!, testmode!,
    Conv, MaxPool
using Flux.Optimise: ADAM, update!
using Flux.Data.MNIST: images, labels
using Flux.Data: DataLoader

using Printf: @sprintf
using Parameters: @with_kw
using Plots

using Yama: create_surface, SurfaceArgs, create_trajectory

pyplot()
Plots.PyPlotBackend()

function make_minibatch(data, labels)
    batch_size = length(data)
    batch = Array{Float32}(undef, size(data[1])..., 1, batch_size)
    @inbounds for i in 1:batch_size
        batch[:, :, :, i] = Float32.(data[i])
    end

    batch, Float32.(onehotbatch(labels, 0:9))
end

@with_kw struct Args
    lr::Float32 = 3e-4
    epochs::Int32 = 20
    batch_size::Int32 = 32
    save_dir::String = "./cnn-checkpoints"
    use_gpu::Bool = true
end

function build_model(imgsize = (28, 28, 1))
    Chain(
        Conv((3, 3), 1 => 16, pad=(1, 1), relu), MaxPool((2, 2)),
        Conv((3, 3), 16 => 32, pad=(1, 1), relu), MaxPool((2, 2)),
        Conv((3, 3), 32 => 32, pad=(1, 1), relu), MaxPool((2, 2)),
        flatten,
        Dense(3 * 3 * 32, 10),
    )
end

function evaluation_step(model, batch, args::SurfaceArgs)
    x, y = make_minibatch(batch...)
    if args.use_gpu
        x, y = x |> gpu, y |> gpu
    end
    logitcrossentropy(model(x), y)
end

function save_model(model, path::String)
    checkpoint_weights = params(model) .|> cpu
    @save path checkpoint_weights
end

function train(; kws...)
    args = Args(; kws...)
    mkpath(args.save_dir)

    model = build_model()
    if args.use_gpu
        model = model |> gpu
    end
    model_parameters = params(model)

    optimizer = ADAM(args.lr)
    train_loader = DataLoader(
        images(), labels(), batchsize=args.batch_size, shuffle=false,
    )
    val_loader = DataLoader(
        images(:test), labels(:test), batchsize=args.batch_size,
    )

    save_model(model, joinpath(args.save_dir, "checkpoint-epoch=00.bson"))

    for epoch in 1:args.epochs
        trainmode!(model)
        for (i, (data, labels)) in enumerate(train_loader)
            x, y = make_minibatch(data, labels)
            if args.use_gpu
                x, y = x |> gpu, y |> gpu
            end
            gradients = gradient(model_parameters) do
                logitcrossentropy(model(x), y)
            end
            update!(optimizer, model_parameters, gradients)
        end

        testmode!(model)
        validation_loss = 0.0f0
        for (i, (data, labels)) in enumerate(val_loader)
            x, y = make_minibatch(data, labels)
            if args.use_gpu
                x, y = x |> gpu, y |> gpu
            end
            validation_loss += logitcrossentropy(model(x), y) |> cpu
        end
        validation_loss = validation_loss / length(val_loader)
        println("Epoch $epoch, Validation loss $(validation_loss)")

        save_model(
            model,
            joinpath(args.save_dir, "checkpoint-epoch=$(@sprintf("%02d", epoch)).bson"),
        )
    end
end

function main_surface()
    model = build_model()
    args = SurfaceArgs(
        batch_size=32,
        x_directions_file="./cnn-checkpoints/checkpoint-epoch=00.bson",
    )
    save_file = "./cnn-surface.bson"
    model_file = "./cnn-checkpoints/checkpoint-epoch=10.bson"

    @load model_file checkpoint_weights
    loadparams!(model, checkpoint_weights)

    # loader = DataLoader(images(:test), labels(:test), batchsize=batch_size)
    loader = DataLoader(images(), labels(), batchsize=args.batch_size)
    coordinates, loss_surface = create_surface(
        model, loader, evaluation_step, args,
    )

    @save save_file coordinates loss_surface
    # @load save_file coordinates loss_surface
    surface(coordinates..., loss_surface, linewidth=0, antialiased=false)
    gui()
end

function get_losses(position, coordinates, loss_surface)::Array{Float32, 1}
    x_coordinates, y_coordinates = coordinates
    x_nearest = argmin(abs.(x_coordinates .- position[1]))
    y_nearest = argmin(abs.(y_coordinates .- position[2]))
    [x_coordinates[x_nearest], y_coordinates[y_nearest], loss_surface[x_nearest, y_nearest]]
end

function main_trajectory()
    checkpoints = readdir("./cnn-checkpoints", join=true)
    checkpoint_file = checkpoints[end]
    checkpoints = checkpoints[1:end - 1]
    @load checkpoint_file checkpoint_weights
    target_weights = params(checkpoint_weights)

    positions = create_trajectory(target_weights, checkpoints)[2:end, :]
    positions ./= maximum(abs.(positions), dims=1)

    surface_file = "./cnn-surface.bson"
    @load surface_file coordinates loss_surface

    positions_losses = Array{Float32, 2}(undef, 3, size(positions, 1))
    for i in 1:size(positions, 1)
        positions_losses[:, i] = get_losses(
            @view(positions[i, :]), coordinates, loss_surface,
        )
    end

    surface(coordinates..., loss_surface, linewidth=0, antialiased=false)
    plot!(
        positions_losses[2, :], positions_losses[1, :], positions_losses[3, :],
        w=2, label="Optimizer path", color="white", linestyle=:dashdot,
    )
    gui()

    # contourf(coordinates..., loss_surface)
    # plot!(positions[:, 2], positions[:, 1], color="white", linestyle=:dashdot, w=2)
    # gui()
end