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

          Solving the HJB and KFE
            And then use the KFE to find the equilibrium interest rate

               Translated from matlab code on Ben Moll's website:
                              http://www.princeton.edu/~moll/HACTproject.htm

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

using LinearAlgebra, SparseArrays, Plots

γ = 2 # parameter from CRRA utility
ρ = 0.05 # the discount rate

z1=.1
z2 =.2
z= [z1 z2]

λ1 = 1.2
λ2 = 1.2

λ = [λ1 λ2]

r0=0.03 #inital guess
rmin= 0.01
rmax=0.04

H = 500 # number of points in the grid space
amin= -0.15
amax= 5

a = LinRange(amin,amax,H)
da = (amax-amin)/(H-1)

aa = [a a]
zz = ones(H,1)*z

maxit = 100 # the maximum number of iterations we allow the finite differencing algorithm
crit = 10^(-6) # our critical value
Δ = 40

dVf=zeros(H,2)
dVb = zeros(H,2)
c= zeros(H,2)

# Create the matrix that captures the systems dynamic from λ process
Aswitch = [-sparse(I,H,H)*λ[1] sparse(I,H,H)*λ[1]; sparse(I,H,H)*λ[2] -sparse(I,H,H)*λ[2]]

# Set up parameters for outer loop
Ir = 40 # number of iterations, values to check
crit_S=10^(-50)
r_grid =LinRange(rmin,rmax,Ir)

# Inital Guesses
r = r0

v0=zeros(H,2)
v0[:,1] = (z[1] .+ max(r,0.01).*a).^(1-γ)/(1-γ)/ρ
v0[:,2] = (z[2] .+ max(r,0.01).*a).^(1-γ)/(1-γ)/ρ

global v0=v0

dist=[]
V_n=[]
g_r=[]
dV_r=[]
V_r=[]
c_r=[]
adot=[]
r_r=[]
rmin_r =[]
rmax_r=[]

S =zeros(Ir,1)

# Loop that finds the equilibrium interest rate
for ir in 1:Ir
    push!(r_r,r)
    push!(rmin_r,rmin)
    push!(rmax_r,rmax)

    if ir>1
        global v0 = V_r[ir-1]
    end

    global v = v0

    # The finite differeing loop
    for n in 1:maxit
        global V=v
        push!(V_n,V)
        # Forward differencing
        dVf[1:H-1,:]=(V[2:H,:]-V[1:H-1,:])/da
        dVf[H,:] = (z .+ r.*amax).^(-γ) # impose state constraint at the max, just in case
        # Backward differencing
        dVb[2:H,:]=(V[2:H,:]-V[1:H-1,:])/da
        dVb[1,:]= (z .+r.*amin).^(-γ)

        # Find consumption and savings with the forward difference
        cf = dVf.^(-1/γ)
        sf = zz +r.*aa -cf

        # Find consumption and savings with the backward difference
        cb= dVb.^(-1/γ)
        sb = zz +r.*aa -cb

        # Find consumption and savings at steady state
        c0 = zz + r.*aa
        dV0 = c0.^(-γ)

        # Now implement the upwind scheme in order to select the best differencing method
        If = sf .> 0
        Ib = sb .< 0
        I0 = (1 .- If - Ib)

        # State constraint at amin is automatically implemented

        global dV_Upwind = dVf.*If + dVb.*Ib + dV0.*I0
        global c = dV_Upwind.^(-1/γ)

        u = c.^(1-γ)/(1-γ)
        # Construct matrix for the evolution of the system
        X = -min.(sb,0)/da
        Y= -max.(sf,0)/da + min.(sb,0)/da
        Z = max.(sf,0)/da

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

        global A = AA + Aswitch

        B = (ρ +1/Δ)*sparse(I,2*H,2*H) -A
        u_stacked = [u[:,1]; u[:,2]]
        V_stacked = [V[:,1]; V[:,2]]

        b = u_stacked + V_stacked/Δ
        V_stacked = B\b # Solves the system of equations

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

        V_change = V - v

        global v = V

        push!(dist, findmax(abs.(V_change))[1])

        if dist[n] < crit
            println("Value Function converged, Iteration=")
            println(n)
            break
        end
    end

    #====================================================
            Fokker-Planck Equation
    ====================================================#

    AT = A' # we need the transpose of our A matrix
    b =zeros(2*H,1)

    # We need to fix one value to prevent the matrix from being singular
    i_fix =1
    b[i_fix] = .1
    row = [zeros(1,i_fix-1) 1 zeros(1,2*H-i_fix)]
    AT[i_fix,:] = row

    # Using this solve the linear system
    gg = AT\b
    g_sum = gg'*ones(2*H,1)*da
    gg = gg./g_sum # we need gg to sum to one

    g = [gg[1:H] gg[H+1:2*H]]

    # Graphs
    push!(g_r, g)
    push!(adot, zz+r.*aa-c)
    push!(V_r, V)
    push!(c_r, c)
    push!(dV_r, dV_Upwind)

    S[ir] = g[:,1]'*a*da + g[:,2]'*a*da

    if S[ir]>crit_S
        println("Excess Supply")
        global rmax =r
        global r = 0.5*(r+rmin)
    elseif S[ir]<-crit_S
        println("Excess Demand")
        global rmin=r
        global r = 0.5*(r+rmax)
    elseif abs(S[ir])<crit_S
        println("Equilibrium Found")
        println(r)
        break
    end

#end of outer loop
end

# Graphs
amax1 = 0.6
amin1 = amin-0.03

plot(a,adot[end][:,1], label="\$ s_1(a)\$", xlabel="Wealth", ylabel="\$s_i(a)\$", xlims=(amin1,amax1))
plot!(a,adot[end][:,2], label="\$ s_2(a)\$")
plot!(LinRange(amin1,amax1,H),zeros(H,1), label="", line=:dash, color=:black)

plot(a,g_r[end][:,1], label="\$ g_1(a)\$", xlabel="Wealth", ylabel="\$g_i(a)\$", xlims=(amin1,amax1), ylims=(0,3))
plot!(a,g_r[end][:,2], label="\$ g_2(a)\$")


text:
# from: https://juliaparallel.github.io/MPI.jl/stable/examples/01-hello/
# run: ~/.julia/bin/mpiexecjl -n 4 julia --project scripts/hello_mpi.jl
using MPI
MPI.Init()

comm = MPI.COMM_WORLD
println("Hello world, I am $(MPI.Comm_rank(comm)) of $(MPI.Comm_size(comm))")
MPI.Barrier(comm)

ext: jl


text:
# globals
global use_cache = true

const cache_loc = joinpath(dirname(@__FILE__), "../.cache")

function cache_init()
    !isdir(cache_loc) && mkdir(cache_loc)
end

"""
Allows the caching mechanism in Descartes to be disabled (`Cache.enabled(false)`)
and enabled(`Cache.enabled(false)`).
"""
function cache_enabled(t::Bool)
    global use_cache = t
end

"""
Get the JLD cache for a primitive.
"""
function cache_location(a::AbstractPrimitive)
    hn = hash(a)
    joinpath(cache_loc, "$hn.jld")
end

"""
Load an object from the primitive cache.
"""
function cache_load(a::AbstractPrimitive, key)
    fl = cache_location(a)
    jldopen(fl, "r") do f
        read(f, key)
    end
end

"""
Add `obj` under `key` to the primitive's cache.
"""
function cache_add(a::AbstractPrimitive, key, obj)
    fl = cache_location(a)
    !isfile(fl) && touch(fl)
    jldopen(fl, "w") do f
        write(f, key, obj)
    end
end

"""
Check if the cache contains the primitive and key
"""
function cache_contains(a::AbstractPrimitive, key)
    fl = cache_location(a)
    if isfile(fl)
        # check if the key is in there
        return jldopen(fl) do f
                   key in names(f)
               end
    else
        return false
    end
end


text:
function FaradayStep!(b, e, m, cb, nx, ny, nz)
    @turbo for k = 1:nz-1
        for j = 1:ny-1
            for i = 1:nx-1
                b[i, j, k, 1] =
                    b[i, j, k, 1] + cb[i, j, k, 1] * (e[i, j, k, 3] - e[i, j+1, k, 3] + e[i, j, k+1, 2] - e[i, j, k, 2]) - m[i, j, k, 1]
                b[i, j, k, 2] =
                    b[i, j, k, 2] + cb[i, j, k, 2] * (e[i, j, k, 1] - e[i, j, k+1, 1] + e[i+1, j, k, 3] - e[i, j, k, 3]) - m[i, j, k, 2]
                b[i, j, k, 3] =
                    b[i, j, k, 3] + cb[i, j, k, 3] * (e[i, j, k, 2] - e[i+1, j, k, 2] + e[i, j+1, k, 1] - e[i, j, k, 1]) - m[i, j, k, 3]
            end
        end
    end
    @turbo for k = 1:nz-1
        for j = 1:ny-1
            for i = nx:nx
                b[i, j, k, 1] =
                    b[i, j, k, 1] + cb[i, j, k, 1] * (e[i, j, k, 3] - e[i, j+1, k, 3] + e[i, j, k+1, 2] - e[i, j, k, 2]) - m[i, j, k, 1]
            end
        end
    end
    @turbo for k = 1:nz-1
        for j = ny:ny
            for i = 1:nx-1
                b[i, j, k, 2] =
                    b[i, j, k, 2] + cb[i, j, k, 2] * (e[i, j, k, 1] - e[i, j, k+1, 1] + e[i+1, j, k, 3] - e[i, j, k, 3]) - m[i, j, k, 2]
            end
        end
    end
    @turbo for k = nz:nz
        for j = 1:ny-1
            for i = 1:nx-1
                b[i, j, k, 3] =
                    b[i, j, k, 3] + cb[i, j, k, 3] * (e[i, j, k, 2] - e[i+1, j, k, 2] + e[i, j+1, k, 1] - e[i, j, k, 1]) - m[i, j, k, 3]
            end
        end
    end
end

            #    ex hx ey hy ez hz
function FaradayStep!(f, m, cb, nx, ny, nz)
    @turbo for k = 1:nz-1
        for j = 1:ny-1
            for i = 1:nx-1
                f[i, j, k, 2] =
                    f[i, j, k, 2] + cb[i, j, k, 1] * (f[i, j, k, 5] - f[i, j+1, k, 5] + f[i, j, k+1, 3] - f[i, j, k, 3]) - m[i, j, k, 1]
                f[i, j, k, 4] =
                    f[i, j, k, 4] + cb[i, j, k, 2] * (f[i, j, k, 1] - f[i, j, k+1, 1] + f[i+1, j, k, 5] - f[i, j, k, 5]) - m[i, j, k, 2]
                f[i, j, k, 6] =
                    f[i, j, k, 6] + cb[i, j, k, 3] * (f[i, j, k, 3] - f[i+1, j, k, 3] + f[i, j+1, k, 1] - f[i, j, k, 1]) - m[i, j, k, 3]
            end
        end
    end
    @turbo for k = 1:nz-1
        for j = 1:ny-1
            for i = nx:nx
                f[i, j, k, 4] =
                f[i, j, k, 4] + cb[i, j, k, 1] * (f[i, j, k, 3] - f[i, j+1, k, 3] + f[i, j, k+1, 2] - f[i, j, k, 2]) - m[i, j, k, 1]
            end
        end
    end
    @turbo for k = 1:nz-1
        for j = ny:ny
            for i = 1:nx-1
                f[i, j, k, 5] =
                f[i, j, k, 5] + cb[i, j, k, 2] * (f[i, j, k, 1] - f[i, j, k+1, 1] + f[i+1, j, k, 3] - f[i, j, k, 3]) - m[i, j, k, 2]
            end
        end
    end
    @turbo for k = nz:nz
        for j = 1:ny-1
            for i = 1:nx-1
                f[i, j, k, 6] =
                f[i, j, k, 6] + cb[i, j, k, 3] * (f[i, j, k, 2] - f[i+1, j, k, 2] + f[i, j+1, k, 1] - f[i, j, k, 1]) - m[i, j, k, 3]
            end
        end
    end
end


text:
# Indexing at a point
@inline function (itp::GriddedInterpolation{T,N})(x::Vararg{Number,N}) where {T,N}
    @boundscheck (checkbounds(Bool, itp, x...) || Base.throw_boundserror(itp, x))
    wis = weightedindexes((value_weights,), itpinfo(itp)..., x)
    coefficients(itp)[wis...]
end
@inline function (itp::GriddedInterpolation)(x::Vararg{UnexpandedIndexTypes})
    itp(to_indices(itp, x)...)
end

@inline function gradient(itp::GriddedInterpolation{T,N}, x::Vararg{Number,N}) where {T,N}
    @boundscheck (checkbounds(Bool, itp, x...) || Base.throw_boundserror(itp, x))
    wis = weightedindexes((value_weights, gradient_weights), itpinfo(itp)..., x)
    SVector(map(inds->coefficients(itp)[inds...], wis))
end

itpinfo(itp::GriddedInterpolation) = (tcollect(itpflag, itp), itp.knots)

weightedindex_parts(fs::F, itpflag::Gridded, ax, x) where F =
    weightedindex_parts(fs, degree(itpflag), ax, x)

roundbounds(x::Integer, knotvec::AbstractVector) = gridded_roundbounds(x, knotvec)
roundbounds(x::Number, knotvec::AbstractVector) = gridded_roundbounds(x, knotvec)
function gridded_roundbounds(x, knotvec::AbstractVector)
    i = find_knot_index(knotvec, x)
    iclamp = max(i, first(axes1(knotvec)))
    inext = min(iclamp+1, last(axes1(knotvec)))
    ifelse(i < iclamp, i+1, ifelse(x - knotvec[iclamp] < knotvec[inext] - x, i, inext))
end

floorbounds(x::Integer, knotvec::AbstractVector) = gridded_floorbounds(x, knotvec)
floorbounds(x, knotvec::AbstractVector) = gridded_floorbounds(x, knotvec)
function gridded_floorbounds(x, knotvec::AbstractVector)
    i = find_knot_index(knotvec, x)
    max(i, first(axes1(knotvec)))
end

@inline find_knot_index(knotv, x) = searchsortedfirst(knotv, x, first(axes1(knotv)), length(knotv), Base.Order.ForwardOrdering()) - 1

function weightedindex_parts(fs::F, deg::Degree, knotvec::AbstractVector, x) where F
    i = find_knot_index(knotvec, x)
    ax1 = axes1(knotvec)
    iclamp = clamp(i, first(ax1), last(ax1)-1)
    weightedindex(fs, deg, knotvec, x, iclamp)
