text:
runbenchmark(::Val{B}; kwargs...) where {B} =
  error("Benchmark $B not found")

# If you call runbenchmarks(stdout; kwargs...) it outputs to the console
function runbenchmarks(io; write_results = OPTIONS["write_results"],
                           test = OPTIONS["test"],
                           blas_num_threads = OPTIONS["blas_num_threads"],
                           blocksparse_num_threads = OPTIONS["blocksparse_num_threads"],
                           maxdims = OPTIONS["maxdims"],
                           cpp_or_julia = OPTIONS["cpp_or_julia"],
                           benchmarks = OPTIONS["benchmarks"],
                           cpp_itensor_version = OPTIONS["cpp_itensor_version"],
                           cpp_itensor_dir = OPTIONS["cpp_itensor_dir"],
                           julia_itensor_version = OPTIONS["julia_itensor_version"],
                           splitblocks = OPTIONS["splitblocks"])
  println(stdout, "Running benchmark with the following options:")
  if io == stdout
    println(stdout, "io = stdout")
  else
    @show io
  end
  @show write_results
  @show test
  @show blas_num_threads
  @show blocksparse_num_threads
  print(stdout, "maxdims = ")
  show(stdout, MIME("text/plain"), maxdims)
  println(stdout)
  @show cpp_or_julia
  @show benchmarks
  @show cpp_itensor_version
  @show julia_itensor_version
  @show splitblocks
  println(stdout)

  if julia_itensor_version ≠ ITensors.version()
    error("Attempting to benchmark with ITensors.jl v$julia_version, but the current version is v$(ITensors.version()). You should fix the ITensors.jl version by running `] add ITensors@$julia_version` at the Julia command prompt.")
  end

  warning_seperator = "X"^70

  println(io)
  println(stdout)
  println(io, warning_seperator)
  println(stdout, warning_seperator)
  if write_results
    println(io, "XXX WARNING: benchmark results are set to be written to disk, may overwrite previous results")
    println(stdout, "XXX WARNING: benchmark results are set to be written to disk, may overwrite previous results")
  else
    println(io, "XXX WARNING: benchmark results are not set to be written to disk.")
    println(stdout, "XXX WARNING: benchmark results are not set to be written to disk.")
  end
  println(io, warning_seperator)
  println(stdout, warning_seperator)
  println(io)
  println(stdout)

  seperator = "#"^70

  outputlevel = 1

  println(io)
  println(io, "Using Julia v$(VERSION)")

  println(io)
  println(io, "Using ITensors.jl v$julia_itensor_version")

  println(io)
  println(io, "Using C++ ITensor v$cpp_itensor_version")

  println(io)
  println(io, "Julia is using BLAS vendor $(BLAS.vendor()). Please double check this is the same as the one you compile C++ ITensor against.")
  println(io, "Using $blas_num_threads BLAS threads")

  # Disable Strided threads
  ITensors.Strided.disable_threads()

  println(io)
  println(io, "Benchmarking $benchmarks")
  println(io)

  println(io, "Bond dimensions set to:")
  show(io, MIME("text/plain"), maxdims)
  println(io)
  println(io)

  for benchmark in benchmarks
    println(io, "X"^70)
    println(io, "X"^70)
    println(io)
    maxdims_benchmark = maxdims[benchmark]
    if test ≠ false
      println(io, "XXX TESTING: only running the benchmarks for the specified set of bond dimensions maxdims[$test] = $(maxdims_benchmark[test]).")
      maxdims_benchmark = maxdims_benchmark[test]
      println(io)
    end
    st = "Running benchmark $benchmark for bond dimensions $maxdims_benchmark, BLAS threads $blas_num_threads, and block sparse threads $blocksparse_num_threads"
    io ≠ stdout && println(stdout, "\n" * st)
    println(io, st)
    println(io)
    # Print the more detailed description of the benchmark
    println(io, descriptions[benchmark])
    println(io)
    valbenchmark = Val(Symbol(benchmark))
    src_dir = joinpath(pkgdir(@__MODULE__), "src")
    data_dir = joinpath(pkgdir(@__MODULE__), "data")
    benchmark_dir = joinpath(src_dir, "bench_$benchmark")
    for maxdim in maxdims_benchmark,
        blas_num_thread in blas_num_threads,
        blocksparse_num_thread in blocksparse_num_threads

      # Path within data/ directory to write results
      writepath = "bench_$benchmark"
      writepath = joinpath(writepath, "blas_num_threads_$(blas_num_thread)")
      writepath = joinpath(writepath, "blocksparse_num_threads_$(blocksparse_num_thread)")
      if splitblocks
        writepath = joinpath(writepath, "splitblocks_$(splitblocks)")
      end

      # Check consistency of threads that are set
      if blocksparse_num_thread == 1
        ITensors.disable_threaded_blocksparse()
      elseif blocksparse_num_thread ≠ Threads.nthreads()
        error("Trying to benchmark with $blocksparse_num_thread number of threads for block sparse operations, but Julia was started with $(Threads.nthreads()) threads. Start Julia with $blocksparse_num_thread threads using `JULIA_NUM_THREADS=$blocksparse_num_thread julia`, `julia --threads=$blocksparse_num_thread`, or `julia -t $blocksparse_num_thread`.")
      else
        ITensors.enable_threaded_blocksparse()
      end

      if blocksparse_num_thread > 1 && blas_num_thread > 1
        error("Trying to benchmark with $blocksparse_num_thread number of threads for block sparse operations and $blas_num_thread BLAS threads. These two types of multithreading conflict with each other, and therefore is not recommended to use.")
      end
      BLAS.set_num_threads(blas_num_thread)

      st = "Maximum bond dimension set to $maxdim, BLAS threads set to $blas_num_thread, block sparse threads set to $blocksparse_num_thread, and splitblocks is set to $splitblocks."
      println(io, "Benchmark located in path $benchmark_dir")
      println(io)
      io ≠ stdout && println(stdout, " " * st)
      println(io, st)
      println(io)
      if isnothing(cpp_or_julia) || cpp_or_julia == "julia"
        st = "Run Julia benchmark $benchmark for bond dimension $maxdim, $blas_num_thread BLAS threads, $blocksparse_num_thread block sparse threads, and splitblocks $splitblocks."
        println(io, seperator)
        io ≠ stdout && println(stdout, " "^2 * st)
        println(io, st)
        println(io)

        # Trigger compilation
        comp_time = runbenchmark(valbenchmark; maxdim = 10, nsweeps = 2,
                                 outputlevel = 0, splitblocks = splitblocks)
        if isnothing(comp_time)
          println(io, "Moving on to the next benchmark.")
          println(stdout, "Moving on to the next benchmark.")
          continue
        end

        local time
        local maxdim_
        output = @capture_out begin
          time = @elapsed maxdim_ = runbenchmark(valbenchmark; maxdim = maxdim, outputlevel = outputlevel, splitblocks = splitblocks)
        end
        println(io, output)

        println(io)
        println(io, "Maximum dimension = $maxdim_")
        st = "Total runtime = $time seconds"
        io ≠ stdout && println(stdout, " "^3 * st)
        println(io, st)
        println(io)

        if write_results
          filepath = joinpath(data_dir, "julia")
          filepath = joinpath(filepath, "v$julia_itensor_version")
          filepath = joinpath(filepath, writepath)
          filepath = joinpath(filepath, "maxdim_$(maxdim).txt")
          println(io, "Writing results to $filepath")
          println(io)
          mkpath(dirname(filepath))
          writedlm(filepath, time)
        end
      end
      if !splitblocks && (isnothing(cpp_or_julia) || cpp_or_julia == "c++")
        st = "Run C++ benchmark $benchmark for bond dimension $maxdim, $blas_num_thread BLAS threads, $blocksparse_num_thread block sparse threads, and splitblocks $splitblocks."
        println(io, seperator)
        io ≠ stdout && println(stdout, " "^2 * st)
        println(io, st)
        println(io)

        # Create the file itensor_dir.mk that has the directory where
        # the version of C++ ITensor lives
        if !isfile(joinpath(cpp_itensor_dir, "itensor", "core.h"))
          error("Attempting to use ITensor v$cpp_itensor_version. The installation location is currently set to the directory $cpp_itensor_dir. If you have C++ ITensor v$cpp_itensor_version installed already, you can point to the location by setting `cpp_itensor_dir` to the full path of the installation. Otherwise, download and compile C++ ITensor v$cpp_itensor_version using the direction in http://www.itensor.org/docs.cgi?page=install&vers=cppv3 (and check out the version with `git checkout v$cpp_itensor_version before building).")
        end

        # Make a random directory to store tempororary build files
        rand_dir = joinpath(src_dir, ".tmp_" * randstring(16))
        mkdir(rand_dir)

        open(joinpath(rand_dir, "itensor_dir.mk"), "w") do io
          write(io, "LIBRARY_DIR=$cpp_itensor_dir")
        end
        # Generate the C++ function we will compile
        open(joinpath(rand_dir, "runbenchmark.cc"), "w") do io
          cpp_file =
          """
          #include "$(joinpath(src_dir, "bench_$benchmark", "runbenchmark.h"))"
          int
          main(int argc, char *argv[])
            {
            auto maxdim = std::stoi(argv[1]);
            auto outputlevel = 0;
            if(argc > 2)
              outputlevel = std::stoi(argv[2]);
            auto maxdim_ = run({"Maxdim = ", maxdim,
                                "OutputLevel = ", outputlevel});
            return 0;
            }
          """
          write(io, cpp_file)
        end
        # Trigger rebuild
        println(io, "Compile the C++ benchmark in temporary directory $rand_dir")
        cp(joinpath(src_dir, "Makefile"), joinpath(rand_dir, "Makefile"))
        output = @capture_out begin
          run(`make --directory=$rand_dir`)
        end
        println(io, output)
        println(io)
        rm(joinpath(rand_dir, "runbenchmark.cc"))
        runbenchmark_exec = joinpath(rand_dir, "runbenchmark")
        open(joinpath(rand_dir, "runbenchmark.sh"), "w") do io
          write(io, """
                    #!/bin/bash
                    MKL_NUM_THREADS=$blas_num_thread OPENBLAS_NUM_THREADS=$blas_num_thread OMP_NUM_THREADS=$blocksparse_num_thread $runbenchmark_exec $maxdim $outputlevel
                    """)
        end
        chmod(joinpath(rand_dir, "runbenchmark.sh"), 0o777)
        local time
        println(io, "Run the C++ benchmark")
        output = @capture_out begin
          time = @elapsed run(`$(joinpath(rand_dir, "runbenchmark.sh"))`)
        end
        println(io, output)

        # Clean up the build files
        println(io, "Cleaning up temporary directory $rand_dir used for building C++ benchmark")
        rm(joinpath(rand_dir, "runbenchmark.sh"))
        rm(joinpath(rand_dir, "runbenchmark"))
        rm(joinpath(rand_dir, "runbenchmark.o"))
        rm(joinpath(rand_dir, "itensor_dir.mk"))
        rm(joinpath(rand_dir, "Makefile"))
        rm(rand_dir)

        println(io)
        println(io, "Maximum dimension = $maxdim")
        st = "Total runtime = $time seconds"
        io ≠ stdout && println(stdout, " "^3 * st)
        println(io, st)
        println(io)

        if write_results
          # TODO: add version number to data file name
          filepath = joinpath(data_dir, "c++")
          filepath = joinpath(filepath, "v$cpp_itensor_version")
          filepath = joinpath(filepath, writepath)
          filepath = joinpath(filepath, "maxdim_$(maxdim).txt")
          println(io, "Writing results to $filepath")
          println(io)
          mkpath(dirname(filepath))
          writedlm(filepath, time)
        end
      end
    end
  end
  return nothing
end

function logfile_name(rootdir, benchmarks, blas_num_threads,
                      blocksparse_num_threads, splitblocks)
  logfile = rootdir
  logfile = joinpath(logfile, "logs")
  logfile = joinpath(logfile, "$(now())")
  logfile *= "_benchmark_$(benchmarks)"
  logfile *= "_blas_num_threads_$(blas_num_threads)"
  logfile *= "_blocksparse_num_threads_$(blocksparse_num_threads)"
  if splitblocks
    logfile *= "_splitblocks_$(splitblocks)"
  end
  logfile *= ".log"
  return logfile
end

# This version outputs to a file
function runbenchmarks(; blas_num_threads = OPTIONS["blas_num_threads"],
                         blocksparse_num_threads = OPTIONS["blocksparse_num_threads"],
                         benchmarks = OPTIONS["benchmarks"],
                         splitblocks = OPTIONS["splitblocks"],
                         logfile = logfile_name(pkgdir(@__MODULE__), benchmarks, blas_num_threads, blocksparse_num_threads, splitblocks),
                         kwargs...)
  println(stdout, "Saving information about the benchmark run to the file $logfile")
  println(stdout)
  open(logfile, "w") do io
    runbenchmarks(io; blas_num_threads = blas_num_threads,
                      blocksparse_num_threads = blocksparse_num_threads,
                      benchmarks = benchmarks,
                      splitblocks = splitblocks, kwargs...)
  end # open
  return nothing
end


text:
using CompScienceMeshes
using Test

vrts = [
    point(0.0, 0.0, 0.0),
    point(1.0, 0.0, 0.0),
    point(0.0, 1.0, 0.0),
]

