text:
module ComradeGalactic

using Comrade
using Reexport
using Distributions
using ForwardDiff
using LinearAlgebra

@reexport using GalacticOptim

export laplace

function GalacticOptim.OptimizationFunction(post::Comrade.Posterior, args...; kwargs...)
    throw("Transform the posterior first using `asflat` or `ascube`")
end

function GalacticOptim.OptimizationFunction(post::Comrade.TransformedPosterior, args...; kwargs...)
    ℓ(x,p) = -logdensityof(post, x)
    return OptimizationFunction(ℓ, args...; kwargs...)
end

"""
    `laplace(prob, opt, args...; kwargs...)`
Compute the Laplace or Quadratic approximation to the prob or posterior.
The `args` and `kwargs` are passed the the GalacticOptim.solve function.

Note the quadratic approximation is in the space of the transformed posterior
not the usual parameter space. This is better for constrained problems where
we may run up against a boundary.
"""
function laplace(prob::OptimizationProblem, opt, args...; kwargs...)
    sol = solve(prob, opt, args...; kwargs...)
    f = Base.Fix2(prob.f, nothing)
    J = ForwardDiff.hessian(f, sol)
    @. J += 1e-5 # add noise to help with positive definitness
    h = J*sol
    return MvNormalCanon(h, Symmetric(J))
end

# """
#     solve(prob, opt, args...; transform=nothing, kwargs...)
# Find the optimum of the problem.

# Optional keyword argument `transform``
# will transform the solution to parameter space, and the return will be
# a tuple with the galactic optim solution in the first element and the optimum
# location in model space in the second argument
# """
# function GalacticOptim.solve(prob::OptimizationProblem, opt, transform::Union{Nothing, HypercubeTransform.AbstractHypercubeTransform, HypercubeTransform.TransformVariables.AbstractTransform, HypercubeTransform.NamedFlatTransform}, args...; kwargs...)
#     sol = solve(prob, opt, args...; kwargs...)
#     if isnothing(transform)
#         return sol
#     else
#         p = HypercubeTransform.transform(transform, sol.u)
#         return sol, p
#     end
# end

end

ext: jl


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: dynamodb_streams
using AWS.Compat
using AWS.UUIDs

"""
    DescribeStream()

Returns information about a stream, including the current status of the stream, its Amazon Resource Name (ARN), the composition of its shards, and its corresponding DynamoDB table.  You can call DescribeStream at a maximum rate of 10 times per second.  Each shard in the stream has a SequenceNumberRange associated with it. If the SequenceNumberRange has a StartingSequenceNumber but no EndingSequenceNumber, then the shard is still open (able to receive more stream records). If both StartingSequenceNumber and EndingSequenceNumber are present, then that shard is closed and can no longer receive more data.

# Required Parameters
- `StreamArn`: The Amazon Resource Name (ARN) for the stream.

# Optional Parameters
- `ExclusiveStartShardId`: The shard ID of the first item that this operation will evaluate. Use the value that was returned for LastEvaluatedShardId in the previous operation.
- `Limit`: The maximum number of shard objects to return. The upper limit is 100.
"""
describe_stream(StreamArn; aws_config::AWSConfig=global_aws_config()) = dynamodb_streams("DescribeStream", Dict{String, Any}("StreamArn"=>StreamArn); aws_config=aws_config)
describe_stream(StreamArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dynamodb_streams("DescribeStream", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("StreamArn"=>StreamArn), args)); aws_config=aws_config)