end

function weightedindex(fs::F, deg::Constant, knotvec, x, iclamp) where F
    pos, δx = positions(deg, knotvec,  x)
    (position=pos, coefs=fmap(fs, deg, δx))
end
function weightedindex(fs::F, deg::Degree, knotvec, x, iclamp) where F
    @inbounds l, u = knotvec[iclamp], knotvec[iclamp+1]
    δx = (x - l)/(u - l)
    (position=iclamp, coefs=rescale_gridded(fs, fmap(fs, deg, δx), u-l))
end

rescale_gridded(fs::F, coefs, Δx) where F =
    (rescale_gridded(fs[1], coefs[1], Δx), rescale_gridded(Base.tail(fs), Base.tail(coefs), Δx)...)
rescale_gridded(::Tuple{}, ::Tuple{}, Δx) = ()
rescale_gridded(::typeof(value_weights), coefs, Δx) = coefs
rescale_gridded(::typeof(gradient_weights), coefs, Δx) = coefs./Δx
rescale_gridded(::typeof(hessian_weights), coefs, Δx) = coefs./Δx.^2

@inline function (itp::GriddedInterpolation{T,N})(x::Vararg{Union{Number,AbstractVector},N}) where {T,N}
    @boundscheck (checkbounds(Bool, itp, x...) || Base.throw_boundserror(itp, x))
    itps = tcollect(itpflag, itp)
    wis = dimension_wis(value_weights, itps, itp.knots, x)
    coefs = coefficients(itp)
    ret = [coefs[i...] for i in Iterators.product(wis...)]
    reshape(ret, shape(wis...))
end

function dimension_wis(f::F, itps, knots, xs) where F
    itpflag, knotvec, x = itps[1], knots[1], xs[1]
    function makewi(y)
        pos, coefs = weightedindex_parts((f,), itpflag, knotvec, y)
        maybe_weightedindex(pos, coefs[1])
    end
    (makewi.(x), dimension_wis(f, Base.tail(itps), Base.tail(knots), Base.tail(xs))...)
end
function dimension_wis(f::F, itps::Tuple{NoInterp,Vararg{Any}}, knots, xs) where F
    (Int.(xs[1]), dimension_wis(f, Base.tail(itps), Base.tail(knots), Base.tail(xs))...)
end
dimension_wis(f, ::Tuple{}, ::Tuple{}, ::Tuple{}) = ()

function getindex_return_type(::Type{GriddedInterpolation{T,N,TCoefs,IT,K}}, argtypes) where {T,N,TCoefs,IT<:DimSpec{Gridded},K}
    Tret = TCoefs
    for a in argtypes
        Tret = Base.promote_op(*, Tret, a)
    end
    Tret
end


text:
# Testing constant non-deterministic input
c = zeros(4); r = 0.1
U = ConstantInput(BallInf(c, r))

#@test length(U) == 1 # constant input => infinitely constant, length not implemented
inputs = next_set(U)
input_remains_constant = inputs.center == c && inputs.radius == r
@test input_remains_constant
for i in 1:5
    inputs = next(U, i)[1]
    input_remains_constant = input_remains_constant && inputs.center == c && inputs.radius == r
end

# Testing time-varying non-deterministic input
c = [zeros(4), 1. + ones(4), 2. + ones(4)]
r = [0.1, 0.2, 0.3]

U = VaryingInput([BallInf(c[i], r[i]) for i in 1:3])

@test length(U) == 3
input_is_varying = true
i = 1
for u in U
    input_is_varying = input_is_varying && u.center == c[i] && u.radius == r[i]
    i += 1
end
@test input_is_varying


text:
import Nemo.setcoeff!, Nemo.exp, Base.start, Base.next, Base.done, Nemo.lift, Hecke.lift
export start, next, done, SetPrimes, psi_lower, psi_upper

#function setcoeff!(g::fmpz_mod_rel_series, i::Int64, a::Nemo.GenRes{Nemo.fmpz})
#  setcoeff!(g, i, lift(a))
#end

function bernstein(h::Int, it::Any, Q = FlintQQ, cl = ceil, a::Int = 776)
  #use 771 and cl = floor to get decent upper bounds
  # more on the choice of 776 and 771 in Dan's paper

  #println("in bernstein with a=$a and cl(0.5) = $(cl(0.5))")
  R,t = PowerSeriesRing(Q, a*h+1, "t")

  #implements https://cr.yp.to/papers/psi.pdf
  # for smoothness for ideals, replace next_prime by the list of the norms of the prime
  # ideals
  # the sums of the coefficients of exp(g) are bounds for psi

  st = start(it)
  p, st = next(it, st)
  g = R(0)
  tp = R(0)
  lpp = Int(cl(Float64(log(p))/log(2)*a))

  function do_single!(g::SeriesElem, pp::Int, np::Int)
    i = pp
    nu = Q(np)*Q(pp)
    while i <= a*h
      A = coeff(g, i)
      A += divexact(nu, Q(i))
      setcoeff!(g, i, A)
      i += pp
    end
  end

  while true
    pp = lpp
    np = 0
    while pp == lpp && !done(it, st)
      np += 1
      p, st = next(it, st)
      lpp = Int(cl(Float64(log(p))/log(2)*a))
      if done(it, st)
        break
      end
    end

    if done(it, st) && pp == lpp
      np += 1
    end

    do_single!(g, pp, np)

    if done(it, st)
      if pp != lpp
        do_single!(g, lpp, 1)
      end
      return g
    end
  end
end

function _exp(a::fmpz_mod_rel_series)
  R = base_ring(parent(a))
  Rx,x = PolynomialRing(R)
  A = Rx()
  for i=0:length(a)-1
    setcoeff!(A, i, coeff(a, i))
  end
  E = Rx()
  ccall((:nmod_poly_exp_series, :libflint), Void, (Ptr{nmod_poly}, Ptr{nmod_poly}, Int64), &E, &A, length(a))
  r = parent(a)()
  for i=0:length(E)-1
    setcoeff!(r, i, coeff(E, i))
  end
  return r
end

immutable SetPrimes{T}
  from::T
  to::T
  mod::T # if set (i.e. >1), only primes p % mod == a are returned
  a::T
  function SetPrimes(f::T, t::T)
    r = new(f, t, T(1), T(0))
    return r
  end
  function SetPrimes(f::T, t::T, mod::T, val::T)
    r = new(f, t, mod, val)
    if gcd(mod, val) != 1
      error("modulus and value need to be coprime")
    end
    return r
  end
end

function SetPrimes{T}(f::T, t::T)
  return SetPrimes{T}(f, t)
end

function SetPrimes{T}(f::T, t::T, mod::T, val::T)
  return SetPrimes{T}(f, t, mod, val)
end

function start{T<: Union{Integer, fmpz}}(A::SetPrimes{T})
  curr = A.from
  if A.mod >1 && curr % A.mod != A.a
    curr += A.mod - (curr % A.mod) + A.a
  end
  while !isprime(curr)
    curr += A.mod
  end
  return curr
end

function next{T<: Union{Integer, fmpz}}(A::SetPrimes{T}, st::T)
  p = st
  if A.mod >1
    m = A.mod
  else
    if p==2
      st = T(3)
      return p, st
    end
    m = T(2)
  end
  st = p+m
  while !isprime(st)
    st += m
  end
  return p, st
end

function done{T <: Union{Integer, fmpz}}(A::SetPrimes{T}, st::T)
  return A.to != -1 && st > A.to
end

eltype{T <: Union{Integer, fmpz}}(::SetPrimes{T}) = T

function lift(R::FmpzRelSeriesRing, f::fmpz_mod_rel_series)
  r = R()
  for i=0:length(f)-1
    setcoeff!(r, i, lift(coeff(f, i)))
  end
  return r
end

function psi_lower(N::fmpz, pr, a::Int=776, cl = ceil)
  p = fmpz(next_prime(2^60))
  n = Int(ceil(log(N)/log(2)))
#  println("precision of $n")
  f = _exp(bernstein(n, pr, ResidueRing(FlintZZ, p), cl, a))
  Rt, t = PowerSeriesRing(FlintZZ, n*a+1, "t")
  f = lift(Rt, f)
  pp = p
  while pp < N
    p = next_prime(p)
#    println("p: $p, pp: $pp N:$N")
    g = _exp(bernstein(n, pr, ResidueRing(FlintZZ, p), cl, a))
    @assert length(g) == length(f)
    for i=0:length(f)
      setcoeff!(f, i, crt(coeff(f, i), pp, lift(coeff(g, i)), p))
    end
    pp *= p
  end
  res = []
  s = 0
  j = 0
  for i=0:n
    while j <= i*a && j < length(f)
      s += coeff(f, j)
      j += 1
    end
    push!(res, s)
    if j >= length(f) break; end
  end
  return res, f  # res[i] <= psi(2^(i-1), B)
end

function psi_lower(N::fmpz, B::Int, a::Int = 776, cl = ceil)
  return psi_lower(fmpz(N), SetPrimes{Int}(2, B), a, cl)
end

function psi_lower(N::Integer, B::Int, a::Int = 776, cl = ceil)
  return psi_lower(fmpz(N), SetPrimes{Int}(2, B), a, cl)
end

function psi_upper(N::fmpz, B::Int, a::Int=771, fl = floor)
  return psi_lower(N, SetPrimes{Int}(2, B), a, fl)
end

function psi_upper(N::Integer, B::Int, a::Int=771, fl = floor)
  return psi_lower(fmpz(N), SetPrimes{Int}(2, B), a, fl)
end

function psi_lower(N::Integer, B::Hecke.NfFactorBase, a::Int=776, cl = ceil)
  lp = sort(fmpz[norm(x) for x=B.ideals])
  return psi_lower(fmpz(N), lp, a, cl)
end

function psi_lower(N::fmpz, B::Hecke.NfFactorBase, a::Int=776, cl = ceil)
  lp = sort(fmpz[norm(x) for x=B.ideals])
  return psi_lower(N, lp, a, cl)
end

function psi_upper(N::Integer, B::Hecke.NfFactorBase, a::Int=771, cl = floor)
  lp = sort(fmpz[norm(x) for x=B.ideals])
  return psi_lower(fmpz(N), lp, a, cl)
end

function psi_upper(N::fmpz, B::Hecke.NfFactorBase, a::Int=771, cl = floor)
  lp = sort(fmpz[norm(x) for x=B.ideals])
  return psi_lower(N, lp, a, cl)
end


text:
module adapters_hive_repo_test

using Test # @test
using Octo.Adapters.Hive # Repo Schema Raw USE

Repo.debug_sql()

Repo.connect(adapter=Octo.Adapters.Hive, host="localhost", port=10000)

struct Employee
end
Schema.model(Employee, table_name="employee", primary_key="eid")

Repo.execute([DROP TABLE IF EXISTS Employee])
Repo.execute(Raw("""
CREATE TABLE IF NOT EXISTS employee (eid Int, name String, salary Double)
"""))

Repo.insert!(Employee, [(1, "Jeremy", 10000.50), (2, "Cloris", 20000.50)])

df = Repo.query(Employee)

@test df[1][Symbol("employee.name")] == "Jeremy"
@test df[2][Symbol("employee.name")] == "Cloris"

@test df[1][Symbol("employee.salary")] == 10000.50
@test df[2][Symbol("employee.salary")] == 20000.50

@test size(df) == (2,)

Repo.disconnect()

end # module adapters_hive_repo_test


text:
import Gaugefields:comm,setvalue!

"""
Struct for DomainwallFermion
"""
struct DomainwallFermion_5D_wing_mpi{NC,WilsonFermion} <: Abstract_DomainwallFermion_5D{NC,WilsonFermion} #<: AbstractFermionfields_5D{NC}
    w::Array{WilsonFermion,1}
    NC::Int64
    NX::Int64
    NY::Int64
    NZ::Int64
    NT::Int64
    L5::Int64
    Dirac_operator::String
    NWilson::Int64
    PEs::NTuple{4,Int64}
    PN::NTuple{4,Int64}
    mpiinit::Bool
    myrank::Int64
    nprocs::Int64
    myrank_xyzt::NTuple{4,Int64}
    mpi::Bool

    function DomainwallFermion_5D_wing_mpi(L5,NC::T,NX::T,NY::T,NZ::T,NT::T,PEs) where T<: Integer

        NV = NX*NY*NZ*NT
        @assert NX % PEs[1] == 0 "NX % PEs[1] should be 0. Now NX = $NX and PEs = $PEs"
        @assert NY % PEs[2] == 0 "NY % PEs[2] should be 0. Now NY = $NY and PEs = $PEs"
        @assert NZ % PEs[3] == 0 "NZ % PEs[3] should be 0. Now NZ = $NZ and PEs = $PEs"
        @assert NT % PEs[4] == 0 "NT % PEs[4] should be 0. Now NT = $NT and PEs = $PEs"

        PN = (NX ÷ PEs[1],
                    NY ÷ PEs[2],
                    NZ ÷ PEs[3],
                    NT ÷ PEs[4],
            )

        nprocs = MPI.Comm_size(comm)
        @assert prod(PEs) == nprocs "num. of MPI process should be prod(PEs). Now nprocs = $nprocs and PEs = $PEs"
        myrank = MPI.Comm_rank(comm)

        myrank_xyzt = get_myrank_xyzt(myrank,PEs)

        x = WilsonFermion_4D_mpi(NC,NX,NY,NZ,NT,PEs)
        xtype = typeof(x)
        w = Array{xtype,1}(undef,L5)
        w[1] = x
        for i=2:L5
            w[i] = similar(x)
        end
        #println(w[2][1,1,1,1,1,1])
        NWilson = length(x)
        Dirac_operator = "Domainwall"
        mpi = true
        mpiinit = true
        return new{NC,xtype}(w,NC,NX,NY,NZ,NT,L5,Dirac_operator,NWilson,
                Tuple(PEs),PN,mpiinit,myrank,nprocs,myrank_xyzt,mpi)
    end

end

function Base.similar(x::DomainwallFermion_5D_wing_mpi{NC,WilsonFermion} ) where {NC,WilsonFermion}
    return DomainwallFermion_5D_wing_mpi(x.L5,NC,x.NX,x.NY,x.NZ,x.NT,x.PEs)
end

#=

