text:
@cache mutable struct DImplicitEulerCache{uType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  uprev2::uType
  tmp::uType
  atmp::uNoUnitsType
end

mutable struct DImplicitEulerConstantCache <: OrdinaryDiffEqConstantCache end

function alg_cache(alg::DImplicitEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{false})
  DImplicitEulerConstantCache()
end

function alg_cache(alg::DImplicitEuler,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,
                   tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Val{true})
  atmp = similar(u,uEltypeNoUnits)
  tmp = similar(u)

  DImplicitEulerCache(u,uprev,uprev2,tmp,atmp)
end


text:
# ---
# title: 54. Spiral Matrix
# id: problem54
# author: Qling
# date: 2021-03-15
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/spiral-matrix/description/>
# hidden: true
# ---
#
# Given an `m x n` `matrix`, return _all elements of the_ `matrix` _in spiral
# order_.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/13/spiral1.jpg)
#
#
#
#     Input: matrix = [[1,2,3],[4,5,6],[7,8,9]]
#     Output: [1,2,3,6,9,8,7,4,5]
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/13/spiral.jpg)
#
#
#
#     Input: matrix = [[1,2,3,4],[5,6,7,8],[9,10,11,12]]
#     Output: [1,2,3,4,8,12,11,10,9,5,6,7]
#
#
#
#
# **Constraints:**
#
#   * `m == matrix.length`
#   * `n == matrix[i].length`
#   * `1 <= m, n <= 10`
#   * `-100 <= matrix[i][j] <= 100`
#
#
## @lc code=start
using LeetCode

function spiral_order(matrix::Vector{Vector{Int}})::Vector{Int}
    res = Int[]
    while !isempty(matrix)
        for each in popfirst!(matrix)
            push!(res, each)
        end

        if !isempty(matrix) && !isempty(matrix[1])
            for row in matrix
                push!(res, pop!(row))
            end
        end

        if !isempty(matrix)
            for each in pop!(matrix)[end:-1:1]
                push!(res, each)
            end
        end

        if !isempty(matrix) && !isempty(matrix[1])
            for row in matrix[end:-1:1]
                push!(res, popfirst!(row))
            end
        end
    end

    return res
end
## @lc code=end

ext: jl


text:
struct ChineseTree
    label::String
    adj::Vector{ChineseTree}
    # prob::Float64
end

function label(s::String)
    t = findfirst(!isletter, s)
    return (t == nothing || t == 1) ? s : s[1:t - 1]
end

isleaf(tree::ChineseTree) = length(tree.adj) == 0
isposn(tree::ChineseTree) = length(tree.adj) == 1 && isleaf(tree.adj[1])

function dfstraverse(tree::ChineseTree, visitor::Function)
    for chd in tree.adj
        dfstraverse(chd, visitor)
    end
    visitor(tree)
    return nothing
end

function ChineseTree(chars::Vector{Char}; l = findfirst(isequal('('), chars), trim = label)
    nchar = length(chars)
    l += 1
    r = l
    while chars[r] != '(' && chars[r] != ')' r += 1 end
    if (chars[r] == ')')
        ss = split(join(chars[l:r - 1]))
        @assert length(ss) == 2
        leaf = ChineseTree(String(ss[2]), ChineseTree[])
        posn = ChineseTree(trim(String(ss[1])), ChineseTree[leaf])
        return posn
    else
        @assert chars[r] == '('
        ss = split(join(chars[l:r - 1]))
        # @assert length(ss) == 1 string(l, "  ", chars)
        if length(ss) == 0
            return ChineseTree(chars; l = r)
        end
        @assert length(ss) == 1
        fa = ChineseTree(trim(String(ss[1])), ChineseTree[])
        # setlabel!(obj, cur, ss[1])
        nlb = 1
        while nlb > 0
            if (chars[r] == '(')
                if nlb == 1 l = r end
                nlb += 1
            elseif chars[r] == ')'
                nlb -= 1
                if nlb == 1
                    push!(fa.adj, ChineseTree(chars; l = l, trim = trim))
                end
            end
            r += 1
        end
        return fa
    end
end

function tochars(s::String)
    cps = Vector{Char}(undef, length(s))
    for (i, cp) in enumerate(s)
        cps[i] = cp
    end
    return cps
end
ChineseTree(ct::String; trim = label) = ChineseTree(tochars(ct); trim = trim)

import Base.size
function size(tree::ChineseTree)
    if isleaf(tree) return (1, 1)
    else
        h = w = 0
        for chd in tree.adj
            ch, cw = size(chd)
            h += ch; w = max(w, cw + 1)
        end
        return (h, w)
    end
end

import Base.display

function display(obj::ChineseTree)
    xlim, ylim = size(obj)
    mat = fill("", xlim, ylim)
    ileaf = 1
    colors = Dict{String, Int}()
    edges = Vector{Tuple{Int, Int, Int, Int}}()
    function dfs(cur::ChineseTree)
        if !haskey(colors, cur.label)
            colors[cur.label] = (1 + (length(colors) + 1) * 10) % 256
        end
        if isleaf(cur)
            mat[ileaf, ylim] = cur.label
            ileaf += 1
            return (ileaf - 1, ylim)
        else
            cxys = Vector{Tuple{Int, Int}}()
            x = xlim; y = ylim;
            for chd in cur.adj
                cx, cy = dfs(chd)
                x = min(x, cx)
                y = min(y, cy - 1)
                push!(cxys, (cx, cy))
            end
            for cxy in cxys push!(edges, (x, y, cxy[1], cxy[2])) end
            mat[x, y] = cur.label
            return (x, y)
        end
    end
    dfs(obj);
    for (lx, ly, rx, ry) in edges
        for x = lx + 1:rx - 1
            new = old = mat[x, ly]
            if old == "" new = "│" end
            if old == "└" new = "├" end
            mat[x, ly] = new
        end
        if lx < rx mat[rx, ly] = "└" end
        mat[rx, ly + 1:ry - 1] .= "─";
    end
    yw = mapreduce(length, max, mat; dims = 1) .+ 2
    for x = 1:xlim, y = 1:ylim
        s = mat[x, y]; ns = length(s)
        # words
        if y == ylim println(" ", s); continue end
        padl = div(yw[y] - ns, 2)
        if (y == ylim - 1)
            padl = yw[y] - ns
        end
        padr = yw[y] - ns - padl
        padlc = padrc = '─'
        if s == "" || s == "│"
            padlc = padrc = ' '
        elseif s == "└" || s == "├" || x == 1 && y == 1
            padlc = ' '
        end
        ns = padlc ^ padl * s * padrc ^ padr
        for i = 1:padl print(padlc) end
        if haskey(colors, s) printstyled(s; color = colors[s]) else print(s) end
        for i = 1:padr print(padrc) end
    end
    return nothing
end

mutable struct Block
    chars::Vector{Char}
    nlb::Int
end
Block() = Block(Vector{Char}(), 0)
import Base.push!

function push!(block::Block, chars::Vector{Char})
    for c in chars
        if c == '('
            block.nlb += 1
        elseif c == ')'
            block.nlb -= 1
        end
        push!(block.chars, c)
    end
end

function ok(block::Block)
    return block.nlb == 0
end

const ChineseTreeBank = Vector{Vector{ChineseTree}}

function parsectb(;ctbrootpath::String = "/Users/dy/ChineseTokenizers.jl/data/ctb8.0", trim = label)
    olddir = pwd()
    newdir = joinpath(realpath(ctbrootpath), "data", "bracketed")
    res = Vector{Vector{ChineseTree}}()
    cd(() -> begin
    for (i, f) in enumerate(readdir("."))
        resf = Vector{ChineseTree}()
        push!(res, resf)
        block = Block()
        open(f, "r") do io
            for line in eachline(io)
                if startswith(line, "(")
                    push!(block, tochars(line))
                    if !ok(block)
                        for line in eachline(io)
                            push!(block, tochars(line))
                            if ok(block) break end
                        end
                    end
                    push!(resf, ChineseTree(block.chars; trim = trim))
                    block = Block()
                end
            end
        end
    end
end, newdir)
    return res
end

# ctb = parsectb("data/ctb8.0")

# using DataFrames
# import Base.stat
# function stat(ctb::ChineseTreeBank)
#     from = Vector{String}(); to = Vector{Vector{String}}()
#     from2 = Vector{String}(); to2 = Vector{String}()
#     function visitor(tree, cur)
#         nchds = length(tree.adj[cur])
#         if nchds > 0
#             if nchds == 1
#                 chd = tree.adj[cur][1]
#                 if length(tree.adj[chd]) == 0
#                     if tree.label[cur] != "-NONE-"
#                         push!(from2, tree.label[cur])
#                         push!(to2, tree.label[chd])
#                     end
#                 elseif tree.label[chd] != "-NONE-"
#                     push!(from, tree.label[cur])
#                     push!(to, tree.label[tree.adj[cur]])
#                 end
#             else
#                 push!(from, tree.label[cur])
#                 push!(to, tree.label[tree.adj[cur]])
#             end
#         end
#     end
#     for vec in ctb
#         for tree in vec
#             dfstraverse(tree, visitor)
#         end
#     end
#     function f(df::DataFrame)
#         return by(df, [:from, :to], tot = :from => length, sort = true)
#     end
#     inn = f(DataFrame(from = from, to = to))
#     pos = f(DataFrame(from = from2, to = to2))
#     inn, pos
# end

function posword(ct::ChineseTree)
    res = Vector{Tuple{String, String}}()
    function visitor(cur)
        if isposn(cur) push!(res, (cur.label, cur.adj[1].label)) end
    end
    dfstraverse(ct, visitor)
    return res
end

function writeposword(filepath::String, ctb::ChineseTreeBank)
    open(filepath, "w") do io
        for ctv in ctb
            for tree in ctv
                res = filter(x -> x[1] != "-NONE-", posword(tree))
                println(io, length(res))
                for r in res
                    println(io, r[1], " ", r[2])
                end
            end
        end
    end
end
# ChineseTreeNode = Union{InnerTreeNode, PosTreeNode, LeafTreeNode}
# ChineseTreeNode(label::String, id::String) == InnerTreeNode(label, id, Int[], 0.0)

function cnf(root::ChineseTree)::ChineseTree
    nchd = length(root.adj)
    if nchd == 0
        return root
    elseif nchd <= 2
        newroot = ChineseTree(root.label, ChineseTree[])
        for i = 1:nchd push!(newroot.adj, cnf(root.adj[i])) end
        return newroot
    else
        newroot = ChineseTree(root.label, ChineseTree[])
        newright = ChineseTree(join(map(x -> x.label, root.adj[2:end]), "+"), root.adj[2:end])
        push!(newroot.adj, cnf(root.adj[1]))
        push!(newroot.adj, cnf(newright))
        return newroot
    end
end

function decnf(root::ChineseTree)
    nodes = decnf2(root)
    return length(nodes) == 1 ? nodes[1] : ChineseTree(join(map(x -> x.label, nodes), "+"), nodes)
end

function decnf2(root::ChineseTree)::Vector{ChineseTree}
    nchd = length(root.adj)
    if nchd == 0 return ChineseTree[root] end
    if nchd == 1
        newroot = ChineseTree(root.label, decnf2(root.adj[1]))
        return ChineseTree[newroot]
    end
    newroot = ChineseTree(root.label, append!(decnf2(root.adj[1]), decnf2(root.adj[2])))
    istmp = in('+', root.label)
    return istmp ? newroot.adj : ChineseTree[newroot]
end

import Base.==
function ==(obj1::ChineseTree, obj2::ChineseTree)
    res = obj1.label == obj2.label
    if !res return false end
    if length(obj1.adj) != length(obj2.adj) return false end
    for i in 1:length(obj1.adj)
        if obj1.adj[i] != obj2.adj[i] return false end
    end
    return true
end


text:
using XVGReader
using Test

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


text:
cubox(p::CudaPtr) = cubox(p.ptr)
cubox(a::CdArray) = cubox(rawpointer(a))
cubox{T}(x::T) = T[x]

function launch(f::CuFunction, grid::CudaDims, block::CudaDims, args::Tuple; shmem_bytes::Int=4, stream=null_stream)
    gx = get_dim_x(grid)
    gy = get_dim_y(grid)
    gz = get_dim_z(grid)

    tx = get_dim_x(block)
    ty = get_dim_y(block)
    tz = get_dim_z(block)

    kernel_args = [cubox(arg) for arg in args]

    checkdrv(ccall((:cuLaunchKernel, libcuda), Cint, (
        Ptr{Void},  # function
        Cuint,  # grid dim x
        Cuint,  # grid dim y
        Cuint,  # grid dim z
        Cuint,  # block dim x
        Cuint,  # block dim y
        Cuint,  # block dim z
        Cuint,  # shared memory bytes,
        Ptr{Void}, # stream
        Ptr{Ptr{Void}}, # kernel parameters,
        Ptr{Ptr{Void}}), # extra parameters
        f.handle, gx, gy, gz, tx, ty, tz, shmem_bytes, stream, kernel_args, 0))