idcs = [
    index(1,2,3)
]

axis = 0.5π * [1,0,0]
G1 = Mesh(vrts, idcs)
G2 = rotate(G1, axis)
G3 = rotate(G1, 2*axis)

G = weld(G1, G2);
G = weld(G1, G2, G3)

@test numcells(G) == numcells(G1) + numcells(G2) + numcells(G3)
@test numvertices(G) == 5

edges = skeleton(G,1)
pairs = cellpairs(G, edges, dropjunctionpair=true)

@test size(pairs) == (2,8)

mesh = meshrectangle(1.0, 1.0, 0.5)
edges = skeleton(mesh,1)
cps = cellpairs(mesh, edges)
@test size(cps) == (2,16)

# width, height = 1.0, 0.5; h = 0.1
# mesh = meshrectangle(width, height, h)
# edges = skeleton(mesh,1)
# cps = cellpairs(mesh, edges, dropjunctionpair=true)

ext: jl


text:
function publish_tweet!(
    state::Tuple{AbstractGraph, AbstractArray}, tweet_list::AbstractArray, agent_idx::Integer,
    tick_nr::Integer=0
)
    graph, agent_list = state
    this_agent = agent_list[agent_idx]
    tweet_opinion = this_agent.opinion + 0.1 * (2 * rand() - 1)
    # upper opinion limit is 1
    if tweet_opinion > 1
        tweet_opinion = 1.0
    # lower opinion limit is -1
    elseif tweet_opinion < -1
        tweet_opinion = -1.0
    end
    tweet = Tweet(tweet_opinion, length(outneighbors(graph, agent_idx)), agent_idx, tick_nr)
    push!(tweet_list, tweet)
    # send tweet to each outneighbor
    for neighbor in outneighbors(graph, agent_idx)
        push!(agent_list[neighbor].feed, tweet)
    end
    return state, tweet_list
end

max_stars_repo_name: digitalemuendigkeit/twitter-abm

max_issues_repo_name: digitalemuendigkeit/twitter-abm

max_forks_repo_name: digitalemuendigkeit/twitter-abm


text:
function cfg_net(
    ;
    agent_count::Int64=100,
    m0::Int64=10,
    growth_rate::Int64=4,
    new_follows::Int64=4,
    initial_follows::Int64=4
    )
    return (agent_count=agent_count,m0=m0,growth_rate=growth_rate,new_follows=new_follows, initial_follows=initial_follows)
end

function cfg_sim(
    ;
    n_iter::Int64=100,
    max_inactive_ticks::Int64=2
    )
    return (n_iter=n_iter,max_inactive_ticks=max_inactive_ticks)
end

function cfg_ot(
    ;
    like::Float64 = 0.2,
    retweet::Float64 = 0.3,
    backfire::Float64 = 0.4,
    check_unease::Float64 = 0.3,
    follow::Float64 = 0.2,
    unfollow::Float64 = 0.4
    )
    return (like=like,retweet=retweet,backfire=backfire,check_unease=check_unease,follow=follow,unfollow=unfollow)
end

function cfg_ag(
    ;
    own_opinion_weight::Float64 = 0.99,
    check_decrease::Float64 = 0.9,
    inclin_interact_lambda::Float64 = log(25),
    unfollow_rate::Float64 = 0.2,
    )
    return (own_opinion_weight=own_opinion_weight,check_decrease=check_decrease,inclin_interact_lambda=inclin_interact_lambda,unfollow_rate=unfollow_rate)
end

function cfg_feed(
    ;
    feed_size::Int64 = 10,
    tweet_decay::Float64 = 0.5
    )
    return (feed_size = feed_size, tweet_decay = tweet_decay)
end

struct Config
    network::NamedTuple{(:agent_count, :m0, :growth_rate, :new_follows, :initial_follows),NTuple{5,Int64}}
    simulation::NamedTuple{(:n_iter, :max_inactive_ticks),NTuple{2,Int64}}
    opinion_treshs::NamedTuple{(:like, :retweet, :backfire, :check_unease, :follow, :unfollow),NTuple{6,Float64}}
    agent_props::NamedTuple{(:own_opinion_weight, :check_decrease, :inclin_interact_lambda, :unfollow_rate),NTuple{4,Float64}}
    feed_props::NamedTuple{(:feed_size, :tweet_decay), <:Tuple{Int64, Float64}}

    function Config(;network = cfg_net(), simulation = cfg_sim(), opinion_treshs = cfg_ot(), agent_props = cfg_ag(), feed_props = cfg_feed())

        new(network, simulation, opinion_treshs, agent_props, feed_props)
    end
end


text:
#####################################################
#nonlinearities

function relu!(inp, outp)
	for i in 1:length(inp)
		outp[i] = max(0, inp[i])
	end
end

function relu!(inp, outp, t)
	for i in 1:length(inp)
		outp[i] = max(0, inp[i] - t[i])
	end
end

function relu_diff!(ax, error)
	for i in 1:length(ax)
		error[i] *= ax[i] > 0
	end
end

function lin!(inp, outp)
	for i in 1:length(inp)
		outp[i] = inp[i]
	end
end

function lin_diff!(ax, error)
end

function sigm!(inp, outp)
	for i in 1:length(inp)
		outp[i] = 1. / (1. + exp(-inp[i]))
	end
end

function sigm_diff!(ax, error)
	for i in 1:length(ax)
		error[i] *= 1. / (1 + exp(-ax[i]))*(1. - 1. / (1 + exp(-ax[i])))
	end
end

function pwl!(inp, outp, t)
	for i in 1:length(inp)
		outp[i] = clamp(inp[i] - t[i],0.,1.) #piece-wise linear
	end
end

function sigm!(inp, outp, t; λ = 10.) #10
	for i in 1:length(inp)
		outp[i] = 1. / (1. + exp(- λ * (inp[i] - t[i])))
	end
end


text:
export circ, circ!
export quadratic!, quadratic

"""
    circ(radius_aperture, x, y)

Respresenting an aperture.
Being 0 if `x^2 + y^2 > radius_aperture^2` and otherwise 1.

# Examples
```julia-repl
julia> circ(1, 0, 1)
1

julia> circ(1, 0.5, 0.5)
1

julia> circ(1, 1, 1)
0
```

"""
function circ(radius_aperture, x, y)
    if x^2 + y ^2 > radius_aperture^2
        return 0
    else
        return 1
    end
end

"""
    circ(radius_aperture, r)

Respresenting an aperture.
Being 0 if `r^2 > radius_aperture^2` and otherwise 1.

# Examples
```julia-repl
julia> circ(1, 0)
1

julia> circ(1, 1)
1

julia> circ(1, 1.01)
0
```
"""
function circ(radius_aperture, r)
    if r^2 > radius_aperture ^2
        return 0
    else
        return 1
    end
end

"""
    circ(arr, L, radius_aperture)

Apply circular aperture of radius `radius_aperture` to an
array which has field width of `L`.
`circ!` is also available.

```jldoctest
julia> circ(ones((5, 5)), 5, 2.5)
5×5 Array{Float64,2}:
 0.0  0.0  1.0  0.0  0.0
 0.0  1.0  1.0  1.0  0.0
 1.0  1.0  1.0  1.0  1.0
 0.0  1.0  1.0  1.0  0.0
 0.0  0.0  1.0  0.0  0.0
```
"""
function circ(arr::AbstractArray, L, radius_aperture)
    radius_arr = Zygote.@ignore PhysicalOptics.rr(size(arr), typeof(arr), L=L)
    pupil = Zygote.@ignore (radius_aperture .≥ radius_arr)
    arr_new = arr .* pupil
    return arr_new
end

function circ!(arr::AbstractArray, L, radius_aperture)
    radius_arr = PhysicalOptics.rr(size(arr), L=L)
    arr .*= (radius_aperture .≥ radius_arr)
    return arr
end

"""
    quadratic(arr, L, diameter, [, diameter=(0,0)])

Apply a quadratic aperture of `diameter` to an
array which has a field width of `L`.
The aperture is shifted by `Δx, Δy` with respect to the center

There is also an in-place version `quadratic!`

```jldoctest
julia> quadratic(ones((4, 4)), 10, 5)
4×4 Array{Float64,2}:
 0.0  0.0  0.0  0.0
 0.0  1.0  1.0  1.0
 0.0  1.0  1.0  1.0
 0.0  1.0  1.0  1.0

julia> quadratic(ones((5, 5)), 10, 5)
5×5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0
 0.0  1.0  1.0  1.0  0.0
 0.0  1.0  1.0  1.0  0.0
 0.0  1.0  1.0  1.0  0.0
 0.0  0.0  0.0  0.0  0.0

julia> quadratic(ones((5, 5)), 5, 2.5, (0.0, 1.25))
5×5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  1.0  1.0
 0.0  0.0  1.0  1.0  1.0
 0.0  0.0  1.0  1.0  1.0
 0.0  0.0  0.0  0.0  0.0
```
"""
function quadratic(arr, L, diameter,  offset::NTuple{2, <: Number}=(0, 0))
    return quadratic!(copy(arr), L, diameter, offset)
end

function quadratic!(arr, L, diameter, offset::NTuple{2, <: Number}=(0,0))
    # need to take abs to compare with diameter
    x = abs.(-offset[2] .+ fftpos(L[2], size(arr, 2)))'
    y = abs.(-offset[1] .+ fftpos(L[1], size(arr, 1)))

    arr .*= (diameter[2] / 2 .≥ x)
    arr .*= (diameter[1] / 2 .≥ y)

    return arr
end

function quadratic!(arr, L::Number, diameter::Number, offset::NTuple{2, <:Number}=(0,0))
    return quadratic!(arr, (L, L), (diameter, diameter), offset)
end

max_stars_repo_name: JuliaPhysics/PhysicalOptics

max_issues_repo_name: JuliaPhysics/PhysicalOptics

max_forks_repo_name: JuliaPhysics/PhysicalOptics


text:
"""Type to record neighbor labels and their counts."""
mutable struct NeighComm
    neigh_pos::Vector{Int}
    neigh_cnt::Vector{Int}
    neigh_last::Int
end

mutable struct NodeStatus
    neigh_pos::Vector{Int}
    neigh_cnt::Vector{Bool}
    neigh_last::Int
end

function range_shuffle!(r, a::AbstractVector)
    @inbounds for i=r:-1:2
        j = StatsBase.randi(i)
        a[i],a[j] = a[j],a[i]
    end
end

"""Return the most frequency label."""
function pre_vote!(g, m, c::NeighComm, u, node_status)
    @inbounds for i=1:c.neigh_last-1
        c.neigh_cnt[c.neigh_pos[i]] = -1
    end
    c.neigh_last = 1
    c.neigh_pos[1] = m[u]
    c.neigh_cnt[c.neigh_pos[1]] = 0
    c.neigh_last = 2
    max_cnt = 0
    for neigh in neighbors(g,u)
        #if in(neigh, all_active_nodes)
        if node_status.neigh_cnt[neigh]
            neigh_comm = m[neigh]
            if c.neigh_cnt[neigh_comm] < 0
                c.neigh_cnt[neigh_comm] = 0
                c.neigh_pos[c.neigh_last] = neigh_comm
                c.neigh_last += 1
            end
            c.neigh_cnt[neigh_comm] += 1
            if c.neigh_cnt[neigh_comm] > max_cnt
                max_cnt = c.neigh_cnt[neigh_comm]
            end
        end
    end
    # ties breaking randomly
    range_shuffle!(c.neigh_last-1, c.neigh_pos)
    for lbl in c.neigh_pos
        if c.neigh_cnt[lbl] == max_cnt
            return lbl
        end
    end
end

function pre_update!(g, m, c::NeighComm, active_nodes, random_order, node_status)
    while !isempty(active_nodes)
        #random_order = collect(active_nodes)
        i = 0
        for u in active_nodes
            i += 1
            random_order[i] = u
        end
        #shuffle!(random_order)
        range_shuffle!(i, random_order)
        for j=1:i
            u = random_order[j]
            old_comm = m[u]
            new_comm = pre_vote!(g, m, c, u, node_status)
            if new_comm != old_comm
                for v in neighbors(g,u)
                    #if in(v, all_active_nodes)
                    if node_status.neigh_cnt[v]
                        push!(active_nodes, v)
                        m[u] = new_comm
                    end
                end
            else
                delete!(active_nodes, u)
            end
        end
    end
end

"""Return the most frequency label."""
function vote!(g, m, c::NeighComm, u)
    @inbounds for i=1:c.neigh_last-1
        c.neigh_cnt[c.neigh_pos[i]] = -1
    end
    c.neigh_last = 1
    c.neigh_pos[1] = m[u]
    c.neigh_cnt[c.neigh_pos[1]] = 0
    c.neigh_last = 2
    max_cnt = 0
    for neigh in neighbors(g,u)
        neigh_comm = m[neigh]
        if c.neigh_cnt[neigh_comm] < 0
            c.neigh_cnt[neigh_comm] = 0
            c.neigh_pos[c.neigh_last] = neigh_comm
            c.neigh_last += 1
        end
        c.neigh_cnt[neigh_comm] += 1
        if c.neigh_cnt[neigh_comm] > max_cnt
            max_cnt = c.neigh_cnt[neigh_comm]
        end
    end
    # ties breaking randomly
    range_shuffle!(c.neigh_last-1, c.neigh_pos)
    for lbl in c.neigh_pos
        if c.neigh_cnt[lbl] == max_cnt
            return lbl
        end
    end