function D5DWx!(xout::DomainwallFermion_5D_wing_mpi{NC,WilsonFermion} ,U::Array{G,1},
    x::DomainwallFermion_5D_wing_mpi{NC,WilsonFermion} ,m,A,L5) where  {NC,WilsonFermion,G <: AbstractGaugefields}

    #temp = temps[4]
    #temp1 = temps[1]
    #temp2 = temps[2]
    clear_fermion!(xout)
    ratio = 1
    #ratio = xout.L5/L5
    if L5 != xout.L5
        @assert L5 % 2 == 0
        irange = Int64[]
        irange_out = Int64[]
        #irange = 1:L5
        #irange_out = (L5+1):xout.L5

        for i5=1:xout.L5
            if i5 <= div(L5,2) || i5 >= xout.L5-div(L5,2)+1
                push!(irange,i5)
            else
                push!(irange_out,i5)
            end

        end

        #for i5 in irange_out
        #    axpy!(1,x.w[i5],xout.w[i5])
        #end
    else
        irange = 1:L5
    end

    for i5 in irange
        j5=i5
        D4x!(xout.w[i5],U,x.w[j5],A,4) #Dw*x
        #Dx!(xout.w[i5],U,x.w[j5],A) #Dw*x
        #Wx!(xout.w[i5],U,x.w[j5],temps) #Dw*x
        #1/(2*A.κ)
        massfactor = -(1/(2*A.κ) + 1)
        set_wing_fermion!(xout.w[i5])
        #add!(ratio,xout.w[i5],ratio,x.w[j5]) #D = x + Ddagw*x
        add!(ratio,xout.w[i5],ratio*massfactor,x.w[j5]) #D = x + Dw*x
        set_wing_fermion!(xout.w[i5])

        #println("xout ",xout.w[i5][1,1,1,1,1,1])

        j5=i5+1
        if 1 <= j5 <= xout.L5
            #-P_- -> - P_+ :gamma_5 of LTK definition
            if xout.L5 != 2
                #mul_1minusγ5x_add!(xout.w[i5],x.w[j5],-1*ratio)
                mul_1plusγ5x_add!(xout.w[i5],x.w[j5],ratio)
                set_wing_fermion!(xout.w[i5])

            end
        end

        j5=i5-1
        if 1 <= j5 <= xout.L5
            #-P_+ -> - P_- :gamma_5 of LTK definition
            if xout.L5 != 2
                #mul_1plusγ5x_add!(xout.w[i5],x.w[j5],-1*ratio)
                mul_1minusγ5x_add!(xout.w[i5],x.w[j5],ratio)
                set_wing_fermion!(xout.w[i5])
            end
        end

        if xout.L5 != 1
            if i5==1
                j5 = xout.L5
                #mul_1plusγ5x_add!(xout.w[i5],x.w[j5],m*ratio)
                mul_1minusγ5x_add!(xout.w[i5],x.w[j5],-m*ratio)
                set_wing_fermion!(xout.w[i5])
            end

            if i5== xout.L5
                j5 = 1
                #mul_1minusγ5x_add!(xout.w[i5],x.w[j5],m*ratio)
                mul_1plusγ5x_add!(xout.w[i5],x.w[j5],-m*ratio)
                set_wing_fermion!(xout.w[i5])
            end
        end

        #println("xout ",xout.w[i5][1,1,1,1,1,1])

    end

    if L5 != xout.L5
        for i5 in irange_out
            axpy!(1,x.w[i5],xout.w[i5])
        end
    end

    set_wing_fermion!(xout)

    return
end

function D5DWdagx!(xout::DomainwallFermion_5D_wing_mpi{NC,WilsonFermion} ,U::Array{G,1},
    x::DomainwallFermion_5D_wing_mpi{NC,WilsonFermion} ,m,A,L5) where  {NC,WilsonFermion,G <: AbstractGaugefields}

    #temp = temps[4]
    #temp1 = temps[1]
    #temp2 = temps[2]
    clear_fermion!(xout)
    ratio = 1
    #ratio = xout.L5/L5

    if L5 != xout.L5
        @assert L5 % 2 == 0
        irange = Int64[]
        irange_out = Int64[]
        #irange = 1:L5
        #irange_out = (L5+1):xout.L5

        for i5=1:xout.L5
            if i5 <= div(L5,2) || i5 >= xout.L5-div(L5,2)+1
                push!(irange,i5)
            else
                push!(irange_out,i5)
            end

        end

        #for i5 in irange_out
        #    axpy!(1,x.w[i5],xout.w[i5])
        #end
    else
        irange = 1:L5
    end

    for i5 in irange
        j5=i5
        #Ddagx!(xout.w[i5],U,x.w[j5],A) #Ddagw*x
        D4dagx!(xout.w[i5],U,x.w[j5],A,4) #Dw*x
        #Wx!(xout.w[i5],U,x.w[j5],temps) #Dw*x
        #1/(2*A.κ)
        massfactor = -(1/(2*A.κ) + 1)
        #println(massfactor)

        #Wdagx!(xout.w[i5],U,x.w[j5],temps) #Ddagw*x
        set_wing_fermion!(xout.w[i5])
        add!(ratio,xout.w[i5],ratio*massfactor,x.w[j5]) #D = x + Dw*x
        #add!(ratio,xout.w[i5],ratio,x.w[j5]) #D = x + Ddagw*x
        set_wing_fermion!(xout.w[i5])

        j5=i5+1
        if 1 <= j5 <= xout.L5
            #-P_-
            if xout.L5 != 2
                #mul_1plusγ5x_add!(xout.w[i5],x.w[j5],-1*ratio)
                mul_1minusγ5x_add!(xout.w[i5],x.w[j5],ratio)
                set_wing_fermion!(xout.w[i5])
            end
        end

        j5=i5-1
        if 1 <= j5 <= xout.L5
            #-P_+
            if xout.L5 != 2
                #mul_1minusγ5x_add!(xout.w[i5],x.w[j5],-1*ratio)
                mul_1plusγ5x_add!(xout.w[i5],x.w[j5],ratio)
                set_wing_fermion!(xout.w[i5])
            end
        end

        if L5 != 1
            if i5==1
                j5 = xout.L5
                #mul_1minusγ5x_add!(xout.w[i5],x.w[j5],m*ratio)
                mul_1plusγ5x_add!(xout.w[i5],x.w[j5],-m*ratio)
                set_wing_fermion!(xout.w[i5])
            end

            if i5==xout.L5
                j5 = 1
                #mul_1plusγ5x_add!(xout.w[i5],x.w[j5],m*ratio)
                mul_1minusγ5x_add!(xout.w[i5],x.w[j5],-m*ratio)
                set_wing_fermion!(xout.w[i5])
            end
        end

    end

    #if L5 != xout.L5
    #    for i5=L5+1:xout.L5
    #        axpy!(1,x.w[i5],xout.w[i5])
    #    end
    #end

    if L5 != xout.L5
        for i5 in irange_out
            axpy!(1,x.w[i5],xout.w[i5])
        end
    end

    set_wing_fermion!(xout)

    return
end

=#

#=
"""
c-------------------------------------------------c
c     Random number function for Gaussian  Noise
    with σ^2 = 1/2
c-------------------------------------------------c
    """
function gauss_distribution_fermion!(x::DomainwallFermion_5D_wing_mpi{NC,NDW}) where {NC,NDW}
    NX = x.NX
    NY = x.NY
    NZ = x.NZ
    NT = x.NT
    n6 = size(x.w[1].f)[6]
    σ = sqrt(1/2)

    L5 = length(x.w)
    for iL = 1:L5
        for ialpha = 1:n6
            for it=1:x.PN[4]
                for iz=1:x.PN[3]
                    for iy=1:x.PN[2]
                        for ix=1:x.PN[1]
                            for ic=1:NC
                                v = σ*randn()+im*σ*randn()
                                setvalue!(x.w[iL],v,ic,ix,iy,iz,it,ialpha)
                                #x[ic,ix,iy,iz,it,ialpha] = σ*randn()+im*σ*randn()
                            end
                        end
                    end
                end
            end
        end
        set_wing_fermion!(x.w[iL])
    end
    return
end

"""
c-------------------------------------------------c
c     Random number function for Gaussian  Noise
    with σ^2 = 1/2
c-------------------------------------------------c
    """
function gauss_distribution_fermion!(x::DomainwallFermion_5D_wing_mpi{NC,NDW},randomfunc,σ) where {NC,NDW}

    NX = x.NX
    NY = x.NY
    NZ = x.NZ
    NT = x.NT
    n6 = size(x.w[1].f)[6]
    #σ = sqrt(1/2)

    L5 = length(x.w)
    for iL = 1:L5

        for mu = 1:n6
            for ic=1:NC
                for it=1:x.PN[4]
                    for iz=1:x.PN[3]
                        for iy=1:x.PN[2]
                            for ix=1:x.PN[1]
                                v1 = sqrt(-log(randomfunc()+1e-10))
                                v2 = 2pi*randomfunc()

                                xr = v1*cos(v2)
                                xi = v1 * sin(v2)

                                v = σ*xr + σ*im*xi
                                setvalue!(x.w[iL],v,ic,ix,iy,iz,it,mu)

                                #x[ic,ix,iy,iz,it,mu] = σ*xr + σ*im*xi
                            end
                        end
                    end
                end
            end
        end

        set_wing_fermion!(x.w[iL])
    end

    return
end

=#


text:
"""
    struct CRMomentumDistributionConfig{T1,T2}
        pmin::T1
        pmax::T1
        Nbins::T2
        bin_width::T1
        mp::T1
        me::T1
        c::T1
        mc_e::T1
        mc_p::T1
    end

Config object to obtain distribution spectrum in momentum space.
"""
struct CRMomentumDistributionConfig
    pmin::Float64
    pmax::Float64
    Nbins::Int64
    bin_width::Float64
    mp::Float64
    me::Float64
    c::Float64
    mc_e::Float64
    mc_p::Float64

    function CRMomentumDistributionConfig(pmin::Real = 0.0, pmax::Real = 0.0,
        Nbins::Integer = 24, mode::Integer = 3)

        # the original version used cgs units
        if mode == 1
            CNST_ME = 9.1095e-28
            CNST_MP = 1.6726e-24
            CNST_C = 2.9979e10
        else
            CNST_ME = 1.0
            CNST_MP = 1.0
            CNST_C = 1.0
        end

        MCe = CNST_ME * CNST_C
        MCp = CNST_MP * CNST_C
        bin_width = log10(pmax / pmin) / Nbins

        new(pmin, pmax, Nbins, bin_width, CNST_MP, CNST_ME, CNST_C, MCe, MCp)

    end
end

"""
    struct CRMomentumDistribution
        bound::Vector{Float64}
        norm::Vector{Float64}
    end

Distribution function reconstruction.
"""
struct CRMomentumDistribution
    bound::Vector{Float64}
    norm::Vector{Float64}
end

"""
    CRMomentumDistribution(CR_N::Vector{<:Real}, CR_S::Vector{<:Real}, CR_C::Real,
                                    pmin::Real, pmax::Real, mc_p::Real)

Function to construct the `CRMomentumDistribution` via properties from a snapshot file.
"""
function CRMomentumDistribution(CR_N::Vector{<:Real}, CR_S::Vector{<:Real}, CR_C::Real,
    pmin::Real, pmax::Real, mc_p::Real, mode::Integer = 3)

    bound, norm = norm_spectrum(CR_N, CR_S, CR_C, pmin, pmax, mc_p, mode)
    CRMomentumDistribution(bound, norm)
end

"""
    struct CR_NormSpectrum
        bin::Vector{Float64}
        energy::Vector{Float64}
    end

Stores the norm spectrum.
"""
struct CR_NormSpectrum
    bin::Vector{Float64}
    norm::Vector{Float64}
end

"""
    struct CR_EnergySpectrum
        bin::Vector{Float64}
        energy::Vector{Float64}
    end

Stores the energy spectrum.
"""
struct CR_EnergySpectrum
    bin::Vector{Float64}
    energy::Vector{Float64}
end

"""
    struct CR_NumSpectrum
        bin::Vector{Float64}
        density::Vector{Float64}
    end

Stores the number density spectrum.
"""
struct CR_NumSpectrum
    bin::Vector{Float64}
    density::Vector{Float64}
end


text:
function make_increasing_vector()
	x0 = -1.0;
	dWageV = [0.1, 3.2];
    wageInter = Param(:x0, "Wage HSG", "wHSG",
        x0, x0,  x0 .* 0.05, x0 .* 10,  true);
    dWage = Param(:dxV, "Wage gradient", "dw",
        dWageV, dWageV, dWageV .* 0.05, dWageV .* 10, true);

	ownId = ObjectId(:test1);
    pvec = ParamVector(ownId,  [wageInter, dWage]);
	return IncreasingVector(ownId, pvec, x0, dWageV)
end

# function increasing_vector_test()
	@testset "IncreasingVector" begin
		iv = make_increasing_vector();
		xV = values(iv);
		@test all(diff(xV) .> 0.0)
		@test xV[1] ≈ iv.x0
		@test xV[end] ≈ iv.x0 + sum(iv.dxV)
	end
# end

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


text:
using Pkg;
Pkg.activate(@__DIR__)
Pkg.instantiate()

using MLJ
using PrettyPrinting
import DataFrames: DataFrame, select!, Not, describe
import Statistics
using Dates
using UrlDownload
using CSV

df = DataFrame(urldownload("https://raw.githubusercontent.com/tlienart/DataScienceTutorialsData.jl/master/data/kc_housing.csv", true))
describe(df)

df.is_renovated = df.yr_renovated .== 0

select!(df, Not([:id, :date, :yr_renovated]))
CSV.write(joinpath(@__DIR__, "..", "house.csv"), df)


text:
module MOITests

using BARON
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

const optimizer = BARON.Optimizer(PrLevel=0)
const caching_optimizer = MOIU.CachingOptimizer(MOIU.Model{Float64}(), BARON.Optimizer(PrLevel=0));

