text:
using StatisticalNonlocality, LinearAlgebra, FFTW, SparseArrays, JLD2
import StatisticalNonlocality: chebyshev, fourier_nodes, fourier_wavenumbers
import StatisticalNonlocality: droprelativezeros!

# transition rate
γ = 1e0  # 1e0 # 1e2
# diffusivity
κ = 1e-2 # 1e0 # 1e-2

# minimal configuration is M = 64, N = 8
N = 8 * 4
M = 8 * 4

# The first step is to build the block operators
a, b = 0, 4π
k = fourier_wavenumbers(N, L = b - a)
x = fourier_nodes(N, a = a, b = b)

ℱ = fft(I + zeros(N, N), 1)
ℱ⁻¹ = ifft(I + zeros(N, N), 1)
Dx = real.(ℱ⁻¹ * Diagonal(im .* k) * ℱ)
Dy = copy(Dx)
x = reshape(x, (N, 1))
y = reshape(x, (1, N))
∂y = kron(Dy, I + zeros(N, N))
∂x = kron(I + zeros(N, N), Dx)
Δ = ∂x^2 + ∂y^2

# convenience functions
advection_operator(ψ, ∂x, ∂y) = Diagonal(∂y * ψ[:]) * ∂x + Diagonal(-∂x * ψ[:]) * ∂y

makelocal(E) = Array(Diagonal(sum(E, dims = 2)[:]))

function localdiffusivity(E, N, M)
    local localdiff = makelocal(E)
    localdiff = [localdiff[i, i] for i = 1:size(E)[1]]
    local localdiff = reshape(localdiff, (N, M ))
    return localdiff
end

function avglocaldiffusivity(E, N, M)
    local localdiff = localdiffusivity(E, N, M)
    local localdiff = sum(localdiff, dims = 1)[:] ./ N # zero'th mode
    return localdiff
end

# the local diffusivity is proportional to 1/γ
function grabdiagonal(A)
    MM = minimum(size(A))
    diagA = zeros(MM)
    for i in 1:MM
        diagA[i] = A[i, i]
    end
    return diagA
end

# random phase wavenumbers
wavemax = 0.5 # change to 3 when automated
𝓀 = collect(-wavemax:0.5:wavemax)
𝓀ˣs = reshape(𝓀, (length(𝓀), 1))
𝓀ʸs = reshape(𝓀, (1, length(𝓀)))
A = @. 0.1 * (𝓀ˣs * 𝓀ˣs + 𝓀ʸs * 𝓀ʸs)^(-11 / 12)
A[A.==Inf] .= 0.0

