text:
module slab
using MPI

export *

# Spectral transformation of three dimensional data aligned
# such that the last component is parallelized across processes

for (name, plan, fft12) in zip((:r2cSlab, :c2cSlab),
                               (:(FFTW.rFFTWPlan{T}), :(FFTW.cFFTWPlan)),
                               (:plan_rfft, :plan_fft))
@eval begin
    immutable ($name){T<:Real}
        # Global shape
        N::Array{Int, 1}
        # Global size of domain
        L::Array{T, 1}
        # Plans
        plan12::($plan)
        plan3::FFTW.cFFTWPlan
        # Work arrays for transformations
        vT::Array{Complex{T}, 3}
        vT_view::Array{Complex{T}, 4}
        v::Array{Complex{T}, 3}
        v_view::Array{Complex{T}, 4}
        v_recv::Array{Complex{T}, 3}
        v_recv_view::Array{Complex{T}, 4}
        # Communicator
        comm::MPI.Comm
        # Amount of data to be send by MPI
        chunk::Int
        num_processes::Int

        # Constructor
        function ($name)(N, L, comm)
            # Verify input
            Nh = N[1]÷2+1
            p = MPI.Comm_size(comm)
            Np = N÷p

            # Allocate work arrays
            vT, v = Array{Complex{T}}(Nh, N[2], Np[3]), Array{Complex{T}}(Nh, Np[2], N[3])
            vT_view, v_view = reshape(vT, (Nh, Np[2], p, Np[3])), reshape(v, (Nh, Np[2], Np[3], p))
            # For MPI.Alltoall! preallocate the receiving buffer
            v_recv = similar(v); v_recv_view = reshape(v_recv, (Nh, Np[2], Np[3], p))

            # Plan Fourier transformations
            A = zeros(T, (N[1], N[2], Np[3]))
            plan12 = ($fft12)(A, (1, 2))
            plan3 = plan_fft!(v, (3, ))
            # Compute the inverse plans
            inv(plan12); inv(plan3)

            chunk = Nh*Np[2]*Np[3]
            # Now we are ready
            new(N, L, plan12, plan3,
                vT, vT_view, v, v_view, v_recv, v_recv_view,
                comm, chunk, p)
        end
    end

    # Constructor
    ($name){T<:Real}(N::Array{Int, 1}, L::Array{T, 1}, comm::Any) = ($name){T}(N, L, comm)

    # Transform real to complex as complex = T o real
    function rfft3{T<:Real}(F::($name){T}, fu::AbstractArray{Complex{T}, 3}, u::AbstractArray{T})
        A_mul_B!(F.vT, F.plan12, u)
        permutedims!(F.v_view, F.vT_view, [1, 2, 4, 3])
        MPI.Alltoall!(F.v_recv_view, F.v_view, F.chunk, F.comm)
        F.plan3*F.v_recv; fu[:] = F.v_recv
    end

    # Transform complex to real as real = T o complex
    function irfft3{T<:Real}(F::($name){T}, u::AbstractArray{T}, fu::AbstractArray{Complex{T}, 3})
        F.plan3.pinv*fu; F.v[:] = fu
        MPI.Alltoall!(F.v_recv_view, F.v_view, F.chunk, F.comm)
        permutedims!(F.vT_view, F.v_recv_view, [1, 2, 4, 3])
        A_mul_B!(u, F.plan12.pinv, F.vT)
    end

    function real_shape{T<:Real}(F::($name){T})
        return (F.N[1], F.N[2], F.N[3]÷F.num_processes)
    end

    function complex_shape{T<:Real}(F::($name){T})
        return (F.N[1]÷2+1, F.N[2]÷F.num_processes, F.N[3])
    end

    end # begin
end # for
end # module


text:
module TestSurfaceFluxes

using Random
using ClimateMachine
ClimateMachine.init()
const ArrayType = ClimateMachine.array_type()

using Test
using NonlinearSolvers
using ClimateMachine.SurfaceFluxes
const SF = SurfaceFluxes
using StaticArrays
import ClimateMachine.MPIStateArrays: array_device
import KernelAbstractions: CPU

using CLIMAParameters: AbstractEarthParameterSet
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

@testset "SurfaceFluxes - FMS Profiles" begin
    FT = Float32

    ## Discretisation altitude z
    z = ArrayType(FT[
        29.432779269303,
        30.0497139076724,
        31.6880000418153,
        34.1873479240475,
    ])
    ## Virtual Dry Static Energy at height z
    pt = ArrayType(FT[
        268.559120403867,
        269.799228886728,
        277.443023238556,
        295.79192777341,
    ])
    ## Surface Pottemp
    pt0 = ArrayType(FT[
        273.42369841804,
        272.551410044203,
        278.638168565727,
        298.133068766049,
    ])
    ## Roughness lengths
    z0 = ArrayType(FT[
        5.86144925739178e-05,
        0.0001,
        0.000641655193293549,
        3.23383768877187e-05,
    ])
    zt = ArrayType(FT[
        3.69403636275411e-05,
        0.0001,
        1.01735489109205e-05,
        7.63933834969505e-05,
    ])
    zq = ArrayType(FT[
        5.72575636226887e-05,
        0.0001,
        5.72575636226887e-05,
        5.72575636226887e-05,
    ])
    ## Speed
    speed = ArrayType(FT[
        2.9693638452068,
        2.43308757772094,
        5.69418282305367,
        9.5608693754561,
    ])
    ## Scale velocity and moisture
    u_star = ArrayType(FT[
        0.109462510724615,
        0.0932942802513508,
        0.223232887323184,
        0.290918439028557,
    ])
    q_star = ArrayType(FT[
        0.000110861442197537,
        9.44983279664197e-05,
        4.17643828631936e-05,
        0.000133135421415819,
    ])
    # No explicit buoyancy terms in ClimateMachine
    b_star = ArrayType([
        0.00690834676781433,
        0.00428178089592372,
        0.00121229800895103,
        0.00262353784027441,
    ])

    for ii in 1:length(u_star)
        # Data at first interior node (x_ave)
        qt_ave = FT(0)
        z_ave = Tuple(z)[ii]
        vdse_ave = Tuple(pt)[ii]
        u_ave = Tuple(speed)[ii]
        x_ave = ArrayType(FT[u_ave, vdse_ave, qt_ave])

        ## Initial guesses for MO parameters
        LMO_init = eps(FT)
        u_star_init = FT(0.1)
        th_star_init = -FT(0.1)
        qt_star_init = -FT(1e-5)
        MO_param_guess =
            ArrayType(FT[LMO_init, u_star_init, th_star_init, qt_star_init])

        # Surface values for variables
        u_sfc = FT(0)
        thv_sfc = Tuple(pt0)[ii]
        qt_sfc = Tuple(q_star)[ii]
        z_sfc = FT(0)
        x_s = ArrayType(FT[u_sfc, thv_sfc, qt_sfc])

        # Roughness
        z_rough = ArrayType(FT[Tuple(z0)[ii], Tuple(zt)[ii], Tuple(zq)[ii]])

        # Constants
        Δz = Tuple(z)[ii]

        args = (
            param_set,
            MO_param_guess,
            x_ave,
            x_s,
            z_rough,
            vdse_ave,
            z_ave / 2,
        )

        ## Assuming surface fluxes are not given
        result = surface_conditions(args..., SF.DGScheme())

        x_star = Array(result.x_star)
        u_star_arr = Array(u_star)
        @test (abs((x_star[1] - u_star_arr[ii]) / u_star_arr[ii])) <= FT(0.15)

        # Keeping FV-based scheme in case needed
        result = surface_conditions(args..., SF.FVScheme())

    end
end

# Shuffing disabled in GPU, Tuple for scalar indexing
cpu_shuffle(array::ArrayType) =
    array_device(array) isa CPU ? shuffle(array) : Tuple(array)

@testset "SurfaceFluxes - Recovery" begin
    FT = Float32

    # Stochastic sampling of parameter space
    for j in 1:20
        # Define MO parameters to be recovered
        u_star = cpu_shuffle(ArrayType(FT[0.05, 0.1, 0.2, 0.3, 0.4]))
        θ_star = cpu_shuffle(ArrayType(FT[0, 50, 100, 200, 300])) # Functions not robust to θ_star < 0

        # Define temperature parameters
        θ_scale = cpu_shuffle(ArrayType(FT[290, 280, 270, 260, 250]))
        θ_s = cpu_shuffle(ArrayType(FT[290, 280, 270, 260, 250]))

        # Define a set of heights
        z = cpu_shuffle(ArrayType(FT[5, 10, 15, 20, 50]))
        z_rough = cpu_shuffle(ArrayType(FT[0.001, 0.01, 0.1, 0.5, 1])) # Must be smaller than z

        # Relative initialization of MO parameters
        LMO_init = cpu_shuffle(ArrayType(FT[-1.0, 0.2, 0.5, 1.0, 2.0]))
        u_star_init = cpu_shuffle(ArrayType(FT[0.2, 0.5, 1.0, 1.5, 2.0])) # Must be positive
        θ_star_init = cpu_shuffle(ArrayType(FT[-1e-6, 1e-6, -1e-6, 1e-6, 1e-6]))

        for ii in 1:length(u_star)
            x_star_given = Tuple(ArrayType(FT[u_star[ii], θ_star[ii]]))
            L = monin_obukhov_length(
                param_set,
                u_star[ii],
                θ_scale[ii],
                -u_star[ii] * θ_star[ii],
            )
            x_s = ArrayType(FT[FT(0), θ_s[ii]])

            u_in = recover_profile(
                param_set,
                z[ii],
                u_star[ii],
                Tuple(x_s)[1],
                z_rough[ii],
                L,
                SF.MomentumTransport(),
                SF.DGScheme(),
            )
            θ_in = recover_profile(
                param_set,
                z[ii],
                θ_star[ii],
                Tuple(x_s)[2],
                z_rough[ii],
                L,
                SF.HeatTransport(),
                SF.DGScheme(),
            )

            x_in = ArrayType(FT[u_in, θ_in])

            MO_param_guess = ArrayType(FT[
                LMO_init[ii] * L,
                u_star_init[ii] * u_star[ii],
                θ_star_init[ii] * θ_star[ii],
            ])
            recovered = surface_conditions(
                param_set,
                MO_param_guess,
                x_in,
                x_s,
                z_rough[ii],
                θ_scale[ii],
                z[ii],
                SF.DGScheme(),
            )

            # Test is recovery under 5% error for all variables
            x_star_recovered = Tuple(recovered.x_star)
            @test abs(x_star_recovered[1] - x_star_given[1]) /
                  (abs(x_star_given[1]) + FT(1e-3)) < FT(0.05)
            @test abs(x_star_recovered[2] - x_star_given[2]) /
                  (abs(x_star_given[2]) + FT(1e-3)) < FT(0.05)
        end
    end
end

include("test_universal_functions.jl")

end # module TestSurfaceFluxes

ext: jl


text:
include("structs.jl")

const R = 8.314
MOLAR_MASS_IGAS = 0.028
MOLAR_MASS_PENTANE = 0.07215
tait_C₅H₁₂ = TaitEoS(35e6, 1e5, 0.2105, 616.18)
ideal_gas = IdealGasEoS(298.0, MOLAR_MASS_IGAS)

max_stars_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_issues_repo_name: sofiabelen/Ideal-Gas-Compressible-Fluid

max_forks_repo_name: sofiabelen/Two-Phase-Filtration


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

###############
#
#  Pre-Cluster
#
#  Clustering will run within GenX automatically if using the TimeDomainReduction parameter
#  in the settings file and if there are no pre-existing clustered inputs in the
#  directory, but you can use this to cluster independently beforehand or to
#  investigate the clustering process and results.
#
#  Run with the following:
#    include("PreCluster.jl")
#
#  Beforehand, set the following parameters in this file:
#    REQUIRED
#    1) inpath -- The directory where the input data is stored
#    OPTIONAL
#    2) v -- True or False, verbosity
#
#  Jack Morris (3/5/2021)
#
###############

working_path = pwd()
src_path = cd(pwd, "..") # pwd() parent
genx_path = cd(pwd, "..\\..") # pwd() grandparent

push!(LOAD_PATH, working_path)
push!(LOAD_PATH, src_path)
push!(LOAD_PATH, genx_path)

using GenX
using DataFrames
using YAML

# Input Data Directory
#inpath = "$working_path/Example_Clustering_SE"
#inpath = "$genx_path/input_data/Inputs/RealSystemExample/ISONE_Trizone_FullTimeseries"
inpath = "$genx_path\\Example_Systems\\SmallNewEngland\\ThreeZones"

# Verbosity and Plot Choices
v = false

# Settings
#settings_path = joinpath(genx_path, "GenX_settings.yml")
settings_path = joinpath(inpath, "Settings")
mysetup = YAML.load(open(joinpath(settings_path, "genx_settings.yml")))
TDRpath = joinpath(inpath, mysetup["TimeDomainReductionFolder"])

if (isfile(TDRpath*"/Load_data.csv")) || (isfile(TDRpath*"/Generators_variability.csv")) || (isfile(TDRpath*"/Fuels_data.csv"))
    println("Data is already clustered. Delete '*_clustered.csv' files and retry: ")
    println(inpath)
elseif mysetup["TimeDomainReduction"] == 0
    println("Trying to pre-cluster the inputs, but the TimeDomainReduction setting is set to 0. Set to 1 and try again.")
    println(settings_path)
else
    Rep_Period, Weights, RMSE, TDRsetup, col_to_zone_map = cluster_inputs(inpath, settings_path, mysetup, v)
    if v
        for res in RMSE
            println(res)
        end
    end
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
# makemask.jl

using ImageGeoms
#include("../src/makemask.jl")

using Test: @test, @testset, @test_throws, @inferred