end

function update!(g, m, c::NeighComm, active_nodes, random_order)
    num_active_nodes = 0
    while !isempty(active_nodes)
        if length(active_nodes) > num_active_nodes
            num_active_nodes = length(active_nodes)
        end
        #random_order = collect(active_nodes)
        i = 0
        for u in active_nodes
            i += 1
            random_order[i] = u
        end
        #shuffle!(random_order)
        range_shuffle!(i, random_order)
        for j=1:i
            u = random_order[j]
            old_comm = m[u]
            new_comm = vote!(g, m, c, u)
            if new_comm != old_comm
                for v in neighbors(g,u)
                    push!(active_nodes, v)
                    m[u] = new_comm
                end
            else
                delete!(active_nodes, u)
            end
        end
    end
    num_active_nodes
end

function lpa_addnode!(g, m, u)
    if !hasnode(g, u)
        addnode!(g, u)
        m[u] = u
    end
end

function lpa_deletenode!(g, m, c, u, active_nodes, random_order, node_status)
    if hasnode(g, u)
        empty!(active_nodes)
        u_comm = m[u]
        delnode!(g, u)
        #delete!(m, u)
        m[u] = 0
        node_status.neigh_last = 0
        for v in nodes(g)
            if m[v] == u_comm
                push!(active_nodes, v)
                node_status.neigh_last += 1
                node_status.neigh_pos[node_status.neigh_last] = v
                node_status.neigh_cnt[v] = true
                m[v] = v
            end
        end
        pre_update!(g, m, c, active_nodes, random_order, node_status)
        for i=1:node_status.neigh_last
            push!(active_nodes, node_status.neigh_pos[i])
            node_status.neigh_cnt[node_status.neigh_pos[i]] = false
        end
        return update!(g, m, c, active_nodes, random_order)
    end
    return 0
end

function lpa_addedge!(g, m, c, u, v, active_nodes, random_order, node_status)
    lpa_addnode!(g, m, u)
    lpa_addnode!(g, m, v)
    if !hasedge(g, u, v)
        u_comm = m[u]
        v_comm = m[v]
        if u_comm != v_comm
            ku_int = 0
            ku_ext = 0
            kv_int = 0
            kv_ext = 0
            for neigh in neighbors(g,u)
                if m[neigh] == u_comm
                    ku_int += 1
                else
                    ku_ext += 1
                end
            end
            for neigh in neighbors(g,v)
                if m[neigh] == v_comm
                    kv_int += 1
                else
                    kv_ext += 1
                end
            end
            if ku_ext+1 >= ku_int || kv_ext+1 >= kv_int
                #empty!(active_nodes)
                #empty!(active_nodes)
                node_status.neigh_last = 0
                for i in nodes(g)
                    i_comm = m[i]
                    if i_comm == u_comm || i_comm == v_comm
                        push!(active_nodes, i)
                        #push!(active_nodes, i)
                        node_status.neigh_last += 1
                        node_status.neigh_pos[node_status.neigh_last] = i
                        node_status.neigh_cnt[i] = true
                        m[i] = i
                    end
                end
                addedge!(g, u, v)
                pre_update!(g, m, c, active_nodes, random_order, node_status)
                #empty!(pre_active_nodes)
                for i=1:node_status.neigh_last
                    push!(active_nodes, node_status.neigh_pos[i])
                    node_status.neigh_cnt[node_status.neigh_pos[i]] = false
                end
                return update!(g, m, c, active_nodes, random_order)
            else
                addedge!(g, u, v)
            end
        else
            addedge!(g, u, v)
        end
    end
    return 0
end

function lpa_deleteedge!(g, m, c, u, v, active_nodes, random_order, node_status)
    if hasedge(g, u, v)
        u_comm = m[u]
        v_comm = m[v]
        if u_comm == v_comm
            empty!(active_nodes)
            node_status.neigh_last = 0
            for i in nodes(g)
                if m[i] == u_comm
                    #push!(pre_active_nodes, i)
                    push!(active_nodes, i)
                    node_status.neigh_last += 1
                    node_status.neigh_pos[node_status.neigh_last] = i
                    node_status.neigh_cnt[i] = true
                    m[i] = i
                end
            end
            deledge!(g, u, v)
            pre_update!(g, m, c, active_nodes, random_order, node_status)
            for i=1:node_status.neigh_last
                push!(active_nodes, node_status.neigh_pos[i])
                node_status.neigh_cnt[node_status.neigh_pos[i]] = false
            end
            return update!(g, m, c, active_nodes, random_order)
        else
            deledge!(g, u, v)
        end
    end
    return 0
end


text:
"""
    TagBot(; destination="TagBot.yml", registry=nothing, dispatch=false)

Adds GitHub release support via [TagBot](https://github.com/JuliaRegistries/TagBot).

## Keyword Arguments
- `destination::AbstractString`: Destination of the workflow file,
  relative to `.github/workflows`.
- `registry::Union{AbstractString, Nothing}`: Custom registry, in the format `owner/repo`.
- `dispatch::Bool`: Whether or not to enable the `dispatch` option.
"""
@with_kw_noshow struct TagBot <: BasicPlugin
    destination::String = "TagBot.yml"
    registry::Union{String, Nothing} = nothing
    dispatch::Bool = false
end

source(::TagBot) = default_file("github", "workflows", "TagBot.yml")
destination(p::TagBot) = joinpath(".github", "workflows", p.destination)
tags(::TagBot) = "<<", ">>"

view(p::TagBot, ::Template, ::AbstractString) = Dict(
    "HAS_DISPATCH" => p.dispatch,
    "REGISTRY" => p.registry,
)


text:
module Primes

function randprime(bits)
    for i in 1:100_000
        x=rand(BigInt(1):BigInt(2)^bits-1)
        if isprime(x)
            return (true, x, i)
        end
    end
    return (false, 0, 0)
end

end

max_stars_repo_name: drichardson/examples

max_issues_repo_name: drichardson/examples

max_forks_repo_name: drichardson/examples


text:
abstract type AbstractDensity <: AbstractMCMC.AbstractModel end
abstract type AbstractApproxPosterior <: AbstractDensity end
#=
unconditional_sample(rng::AbstractRNG,density::AbstractDensity) = error("define a method to cover unconditional_sample(rng::AbstractRNG,density::"*repr(typeof(density))*").")
loglike(density::AbstractDensity,sample) = error("define a method to cover logpdf(density::"*repr(typeof(density))*",sample). must return a named tuple (logprior = ?, loglikelihood = ?)")
length(density::AbstractDensity) = error("define a method to cover length(density::"*repr(typeof(density))*").")
accept(density::AbstractDensity,rng::AbstractRNG,old_ld,new_ld) = error("define a method to cover accept(density::"*repr(typeof(density))*",rng::AbstractRNG,old_ld,new_ld). must return boolean to accept or reject a transition from old_ld → new_ld")
=#

struct Particle{Xt}
    x::Xt
    Particle(x::T) where {T} = new{T}(x)
end

op(f, a::Particle, b::Particle) = Particle(op(f, a.x, b.x))
op(f, a::Particle, b::Number) = Particle(op(f, a.x, b))
op(f, a::Number, b::Particle) = Particle(op(f, a, b.x))
op(f, a::Number, b::Number) = f(a, b)
op(f, a, b) = op.(Ref(f), a, b)

op(f, a::Particle) = Particle(op(f, a.x))
op(f, a::Number) = f(a)
op(f, a) = op.(Ref(f), a)

op(f, args...) = foldl((x, y) -> op(f, x, y), args)

push_p(density::AbstractDensity, p::Particle) = p
push_p(density::AbstractApproxPosterior, p::Particle) = Particle(push_p(density.prior, p.x))
push_p(density::Factored, p) = push_p.(density.p, p)
push_p(density::Distribution, p) = push_p.(Ref(density), p)
push_p(density::ContinuousDistribution, p::Number) = float(p)
push_p(density::DiscreteDistribution, p::Number) = round(Int, p)

unconditional_sample(rng::AbstractRNG, density::AbstractApproxPosterior) =
    Particle(rand(rng, density.prior))

length(density::AbstractApproxPosterior) = length(density.prior)

struct ApproxKernelizedPosterior{P<:Distribution,C,S<:Real} <: AbstractApproxPosterior
    prior::P
    cost::C
    scale::S
    ApproxKernelizedPosterior(
        prior::T1,
        cost::T2,
        target_average_cost::T3,
    ) where {T1,T2,T3} = new{T1,T2,T3}(prior, cost, target_average_cost)
end

function loglike(density::ApproxKernelizedPosterior, sample::Particle)
    lp = logpdf(density.prior, sample.x)
    ll = lp
    if isfinite(lp)
        ll = -0.5 * abs2(density.cost(sample.x) / density.scale)
    end
    (logprior = lp, loglikelihood = ll)
end

is_valid_logdensity(density::ApproxKernelizedPosterior, ld) = isfinite(sum(ld))

function accept(
    density::ApproxKernelizedPosterior,
    rng::AbstractRNG,
    old_ld,
    new_ld,
    ld_correction,
)
    isfinite(ld_correction) || error("ld_correction is invalid")
    is_valid_logdensity(density, old_ld) || error("starting sample invalid.")
    is_valid_logdensity(density, new_ld) || return false

    lW = ld_correction + sum(new_ld) - sum(old_ld)
    return -randexp(rng) <= lW
end
struct ApproxPosterior{P<:Distribution,C,S<:Real} <: AbstractApproxPosterior
    prior::P
    cost::C
    maxcost::S
    ApproxPosterior(prior::T1, cost::T2, max_cost::T3) where {T1,T2,T3} =
        new{T1,T2,T3}(prior, cost, max_cost)
end

function loglike(density::ApproxPosterior, sample::Particle)
    lp = logpdf(density.prior, sample.x)
    cs = -lp
    if isfinite(lp)
        cs = density.cost(sample.x)
    end
    (logprior = lp, cost = cs)
end

is_valid_logdensity(density::ApproxPosterior, ld) =
    isfinite(ld.cost) && isfinite(ld.logprior)

function accept(density::ApproxPosterior, rng::AbstractRNG, old_ld, new_ld, ld_correction)
    isfinite(ld_correction) || error("ld_correction is invalid")
    is_valid_logdensity(density, old_ld) || error("starting sample invalid.")
    is_valid_logdensity(density, new_ld) || return false

    lW = ld_correction + new_ld.logprior - old_ld.logprior
    lW2 = max(density.maxcost, old_ld.cost) - new_ld.cost
    (-randexp(rng) <= lW) && lW2 >= 0
end
struct CommonLogDensity{N,A,B} <: AbstractDensity
    sample_init::A
    lπ::B
    CommonLogDensity(nparameters::Int, sample_init::A, lπ::B) where {A,B} =
        new{nparameters,A,B}(sample_init, lπ)
end

unconditional_sample(rng::AbstractRNG, density::CommonLogDensity) =
    Particle(density.sample_init(rng))

length(density::CommonLogDensity{N}) where {N} = N

function loglike(density::CommonLogDensity, sample::Particle)
    density.lπ(sample.x)
end

is_valid_logdensity(density::CommonLogDensity, ld) = isfinite(ld)

function accept(density::CommonLogDensity, rng::AbstractRNG, old_ld, new_ld, ld_correction)
    isfinite(ld_correction) || error("ld_correction is invalid")
    is_valid_logdensity(density, old_ld) || error("starting sample invalid.")
    is_valid_logdensity(density, new_ld) || return false
    return -randexp(rng) <= ld_correction + new_ld - old_ld
end

"""
    ApproxKernelizedPosterior(
        prior::Distribution,
        cost::Function,
        target_average_cost::Real
    )
this function will return a type which can be used in the `sample` function as an ABC density,
this type works by assuming Gaussianly distributed errors 𝒩(0,ϵ), ϵ is specified in the variable `target_average_cost`.
"""
ApproxKernelizedPosterior
"""
    ApproxPosterior(
        prior::Distribution,
        cost::Function,
        max_cost::Real
    )
this function will return a type which can be used in the `sample` function as an ABC density,
this type works by assuming uniformly distributed errors in [-ϵ,ϵ], ϵ is specified in the variable `max_cost`.
"""
ApproxPosterior

"""
    CommonLogDensity(nparameters, sample_init, lπ)
this function will return a type for performing classical MCMC via the `sample` function.

`nparameters`: total number of parameters per sample.

`sample_init`: function which accepts an `RNG::AbstractRNG` and returns a sample for `lπ`.

`lπ`: function which accepts a sample, and returns a log-density float value.
"""
CommonLogDensity

export CommonLogDensity, ApproxPosterior, ApproxKernelizedPosterior


text:
module Minijyro

using Distributions

export sample!,
    MinijyroModel,
    AbstractHandler,
    TraceHandler,
    LogJointHandler,
    ConditionHandler,
    ReplayHandler,
    EscapeHandler,
    EscapeException,
    queue,
    @jyro,
    handle!,
    handle,
    discrete_enumeration,
    hmc,
    nuts,
    enter!,
    exit!