tic = Base.time()
for (i, 𝓀ˣ) in enumerate(𝓀ˣs), (j, 𝓀ʸ) in enumerate(𝓀ʸs)
    filename = "lots_of_diffusivities_kx_" * string(𝓀ˣ) * "ky_" * string(𝓀ʸ) * ".jld2"
    println("-----------")
    println("currently on wavenumber 𝓀ˣ=", 𝓀ˣ, ", 𝓀ʸ=", 𝓀ʸ)
    println("The amplitude is ", A[i, j])
    println("-----------")
    U₀ = A[i, j]

    ψ¹ = @. U₀ * sin(𝓀ˣ * x + 𝓀ʸ * y)
    ψ² = @. U₀ * cos(𝓀ˣ * x + 𝓀ʸ * y)
    ψ³ = @. -U₀ * sin(𝓀ˣ * x + 𝓀ʸ * y)
    ψ⁴ = @. -U₀ * cos(𝓀ˣ * x + 𝓀ʸ * y)

    A¹ = advection_operator(ψ¹, ∂x, ∂y)
    A² = advection_operator(ψ², ∂x, ∂y)
    A³ = advection_operator(ψ³, ∂x, ∂y)
    A⁴ = advection_operator(ψ⁴, ∂x, ∂y)

    L = [(γ*I-κ*Δ) 0.0*I 0.5*A¹; 0.0*I (γ*I-κ*Δ) -0.5*A²; A¹ -A² (2*γ*I-κ*Δ)]
    G = L \ I

    # Create the block linear operator for components of the diffusivity tensor
    u¹ = Diagonal(∂y * ψ¹[:])
    u² = Diagonal(∂y * ψ²[:])
    v¹ = Diagonal(-∂x * ψ¹[:])
    v² = Diagonal(-∂x * ψ²[:])

    U = -0.5 * [u¹; u²; 0 * I]
    Uᵀ = [u¹ u² 0 * I]
    V = -0.5 * [v¹; v²; 0 * I]
    Vᵀ = [v¹ v² 0 * I]

    EF¹¹ = Uᵀ * G * U # flux in x due to gradients in x
    EF¹² = Uᵀ * G * V # flux in x due to gradients in z
    EF²¹ = Vᵀ * G * U # flux in z due to gradients in x
    EF²² = Vᵀ * G * V # flux in z due to gradients in z

    E = [EF¹¹ EF¹²; EF²¹ EF²²]

    λE = eigvals(E)
    bools = real.(λE) .≤ eps(1e3 * maximum(abs.(λE)))
    println(sum(bools) == length(λE))

    # OUTPUT
    # always overwrite
    data_directory = "data"
    mkpath(data_directory)
    filepath = data_directory * "/" * filename
    if isfile(filepath)
        rm(filepath)
    end
    ##
    file = jldopen(filepath, "a+")
    # diffusivity
    groupname = "diffusivity"
    JLD2.Group(file, groupname)
    file[groupname]["K11"] = EF¹¹
    file[groupname]["K12"] = EF¹²
    file[groupname]["K21"] = EF²¹
    file[groupname]["K22"] = EF²²
    # parameters
    groupname = "parameters"
    JLD2.Group(file, groupname)
    file[groupname]["γ"] = γ
    file[groupname]["κ"] = κ
    # grid
    groupname = "grid"
    JLD2.Group(file, groupname)
    file[groupname]["x"] = x
    file[groupname]["y"] = y
    # transition matrix
    groupname = "transition"
    JLD2.Group(file, groupname)

    T = [-1 1 0 0; 0 -1 1 0; 0 0 -1 1; 1 0 0 -1]
    sT = (T + T') / 2

    file[groupname]["T"] = γ * sT
    # streamfunction
    groupname = "streamfunction"
    JLD2.Group(file, groupname)
    file[groupname]["ψ¹"] = ψ¹
    file[groupname]["ψ²"] = ψ²
    file[groupname]["ψ³"] = ψ³
    file[groupname]["ψ⁴"] = ψ⁴
    # velocities
    groupname = "velocities"
    JLD2.Group(file, groupname)
    file[groupname]["u¹"] = reshape(grabdiagonal(u¹), size(ψ¹))
    file[groupname]["u²"] = reshape(grabdiagonal(u²), size(ψ¹))
    file[groupname]["v¹"] = reshape(grabdiagonal(v¹), size(ψ¹))
    file[groupname]["v²"] = reshape(grabdiagonal(v²), size(ψ¹))

    # local diffusivities
    groupname = "localdiffusivity"
    JLD2.Group(file, groupname)
    file[groupname]["κ11"] = localdiffusivity(EF¹¹, N, N)
    file[groupname]["κ12"] = localdiffusivity(EF¹², N, N)
    file[groupname]["κ21"] = localdiffusivity(EF²¹, N, N)
    file[groupname]["κ22"] = localdiffusivity(EF²², N, N)

    # zonally averaged local diffusivities
    groupname = "averagelocaldiffusivity"
    JLD2.Group(file, groupname)
    file[groupname]["κ11"] = avglocaldiffusivity(EF¹¹, N, N)
    file[groupname]["κ12"] = avglocaldiffusivity(EF¹², N, N)
    file[groupname]["κ21"] = avglocaldiffusivity(EF²¹, N, N)
    file[groupname]["κ22"] = avglocaldiffusivity(EF²², N, N)

    close(file)

end
toc = Base.time()

println("The time for the simulation is ", (toc - tic)/60, " minutes ")

max_stars_repo_name: sandreza/StatisticalNonlocality

max_issues_repo_name: sandreza/StatisticalNonlocality

max_forks_repo_name: sandreza/StatisticalNonlocality


text:
import BinDeps

rm′(f) = (isdir(f) || isfile(f)) && rm(f, recursive = true)

const version = "4.0.4"

folder() = normpath(joinpath(@__FILE__, "../../../deps"))

@static if Sys.isapple()
    uninstall() = map(rm′, filter(x -> !endswith(x, "build.jl"), readdir(folder())))
else
    uninstall() = rm′(joinpath(folder(), "atom"))
end

isinstalled() = Sys.isapple() ?
    isfile(joinpath(folder(), "version")) :
    isdir(joinpath(folder(), "atom"))

function install()
  dir = folder()
  if Sys.isapple()
    _icons = normpath(joinpath(@__FILE__, "../../../res/julia-icns.icns"))
  end
  !isdir(dir) && mkpath(dir)
  uninstall()
  cd(dir) do
    download(x) = run(BinDeps.download_cmd(x, basename(x)))

    download("http://junolab.s3.amazonaws.com/blink/julia.png")

    if Sys.isapple()
      file = "electron-v$version-darwin-x64.zip"
      download("https://github.com/electron/electron/releases/download/v$version/$file")
      run(`unzip -q $file`)
      rm(file)
      if isdir("Julia.app") # Issue # 178
          rm("Julia.app", recursive=true)
      end
      run(`mv Electron.app Julia.app`)
      run(`mv Julia.app/Contents/MacOS/Electron Julia.app/Contents/MacOS/Julia`)
      run(`sed -i.bak 's/Electron/Julia/' Julia.app/Contents/Info.plist`)
      run(`cp $_icons Julia.app/Contents/Resources/electron.icns`)
      run(`touch Julia.app`)  # Apparently this is necessary to tell the OS to double-check for the new icons.
    end

    if Sys.iswindows()
      arch = Int == Int64 ? "x64" : "ia32"
      file = "electron-v$version-win32-$arch.zip"
      download("https://github.com/electron/electron/releases/download/v$version/$file")
      try
        run(BinDeps.unpack_cmd(file, "atom", ".zip", nothing))
      catch e
        if e isa Base.IOError && e.code == Base.UV_ENOENT # No such file or directory
          error("Unable to find 7z.exe. Try running `make win-extras` if you built Julia from source.")
        else
          rethrow(e)
        end
      finally
        rm(file)
      end
    end

    if Sys.islinux()
      arch = Int == Int64 ? "x64" : "ia32"
      file = "electron-v$version-linux-$arch.zip"
      download("https://github.com/electron/electron/releases/download/v$version/$file")
      run(`unzip -q $file -d atom`)
      rm(file)
    end
  end
end

ext: jl


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Random
using SuiteSparse, LinearAlgebra, SparseArrays

include("umfpack.jl")
include("cholmod.jl")
include("spqr.jl")


text:
#---------------------------------------------------------------------#
#This code computes the Legendre Polynomials and its 1st and 2nd derivatives
#Written by F.X. Giraldo on 4/2008
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey, CA 93943-5216
#---------------------------------------------------------------------#
function legendre_poly(p,x)

    L1=0;L1_1=0;L1_2=0
    L0=1;L0_1=0;L0_2=0

    for i=1:p
        L2=L1;L2_1=L1_1;L2_2=L1_2
        L1=L0;L1_1=L0_1;L1_2=L0_2
        a=(2*i-1)/i
        b=(i-1)/i
        L0=a*x*L1 - b*L2
        L0_1=a*(L1 + x*L1_1) - b*L2_1
        L0_2=a*(2*L1_1 + x*L1_2) - b*L2_2
    end
    return (L0,L0_1,L0_2)
end


text:
# This file was generated, do not modify it. # hide
minrecall = callbacks3.hyperparameters
optimize!(G1, minrecall, queries=equeries)
optimize!(G2, minrecall, queries=equeries)
optimize!(G3, minrecall, queries=equeries)

@time I1, D1 = searchbatch(G1, queries, k)
@time I2, D2 = searchbatch(G2, queries, k)
@time I3, D3 = searchbatch(G3, queries, k)
nothing # hide

@show macrorecall(goldI, I1)
@show macrorecall(goldI, I2)
@show macrorecall(goldI, I3)

max_stars_repo_name: sadit/SimilaritySearchDemos

max_issues_repo_name: sadit/SimilaritySearchDemos

max_forks_repo_name: sadit/SimilaritySearchDemos


text:
using Test
using Dash
using JSON2
using SnoopCompile
using HTTP
@testset "dev" begin
external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"]

styles = (
    pre = (
        border = "thin lightgrey solid",
        overflowX = "scroll"
    ),
)

app = dash(external_stylesheets=external_stylesheets)
app.layout = html_div() do
        html_div(className="three columns") do
        end,
        html_div(className="three columns") do
        end,
        html_div(className="three columns") do
        end,
        html_div(className="three columns") do
        end
    end

function pretty_json(t)
    io = IOBuffer()
    JSON2.pretty(io, JSON2.write(t), offset = 2)
    return String(take!(io))
end

#inf_timing = @snoopi tmin = 0.1 run_server(app, "0.0.0.0", 8080)
precompile(Dash.make_handler, (Dash.DashApp,))
println("====")
@time handler = Dash.make_handler(app)
@time handler = Dash.make_handler(app)
println("====")
inf_timing = @snoopi HTTP.handle(handler, HTTP.Request("GET", "/"))
t = 0.
for itm in inf_timing
    @show itm
    t += itm[1]
end
@show t
pc = SnoopCompile.parcel(inf_timing)
@show pc
SnoopCompile.write("venv/precompile", pc)
end


text:
using LinearAlgebra
include("gates.jl")

# common states
ket_0 = Complex[1; 0]
ket_1 = Complex[0; 1]
ket_plus = Complex[1; 1] ./ sqrt(2)
ket_minus = Complex[1; -1] ./ sqrt(2)

# TODO
# struct Op
#     U::Matrix{Complex}
#     isUnitary::Bool
# end

"""
    conj_transpose(s)

conj_transpose is a clearer name for adjoint for students.
it is identical to `s'`
"""
function conj_transpose(s::Union{AbstractVecOrMat})
    conj.(transpose(s))
end

"""
    commutator(A, B)

commutator(A, B) := A * B - B * A
Takes in two operators, if the result is zero, the operators commute.
"""
function commutator(A, B)
    A * B - B * A
end

"""
    ops_commute(A, B)

ops_commute(A, B) := commutator(A,B) ≈ fill(0., (dim, dim))
for ops to commute they must be square and equal dim, so dim
    is arbitrarily size(A, 1)
"""
function ops_commute(A, B)
    if size(A) != size(B)
        return false
    end
    dim = size(A, 1)
    commutator(A,B) ≈ fill(0., (dim, dim))
end

    """
    is_unitary(s)

A unitary matrix U is defined by U†U = UU† = I.
returns true if `s` is unitary.
Note `s` must be square.
"""
# TODO: fix 1/sqrt(2) and normalization
function is_unitary(s)
    LinearAlgebra.checksquare(s)
    s * conj_transpose(s) ≈ Id
end

"""
    is_hermitian(s)

A hermitian matrix U is defined by U = †U.
returns true if s is hermitian.
"""
# throws error if not hermitian
function is_hermitian(s)
    LinearAlgebra.checksquare(s);
    s ≈ conj_transpose(s)
end

"""
    is_normed(s)

if sum(s .^ 2)) ≈ 1., returns true
"""
# throws error if s is not normalized
function is_normed(s)
    mag = convert.(Float64, sum(s .^ 2))
    @assert mag ≈ 1.
end

"""
    is_linear(s)

"""
function is_linear()

end

"""
    inner(a, b)

given two equal dimension column vectors, a and b,
returns the inner product <a|b> where <a| == conj_transpose(a)
this is a duplicate function of dot(a, b) or (a ⋅ b).
"""
# use ⋅ ?
function inner(a, b)
    conj_transpose(a) * b
end

"""
    outer(a, b)

given two equal dimension column vectors, |a> and |b>,
returns the inner product |a><b| where <b| == conj_transpose(b)
note, this function does NOT return the same as kron(a, b).
outer returns a matrix, whereas kron returns a vector
"""
# todo splat n ops
# add assert equal dims
function outer(a, b)
    a * conj_transpose(b)
end

"""
    Cgate(U)

Cgate(U) creates a controlled U gate, given a unitary U.
"""
Cgate(U) = kron(outer(ket_0, ket_0), Id) + kron(outer(ket_1, ket_1), U)

"""
    bell_states()

returns the four bell states as vectors in order of binary increasing
qubit input states (00, 01, 10, 11)
"""
function bell_states()
    # # bell states
    # 00: 00 + 11
    # 10: 00 - 11
    # 01: 01 + 10
    # 11: 01 - 10
    β_00 = Cgate(X) * kron(H * ket_0, ket_0)
    β_01 = Cgate(X) * kron(H * ket_0, ket_1)
    β_10 = Cgate(X) * kron(H * ket_1, ket_0)
    β_11 = Cgate(X) * kron(H * ket_1, ket_1)
    β_00, β_01, β_10, β_11
end

"""
    cauchy_schwarz(u, v)

given arbitrary column vectors `|u⟩, |v⟩`, we want to show that,

⟨u|u⟩⟨v|v⟩ ≥ |⟨u|v⟩|²
"""

function cauchy_schwarz(u, v)

end
# function eulerform(U)
    # thm
    # U single qubit unitary ⟹ ∃ α, β, γ, δ ∈ ℝ, st
    # U = exp(iα) * [exp(-i\)]

# end

max_stars_repo_name: anandijain/MyQuantum

max_issues_repo_name: anandijain/MyQuantum

max_forks_repo_name: anandijain/MyQuantum


text:
using Test, MPI
using LinearAlgebra

using ClimateMachine
using ClimateMachine.MPIStateArrays

ClimateMachine.init()
const ArrayType = ClimateMachine.array_type()
const mpicomm = MPI.COMM_WORLD

mpisize = MPI.Comm_size(mpicomm)
mpirank = MPI.Comm_rank(mpicomm)

@testset "MPIStateArray reductions" begin

    localsize = (4, 6, 8)

    A = Array{Float32}(reshape(1:prod(localsize), localsize))
    globalA = vcat([A for _ in 1:mpisize]...)

    QA = MPIStateArray{Float32}(mpicomm, ArrayType, localsize...)
    QA .= A

    @test norm(QA, 1) ≈ norm(globalA, 1)
    @test norm(QA) ≈ norm(globalA)
    @test norm(QA, Inf) ≈ norm(globalA, Inf)

    @test norm(QA; dims = (1, 3)) ≈ mapslices(norm, globalA; dims = (1, 3))
    @test norm(QA, 1; dims = (1, 3)) ≈
          mapslices(S -> norm(S, 1), globalA, dims = (1, 3))
    @test norm(QA, Inf; dims = (1, 3)) ≈
          mapslices(S -> norm(S, Inf), globalA, dims = (1, 3))

    B = Array{Float32}(reshape(reverse(1:prod(localsize)), localsize))
    globalB = vcat([B for _ in 1:mpisize]...)

    QB = similar(QA)
    QB .= B

    @test isapprox(euclidean_distance(QA, QB), norm(globalA .- globalB))
    @test isapprox(dot(QA, QB), dot(globalA, globalB))

    C = fill(Float32(mpirank + 1), localsize)
    globalC = vcat([fill(i, localsize) for i in 1:mpisize]...)
    QC = similar(QA)
    QC .= C

    @test sum(QC) == sum(globalC)
    @test Array(sum(QC; dims = (1, 3))) == sum(globalC; dims = (1, 3))
    @test maximum(QC) == maximum(globalC)
    @test Array(maximum(QC; dims = (1, 3))) == maximum(globalC; dims = (1, 3))
    @test minimum(QC) == minimum(globalC)
    @test Array(minimum(QC; dims = (1, 3))) == minimum(globalC; dims = (1, 3))
end


text:
module Get_products_test

include()

function createDict()
    io = getProductsPages( authenticate("davidefavaro","Tirocinio"), 200 )
    dict = getPageProducts( io )
    return dict
end

end # module

max_stars_repo_name: DavideFavaro/Tirocinio

max_issues_repo_name: DavideFavaro/Tirocinio

max_forks_repo_name: DavideFavaro/Tirocinio


text:
using Combinatorics
using InteractiveUtils

"""
    ops::Tuple{String}
    delta::TransferConfig{1}
    parity::Tuple{Int}
    clusters::Tuple{Cluster}
    ints::Array{Float64}
    cache::Dict

input:
- delta = list of change of Na,Nb,state
			e.g., [(-1,-1),(1,1),(0,0)] means alpha and beta transition
			from cluster 1 to 2, cluster 3 is fock diagonal
- ops   = list of operators
			e.g., ["ab","AB",""]

- ints  = tensor containing the integrals for this block
			e.g., ndarray([p,q,r,s]) where p,q are in 1 and r,s are in 2

- data contained in object
		active: list of clusters which have non-identity operators
			this includes fock-diagonal couplings,
			e.g., ["Aa","","Bb"] would have active = [0,2]
- parity: does each operator have even or odd number of second quantized operators
"""
abstract type ClusteredTerm end

struct ClusteredTerm1B <: ClusteredTerm
    ops::Tuple{String}
    delta::TransferConfig{1}
    parity::Tuple{Int}
    clusters::Tuple{Cluster}
    ints::Array{Float64,1}
    cache::Dict
end

struct ClusteredTerm2B <: ClusteredTerm
    ops::Tuple{String,String}
    #delta::Tuple{Tuple{Int16,Int16},Tuple{Int16,Int16}}
    delta::TransferConfig{2}
    parity::Tuple{Int,Int}
    #active::Vector{Int16}
    clusters::Tuple{Cluster,Cluster}
    ints::Array{Float64,2}
    cache::Dict
end

struct ClusteredTerm3B <: ClusteredTerm
    ops::Tuple{String,String,String}
    delta::TransferConfig{3}
    parity::Tuple{Int,Int,Int}
    #active::Vector{Int16}
    clusters::Tuple{Cluster,Cluster,Cluster}
    ints::Array{Float64,3}
    cache::Dict
end

struct ClusteredTerm4B <: ClusteredTerm
    ops::Tuple{String,String,String,String}
    delta::TransferConfig{4}
    parity::Tuple{Int,Int,Int,Int}
    clusters::Tuple{Cluster,Cluster,Cluster,Cluster}
    ints::Array{Float64,4}
    cache::Dict
end

#function ClusteredTerm(ops, delta::Vector{Tuple{Int}}, clusters, ints)
#end

function Base.display(t::ClusteredTerm1B)
    @printf( " 1B: %2i          :", t.clusters[1].idx)
    println(t.ops)
end
function Base.display(t::ClusteredTerm2B)
    @printf( " 2B: %2i %2i       :", t.clusters[1].idx, t.clusters[2].idx)
    println(t.ops)
end
function Base.display(t::ClusteredTerm3B)
    @printf( " 3B: %2i %2i %2i    :", t.clusters[1].idx, t.clusters[2].idx, t.clusters[3].idx)
    println(t.ops)
end
function Base.display(t::ClusteredTerm4B)
    @printf( " 4B: %2i %2i %2i %2i :", t.clusters[1].idx, t.clusters[2].idx, t.clusters[3].idx, t.clusters[4].idx)
    println(t.ops)
end

#######################################################################################################333
"""
    trans::Dict{TransferConfig,Vector{ClusteredTerm}}
"""
struct ClusteredOperator
    trans::Dict{TransferConfig,Vector{ClusteredTerm}}
end
ClusteredOperator() = ClusteredOperator(Dict{TransferConfig,Vector{ClusteredTerm}}())
function flush_cache(clustered_ham::ClusteredOperator)
    for (ftrans, terms) in clustered_ham.trans
        for term in terms
            for cached in keys(term.cache)
                delete!(term.cache, cached)
            end
            #display(term.cache)
        end
    end
end
Base.haskey(co::ClusteredOperator, i) = return haskey(co.trans,i)
Base.setindex!(co::ClusteredOperator, i, j) = return co.trans[j] = i
Base.getindex(co::ClusteredOperator, i) = return co.trans[i]
Base.iterate(co::ClusteredOperator, state=1) = iterate(co.trans, state)
flush_cache(h::Dict{TransferConfig,Vector{ClusteredTerm}}) = flush_cache(ClusteredOperator(h))
mem_used_by_cache(h::Dict{TransferConfig,Vector{ClusteredTerm}}) = mem_used_by_cache(ClusteredOperator(h))
Base.convert(::Type{ClusteredOperator}, input::Dict{TransferConfig,Vector{ClusteredTerm}}) = ClusteredOperator(input)
@inline Base.length(co::ClusteredOperator) = return length(co.trans)
function Base.display(co::ClusteredOperator)
    for (ftrans, terms) in co
        display(ftrans)
        for term in terms
            display(term)
        end
    end
end

function mem_used_by_cache(h::ClusteredOperator)
    mem = 0.0
    for (ftrans, terms) in h.trans
        for term in terms
            mem += sizeof(term.cache)
        end
    end
    return mem
end

function bubble_sort(inp)
    #={{{=#
    cmpcount, swapcount = 0, 0
    blist = copy(inp)
    bperm = collect(1:length(inp))
    for j in 1:length(blist)
        for i in 1:(length(blist)-j)
            cmpcount += 1
            if blist[i] > blist[i+1]
                swapcount += 1
                blist[i], blist[i+1] = blist[i+1], blist[i]
                bperm[i], bperm[i+1] = bperm[i+1], bperm[i]
            end
        end
    end
    return bperm, swapcount
#=}}}=#
end

"""
    extract_terms(ints::InCoreInts, clusters)

Extract all ClusteredTerm types from a given 1e integral tensor
and a list of clusters
returns `terms::Dict{TransferConfig,Vector{ClusteredTerm}}`
"""
function extract_ClusteredTerms(ints::InCoreInts, clusters)
    norb = 0
    for ci in clusters
        norb += length(ci)
    end
    length(size(ints.h1)) == 2 || throw(Exception)
    size(ints.h1,1) == norb || throw(Exception)
    size(ints.h1,2) == norb || throw(Exception)

    #terms = Dict{TransferConfig,Vector{ClusteredTerm}}()
    terms = ClusteredOperator()
    #terms = Dict{Vector{Tuple{Int16,Int16}},Vector{ClusteredTerm}}()
    #terms = Dict{Tuple,Vector{ClusteredTerm}}()
    n_clusters = length(clusters)
    ops_a = Array{String}(undef,n_clusters)
    ops_b = Array{String}(undef,n_clusters)
    fill!(ops_a,"")
    fill!(ops_b,"")

    zero_fock = TransferConfig([(0,0) for i in clusters])
    #zero_fock::Vector{Tuple{Int16,Int16}} = [(0,0) for i in clusters]
    #zero_fock = Tuple([(0,0) for i in clusters])
    terms[zero_fock] = Vector{ClusteredTerm}()

    # 1-body terms
    if true
        for ci in clusters
#={{{=#
            # instead of forming p'q and p'q'sr just precontract and keep them in
            # ClusterOps
            term = ClusteredTerm1B(("H",), ((0,0),), (0,), (ci,), zeros(1),Dict())
            push!(terms[zero_fock],term)
#=}}}=#
        end
    end

    # 2-body 1-electron terms
    if true
        for ci in clusters
            for cj in clusters
                #={{{=#
                i = ci.idx
                j = cj.idx

                i < j || continue

                spin_cases =[["A","a"],
                             ["B","b"]
                            ]

                fock_cases =[[(1,0),(-1,0)],
                             [(0,1),(0,-1)]
                            ]

                termstr = []
                append!(termstr,unique(permutations([ci,cj])))
                #append!(termstr,unique(permutations([ci,cj,cj,cj])))
                #append!(termstr,unique(permutations([ci,ci,ci,cj])))

                #
                #   (pr|qs) p'q'sr
                #
                for term in termstr

                    #
                    #   find permutations and sign needed to sort the indices
                    #   such that clusters increase from left to right
                    perm, countswap = bubble_sort(term)
                    perm == sortperm(term, alg=MergeSort)|| throw(Exception)

                    permsign = 1
                    if countswap%2 != 0
                        permsign = -1
                    end

                    #vprqs = view(ints.h2,term[1].orb_list, term[2].orb_list, term[3].orb_list, term[4].orb_list)
                    hpq = view(ints.h1,term[1].orb_list, term[2].orb_list)

                    #
                    # now align (pqsr) ints so that they align with indices from operators after sorting
                    # in this ordering, one can simply contract by sum(v .* d)
                    h = permsign .* permutedims(hpq,perm)

                    for sidx in 1:length(spin_cases)
                        oper = spin_cases[sidx][perm]
                        fock = fock_cases[sidx][perm]
                        oper1 = ""
                        oper2 = ""
                        fock1 = [0,0]
                        fock2 = [0,0]
                        for cidx in 1:length(term[perm])
                            if term[perm][cidx] == ci
                                oper1 *= oper[cidx]
                                fock1 .+= fock[cidx]
                            elseif term[perm][cidx] == cj
                                oper2 *= oper[cidx]
                                fock2 .+= fock[cidx]
                            else
                                throw(Exception)
                            end
                        end

                        parity1 = 0
                        parity2 = 0
                        if length(oper1)%2 != 0
                            parity1 = 1
                        end
                        if length(oper2)%2 != 0
                            parity2 = 1
                        end
                        parity = (parity1, parity2)

                        clusteredterm = ClusteredTerm2B((oper1,oper2), (Tuple(fock1),Tuple(fock2)), parity, (ci, cj), h, Dict())
                        #display(clusteredterm)
                        focktrans = replace(zero_fock, (ci.idx, cj.idx), (fock1, fock2))
#                        focktrans = [zero_fock...]
#                        focktrans[ci.idx] = Tuple(fock1)
#                        focktrans[cj.idx] = Tuple(fock2)
#                        focktrans = TransferConfig(focktrans)
                        if haskey(terms,focktrans)
                            push!(terms[focktrans], clusteredterm)
                        else
                            terms[focktrans] = [clusteredterm]
                        end

                    end
                end
            end
        end
        #=}}}=#
    end

    # 2-body 2-electron terms
    if true
        for ci in clusters
            for cj in clusters
                #={{{=#
                i = ci.idx
                j = cj.idx

                i < j || continue

                spin_cases =[["A","A","a","a"],
                             ["B","B","b","b"],
                             ["A","B","b","a"],
                             ["B","A","a","b"]
                            ]

                fock_cases =[[(1,0),(1,0),(-1,0),(-1,0)],
                             [(0,1),(0,1),(0,-1),(0,-1)],
                             [(1,0),(0,1),(0,-1),(-1,0)],
                             [(0,1),(1,0),(-1,0),(0,-1)]
                            ]

                termstr = []
                append!(termstr,unique(permutations([ci,ci,cj,cj])))
                append!(termstr,unique(permutations([ci,cj,cj,cj])))
                append!(termstr,unique(permutations([ci,ci,ci,cj])))

                #
                #   (pr|qs) p'q'sr
                #
                for term in termstr

                    #
                    #   find permutations and sign needed to sort the indices
                    #   such that clusters increase from left to right
                    perm, countswap = bubble_sort(term)
                    perm == sortperm(term, alg=MergeSort)|| error("problem with bubble_sort")

                    permsign = 1
                    if countswap%2 != 0
                        permsign = -1
                    end

                    vprqs = view(ints.h2,term[1].orb_list, term[4].orb_list, term[2].orb_list, term[3].orb_list)
                    #
                    # align (prqs) ints so that they align with indices from operators before sorting
                    vpqsr = permutedims(vprqs,[1,3,4,2])

                    #
                    # now align (pqsr) ints so that they align with indices from operators after sorting
                    # in this ordering, one can simply contract by sum(v .* d)
                    v = (.5 * permsign) .* permutedims(vpqsr,perm)

                    #
                    # now reshape ints so they can be contracted like gamma(pqr) V(pqr,s) gamma(s)
                    newshape = [1,1]

                    for cidx in 1:length(term[perm])
                        if term[perm][cidx] == ci
                            newshape[1] *= size(v,cidx)
                        elseif term[perm][cidx] == cj
                            newshape[2] *= size(v,cidx)
                        else
                            throw(Exception)
                        end
                    end

                    for sidx in 1:length(spin_cases)
                        oper = spin_cases[sidx][perm]
                        fock = fock_cases[sidx][perm]
                        oper1 = ""
                        oper2 = ""
                        fock1 = [0,0]
                        fock2 = [0,0]
                        for cidx in 1:length(term[perm])
                            if term[perm][cidx] == ci
                                oper1 *= oper[cidx]
                                fock1 .+= fock[cidx]
                            elseif term[perm][cidx] == cj
                                oper2 *= oper[cidx]
                                fock2 .+= fock[cidx]
                            else
                                throw(Exception)
                            end
                        end
                        vcurr = deepcopy(v)

                        if true
                            if oper1 == "BA"
                                oper1 = "AB"
                                vcurr = -permutedims(vcurr,[2,1,3,4])
                            elseif oper1 == "ba"
                                oper1 = "ab"
                                vcurr = -permutedims(vcurr,[2,1,3,4])
                            elseif oper1 == "BAa"
                                oper1 = "ABa"
                                vcurr = -permutedims(vcurr,[2,1,3,4])
                            elseif oper1 == "BAb"
                                oper1 = "ABb"
                                vcurr = -permutedims(vcurr,[2,1,3,4])
                            elseif oper1 == "Bab"
                                oper1 = "Bba"
                                vcurr = -permutedims(vcurr,[1,3,2,4])
                            elseif oper1 == "Aab"
                                oper1 = "Aba"
                                vcurr = -permutedims(vcurr,[1,3,2,4])
                            end

                            if oper2 == "BA"
                                oper2 = "AB"
                                vcurr = -permutedims(vcurr,[1,2,4,3])
                            elseif oper2 == "ba"
                                oper2 = "ab"
                                vcurr = -permutedims(vcurr,[1,2,4,3])
                            elseif oper2 == "BAa"
                                oper2 = "ABa"
                                vcurr = -permutedims(vcurr,[1,3,2,4])
                            elseif oper2 == "BAb"
                                oper2 = "ABb"
                                vcurr = -permutedims(vcurr,[1,3,2,4])
                            elseif oper2 == "Bab"
                                oper2 = "Bba"
                                vcurr = -permutedims(vcurr,[1,2,4,3])
                            elseif oper2 == "Aab"
                                oper2 = "Aba"
                                vcurr = -permutedims(vcurr,[1,2,4,3])
                            end
                        end

                        vcurr = copy(reshape(vcurr,newshape...))

                        parity1 = 0
                        parity2 = 0
                        if length(oper1)%2 != 0
                            parity1 = 1
                        end
                        if length(oper2)%2 != 0
                            parity2 = 1
                        end
                        parity = (parity1, parity2)

                        clusteredterm = ClusteredTerm2B((oper1,oper2), (Tuple(fock1),Tuple(fock2)), parity, (ci, cj), vcurr, Dict())
                        #display(clusteredterm)
                        focktrans = replace(zero_fock, (ci.idx, cj.idx), (fock1, fock2))
                        if haskey(terms,focktrans)
                            push!(terms[focktrans], clusteredterm)
                        else
                            terms[focktrans] = [clusteredterm]
                        end

                    end

                    #all(isapprox(vpqsr,0.0)) || continue
                end

                #=}}}=#
            end
        end
    end

    # 3-body 2-electron terms
    if true
        for ci in clusters
            for cj in clusters
                for ck in clusters
                    #={{{=#
                    i = ci.idx
                    j = cj.idx
                    k = ck.idx

                    i < j < k || continue

                    spin_cases =[["A","A","a","a"],
                                 ["B","B","b","b"],
                                 ["A","B","b","a"],
                                 ["B","A","a","b"]
                                ]

                    fock_cases =[[(1,0),(1,0),(-1,0),(-1,0)],
                                 [(0,1),(0,1),(0,-1),(0,-1)],
                                 [(1,0),(0,1),(0,-1),(-1,0)],
                                 [(0,1),(1,0),(-1,0),(0,-1)]
                                ]

                    termstr = []
                    append!(termstr,unique(permutations([ci,ci,cj,ck])))
                    append!(termstr,unique(permutations([ci,cj,cj,ck])))
                    append!(termstr,unique(permutations([ci,cj,ck,ck])))

                    #
                    #   (pr|qs) p'q'sr
                    #
                    for term in termstr

                        #
                        #   find permutations and sign needed to sort the indices
                        #   such that clusters increase from left to right
                        perm, countswap = bubble_sort(term)
                        perm == sortperm(term, alg=MergeSort)|| error("problem with bubble_sort")

                        permsign = 1
                        if countswap%2 != 0
                            permsign = -1
                        end

                        vprqs = view(ints.h2,term[1].orb_list, term[4].orb_list, term[2].orb_list, term[3].orb_list)
                        #
                        # align (prqs) ints so that they align with indices from operators before sorting
                        vpqsr = permutedims(vprqs,[1,3,4,2])

                        #
                        # now align (pqsr) ints so that they align with indices from operators after sorting
                        # in this ordering, one can simply contract by sum(v .* d)
                        v = (.5 * permsign) .* permutedims(vpqsr,perm)

                        #
                        # now reshape ints so they can be contracted like gamma(pqr) V(pqr,s) gamma(s)
                        newshape = [1,1,1]

                        for cidx in 1:length(term[perm])
                            if term[perm][cidx] == ci
                                newshape[1] *= size(v,cidx)
                            elseif term[perm][cidx] == cj
                                newshape[2] *= size(v,cidx)
                            elseif term[perm][cidx] == ck
                                newshape[3] *= size(v,cidx)
                            else
                                throw(Exception)
                            end
                        end

                        for sidx in 1:length(spin_cases)
                            oper = spin_cases[sidx][perm]
                            fock = fock_cases[sidx][perm]
                            oper1 = ""
                            oper2 = ""
                            oper3 = ""
                            fock1 = [0,0]
                            fock2 = [0,0]
                            fock3 = [0,0]
                            for cidx in 1:length(term[perm])
                                if term[perm][cidx] == ci
                                    oper1 *= oper[cidx]
                                    fock1 .+= fock[cidx]
                                elseif term[perm][cidx] == cj
                                    oper2 *= oper[cidx]
                                    fock2 .+= fock[cidx]
                                elseif term[perm][cidx] == ck
                                    oper3 *= oper[cidx]
                                    fock3 .+= fock[cidx]
                                else
                                    throw(Exception)
                                end
                            end
                            vcurr = deepcopy(v)

                            if true
                                if oper1 == "BA"
                                    oper1 = "AB"
                                    vcurr = -permutedims(vcurr,[2,1,3,4])
                                elseif oper1 == "ba"
                                    oper1 = "ab"
                                    vcurr = -permutedims(vcurr,[2,1,3,4])
                                end

                                if oper2 == "BA"
                                    oper2 = "AB"
                                    vcurr = -permutedims(vcurr,[1,3,2,4])
                                elseif oper2 == "ba"
                                    oper2 = "ab"
                                    vcurr = -permutedims(vcurr,[1,3,2,4])
                                end

                                if oper3 == "BA"
                                    oper3 = "AB"
                                    vcurr = -permutedims(vcurr,[1,2,4,3])
                                elseif oper3 == "ba"
                                    oper3 = "ab"
                                    vcurr = -permutedims(vcurr,[1,2,4,3])
                                end
                            end

                            vcurr = copy(reshape(vcurr,newshape...))

                            #core,factors = tucker_decompose(vcurr)

                            parity1 = 0
                            parity2 = 0
                            parity3 = 0
                            if length(oper1)%2 != 0
                                parity1 = 1
                            end
                            if length(oper2)%2 != 0
                                parity2 = 1
                            end
                            if length(oper3)%2 != 0
                                parity3 = 1
                            end
                            parity = (parity1, parity2, parity3)

                            clusteredterm = ClusteredTerm3B((oper1,oper2,oper3), (Tuple(fock1),Tuple(fock2),Tuple(fock3)), parity, (ci, cj, ck), vcurr, Dict())
                            #display(clusteredterm)
                            focktrans = replace(zero_fock, (ci.idx, cj.idx, ck.idx), (fock1, fock2, fock3))
                            if haskey(terms,focktrans)
                                push!(terms[focktrans], clusteredterm)
                            else
                                terms[focktrans] = [clusteredterm]
                            end
                        end
                        #all(isapprox(vpqsr,0.0)) || continue
                    end

                    #=}}}=#
                end
            end
        end
    end

    # 4-body 2-electron terms
    if true
        for ci in clusters
            for cj in clusters
                for ck in clusters
                    for cl in clusters
                        #={{{=#
                        i = ci.idx
                        j = cj.idx
                        k = ck.idx
                        l = cl.idx

                        i < j < k < l|| continue

                        spin_cases =[["A","A","a","a"],
                                     ["B","B","b","b"],
                                     ["A","B","b","a"],
                                     ["B","A","a","b"]
                                    ]

                        fock_cases =[[(1,0),(1,0),(-1,0),(-1,0)],
                                     [(0,1),(0,1),(0,-1),(0,-1)],
                                     [(1,0),(0,1),(0,-1),(-1,0)],
                                     [(0,1),(1,0),(-1,0),(0,-1)]
                                    ]

                        termstr = []
                        append!(termstr,unique(permutations([ci,cj,ck,cl])))

                        #
                        #   (pr|qs) p'q'sr
                        #
                        for term in termstr

                            #
                            #   find permutations and sign needed to sort the indices
                            #   such that clusters increase from left to right
                            perm, countswap = bubble_sort(term)
                            perm == sortperm(term, alg=MergeSort)|| error("problem with bubble_sort")

                            permsign = 1
                            if countswap%2 != 0
                                permsign = -1
                            end

                            vprqs = view(ints.h2,term[1].orb_list, term[4].orb_list, term[2].orb_list, term[3].orb_list)
                            #
                            # align (prqs) ints so that they align with indices from operators before sorting
                            vpqsr = permutedims(vprqs,[1,3,4,2])

                            #
                            # now align (pqsr) ints so that they align with indices from operators after sorting
                            # in this ordering, one can simply contract by sum(v .* d)
                            v = (.5 * permsign) .* permutedims(vpqsr,perm)

                            #
                            # no reshape needed

                            for sidx in 1:length(spin_cases)
                                oper = spin_cases[sidx][perm]
                                fock = fock_cases[sidx][perm]
                                oper1 = ""
                                oper2 = ""
                                oper3 = ""
                                oper4 = ""
                                fock1 = [0,0]
                                fock2 = [0,0]
                                fock3 = [0,0]
                                fock4 = [0,0]
                                for cidx in 1:length(term[perm])
                                    if term[perm][cidx] == ci
                                        oper1 *= oper[cidx]
                                        fock1 .+= fock[cidx]
                                    elseif term[perm][cidx] == cj
                                        oper2 *= oper[cidx]
                                        fock2 .+= fock[cidx]
                                    elseif term[perm][cidx] == ck
                                        oper3 *= oper[cidx]
                                        fock3 .+= fock[cidx]
                                    elseif term[perm][cidx] == cl
                                        oper4 *= oper[cidx]
                                        fock4 .+= fock[cidx]
                                    else
                                        throw(Exception)
                                    end
                                end

                                parity1 = 0
                                parity2 = 0
                                parity3 = 0
                                parity4 = 0
                                if length(oper1)%2 != 0
                                    parity1 = 1
                                end
                                if length(oper2)%2 != 0
                                    parity2 = 1
                                end
                                if length(oper3)%2 != 0
                                    parity3 = 1
                                end
                                if length(oper4)%2 != 0
                                    parity4 = 1
                                end
                                parity = (parity1, parity2, parity3, parity4)

                                clusteredterm = ClusteredTerm4B((oper1,oper2,oper3,oper4), (Tuple(fock1),Tuple(fock2),Tuple(fock3),Tuple(fock4)), parity, (ci, cj, ck, cl), v, Dict())
                                focktrans = replace(zero_fock, (ci.idx, cj.idx, ck.idx, cl.idx), (fock1, fock2, fock3, fock4))
                                if haskey(terms,focktrans)
                                    push!(terms[focktrans], clusteredterm)
                                else
                                    terms[focktrans] = [clusteredterm]
                                end
                            end
                        end
                        #=}}}=#
                    end
                end
            end
        end
    end

    unique!(terms)

    return terms
end

"""
    extract_S2(clusters)

Form a clustered operator type for the S^2 operator
"""
function extract_S2(clusters)
            #={{{=#

    norb = 0
    for ci in clusters
        norb += length(ci)
    end

    terms = ClusteredOperator()
    n_clusters = length(clusters)
    ops_a = Array{String}(undef,n_clusters)
    ops_b = Array{String}(undef,n_clusters)
    fill!(ops_a,"")
    fill!(ops_b,"")

    zero_fock = TransferConfig([(0,0) for i in clusters])
    terms[zero_fock] = Vector{ClusteredTerm}()
    for ci in clusters
        fock1 = (0,0)
        clusteredterm = ClusteredTerm1B(("S2",), (fock1,), (0,), (ci, ), ones(1), Dict())

        focktrans = replace(zero_fock, (ci.idx,), (fock1,))
        if haskey(terms,focktrans)
            push!(terms[focktrans], clusteredterm)
        else
            terms[focktrans] = [clusteredterm]
        end

        for cj in clusters
            i = ci.idx
            j = cj.idx

            i < j || continue

            fock1 = (1,-1)
            fock2 = (-1,1)
            clusteredterm = ClusteredTerm2B(("S+","S-"), (fock1,fock2), (0,0), (ci, cj), ones(length(ci),length(cj)), Dict())

            focktrans = replace(zero_fock, (ci.idx, cj.idx), (fock1, fock2))
            if haskey(terms,focktrans)
                push!(terms[focktrans], clusteredterm)
            else
                terms[focktrans] = [clusteredterm]
            end

            fock1 = (-1,1)
            fock2 = (1,-1)
            clusteredterm = ClusteredTerm2B(("S-","S+"), (fock1,fock2), (0,0), (ci, cj), ones(length(ci),length(cj)), Dict())

            focktrans = replace(zero_fock, (ci.idx, cj.idx), (fock1, fock2))
            if haskey(terms,focktrans)
                push!(terms[focktrans], clusteredterm)
            else
                terms[focktrans] = [clusteredterm]
            end

            fock1 = (0,0)
            fock2 = (0,0)
            clusteredterm = ClusteredTerm2B(("Sz","Sz"), (fock1,fock2), (0,0), (ci, cj), 2*ones(1,1), Dict())

            focktrans = replace(zero_fock, (ci.idx, cj.idx), (fock1, fock2))
            if haskey(terms,focktrans)
                push!(terms[focktrans], clusteredterm)
            else
                terms[focktrans] = [clusteredterm]
            end
        end
    end
    return terms
end
#=}}}=#

"""
    unique!(clustered_ham::ClusteredOperator)

combine terms to keep only unique operators
"""
function unique!(clustered_ham::ClusteredOperator)
#function unique!(clustered_ham::Dict{TransferConfig,Vector{ClusteredTerm}})
#={{{=#
    println(" Remove duplicates")
    #
    # first just remove duplicates
    nstart = 0
    nfinal = 0
#    for (ftrans, terms) in clustered_ham
#        tmp = deepcopy(terms)
#        for term in terms
#            println(term.ops)
##            swap = false
##            idx = 0
##            for (opidx,op) in enumerate(term.ops)
##                if op == "BA"
##                    swap = true
##                    idx = opidx
##                end
##            end
##            if swap
##                println(size(term.ints))
##                v = term.ints
##                v = copy(reshape(v, (size(v,1), size(v,2), size(v,3), size(v,4))))
##                v = -permutedims(v,[2,1,3,4])
##                v = reshape(v, (size(v,1)*size(v,2)*size(v,3), size(v,4)))
##
##                newterm = ClusteredTerm2B(("ABb","a"), term.delta, term.clusters, term.ints)
##                push!(tmp,newterm)
##            end
#
#            if term.ops == ("BAb","a")
#
#                println(size(term.ints))
#                v = term.ints
#                v = copy(reshape(v, (size(v,1), size(v,2), size(v,3), size(v,4))))
#                v = -permutedims(v,[2,1,3,4])
#                v = reshape(v, (size(v,1)*size(v,2)*size(v,3), size(v,4)))
#
#                newterm = ClusteredTerm2B(("ABb","a"), term.delta, term.clusters, term.ints)
#                push!(tmp,newterm)
#            else
#                push!(tmp,term)
#            end
#        end
#
#        clustered_ham[ftrans] = tmp
#    end

#        tmp = deepcopy(terms)
#        for term in tmp
#            for op in term.ops
#                if op == "BAb"
#
#                    fconfig =
#                    newterm = ClusteredTerm2B((oper1,oper2,oper3,oper4), ftrans, term.clusters, v)

    for (ftrans, terms) in clustered_ham
        unique = Dict()
        for term in terms
            nstart += 1
            keystr = ""
            for (i,j) in zip(term.ops,term.clusters)
                keystr *= string(i,"(",j.idx,")")
            end
            if haskey(unique,keystr)
                unique[keystr].ints .+= term.ints
            else
                unique[keystr] = deepcopy(term)
            end
        end
        clustered_ham[ftrans] = Vector{ClusteredTerm}()
        for (keystr, term) in unique
            push!(clustered_ham[ftrans],term)
        end
        nfinal += length(clustered_ham[ftrans])
    end

    @printf(" Number of terms reduced from %5i to %5i\n", nstart, nfinal)
#=}}}=#
end

function check_term(term::ClusteredTerm1B,
                            fock_bra::FockConfig, bra::T,
                            fock_ket::FockConfig, ket::T) where T<:Union{ClusterConfig, TuckerConfig}
    length(fock_bra) == length(fock_ket) || throw(Exception)
    length(bra) == length(ket) || throw(Exception)
    n_clusters = length(bra)
    #
    # make sure inactive clusters are diagonal
    for ci in 1:n_clusters
        ci != term.clusters[1].idx || continue

        fock_bra[ci] == fock_ket[ci] || return false
        bra[ci] == ket[ci] || return false
    end
    fock_bra[term.clusters[1].idx] == fock_ket[term.clusters[1].idx] .+ term.delta[1] || return false
    return true
end

function check_term(term::ClusteredTerm2B,
                            fock_bra::FockConfig, bra::T,
                            fock_ket::FockConfig, ket::T) where T<:Union{ClusterConfig, TuckerConfig}
    length(fock_bra) == length(fock_ket) || throw(Exception)
    length(bra) == length(ket) || throw(Exception)
    n_clusters = length(bra)
    #
    #
    #
    # make sure inactive clusters are diagonal
    for ci in 1:n_clusters
        ci != term.clusters[1].idx || continue
        ci != term.clusters[2].idx || continue

        fock_bra[ci] == fock_ket[ci] || return false
        bra[ci] == ket[ci] || return false
    end
    fock_bra[term.clusters[1].idx] == fock_ket[term.clusters[1].idx] .+ term.delta[1] || return false
    fock_bra[term.clusters[2].idx] == fock_ket[term.clusters[2].idx] .+ term.delta[2] || return false
    return true
end

function check_term(term::ClusteredTerm3B,
                            fock_bra::FockConfig, bra::T,
                            fock_ket::FockConfig, ket::T) where T<:Union{ClusterConfig, TuckerConfig}
    length(fock_bra) == length(fock_ket) || throw(Exception)
    length(bra) == length(ket) || throw(Exception)
    n_clusters = length(bra)
    #
    # make sure inactive clusters are diagonal
    for ci in 1:n_clusters
        ci != term.clusters[1].idx || continue
        ci != term.clusters[2].idx || continue
        ci != term.clusters[3].idx || continue

        fock_bra[ci] == fock_ket[ci] || return false
        bra[ci] == ket[ci] || return false
    end
    fock_bra[term.clusters[1].idx] == fock_ket[term.clusters[1].idx] .+ term.delta[1] || return false
    fock_bra[term.clusters[2].idx] == fock_ket[term.clusters[2].idx] .+ term.delta[2] || return false
    fock_bra[term.clusters[3].idx] == fock_ket[term.clusters[3].idx] .+ term.delta[3] || return false
    return true
end

function check_term(term::ClusteredTerm4B,
                            fock_bra::FockConfig, bra::T,
                            fock_ket::FockConfig, ket::T) where T<:Union{ClusterConfig, TuckerConfig}
    length(fock_bra) == length(fock_ket) || throw(Exception)
    length(bra) == length(ket) || throw(Exception)
    n_clusters = length(bra)
    #
    # make sure inactive clusters are diagonal
    for ci in 1:n_clusters
        ci != term.clusters[1].idx || continue
        ci != term.clusters[2].idx || continue
        ci != term.clusters[3].idx || continue
        ci != term.clusters[4].idx || continue

        fock_bra[ci] == fock_ket[ci] || return false
        bra[ci] == ket[ci] || return false
    end
    fock_bra[term.clusters[1].idx] == fock_ket[term.clusters[1].idx] .+ term.delta[1] || return false
    fock_bra[term.clusters[2].idx] == fock_ket[term.clusters[2].idx] .+ term.delta[2] || return false
    fock_bra[term.clusters[3].idx] == fock_ket[term.clusters[3].idx] .+ term.delta[3] || return false
    fock_bra[term.clusters[4].idx] == fock_ket[term.clusters[4].idx] .+ term.delta[4] || return false
    return true
end

function compute_terms_state_sign(term::ClusteredTerm, fock_ket::FockConfig)
    #
    # determine sign from rearranging clusters if odd number of operators
    state_sign = 1
    for (oi,o) in enumerate(term.ops)
        if term.parity[oi] == 1  #only count electrons if operator is odd
            n_elec_hopped = 0
            for ci in 1:term.clusters[oi].idx-1
                n_elec_hopped += fock_ket[ci][1] + fock_ket[ci][2]
            end
            if n_elec_hopped % 2 != 0
                state_sign = -state_sign
            end
        end
    end
    return state_sign
end

function print_fock_sectors(sector::Vector{Tuple{T,T}}) where T<:Integer
    print("  ")
    for ci in sector
        @printf("(%iα,%iβ)", ci[1],ci[2])
    end
    println()
end

"""
    extract_1body_operator(clustered_ham::ClusteredOperator; op_string="H")

Extract a 1-body operator for use in perturbation theory
- `op_string`: either H or Hcmf
"""
function extract_1body_operator(clustered_ham::ClusteredOperator; op_string="H")
    out = ClusteredOperator()
    for (ftrans, terms) in clustered_ham
        for term in terms
            if term isa ClusteredTerm1B

                term2 = ClusteredTerm1B((op_string,), term.delta, term.parity, term.clusters, term.ints, term.cache)
                if haskey(out, ftrans)
                    push!(out[ftrans], term2)
                else
                    out[ftrans] = [term2]
                end
            end
        end
    end
    return out
end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
using Test
import Julie

@testset "julie_test" begin
    @test #Something
end

max_stars_repo_name: ErikQQY/Julie

max_issues_repo_name: ErikQQY/Julie

max_forks_repo_name: ErikQQY/Julie


text:
function f_bad()
    x = 0
    Threads.@threads for i in 1:10^7
        x += 1
    end
    return x
end

function f_atomic()
    x = Threads.Atomic{Int}(0)
    Threads.@threads for i in 1:10^7
        Threads.atomic_add!(x, 1)
    end
    return x[]
end

function f_spin()
    l = Threads.SpinLock()
    x = 0
    Threads.@threads for i in 1:10^7
        Threads.lock(l)
        x += 1
        Threads.unlock(l)
    end
    return x
end

function f_mutex()
    l = Threads.Mutex()
    x = 0
    Threads.@threads for i in 1:10^7
        Threads.lock(l)
        x += 1
        Threads.unlock(l)
    end
    return x
end

for f in [f_bad, f_atomic, f_spin, f_mutex]
    println(f)
    @time println(f())
    @time println(f())
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
export Parameters

immutable Parameters
    mu::Integer
    lambda::Integer
    mutrate::AbstractFloat
    targetfitness::AbstractFloat

    numinputs::Integer
    numoutputs::Integer

    numperlevel::Integer
    numlevels::Integer
    numlevelsback::Integer

    funcs::Vector{Func}
    fitfunc::Function
end

function Parameters(numinputs, numoutputs, numperlevel, numlevels, numlevelsback)
    mu = 1
    lambda = 4
    mutrate = 0.05
    targetfitness = 1.0
    funcs = default_funcs()
    fitfunc = hamming_max

    return Parameters(mu, lambda, mutrate, targetfitness, numinputs, numoutputs, numperlevel, numlevels, numlevelsback, funcs, fitfunc)
end


text:
#=
Part of this work is derived from astropy/photutils. The relevant derivations
are considered under a BSD 3-clause license. =#

"""
    RectangularAperture(x, y, w, h, θ=0)
    RectangularAperture(position, w, h, θ=0)

A rectangular aperture.

A rectangular aperture with width `w`, height `h`, and position angle `θ` in degrees.

# Examples
```jldoctest
julia> ap = RectangularAperture(0, 0, 10, 4, 0)
5×11 RectangularAperture{Int64} with indices -2:2×-5:5:
 0.25  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.25
 0.5   1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.5
 0.5   1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.5
 0.5   1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  1.0  0.5
 0.25  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.25
```
"""
struct RectangularAperture{T <: Number} <: AbstractAperture{T}
    x::T
    y::T
    w::T
    h::T
    theta::T
    function RectangularAperture(x::T, y::T, w::T, h::T, θ::T) where T <: Number
        (w < 0  || h < 0) && error("Invalid side lengths w=$w, h=$h. Both must be positive.")
        new{T}(x, y, w, h, mod(θ, 360))
    end
end

RectangularAperture(center, w, h, θ) = RectangularAperture(center..., w, h, θ)
RectangularAperture(x, y, w, h, θ) = RectangularAperture(promote(x, y, w, h, θ)...)

function Base.show(io::IO, ap::RectangularAperture)
    print(io, "RectangularAperture($(ap.x), $(ap.y), w=$(ap.w), h=$(ap.h), θ=$(ap.theta)°)")
end

function bounds(ap::RectangularAperture)
    w2 = ap.w / 2
    h2 = ap.h / 2
    sint, cost = sincos(deg2rad(ap.theta))

    dx1 = abs(w2 * cost - h2 * sint)
    dy1 = abs(w2 * sint + h2 * cost)
    dx2 = abs(w2 * cost + h2 * sint)
    dy2 = abs(w2 * sint - h2 * cost)

    dx = max(dx1, dx2)
    dy = max(dy1, dy2)

    xmin = ceil(Int, ap.x - dx - 0.5)
    ymin = ceil(Int, ap.y - dy - 0.5)
    xmax = ceil(Int, ap.x + dx - 0.5)
    ymax = ceil(Int, ap.y + dy - 0.5)
    return xmin, xmax, ymin, ymax
end

function overlap(ap::RectangularAperture, i, j)
    y = i - ap.y
    x = j - ap.x
    flags = (
        inside_rectangle(x - 0.5, y - 0.5, ap.w, ap.h, ap.theta),
        inside_rectangle(x - 0.5, y + 0.5, ap.w, ap.h, ap.theta),
        inside_rectangle(x + 0.5, y - 0.5, ap.w, ap.h, ap.theta),
        inside_rectangle(x + 0.5, y + 0.5, ap.w, ap.h, ap.theta)
    )
    all(flags) && return Inside
    all(!, flags) && return Outside

    return Partial
end

partial(ap::RectangularAperture, x, y) = rectangular_overlap_exact(x - 0.5, y - 0.5, x + 0.5, y + 0.5, ap.w, ap.h, ap.theta)
partial(sub_ap::Subpixel{T,<:RectangularAperture}, x, y) where {T} = rectangular_overlap_single_subpixel(x - 0.5, y - 0.5, x + 0.5, y + 0.5, sub_ap.ap.w, sub_ap.ap.h, sub_ap.ap.theta, sub_ap.N)

#######################################################

"""
    RectangularAnnulus(x, y, w_in, w_out, h_out, θ=0)
    RectangularAnnulus(position, w_in, w_out, h_out, θ=0)

A rectangular annulus with inner width `w_in`, outer width `w_out`, outer height `h_out`, and position angle `θ` in degrees. `h_in` is automatically calculated from `w_in / w_out * h_out`. Note that `w_out ≥ w_in > 0`.

# Examples
```jldoctest
julia> ap = RectangularAnnulus(0, 0, 5, 10, 8, 45)
13×13 RectangularAnnulus{Float64} with indices -6:6×-6:6:
 0.0       0.0       0.0       0.0       …  0.0       0.0       0.0
 0.0       0.0       0.0       0.0          0.0       0.0       0.0
 0.0       0.0       0.0       0.568542     0.0       0.0       0.0
 0.0       0.0       0.568542  1.0          0.0       0.0       0.0
 0.0       0.568542  1.0       1.0          0.5       0.0       0.0
 0.528175  1.0       1.0       1.0       …  0.941125  0.0       0.0
 0.0       0.941125  1.0       0.678301     1.0       0.941125  0.0
 0.0       0.0       0.941125  1.0          1.0       1.0       0.528175
 0.0       0.0       0.5       0.941125     1.0       0.568542  0.0
 0.0       0.0       0.0       0.5          0.568542  0.0       0.0
 0.0       0.0       0.0       0.0       …  0.0       0.0       0.0
 0.0       0.0       0.0       0.0          0.0       0.0       0.0
 0.0       0.0       0.0       0.0          0.0       0.0       0.0
```
"""
struct RectangularAnnulus{T <: Number} <: AbstractAperture{T}
    x::T
    y::T
    w_in::T
    w_out::T
    h_in::T
    h_out::T
    theta::T
end

RectangularAnnulus(center, w_in, w_out, h_out, θ=0) = RectangularAnnulus(center..., w_in, w_out, h_out, θ)
RectangularAnnulus(x::Number, y::Number, w_in, w_out, h_out, θ=0) = RectangularAnnulus(promote(x, y, w_in, w_out, h_out * w_in / w_out, h_out, θ)...)
# RectangularAnnulus(x, y, w_in, w_out, h_out) = RectangularAnnulus(promote(x, y, w_in, w_out, h_out, 0)...)

function Base.show(io::IO, ap::RectangularAnnulus)
    print(io, "RectangularAnnulus($(ap.x), $(ap.y), w_in=$(ap.w_in), w_out=$(ap.w_out), h_in=$(ap.h_in), h_out=$(ap.h_out), θ=$(ap.theta)°)")
end

function overlap(ap::RectangularAnnulus, i, j)
    y = i - ap.y
    x = j - ap.x
    flags_out = (
        inside_rectangle(x - 0.5, y - 0.5, ap.w_out, ap.h_out, ap.theta),
        inside_rectangle(x - 0.5, y + 0.5, ap.w_out, ap.h_out, ap.theta),
        inside_rectangle(x + 0.5, y - 0.5, ap.w_out, ap.h_out, ap.theta),
        inside_rectangle(x + 0.5, y + 0.5, ap.w_out, ap.h_out, ap.theta)
    )

    flags_in = (
        inside_rectangle(x - 0.5, y - 0.5, ap.w_in, ap.h_in, ap.theta),
        inside_rectangle(x - 0.5, y + 0.5, ap.w_in, ap.h_in, ap.theta),
        inside_rectangle(x + 0.5, y - 0.5, ap.w_in, ap.h_in, ap.theta),
        inside_rectangle(x + 0.5, y + 0.5, ap.w_in, ap.h_in, ap.theta)
    )

   all(flags_out) && all(!, flags_in) && return Inside
   all(flags_in) || all(!, flags_out) && return Outside

    return Partial
end

function bounds(ap::RectangularAnnulus)
    w2 = ap.w_out / 2
    h2 = ap.h_out / 2
    sint, cost = sincos(deg2rad(ap.theta))

    dx1 = abs(w2 * cost - h2 * sint)
    dy1 = abs(w2 * sint + h2 * cost)
    dx2 = abs(w2 * cost + h2 * sint)
    dy2 = abs(w2 * sint - h2 * cost)

    dx = max(dx1, dx2)
    dy = max(dy1, dy2)

    xmin = ceil(Int, ap.x - dx - 0.5)
    ymin = ceil(Int, ap.y - dy - 0.5)
    xmax = ceil(Int, ap.x + dx - 0.5)
    ymax = ceil(Int, ap.y + dy - 0.5)
    return (xmin, xmax, ymin, ymax)
end

function partial(ap::RectangularAnnulus, x, y)
    f1 = rectangular_overlap_exact(x - 0.5, y - 0.5, x + 0.5, y + 0.5, ap.w_out, ap.h_out, ap.theta)
    f2 = rectangular_overlap_exact(x - 0.5, y - 0.5, x + 0.5, y + 0.5, ap.w_in, ap.h_in, ap.theta)
    return f1 - f2
end

function partial(sub_ap::Subpixel{T,<:RectangularAnnulus}, x, y) where T
    ap = sub_ap.ap
    f1 = rectangular_overlap_single_subpixel(x - 0.5, y - 0.5, x + 0.5, y + 0.5, ap.w_out, ap.h_out, ap.theta, sub_ap.N)
    f2 = rectangular_overlap_single_subpixel(x - 0.5, y - 0.5, x + 0.5, y + 0.5, ap.w_in, ap.h_in, ap.theta, sub_ap.N)
    return f1 - f2
end


text:
import PAPI

const NUM_EVENTS = 2
const THRESHOLD = 10_000

function computation_mult()
    tmp = 1.0
    for i=1:THRESHOLD
        tmp *= i
    end
    return tmp
end

function computation_add()
    tmp = 0
    for i=1:THRESHOLD
        tmp += i
    end
    return tmp
end

function main()

    # Precompile functions so that we are also not
    # measuring compilation overhead
    precompile(computation_add, ())
    precompile(computation_mult,())

    cs = PAPI.EventSet([PAPI.TOT_INS, PAPI.TOT_CYC])

    info("There are $(PAPI.num_counters()) counters on this system")

    # Initialize the PAPI library and start counting
    # the events named in the events array.
    # This function implicitly stops and initializes
    # any counters running as a result of a previous call
    # to PAPI.start_counters()

    info("Counters started")
    PAPI.start_counters(cs)

    computation_add()

    values = PAPI.read_counters!(cs)

    @printf("The total instructions executed for addition are %lld \n",values[1]);
    @printf("The total cycles used are %lld \n", values[2] );
end

main()


text:
instructions = String[]
global line = chomp(readline());
while line != ""
    global line
    push!(instructions, line)
    global line = chomp(readline());
end

called = Int[]
global inst = 0
global accumulator = 0
while !(inst in called)
    push!(called, inst)
    elems = split(instructions[inst + 1])
    op = elems[1]
    amount = parse(Int, elems[2])
    if op == "acc"
        global accumulator += amount
    end

    if op == "jmp"
        global inst += amount
    else
        global inst += 1
    end
end

println(accumulator)

max_stars_repo_name: noeppi-noeppi/aoc

max_issues_repo_name: noeppi-noeppi/aoc

max_forks_repo_name: noeppi-noeppi/aoc


text:
# Simple test of fault tolerance functions

using Test
using Main.CGP
MyFunc = Main.CGP.MyFunc
funcs = default_funcs(3)
Not = Main.CGP.Not  # Not() only inverts the relevant bits
c = build_chromosome(
  (1,2,3),
  ((NAND,Integer[3, 1]),(NOR,Integer[2, 1]),(OR,Integer[4, 5])),
  (5,6))

context = construct_context(c.params.numinputs)
outv = output_values(c)

@testset "test correctness of execute_chromomsome_ft() and fault_tolerance_fitness()" begin
# Compute output for c
n1 = context[1]
n2 = context[2]
n3 = context[3]
n4 = Nand(n3,n1)
n5 = Nor(n2,n1)
n6 = n4 | n5
n7 = n5
n8 = n6
@test outv == [n7,n8]

# compute output for c with gate 6 negated
p6 = Not(n6)
p7 = n5
p8 = p6
e6 = execute_chromosome_ft( c, context, 6 );
@test  e6 == [p7,p8]

# compute output for c with gate 5 negated
q5 = Not(n5)
q6 = n4 | q5
q7 = q5
q8 = q6
e5 = execute_chromosome_ft( c, context, 5 );
@test  e5 == [q7,q8]

# compute output for c with gate 4 negated
r4 = Not(n4)
r5 = Nor(n2,n1)
r6 = r4 | r5
r7 = r5
r8 = r6
e4 = execute_chromosome_ft( c, context, 4 );
@test  e4 == [r7,r8]

# compute output for c with gate 3 negated
s3 = Not(n3)
s4 = Nand(s3,n1)
s5 = Nor(n2,n1)
s6 = s4 | s5
s7 = s5
s8 = s6
e3 = execute_chromosome_ft( c, context, 3 );
@test  e3 == [s7,s8]

funcs = default_funcs(2)  # Reset for chromosome c1 which has 2 inputs
# Flipping the output of gate 3 of c1 changes the output from 0x00 to 0x0f
c1 = build_chromosome( (1,2), ((ZERO,Integer[]),), (3,))
println("fault_tolerance_fitness(c1): ",fault_tolerance_fitness(c1))
@test isapprox(fault_tolerance_fitness(c1),0.0)

funcs = default_funcs(3)  # Reset for chromosome nulc which has 3 inputs
# Flipping the output of gate 4 (the first interior node) has no effect on the output
# Flipping the output of gate 5 (the second interior node)  flips the first output but not the second
# Flipping the output of gate 6 (the third interior node)  flips the second output but not the first
nulc = build_chromosome(
  (1,2,3),
  ((NAND,Integer[3, 1]),(ZERO,Integer[]),(ONE,Integer[])),
  (5,6))

@test isapprox(fault_tolerance_fitness(nulc),2/3)
end  # @testset


text:
# Add manual additions to the emoji table here, they will take precedence over
# the automatic entries
# It is a vector of Pair{String,String}

const manual = []


text:
function iliev(k)
#
#  generalization of Iliev example:
#
#    (x-1)^k (x-2)^2k (x-3)^4k
#
    p = poly([ones(k);2ones(2k);3ones(3k)])
    z = [1.0 k; 2 2k; 3 4k]
    p, PolyZeros(z)
end


text:
if !isfile(joinpath(@__DIR__, "already_showed"))
    print_with_color(Base.info_color(), STDERR,
    """
    PkgBenchmark has been completely rewritten. Please see https://github.com/JuliaCI/PkgBenchmark.jl/
    for updated documentation and examples. Code written for previous versions of PkgBenchmark is
    unlikely to still work.
    """)
    touch("already_showed")
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_static.jl --help` to see a usage message.
using BinaryBuilder, Pkg

include("config.jl")

# Collection of sources required to complete build
sources = [ArchiveSource(
    "https://github.com/chemfiles/chemfiles/archive/$version.tar.gz", chemfiles_sha256
)]

# Bash recipe for building across all platforms
script = """
cd \${WORKSPACE}/srcdir/chemfiles-*/
mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=\${prefix} \\
      -DCMAKE_TOOLCHAIN_FILE=\${CMAKE_TARGET_TOOLCHAIN} \\
      -DCMAKE_POSITION_INDEPENDENT_CODE=ON \\
      -DCMAKE_BUILD_TYPE=Release \\
      -DBUILD_SHARED_LIBS=$(BUILD_SHARED_LIBS ? "ON" : "OFF") \\
      ..
make -j\${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
# Also build for M1 macs
push!(platforms, Platform("aarch64", "macos"))

# Build the tarballs
build_tarballs(
    ARGS,
    name,
    version,
    sources,
    script,
    platforms,
    Product[],
    Dependency[];
    verbose=true,
    julia_compat="1.6"
)

for path in readdir("products/")
    if BUILD_SHARED_LIBS && startswith(path, "Chemfiles.")
        new_path = replace(path, "Chemfiles." => "chemfiles-shared.")
    elseif !BUILD_SHARED_LIBS && startswith(path, "Chemfiles.")
        new_path = replace(path, "Chemfiles." => "chemfiles-static.")
    else
        continue
    end
    Base.Filesystem.mv(joinpath("products/", path), joinpath("products/", new_path))
end

println("====================================================\n")
println("Done builing chemfiles for all platforms,")
println("make sure to upload everything in products/")

max_stars_repo_name: chemfiles/chemfiles-prebuilt

max_issues_repo_name: chemfiles/chemfiles-prebuilt

max_forks_repo_name: chemfiles/chemfiles-prebuilt


text:
abstract type AbstractTrajectory{n,m,T} <: AbstractVector{T} end

has_terminal_control(Z::AbstractTrajectory) = !RobotDynamics.is_terminal(Z[end])
traj_size(Z::AbstractTrajectory{n,m}) where {n,m} = n,m,length(Z)
num_vars(Z::AbstractTrajectory) = num_vars(traj_size(Z)..., has_terminal_control(Z))
eachcontrol(Z::AbstractTrajectory) = has_terminal_control(Z) ? Base.OneTo(length(Z)) : Base.OneTo(length(Z)-1)
state_dim(Z::AbstractTrajectory{n}) where {n} = n
control_dim(Z::AbstractTrajectory{<:Any,m}) where {m} = m

function num_vars(n::Int, m::Int, N::Int, equal::Bool=false)
    Nu = equal ? N : N-1
    return N*n + Nu*m
end

function Base.copyto!(dest::AbstractTrajectory, src::AbstractTrajectory)
	@assert traj_size(dest) == traj_size(src)
	for k = 1:length(src)
		dest[k] = src[k]
	end
	return src
end

@inline states(Z::AbstractTrajectory) = state.(Z)
function controls(Z::AbstractTrajectory)
	return [control(Z[k]) for k in eachcontrol(Z) ]
end

states(Z::AbstractTrajectory, inds::AbstractVector{<:Integer}) = [states(Z, i) for i in inds]
controls(Z::AbstractTrajectory, inds::AbstractVector{<:Integer}) = [controls(Z, i) for i in inds]

states(Z::AbstractTrajectory, ind::Integer) = [state(z)[ind] for z in Z]
controls(Z::AbstractTrajectory, ind::Integer) = [control(z)[ind] for z in Z]

get_data(Z::AbstractTrajectory) = get_data.(Z)

function Base.isapprox(Z1::AbstractTrajectory, Z2::AbstractTrajectory)
    all(zs->zs[1] ≈ zs[2], zip(Z1,Z2))
end

function set_dt!(Z::AbstractTrajectory, dt::Real)
    t = Z[1].t
    for z in Z
        z.dt = dt
        z.t = t
        t += dt
    end
    return t
end

"""
    Traj{n,m,T,KP}

A vector of `AbstractKnotPoint`s of type `KP` with state dimension `n`,
control dimension `m`, and value type `T`

Supports iteration and indexing.

# Constructors
    Traj(n, m, dt, N, equal=false)
    Traj(x, u, dt, N, equal=false)
    Traj(X, U, dt, t)
    Traj(X, U, dt)
"""
struct Traj{n,m,T,KP} <: AbstractTrajectory{n,m,T}
	data::Vector{KP}
	function Traj(Z::Vector{<:AbstractKnotPoint{n,m,<:Any,T}}) where {n,m,T}
		new{n,m,T,eltype(Z)}(Z)
	end
end

# AbstractArray interface
@inline Base.iterate(Z::Traj, k::Int) = iterate(Z.data, k)
Base.IteratorSize(Z::Traj) = Base.HasLength()
Base.IteratorEltype(Z::Traj) = Base.IteratorEltype(Z.data)
@inline Base.eltype(Z::Traj) = eltype(Z.data)
@inline Base.length(Z::Traj) = length(Z.data)
@inline Base.size(Z::Traj) = size(Z.data)
@inline Base.getindex(Z::Traj, i) = Z.data[i]
@inline Base.setindex!(Z::Traj, v, i) = Z.data[i] = v
@inline Base.firstindex(Z::Traj) = 1
@inline Base.lastindex(Z::Traj) = lastindex(Z.data)
Base.IndexStyle(::Traj) = IndexLinear()

Traj(Z::Traj) = Z

function Base.copy(Z::AbstractTrajectory{Nx,Nu}) where {Nx,Nu}
    Traj([KnotPoint{Nx,Nu}(copy(z.z), z.t, z.dt) for z in Z])
end

function Traj(n::Int, m::Int, dt::AbstractFloat, N::Int; equal=false)
    x = NaN*@SVector ones(n)
    u = @SVector zeros(m)
    Traj(x,u,dt,N; equal=equal)
end

function Traj(x::SVector, u::SVector, dt::AbstractFloat, N::Int; equal=false)
    equal ? uN = N : uN = N-1
    Z = [KnotPoint(x,u,(k-1)*dt,dt) for k = 1:uN]
    if !equal
        m = length(u)
        push!(Z, KnotPoint(x,u*0,(N-1)*dt,0.0))
    end
    return Traj(Z)
end

function Traj(X::Vector, U::Vector, dt::Vector, t=cumsum(dt) .- dt[1])
    n,m = length(X[1]), length(U[1])
    Z = [KnotPoint{n,m}(n,m,[X[k]; U[k]], t[k], dt[k]) for k = 1:length(U)]
    if length(U) == length(X)-1
        push!(Z, KnotPoint{n,m}(n,m,[X[end]; U[1]*0],t[end],0.0))
    end
    return Traj(Z)
end

function setstates!(Z::Traj, X)
    for k in eachindex(Z)
		setstate!(Z[k], X[k])
    end
end

function setstates!(Z::Traj, X::AbstractMatrix)
    for k in eachindex(Z)
		setstate!(Z[k], X[:,k])
    end
end

function setcontrols!(Z::AbstractTrajectory, U)
    for k in 1:length(Z)-1
		setcontrol!(Z[k], U[k])
    end
end

function setcontrols!(Z::AbstractTrajectory, U::AbstractMatrix)
    for k in 1:length(Z)-1
		setcontrol!(Z[k], U[:,k])
    end
end

function setcontrols!(Z::AbstractTrajectory, u::SVector)
    for k in 1:length(Z)-1
		setcontrol!(Z[k], u)
    end
end

function settimes!(Z::AbstractTrajectory, ts)
    for k in eachindex(ts)
        Z[k].t = ts[k]
        k < length(ts) && (Z[k].dt = ts[k+1] - ts[k])
    end
end

function gettimes(Z::Traj)
    [z.t for z in Z]
end

function shift_fill!(Z::Traj, n=1)
    N = length(Z)
    isterm = is_terminal(Z[end])
    for k = 1+n:N
        Z[k-n] = copy(Z[k])
    end
    xf = state(Z[N-n])
    uf = control(Z[N-n])
    dt = Z[N-n-1].dt
    for k = N-n:N
        setstate!(Z[k], xf)
        Z[k].t = Z[k-1].t + dt
        if k == N && is_terminal(Z[k])
            Z[k].dt = 0
        else
            setcontrol!(Z[k], uf)
            Z[k].dt = dt
        end
    end
end

function Base.copyto!(Z::Traj, Z0::Traj)
	@assert length(Z) == length(Z0)
	for k in eachindex(Z)
		copyto!(Z[k].z, Z0[k].z)
	end
end

function Base.copyto!(Z::Union{Vector{<:KnotPoint},Traj{<:Any,<:Any,<:Any,<:KnotPoint}}, Z0::Traj)
	@assert length(Z) == length(Z0)
	for k in eachindex(Z)
		Z[k].z = Z0[k].z
	end
end

#~~~~~~~~~~~~~~~~~~~~~~~~~~ FUNCTIONS ON TRAJECTORIES ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

state_diff_jacobian!(model::AbstractModel, G, Z::AbstractTrajectory) =
    state_diff_jacobian!(statevectortype(model), model, G, Z)
state_diff_jacobian!(::EuclideanState, model::AbstractModel, G, Z::AbstractTrajectory) = nothing

function state_diff_jacobian!(::RotationState, model::AbstractModel, G, Z::AbstractTrajectory)
	for k in eachindex(Z)
		G[k] .= 0
		state_diff_jacobian!(RotationState(), model, G[k], Z[k])
	end
end

function rollout!(sig::FunctionSignature, model::DiscreteDynamics, Z::AbstractTrajectory, x0=state(Z[1]))
    setstate!(Z[1], x0)
    for k = 2:length(Z)
        RobotDynamics.propagate_dynamics!(sig, model, Z[k], Z[k-1])
    end
end


text:
using Gadfly
using Cairo
using DataFrames

K2eV = 8.61733e-5 #K->eV

plottheme = Gadfly.Theme(
    highlight_width=0pt,
    default_point_size=2pt,
    line_width=2pt,
    minor_label_font_size=13pt,
    minor_label_color=colorant"black",
    major_label_font_size=15pt,
    major_label_color=colorant"black",
    key_label_color=colorant"black",
    key_title_color=colorant"black",
    key_title_font_size=12pt,
    key_label_font_size=10pt,
    grid_color=colorant"black")

function b2(e,t)
  return 1.0e12*exp(-e/t)
end

function τ(e,t)
  rn = rand(1)
  return (-1.0/b2(e,t))*log(rn[1])
end

df100 = DataFrame()
df100[:temp] = 1:1:1.0e3
df100[:time] = 0.0
df100[:ed] = 100.0
df100[:str] = "0.008"

df200 = DataFrame()
df200[:temp] = 1:1:1.0e3
df200[:time] = 0.0
df200[:ed] = 200.0
df200[:str] = "0.01"

df300 = DataFrame()
df300[:temp] = 1:1:1.0e3
df300[:time] = 0.0
df300[:ed] = 300.0
df300[:str] = "0.02"

df400 = DataFrame()
df400[:temp] = 1:1:1.0e3
df400[:time] = 0.0
df400[:ed] = 400.0
df400[:str] = "0.03"

df500 = DataFrame()
df500[:temp] = 1:1:1.0e3
df500[:time] = 0.0
df500[:ed] = 500.0
df500[:str] = "0.04"

df = DataFrame()
df = vcat(df, df100, df200, df300, df400, df500)

for i in 1:size(df[:temp],1)
  if i%1000 == 0 println("$i of $(size(df,1))") end
  df[:time][i] = τ(df[:ed][i],df[:temp][i])
  df[:ed][i] = df[:ed][i]*K2eV
end

p1 = plot(
          df,
          x="temp",
          y="time",
          color="str",
          Geom.smooth,
          Scale.color_discrete,
          Scale.y_log10,
          Scale.x_log10,
          Guide.xlabel("Temperature (K)"),
          Guide.ylabel("Hopping time (s)"),
          Guide.colorkey("Eb(eV)"),
         )
push!(p1,plottheme)

draw(PDF("f0.pdf",6inch,4.5inch),p1)

println("Ending script")


text:
# Orthogonality tests

using GaussQuadrature
using Printf

T = BigFloat
println("\nFloating point type data type is ", T)
println("\teps = ", eps(T))

xlegendre(n, alpha, beta, endpt) = legendre(T, n, endpt)
xchebyshev_first(n, alpha, beta, endpt) = chebyshev(T, n, 1, endpt)
xchebyshev_second(n, alpha, beta, endpt) = chebyshev(T, n, 2, endpt)
xhermite(n, alpha, beta, endpt) = hermite(T, n)
xlaguerre(n, alpha, beta, endpt) = laguerre(n, alpha, endpt)
xlogweight(n, alpha, beta , endpt) = logweight(T, n,
                                     round(Int64,alpha), endpt)

xlegendre_coefs(n, alpha, beta) = legendre_coefs(T, n)
xchebyshev_first_coefs(n, alpha, beta) = chebyshev_coefs(T, n, 1)
xchebyshev_second_coefs(n, alpha, beta) = chebyshev_coefs(T, n, 2)
xhermite_coefs(n, alpha, beta) = hermite_coefs(T, n)
xlaguerre_coefs(n, alpha, beta) = laguerre_coefs(n, alpha)
xlogweight_coefs(n, alpha, beta) = logweight_coefs(T, n,
                                   round(Int64,alpha))

rule = [xlegendre, xchebyshev_first, xchebyshev_second,
        xhermite, jacobi, xlaguerre, xlogweight]

coefs = [ xlegendre_coefs, xchebyshev_first_coefs,
          xchebyshev_second_coefs, xhermite_coefs, jacobi_coefs,
          xlaguerre_coefs, xlogweight_coefs ]

name = ["Legendre", "Chebyshev (first)", "Chebyshev (second)",
         "Hermite", "Jacobi", "Laguerre", "Log weight" ]

alpha = one(T)
beta  = one(T)/2
xargs = [(),(1),(2), (), (alpha, beta), (alpha)]

function discrepancy(n, dop, w, p)
    d = 0.0
    for k = 0:dop
        for j = 0:min(k-1, dop-k)
            s = 0.0
            for l = 1:n
                s += w[l] * p[l,j+1] * p[l,k+1]
            end
            d = max(d, abs(s))
        end
        if 2*k <= dop
            s = 0.0
            for l = 1:n
                s += w[l] * p[l,k+1]^2
            end
            d = max(d, abs(s-1.0))
        end
    end
    return d
end

function dop(n, endpt)
    if endpt == neither
        return 2*n - 1
    elseif endpt == left || endpt == right
        return 2*n - 2
    else
        return 2*n - 3
    end
end

function test_rule(descr, nmin, nmax, rule, coefs, name, endpt)
    println("\nTesting ", descr, " rules with ", nmin, " to ",
            nmax, " points")
    for case in zip(rule, coefs, name)
        rulefunc, coeffunc, rulename = case
        maxd = 0.0
        for n = nmin:nmax
            x, w = rulefunc(n, alpha, beta, endpt)
            a, b = coeffunc(2n, alpha, beta)
            p = orthonormal_poly(x, a, b)
            d = discrepancy(n, dop(n, endpt), w, p)
            maxd = max(d, maxd)
        end
        @printf("\t%s: %12.4e\n", rpad(rulename, 20, ' '), maxd)
    end
end

test_rule("plain Gauss", 1, 5, rule, coefs, name, neither)

filter!(x->(x!=xhermite), rule)
filter!(x->(x!=xhermite_coefs), coefs)
filter!(x->(x!="Hermite"), name)

test_rule("left Radau", 1, 5, rule, coefs, name, left)

filter!(x->(x!=xlaguerre), rule)
filter!(x->(x!=xlaguerre_coefs), coefs)
filter!(x->(x!="Laguerre"), name)

test_rule("right Radau", 1, 5, rule, coefs, name, right)
test_rule("Lobatto", 2, 5, rule, coefs, name, both)


text:
# algorithm

abstract type Sinkhorn end

# solver

struct SinkhornSolver{A<:Sinkhorn,M,N,CT,E<:Real,T<:Real,R<:Real,C1,C2}
    source::M
    target::N
    C::CT
    eps::E
    alg::A
    atol::T
    rtol::R
    maxiter::Int
    check_convergence::Int
    cache::C1
    convergence_cache::C2
end

function build_solver(
    μ::AbstractVecOrMat,
    ν::AbstractVecOrMat,
    C::AbstractMatrix,
    ε::Real,
    alg::Sinkhorn;
    atol=nothing,
    rtol=nothing,
    check_convergence=10,
    maxiter::Int=1_000,
)
    # check that source and target marginals have the correct size and are balanced
    checksize(μ, ν, C)
    size2 = checksize2(μ, ν)
    checkbalanced(μ, ν)

    # compute type
    T = float(Base.promote_eltype(μ, ν, one(eltype(C)) / ε))

    # build caches
    cache = build_cache(T, alg, size2, μ, ν, C, ε)
    convergence_cache = build_convergence_cache(T, size2, μ)

    # set tolerances
    _atol = atol === nothing ? 0 : atol
    _rtol = rtol === nothing ? (_atol > zero(_atol) ? zero(T) : sqrt(eps(T))) : rtol

    # create solver
    solver = SinkhornSolver(
        μ, ν, C, ε, alg, _atol, _rtol, maxiter, check_convergence, cache, convergence_cache
    )

    return solver
end

# convergence caches

struct SinkhornConvergenceCache{U,S<:Real}
    tmp::U
    norm_source::S
end

struct SinkhornBatchConvergenceCache{U,S,C}
    tmp::U
    tmp2::U
    norm_source::S
    norm_uKv::U
    norm_diff::U
    isconverged::C
end

function build_convergence_cache(::Type{T}, ::Tuple{}, μ::AbstractVector) where {T}
    tmp = similar(μ, T)
    norm_μ = sum(abs, μ)
    return SinkhornConvergenceCache(tmp, norm_μ)
end

function build_convergence_cache(
    ::Type{T}, size2::Tuple{Int}, μ::AbstractVecOrMat
) where {T}
    tmp = similar(μ, T, size(μ, 1), size2...)
    tmp2 = similar(tmp)
    norm_μ = μ isa AbstractVector ? sum(abs, μ) : sum(abs, μ; dims=1)
    norm_uKv = similar(tmp, 1, size2...)
    norm_diff = similar(tmp, 1, size2...)
    isconverged = similar(tmp, Bool, 1, size2...)
    return SinkhornBatchConvergenceCache(
        tmp, tmp2, norm_μ, norm_uKv, norm_diff, isconverged
    )
end

# Sinkhorn algorithm steps (see solve!)
function init_step!(solver::SinkhornSolver)
    return A_batched_mul_B!(solver.cache.Kv, solver.cache.K, solver.cache.v)
end

function step!(solver::SinkhornSolver, iter::Int)
    μ = solver.source
    ν = solver.target
    cache = solver.cache
    u = cache.u
    v = cache.v
    Kv = cache.Kv
    K = cache.K

    u .= μ ./ Kv
    At_batched_mul_B!(v, K, u)
    v .= ν ./ v
    return A_batched_mul_B!(Kv, K, v)
end

function check_convergence(solver::SinkhornSolver)
    return OptimalTransport.check_convergence(
        solver.source,
        solver.cache.u,
        solver.cache.Kv,
        solver.convergence_cache,
        solver.atol,
        solver.rtol,
    )
end

# API

"""
    sinkhorn(
        μ, ν, C, ε, alg=SinkhornGibbs();
        atol=0, rtol=atol > 0 ? 0 : √eps, check_convergence=10, maxiter=1_000,
    )

Compute the optimal transport plan for the entropically regularized optimal transport
problem with source and target marginals `μ` and `ν`, cost matrix `C` of size
`(length(μ), length(ν))`, and entropic regularization parameter `ε`.

The optimal transport plan `γ` is of the same size as `C` and solves
```math
\\inf_{\\gamma \\in \\Pi(\\mu, \\nu)} \\langle \\gamma, C \\rangle
+ \\varepsilon \\Omega(\\gamma),
```
where ``\\Omega(\\gamma) = \\sum_{i,j} \\gamma_{i,j} \\log \\gamma_{i,j}`` is the entropic
regularization term.

Every `check_convergence` steps it is assessed if the algorithm is converged by checking if
the iterate of the transport plan `G` satisfies
```julia
isapprox(sum(G; dims=2), μ; atol=atol, rtol=rtol, norm=x -> norm(x, 1))
```
The default `rtol` depends on the types of `μ`, `ν`, and `C`. After `maxiter` iterations,
the computation is stopped.

Batch computations for multiple histograms with a common cost matrix `C` can be performed by
passing `μ` or `ν` as matrices whose columns correspond to histograms. It is required that
the number of source and target marginals is equal or that a single source or single target
marginal is provided (either as matrix or as vector). The optimal transport plans are
returned as three-dimensional array where `γ[:, :, i]` is the optimal transport plan for the
`i`th pair of source and target marginals.

See also: [`sinkhorn2`](@ref)
"""
function sinkhorn(μ, ν, C, ε, alg::Sinkhorn; kwargs...)
    # build solver
    solver = build_solver(μ, ν, C, ε, alg; kwargs...)

    # perform Sinkhorn algorithm
    solve!(solver)

    # compute optimal transport plan
    γ = plan(solver)

    return γ
end

"""
    sinkhorn2(
        μ, ν, C, ε, alg=SinkhornGibbs(); regularization=false, plan=nothing, kwargs...
    )

Solve the entropically regularized optimal transport problem with source and target
marginals `μ` and `ν`, cost matrix `C` of size `(length(μ), length(ν))`, and entropic
regularization parameter `ε`, and return the optimal cost.

A pre-computed optimal transport `plan` may be provided. The other keyword arguments
supported here are the same as those in the [`sinkhorn`](@ref) function.

!!! note
    As the `sinkhorn2` function in the Python Optimal Transport package, this function
    returns the optimal transport cost without the regularization term. The cost
    with the regularization term can be computed by setting `regularization=true`.

See also: [`sinkhorn`](@ref)
"""
function sinkhorn2(μ, ν, C, ε, alg::Sinkhorn; regularization=false, plan=nothing, kwargs...)
    γ = if plan === nothing
        sinkhorn(μ, ν, C, ε, alg; kwargs...)
    else
        # check dimensions
        checksize(μ, ν, C)
        size(plan) == size(C) || error(
            "optimal transport plan `plan` and cost matrix `C` must be of the same size",
        )
        plan
    end
    cost = if regularization
        dot_matwise(γ, C) .+
        ε * reshape(sum(LogExpFunctions.xlogx, γ; dims=(1, 2)), size(γ)[3:end])
    else
        dot_matwise(γ, C)
    end

    return cost
end


text:
##########################
# FunctionCell Interface #
##########################
"""
	networklearner(Adj, fl_train, fl_exec, fr_train, ft_exec [;kwargs])

Constructs an untrained cell that when piped data inside, returns a network learner trained
function cell based on the input data and labels.

# Arguments
  * `x::AbstractMatrix` training data (used by `fl_train`, `fl_exec`; if `use_local_data==true`, it is also used by `fr_train`)
  * `Adj::Vector{AbstractAdjacency}` a vector containing the observation relational structures (adjacency objects)
  * `fl_train` local model training 'function'; can be anything that supports the call `fl_train((X,y))`
  * `fl_exec` local model prediction 'function'; can be anything that supports the call `fl_exec(Ml,X)` where `Ml = fl_train((X,y))`
  * `fr_train` relational model training `function`; can be anything that suports the call `fr_train((Xr,y))`
  * `fr_exec` relational model prediction `function`; can be anything that suports the call `fr_exec(Mr,Xr)` where `Mr = fr_train((Xr,y))`
  and `Xr` is a dataset of relational variables generated by the relational learner using the results of the local model prediction
  function and adjacency structures.

# Keyword arguments
  * `priors::Vector{Float64}` class priors (if applicable)
  * `learner::Symbol` relational learner (i.e. variable generator); available options `:rn`, `:wrn`, `:bayesrn` and `:cdrn` (default `:wrn`)
  * `inference::Symbol` collective inference method; available options `:rl`, `:ic` and `:gs` (default `:rl`)
  * `normalize::Bool` whether to normalize the relational variables per-observation to the L1 norm (default `true`)
  * `use_local_data::Bool` whether the relational model should use the local data provided (i.e. in `X`) (default `true`)
  * `f_targets::Function` function that extracts targets from estimates generated by the local/relational models
  (default `f_targets = x->MLDataPattern.targets(indmax,x)`)
  * `obsdim::Int` observation dimension (default `2`)
  * `tol::Float64` maximum admissible mean estimate error for collective inference convergence (default `1e-6`)
  * `κ::Float64` relaxation labeling starting constant, used if `learner == :rl` (default `1.0`)
  * `α::Float64` relaxation labeling decay constant, used if `learner == :rl` (default `0.99`)
  * `maxiter::Int` maximum number of iterations for collective inference (default `100`)
  * `bratio::Float64` percentage of iterations i.e. `maxiter` used for Gibbs sampling burn-in (default `0.1`)

Read the `NetworkLearning.jl` documentation for more information.
"""
networklearner(Adj, fl_train, fl_exec, fr_train, fr_exec; kwargs...) =
	FunctionCell(networklearner, (Adj, fl_train, fl_exec, fr_train, fr_exec), ModelProperties(), "Network Learner"; kwargs...)

############################
# DataCell/Array Interface #
############################
"""
	networklearner(x, Adj, fl_train, fl_exec, fr_train, ft_exec [;kwargs])

Trains a network learner model that using the data `x`.
"""
# Training
networklearner(x::T where T<:CellDataL, Adj, fl_train, fl_exec, fr_train, fr_exec; kwargs...) =
	networklearner((getx!(x), gety(x)), Adj, fl_train, fl_exec, fr_train, fr_exec; kwargs...)
networklearner(x::Tuple{T,S} where T<:AbstractVector where S<:AbstractVector, Adj, fl_train, fl_exec, fr_train, fr_exec; kwargs...) =
	networklearner((mat(x[1], LearnBase.ObsDim.Constant{2}()), x[2]), Adj, fl_train, fl_exec, fr_train, fr_exec; kwargs...)
networklearner(x::Tuple{T,S} where T<:AbstractMatrix where S<:AbstractVector, Adj, fl_train, fl_exec, fr_train, fr_exec; kwargs...) = begin

	@assert nobs(x[1]) == nobs(x[2]) "[networklearner] Expected $(nobs(x[1])) labels/values, got $(nobs(x[2]))."

	# Transform labels first
	enc = labelencn(x[2])
	yenc = label2ind.(x[2],enc)

	# Train model
	nldata = NetworkLearning.fit(NetworkLearning.NetworkLearnerObs, getobs(x[1]), yenc,
			      		Adj, fl_train, fl_exec, fr_train, fr_exec; kwargs...)

	# Build model properties
	modelprops = ModelProperties(nvars(x[1]), length(enc.label), enc)

	FunctionCell(networklearner, Model(nldata, modelprops), "Network Learner")

end

# Execution
networklearner(x::T where T<:CellData, model::Model{<:NetworkLearning.AbstractNetworkLearner}) =
	datacell(networklearner(getx!(x), model), gety(x))

networklearner(x::T where T<:AbstractVector, model::Model{<:NetworkLearning.AbstractNetworkLearner}) =
	networklearner(mat(x, LearnBase.ObsDim.Constant{2}()), model)

networklearner(x::T where T<:AbstractMatrix, model::Model{<:NetworkLearning.AbstractNetworkLearner}) =
	NetworkLearning.predict(model.data, getobs(x))


text:
using Test, Sockets, Erjulix, ErlangTerm

localhost = Sockets.localhost
erl = UDPSocket()     # erlang test socket
ep = Erjulix.getPort(1000)
@test bind(erl, localhost, ep)

pport = Erjulix.getPort(1000)
println("pServer setup")
ps = pServer(pport)

@test ps.state == :runnable
send(erl, localhost, pport, serialize(:srv))
hp, pkg = recvfrom(erl)
msg = deserialize(pkg)
@test hp.host == localhost
@test first(msg) == :ok
@test last(msg) == repr(Erjulix._ESM[end])
println("EvalServer at $(hp.host):$(hp.port)")

# test eval
println("test :eval")
send(erl, hp.host, hp.port, serialize((:eval, "sum(1:10)")))
println("msg sent!")
sleep(0.5)
@test Erjulix._ESM[end]._eServer.state == :runnable
println("server ok")
pkg = recv(erl)
println("got message from server")
msg = deserialize(pkg)
@test msg == (:ok, 55)

# test call
println("test :call")
send(erl, hp.host, hp.port, serialize((:call, :sum, [collect(11:20)])))
sleep(0.2)
@test Erjulix._ESM[end]._eServer.state == :runnable
pkg = recv(erl)
msg = deserialize(pkg)
@test msg == (:ok, 155)

# test set
println("test :set")
send(erl, hp.host, hp.port, serialize((:set, :a, collect(21:30))))
sleep(0.2)
@test Erjulix._ESM[end]._eServer.state == :runnable
pkg = recv(erl)
msg = deserialize(pkg)
@test msg == (:ok, :nil)
@test Erjulix._ESM[end].a == collect(21:30)

# test exit
println("test :exit")
send(erl, hp.host, hp.port, serialize(:exit))
sleep(0.5)
@test Erjulix._ESM[end]._eServer.state == :done

send(erl, localhost, pport, serialize(:exit))
sleep(0.5)
@test ps.state == :done

max_stars_repo_name: pbayer/erjulix

max_issues_repo_name: pbayer/erjulix

max_forks_repo_name: pbayer/erjulix