@testset "makemask" begin
    @test MaskAll <: Mask
    @test MaskCircle <: Mask
    @test MaskEllipse <: Mask
    @test MaskAllButEdge <: Mask

    i2 = ImageGeom( dims=(3,4) )
    i3 = ImageGeom( dims=(3,4,5) )

    m = @inferred makemask(i2, MaskAll())
    @test m == trues(i2)

    m = @inferred makemask(i3, MaskAll())
    @test m == trues(i3)

    m = @inferred makemask(i2, MaskAllButEdge())
    @test m == Bool[0 0 0 0; 0 1 1 0; 0 0 0 0]

    m = @inferred makemask(i3, MaskAllButEdge())
    @test sum(m) == 2 * size(i3,3)

    i2 = ImageGeom( dims=(8,6) )
    i3 = ImageGeom( dims=(8,6,4) )

    m2 = @inferred makemask(i2, MaskCircle())
    @test sum(m2) == 12

    m3 = @inferred makemask(i3, MaskCircle())
    @test sum(m3) == 12 * size(i3,3)
    @test m3[:,:,end] == m2

    m2 = @inferred makemask(i2, MaskEllipse())
    @test sum(m2) == 20

    m3 = @inferred makemask(i3, MaskEllipse())
    @test sum(m3) == 20 * size(i3,3)

    i1 = ImageGeom( dims=(8,) )
    @test_throws MethodError makemask(i1, MaskAllButEdge())

    i2 = @inferred ImageGeom( MaskCircle() ; dims=(8,6) )
    @test i2.mask == makemask(i2, MaskCircle())
end


text:
@inline function start_worker_processes(num_workers::Integer;
                                        cookie::String,
                                        delay,
                                        julia_worker_timeout::Integer)
    output_directories = Vector{String}(undef, num_workers)
    processes = Vector{Base.Process}(undef, num_workers)
    output_strings = Vector{String}(undef, num_workers)
    for i = 1:num_workers
        output_directories[i] = mktempdir(; cleanup = true)
    end
    for i = 1:num_workers
        output_directory = output_directories[i]
        stdout_file = joinpath(output_directory, "stdout.txt")
        stderr_file = joinpath(output_directory, "stderr.txt")
        p = withenv("JULIA_WORKER_TIMEOUT" => "$(julia_worker_timeout)") do
            cmd = `julia --worker=$(cookie)`
            my_pipeline = pipeline(cmd, stdout = stdout_file, stderr = stderr_file)
            run(my_pipeline; wait = false)
        end
        processes[i] = p
    end
    sleep(delay)
    for i = 1:num_workers
        output_directory = output_directories[i]
        stdout_file = joinpath(output_directory, "stdout.txt")
        output = read(stdout_file, String)::String
        output_strings[i] = output
    end
    return processes, output_strings
end


text:
@testset "Function: nrbline" begin
  line = nrbline();
  @test line.coefs == [0.0 1.0; 0.0 0.0; 0.0 0.0; 1.0 1.0];
  @test line.knots == vec([0.0 0.0 1.0 1.0]);
  @test line.number == 2;
  @test line.order == 2;
  @test typeof(line)<:NURBS1D;

  line = nrbline([-1.0;-1.0;-3.0],[3.0;1.0;1.0]);
  @test line.coefs == [-1.0 3.0;-1.0 1.0;-3.0 1.0; 1.0 1.0];
  @test line.knots == vec([0.0 0.0 1.0 1.0]);
  @test line.number == 2;
  @test line.order == 2;
  @test typeof(line)<:NURBS1D;

  @test_throws ArgumentError nrbline([-1.0;-1.0;-3.0],[3.0;1.0;1.0;2.0]);
  @test_throws ArgumentError nrbline([-1.0;-1.0;-3.0],[3.0;1.0]);
end # testset


text:
@testset "CoNLL-U" begin

    trees = collect(test_treebank("english.conllu"))

    @test length(trees) == 2
    @test length.(trees) == [7, 19]

    for C in [ArcStandard(), ArcEager(), ListBasedNonProjective()], tree in trees
        tokens = form.(tree)
        oracle = StaticOracle(C)
    end

    for C in [ArcEager(), ArcHybrid()]
        oracle = StaticOracle(C, arc = untyped)
        for (cfg, t) in DependencyTrees.xys(oracle, trees)
            @test DT.args(t) == ()
        end
    end

    @test DependencyTrees.noval(CoNLLU).head == -1

    # make sure the errors get thrown correctly
    @test_throws MultiWordTokenError CoNLLU("18-19	cannot	_	_	_	_	_	_	_	SpaceAfter=No")
    @test_throws EmptyTokenError CoNLLU("0.1	nothing	_	_	_	_	_	_	_	_")

    c = CoNLLU("1	Distribution	distribution	NOUN	S	Number=Sing	7	nsubj	_	_")
    @test c.feats == ["Number=Sing"]

    @test_throws Exception CoNLLU("1	2	3")

    sent = """
1	They	they	PRON	PRP	Case=Nom|Number=Plur	2	nsubj	2:nsubj|4:nsubj	_
2	buy	buy	VERB	VBP	Number=Plur|Person=3|Tense=Pres	0	root	0:root	_
3	and	and	CONJ	CC	_	4	cc	4:cc	_
4	sell	sell	VERB	VBP	Number=Plur|Person=3|Tense=Pres	2	conj	0:root|2:conj	_
5	books	book	NOUN	NNS	Number=Plur	2	obj	2:obj|4:obj	_
6	.	.	PUNCT	.	_	2	punct	2:punct	_"""

    graph = DependencyTree{CoNLLU}(sent)
    for d in graph.tokens
        @test untyped(d) == ()
        @test typed(d) == (d.deprel,)
    end
    @test startswith(showstr(graph), "DependencyTree{CoNLLU}\n1\tThey")
    for (i, tok) in enumerate(tokens(graph))
        @test startswith(showstr(tok), string(tok.id))
    end

    c = CoNLLU("1	They	they	PRON	PRP	Case=Nom|Number=Plur	2	nsubj	2:nsubj|4:nsubj	_")
    @test length(c.deps) == 2

    sent = """
1	Sue	Sue	_	_	_	2	_	_	_
2	likes	like	_	_	_	0	_	_	_
3	coffee	coffee	_	_	_	2	_	_	_
4	and	and	_	_	_	2	_	_	_
5	Bill	Bill	_	_	_	2	_	_	_
5.1	likes	like	_	_	_	_	_	_	_
6	tea	tea	_	_	_	2	_	_	_
"""

    graph = DependencyTree{CoNLLU}(sent)
    @test length(graph.emptytokens) == 1
    @test length(graph) == 6

    sent = """
1-2	vámonos	_	_	_	_	_	_	_	_
1	vamos	ir	_	_	_	0	_	_	_
2	nos	nosotros	_	_	_	1	_	_	_
3-4	al	_	_	_	_	_	_	_	_
3	a	a	_	_	_	5	_	_	_
4	el	el	_	_	_	5	_	_	_
5	mar	mar	_	_	_	1	_	_	_
"""

    graph = DependencyTree{CoNLLU}(sent)
    @test length(graph.emptytokens) == 0
    @test length(graph.mwts) == 2
    @test length(graph) == 5
end


text:
using Pkg; Pkg.activate(@__DIR__)
using StructArrays
import JuKeBOX as Jk
using GLMakie
using Makie

using JuKeBOX
using ForwardDiff
using ComradeBase

x = range(-15.0, 15.0, length=128)
y = range(-15.0, 15.0, length=128)

ix,iy = Tuple(ind)

g, acc = bam(2, spin, 1.0, 6.0, 5.0, 0.9, π/2, π/2)
o = Observer(1.0, inc)

s = Jk.SimpleModel(acc, g, o)

function f(spin, inc, x, y)
    g, acc = bam(1, spin, 1.0, 6.0, 5.0, 0.9, π/2, π/2)
    o = Jk.Observer(1.0, inc)
    first(raytrace(x, y, g, o, acc))
end