include("minijyro_types.jl")
include("dsl.jl")
include("handlers.jl")
include("inference.jl")

"""
    sample!(args...)

Allows sampling from probability distributions with side effects specified by
effect handlers in `handlers_stack`.

# Arguments
- `trace::Dict`: dict that can be used to store things in a global state
- `handlers_stack::Array{AbstractHandler,1}`: stack of effect handlers
- `name::Any`: unique name for the sample site
- `dist::Distributions.Distribution`: distribution to sample from
"""
function sample!(
    trace::Dict,
    handlers_stack::Array{AbstractHandler,1},
    name::Any,
    dist::Distributions.Distribution
)
    if length(handlers_stack) == 0
        return rand(dist)
    end

    # NOTE: In this simplified case we could replace :args with :dist because
    #       our "effect" is always sampling from a distribution.
    # NOTE: This could be a custom type but is a Dict for simplicity.
    initial_msg = Dict(
        :fn => rand,
        :args => (dist, ),
        :name => name,
        :value => nothing,
        :is_observed => false,
        :done => false, # Use if we still want other handlers to be active.
        :stop => false, # Use if we want handlers below on the stack not to be active.
        :continuation => nothing # Function to be called after effect handlers have been applied.
    )
    msg = apply_stack!(trace, handlers_stack, initial_msg)
    return msg[:value]
end

"""
    apply_stack!(trace::Dict, handlers_stack::Array{AbstractHandler,1}, msg::Dict)

Apply all the effect handlers from `handlers_stack` to the effect described in
`msg`.

# Arguments
- `trace::Dict`: dict which can be used to store things in a global state
- `handlers_stack::Array{AbstractHandler,1}`: stack of effect handlers
- `msg::Dict`: dict describing the effectful computation

# Returns
- `Dict`: message after effect handlers are applied
"""
function apply_stack!(
    trace::Dict,
    handlers_stack::Array{AbstractHandler,1},
    msg::Dict
)
    @assert length(handlers_stack) > 0

    handler_counter = 0
    # Loop through handlers from top of the stack to the bottom.
    for handler in handlers_stack[end:-1:1]
        handler_counter += 1
        process_message!(trace, handler, msg)
        if msg[:stop]
            break
        end
    end

    if !(msg[:value] != nothing || msg[:done])
        msg[:value] = msg[:fn](msg[:args]...)
    end

    # Loop through handlers from bottom of the stack to the top.
    # If we exited the first loop early then we will start looping from the
    # handler which caused the loop to exit.
    for handler in handlers_stack[end-handler_counter+1:end]
        postprocess_message!(trace, handler, msg)
    end

    if msg[:continuation] != nothing
        msg[:continuation](trace, msg)
    end

    return msg
end
end # module


text:
function IterativeLinearSolver(arg0::IterationManager)
    return IterativeLinearSolver((IterationManager,), arg0)
end

function IterativeLinearSolver(arg0::jint)
    return IterativeLinearSolver((jint,), arg0)
end

function solve(obj::IterativeLinearSolver, arg0::RealLinearOperator, arg1::RealVector)
    return jcall(obj, "solve", RealVector, (RealLinearOperator, RealVector), arg0, arg1)
end

function solve(obj::IterativeLinearSolver, arg0::RealLinearOperator, arg1::RealVector, arg2::RealVector)
    return jcall(obj, "solve", RealVector, (RealLinearOperator, RealVector, RealVector), arg0, arg1, arg2)
end

function solve_in_place(obj::IterativeLinearSolver, arg0::RealLinearOperator, arg1::RealVector, arg2::RealVector)
    return jcall(obj, "solveInPlace", RealVector, (RealLinearOperator, RealVector, RealVector), arg0, arg1, arg2)
end


text:
using DelayDiffEq, DiffEqProblemLibrary.DDEProblemLibrary
using LinearAlgebra, Test

DDEProblemLibrary.importddeproblems()
const PROB_WALTMAN = DDEProblemLibrary.prob_dde_RADAR5_waltman_5
const PROB_KWARGS = (reltol = 1e-6, abstol = [1e-21, 1e-21, 1e-21, 1e-21, 1e-9, 1e-9])

# solution at final time point T = 300 obtained from RADAR5
# with relative tolerance 1e-6 and componentwise absolute tolerances
# 1e-21, 1e-21, 1e-21, 1e-21, 1e-9, and 1e-9
const RADAR5_SOL = [6.154488183e-16, 3.377120916e-7, 4.22140331e-7,
                    2.142554562e-6, 299.9999999, 299.6430338]

function test_waltman_sol(sol)
    @test sol.retcode === :Success
    @test sol.t[end] == 300

    # compare solution at the final time point with RADAR5
    for i in 1:6
        @test sol.u[end][i] ≈ RADAR5_SOL[i] rtol = 1e-3 atol = 1e-17
    end
end

# standard factorization
sol1 = solve(PROB_WALTMAN, MethodOfSteps(Rosenbrock23()); PROB_KWARGS...)
test_waltman_sol(sol1)

# in-place LU factorization
sol2 = solve(PROB_WALTMAN, MethodOfSteps(Rosenbrock23(linsolve=LinSolveFactorize(lu!)));
             PROB_KWARGS...)
test_waltman_sol(sol2)

# out-of-place LU factorization
sol3 = solve(PROB_WALTMAN, MethodOfSteps(Rosenbrock23(linsolve=LinSolveFactorize(lu)));
             PROB_KWARGS...)
test_waltman_sol(sol3)

# compare in-place and out-of-place LU factorization
@test sol2.t == sol3.t
@test sol2.u == sol3.u


text:
@rule typeof(+)(:in1, Marginalisation) (m_out::PointMass, m_in2::PointMass) = PointMass(mean(m_out) - mean(m_in2))

@rule typeof(+)(:in1, Marginalisation) (m_out::UnivariateNormalDistributionsFamily, m_in2::PointMass) = NormalMeanVariance(mean(m_out) - mean(m_in2), var(m_out))
@rule typeof(+)(:in1, Marginalisation) (m_out::PointMass, m_in2::UnivariateNormalDistributionsFamily) = NormalMeanVariance(mean(m_out) - mean(m_in2), var(m_in2))

@rule typeof(+)(:in1, Marginalisation) (m_out::MultivariateNormalDistributionsFamily, m_in2::PointMass) = MvNormalMeanCovariance(mean(m_out) - mean(m_in2), cov(m_out))
@rule typeof(+)(:in1, Marginalisation) (m_out::PointMass, m_in2::MultivariateNormalDistributionsFamily) = MvNormalMeanCovariance(mean(m_out) - mean(m_in2), cov(m_in2))

# TODO: Do we need it?
# Specific case for precision parametrisation
@rule typeof(+)(:in1, Marginalisation) (m_out::NormalMeanPrecision, m_in2::PointMass) = NormalMeanPrecision(mean(m_out) - mean(m_in2), precision(m_out))
@rule typeof(+)(:in1, Marginalisation) (m_out::PointMass, m_in2::NormalMeanPrecision) = NormalMeanPrecision(mean(m_out) - mean(m_in2), precision(m_in2))

@rule typeof(+)(:in1, Marginalisation) (m_out::NormalMeanPrecision, m_in2::NormalMeanPrecision) = begin
    mout, pout = mean_precision(m_out)
    min2, pin2 = mean_precision(m_in2)
    return NormalMeanPrecision(mout - min2, pout * pin2 / (pout + pin2))
end

@rule typeof(+)(:in1, Marginalisation) (m_out::UnivariateNormalDistributionsFamily, m_in2::UnivariateNormalDistributionsFamily) = begin
    mout, vout = mean_var(m_out)
    min2, vin2 = mean_var(m_in2)
    return NormalMeanVariance(mout - min2, vout + vin2)
end

@rule typeof(+)(:in1, Marginalisation) (m_out::MultivariateNormalDistributionsFamily, m_in2::MultivariateNormalDistributionsFamily) = begin
    mout, vout = mean_cov(m_out)
    min2, vin2 = mean_cov(m_in2)
    return MvNormalMeanCovariance(mout - min2, vout + vin2)
end


text:
function zad41(A,key)
    br = 0
    for i=1:length(A)
        if(A[i] == key)
            br +=1
        end
        if(A[i] != key && br !=0)
            break
        end
    end
    return br
end


text:
# # plot(predl1[1:50])
# # plot!(yy_train[1:50])
# # 10,3,4,5,20
# i1=findall(abs.(predl1-yy_train).>0.1)[1]
# # plot(0:9,predl1_train_all[i1,:])
# yy_train_all[i1]
# yy_train[i1]
# predl1[i1]
# using Images
# train_tensor = MNIST.traintensor()
# train_images = MNIST.convert2image(1.0 .-train_tensor)
# train_images[:,:,i1]

using Images
test_tensor = MNIST.testtensor()
test_images = MNIST.convert2image(1.0 .-test_tensor)
i1=findall(abs.(yy_test_pred-yy_test_all).>0.1)
nn=9
test_images[:,:,i1[nn]]
yy_test_all[i1[nn]]
yy_test_pred[i1[nn]]

# x1=Vector{Float64}(undef,0)
# for i=1:28
#     global x1
#     x1=[x1;float.(train_images[i,1:28,1])]
# end
#
# maximum(abs.(xx_train[1,:].-x1))
#
# plot()
# for i=1:28
#     display(plot!(xx_train[1,(i-1)*28+1:i*28]))
# end
#
#
#
#
#
# train_images[:,:,1]
# x1=Array{Float64}(undef,0,14)
# for i=1:2:28
#     global x1
#     x1=[x1;float.(train_images[i,1:2:end,1])']
# end
# Gray.(x1)
#
#
# im1=train_images[:,:,1]
# im1f=float.(im1.+.1)
# img=Gray.(im1f)


text:
import Yota: ungetindex, ungetindex!

@testset "helpers: getindex" begin
    x = rand(4)
    dx = zero(x)
    dy = 1.0
    @test ungetindex!(dx, x, dy, 2) == [0.0, 1, 0, 0]
    @test ungetindex(x, dy, 2) == [0.0, 1, 0, 0]

    x = rand(4)
    dx = zero(x)
    dy = [1, 1, 1]
    @test ungetindex!(dx, x, dy, [1, 3, 1]) == [2.0, 0, 1, 0]
    @test ungetindex(x, dy, [1, 3, 1]) == [2.0, 0, 1, 0]

    x = rand(4, 5)
    dx = zero(x)
    dy = ones(4, 3)
    expected = [2 0 1 0 0;
                2 0 1 0 0;
                2 0 1 0 0;
                2 0 1 0 0.0]
    @test ungetindex!(dx, x, dy, :, [1, 3, 1]) == expected
    @test ungetindex(x, dy, :, [1, 3, 1]) == expected

    ## additional smoke tests

    x = rand(3, 4)
    I = (1:2, 2:3)
    y = x[I...]
    dy = ones(size(y))
    dx = zero(x)
    ungetindex!(dx, x, dy, I...)

    x = rand(3, 4)
    I = (1, 2:3)
    y = x[I...]
    dy = ones(size(y))
    dx = zero(x)
    ungetindex!(dx, x, dy, I...)

    x = rand(3, 4)
    I = (1, :)
    y = x[I...]
    dy = ones(size(y))
    dx = zero(x)
    ungetindex!(dx, x, dy, I...)

    x = rand(3, 4)
    I = (1, [1, 3])
    y = x[I...]
    dy = ones(size(y))
    dx = zero(x)
    ungetindex!(dx, x, dy, I...)

    # single-element Cartesian index
    x = rand(3, 4)
    I = (1, 2)
    y = x[I...]
    dy = 1.0
    dx = zero(x)
    ungetindex!(dx, x, dy, I...)

    # single-element linear index
    x = rand(3, 4)
    I = (5,)
    y = x[I...]
    dy = 1.0
    dx = zero(x)
    ungetindex!(dx, x, dy, I...)

    # CUDA
    if CUDA.functional()
        x = rand(4) |> cu
        dx = zero(x) |> cu
        dy = 1.0f0
        @test ungetindex!(dx, x, dy, 2) == cu([0.0, 1, 0, 0])
        @test ungetindex(x, dy, 2) == cu([0.0, 1, 0, 0])

        x = rand(4, 5) |> cu
        dx = zero(x) |> cu
        dy = ones(4, 3) |> cu
        @test ungetindex!(dx, x, dy, :, [1, 3, 1]) == cu(expected)
        @test ungetindex(x, dy, :, [1, 3, 1]) == cu(expected)
    end

end


text:
using BandedMatrices, LinearAlgebra, ArrayLayouts, FillArrays, Test, Base64
import BandedMatrices: banded_mul!, isbanded, AbstractBandedLayout, BandedStyle,
                        rowsupport, colsupport, _BandedMatrix, BandedColumns
import ArrayLayouts: OnesLayout, UnknownLayout

struct PseudoBandedMatrix{T} <: AbstractMatrix{T}
    data::Array{T}
    l::Int
    u::Int
end

Base.size(A::PseudoBandedMatrix) = size(A.data)
function Base.getindex(A::PseudoBandedMatrix, j::Int, k::Int)
    l, u = bandwidths(A)
    if -l ≤ k-j ≤ u
        A.data[j, k]
    else
        zero(eltype(A.data))
    end