@testset "Unit" begin
    config = MOIT.TestConfig(atol=1e-5, rtol=1e-4, infeas_certificates=true, duals=false)
    # A number of test cases are excluded because loadfromstring! works only
    # if the solver supports variable and constraint names.
    exclude = ["delete_variable", # Deleting not supported.
               "delete_variables", # Deleting not supported.
               "solve_zero_one_with_bounds_1", # loadfromstring!
               "solve_zero_one_with_bounds_2", # loadfromstring!
               "solve_zero_one_with_bounds_3", # loadfromstring!
               "getconstraint", # Constraint names not suported.
               "solve_with_upperbound", # loadfromstring!
               "solve_with_lowerbound", # loadfromstring!
               "solve_integer_edge_cases", # loadfromstring!
               "solve_affine_lessthan", # loadfromstring!
               "solve_affine_greaterthan", # loadfromstring!
               "solve_affine_equalto", # loadfromstring!
               "solve_affine_interval", # loadfromstring!
               "solve_duplicate_terms_vector_affine", # Vector variables
               "solve_blank_obj", # loadfromstring!
               "solve_affine_deletion_edge_cases", # Deleting not supported.
               "number_threads", # NumberOfThreads not supported
               "delete_nonnegative_variables", # get ConstraintFunction n/a.
               "update_dimension_nonnegative_variables", # get ConstraintFunction n/a.
               "delete_soc_variables", # VectorOfVar. in SOC not supported
               "solve_result_index", # DualObjectiveValue not supported
               "time_limit_sec", # Supported by Optimizer, but not by MOIU.Model
               "silent",
               "raw_status_string",
               "solve_qp_edge_cases",
               "solve_objbound_edge_cases"
               ]
    MOIT.unittest(caching_optimizer, config, exclude)
end

MOI.empty!(optimizer)

@testset "MOI Continuous Linear" begin
    config = MOIT.TestConfig(atol=1e-5, rtol=1e-4, infeas_certificates=false, duals=false)
    excluded = String[
        "linear7", # vector constraints
        "linear8b", # certificate provided in this case (result count is 1)
        "linear8c", # should be unbounded below, returns "Preprocessing found feasible solution with value -.200000000000E+052"
        "linear15", # vector constraints
        "partial_start" # TODO
    ]
    # MOIT.partial_start_test(optimizer, config)
    MOIT.contlineartest(caching_optimizer, config, excluded)
    MOIT.linear8btest(caching_optimizer, MOIT.TestConfig(atol=1e-5, rtol=1e-4, infeas_certificates=true, duals=false))
end

MOI.empty!(optimizer)

@testset "MOI Integer Linear" begin
    config = MOIT.TestConfig(atol=1e-5, rtol=1e-4, infeas_certificates=false, duals=false)
    excluded = String[
        "int2", # SOS1
        "indicator1", # ACTIVATE_ON_ONE
        "indicator2", # ACTIVATE_ON_ONE
        "indicator3", # ACTIVATE_ON_ONE
        "indicator4", # ACTIVATE_ON_ONE
    ]
    MOIT.intlineartest(caching_optimizer, config, excluded)
end

MOI.empty!(optimizer)

@testset "MOI Continuous Quadratic" begin
    # TODO: rather high tolerances
    config = MOIT.TestConfig(atol=1e-3, rtol=1e-3, infeas_certificates=false, duals=false)
    excluded = String[
        "qcp1", # vector constraints
    ]
    MOIT.contquadratictest(caching_optimizer, config, excluded)
end

MOI.empty!(optimizer)
bridged = MOIB.full_bridge_optimizer(optimizer, Float64)

@testset "MOI Nonlinear" begin
    config = MOIT.TestConfig(atol=1e-5, rtol=1e-4, infeas_certificates=false, duals=false)
    excluded = String[
        "nlp_objective_and_moi_objective",
    ]
    MOIT.nlptest(bridged, config, excluded)
end

end # module


text:
module SimpleGraphics

__precompile__(true)

export 	toSVG, Show, Color, GraphicsElement, VertexList, GraphicsAttributes, Blank,
				Group, Rectangle, Circle, Polygon, translate, translate!,
	   		rotate, rotate!, transform, transform!,
				centreOfGravity, TextElement, maxExtents, Export

abstract type GraphicsElement end
const GraphicsAttributes = Dict{String,Any}
const VertexList = Array{Float32,2}

# =======================================================
# Blank: a NOP GraphicsElement
# =======================================================

struct Blank <: GraphicsElement
end

# =======================================================
# Groups
# =======================================================

struct Group <: GraphicsElement
	elements::Array{GraphicsElement,1}
	attr::GraphicsAttributes
	clipPath::GraphicsElement

	function Group(elements=Array{GraphicsElement,1}([]),attr=GraphicsAttributes([]),cp=Blank())
		new(elements,attr,cp)
	end
end

function push!(g::Group,e::GraphicsElement)
	 Base.push!(g.elements,e)
end

# =======================================================
# Rectangles
# =======================================================

#  <rect x="50" y="20" width="150" height="150"  style="fill:blue;stroke:pink;stroke-width:5;fill-opacity:0.1;stroke-opacity:0.9" />

struct Rectangle <: GraphicsElement
	vertices::VertexList
	attr::GraphicsAttributes

	function Rectangle(x::Real,y::Real,w::Real,h::Real,attr=GraphicsAttributes([]))
			r=new([x y ; x+w y+h],attr)
	end

end

# =======================================================
# Polygons and Polylines
# =======================================================

struct
	Polygon <: GraphicsElement
	vertices::VertexList
	attr::GraphicsAttributes
	Polygon(list::VertexList,attr=GraphicsAttributes([]))=new(list,attr)
end

# =======================================================
# Ellipses and Circles
# =======================================================

mutable struct Circle <: GraphicsElement
	vertices::VertexList
	attr::GraphicsAttributes
	Circle(list::VertexList) = new(list,GraphicsAttributes([]) )
	Circle(x::Number,y::Number,r::Number) = new(VertexList([x y ; x+r y]),GraphicsAttributes([]))
end

# =======================================================
# Text elements
# =======================================================

struct TextElement <: GraphicsElement
	vertices::VertexList
	text::String
	orientation::Float32
	attr::GraphicsAttributes
	TextElement(v::VertexList,s::String,attr=GraphicsAttributes([]),o=0.0) =
		new(v,s,o,attr)
end

# =======================================================
# SVG Output
# =======================================================

# The following function maps vertices onto the SVG coordinate system

function mapSVG(v::VertexList)
	 metric = Float32[1 0 ; 0 -1]
	 return v*metric
end

function toSVG(g::Group)

	s=""
	clipName="$(round(10^7*rand()))"

	if !isa(g.clipPath,Blank)
		s="$(s)<clipPath id=\"$(clipName)\">"
		s="$(s) $(toSVG(g.clipPath))"
		s="$(s)</clipPath>"
		s="$(s)<g clip-path=\"url(#$(clipName))\""
	else
		s="$(s)<g "
	end

	for k in keys(g.attr)
		s="$s $k=\"$(g.attr[k])\""
	end
	s="$s >"

	for p in g.elements
		s="$s $(toSVG(p))\n"
	end
	s="$s </g>\n"
end

function toSVG(p::Rectangle)
	vert=mapSVG(p.vertices)
	x=vert[1,1]
	y=vert[1,2]
	w=round(vert[2,1]-vert[1,1],3)
	h=-round(vert[2,2]-vert[1,2],3)
	s="<rect x=\"$x\" y=\"$(y-h)\" width=\"$w\" height=\"$h\" "
	for k in keys(p.attr)
		s="$s $k=\"$(p.attr[k])\""
	end
	s = "$s />"
end

# first point: centre
# second point: on the periphery
function toSVG(c::Circle)
	v=mapSVG(c.vertices)
	s="<circle cx=\"$(round(v[1,1],3))\" cy=\"$(round(v[1,2],3))\" "
	r=norm(v[2,:]-v[1,:])
	s="$s r=\"$(round(r,3))\""
	for k in keys(c.attr)
		s="$s $k=\"$(c.attr[k])\""
	end
	s="$s />"
	s
end

function toSVG(t::TextElement)
	v=mapSVG(t.vertices)
	s="<text transform=\"rotate($(t.orientation),$(round(v[1,1],3)),$(round(v[1,2],3))) translate($(round(v[1,1],3)),$(round(v[1,2],3)))\""
	for k in keys(t.attr)
		s="$s $k=\"$(t.attr[k])\""
	end
	s="$s>"
	s="$s $(t.text) </text>"
end

function toSVG(p::Polygon)
	if get(p.attr,"closed",true)
		# s="<polyline fill=\"none\" stroke=\"blue\" stroke-width=\"0.1\" "
		s="<polyline "
	else
		s="<polyline "
	end
	for k in keys(p.attr)
		s="$s $k=\"$(p.attr[k])\""
	end
	s="$s points=\""
	n,l=size(p.vertices)
	v=mapSVG(p.vertices)
	for k=1:n
	   s="$s $(round(v[k,1],3)),$(round(v[k,2],3)) "
	end
	s="$s \""
	s = "$s />"
end

function toSVG(g::Blank)
	return ""
end

# =======================================================
# Generic Transformations
# =======================================================

function translate!(p::GraphicsElement,dx::Real,dy::Real)
	p.vertices .+= [dx dy]
end

function translate!(g::Group,dx::Real,dy::Real)
	for e in g.elements
		translate!(e,dx,dy)
	end
end

function translate(p::GraphicsElement,dx::Real,dy::Real)
	pn=deepcopy(p)
	translate!(pn,dx,dy)
	pn
end

function rotate!(p::GraphicsElement,θ::Real,axis=[0. 0.])
	rmatrix=[cos(θ) sin(θ) ; -sin(θ) cos(θ)]
	p.vertices .-= axis
	p.vertices *= rmatrix
	p.vertices .+= axis
end

function rotate!(g::Group,θ::Real,axis=[0. 0.])
	for e in g.elements
		rotate!(e,θ,axis)
	end
end

function transform!(p::GraphicsElement,Q=Float32[1 0; 0 1])
	p.vertices *= Q
end

function transform!(g::Group,Q=Float32[1 0; 0 1])
	for e in g.elements
		transform!(e,Q)
	end
end

function transform(p::GraphicsElement,Q=Float32[1 0; 0 1])
	pn=deepcopy(p)
	transform!(p,Q)
	return p
end

function rotate(p::GraphicsElement,θ::Real,axis=[0. 0.])
	p1=deepcopy(p)
	rotate!(p1,θ,axis)
	p1
end

function centreOfGravity(p::GraphicsElement)
	cg=sum(p.vertices,1)/(size(p.vertices)[1])
end

# ================================================
# Output
# ================================================

"""
`maxExtents(g,def)` returns a 2x2 array with the bottom left and top right points of the
bounding box fitting around `g` _and_ `def`, where `def` is a 2x2 array
"""
function maxExtents(v::VertexList,def=VertexList([0 0 ; 0 0]))
	return VertexList([minimum([v;def] ,1) ; maximum([v;def],1)])
end

function maxExtents(g::GraphicsElement,def=VertexList([0 0 ; 0 0]))
	v= (:vertices in fieldnames(g)) ? g.vertices : VertexList([])
	return maxExtents(v,def)
end

function maxExtents(g::Group,def=VertexList([0 0 ; 0 0]))
	 v=def
	 for e in g.elements
		 v=[v; maxExtents(e,def)]
	 end
	 return maxExtents(v,def)
 end

function Show(g::GraphicsElement)
	bbox=maxExtents(g)
  display("image/svg+xml",svgWrap(toSVG(g),ext=1.1*mapSVG(bbox)))
end

function Export(f::IOStream,g::GraphicsElement)
	bbox=maxExtents(g)
	write(f,svgWrap(toSVG(g),ext=1.1*mapSVG(bbox)))
end

function Export(s::String,g::GraphicsElement)
	f=open(s,"w")
	Export(f,g)
	close(f)
end

function svgWrap(s::String;width="15cm",height="15cm",ext=Float32[-500 -500; 500 500])

  return "<?xml version=\"1.0\" standalone=\"no\"?>" *
         "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\" \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">" *
         "<svg width=\"$width\" height=\"$height\" viewBox=\"$(ext[1,1]) $(ext[2,2]) $(ext[2,1]-ext[1,1]) $(-ext[2,2]+ext[1,2]) \" xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\">"*
         s *
         "</svg>"
end

Base.push!(g::SimpleGraphics.Group,e::SimpleGraphics.GraphicsElement)=SimpleGraphics.push!(g,e)

end # module SimpleGraphics

max_stars_repo_name: marcel-utz/NMR

max_issues_repo_name: marcel-utz/NMR


text:
module BenchGetSetIndex

using BenchmarkTools
using ITensors

suite = BenchmarkGroup()