end

# Note this is asynchronous wrt the host, but you can wait on the stream
function cudasleep(secs; dev::Integer=device(), stream=null_stream)
    device(dev)
    rate = attribute(dev, rt.cudaDevAttrClockRate)
    func = ptxdict[(dev, "clock_block")]
    twatchdog = 1.95  # watchdog timer kicks in after 2 secs
    while secs > 0
        tics = int64(1000*rate*min(twatchdog, secs))  # rate is in kHz
        secs -= twatchdog
        launch(func, 1, 1, (tics,), stream=stream)
    end
end


text:
"""
  outer!(z::Vector{Float}, a::Vector{Float}, b::Vector{Float})

Return z - Vector{Float} , vectorsed outer/kroneker product o vectors a and b
Auxiliary function for rawmoment
"""
function outer!(z::Vector{T}, a::Vector{T}, b::Vector{T}) where T<: AbstractFloat
    sa = size(a,1)
    sb = size(b,1)
    m = 1
    for i = 1:sb for j = 1:sa
        @inbounds z[m] = a[j]*b[i]
        m += 1
    end
    end
end

"""
  updvec!(A::Vector{Float}, B::Vector{Float})

Returns updated Vector{Float} A, by adding elementwisely Vector{Float} B
Auxiliary function for rawmoment
"""
function updvec!(A::Vector{T}, B::Vector{T}) where T<: AbstractFloat
  n = size(A, 1)
  for i=1:n
    @inbounds A[i] += B[i]
  end
  return A
end

"""
  rawmoment(X::Matrix{T}, m::Int = 4)

Simmilar to raw_moments_upto_p in R, does not expoloit tensor's symmetry
pyramid structures and blocks

Returns Array{Float, m}, the m'th moment's tensor
"""
function rawmoment(X::Matrix{T}, m::Int = 4) where T<: AbstractFloat
  t,n = size(X)
  if m == 1
    return mean(X, dims=1)[1,:]
  else
    z = [map(i -> zeros(T, n^i), 1:m)...]
    y = zeros(T, n^m)
    for i in 1:t
      xi = X[i, :]
      z[1] = xi
      for j in 2:m
        outer!(z[j], xi, z[j-1])
      end
      updvec!(y, z[m])
    end
  end
  reshape(y/t, fill(n, m)...)
end

"""
  raw_moments_upto_k(X::Matrix, k::Int = 4)

Returns [Array{Float, 1}, ..., Array{Float, k}] noncentral moment tensors of
order 1, ..., k
"""
raw_moments_upto_k(X::Matrix{T}, k::Int = 4) where T<: AbstractFloat =
  [rawmoment(X, i) for i in 1:k]

"""
  cumulants_from_moments(raw::Vector{Array{Float, i}, m = 1:k})

Returns [Array{Float, 1}, ..., Array{Float, k}] cumulant tensors of order 1, ..., k

Uses relation between cumulants and multivariate moments from e.g.
c_{ijkl} = m_{ijkl} - m_{ijk} m_{l} [4] - m_{ij} m_{kl} [3] + 2 m_{ij} m_{k} m_{l}
-6 m_{i} m_{j} m_{k} m_{l}
"""
function cumulants_from_moments(raw::Vector{Array{T}}) where T<: AbstractFloat
  k = length(raw)
  cumarr = Array{Array{Float64}}(undef, k)
  for j in 1:k
    dimr = size(raw[j])
    cumu = zeros(Float64, dimr)
    ldim = length(dimr)
    spp = collect(partitions(1:ldim))
    qpp = [(-1)^i*factorial(i) for i in 0:(ldim-1)]
    sppl = [map(length, spp[i]) for i in 1:length(spp)]
    for i in 1:prod(dimr)
      @inbounds ind = Tuple(CartesianIndices(dimr)[i])
      @inbounds cumu[ind...] = onecumulant(ind, raw, spp, sppl, qpp)
    end
    cumarr[j] = cumu
  end
  cumarr
end

"""
  onecumulant(ind::Tuple, raw::Vector{Array}, spp::Vector, sppl::Vector{Vector}, dpp::Vector)

raw - vector of moment's tensors, spp - vector of partitions, sppl - vector of sizes of partitions
dpp - vector of a factor for each product of moments (a factorial factor).

Returns Array{Float, n} the n'th cumulant tensor

"""
function onecumulant(ind::Tuple, raw::Vector{Array{T}}, spp::Vector,
          sppl::Vector{Vector{Int}}, dpp::Vector{Int}) where T<: AbstractFloat
  ret = zero(T)
  for i in 1:length(spp)
    part = spp[i]
    beln = length(part)
    k = sppl[i]
    temp = one(T)
    for r in 1:beln
      temp *= raw[k[r]][ind[part[r]]...]
    end
    ret += dpp[beln]*temp
  end
  ret
end

"""
  cumulatsfrommoments(x::Matrix{Float}, k::Int)

Returns a vector of 1,2, .., k dims Arrays{Float} of cumulant tensors
"""
mom2cums(x::Matrix{T}, k::Int) where T<: AbstractFloat =
  cumulants_from_moments(raw_moments_upto_k(x, k))