end
function Base.setindex!(A::PseudoBandedMatrix, v, j::Int, k::Int)
    l, u = bandwidths(A)
    if -l ≤ k-j ≤ u
        A.data[j, k] = v
    else
        error("out of band.")
    end
end

struct PseudoBandedLayout <: AbstractBandedLayout end
Base.BroadcastStyle(::Type{<:PseudoBandedMatrix}) = BandedStyle()
BandedMatrices.MemoryLayout(::Type{<:PseudoBandedMatrix}) = PseudoBandedLayout()
BandedMatrices.isbanded(::PseudoBandedMatrix) = true
BandedMatrices.bandwidths(A::PseudoBandedMatrix) = (A.l , A.u)
BandedMatrices.inbands_getindex(A::PseudoBandedMatrix, j::Int, k::Int) = A.data[j, k]
BandedMatrices.inbands_setindex!(A::PseudoBandedMatrix, v, j::Int, k::Int) = setindex!(A.data, v, j, k)
LinearAlgebra.fill!(A::PseudoBandedMatrix, v) = fill!(A.data,v)

@testset "banded matrix interface" begin
    @testset "Zeros" begin
        @test isbanded(Zeros(5,6))
        @test bandwidths(Zeros(5,6)) == (0,0)
        @test BandedMatrices.inbands_getindex(Zeros(5,6), 1,2) == 0
    end

    @testset "Eye" begin
        @test isbanded(Eye(5))
        @test bandwidths(Eye(5)) == (0,0)
        @test BandedMatrices.inbands_getindex(Eye(5), 1,1) == 1

        V = view(Eye(5),2:3,:)
        @test bandwidths(V) == (-1,1)
        @test MemoryLayout(typeof(V)) isa BandedColumns{OnesLayout}
        @test BandedMatrix(V) == [0 1 0 0 0; 0 0 1 0 0]

        B = brand(5,5,1,2)
        @test Eye(5) * B == B
        @test B * Eye(5) == B
        @test muladd!(2.0, Eye(5), B, 0.0, zeros(5,5)) == 2B
        @test muladd!(2.0, B, Eye(5), 0.0, zeros(5,5)) == 2B
    end

    @testset "Diagonal" begin
        A = Diagonal(ones(5,5))
        @test isbanded(A)
        @test bandwidths(A) == (0,0)
        @test BandedMatrices.inbands_getindex(A, 1,1) == 1
        BandedMatrices.inbands_setindex!(A, 2, 1,1)
        @test A[1,1] == 2
        @test A[1,2] == 0
        @test BandedMatrices.@inbands(A[1,2]) == 2

        @test MemoryLayout(view(A, 1:3,2:4)) isa BandedColumns{DenseColumnMajor}
        @test MemoryLayout(view(A, [1,2,3],2:4)) isa UnknownLayout
    end

    A = SymTridiagonal([1,2,3],[4,5])
    @test isbanded(A)
    @test bandwidths(A) == (1,1)
    @test BandedMatrices.inbands_getindex(A, 1,1) == 1
    BandedMatrices.inbands_setindex!(A, 2, 1,1)
    @test A[1,1] == 2

    A = PseudoBandedMatrix(rand(5, 4), 2, 2)
    B = rand(5, 4)
    C = copy(B)

    @test Matrix(B .= 2.0 .* A .+ B) ≈ 2*Matrix(A) + C

    A = PseudoBandedMatrix(rand(5, 4), 1, 2)
    B = (z -> exp(z)-1).(A)
    @test B isa BandedMatrix
    @test bandwidths(B) == bandwidths(A)
    @test B == (z -> exp(z)-1).(Matrix(A))

    A = PseudoBandedMatrix(rand(5, 4), 1, 2)
    B = A .* 2
    @test B isa BandedMatrix
    @test bandwidths(B) == bandwidths(A)
    @test B == 2Matrix(A) == (2 .* A)

    A = PseudoBandedMatrix(rand(5, 4), 1, 2)
    B = PseudoBandedMatrix(rand(5, 4), 2, 1)

    @test A .+ B isa BandedMatrix
    @test bandwidths(A .+ B) == (2,2)
    @test A .+ B == Matrix(A) + Matrix(B)

    A = PseudoBandedMatrix(rand(5, 4), 1, 2)
    B = PseudoBandedMatrix(rand(5, 4), 2, 3)
    C = deepcopy(B)
    @test Matrix(C .= 2.0 .* A .+ C) ≈ 2*Matrix(A) + B ≈ 2*A + Matrix(B) ≈ 2*Matrix(A) + Matrix(B) ≈ 2*A + B

    y = rand(4)
    z = zeros(5)
    muladd!(1.0, A, y, 0.0, z)
    @test z ≈ A*y ≈ Matrix(A)*y

    @test bandwidths(BandedMatrix(A)) ==
            bandwidths(BandedMatrix{Float64}(A)) ==
            bandwidths(BandedMatrix{Float64,Matrix{Float64}}(A)) ==
            bandwidths(convert(BandedMatrix{Float64}, A)) ==
            bandwidths(convert(BandedMatrix{Float64,Matrix{Float64}},A)) ==
            bandwidths(A)
end

@testset "Bi/Tri/Diagonal" begin
    D = Diagonal(randn(5))
    @test bandwidths(D) == (0,0)

    A =  brand(5,5,1,1)
    @test D*A isa BandedMatrix
    @test D*A == D*Matrix(A)

    @test A*D isa BandedMatrix
    @test A*D == Matrix(A)*D

    M = MulAdd(2.0, A, D, 1.0, A)
    @test copyto!(similar(M), M) ≈ 2A*D + A
    M = MulAdd(2.0, D, A, 1.0, A)
    @test copyto!(similar(M), M) ≈ 2D*A + A

    D = Diagonal(rand(Int,10))
    B = brand(10,10,-1,2)
    J = SymTridiagonal(randn(10), randn(9))
    T = Tridiagonal(randn(9), randn(10), randn(9))
    Bl = Bidiagonal(randn(10), randn(9), :L)
    Bu = Bidiagonal(randn(10), randn(9), :U)

    @test bandwidths(D) == (0,0)
    @test bandwidths(J) == bandwidths(T) == (1,1)
    @test bandwidths(Bu) == (0,1)
    @test bandwidths(Bl) == (1,0)

    @test isbanded(D)
    @test isbanded(J)
    @test isbanded(T)
    @test isbanded(Bu)
    @test isbanded(Bl)

    @test BandedMatrix(J) == J
    @test BandedMatrix(T) == T
    @test BandedMatrix(Bu) == Bu
    @test BandedMatrix(Bl) == Bl

    @test Base.BroadcastStyle(Base.BroadcastStyle(typeof(B)), Base.BroadcastStyle(typeof(D))) ==
        Base.BroadcastStyle(Base.BroadcastStyle(typeof(D)), Base.BroadcastStyle(typeof(B))) ==
        Base.BroadcastStyle(Base.BroadcastStyle(typeof(B)), Base.BroadcastStyle(typeof(T))) ==
        Base.BroadcastStyle(Base.BroadcastStyle(typeof(T)), Base.BroadcastStyle(typeof(B))) ==
        Base.BroadcastStyle(Base.BroadcastStyle(typeof(B)), Base.BroadcastStyle(typeof(J))) ==
        Base.BroadcastStyle(Base.BroadcastStyle(typeof(J)), Base.BroadcastStyle(typeof(B))) ==
        Base.BroadcastStyle(Base.BroadcastStyle(typeof(Bl)), Base.BroadcastStyle(typeof(B))) ==
        Base.BroadcastStyle(Base.BroadcastStyle(typeof(Bu)), Base.BroadcastStyle(typeof(B))) ==
            BandedStyle()

    A = B .+ D
    Ã = D .+ B
    @test A isa BandedMatrix
    @test Ã isa BandedMatrix
    @test bandwidths(A) == bandwidths(Ã) == (0,2)
    @test Ã == A == Matrix(B) + Matrix(D)

    A = B .+ J
    Ã = J .+ B
    @test A isa BandedMatrix
    @test Ã isa BandedMatrix
    @test bandwidths(A) == bandwidths(Ã) == (1,2)
    @test Ã == A == Matrix(B) + Matrix(J)

    A = B .+ T
    Ã = T .+ B
    @test A isa BandedMatrix
    @test Ã isa BandedMatrix
    @test bandwidths(A) == bandwidths(Ã) == (1,2)
    @test Ã == A == Matrix(B) + Matrix(T)

    A = B .+ Bl
    Ã = Bl .+ B
    @test A isa BandedMatrix
    @test Ã isa BandedMatrix
    @test bandwidths(A) == bandwidths(Ã) == (1,2)
    @test Ã == A == Matrix(B) + Matrix(Bl)

    A = B .+ Bu
    Ã = Bu .+ B
    @test A isa BandedMatrix
    @test Ã isa BandedMatrix
    @test bandwidths(A) == bandwidths(Ã) == (0,2)
    @test Ã == A == Matrix(B) + Matrix(Bu)

    @test layout_getindex(D,1:10,1:10) isa BandedMatrix
    @test layout_getindex(B,1:10,1:10) isa BandedMatrix
    @test layout_getindex(Bu,1:10,1:10) isa BandedMatrix
    @test layout_getindex(J,1:10,1:10) isa BandedMatrix
    @test layout_getindex(T,1:10,1:10) isa BandedMatrix
end

@testset "rot180" begin
    A = brand(5,5,1,2)
    R = rot180(A)
    @test bandwidths(R) == (2,1)
    @test R == rot180(Matrix(A))

    A = brand(5,4,1,2)
    R = rot180(A)
    @test bandwidths(R) == (3,0)
    @test R == rot180(Matrix(A))

    A = brand(5,6,1,-1)
    R = rot180(A)
    @test bandwidths(R) == (-2,2)
    @test R == rot180(Matrix(A))

    A = brand(6,5,-1,1)
    R = rot180(A)
    @test bandwidths(R) == (2,-2)
    @test R == rot180(Matrix(A))
end


text:
function msa(network::AbstractNetwork,
              trips::AbstractMatrix{T},
              costfn::CostFunction;
              basedon=:link,
              errtol=1e-4) where {T<:Real, U<:Integer}
    # initialize
    flows = allornothing(network, trips, costfn; basedon=basedon)

    # start iteration
    err = 1.
    iterno = 2 # first iteration @ initialization
    while err > errtol
        ## find target solution
        linkcosts = costfn(flows).costs
        shortflows = allornothing(network, trips, linkcosts; basedon=basedon)

        ## find stepsize
        μ = 1/iterno

        ## calculate new link travel times
        flows += μ * (shortflows - flows)

        ## calculate error
        err = sum(flows .* linkcosts)/sum(shortflows .* linkcosts) - 1.
        iterno += 1
    end
    return (flows, err)
end


text:
using Plots
using PrettyTables

include("rfiod.jl")
using Main.RFIOD

io = stdout
# Uncomment below to save to file instead
# io = open("output.txt", "w")

F₁(x) = 2x - 3sin(x) + 5
F₁′(x) = 2 - 3cos(x)
F₁root = -2.8832368725582835

F₂(x) = x^3 - 8.5x^2 + 20x - 8
F₂′(x) = 3x^2 - 17x + 20
F₂root = 4

# -- Question 1 -- #

q1plot = plot(F₁, -5:0.1:5, framestyle=:origin, legend = false, title = "F₁(x) = 2x - 3sin(x) + 5 has a single real root in [-5,5]")
# F₁′(x) > 0 outside of this interval with F₁(-5) < 0, F₁(5) > 0, so no roots outside

# -- Binary Search Programming Task -- #
println(io, "-- Binary Search Programming Task --")

# Initial interval [-5, 5]
BSresult1 = binarysearch(F₁, -5, 5)
# Initial interval [-5, 5], but with inputs swapped
BSresult2 = binarysearch(F₁, 5, -5)
# Initial interval [-2.88323688, -2.88323687]
BSresult3 = binarysearch(F₁, -2.88323688, -2.88323687)
# Initial interval [-100, 100]
BSresult4 = binarysearch(F₁, -100, 100)

println(io, "Binary search with initial interval [a,b]")
# Table created for display
tabledata = hcat(
    [-5, 5, -2.88323688, -100],
    [5, -5, -2.88323687, 100],
    Plots.unzip([BSresult1, BSresult2, BSresult3, BSresult4])...
)
pretty_table(io, tabledata, header = ["a", "b", "root", "iterations"])

# Initial interval [-1, 1], which does not contain a root
BSresult5 = try binarysearch(F₁, -1, 1) catch ErrorException "Error, root not found" end
println(io, "Binary search with the invalid initial interval [-1, 1] gives:")
show(io, BSresult5)

# -- Question 3 -- #
println(io, "\n\n-- Question 3 --")

f(k) = x -> x - F₁(x)/(2+k)