function img(spin, inc)
    p = [spin, inc]
    x = range(-15.0, 15.0, length=128)
    y = range(-15.0, 15.0, length=128)
    println(x)
    broadcast(x,y') do xx,yy
        g(x) = f(x[1], x[2], xx, yy)
        ForwardDiff.gradient(g, p)
    end
end

pimg = f.(spin, inc, x, y')
heatmap(pimg)

spin = rand()
inc = rand()*π/2

dimg = img(spin, inc)
incderiv = last.(dimg)
ind = sum(isnan.(incderiv))
heatmap(incderiv)
aderiv = first.(dimg)
heatmap(aderiv)

fov = 30.0
alpha = range(-15.0, 15.0, length=124)
beta  = range(-15.0, 15.0, length=124)
polim = StructArray((I=zeros(length(alpha), length(beta)), Q = zeros(length(alpha), length(beta)), U = zeros(length(alpha), length(beta))))

inc = Observable(0.01)
spin = Observable(0.5)
g, acc = Jk.bam(2, 0.25, 1.0, 6.0, 5.0, 0.9, π/2, π/2)
l1 = on(inc) do val
    obs = Jk.Observer(1.0, val)
    g = Jk.Kerr(spin.val)
    t = @elapsed Jk.traceimg!(polim, alpha, beta, g, obs, acc)
    println("I took $(t*1000.0) ms to run")
    io[] = polim.I'
end

l2 = on(spin) do val
    obs = Jk.Observer(1.0, inc.val)
    g = Jk.Kerr(val)
    t = @elapsed Jk.traceimg!(polim, alpha, beta, g, obs, acc)
    println("I took $(t*1000.0) ms to run")
    println("Using spin $val and inc $(inc.val)")
    io[] = polim.I'
end

io = Observable(polim.U')

fig, ax = image(alpha, beta, io, colormap=:afmhot,
                axis=(aspect=1, xlabel="RA (μas)", ylabel="DEC (μas)"),
                #title="Inclination = $inc"
                )
Makie.deactivate_interaction!(ax, :rectanglezoom)
spoint = select_point(ax.scene)

on(spoint) do z
    x,y = z
    inc[] =  (fov/2+y)/fov*π/2
    spin[] = x/(fov/2)
end

function test(a, θ)
    g, acc = bam(1, a, 1.5, 5.0, 1.0, 0.9, π/2, π/3)
    o = Observer(1.0, θ)
    #s = Jk.SimpleModel(acc, g, o)
    #sI = ComradeBase.SingleStokes(s, :I)
    v = raytrace(5.0, 2.0, g, o, acc)
    return norm(StokesVector(v[1],v[2],v[3],v[4]))
end

function foo(u)
    x = u[1]
    y = u[2]
    v = ComradeBase.StokesVector((x^2, y^2, (x-y)^2, x^2+y^2))
    return norm(v)
end


text:
export blockCG

function blockCG{T1,T2}(A::SparseMatrixCSC{T1,Int},B::Array{T2,2}; kwargs...)
	X = zeros(promote_type(T1,T2),size(A,2),size(B,2)) # pre-allocate
	return blockCG(V -> At_mul_B!(1.0,A,V,0.0,X),B;kwargs...) # multiply with transpose of A for efficiency
end

blockCG(A,B::Array;kwargs...) = blockCG(X -> A*X,B;kwargs...)

"""
blockCG(A,B,X)

Preconditioned Conjugate Gradient Method for solving

A * X = B

Input:

	A           - matrix or function computing A*x
	B           - array of right hand sides
	tol         - error tolerance
	maxIter     - maximum number of iterations
	M           - preconditioner, a function that computes M\\x
	X           - array of starting guesses (will be overwritten)
	out         - flag for output (-1: no output, 0: only errors, 1: final status, 2: residual norm at each iteration)
	ortho       - flag for re-orthogonalization (default: false)
	pinvTol     - tolerance for pseudoinverse (default: eps(T)*size(B,1))
	storeInterm - flag for storing iterates (default: false)

Output:

  X       - approximate solutions (2D array if !storeInterm, 3D array else)
  flag    - exit flag ( 0 : desired tolerance achieved,
  					 -1 : maxIter reached without converging
  					 -9 : right hand side was zero)
  err     - norm of relative residual, i.e., norm(A*X-B)/norm(B)
  iter    - number of iterations
  resvec  - norm of relative residual at each iteration

Reference:

	O'Leary, D. P. (1980). The block conjugate gradient algorithm and related methods.
	Linear Algebra and Its Applications, 29, 293–322. http://doi.org/10.1016/0024-3795(80)90247-5

"""
function blockCG{T<:AbstractFloat}(A::Function,B::Array{T};X=zeros(T,size(B)),M::Function=identity,maxIter=20,tol=1e-2,ortho::Bool=false,pinvTol =eps(T)*size(B,1),out::Int=0,storeInterm::Bool=false)

if norm(B)==0; return zeros(eltype(B),size(B)),-9,0.0,0,[0.0]; end

R = copy(B)
Z = M(R)
if ortho
    P, = mgs(Z);
else
	P = copy(Z); #P must be different memory than Z and R.
end

n, nrhs = size(X)
nB      = computeNorm(B)
resmat  = zeros(maxIter,nrhs)

if storeInterm
    Xout = zeros(n,nrhs,maxIter)	# allocate space for intermediates
end

if out==2
	println("=== blockcg ===")
	println(@sprintf("%4s\t%7s","iter","max(relres)"))
end

# pre-allocate
PTQ = zeros(T,nrhs,nrhs)
PTR = zeros(T,nrhs,nrhs)
QTZ = zeros(T,nrhs,nrhs)

flag = -1;
iter = 0
for iter=1:maxIter
	Q = A(P)
    # PTQ   = P'*Q
	BLAS.gemm!('T','N',1.0,P,Q,0.0,PTQ)

    # Alpha = (PTQ)\(P'*R);
    BLAS.gemm!('T','N',1.0,P,R,0.0,PTR)
	pinvPTQ = getPinv!(PTQ,pinvTol)
    Alpha = pinvPTQ*PTR

    # X     += P*Alpha
    BLAS.gemm!( 'N','N',1.0, P, Alpha, 1.0, X)
	if storeInterm; Xout[:,:,iter] = X; end
    # R     -= Q*Alpha
    BLAS.gemm!('N','N',-1.0,Q,Alpha,1.0,R)

    resmat[iter,:] = computeNorm(R)./nB
    if out==2
		println(@sprintf("%3d\t%1.2e",iter,maximum(resmat[iter,:])))
	end
    if maximum(resmat[iter,:]) < tol
	    flag = 0;
        break;
    end

    Z     = M(R)
    #Beta  = -(PTQ)\(Q'*Z);
	BLAS.gemm!('T','N',1.0,Q,Z,0.0,QTZ)
    Beta  = -pinvPTQ*QTZ

	# Z might be just R here - don't overwite it! Q is not needed.
	Q[:] = Z;
	BLAS.gemm!('N','N',1.0,P,Beta,1.0,Q);
	P[:] = Q;

    if ortho
        P,     =  mgs!(P)
 	end
end
if out>=0
	if flag==-1
		println(@sprintf("blockCG iterated maxIter (=%d) times but reached only residual norm %1.2e instead of tol=%1.2e.",																							maxIter,maximum(resmat[iter,:]),tol))
	elseif flag==0 && out>=1
		println(@sprintf("blockCG achieved desired tolerance at iteration %d. Residual norm is %1.2e.",iter,maximum(resmat[iter,:])))
	end
end
if storeInterm
    return Xout[:,:,1:iter],flag,resmat[iter,:],iter,resmat[1:iter,:]
else
    return X,flag,resmat[iter,:],iter,resmat[1:iter,:]
end
end

function computeNorm(R)
	n,nrhs = size(R)
	res    = zeros(nrhs)
	for k=1:nrhs
		for i=1:n
			res[k]+=R[i,k]*R[i,k]
		end
	end
	return sqrt(res)
end

function getPinv!(A,pinvTol)
	SVD         = svdfact!(A)
	Sinv        = zeros(length(SVD.S))
    index       = SVD.S .> pinvTol*maximum(SVD.S)
    Sinv[index] = 1.0./ SVD.S[index]
    Sinv[find(!isfinite(Sinv))] = 0.0
    return SVD.Vt'scale(Sinv, SVD.U')
end


text:
"""
    RecurrentQNetworkArchitecture
    specify an architecture with the following form:
    [CONV]-[FC]-[LSTM]-[FC]
"""
@with_kw mutable struct RecurrentQNetworkArchitecture
   fc_in::Vector{Int64} = Vector{Int64}[]
   convs::Vector{Tuple{Int64, Vector{Int64}, Int64}} = Vector{Tuple{Int64, Tuple{Int64, Int64}, Int64}}[]
   fc_out::Vector{Int64} = Vector{Int64}[]
   lstm_size::Int64 = 32
end

"""
    DeepRecurrentQLearningSolver
Deep Q learning with a recurrent module to solve POMDPs

## Fields:
- `arch::RecurrentQNetworkArchitecture` Specify the architecture of the Q network default = QNetworkArchitecture(conv=[], fc=[])
- `lr::Float64` learning rate default = 0.005
- `max_steps::Int64` total number of training step default = 1000
- `target_update_freq::Int64` frequency at which the target network is updated default = 500
- `batch_size::Int64` batch size sampled from the replay buffer default = 32
- `trace_length::Int64` trajectory length used to train the LSTM
- `train_freq::Int64` frequency at which the active network is updated default  = 4
- `log_freq::Int64` frequency at which to logg info default = 100
- `eval_freq::Int64` frequency at which to eval the network default = 100
- `num_ep_eval::Int64` number of episodes to evaluate the policy default = 100
- `eps_fraction::Float64` fraction of the training set used to explore default = 0.5
- `eps_end::Float64` value of epsilon at the end of the exploration phase default = 0.01
- `double_q::Bool` double q learning udpate default = true
- `dueling::Bool` dueling structure for the q network default = true
- `prioritized_replay::Bool` enable prioritized experience replay default = true
- `prioritized_replay_alpha::Float64` default = 0.6
- `prioritized_replay_epsilon::Float64` default = 1e-6
- `prioritized_replay_beta::Float64` default = 0.4
- `buffer_size::Int64` size of the experience replay buffer default = 1000
- `max_episode_length::Int64` maximum length of a training episode default = 100
- `train_start::Int64` number of steps used to fill in the replay buffer initially default = 200
- `grad_clip::Bool` enables gradient clipping default = true
- `clip_val::Float64` maximum value for the grad norm default = 10.0
- `rng::AbstractRNG` random number generator default = MersenneTwister(0)
- `verbose::Bool` default = true
"""
@with_kw mutable struct DeepRecurrentQLearningSolver
    arch::RecurrentQNetworkArchitecture = RecurrentQNetworkArchitecture()
    lr::Float64 = 0.001
    max_steps::Int64 = 1000
    target_update_freq::Int64 = 500
    batch_size::Int64 = 32
    trace_length = 6
    train_freq::Int64  = 4
    log_freq::Int64 = 100
    eval_freq::Int64 = 100
    num_ep_eval::Int64 = 100
    eps_fraction::Float64 = 0.5
    eps_end::Float64 = 0.01
    double_q::Bool = true
    dueling::Bool = true
    buffer_size::Int64 = 10000
    max_episode_length::Int64 = 100
    train_start::Int64 = 200
    grad_clip::Bool = true
    clip_val::Float64 = 10.0
    rng::AbstractRNG = MersenneTwister(0)
    logdir::String = "log"
    save_freq::Int64 = 10000
    verbose::Bool = true
end

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
include("./src/MiniPascal.jl")
using .MiniPascal
const m = MiniPascal

if length(ARGS) != 1
  println("USAGE: julia compile.jl sourcefile")
  exit(1)
end

filename = ARGS[1]

if !isfile(filename)
  println("File $(filename) not found.")
  exit(1)
end

source = open(ARGS[1]) do file
  read(file, String)
end

source = try
  ascii(source)
catch e
  println("The source code contains non-ascii characters.")
  exit(1)
end

try
  io = IOBuffer()
  m.generate(source, io)
  code = String(take!(io))
  filename = "out.ll"
  file = open(filename, "w")
  println(file, code)
  close(file)
catch e
  showerror(stdout, e)
end

println()
exit(0)

max_stars_repo_name: juhamyllari/minipascal-julia

max_issues_repo_name: juhamyllari/minipascal-julia

max_forks_repo_name: juhamyllari/minipascal-julia


text:
using InfiniteLinearAlgebra, FillArrays, LazyArrays, ArrayLayouts, Test

@testset "Compact" begin
    A = ApplyMatrix(hvcat, 2, randn(5,5), Zeros(5,∞), Zeros(∞,5), Zeros(∞,∞))
    b = [randn(10); zeros(∞)];
    @test ((I + A) \ b)[1:10] ≈ (I+A)[1:10,1:10] \ b[1:10]

    C = zeros(∞,∞);
    C[1:5,1:5] .= randn.()
    @test_skip ((I + C) \ b)[1:10] ≈ (I+C)[1:10,1:10] \ b[1:10]
end


text:
using Documenter
using VectorDataUtils

DocMeta.setdocmeta!(VectorDataUtils, :DocTestSetup, :(using VectorDataUtils); recursive=true)
doctest(VectorDataUtils)

makedocs(
    sitename = "VectorDataUtils",
    format = Documenter.HTML(),
    modules = [VectorDataUtils],
    pages = [
        "Home" => "index.md",
        "Contents" => [
            "norms.md",
            "distances.md",
            "running_stats.md",
            "kmeans.md",
            "utilities.md"
        ],
    ],
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/vidhyasaharan/VectorDataUtils.jl.git"
)


text:
using Test
using Comonicon
using Comonicon.Options

@testset "load options" begin
    options = read_options(pkgdir(Comonicon, "test"))
    display(options)
    @test options == Options.Comonicon(;
        name = "foo",
        install = Options.Install(; compile = "min"),
        sysimg = Options.SysImg(;
            cpu_target = "native",
            precompile = Options.Precompile(; execution_file = ["deps/precopmile.jl"]),
        ),
        download = Options.Download(; user = "Roger-luo", repo = "Foo.jl"),
        application = Options.Application(;
            assets = Options.Asset[asset"PkgTemplate: templates", asset"assets/images"],
            incremental = true,
            filter_stdlibs = false,
        ),
    )
end

@testset "has_comonicon_toml" begin
    @test Options.has_comonicon_toml(Main) == false
end


text:
# module escape_dimer_one_hit
# export escape_exit_num_dimer
using DifferentialEquations
using Roots
using PolynomialRoots

using Printf
# using MATLAB
# using ProgressMeter

include("leap_frog_definitions.jl")
max_hit_q1_forward_one = 1
H=.25

condition_max_hits_forward_one(u, t, integrator) = u[3] == max_hit_q1_forward_one
affect_stop!(integrator) = terminate!(integrator)
function condition_hits_PSS_q1(u, t, integrator) # Event when event_f(u,t) == 0
    u[1]
end
function affect_update_iterator_q1!(integrator)
    q2=integrator.u[2]
    p2=integrator.u[5]
    p1=integrator.u[4]
    P1=P1_poly(q2,p2, H)
    if !isempty(P1) && abs(p1-P1)<1e-5
        integrator.u[3] = integrator.u[3] + 1
    end
end

callback_max_hits_forward_one = DiscreteCallback(condition_max_hits_forward_one, affect_stop!)
callback_hits_PSS_q1_forward_one = ContinuousCallback(condition_hits_PSS_q1, nothing,  affect_update_iterator_q1!, rootfind = true)
cb_forward_one = CallbackSet(callback_hits_PSS_q1_forward_one, callback_max_hits_forward_one)
function forward_one(Q2, P2, t_end, H)

    P1 = P1_poly(Q2, P2, H)
    if ~isempty(P1)
        q0, p0 = [zeros(3) for i in 1:2]
        q0[1] = 0
        q0[2] = Q2
        q0[3] = 0
        p0[1] = P1
        p0[2] = P2
        p0[3] = 0
       #constructor for ODE
        prob = HamiltonianProblem{true}(Hamiltonian_Dimer, q0, p0, (0., t_end));
        #solve ode , save_everystep=false is important to prevent sol to include all points, not just event points
        # sol = solve(prob,  Vern9(), maxiters = 1e20, reltol = 1e-13, abstol = 1e-15, callback = cb_forward_one, save_start = false, save_end = true, save_everystep = false)
        # sol=solve(prob,  RK4(), maxiters = 1e20, callback=cb_forward_one,save_start=false,save_end=true,save_everystep=false)
        # sol=solve(prob, RK4(),maxiters=1e20,reltol=1e-8,abstol=1e-10,callback=cb,save_start=true,save_end=true,save_everystep=false)
        sol = solve(prob, Vern9(), maxiters = 1e20, reltol = 1e-13, abstol = 1e-15, callback = cb_forward_one, save_start = true, save_end = false, save_everystep = false)
        # uf = zeros(6)
        # uf[1] = sol[1,3]
        # uf[2] = sol[2,3]
        # uf[3] = sol[3,3]
        # uf[4] = sol[4,3]
        # uf[5] = sol[5,3]
        # uf[6] = sol[6,3]
        # q1 = uf[1]
        # q2 = uf[2]
        # p1 = uf[4]
        # p2 = uf[5]
        if  sol[3,end] == 1
            uf = zeros(6)
            uf[1] = sol[1,3]
            uf[2] = sol[2,3]
            uf[3] = sol[3,3]
            uf[4] = sol[4,3]
            uf[5] = sol[5,3]
            uf[6] = sol[6,3]
            q1 = uf[1]
            q2 = uf[2]
            p1 = uf[4]
            p2 = uf[5]
            q = [q1, q2]
            p = [p1, p2]
            diff = abs(Hamiltonian_Dimer(q, p, 0) - H)
            if diff < 1e-1
                return q1, q2, p1, p2
            end
        end
        # if dH<1e-5
        #     if is_it(uf,H,tol_dist)
        #         return sol.u[end][5]
        #     else
        #         return -1
        #     end
        # end
    end
end
function unzipper2(array_of_points)
    N = length(array_of_points)
    Q = zeros(N)
    P = zeros(N)
    for (index, point) in enumerate(array_of_points)
        Q[index] = point[2]
        P[index] = point[4]
    end
    return Q, P
end
function unzipper3(array_of_points)
    N = length(array_of_points)
    Q0 = zeros(N)
    P0 = zeros(N)
    Q1 = zeros(N)
    P1 = zeros(N)
    for (index, point) in enumerate(array_of_points)
        Q0[index] = point[2]
        P0[index] = point[4]
        Q1[index] = point[6]
        P1[index] = point[8]
    end
    return Q0, P0, Q1, P1
end

function unzipper3(array_of_points)
    N = length(array_of_points)
    Q0 = zeros(N)
    P0 = zeros(N)
    Q1 = zeros(N)
    P1 = zeros(N)
    for (index, point) in enumerate(array_of_points)
        Q0[index] = point[2]
        P0[index] = point[4]
        Q1[index] = point[6]
        P1[index] = point[8]
    end
    return Q0, P0, Q1, P1
end
function unzipper_quad(array_of_points)
    N = length(array_of_points)
    Q = zeros(N)
    P = zeros(N)
    for (index, point) in enumerate(array_of_points)
        Q[index] = point[2]
        P[index] = point[4]
    end
    full_Q = vcat(Q, Q, -Q, -Q)
    full_P = vcat(P, -P, P, -P)
    return full_Q, full_P
end
function unzipper_half(array_of_points)
    N = length(array_of_points)
    Q = zeros(N)
    P = zeros(N)
    for (index, point) in enumerate(array_of_points)
        Q[index] = point[2]
        P[index] = point[4]
    end
    full_Q = vcat(Q, -Q)
    full_P = vcat(P, -P)
    return full_Q, full_P
end


text:
# Example of vOptSpecific, for the 2OSP
# July 2017

using vOptSpecific

n =4
p = [2, 4, 3, 1]
d = [1, 2, 4, 6]

id = set2OSP(n, p, d)
z1, z2 , S = vSolve(id) #Default solver is OSP_Van_Wassenhove1980()

for i = 1:length(S)
    println("(", z1[i], ", ", z2[i], ") : ", S[i])
end

max_stars_repo_name: vOptSolver/vOptSpecific

max_issues_repo_name: vOptSolver/vOptSpecific

max_forks_repo_name: vOptSolver/vOptSpecific


text:
f(x, y, z) = exp.(x .* y .* z) .* sin.(π .* x) .* sin.(π .* y) .* sin.(π .* z)

function advection_dominated(;N = 50, β = 1000.0)
    # Problem: Δu + βuₓ = f
    # u = 0 on the boundaries
    # f(x, y, z) = exp(xyz) sin(πx) sin(πy) sin(πz)
    # 2nd order central differences (shows serious wiggles)

    # Total number of unknowns
    n = N^3

    # Mesh width
    h = 1.0 / (N + 1)

    # Interior points only
    xs = range(0, stop=1, length=N + 2)[2 : N + 1]

    # The Laplacian
    Δ = laplace_matrix(Float64, N, 3) ./ -h^2

    # And the dx bit.
    ∂x_1d = spdiagm(-1 => fill(-β / 2h, N - 1), 1 => fill(β / 2h, N - 1))
    ∂x = kron(speye(N^2), ∂x_1d)

    # Final matrix and rhs.
    A = Δ + ∂x
    b = reshape([f(x, y, z) for x ∈ xs, y ∈ xs, z ∈ xs], n)

    A, b
end

function laplace_matrix(::Type{T}, n, dims) where T
    D = second_order_central_diff(T, n)
    A = copy(D)

    for idx = 2 : dims
        A = kron(A, speye(n)) + kron(speye(size(A, 1)), D)
    end

    A
end

second_order_central_diff(::Type{T}, dim) where {T} = convert(
    SparseMatrixCSC{T, Int},
    SymTridiagonal(fill(2 * one(T), dim), fill(-one(T), dim - 1))
)


text:
using Oxygen

println("Testing Atom")

@testset "OxygenAtom" begin
    carbon = OxygenAtom(6, "C", 12.011, Oxygen.Electrons(6), 2.55)
    @test typeof(carbon.number) == Int
    @test typeof(carbon.symbol) == String
    @test typeof(carbon.mass) == Float64
    @test typeof(carbon.orbitals) == Oxygen.Electrons
    @test typeof(carbon.electronegativity) == Float64
    @test carbon == OxygenAtom(6)
    @test carbon == OxygenAtom("C")
end

@testset "Printing" begin
    carbon = OxygenAtom(6, "C", 12.011, Oxygen.Electrons(6), 2.55)
    @test repr(carbon) == "\"C\" atom (number 6)"
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct Subscription <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    subscriptionId::Any # spec type: Union{ Nothing, String } # spec name: subscriptionId
    displayName::Any # spec type: Union{ Nothing, String } # spec name: displayName
    tenantId::Any # spec type: Union{ Nothing, String } # spec name: tenantId
    state::Any # spec type: Union{ Nothing, String } # spec name: state
    subscriptionPolicies::Any # spec type: Union{ Nothing, SubscriptionPolicies } # spec name: subscriptionPolicies
    authorizationSource::Any # spec type: Union{ Nothing, String } # spec name: authorizationSource
    managedByTenants::Any # spec type: Union{ Nothing, Vector{ManagedByTenant} } # spec name: managedByTenants
    tags::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: tags

    function Subscription(;id=nothing, subscriptionId=nothing, displayName=nothing, tenantId=nothing, state=nothing, subscriptionPolicies=nothing, authorizationSource=nothing, managedByTenants=nothing, tags=nothing)
        o = new()
        validate_property(Subscription, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(Subscription, Symbol("subscriptionId"), subscriptionId)
        setfield!(o, Symbol("subscriptionId"), subscriptionId)
        validate_property(Subscription, Symbol("displayName"), displayName)
        setfield!(o, Symbol("displayName"), displayName)
        validate_property(Subscription, Symbol("tenantId"), tenantId)
        setfield!(o, Symbol("tenantId"), tenantId)
        validate_property(Subscription, Symbol("state"), state)
        setfield!(o, Symbol("state"), state)
        validate_property(Subscription, Symbol("subscriptionPolicies"), subscriptionPolicies)
        setfield!(o, Symbol("subscriptionPolicies"), subscriptionPolicies)
        validate_property(Subscription, Symbol("authorizationSource"), authorizationSource)
        setfield!(o, Symbol("authorizationSource"), authorizationSource)
        validate_property(Subscription, Symbol("managedByTenants"), managedByTenants)
        setfield!(o, Symbol("managedByTenants"), managedByTenants)
        validate_property(Subscription, Symbol("tags"), tags)
        setfield!(o, Symbol("tags"), tags)
        o
    end
end # type Subscription

const _property_map_Subscription = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"), Symbol("subscriptionId")=>Symbol("subscriptionId"), Symbol("displayName")=>Symbol("displayName"), Symbol("tenantId")=>Symbol("tenantId"), Symbol("state")=>Symbol("state"), Symbol("subscriptionPolicies")=>Symbol("subscriptionPolicies"), Symbol("authorizationSource")=>Symbol("authorizationSource"), Symbol("managedByTenants")=>Symbol("managedByTenants"), Symbol("tags")=>Symbol("tags"))
const _property_types_Subscription = Dict{Symbol,String}(Symbol("id")=>"String", Symbol("subscriptionId")=>"String", Symbol("displayName")=>"String", Symbol("tenantId")=>"String", Symbol("state")=>"String", Symbol("subscriptionPolicies")=>"SubscriptionPolicies", Symbol("authorizationSource")=>"String", Symbol("managedByTenants")=>"Vector{ManagedByTenant}", Symbol("tags")=>"Dict{String, String}")
Base.propertynames(::Type{ Subscription }) = collect(keys(_property_map_Subscription))
Swagger.property_type(::Type{ Subscription }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_Subscription[name]))}
Swagger.field_name(::Type{ Subscription }, property_name::Symbol) =  _property_map_Subscription[property_name]

const _allowed_Subscription_state = ["Enabled", "Warned", "PastDue", "Disabled", "Deleted"]

function check_required(o::Subscription)
    true
end

function validate_property(::Type{ Subscription }, name::Symbol, val)
    if name === Symbol("state")
        Swagger.validate_param(name, "Subscription", :enum, val, _allowed_Subscription_state)
    end
end


text:
using Test
using Random
using LinearAlgebra
using Uchiyama

include("test_periodic_hard_spheres.jl")

@testset "step" begin

   n = 2
   ϵ = 0.1
   q = [[0.5,0.2], [0.5,0.8]]
   v = [[0,1], [0,-1]]

   rng = MersenneTwister(42)
   p = SquareParticles(rng, n, ϵ)
   for i in 1:n
       p.q[i] = q[i]
       p.v[i] = v[i]
   end

   collisions = PeriodicCollisions(p)

   @test step!(p, collisions) ≈ 0.2

   @test p.q[1] ≈ [0.5,0.4]
   @test p.v[1] ≈ [0,-1]
   @test p.q[2] ≈ [0.5,0.6]
   @test p.v[2] ≈ [0, 1]

   @test step!(p, collisions) ≈ 0.3

   @test p.q[1] ≈ [0.5,0.1]
   @test p.v[1] ≈ [0,1]
   @test p.q[2] ≈ [0.5,0.9]
   @test p.v[2] ≈ [0,-1]

end

@testset "overlap" begin

   import Uchiyama: overlap

   q  = [0.5, 0.5]

   p  = [0.6, 0.6]
   ϵ  = 0.1
   @test overlap( p, q, ϵ )
   p  = [0.65, 0.65]
   @test overlap( p, q, ϵ )
   p  = [0.65, 0.5]
   @test overlap( p, q, ϵ )
   p  = [0.5, 0.65]
   @test overlap( p, q, ϵ )
   p  = [0.8, 0.8]
   @test !overlap( p, q, ϵ )
   p  = [0.1, 0.1]
   @test !overlap( p, q, ϵ )
   p  = [0.69, 0.69]
   @test overlap( p, q, ϵ )
   p  = [0.41, 0.4]
   @test overlap( p, q, ϵ )

end

@testset "init_particles" begin

     n = 40
     ϵ = 0.01
     rng = MersenneTwister(1234)

     sq = SquareParticles(rng, n, ϵ)
     @test true
     hs = HardSpheres(rng, n, ϵ)
     @test true

end

include("test_tempscoll.jl")
include("test_fantome.jl")


text:
"""
    find_biclique_cover(skeleton, faces)

Given a skeleton of a planar graph and its faces, finds a biclique cover for the complement
of the respective finite element graph with a divide-and-conquer approach: construct a biclique
from a separator for the current finite element graph and call the function recursively on the
induced subskeletons and subfaces of A ∪ C and B ∪ C.
"""
function find_biclique_cover(skeleton::LabeledGraph{T}, faces::Set{Vector{T}})::Set{Pair{Set{T}, Set{T}}} where {T}
    face_pairs = _find_face_pairs(faces)
    feg = ClutteredEnvPathOpt._find_finite_element_graph(skeleton, face_pairs)

    if ClutteredEnvPathOpt.lg_is_complete(feg)
        return Set{Pair{Set{T}, Set{T}}}()
    end

    (C, A, B) = _find_feg_separator_lt_no_empty(skeleton, face_pairs)
    # println("Set A: $A\nSet B: $B\nSet C: $C")

    if !isempty(A) && !isempty(B)
        node = Set([A => B])
    else
        node = Set{Pair{Set{T}, Set{T}}}()
        println("A or B are empty.\nSkeleton vertices are $(skeleton.labels)")
    end

    if !isempty(union(A,C))
        skeleton_ac, faces_ac = _find_skeleton_faces(union(A, C), skeleton, faces)
        if !isempty(faces_ac)
            left = find_biclique_cover(skeleton_ac, faces_ac)
        else
            left = Set{Pair{Set{T}, Set{T}}}()
        end
    else
        left = Set{Pair{Set{T}, Set{T}}}()
    end

    if !isempty(union(B,C))
        skeleton_bc, faces_bc = _find_skeleton_faces(union(B, C), skeleton, faces)
        if !isempty(faces_bc)
            right = find_biclique_cover(skeleton_bc, faces_bc)
        else
            right = Set{Pair{Set{T}, Set{T}}}()
        end
    else
        right = Set{Pair{Set{T}, Set{T}}}()
    end

    return union(node, left, right)
end

"""
    find_biclique_cover_debug(skeleton, faces, points, obstacles, file_name)

Same as find_biclique_cover, with the added features of plotting the new skeletons and
recording separator information at each call.
"""
function find_biclique_cover_debug(skeleton::LabeledGraph{T}, faces::Set{Vector{T}}, points::Vector{Any}, obstacles, file_name)::Set{Pair{Set{T}, Set{T}}} where {T}
    ClutteredEnvPathOpt.plot_edges(skeleton, points, obstacles)
    f = open(file_name, "a")
    for (i,face) in enumerate(faces)
        write(f, "Face $i: $face\n")
        # println("Face $i: $face")
    end
    flush(f)
    close(f)

    face_pairs = _find_face_pairs(faces)
    feg = ClutteredEnvPathOpt._find_finite_element_graph(skeleton, face_pairs)

    if ClutteredEnvPathOpt.lg_is_complete(feg)
        f = open(file_name, "a")
        write(f, "FEG is complete.\n")
        flush(f)
        close(f)

        # println("FEG is complete.")
        return Set{Pair{Set{T}, Set{T}}}()
    end

    (C, A, B) = _find_feg_separator_lt_no_empty(skeleton, face_pairs)

    if !isempty(A) && !isempty(B)
        node = Set([A => B])
    else
        node = Set{Pair{Set{T}, Set{T}}}()
        f = open(file_name, "a")
        write(f, "A or B is empty.\n")
        flush(f)
        close(f)

        # println("A or B is empty.")
    end

    if !isempty(union(A,C))
        skeleton_ac, faces_ac = _find_skeleton_faces(union(A, C), skeleton, faces)
        if !isempty(faces_ac)  # I think this will always be nonempty of A ∪ C is nonempty
            f = open(file_name, "a")
            write(f, "\nPlot \nSet A: $A\nSet B: $B\nSet C: $C\nA ∪ C\n")
            flush(f)
            close(f)

            # println("\nSet A: $A\nSet B: $B\nSet C: $C\nA ∪ C")
            left = find_biclique_cover_debug(skeleton_ac, faces_ac, points, obstacles, file_name)
        else
            left = Set{Pair{Set{T}, Set{T}}}()
            f = open(file_name, "a")
            write(f, "faces_ac is empty.\n")
            flush(f)
            close(f)

            # println("faces_ac is empty.")
        end
    else
        left = Set{Pair{Set{T}, Set{T}}}()
        f = open(file_name, "a")
        write(f, "A ∪ C is empty.\n")
        flush(f)
        close(f)

        # println("A ∪ C is empty.")
    end

    if !isempty(union(B,C))
        skeleton_bc, faces_bc = _find_skeleton_faces(union(B, C), skeleton, faces)
        if !isempty(faces_bc)
            f = open(file_name, "a")
            write(f, "\nPlot \nSet A: $A\nSet B: $B\nSet C: $C\nB ∪ C\n")
            flush(f)
            close(f)

            # println("\nSet A: $A\nSet B: $B\nSet C: $C\nB ∪ C")
            right = find_biclique_cover_debug(skeleton_bc, faces_bc, points, obstacles, file_name)
        else
            right = Set{Pair{Set{T}, Set{T}}}()
            f = open(file_name, "a")
            write(f, "faces_bc is empty.\n")
            flush(f)
            close(f)

            # println("faces_bc is empty.")
        end
    else
        right = Set{Pair{Set{T}, Set{T}}}()
        f = open(file_name, "a")
        write(f, "B ∪ C is empty.\n")
        flush(f)
        close(f)

        # println("B ∪ C is empty.")
    end

    return union(node, left, right)
end

function find_biclique_cover_one_iter(skeleton::LabeledGraph{T}, faces::Set{Vector{T}}) where {T}
    face_pairs = _find_face_pairs(faces)
    feg = ClutteredEnvPathOpt._find_finite_element_graph(skeleton, face_pairs)

    if ClutteredEnvPathOpt.lg_is_complete(feg)
        return Set{Pair{Set{T}, Set{T}}}()
    end

    (C, A, B) = _find_feg_separator_lt_no_empty(skeleton, face_pairs)
    # println("Set A: $A\nSet B: $B\nSet C: $C")

    # if !isempty(A) && !isempty(B)
    #     node = Set([A => B])
    # else
    #     node = Set{Pair{Set{T}, Set{T}}}()
    #     println("A or B are empty.\nSkeleton vertices are $(skeleton.labels)")
    # end

    if !isempty(union(A,C))
        skeleton_ac, faces_ac = _find_skeleton_faces(union(A, C), skeleton, faces)
    #     if !isempty(faces_ac)
    #         left = find_biclique_cover(skeleton_ac, faces_ac)
    #     else
    #         left = Set{Pair{Set{T}, Set{T}}}()
    #     end
    # else
    #     left = Set{Pair{Set{T}, Set{T}}}()
    end

    if !isempty(union(B,C))
        skeleton_bc, faces_bc = _find_skeleton_faces(union(B, C), skeleton, faces)
    #     if !isempty(faces_bc)
    #         right = find_biclique_cover(skeleton_bc, faces_bc)
    #     else
    #         right = Set{Pair{Set{T}, Set{T}}}()
    #     end
    # else
    #     right = Set{Pair{Set{T}, Set{T}}}()
    end

    # return union(node, left, right)
    return (C,A,B), skeleton_ac, faces_ac, skeleton_bc, faces_bc
end

"""
    biclique_merger(cover, feg)

Given a biclique cover, merges bicliques together when possible for a smaller cover.
"""
function biclique_merger(cover::Set{Pair{Set{T}, Set{T}}}, feg::LabeledGraph{T})::Set{Pair{Set{T},Set{T}}} where {T}
    # Sort bicliques by size so that we attempt merging smaller bicliques first
    cover = sort(collect(cover), by=(s->length(s.first)+length(s.second)))

    # Loop until we are unable to find any merges
    merge_found = true
    while merge_found
        for i = 1:length(cover)
            for j = i+1:length(cover)
                # TODO: Check if this is right
                # if (!isempty(intersect(cover[i].first, cover[j].first)) && !isempty(intersect(cover[i].first, cover[j].second))) ||
                #     (!isempty(intersect(cover[i].second, cover[j].first)) && !isempty(intersect(cover[i].second, cover[j].second)))
                #     continue
                # end

                potential_merge_found = false
                merge_found = false
                A1_adj_A2 = false
                A1_adj_B2 = false
                B1_adj_A2 = false
                B1_adj_B2 = false
                for u in cover[i].first
                    for v in cover[j].first
                        if LightGraphs.has_edge(feg.graph, feg.labels[u], feg.labels[v])
                            A1_adj_A2 = true
                            break
                        end
                    end
                    if A1_adj_A2
                        break
                    end
                end

                for u in cover[i].first
                    for v in cover[j].second
                        if LightGraphs.has_edge(feg.graph, feg.labels[u], feg.labels[v])
                            A1_adj_B2 = true
                            break
                        end
                    end
                    if A1_adj_B2
                        break
                    end
                end

                if A1_adj_A2 && A1_adj_B2
                    continue  # cannot merge the bicliques in this case
                end

                if !A1_adj_A2
                    for u in cover[i].second
                        for v in cover[j].second
                            if LightGraphs.has_edge(feg.graph, feg.labels[u], feg.labels[v])
                                B1_adj_B2 = true
                                break
                            end
                        end
                        if B1_adj_B2
                            break
                        end
                    end
                end

                if !A1_adj_B2
                    for u in cover[i].second
                        for v in cover[j].first
                            if LightGraphs.has_edge(feg.graph, feg.labels[u], feg.labels[v])
                                B1_adj_A2 = true
                                break
                            end
                        end
                        if B1_adj_A2
                            break
                        end
                    end
                end

                # Two possible merges. If both are possible, will need to see if it matters which merge we choose
                A = Set([])
                B = Set([])
                if !A1_adj_A2 && !B1_adj_B2
                    # A = A1 ∪ B2
                    A = union(cover[i].first, cover[j].second)
                    # B = B1 ∪ A2
                    B = union(cover[i].second, cover[j].first)

                    potential_merge_found = true
                elseif !A1_adj_B2 && !B1_adj_A2
                    # A = A1 ∪ A2
                    A = union(cover[i].first, cover[j].first)
                    # B = B1 ∪ B2
                    B = union(cover[i].second, cover[j].second)

                    potential_merge_found = true
                end

                # Second condition is to prevent vertices from being in both A and B
                if potential_merge_found && isempty(intersect(A,B)) # don't need second condition if above is correct
                    merge_found = true
                    deleteat!(cover, j)
                    deleteat!(cover, i)
                    push!(cover, A => B)
                    sort!(cover, by=(s->length(s.first)+length(s.second)))
                    break
                end

            end
            if merge_found
                break
            end
        end

    end

    return Set(cover)
end

struct BicliqueCoverTree
    biclique
    left
    right
end

function find_biclique_cover_as_tree(skeleton::LabeledGraph{T}, faces::Set{Vector{T}}) where {T}
    face_pairs = ClutteredEnvPathOpt._find_face_pairs(faces)
    feg = ClutteredEnvPathOpt._find_finite_element_graph(skeleton, face_pairs)

    if ClutteredEnvPathOpt.lg_is_complete(feg)
        return nothing
    end

    (C, A, B) = ClutteredEnvPathOpt._find_feg_separator_lt_no_empty(skeleton, face_pairs)

    if !isempty(A) && !isempty(B)
        node = A => B
    else
        node = nothing
    end
    if !isempty(union(A,C))
        skeleton_ac, faces_ac = ClutteredEnvPathOpt._find_skeleton_faces(union(A, C), skeleton, faces)
        if !isempty(faces_ac)
            left = ClutteredEnvPathOpt.find_biclique_cover_as_tree(skeleton_ac, faces_ac)
        else
            left = nothing
        end
    else
        left = nothing
    end
    if !isempty(union(B,C))
        skeleton_bc, faces_bc = ClutteredEnvPathOpt._find_skeleton_faces(union(B, C), skeleton, faces)
        if !isempty(faces_bc)
            right = ClutteredEnvPathOpt.find_biclique_cover_as_tree(skeleton_bc, faces_bc)
        else
            right = nothing
        end
    else
        right = nothing
    end

    return BicliqueCoverTree(node, left, right)
end

function flatten(tree)
    if (tree === nothing)
        return []
    end

    return union([tree.biclique], flatten(tree.left), flatten(tree.right))
end

function getEdges(tree, cover)
    res = ""

    if (tree === nothing)
        return res
    end

    node_num = findfirst(x -> x === tree.biclique, cover)

    if (tree.left !== nothing)
        left_num = findfirst(x -> x === tree.left.biclique, cover)
        res = string(res, node_num, " -> ", left_num, " ", getEdges(tree.left, cover))
    end

    if (tree.right !== nothing)
        right_num = findfirst(x -> x === tree.right.biclique, cover)
        res = string(res, node_num, " -> ", right_num, " ", getEdges(tree.right, cover))
    end

    return res
end

function tree2digraph(tree)
    cover = collect(flatten(tree))

    edges = getEdges(tree, cover)

    return (cover, string("Digraph G { ", edges, " }"))
end

"""
    _find_skeleton_faces(vertices, old_skeleton, old_faces)

Find the skeleton and faces of a finite element subgraphgraph given a subset of
vertices, a finite element graph's skeleton and the faces of its planar
embedding. Returns a (skeleton, face vector) tuple.

Note: ensure the faces are passed in as a set of clockwise-ordered vertices
"""
function _find_skeleton_faces(vertices::Set{T}, old_skeleton::LabeledGraph{T}, old_faces::Set{Vector{T}})::Tuple{LabeledGraph{T}, Set{Vector{T}}} where {T}
    skeleton = copy(old_skeleton)
    for vertex in keys(old_skeleton.labels)
        if !in(vertex, vertices)
            ClutteredEnvPathOpt.lg_rem_vertex!(skeleton, vertex)
        end
    end

    # Remove obsolete vertices from each face and return faces with at least 1 vertex
    faces = @pipe map(face -> filter(vertex -> in(vertex, vertices), face), collect(old_faces)) |>
                filter(face -> length(face) > 0, _)

    # If a new face is a subset of another new face, don't include it
    new_faces = []
    removed_indices = []
    for i = 1:length(faces)
        indicator = 0
        for j = 1:length(faces)
            if j == i || j in removed_indices
                continue
            end
            if intersect(Set(faces[i]), Set(faces[j])) == Set(faces[i]) && length(faces[j]) >= length(faces[i])
                indicator += 1
                break
            end
        end
        if indicator != 0
            push!(removed_indices, i)
        else
            push!(new_faces, faces[i])
        end
    end

    # Close the new subface
    for face in new_faces
        if length(face) > 1
            for i = 2:length(face)
                ClutteredEnvPathOpt.lg_add_edge!(skeleton, face[i-1], face[i])
            end
            ClutteredEnvPathOpt.lg_add_edge!(skeleton, face[end], face[1])
        end
    end

    return skeleton, Set(new_faces)
end

"""
    _find_face_pairs(faces)

Convert a set of list of vertices included in a face to a set of sets of edges
comprising a face.
"""
function _find_face_pairs(faces::Set{Vector{T}})::Set{Set{Pair{T, T}}} where {T}
    face_pairs = Set{Set{Pair{T, T}}}()
    for face in faces
        pairs = Set{Pair{T, T}}()

        for i in 2:length(face)
            push!(pairs, face[i - 1] => face[i])
        end

        push!(pairs, face[end] => face[1])

        push!(face_pairs, pairs)
    end

    return face_pairs
end

"""
    _find_feg_separator_lt_no_empty(skeleton, face_pairs)

Find the separator of a finite element graph, repeating the process if for some
given root either A or B is empty. Returns the separator for which |A|*|B| is largest.
"""
function _find_feg_separator_lt_no_empty(skeleton::LabeledGraph{T}, face_pairs::Set{Set{Pair{T, T}}})::Tuple{Set{T}, Set{T}, Set{T}} where {T}
    feg = ClutteredEnvPathOpt._find_finite_element_graph(skeleton, face_pairs)
    # for root in keys(skeleton.labels)
    #     (C, A, B) = find_feg_separator_lt(skeleton, face_pairs, root)
    #     (C, A, B) = pp_expell(feg, C, A, B)

    #     if !isempty(A) && !isempty(B)
    #         return (C, A, B)
    #     end
    # end

    # return (Set([]), Set([]), Set([]))
    #####################################################################################################

    separators = @pipe map(root -> find_feg_separator_lt(skeleton, face_pairs, root), collect(keys(skeleton.labels))) |>
        map(sep -> pp_expell(feg, sep[1], sep[2], sep[3]), _)

    balances = map(separator -> begin
        if isempty(separator[2]) || isempty(separator[3]) # A or B empty is no bueno
            return -1.0
        elseif isempty(separator[1]) # C empty, A and B are not, prioritize this case
            return 1.0 + min(length(separator[2]) / length(separator[3]), length(separator[3]) / length(separator[2]))
        else # None are empty
            return 1.0  # the separator that finds the most edges will be chosen
            # return min(length(separator[2]) / length(separator[3]), length(separator[3]) / length(separator[2]))
        end
    end, separators)

    max_balance = max(balances...)
    # println("Max Balance: $max_balance")

    if max_balance < -0.5
        # If at least one of A or B is empty, the subskeleton would be the same and we'll have infinite recursion
        # In such a case, our SEPARATOR ALGORITHM and PP_EXPELL cannot produce two nonempty A and B
        # Return empty sets since there'll be no biclique
        return (Set([]), Set([]), Set([]))
    else
        # Return the biclique that makes the most edges (e.g. |A|*|B| is largest)
        select = filter(b -> balances[b] == max_balance, 1:length(balances))
        max_edges = 0
        best_index = 0
        for i in select
            if length(separators[i][2])*length(separators[i][3]) > max_edges
                max_edges = length(separators[i][2])*length(separators[i][3])
                best_index = i
            end
        end
        return separators[best_index]
    end

end

"""
    _is_valid_biclique_cover(lg, cover)

Tests whether or not a cover is a valid biclique cover of a given graph.
"""
function _is_valid_biclique_cover(lg::LabeledGraph{T}, cover::Set{Pair{Set{T}, Set{T}}})::Bool where {T}
    e_bar = LightGraphs.edges(LightGraphs.complement(lg.graph))

    edges = @pipe map(pair -> _cartesian_product(pair.first, pair.second), collect(cover)) |> reduce(union!, _, init=Set{Pair{T, T}}())
    test = LightGraphs.SimpleGraph(LightGraphs.nv(lg.graph))

    for edge in edges
        # edge.first/second will be our nodes, so we map them to the nodes they are in the graph
        LightGraphs.add_edge!(test, lg.labels[edge.first], lg.labels[edge.second])
    end

    return LightGraphs.edges(test) == e_bar
end

"""
    _is_valid_biclique_cover_diff(lg, cover)

Tests whether or not a cover is a valid biclique cover of a given graph. Returns
their difference in cardinality (positive => missing edges, negative => surplus of edges)
and a LabeledGraph containing any missing edges.
"""
function _is_valid_biclique_cover_diff(lg::LabeledGraph{T}, cover::Set{Pair{Set{T}, Set{T}}})::Tuple{Bool,Number,Set{LightGraphs.SimpleGraphs.SimpleEdge{T}},LabeledGraph{T}} where {T}
    e_bar = LightGraphs.edges(LightGraphs.complement(lg.graph))

    edges = @pipe map(pair -> _cartesian_product(pair.first, pair.second), collect(cover)) |> reduce(union!, _, init=Set{Pair{T, T}}())
    test = LightGraphs.SimpleGraph(LightGraphs.nv(lg.graph))
    rev = _reverse_labels(lg.labels)

    for edge in edges
        # edge.first/second will be our nodes, so we map them to the nodes they are in the graph
        LightGraphs.add_edge!(test, lg.labels[edge.first], lg.labels[edge.second])
    end

    # Create a LabeledGraph containing only the missing edges
    set_diff = setdiff(Set(e_bar), Set(LightGraphs.edges(test)))
    missing_e_graph = LightGraphs.SimpleGraph(LightGraphs.nv(lg.graph))
    missing_e_lg = LabeledGraph{Int}(missing_e_graph, lg.labels)
    for edge in set_diff
        lg_add_edge!(missing_e_lg, rev[edge.src], rev[edge.dst])
    end

    return LightGraphs.edges(test) == e_bar, length(e_bar) - LightGraphs.ne(test), set_diff, missing_e_lg
end

"""
    _cartesian_product(A, B)

Returns the cartesian product of two sets
"""
function _cartesian_product(A::Set{V}, B::Set{W})::Set{Pair{V, W}} where {V, W}
    res = Set{Pair{V, W}}()

    for a in A
        for b in B
            push!(res, Pair(a, b))
        end
    end

    return res
end

"""
    _find_feg_separator_lt_no_empty_visual(skeleton, face_pairs, locations)
Finds (C, A, B) and plots the edges found between A and B.
"""
function _find_feg_separator_lt_no_empty_visual(skeleton::LabeledGraph{T}, face_pairs::Set{Set{Pair{T, T}}}, points::Vector{Any})::Tuple{Set{T}, Set{T}, Set{T}} where {T}
    (C,A,B) = _find_feg_separator_lt_no_empty(skeleton, face_pairs)

    if isempty(A) || isempty(B)
        return (C,A,B)
    end

    # Create edges between A and B
    new_edges = ClutteredEnvPathOpt._cartesian_product(A, B)

    Plots.plot()
    Plots.scatter!(map(point -> point.first, points), map(point -> point.second, points))

    for a in A
        Plots.scatter!([points[a].first], [points[a].second], color=:red)
    end
    for b in B
        Plots.scatter!([points[b].first], [points[b].second], color=:yellow)
    end

    for edge in new_edges
        (x1,y1) = points[edge.first]
        (x2,y2) = points[edge.second]
        Plots.plot!([x1, x2], [y1, y2], color=:green, linestyle=:dash)
    end

    display(Plots.plot!(legend=:false,xlims=(-.05,1.05),ylims=(-0.05,1.05), title="Biclique"))

    return (C,A,B)
end

"""
    _find_biclique_cover_visual(skeleton, faces, locations)
Finds a biclique cover and plots the edges of each biclique.
"""
function _find_biclique_cover_visual(skeleton::LabeledGraph{T}, faces::Set{Vector{T}}, points::Vector{Any})::Set{Pair{Set{T}, Set{T}}} where {T}
    face_pairs = _find_face_pairs(faces)
    feg = ClutteredEnvPathOpt._find_finite_element_graph(skeleton, face_pairs)

    if ClutteredEnvPathOpt.lg_is_complete(feg)
        return Set{Pair{Set{T}, Set{T}}}()
    end

    (C, A, B) = _find_feg_separator_lt_no_empty_visual(skeleton, face_pairs, points)
    # println("Set A: $A\nSet B: $B\nSet C: $C")

    if !isempty(A) && !isempty(B)
        node = Set([A => B])
    else
        node = Set{Pair{Set{T}, Set{T}}}()
    end
    if !isempty(union(A,C))
        skeleton_ac, faces_ac = _find_skeleton_faces(union(A, C), skeleton, faces)
        if !isempty(faces_ac)
            # println("\nSet A: $A\nSet B: $B\nSet C: $C\nA ∪ C")
            left = _find_biclique_cover_visual(skeleton_ac, faces_ac, points)
        else
            left = Set{Pair{Set{T}, Set{T}}}()
        end
    else
        left = Set{Pair{Set{T}, Set{T}}}()
    end
    if !isempty(union(B,C))
        skeleton_bc, faces_bc = _find_skeleton_faces(union(B, C), skeleton, faces)
        if !isempty(faces_bc)
            # println("\nSet A: $A\nSet B: $B\nSet C: $C\nB ∪ C")
            right = _find_biclique_cover_visual(skeleton_bc, faces_bc, points)
        else
            right = Set{Pair{Set{T}, Set{T}}}()
        end
    else
        right = Set{Pair{Set{T}, Set{T}}}()
    end

    return union(node, left, right)
end


text:
""" Loads Iris dataset for tests. """
function load_irisdata()
    path = "data/iris.data"
    csv_data = CSV.File(path, header=false)
    iris_names = ["sepal_len", "sepal_wid", "petal_len", "petal_wid", "class"]
    df = DataFrame(csv_data)
    rename!(df, iris_names)
    dropmissing!(df)
    return df
end

""" Turns classes in DataFrame (last column) into binary classes. """
function binarize(df = load_irisdata())
    df_bin = copy(df)
    X = df_bin[:,1:end-1]
    Y = df_bin[:, end]
    classes = unique(Y)
    for i = 1:length(Y)
        if Y[i] != classes[1]
            Y[i] = classes[2]
        end
    end
    df_bin[:, end] = Y
    return df_bin
end


text:
function _backend_services_create_jwt(config::BackendServicesConfig)
    jti = Random.randstring(150)

    now = TimeZones.now(TimeZones.localzone())
    expiration_time = now + Dates.Minute(4)
    expiration_time_seconds_since_epoch_utc = integer_seconds_since_the_epoch_utc(expiration_time)

    jwt_payload_claims_dict = Dict(
        "iss" => config.iss,
        "sub" => config.sub,
        "aud" => config.token_endpoint,
        "jti" => jti,
        "exp" => expiration_time_seconds_since_epoch_utc,
    )
    jwt = JSONWebTokens.encode(config.private_key, jwt_payload_claims_dict)

    return jwt
end

"""
    backend_services(config::BackendServicesConfig)
"""
function backend_services(config::BackendServicesConfig)
    jwt = _backend_services_create_jwt(config)

    body_params = Dict{String, String}()
    body_params["grant_type"] = "client_credentials"
    body_params["client_assertion_type"] = "urn:ietf:params:oauth:client-assertion-type:jwt-bearer"
    body_params["client_assertion"] = jwt

    if config.scope !== nothing
        body_params["scope"] = config.scope
    end

    _response = HTTP.request(
        "POST",
        config.token_endpoint;
        headers = Dict("Content-Type" => "application/x-www-form-urlencoded"),
        body = URIs.escapeuri(body_params),
    )

    access_token_response = JSON3.read(String(_response.body))
    access_token = access_token_response.access_token

    access_token_is_jwt, access_token_jwt_decoded = try_decode_jwt(access_token)

    backend_services_result = BackendServicesResult(;
        access_token             = access_token,
        access_token_is_jwt      = access_token_is_jwt,
        access_token_jwt_decoded = access_token_jwt_decoded,
        access_token_response    = access_token_response,
    )

    return backend_services_result
end


text:
"""
Test problem to give the chord method problems.
The line search will fail in the middle of the iteration
and demand a recompute of the derivative.
"""
function linatan(x)
    alpha = 0.01
    return (1.0 + alpha * x) * atan(x)
end


text:
using Brillouin, Test
using LinearAlgebra
using StaticArrays

@testset "WignerSeitz" begin
    # --- `wignerseitz` ---
    # hexagonal lattice (example from space group 147)
    Rs = convert(SVector{3,SVector{3,Float64}}, [[1.0,0,0], [-.5,√3/2,0], [0,0,1.25]])      # direct basis
    Gs = convert(SVector{3,SVector{3,Float64}}, [[2π,2π/√3,0], [0,4π/√3,0], [0,0,2π/1.25]]) # reciprocal basis
    cell = wignerseitz(Gs)
    @test length(cell) == 8
    @test length(vertices(cell)) == 12
    @test basis(cell) == Gs

    # test that everything works the same if we use ordinary vectors instead of SVectors
    @test wignerseitz(convert(Vector{SVector{3,Float64}}, Gs)) == cell
    @test wignerseitz(convert(Vector{Vector{Float64}}, Gs))    == cell

    # ------------------------------------------------------------------------------------ #
    # for CI and testing generally, we cannot depend on a specific sorting/ordering
    # from `wignerseitz` (and ultimately, `PySpatial.Voronoi` or `PySpatial.ConvexHull`),
    # so to be robust to this, we compare output after sorting to a canonical order
    function canonicalize!(c::Cell)
        vs = vertices(c)
        fs = faces(c)
        # sort the vertices in some canonical order (we also use `round` to be robust to
        # floating point differences across architectures)
        I = sortperm(vs, by=v->round.(v, digits=5))
        permute!(vs, I)
        # relabel the face-indices to match the new vertex sorting
        replace!.(fs, (I .=> eachindex(I))...)
        # sort the faces in some canonical order; note that each face
        # must themselves also be sorted (because they are cyclical)
        map!(circshift_minimal_first, fs, fs)
        sort!(fs)

        return c
    end

    function circshift_minimal_first(v::Vector)
        _, idx = findmin(v)
        return circshift(v, length(v)-idx+1)
    end

    # actually canonicalize `cell`
    canonicalize!(cell)

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

    # test show/display (canonicalization is crucial here)
    cell_show_reference = """
    Cell{3} (8 faces, 12 vertices):
      verts: [-0.666667, 0.333333, -0.5]
             [-0.666667, 0.333333, 0.5]
             [-0.333333, -0.333333, -0.5]
             [-0.333333, -0.333333, 0.5]
             [-0.333333, 0.666667, -0.5]
             [-0.333333, 0.666667, 0.5]
             [0.333333, -0.666667, -0.5]
             [0.333333, -0.666667, 0.5]
             [0.333333, 0.333333, -0.5]
             [0.333333, 0.333333, 0.5]
             [0.666667, -0.333333, -0.5]
             [0.666667, -0.333333, 0.5]
      faces: [1, 2, 6, 5]
             [1, 3, 4, 2]
             [1, 5, 9, 11, 7, 3]
             [2, 4, 8, 12, 10, 6]
             [3, 7, 8, 4]
             [5, 6, 10, 9]
             [7, 11, 12, 8]
             [9, 10, 12, 11]
      basis: [6.283185, 3.627599, 0.0]
             [0.0, 7.255197, 0.0]
             [0.0, 0.0, 5.026548]"""

    cell_show_referenceᶜ = """
    Cell{3} (8 faces, 12 vertices):
      verts: [-4.18879, 0.0, -2.513274]
             [-4.18879, 0.0, 2.513274]
             [-2.094395, -3.627599, -2.513274]
             [-2.094395, -3.627599, 2.513274]
             [-2.094395, 3.627599, -2.513274]
             [-2.094395, 3.627599, 2.513274]
             [2.094395, -3.627599, -2.513274]
             [2.094395, -3.627599, 2.513274]
             [2.094395, 3.627599, -2.513274]
             [2.094395, 3.627599, 2.513274]
             [4.18879, 0.0, -2.513274]
             [4.18879, 0.0, 2.513274]
      faces: [1, 2, 6, 5]
             [1, 3, 4, 2]
             [1, 5, 9, 11, 7, 3]
             [2, 4, 8, 12, 10, 6]
             [3, 7, 8, 4]
             [5, 6, 10, 9]
             [7, 11, 12, 8]
             [9, 10, 12, 11]
      basis: [6.283185, 3.627599, 0.0]
             [0.0, 7.255197, 0.0]
             [0.0, 0.0, 5.026548]""" # cartesian representation
    test_show(cell, cell_show_reference)
    test_show(Brillouin.cartesianize(cell), cell_show_referenceᶜ)

    # cartesianize/latticize(!)
    @test Brillouin.latticize!(Brillouin.cartesianize!(deepcopy(cell))) ≈ cell
    cellᶜ = Brillouin.cartesianize!(deepcopy(cell))
    @test Brillouin.cartesianize!(Brillouin.latticize(cellᶜ, basis(cell)[[3,2,1]])) ≈ cellᶜ

    # test indexing/iteration of Cell struct
    @test cell[1] ≈ [[-2/3,  1/3, -1/2],
                     [-2/3,  1/3, 1/2],
                     [-1/3,  2/3, 1/2],
                     [-1/3,  2/3, -1/2]]
    @test cell[4] ≈ [[-2/3,  1/3, 1/2],
                     [-1/3, -1/3, 1/2],
                     [ 1/3, -2/3, 1/2],
                     [ 2/3, -1/3, 1/2],
                     [ 1/3,  1/3, 1/2],
                     [-1/3,  2/3, 1/2]]
    @test collect(cell) == [cell[i] for i in eachindex(cell)] # iteration vs. indexing

    # test error on out-of-bounds indexing error
    @test_throws BoundsError cell[9]

    # test :triangles output
    cell′ = wignerseitz(Gs; merge = true)
    @test all(v -> length(v) == 3, vertices(cell′))

    # ------------------------------------------------------------------------------------ #
    # 2D test case
    Rs_2d = convert(SVector{2,SVector{2,Float64}}, [[1.0,0], [-.5,√3/2]])
    cell_2d = wignerseitz(Rs_2d)
    @test length(vertices(cell_2d)) == 6
    @test length(only(faces(cell_2d))) == 6
    @test basis(cell_2d) == Rs_2d

    # 1D test case
    Rs_1d = convert(SVector{1,SVector{1,Float64}}, [[1.5]])
    cell_1d = wignerseitz(Rs_1d)
    @test vertices(cell_1d) == [[-0.5], [0.5]]
    @test only(faces(cell_1d)) == [1,2]
    @test basis(cell_1d) == Rs_1d
    @test_throws DimensionMismatch wignerseitz([SVector(1.0), SVector(1.0)])
end

@testset "reduce_to_wignerseitz" begin
    Rs = convert(SVector{3,SVector{3,Float64}}, [[1.0,0,0], [-.5,√3/2,0], [0,0,1.25]])
    cell = wignerseitz(Rs)
    vs = vertices(cell)

    # we don't guarantee anything about points exactly on the border of the Wigner-Seitz
    # both anything else should be solid: hence, check that points slightly inside stay
    # put after reduction, while points outside definitely change
    @test all(r->reduce_to_wignerseitz(r, Rs) ≈ r, .999*vs)
    @test all(r->norm(reduce_to_wignerseitz(r, Rs) - r) > norm(r)/1.1, 1.001*vs)

    # pick a point inside the cell, then shuffle it out by a few lattice vectors; check
    # that it is shuffled back in
    r = SVector(0.12, 0.32, -0.24)
    r′ = r + SVector(3, -4, 8)
    @test reduce_to_wignerseitz(r, Rs) ≈ r
    @test reduce_to_wignerseitz(r′, Rs) ≈ r

    # test a point in the rectangular unit cell that isn't in the Wigner-Seitz cell
    ruc = SVector(-0.49, 0.49, 0.49)
    @test !(reduce_to_wignerseitz(ruc, Rs) ≈ ruc)

    # test utility accessors
    r′′ = convert(Vector{Float64}, r)
    @test reduce_to_wignerseitz(r′′, Rs) ≈ r
    @test_throws DimensionMismatch reduce_to_wignerseitz(r′′[1:2], Rs)
    @test_throws DimensionMismatch reduce_to_wignerseitz(r′′, [Rs[1][1:3], Rs[2][1:3], Rs[3][1:2]])
end


text:
#!/usr/bin/julia

# TODO allow specification of log name
# TODO auto log name by date

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    "script"
        help = "the script to be run"
        required = true
    "name"
        help = "name for the simulation"
        required = false
end

args = parse_args(ARGS, s)
name = args["name"] == nothing ? args["script"] : args["name"]
log = "simlog.md" # TODO make arg
datestring = Dates.format(Dates.now(), "u d HH:MM") # TODO make arg

tmp = tempname()
script_contents = readall(args["script"])
run(pipeline(`julia $(args["script"])`, `tee $tmp`)) # TODO make cross platform

logstring = """
# [$datestring] $name

## Input
```julia
$script_contents
```
## Output
```
$(readall(tmp))
```
"""

open(log, "a") do f
    write(f, logstring)
end


text:
"""
    struct KmeansResult{T<:AbstractFloat}
        K::Int
        centroids::Vector{Vector{T}}
        cluster::Vector{Int}
        withinss::T
        iter::Int
    end

Object resulting from kmeans algorithm that contains the number of clusters, centroids, clusters prediction,
total-variance-within-cluster and number of iterations until convergence.
"""
struct KmeansResult{T<:AbstractFloat}
    K::Int
    centroids::Vector{Vector{T}}
    cluster::Vector{Int}
    withinss::T
    iter::Int
end

function Base.print(io::IO, model::KmeansResult{T}) where {T<:AbstractFloat}
    p = ["     $(v)\n" for v in model.centroids]

    print(IOContext(io, :limit => true), "KmeansResult{$T}:
 K = $(model.K)
 centroids = [\n", p..., " ]
 cluster = ", model.cluster, "
 within-cluster sum of squares = $(model.withinss)
 iterations = $(model.iter)")
end

Base.show(io::IO, model::KmeansResult) = print(io, model)

"""
    ClusterAnalysis.euclidean(a::AbstractVector, b::AbstractVector)

Calculate euclidean distance from two vectors. √∑(aᵢ - bᵢ)².

# Arguments (positional)
- `a`: First vector.
- `b`: Second vector.

# Example
```julia
julia> using ClusterAnalysis

julia> a = rand(100); b = rand(100);

julia> ClusterAnalysis.euclidean(a, b)
3.8625780213774954
```
"""
function euclidean(a::AbstractVector{T},
                   b::AbstractVector{T}) where {T<:AbstractFloat}
    @assert length(a) == length(b)

    # euclidean(a, b) = √∑(aᵢ- bᵢ)²
    s = zero(T)
    @simd for i in eachindex(a)
        @inbounds s += (a[i] - b[i])^2
    end
    return √s
end

"""
    ClusterAnalysis.squared_error(data::AbstractMatrix)
    ClusterAnalysis.squared_error(col::AbstractVector)

Function that evaluate the kmeans, using the Sum of Squared Error (SSE).

# Arguments (positional)
- `data` or `col`: Matrix of data observations or a Vector which represents one column of data.

# Example
```julia
julia> using ClusterAnalysis

julia> a = rand(100, 4);

julia> ClusterAnalysis.squared_error(a)
34.71086095943974

julia> ClusterAnalysis.squared_error(a[:, 1])
10.06029322934825
```
"""
function squared_error(data::AbstractMatrix{T}) where {T<:AbstractFloat}
    error = zero(T)
    @simd for i in axes(data, 2)
        error += squared_error(view(data, :, i))
    end
    return error
end

function squared_error(col::AbstractVector{T}) where {T<:AbstractFloat}
    μ = mean(col)
    error = zero(T)
    @simd for i in eachindex(col)
        @inbounds error += (col[i] - μ)^2
    end
    return error
end

"""
    ClusterAnalysis.totalwithinss(data::AbstractMatrix, K::Int, cluster::Vector)

Calculate the total-variance-within-cluster using the `squared_error()` function.

# Arguments (positional)
- `data`: Matrix of data observations.
- `K`: number of clusters.
- `cluster`: Vector of cluster for each data observation.

# Example
```julia
julia> using ClusterAnalysis
julia> using CSV, DataFrames

julia> iris = CSV.read(joinpath(pwd(), "path/to/iris.csv"), DataFrame);
julia> df = iris[:, 1:end-1];
julia> model = kmeans(df, 3);

julia> ClusterAnalysis.totalwithinss(Matrix(df), model.K, model.cluster)
78.85144142614601
```
"""
function totalwithinss(data::AbstractMatrix{T}, K::Int, cluster::AbstractVector{Int}) where {T<:AbstractFloat}
    # evaluate total-variance-within-clusters
    error = zero(T)
    @simd for k in 1:K
        error += squared_error(data[cluster .== k, :])
    end
    return error
end

"""
    kmeans(table, K::Int; nstart::Int = 1, maxiter::Int = 10, init::Symbol = :kmpp)
    kmeans(data::AbstractMatrix, K::Int; nstart::Int = 1, maxiter::Int = 10, init::Symbol = :kmpp)

Classify all data observations in k clusters by minimizing the total-variance-within each cluster.

# Arguments (positional)
- `table` or `data`: table or Matrix of data observations.
- `K`: number of clusters.

## Keyword
- `nstart`: number of starts.
- `maxiter`: number of maximum iterations.
- `init`: centroids inicialization algorithm - `:kmpp` (default) or `:random`.

# Example
```julia
julia> using ClusterAnalysis
julia> using CSV, DataFrames

julia> iris = CSV.read(joinpath(pwd(), "path/to/iris.csv"), DataFrame);
julia> df = iris[:, 1:end-1];

julia> model = kmeans(df, 3)
KmeansResult{Float64}:
 K = 3
 centroids = [
     [5.932307692307693, 2.755384615384615, 4.42923076923077, 1.4384615384615382]
     [5.006, 3.4279999999999995, 1.462, 0.24599999999999997]
     [6.874285714285714, 3.088571428571429, 5.791428571428571, 2.117142857142857]
 ]
 cluster = [2, 2, 2, 2, 2, 2, 2, 2, 2, 2  …  3, 3, 1, 3, 3, 3, 1, 3, 3, 1]
 within-cluster sum of squares = 78.85144142614601
 iterations = 7
```

# Pseudo-code of the algorithm:
* Repeat `nstart` times:
    1. Initialize `K` clusters centroids using KMeans++ algorithm or random init.
    2. Estimate clusters.
    3. Repeat `maxiter` times:
        * Update centroids using the mean().
        * Reestimates the clusters.
        * Calculate the total-variance-within-cluster.
        * Evaluate the stop rule.
* Keep the best result (minimum total-variance-within-cluster) of all `nstart` executions.

For more detailed explanation of the algorithm, check the
[`Algorithm's Overview of KMeans`](https://github.com/AugustoCL/ClusterAnalysis.jl/blob/main/algo_overview/kmeans_overview.md).
"""
function kmeans(table, K::Int; kwargs...)
    Tables.istable(table) ? (data = Tables.matrix(table)) : throw(ArgumentError("The table argument passed does not implement the Tables.jl interface."))
    return kmeans(data, K; kwargs...)
end

kmeans(data::AbstractMatrix{T}, K::Int; kwargs...) where {T} = kmeans(convert(Matrix{Float64}, data), K; kwargs...)

function kmeans(data::AbstractMatrix{T}, K::Int;
                nstart::Int = 1,
                maxiter::Int = 10,
                init::Symbol = :kmpp) where {T<:AbstractFloat}

    nl = size(data, 1)

    centroids = Vector{Vector{T}}(undef, K)
    cluster = Vector{Int}(undef, nl)
    withinss = Inf
    iter = 0

    # run multiple kmeans to get the best result
    for _ in 1:nstart

        new_centroids, new_cluster, new_withinss, new_iter = _kmeans(data, K, maxiter, init)

        if new_withinss < withinss
            centroids .= new_centroids
            cluster .= new_cluster
            withinss = new_withinss
            iter = new_iter
        end
    end

    return KmeansResult(K, centroids, cluster, withinss, iter)
end

function _kmeans(data::AbstractMatrix{T}, K::Int, maxiter::Int, init::Symbol) where {T<:AbstractFloat}

    nl = size(data, 1)

    # generate random centroids
    centroids = _initialize_centroids(data, K, init)

    # first clusters estimate
    cluster = Vector{Int}(undef, nl)
    for (i, obs) in enumerate(eachrow(data))
        dist = [euclidean(obs, c) for c in centroids]
        @inbounds cluster[i] = argmin(dist)
    end

    # first evaluation of total-variance-within-cluster
    withinss = totalwithinss(data, K, cluster)

    # variables to update during the iterations
    new_centroids = copy(centroids)
    new_cluster = copy(cluster)
    iter = 1
    norms = norm.(centroids)

    # start kmeans iterations until maxiter or convergence
    for _ in 2:maxiter

        # update new_centroids using the mean
        @simd for k in 1:K             # mean.(eachcol(data[new_cluster .== k, :]))
            @inbounds new_centroids[k] = vec(mean(view(data, new_cluster .== k, :), dims = 1))
        end

        # estimate cluster to all observations
        for (i, obs) in enumerate(eachrow(data))
            dist = [euclidean(obs, c) for c in new_centroids]
            @inbounds new_cluster[i] = argmin(dist)
        end

        # update iter, withinss-variance and calculate centroid norms
        new_withinss = totalwithinss(data, K, new_cluster)
        new_norms = norm.(new_centroids)
        iter += 1

        # convergence rule
        norm(norms - new_norms) ≈ 0 && break

        # update centroid norms
        norms .= new_norms

        # update centroids, cluster and whithinss
        if new_withinss < withinss
            centroids .= new_centroids
            cluster .= new_cluster
            withinss = new_withinss
        end

    end

    return centroids, cluster, withinss, iter
end

function _initialize_centroids(data::AbstractMatrix{T}, K::Int, init::Symbol) where {T<:AbstractFloat}
    nl = size(data, 1)

    if init == :random
        indexes = rand(1:nl, K)
        centroids = Vector{T}[data[i, :] for i in indexes]
    elseif init == :kmpp
        centroids = Vector{Vector{T}}(undef, K)
        centroids[1] = data[rand(1:nl), :]

        # distance vector for each observation
        dists = Vector{T}(undef, nl)

        # get each new centroid by the furthest observation (maximum distance)
        for k in 2:K

            # for each observation get the nearest centroid by the minimum distance
            for (i, row) in enumerate(eachrow(data))
                dist_c = [euclidean(row, c) for c in @view centroids[1:(k-1)]]
                @inbounds dists[i] = minimum(dist_c)
            end

            # new centroid by the furthest observation
            @inbounds centroids[k] = data[argmax(dists), :]
        end
    else
        throw(ArgumentError("The symbol :$init is not a valid argument. Use :random or :kmpp."))
    end

    return centroids
end


text:
module Intervals

import Base:
    +, -, \,
    <, >, ==, !=, ≤,
    ∈, ∩, ∪, ⊆, ⊇, ~,
    max, min,
    isless, isequal, isdisjoint,
    iterate, isempty, length, copy, getindex, to_index,
    show

export Interval, IntervalSet
export containing, partition

# -------------------------------------------------------------------------
# atomic interval

# represents an atomic (closed, open) interval
struct Interval{T<:Real}
    lo :: T
    hi :: T
    Interval{T}(lo,hi) where T <: Real = hi <= lo ? throw(ArgumentError("invalid interval hi=$hi ≤ low=$lo")) : new(lo, hi)
end

# ---------------------------------
# constructors

Interval(lo::T, hi::T) where T <: Real = Interval{T}(lo,hi)
Interval(x::T) where T <: Real = Interval{T}(x,x)
Interval(t::Tuple) = Interval(t...)

copy(a::Interval) = Interval(a.lo, a.hi)

# ---------------------------------
# repl properties

show(io::IO, a::Interval) = print(io, "[$(a.lo),$(a.hi))")

# ---------------------------------
# operations

# collection/

length(a::Interval) = a.hi - a.lo
to_index(a::Interval{T}) where T <: Integer = a.lo:(a.hi-1)

# logical

==(a::Interval, b::Interval) = a.lo == b.lo && a.hi == b.hi
≤(a::Interval, b::Interval)  = a.lo ≤ b.lo && a.hi ≤ b.hi
<(a::Interval, b::Interval)  = a.lo < b.lo && a.hi < b.hi

isless(a::Interval, b::Interval)  = a < b
isequal(a::Interval, b::Interval) = a == b

# set relationships
∈(x::T, a::Interval{T}) where T <: Real = a.lo ≤ x ≤ a.hi

⊆(a::Interval, b::Interval) = b.lo ≤ a.lo && a.hi ≤ b.hi
⊂(a::Interval, b::Interval) = (a==b) ? false : a ⊆ b
⊇(a::Interval, b::Interval) = b ⊆ a
⊃(a::Interval, b::Interval) = b ⊂ a

# arithmetic
+(a::Interval) = a
-(a::Interval) = Interval(-a.lo, -a.hi)

+(a::Interval{T}, x::T) where T <: Real = Interval{T}(a.lo+x, a.hi+x)
-(a::Interval{T}, x::T) where T <: Real = Interval{T}(a.lo-x, a.hi-x)

+(x::T, a::Interval{T}) where T <: Real = +a+x
-(x::T, a::Interval{T}) where T <: Real = -a+x

isdisjoint(a::Interval, b::Interval) = b.hi < a.lo || a.hi < b.lo

function ∩(a::Interval, b::Interval)
    try
        return Interval(max(a.lo,b.lo), min(a.hi,b.hi))
    catch err
        if isa(err, ArgumentError)
            return nothing # intervals were disjoint
        end

        throw(err) # unexpected error
    end
end

function \(a::Interval, b::Interval)
    ab = a ∩ b

    isnothing(ab) && return [a]                         # disjoint
    ab == a       && return typeof(a)[]                 # a ⊆ b
    a.lo == ab.lo && return [Interval(ab.hi, a.hi)]     # a ≥ b
    a.hi == ab.hi && return [Interval(a.lo, ab.lo)]     # a ≤ b

    return [Interval(a.lo, b.lo), Interval(b.hi, a.hi)] # a ⊃ b
end

function partition(a::Interval, b::Interval)
	ab = a ∩ b
	isnothing(ab) && return [a] # disjoint
	ab == a		  && return [a] # a ⊆ b
	a.lo == ab.lo && return [ab, Interval(ab.hi, a.hi)]
	a.hi == ab.hi && return [Interval(a.lo, ab.lo), ab]

	return [Interval(a.lo, ab.lo), ab, Interval(ab.hi, a.hi)] # a ⊃ b
end

# -------------------------------------------------------------------------
# interval sets

# ---------------------------------
# helpers

# NOTE: assumes (closed, open) intervals
function merge!(I::Array{Interval{T}}) where T <: Real
    i = 1
    while i < length(I)
        if I[i].hi < I[i+1].lo
            i += 1
            continue
        end

        I[i] = Interval(I[i].lo, I[i+1].hi)
        deleteat!(I, i+1)
    end
end

# ---------------------------------
# typedef

struct IntervalSet{T<:Real}
    Is  :: Array{Interval{T},1}
    min :: T
    max :: T

    function IntervalSet{T}(Is, min, max) where T <: Real
        # ensures atomic intervals are disjoint
        𝕀 = Interval{T}[]
        for I in Is
            n = length(𝕀)
            W = [I]
            for i ∈ 1:n
                W = [x for w in W for x in w \ 𝕀[i]]
            end

            append!(𝕀, W)
        end

        sort!(𝕀)
        merge!(𝕀)

        return new(𝕀, min, max)
    end
end

# ---------------------------------
# constructors

IntervalSet(min::T, max::T, itr)                    where T <: Real    = IntervalSet{T}(itr, min, max)
IntervalSet(min::T, max::T, Is::Array{Interval{T}}) where T <: Real    = IntervalSet{T}(Is, min, max)
IntervalSet(min::T, max::T, I::Interval{T})         where T <: Real    = IntervalSet{T}([I], min, max)

IntervalSet(min::T, max::T)                         where T <: Real    = IntervalSet{T}(Interval{T}[], min, max)
IntervalSet(Is::Array{Interval{T}})                 where T <: Real    = IntervalSet{T}(Is, typemin(T), typemax(T))
IntervalSet(I::Interval{T})                         where T <: Real    = IntervalSet{T}([I], typemin(T), typemax(T))

IntervalSet(min, max, Is::Tuple{T,T}...)            where T <: Real    = IntervalSet{T}((Interval(I) for I in Is), min, max)
IntervalSet(Is::Tuple{T,T}...)                      where T <: Real    = IntervalSet{T}((Interval(I) for I in Is), typemin(T), typemax(T))

copy(I::IntervalSet{T}) where T <: Real = IntervalSet{T}(copy(I.Is), I.min, I.max)

# ---------------------------------
# repl properties

show(io::IO, Is::IntervalSet) = print(io, Is.Is)

# ---------------------------------
# collection properties

length(I::IntervalSet)  = length(I.Is)

iterate(I::IntervalSet)        = length(I) > 0 ? (I.Is[1], 2) : nothing
iterate(I::IntervalSet, state) = (state <= length(I)) ? (I.Is[state], state+1) : nothing

getindex(I::IntervalSet, i) = getindex(I.Is, i)
getindex(A::T, I::IntervalSet) where T <: AbstractArray = [a for i in I for a in A[i]]

min(I::IntervalSet) = I.Is[1].lo
max(I::IntervalSet) = I.Is[end].hi

# ---------------------------------
# operations

# logical
isempty(I::IntervalSet) = length(I) == 0

# arithmetic

function ~(I::IntervalSet)
    starts = (i == 0 ? I.min : I[i].hi for i in 0:length(I))
    stops  = ((i<=length(I)) ? I[i].lo : I.max for i in 1:(length(I)+1))

    return IntervalSet(I.min, I.max, (Interval(start, stop) for (start, stop) in zip(starts,stops) if start < stop))
end

\(I::IntervalSet{T}, j::Interval{T})    where T <: Real = IntervalSet{T}((x for i ∈ I for x in i \ j), I.min, I.max)
\(I::IntervalSet{T}, J::IntervalSet{T}) where T <: Real = foldl(\, J; init=I)

# NOTE: extra generator voodoo avoids creating a list with nothings and then filtering them out
∩(I::IntervalSet{T}, j::Interval{T}) where T <: Real =
    IntervalSet{T}(
        ( x for i ∈ I for x in [i ∩ j] if !isnothing(x) ),
        I.min,
        I.max
    )
∩(j::Interval, I::IntervalSet) = I ∩ j

∩(I::IntervalSet{T}, J::IntervalSet{T}) where T <: Real =
    IntervalSet{T}(
       sort([ x for i ∈ I for j ∈ J for x in [i ∩ j] if !isnothing(x) ]),
       min(I.min,J.min),
       max(I.max,J.max)
    )

function ∪(I::IntervalSet, J::IntervalSet)
    IJ = copy(I)
    K  = J \ I
    for k ∈ K
        ι = searchsortedfirst(IJ.Is, k)
        insert!(IJ.Is, ι, k)
        merge!(IJ.Is)
    end

    return IJ
end

∪(I::IntervalSet, j::Interval) = I ∪ IntervalSet(j)

isdisjoint(A::IntervalSet, b::Interval) = all(isdisjoint(a,b) for a in A)

function containing(A::IntervalSet{T}, x::T) where T <: Real
    for a in A
        if x ∈ a
            return a
        end
    end

    return nothing
end

function containing(A::IntervalSet, b::Interval)
    for a in A
        if b ⊆ a
            return a
        end
    end

    return nothing
end

# -------------------------------------------------------------------------
# tests

function test()
    i = Interval(1, 10)
    j = Interval(0, 11)

    @show i
    @show j

    @show i \ j
    @show j \ i

    I = IntervalSet(0, 50, (0, 4), (2, 5), (7, 8), (10,15), (17,19))
    J = IntervalSet(0, 50, (3, 4), (6, 8), (11,14), (16,18))

    @show I
    @show J

    @show I \ i
    @show I \ j

    @show I \ J

    @show I ∩ j

    @show I ∪ J

    @show ~I

    return nothing
end

end

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
__precompile__()

module Polyhedra

using Compat
using Compat.SparseArrays
using Compat.LinearAlgebra

using GeometryTypes

using MultivariatePolynomials

export PolyhedraLibrary, Polyhedron, FullDim

abstract type PolyhedraLibrary end
abstract type Polyhedron{N,T} <: GeometryPrimitive{N,T} end

using StaticArrays
using StaticArrays.FixedSizeArrays: FixedVector

import MathProgBase
const MPB = MathProgBase
const MPBSI = MPB.SolverInterface

# Similar to StaticArrays.Size
struct FullDim{N}
    function FullDim{N}() where N
        new{N::Int}()
    end
end
fulldim(::FullDim{N}) where N = N
Base.:+(d1::FullDim{N1}, d2::FullDim{N2}) where {N1, N2} = FullDim{N1+N2}()
FullDim(v::AbstractVector) = FullDim{length(v)}()
FullDim(::Union{StaticArrays.SVector{N}, Type{<:StaticArrays.SVector{N}}}) where N = FullDim{N}()
if VERSION >= v"0.7-"
    Base.broadcastable(fd::FullDim) = Ref(fd)
end

MultivariatePolynomials.coefficienttype(::Union{AbstractVector{T}, Type{<:AbstractVector{T}}}) where T = T
similar_type(::Type{<:Vector}, ::FullDim, ::Type{T}) where T = Vector{T}
similar_type(::Type{SparseVector{S, IT}}, ::FullDim, ::Type{T}) where {S, IT, T} = SparseVector{T, IT}
# We define the following method to avoid the fallback to call FullDim(Vector{...})
similar_type(::Type{<:Vector}, ::Type{T}) where T = Vector{T}
similar_type(::Type{SparseVector{S, IT}}, ::Type{T}) where {S, IT, T} = SparseVector{T, IT}
function similar_type(::Type{SAT}, ::FullDim{D}, ::Type{Tout}) where {SAT <: StaticArrays.SVector, D, Tout}
    StaticArrays.similar_type(SAT, Tout, StaticArrays.Size(D))
end

similar_type(::Type{<:Point}, ::FullDim{N}, ::Type{T}) where {N, T} = Point{N,T}
similar_type(::Type{<:Vec}, ::FullDim{N}, ::Type{T}) where {N, T} = Vec{N,T}

emptymatrix(::Type{MT}, m, n) where {MT<:AbstractMatrix} = MT(undef, m, n)
emptymatrix(::Type{SparseMatrixCSC{T, Int}}, m, n) where T = spzeros(T, m, n)
similar_type(::Type{<:Matrix}, ::Type{T}) where T = Matrix{T}
similar_type(::Type{SparseMatrixCSC{S, I}}, ::Type{T}) where {S, I, T} = SparseMatrixCSC{T, I}

# Interface/Definitions
include("elements.jl")
include("comp.jl")
include("representation.jl")
include("indices.jl")
include("incidence.jl")
include("iterators.jl")
include("polyhedron.jl")

# For retro-compatibility with CDD and LRS, remove in v0.3.4
vvectortype(RepT::Type{<:VRep}) = vectortype(RepT)
hvectortype(RepT::Type{<:HRep}) = vectortype(RepT)
# TODO Only define vectortype for the type for Polyhedron v0.4
function vectortype(RepT::Type{<:Polyhedron})
    @assert hvectortype(RepT) == vvectortype(RepT)
    hvectortype(RepT)
end
vectortype(HRepT::Type{<:HRepresentation}) = hvectortype(HRepT)
vectortype(VRepT::Type{<:VRepresentation}) = vvectortype(VRepT)
vectortype(p::Rep) = vectortype(typeof(p))

vectortype(::Type{<:AbstractSparseArray{T}}) where T = SparseVector{T, Int}
vectortype(::Type{<:AbstractMatrix{T}}) where T = Vector{T}

const arraytype = vectortype

hmatrixtype(RepT::Type{<:HRep}, T::Type) = matrixtype(similar_type(hvectortype(RepT), T))
vmatrixtype(RepT::Type{<:VRep}, T::Type) = matrixtype(similar_type(vvectortype(RepT), T))
matrixtype(::Type{<:AbstractVector{T}}) where T = Matrix{T}
matrixtype(::Type{<:AbstractSparseVector{T}}) where T = SparseMatrixCSC{T, Int}

function similar_type(::Type{ET}, ::Type{Tout}) where {Tout, ET<:Union{HRepElement, VRepElement, Rep}}
    similar_type(ET, FullDim(ET), Tout)
end
function similar_type(::Type{ET}, d::FullDim) where {ET<:Union{HRepElement, VRepElement, Rep}}
    similar_type(ET, d, coefficienttype(ET))
end

# Operations
include("repop.jl")
include("center.jl")
include("repelemop.jl")
include("aff.jl")
include("redundancy.jl")
include("projection.jl")

# Implementations of representations
include("vecrep.jl")
include("mixedrep.jl")
include("lphrep.jl")
include("jump.jl")
include("matrep.jl")
include("liftedrep.jl")
include("doubledescription.jl") # FIXME move it after projection.jl once it stops depending on LiftedRep
include("interval.jl") # 1D polyhedron
include("simplepolyhedron.jl")

# Optimization
include("opt.jl")
include("lpqp_to_polyhedra.jl")
include("polyhedra_to_lpqp.jl")
include("vrepsolver.jl")
include("default.jl")

# Visualization
include("show.jl")
include("recipe.jl")
include("decompose.jl")

end # module


text:
export ϵ₀, μ₀, c₀, η₀
export Polarization, TM, TE
export Direction, x̂, ŷ, ẑ
export Point

const ϵ₀ = 8.85418782e-12;
const μ₀ = 1.25663706e-6;
const c₀ = sqrt(1/ϵ₀/μ₀);
const η₀ = sqrt(μ₀/ϵ₀);

const Float  =  Float64;
const Tuple2 =  NTuple{2};

const XX = 1
const YY = 2
const ZZ = 3

@enum Direction x̂=1 ŷ=2 ẑ=3
@enum DerivativeDirection Forward=1 Backward=2
@enum Polarization TM=1 TE=2
@enum FDFDMatSymmetry CSym=1 CNSym=2

struct Point <: FieldVector{2, Float64}
    x::Float64
    y::Float64
end