"""
    GetRecords()

Retrieves the stream records from a given shard. Specify a shard iterator using the ShardIterator parameter. The shard iterator specifies the position in the shard from which you want to start reading stream records sequentially. If there are no stream records available in the portion of the shard that the iterator points to, GetRecords returns an empty list. Note that it might take multiple calls to get to a portion of the shard that contains stream records.   GetRecords can retrieve a maximum of 1 MB of data or 1000 stream records, whichever comes first.

# Required Parameters
- `ShardIterator`: A shard iterator that was retrieved from a previous GetShardIterator operation. This iterator can be used to access the stream records in this shard.

# Optional Parameters
- `Limit`: The maximum number of records to return from the shard. The upper limit is 1000.
"""
get_records(ShardIterator; aws_config::AWSConfig=global_aws_config()) = dynamodb_streams("GetRecords", Dict{String, Any}("ShardIterator"=>ShardIterator); aws_config=aws_config)
get_records(ShardIterator, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dynamodb_streams("GetRecords", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ShardIterator"=>ShardIterator), args)); aws_config=aws_config)

"""
    GetShardIterator()

Returns a shard iterator. A shard iterator provides information about how to retrieve the stream records from within a shard. Use the shard iterator in a subsequent GetRecords request to read the stream records from the shard.  A shard iterator expires 15 minutes after it is returned to the requester.

# Required Parameters
- `ShardId`: The identifier of the shard. The iterator will be returned for this shard ID.
- `ShardIteratorType`: Determines how the shard iterator is used to start reading stream records from the shard:    AT_SEQUENCE_NUMBER - Start reading exactly from the position denoted by a specific sequence number.    AFTER_SEQUENCE_NUMBER - Start reading right after the position denoted by a specific sequence number.    TRIM_HORIZON - Start reading at the last (untrimmed) stream record, which is the oldest record in the shard. In DynamoDB Streams, there is a 24 hour limit on data retention. Stream records whose age exceeds this limit are subject to removal (trimming) from the stream.    LATEST - Start reading just after the most recent stream record in the shard, so that you always read the most recent data in the shard.
- `StreamArn`: The Amazon Resource Name (ARN) for the stream.

# Optional Parameters
- `SequenceNumber`: The sequence number of a stream record in the shard from which to start reading.
"""
get_shard_iterator(ShardId, ShardIteratorType, StreamArn; aws_config::AWSConfig=global_aws_config()) = dynamodb_streams("GetShardIterator", Dict{String, Any}("ShardId"=>ShardId, "ShardIteratorType"=>ShardIteratorType, "StreamArn"=>StreamArn); aws_config=aws_config)
get_shard_iterator(ShardId, ShardIteratorType, StreamArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dynamodb_streams("GetShardIterator", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ShardId"=>ShardId, "ShardIteratorType"=>ShardIteratorType, "StreamArn"=>StreamArn), args)); aws_config=aws_config)

"""
    ListStreams()

Returns an array of stream ARNs associated with the current account and endpoint. If the TableName parameter is present, then ListStreams will return only the streams ARNs for that table.  You can call ListStreams at a maximum rate of 5 times per second.

# Optional Parameters
- `ExclusiveStartStreamArn`: The ARN (Amazon Resource Name) of the first item that this operation will evaluate. Use the value that was returned for LastEvaluatedStreamArn in the previous operation.
- `Limit`: The maximum number of streams to return. The upper limit is 100.
- `TableName`: If this parameter is provided, then only the streams associated with this table name are returned.
"""
list_streams(; aws_config::AWSConfig=global_aws_config()) = dynamodb_streams("ListStreams"; aws_config=aws_config)
list_streams(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = dynamodb_streams("ListStreams", args; aws_config=aws_config)


text:
using SubmodularMaximization
using POMDPs
using MCTS
using PyPlot
using Base.Iterators
using Base.Threads
using JLD2
using Statistics

close()

data_file = "./data/horizon_entropy_data.jld2"
reprocess = false

steps = 100
num_robots = 4

trials = 1:20

horizons = 1:6

# We will drop a fraction of each trial so that the filters have time to
# converge to steady states
trial_steps = 20:steps

# Note: we will run trials in threads so the solvers do not have to be threaded
solvers = [solve_sequential, solve_myopic]
solver_strings = map(string, solvers)

function get_data()
  data_mutex = Mutex()
  data = Dict()

  if !isfile(data_file) || reprocess
    tests = product(solvers, horizons, trials)
    test_partition = partition(tests, nthreads())

    for (ii, block) in enumerate(test_partition)
      println("Block num: ", ii, " of ", length(test_partition))

      @time @threads for (solver, horizon, trial) in block

        configs = MultiRobotTargetTrackingConfigs(num_robots,
                                                  horizon=horizon,
                                                 )

        trial_data = target_tracking_experiment(steps=steps,
                                                num_robots=num_robots,
                                                configs=configs,
                                                solver=solver
                                               )

        lock(data_mutex)
        data[string(solver), horizon, trial] =
          (data=trial_data, configs=configs)
        unlock(data_mutex)
      end
    end

    @save data_file data
  else
    @load data_file data
  end

  data
end

@time data = get_data()

entropies = Dict(key=>map(x->entropy(x.target_filters), value.data)
                 for (key, value) in data)

concatenated_entropy = map(product(solver_strings, horizons)) do (solver,
                                                                  horizon)
  vcat(map(trials) do trial
    entropies[solver, horizon, trial][trial_steps]
  end...)
end[:]
titles = map(product(solver_strings, horizons)) do (solver, horizon)
  string(solver, " ", horizon)
end[:]

boxplot(concatenated_entropy, notch=false, vert=false)
title("Entropy distribution")
xlabel("Entropy (bits)")
yticks(1:length(titles), titles)

save_fig("fig", "entropy_by_solver_horizon")

means = map(mean, concatenated_entropy)
for (title, mean) in zip(titles, means)
  println(title, "(mean): ", mean)
end

println()

println("Sequential entropy is so much lower (better) than myopic:")
gaps =  means[2:2:end] - means[1:2:end]
for (horizon, gap) in zip(horizons, gaps)
  println("Horizon: ", horizon, " Gap: ", gap)
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
# declare 3 global variables
x = y = 0
#in test function a and x are local
function test()
    local x = 1
    println("x = $x")
    global y = 1
end
test()
println("x = $x")
println("y = $y")

max_stars_repo_name: elucian/julia

max_issues_repo_name: elucian/julia

max_forks_repo_name: elucian/julia


text:
#NOTE: May also want to check multiple years
const YEAR_PAIRS = [(2009,2010),(2013,2014)]

function singleYearPlots(data::NCCSData; yearPairs::Vector{Tuple{Int,Int}} = YEAR_PAIRS,
  outputPath::String = GRAPH_PATH)
  local plotNames::Vector{String} = Vector{String}()
  local plots::Vector{PlotContainer} = Vector{PlotContainer}()
  local subdf::SubDataFrame

  for p ∈ YEAR_PAIRS
    @assert p[2] == p[1] + 1 #only support one year difference for now
    subdf = view(data.df, data.df[:fisyr] .== p[2], :)
    subdf = view(subdf, (!ismissing).(subdf[:lpreturnbyyr]), :)
    subdf = view(subdf, (!ismissing).(subdf[:preturnbyyr]), :)

    push!(plotNames, "Performance Percentile in $(p[2]) v. $(p[1])")
    push!(plots, plot(subdf, x=:lpreturnbyyr, y=:preturnbyyr,
      #size=:lnetassets,
      Guide.xlabel("Return in $(p[1])"),Guide.ylabel("Return in $(p[2])"),
      Geom.point,
      Coord.cartesian(xmin = 0.0, xmax=1.0, ymin=0.0, ymax = 1.0)))
  end

  for i ∈ 1:length(plots) #write the graphs
    draw(PNG("$outputPath\\$(plotNames[i]).png", 9inch, 7inch), plots[i])
    #savefig(plots[i], "$outputPath\\$(plotNames[i])$(graphsuffix).$ext")
    #println("Graph $(plotNames[i]) written.")
  end

  return nothing
end

function persistancePlots(data::NCCSData)

  singleYearPlots(data)

  println("completed persistance plots")
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
#=
AdvancedPrim:
- Julia version:
- Author: samuel
- Date: 2021-09-10
=#

using SuiteSparseGraphBLAS
using SparseArrays
using LinearAlgebra

function stoppingCondition(m, n)

infinito = true

for i = 1:n
        if(m[i] == 0.0)
            return false
    end
end
    return infinito
end

#A is the input matrix, n is the number of nodes, m is the source matrix. Output: minimum spanning tree and minimum spanning tree cost
function a_prim(A, n, m)

    d = GBVector{Float64}(n)

    weight = 0.0

    d = A[1,:]
    mst = GBVector{Float64}(n) #minimum spanning tree

    print("INITIAL WEIGHT: ")
    print(weight)
    print("\n\n")
    while(stoppingCondition(m, n) == false)

        u = argmin(m'+d)
        m[u[2]] = Inf
        push!(mst, d[u[2]])
        weight = weight + d[u[2]]
        print("WEIGHT: ")
        print(weight)
        d = emul(d, A[u[1],:],  BinaryOps.MIN)
        print("  ITERATION FINISHED")
        print("\n\n\n")
    end

    print("COMPLETED")
    print(" ALL THE NODES HAVE BEEN PROCESSED (INF)")
    print("\n\n\n")
    print(m')
    print("\n\n\n")
    print("MINIMUM SPANNING TREE EDGES")
    print("\n\n\n")
    print(mst)
    return weight, mst

end


text:
# Do two-sided equal variances T tests on data from a CSV file.
using CSV
using DataFrames
using HypothesisTests

# Returns a dafaframe by reading a csv file
# make sure that comment lines are preceded by #
function import_csv( file::String )
  CSV.read(file,comment="#")
end

function two_sided_t_test_pvalue( value_field::CSV.Column{Int64,Int64}, discriminator_field::CSV.Column{Bool,Bool} )
  pvalue( EqualVarianceTTest( value_field[ discriminator_field ], value_field[ .~discriminator_field ] ))
end

function two_sided_t_test( value_field::CSV.Column{Int64,Int64}, discriminator_field::CSV.Column{Bool,Bool} )
   EqualVarianceTTest( value_field[ discriminator_field ], value_field[ .~discriminator_field ] )
end

function two_sided_t_test_pvalue( value_field::CSV.Column{Float64,Float64}, discriminator_field::CSV.Column{Bool,Bool} )
  pvalue( EqualVarianceTTest( value_field[ discriminator_field ], value_field[ .~discriminator_field ] ))
end

function two_sided_t_test( value_field::CSV.Column{Float64,Float64}, discriminator_field::CSV.Column{Bool,Bool} )
   EqualVarianceTTest( value_field[ discriminator_field ], value_field[ .~discriminator_field ] )
end

# Do specific tests for nactive, complexity, degeneracy, and steps using fault_tol as the discriminator field
function nactive_test( file::String )
  df = import_csv( file )
  pvalue( EqualVarianceTTest( df.nactive[ df.fault_tol ], df.nactive[ .~df.fault_tol ] ) )
end

function complexity_test( file::String )
  df = import_csv( file )
  pvalue( EqualVarianceTTest( df.complexity[ df.fault_tol ], df.complexity[ .~df.fault_tol ] ) )
end

function degeneracy_test( file::String )
  df = import_csv( file )
  pvalue( EqualVarianceTTest( df.degeneracy[ df.fault_tol ], df.degeneracy[ .~df.fault_tol ] ) )
end

function steps_test( file::String )
  df = import_csv( file )
  pvalue( EqualVarianceTTest( df.steps[ df.fault_tol ], df.steps[ .~df.fault_tol ] ) )
end


text:
using Reexport

@reexport module JointDistanceDistribution
    include("joint_distance_distribution.jl")
end


text:
"""
Creation and manipulation of ∇×∇× differential operator
"""
module Curlcurls

export Curlcurl

files = (
    "1D/Curlcurls1D.jl",
    "2D/Symmetric/Curlcurls2D.jl",
    # "3D/Curlcurls3D.jl"
    )

using ..Domains
using ..Lattices
using ..Laplacians
using ..Points
using SparseArrays

import LinearAlgebra: I
import ..Symmetric, ..Unsymmetric

struct Curlcurl{N,CLASS}
    cc0::SparseMatrixCSC{ComplexF64,Int} # constant
    cc1::SparseMatrixCSC{ComplexF64,Int} # linear in kx,ky,kz
    cc2::SparseMatrixCSC{ComplexF64,Int} # quadratic in kx,ky,kz
    complex_scaling::SparseMatrixCSC{ComplexF64,Int}
end

foreach(include,files)

function Base.conj!(cc::Curlcurl{N}) where N
    conj!(cc.cc0)
    conj!(cc.cc1)
    conj!(cc.cc2)
    return cc
end
Base.conj(cc::Curlcurl{N}) where N = Curlcurl{N}(conj(cc.cc0),conj(cc.cc1),conj(cc.cc2),conj(cc.complex_scaling))

################################################################################
# Pretty Printing

import ..PRINTED_COLOR_DARK
import ..PRINTED_COLOR_VARIABLE

function Base.show(io::IO,cc::Curlcurl{N,CLASS}) where {N,CLASS}
    print(io,"$(N)D ")
    CLASS<:Symmetric ? print(io,"Symmetric ") : nothing
    CLASS<:Unsymmetric ? print(io,"Unsymmetric ") : nothing
    printstyled(io,"Curlcurl",color=PRINTED_COLOR_DARK)
    print(io,"(call with ")
    printstyled(io,"ky",color=PRINTED_COLOR_VARIABLE)
    print(io,", ")
    printstyled(io,"kz",color=PRINTED_COLOR_VARIABLE)
    print(io,")")
end

end


text:
import Base: print

function Base.show(io::IO, p::Polynomial)
    # if p.homogenized
    #     vars = ["x$i" for i=0:nvariables(p)-1]
    # else
    #     vars = ["x$i" for i=1:nvariables(p)]
    # end
    print_poly(io, p, variables(p))
end

# function Base.show(io::IO, P::PolynomialSystem)
#     for p in P.polys
#         print_poly(io, p, P.vars)
#         print(io, "\n")
#     end
# end

#helpers

function print_poly(io::IO, p::Polynomial{T}, vars) where T
    first = true
    exps = exponents(p)
    cfs = coefficients(p)

    m, n = size(exps)

    for i=1:n
        exp = exps[:, i]
        coeff = cfs[i]

        if (!first && show_plus(coeff))
            print(io, "+")
        end
        first = false

        if (coeff != 1 && coeff != -1) || exp == zeros(Int, m)
            show_coeff(io, coeff)
        elseif coeff == -1
            print(io, "-")
        end

        for (var, power) in zip(vars, exp)
            if power == 1
                print(io, "$(pretty_var(var))")
            elseif power > 1
                print(io, "$(pretty_var(var))$(pretty_power(power))")
            end
        end
    end

    if first
        print(io, zero(T))
    end
end

function unicode_subscript(i)
    if i == 0
        "\u2080"
    elseif i == 1
        "\u2081"
    elseif i == 2
        "\u2082"
    elseif i == 3
        "\u2083"
    elseif i == 4
        "\u2084"
    elseif i == 5
        "\u2085"
    elseif i == 6
        "\u2086"
    elseif i == 7
        "\u2087"
    elseif i == 8
        "\u2088"
    elseif i == 9
        "\u2089"
    end
end

function unicode_superscript(i)
    if i == 0
        "\u2070"
    elseif i == 1
        "\u00B9"
    elseif i == 2
        "\u00B2"
    elseif i == 3
        "\u00B3"
    elseif i == 4
        "\u2074"
    elseif i == 5
        "\u2075"
    elseif i == 6
        "\u2076"
    elseif i == 7
        "\u2077"
    elseif i == 8
        "\u2078"
    elseif i == 9
        "\u2079"
    end
end

pretty_power(pow::Int) = join(map(unicode_superscript, reverse(digits(pow))))

function pretty_var(var::String)
    m = match(r"([a-zA-Z]+)(?:_*)(\d+)", var)
    if m === nothing
        var
    else
        base = string(m.captures[1])
        index = parse(Int, m.captures[2])
        base * join(map(unicode_subscript, reverse(digits(index))))
    end
end
pretty_var(var) = pretty_var(string(var))

# helpers
show_plus(x::Real) = x >= 0
show_plus(x::Complex) = x != -1

show_coeff(io::IO, x::Real) = print(io, x)
function show_coeff(io::IO, x::Complex)
    if imag(x) == 0.0
        print(io, convert(Float64, x))
    else
        print(io, "($(x))")
    end
end


text:
# Autogenerated wrapper script for FGlT_jll for aarch64-apple-darwin
export libfglt

using CompilerSupportLibraries_jll
using cilkrts_jll
JLLWrappers.@generate_wrapper_header("FGlT")
JLLWrappers.@declare_library_product(libfglt, "@rpath/libfglt.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, cilkrts_jll)
    JLLWrappers.@init_library_product(
        libfglt,
        "lib/libfglt.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export Event, Events, handle_event_point

struct Event{T}
    point::Point{T}
    segments::Set{Segment{T}}
end

Event(point::Point{T}) where {T<:Number} = Event(point, Set{Segment{T}}())

Event(x::Float64, y::Float64) = Event(Point(x, y))
Event(x::Number, y::Number) = Event(convert(Float64, x), convert(Float64, y))

Base.:(==)(e1::Event, e2::Event) = e1.point == e2.point
Base.isless(e1::Event, e2::Event) = e1.point < e2.point


text:
using CPLEX, JuMP, MathProgBase

mod = Model(solver=CplexSolver(CPX_PARAM_PRELINEAR=0, CPX_PARAM_PREIND=0, CPX_PARAM_ADVIND=0, CPX_PARAM_MIPSEARCH=1,CPX_PARAM_MIPCBREDLP=0,CPX_PARAM_CUTSFACTOR=1))
m_internal = MathProgBase.model(CplexSolver())

MathProgBase.loadproblem!(m_internal, "location.lp")

# grab MathProgBase data
c = MathProgBase.getobj(m_internal)
A = MathProgBase.getconstrmatrix(m_internal)
m, n = size(A)
xlb = MathProgBase.getvarLB(m_internal)
xub = MathProgBase.getvarUB(m_internal)
l = MathProgBase.getconstrLB(m_internal)
u = MathProgBase.getconstrUB(m_internal)
vtypes = MathProgBase.getvartype(m_internal)

# populate JuMP model with data from internal model
@defVar(mod, x[1:n])
for i in 1:n
    setLower(x[i], xlb[i])
    setUpper(x[i], xub[i])
    (vtypes[i] == :Bin || vtypes[i] == :Int) ? mod.colCat[x[i].col] = :Int : nothing # change vartype to integer when appropriate
end
At = A' # transpose to get useful row-wise sparse representation
for i in 1:At.n
    @addConstraint( mod, l[i] <= sum{ At.nzval[idx]*x[At.rowval[idx]], idx = At.colptr[i]:(At.colptr[i+1]-1) } <= u[i] )
end
@setObjective(mod, Min, sum{ c[i]*x[i], i=1:n })

function callback(cb)
    xval = MathProgBase.cbgetlpsolution(cb)
    objval = cbgetnodeobjval(cb)
    best = MathProgBase.cbgetobj(cb)
    feas = cbgetintfeas(cb)

    if  MathProgBase.cbgetexplorednodes(cb) < 2
        maxinf = maxobj = bestj = -Inf
        for j in 1:MathProgBase.numvar(mod)
            if feas[j] == 1
                xj_inf = xval[j] - floor(xval[j])
                if xj_inf > 0.5
                    xj_inf = 1.0 - xj_inf
                end
                if (xj_inf >= maxinf) && (xj_inf > maxinf || abs(c[j]) >= maxobj)
                    bestj = j
                    maxinf = xj_inf
                    maxobj = abs(c[j])
                end
            end
        end
        if bestj >= 1
            xj_lo = floor(xval[bestj])
            #println("branching on $bestj")
            addBranch(cb, x[bestj] >= xj_lo + 1, objval)
            addBranch(cb, x[bestj] <= xj_lo,     objval)
        end
    end
end

solve(mod; load_model_only=true)
setBranchCallback(mod, callback)
solve(mod)


text:
using CircularList
using Test

@testset "CircularList.jl" begin

    # general/peek
    CL = circularlist(0)
    @test length(CL) == 1
    @test current(CL) == previous(CL)
    @test current(CL) == next(CL)

    # insert!
    start = current(CL)
    insert!(CL, 1)
    insert!(CL, 2)
    @test length(CL) == 3
    @test start.next.data == 1
    @test start.next.next.data == 2
    @test start.next.next.next.data == 0

    # delete!
    @test current(CL).data == 2
    delete!(CL)
    @test length(CL) == 2
    @test current(CL).data == 1
    @test current(CL).next.data == 0
    @test current(CL).next.next.data == 1
    @test current(CL).prev.data == 0
    @test current(CL).prev.prev.data == 1

    # auto resize feature
    CL = circularlist("str_0", capacity = 5)
    for i in 1:100
        insert!(CL, "str_$i")
    end
    @test length(CL) == 101
    @test CL.capacity == 160   # 5 * 2 * 2 * 2 * 2 * 2 = 160

    # forward -- testing next links
    let curr = current(CL)
        for i in 1:length(CL)
            forward!(CL)
        end
        @test current(CL) == curr
    end

    # backward -- testing prev links
    let curr = current(CL)
        for i in 1:length(CL)
            backward!(CL)
        end
        @test current(CL) == curr
    end

    CL = circularlist(collect(1:10))
    @test eltype(CL) == Int
    @test current(CL).data == 1
    @test previous(CL).data == 10
    @test next(CL).data == 2

    # test iteration
    @test [x for x in CL] == collect(1:10)
    @test sum(x for x in CL) == 55

    # test head/tail
    @test head(CL).data == 1
    @test tail(CL).data == 10

    # mixing delete! & insert! to ensure integrity
    delete!(CL); delete!(CL); delete!(CL)
    insert!(CL,11); insert!(CL,12); insert!(CL,13)
    @test length(CL) == 10
    @test size(CL) == (10,)

    # dummy tests to pass coverge
    println("dummy tests")
    show(CL)
    println()
    show(current(CL))
    println()
    println()
end


text:
# This file was generated, do not modify it. # hide
using MLJ
using PrettyPrinting
import DataFrames: DataFrame, select!, Not, describe
import Statistics
using Dates
using PyPlot
ioff() # hide
using UrlDownload

MLJ.color_off() # hide

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


text:
# -------- Opening a window ---------------
# Forward reference for @cfunction
function windowEventWatcher end
const window_event_watcher_cfunc = Ref(Ptr{Nothing}(0))

const window_paused = Threads.Atomic{UInt8}(0) # Whether or not the game should be running (if lost focus)

function makeWinRenderer(title = "GameZero Julia", w=400,h=400 )
    global winWidth, winHeight, winWidth_highDPI, winHeight_highDPI

    win = SDL2.CreateWindow(title,
        Int32(SDL2.WINDOWPOS_CENTERED()), Int32(SDL2.WINDOWPOS_CENTERED()), Int32(w), Int32(h),
        UInt32(SDL2.WINDOW_ALLOW_HIGHDPI|SDL2.WINDOW_OPENGL|SDL2.WINDOW_SHOWN));
        SDL2.SetWindowMinimumSize(win, Int32(w), Int32(h))
    window_event_watcher_cfunc[] = @cfunction(windowEventWatcher, Cint, (Ptr{Nothing}, Ptr{SDL2.Event}))
    SDL2.AddEventWatch(window_event_watcher_cfunc[], win);

    renderer = SDL2.CreateRenderer(win, Int32(-1), UInt32(SDL2.RENDERER_ACCELERATED | SDL2.RENDERER_PRESENTVSYNC))
    SDL2.SetRenderDrawBlendMode(renderer, UInt32(SDL2.BLENDMODE_BLEND))
    return win,renderer
end

# This function handles all window events.
# We currently do no allow window resizes
function windowEventWatcher(data_ptr::Ptr{Cvoid}, event_ptr::Ptr{SDL2.Event})::Cint
    global winWidth, winHeight, cam, window_paused, renderer, win
    ev = unsafe_load(event_ptr, 1)
    ee = ev._Event
    t = UInt32(ee[4]) << 24 | UInt32(ee[3]) << 16 | UInt32(ee[2]) << 8 | ee[1]
    t = SDL2.Event(t)
    if (t == SDL2.WindowEvent)
        event = unsafe_load( Ptr{SDL2.WindowEvent}(pointer_from_objref(ev)) )
        winevent = event.event;  # confusing, but that's what the field is called.
        if (winevent == SDL2.WINDOWEVENT_FOCUS_LOST || winevent == SDL2.WINDOWEVENT_HIDDEN || winevent == SDL2.WINDOWEVENT_MINIMIZED)
            # Stop game playing when out of focus
                window_paused[] = 1
            #end
        elseif (winevent == SDL2.WINDOWEVENT_FOCUS_GAINED || winevent == SDL2.WINDOWEVENT_SHOWN)
            window_paused[] = 0
        end
    end
    return 0
end

function getWindowSize(win)
    w,h,w_highDPI,h_highDPI = Int32[0],Int32[0],Int32[0],Int32[0]
    SDL2.GetWindowSize(win, w, h)
    SDL2.GL_GetDrawableSize(win, w_highDPI, h_highDPI)
    return w[],h[],w_highDPI[],h_highDPI[]
end


text:
using Distributed
using LinearAlgebra

BLAS.set_num_threads(1)

# Apply a sum of distributed matrices and reduce
# after applying to a vector.

# Number of matrices in our sum.
nmats = 2

# Specify the number of processors Julia will use.
# For now we just say it is the same as the number of
# matrices in our sum.
addprocs(nmats)

# Size of out matrices.
@everywhere n = 4

# Set of matrices in our sum, each created on a
# separate worker.
As = [@spawnat(workers()[nmat], randn(n, n)) for nmat in 1:nmats]

# Vector we will apply to our matrices, distributed to
# every worker.
@everywhere v = randn(n)

# Distributed matrix-vector multiplication for each distributed
# matrix in our sum, then reduced over with `+`.
Av = @distributed (+) for n in 1:nmats
  # Need to call `fetch` to instantiate the matrix
  # on the worker where it was defined above.
  fetch(As[n]) * v
end

# Test by fetching and summing locally
A = sum(fetch.(As))
@show norm(Av - A * v)

max_stars_repo_name: mtfishman/ITensorParallel


text:
#Subject
"""
    setkelauto!(data::T, kelauto::Bool) where T <: PKSubject

Set range for elimination parameters calculation for subject.

* `data`     - PK subject;
* `kelauto`  - value.
"""
function setkelauto!(data::T, kelauto::Bool) where T <: PKSubject
    if !kelauto
        if !(data.kelrange.kelend > data.kelrange.kelstart > 0) error("Start point: $(data.kelrange.kelstart) end point: $(data.kelrange.kelend), check that data.kelrange.kelend > data.kelrange.kelstart > 0") end
    end
    data.kelauto = kelauto
    data
end

#DS ind Int
"""
    setkelauto!(data::DataSet{T}, kelauto::Bool, ind::Int) where T <: PKSubject
"""
function setkelauto!(data::DataSet{T}, kelauto::Bool, ind::Int) where T <: PKSubject
    setkelauto!(data[ind], kelauto)
    data
end
#DS iter Int
"""
    setkelauto!(data::DataSet{T}, kelauto::Bool, inds::Union{Vector{Int}, UnitRange{Int}, Tuple{Vararg{Int}}}) where T <: PKSubject
"""
function setkelauto!(data::DataSet{T}, kelauto::Bool, inds::Union{Vector{Int}, UnitRange{Int}, Tuple{Vararg{Int}}}) where T <: PKSubject
    for i in inds
        setkelauto!(data[i], kelauto)
    end
    data
end
#DS all
"""
    setkelauto!(data::DataSet{T}, kelauto::Bool) where T <: PKSubject
"""
function setkelauto!(data::DataSet{T}, kelauto::Bool) where T <: PKSubject
    for i = 1:length(data)
        setkelauto!(data[i], kelauto)
    end
    data
end
#DS Dict
"""
    setkelauto!(data::DataSet{T}, kelauto::Bool, sort::Dict) where T <: PKSubject
"""
function setkelauto!(data::DataSet{T}, kelauto::Bool, sort::Dict) where T <: PKSubject
    for i = 1:length(data)
        if sort ⊆ data[i].id setkelauto!(data[i], kelauto) end
    end
    data
end

#GET subj
"""
    getkelauto!(data::T) where T <: PKSubject
"""
function getkelauto(data::T) where T <: PKSubject
    data.kelauto
end


text:
module O_structs

export u_type                     # structure of uncertainty parameters
export ms_type                    # structure of master problem sets
export mp_type                    # structure of master problem parameters
export ps_type                    # structure of subproblem sets
export pp_type                    # structure of subproblem parameters

mutable struct u_type
       ni::Int64                  # number of subproblems
       nx::Int64                  # size of vector x
      nx0::Int64                  # size of master vector decisions x0
       nh::Int64                  # size of vector h
       nc::Int64                  # size of vector c
        h::Array{Float64,2}       # array of vector h
        c::Array{Float64,2}       # array of vector c
end

mutable struct ms_type
       P::UnitRange{Int64}       # set of technologies
      I0::UnitRange{Int64}       # set of "investment" nodes
       I::UnitRange{Int64}       # set of "operational" nodes
end

mutable struct mp_type
       κ::Float64                 # years of operational problem (yr)
       xh::Array{Float64,2}       # historical installed capacity (GW)
       xm::Array{Float64,1}       # maximum accumulated capacity (GW)
       ci::Array{Float64,2}       # investment cost (£/GW)
       cf::Array{Float64,1}       # fix OM cost (£/GWyr)
       π0::Array{Float64,1}       # probability associated to "investment" node i0 (-)
        π::Array{Float64,1}       # probability associated to "operational" node i (-)
      map::Vector{Array{Int64,1}} # map "operational" node i → "investment" node i0
end

mutable struct ps_type
        P::UnitRange{Int64}       # set of technologies
        G::UnitRange{Int64}       # set of conventional technologies
        B::UnitRange{Int64}       # set of storage technologies
        R::UnitRange{Int64}       # set of renewable technologies
        S::UnitRange{Int64}       # set of seasons
        H::UnitRange{Int64}       # hours within a season
       G0::UnitRange{Int64}       # set of conventional technologies exept nuclear
       gn::Int64                  # nuclear index
       b0::Int64                  # index 0 in x of storage technologies
       r0::Int64                  # index 0 in x of renewables technologies
end

mutable struct pp_type
      cvg::Array{Float64,1}       # var OM cost (conv. tech.) (£/MWh)
      cfg::Array{Float64,1}       # fuel cost (conv. tech.) (£/MWh)
       eg::Array{Float64,1}       # CO2 emissions (conv. tech.) (tCO₂/MWh)
       ηg::Array{Float64,1}       # efficiency (conv. tech.) (-)
       rg::Array{Float64,1}       # ramping limitations (conv. tech.) (MW)
       ηb::Array{Float64,1}       # efficiency (storage tech.) (-)
       pb::Array{Float64,1}       # ch/disch power (storage tech.) (MW/MW)
       pr::Array{Float64,3}       # renewable energy prod (re. tech.) (MW/MW)
       pd::Array{Float64,2}       # energy demand (MW)
       νd::Float64                # scaling demand (-)
       cs::Float64                # shedding demand cost (£/MW)
       co::Float64                # CO₂ cost (conv. tech.) (£/tCO₂)
      lco::Float64                # CO₂ yearly limit (tCO₂)
      ρco::Float64                # scaling CO₂ limit (-)
        α::Float64                # seasonal weight (-)
end

end


text:
using Random, JuMP
using NLPModels
using NLPModelsJuMP
using MadNLP
using NLPModelsIpopt, Ipopt
using BlockNLPModels
using Test
using LinearAlgebra
using SparseArrays

include("simple_block_model.jl")


text:
using Combinatorics

term = "RCode"
i = 0
pcnt = factorial(length(term))
print("All the permutations of ", term, " (", pcnt, "):\n    ")
for p in permutations(split(term, ""))
    global i
    print(join(p), " ")
    i += 1
    i %= 12
    i != 0 || print("\n    ")
end
println()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    evalcomplexLorentzian(u, α, β::T, λ, Ω)::Complex{T} where T <: Real

Evaluates the complex Lorentzian function. Returns α*exp(im*β)/(λ + im*(2*π*u - Ω)).
"""
function evalcomplexLorentzian(u, α, β::T, λ, Ω)::Complex{T} where T <: Real
    return α*exp(im*β)/(λ + im*(2*π*u - Ω))
end

"""
    gettimerange(N::Int, fs::T) where T

Returns the time stamps for a sequence, starting at time 0. Returns zero(T):Ts:(N-1)*Ts, Ts = 1/fs.
"""
function gettimerange(N::Int, fs::T) where T
    Ts::T = 1/fs

    return zero(T):Ts:(N-1)*Ts
end

"""
    getwraparoundDFTfreqs(N::Int, fs::T, ν_begin::T) where T

Get the frequency indices for a sequence of length `N` that is sampled at `fs` Hz, with a starting frequency in `ν_begin` Hz.

...
# Inputs:

# Outputs, in the order returned.

- `U_DFT::LinRange{T,Int}`: frequency indices (in Hz) for `fft(x)`, where `length(x)` is `N` and values in `x` were sampled at `fs`.
- `U_y::Vector{T}`: frequency indices (in Hz) for the shited version of `fft(x)` that has the first frequency larger than `ν_begin` as the first entry. When fitting 1D 1H NMR data in the frequency domain, `U_y` is used instead of `U_DFT` with `ν_begin` being the 0 ppm frequency in Hz. One can then convert `U_y`, which is in Hz, to ppm, by using a conversion function appropriate for this experiment. See the `hz2ppmfunc` output from `loadspectrum()`.
- `U_inds::Vector{Int}`: `ff(x)[U_inds]` would give the correct spectrum values that correspond to the frequency indices in `U_y`.

Note that `U_y` doesn't contain the same numerical values as `U_DFT[U_inds]` because the latter is not neccessarily a monotonically increasing sequence, but `U_y` is. They both describe the same discrete-time frequencies though.
, , , , ,

# Examples
See load_experiment.jl in the /examples folder for an example.
"""
function getwraparoundDFTfreqs(N::Int, fs::T, ν_begin::T) where T

    U0 = getDFTfreqrange(N, fs)
    out, inds = wrapfreqrange(U0, ν_begin, fs)

    return U0, out, inds
end

# Assumes U0 is sorted in ascending order.
function wrapfreqrange(U0, ν_begin::T, fs::T) where T <: Real

    N = length(U0)
    ind = findfirst(xx->(xx>ν_begin), U0)

    # TODO handle these exceptions with more grace.
    @assert typeof(ind) == Int
    @assert ind <= N

    out = Vector{T}(undef, N)
    #out[1:ind] = U0[1:ind] .+ fs
    #out[ind+1:end] = U0[ind+1:end]

    M = N-ind
    out[1:M] = U0[ind+1:end]
    out[M+1:end] = U0[1:ind] .+ fs

    inds = collect(1:N)
    inds[1:M] = collect(ind+1:N)
    inds[M+1:end] = collect(1:ind)

    return out, inds
end

function getDFTfreqrange(N::Int, fs::T)::LinRange{T} where T
    a = zero(T)
    b = fs-fs/N

    return LinRange(a, b, N)
end

  # case 1D.
function computeDTFTch3eq29(h, u::T, Λ)::Complex{T} where T <: Real

    # debug_array = zeros(Complex{T}, length(Λ))
    # debug_array2 = zeros(Complex{T}, length(Λ))
    # fs = 1/(Λ[2]-Λ[1])

    running_sum = zero(T)
    for i = 1:length(Λ)
        x = Λ[i]

        running_sum += h[i]*exp(-im*2*π*u*x)

        # # debug.
        # debug_array[i] = h[i]*exp(-im*2*π*u*x)
        #
        # #debug_array2[i] = h[i]*exp(-im*2*π*u*x +im*35*2*pi)
        # debug_array2[i] = h[i]*exp(-im*2*π*(u+3*fs)*x)
    end

    #println("debug: norm(debug_array -debug_array2 ) = ", norm(debug_array -debug_array2 ) )
    return running_sum
end


text:
module TestBroadcasting

using Test
using InitialValues

OPS = [*, +, |, &, min, max, Base.add_sum, Base.mul_prod]

@testset for op in OPS
    @test op.(InitialValue(op), 1:3) == 1:3
end

@testset "binary op" begin
    @test InitialValue(+) .+ (1:3) == 1:3
    @test InitialValue(*) .* (1:3) == 1:3
    @test InitialValue(|) .| (1:3) == 1:3
    @test InitialValue(&) .& (1:3) == 1:3
end

end  # module


text:
# Ported from: https://github.com/juj/RectangleBinPack/blob/master/GuillotineBinPack.cpp

@enum FreeRectChoiceHeuristic begin
    RectBestAreaFit #< -BAF
    RectBestShortSideFit #< -BSSF
    RectBestLongSideFit #< -BLSF
    RectWorstAreaFit #< -WAF
    RectWorstShortSideFit #< -WSSF
    RectWorstLongSideFit #< -WLSF
end

@enum GuillotineSplitHeuristic begin
    SplitShorterLeftoverAxis #< -SLAS
    SplitLongerLeftoverAxis #< -LLAS
    SplitMinimizeArea #< -MINAS, Try to make a single big rectangle at the expense of making the other small.
    SplitMaximizeArea #< -MAXAS, Try to make both remaining rectangles as even-sized as possible.
    SplitShorterAxis #< -SAS
    SplitLongerAxis #< -LAS
end

struct GuillotinePacker{T}
    bin_width::Int
    bin_height::Int
    merge::Bool
    rect_choice::FreeRectChoiceHeuristic
    split_method::GuillotineSplitHeuristic
    free_rectangles::Vector{Rect2D{T}}
    used_rectangles::Vector{Rect2D{T}}

end

function GuillotinePacker(width::T, height::T; merge::Bool=true,
                rect_choice::FreeRectChoiceHeuristic=RectBestAreaFit,
                split_method::GuillotineSplitHeuristic=SplitMinimizeArea) where T

    packer = GuillotinePacker(width, height, merge, rect_choice, split_method,
                              Rect2D{T}[], Rect2D{T}[])
    push!(packer.free_rectangles, Rect2D(T(0), T(0), width, height))
    return packer
end

function Base.push!(packer::GuillotinePacker, rects::AbstractVector{Rect2D{T}}) where T
    rects = copy(rects)
    free_rectangles = packer.free_rectangles
    used_rectangles = packer.used_rectangles
    rect_choice = packer.rect_choice
    split_method = packer.split_method
    # Remember variables about the best packing choice we have made so far during the iteration process.
    best_free_rect = 0
    best_rect = 0
    best_flipped = false
    # Pack rectangles one at a time until we have cleared the rects array of all rectangles.
    # rects will get destroyed in the process.
    while length(rects) > 0
        # Stores the penalty score of the best rectangle placement - bigger=worse, smaller=better.
        best_score = typemax(Int)

        for i in 1:length(free_rectangles)
            force_jump = false
            rect_i = free_rectangles[i]
            for j in 1:length(rects)
                # If this rectangle is a perfect match, we pick it instantly.
                if all(widths(rects[j]) .== widths(rect_i))
                    best_free_rect = i
                    best_rect = j
                    best_flipped = false
                    best_score = typemin(Int)
                    force_jump = true # Force a jump out of the outer loop as well - we got an instant fit.
                    break
                # If flipping this rectangle is a perfect match, pick that then.
                elseif (height(rects[j]) == width(rect_i) && width(rects[j]) == height(rect_i))
                    best_free_rect = i
                    best_rect = j
                    best_flipped = true
                    best_score = typemin(Int)
                    force_jump = true # Force a jump out of the outer loop as well - we got an instant fit.
                    break
                # Try if we can fit the rectangle upright.
                elseif all(widths(rects[j]) .<= widths(rect_i))
                    score = score_by_heuristic(width(rects[j]), height(rects[j]), rect_i, rect_choice)
                    if (score < best_score)
                        best_free_rect = i
                        best_rect = j
                        best_flipped = false
                        best_score = score
                    end
                # If not, then perhaps flipping sideways will make it fit?
                elseif (height(rects[j]) <= width(rect_i) && width(rects[j]) <= height(rect_i))
                    score = score_by_heuristic(height(rects[j]), width(rects[j]), rect_i, rect_choice)
                    if (score < best_score)
                        best_free_rect = i
                        best_rect = j
                        best_flipped = true
                        best_score = score
                    end
                end
            end
            force_jump && break
        end

        # If we didn't manage to find any rectangle to pack, abort.
        (best_score == typemax(Int)) && return nothing

        # Otherwise, we're good to go and do the actual packing.
        xy = minimum(free_rectangles[best_free_rect])
        if best_flipped
            new_node = Rect2D(xy, reverse(widths(rects[best_rect])))
        else
            new_node = Rect2D(xy, widths(rects[best_rect]))
        end

        # Remove the free space we lost in the bin.
        split_freerect_by_heuristic(packer, free_rectangles[best_free_rect], new_node, split_method)
        splice!(free_rectangles, best_free_rect)

        # Remove the rectangle we just packed from the input list.
        splice!(rects, best_rect)

        # Perform a Rectangle Merge step if desired.
        if packer.merge
            MergeFreeList(packer)
        end

        # Remember the new used rectangle.
        push!(used_rectangles, new_node)
    end
end

function Base.push!(packer::GuillotinePacker, rect::Rect2D)
    w, h = widths(rect)
    free_rectangles = packer.free_rectangles
    used_rectangles = packer.used_rectangles
    # Find where to put the new rectangle.
    new_rect, free_node_idx = FindPositionForNewNode(packer, w, h, packer.rect_choice)

    # Abort if we didn't have enough space in the bin.
    if (height(new_rect) == 0)
        return nothing
    end

    # Remove the space that was just consumed by the new rectangle.
    split_freerect_by_heuristic(packer, free_rectangles[free_node_idx], new_rect, packer.split_method)
    splice!(free_rectangles, free_node_idx)

    # Perform a Rectangle Merge step if desired.
    if packer.merge
        MergeFreeList(packer)
    end

    # Remember the new used rectangle.
    push!(used_rectangles, new_rect)

    return new_rect
end

#/ Computes the ratio of used surface area to the total bin area.
function Occupancy(packer::GuillotinePacker)
    used_rectangles = packer.used_rectangles
    #TODO The occupancy rate could be cached/tracked incrementally instead
    # of looping through the list of packed rectangles here.
    surface_area = sum((prod(widths(rect)) for rect in used_rectangles))
    return surface_area / (packer.bin_width * packer.bin_height)

end

#/ Returns the heuristic score value for placing a rectangle of size width*height into free_rect. Does not try to rotate.
function score_by_heuristic(width, height, free_rect::Rect2D, rect_choice::FreeRectChoiceHeuristic)
    rect_choice == RectBestAreaFit && return ScoreBestAreaFit(width, height, free_rect)
    rect_choice == RectBestShortSideFit && return ScoreBestShortSideFit(width, height, free_rect)
    rect_choice == RectBestLongSideFit && return ScoreBestLongSideFit(width, height, free_rect)
    rect_choice == RectWorstAreaFit && return ScoreWorstAreaFit(width, height, free_rect)
    rect_choice == RectWorstShortSideFit && return ScoreWorstShortSideFit(width, height, free_rect)
    rect_choice == RectWorstLongSideFit && return ScoreWorstLongSideFit(width, height, free_rect)
end

function ScoreBestAreaFit(w, h, free_rect::Rect2D)
    return width(free_rect) * height(free_rect) - w * h
end

function ScoreBestShortSideFit(w, h, free_rect::Rect2D)
    leftoverHoriz = abs(width(free_rect) - w)
    leftoverVert = abs(height(free_rect) - h)
    leftover = min(leftoverHoriz, leftoverVert)
    return leftover
end

function ScoreBestLongSideFit(w, h, free_rect::Rect2D)
    leftoverHoriz = abs(width(free_rect) - w)
    leftoverVert = abs(height(free_rect) - h)
    leftover = max(leftoverHoriz, leftoverVert)
    return leftover

end

function ScoreWorstAreaFit(width, height, free_rect::Rect2D)
    return -ScoreBestAreaFit(width, height, free_rect)
end

function ScoreWorstShortSideFit(width, height, free_rect::Rect2D)
    return -ScoreBestShortSideFit(width, height, free_rect)
end

function ScoreWorstLongSideFit(width, height, free_rect::Rect2D)
    return -ScoreBestLongSideFit(width, height, free_rect)
end

function FindPositionForNewNode(packer::GuillotinePacker, w, h, rect_choice::FreeRectChoiceHeuristic)
    free_rectangles = packer.free_rectangles
    best_node = Rect2D(0,0,0,0)
    free_node_idx = 0
    best_score = typemax(Int)
    #/ Try each free rectangle to find the best one for placement.
    for i in 1:length(free_rectangles)
        rect_i = free_rectangles[i]
        # If this is a perfect fit upright, choose it immediately.
        if (w == width(rect_i) && h == height(rect_i))
            best_node = Rect2D(minimum(rect_i), Vec(w, h))
            free_node_idx = i
            break
        # If this is a perfect fit sideways, choose it.
        elseif (h == width(rect_i) && w == height(rect_i))
            best_node = Rect2D(minimum(rect_i), Vec(h, w))
            best_score = typemin(Int)
            free_node_idx = i
            break
            # Does the rectangle fit upright?
        elseif (w <= width(rect_i) && h <= height(rect_i))
            score = score_by_heuristic(w, h, rect_i, rect_choice)
            if (score < best_score)
                best_node = Rect2D(minimum(rect_i), Vec(w, h))
                best_score = score
                free_node_idx = i
            end
            # Does the rectangle fit sideways?
        elseif (h <= width(rect_i) && w <= height(rect_i))
            score = score_by_heuristic(h, w, rect_i, rect_choice)
            if (score < best_score)
                best_node = Rect2D(minimum(rect_i), Vec(h, w))
                best_score = score
                free_node_idx = i
            end
        end
    end
    return best_node, free_node_idx
end

function split_freerect_by_heuristic(packer::GuillotinePacker, free_rect::Rect2D, placed_rect::Rect2D, method::GuillotineSplitHeuristic)
    # Compute the lengths of the leftover area.
    w, h = widths(free_rect) .- widths(placed_rect)

    # Placing placed_rect into free_rect results in an L-shaped free area, which must be split into
    # two disjorectangles. This can be achieved with by splitting the L-shape using a single line.
    # We have two choices: horizontal or vertical.

    # Use the given heuristic to decide which choice to make.
    split_horizontal = false
    if method == SplitShorterLeftoverAxis
        # Split along the shorter leftover axis.
        split_horizontal = (w <= h)
    elseif method == SplitLongerLeftoverAxis
        # Split along the longer leftover axis.
        split_horizontal = (w > h)
    elseif method == SplitMinimizeArea
        # Maximize the larger area == minimize the smaller area.
        # Tries to make the single bigger rectangle.
        split_horizontal = (width(placed_rect) * h > w * height(placed_rect))
    elseif method == SplitMaximizeArea
        # Maximize the smaller area == minimize the larger area.
        # Tries to make the rectangles more even-sized.
        split_horizontal = (width(placed_rect) * h <= w * height(placed_rect))
    elseif method == SplitShorterAxis
        # Split along the shorter total axis.
        split_horizontal = (width(free_rect) <= height(free_rect))
    elseif method == SplitLongerAxis
        # Split along the longer total axis.
        split_horizontal = (width(free_rect) > height(free_rect))
    else
        error("Method unknown: $(method)")
    end
    # Perform the actual split.
    split_freerect_along_axis(packer, free_rect, placed_rect, split_horizontal)
end

#/ This function will add the two generated rectangles into the free_rectangles array. The caller is expected to
#/ remove the original rectangle from the free_rectangles array after that.
function split_freerect_along_axis(packer::GuillotinePacker, free_rect::Rect2D, placed_rect::Rect2D, split_horizontal::Bool)
    free_rectangles = packer.free_rectangles
    # Form the two new rectangles.
    x, y = minimum(free_rect)
    if split_horizontal
        bottom = Rect2D(x, y + height(placed_rect), width(free_rect), height(free_rect) - height(placed_rect))
        right = Rect2D(x + width(placed_rect), y, width(free_rect) - width(placed_rect), height(placed_rect))
    else # Split vertically
        bottom = Rect2D(x, y + height(placed_rect), width(placed_rect), height(free_rect) - height(placed_rect))
        right = Rect2D(x + width(placed_rect), y, width(free_rect) - width(placed_rect), height(free_rect))
    end

    # Add the new rectangles into the free rectangle pool if they weren't degenerate.
    if (width(bottom) > 0 && height(bottom) > 0)
        push!(free_rectangles, bottom)
    end
    if (width(right) > 0 && height(right) > 0)
        push!(free_rectangles, right)
    end
end

function MergeFreeList(packer::GuillotinePacker)
    free_rectangles = packer.free_rectangles
    # Do a Theta(n^2) loop to see if any pair of free rectangles could me merged into one.
    # Note that we miss any opportunities to merge three rectangles into one. (should call this function again to detect that)
    i = 0; j = 0
    while i < length(free_rectangles)
        i += 1; j = i + 1
        rect_i = free_rectangles[i]
        while j < length(free_rectangles)
            j += 1
            if (width(rect_i) == width(free_rectangles[j]) && minimum(rect_i)[1] == minimum(free_rectangles[j])[1])
                if (minimum(rect_i)[2] == maximum(free_rectangles[j])[2])
                    new_y = minimum(rect_i)[2] - height(free_rectangles[j])
                    new_h = height(rect_i) + height(free_rectangles[j])
                    free_rectangles[i] = Rect2D(minimum(rect_i)[1], new_y, Vec(width(rect_i), new_h))
                    splice!(free_rectangles, j)
                    j -= 1
                elseif (maximum(rect_i)[2] == minimum(free_rectangles[j])[2])
                    new_h = height(rect_i) + height(free_rectangles[j])
                    free_rectangles[i] = Rect2D(minimum(rect_i), Vec(width(rect_i), new_h))
                    splice!(free_rectangles, j)
                    j -= 1
                end
            elseif (height(rect_i) == height(free_rectangles[j]) && minimum(rect_i)[2] == minimum(free_rectangles[j])[2])
                if (minimum(rect_i)[1] == minimum(free_rectangles[j])[1] + width(free_rectangles[j]))
                    new_x = minimum(rect_i)[1] - width(free_rectangles[j])
                    new_w = width(rect_i) + width(free_rectangles[j])
                    free_rectangles[i] = Rect2D(new_x, minimum(rect_i)[2], new_w, height(rect_i))
                    splice!(free_rectangles, j)
                    j -= 1
                elseif maximum(rect_i)[1] == minimum(free_rectangles[j])[1]
                    new_w = Vec(width(rect_i) + width(free_rectangles[j]), height(rect_i))
                    free_rectangles[i] = Rect2D(minimum(rect_i), new_w)
                    splice!(free_rectangles, j)
                    j -= 1
                end
            end
        end
    end
end


text:
"""
    MinMax()
MinMax Scaler. This scaler scale data as:
``\\tilde{\\boldsymbol{x}} = \\frac{\\boldsymbol{x}-min(\\boldsymbol{x})}{max(\\boldsymbol{x})-min(\\boldsymbol{x})}``

# Example
```jldoctest preprocessing
julia> scaler = MinMax()
MinMax(Float64[])

julia> fit!(scaler, x)
2×2 Matrix{Float64}:
  1.39855   0.954091
 -1.99789  -3.27656

julia> transform(scaler, x)
20×2 Matrix{Float64}:
 0.597368  0.879853
 0.591301  0.826179
 0.88601   0.942311
 0.992432  1.0
 0.324455  0.710522
 0.877757  0.886514
 1.0       0.927088
 0.199996  0.561398
 0.592578  0.831614
 0.675601  0.888666
 0.488471  0.749707
 0.715552  0.866175
 0.841559  0.908526
 0.442398  0.807779
 0.871787  0.905139
 0.179268  0.524104
 0.72641   0.944478
 0.607941  0.895508
 0.153707  0.43407
 0.0       0.0
```
"""
mutable struct MinMax
    p::AbstractVecOrMat
    MinMax() = new(Array{Float64}(undef, 0))
end

function fit!(scaler::MinMax, x; dims=1)
    scaler.p = vcat(maximum(x, dims=dims), minimum(x, dims=dims))
end

mms(x, ma, mi) = @. (x-mi) / (ma - mi)

function transform(scaler::MinMax, x; dims=1)
    p = scaler.p
    check_size(x, p, dims)
    y = similar(x, Float32)
    if dims == 1
        for i in 1 : size(x, 2)
            y[:, i] = mms(x[:, i], p[:, i]...)
        end
    elseif dims == 2
        for i in 1 : size(x, 1)
            y[i, :] = mms(x[i, :], p[:, i]...)
        end
    end
    return y
end

function fit_transform!(scaler::MinMax, x; dims=1)
    fit!(scaler, x, dims=dims)
    transform(scaler, x, dims=dims)
end

imms(x, ma, mi) = @. x*(ma-mi)+mi

function inv_transform(scaler::MinMax, x; dims=1)
    p = scaler.p
    check_size(x, p, dims)
    y = similar(x, Float32)
    if dims == 1
        for i in 1 : size(x, 2)
            y[:, i] = imms(x[:, i], p[:, i]...)
        end
    elseif dims == 2
        for i in 1 : size(x, 1)
            y[i, :] = imms(x[i, :], p[:, i]...)
        end
    end
    return y
end


text: include("./KroneckerProducts.jl")


text:
@testset "overloads" begin
    rng = MersenneTwister(123456)

    k1 = LinearKernel()
    k2 = SqExponentialKernel()
    k3 = RationalQuadraticKernel()

    for (op, T) in ((+, KernelSum), (*, KernelProduct), (⊗, KernelTensorProduct))
        if T === KernelTensorProduct
            v2_1 = rand(rng, 2)
            v2_2 = rand(rng, 2)
            v3_1 = rand(rng, 3)
            v3_2 = rand(rng, 3)
            v4_1 = rand(rng, 4)
            v4_2 = rand(rng, 4)
        else
            v2_1 = v3_1 = v4_1 = rand(rng, 3)
            v2_2 = v3_2 = v4_2 = rand(rng, 3)
        end
        k = T(k1, k2)

        @test op(k1, k2)(v2_1, v2_2) == k(v2_1, v2_2)
        @test op(k, k3)(v3_1, v3_2) == T((k1, k2, k3))(v3_1, v3_2)
        @test op(k3, k)(v3_1, v3_2) == T((k3, k1, k2))(v3_1, v3_2)
        @test op(k, k)(v4_1, v4_2) == T((k1, k2, k1, k2))(v4_1, v4_2)
        @test op(k1, k2) == T([k1, k2]) == T((k1, k2))

        @test op(T([k1, k2]), T([k2, k1])).kernels == [k1, k2, k2, k1]
        @test op(T([k1, k2]), k3).kernels == [k1, k2, k3]
        @test op(k3, T([k1, k2])).kernels == [k3, k1, k2]

        @test op(T((k1, k2)), T((k2, k1))).kernels == (k1, k2, k2, k1)
        @test op(T((k1, k2)), k3).kernels == (k1, k2, k3)
        @test op(k3, T((k1, k2))).kernels == (k3, k1, k2)
    end
end


text:
"""
    to_symbol(input::String)

Cleans up `input` string and converts it into a symbol
"""
function to_symbol(input::String)
    fixed = replace(input, r"[^\w\ \-\_]"=>"")
    fixed = replace(fixed, r"[\ \-\_]+"=>"_")
    Symbol(replace(fixed, r"^[\d]"=>s"_\g<0>"))
end

"""
    myendian()

Returns the TIFF endian byte order expected if its endianness matches the host's
"""
function myendian()
    if ENDIAN_BOM == 0x04030201
        return 0x4949
    elseif ENDIAN_BOM == 0x01020304
        return 0x4d4d
    end
end

"""
    check_bswap(io::Union{Stream, IOStream})

Check endianness of TIFF file to see if we need to swap bytes
"""
function check_bswap(io::Union{Stream, IOStream})
    seekstart(io)
    endianness = read(io, UInt16)
    # check if we need to swap byte order
    need_bswap = endianness != myendian()

    tiff_version = Array{UInt8}(undef, 2)
    read!(io, tiff_version)
    (tiff_version != [0x2a, 0x00] && tiff_version != [0x00, 0x2a]) && error("Big-TIFF files aren't supported yet")

    return need_bswap
end

"""
Extract the name of the file backing a stream
"""
extract_filename(io::IOStream) = split(io.name, " ")[2][1:end-1]
extract_filename(io::Stream) = io.filename


text:
using Test, YAAD, YAAD.TestUtils
using LinearAlgebra

@testset "linalg.jl" begin

@test gradcheck(tr, Variable(rand(10, 10)))
@test gradcheck(*, Variable(rand(10, 10)), Variable(rand(10, 10)))
@test gradcheck(transpose, Variable(rand(10, 10)))

end


text:
# ebi_blast.jl - EBI BLAST REST Interface
# ...


text:
struct Interval{K}
    lo::K
    hi::K
    function Interval{K}(l::K, h::K) where {K}
        @assert !(h < l) "Invalid interval"
        new{K}(l, h)
    end
end

function Base.getindex(itv::Interval, i::Int)
    if i == 1
        return itv.lo
    elseif i == 2
        return itv.hi
    else
        error("Invalid index for Interval")
    end
end

Interval(l::K, h::K) where K = Interval{K}(l, h)

Base.convert(::Type{Interval{K}}, t::Tuple{K, K}) where K = Interval(t...)

Base.promote_rule(::Type{Interval{T}}, ::Type{Interval{S}}) where {T, S} =
    Interval{promote_type(T, S)}

Base.show(io::IO, i::Interval) = Base.print(io, "($(i.lo), $(i.hi))")

function overlaps(i1::Union{Tuple{K, K}, Interval{K}}, i2::Interval{K}) where K
    (i1, i2) = i1[1] < i2[1] ? (i1, i2) : (i2, i1)
    return !(i2[1] < i1[1]) && !(i1[2] < i2[1])
end

Base.isless(i1::Interval{K1}, i2::Interval{K2}) where {K1, K2} =
    Base.isless(promote(i1, i2)...)

Base.isless(i1::Interval{K}, i2::Interval{K}) where {K} =
    (i1.lo < i2.lo) || (!(i1.lo < i2.lo || i2.lo < i1.lo) && (i1.hi < i2.hi))

mutable struct IntervalKey{K}
    i::Interval{K}
    submax::K
    IntervalKey{K}(i::Interval{K}) where {K} = new{K}(i, i.hi)
end

Base.show(io::IO, k::IntervalKey) =
    print(io, "($(k.i.lo), $(k.i.hi), $(k.submax))")

IntervalKey(lo::K, hi::K) where {K} = IntervalKey{K}(Interval(lo, hi))
IntervalKey(i::Interval{K}) where {K} = IntervalKey(i.lo, i.hi)

const IntervalTree{K, V} = RBTree{IntervalKey{K}, V}
const IntervalNode{K, V} = RBNode{IntervalKey{K}, V}

Base.isless(i1::IntervalKey{K}, i2::IntervalKey{K}) where {K} =
    Base.isless(i1.i, i2.i)

function update_submax!(t::IntervalTree{K, V},
                        n::IntervalNode{K, V}) where {K, V}
    isnil(t, n)   && return
    isnil(t, n.l) && isnil(t, n.r) && (n.k.submax = n.k.i.hi; return)
    isnil(t, n.l) && (n.k.submax = max(n.k.i.hi, n.r.k.submax); return)
    isnil(t, n.r) && (n.k.submax = max(n.k.i.hi, n.l.k.submax); return)
    (n.k.submax = max(n.k.i.hi, n.l.k.submax, n.r.k.submax); return)
end

#=
Solution similar to as suggested here:

https://algs4.cs.princeton.edu/93intersection/
=#

@inline function _intersect(t::IntervalTree{K, V},
                            x::IntervalNode{K, V},
                            i::Union{Tuple{K, K}, Interval{K}},
                            nodes::Vector{IntervalNode{K, V}}) where {K, V}
    isnil(t, x) && return false
    lfound = !isnil(t, x.l) && !(x.l.k.submax < i[1]) &&
        _intersect(t, x.l, i, nodes)
    found = overlaps(i, x.k.i)
    found && push!(nodes, x)
    rfound = (lfound || isnil(t, x.l) || x.l.k.submax < i[1]) &&
        _intersect(t, x.r, i, nodes)
    return lfound || found || rfound
end

@inline function Base.intersect(t::IntervalTree{K, V},
                                i::Interval{K}) where {K, V}
    nodes = Vector{IntervalNode{K, V}}()
    _intersect(t, t.root, i, nodes)
    retval = similar(nodes, Pair{Interval{K}, V})
    for i = 1:length(nodes)
        x = nodes[i]
        retval[i] = x.k.i => x.v
    end
    return retval
end

@inline function Base.intersect(t::IntervalTree{K, V},
                                i::Tuple{K, K}) where {K, V}
    nodes = Vector{IntervalNode{K, V}}()
    _intersect(t, t.root, i, nodes)
    length(nodes) == 0 && return nothing
    x = nodes[1]
    return x.k.i => x.v
end

@inline function intersects(t::IntervalTree{K, V}, i::Interval{K}) where {K, V}
    x = t.root
    while !isnil(t, x) && !overlaps(i, x.k.i)
        x = isnil(t, x.l) || (x.l.k.submax < i.lo) ? x.r : x.l
    end
    return !isnil(t, x)
end

_search(t::IntervalTree{K, V}, x::IntervalNode{K, V},
        i::Interval{K}) where {K, V} = _search(t, x, IntervalKey(i))

@inline function Base.get(t::IntervalTree{K, V},
                          i::Interval{K}, v::V) where {K, V}
    isempty(t) && return v
    n, d = _search(t, t.root, IntervalKey(i))
    d != 0 && return v
    return n.v
end

@inline function Base.get!(t::IntervalTree{K, V},
                           i::Interval{K}, v::V) where {K, V}
    if isempty(t)
        insert!(t, i, v)
        return v
    end
    n, d = _search(t, t.root, IntervalKey(i))
    d == 0 && return n.v
    insert!(t, i, v)
    return v
end

@inline function Base.setindex!(t::IntervalTree{K, V},
                                v::V, i::Interval{K}) where {K, V}
    if isempty(t)
        insert!(t, i, v)
        return v
    end
    n, d = _search(t, t.root, IntervalKey(i))
    d == 0 && delete!(t, i)
    insert!(t, i, v)
    return v
end

@inline function Base.getindex(t::IntervalTree{K, V},
                               i::Interval{K}) where {K, V}
    n, d = _search(t, t.root, IntervalKey(i))
    d != 0 && error("Index $i not found.")
    return n.v
end

@inline function left_rotate!(t::IntervalTree{K, V},
                              x::IntervalNode{K, V}) where {K, V}
    y = x.r
    x.r = y.l
    !isnil(t, y.l) && (y.l.p = x)
    y.p = x.p
    if isnil(t, x.p)
        t.root = y
    elseif x === x.p.l
        x.p.l = y
    else
        x.p.r = y
    end
    y.l = x
    x.p = y

    update_submax!(t, x)
    update_submax!(t, y)
    return
end

@inline function right_rotate!(t::IntervalTree{K, V},
                               y::IntervalNode{K, V}) where {K, V}
    x = y.l
    y.l = x.r
    !isnil(t, x.r) && (x.r.p = y)
    x.p = y.p
    if isnil(t, y.p)
        t.root = x
    elseif y === y.p.r
        y.p.r = x
    else
        y.p.l = x
    end
    x.r = y
    y.p = x

    update_submax!(t, y)
    update_submax!(t, x)
    return
end

function Base.insert!(t::IntervalTree{K, V},
                      k::IntervalKey{K},
                      v::V) where {K, V}
    z = RBNode(t, k, v)
    y = t.nil
    x = t.root

    while x !== t.nil
        y = x
        x = k < x.k ? x.l :
            !t.unique || (x.k < k) ? x.r :
            t.unique && error("Key $k already exists.")
    end
    z.p = y
    if y === t.nil
        t.root = z
    elseif k < y.k
        y.l = z
    else
        y.r = z
    end
    z.red = true
    zz = z
    # Ensure new node range is part of the tree subtree max calculations
    while !isnil(t, zz)
        update_submax!(t, zz)
        zz = zz.p
    end
    _insert_fixup!(t, z)
    t.n += 1
    return t
end

function Base.insert!(t::IntervalTree{K, V}, lo::K, hi::K, v::V) where {K, V}
    ik = IntervalKey(lo, hi)
    Base.insert!(t, ik, v)
    return t
end

function Base.insert!(t::IntervalTree{K, V}, i::Interval{K}, v::V) where {K, V}
    ik = IntervalKey(i)
    Base.insert!(t, ik, v)
    return t
end

function Base.iterate(it::Iterator{IntervalKey{K}, V,
                                   IntervalTree{K, V}},
                      n::IntervalNode{K, V}) where {K, V}
    n === it.to && return nothing
    return ((n.k.i => n.v), _successor(n, it.tree))
end

Base.eltype(it::Iterator{IntervalKey{K}, V,
                         IntervalTree{K, V}}) where {K, V} = Pair{Interval{K}, V}

function Base.minimum(t::IntervalTree)
    isempty(t) && error("Empty tree cannot have a minimum")
    n = _minimum(t.root, t)
    return n.k.i => n.v
end

function Base.maximum(t::IntervalTree)
    isempty(t) && error("Empty tree cannot have a maximum")
    n = _maximum(t.root, t)
    return n.k.i => n.v
end

@inline function _delete!(t::IntervalTree{K, V},
                          z::IntervalNode{K, V}) where {K, V}
    y = z
    y_is_red = y.red
    if z.l === t.nil
        x = z.r
        _transplant!(t, z, z.r)
    elseif z.r === t.nil
        x = z.l
        _transplant!(t, z, z.l)
    else
        y = _minimum(z.r, t)
        y_is_red = y.red
        x = y.r
        if y.p === z
            x.p = y
        else
            _transplant!(t, y, y.r)
            y.r = z.r
            y.r.p = y
        end
        _transplant!(t, z, y)
        y.l = z.l
        y.l.p = y
        y.red = z.red
    end
    # Ensure that the delete of node is reflected up.
    tx = x
    while true
        update_submax!(t, tx)
        tx = tx.p
        isnil(t, tx) && break
    end
    if !y_is_red
        _delete_fixup!(t, x)
    end
    return z
end

@inline function Base.delete!(t::IntervalTree{K, V}, i::Interval{K}) where {K, V}
    isempty(t) && error("Cannot delete from empty tree")
    n, d = _search(t, t.root, i)
    d != 0 && return nothing
    _delete!(t, n)
    t.n -= 1
    return n.k.i => n.v
end


text:
"""
    dist_matrix(distr1, distr2; dist_met = kmer_seeded_edit_dist)

distances[i, j] is distance from distr1[i] to distr2[j]
"""
function dist_matrix(distr1, distr2; dist_met = kmer_seeded_edit_dist)
    if distr1 == distr2
        return symmetric_dist_matrix(distr1, dist_met = dist_met)
    end

    distances = zeros(length(distr1), length(distr2))
    for i in 1:length(distr1)
        for j in 1:length(distr2)
            distances[i, j] = dist_met(distr1[i], distr2[j])
        end
    end
    return distances
end

"""
    dist_matrix_mt(distr1, distr2; dist_met = kmer_seeded_edit_dist)

dist_matrix but multithreaded instead
"""
function dist_matrix_mt(distr1, distr2; dist_met = kmer_seeded_edit_dist)
    if distr1 == distr2
        return symmetric_dist_matrix_mt(distr1, dist_met = dist_met)
    end

    distances = zeros(length(distr1), length(distr2))
    @sync @parallel for i in 1:length(distr1)
        for j in 1:length(distr2)
            distances[i, j] = dist_met(distr1[i], distr2[j])
        end
    end
    return distances
end

"""
    symmetric_dist_matrix(distr; dist_met = kmer_seeded_edit_dist)

similar to dist_matrix but where distr1 == distr2. Automatically called by dist_matrix
"""
function symmetric_dist_matrix(distr; dist_met = kmer_seeded_edit_dist)
    distances = zeros(length(distr), length(distr))
    for i in 1:length(distr)
        for j in i+1:length(distr)
            distances[i, j] = distances[j, i]  = dist_met(distr[i], distr[j])
        end
    end
    return distances
end

"""
symmetric_dist_matrix_mt(distr; dist_met = kmer_seeded_edit_dist)

similar to dist_matrix_mt but where distr1 == distr2. Automatically called by dist_matrix_mt
"""
function symmetric_dist_matrix_mt(distr; dist_met = kmer_seeded_edit_dist)
    distances = zeros(length(distr), length(distr))
    @sync @parallel for i in 1:length(distr)
        for j in i+1:length(distr)
            distances[i, j] = distances[j, i]  = dist_met(distr[i], distr[j])
        end
    end
    return distances
end

"""
smd_mf(distances_matrix::Array{Float64,2};
                freq1::Vector{Float64}=Float64[],
                freq2::Vector{Float64}=Float64[],
                unbounded_first::Bool=false,
                unbounded_second::Bool=false)

Mutation distance from one population to another
"""
function smd_mf(distances_matrix::Array{Float64,2};
                freq1::Vector{Float64}=Float64[],
                freq2::Vector{Float64}=Float64[],
                unbounded_first::Bool=false,
                unbounded_second::Bool=false)
    if (size(distances_matrix)[1] != length(freq1))
        freq1 = ones(size(distances_matrix)[1])
    end
    if (size(distances_matrix)[2] != length(freq2))
        freq2 = ones(size(distances_matrix)[2])
    end

    # init solver
    m = Model(solver = ClpSolver(SolveType=5))

    # flow is a matrix with how much of each distr1 elem maps to how
    # much of each distr2 elem
    @variable(m, flow[1:length(freq1), 1:length(freq2)] >= 0 )

    # constraint: The occurrence of a element "i" in flow is equal to
    # its frequency (normalized)
    if !unbounded_first
        for i = 1:length(freq1)
            @constraint(m, sum(flow[i, :]) == freq1[i] / sum(freq1))
        end
    end

    # constraint: The occurence of a element "j" in flow is equal to
    # its frequency (normalized)
    if !unbounded_second
        for j = 1:length(freq2)
            @constraint(m, sum(flow[:, j]) == freq2[j] / sum(freq2))
        end
    end

    # minimize sum of distance flows
    @objective(m, Min, sum(flow .* distances_matrix))

    status = solve(m)
    return (getobjectivevalue(m), (getvalue(flow)))
end

"""
smd_sum(distances_matrix::Array{Float64,2};
                freq1::Vector{Float64}=Float64[],
                freq2::Vector{Float64}=Float64[])
"""
function smd_sum(distances_matrix::Array{Float64,2};
                freq1::Vector{Float64}=Float64[],
                freq2::Vector{Float64}=Float64[])
	return smd_mf(distances_matrix, freq1, freq2, unbounded_first = true)[2]
		+ smd_mf(distances_matrix, freq1, freq2, unbounded_second = true)[2]
end

"""
smd_distance_wrapper(distmat::Array{Float64,2}, inds1, inds2)

Wrapper for the smd_mf call
"""
function smd_distance_wrapper(distmat::Array{Float64,2}, inds1, inds2)
    return smd_mf(distmat[inds1, :][:, inds2])[1]
end

"""
fst_wrapper(distmat, inds1, inds2)

Wrapper function
"""
function fst_wrapper(distmat, inds1, inds2)
    l1 = length(inds1)
    l2 = length(inds2)

    within = (sum([distmat[x,y] for x in inds1, y in inds1])/(l1  - 1 ) + (sum([distmat[x,y] for x in inds2, y in inds2])/(l2 - 1)))/(l1 + l2)
    between = mean([distmat[x,y] for x in inds1, y in inds2])
    return (between - within)/between
end

"""
permutation_test(distmat::Array{Float64,2}; l1 = nothing, l2 = nothing, tests=10000, dist_func = smd_distance_wrapper, randvariation=true)
"""
function permutation_test(distmat::Array{Float64,2}; l1 = nothing, l2 = nothing, tests=10000, dist_func = smd_distance_wrapper, randvariation=true)
    #Get the upper right quadrant of the huge matrix

    if l1 == nothing || l2 == nothing
        l1 = trunc(Int, (size(distmat)[1])/2)
        l2 = ceil(Int, (size(distmat)[1])/2)
    end

    group1_selects = 1:l1
    group2_selects = group1_selects + l2

    baseline = dist_func(distmat, group1_selects, group2_selects)

    array_of_distances = []
    #array_of_distances = zeros(Float64, tests)

    for i in 1:tests
        # permuted = shuffle(1:(l1+l2))
        # group1_selects = permuted[1:l1]
        # group2_selects = permuted[l1+1:end]
        group1_selects = shuffle(1:(l1+l2))[1:l1]
        group2_selects = [x for x in 1:(l1+l2) if !(x in group1_selects)]

        dist = dist_func(distmat, group1_selects, group2_selects)
        push!(array_of_distances, dist)
        #array_of_distances[i] = dist
    end

    if randvariation
        baseline += randn() * 0.00000001
        dist_func = [x + (randn() * 0.00000001) for x in array_of_distances]
    end
    return array_of_distances, baseline
    #return length(array_of_distances .> baseline) > 0 ? sum(array_of_distances .> baseline)/length(array_of_distances) : 1
end


text:
using Remark, FileWatching

while true
    Remark.slideshow(@__DIR__; options = Dict("ratio" => "16:9"), credit=false,
    title = "How similar do two strings look? Visual distances in Julia")
    @info "Rebuilt"
    FileWatching.watch_folder(joinpath(@__DIR__, "src"))
end


text:
using JuliaFormatter

handle("format") do data
  @destruct [
    text,
    indent || 4,
    margin || 92,
    always_for_in || false,
    whitespace_typedefs || false,
    whitespace_ops_in_indices || false,
    remove_extra_newlines || false,
    import_to_using || false,
    pipe_to_function_call || false,
    short_to_long_function_def || false,
    always_use_return || false,
    use_YAS_style || false
  ] = data

  style = (use_YAS_style && isdefined(JuliaFormatter, :YASStyle)) ?
    JuliaFormatter.YASStyle() :
    JuliaFormatter.DefaultStyle()

  return Dict(:formattedtext => format_text′(
    text;
    indent = indent,
    margin = margin,
    always_for_in = always_for_in,
    whitespace_typedefs = whitespace_typedefs,
    whitespace_ops_in_indices = whitespace_ops_in_indices,
    remove_extra_newlines = remove_extra_newlines,
    import_to_using = import_to_using,
    pipe_to_function_call = pipe_to_function_call,
    short_to_long_function_def = short_to_long_function_def,
    always_use_return = always_use_return,
    style = style
  ))
end

# HACK: extract keyword arguments of `format_text`; `Base.kwarg_decl` isn't available as of v1.0
const FORMAT_TEXT_KWARGS = let
  ms = collect(methods(format_text))
  filter!(m->m.module==JuliaFormatter, ms)
  m = match(r";(.+)\)", string(first(ms)))
  m === nothing ? Symbol[] : Symbol.(strip.(split(m.captures[1]), Ref((' ', ','))))
end
function format_text′(text; kwargs...)
  # only pass valid keyword arguments to `format_text`
  valid_kwargs_dict = filter(kwargs) do (k, _)
    @static isempty(FORMAT_TEXT_KWARGS) ? false : in(k,FORMAT_TEXT_KWARGS)
  end
  ks = (collect(keys(valid_kwargs_dict))...,)
  vs = (collect(values(valid_kwargs_dict))...,)
  valid_kwargs_nt = NamedTuple{ks}(vs)
  valid_kwargs = pairs(valid_kwargs_nt)

  format_text(text; valid_kwargs...)
end


text:
using CPIMC, CPIMC.PlaneWaves, CPIMC.UniformElectronGas, CPIMC.DefaultUpdates, DataStructures
import CPIMC: orbs, adjacent_kinks_affecting_orbs, kinks_affecting_orbs, prev, next, next_affecting, prev_affecting, τ_prev_affecting, τ_next_affecting, τ_borders, isunaffected, isunaffected_in_interval, time_ordered_orbs, occupations_at, longest_type_1_chain_length, right_type_1_count, kinks_from_periodic_interval, times_from_periodic_interval, Δ, Woffdiag_element, ΔWoffdiag_element, ΔWdiag_element, ΔW_diag, add_orbs, add_orbs!, drop_orbs, drop_orbs!, drop_kinks, drop_kinks!, add_kinks, add_kinks!

S = sphere_with_same_spin(PlaneWave((0,0,0)),dk=1)
a = PlaneWave((-2,0,0))
b = PlaneWave((3,0,0))
c = PlaneWave((0,0,0))
d = PlaneWave((1,0,0))
e = PlaneWave((5,9,9))
f = PlaneWave((1,1,1))

g = PlaneWave(a.vec + e.vec - f.vec, Up)
h = PlaneWave(c.vec + d.vec - g.vec, Up)

sd = Kinks( ImgTime(0.2) => T4(a,b,c,d),
       ImgTime(0.5) => T4(c,d,a,b),
       ImgTime(0.6) => T4(b,a,d,c),
       ImgTime(0.8) => T4(d,c,b,a) )

conf = Configuration(sphere(PlaneWave((0,0,0),Up),dk=1),sd)

@testset "orbs" begin
    @test orbs(T4(a,b,c,d))[1] == a
    @test orbs(T4(a,b,c,d))[2] == b
    @test orbs(T4(a,b,c,d))[3] == c
    @test orbs(T4(a,b,c,d))[4] == d
    @test orbs(T4(a,b,c,d)) == (a,b,c,d)
    @test orbs(T2(a,b)) == (a,b)
end

@testset "adjacent_kinks_affecting_orbs" begin
    @test adjacent_kinks_affecting_orbs(conf, Set([a]), ImgTime(0.0)) == (ImgTime(0.8) => T4(d,c,b,a), ImgTime(0.2) => T4(a,b,c,d))
    @test adjacent_kinks_affecting_orbs(conf, Set([b]), ImgTime(0.1)) == (ImgTime(0.8) => T4(d,c,b,a), ImgTime(0.2) => T4(a,b,c,d))
    @test adjacent_kinks_affecting_orbs(conf, Set([c]), ImgTime(0.2)) == (ImgTime(0.8) => T4(d,c,b,a), ImgTime(0.5) => T4(c,d,a,b))
    @test adjacent_kinks_affecting_orbs(conf, Set([d]), ImgTime(0.3)) == (ImgTime(0.2) => T4(a,b,c,d), ImgTime(0.5) => T4(c,d,a,b))
    @test adjacent_kinks_affecting_orbs(conf, Set([a,b]), ImgTime(0.4)) == (ImgTime(0.2) => T4(a,b,c,d), ImgTime(0.5) => T4(c,d,a,b))
    @test adjacent_kinks_affecting_orbs(conf, Set([b,e]), ImgTime(0.5)) == (ImgTime(0.2) => T4(a,b,c,d), ImgTime(0.6) => T4(b,a,d,c))
    @test adjacent_kinks_affecting_orbs(conf, Set([c,d]), ImgTime(0.6)) == (ImgTime(0.5) => T4(c,d,a,b), ImgTime(0.8) => T4(d,c,b,a))
    @test adjacent_kinks_affecting_orbs(conf, Set([d,a]), ImgTime(0.7)) == (ImgTime(0.6) => T4(b,a,d,c), ImgTime(0.8) => T4(d,c,b,a))
    @test adjacent_kinks_affecting_orbs(conf, Set([a,c,d]), ImgTime(0.8)) == (ImgTime(0.6) => T4(b,a,d,c), ImgTime(0.2) => T4(a,b,c,d))
    @test adjacent_kinks_affecting_orbs(conf, Set([b,c,d]), ImgTime(0.9)) == (ImgTime(0.8) => T4(d,c,b,a), ImgTime(0.2) => T4(a,b,c,d))
    @test adjacent_kinks_affecting_orbs(conf, Set([e]), ImgTime(0.9)) == (nothing,nothing)
end

@testset "kinks_affecting_orbs" begin
    @test kinks_affecting_orbs(conf, Set([a])) == sd
    @test kinks_affecting_orbs(conf, Set([b])) == sd
    @test kinks_affecting_orbs(conf, Set([c])) == sd
    @test kinks_affecting_orbs(conf, Set([d])) == sd
    @test kinks_affecting_orbs(conf, Set([e])) == Kinks{PlaneWave}()
end

@testset "prev" begin
    @test prev(conf.kinks, ImgTime(0.0)) == 4
    @test prev(conf.kinks, ImgTime(0.5)) == 1
    @test prev(conf.kinks, ImgTime(0.1)) == 4
    @test prev(conf.kinks, ImgTime(0.7)) == 3
end

@testset "next" begin
    @test next(conf.kinks, ImgTime(0.1)) == 1
    @test next(conf.kinks, ImgTime(0.5)) == 3
    @test next(conf.kinks, ImgTime(0.0)) == 1
    @test next(conf.kinks, ImgTime(0.8)) == 1
end

@testset "prev_affecting" begin
    @test prev_affecting(conf.kinks, (a,), ImgTime(0.0)) == 4
    @test prev_affecting(conf.kinks, (a,e), ImgTime(0.5)) == 1
    @test prev_affecting(conf.kinks, (b,d), ImgTime(0.1)) == 4
    @test prev_affecting(conf.kinks, (b,d), ImgTime(0.1)) == 4
    @test prev_affecting(conf.kinks, (e,f), ImgTime(0.1)) == 0
end

@testset "next_affecting" begin
    @test next_affecting(conf.kinks, (a,), ImgTime(0.0)) == 1
    @test next_affecting(conf.kinks, (a,e), ImgTime(0.5)) == 3
    @test next_affecting(conf.kinks, (b,d), ImgTime(0.1)) == 1
    @test next_affecting(conf.kinks, (b,d), ImgTime(0.8)) == 1
    @test next_affecting(conf.kinks, (e,f), ImgTime(0.1)) == 0
end

@testset "τ_prev_affecting" begin
    @test τ_prev_affecting(conf.kinks, Set([a]), ImgTime(0.0)) == ImgTime(0.8)
    @test τ_prev_affecting(conf.kinks, Set([a,e]), ImgTime(0.5)) == ImgTime(0.2)
    @test τ_prev_affecting(conf.kinks, Set([b,d]), ImgTime(0.1)) == ImgTime(0.8)
end

@testset "τ_next_affecting" begin
    @test τ_next_affecting(conf.kinks, Set([a]), ImgTime(0.0)) == ImgTime(0.2)
    @test τ_next_affecting(conf.kinks, Set([a,e]), ImgTime(0.5)) == ImgTime(0.6)
    @test τ_next_affecting(conf.kinks, Set([b,d]), ImgTime(0.1)) == ImgTime(0.2)
end

@testset "isunaffected" begin
    @test !isunaffected(conf.kinks, a)
    @test !isunaffected(conf.kinks, b)
    @test !isunaffected(conf.kinks, c)
    @test !isunaffected(conf.kinks, d)
    @test isunaffected(conf.kinks, e)
end

@testset "isunaffected_in_interval" begin
    @test !isunaffected_in_interval(conf.kinks, a, ImgTime(0.1), ImgTime(0.7))
    @test isunaffected_in_interval(conf.kinks, a, ImgTime(0.05), ImgTime(0.15))
    @test isunaffected_in_interval(conf.kinks, a, ImgTime(0.9), ImgTime(0.15))
    @test !isunaffected_in_interval(conf.kinks, b, ImgTime(0.1), ImgTime(0.9))
    @test isunaffected_in_interval(conf.kinks, b, ImgTime(0.9), ImgTime(0.1))
    @test !isunaffected_in_interval(conf.kinks, c, ImgTime(0.1), ImgTime(0.7))
    @test !isunaffected_in_interval(conf.kinks, d, ImgTime(0.1), ImgTime(0.7))
    @test isunaffected_in_interval(conf.kinks, e, ImgTime(0.1), ImgTime(0.7))
end

@testset "time_ordered_orbs(::SortedDict{ImgTime,<:Kink})" begin
    @test time_ordered_orbs(conf.kinks)[1] == a
    @test time_ordered_orbs(conf.kinks)[4] == d
    @test time_ordered_orbs(conf.kinks)[end] == a
    @test time_ordered_orbs(conf.kinks)[end-4] == c
end

@testset "kinks_from_periodic_interval(ck::SortedDict{ImgTime,<:Kink}, τ1, τ2)" begin
    ImgTime(0.2) => T4(a,b,c,d), ImgTime(0.5) => T4(c,d,a,b), ImgTime(0.6) => T4(b,a,d,c), ImgTime(0.8) => T4(d,c,b,a)
    @test kinks_from_periodic_interval(sd, 0.1, 0.6) == Kinks(ImgTime(0.2) => T4(a,b,c,d), ImgTime(0.5) => T4(c,d,a,b))
    @test kinks_from_periodic_interval(sd, 0.7, 0.4) == Kinks(ImgTime(0.2) => T4(a,b,c,d), ImgTime(0.8) => T4(d,c,b,a))
    @test kinks_from_periodic_interval(sd, 0.5, 0.7) == Kinks(ImgTime(0.6) => T4(b,a,d,c))

    @test isempty(kinks_from_periodic_interval(sd, 0.3, 0.4))
    @test kinks_from_periodic_interval(sd, 0.4, 0.3) == sd

    @test isempty(kinks_from_periodic_interval(sd, 0.5, 0.5))
    @test isempty(kinks_from_periodic_interval(Kinks{PlaneWave{3}}(), 0.1, 0.9))
end

@testset "times_from_periodic_interval" begin
    t1 = ImgTime(0.1)
    t2 = ImgTime(0.3)
    @test times_from_periodic_interval(sd, t1, t2) == [ImgTime(0.2)]
    @test times_from_periodic_interval(sd, t2, t1) == [ImgTime(0.5), ImgTime(0.6), ImgTime(0.8)]

    # test for times with no kinks in between
    t1 = ImgTime(0.3)
    t2 = ImgTime(0.4)
    @test times_from_periodic_interval(sd, t1, t2) == []
end

@testset "Δ(τ1::ImgTime,τ2::ImgTime)" begin
    @test Δ(ImgTime(0.3), ImgTime(0.5)) == ImgTime(0.5) - ImgTime(0.3)
    @test Δ(ImgTime(0.8), ImgTime(0.2)) == ImgTime(1) + ImgTime(0.2) - ImgTime(0.8)
    @test iszero( Δ(ImgTime(0.8),ImgTime(0.8)) )

    @test float(Δ(ImgTime(0.3), ImgTime(0.5))) ≈ float(ImgTime(0.2))
end

@testset "Type_1_investigation" begin
    g = PlaneWave(a.vec + e.vec - f.vec, Up)
    h = PlaneWave(c.vec + d.vec - g.vec, Up)

    Type_1_chain = Kinks( ImgTime(0.2) => T4(a,b,c,d),
                          ImgTime(0.5) => T4(f,g,e,a),
                          ImgTime(0.6) => T4(c,d,h,g),
                          ImgTime(0.8) => T4(e,h,b,f) )

    @test (a.vec + b.vec - c.vec - d.vec) == PlaneWave((0,0,0)).vec
    @test (f.vec + g.vec - e.vec - a.vec) == PlaneWave((0,0,0)).vec
    @test (c.vec + d.vec - h.vec - g.vec) == PlaneWave((0,0,0)).vec
    @test (e.vec + h.vec - b.vec - f.vec) == PlaneWave((0,0,0)).vec

    occs = setdiff!(union!(sphere(PlaneWave((0,0,0),Up),dk=1), Set([e,h])), Set([g,f]))
    conf_Type_1 = Configuration(occs,Type_1_chain)
    @test (occupations_at(conf_Type_1, ImgTime(0.9)) == occs)
    @test longest_type_1_chain_length(conf_Type_1.kinks) == 4
    @test right_type_1_count(conf_Type_1.kinks) == 4

    add_kinks!(Type_1_chain, (ImgTime(0.52) => T4(e,a,g,f), ImgTime(0.54) => T4(g,f,a,e)))

    conf_Type_1 = Configuration(occs,Type_1_chain)

    @test length(conf_Type_1.kinks) == 6

    @test longest_type_1_chain_length(conf_Type_1.kinks) == 4
    @test right_type_1_count(conf_Type_1.kinks) == 4

    add_kinks!(Type_1_chain, (ImgTime(0.82) => T4(b,a,d,c), ImgTime(0.84) => T4(c,d,a,b)))

    conf_Type_1 = Configuration(occs,Type_1_chain)
    @test longest_type_1_chain_length(conf_Type_1.kinks) == 3
    @test right_type_1_count(conf_Type_1.kinks) == 4

end

@testset "ΔWoffdiag_element" begin

    mod = UEG()

    ens = CEnsemble(2.0, 5.680898543560106, 7)# θ: 0.125, λ: 0.09945178864947428

    t1 = ImgTime(0.3)
    t2 = ImgTime(0.7)
    τ3 = ImgTime(0.9)

    @test ΔWoffdiag_element(mod, ens, kinks_from_periodic_interval(sd, ImgTime(0.3), ImgTime(0.9)), kinks_from_periodic_interval(sd, ImgTime(0.3), ImgTime(0.7))) ≈ Woffdiag_element(mod, ens, d,c,b,a)
end

@testset "ΔWdiag_element(::Model, ::Ensemble, ::Configuration, i, j, k, l, τ1, τ2)" begin

    mod = UEG()

    i = PlaneWave((0,-4,0))
    j = PlaneWave((0,3,1))
    k = PlaneWave((0,0,1))
    l = PlaneWave((0,-1,0))

    τ1 = ImgTime(0.1)
    τ2 = ImgTime(0.3)
    τ3 = ImgTime(0.4)

    λ = 0.8
    β = 0.02
    N = length(conf.occupations)

    @test ΔWdiag_element(mod, CEnsemble(λ, β, N), conf, i, j, k, l, τ2, τ3) ≈ ΔW_diag(mod, i, j, k, l, occupations_at(conf,τ2)) * (τ3 - τ2) * λ

    @test ΔWdiag_element(mod, CEnsemble(λ, β, N), conf, i, j, k, l, τ1, τ3) ≈ ( ΔW_diag(mod, i, j, k, l, occupations_at(conf,τ1)) * (ImgTime(0.2) - τ1)
                                                                                + ΔW_diag(mod, i, j, k, l, occupations_at(conf,ImgTime(0.2))) * (τ3 - ImgTime(0.2))
                                                                                ) * λ

    @test ΔWdiag_element(mod, CEnsemble(λ, β, N), conf, i, j, k, l, τ3, τ1) ≈ ( ΔW_diag(mod, i, j, k, l, occupations_at(conf,τ3)) * (ImgTime(0.5) - τ3)
                                                                                + ΔW_diag(mod, i, j, k, l, occupations_at(conf,ImgTime(0.5))) * (ImgTime(0.6) - ImgTime(0.5))
                                                                                + ΔW_diag(mod, i, j, k, l, occupations_at(conf,ImgTime(0.6))) * (ImgTime(0.8) - ImgTime(0.6))
                                                                                + ΔW_diag(mod, i, j, k, l, occupations_at(conf,ImgTime(0.8))) * (ImgTime(1) + τ1 - ImgTime(0.8))
                                                                                ) * λ
end

@testset "τ_borders" begin
    @test τ_borders(conf, Set([a]), ImgTime(0.0)) == (ImgTime(0.8), ImgTime(0.2))
    @test τ_borders(conf, Set([a,e]), ImgTime(0.5)) == (ImgTime(0.2), ImgTime(0.6))
    @test τ_borders(conf, Set([b,d]), ImgTime(0.1)) == (ImgTime(0.8), ImgTime(0.2))
    @test τ_borders(conf, Set([e]), ImgTime(0.9)) == (ImgTime(0.0), ImgTime(1.0))

    conf1 = Configuration(sphere(PlaneWave((0,0,0),Up),dk=1))
    m = UEG()
    ens = CEnsemble(2,2,7)
    for _ in (1:5)
        dv, Δ = add_type_B!(m, ens, conf1)
        if dv == 0
            apply_back_step!(conf1, Δ)
        end
        dv, Δ = add_type_C!(m, ens, conf1)
        if dv == 0
            apply_back_step!(conf1, Δ)
        end
        dv, Δ = add_type_D!(m, ens, conf1)
        if dv == 0
            apply_back_step!(conf1, Δ)
        end
        dv, Δ = add_type_E!(m, ens, conf1)
        if dv == 0
            apply_back_step!(conf1, Δ)
        end
    end

    for _ in 1:10
        kink = rand(conf1.kinks)

        @test first(τ_borders(conf1, orbs(last(kink)), first(kink))) == τ_prev_affecting(conf1.kinks, orbs(last(kink)), first(kink))
        @test last(τ_borders(conf1, orbs(last(kink)), first(kink))) == τ_next_affecting(conf1.kinks, orbs(last(kink)), first(kink))
    end
end

@testset "add_orbs" begin
    c2 = Set{PlaneWave{3}}()
    @test add_orbs(c2,(a,)) == Set([a])
    @test add_orbs(c2,(a,b,c,)) == Set([a,b,c])
    @test add_orbs(c2,nothing) == c2
end

@testset "add_orbs!" begin
    c2 = Set{PlaneWave{3}}()
    add_orbs!(c2,(a,))
    @test c2 == Set([a])
    add_orbs!(c2,(b,c))
    @test c2 == Set([a,b,c])
    add_orbs!(c2,nothing)
    @test c2 == Set([a,b,c])
end

@testset "drop_orbs" begin
    @test drop_orbs(Set([a,b,c]), (a,b)) == Set([c])
    @test drop_orbs(Set([a,b,c]), (c,)) == Set([a,b])
    @test drop_orbs(Set([a,b,c]), nothing) == Set([a,b,c])
end

@testset "drop_orbs!" begin
    c2 = Set([a,b,c])
    drop_orbs!(c2, nothing)
    @test c2 == Set([a,b,c])
    drop_orbs!(c2, (a,))
    @test c2 == Set([b,c])
    drop_orbs!(c2, (b,c))
    @test c2 == Set{PlaneWave{3}}([])
end

@testset "add_kinks" begin
    @test add_kinks(Kinks( ImgTime(0.1) => T2(a,b), ImgTime(0.2) => T2(b,a), ImgTime(0.3) => T2(a,c) )
,(ImgTime(0.25) => T2(b,d),)) == Kinks( ImgTime(0.1) => T2(a,b), ImgTime(0.2) => T2(b,a), ImgTime(0.25) => T2(b,d), ImgTime(0.3) => T2(a,c) )
    @test add_kinks(Kinks( ImgTime(0.9) => T2(a,b), ImgTime(0.95) => T2(b,a) ), (ImgTime(0.1) => T2(a,b), ImgTime(0.2) => T2(b,a), ImgTime(0.3) => T2(a,c))) == Kinks( ImgTime(0.1) => T2(a,b), ImgTime(0.2) => T2(b,a), ImgTime(0.3) => T2(a,c), ImgTime(0.9) => T2(a,b), ImgTime(0.95) => T2(b,a) )
    @test add_kinks(sd, nothing) == sd
end

@testset "add_kinks!" begin
    c2 = Kinks( ImgTime(0.1) => T2(a,b), ImgTime(0.2) => T2(b,a), ImgTime(0.3) => T2(a,c) )
    add_kinks!(c2,(ImgTime(0.4) => T2(d,e),))
    @test c2 == Kinks( ImgTime(0.1) => T2(a,b), ImgTime(0.2) => T2(b,a), ImgTime(0.3) => T2(a,c), ImgTime(0.4) => T2(d,e))
    add_kinks!(c2,(ImgTime(0.9) => T2(a,b), ImgTime(0.95) => T2(b,a)))
    @test c2 == Kinks( ImgTime(0.1) => T2(a,b), ImgTime(0.2) => T2(b,a), ImgTime(0.3) => T2(a,c), ImgTime(0.4) => T2(d,e), ImgTime(0.9) => T2(a,b), ImgTime(0.95) => T2(b,a))
    add_kinks!(c2,nothing)
    @test c2 == Kinks( ImgTime(0.1) => T2(a,b), ImgTime(0.2) => T2(b,a), ImgTime(0.3) => T2(a,c), ImgTime(0.4) => T2(d,e), ImgTime(0.9) => T2(a,b), ImgTime(0.95) => T2(b,a))

end

@testset "drop_kinks!" begin
    c2 = Kinks( ImgTime(0.1) => T2(a,b), ImgTime(0.2) => T2(b,a), ImgTime(0.3) => T2(a,c), ImgTime(0.4) => T2(d,e), ImgTime(0.9) => T2(a,b), ImgTime(0.95) => T2(b,a))
    drop_kinks!(c2, (ImgTime(0.1) => T2(a,b),))
    @test c2 == Kinks( ImgTime(0.2) => T2(b,a), ImgTime(0.3) => T2(a,c), ImgTime(0.4) => T2(d,e), ImgTime(0.9) => T2(a,b), ImgTime(0.95) => T2(b,a))
    drop_kinks!(c2, (ImgTime(0.4) => T2(d,e), ImgTime(0.9) => T2(a,b)))
    @test c2 == Kinks( ImgTime(0.2) => T2(b,a), ImgTime(0.3) => T2(a,c), ImgTime(0.95) => T2(b,a))
    drop_kinks!(c2, nothing)
    @test c2 == Kinks( ImgTime(0.2) => T2(b,a), ImgTime(0.3) => T2(a,c), ImgTime(0.95) => T2(b,a))
end

@testset "drop_kinks" begin
    c2 = Kinks( ImgTime(0.1) => T2(a,b), ImgTime(0.2) => T2(b,a), ImgTime(0.3) => T2(a,c), ImgTime(0.4) => T2(d,e), ImgTime(0.9) => T2(a,b), ImgTime(0.95) => T2(b,a))
    @test drop_kinks(c2, (ImgTime(0.1) => T2(a,b),)) == Kinks( ImgTime(0.2) => T2(b,a), ImgTime(0.3) => T2(a,c), ImgTime(0.4) => T2(d,e), ImgTime(0.9) => T2(a,b), ImgTime(0.95) => T2(b,a))
    @test drop_kinks(c2, (ImgTime(0.1) => T2(a,b), ImgTime(0.4) => T2(d,e), ImgTime(0.9) => T2(a,b))) == Kinks( ImgTime(0.2) => T2(b,a), ImgTime(0.3) => T2(a,c), ImgTime(0.95) => T2(b,a))
    @test drop_kinks(sd, nothing) == sd
end


text:
# ---
# title: 1315. Sum of Nodes with Even-Valued Grandparent
# id: problem1315
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/sum-of-nodes-with-even-valued-grandparent/description/>
# hidden: true
# ---
#
# Given a binary tree, return the sum of values of nodes with even-valued
# grandparent.  (A _grandparent_ of a node is the parent of its parent, if it
# exists.)
#
# If there are no nodes with an even-valued grandparent, return `0`.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2019/07/24/1473_ex1.png)**
#
#
#
#     Input: root = [6,7,8,2,7,1,3,9,null,1,4,null,null,null,5]
#     Output: 18
#     Explanation: The red nodes are the nodes with even-value grandparent while the blue nodes are the even-value grandparents.
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is between `1` and `10^4`.
#   * The value of nodes is between `1` and `100`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
push!(LOAD_PATH, "./solutions")

using Test

using ArithmeticKoans
@testset "Arithmetic Operators" begin
    @test ArithmeticKoans.sum(1, 2) == 3
    @test ArithmeticKoans.difference(8, 1) == 7
    @test ArithmeticKoans.multiplication(10, 2) == 20
    @test ArithmeticKoans.division(35, 5) == 7.0
    @test ArithmeticKoans.integer_division(35, 5) == 7
    @test ArithmeticKoans.power(2, 2) == 4
    @test ArithmeticKoans.modulo(12, 10) == 2
end

using LogicalKoans
@testset "Logical Operators" begin
    @test LogicalKoans.not(true) == false
    @test LogicalKoans.not(false) == true
    @test LogicalKoans.equals(true, true) == true
    @test LogicalKoans.equals(false, false) == true
    @test LogicalKoans.not_equals(true, false) == true
    @test LogicalKoans.less_than(1, 10) == true
    @test LogicalKoans.greater_than(10, 1) == true
    @test LogicalKoans.less_than_or_equal_to(2, 2) == true
    @test LogicalKoans.less_than_or_equal_to(1, 10) == true
    @test LogicalKoans.less_than_or_equal_to(10, 1) == false
    @test LogicalKoans.greater_than_or_equal_to(2, 2) == true
    @test LogicalKoans.greater_than_or_equal_to(10, 1) == true
    @test LogicalKoans.greater_than_or_equal_to(1, 10) == false
end

using BitwiseKoans
@testset "Bitwise Operators" begin
    @test BitwiseKoans.not(2) == -3
    @test BitwiseKoans.and(3, 5) == 1
    @test BitwiseKoans.or(2, 4) == 6
    @test BitwiseKoans.xor_koan(2, 4) == 6
    @test BitwiseKoans.right_shift(2, 1) == 1
    @test BitwiseKoans.left_shift(2, 1) == 4
end

# using ComplexNumberKoans
# @testset "Complex Numbers" begin
# end

# using StringKoans
# @testset "Strings" begin
# end

# using ConversionKoans
# @testset "Type Conversions" begin
# end

# using ArrayKoans
# @testset "Arrays" begin
# end

# using VectorOperationKoans
# @testset "Vector Operators" begin
# end

# using MultiDimensionalArrayKoans
# @testset "Multi-dimensional Arrays" begin
# end

# using IndexingKoans
# @testset "Indexing" begin
# end

# using IterationKoans
# @testset "Iteration" begin
# end

# using ModuleKoans
# @testset "Modules" begin
# end

# using StructKoans
# @testset "Structs" begin
# end

# using UnionKoans
# @testset "Unions" begin
# end

# using FunctionKoans
# @testset "Functions" begin
# end

using HighOrderFunctionKoans
@testset "High Order Functions" begin
    f = (x, y) -> x + y
    @test HighOrderFunctionKoans.round_list([1.5, 2.3, 4.7]) == [2.0, 2.0, 5.0]
    @test HighOrderFunctionKoans.apply_fn(f, 7) == 12
    @test HighOrderFunctionKoans.add_n([1, 2, 4, 5],  3) == [4, 5, 7, 8]
    @test HighOrderFunctionKoans.remove_multiples([10, 12, 15, 18, 20], 3) == [10, 20]
end

# using MethodKoans
# @testset "Methods" begin
# end

using DestructuringKoans
@testset "Destructuring" begin
    @test DestructuringKoans.add_separately(1+2im, 3-7im) == (4, -5)
    @test DestructuringKoans.add_coordinates((2, 3), (4, -1)) == (6, 2)
end

using MultipleDispatchKoans
@testset "Multiple Dispatch" begin
    dog = MultipleDispatchKoans.Dog("goodboi")
    cat = MultipleDispatchKoans.Cat("Sir Catcelot")
    @test MultipleDispatchKoans.hello() == "hello stranger"
    @test MultipleDispatchKoans.hello("julia") == "hello julia"
    @test MultipleDispatchKoans.pet_talk(dog) == "$(dog.name): bark"
    @test MultipleDispatchKoans.pet_talk(cat) == "$(cat.name): meow"
end

# using ValueTypeKoans
# @testset "Value Types" begin
# end

using PythonInterop
@testset "Python Interop" begin
    @test PythonInterop.py_sinpi(0) == 0
    @test round(PythonInterop.py_sinpi(1/pi), Base.Rounding.RoundToZero, digits=2) == 0.84

    pyfib = PythonInterop.get_pyfib()
    @test PythonInterop.jlfib(1, PythonInterop.jlfib) == 1
    @test PythonInterop.jlfib(5, PythonInterop.jlfib) == 5
    @test PythonInterop.jlfib(5, pyfib) == 5
    @test PythonInterop.jlfib(7, pyfib) == 13
end

max_stars_repo_name: lambdaclass/julialings

max_issues_repo_name: lambdaclass/julialings

max_forks_repo_name: lambdaclass/julialings


text:
# ---------------------------------------------------------------
struct STST
    buffs::Vector{Vector{Float64}}
    nbuffs::Int
    bufflen::Int
    absmax::Vector
    state::Vector
    function STST(nbuffs, bufflen)
        @assert nbuffs > 0
        @assert bufflen > 0
        buffs = [zeros(bufflen) for _ in 1:nbuffs]
        absmax = zeros(bufflen)
        curridx = 1
        firstfill = 0 # false
        state = [curridx, firstfill]
        new(buffs, nbuffs, bufflen, absmax, state)
    end
end

function Base.push!(stst::STST, vs...)
    curridx, firstfill = stst.state

    for i in eachindex(stst.buffs)
        stst.buffs[i][curridx] = vs[i]
        stst.absmax[i] = max(abs(vs[i]), stst.absmax[i])
    end

    curridx = (curridx >= stst.bufflen) ?
        (firstfill = 1) : (curridx + 1)
    stst.state .= [curridx, firstfill]
end

function is_steady(stst, ststth)
    _ , firstfill = stst.state
    (firstfill == 0) && return false

    for i in eachindex(stst.buffs)
        max_ = maximum(stst.buffs[i])
        min_ = minimum(stst.buffs[i])
        norm_ = stst.absmax[i]
        norm_ = iszero(norm_) ? 1.0 : norm_
        var_ = (max_ - min_) / norm_
        (var_ > ststth) && return false
    end
    return true
end


text:
# hello world
println("Hello World")

hello_world = "hello world"
println(hello_world)
println(typeof(hello_world))

# sum
a = 10 + 7
println(a)

# diff
a = 10 - 7
println(a)

# product
a = 10 * 7
println(a)

# quotient
a = 10/7
println(a)

# quotient
a = 7\10
println(a)

# power
a = 10^7
println(a)

# modulus
a = 10%7
println(a)