"""
    q3FPI(k, x₀, Nₘₐₓ, [ε])

Run the fixed point iteration algorithm used in Question 3, with parameter `k` and starting value `x₀`.

Terminate if consecutive iterations differ by less than `ε`, or after `Nₘₐₓ` iterations.

Outputs a table to the REPL with the results, and returns a staircase/cobweb diagram.

# Examples

```julia-repl
julia> q3FPI(3,-2,10)
Fixed point iteration with k = 3, Nₘₐₓ = 10
┌─────┬──────────┬────────────┬───────────┐
│   n │       xₙ │         ϵₙ │   ϵₙ/ϵₙ₋₁ │
├─────┼──────────┼────────────┼───────────┤
│ 0.0 │     -2.0 │   0.883237 │       0.0 │
│ 1.0 │ -2.74558 │   0.137658 │  0.155857 │
│ 2.0 │ -2.87879 │ 0.00444334 │  0.032278 │
│ 3.0 │ -2.88315 │ 9.00029e-5 │ 0.0202557 │
│ 4.0 │ -2.88324 │ 1.79286e-6 │ 0.0199201 │
│ 5.0 │ -2.88324 │ 3.57019e-8 │ 0.0199133 │
└─────┴──────────┴────────────┴───────────┘
```

"""
function q3FPI(k, x₀, Nₘₐₓ, ε = 1e-5)

    result = fixedpointiteration(f(k), x₀, Nₘₐₓ, ε)
    n = length(result) - 1

    errors = result .- F₁root
    errorratios = [0, (errors[2:(n+1)] ./ errors[1:n])...]
    tabledata = hcat(
        collect(0:n),
        result,
        errors,
        errorratios
    )
    println(io, "Fixed point iteration with k = $k, Nₘₐₓ = $Nₘₐₓ")
    pretty_table(io, tabledata, header = ["n", "xₙ", "ϵₙ", "ϵₙ/ϵₙ₋₁"])

    FPIplot = staircasediagram(f(k), result)
    plot!(FPIplot, title = "k = $k")
    return FPIplot

end

# k = 0, ε = 10⁻⁵, x₀ = -2, Nₘₐₓ = 10
q3plot1 = q3FPI(0, -2, 10)

# k = 4 gives monotonic convergence as f(8)'(x*) > 0
q3plot2 = q3FPI(4, -2, 20)

# k = 2 gives oscillatory convergence as f(1)'(x*) < 0
q3plot3 = q3FPI(2, -2, 20)

# k = 16, Nₘₐₓ = 50
q3FPI(16, -2, 50)

# -- Question 4 -- #
println(io, "\n-- Question 4 --")

g(x) = x - F₂(x)/20

q4result = fixedpointiteration(g, 5, 1000, 1e-5)
# n and xₙ take the final values given below
n = length(q4result) - 1
xₙ = q4result[end]
println(io, "Fixed point iteration to find the double root of
F₂(x) = x³ - 8.5x² + 20x - 8 results in n = $n iterations
with estimated root xₙ = $xₙ")

# -- Question 5 -- #
println(io, "\n-- Question 5 --")

"""
    q5NR(f, x₀, Nₘₐₓ, root, order)

Run the Newton-Raphson algorithm used in Question 5 on the function `f(x) = x - F(x)/F'(x)` starting at `x₀`

Terminate the fixed point iteration after at most Nₘₐₓ iterations.

Outputs a table with errors to `root` and error ratios of order `order`, and returns the vector of iterations.

# Examples

```julia-repl
julia> q5NR(x -> x - sin(x)/cos(x), 3, 10, π, 2)
┌─────┬─────────┬──────────────┬──────────────┐
│   n │      xₙ │           ϵₙ │     ϵₙ/ϵₙ₋₁² │
├─────┼─────────┼──────────────┼──────────────┤
│ 0.0 │     3.0 │    -0.141593 │          0.0 │
│ 1.0 │ 3.14255 │  0.000953889 │    0.0475791 │
│ 2.0 │ 3.14159 │ -2.89316e-10 │ -0.000317963 │
│ 3.0 │ 3.14159 │          0.0 │          0.0 │
│ 4.0 │ 3.14159 │          0.0 │          NaN │
└─────┴─────────┴──────────────┴──────────────┘
5-element Vector{Float64}:
 3.0
 3.142546543074278
 3.141592653300477
 3.141592653589793
 3.141592653589793
```
"""
function q5NR(f, x₀, Nₘₐₓ, root, order)
    result = fixedpointiteration(f, x₀, Nₘₐₓ)
    n = length(result) - 1
    errors = result .- root
    errorratios = [0, (errors[2:(n+1)] ./ errors[1:n].^order)...]
    tabledata = hcat(
        collect(0:n),
        result,
        errors,
        errorratios
    )
    pretty_table(io, tabledata, header = ["n", "xₙ", "ϵₙ", "ϵₙ/ϵₙ₋₁" * (order == 2 ? "²" : "")])
    return result
end

"""
    q5NRgraph(F, result, nlines)
Create a graph from the first `nvalues` iterations of `result`, coming from the Newton-Raphson algorithm on `F`.

# Examples

```julia-repl
julia> q5NRgraph(sin, q5NR(x -> x - sin(x)/cos(x), 3, 10, π, 2), 4)
```
"""
function q5NRgraph(F, result, nlines)
    NRplot = scatter((result[1],F(result[1])), markercolor = :black, framestyle = :origin,
        title = "Newton-Raphson with x₀ = $(result[1])", legend = false)
    for n = 1:nlines
        plot!(NRplot, [result[n],result[n+1],result[n+1]],
            [F(result[n]), 0, F(result[n+1])], linecolor = :black)
    end
    plot!(NRplot, F, linecolor = :orange)
    return NRplot
end

f₁(x) = x - F₁(x)/F₁′(x)
f₂(x) = x - F₂(x)/F₂′(x)

# x₀ = 0 for F₁ gives divergence
println(io, "Newton-Raphson method for \nF₁(x) = 2x - 3sin(x) + 5 with x₀ = 0")
q5result1 = q5NR(f₁, 0, 10, F₁root, 2)
q5plot1 = q5NRgraph(F₁, q5result1, 4)

# x₀ = -2 for F₁ gives convergence
println(io, "Newton-Raphson method for \nF₁(x) = 2x - 3sin(x) + 5 with x₀ = -2")
q5result2 = q5NR(f₁, -2, 10, F₁root, 2)
q5plot2 = q5NRgraph(F₁, q5result2, 2)

# x₀ = 5 for F₂ gives convergence
println(io, "Newton-Raphson method for \nF₂(x) = x³ - 8.5x² + 20x - 8 with x₀ = 5")
q5result3 = q5NR(f₂, 5, 50, F₂root, 1)
q5plot3 = q5NRgraph(F₂, q5result3, 10)

# Saves images automatically, uncomment to use
#=

# File extension e.g. ".png", ".pdf"
ext = ".pdf"

# Makes a new directory for the images (or ignores the error if one exists already)
try
    mkdir("images")
catch IOError
    nothing
end

# Saves the figures created above
try
    savefig(q1plot,"images\\q1-singleroot$ext")
    savefig(q3plot1,"images\\q3-k0-graph$ext")
    savefig(q3plot2,"images\\q3-k4-graph$ext")
    savefig(q3plot3,"images\\q3-k2-graph$ext")
    savefig(q5plot1,"images\\q5-f1divergence-graph$ext")
    savefig(q5plot2,"images\\q5-f1convergence-graph$ext")
    savefig(q5plot3,"images\\q5-f2convergence-graph$ext")
catch IOError
    println("Failed to save images")
end

=#

io == stdout || close(io)
nothing


text:
module JungleHelperCassetteCustomPreviewMusic

using ..Ahorn, Maple

@pardef CassetteCustomPreviewMusic(x1::Integer, y1::Integer, x2::Integer=x1, y2::Integer=y1, musicEvent::String="event:/game/general/cassette_preview", musicParamName::String="remix", musicParamValue::Number=1.0) =
    Entity("JungleHelper/CassetteCustomPreviewMusic", x=x1, y=y1, nodes=Tuple{Int, Int}[(0, 0), (x2, y2)], musicEvent=musicEvent, musicParamName=musicParamName, musicParamValue=musicParamValue)

const placements = Ahorn.PlacementDict(
    "Cassette (Custom Preview Music) (Jungle Helper)" => Ahorn.EntityPlacement(
        CassetteCustomPreviewMusic,
        "point",
        Dict{String, Any}(),
        function(entity)
            entity.data["nodes"] = [
                (Int(entity.data["x"]) + 32, Int(entity.data["y"])),
                (Int(entity.data["x"]) + 64, Int(entity.data["y"]))
            ]
        end
    ),
)

Ahorn.nodeLimits(entity::CassetteCustomPreviewMusic) = 2, 2

sprite = "collectables/cassette/idle00.png"

function Ahorn.selection(entity::CassetteCustomPreviewMusic)
    x, y = Ahorn.position(entity)
    controllX, controllY = Int.(entity.data["nodes"][1])
    endX, endY = Int.(entity.data["nodes"][2])

    return [
        Ahorn.getSpriteRectangle(sprite, x, y),
        Ahorn.getSpriteRectangle(sprite, controllX, controllY),
        Ahorn.getSpriteRectangle(sprite, endX, endY)
    ]
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::CassetteCustomPreviewMusic)
    px, py = Ahorn.position(entity)
    nodes = entity.data["nodes"]

    for node in nodes
        nx, ny = Int.(node)

        Ahorn.drawArrow(ctx, px, py, nx, ny, Ahorn.colors.selection_selected_fc, headLength=6)
        Ahorn.drawSprite(ctx, sprite, nx, ny)
        px, py = nx, ny
    end
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::CassetteCustomPreviewMusic, room::Maple.Room) = Ahorn.drawSprite(ctx, sprite, 0, 0)

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
using Dates
using ArgParse
using Printf
using Random
using GaussianProcesses

# before_all = Dates.now()

include("../../btg.jl")
path = "../../datasets/"
include(path * "load_bitcoin.jl")

data, label = bitcoin(path)

max_label = maximum(label)
label = label ./ max_label

xs = copy(data)
xs = xs/maximum(data)

###
###  Manually enter parsed args entries so we can run this file as a script
###
##
parsed_args = Dict()
#push!(parsed_args, "randseed" => 1234)
push!(parsed_args, "randseed" => 120)
push!(parsed_args, "lmax" => 0.3)
push!(parsed_args, "lmin" => 0.2)
push!(parsed_args, "percent_train" => 0.1)
push!(parsed_args, "p" => 1)

randseed = parsed_args["randseed"]; rng = MersenneTwister(randseed)
Random.seed!(randseed);

function covariate_fun(x, p)
    n = size(x, 1)
    d = size(x, 2)
    if p == 1
        return ones(n, 1)
    elseif p == 1 + d
        return hcat(ones(n), x)
    else
        throw(ArgumentError("Only support constant or linear convariate."))
    end
end

function zero_fun(x, p)
    n = size(x, 1)
    d = size(x, 2)
    if p == 1
        return zeros(n, 1)
    elseif p == 1 + d
        return hcat(zeros(n), x)
    else
        throw(ArgumentError("Only support constant or linear convariate."))
    end
end

# training set
p_train = parsed_args["percent_train"]
n_train = (Int64)(round(p_train * length(data)))
n_test = length(data) - n_train
"""
Input percent train
"""
function sample(inds, percent)
    num = (Int64)(round(percent * length(data)))
    #train = sort(rand(1:490, (num,)))
    ind_train = sort(randperm(length(data))[1:num])
    ind_test = copy(inds)
    deleteat!(ind_test, ind_train)
    #work with indices up until this point
    x_train = xs[ind_train]
    x_test = xs[ind_test]
    return ind_train, x_train, ind_test, x_test
end

ind_train, x_train, ind_test, x_test  = sample(data, p_train) #training and testing pts
train_label = label[ind_train]
test_label = label[ind_test]
y_train = train_label

x = reshape(x_train, n_train, 1)
#Fx = covariate_fun(x, parsed_args["p"])
Fx = zero_fun(x, parsed_args["p"])
trainingData0 = trainingData(x, Fx, train_label)

d = getDimension(trainingData0); n = getNumPts(trainingData0); p = getCovDimension(trainingData0)

x_test = reshape(x_test, n_test, 1)
#Fx_test = covariate_fun(x_test, parsed_args["p"])
Fx_test = zero_fun(x_test, parsed_args["p"])
y_test = test_label

elapsedmin = 0

#parameter setting
# myquadtype = parsed_args["sparse"] ? ["SparseCarlo", "SparseCarlo"] : ["QuasiMonteCarlo", "QuasiMonteCarlo"]
myquadtype = ["Gaussian", "Gaussian"]
#rangeλ = [-3.0 3]
rangeλ = reshape([-0.2], 1, 1)
#rangeλ = reshape([-0.3 -0.1], 1, 2)
lmin = parsed_args["lmin"]
lmax = parsed_args["lmax"]
#rangeθ = [1/lmax^2 1/lmin^2]
#rangeθ = reshape([300.0 450.0], 1, 2)
rangeθ = reshape([350.0], 1, 1)

#@info "rangeθ, rangel:" rangeθ, [lmin lmax]
# rangeθ = [0.111 25]
# build btg model

function load_bitcoin_btg()
    return btg(trainingData0, rangeθ, rangeλ; priorθ = inverseUniform(rangeθ), quadtype = myquadtype)
end

function getseed()
    return randseed
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
#######################################################################
# Transform system
#######################################################################