i = Index(100)
A = randomITensor(i', i)

suite["setindex!"] = @benchmarkable $A[1, 1] = 1.0
suite["setindex! end"] = @benchmarkable $A[1, end] = 1.0
suite["setindex! indval"] = @benchmarkable $A[$(i') => 1, $i => 1] = 1.0
suite["setindex! indval end"] = @benchmarkable $A[$(i') => 1, $i => end] = 1.0
suite["getindex"] = @benchmarkable $A[1, 1]
suite["getindex end"] = @benchmarkable $A[1, end]
suite["getindex indval"] = @benchmarkable $A[$(i') => 1, $i => 1]
suite["getindex indval end"] = @benchmarkable $A[$(i') => 1, $i => end]
end

BenchGetSetIndex.suite


text:
n = rand(1:10)
@test matrixdepot("parter", n) == matrixdepot("parter", Float64, n)

A = matrixdepot("cauchy", [1.:n;] + 0.5, -[1.:n;])

@test_approx_eq matrixdepot("parter", n) A
println("'parter' passed test...")


text:
# Plot the total constraint violation of the course of the solver

using TrajOptSOCPs, Plots

function plotConstraintViolation(cM::constraintManager, trajList,
                                 penalty::Float64 = 1.0)
    cVList = getConstraintViolationList(cM, trajList, penalty)

    pltCV = plot(safeNorm(cVList), markershape = :square, yaxis = :log)
    title!("Inequality Constraint Violation")
    ylabel!("Constraint Violation at Penalty of $penalty")
    xlabel!("Step")

    pltCV2 = plot()

    if typeof(cM) == constraintManager_Dynamics
        cVListDyn = [safeNorm(norm(evalAffineEq(cM, traj))) for traj in trajList]
        plot!(cVListDyn, markershape = :square, markercolor = :darkred,
                    linecolor = :darkred, yaxis = :log)
        title!("Dynamics Constraint Violation")
        ylabel!("Constraint Violation at Penalty of $penalty")
        xlabel!("Step")
    end

    return pltCV, pltCV2
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
abstract type AbstractBorderedLinearSolver <: AbstractLinearSolver end

# call for using BorderedArray input, specific to Arclength Continuation
(lbs::AbstractBorderedLinearSolver)(J, dR::vec1, dz::vec2, R, n::T, theta::T; shift::Ts = nothing) where {vec1, vec2, T, Ts} = (lbs)(J, dR, dz.u, dz.p, R, n, theta / length(dz.u), one(T) - theta; shift = shift)

####################################################################################################
@with_kw struct BorderingBLS{S <: AbstractLinearSolver, Ttol} <: AbstractBorderedLinearSolver
	solver::S = DefaultLS()
	tol::Ttol = 1e-12
	checkPrecision::Bool = false
end

# dummy constructor to simplify user passing options to continuation
BorderingBLS(ls::AbstractLinearSolver) = BorderingBLS(solver = ls)

# solve in dX, dl
#          J  * dX +       dR   * dl = R
# xiu * dz.u' * dX + xip * dz.p * dl = n
function (lbs::BorderingBLS{S, Ttol})(  J, dR,
								dzu, dzp::T, R, n::T,
								xiu::T = T(1), xip::T = T(1); shift::Ts = nothing)  where {T, S, Ts, Ttol}
	# the following parameters are used for the pseudo arc length continuation
	# xiu = theta / length(dz.u)
	# xip = one(T) - theta

	# we make this branching to avoid applying a zero shift
	if isnothing(shift)
		x1, x2, _, (it1, it2) = lbs.solver(J, R, dR)
	else
		x1, x2, _, (it1, it2) = lbs.solver(J, R, dR; a₀ = shift)
	end

	dl = (n - dot(dzu, x1) * xiu) / (dzp * xip - dot(dzu, x2) * xiu)

	# dX = x1 .- dl .* x2
	axpy!(-dl, x2, x1)

	# we check the precision of the solution by the bordering algorithm
	# mainly for debugging purposes
	if lbs.checkPrecision
		# at this point, x2 is not used anymore, we can use it for computing the residual
		# hence x2 = J*x1 + dl*dR - R
		x2 = apply(J, x1)
		axpy!(dl, dR, x2)
		axpy!(-1, R, x2)

		printstyled(color=:red,"--> res = ", ( norm(x2), abs(n - xip*dzp*dl -xiu* dot(dzu, x1))), "\n")
	end

	return x1, dl, true, (it1, it2)
end
####################################################################################################
# this interface should work for Sparse Matrices as well as for Matrices
struct MatrixBLS{S <: AbstractLinearSolver} <: AbstractBorderedLinearSolver
	solver::S
end

# dummy constructor to simplify user passing options to continuation
MatrixBLS() = MatrixBLS(DefaultLS())

# case of a scalar additional linear equation
function (lbs::MatrixBLS)(J, dR,
						dzu, dzp::T, R::vectype, n::T,
						xiu::T = T(1), xip::T = T(1); shift::Ts = nothing)  where {T <: Number, vectype <: AbstractVector, S, Ts}

	if isnothing(shift)
		A = J
	else
		A = J + shift * I
	end

	A = hcat(A, dR)
	A = vcat(A, vcat(vec(dzu) .* xiu, dzp * xip)')

	# solve the equations and return the result
	rhs = vcat(R, n)
	res = A \ rhs
	return res[1:end-1], res[end], true, 1
end

# # case of a scalar additional linear equation
# function (lbs::MatrixBLS)(J, dR,
# 						dzu, dzp, R::vectype, n::Tv,
# 						xiu::T = T(1), xip::T = T(1); shift::Ts = nothing)  where {Tv, T <: Number, vectype <: AbstractVector, S, Ts}
# 	@warn "Experimental!! Solving bordered linear system with non scalar bordered equation"
# 	N = size(J, 1)
# 	n1, n2 = size(dzp)
#
# 	rhs = vcat(R, n)
# 	A = similar(J, N + n1, N + n2)
#
# 	if shift == nothing
# 		A[1:N, 1:N] .= J
# 	else
# 		A[1:N, 1:N] .= J + shift * I
# 	end
#
# 	A[1:N, end-n1+1:end] .= dR
# 	A[end-n1+1:end, 1:N] .= (dzu) .* xiu
# 	A[end-n1+1:end, end-n1+1:end]  .= dzp .* xip
# 	res = A \ rhs
# 	return res[1:end-n1], res[end-n1+1:end], true, 1
# end
####################################################################################################
# composite type to save the bordered linear system with expression
# [ J	a]
# [b'	c]
# It then solved using Matrix Free algorithm applied to the full operator and not just J as for MatrixFreeBLS
#
struct MatrixFreeBLSmap{Tj, Ta, Tb, Tc}
	J::Tj
	a::Ta
	b::Tb
	c::Tc
end

function (lbmap::MatrixFreeBLSmap{Tj, Ta, Tb, Tc})(x::BorderedArray{Ta, Tc}) where {Tj, Ta, Tb, Tc <: Number}
	# This implements the case where Tc is a number, ie there is one scalar constraint in the
	# bordered linear system
	out = similar(x)
	copyto!(out.u, apply(lbmap.J, x.u))
	axpy!(x.p, lbmap.a, out.u)
	out.p = dot(lbmap.b, x.u)  + lbmap.c  * x.p
	return out
end

function (lbmap::MatrixFreeBLSmap{Tj, Ta, Tb, Tc})(x::AbstractArray) where {Tj, Ta <: AbstractArray, Tb, Tc <: Number}
	# This implements the case where Tc is a number, ie there is one scalar constraint in the
	# bordered linear system
	out = similar(x)
	xu = @view x[1:end-1]
	xp = x[end]
	# copyto!(out.u, apply(lbmap.J, x.u))
	out[1:end-1] .= @views apply(lbmap.J, xu) .+ xp .* lbmap.a
	out[end] = @views dot(lbmap.b, xu)  + lbmap.c  * xp
	return out
end

struct MatrixFreeBLS{S} <: AbstractBorderedLinearSolver
	solver::S
end

# dummy constructor to simplify user passing options to continuation
MatrixFreeBLS() = MatrixFreeBLS(DefaultLS())

extractVector(x::AbstractVector) = @view x[1:end-1]
extractVector(x::BorderedArray) = x.u

extractParameter(x::AbstractVector) = x[end]
extractParameter(x::BorderedArray) = x.p

# We restrict to bordered systems where the added component is scalar
function (lbs::MatrixFreeBLS{S})(J, 		dR,
								dzu, 	dzp::T, R, n::T,
								xiu::T = T(1), xip::T = T(1); shift::Ts = nothing) where {T <: Number, S, Ts}
	linearmap = MatrixFreeBLSmap(J, dR, rmul!(copy(dzu), xiu), dzp * xip)
	if S <: GMRESIterativeSolvers
		rhs = vcat(R, n)
	else
		rhs = BorderedArray(copy(R), n)
	end
	sol, cv, it = lbs.solver(linearmap, rhs)
	return extractVector(sol), extractParameter(sol), cv, it
end
####################################################################################################
# Nested algorithm for solving the bordered linear system


text:
@testset "TestTTest" begin

    # OneSampleTTest tests
    # Loop over n, d, p, s and one_sided
    # TODO: Test one-sided vs two-sided
    params = (
        (33.36712895331925, 0.5, 0.80, 0.05, false),
        (43.995480906737704, 0.5, 0.90, 0.05, false),
        (53.940614296961286, 0.5, 0.95, 0.05, false),
        (16, 0.7494656445134399, 0.8, 0.05, false),
        (32, 0.5112738234562735, 0.8, 0.05, false),
        (64, 0.35565012796825096, 0.8, 0.05, false),
        (16, 0.5, 0.46487, 0.05, false),
        (32, 0.5, 0.782276, 0.05, false),
        (64, 0.5, 0.9760556, 0.05, false),
    )

    for (n, d, p, s, one_sided) in params
        p_est = power(OneSampleTTest, n, d, 1.0, s, one_sided)
        n_est = sample_size(OneSampleTTest, d, 1.0, p, s, one_sided)
        d_est = effect_size(OneSampleTTest, n, p, s, one_sided)

        @test abs(p_est - p) < 1e-4
        @test abs(n_est - n) < 1e-4
        @test abs(d_est - d) < 1e-4
    end

    # EqualVarianceTTest tests
    # Loop over n, d, p, s and one_sided
    # TODO: Test one-sided vs two-sided
    params = (
        (63.76561019093564, 0.5, 0.80, 0.05, false),
        (85.03128413735288, 0.5, 0.90, 0.05, false),
        (104.92794298436287, 0.5, 0.95, 0.05, false),
        (16, 1.0236637440012688, 0.8, 0.05, false),
        (32, 0.7114799360863416, 0.8, 0.05, false),
        (64, 0.4990691779657755, 0.8, 0.05, false),
        (16, 0.5, 0.2777445, 0.05, false),
        (32, 0.5, 0.5036382, 0.05, false),
        (64, 0.5, 0.8014596, 0.05, false),
    )

    for (n, d, p, s, one_sided) in params
        p_est = power(EqualVarianceTTest, n, d, 1.0, s, one_sided)
        n_est = sample_size(EqualVarianceTTest, d, 1.0, p, s, one_sided)
        d_est = effect_size(EqualVarianceTTest, n, p, s, one_sided)

        @test abs(p_est - p) < 1e-4
        @test abs(n_est - n) < 1e-4
        @test abs(d_est - d) < 1e-4
    end
end


text:
function dataframe(response::Requests.Response)

    buffer = PipeBuffer() # open a buffer in which to dump data
    df     = DataFrame()  # init empty DataFrame

    try
        # Write the data to the buffer
        write(buffer, response.data)

        # Use DataFrame's readtable to read the data directly from buffer
        df = readtable(buffer)

        # Convert dates to DateTime object
        df[:Date] = Date{ISOCalendar}[date(d) for d in df[:Date]]
    finally
        close(buffer)
    end

    return df
end


text:
using RegularizedLeastSquares, LinearAlgebra
# Packages for testing only
using Random, Test
using SparsityOperators: fft

include("Kaczmarz.jl")
include("testProxMaps.jl")
include("testSolvers.jl")


text:
using GridUtilities
using Test
##using TestSetExtensions

#@test isempty(detect_ambiguities(ViscousFlow))
include("histories.jl")

#@testset ExtendedTestSet "All tests" begin
#    @includetests ARGS
#end

#if isempty(ARGS)
#    include("../docs/make.jl")
#end


text:
#####################################################################
##### GENERIC FUNCTIONS                                         #####
#####################################################################

isbinary(y) = sort(unique(y)) == [0, 1]

#####################################################################
#####                     QUERY PARSING                          ####
#####################################################################

function parse_queries(queryfile::String)
    config = TOML.parsefile(queryfile)
    queries = Query[]
    for (queryname, querydict) in config
        if lowercase(queryname) ∉ ("threshold", "snps")
            rsids = keys(querydict)
            names = Tuple(Symbol(x) for x in rsids)
            rsid_vals = [split(filter(x->!isspace(x), querydict[rsid]), "->") for rsid in rsids]
            case = NamedTuple{names}(rsid_val[2] for rsid_val in rsid_vals)
            control = NamedTuple{names}(rsid_val[1] for rsid_val in rsid_vals)
            push!(queries,  Query(case=case, control=control, name=queryname))
        end
    end
    return queries
end

#####################################################################
#####                 PHENOTYPE PARSING                          ####
#####################################################################

phenotypesnames(phenotype_listfile::Nothing) = nothing

function phenotypesnames(phenotype_listfile::String)
    phenotypeslist = open(readlines, phenotype_listfile)
    return vcat(["SAMPLE_ID"], phenotypeslist)
end

function load_phenotypes(phenotypes_datafile, phenotypes_listfile)
    columns = phenotypesnames(phenotypes_listfile)
    return CSV.File(phenotypes_datafile, select=columns) |> DataFrame
end

#####################################################################
#####                 CV ADAPTIVE FOLDS                          ####
#####################################################################

countuniques(v::AbstractVector) = [count(==(u), v) for u in unique(v)]
countuniques(table) =
    countuniques([values(x) for x in Tables.namedtupleiterator(table)])

"""
    AdaptiveCV(cv::Union{CV, StratifiedCV})

Implements the rule of thum given here: https://www.youtube.com/watch?v=WYnjja8DKPg&t=4s
"""
mutable struct AdaptiveCV <: MLJBase.ResamplingStrategy
    cv::Union{CV, StratifiedCV}
end

function MLJBase.train_test_pairs(cv::AdaptiveCV, rows, y)
    # Compute n-eff
    n = nrows(y)
    neff =
        if scitype(first(y)) == MLJBase.Continuous
            n
        else
            counts = countuniques(y)
            nrare = minimum(counts)
            min(n, 5*nrare)
        end

    # Compute number of folds
    nfolds =
        if neff < 30
            neff
        elseif neff < 500
            20
        elseif neff < 5000
            10
        elseif neff < 10_000
            5
        else
            3
        end

    # Update the underlying n_folds
    adapted_cv = typeof(cv.cv)(nfolds=nfolds, shuffle=cv.cv.shuffle, rng=cv.cv.rng)

    return MLJBase.train_test_pairs(adapted_cv, rows, y)
end


text:
# SCRIPT

function main()
    for a in 1:998
        for b in (a + 1):(1000 - a)
            c = 1000 - (a + b);
            if (a * a + b * b == c * c)
                println(a, " ", b, " ", c)
                exit()
            end
        end
    end
end

main()

max_stars_repo_name: shanedrabing/project-euler

max_issues_repo_name: shanedrabing/project-euler

max_forks_repo_name: shanedrabing/project-euler


text:
using Base.Test
import Utils

@test Utils.sumt(Int, i for i = 1:10 if i < 0) === 0
@test Utils.prodt(Int, i for i = 1:10 if i < 0) === 1

@test Utils.simplexvol([[0., 0.], [1., 0.], [0., 1.]]) == 0.5

@test Utils.lbinomial(10, 5) ≈ 5.5294290875114233067303074602045842750830661139479


text:
"""
Buffer backed by memory of type `M`.

Offsets, strides and sizes are always expressed in bytes.
"""
abstract type Buffer{M<:Memory} <: LavaAbstraction end

memory_type(::Type{<:Buffer{M}}) where {M} = M
memory_type(t) = memory_type(typeof(t))

vk_handle_type(::Type{<:Buffer}) = Vk.Buffer

Vk.bind_buffer_memory(buffer::Buffer, memory::Memory) = Vk.bind_buffer_memory(device(buffer), buffer, memory, offset(memory))

# abstract type DataOperation end

# abstract type Transfer <: DataOperation end
# abstract type VideoDecode <: DataOperation end
# abstract type VideoEncode <: DataOperation end

# abstract type StorageType end

# abstract type Vertex <: StorageType end
# abstract type Index{T} <: StorageType end

# abstract type SparsityType end

# abstract type SparseResidency <: SparsityType end
# abstract type SparseBinding <: SparsityType end

# abstract type Sparse{T<:SparsityType,O<:LavaAbstraction} <: StorageType end

abstract type DenseBuffer{M<:DenseMemory} <: Buffer{M} end

isallocated(buffer::DenseBuffer) = isdefined(buffer.memory, 1)
memory(buffer::DenseBuffer) = buffer.memory[]

struct BufferBlock{M<:DenseMemory} <: DenseBuffer{M}
    handle::Vk.Buffer
    size::Int
    usage::Vk.BufferUsageFlag
    queue_family_indices::Vector{Int8}
    sharing_mode::Vk.SharingMode
    memory::Ref{M}
end

Base.size(buffer::BufferBlock) = buffer.size

device_address(buffer::BufferBlock) = Vk.get_buffer_device_address(device(buffer), Vk.BufferDeviceAddressInfo(handle(buffer)))

function BufferBlock(device, size; usage = Vk.BufferUsageFlag(0), queue_family_indices = queue_family_indices(device), sharing_mode = Vk.SHARING_MODE_EXCLUSIVE, memory_type = MemoryBlock, allocate = false)
    usage |= Vk.BUFFER_USAGE_SHADER_DEVICE_ADDRESS_BIT
    info = Vk.BufferCreateInfo(size, usage, sharing_mode, queue_family_indices)
    handle = unwrap(create(BufferBlock, device, info))
    buffer = BufferBlock(handle, size, usage, convert(Vector{Int8}, queue_family_indices), sharing_mode, Ref{memory_type}())
end

function Base.show(io::IO, block::BufferBlock)
    print(io, BufferBlock, "($(Base.format_bytes(size(block))), $(block.usage)")
    if !isallocated(block)
        print(io, ", unallocated")
    end
    print(io, ')')
end

function buffer(device, data, device_local::Val{false}; kwargs...)
    buffer = BufferBlock(device, sizeof(data); kwargs...)
    allocate!(buffer, MEMORY_DOMAIN_HOST)
    copyto!(buffer, data; device)
    buffer
end

"""
Create a new buffer and fill it with data.

Note that as the buffer is device-local, a transfer via staging buffer is required.
In this case, a tuple of `(buffer, state)` is returned.
`state` is an [`ExecutionState`](@ref) whose bound resources **must** be
preserved until the transfer has been completed. See the related documentation for more details.

!!! warning
    Failure to comply will result in undefined behavior, potentially leading to crashes or invalid data.
"""
function buffer(device, data, device_local::Val{true}; usage = Vk.BufferUsageFlag(0), kwargs...)
    usage |= Vk.BUFFER_USAGE_TRANSFER_DST_BIT
    buffer = BufferBlock(device, sizeof(data); usage, kwargs...)
    allocate!(buffer, MEMORY_DOMAIN_DEVICE)
    state = copyto!(buffer, data; device)
    buffer, state
end

"""
Create a new buffer and fill it with data.

By default, the buffer will be host-visible.

You can use `buffer(device, data, Val(true))` to make it device-local.
Note that you will have to preserve a temporary staging
buffer while the command is executing, as indicated in the documentation.
"""
buffer(device, data; kwargs...) = buffer(device, data, Val(false); kwargs...)

function Base.similar(buffer::BufferBlock{T}; memory_domain = nothing, usage = buffer.usage) where {T}
    similar = BufferBlock(device(buffer), size(buffer); usage, buffer.queue_family_indices, buffer.sharing_mode, memory_type = T)
    if isallocated(buffer)
        memory_domain = @something(memory_domain, memory(buffer).domain)
        allocate!(similar, memory_domain)
    end
    similar
end

function Base.copyto!(buffer::BufferBlock, data; device = nothing)
    mem = memory(buffer)
    if Vk.MEMORY_PROPERTY_DEVICE_LOCAL_BIT in properties(mem)
        device::Device
        tmp = similar(buffer; memory_domain = MEMORY_DOMAIN_HOST, usage = Vk.BUFFER_USAGE_TRANSFER_SRC_BIT)
        copyto!(tmp, data)
        transfer(device, tmp, buffer; semaphore = Vk.SemaphoreSubmitInfoKHR(Vk.Semaphore(device), 0, 0), free_src = true, signal_fence = true)
    elseif Vk.MEMORY_PROPERTY_HOST_VISIBLE_BIT in properties(mem)
        map(mem) do ptr
            ptrcopy!(ptr, data)
        end
    else
        error("Buffer not visible neither to device nor to host (memory properties: $(properties(mem))).")
    end
end

function ptrcopy!(ptr, data::DenseArray{T}) where {T}
    GC.@preserve data unsafe_copyto!(Ptr{T}(ptr), pointer(data), length(data))
end

function ptrcopy!(ptr, data::T) where {T}
    ref = Ref{T}(data)
    GC.@preserve ref unsafe_copyto!(Ptr{T}(ptr), Base.unsafe_convert(Ptr{T}, ref), 1)
end

function ptrcopy!(ptr, data::String)
    GC.@preserve data unsafe_copyto!(Ptr{UInt8}(ptr), Base.unsafe_convert(Ptr{UInt8}, data), sizeof(data))
end

function transfer(device, src::Buffer, dst::Buffer; command_buffer = request_command_buffer(device, Vk.QUEUE_TRANSFER_BIT), signal_fence = false, semaphore = nothing, free_src = false)
    @assert size(src) == size(dst)

    Vk.cmd_copy_buffer(command_buffer, src, dst, [Vk.BufferCopy(offset(src), offset(dst), size(src))])

    signal_semaphores = []
    !isnothing(semaphore) && push!(signal_semaphores, semaphore)
    info = Vk.SubmitInfo2KHR([], [Vk.CommandBufferSubmitInfoKHR(command_buffer)], signal_semaphores)
    if free_src
        submit(device, command_buffer.queue_family_index, info; signal_fence, semaphore, free_after_completion = [Ref(src)])
    else
        submit(device, command_buffer.queue_family_index, info; signal_fence, semaphore, release_after_completion = [Ref(src)])
    end
end

Base.collect(buffer::Buffer, device = nothing) = collect(memory(buffer), size(buffer), device)
Base.collect(::Type{T}, buffer::Buffer, device = nothing) where {T} = reinterpret(T, collect(buffer, device))

struct SubBuffer{B<:DenseBuffer} <: Buffer{SubMemory}
    buffer::B
    offset::Int
    stride::Int
    size::Int
end

@forward SubBuffer.buffer handle

Base.size(sub::SubBuffer) = sub.size

offset(buffer::Buffer) = 0
offset(buffer::SubBuffer) = buffer.offset

stride(buffer::Buffer) = 0
stride(buffer::SubBuffer) = buffer.stride

memory(sub::SubBuffer) = @view memory(sub.buffer)[offset(sub):(size(sub) - offset(sub))]

device_address(sub::SubBuffer) = device_address(sub.buffer) + UInt64(sub.offset)

function Base.view(buffer::DenseBuffer, range::StepRange)
    range_check(buffer, range)
    SubBuffer(buffer, range.start, range.step, range.stop - range.start)
end

function Base.view(buffer::DenseBuffer, range::UnitRange)
    range_check(buffer, range)
    SubBuffer(buffer, range.start, 0, range.stop - range.start)
end

Base.firstindex(buffer::Buffer) = offset(buffer)
Base.lastindex(buffer::Buffer) = size(buffer)

"""
Allocate a `MemoryBlock` and bind it to the provided buffer.
"""
function allocate!(buffer::DenseBuffer, domain::MemoryDomain)
    _device = device(buffer)
    reqs = Vk.get_buffer_memory_requirements(_device, buffer)
    memory = MemoryBlock(_device, reqs.size, reqs.memory_type_bits, domain)
    bind!(buffer, memory)
end

function bind!(buffer::BufferBlock, memory::Memory)
    unwrap(Vk.bind_buffer_memory(buffer, memory))
    buffer.memory[] = memory
    memory.is_bound[] = true
    buffer
end


text:
using Test, Pkg
using MPI

if haskey(Pkg.installed(), "CuArrays")
    using CuArrays
    ArrayType = CuArray
else
    ArrayType = Array
end

MPI.Init()

comm = MPI.COMM_WORLD
size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)

for T in Base.uniontypes(MPI.MPIDatatype)

    # Allocating version
    a = ArrayType(fill(T(rank), size))
    b = MPI.Alltoall(a, 1, comm)
    @test b isa ArrayType{T}
    @test b == ArrayType{T}(0:size-1)

    # Non Allocating version
    a = ArrayType(fill(T(rank),size))
    b = ArrayType{T}(undef, size*1)
    MPI.Alltoall!(a, b, 1, comm)
    @test b == ArrayType{T}(0:size-1)

    # IN_PLACE version
    a = ArrayType{T}(fill(T(rank),size))
    MPI.Alltoall!(MPI.IN_PLACE, a, 1, comm)
    @test a == ArrayType{T}(0:size-1)
end

MPI.Finalize()
@test MPI.Finalized()


text:
using DifferentialEquations, StatsBase, StatsPlots
using DifferentialEvolutionMCMC, Distributions, Random, StatsBase, LabelledArrays

Random.seed!(42)

function lotka_volterra(du, u, p, t)
    x, y = u
    α, β, γ, δ  = p
    du[1] = (α - β * y) * x
    du[2] = (δ * x - γ) * y
end

# true parameters
p = [1.5, 1.0, 3.0, 1.0]
# initial state values
u0 = [1.0,1.0]
problem  = ODEProblem(lotka_volterra, u0, (0.0,10.0), p)
sol = solve(problem, Tsit5())
plot(sol)

sol1 = solve(problem, Tsit5(),saveat=0.1)
data = Array(sol1) + 0.5 * randn(size(Array(sol1)))
plot(sol1, alpha = 0.3, legend = false); scatter!(sol1.t, data')

function prior_loglike(α, β, γ, δ, σ)
    LL = 0.0
    LL += logpdf(truncated(Normal(1.5, 0.5), 0.5 ,2.5), α)
    LL += logpdf(truncated(Normal(1.2, 0.5), 0, 2), β)
    LL += logpdf(truncated(Normal(3.0, 0.5), 1, 4), γ)
    LL += logpdf(truncated(Normal(1.0, 0.5), 0, 2), δ)
    LL += logpdf(InverseGamma(2, 3), σ)
    return LL
end

function sample_prior()
    α = rand(truncated(Normal(1.5, 0.5), 0.5 ,2.5))
    β = rand(truncated(Normal(1.2, 0.5), 0, 2))
    γ = rand(truncated(Normal(3.0, 0.5), 1, 4))
    δ = rand(truncated(Normal(1.0, 0.5), 0, 2))
    σ = rand(InverseGamma(2, 3))
    return [α, β, γ, δ, σ]
end

bounds = (
    (.5,2.5),
    (0, 2),
    (1,4),
    (0,2),
    (0,Inf)
)

names = (:α,:β,:γ,:δ,:σ)

# Log Likelihood function
function loglike(data, problem, θ...)
    prob = remake(problem, p=θ)
    σ = θ[end]
    predicted = solve(prob, Tsit5(), saveat=0.1)
    LL = 0.0
    for i in 1:length(predicted)
        LL += logpdf(MvNormal(predicted[i], σ), data[:,i])
    end
    return LL
end

model = DEModel(
    problem;
    loglike,
    prior_loglike,
    sample_prior,
    data,
    names
)

de = DE(;sample_prior, bounds, burnin=1000, Np=12, n_groups=3)
n_iter = 3000
chains = sample(model, de, MCMCThreads(), n_iter, progress=true)

chain_array = Array(chains)
pl = plot(grid=false)
for k in 1:300
    p = chain_array[rand(1:2000),:]
    resol = solve(remake(problem, p=p[1:4]),Tsit5(), saveat=0.1)
    sim_data = Array(resol) + p[end] * randn(size(Array(resol)))
    plot!(resol.t, sim_data', alpha=0.1, color = "#BBBBBB", legend = false)
end
scatter!(sol1.t, data');
plot!(sol1, w=1, legend = false)


text:
const PROJECT_DIR = joinpath(@__DIR__, "..")

using Pkg
Pkg.activate(PROJECT_DIR)

using YAML
using MolecularGraph
using MolecularGraph.Graph
using MolecularGraph.Util

const INPUT_DIR = joinpath(PROJECT_DIR, "./assets/raw/functionalgroup")
const INPUT_SA_DIR = joinpath(PROJECT_DIR, "./assets/raw/structuralalerts")
const OUTPUT_FILE = joinpath(PROJECT_DIR, "./assets/const/default_query_relations.yaml")

function run()
    # read file directories
    paths = [joinpath(INPUT_DIR, f) for f in readdir(INPUT_DIR)]
    # append!(paths, [joinpath(INPUT_SA_DIR, f) for f in readdir(INPUT_SA_DIR)])

    # read list of functional groups from .yaml files
    fgrecords = []
    for p in paths
        basename(p) == ".DS_Store" && continue
        for rcd in YAML.load(open(p))
            m = smartstomol(rcd["query"])
            if basename(p) == "PAINS.yaml"
                m = removehydrogens(m)
                inferatomaromaticity!(m)
            end
            convertnotquery!(m)
            rcd["qmol"] = m
            rcd["source"] = split(basename(p), ".")[1]
            push!(fgrecords, rcd)
        end
    end
    # Detect query relationships
    isaedges = Set{Tuple{Int,Int}}()
    hasedges = Set{Tuple{Int,Int}}()
    for u in 1:length(fgrecords)
        @debug "---" fgrecords[u]["key"]
        u % 50 == 0 && println(u, " records processed...")
        for v in 1:length(fgrecords)
            u == v && continue
            # isa
            if hasexactmatch(fgrecords[u]["qmol"], fgrecords[v]["qmol"])
                if (v, u) in isaedges
                    pop!(isaedges, (v, u))  # duplicate
                    if !haskey(fgrecords[u], "aliases")
                        fgrecords[u]["aliases"] = []
                    end
                    if !haskey(fgrecords[v], "aliases")
                        fgrecords[v]["aliases"] = []
                    end
                    push!(fgrecords[u]["aliases"], fgrecords[v]["key"])
                    push!(fgrecords[v]["aliases"], fgrecords[u]["key"])
                    @debug "dupes" fgrecords[u]["key"] fgrecords[v]["key"]
                else
                    push!(isaedges, (u, v))
                    @debug "isa" fgrecords[u]["key"] fgrecords[v]["key"]
                end
            end
            # has
            if hassubstructmatch(
                    fgrecords[u]["qmol"], fgrecords[v]["qmol"], atommatcher=recursiveatommatch)
                if (v, u) in hasedges
                    pop!(hasedges, (v, u))  # duplicate
                else
                    push!(hasedges, (u, v))
                    @debug "has" fgrecords[u]["key"] fgrecords[v]["key"]
                end
            end
        end
    end
    # Generate relationship graph
    isagraph = plaindigraph(length(fgrecords), isaedges)
    redisaedges = transitive_reduction(isagraph)
    hasgraph = plaindigraph(length(fgrecords), hasedges)
    redhasedges = transitive_reduction(hasgraph)
    # Add relationship to records
    for e in redisaedges
        (u, v) = getedge(isagraph, e)
        if !haskey(fgrecords[u], "isa")
            fgrecords[u]["isa"] = []
        end
        push!(fgrecords[u]["isa"], fgrecords[v]["key"])
    end
    for e in redhasedges
        (u, v) = getedge(hasgraph, e)
        (u, v) in isaedges && continue
        if !haskey(fgrecords[u], "has")
            fgrecords[u]["has"] = []
        end
        push!(fgrecords[u]["has"], fgrecords[v]["key"])
    end
    for rcd in fgrecords
        delete!(rcd, "qmol")
    end
    # Write records to the file
    YAML.write_file(OUTPUT_FILE, fgrecords)
end

run()


text:
"""
This does 0-1 (each item can be chosen only once) knapsack :
pack capacity = capacity
weight of each item = weights
value of each item = values
dp array is what the function works on
It returns the ans (dp[capacity])
```
julia> dp=zeros(Int,30)
julia> ZeroOnePack!(20,[1,3,11],[2,5,30],dp)
37
```
"""
function ZeroOnePack!(capacity::N,weights::V,values::V,dp::V)where N<:Number where V<:AbstractVector
    for i in 1:length(weights)
        j=capacity
        while j>weights[i] # reversed loop
            dp[j]=max(dp[j],dp[j-weights[i]]+values[i])
            j-=1
        end
        dp[weights[i]]=max(dp[weights[i]],values[i]) # dp[j]=max(dp[j],dp[0]+values[i])
    end
    return dp[capacity]
end

"""
This does complete/infinite (each item can be chosen infinite times) knapsack :
pack capacity = capacity
weight of each item = weights
value of each item = values
dp array is what the function works on
It returns the ans (dp[capacity])
```
julia> dp=zeros(Int,30)
julia> CompletePack!(20,[1,2,9],[1,3,20],dp)
43
```
"""
function CompletePack!(capacity::N,weights::V,values::V,dp::V)where N<:Number where V<:AbstractVector
    for i in 1:length(weights)
        dp[weights[i]]=max(dp[weights[i]],values[i]) # dp[j]=max(dp[j],dp[0]+values[i])
        for j in weights[i]+1:capacity
            dp[j]=max(dp[j],dp[j-weights[i]]+values[i])
        end
    end
    return dp[capacity]
end

max_stars_repo_name: arubhardwaj/Julia

max_issues_repo_name: arubhardwaj/Julia

max_forks_repo_name: arubhardwaj/Julia


text:
### A Pluto.jl notebook ###
# v0.18.1

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 49cb409b-e564-47aa-9dae-9bc5bffa991d
using UUIDs

# ╔═╡ 729ae3bb-79c2-4fcd-8645-7e0071365537
md"""
# PlutoHooks.jl

Bring your notebook to life! This is an abstraction based on [React.js Hooks](https://reactjs.org/docs/hooks-intro.html) to implement "react-like" features in [Pluto.jl](https://plutojl.org). It allows cells to carry information and processes between updates, and even update itself.

There is a lot you can do with this, but some examples:
- Run a process and relay it's output to the rest of your notebook.
- Watch a file and reload the content when it changes.
- Do a computation on separate thread while the rest of notebook continue running.

You need to use Pluto version >= 0.17.2.
"""

# ╔═╡ bc0e4219-a40b-46f5-adb2-f164d8a9bbdb
"""
	@use_memo(deps::Vector{Any}) do
		# Expensive computation/loading
	end

Does a computation only when the deps array has changed.
This is useful for heavy computations as well as resource fetches like file reading or fetching from the network.

```julia
# Only read a file once
@use_memo([filename]) do
	read(filename)
end
```

```julia
@use_memo([a, b]) do
	a + b # But they're like really big numbers
end
```
"""
macro use_memo(f, deps)
	quote
		ref = @use_ref(nothing)
		if @use_did_deps_change($(esc(deps)))
			ref[] = $(esc(f))()
		end
		ref[]
	end
end

# ╔═╡ 0f632b57-ea01-482b-b93e-d69f962a6d92
md"""
## Not really hooks but internally very hook-ish

These are all for making sure you have some level of Pluto-ness active. These are made to work outside of Pluto as well, but obviously give you the opposite results :P
"""

# ╔═╡ 8c2e9cad-eb63-4af5-8b52-629e8d3439bd
"""
	is_running_in_pluto_process()

This doesn't mean we're in a Pluto cell, e.g. can use @bind and hooks goodies.
It only means PlutoRunner is available (and at a version that technically supports hooks)
"""
function is_running_in_pluto_process()
	isdefined(Main, :PlutoRunner) &&
	# Also making sure my favorite goodies are present
	isdefined(Main.PlutoRunner, :GiveMeCellID) &&
	isdefined(Main.PlutoRunner, :GiveMeRerunCellFunction) &&
	isdefined(Main.PlutoRunner, :GiveMeRegisterCleanupFunction)
end

# ╔═╡ df0645b5-094a-45b9-b72a-ab7ef9901fa1
"""
	is_inside_pluto(mod::Module)

This can be useful to implement the behavior for when the hook is called outside Pluto but in the case where Pluto **can** be loaded.
"""
function is_inside_pluto(mod::Module)
	# Note: this could be moved to AbstractPlutoDingejtes
	startswith(string(nameof(mod)), "workspace#") &&
		isdefined(mod, Symbol("@bind"))
end

# ╔═╡ c82c8aa9-46a9-4110-88af-8638625222e3
"""
	@use_ref(initial_value::Any)::Ref{Any}

Creates a Ref that is stable over multiple **implicit** runs of the cell. Implicit run meaning a variable (or a bond) used in this cell is updated, causing this cell to re-run. When you press shift+enter in Pluto however, this ref will reset.

This is useful to keep state around between runs.
"""
macro use_ref(initial_value=nothing)
	if !is_inside_pluto(__module__)
        ref = Ref{Any}()
		return quote
            ref = $(ref)
            ref[] = $(esc(initial_value))
            ref
        end
	end

	ref_ref = Ref(Ref{Any}())

	quote
		# use_did_deps_change with empty array so it only refreshes
		# initially and on cell refresh.
		# You might wonder: But Michiel, this already refreshes on cell refresh,
		# 	because the macro will rerun!
		# I hear you. But I have bigger plans..... 😈
		if @use_did_deps_change([])
			$ref_ref[] = Ref{Any}($(esc(initial_value)))
		end
		$ref_ref[]
	end
end

# ╔═╡ 1df0a586-3692-11ec-0171-0b48a4a1c4bd
"""
	state, set_state = @use_state(initial_value::Any)

Returns a tuple for an update-able value. `state` will be whatever value you put in, and `set_state` is a function that you can call with a value, and it will set `state` to that value, and re-run the cell. Useful in combination with [`@use_effect`](@ref):

```julia
web_response = let
	state, set_state = @use_state(nothing)

	@use_effect([]) do
		schedule(Task() do
			response = HTTP.fetch("idk-what-api-HTTP.jl-has")
			set_state(response)
		end)
	end

	state
end
```

Be careful to have your [`@use_effect`](@ref) not rely on `state`, because it will most likely not have a reference to the latest state, but to the `state` at the moment that the [`@use_effect`](@ref) first ran.

To circumvent the most common case where this is a problem, you can pass a function to `set_state` that recieves the previous state as an argument:

```julia
counter = begin
	state, set_state = @use_state(0)

	@use_effect([]) do
		schedule(Task() do
			while true
				sleep(1)
				set_state(function(previous_state)
					previous_state + 1
				end)
			end
		end)

		# In the real world this should also return a cleanup function,
		# More on that in the docs for @use_effect
	end

	state
end
```
"""
macro use_state(initial_value)
	if !is_inside_pluto(__module__)
		return quote
			($(esc(initial_value)), x -> nothing)
		end
	end

	quote
		rerun_cell_fn = @give_me_rerun_cell_function()
		state_ref = @use_ref($(esc(initial_value)))

		# But there are no deps! True, but this takes care of initialization,
		# and the case that @use_deps creates, where we switch the cell_id around.
		# if @use_did_deps_change([])
		# 	state_ref[] = $(esc(initial_value))
		# end

		# TODO Make set_state throw when used after a cell is disposed
		# .... (so this would require @use_effect)
		# .... Reason I want this is because it will help a bunch in spotting
		# .... tasks that don't get killed, stuff like that.
		set_state = (new) -> begin
			new_value = if hasmethod(new, Tuple{typeof(new)})
				new(state_ref[])
			else
				new
			end

			state_ref[] = new_value
			rerun_cell_fn()
		end
		(state_ref[], set_state)
	end
end

# ╔═╡ cd048a16-37f5-455e-8b6a-c098d5f83b96
"""
	@use_deps(deps::Vector) do
		# ... others hooks ...
	end

Experimental function to wrap a bunch of macros in a fake cell that fully refreshes when the deps provided change. This is useful if you make a macro that wraps a bunch of Pluto Hooks, and you just want to refresh the whole block when something changes. This also clears [`@use_ref`](@ref)'s and [`@use_state`](@ref)'s, even though these don't even have a deps argument.

Not entirely sure how much this is necessary (or if I'm missing something obvious that doesn't make it necessary).

Also, this name does **not** spark joy.
"""
macro use_deps(fn_expr, deps)
	if !is_inside_pluto(__module__)
		return quote
			$(esc(deps))

			$(esc(fn_expr))()
		end
	end

	cell_id_ref = Ref{UUID}(uuid4())

	quote
		if @use_did_deps_change($(esc(deps)))
			$cell_id_ref[] = uuid4()
		end

		with_cell_id($(esc(fn_expr)), $cell_id_ref[])
	end
end

# ╔═╡ 89b3f807-2e24-4454-8f4c-b2a98aee571e
"""
	@use_effect(deps::Vector{Any}) do
		# Side effect
		return function()
			# Cleanup
		end
	end

Used to run a side effects that can create long running processes. A good example of this would be a HTTP server, or a task that runs an async process. Maybe it's a veeeery long running HTTP request, possibly a websocket connection to an API. 🌈

Likely to be used with [`@use_state`](@ref), as without it it's kinda useless. You want to get the values you fetch in the `@use_effect` back into the notebook, and that is what `@use_state` is for.

The function returned from `@use_effect` is called whenever the process is supposed to be stopped. This is either when the deps change, the cell is explicitly re-run or the cell is deleted. Make sure you write good cleanup functions! It's often seen as an afterthought, but it can make your notebook experience so much better.

```julia
# Ping me if you have a better real world example that uses deps
# Also don't copy this verbatim, we'll have `@use_task` that is smarter
#   with it's cleanup and returns the task state!
@use_effect([log_prefix])
	task = schedule(Task() do
		while true
			sleep(1)
			@info "
		end
	end)

	return function()
		Base.schedule(task, InterruptException(), error=true)
	end
end
```
"""
macro use_effect(f, deps)
	if !is_inside_pluto(__module__)
		return quote
			$(esc(deps))
			$(esc(f))()
		end
	end

	# For some reason the `cleanup_ref` using @use_ref or assigned outside the
	# 	`register_cleanup_fn() do ... end` (and not interpolated directly into it)
	# 	is `nothing` when the cleanup function actually ran...
	# Honestly, no idea how that works, like... `cleanup_ref[]` can be nothing sure,
	# 	but `cleanup_ref` self can be `nothing`???
	cleanup_ref = Ref{Function}(() -> nothing)
	quote
		cleanup_ref = $(cleanup_ref)

		register_cleanup_fn = @give_me_register_cleanup_function()
		register_cleanup_fn() do
			$(cleanup_ref)[]()
		end

		if @use_did_deps_change($(esc(deps)))
			cleanup_ref[]()

			local cleanup_func = $(esc(f))()
			if cleanup_func isa Function
				cleanup_ref[] = cleanup_func
			end
		end

		nothing
	end
end

# ╔═╡ 3f632c14-5f25-4426-8bff-fd315db55db5
export @use_ref, @use_state, @use_memo, @use_effect, @use_deps

# ╔═╡ c461f6da-a252-4cb4-b510-a4df5ab85065
"""
	@use_did_deps_change(deps::Vector{Any})

The most base-level `use_xxx` macro that we have, and I hope we can make it so you don't actually need this in your own code. It will, when called with deps, return `true` if the deps imply a refresh.

It will always return `true` when run the first time.
After that it will.

1. `deps=nothing` will return `true`
2. `deps=[]` will return `false`
3. `deps=[something_else...]` will return true when the deps are different than they were before
"""
macro use_did_deps_change(deps)
	if !is_inside_pluto(__module__)
		return quote
			$(esc(deps))
			true # Simulates the first run
		end
	end

	# Can't use @use_ref because this is used by @use_ref
	initialized_ref = Ref(false)
	last_deps_ref = Ref{Any}(nothing)
	last_cell_id_ref = Ref{Any}(nothing)

	quote
		initialized_ref = $(initialized_ref)
		last_deps_ref = $(last_deps_ref)
		last_cell_id_ref = $(last_cell_id_ref)

		current_deps = $(esc(deps))
		current_cell_id = @give_me_the_pluto_cell_id()

		if initialized_ref[] == false
			initialized_ref[] = true
			last_deps_ref[] = current_deps
			last_cell_id_ref[] = current_cell_id
			true
		else
			# No dependencies? Always re-render!
			if current_deps === nothing
				true
			elseif (
				# There is a problem here with either cell_id or one of the deps
				# being missing... >_> Not sure what would be good here,
				# === would fix missing, but would also make all comparisons strict.
				# Explicitly checking for missing... ? 🤮
				last_deps_ref[] == current_deps &&
				last_cell_id_ref[] == current_cell_id
			)
				false
			else
				last_deps_ref[] = current_deps
				last_cell_id_ref[] = current_cell_id
				true
			end
		end
	end
end

# ╔═╡ 84736507-7ea9-4b4b-9b70-b1e9b4b33cde
md"""
### Until I get the PlutoTest PR out
"""

# ╔═╡ 118991d7-f470-4775-ac44-4638f4989d58
md"""
## PlutoRunner-based internals

These are, I hope, the only parts that need to explicitly reference PlutoRunner.
Each of these inserts a reference to a special PlutoRunner object into the resulting expression, and that special object will be caught by PlutoRunner while evaluating the cell, and replaced with the actual value.

It seems a bit over-engineered, and I guess it is, BUT, this makes it possible to have a very strict sense of what cell is actually running what function. Also it allows other macros (specifically [`@use_deps`](@ref)) to insert it's own values instead of Plutos, thus kinda creating a cell-in-a-cell 😏
"""

# ╔═╡ 405fb702-cf4a-4d34-b8ed-d3258a61256b
const overwritten_cell_id = Ref{Union{Nothing,UUID}}(nothing)

# ╔═╡ 39aa6082-40ca-40c3-a2c0-4b6221edda32
"""
	@give_me_the_pluto_cell_id()

> ⚠ Don't use this directly!! if you think you need it, you might actually need [`@use_did_deps_change([])`](@ref) but even that is unlikely.

Used inside a Pluto cell this will resolve to the current cell UUID.
Outside a Pluto cell it will throw an error.
"""
macro give_me_the_pluto_cell_id()
	if is_running_in_pluto_process()
		:(something(overwritten_cell_id[], dont_be_pluto_special_value($(Main.PlutoRunner.GiveMeCellID()))))
	else
		:(throw(NotRunningInPlutoCellException()))
	end
end

# ╔═╡ d9d14e60-0c91-4eec-ba28-82cf1ebc115f
"""
	@use_is_pluto_cell()

Returns whether or not this expression is running inside a Pluto cell.
This goes further than checking if the process we're running in is started using Pluto, this actually checks if this code is part of the code that gets evaluated within a cell. Meant to be used directly in a Pluto cell, or returned from a macro.

This is nestable, so you can use `@use_is_pluto_cell()` inside your macro own and, as long as that macro is used in a Pluto cell directly, it will return true.

Using this inside a function will return whether or not that function is defined in a Pluto cell. If you then call that function from a script, it will still return true:

```julia
# 49cb409b-e564-47aa-9dae-9bc5bffa991d
function wrong_use_of_use_is_pluto_cell()
	return @use_is_pluto_cell()
end

# 49cb409b-e564-47aa-9dae-9bc5bffa991d
# eval circumvents Pluto-ness
eval(quote
	@use_is_pluto_cell() // false
	wrong_use_of_use_is_pluto_cell() // true
end)
```
"""
macro use_is_pluto_cell()
	# Normally you don't need this,
	# but for some reason skip_as_script seems to want it still
	var"@give_me_the_pluto_cell_id"

	give_me_cell_id = is_running_in_pluto_process() ?
		Main.PlutoRunner.GiveMeCellID() :
		nothing

	quote
		if (
			is_running_in_pluto_process() &&
			$(give_me_cell_id) != Main.PlutoRunner.GiveMeCellID()
		)
			true
		else
			false
		end
	end
end

# ╔═╡ cce13aec-7cf0-450c-bc93-bcc4e2a70dfe
"""
	@skip_as_script expr

Only run the expression if you're running inside a pluto cell. Small wrapper around [`@use_is_pluto_cell`](@ref).

"""
macro skip_as_script(expr)
	var"@use_is_pluto_cell"

	quote
		if @use_is_pluto_cell()
			$(esc(expr))
		else
			nothing
		end
	end
end

# ╔═╡ 71963fa5-82f0-4c8d-9368-0d6ba317f59e
# Notice that even though we run it in this cell's module,
# that doesn't count as "being in Pluto" enough.
@skip_as_script let
	is_pluto_cell = eval(quote
		@use_is_pluto_cell()
	end)
	if is_pluto_cell
		error("❌ eval() thinks it is a Pluto cell! What!!")
	else
		md"✅ Nice, eval() is indeed not the Pluto cell"
	end
end

# ╔═╡ ec74d9b7-b2ff-4758-a305-c3f30509a786
"""
	@only_as_script expr

Only run the expression if you're **not** running inside a pluto cell. Small wrapper around [`@use_is_pluto_cell`](@ref).

"""
macro only_as_script(expr)
	var"@use_is_pluto_cell"

	quote
		if @use_is_pluto_cell()
			nothing
		else
			$(esc(expr))
		end
	end
end

# ╔═╡ 92cfc989-5862-4314-ae1b-9cbfc4b42b40
export @use_is_pluto_cell, @skip_as_script, @only_as_script

# ╔═╡ 014d0172-3425-4429-b8d6-1d195bc60a66
@skip_as_script let
	if @use_is_pluto_cell()
		md"✅ Nice, we are indeed running in Pluto"
	else
		error("❌ Uhhhhhh")
	end
end

# ╔═╡ 3d2516f8-569e-40e4-b1dd-9f024f9266e4
"""
	@give_me_rerun_cell_function()

> ⚠ Don't use this directly!! if you think you need it, you need [`@use_state`](@ref).

Used inside a Pluto cell this will resolve to a function that, when called, will cause the cell to be re-run (in turn re-running all dependent cells).
Outside a Pluto cell it will throw an error.
"""
macro give_me_rerun_cell_function()
	if is_running_in_pluto_process()
		:(dont_be_pluto_special_value($(Main.PlutoRunner.GiveMeRerunCellFunction())))
	else
		:(throw(NotRunningInPlutoCellException()))
	end
end

# ╔═╡ cf55239c-526b-48fe-933e-9e8d56161fd6
"""
	@give_me_register_cleanup_function()

> ⚠ Don't use this directly!! if you think you need it, you need [`@use_effect`](@ref).

Used inside a Pluto cell this will resolve to a function that call be called with yet another function, and then will call that function when the cell gets explicitly re-run. ("Explicitly re-run" meaning all `@use_ref`s get cleared, for example).
Outside a Pluto cell it will throw an error.
"""
macro give_me_register_cleanup_function()
	if is_running_in_pluto_process()
		:(dont_be_pluto_special_value(
			$(Main.PlutoRunner.GiveMeRegisterCleanupFunction())
		))
	else
		:(throw(NotRunningInPlutoCellException()))
	end
end

# ╔═╡ 86a2f051-c554-4a1c-baee-8d01653c15be
"""
	with_cell_id(f, cell_id)

> ⚠ Don't use this directly!! if you think you need it, you need [`@use_deps`](@ref).

Used inside a cell to get the "proxy" cell id. This could be the real one but also a fake one in case your hook is called from another hook.
"""
function with_cell_id(f::Function, cell_id)
	previous_cell_id = overwritten_cell_id[]
	overwritten_cell_id[] = cell_id
	try
		f()
	finally
		overwritten_cell_id[] = previous_cell_id
		nothing
	end
end

# ╔═╡ b36e130e-578b-42cb-8e3a-763f6b97108d
md"""
### Very cool small helpers

These are just to make [`@give_me_the_pluto_cell_id`](@ref), [`@give_me_rerun_cell_function`](@ref) and [`@give_me_register_cleanup_function`](@ref) throw whenever you're not in Pluto.

One more reason to not call these directly.
"""

# ╔═╡ ff97bcce-1d29-469e-a4be-5dc902676057
Base.@kwdef struct NotRunningInPlutoCellException <: Exception end

# ╔═╡ 78d28d07-5912-4306-ad95-ad245797889f
function Base.showerror(io::IO, expr::NotRunningInPlutoCellException)
	print(io, "NotRunningInPlutoCell: Expected to run in a Pluto cell, but wasn't! We'll try to get these hooks to work transparently when switching from Pluto to a script.. but not yet, so just as a precaution: this error!")
end

# ╔═╡ 1b8d6be4-5ba4-42a8-9276-9ef687a8a7a3
if is_running_in_pluto_process()
	function dont_be_pluto_special_value(x::Main.PlutoRunner.SpecialPlutoExprValue)
		throw(NotRunningInPlutoCellException())
	end
end

# ╔═╡ f168c077-59c7-413b-a0ac-c0fd72781b72
dont_be_pluto_special_value(x::Any) = x

# ╔═╡ 9ec6b9c5-6bc1-4033-ab93-072f783184e9
md"""
### Until I get the PlutoTest PR out
"""

# ╔═╡ fd653af3-be53-4ddd-b69d-3967ef6d588a
md"#### `@give_me_the_pluto_cell_id()`"

# ╔═╡ b25ccaf1-cf46-4eea-a4d9-16c68cf56fad
@skip_as_script try
	eval(quote
		@give_me_the_pluto_cell_id()
	end)
	error("❌ This should throw a NotRunningInPlutoCellException.. but didn't!")
catch e
	if e isa NotRunningInPlutoCellException
		md"✅ Nice, we got an exception like we should"
	else
		rethrow(e)
	end
end

# ╔═╡ e5905d1e-33ec-47fb-9f98-ead82eb03be8
@skip_as_script begin
	cell_id = @give_me_the_pluto_cell_id()
	if cell_id isa UUID
		md"✅ Nice, we got the cell UUID"
	else
		error("❌ What the? Got a typeof($(typeof(cell_id)))")
	end
end

# ╔═╡ 274c2be6-6075-45cf-b28a-862c8bf64bd4
md"""
## Examples/Experiments

Ideally, these functions would be in their own package (so they can update without PlutoHooks updating), but for now we keep them here to show of and test the stuff above.
"""

# ╔═╡ 90f051be-4384-4383-9a56-2aa584687dc3
macro use_reducer(fn, deps=nothing)
	quote
		ref = @use_ref(nothing)

		current_value = ref[]
		if @use_did_deps_change($(esc(deps)))
			next_value = $(esc(fn))(current_value)
			ref[] = next_value
		end

		ref[]
	end
end

# ╔═╡ c8c560bf-3ef6-492f-933e-21c898fb2db6
md"### `@use_task`"

# ╔═╡ 9ec99592-955a-41bd-935a-b34f37bb5977
macro use_task(f, deps)
	quote
		error("@use_task was moved to PlutoLinks.jl")
	end
end

# ╔═╡ 56f2ff19-c6e8-4858-8e6a-3b790fae7ecb
md"### `@use_file(filename)`"

# ╔═╡ e240b167-560c-4dd7-9801-30467d8758be
macro use_file_change(filename)
	quote
		error("@use_file_change was moved to PlutoLinks.jl")
	end
end

# ╔═╡ 461231e8-4958-46b9-88cb-538f9151a4b0
macro use_file(filename)
	quote
		error("@use_file was moved to PlutoLinks.jl")
	end
end

# ╔═╡ 9af74baf-6571-4a0c-b0c0-989472f18f7a
md"### `@ingredients(julia_file_path)`"

# ╔═╡ d84f47ba-7c18-4d6c-952c-c9a5748a51f8
macro ingredients(filename)
	quote
		error("@ingredients was moved to PlutoLinks.jl")
	end
end

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[Random]]
deps = ["SHA", "Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
"""

# ╔═╡ Cell order:
# ╟─729ae3bb-79c2-4fcd-8645-7e0071365537
# ╠═49cb409b-e564-47aa-9dae-9bc5bffa991d
# ╠═3f632c14-5f25-4426-8bff-fd315db55db5
# ╠═92cfc989-5862-4314-ae1b-9cbfc4b42b40
# ╟─c82c8aa9-46a9-4110-88af-8638625222e3
# ╟─1df0a586-3692-11ec-0171-0b48a4a1c4bd
# ╟─cd048a16-37f5-455e-8b6a-c098d5f83b96
# ╟─89b3f807-2e24-4454-8f4c-b2a98aee571e
# ╟─bc0e4219-a40b-46f5-adb2-f164d8a9bbdb
# ╟─c461f6da-a252-4cb4-b510-a4df5ab85065
# ╟─0f632b57-ea01-482b-b93e-d69f962a6d92
# ╟─d9d14e60-0c91-4eec-ba28-82cf1ebc115f
# ╟─cce13aec-7cf0-450c-bc93-bcc4e2a70dfe
# ╟─ec74d9b7-b2ff-4758-a305-c3f30509a786
# ╟─8c2e9cad-eb63-4af5-8b52-629e8d3439bd
# ╟─df0645b5-094a-45b9-b72a-ab7ef9901fa1
# ╟─84736507-7ea9-4b4b-9b70-b1e9b4b33cde
# ╟─014d0172-3425-4429-b8d6-1d195bc60a66
# ╟─71963fa5-82f0-4c8d-9368-0d6ba317f59e
# ╟─118991d7-f470-4775-ac44-4638f4989d58
# ╟─405fb702-cf4a-4d34-b8ed-d3258a61256b
# ╟─39aa6082-40ca-40c3-a2c0-4b6221edda32
# ╟─3d2516f8-569e-40e4-b1dd-9f024f9266e4
# ╟─cf55239c-526b-48fe-933e-9e8d56161fd6
# ╟─86a2f051-c554-4a1c-baee-8d01653c15be
# ╟─b36e130e-578b-42cb-8e3a-763f6b97108d
# ╠═ff97bcce-1d29-469e-a4be-5dc902676057
# ╟─78d28d07-5912-4306-ad95-ad245797889f
# ╟─1b8d6be4-5ba4-42a8-9276-9ef687a8a7a3
# ╟─f168c077-59c7-413b-a0ac-c0fd72781b72
# ╟─9ec6b9c5-6bc1-4033-ab93-072f783184e9
# ╟─fd653af3-be53-4ddd-b69d-3967ef6d588a
# ╟─b25ccaf1-cf46-4eea-a4d9-16c68cf56fad
# ╟─e5905d1e-33ec-47fb-9f98-ead82eb03be8
# ╟─274c2be6-6075-45cf-b28a-862c8bf64bd4
# ╟─90f051be-4384-4383-9a56-2aa584687dc3
# ╟─c8c560bf-3ef6-492f-933e-21c898fb2db6
# ╠═9ec99592-955a-41bd-935a-b34f37bb5977
# ╟─56f2ff19-c6e8-4858-8e6a-3b790fae7ecb
# ╠═e240b167-560c-4dd7-9801-30467d8758be
# ╠═461231e8-4958-46b9-88cb-538f9151a4b0
# ╟─9af74baf-6571-4a0c-b0c0-989472f18f7a
# ╠═d84f47ba-7c18-4d6c-952c-c9a5748a51f8
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
module Baskets

using ..CurrenciesBase
import Base: +, -, *, /, ==

export Basket

include("impl.jl")
include("type.jl")
include("mixed.jl")

end