text:
# brick wall scanning for a single MERA layer with treatment to the tail
function correlation_matrix_to_gmps_brickwall_tailed(
  Λ0::AbstractMatrix{ElT},
  inds::Vector{Int};
  eigval_cutoff::Float64=1e-8,
  maxblocksize::Int=size(Λ0, 1),
) where {ElT<:Number}
  Λ = Hermitian(Λ0)
  N = size(Λ, 1)
  V = Circuit{ElT}([])
  #ns = Vector{real(ElT)}(undef, 2*N)
  err_tot = 0.0
  indsnext = Int[]
  relinds = Int[]
  for i in 1:N
    if i % 2 == 0
      append!(indsnext, inds[i])
      append!(relinds, i)
      continue
    end
    blocksize = 0
    n = 0.0
    err = 0.0
    p = Int[]
    uB = 0.0
    # find the block whose lowest eigenvalue is within torelence
    for blocksize in 1:maxblocksize
      j = min(i + blocksize, N)
      ΛB = deepcopy(Λ[i:j, i:j]) #@view Λ[i:j, i:j] # \LambdaB is still part of Lambda
      nB, uB = eigen(Hermitian(ΛB))
      # sort by -(n * log(n) + (1 - n) * log(1 - n)) in ascending order
      p = sortperm(nB; by=entropy)
      n = nB[p[1]]
      err = min(n, 1 - n)
      err ≤ eigval_cutoff && break
    end
    # keep the node if the err cannot be reduced
    if i + maxblocksize >= N && err > eigval_cutoff
      append!(indsnext, inds[i])
      append!(relinds, i)
      continue
    end
    err_tot += err
    #ns[i] = n # eigenvalue
    v = deepcopy(uB[:, p[1]]) #@view uB[:, p[1]] # eigenvector of the correlation matrix
    g, _ = givens_rotations(v) # convert eigenvector into givens rotation
    shift!(g, i - 1) # shift rotation location
    # In-place version of:
    # V = g * V
    lmul!(g, V)
    #@show g
    Λ = Hermitian(g * Λ * g') #isolate current site i
  end
  return Λ, V, indsnext, relinds
end

# shift givens rotation indexes according to the inds
function shiftByInds!(G::Circuit, inds::Vector{Int})
  for (n, g) in enumerate(G.rotations)
    G.rotations[n] = Givens(inds[g.i1], inds[g.i2], g.c, g.s)
  end
  return G
end

"""
    correlation_matrix_to_gmera(Λ::AbstractMatrix{ElT}; eigval_cutoff::Float64 = 1e-8, maxblocksize::Int = size(Λ0, 1))
Diagonalize a correlation matrix through MERA layers,
output gates and eigenvalues of the correlation matrix
"""
# Combine gates for each MERA layer
function correlation_matrix_to_gmera(
  Λ0::AbstractMatrix{ElT}; eigval_cutoff::Float64=1e-8, maxblocksize::Int=size(Λ0, 1)
) where {ElT<:Number}
  Λ = Hermitian(Λ0)
  N = size(Λ, 1)
  Nnew = N - 1
  inds = collect(1:N)
  V = Circuit{ElT}([])
  Λtemp = deepcopy(Λ)
  layer = 0                                 # layer label of MERA
  while N > Nnew # conditioned on the reduction of nodes
    N = Nnew
    # indsnext: next layer indexes with original matrix labels
    # relinds: next layer indexes with labels from the last layer
    Λr, C, indsnext, relinds = correlation_matrix_to_gmps_brickwall_tailed(
      Λtemp, inds; eigval_cutoff=eigval_cutoff, maxblocksize=maxblocksize
    )
    shiftByInds!(C, inds)    # shift the index back to the original matrix
    inds = indsnext
    Λtemp = deepcopy(Λr[relinds, relinds]) # project to even site for next layer based on keeping indexes relinds
    Nnew = size(Λtemp, 1)
    lmul!(C, V)                           # add vector of givens rotation C into the larger vector V
    #V = C * V
    layer += 1
    #Λ = ITensors.Hermitian(C * Λ * C')
  end
  # gmps for the final layer
  Λr, C = correlation_matrix_to_gmps(
    Λtemp; eigval_cutoff=eigval_cutoff, maxblocksize=maxblocksize
  )
  shiftByInds!(C, inds)
  lmul!(C, V)
  Λ = V * Λ0 * V'
  ns = real.(diag(Λ))
  return ns, V
end

# output the MERA gates and eigenvalues of correlation matrix from WF
function slater_determinant_to_gmera(Φ::AbstractMatrix; kwargs...)
  return correlation_matrix_to_gmera(conj(Φ) * transpose(Φ); kwargs...)
end

# ouput the MPS based on the MERA gates
function correlation_matrix_to_mera(
  s::Vector{<:Index},
  Λ::AbstractMatrix;
  eigval_cutoff::Float64=1e-8,
  maxblocksize::Int=size(Λ, 1),
  kwargs...,
)
  @assert size(Λ, 1) == size(Λ, 2)
  ns, C = correlation_matrix_to_gmera(
    Λ; eigval_cutoff=eigval_cutoff, maxblocksize=maxblocksize
  )
  if all(hastags("Fermion"), s)
    U = [ITensor(s, g) for g in reverse(C.rotations)]
    ψ = MPS(s, n -> round(Int, ns[n]) + 1, U; kwargs...)
  elseif all(hastags("Electron"), s)
    isodd(length(s)) && error(
      "For Electron type, must have even number of sites of alternating up and down spins.",
    )
    N = length(s)
    if isspinful(s)
      error(
        "correlation_matrix_to_mps(Λ::AbstractMatrix) currently only supports spinless Fermions or Electrons that do not conserve Sz. Use correlation_matrix_to_mps(Λ_up::AbstractMatrix, Λ_dn::AbstractMatrix) to use spinful Fermions/Electrons.",
      )
    else
      sf = siteinds("Fermion", 2 * N; conserve_qns=true)
    end
    U = [ITensor(sf, g) for g in reverse(C.rotations)]
    ψf = MPS(sf, n -> round(Int, ns[n]) + 1, U; kwargs...)
    ψ = MPS(N)
    for n in 1:N
      i, j = 2 * n - 1, 2 * n
      C = combiner(sf[i], sf[j])
      c = combinedind(C)
      ψ[n] = ψf[i] * ψf[j] * C
      ψ[n] *= δ(dag(c), s[n])
    end
  else
    error("All sites must be Fermion or Electron type.")
  end
  return ψ
end

function slater_determinant_to_mera(s::Vector{<:Index}, Φ::AbstractMatrix; kwargs...)
  return correlation_matrix_to_mera(s, conj(Φ) * transpose(Φ); kwargs...)
end

# G the circuit from the gates, N is the total number of sites
function UmatFromGates(G::Circuit, N::Int)
  U = Matrix{Float64}(I, N, N)
  n = size(G.rotations, 1)
  for k in 1:n
    rot = G.rotations[k]
    U = rot * U
  end
  return U
end

# compute the energy of the state based on the gates
function EfromGates(H::Matrix{<:Number}, U::Matrix{<:Number})
  Htemp = U * H * U'
  Etot = 0
  N = size(U, 1)
  for i in 1:N
    if Htemp[i, i] < 0.0
      Etot += Htemp[i, i]
    end
  end
  return Etot
end


text:
@testset "Friction factor" begin

  @test_throws DomainError FrictionFactor(ReynoldsNumber(-10),RoughnessRatio(0.001))

  @test_throws DomainError FrictionFactor(ReynoldsNumber(10),RoughnessRatio(-0.001))

  @test FrictionFactor(ReynoldsNumber(1e5),RoughnessRatio(0.002)) < 0.035

end

max_stars_repo_name: UCLAMAEThreads/ThermofluidQuantities

max_issues_repo_name: UCLAMAEThreads/ThermofluidQuantities

max_forks_repo_name: UCLAMAEThreads/ThermofluidQuantities


text:
# Implements Armstrong's CitZAF algorithm as presented in the Green Book (Heinrich, 1991)

struct CitZAF <: MatrixCorrection
    material::Material
    E0::Float64  # Beam energy
    subshell::Union{Nothing,AtomicSubShell}
    Ea::Float64 # Sub-shell energy
    γ0::Float64
    q::Float64
    α::Float64
    β::Float64

    function CitZAF(mat::Material, ashell::Union{AtomicSubShell,Nothing}, ea::Float64, e0::Float64)
        function ϕ0L(u0, ηbar)
            a = 3.43378 + (-10.7872 + (10.97628 - 3.62286 / u0) / u0) / u0
            b = -0.59299 + (21.55329 + (-30.55248 + 9.59218 / u0) / u0) / u0
            return 1.0 + (ηbar / (1.0 + ηbar)) * (a + b * log(1.0 + ηbar))
        end
        function ηbarLS(z, e0k)
            h1 = 1.0e-4 * (-52.3791 + z * (150.48371 + z * (-1.67373 + z * 0.00716)))
            h2 = 1.0e-4 * (-1112.8 + z * (30.289 - z * 0.15498))
            return h1 * (1.0 + h2 * log(e0k / 20.0))
        end
        e0k, eck, u0 = 0.001 * e0, 0.001 * ea, e0 / ea
        @assert u0 >= 1.0
        zbar =
            sum(value(mat[elm]) * z(elm) / a(elm, mat) for elm in keys(mat)) /
            sum(value(mat[elm]) / a(elm, mat) for elm in keys(mat))
        abar = sum(value(mat[elm]) for elm in keys(mat)) / sum(value(mat[elm]) / a(elm, mat) for elm in keys(mat))
        γ0 = (5π * u0) / ((u0 - 1.0) * log(u0)) * (log(u0) - 5.0 + 5.0 * u0^-0.2)
        ηbar = mapreduce(elm -> mat[elm] * ηbarLS(z(elm), e0k), +, keys(mat))
        ϕ0 = ϕ0L(u0, ηbar) # 1.0 + 2.8*(1.0 -  0.9/u0)*ηbar
        q = (γ0 - ϕ0) / γ0
        α = 2.97e5 * (zbar^1.05 / (abar * e0k^1.25)) * sqrt(log(1.166 * e0 / J(XPP, mat)) / (e0k - eck))
        β = (8.5e5 * zbar^2) / (abar * e0k^2 * (γ0 - 1.0))
        return new(mat, e0, ashell, ea, γ0, q, α, β)
    end
end

ϕ(cz::CitZAF, ρz) = cz.γ0 * exp(-(cz.α * ρz)^2) * (1.0 - cz.q * exp(-cz.β * ρz))

F(cz::CitZAF) = ((cz.α - cz.q * cz.α + cz.β) * cz.γ0) / (cz.α * (cz.α + cz.β))

function Fχ(cz::CitZAF, θtoa::Float64)
    @assert isnothing(cz.subshell)  "Use only for continuum correction"
    χm = χ(material(cz), cz.Ea, θtoa)
    return cz.γ0 * (1.0 / (cz.α + χm) - cz.q / (cz.α + cz.β + χm))
end

function Fχ(cz::CitZAF, xray::CharXRay, θtoa::Float64)
    @assert !isnothing(cz.subshell) "Use only for characteristic correction"
    @assert inner(xray) == cz.subshell
    χm = χ(material(cz), xray, θtoa)
    return cz.γ0 * (1.0 / (cz.α + χm) - cz.q / (cz.α + cz.β + χm))
end

function Fχp(cz::CitZAF, xray::CharXRay, θtoa::Float64, τ::Float64)
    @assert isnothing(cz.subshell) || (inner(xray) == cz.subshell)
    χm = χ(material(cz), xray, θtoa)
    return ((1.0 - exp(-τ * (cz.α + χm))) * cz.γ0) / (cz.α + χm) +
           ((-1.0 + exp(-τ * (cz.α + cz.β + χm))) * cz.q * cz.γ0) / (cz.α + cz.β + χm)
end

matrixcorrection(::Type{CitZAF}, mat::Material, ashell::AtomicSubShell, e0::Float64) = CitZAF(mat, ashell, energy(ashell), e0)

continuumcorrection(::Type{CitZAF}, mat::Material, ea::Float64, e0::Float64) = CitZAF(mat, nothing, ea, e0)

Base.range(::Type{CitZAF}, mat::Material, e0::Float64, inclDensity=true) = range(XPP, mat, e0, inclDensity)

NeXLCore.minproperties(::Type{CitZAF}) = (:BeamEnergy, :TakeOffAngle)

max_stars_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_issues_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_forks_repo_name: NicholasWMRitchie/NeXLMatrixCorrection


text:
function get_fitness(ea::Function, ctype::DataType)
    try
        maxfit, best = ea(ctype, nin, nout, fit)
        return -maxfit
    catch
        Logging.info(@sprintf("Error: %s\nE:%s %s %s", string(catch_stacktrace()),
                              string(ea), string(ctype), CGP.Config.to_string()))
        return Inf
    end
end

function oneplus_config(cfg::Array{Float64}, ctype::DataType)
    cfg = mod.(cfg, 1.0)
    mut = CGP.Config.index_in(CGP.mutations, cfg[1])
    CGP.Config.init(Dict("mutate_method" => string("\"", mut, "\""),
                         "recurrency" => cfg[2],
                         "input_mutation_rate" => cfg[3],
                         "output_mutation_rate" => cfg[4],
                         "node_mutation_rate" => cfg[5],
                         "add_node_rate" => cfg[6],
                         "delete_node_rate" => cfg[7],
                         "add_mutation_rate" => cfg[8],
                         "delete_mutation_rate" => cfg[9]))
    get_fitness(oneplus, ctype)
end

function ga_config(cfg::Array{Float64}, ctype::DataType)
    cfg = mod.(cfg, 1.0)
    mut = CGP.Config.index_in(CGP.mutations, cfg[1])
    cross = CGP.Config.index_in(CGP.crossovers, cfg[2])
    CGP.Config.init(Dict("mutate_method" => string("\"", mut, "\""),
                         "crossover_method" => string("\"", cross, "\""),
                         "recurrency" => cfg[3],
                         "input_mutation_rate" => cfg[4],
                         "output_mutation_rate" => cfg[5],
                         "node_mutation_rate" => cfg[6],
                         "add_node_rate" => cfg[7],
                         "delete_node_rate" => cfg[8],
                         "add_mutation_rate" => cfg[9],
                         "delete_mutation_rate" => cfg[10],
                         "ga_elitism_rate" => cfg[11],
                         "ga_crossover_rate" => cfg[12],
                         "ga_mutation_rate" => cfg[13]))
    get_fitness(GA, ctype)
end

function neat_config(cfg::Array{Float64}, ctype::DataType)
    cfg = mod.(cfg, 1.0)
    mut = CGP.Config.index_in(CGP.mutations, cfg[1])
    cross = CGP.Config.index_in(CGP.crossovers, cfg[2])
    dist = CGP.Config.index_in(CGP.distances, cfg[3])
    CGP.Config.init(Dict("mutate_method" => string("\"", mut, "\""),
                         "crossover_method" => string("\"", cross, "\""),
                         "distance_method" => string("\"", dist, "\""),
                         "recurrency" => cfg[4],
                         "input_mutation_rate" => cfg[5],
                         "output_mutation_rate" => cfg[6],
                         "node_mutation_rate" => cfg[7],
                         "add_node_rate" => cfg[8],
                         "delete_node_rate" => cfg[9],
                         "add_mutation_rate" => cfg[10],
                         "delete_mutation_rate" => cfg[11],
                         "speciation_thresh" => cfg[12],
                         "ga_crossover_rate" => cfg[13],
                         "ga_mutation_rate" => cfg[14]))
    get_fitness(cgpneat, ctype)
end

function param_sweep()
    for ctype in CGP.chromosomes
        pure_cmaes(x->oneplus_config(x, ctype), rand(9), 0.1*ones(9);
                   lambda=10, stopeval=100)
        pure_cmaes(x->ga_config(x, ctype), rand(13), 0.1*ones(13);
                   lambda=10, stopeval=100)
        pure_cmaes(x->neat_config(x, ctype), rand(14), 0.1*ones(14);
                   lambda=10, stopeval=100)
    end
end


text:
module TestInOrder
using Test

include("../src/Tree/inorder.jl")

@testset "Inorder                        " begin
    @test binary_inorder(1) == [0]
    @test binary_inorder(2) == [1, 0, 2]
    @test binary_inorder(3) == [3, 1, 4, 0, 5, 2, 6]
    @test binary_inorder(4) == [7, 3, 8, 1, 9, 4, 10, 0, 11, 5, 12, 2, 13, 6, 14]
end

end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
##  SCRIPT NAME:   utilityH.jl
##  DATE:          March 2015
#
#   This script calculates the felicity the husband enjoys for every value
#   of his private consumption and the family level public consumption.
#
#
##  ALEXANDROS THELOUDIS, UCL

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

function utilityH(cH, K)

  utility = alphaCH*(cH^(1-gammaCH))/(1-gammaCH) + (1-alphaCH)*(K^(1-gammaKH))/(1-gammaKH)

  return utility

end

max_stars_repo_name: floswald/ucl-econ-julia

max_issues_repo_name: floswald/ucl-econ-julia

max_forks_repo_name: floswald/ucl-econ-julia


text:
export deconvolution

"""
    deconvolution(measured, psf; <keyword arguments>)

Computes the deconvolution of `measured` and `psf`. Return parameter is a tuple with
two elements. The first entry is the deconvolved image. The second return parameter
is the output of the optimization of Optim.jl

Multiple keyword arguments can be specified for different loss functions,
regularizers and mappings.

# Arguments
- `loss=Poisson()`: the loss function taking a vector the same shape as measured.
- `regularizer=GR()`: A regularizer function, same form as `loss`.
- `λ=0.05`: A float indicating the total weighting of the regularizer with
    respect to the global loss function
- `background=0`: A float indicating a background intensity level.
- `mapping=Non_negative()`: Applies a mapping of the optimizer weight. Default is a
              parabola which achieves a non-negativity constraint.
- `iterations=nothing`: Specifies a number of iterations after the optimization.
    definitely should stop. By default 20 iterations will be selected by generic_invert.jl,
    if `nothing` is provided.
- `conv_dims`: A tuple indicating over which dimensions the convolution should happen.
               per default `conv_dims=1:ndims(psf)`
- `plan_fft=true`: Boolean whether plan_fft is used. Gives a slight speed improvement.
- `padding=0`: an float indicating the amount (fraction of the size in that dimension)
    of padded regions around the reconstruction. Prevents wrap around effects of the FFT.
    A array with `size(arr)=(400, 400)` with `padding=0.05` would result in reconstruction size of
    `(440, 440)`. However, if padding is >= 0.0, we only return the reconstruction cropped to the original size.
    For negative paddings, the absolute value is used, but the result maintains the padded size.
    `padding=0` disables any padding.
- `opt_package=Opt_Optim`: decides which backend for the optimizer is used.
- `opt=LBFGS()`: The chosen optimizer which must fit to `opt_package`
- `opt_options=nothing`: Can be a options file required by Optim.jl. Will overwrite iterations.
- `initial=mean(measured)`: defines a value (or array) with the initial guess. This will be pulled through the inverse mapping function
                     and extended with a mean value (if border regions are used).
- `debug_f=nothing`: A debug function which must take a single argument, the current reconstruction.

# Example
```julia-repl
julia> using DeconvOptim, TestImages, Colors, Noise;

julia> img = Float32.(testimage("resolution_test_512"));

julia> psf = Float32.(generate_psf(size(img), 30));

julia> img_b = conv(img, psf);

julia> img_n = poisson(img_b, 300);

julia> res, o = deconvolution(img_n, psf);
```
"""
function deconvolution(measured::AbstractArray{T, N}, psf;
        loss=Poisson(),
        regularizer=GR(),
        λ=T(0.05),
        background=zero(T),
        mapping=Non_negative(),
        iterations=nothing,
        conv_dims = ntuple(+, ndims(psf)),
        padding=0.00,
        opt_options=nothing,
        opt=LBFGS(linesearch=BackTracking()),
        initial=mean(measured),
        debug_f=nothing,
        opt_package=Opt_Optim) where {T, N}

    # rec0 will be an array storing the final reconstruction
    # we choose it larger than the measured array to reduce
    # wrap around artifacts of the Fourier Transform
    # we create a array size_padded which stores a new array size
    # our reconstruction array will be larger than measured
    # to prevent wrap around artifacts
    size_padded = []
    for i = 1:ndims(measured)
        # if the size of the i-th dimension is 1
        # don't do any padding because there won't be no
        # convolution in that dimension
        if  size(measured)[i] == 1
            push!(size_padded, 1)
        else
            # only pad, if padding is true
            if ~iszero(padding)
                # 2 * ensures symmetric padding
                # minimum padding is 2 (4 in total) on each side
                x = next_fast_fft_size(max(4, 2 * round(Int, size(measured)[i] * abs(padding))))
            else
                x = 0
            end
            push!(size_padded, size(measured)[i] + x)
        end
    end

    # we divide by the mean to normalize
    rescaling = mean(measured)
    measured = measured ./ rescaling
    initial = initial ./ rescaling

    # create rec0 which will be the initial guess for the reconstruction
    rec0 = similar(measured, (size_padded)...)
    fill!(rec0, one(eltype(measured)))

    # alternative rec0_center, unused at the moment
    #rec0_center = m_invf(abs.(conv(measured, psf, conv_dims)))
    #
    # take the mean as the initial guess
    # therefore has the same total energy at the initial guess as
    # measured
    csize = isa(initial,AbstractArray) ? size(initial) : size(measured)
    one_arr = similar(measured, size(measured))
    fill!(one_arr, mean(measured))
    center_set!(rec0, one_arr .* initial)
    mf, mf_inv = get_mapping(mapping)
    rec0 = mf_inv(rec0)
    # psf_n is the psf with the same size as rec0 but only in that dimensions
    # that were supported by the initial psf. Broadcasting of psf with less
    # dimensions is still supported
    # we put the small psf into the new one
    # it is important to pad the PSF instead of the OTF

    psf_new_size = Array{Int}(undef, 0)
    for i = 1:ndims(psf)
        push!(psf_new_size, size(rec0)[i])
    end

    psf_new_size = tuple(psf_new_size...)
    psf_n = similar(rec0, psf_new_size)
    fill!(psf_n, zero(eltype(rec0)))
    psf_n = center_set!(psf_n, fftshift(psf))
    psf = ifftshift(psf_n)

    # the psf should be normalized to 1
    psf ./= sum(psf)

    otf, conv_temp = plan_conv(rec0, psf, conv_dims)

    # forward model is a convolution
    # due to numerics, we need to clip at 0
    # analytically it's a convolution psf ≥ 0 and image ≥ 0
    # so it must be conv(psf, image) ≥ 0
    forward(x) = let
        if iszero(background)
            center_extract((conv_aux(conv_temp, x, otf)), size(measured))
        else
            center_extract((conv_aux(conv_temp, x, otf) .+ background), size(measured))
        end
    end
    # pass to more general optimization
    res_out, res = invert(measured, rec0, forward;
                          iterations=iterations, λ=λ,
                          regularizer=regularizer,
                          opt=opt,
                          opt_options=opt_options,
                          mapping=mapping,
                          loss=loss,
                          debug_f=debug_f, opt_package=opt_package)

    res_out .*= rescaling
    # since we do some padding we need to extract the center part
    # for negative paddings, keep the large size.
    if padding > 0.0
        res_out = center_extract(res_out, size(measured))
    end
    return res_out, res
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %% [markdown]
# See https://discourse.julialang.org/t/whats-wrong-with-my-chi-squared-goodness-of-fits-tests/64334

# %%
using Distributions, StatsPlots
P = plot(DiscreteUniform(0, 100); label="DiscreteUniform(0, 100)", xtick=0:10:100, ylim=(0, 0.013))
Q = plot(Binomial(100, 0.7); label="Binomial(100, 0.7)", legend=:topleft, xtick=0:10:100)
plot(P, Q; size=(800, 250))

# %%
using Distributions, Random, StatsPlots

# The test definitions.....
function computeDensity(data, supp)
    counts =  [count(i -> i==s,data) for s in supp]
    if length(data) > sum(counts)
        error("There are some data not in the support !")
    end
    return counts
end

"""
Modified version of `goodnessOfFitDiscrete` function in
https://discourse.julialang.org/t/whats-wrong-with-my-chi-squared-goodness-of-fits-tests/64334

* The old argument `support` is deleted.
* The new local variable `supp` is defined to be `support(f₀)`
"""
function goodnessOfFitDiscrete(data, f₀; compressedData=true, α=0.05, d=0)
    supp       = support(f₀)
    if !compressedData
        data   = computeDensity(data, supp)
    end
    K          = length(supp)
    N          = sum(data)
    p̂          = data ./ N
    df         = K - d - 1
    p0         = pdf.(f₀,supp)
    T          = N * sum((p̂[k] - p0[k])^2/p0[k] for k in 1:K)
    χDist      = Chisq(df)
    rejectedH₀ = T > quantile(χDist, 1-α)
    p_value    = 1 - cdf(χDist, T)
    return (testValue=T, threashold=quantile(χDist,1-α), rejectedH₀=rejectedH₀, p_value=p_value)
end

function repeat_tests(f₀; datasize = 10000, repetitions = 10000, α = 0.05, d = 0)
    testValue = zeros(repetitions)
    rejectedH₀ = falses(repetitions)
    data = rand(f₀, datasize)
    for rep in 1:repetitions
        rand!(f₀, data)
        out = goodnessOfFitDiscrete(data, f₀; compressedData=false, α, d)
        testValue[rep] = out.testValue
        rejectedH₀[rep] = out.rejectedH₀
    end
    α_real = sum(rejectedH₀)/repetitions
    (; α_real, testValue, rejectedH₀)
end

function plot_testValue_dist(f₀; datasize = 10000, repetitions = 10000, α = 0.05, d=0)
    α_real, testValue, rejectedH₀ = repeat_tests(f₀; datasize, repetitions, α, d)

    title = "$f₀   (real α = $α_real / nominal α = $α)"
    plot(; title, titlefontsize=10)

    xlim=(0, 2length(support(f₀)))
    a, b = xlim
    histogram!(testValue[a .≤ testValue .≤ b]; norm=true, alpha=0.3, label="testValue", xlim)

    df = length(support(f₀)) - d - 1
    plot!(Chisq(df), a, b; label="Chisq(df = $df)", ls=:dash, lw=2)
end

# %%
plot_testValue_dist(DiscreteUniform(0, 100))

# %%
plot_testValue_dist(Binomial(100, 0.7))

# %%
plot_testValue_dist(Binomial(100, 0.7); d = 40)

# %%

max_stars_repo_path: 0010/On chi-squared goodness of fits tests.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
module LeastSquares

export curve_fit, estimate_errors, estimate_covar
export lsfit, @afuns

#using Optim
using Calculus
##TODO: not clear I need Distributions
#import Distributions: quantile
import Optim: OptimizationState, MultivariateOptimizationResults, OptimizationTrace
include("levenberg_marquardt.jl")
include("curve_fit.jl")
include("lsfit.jl")

end # module

max_stars_repo_name: gabrielgellner/LeastSquares

max_issues_repo_name: gabrielgellner/LeastSquares

max_forks_repo_name: gabrielgellner/LeastSquares


text:
@testset "Quadratic" begin
    using STMO.Quadratic

    p, q, r = 4.0, 3.0, 1.0

    @test solve_quadratic(p, q, r) ≈ -q / p
    @test_throws AssertionError solve_quadratic(-3, 1)

    P = [4 1; 1 2.0]
    q = [-3.0, 1.0]

    @test solve_quadratic(P, q) ≈ - P \ q
    @test norm(P * solve_quadratic(P, q) + q) < 1e-7

    x0 = zeros(2)
    @test gradient_descent(P, q, copy(x0), ϵ=1e-12) ≈ solve_quadratic(P, q)
    tracker = PathTrack(x0)

    @test tracker isa PathTrack
    gradient_descent(P, q, copy(x0), tracker=tracker)

    @test nsteps(tracker) > 1

    @test gradient_descent(P, q, copy(x0), ϵ=1e-12, β=0.2) ≈ solve_quadratic(P, q)
end


text:
import FITSIO
import JLD
using DataStructures

import ..Log
using ..Model
import ..SDSSIO
import ..Infer
import ..SDSSIO: RunCamcolField
import ..PSF

using ..DeterministicVI
using ..DeterministicVI.ConstraintTransforms: ConstraintBatch, DEFAULT_CHUNK
using ..DeterministicVI.ElboMaximize: ElboConfig, maximize!

const PeakFlops = false
gdlog = nothing

# 20 minute threshold
const BoxMaxThreshold = convert(UInt64, 20*60*1e9)::UInt64
const OneMinute = 0x0000000df8475800 # in nanoseconds

mutable struct BoxKillSwitch <: Function
    started::UInt64
    numfin_threshold::Int64
    numfin::Atomic{Int64}
    finished::Vector{UInt64}
    lastfin::Atomic{Int64}
    killed::Bool
    messaged::Bool

    BoxKillSwitch() = new(0, ceil(Int64, nthreads()/2), Atomic{Int64}(0),
                          fill(typemax(UInt64), nthreads()),
                          Atomic{Int64}(0), false, false)
end

if PeakFlops
function (f::BoxKillSwitch)(x)
    now_ns = time_ns()
    if f.started == 0
        # benign race here
        f.started = now_ns
        return false
    end
    intv = now_ns - f.started
    if intv > BoxMaxThreshold
        # and here
        return f.killed = true
    end
    if div(intv, OneMinute) > 0
        f.started = now_ns
        n_active = 0
        n_inactive = 0
        for elbo_vars in DeterministicVI.ElboMaximize.ELBO_VARS_POOL
            n_active += elbo_vars.active_pixel_counter[]
            n_inactive += elbo_vars.inactive_pixel_counter[]
        end
        message(gdlog, "pixel visits: ($n_active,$n_inactive)")
    end
    return f.killed
end
else # !PeakFlops
function (f::BoxKillSwitch)(x)
    now_ns = time_ns()
    if f.started == 0
        # benign race here
        f.started = now_ns
        return false
    end
    intv = now_ns - f.started
    if intv > BoxMaxThreshold
        # and here
        return f.killed = true
    end
#=
    if f.numfin[] >= f.numfin_threshold
        lastfinat = f.finished[f.lastfin[]]
        if now_ns - lastfinat > floor(Int64, (lastfinat - f.started) * 0.25)
            #f.killed = true
            if !f.messaged
                f.messaged = true
                ttms = [f.finished[i] == typemax(UInt64) ? "" :
                            @sprintf("%d:%.3f ", i, (f.finished[i] - f.started) / 1e6)
                        for i in 1:length(f.finished)]
                Log.message("imbalance threshold: $(ttms...)ms")
            end
        end
=#
    return f.killed
end
end

function ks_reset_finished(ks::BoxKillSwitch)
    ks.numfin[] = 0
    fill!(ks.finished, typemax(UInt64))
    ks.lastfin[] = 0
    ks.messaged = false
end

function union_find!(i, components_tree)
    root = i
    while components_tree[i] != i
        i = components_tree[i]
    end

    while root != i
        root2 = components_tree[root]
        components_tree[root] = i
        root = root2
    end
    return i
end

"""
Computes connected components given a conflict graph, and the range of sources
to process. Writes to the components argument.

- source_start - start source
- end_source - end source (inclusively processed)
- sources - the actual source ids to find connected components for
    (indexed by source_start and source_end)
- neighbor_map - conflict graph of sources
- components_tree - the components array to do union find on.
- components - Dict{Int64, Vector{Int64}} dictionary from cc_id to
              target_source index in that component
- src_to_target_index - a mapping from light source id -> index
            within the target_sources array. This is required
            since we need to write the index within target sources
            to the output, rather than the source id itself.
Returns:
- Nothing
"""
function compute_connected_components(source_start, source_end, sources, neighbor_map,
                                      components_tree, components, src_to_target_index)
    # Construct reverse map from source_id -> index in source. Only do this for
    # the specified portion we are computing connected components.
    source_to_index = Dict{Int64, Int64}()
    index_to_source = Dict{Int64, Int64}()
    for i = source_start:source_end
        source_to_index[sources[i]] = i-source_start+1
        index_to_source[i-source_start+1] = sources[i]
    end

    for i = source_start:source_end
        components_tree[i-source_start+1] = i-source_start+1
    end

    # Run union find algorithm to find connected components.
    for i = source_start:source_end
        target = union_find!(i-source_start+1, components_tree)
        for neighbor in neighbor_map[index_to_source[i-source_start+1]]
            if haskey(source_to_index, neighbor)
                neighbor_idx = source_to_index[neighbor]
                component_source = union_find!(neighbor_idx, components_tree)
                components_tree[component_source] = target
            end

        end
    end

    # Write to components dictionary.
    for i = source_start:source_end
        index = union_find!(i-source_start+1, components_tree)
        if !haskey(components, index)
            components[index] = Vector{Int64}()
        end
        push!(components[index], src_to_target_index[index_to_source[i-source_start+1]])
    end
end

"""
Partitions sources via the cyclades algorithm.

- n_threads - number of threads to which to distribute sources
- target_sources - array of target sources. Elements should match keys of neighbor_map.
- neighbor_map - graph of connections of sources

Returns:
- An array of vectors representing the workload of each thread ([thread][batch][sources(indices)])
"""
function partition_cyclades(n_threads, target_sources, neighbor_map; batch_size=60)
    #Log.info("Starting Cyclades partitioning...")
    tic()

    n_sources = length(target_sources)
    n_total_batches = convert(Int64, ceil(n_sources / batch_size))

    # Construct src_to_indx map
    src_to_indx = Dict{Int64, Int64}()
    for i=1:n_sources
        src_to_indx[target_sources[i]] = i
    end

    # The final workload distribution.
    thread_sources_assignment = Vector{Vector{Vector{Int64}}}(n_threads)
    for thread = 1:n_threads
        thread_sources_assignment[thread] = Vector{Vector{Int64}}(n_total_batches)
        for batch = 1:n_total_batches
            thread_sources_assignment[thread][batch] = Vector{Int64}()
        end
    end

    # First shuffle the sources. Note Cyclades is serially equivalent
    # to this permutation of sources.
    sources = [x for x in keys(neighbor_map)]
    shuffle!(sources)

    # Process batch_size number of sources at a time.
    # TODO max - parallelize everything below (particularly CC computation)

    # The components tree is for union-find.
    components_tree = [Vector{Int64}(batch_size) for i=1:n_threads]

    # We have n_total_batches components, where each component is a dictionary
    # from the component id, to a list of sources in that component
    components = [Dict{Int64, Vector{Int64}}() for i=1:n_total_batches]
    sources_in_components = 0
    for (cur_batch, source_index) in enumerate(collect(1:batch_size:n_sources))
        source_start = source_index
        source_end = min(source_index+batch_size-1, n_sources)
        @assert source_end - source_start + 1 <= batch_size

        # TODO max - parallelize this
        compute_connected_components(source_start, source_end, sources, neighbor_map,
                                     components_tree[1],
                                     components[cur_batch],
                                     src_to_indx)
    end

    assigned_sources = 0
    # Load balance the connected components within each batch into thread_sources_assignment.
    for (cur_batch, cur_batch_component) in enumerate(components)
        # Priority queue for load balancing.
        pqueue = PriorityQueue(Base.Order.Forward, (i, 0) for i=1:n_threads)

        # Assign non-conflicting group of sources to different threads
        for (component_group_id, sources_of_component) in cur_batch_component
            least_loaded_thread = peek(pqueue)[1]
            for source_to_assign in sources_of_component
                push!(thread_sources_assignment[least_loaded_thread][cur_batch], source_to_assign)
                assigned_sources += 1
            end
            pqueue[least_loaded_thread] += length(sources_of_component)
        end
    end

    #Log.info("Cyclades - Assigned sources: $(assigned_sources) vs correct number of sources: $(n_sources)")
    @assert assigned_sources == n_sources
    #Log.info("Cyclades - Number of batches: $(n_total_batches)")
    #Log.info("Finished Cyclades partitioning.  Elapsed time: $(toq()) seconds")

    for cur_batch = 1:length(collect(1:batch_size:n_sources))
        load_balance_for_batch = [length(thread_sources_assignment[t][cur_batch]) for t=1:n_threads]
        #Log.info("Load balance for batch $(cur_batch) - $(load_balance_for_batch)")
    end

    thread_sources_assignment
end

"""
Partitions sources via the cyclades algorithm.

- target_sources - array of target sources. Elements should match keys of neighbor_map.
- neighbor_map - graph of connections of sources

Returns:
- An array of vectors representing the workload of each thread ([thread][batch][sources(indices)])
"""
function partition_cyclades_dynamic(target_sources, neighbor_map; batch_size=60)
    #Log.info("Starting Cyclades partitioning...")
    tic()

    n_sources = length(target_sources)
    n_total_batches = convert(Int64, ceil(n_sources / batch_size))

    # Construct src_to_indx map
    src_to_indx = Dict{Int64, Int64}()
    for i=1:n_sources
        src_to_indx[target_sources[i]] = i
    end

    # The final workload distribution.
    thread_sources_assignment = Vector{Vector{Vector{Int64}}}(n_total_batches)
    for batch = 1:n_total_batches
        thread_sources_assignment[batch] = Vector{Vector{Int64}}()
    end

    # First shuffle the sources. Note Cyclades is serially equivalent
    # to this permutation of sources.
    sources = [x for x in keys(neighbor_map)]
    shuffle!(sources)

    # Process batch_size number of sources at a time.
    # TODO max - parallelize everything below (particularly CC computation)

    # The components tree is for union-find.
    components_tree = Vector{Int64}(batch_size)

    # We have n_total_batches components, where each component is a dictionary
    # from the component id, to a list of sources in that component
    components = [Dict{Int64, Vector{Int64}}() for i=1:n_total_batches]
    sources_in_components = 0
    for (cur_batch, source_index) in enumerate(collect(1:batch_size:n_sources))
        source_start = source_index
        source_end = min(source_index+batch_size-1, n_sources)
        @assert source_end - source_start + 1 <= batch_size

        # TODO max - parallelize this
        compute_connected_components(source_start, source_end, sources, neighbor_map,
                                     components_tree,
                                     components[cur_batch],
                                     src_to_indx)
    end

    assigned_sources = 0
    # Add sources to sources assignment
    for (cur_batch, cur_batch_component) in enumerate(components)
        for (component_group_id, sources_of_component) in cur_batch_component
            push!(thread_sources_assignment[cur_batch], copy(sources_of_component))
        assigned_sources += length(sources_of_component)
        end
    end

    #Log.info("Cyclades - Assigned sources: $(assigned_sources) vs correct number of sources: $(n_sources)")
    @assert assigned_sources == n_sources
    #Log.info("Cyclades - Number of batches: $(n_total_batches)")
    #Log.info("Finished Cyclades partitioning.  Elapsed time: $(toq()) seconds")

    #Log.info("Assigned sources: $(thread_sources_assignment)")

    thread_sources_assignment
end

"""
Partitions sources across threads equally.

- n_threads - number of threads to which to distribute processing the sources.
- n_sources - the number of total sources to process.

Returns:
- An array of vectors representing the workload of each thread ([thread][batch][sources]).
  Note for partition equally, there is only 1 batch.

"""
function partition_equally(n_threads, n_sources)
    #Log.info("Starting basic source partitioning...")
    tic()
    n_sources_per_thread = floor(Int64, n_sources / n_threads)
    thread_sources_assignment = Vector{Vector{Vector{Int64}}}(n_threads)
    n_sources_assigned = 0
    for thread = 1:n_threads
        start_source = (thread-1) * n_sources_per_thread
        end_source = thread * n_sources_per_thread
        if thread == n_threads
            end_source = n_sources
        end
        # Only 1 batch for partitioning equally
        thread_sources_assignment[thread] = Vector{Vector{Int64}}(1)
        thread_sources_assignment[thread][1] = Vector{Int64}(end_source-start_source)
        for source_assignment = start_source:end_source-1
            thread_sources_assignment[thread][1][source_assignment-start_source+1] = source_assignment+1
            n_sources_assigned += 1
        end
    end
    @assert n_sources_assigned == n_sources
    #Log.info("Finished basic source partitioning. Elapsed time: $(toq()) seconds")
    thread_sources_assignment
end

"""
Return a Cyclades partitioning or an equal partitioning of the target
sources.
"""
function partition_box(npartitions::Int, target_sources::Vector{Int},
                       neighbor_map::Vector{Vector{Int}}, ea_vec;
                       cyclades_partition=true,
                       batch_size=7000)
    if cyclades_partition
        cyclades_neighbor_map = Dict{Int64,Vector{Int64}}()
        for (index, neighbors) in enumerate(neighbor_map)
            cyclades_neighbor_map[target_sources[index]] = neighbors
        end
        #return partition_cyclades(npartitions, target_sources,
        #                          cyclades_neighbor_map,
        #                          batch_size=batch_size)
        #return partition_cyclades_dynamic(target_sources,
        #                                  cyclades_neighbor_map,
        #                                  batch_size=batch_size)
	return partition_cyclades_dynamic_auto_batchsize(target_sources, cyclades_neighbor_map, ea_vec)
    else
        return partition_equally(npartitions, length(target_sources))
    end
end

function estimate_time(patches)
    sum(sum(p.active_pixel_bitmap) for p in patches)
end

function load_balance_across_threads(n_threads, times)
    ts = [0 for i=1:n_threads]
    for t in times
        minimum,index_min = findmin(ts)
        ts[index_min] += t
    end
    return ts
end

"""
Partitions sources via the cyclades algorithm. Finds the batch size which returns most balanced CC distribution.
- target_sources - array of target sources. Elements should match keys of neighbor_map.
- neighbor_map - graph of connections of sources
"""
function partition_cyclades_dynamic_auto_batchsize(target_sources, neighbor_map, ea_vec)
    n_threads = nthreads()

    # Sample batch sizes at intervals
    n_to_sample = 100
    stepsize = max(1, trunc(Int, length(target_sources) / n_to_sample))

    # Best load imbalance of the batch sizes to test
    best_score = Inf
    best_result = Inf
    best_batch_size = -Inf

    for batch_size_to_use = 1 : stepsize : length(target_sources)+1
        ccs = partition_cyclades_dynamic(target_sources, neighbor_map, batch_size=batch_size_to_use)
        score = 0
        for batch in ccs
            # Find average load imbalance within the batch as a percentage
	    times = [sum([estimate_time(ea_vec[source_index].patches) for source_index in component]) for component in batch]
            times = load_balance_across_threads(n_threads, times)
            estimated_imbalance = mean(maximum(times) - times)
	    score += estimated_imbalance
        end
        if score <= best_score
            best_result = ccs
            best_score = score
            best_batch_size = batch_size_to_use
        end
    end
    for batch in best_result
        sizes = [length(x) for x in batch]
    end
    best_result
end

"""
Pre-allocate elbo args variational parameters; configurations need to
be persisted across calls to `maximize!()` so that location constraints
don't shift from their initial position.
"""
function setup_vecs(n_sources::Int, target_sources::Vector{Int},
                    catalog::Vector{CatalogEntry})
    ts_vp = Dict{Int64, Array{Float64}}()
    for ts in target_sources
        cat = catalog[ts]
        ts_vp[ts] = generic_init_source(cat.pos)
    end

    ea_vec = Vector{ElboArgs}(n_sources)
    vp_vec = Vector{VariationalParams{Float64}}(n_sources)
    cfg_vec = Vector{ElboConfig{DEFAULT_CHUNK,Float64}}(n_sources)

    return ea_vec, vp_vec, cfg_vec, ts_vp
end

"""
Like one_node_infer, uses multiple threads on one node to fit the Celeste
model over numerous iterations.

catalog - the catalog of light sources
target_sources - light sources to optimize
neighbor_map - ligh_source index -> neighbor light_source id

cyclades_partition - use the cyclades algorithm to partition into non conflicting batches for updates.
batch_size - size of a single batch of sources for updates
within_batch_shuffling - whether or not to process sources within a batch randomly
joint_inference_terminate - whether to terminate once sources seem to be stable
n_iters - number of iterations to optimize. 1 iteration optimizes a full pass over target
          sources if optimize_fixed_iters=true.

Returns:

- Vector of OptimizedSource results
"""
function one_node_joint_infer(config::Config, catalog, target_sources, neighbor_map, images;
                              cyclades_partition::Bool=true,
                              batch_size::Int=7000,
                              within_batch_shuffling::Bool=true,
                              n_iters::Int=3,
                              timing=InferTiming())
    # Seed random number generator to ensure the same results per run.
    srand(42)

    n_threads = nthreads()

    # Partition the sources
    n_sources = length(target_sources)
    Log.info("Optimizing $(n_sources) sources")

    #Log.info(batched_connected_components)

    Log.info("Done assigning sources to threads for processing")

    ea_vec, vp_vec, cfg_vec, target_source_variational_params =
            setup_vecs(n_sources, target_sources, catalog)

    ks = BoxKillSwitch()

    # Initialize elboargs in parallel
    tic()
    thread_initialize_sources_assignment::Vector{Vector{Vector{Int64}}} = partition_equally(n_threads, n_sources)

    initialize_elboargs_sources!(config, ea_vec, vp_vec, cfg_vec, thread_initialize_sources_assignment,
                                 catalog, target_sources, neighbor_map, images,
                                 target_source_variational_params;
                                 termination_callback=ks)
    timing.init_elbo = toq()

    #thread_sources_assignment = partition_box(n_threads, target_sources,
    #                                  neighbor_map;
    #                                  cyclades_partition=cyclades_partition,
    #                                  batch_size=batch_size)
    batched_connected_components = partition_box(n_threads, target_sources,
                                                 neighbor_map, ea_vec;
                                                 cyclades_partition=cyclades_partition,
                                                 batch_size=batch_size)

    # Process sources in parallel
    tic()

    #process_sources!(images, ea_vec, vp_vec, cfg_vec,
    #                 thread_sources_assignment,
    #                 n_iters, within_batch_shuffling)
    process_sources_dynamic!(images, ea_vec, vp_vec, cfg_vec,
                             batched_connected_components,
                             n_iters, within_batch_shuffling;
                             kill_switch=ks)

    timing.opt_srcs = toq()
    timing.num_srcs = n_sources

    # Return add results to vector
    results = OptimizedSource[]

    if ks.killed
        Log.message("optimization terminated due to timeout")
    else
        for i = 1:n_sources
            entry = catalog[target_sources[i]]
            result = OptimizedSource(entry.pos[1],
                                     entry.pos[2],
                                     vp_vec[i][1])
            push!(results, result)
        end
    end

    show_pixels_processed()

    results
end

# legacy wrapper
function one_node_joint_infer(catalog, target_sources, neighbor_map, images;
                              cyclades_partition::Bool=true,
                              batch_size::Int=7000,
                              within_batch_shuffling::Bool=true,
                              n_iters::Int=3,
                              timing=InferTiming())
    one_node_joint_infer(
        Config(),
        catalog,
        target_sources,
        neighbor_map,
        images,
        cyclades_partition=cyclades_partition,
        batch_size=batch_size,
        within_batch_shuffling=within_batch_shuffling,
        n_iters=n_iters,
        timing=timing
    )
end

function initialize_elboargs_sources!(config::Config, ea_vec, vp_vec, cfg_vec,
                                      thread_initialize_sources_assignment,
                                      catalog, target_sources, neighbor_map, images,
                                      target_source_variational_params;
                                      termination_callback=nothing)
    Threads.@threads for i in 1:nthreads()
        try
            for batch in 1:length(thread_initialize_sources_assignment[i])
                for source_index in thread_initialize_sources_assignment[i][batch]
                    init_elboargs(config, source_index, catalog, target_sources,
                                  neighbor_map, images, ea_vec, vp_vec, cfg_vec,
                                  target_source_variational_params;
                                  termination_callback=termination_callback)
                end
            end
        catch ex
            Log.exception(ex)
            rethrow()
        end
    end
end

"""
Initialize elbo args for the specified target source.
"""
function init_elboargs(config::Config,
                       ts::Int,
                       catalog::Vector{CatalogEntry},
                       target_sources::Vector{Int},
                       neighbor_map::Vector{Vector{Int}},
                       images::Vector{Image},
                       ea_vec::Vector{ElboArgs},
                       vp_vec::Vector{VariationalParams{Float64}},
                       cfg_vec::Vector{ElboConfig{DEFAULT_CHUNK,Float64}},
                       ts_vp::Dict{Int64,Array{Float64}};
                       termination_callback=nothing)
    try
        entry_id = target_sources[ts]
        entry = catalog[entry_id]
        neighbor_ids = neighbor_map[ts]
        neighbors = catalog[neighbor_ids]
        cat_local = vcat([entry], neighbors)
        ids_local = vcat([entry_id], neighbor_ids)

        patches = Infer.get_sky_patches(images, cat_local)
        Infer.load_active_pixels!(config, images, patches)
        # Load vp with shared target source params, and also vp
        # that doesn't share target source params
        vp = Vector{Float64}[haskey(ts_vp, x) ?
                             ts_vp[x] :
                             catalog_init_source(catalog[x])
                             for x in ids_local]
        ea = ElboArgs(images, patches, [1])

        ea_vec[ts] = ea
        vp_vec[ts] = vp
        cfg_vec[ts] = ElboConfig(ea, vp;
                termination_callback=termination_callback)
    catch exc
        if is_production_run || nthreads() > 1
            Log.exception(exc)
        else
            rethrow()
        end
    end
end

function process_sources!(images::Vector{Model.Image},
                          ea_vec::Vector{ElboArgs},
                          vp_vec::Vector{VariationalParams{Float64}},
                          cfg_vec::Vector{ElboConfig{DEFAULT_CHUNK,Float64}},
                          thread_sources_assignment::Vector{Vector{Vector{Int64}}},
                          n_iters::Int,
                          within_batch_shuffling::Bool)
    n_threads::Int = nthreads()
    n_batches::Int = length(thread_sources_assignment[1])
    for iter in 1:n_iters
        # Process every batch of every iteration. We do the batches on the outside
        # Since there is an implicit barrier after the inner threaded for loop below.
        # We want this barrier because there may be conflict _between_ Cyclades batches.
        for batch in 1:n_batches
            process_sources_elapsed_times::Vector{Float64} = Vector{Float64}(n_threads)
            # Process every batch of every iteration with n_threads
            Threads.@threads for i in 1:n_threads
                try
                    tic()
                    process_sources_kernel!(ea_vec, vp_vec, cfg_vec,
                                            thread_sources_assignment[i::Int][batch::Int],
                                            within_batch_shuffling)
                    process_sources_elapsed_times[i::Int] = toq()
                catch exc
                    Log.exception(exc)
                    rethrow()
                end
            end
            Log.info("Batch $(batch) - $(process_sources_elapsed_times)")
        end
    end
end

function process_sources_dynamic!(images::Vector{Model.Image},
                                  ea_vec::Vector{ElboArgs},
                                  vp_vec::Vector{VariationalParams{Float64}},
                                  cfg_vec::Vector{ElboConfig{DEFAULT_CHUNK,Float64}},
                                  thread_sources_assignment::Vector{Vector{Vector{Int64}}},
                                  n_iters::Int,
                                  within_batch_shuffling::Bool;
                                  kill_switch=nothing)
    Log.info("Processing with dynamic connected components load balancing")

    n_threads::Int = nthreads()
    n_batches::Int = length(thread_sources_assignment)

    l = SpinLock()

    total_idle_time = 0
    total_sum_of_thread_times = 0

    for iter in 1:n_iters

        # Process every batch of every iteration. We do the batches on the outside
        # Since there is an implicit barrier after the inner threaded for loop below.
        # We want this barrier because there may be conflict _between_ Cyclades batches.
        for batch in 1:n_batches

            connected_components_index = 1
            process_sources_elapsed_times::Vector{Float64} = Vector{Float64}(n_threads)

            # Process every batch of every iteration with n_threads
            Threads.@threads for i in 1:n_threads
                try
                    tic()
                    while true
                        cc_index = -1

                        #ccall(:jl_,Void,(Any,), "this is thread number $(Threads.threadid()) $(cc_index)")

                        lock(l)
                        cc_index = connected_components_index
                        connected_components_index += 1
                        unlock(l)

                        if cc_index > length(thread_sources_assignment[batch::Int])
                            break
                        end

                        #ccall(:jl_,Void,(Any,), "this is thread number $(Threads.threadid()) $(cc_index) popped")

                        process_sources_kernel!(ea_vec, vp_vec, cfg_vec,
                                                thread_sources_assignment[batch::Int][cc_index],
                                                within_batch_shuffling)
                        if kill_switch != nothing && kill_switch.killed
                            break
                        end
                    end
                    process_sources_elapsed_times[i::Int] = toq()
                catch exc
                    Log.exception(exc)
                    rethrow()
                end
                if kill_switch != nothing
                    tid = threadid()
                    kill_switch.finished[tid] = time_ns()
                    kill_switch.lastfin[] = tid
                    atomic_add!(kill_switch.numfin, 1)
                end
            end
            Log.info("Batch $(batch) - $(process_sources_elapsed_times)")
	    avg_thread_idle_time = mean(maximum(process_sources_elapsed_times) - process_sources_elapsed_times)
	    maximum_thread_time = maximum(process_sources_elapsed_times)
            idle_percent = 100.0 * (avg_thread_idle_time / maximum_thread_time)
            Log.info("Batch $(batch) avg threads idle: $(round(Int, idle_percent))% ($(avg_thread_idle_time) / $(maximum_thread_time))")
	    total_idle_time += sum(maximum(process_sources_elapsed_times) - process_sources_elapsed_times)
            total_sum_of_thread_times += sum(process_sources_elapsed_times)
            if kill_switch != nothing
                if kill_switch.killed
                    break
                else
                    ks_reset_finished(kill_switch)
                end
            end
        end
        if kill_switch != nothing && kill_switch.killed
            break
        end
    end
    Log.info("Total idle time: $(round(Int, total_idle_time)), Total sum of threads times: $(round(Int, total_sum_of_thread_times))")
end

# Process partition of sources. Multiple threads call this function in parallel.
function process_sources_kernel!(ea_vec::Vector{ElboArgs},
                                 vp_vec::Vector{VariationalParams{Float64}},
                                 cfg_vec::Vector{ElboConfig{DEFAULT_CHUNK,Float64}},
                                 source_assignment::Vector{Int64},
                                 within_batch_shuffling::Bool)
    try
        # Shuffle the source assignments within each batch of each process.
        # This is disabled by default because it ruins the deterministic outcome
        # required by the test cases.
        if within_batch_shuffling
            shuffle!(source_assignment)
        end

	#tic()
        for i in source_assignment
            maximize!(ea_vec[i], vp_vec[i], cfg_vec[i])
            if cfg_vec[i].optim_options.callback.killed
                break
            end
        end
	#ccall(:jl_,Void,(Any,), "this is thread number $(Threads.threadid()) took $(toq()) seconds to process batch with $(length(source_assignment)) sources")
    catch ex
        if is_production_run || nthreads() > 1
            Log.exception(ex)
        else
            rethrow(ex)
        end
    end
end

function get_pixels_processed()
    n_active = 0
    n_inactive = 0
    for elbo_vars in DeterministicVI.ElboMaximize.ELBO_VARS_POOL
        n_active += elbo_vars.active_pixel_counter[]
        n_inactive += elbo_vars.inactive_pixel_counter[]
        elbo_vars.active_pixel_counter[] = 0
        elbo_vars.inactive_pixel_counter[] = 0
    end
    return (n_active, n_inactive)
end

# Count and display the number of active and inactive pixels processed.
function show_pixels_processed()
    n_active, n_inactive = get_pixels_processed()
    Log.message("$(Time(now())): (active,inactive) pixels processed: \($n_active,$n_inactive\)")
end


text:
module XGPaint

using Interpolations
using QuadGK
using Roots
using Cosmology
using Unitful, UnitfulAstro
using Parameters
using Random
using Healpix
using PhysicalConstants.CODATA2018
using StatsBase
import ThreadsX
import Distributions

include("./pixell.jl")
include("./model.jl")
include("./util.jl")
include("./cib.jl")
include("./radio.jl")

export get_cosmology, read_halo_catalog_hdf5
export Radio_Sehgal2009, CIB_Planck2013
export paint!, generate_sources
export  paint_with_mask!, estimate_fluxes!, make_multifrequency_dataset
export fixsen_model, BlackBody, estimate_normalization
export read_from_hdf5, read_fluxes_hdf5, write_to_hdf5
end # module


text:
function SeisSmooth1(A;L=7)
	a = isodd(L) ? 0 : 1
	kernel = ones(eltype(A),L)/L;
	A = conv(A,kernel)
	A = A[floor(Int,L/2) + 1:end-floor(Int,L/2) + a]
	return A
end


text:
# Tests for Cost DDF DEA Models
@testset "CostDDFDEAModel" begin

    # Test using Book data
    X = [2 2; 1 4; 4 1; 4 3; 5 5; 6 1; 2 5; 1.6 8];
    Y = [1; 1; 1; 1; 1; 1; 1; 1];
    W = [1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 1];

    # With Gx = :Ones
    deacostddfones = deacostddf(X, Y, W, Gx = :Ones)

    @test typeof(deacostddfones) == CostDDFDEAModel

    @test nobs(deacostddfones) == 8
    @test ninputs(deacostddfones) == 2
    @test noutputs(deacostddfones) == 1
    @test ismonetary(deacostddfones) == false

    @test efficiency(deacostddfones) == efficiency(deacostddfones, :Economic)
    @test efficiency(deacostddfones, :Economic)   ≈ [0.0; 0.5; 0.5; 1.5; 3.0; 1.5; 1.5; 2.8]
    @test efficiency(deacostddfones, :Technical)  ≈ [0.0; 0.0; 0.0; 4/3; 3; 0.0; 1.0; 0.6]
    @test efficiency(deacostddfones, :Allocative) ≈ [0.0; 0.5; 0.5; 1/6; 0.0; 1.5; 0.5; 2.2]
    @test normfactor(deacostddfones) ≈ [2.0; 2.0; 2.0; 2.0; 2.0; 2.0; 2.0; 2.0]

    deacostddfonesmon = deacostddf(X, Y, W, Gx = :Ones, monetary = true)

    @test ismonetary(deacostddfonesmon) == true

    @test efficiency(deacostddfonesmon) ≈ efficiency(deacostddfones) .* normfactor(deacostddfones)
    @test efficiency(deacostddfonesmon, :Technical) ≈ efficiency(deacostddfones, :Technical) .* normfactor(deacostddfones)
    @test efficiency(deacostddfonesmon, :Allocative) ≈ efficiency(deacostddfones, :Allocative) .* normfactor(deacostddfones)

    # With Gx = :Ones and CRS (Same result in this example)
    deacostddfonescrs = deacostddf(X, Y, W, Gx = :Ones, rts = :CRS)

    @test efficiency(deacostddfonescrs, :Economic)   ≈ [0.0; 0.5; 0.5; 1.5; 3.0; 1.5; 1.5; 2.8]
    @test efficiency(deacostddfonescrs, :Technical)  ≈ [0.0; 0.0; 0.0; 4/3; 3; 0.0; 1.0; 0.6]
    @test efficiency(deacostddfonescrs, :Allocative) ≈ [0.0; 0.5; 0.5; 1/6; 0.0; 1.5; 0.5; 2.2]

    # With Gx = :Observed
    deacostddfobs = deacostddf(X, Y, W, Gx = :Observed)

    @test efficiency(deacostddfobs, :Economic)   ≈ [0.0; 0.2; 0.2; 0.428571; 0.6; 0.428571; 0.428571; 0.583333] atol = 1e-5
    @test efficiency(deacostddfobs, :Technical)  ≈ [0.0; 0.0; 0.0; 0.4; 0.6; 0.0; 1/3; 6/16] atol = 1e-5
    @test efficiency(deacostddfobs, :Allocative) ≈ [0.0; 0.2; 0.2; 0.028571; 0.0; 0.428571; 0.095238; 0.208333] atol = 1e-5

    # With Gx = :Mean
    deacostddfmean = deacostddf(X, Y, W, Gx = :Mean)

    @test efficiency(deacostddfmean, :Economic)   ≈ [0.0; 0.146520; 0.146520; 0.439560; 0.879121; 0.439560; 0.439560; 0.820513] atol = 1e-5
    @test efficiency(deacostddfmean, :Technical)  ≈ [0.0; 0.0; 0.0; 0.382775; 0.861244; 0.0; 0.299252; 0.1875] atol = 1e-5
    @test efficiency(deacostddfmean, :Allocative) ≈ [0.0; 0.146520; 0.146520; 0.056785; 0.017877; 0.439560; 0.140309; 0.633013] atol = 1e-5

    # With Gx = :Monetary
    deacostddfmon = deacostddf(X, Y, W, Gx = :Monetary)

    @test efficiency(deacostddfmon, :Economic)   ≈ [0.0; 1.0; 1.0; 3.0; 6.0; 3.0; 3.0; 5.6]
    @test efficiency(deacostddfmon, :Technical)  ≈ [0.0; 0.0; 0.0; 16/6; 6.0; 0.0; 2.0; 1.2]
    @test efficiency(deacostddfmon, :Allocative) ≈ [0.0; 1.0; 1.0; 1/3; 0.0; 3.0; 1.0; 4.4]

    # With custom direction
    deacostddfcustom = deacostddf(X, Y, W, Gx = ones(size(X)))
    @test efficiency(deacostddfcustom) ≈ efficiency(deacostddfones)

    # Print
    show(IOBuffer(), deacostddfones)

    # Test errors
    @test_throws DimensionMismatch deacostddf([1; 2 ; 3], [4 ; 5], [1; 1; 1] , Gx = :Ones) #  Different number of observations
    @test_throws DimensionMismatch deacostddf([1; 2; 3], [4; 5; 6], [1; 2; 3; 4], Gx = :Ones) # Different number of observation in prices
    @test_throws DimensionMismatch deacostddf([1 1; 2 2; 3 3 ], [4; 5; 6], [1 1 1; 2 2 2; 3 3 3], Gx = :Ones) # Different number of input prices and inputs
    @test_throws DimensionMismatch deacostddf([1; 2; 3] , [4; 5; 6], [1; 1; 1], Gx = [1 1; 2 2; 3 3]) # Different size of inputs direction
    @test_throws ArgumentError deacostddf([1; 2; 3], [4; 5; 6], [1; 2; 3], Gx = :Ones, rts = :Error) # Invalid returns to scale
    @test_throws ArgumentError deacostddf([1; 2; 3], [1; 2; 3], [1; 1; 1], Gx = :Error) # Invalid inuts direction

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 7233fedb-3123-4040-97d9-2cd3cd285a67
begin
	using Pkg
	Pkg.activate(mktempdir())
	Pkg.Registry.update()
	Pkg.add("Yao")
	Pkg.add("YaoPlots")
end

# ╔═╡ 06657244-3d4d-11eb-3a7c-f79844513992
using Yao, YaoPlots

# ╔═╡ e0cfd952-3d49-11eb-195b-b34f3836e322
md"# Deutsch Algorithm"

# ╔═╡ 06c7fa22-3d4a-11eb-2f15-b7ef46b75be3
md"Deutsch's Algorithm was the first algorithm to prove that quantum computers can perform some tasks faster than classical ones, although the speedup in this case is a trivial one.

The problem looks somewhat like this :-

	Consider 4 functions 𝑓₀, 𝑓₁, 𝑓₂ and 𝑓₃. Each of them can take 0 or 1 as input.

	𝑓₀(0) = 0, 𝑓₀(1) = 0 -> Constant function
	𝑓₁(0) = 0, 𝑓₁(1) = 1 -> Balanced function
	𝑓₂(0) = 1, 𝑓₂(1) = 0 -> Balanced function
	𝑓₃(0) = 1, 𝑓₃(1) = 1 -> Constant function

	Given one of these four functions at random, how many evaluations must we make to confirm whether the given function is balanced or constant.

Classically, it takes at least two evaluations to confirm whether the given function is balanced or constant."

# ╔═╡ 063623ea-3d4a-11eb-3471-ef67e8f31249
md"Using quantum computers, we can confirm whether the given function is balanced or not using one evaluation.

First we construct gates that correspond to the 4 functions. Here, let's consider a circuit, which takes the qubits `` |x〉 `` and `` |y〉, `` and returns the qubits `` |x〉 `` and ``|y⊕𝑓_i(x)〉`` respectively, where `` i `` can be any random number between `` 0 `` to `` 3 ``. Let's call this circuit `` F_i `` ."

# ╔═╡ f95221b4-4452-11eb-3cf6-d330a5c5aa11
YaoPlots.plot(put(2, 1:2 => label(chain(2), "𝐹ᵢ")))

# ╔═╡ 72324dca-4453-11eb-1c26-2b5c1031570b
md"Below is my implementation of 𝐹ᵢ"

# ╔═╡ f1f2f676-3d4c-11eb-307b-f7a7408372e0
begin
	𝐹₀ = chain(2)
	plot(𝐹₀)
end

# ╔═╡ 6514718e-3d4d-11eb-088b-9794ec4e7bc4
begin
	𝐹₁ = chain(2, control(1, 2=>X))
	plot(𝐹₁)
end

# ╔═╡ 7d3fd99c-3d4d-11eb-0100-67a39154fef1
begin
	𝐹₂ = chain(2, control(1, 2=>X), put(2=>X))
	plot(𝐹₂)
end

# ╔═╡ 92024f36-3d4d-11eb-3755-1f4538ed6003
begin
	𝐹₃ = chain(2, put(2=>X))
	plot(𝐹₃)
end

# ╔═╡ 2ef8834e-3d50-11eb-233c-fd44548a80b9
circuit = [𝐹₀ 𝐹₁ 𝐹₂ 𝐹₃]

# ╔═╡ 505429aa-3d50-11eb-333f-a75bc433b9ee
r = rand(0:3)

# ╔═╡ 938c4318-3d50-11eb-3060-31bc11556300
md" The new question now is :-

	Given one of the four circuits at random, how many evaluations would it take to find whether the underlying function is balanced or constant?"

# ╔═╡ 559aa08a-3d56-11eb-2103-dd4114a41999
md"If, in the above circuits we pass `` |0〉 `` or `` |1〉 ``, it'll be the same as classical computers -  it'll take 2 evaluations.
To achieve this, we'll pass a superposition of `` |0〉 `` and `` |1〉 ``, to the `` F_i `` circuit, and then we'll pass an H gate to the first qubit and measure it."

# ╔═╡ c1ca5030-4453-11eb-1814-dd9862cb5aa1
md"The final circuit for Deutsch Algorithm where the input is `` |01〉 `` looks somewhat like this"

# ╔═╡ f4e5006e-4453-11eb-3a08-a35438601cf1
plot(chain(2, repeat(H,1:2), put(1:2 => label(chain(2), "𝐹ᵢ")), put(1=>H), Measure(2, locs=1:1)))

# ╔═╡ 5dec54cc-4454-11eb-33fe-1b90858ea78a
md"Below is the final implementation, using the `` F_i `` we constructed before."

# ╔═╡ 5d52c8a8-3d50-11eb-1576-cb5f3485262b
begin
	input = ArrayReg(bit"10") #Remember, the circuit takes the qubits in reverse order
	result = input |> chain(2, repeat(H,1:2), put(1:2=>circuit[r+1]), put(1=>H)) |> r->measure(r,1)
end

# ╔═╡ a04dbca6-3d50-11eb-24c2-5dbef3190e76
result[1] == bit"1" ? "Balanced" : "Constant"

# ╔═╡ 8f5ac082-3d58-11eb-3952-8b8a74dc4e26
md"As you can see from the value of r, the result is correct."

# ╔═╡ f3568154-3ed2-11eb-125d-33a32b2030a7
md"# Deutsch–Jozsa Algorithm"

# ╔═╡ f5b13e08-3ed2-11eb-0790-d36bf95fe54f
md"The Deutsch-Jozsa Algorithm is the _general_ version of the Deutsch Algorithm. The problem is : --

	The functions are now of n variables. The inputs to each of these n variables can either be 0 or 1. So can be the output. The function can either be constant, where all the inputs get sent to 0 or all the inputs get sent to 1, or balanced, where half the inputs get sent to 0 and the rest get sent to 1.
	To illustrate this, let's take an example of 4 qubits. Which means 2⁴ possible inputs, as each input can either be 0 or 1.
	Let's check each input -
	(0,0,0,0)	(1,1,1,1)
	(0,0,0,1)	(1,1,1,0)
	(0,0,1,0)	(1,1,0,1)
	(0,1,0,0)	(1,0,1,1)
	(1,0,0,0)	(0,1,1,1)
	(0,0,1,1)	(1,1,0,0)
	(1,1,0,0)	(0,0,1,1)
	(0,1,1,0)	(1,0,0,1)
	So if f(0,0,0,0) = 0, f(all 2⁴ combinations) = 0
	or f(0,0,0,0) = 1, f(all 2⁴ combinations) = 1
	The function is constant. Else, it's balanced."

# ╔═╡ 0a3fe586-3ed3-11eb-2c74-010072b424b2
md"Deutsch Algorithm is a case of Deutsch-Josza Algorithm, where the input is 1 qubit. Drawing inspiration from our above circuit, let's make a circuit for n inputs."

# ╔═╡ 2bb45ad8-3ed5-11eb-1647-8f67da7adba6
n = 3

# ╔═╡ 0c3fe0ac-3ed3-11eb-29f4-3f693e893928
begin
	ℱ₀ = chain(n+1)
	plot(ℱ₀)
end

# ╔═╡ 61492872-3ed5-11eb-2830-d33021e6bda6
begin
	ℱ₁ = chain(n+1, [control(k, n+1=>X) for k in 1:n])
	plot(ℱ₁)
end

# ╔═╡ 60951556-3ed5-11eb-0038-2576cbddd974
begin
	ℱ₂ = chain(n+1, chain(n+1, [control(k, n+1=>X) for k in 1:n]), put(n+1=>X))
	plot(ℱ₂)
end

# ╔═╡ 5f80f308-3ed5-11eb-3359-fddaa390f5ae
begin
	ℱ₃ = chain(n+1, put(n+1=>X))
	plot(ℱ₃)
end

# ╔═╡ ff428a74-3ed4-11eb-093d-85e41e536ad3
ℱ = [ℱ₀ ℱ₁ ℱ₂ ℱ₃]

# ╔═╡ d605a88a-3ed4-11eb-2189-7d6bb267dd25
newr = rand(0:3)

# ╔═╡ e8972db6-3ed4-11eb-2e01-0139eff96740
begin
	i = join(ArrayReg(bit"1"), zero_state(n))
	output = i |> chain(n+1, repeat(H, 1:n+1), put(1:n+1=>ℱ[newr+1]), repeat(H,1:n)) |> r->measure(r,1:n)
end

# ╔═╡ 8e9d0e2e-3ed8-11eb-326e-21f2893dadc5
md"If output of the measured qubits is $ $(\"0\"^n) $, then the function is constant, else if, the output is $ $(\"1\"^n) $, the function is balanced"

# ╔═╡ cc44243a-3ed8-11eb-2ef8-d98a6700769b
output == measure(zero_state(n)) ? "Constant" : "Balanced"

# ╔═╡ 0ed997be-4157-11eb-3f71-a16557d72c97
md"Please note, you can change the value of n above to get the Deutsch-Josza Algorithm for different inputs. `` n = 1 `` will give the Deutsch Algorithm."

# ╔═╡ Cell order:
# ╠═7233fedb-3123-4040-97d9-2cd3cd285a67
# ╟─e0cfd952-3d49-11eb-195b-b34f3836e322
# ╟─06c7fa22-3d4a-11eb-2f15-b7ef46b75be3
# ╟─063623ea-3d4a-11eb-3471-ef67e8f31249
# ╠═06657244-3d4d-11eb-3a7c-f79844513992
# ╟─f95221b4-4452-11eb-3cf6-d330a5c5aa11
# ╟─72324dca-4453-11eb-1c26-2b5c1031570b
# ╠═f1f2f676-3d4c-11eb-307b-f7a7408372e0
# ╠═6514718e-3d4d-11eb-088b-9794ec4e7bc4
# ╠═7d3fd99c-3d4d-11eb-0100-67a39154fef1
# ╠═92024f36-3d4d-11eb-3755-1f4538ed6003
# ╠═2ef8834e-3d50-11eb-233c-fd44548a80b9
# ╠═505429aa-3d50-11eb-333f-a75bc433b9ee
# ╟─938c4318-3d50-11eb-3060-31bc11556300
# ╟─559aa08a-3d56-11eb-2103-dd4114a41999
# ╟─c1ca5030-4453-11eb-1814-dd9862cb5aa1
# ╟─f4e5006e-4453-11eb-3a08-a35438601cf1
# ╟─5dec54cc-4454-11eb-33fe-1b90858ea78a
# ╠═5d52c8a8-3d50-11eb-1576-cb5f3485262b
# ╠═a04dbca6-3d50-11eb-24c2-5dbef3190e76
# ╟─8f5ac082-3d58-11eb-3952-8b8a74dc4e26
# ╟─f3568154-3ed2-11eb-125d-33a32b2030a7
# ╟─f5b13e08-3ed2-11eb-0790-d36bf95fe54f
# ╟─0a3fe586-3ed3-11eb-2c74-010072b424b2
# ╠═2bb45ad8-3ed5-11eb-1647-8f67da7adba6
# ╠═0c3fe0ac-3ed3-11eb-29f4-3f693e893928
# ╠═61492872-3ed5-11eb-2830-d33021e6bda6
# ╠═60951556-3ed5-11eb-0038-2576cbddd974
# ╠═5f80f308-3ed5-11eb-3359-fddaa390f5ae
# ╠═ff428a74-3ed4-11eb-093d-85e41e536ad3
# ╠═d605a88a-3ed4-11eb-2189-7d6bb267dd25
# ╠═e8972db6-3ed4-11eb-2e01-0139eff96740
# ╟─8e9d0e2e-3ed8-11eb-326e-21f2893dadc5
# ╠═cc44243a-3ed8-11eb-2ef8-d98a6700769b
# ╟─0ed997be-4157-11eb-3f71-a16557d72c97


text:
using Genie.Configuration, Logging

const config = Settings(
    server_port=8000,
    server_host="127.0.0.1",
    log_level=Logging.Debug,
    log_to_file=true,
    server_handle_static_files=true
)

ENV["JULIA_REVISE"] = "off"

max_stars_repo_name: MLH-Fellowship/ZoomCord

max_issues_repo_name: MLH-Fellowship/ZoomCord

max_forks_repo_name: MLH-Fellowship/ZoomCord


text:
__precompile__()

module MyGraph

using DataStructures, LightGraphs

import Base: ==, isless, getindex, setindex!, isempty, isequal, eltype

import LightGraphs: AbstractGraph, edgetype,has_edge,has_vertex,inneighbors,ne,
                    nv,outneighbors,vertices,is_directed

export  AbstractDictGraph,
        DictGraph

abstract type AbstractDictGraph{V,E} <: AbstractGraph{V} end

mutable struct DictGraph{V,E} <: AbstractDictGraph{V,E}
    is_directed::Bool
    is_weighted::Bool
    numvertices::Int
    numedges::Int
    vvwr::Dict{V,Dict{V,E}} #vertex => {(vertex,weight),...} relation

    function DictGraph{V,E}(;isdirected = false, isweighted = false, vvwr = Dict{V, Dict{V,E}}()) where {V,E}
        new{V,E}(isdirected, isweighted, 0, 0, vvwr)
    end
end

function DictGraph(G::DictGraph)
    DictGraph(G.isdirected, G.isweighted, G.numvertices, G.numedges, Dict(G.vvwr))
end

#Basic operations over graph.
vertices(G::AbstractDictGraph) = collect(keys(G.vvwr))

function add_vertex!(G::AbstractDictGraph{V,E}, vertices::Vararg{V,N} where N) :: Nothing where {V,E}
    for vertex in vertices
        (!haskey(G.vvwr, vertex) && push!(G.vvwr, vertex => Dict{V,E}())) && G.numvertices += 1
    end
    nothing
end

function add_vertex(G::AbstractDictGraph{V,E}, vertices::Vararg{UInt, N} where N)::DictGraph where {V,E}
    T::DictGraph = DictGraph(G)
    add_vertex!(T,vertices...)
    T
end

function rem_vertex!(G::AbstractDictGraph{V,E}, vertices::Vararg{V,N} where N)::Nothing where {V,E}
    function remvertex!(G::AbstractDictGraph{V,E}, vertex::V) :: Nothing where {V,E}
        try
            G[v];
            for k in vertices(G)
                k != v && delete!(G[k],v)
            end

            delete!(G.vvwr, v)
            G.numvertices -= 1
        catch e
            if isa(e, KeyError)
                println("vertex $v does not belongs to this graph.")
            else
                error("G: unknown error.")
            end
        end # try
    end
    foreach(vertex -> remvertex!(G,vertex), vertices)
end

function rem_vertex(G::DictGraph, vertices::Vararg{UInt,N} where N)::DictGraph
    T = DictGraph(G)
    rem_vertex!(T,vertices...)
    T
end

function addedge!(G::DictGraph, source::Int, target::Int; weight::Real = 1.0) :: Nothing
    add_vertex!(G, source)
    add_vertex!(G, target)
    push!(G[source], target => weight)
    !G.isdirected && push!(G[target], source => weight)
    G.numedges += 1
end

function addedge!(G::DictGraph, edge::Tuple{UInt,UInt,Real})::Nothing
    addedge!(G,edge[1],edge[2]; weight = edge[3])
end

function addedge!(G::DictGraph, edges::Vararg{Tuple{UInt,UInt,Real}, N} where N)::Nothing
    foreach(edge -> addedge!(G,edge), edges)
end

function addedge(G::DictGraph, source::Int, target::Int; weight::Real = 1.0)::DictGraph
    T = DictGraph(G);
    addedge!(T, v, u; weight = weight);
    T
end

function addedge(G::DictGraph, edges::Vararg{Tuple{UInt,UInt,Real}, N} where N)::DictGraph
    T = DictGraph(G)
    addedge!(T,edges...)
    T
end

function getindex(G::DictGraph, i::Int)
    return G.vvwr[i]
end

function getedgeweight(G::DictGraph, source::Int, target::Int)::Real
    !G.isweighted || error("Graph is not weghted.")
    haskey(G.vvwr, source) || error("vertex $source does not belongs to this graph.")
    haskey(G.vvwr, source) || error("vertex $target does not belongs to this graph.")
    G[source][target]
end

function lowestweightedge(G::DictGraph, vertex::Int, only::Vararg{Int,N} where N)::Pair{Int, Real}
    !G.isweighted || error("Graph $(G.label) is not weghted.")
    haskey(G.vvwr) || error("the vertex $vertex does not belongs to this graph.")
    E = Set(only)
    foreach(e -> haskey(G[e]) || error("the vertex $vertex does not belongs to this graph."), S)
    t = Pair(-1, Inf32);
    for kv in G[vertex]
         kv[1] ∈ E && kv[2] < t[2] && (t = kv)
    end
    t
end

function minweightedge(G::DictGraph)::Pair{Real, Tuple{Int,Int}}
    D = Array{Pair{Real, Tuple{Int,Int}}}()
    for v in vertices(G)
        x = minweightedgefromvertex(G,v)
        push!(D, x[2] => (v, x[1]))
    end
    min(D)
end

include("connectedgraph.jl")
#algorihms
include("algorithms/dijkstra.jl")
include("algorithms/DFS.jl")
include("algorithms/BFS.jl")
include("algorithms/printpath.jl")
include("algorithms/kruskal.jl")
include("algorithms/primjarnik.jl")
include("algorithms/alg1.jl")

end #module → ↔


text:
using Serialization;
using Images, ImageTransformations, Interpolations;

dim = 10
view = 3

pixels = 100.0 * randn(Float32, dim, dim)

println("; dims: ", size(pixels))

@time res = Float32.(imresize(pixels, (view, view)))

# println(res, "; dims: ", size(res))
display(res)

data = rand(dim, dim)
mask = map(x -> (x > 0.5), data)
println(mask, typeof(mask))

res = Bool.(imresize(mask, (view, view), method = Constant()))
println(res, typeof(res))
display(res)

view_pixels = deserialize("/Users/chris/view_pixels.bin")
view_mask = deserialize("/Users/chris/view_mask.bin")

view_width = Int64(372)
view_height = Int64(372)

display(view_pixels[1:5, 1:5])
display(view_mask[1:5, 1:5])

view_pixels = Float32.(imresize(view_pixels, (view_width, view_height)),)
view_mask =
    Bool.(
        imresize(
            view_mask,
            (Integer(view_width), Integer(view_height)),
            method = Constant(),
        ),
    ) # use Nearest-Neighbours for the mask