function transform!(sys::System, f::Function; sublats = eachindex(sys.lattice.sublats))
    for s in sublats
        sind = sublatindex(sys.sysinfo, s)
        transform!(sys.lattice.sublats[sind], f)
    end
    sys.lattice.bravais = transform(sys.lattice.bravais, f)
    return sys
end

transform(sys::System, f; kw...) = transform!(deepcopy(sys), f; kw...)

transform!(f::Function; kw...) = x -> transform!(x, f; kw...)
transform(f::Function; kw...) = x -> transform(x, f; kw...)

#######################################################################
# System's Hamiltonian and velocity
#######################################################################

function hamiltonian(sys::System{E,L,T};
                     k = zero(SVector{E,T}), ϕn = blochphases(k, sys),
                     avoidcopy = false) where {E,L,T}
    L == 0 || insertblochphases!(sys.hamiltonian, SVector{L,T}(ϕn))
    return avoidcopy ? sys.hamiltonian.matrix : copy(sys.hamiltonian.matrix)
end

hamiltonian(; kw...) = sys -> hamiltonian(sys; kw...)

velocity(sys::System{E,L}; kw...) where {E,L} = ntuple(i -> velocity(sys, i; kw...), Val(L))
velocity(sys::System{E,L}, axis::Int; kw...) where {E,L} =
    velocity(sys, SMatrix{L,L,Int}(I)[:, axis]; kw...)

function velocity(sys::System{E,L,T}, dϕaxis::Union{Tuple,SVector};
                  k = zero(SVector{E,T}), ϕn = blochphases(k, sys),
                  avoidcopy = false) where {E,L,T}
    L == 0 || insertblochphases!(sys.velocity, SVector{L,T}(ϕn), SVector{L,T}(dϕaxis))
    return avoidcopy ? sys.velocity.matrix : copy(sys.velocity.matrix)
end

function blochphases(k, sys::System{E,L,T}) where {E,L,T}
	length(k) == E || throw(DimensionMismatch(
            "The dimension of the Bloch vector `k` should match the embedding dimension $E"))
	return transpose(bravaismatrix(sys)) * SVector{E,T}(k) / (2pi)
end

#######################################################################
# Operator builder tools
#######################################################################

struct IJV{Tv}
    I::Vector{Int}
    J::Vector{Int}
    V::Vector{Tv}
end

IJV{Tv}() where {Tv} = IJV(Int[], Int[], Tv[])

Base.resize!(ijv::IJV, n) = (resize!(ijv.I, n); resize!(ijv.J, n); resize!(ijv.V, n))
Base.append!(ijv::IJV, ijv2::IJV) =
    (append!(ijv.I, ijv2.I); append!(ijv.J, ijv2.J); append!(ijv.V, ijv2.V);)

function Base.push!(ijv::IJV, s::SMatrix{N,M}, (rowoffset, coloffset)) where {N,M}
    for j in 1:M, i in 1:N
        push!(ijv.I, rowoffset + i)
        push!(ijv.J, coloffset + j)
        push!(ijv.V, s[i,j])
    end
    return nothing
end

Base.isempty(ijv::IJV) = isempty(ijv.I)

struct IJVN{Tv,L}
    IJV::IJV{Tv}
    ndist::SVector{L,Int}
end

IJVN(b::Block) = IJVN(IJV(findnz(b.matrix)...), b.ndist)

struct IJVbuilder{Tv,T,E,L,S<:SystemInfo,EL}
    sysinfo::S
    lattice::Lattice{E,L,T,EL}
    IJVNs::Vector{IJVN{Tv,L}}
    kdtrees::Vector{KDTree{SVector{E,T},Euclidean,T}}
end

IJVbuilder(lat::Lattice{E,L,T}, sysinfo::SystemInfo{Tv}, ) where {E,L,T,Tv} =
    IJVbuilder(sysinfo, lat, IJVN{Tv,L}[],
        Vector{KDTree{SVector{E,T},Euclidean,T}}(undef, length(lat.sublats)))

function Base.push!(builder::IJVbuilder, ijvn::IJVN)
    n = findfirst(bijvn -> bijvn.ndist == ijvn.ndist, builder.IJVNs)
    if n === nothing
        push!(builder.IJVNs, ijvn)
    else
        append!(builder.IJVNs[n].IJV, ijvn.IJV)
    end
    return builder
end

# This violates DRY with the above, but saves some time and allocations (the above is more critical)
function get_or_add_IJVN(builder::IJVbuilder{Tv}, ndist) where {Tv}
    n = findfirst(ijv -> ijv.ndist == ndist, builder.IJVNs)
    if n === nothing
        ijvn = IJVN(IJV{Tv}(), ndist)
        push!(builder.IJVNs, ijvn)
        return ijvn
    else
        return builder.IJVNs[n]
    end
end

Block(ijvn::IJVN, dimh, sysinfo::SystemInfo) =
    Block(ijvn.ndist, sparse(ijvn.IJV.I, ijvn.IJV.J, ijvn.IJV.V, dimh, dimh), sysinfo)

#######################################################################
# Operator builder
#######################################################################

Operator(lat, sysinfo) = Operator(IJVbuilder(lat, sysinfo))

function Operator(builder::IJVbuilder)
    applyterms!(builder, builder.sysinfo.sampledterms...)
    matrix, intra, inters, boundary = assembleblocks(builder)
    ishermitian(matrix) || @warn "Non-hermitian Hamiltonian!"
    return Operator(matrix, intra, inters, boundary)
end

applyterms!(builder, term, terms...) = applyterms!(applyterm!(builder, term), terms...)
applyterms!(builder) = builder

function applyterm!(builder::IJVbuilder{Tv,T,E,L}, (term, sample, sublats)) where {Tv,T,E,L}
    L == 0 || isonsite(term) || checkinfinite(term)
    ndistiter = ndists(term, Val(L))
    sublatsiter = filteredsublats(sublats, builder.sysinfo)
    ijv = IJV{Tv}()
    ijvc = IJV{Tv}()
    for ndist in ndistiter
        foundndist = false
        addedconjugate = false
        for (s1, s2) in sublatsiter
            isvalidlink(term, (s1, s2), ndist) || continue
            dist = builder.lattice.bravais.matrix * ndist
            for (jsource, rsource) in enumerate(builder.lattice.sublats[s2].sites)
                for itarget in _targets(term, rsource, dist, jsource, (s1, s2), builder)
                    isselfhopping(term, (s1, s2), (itarget, jsource), ndist) && continue
                    rtarget = builder.lattice.sublats[s1].sites[itarget] + dist
                    r, dr   = _rdr(rsource, rtarget)
                    smatrix = term(r, dr)
                    iszero(smatrix) && continue
                    rowoffset = builder.sysinfo.offsets[s1] + (itarget - 1) * builder.sysinfo.norbitals[s1]
                    coloffset = builder.sysinfo.offsets[s2] + (jsource - 1) * builder.sysinfo.norbitals[s2]
                    push!(ijv, smatrix, (rowoffset, coloffset))
                    addconjugate = needsconjugate(term, (s1, s2), ndist)
                    addconjugate && push!(ijvc, smatrix', (coloffset, rowoffset))
                    foundndist = true
                    addedconjugate = addedconjugate || addconjugate
                end
            end
        end
        if foundndist
            acceptcell!(ndistiter, ndist)
            if addedconjugate   # this implies ijv is not empty
                push!(builder, IJVN(ijv, ndist))
                push!(builder, IJVN(ijvc, negSVector(ndist)))
                ijv = IJV{Tv}()
                ijvc = IJV{Tv}()
            elseif !isempty(ijv)
                push!(builder, IJVN(ijv, ndist))
                ijv = IJV{Tv}()
            end
        end
    end
    return builder
end

ndists(::Hopping{F,S,Missing}, ::Val{L}) where {F,S,L} = BoxIterator(zero(SVector{L,Int}))
ndists(h::Hopping{F,S,D}, ::Val{L}) where {F,S,D,L} = h.ndists
ndists(::Onsite, ::Val{L}) where {L} = (zero(SVector{L,Int}),)

filteredsublats(::Missing, sysinfo) = allorderedpairs(values(sysinfo.namesdict))
filteredsublats(ss::NTuple{N,Tuple{Int,Int}}, sysinfo) where {N} = unique!(collect(ss))

isselfhopping(term, (s1, s2), (itarget, jsource), ndist) =
    itarget == jsource && s1 == s2 && !isonsite(term) && iszero(ndist)

isvalidlink(term, (s1, s2), ndist) = needsconjugate(term, (s1, s2), ndist) || (iszero(ndist) && s1 == s2)

needsconjugate(::Onsite, args...) = false
needsconjugate(::Hopping{F,Missing,Missing}, (s1, s2), ndist) where {F} = s1 != s2 || ispositive(ndist)
needsconjugate(::Hopping{F,S,Missing}, (s1, s2), ndist) where {F,S} = s1 != s2 || ispositive(ndist)
needsconjugate(::Hopping{F,Missing}, (s1, s2), ndist) where {F} = s1 != s2 || !iszero(ndist)
needsconjugate(::Hopping, (s1, s2), ndist) = s1 != s2 || !iszero(ndist)

function ispositive(ndist)
    result = false
    for i in ndist
        i == 0 || (result = i > 0; break)
    end
    return result
end

@inline checkinfinite(modelterm) = modelterm.range === missing && modelterm.ndist === missing &&
    throw(ErrorException("Tried to implement an infinite-range hopping on an unbounded lattice"))

    function checkblockdims(::SMatrix{N,M}, sdata, sublatsiter) where {N,M}
    @boundscheck for (s1, s2) in sublatsiter
        (N == sdata.norbitals[s1]) && (M == sdata.norbitals[s2]) ||
            throw(DimensionMismatch("Hamiltonian term of dimension ($N,$M) does not match orbital ($(sdata.norbitals[s1]), $(sdata.norbitals[s2]))."))
    end
end

function _targets(term::Hopping{F,S,D,R}, r2, dist, j, (s1, s2), builder) where {F,S,D,R<:Real}
    isassigned(builder.kdtrees, s1) || (builder.kdtrees[s1] = KDTree(builder.lattice.sublats[s1].sites))
    return inrange(builder.kdtrees[s1], r2 - dist, term.range)
end

_targets(term::Hopping{F,S,D,R}, r2, dist, j, (s1, s2), builder) where {F,S,D,R<:Missing} =
    eachindex(builder.sublats[s1].sites)
_targets(term::Onsite, r2, dist, j, (s1, s2), builder) = (j,)

function assembleblocks(builder::IJVbuilder{Tv,T,E,L}) where {Tv,T,E,L}
    dimh = sum(builder.sysinfo.dims)
    intraIJVN = get_or_add_IJVN(builder, zero(SVector{L,Int}))
    intra = Block(intraIJVN, dimh, builder.sysinfo)
    inters = Block{Tv,L}[]
    for ijvn in builder.IJVNs
        iszero(ijvn.ndist) || push!(inters, Block(ijvn, dimh, builder.sysinfo))
    end
    ijv = intraIJVN.IJV
    foreach(ijvn -> iszero(ijvn.ndist) || appendIJV!(ijv, ijvn.IJV), builder.IJVNs)
    matrix = sparse(ijv.I, ijv.J, ijv.V, dimh, dimh)
    # matrix = +(intra.matrix, [inter.matrix for inter in inters]...)  ## This is a bit slower
    boundary = extractboundary(matrix, intra, inters)
    return matrix, intra, inters, boundary
end

appendIJV!(ijv, ijv2::IJV) = (append!(ijv.I, ijv2.I); append!(ijv.J, ijv2.J); append!(ijv.V, ijv2.V))

function extractboundary(matrix, intra::Block{Tv}, inters) where {Tv}
    rowsintra = rowvals(intra.matrix)
    rowsmatrix = rowvals(matrix)
    dimh = size(matrix, 1)
    boundary = Tuple{Int,Int}[]
    for col = 1:dimh
        for inter in inters
            rowsn = rowvals(inter.matrix)
            for ptr in nzrange(inter.matrix, col)
                ptrintra = ptrmatrix = 0
                row = rowsn[ptr]
                for ptrm in nzrange(matrix, col)
                    rowsmatrix[ptrm] == row && (ptrmatrix = ptrm; break)
                end
                iszero(ptrmatrix) && throw(ErrorException(
                    "Unexpected: found element in intercell harmonic not present in work matrix"))
                for ptri in nzrange(intra.matrix, col)
                    rowsintra[ptri] == row && (ptrintra = ptri; break)
                end
                push!(boundary, (ptrmatrix, ptrintra))
                # ptrintra == 0 if inter element is only present in ptrmatrix
            end
        end
    end
    return unique!(sort!(boundary))
end

#######################################################################
# combine systems
#######################################################################

combine(ss...; kw...) = _combine(collectfirst(ss...)...; kw...)
combine(model::Model) = sys -> combine(sys, model)

_combine(systems, (m,)::Tuple{Model}; kw...) = _combine(systems, promote_model(m, systems...); kw...)

function _combine(systems::NTuple{N,System}, model::Model; kw...) where {N}
    sublats = _combinesublats(systems...; kw...)
    bravais = first(systems).lattice.bravais
    lattice = _lattice(bravais, sublats)
    sampledterms = _getsampledterms(systems...)
    sysinfo = SystemInfo(lattice, model, sampledterms...)
    sysinfo_withoutnew = SystemInfo(lattice, model)
    _copyoffsets!(sysinfo_withoutnew, sysinfo)
    builder = IJVbuilder(lattice, sysinfo_withoutnew)
    nsublats = 1
    for system in systems
        offset = sysinfo.offsets[nsublats]
        _combine_block(builder, system.hamiltonian.intra, offset)
        foreach(block -> _combine_block(builder, block, offset), system.hamiltonian.inters)
        nsublats += length(system.lattice.sublats)
    end
    hamiltonian = Operator(builder)
    velocity = boundaryoperator(hamiltonian)
    return System(lattice, hamiltonian, velocity, sysinfo)
end

_getsampledterms(system, systems...) = (system.sysinfo.sampledterms..., _getsampledterms(systems...)...)
_getsampledterms() = ()

function _copyoffsets!(sysinfo1, sysinfo2)
    copy!(sysinfo1.norbitals, sysinfo2.norbitals)
    copy!(sysinfo1.dims, sysinfo2.dims)
    copy!(sysinfo1.offsets, sysinfo2.offsets)
end

function _combinesublats(system::System{Tv,T,E}, systems...; checkbravais = true) where {Tv,T,E}
    !checkbravais || is_bravais_compatible(system, systems...) || throw(ErrorException("Systems with incompatible Bravais matrices"))
    sublats = deepcopy(system.lattice.sublats)
    foreach(sys -> append!(sublats, sys.lattice.sublats), systems)
    return sublats
end

_combine(systems, ::Tuple{}; kw...) = _combine(systems, Model(); kw...)

function _combine_block(builder, block, offset)
    ijvn = IJVN(block)
    ijvn.IJV.I .+= offset
    ijvn.IJV.J .+= offset
    push!(builder, ijvn)
end

is_bravais_compatible() = true
is_bravais_compatible(sys::System, ss::System...) = all(s -> isequal(sys.lattice.bravais, s.lattice.bravais), ss)

function Base.isequal(b1::Bravais{E,L}, b2::Bravais{E,L}) where {E,L}
    vs1 = ntuple(i -> b1.matrix[:, i], Val(L))
    vs2 = ntuple(i -> b2.matrix[:, i], Val(L))
    for v2 in vs2
        found = false
        for v1 in vs1
            (isapprox(v1, v2) || isapprox(v1, -v2)) && (found = true; break)
        end
        !found && return false
    end
    return true
end

#######################################################################
# Grow system
#######################################################################

struct BlockBuilder{Tv,L}
    ndist::SVector{L,Int}
    matrixbuilder::SparseMatrixBuilder{Tv}
end

BlockBuilder{Tv,L}(ndist::SVector, n::Int) where {Tv,L} = BlockBuilder{Tv,L}(ndist, SparseMatrixBuilder{Tv}(n, n))
BlockBuilder{Tv,L}(n::Int) where {Tv,L} = BlockBuilder{Tv,L}(zero(SVector{L,Int}), SparseMatrixBuilder{Tv}(n, n))

Block(bb::BlockBuilder, sysinfo) = Block(bb.ndist, sparse(bb.matrixbuilder), sysinfo)

struct OperatorBuilder{Tv,L}
    intra::BlockBuilder{Tv,L}
    inters::Vector{BlockBuilder{Tv,L}}
end

OperatorBuilder{Tv,L}(n::Int) where {Tv,L} =
    OperatorBuilder(BlockBuilder{Tv,L}(n), BlockBuilder{Tv,L}[])

    dim(o::OperatorBuilder) = size(o.intra.matrixbuilder, 1)

function Operator(b::OperatorBuilder, sysinfo)
    intra = Block(b.intra, sysinfo)
    inters = [Block(inter, sysinfo) for inter in b.inters]
    matrix = +(intra.matrix, [inter.matrix for inter in inters]...)
    boundary = extractboundary(matrix, intra, inters)
    return Operator(matrix, intra, inters, boundary)
end

_truefunc(r) = true

grow(; kw...) = sys -> grow(sys; kw...)
grow(sys::System{E,L}; supercell = SMatrix{L,0,Int}(), region = _truefunc) where {E,L} =
    grow(sys, supercell, region)
grow(sys::System{E,L}, supercell::Integer, region) where {E,L} =
    grow(sys, SMatrix{L,L,Int}(supercell*I), region)
grow(sys::System{E,L}, supercell::NTuple{L,Integer}, region) where {E,L} =
    grow(sys, SMatrix{L,L,Int}(Diagonal(SVector{L,Int}(supercell))), region)
grow(sys::System{E,L}, supercell::NTuple{N,NTuple{L,Integer}}, region) where {E,L,N} =
    grow(sys, toSMatrix(supercell...), region)
grow(sys::System, supercell, region) =
    throw(DimensionMismatch("Possible mismatch between `supercell` and system, or `region` not a function."))

function grow(sys::System{E,L}, supercell::SMatrix{L,L2,Int}, region::Function) where {E,L,L2}
    L2 < L && region == _truefunc && error("Unbounded fill region for $(L-L2) dimensions")
    sublats, sitemaps = _growsublats(sys, supercell, region)
    bravais = _growbravais(sys, supercell)
    sysinfo = _growsysinfo(sys.sysinfo, sublats)
    hamiltonian = _growhamiltonian(sys, supercell, sitemaps, sysinfo)
    velocity = boundaryoperator(hamiltonian)
    return System(Lattice(sublats, bravais), hamiltonian, velocity, sysinfo)
end

_growbravais(sys::System{Tv,T,E,L}, supercell) where {Tv,T,E,L} =
    Bravais(bravaismatrix(sys) * supercell)

function _growsysinfo(oldsysinfo::SystemInfo, sublats)
    sysinfo = deepcopy(oldsysinfo)
    for (i, sublat) in enumerate(sublats)
        sysinfo.nsites[i] = length(sublat.sites)
        sysinfo.dims[i] = sysinfo.norbitals[i] * sysinfo.nsites[i]
        sysinfo.offsets[i + 1] = sysinfo.offsets[i] + sysinfo.dims[i]
    end
    return sysinfo
end

const TOOMANYITERS = 10^8

function _growsublats(sys::System{E,L,T}, supercell, region) where {E,L,T}
    bravais = bravaismatrix(sys)
    iter = BoxIterator(zero(SVector{L,Int}))
    hasnoshadow = _hasnoshadow(supercell)
    firstfound = false
    counter = 0
    for ndist in iter   # We first compute the bounding box
        found = false
        counter += 1; counter == TOOMANYITERS && @warn "Region fill seems non-covergent, check `region`"
        for sublat in sys.lattice.sublats, site in sublat.sites
            r = bravais * ndist + site
            found = hasnoshadow(ndist) && region(r)
            if found || !firstfound
                acceptcell!(iter, ndist)
                firstfound = found
                break
            end
        end
    end
    bbox = boundingboxiter(iter)
    ranges = UnitRange.(bbox...)

    # sitemaps[sublat][oldsiteindex, cell...] is newsiteindex for a given sublat
    sitemaps = OffsetArray{Int,L+1,Array{Int,L+1}}[zeros(Int, 1:length(sublat.sites),
        ranges...) for sublat in sys.lattice.sublats]

    sublats = [empty(sublat) for sublat in sys.lattice.sublats]
    for (s, sublat) in enumerate(sys.lattice.sublats)
        counter = 0
        for cell in CartesianIndices(ranges) , (i, site) in enumerate(sublat.sites)
            ndist = SVector{L,Int}(Tuple(cell))
            r = site + bravais * ndist
            if hasnoshadow(ndist) && region(r)
                push!(sublats[s].sites, r)
                counter += 1
                sitemaps[s][i, Tuple(cell)...] = counter
            end
        end
    end
    return sublats, sitemaps
end

pinvint(s::SMatrix{N,0}) where {N} = (SMatrix{0,0,Int}(), 0)

function pinvint(s::SMatrix{N,M}) where {N,M}
    qrfact = qr(s)
    pinverse = inv(qrfact.R) * qrfact.Q'
    n = det(qrfact.R)^2
    iszero(n) && throw(ErrorException("Supercell is singular"))
    return round.(Int, n * inv(qrfact.R) * qrfact.Q'), round(Int, n)
end

# This is true whenever old ndist is perpendicular to new lattice
_hasnoshadow(supercell) = let invs = pinvint(supercell); ndist -> iszero(_newndist(ndist, invs)); end

_newndist(oldndist, (pinvs, n)) = fld.(pinvs * oldndist, n)
_newndist(oldndist, (pinvs, n)::Tuple{<:SMatrix{0,0},Int}) = SVector{0,Int}()

_wrappedndist(ndist, s, invs) = (nn = _newndist(ndist, invs); (nn, ndist - s * nn))

function _growhamiltonian(sys::System{E,L,T,Tv}, supercell::SMatrix{L,L2}, sitemaps, newsysinfo) where {E,L,T,Tv,L2}
    blocks = sort!(append!([sys.hamiltonian.intra], sys.hamiltonian.inters), by = block -> reverse(block.ndist))
    dimh = sum(s -> sys.sysinfo.norbitals[s] * maximum(sitemaps[s]), eachindex(sys.lattice.sublats))
    invsupercell = pinvint(supercell)
    opbuilder = OperatorBuilder{Tv,L2}(dimh)
    colsdone = 0
    for (s2, sublat2) in enumerate(sys.lattice.sublats)
        cartesian_s2 = CartesianIndices(sitemaps[s2])
        norb2 = sys.sysinfo.norbitals[s2]
        offset2 = sys.sysinfo.offsets[s2]
        for (j, newsitesrc) in enumerate(sitemaps[s2])
            iszero(newsitesrc) && continue  # skip column, source not in region
            ndistsrc = SVector{L,Int}(Base.tail(Tuple(cartesian_s2[j])))
            sitesrc = first(Tuple(cartesian_s2[j]))
            for orb2 in 1:norb2
                col = offset2 + (sitesrc - 1) * norb2 + orb2
                for block in blocks
                    ptrs = nzrange(block.matrix, col)
                    isempty(ptrs) && continue
                    rows = rowvals(block.matrix)
                    vals = nonzeros(block.matrix)
                    ndist = ndistsrc + block.ndist
                    newndist, wrappedndist = _wrappedndist(ndist, supercell, invsupercell)
                    newblock = get_or_add_blockbuilder(opbuilder, newndist, colsdone)
                    for ptr in ptrs
                        site, orb, s1 = tosite(rows[ptr], sys.sysinfo)
                        checkbounds(Bool, sitemaps[s1], site, Tuple(wrappedndist)...) || continue
                        newsitedest = sitemaps[s1][site, Tuple(wrappedndist)...]
                        iszero(newsitedest) && continue
                        row = torow(newsitedest, s1, newsysinfo) + orb - 1
                        pushtocolumn!(newblock.matrixbuilder, row, vals[ptr])
                    end
                end
                colsdone += 1
                finalizecolumn!(opbuilder.intra.matrixbuilder)
                foreach(block -> finalizecolumn!(block.matrixbuilder), opbuilder.inters)
            end
        end
    end
    return Operator(opbuilder, newsysinfo)
end

function get_or_add_blockbuilder(op::OperatorBuilder{Tv,L}, ndist, colsdone) where {Tv,L}
    if iszero(ndist)
        return op.intra
    else
        k = findfirstblock(op.inters, ndist)  # better than findfirst (less runtime variance)
        if iszero(k)    # not found, add a new block builder
            block = BlockBuilder{Tv,L}(ndist, dim(op))
            finalizecolumn!(block.matrixbuilder, colsdone)
            push!(op.inters, block)
            return block
        else
            return op.inters[k]
        end
    end
end
function findfirstblock(blocks, ndist)
    for n in eachindex(blocks)
        blocks[n].ndist == ndist && return n
    end
    return 0
end

#######################################################################
# bound
#######################################################################

bound(;kw...) = sys -> bound(sys; kw...)
bound(sys::System{E,L}; except = ()) where {E,L,L2} = _bound(sys, toSVector(Int, except))

function _bound(sys::System{E,L,T,Tv}, exceptaxes::SVector{L2,Int}) where {E,L,L2,T,Tv}
    sublats = deepcopy(sys.lattice.sublats)
    bravais = Bravais(sys.lattice.bravais.matrix[:, exceptaxes])
    lattice = Lattice(sublats, bravais)
    sysinfo = sys.sysinfo
    matrix = deepcopy(sys.hamiltonian.matrix)
    supercell = SMatrix{L,L,Int}(I)[:,exceptaxes]
    intra =   _projectblock(sys.hamiltonian.intra, supercell)
    inters = Block{Tv,L2}[_projectblock(inter, supercell) for inter in sys.hamiltonian.inters
              if _isselfshadow(inter.ndist, exceptaxes)]
    boundary = extractboundary(matrix, intra, inters)
    ham = Operator(matrix, intra, inters, boundary)
    vel = boundaryoperator(ham)
    return System(lattice, ham, vel, sysinfo)
end

# This is true if the part of ndist orthogonal to exceptaxes is zero
_isselfshadow(ndist::SVector{L}, exceptaxes) where {L} =
    all(i -> iszero(ndist[i]) || i in exceptaxes, 1:L)

function _projectblock(block, scell)
    Block(scell' * block.ndist, copy(block.matrix), block.nlinks)
end