text:
using CategoricalArrays
using DataFrames
using PooledArrays
using Random

fullgc() = (GC.gc(true); GC.gc(true); GC.gc(true); GC.gc(true))

function run_bench(llen, rlen, arr_type, dup_mode, sort_mode, num_keys, join_type)

    @assert llen > 1000
    @assert rlen > 2000

    @assert arr_type in ["int", "pool", "cat", "str"]
    @assert dup_mode in ["uniq", "dup", "manydup"]
    @assert sort_mode in ["sort", "rand"]
    @assert num_keys in ["1", "2"]
    @assert join_type in ["inner", "left", "right", "outer", "semi", "anti"]

    pad = maximum(length.(string.((llen, rlen))))

    if arr_type == "int"
        if dup_mode == "uniq"
            col1 = [1:llen;]
            col2 = [1:rlen;]
        elseif dup_mode == "dup"
            col1 = repeat(1:llen ÷ 2, inner=2)
            col2 = repeat(1:rlen ÷ 2, inner=2)
        else
            @assert dup_mode == "manydup"
            col1 = repeat(1:llen ÷ 20, inner=20)
            col2 = repeat(1:rlen ÷ 20, inner=20)
        end
    elseif arr_type == "pool"
        if dup_mode == "dup"
            col1 = PooledArray(repeat(string.(1:llen ÷ 2, pad=pad), inner=2))
            col2 = PooledArray(repeat(string.(1:rlen ÷ 2, pad=pad), inner=2))
        else
            @assert dup_mode == "manydup"
            col1 = PooledArray(repeat(string.(1:llen ÷ 20, pad=pad), inner=20))
            col2 = PooledArray(repeat(string.(1:rlen ÷ 20, pad=pad), inner=20))
        end
    elseif arr_type == "cat"
        if dup_mode == "dup"
            col1 = categorical(repeat(string.(1:llen ÷ 2, pad=pad), inner=2))
            col2 = categorical(repeat(string.(1:rlen ÷ 2, pad=pad), inner=2))
        else
            @assert dup_mode == "manydup"
            col1 = categorical(repeat(string.(1:llen ÷ 20, pad=pad), inner=20))
            col2 = categorical(repeat(string.(1:rlen ÷ 20, pad=pad), inner=20))
        end
    else
        @assert arr_type == "str"
        if dup_mode == "uniq"
            col1 = string.(1:llen, pad=pad)
            col2 = string.(1:rlen, pad=pad)
        elseif dup_mode == "dup"
            col1 = repeat(string.(1:llen ÷ 2, pad=pad), inner=2)
            col2 = repeat(string.(1:rlen ÷ 2, pad=pad), inner=2)
        else
            @assert dup_mode == "manydup"
            col1 = repeat(string.(1:llen ÷ 20, pad=pad), inner=20)
            col2 = repeat(string.(1:rlen ÷ 20, pad=pad), inner=20)
        end
    end

    Random.seed!(1234)

    if sort_mode == "rand"
        shuffle!(col1)
        shuffle!(col2)
    else
        @assert sort_mode == "sort"
    end

    joinfun = Dict("inner" => innerjoin, "left" => leftjoin,
                    "right" => rightjoin, "outer" => outerjoin,
                    "semi" => semijoin, "anti" => antijoin)[join_type]

    if num_keys == "1"
        df1 = DataFrame(id1=col1)
        df2 = DataFrame(id1=col2)
        joinfun(df1[1:1000, :], df2[1:2000, :], on=:id1)
        joinfun(df2[1:2000, :], df1[1:1000, :], on=:id1)
        fullgc()
        t1 = @elapsed joinfun(df1, df2, on=:id1)
        fullgc()
        t2 = @elapsed joinfun(df2, df1, on=:id1)
    else
        @assert num_keys == "2"
        df1 = DataFrame(id1=col1, id2=col1)
        df2 = DataFrame(id1=col2, id2=col2)
        joinfun(df1[1:1000, :], df2[1:2000, :], on=[:id1, :id2])
        joinfun(df2[1:2000, :], df1[1:1000, :], on=[:id1, :id2])
        fullgc()
        t1 = @elapsed joinfun(df1, df2, on=[:id1, :id2])
        fullgc()
        t2 = @elapsed joinfun(df2, df1, on=[:id1, :id2])
    end

    return (t1 + t2) / 2.0
end

if abspath(PROGRAM_FILE) == @__FILE__
    @info "$ARGS"

    llen = parse(Int, ARGS[1])
    rlen = parse(Int, ARGS[2])
    t = run_bench(llen, rlen, ARGS[3], ARGS[4], ARGS[5], ARGS[6], ARGS[7])
    @info("Timing: $t sec")
end


text:
#!/usr/bin/env julia

using Luxor, Colors

function draw_luxor_demo(fname)
    Drawing(1200, 1400, fname)

    origin() # move 0/0 to center
    background("purple")

    setopacity(0.7)     # opacity from 0 to 1
    sethue(0.3,0.7,0.9) # sethue sets the color but doesn't change the opacity
    setline(20)         # line width
    # graphics functions use :fill, :stroke, :fillstroke, :clip, or leave blank
    rect(-400,-400,800,800, :fill)
    randomhue()
    circle(0, 0, 460, :stroke)

    # clip the following graphics to a circle positioned above the x axis
    circle(0,-200,400,:clip)
    sethue("gold")
    setopacity(0.7)
    setline(10)

    # simple line drawing
    for i in 0:pi/36:2pi - pi/36
        move(0, 0)
        line(cos(i) * 600, sin(i) * 600 )
        stroke()
    end

    # finish clipping
    clipreset()

    # here using Mac OS X fonts
    fontsize(60)
    setcolor"turquoise"
    fontface("Optima-ExtraBlack")
    textwidth = textextents("Luxor")[5]
    # move the text by half the width
    textcentred("Luxor", 0, currentdrawing.height/2 - 400)

    # text on curve
    fontsize(18)
    fontface("Avenir-Black")

    textcurve("THIS IS TEXT ON A CURVE " ^ 14, 0, 550, O)
    finish()
    println("finished test: output in $(fname)")
end

draw_luxor_demo("/tmp/luxor-test1.png") # PNG or PDF file name

ext: jl


text:
using HDF5 # To have access to .hf5 files
using Plots # To be able to plot data

########################################
# Function that read the .hf5 file
########################################

namefile = "../../data/Dump_Diffusion_Coefficients.hf5"

"""
    openData(namefile)

Recovers the data from the .hf5 file `namefile` and returns it as a 4-uple.
In the order nbj, aMeasure, tabj, tabDRRjj.
"""
function openData(namefile)
    file = h5open(namefile, "r")
    nbj = read(file,"nbj")
    nba = read(file,"nba")
    tabaj = read(file,"tabaj")
    tabj = read(file,"tabj")
    taba = read(file,"taba")
    tabDRRjj = read(file,"tabDRRjj")
    close(file)
    return nbj, nba, tabaj, tabj, taba, tabDRRjj
end

########################################
# Getting the data
########################################
println("Recovering plot data...")
nbj, nba, tabaj, tabj, taba, tabDRRjj = openData(namefile)

########################################
# Converting the data into log-scaling
########################################
println("Converting data into log-log scaling...")
for ij=1:nbj
    tabj[ij] = log10(tabj[ij])
end

for ia=1:nba
    taba[ia] = log10(taba[ia])
end

for iGrid=1:nba*nbj
    if (tabDRRjj[iGrid] == 0) # Within loss-cone, set an arbitrary low value for null
        tabDRRjj[iGrid] = -10
    else
        tabDRRjj[iGrid] = log10(tabDRRjj[iGrid])
    end
end

########################################
# Plotting the data
########################################
println("Plotting the data...")
p = Plots.contourf(tabj,taba,tabDRRjj)
Plots.savefig(p,"../../graphs/Julia/Djj.png") # Saves the figure
Plots.display(p)
readline()


text:
@inline function Base.getproperty(ϕ::GridDeformation, name::Symbol)
    if name === :u
        return getfield(ϕ, :u)
    elseif name === :nodes
        return getfield(ϕ, :nodes)
    elseif name === :knots
        Base.depwarn("the field name of `GridDeformation` has changed from `knots` to `nodes`.", :getproperty)
        return getfield(ϕ, :nodes)
    end
    error("field ", name, " is not available")
end

function GridDeformation(u::AbstractArray{FV,N},
                         sz::NTuple{N,L}) where {FV<:SVector,N,L<:Integer}
    Base.depwarn("`GridDeformation(u, size(fixed))` is deprecated, use `GridDeformation(u, axes(fixed))` instead.", :GridDeformation)
    return GridDeformation(u, map(Base.OneTo, sz))
end

# The above causes an ambiguity, resolve it
function GridDeformation(u::AbstractArray{FV,0}, nodes::Tuple{}) where FV<:SVector
    error("node tuple cannot be empty")
end

function GridDeformation(u, nodes::AbstractVector{<:Integer})
    Base.depwarn("`GridDeformation(u, [size(fixed)...])` is deprecated, pass the axes of `fixed` instead.", :GridDeformation)
    GridDeformation(u, map(Base.OneTo, (nodes...,)))
end

function tform2deformation(tform::AffineMap{M,V}, arraysize::DimsLike, gridsize) where {M,V}
    Base.depwarn("""
    `tform2deformation(tform, arraysize, gridsize)` is deprecated.
    Formerly, `tform` was defined so as to apply to `centered(img)`, which shifts
    the axes of `img` to put 0 at the midpoint of `img`.
    Now you should call this as `tform2deformation(tform, axs, gridsize)`, where
    `axs` represents the desired domain of `tform`.
    If you are transitioning old code, either use `axs = axes(fixed)` if `fixed`
    is already centered, or `axs = centeraxes(axes(fixed))` if not.
    """, :tform2deformation)
    axs = map((arraysize...,)) do sz
        halfsz = sz ÷ 2
        IdentityUnitRange(1-halfsz:sz-halfsz)
    end
    return tform2deformation(tform, axs, (gridsize...,))
end

import Base: getindex
@deprecate getindex(ϕ::GridDeformation{T,N,A}, xs::Vararg{Number,N}) where {T,N,A<:AbstractInterpolation} ϕ(xs...)

Base.@deprecate_binding eachknot eachnode
Base.@deprecate_binding knotgrid nodegrid


text:
using HDF5

function save_params(file::HDF5File, rbm::RBM, name::AbstractString)
    write(file, "$(name)___weight", rbm.W')
    write(file, "$(name)___vbias", rbm.vbias)
    write(file, "$(name)___bias", rbm.hbias)
end

function save_params(file::HDF5File, rbm::RBM, name::AbstractString)
    write(file, "$(name)___weight", rbm.W')
    write(file, "$(name)___vbias", rbm.vbias)
    write(file, "$(name)___bias", rbm.hbias)
end

function save_params(filename::AbstractString,rbm::RBM,name::AbstractString)
    h5open(filename,"w") do file
        save_params(file,rbm,name)
    end
end

function append_params(filename::AbstractString,rbm::RBM,name::AbstractString)
    h5open(filename,"r+") do file
        save_params(file,rbm,name)
    end
end

function load_params(file::HDF5File, rbm::RBM, name::AbstractString)
    rbm.W = read(file, "$(name)___weight")'
    rbm.vbias = read(file, "$(name)___vbias")
    rbm.hbias = read(file, "$(name)___bias")
end

function save_params(file::HDF5File, net::Net)
    for i=1:length(net)
        save_params(file, net[i], getname(net, i))
    end
end
save_params(path::AbstractString, net::Net) = h5open(path, "w") do h5
    save_params(h5, net)
end

function load_params(file, net::Net)
    for i=1:length(net)
        load_params(file, net[i], getname(net, i))
    end
end
load_params(path::AbstractString, net::Net) = h5open(path) do h5
    load_params(h5, net)
end

function SaveMonitorHDF5(mon::Monitor,filename::AbstractString)
    h5open(filename , "w") do file
        # write(file, "LastIndex", mon.LastIndex)
        # write(file, "UseValidation", mon.UseValidation)
        write(file, "MonitorEvery", mon.MonitorEvery)
        # write(file, "MonitorVisual", mon.MonitorVisual)
        # write(file, "MonitorText", mon.MonitorText)
        write(file, "Epochs", mon.Epochs)
        write(file, "LearnRate", mon.LearnRate)
        write(file, "Momentum", mon.Momentum)
        write(file, "PseudoLikelihood", mon.PseudoLikelihood)
        write(file, "TAPLikelihood", mon.TAPLikelihood)
        write(file, "ValidationPseudoLikelihood", mon.ValidationPseudoLikelihood)
        write(file, "ValidationTAPLikelihood", mon.ValidationTAPLikelihood)
        write(file, "ReconError", mon.ReconError)
        write(file, "ValidationReconError", mon.ValidationReconError)
        write(file, "BatchTime_µs", mon.BatchTime_µs)
        # write(file, "FigureHandle", mon.FigureHandle)
    end
end


text:
using Distributions
using DataStructures

f(x) = x^2 - 5 * x + 6

function filter_best_solutions(solutions; pa=0.25)
    new_best_solutions=SortedDict(Dict{Float64, Float64}())
    println("***FILTER RCVD : ", length(solutions))
    println("Filterting : ", floor((length(solutions) - length(solutions)*pa)))
    let
        i = 1
        for (key, value) in solutions
            new_best_solutions[key] = value
            i = i + 1
            if i >= (length(solutions) - length(solutions)*pa)
                return merge(new_best_solutions, initialize_population(floor((length(solutions) - length(solutions)*pa))))
            end
        end
    end
end

function get_best_solution(soln)
    best = soln[1]
    for s in soln
        if f(s) < f(best)
            best = s
        end
    end
    return best
end

heaviside(x::AbstractFloat) = ifelse(x < 0, zero(x), ifelse(x > 0, one(x), oftype(x,0.5)))

function generate_new_population(population; pa=0.25)
    new_population=SortedDict(Dict{Float64, Float64}())
    for (key, value) in population
        new_value = value + heaviside(pa - rand()) * (rand(population)[2] - rand(population)[2])
        new_population[f(new_value)] = new_value
    end
    return new_population
end

function initialize_population(total)
    population = SortedDict(Dict{Float64,Float64}())
    for i=1:total
        soln = rand(Float64)
        population[f(soln)] = soln
    end
    return population
end

function cuckoo_search(f;total_cuckoos=25,pa = 0.25,α=0.01, λ=0.5, MaxGenerations = 15)
    cuckoos = initialize_population(total_cuckoos)
    println("Initial population of Cuckoos: ", length(cuckoos))
    L = Levy(λ)
    best = nothing
    t=1
    while t < MaxGenerations
        r_cuckoo = rand(cuckoos)
        new_best = r_cuckoo[2] + α * rand(L)
        if (f(new_best) < r_cuckoo[2])
            best = new_best
        end
        #
        cuckoos = filter_best_solutions(cuckoos, pa=pa)
        println("After FILTER Total Cuckoos: ", length(cuckoos))
        cuckoos = generate_new_population(cuckoos, pa=pa)
        best = first(cuckoos)[2]
        println("Current Best: ", best)
        println("Total Cuckoos: ", length(cuckoos))
        t = t +1
    end
    return best
end

global_minimum=cuckoo_search(f)

println("Final Global Minimum: ", global_minimum)


text:
#########################
### Low rank matrices ###
#########################

struct LowRankMat
    eigenvalues::Vector{Arb}
    leftevs::Vector{ArbRefMatrix}
    rightevs::Vector{ArbRefMatrix}
end
function LowRankMat(eigenvalues, eigenvectors)
    LowRankMat(eigenvalues, eigenvectors, eigenvectors)
end

Base.getindex(A::LowRankMat,i,j) = sum(A.eigenvalues[r] * A.leftevs[r][j] * A.rightevs[r][i] for r=1:length(A.eigenvalues))

function Base.size(m::LowRankMat, i)
    length(m.leftevs[1])
end

function convert_to_prec(A::LowRankMat,prec=precision(BigFloat))
    evs = [Arb(ev,prec=prec) for ev in A.eigenvalues]
    leftvecs = [ArbRefMatrix(vec,prec=prec) for vec in A.leftevs]
    rightvecs = [ArbRefMatrix(vec,prec=prec) for vec in A.rightevs]
    for i=1:length(evs)
        Arblib.get_mid!(evs[i],evs[i])
        Arblib.get_mid!(leftvecs[i],leftvecs[i])
        Arblib.get_mid!(rightvecs[i],rightvecs[i])
    end
    return LowRankMat(evs,leftvecs,rightvecs)
end

######################################
### Completely sampled polynomials ###
######################################

struct SampledMPolyRing{T} <: AbstractAlgebra.Ring
    base_ring
end

SampledMPolyRing(ring::T) where T<:AbstractAlgebra.Ring = SampledMPolyRing{elem_type(ring)}(ring)

==(a::SampledMPolyRing{T}, b::SampledMPolyRing{T}) where T = true # a.base_ring == b.base_ring TODO: this is a workaround

struct SampledMPolyElem <: AbstractAlgebra.RingElem
    evaluations::Dict
end

function SampledMPolyElem(p::MPolyElem, samples)
    SampledMPolyElem(Dict(samples[s] => p(samples[s]...) for s in eachindex(samples)))
end

function SampledMPolyElem(samples, values)
    SampledMPolyElem(Dict(samples[s] => values[s] for s in eachindex(samples)))
end

function evaluate(p::SampledMPolyElem, v...)
    p.evaluations[collect(v)]
end

function (p::SampledMPolyElem)(v...)
    evaluate(p, v...)
end

# addition
function Base.:+(p::SampledMPolyElem, q::SampledMPolyElem)
    @assert keys(p.evaluations) == keys(q.evaluations)
    SampledMPolyElem(Dict(sample => p.evaluations[sample] + q.evaluations[sample] for sample in keys(p.evaluations)))
end

function Base.:+(p::MPolyElem, q::SampledMPolyElem)
    r = copy(q)
    for s in keys(r.evaluations)
        r.evaluations[s] += p(s...)
    end
    r
end

function Base.:+(q::SampledMPolyElem, p::MPolyElem)
    p + q
end

# subtraction
function Base.:-(q::SampledMPolyElem)
    SampledMPolyElem(Dict(sample => -q.evaluations[sample] for sample in keys(q.evaluations)))
end

function Base.:-(p::T, q::SampledMPolyElem) where T <: Union{SampledMPolyElem, MPolyElem}
    p + (-q)
end

function Base.:-(p::SampledMPolyElem, q::MPolyElem)
    p + (-q)
end

# multiplication
function Base.:*(p::MPolyElem, q::SampledMPolyElem)
    r = copy(q)
    for s in keys(r.evaluations)
        r.evaluations[s] *= p(s...)
    end
    r
end

function Base.:*(p, q::SampledMPolyElem)
    r = copy(q)
    for s in keys(r.evaluations)
        r.evaluations[s] *= p(s...)
    end
    r
end

function Base.:*(q::SampledMPolyElem, p::MPolyElem)
    p * q
end

function Base.:*(p::SampledMPolyElem, q::SampledMPolyElem)
    @assert keys(p.evaluations) == keys(q.evaluations)
    r = copy(q)
    for s in keys(r.evaluations)
        r.evaluations[s] *= p(s...)
    end
    r
end

function Base.copy(q::SampledMPolyElem)
    SampledMPolyElem(copy(q.evaluations))
end

"""
    approximatefekete(basis, samples) -> basis, samples

Compute approximate fekete points based on samples and a corresponding orthogonal basis.
The basis consists of sampled polynomials, sampled at `samples`

This preserves a degree ordering of `basis` if present.
"""
function approximatefekete(basis, samples;show_det=false,s=3)
    V, P, samples = approximate_fekete(samples, basis,show_det=show_det,s=s)
    [SampledMPolyElem(samples, V[:,p]) for p in eachindex(basis)], samples
end

###################################
### Low rank matrix polynomials ###
###################################

polys = Union{MPolyElem, SampledMPolyElem}

"""
    LowRankMat(eigenvalues::Vector, rightevs::Vector{Vector}[, leftevs::Vector{Vector}])

The matrix ``∑_i λ_i v_i w_i^T``.

If `leftevs` is not specified, use `leftevs = rightevs`.
"""
struct LowRankMatPol
    eigenvalues::Vector{polys}
    leftevs::Vector{Vector{polys}}
    rightevs::Vector{Vector{polys}}
end

function LowRankMatPol(eigenvalues, eigenvectors)
    LowRankMatPol(eigenvalues, eigenvectors, eigenvectors)
end
LinearAlgebra.transpose(A::Union{LowRankMatPol,LowRankMat}) = typeof(A)(A.eigenvalues,A.rightevs,A.leftevs)
Base.getindex(A::LowRankMatPol,i,j) = sum(A.eigenvalues[r] * A.leftevs[r][j] * A.rightevs[r][i] for r=1:length(A.eigenvalues))

function Base.size(m::LowRankMatPol, i)
    length(m.leftevs[1])
end

function AbstractAlgebra.evaluate(p::LowRankMatPol, sample...; scaling = 1, prec=precision(BigFloat))
    LowRankMat([Arb(scaling*v(sample...), prec=prec) for v in p.eigenvalues],
               [ArbRefMatrix([v(sample...) for v in w], prec=prec) for w in p.leftevs],
               [ArbRefMatrix([v(sample...) for v in w], prec=prec) for w in p.rightevs])
end

#############################
### Low rank SOS problems ###
#############################

"""
    Block(l::Any[, r::Int, s::Int])

Specifies a block corresponding to the positive semidefinite variable `l`.

Specifying `r,s` makes the Block correspond to the `r,s` subblock of the variable `l`.
"""
struct Block
    l::Any
    r::Int
    s::Int
end

function Block(l)
    Block(l, 1, 1)
end

"""
    Constraint(constant, matrixcoeff, freecoeff, samples[, scalings])

Models a polynomial constaint of the form
```math
    c(x) = ∑_l ⟨ Y^l_{r,s}, A^l_{r,s}(x) ⟩ + ∑_i y_i B_i(x)
```
with samples, where `r,s` are defined by the `Block` structure with `l` as first argument

Arguments:
  - `constant::MPolyElem`
  - `matrixcoeff::Dict{Block, LowRankMatPol}`
  - `freecoeff::Dict{Any, MPolyElem}`
  - `samples::Vector{Vector}`
  - `scalings::Vector`
"""
struct Constraint
    constant::MPolyElem
    matrixcoeff::Dict{Block,LowRankMatPol}
    freecoeff::Dict{Any,MPolyElem}
    samples::Vector{Vector}
    scalings::Vector
end
Constraint(constant,matrixcoeff,freecoeff,samples) = Constraint(constant,matrixcoeff,freecoeff,samples,[1 for s in samples])

"""
    Objective(constant, matrixcoeff::Dict{Block, Matrix}, freecoeff::Dict)

The objective for the LowRankPolProblem
"""
struct Objective
    constant
    matrixcoeff::Dict{Block,Matrix}
    freecoeff::Dict
end

"""
    LowRankPolProblem(maximize, objective, constraints)

Combine the objective and constraints into a low-rank polynomial problem
"""
struct LowRankPolProblem #Maybe we  need to rename this to e.g. LowRankPolProblem
    maximize::Bool
    objective::Objective
    constraints::Vector{Constraint}
end

##############################
### Clustered low rank SDP ###
##############################

struct ClusteredLowRankSDP
    maximize::Bool
    constant::Arb #constant is part of the objective, so maybe we should have constant, C, b, A, B, c or something like that.
    A::Vector{Vector{Matrix{Dict{Int,LowRankMat}}}}         # A[j][l][r,s][p]
    B::Vector{ArbRefMatrix}                                 # B[j]
    c::Vector{ArbRefMatrix}                                 # c[j]
    C::BlockDiagonal{ArbRef,BlockDiagonal{ArbRef,ArbRefMatrix}} # C[j][l]
    b::ArbRefMatrix                                         #b
    matrix_coeff_names::Vector{Vector{Any}}                 # the (user defined) names of the blocks
    free_coeff_names::Vector{Any}                           # the (user defined) names of the free variables
end

function ClusteredLowRankSDP(maximize,constant,A,B,c,C,b)
    matrixcoeff_names = [[(j,l) for l in eachindex(A[j])] for j in eachindex(A)]
    freecoeff_names = collect(1:size(b,1))
    return ClusteredLowRankSDP(maximize,constant,A,B,c,C,b,matrixcoeff_names,freecoeff_names)
end

"""
    ClusteredLowRankSDP(sos::LowRankPolProblem; as_free::Vector, prec, verbose])

Define a ClusteredLowRankSDP based on the LowRankPolProblem sos.

The PSD variables defined by the keys in the vector `as_free` will be modelled as extra free variables,
with extra constraints to ensure that they are equal to the entries of the PSD variables.
Remaining keyword arguments:
  - `prec` (default: precision(BigFloat)) - the precision of the result
  - `verbose` (default: false) -  print progress to the standard output
"""
function ClusteredLowRankSDP(sos::LowRankPolProblem; as_free = [], prec=precision(BigFloat),verbose = false)
    # the blocks in as_free will be modelled as extra variables
    if length(as_free) > 0
        if verbose
            println("Modelling the selected PSD variables as free variables...")
        end
        #We do not modify the sos problem
        sos = deepcopy(sos)
        for l in as_free
            #For each constraint, for each Block b with l==b.l, we need to
            #   1) add the entries of the LowRankMatPol to the freecoeff dict with key b,i,j
            #   2) remove the block from the constraint (from the dict matrixcoeff)
            #For each block we use, we need to:
            #   1) add a new constraint equating the free variable with the entry of the block
            m = 0
            n = 0
            for constraint in sos.constraints
                for block in collect(keys(constraint.matrixcoeff))
                    if block.l == l #we do this l now
                        if n == 0
                            n = size(constraint.matrixcoeff[block],1)
                        elseif n != size(constraint.matrixcoeff[block],1)
                            error("The blocks are not of equal sizes")
                        end
                        for i=1:size(constraint.matrixcoeff[block],1)
                            for j=1:size(constraint.matrixcoeff[block],2)
                                if block.r == block.s && i >= j
                                    #diagonal (r,s) block, so we only need to do things when i >= j
                                    if i == j # just the entry
                                        constraint.freecoeff[(l,(block.r-1)*n+i,(block.s-1)*n+j)] = constraint.matrixcoeff[block][i,j]
                                    else #i > j, so we need to take both sides of the diagonal into account
                                        constraint.freecoeff[(l,(block.r-1)*n+i,(block.s-1)*n+j)] = 2*constraint.matrixcoeff[block][i,j]
                                    end
                                elseif block.r > block.s #the other case is r < s, but because it is symmetric we do that now too. (hence the factor 2)
                                    constraint.freecoeff[(l,(block.r-1)*n+i,(block.s-1)*n+j)] = 2 * constraint.matrixcoeff[block][i,j]
                                end
                            end
                        end
                        m = max(block.r,block.s,m)
                        delete!(constraint.matrixcoeff,block)
                    end
                end
            end

            # add equality constraints for the new free variables and the new blocks
            R, (x,) = PolynomialRing(RealField,["x"])
            for i=1:n*m
                for j=1:i
                    if i==j
                        push!(sos.constraints,Constraint(R(0),
                            Dict(Block(l,i,j)=>LowRankMatPol([R(1)],[[R(1)]])), # matrix variable
                            Dict((l,i,j)=>R(-1)), # free variable
                            [[0]],
                            ))
                    else # i != j, so both sides get 1/2
                        push!(sos.constraints,Constraint(R(0),
                            Dict(Block(l,i,j) => LowRankMatPol([R(1)],[[R(1)]]), Block(l,j,i) => LowRankMatPol([R(1)],[[R(1)]])),
                            Dict((l,i,j)=>R(-2)),
                            [[0]],
                        ))
                    end
                end
            end

            # change the objective
            # Two options:
            # 1) hvcat the blocks and set an objective for the new subblocks (PSD side)
            # 2) set an objective for the free variables.
            # 3) set both to 1/2 * objective
            #This is an implementation of (1)
            new_blocks = Dict() #because we delete the blocks afterwards, so we need to store the new blocks separately and put them in the objective afterwards
            for block in collect(keys(sos.objective.matrixcoeff)) # now this is not modified during the loop
                if block.l == l && block.r >= block.s # we only do the lower half, and we make it symmetric
                    for i=1:n
                        for j=1:(block.r == block.s ? i : n)
                            if (block.r-1)*n+i == (block.s-1)*n + j # on the diagonal
                                new_blocks[Block(l, (block.r-1)*n+i, (block.s-1)*n+j)] = hcat(sos.objective.matrixcoeff[block][i,j])
                            else #not on the diagonal of the full matrix, so we need to take both sides into account
                                new_blocks[Block(l, (block.r-1)*n+i, (block.s-1)*n+j)] = hcat(sos.objective.matrixcoeff[block][i,j])
                                new_blocks[Block(l, (block.s-1)*n+j, (block.r-1)*n+i)] = hcat(sos.objective.matrixcoeff[block][i,j])
                            end
                        end
                    end
                end
                if block.l == l
                    delete!(sos.objective.matrixcoeff,block)
                end
            end
            for (block, value) in new_blocks
                sos.objective.matrixcoeff[block] = value
            end

        end
    end
    if verbose
        println("Forming the clusters...")
    end
    # clusters will be a vector of vectors containing block labels (so block.l)
    #NOTE: if a constraint doesn't have PSD variables, this creates an empty cluster.
    clusters = Vector{Any}[]
    for constraintindex in eachindex(sos.constraints)
        constraint = sos.constraints[constraintindex]
        clusterset = Set{Int}()
        for k in keys(constraint.matrixcoeff)
            for i in eachindex(clusters)
                if k.l in clusters[i]
                    push!(clusterset, i)
                end
            end
        end
        clusterset = sort(collect(clusterset))
        if isempty(clusterset)
            if length(keys(constraint.matrixcoeff)) > 0 #only create a new cluster for constraints with matrix blocks
                push!(clusters, [k.l for k in keys(constraint.matrixcoeff)])
            end
        else
            #remove the clusters which we need to combine from 'clusters', and merge them
            u = union(splice!(clusters, clusterset)...)
            append!(u, [k.l for k in keys(constraint.matrixcoeff)])
            push!(clusters, u)
        end
    end

    #NOTE: if a constraint doesn't use PSD variables, this will not assign it to a cluster.
    # find out which constraints belong to which cluster
    # we build a vector of vectors containing constraint indices
    cluster_constraint_index = [Int[] for _ in clusters]
    remaining = Int[]
    for constraintindex in eachindex(sos.constraints)
        constraint = sos.constraints[constraintindex]
        k = first(keys(constraint.matrixcoeff))
        found_cluster = false
        for clusterindex in eachindex(clusters)
            if k.l in clusters[clusterindex]
                push!(cluster_constraint_index[clusterindex], constraintindex)
                found_cluster = true
                break
            end
        end
        if !found_cluster
            #this is a constraint without cluster
            push!(remaining,constraintindex)
        end
    end
    if length(remaining)>0
        # there is at least on constraint which doesn't use PSD variables, we put it in a separate cluster (otherwise we would forget about it)
        #NOTE: this is untested, since we did not encounter this in any problem
        push!(cluster_constraint_index, remaining)
    end
    if verbose
        println("Sampling the low rank polynomial matrices...")
    end
    A = []
    C = []
    matrix_coeff_names = []
    for clusterindex in eachindex(cluster_constraint_index)
        constraintindices = cluster_constraint_index[clusterindex]
        nsubblocks = Dict{Any,Int}()
        subblocksizes = Dict{Any,Int}()
        for constraintindex in constraintindices
            constraint = sos.constraints[constraintindex]
            for block in keys(constraint.matrixcoeff)
                subblocksizes[block.l] = max(size(constraint.matrixcoeff[block], 1),get(subblocksizes,block.l,0))
                nsubblocks[block.l] = max(block.r, block.s, get(nsubblocks, block.l, 0))
            end
        end
        v = [[Dict{Int,LowRankMat}() for _=1:m, _=1:m] for (k, m) in nsubblocks]
		subblock_keys = collect(keys(nsubblocks))

        i = 1
        for constraintindex in constraintindices
            constraint = sos.constraints[constraintindex]
            for (idx,sample) in enumerate(constraint.samples) #This can take quite long for a large amount of samples & a large basis
                for block in keys(constraint.matrixcoeff)
                    v[findfirst(isequal(block.l), subblock_keys)][block.r,block.s][i] = evaluate(constraint.matrixcoeff[block], sample...,scaling = constraint.scalings[idx],prec=prec)
                    # We can set the s,r block here to the transpose of the r,s block. Or to 1/2 * the transpose if we would want that.
                end
                i += 1
            end
        end

        push!(A, v)

        M = [[ArbRefMatrix(subblocksizes[l], subblocksizes[l],prec=prec) for r=1:m, s=1:m] for (l, m) in nsubblocks]
        for block in keys(sos.objective.matrixcoeff)
            if block.l in clusters[clusterindex]
                M[findfirst(isequal(block.l), subblock_keys)][block.r, block.s] = ArbRefMatrix(sos.objective.matrixcoeff[block],prec=prec)
            end
        end
        push!(C, [ArbRefMatrix(hvcat(size(S,2), [S[r,s] for r=1:size(S,1) for s=1:size(S,2)]...),prec=prec) for S in M])
        push!(matrix_coeff_names, subblock_keys)
    end
    C_block = BlockDiagonal([BlockDiagonal([b for b in blocks]) for blocks in C])

    if verbose
        println("Sampling the coefficients of the free variables...")
    end

    freecoefflabels = Set()
    for constraint in sos.constraints
        union!(freecoefflabels, keys(constraint.freecoeff))
    end
    freecoefflabels = collect(freecoefflabels) #We want an ordering for B and b

    # B is a vector of the matrices B^j
    B = ArbRefMatrix[]
    for constraintindices in cluster_constraint_index
        Bj = []
        for constraintindex in constraintindices
            constraint_B = ArbRefMatrix(length(sos.constraints[constraintindex].samples),length(freecoefflabels),prec=prec)
            for s_idx in eachindex(sos.constraints[constraintindex].samples)
                s = sos.constraints[constraintindex].samples[s_idx]
                scaling = sos.constraints[constraintindex].scalings[s_idx]
                for (k,v) in sos.constraints[constraintindex].freecoeff
                    constraint_B[s_idx,findfirst(isequal(k),freecoefflabels)] = Arb(scaling*v(s...))
                end
            end
            push!(Bj,constraint_B)
        end
        push!(B,vcat(Bj...))
    end
    if verbose
        println("Sampling the constants...")
    end
    c = [ArbRefMatrix([sos.constraints[constraintindex].scalings[s_idx]*sos.constraints[constraintindex].constant(sample...) for constraintindex in constraintindices for (s_idx,sample) in enumerate(sos.constraints[constraintindex].samples)],prec=prec) for constraintindices in cluster_constraint_index]

    b = ArbRefMatrix(length(freecoefflabels),1,prec=prec)
    for (k,v) in sos.objective.freecoeff
        b[findfirst(isequal(k), freecoefflabels)] = Arb(v)
    end
    ClusteredLowRankSDP(sos.maximize, Arb(sos.objective.constant,prec=prec), A, B, c, C_block, b, matrix_coeff_names,freecoefflabels)
end

function convert_to_prec(sdp::ClusteredLowRankSDP,prec=precision(BigFloat))
    # We convert everything to the new precision
    new_A = similar(sdp.A)
    new_B = similar(sdp.B)
    new_C = similar(sdp.C)
    new_c = similar(sdp.c)
    new_constant = Arb(sdp.constant,prec=prec)
    for j in eachindex(sdp.A)
        new_A[j] = similar(sdp.A[j])
        for l in eachindex(sdp.A[j])
            new_A[j][l] = similar(sdp.A[j][l])
            for r=1:size(sdp.A[j][l],1)
                for s=1:size(sdp.A[j][l],2)
                    new_A[j][l][r,s] = Dict{Int,LowRankMat}()

                    for p in keys(sdp.A[j][l][r,s])
                        new_A[j][l][r,s][p] = convert_to_prec(sdp.A[j][l][r,s][p], prec)
                    end
                end
            end
            new_C.blocks[j].blocks[l] = ArbRefMatrix(sdp.C.blocks[j].blocks[l],prec=prec)
            Arblib.get_mid!(new_C.blocks[j].blocks[l], new_C.blocks[j].blocks[l])
        end
        new_B[j] = ArbRefMatrix(sdp.B[j], prec=prec)
        Arblib.get_mid!(new_B[j], new_B[j])

        new_c[j] = ArbRefMatrix(sdp.c[j], prec=prec)
        Arblib.get_mid!(new_c[j], new_c[j])
    end

    new_b = ArbRefMatrix(sdp.b, prec=prec)
    Arblib.get_mid!(new_b, new_b)

    return ClusteredLowRankSDP(sdp.maximize,new_constant,new_A,new_B,new_c,new_C,new_b,sdp.matrix_coeff_names,sdp.free_coeff_names)
end

######################
### Solver results ###
######################

# Define the possible status
abstract type Status end
struct Optimal <: Status end
struct NearOptimal <: Status end
struct PrimalFeasible <: Status end
struct DualFeasible <: Status end
struct Feasible <: Status end
struct NotConverged <: Status end

optimal(::Status) = false
optimal(::Optimal) = true

#The name to make clear what results these are. Distinguishing it from e.g. SDPA-GMP if someone uses both?
struct CLRSResults
    x::ArbRefMatrix
    X::BlockDiagonal{ArbRef,BlockDiagonal{ArbRef,ArbRefMatrix}}
    y::ArbRefMatrix
    Y::BlockDiagonal{ArbRef,BlockDiagonal{ArbRef,ArbRefMatrix}}
    primal_objective::Arb
    dual_objective::Arb
    matrix_coeff_names::Vector{Vector{Any}}
    free_coeff_names::Vector{Any}
    # status::Status #optimal, near optimal, primal feasible, dual feasible, feasible (p&d),not converged
end

Base.show(io::IO, x::Optimal) = @printf(io, "pdOpt")
Base.show(io::IO, x::NearOptimal) = @printf(io,"NearOpt")
Base.show(io::IO, x::Feasible) = @printf(io,"pdFeas")
Base.show(io::IO, x::PrimalFeasible) = @printf(io,"pFeas")
Base.show(io::IO, x::DualFeasible) = @printf(io,"dFeas")
Base.show(io::IO, x::Status) = @printf(io,"NOINFO")


text:
#
function f1()
    sleep(0.01)  # sleep for 10ms
    return 0
end
function f2()
    sleep(0.005)  # sleep for 1ms
    return 0
end

function g1()
    sleep(0.001)  # sleep for 1ms
end
function g2()
    sleep(0.05)  # sleep for 10ms
end

# Simulate the outer loop
function h1()
    for i in 1:100
        f1()
        g1()
    end
    println("100 iterations")
    return 0
end

function h2()
    for i in 1:20
        f2()
        g2()
    end
    println("20 iterations")
    return 0
end


text:
using Plots
using LinearAlgebra
using Polynomials

using SparseArrays
using Arpack

using ForwardDiff

using BenchmarkTools

#using SymPy
#komvos1 komvos2 komvos3
xy = [ 0 0;
       0 1]; #orismos komvon san xyz
nodes = size(xy)[2]; #arithmos komvon

#komvos1_start komvos1_end
numel = [1 2;
]; #orismos ton stixion os grammes me toys akraious komvous
N = size(numel)[1]; #arithmos stoixion
dofs = 3*nodes;
#orismos elastikotitas stixion
a = 20e-2;
E = 220e9.*ones(N,1); I = a^4 .*ones(N,1)./12; A = a^2 .*ones(N,1) #is square with axa
rho = 7.8e3;

#orismos ton desmeumenon komvon
BC = zeros(dofs,1)
#ignored = 1:3:dofs;#ignore turn of axis=ignore freedom
BC[1:3] .= 1; #BC[ignored].=1;
fixedofs = findall(x->x==1,BC[:])
freedofs = findall(x->x==0,BC[:])
#boundary conditions
U = zeros(dofs+3*N,1);
F = zeros(dofs+3*N,1); F[[4 5]] .= 10e3
#euresi olon ton deikton pou tha prospelasoun ton K
alldofs = zeros(N,6);
#orismos olikou mitroou stivarotitas
global K = sparse([],[],Float64[],dofs+3*N,dofs+3*N);
global M = sparse([],[],Float64[],dofs+3*N,dofs+3*N);
ignored = []

#build hermite interpolation
function Lagrange(x,xi,xj)
           p = poly(xj[:])
           p = p/p(xi)
           y = p(x)
           return y
end

function H(x,xi,xj)
    Lval(y) = Lagrange(y,xi,xj)
    Ldot = ForwardDiff.derivative(Lval,xi)
    H = (1-2*(x-xi)*Ldot)*Lagrange(x,xi,xj)^2
    return H
end
function Hbar(x,xi,xj)
    Hbar = (x-xi)*Lagrange(x,xi,xj)^2
    return Hbar
end
#define integration functions
function gauss(n)
    w = zeros(n,1)
    global coeffs = zeros(n+1,1);
    for k = Int64.(0:floor(n/2))
        coeffs[2*k+1] = (-1)^k*binomial(n,k)*binomial(2*n-2*k,n)/(2^n);
    end
    coeffs = Poly(reverse(coeffs,dims=1)[:])
    points = roots(coeffs);

    for i=1:n
        L = poly(points[[1:i-1; i+1:end]]);
        L = L/polyval(L,points[i]);
        L = polyint(L);
        w[i] = polyval(L,1)-polyval(L,-1);

    end
    #I = sum(w(i).*f(points));
    return w,points
end
gaussw,gausspoints = gauss(6) #used to integrate by gaussintegral

function gaussintegral(x1,x2,w,points,f)
    g(y) = 0.5*(x2-x1)*f.((x2-x1)*y/2 +(x2+x1)/2)
    I = sum(w.*g.(points))
    return I
end
#global K = zeros(dofs,dofs);
Nt(x) = [Lagrange(x,0,[0.5*L L]);
         Lagrange(x,0.5*L,[0 L]);
         Lagrange(x,L,[0 0.5*L])]
Ntd(x) = ForwardDiff.derivative(Nt,x)

Nb(x) = [H(x,0,[0.5*L L]);
         Hbar(x,0,[0.5*L L]);
         H(x,0.5*L,[0 L]);
         Hbar(x,0.5*L,[0 L]);
         H(x,L,[0 0.5*L]);
         Hbar(x,L,[0 0.5*L]) ]
Nbd(x) = ForwardDiff.derivative(Nb,x)
smolKt(x) = Ntd(x)*Ntd(x)'; smolKb(x) = Nb(x)*Nb(x)';
smolMt(x) = Nt(x)*Nt(x)'; smolMb(x) = Nbd(x)*Nbd(x)';

#start main loop to calculate stiffness matrix
for i=1:N
    i1 = numel[i,1]
    i2 = numel[i,2]
    #calclulate directions
    x1 = xy[1,i1]; x2 = xy[1,i2]; dx = x2-x1;
    y1 = xy[2,i1]; y2 = xy[2,i2]; dy = y2-y1;
    L = sqrt(dx^2+dy^2);
    nx = dx/L; ny = dy/L;
    #direction matrix
    c = nx; s = ny;
    Bbend = [-s c 0  0 0 0  0 0 0;
              0 0 1  0 0 0  0 0 0;
              0 0 0 -s c 0  0 0 0;
              0 0 0  0 0 1  0 0 0;
              0 0 0  0 0 0 -s c 0;
              0 0 0  0 0 0  0 0 1];
    Btens = [c s 0 0 0 0;
             0 0 c s 0 0;
             0 0 0 0 c s];
    kbend = E[i]*I[i]/L^3;

    #insert tensile
    #=
    ktens = E[i]*A[i]/L
    global Ketens = (ktens/3)*[ 7 -8  1;
                               -8 16 -8;
                                1 -8  7];
    Ketens = Btens'*Ketens*Btens
    global Metens = (rho*A[i]*L/15)*[2  1 -1/2;
                                     1  8  1;
                                    -1/2  1  2];
    Metens = Btens'*Metens*Btens;

    Kebend = (kbend/35/(12+13*L)^2).*[4*(25776+55560*L+30695*L^2) 2*L*(80592+173240*L+95025*L^2)   -512*(927+2040*L+1085*L^2) 512*L*(537+1175*L+665*L^2) -4*(56016+127800*L+75985*L^2) 2*L*(18384+42680*L+25585*L^2);
                     0  L^2*4*(118080+25400*L+13835*L^2) -256*L*(492+1040*L+555*L^2) 128*L^2*(354+770*L+435*L^2) -2*L*(17616+40120*L+23985*L^2) 2*L^2*(2832+6600*L+3995*L^2);
                     0  0     8192*(126+270*L+145*L^2)  1024*L*(12+40*L+25*L^2) -512*(1044+2280*L+1235*L^2) 256*L*(516+1120+605*L^2);
                     0  0     0   1024*L^2*(183+405*L+230*L^2) -512*L*(561+1255*L+715*L^2) 128*L^2*(378+850*L+485*L^2);
                     0  0     0     0   4*(189648+419640*L+234065*L^2) -2*L*(84432+186040*L+103025*L^2);
                     0  0     0     0  0   L^2*4*(12192+26680*L+14635*L^2)];
    Kebend = Symmetric(Kebend);
    Kebend = Bbend'*Kebend*Bbend;
    Mebend = rho*A[i]*(L/(12+13*L)^2)*[(337248+750276*L+417883*L^2)/15015   L*(57552+132764*L+76767*L^2)/45045  8*(18408+41602*L+23423*L^2)/30030     -2*L*(31848+60578*L+27535*L^2)/45045 (81168+181304*L+10061*L^2)/30030 -L*(48720+104128*L+54703*L^2)/45045;
                                       0   L^2*(4236+10167*L+6131*L^2)/45045 8*L*(5772+13573*L+7907*L^2)/45045 -2*L^2*(764+1381*L+556*L^2)/15015 L*(58704+139712*L+81263*L^2)/180180 -L^2*(1840+4208*L+2343*L^2)/60060;
                                       0     0   1024*(858+1846+993*L^2)/15015    128*L*(156+556*L+415*L^2)/45045 8*(15912+32706*L+16783*L^2)/15015 -8*L*(4524+9125*L+4587*L^2)/45045;
                                       0     0    0     128*L^2*(471+1029*L+568*L^2)/45045 2*L*(42456+98782*L+56745*L^2)/45045 -2*L^2*(1076+2515*L+1441*L^2)/15015;
                                       0     0    0    0     (316032+673868*L+359463*L^2)/15015 -L*(49440+103456*L+54197*L^2)/45045;
                                       0     0    0    0    0    L^2*(300+615*L+316*L^2)/4095];
    Mebend = Symmetric(Mebend);
    Mebend = Bbend'*Mebend*Bbend;
    =#
    ktens = E[i]*A[i]
    kbend = E[i]*I[i]

    Ketens = ktens*gaussintegral(0,L,gaussw,gausspoints,smolKt)
    Ketens = Btens'*Ketens*Btens
    Metens = rho[i]*A[i]*gaussintegral(0,L,gaussw,gausspoints,smolMt)
    Metens = Btens'*Metens*Btens

    Kebend = kbend.*gaussintegral(0,L,gaussw,gausspoints,smolKb)
    Kebend = Bbend'*Kebend*Bbend
    Mebend = rho[i]*A[i].*gaussintegral(0,L,gaussw,gausspoints,smolMb)
    Mebend = Bbend'*Mebend*Bbend

    tensiledofs = [1 2 4 5 7 8];
    global Ke = Kebend .+0;
    Ke[tensiledofs,tensiledofs] = Ke[tensiledofs,tensiledofs][1,:,1,:]  + Ketens;
        #println(size(Ke),"\n")
    ending = dofs+3*N
    #edofs = [3*i1-2 3*i1-1 3*i1 ending-3*i+3 ending-3*i+2 ending-3*i+1 3*i2-2 3*i2-1 3*i2]
    edofs = [1 2 3  7 8 9  4 5 6]
    push!(ignored,ending-3*i+3, ending-3*i+2, ending-3*i+1)
    K[edofs,edofs] = K[edofs,edofs][1,:,1,:] + Ke;
    #println(K,"\n")
    Me = Mebend .+0;
    Me[tensiledofs,tensiledofs] = Me[tensiledofs,tensiledofs][1,:,1,:]  + Metens;
    M[edofs,edofs] = M[edofs,edofs][1,:,1,:] + Me;
end
Kfinal = K[vcat(freedofs,ignored),vcat(freedofs,ignored)]
Mfinal = M[vcat(freedofs,ignored),vcat(freedofs,ignored)]
rhs = F
rhs = rhs[vcat(freedofs,ignored)]
rhs = rhs - K[vcat(freedofs,ignored),fixedofs]*U[fixedofs]

Ufree = Kfinal\rhs
U[freedofs] = Ufree[freedofs]
#Ufr = reshape(U,(3,nodes))
show(U[:])
#println(varinfo())
f = inv(Matrix(Mfinal))*Matrix(Kfinal)
omega = sqrt.(abs.(eigvals(f)))/(2*pi)
st = eigvecs(f)
plot(omega)

#calculate

max_stars_repo_path: axnai julia/kampsi-2d-high-degree.jl

max_stars_repo_name: billyon/first

max_issues_repo_path: axnai julia/kampsi-2d-high-degree.jl

max_issues_repo_name: billyon/first

max_forks_repo_path: axnai julia/kampsi-2d-high-degree.jl

max_forks_repo_name: billyon/first


text:
export have_method,invoke_callback, @rule, @inline_rule

classify_bool(seq, pred) = begin
    true_elems = eltype(seq)[]
    false_elems = eltype(seq)[]

    for elem in seq
        if pred(elem)
            push!(true_elems,elem)
        else
            push!(false_elems,elem)
        end
    end
    return true_elems, false_elems
end

"""
classify(sequence, key=nothing,value=nothing)

Return a dictionary where the items in sequence
are classified into groups according to the key()
function, with corresponding values calculated
using the value() function. If neither key nor
value are provided, items become their own
values.

This is e.g. used to classify patterns into types,
either literal strings or regular expressions.

"""
classify(seq;key=nothing,value=nothing) = begin
    d = Dict()
    for item in seq
        if !(key==nothing) k = key(item) else k = item end
        if !(value==nothing) v = value(item) else v = item end
        if haskey(d,k)
            push!(d[k],v)
        else
            d[k] = [v]
        end
    end
    return d
end

"""
    bfs(initial,expand)

Call function `expand` on elements of `initial` and
any newly generated elements until
no new elements are generated. A sequence of elements
are returned via a channel.
"""
#TODO: use proper Queue from DataStructures
bfs(initial,expand) = Channel() do q_chan
    # Python: open_q = deque(list(initial))
    open_q = collect(initial)
    visited = Set(open_q)
    # println("Queue start: $visited")
    while length(open_q) > 0
        node = popfirst!(open_q)
        #println("\nReturning $node")
        put!(q_chan,node)
        for next_node in expand(node)
            if !(next_node in visited)
                push!(visited,next_node)
                push!(open_q,next_node)
            end
        end
    end
    #println("BFS has finished!!")
end

# I know there's probably a fancy way to do this in Julia but this
# is basically a copy of the Python version for now
raw"""
    Accepts a list of alternatives, and enumerates all their possible concatinations.

    Examples:
        >>> combine_alternatives([range(2), [4,5]])
        [[0, 4], [0, 5], [1, 4], [1, 5]]

        >>> combine_alternatives(["abc", "xy", '$'])
        [['a', 'x', '$'], ['a', 'y', '$'], ['b', 'x', '$'], ['b', 'y', '$'], ['c', 'x', '$'], ['c', 'y', '$']]

        >>> combine_alternatives([])
        [[]]
"""
combine_alternatives(lists) = begin
    if isempty(lists)
        return [[]]
    end
    @assert all(l->!isempty(l),lists)
    init = [[x] for x in lists[1]]
    return reduce((x,y) -> [vcat(i,j) for i in x for j in y],ll[2:end],init=init)
end

# The Python version
# returns (min, max) possible width of matching strings. We only
# have access to the minimum width through PCRE, so we use that.
# Note that this relies on implementation details of the pcre2
# access in Julia base (not part of the published API) so may
# fail in future.  We do a naive guess for the maximum length.
get_regexp_width(regexp) = begin
    minlength = 0
    try
        r = Regex(regexp)
        minlength = Base.PCRE.info(r.regex,Base.PCRE.INFO_MINLENGTH,Int32)
    catch
        throw(error("Cannot parse $regexp as regexp"))
    end
    # A naive guess. Character ranges are not caught, neither are repeats
    if (occursin("+",regexp) && !occursin("\\+",regexp))|| (
        occursin("*",regexp) && !occursin("\\*",regexp))
        maxlength = 1000
        #println("$regexp has a big maxlength")
    else
        maxlength = minlength
    end
    return (minlength,maxlength)

end

## ============Only in the Julia version ===========
##
## Regular expression escape; copied from Python version

escape_re_string(s::String) = begin
    special_chars = "()[]{}?*+-|^\$\\.&~# \t\n\r\v\f"
    replace(s, Set(special_chars) => x -> "\\" * x)
end

## Find the largest index of matching captures

lastmatch(m::RegexMatch) = begin
    found_match = 0
    for i in 1:length(m.captures)
        if m.captures[i] != nothing
            found_match = i
        end
    end
    return found_match
end

## Partial function

partial(f,a...) = (b...) -> return f(a...,b...)

## Defaultdict: initialise with a Type that produces the
## default value when the key is missing. Mimic use of defaultdict
## in Lark

struct DefaultDict{K,V} <: AbstractDict{K,V}
    d::Dict{K,V}
end

DefaultDict{K,V}() where {K,V} = begin
    DefaultDict(Dict{K,V}())
end

Base.getindex(d::DefaultDict{K,V},k) where {K,V} = begin
    if !haskey(d,k)
        d.d[k] = V()
    end
    return d.d[k]
end

Base.setindex!(d::DefaultDict,v,k) = d.d[k] = v
Base.haskey(d::DefaultDict,k) = haskey(d.d,k)
Base.keys(d::DefaultDict) = keys(d.d)
Base.values(d::DefaultDict) = values(d.d)
Base.iterate(d::DefaultDict) = Base.iterate(d.d)
Base.iterate(d::DefaultDict,s) = Base.iterate(d.d,s)
Base.length(d::DefaultDict) = Base.length(d.d)

"""
    @rule s

`s` is a function `rule_name(t,args)` where `t` is an instance of a subtype of
`Transformer` or `Visitor`. `args` is an array holding values for each node of
the grammar production.
"""
macro rule(s)
    if s.head != :(=) || s.args[1].head != :call
        error("A rule must be a function definition")
    end
    rule_name = QuoteNode(s.args[1].args[1])
    if s.args[1].args[2].head != :(::)
        error("Type must be included in the first argument to define a rule: $s")
    end
    rule_type = s.args[1].args[2].args[2] # the type name
    @debug "Rule name: $rule_name, Rule type $rule_type"
    esc(quote
        Lerche.transformer_func($(s.args[1].args[2]),::Val{$rule_name},meta::Lerche.Meta,$(s.args[1].args[3])) = $(s.args[2])
    end )
end

"""
    @inline_rule s

`s` is a function of the form `rule_name(t,args...)` where `rule_name` is the
name of a grammar rule and `args` are values for the individual elements of that
rule's grammar production. This rule is called by `transform(t,tree)`
and `visit(v,tree)`. `t` is a subtype of `Transformer` and v is a
subtype of `Visitor`.
"""
macro inline_rule(s)
    if s.head != :(=) || s.args[1].head != :call
        error("A rule must be a function definition")
    end
    rule_name = QuoteNode(s.args[1].args[1])
    rule_type = s.args[1].args[2].args[2] # the type name
    @debug "Inline rule name: $rule_name, Rule type $rule_type"
    @debug "Args $(s.args[1].args)"
    esc(quote
        Lerche.transformer_func(x::$rule_type,y::Val{$rule_name},meta::Lerche.Meta,z::Array) = Lerche.transformer_func(x,y,z...)
        Lerche.transformer_func($(s.args[1].args[2]),::Val{$rule_name},$(s.args[1].args[3:end]...)) = $(s.args[2])
    end)
end

"""
    grab!(s::Stack,n)

Pop! the top `n` elements of `s`, returning them as a list in reverse order of removal
"""
grab!(s::Stack{T},n) where T = begin
    output = Array{T,1}(undef,n)
    for i in 1:n
        output[n-i+1] = pop!(s)
    end
    return output
end


text:
using DataFrames, Query, NamedTuples, TypedTables

df1 = DataFrame(a=[1,2,3], b=[1.,2.,3.])
df2 = @Table(c=[2.,4.,2.], d=["John", "Jim","Sally"])

x = @from i in df1 begin
    @join j in df2 on get(i.a) equals convert(Int,j.c) into k
    @where i.a>1
    @select {t1=i,t2=k}
    @collect DataFrame
end

println(x)


text:
using ArgParse2
using Test

function Base.redirect_stdout(f::Function, buf::Base.GenericIOBuffer)
    old_stdout = stdout
    try
        rd, = redirect_stdout()
        ret = f()
        Libc.flush_cstdio()
        flush(stdout)
        write(buf, readavailable(rd))
        return ret
    finally
        redirect_stdout(old_stdout)
    end
end

@testset "Smoke Test" begin
    parser = ArgumentParser(prog = "PROG", description="Smoke Test", epilog="Just testing smoke")
    add_argument! = argument_adder(parser)

    add_argument!("-f", "--foo")
    add_argument!("bar")

    @test parse_args(parser, ["BAR"]) === (bar = "BAR", foo = nothing)
    @test parse_args(parser, ["BAR", "--foo", "FOO"]) === (bar = "BAR", foo = "FOO")
    @test parse_args(parser, ["BAR", "--foo=FOO"]) === (bar = "BAR", foo = "FOO")
    @test_throws ArgumentError parse_args(parser, ["--foo", "FOO"])

    io = IOBuffer()

    redirect_stdout(io) do
        show_help(parser, exit_when_done=false)
        show_usage(parser)
    end
end

@testset "Action" begin
    @testset "Store1" begin
        parser = ArgumentParser()
        add_argument! = argument_adder(parser)

        add_argument!("--foo")

        @test parse_args(parser, split("--foo 1")) === (foo = "1",)
    end

    @testset "Store2" begin
        parser = ArgumentParser()
        add_argument! = argument_adder(parser)
        add_argument!("bar", nargs=2, type=Char)
        add_argument!("--foo", nargs=2)

        @test parse_args(parser, split("a b --foo d e")) == (bar = ['a', 'b'], foo = ["d", "e"])
    end

    @testset "Store Constant" begin
        parser = ArgumentParser()
        add_argument! = argument_adder(parser)

        add_argument!("--foo", action = "store_const", constant = 42)

        @test parse_args(parser, ["--foo"]) === (foo = 42,)
    end

    @testset "Store true/false" begin
        parser = ArgumentParser()
        add_argument! = argument_adder(parser)

        add_argument!("--foo", action = "store_true")
        add_argument!("--bar", action = "store_false")
        add_argument!("--baz", action = "store_false")

        @test parse_args(parser, split("--foo --bar")) === (foo = true, bar = false, baz = true)
    end

    @testset "Append" begin
        parser = ArgumentParser()
        add_argument! = argument_adder(parser)

        add_argument!("--foo", action = "append", default=["0"])

        @test parse_args(parser, split("--foo 1 --foo 2")) == (foo = ["0", "1", "2"],)
    end

    @testset "Append Constant" begin
        parser = ArgumentParser()
        add_argument! = argument_adder(parser)

        add_argument!("--str", dest = "types", action = "append_const", constant = String)
        add_argument!("--int", dest = "types", action = "append_const", constant = Int)

        @test parse_args(parser, split("--str --int")) == (types = [String, Int],)
    end

    @testset "Count" begin
        parser = ArgumentParser()
        add_argument! = argument_adder(parser)

        add_argument!("--verbose", "-v", action = "count", default = 0)

        @test parse_args(parser, ["-vvv"]) === (verbose = 3,)
    end

    @testset "Types" begin
        parser = ArgumentParser()
        add_argument! = argument_adder(parser)

        add_argument!("foo", type=Char)

        @test parse_args(parser, ["a"]) === (foo='a',)
        @test_throws ArgumentError parse_args(parser, ["abc"])
    end
end

@testset "Help" begin
    @testset "No help" begin
        parser = ArgumentParser(add_help = false)
        io = IOBuffer()
        show_help(io, parser, exit_when_done=false)
        output = String(take!(io)) |> rstrip

        @test output == "Usage: PROGRAM"
    end

    @testset "With help" begin
        parser = ArgumentParser()
        io = IOBuffer()
        show_help(io, parser, exit_when_done=false)
        output = String(take!(io)) |> rstrip |> x -> split(x, "\n") |> last |> lstrip

        @test output == "-h, --help   $(ArgParse2.HELP_TEXT)"
    end

    @testset "No Prog" begin
        parser = ArgumentParser(description="Test missing prog")
        add_argument!(parser, "foo")
        io = IOBuffer()
        show_help(io, parser, exit_when_done=false)
        output = String(take!(io)) |> x -> split(x, "\n") |> first

        @test output == "Usage: PROGRAM [-h] foo"
    end

    @testset "Prog" begin
        parser = ArgumentParser(prog="PROG", add_help = false)
        io = IOBuffer()
        show_help(io, parser, exit_when_done=false)
        output = String(take!(io)) |> rstrip

        @test output == "Usage: PROG"
    end

    @testset "Description" begin
        parser = ArgumentParser(add_help = false,
            description = """Frob the knob in the best way possible.
            Frobbing is an important step to take before fooing a bar,
            and is used in many fields when studying isoceles triangles
            in the wild.""")
        io = IOBuffer()
        show_help(io, parser, exit_when_done=false)
        output = String(take!(io)) |> rstrip

        @test output == ("""
        Usage: PROGRAM

        Frob the knob in the best way possible. Frobbing is an important step
        to take before fooing a bar, and is used in many fields when studying
        isoceles triangles in the wild.
        """ |> rstrip)
    end

    @testset "Usage" begin
        parser = ArgumentParser(usage = "PROG [options]", add_help=false)
        io = IOBuffer()
        show_usage(io, parser)
        output = String(take!(io)) |> rstrip

        @test output == "Usage: PROG [options]"
    end

    @testset "Multiple Arguments" begin
        parser = ArgumentParser()
        add_argument! = argument_adder(parser)
        add_argument!("foo", nargs=2)
        add_argument!("-b", nargs="?")
        add_argument!("-a", nargs="*")
        add_argument!("-r", nargs="+")

        io = IOBuffer()
        show_help(io, parser, exit_when_done=false)
        output=String(take!(io)) |> rstrip

        @test output == rstrip("""
        Usage: PROGRAM [-h] [-b [B]] [-a [A [A ...]]] [-r R [R ...]] foo foo

        Positional arguments:
         foo          foo help

        Optional arguments:
         -h, --help   show this help message and exit
         -b [B]       b help
         -a [A [A ...]]
                      a help
         -r R [R ...]
                      r help
        """)
    end
end

@testset "Arguments" begin
    @testset "Negative numerical arguments (1)" begin
        parser = ArgumentParser(prog="PROG")
        add_argument! = argument_adder(parser)

        add_argument!("-x")
        add_argument!("foo", nargs='?')

        # no negative number options, so -1 is a positional argument
        @test parse_args(parser, ["-x", "-1"]) === (foo=nothing, x="-1")

        # no negative number options, so -1 and -5 are positional arguments
        @test parse_args(parser, ["-x", "-1", "-5"]) === (foo="-5", x="-1")

    end

    @testset "Negative numerical arguments (1)" begin
        parser = ArgumentParser(prog="PROG")
        add_argument! = argument_adder(parser)
        add_argument!("-1", dest="one")
        add_argument!("foo", nargs='?')

        # negative number options present, so -1 is an option
        @test parse_args(parser, ["-1", "X"]) === (foo=nothing, one="X")

        # negative number options present, so -2 is an option
        @test_throws KeyError parse_args(parser, ["-2"])

        # negative number options present, so both -1s are options
        @test_throws ArgumentError parse_args(parser, ["-1", "-1"])
    end

    @testset "Multiple Arguments" begin
        @testset "?" begin
            parser = ArgumentParser()
            add_argument! = argument_adder(parser)
            add_argument!("bar", nargs="?")
            add_argument!("-f", nargs="?")

            @test parse_args(parser, []) === (bar=nothing, f=nothing)
            @test parse_args(parser, ["-f"]) === (bar=nothing, f=nothing)
            @test parse_args(parser, ["-f", "-1"]) === (bar=nothing, f="-1")
            @test parse_args(parser, ["-f", "-1", "-2"]) === (bar="-2", f="-1")
            @test parse_args(parser, ["a", "-f", "-1"]) === (bar="a", f="-1")
        end

        @testset "*" begin
            parser = ArgumentParser()
            add_argument! = argument_adder(parser)
            add_argument!("bar", nargs="*")
            add_argument!("-f", nargs="*")

            # Python's behavior
            @test_broken parse_args(parser, []) == (bar=nothing, f=nothing)
            @test_broken parse_args(parser, ["-f"]) == (bar=nothing, f=nothing)
            @test_broken parse_args(parser, ["-f", "-1"]) == (bar=nothing, f=["-1"])
            @test_broken parse_args(parser, ["-f", "-1", "-2"]) == (bar=nothing, f=["-1", "-2"])

            @test parse_args(parser, []) == (bar=[], f=[])
            @test parse_args(parser, ["-f"]) == (bar=[], f=[])
            @test parse_args(parser, ["-f", "-1"]) == (bar=[], f=["-1"])
            @test parse_args(parser, ["-f", "-1", "-2"]) == (bar=[], f=["-1", "-2"])
            @test parse_args(parser, ["a", "-f", "-1", "-2"]) == (bar=["a"], f=["-1", "-2"])
            @test parse_args(parser, ["a", "b", "c", "-f", "-1", "-2"]) == (bar=["a", "b", "c"], f=["-1", "-2"])
        end

        @testset "+ (1)" begin
            parser = ArgumentParser()
            add_argument! = argument_adder(parser)
            add_argument!("-f", "--foo", nargs="+")

            # Python's behavior
            @test_broken parse_args(parser, ["a"]) == (foo=nothing)

            @test parse_args(parser, []) == (foo=[],)
            @test_throws ArgumentError parse_args(parser, ["-f"])
        end

        @testset "+ (2)" begin
            parser = ArgumentParser()
            add_argument! = argument_adder(parser)
            add_argument!("bar", nargs="+")
            add_argument!("-f", "--foo", nargs="+")

            @test_throws ArgumentError parse_args(parser, [])
            # Python's behavior
            @test_broken parse_args(parser, ["a"]) == (bar=["a"], foo=nothing)

            @test parse_args(parser, ["a"]) == (bar=["a"], foo=[])

            @test_throws ArgumentError parse_args(parser, ["a", "-f"])

            @test parse_args(parser, ["a", "-f", "-1"]) == (bar=["a"], foo=["-1"],)
            @test parse_args(parser, ["a", "b", "c", "-f", "-1", "-2"]) == (bar=["a", "b", "c"], foo=["-1", "-2"],)
        end

        @testset "?*" begin
            parser = ArgumentParser()
            add_argument! = argument_adder(parser)
            add_argument!("foo", nargs="?")
            add_argument!("bar", nargs="*")

            @test parse_args(parser, []) == (foo=nothing, bar=[])
            @test parse_args(parser, ["a"]) == (foo="a", bar=[])
            @test parse_args(parser, ["a", "b", "c"]) == (foo="a", bar=["b", "c"])
        end
    end

    @testset "Choices" begin
        @testset "Array" begin
            parser = ArgumentParser()
            add_argument! = argument_adder(parser)
            add_argument!("move", choices=["rock", "paper", "scissors"])

            @test parse_args(parser, ["rock"]) === (move="rock",)
            @test_throws ArgumentError parse_args(parser, ["fire"])
        end

        @testset "Range" begin
            parser = ArgumentParser()
            add_argument! = argument_adder(parser)
            add_argument!("count", choices=1:3)

            @test parse_args(parser, ["1"]) === (count=1,)
            @test_throws ArgumentError parse_args(parser, ["4"])
        end
    end

    @testset "Dest1" begin
        parser = ArgumentParser()
        add_argument! = argument_adder(parser)
        add_argument!("-a", action="append", dest="dest", default=["a"])
        add_argument!("-b", action="append", dest="dest", default=["a"])
        @test parse_args(parser, []) == (dest=["a"],)
        @test parse_args(parser, split("-a c -b d")) == (dest=["c", "d"],)
    end

    @testset "Dest2" begin
        parser = ArgumentParser()
        add_argument! = argument_adder(parser)
        add_argument!("-a", action="append", dest="dest", default=["a"])
        add_argument!("-b", action="append", dest="dest", default=["b"])
        @test_throws ArgumentError parse_args(parser, [])
    end

    @testset "Required Options" begin
        parser = ArgumentParser()
        add_argument! = argument_adder(parser)
        add_argument!("-f", required=true)

        @test_throws ArgumentError parse_args(parser, [])
        @test parse_args(parser, ["-f", "foo"]) === (f="foo",)
    end
end

@testset "Construction Errors" begin
    parser = ArgumentParser()
    add_argument! = argument_adder(parser)
    @test_throws ArgumentError add_argument!("---foo")
    @test_throws ArgumentError add_argument!("-f", nargs="=")

    @test_throws ArgumentError add_argument!("-t", action="store_true", nargs=1)
    @test_throws ArgumentError add_argument!("-f", action="store_false", nargs=1)
    @test_throws ArgumentError add_argument!("-c", action="store_const", nargs=1, constant=1)
    @test_throws ArgumentError add_argument!("-a", action="append_const", nargs=1, constant=1)
    @test_throws ArgumentError add_argument!("-c", action="count", nargs=1)

    @test_throws ArgumentError add_argument!("-c", action="store_const")
    @test_throws ArgumentError add_argument!("-c", action="append_const")

    @test_throws ArgumentError add_argument!("-a", nargs="*", default="aaa")
    @test_throws ArgumentError add_argument!("-c", choices=["a","b"], default="c")
    @test_throws ArgumentError add_argument!("-c", choices=["a","b"], default=["c"])

    @test_throws ArgumentError add_argument!("-a", action="nothing")
end

@testset "Argument/flag Formatting" begin
    @test ArgParse2.format_flag("--flag", 1, "FLAG") === "--flag FLAG"
    @test ArgParse2.format_flag("--flag", 2, "FLAG") === "--flag FLAG FLAG"
    @test ArgParse2.format_flag("--flag", :?, "FLAG") === "--flag [FLAG]"
    @test ArgParse2.format_flag("--flag", :*, "FLAG") === "--flag [FLAG [FLAG ...]]"
    @test ArgParse2.format_flag("--flag", :+, "FLAG") === "--flag FLAG [FLAG ...]"

    @test_throws ArgumentError ArgParse2.format_arg_name("name", :(=))
end

@testset "Misc" begin
    @test ArgParse2.to_symbol(1) === 1
    for action in [:store_true, :store_false, :store_const, :append_const, :count]
        @test_throws ArgumentError ArgParse2.validate_args(
            action, nothing, nothing, nothing, nothing, ["true", "false"], false
        )
    end
    @test ArgParse2.validate_action("hi") === nothing
    @test ArgParse2.validate_nargs(1) === nothing
    @test ArgParse2.get_nargs(nothing, Bool, :store_true) === 0
    @test ArgParse2.get_nargs(nothing, Vector{String}, nothing) === :*
    @test ArgParse2.get_nargs(1, String, nothing) === 1
    @test ArgParse2.get_nargs(:*, String, nothing) === :*
    @test ArgParse2.get_nargs(:+, String, nothing) === :+
    @test ArgParse2.parse_item(String, "somestring", :arg) == "somestring"
    @test_throws ArgumentError ArgParse2.process_zero_arg_flag(nothing, :bad_action, nothing, nothing, nothing, nothing)

    parser = ArgumentParser(description="Test missing prog")
    @test_logs (:warn, r"add_argument\(...\) is deprecated.*") add_argument(parser, "foo")
    @test parse_args(parser, ["bar"]) === (foo = "bar",)
end

nothing


text:
const AnalysisData = ImmutableDict{Type{<:AbstractAnalysis}, Any}

# parametrize metadata by M
mutable struct EClass
    id::Int64
    nodes::Vector{ENode}
    parents::Vector{Pair{ENode, Int64}}
    data::Union{Nothing, AnalysisData}
    # data::M
end
EClass(id) = EClass(id, ENode[], Pair{ENode, Int64}[], nothing)
EClass(id, nodes, parents) = EClass(id, nodes, parents, nothing)

# Interface for indexing EClass
Base.getindex(a::EClass, i) = a.nodes[i]
Base.setindex!(a::EClass, v, i) = setindex!(a.nodes, v, i)
Base.firstindex(a::EClass) = firstindex(a.nodes)
Base.lastindex(a::EClass) = lastindex(a.nodes)

# Interface for iterating EClass
Base.iterate(a::EClass) = iterate(a.nodes)
Base.iterate(a::EClass, state) = iterate(a.nodes, state)

# Showing
function Base.show(io::IO, a::EClass)
    print(io, "EClass $(a.id) (")
    print(io, collect(a.nodes))
    if a.data === nothing
        print(io, ")")
        return
    end
    print(io, ", analysis = {")
    for (k, v) ∈ a.data
        print(io, "$k => $v, ")
    end
    print(io, "})")
end
#
# function addparent!(a::EClass, n::ENode, p::EClass)
#     a.parents[n] = p
# end

function addparent!(a::EClass, n::ENode, id::Int64)
    # if (n => id) ∉ a.parents
        push!(a.parents, (n => id))
    # end
    # a.parents[n] = id
end

function Base.union(to::EClass, from::EClass)
    EClass(to.id, vcat(from.nodes, to.nodes),
        vcat(from.parents, to.parents),
        if to.data !== nothing && from.data !== nothing
            join_analysis_data(to.data, from.data)
        elseif to.data === nothing
            from.data
        else nothing end)
end

function Base.union!(to::EClass, from::EClass)
    append!(to.nodes, from.nodes)
    append!(to.parents, from.parents)
    if to.data !== nothing && from.data !== nothing
        # merge!(to.data, from.data)
        # to.data = join_analysis_data(to.data, from.data)
        to.data = join_analysis_data(to.data, from.data)
    elseif to.data === nothing
        to.data = from.data
    end
    return to
end

function join_analysis_data(d::AnalysisData, dsrc::AnalysisData)
    for (an, val_b) in dsrc
        if haskey(d, an)
            val_a = d[an]
            nv = join(an, val_a, val_b)
            # d[an] = nv
            # WARNING immutable version
            d = Base.ImmutableDict(d,an=>nv)
        end
    end
    return d
end

# Thanks to Shashi Gowda
function hasdata(a::EClass, x::Type{<:AbstractAnalysis})
    a.data === nothing && (return false)
    return haskey(a.data, x)
end

function getdata(a::EClass, x::Type{<:AbstractAnalysis})
    !hasdata(a, x) && error("EClass $a does not contain analysis data for $x")
    return a.data[x]
end

function getdata(a::EClass, x::Type{<:AbstractAnalysis}, default)
    hasdata(a, x) ? a.data[x] : default
end

function setdata!(a::EClass, x::Type{<:AbstractAnalysis}, value)
    # lazy allocation
    a.data === nothing && (a.data = AnalysisData())
    # a.data[x] = value
    a.data = AnalysisData(a.data, x, value)
end

function funs(a::EClass)
    map(a.nodes) do x
        x.head
    end
end

function funs_arity(a::EClass)
    map(a.nodes) do x
        (x.head, arity(x))
    end
end


text:
## A manipulate like widget
## this follows the interface of RStudio's manipulate, which is likely following that of
## Mathematica's

## This is *not* the manipulate of Interact.jl!

using Winston ## Winston *must* load before PySide, this is flaky
using PySide

## We render output into a label. Here are some different renderers...
function render(parent, p::String)
    setText(parent, p)
end
function render(parent::Label, p::Winston.FramedPlot)
    nm = tempname() * ".png"
    file(p, nm)
    setPixmap(parent, Qt[:QPixmap](nm))
end

function render(parent::Label, p::Winston.FramedPlot)
    nm = tempname() * ".png"
    file(p, nm)
    setPixmap(parent, Qt[:QPixmap](nm))
end

## context to store dynamic values
module ManipulateContext
using Winston
end

abstract ManipulateWidget

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

type SliderWidget <: ManipulateWidget
    w
    nm
    label
    initial
    rng
end

slider(nm::String, label::String, rng::Range1, initial::Integer) = SliderWidget(nothing, nm, label, initial, rng)
slider(nm::String, label::String, rng::Range1) = slider(nm, label, rng, minimum(rng))
slider(nm::String,  rng::Union(Range, Range1)) = slider(nm, nm, rng, minimum(rng))

function make_widget(parent, widget::SliderWidget)
    sl = Slider(parent, "Horizontal", widget.rng)
    set_value(sl, widget.initial)
    widget.w = sl
end

type PickerWidget <: ManipulateWidget
    w
    nm
    label
    initial
    vals
end

function make_widget(parent, widget::PickerWidget)
    cb = ComboBox(parent)
    set_items(cb, widget.vals)
    set_value(cb, widget.initial)
    widget.w = cb
end

picker{T <: String}(nm::String, label::String, vals::Vector{T}, initial) = PickerWidget(nothing, nm, label, initial, vals)
picker{T <: String}(nm::String, label::String, vals::Vector{T}) = picker(nm, label, vals, vals[1])
picker{T <: String}(nm::String, vals::Vector{T}) = picker(nm, nm, vals)
picker(nm::String, label::String, vals::Dict, initial) = PickerWidget(nm, label, vals, initial)
picker(nm::String, label::String, vals::Dict) = PickerWidget(nm, label, vals, [string(k) for (k,v) in vals][1])
picker(nm::String, vals::Dict) = picker(nm, nm, vals)

type CheckboxWidget <: ManipulateWidget
    w
    nm
    label
    initial
end
function make_widget(parent, widget::CheckboxWidget)
    cb = CheckBox(parent)
    set_items(cb, widget.label)
    set_value(cb, widget.initial)
    widget.w = cb
end

checkbox(nm::String, label::String, initial::Bool) = CheckboxWidget(nothing, nm, label, initial)
checkbox(nm::String, label::String) = checkbox(nm, label, false)

type ButtonWidget <: ManipulateWidget
    w
    label
    nm
end
function make_widget(parent, widget::ButtonWidget)
    btn = Button(parent)
    set_value(btn, widget.label)
    widget.w = btn
end

button(label::String) = ButtonWidget(nothing, label, nothing)

## Add text widget to gather one-line of text
type EntryWidget <: ManipulateWidget
    w
    nm
    label
    initial
end
function make_widget(parent, widget::EntryWidget)
    e = LineEdit(parent)
    set_value(e, widget.initial)
    widget.w = e
end

entry(nm::String, label::String, initial::String) = EntryWidget(nothing, nm, label, initial)
entry(nm::String, initial::String) = EntryWidget(nm, nm, initial)
entry(nm::String) = EntryWidget(nm, nm, "")

## Expression returns a plot object. Use names as values
function manipulate(ex::Union(Symbol,Expr), w::Widget, controls...)

    lyt = HBoxLayout(w)
    setLayout(w, lyt)

    pg = Splitter(w)
    addWidget(lyt, pg)

    control_pane = Widget(pg)
    graph =Label(pg)               # This would be QtSvg.QSvgWidget if using Gadfly
    addWidget(pg, control_pane)
    addWidget(pg, graph)

    form_lyt = FormLayout(control_pane)
    setLayout(control_pane, form_lyt)

    ## create, layout widgets
    for widget in controls
        make_widget(control_pane, widget)
        addRow(form_lyt, widget.label, widget.w)
    end

    get_values() = [get_value(i.w) for i in  controls]
    get_nms() = map(u -> u.nm, controls)
    function get_vals()
        d = Dict()                      # return Dict of values
        vals = get_values(); keys = get_nms()
        for i in 1:length(vals)
            if !isa(keys[i], Nothing)
                d[keys[i]] = vals[i]
            end
        end
        d
    end

    function dict_to_module(d::Dict) ## stuff values into Manipulate Context
        for (k,v) in d
            eval(ManipulateContext, :($(symbol(k)) = $v))
        end
    end

    function make_graphic(x...)
        d = get_vals()
        dict_to_module(d)
        p = eval(ManipulateContext, ex)
        render(graph, p)
    end

    map(u -> change_slot(u.w, make_graphic), controls)
    make_graphic()
    controls
end

### How this is employed:
ex = quote
    x = linspace( 0, n * pi, 100 )
    c = cos(x)
    s = sin(x)
    p = FramedPlot()
    setattr(p, "title", title)
    if fillbetween
        add(p, FillBetween(x, c, x, s) )
    end
    add(p, Curve(x, c, "color", color) )
    add(p, Curve(x, s, "color", "blue") )
    p
end

w = Widget()
obj = manipulate(ex, w,
                 slider("n", "[0, n*pi]", 1:10)
                 ,entry("title", "Title", "title")
                 ,checkbox("fillbetween", "Fill between?", true)
                 ,picker("color", "Cos color", ["red", "green", "yellow"])
                 ,button("update")
                 )
raise(w)


text:
# Autogenerated wrapper script for ffts_jll for aarch64-linux-musl
export libffts

JLLWrappers.@generate_wrapper_header("ffts")
JLLWrappers.@declare_library_product(libffts, "libffts.so.0.9.0")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libffts,
        "lib/libffts.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using ITensors, ITensorsGPU

function my_polar( T::ITensor, inds...; kwargs... )
    U, S, V, spec, u, v = svd(T, inds; kwargs...)
    replaceindex!(U, u, v)
    return U*V
end

function initQs( A::PEPS, col::Int, next_col::Int; kwargs...)
    Ny, Nx = size(A)
    maxdim::Int = get(kwargs, :maxdim, 1)
    Q         = MPO(Ny, deepcopy(A[:, col]), 0, Ny+1)
    prev_col  = next_col > col ? col - 1 : col + 1
    A_r_inds  = [commonindex(A[row, col], A[row, next_col]) for row in 1:Ny]
    QR_inds   = [Index(ITensors.dim(A_r_inds[row]), "Site,QR,c$col,r$row") for row in 1:Ny]
    A_up_inds = [commonindex(A[row, col], A[row+1, col]) for row in 1:Ny-1]
    Q_up_inds = [Index(ITensors.dim(A_up_inds[row]), "Link,u,Qup$row") for row in 1:Ny-1]
    next_col_inds = [commonindex(A[row, col], A[row, next_col]) for row in 1:Ny]
    prev_col_inds = 0 < prev_col < Nx ? [commonindex(A[row, col], A[row, prev_col]) for row in 1:Ny] : Vector{Index}(undef, Ny)
    for row in 1:Ny
        row < Ny && replaceindex!(Q[row], A_up_inds[row], Q_up_inds[row])
        row > 1  && replaceindex!(Q[row], A_up_inds[row-1], Q_up_inds[row-1])
        replaceindex!(Q[row], next_col_inds[row], QR_inds[row])
    end
    return Q, QR_inds, next_col_inds
end

function gaugeQR(A::PEPS, col::Int, side::Symbol; kwargs...)
    overlap_cutoff::Real = get(kwargs, :overlap_cutoff, 1e-4)
    max_gauge_iter::Int  = get(kwargs, :max_gauge_iter, 100)
    Ny, Nx = size(A)
    is_gpu = !(data(store(A[1,1])) isa Array)
    prev_col_inds = Vector{Index}(undef, Ny)
    next_col = side == :left ? col - 1 : col + 1
    prev_col = side == :left ? col + 1 : col - 1
    Q, QR_inds, next_col_inds = initQs(A, col, next_col; kwargs...)
    left_edge     = col == 1
    right_edge    = col == Nx
    ratio_history = Vector{Float64}()
    ratio         = 0.0
    best_overlap  = 0.0
    Ampo          = MPO(Ny)
    best_Q        = MPO(Ny)
    best_R        = MPO(Ny)
    iter          = 1
    dummy_nexts   = [Index(ITensors.dim(next_col_inds[row]), "DM,Site,r$row") for row in 1:Ny]
    cmb_l         = Vector{ITensor}(undef, Ny)
    iter          = 0
    while best_overlap < overlap_cutoff
        thisTerm  = Vector{ITensor}(undef, Ny)
        thisfTerm = Vector{ITensor}(undef, Ny)
        @timeit "build envs" begin
            for row in 1:Ny
                Ap = dag(deepcopy(A[row, col]))'
                Ap = setprime(Ap, 0, next_col_inds[row]')
                Qp = dag(deepcopy(Q[row]))'
                Qp = setprime(Qp, 0, QR_inds[row]')
                thisTerm[row]  = A[row, col] * Ap * Qp
                thisfTerm[row] = thisTerm[row] * Q[row]
            end
            Envs = thisTerm
            fF   = cumprod(thisfTerm)
            rF   = reverse(cumprod(reverse(thisfTerm)))
            for row in 1:Ny
                if row > 1
                    Envs[row] *= fF[row - 1]
                end
                if row < Ny
                    Envs[row] *= rF[row + 1]
                end
            end
        end
        @timeit "polar decomp" begin
            for row in 1:Ny
                if row < Ny
                    Q_ = my_polar(Envs[row], QR_inds[row], commonindex(Q[row], Q[row+1]); kwargs...)
                    Q[row] = deepcopy(noprime(Q_))
                else
                    Q_ = my_polar(Envs[row], QR_inds[row]; kwargs...)
                    Q[row] = deepcopy(noprime(Q_))
                end
                AQinds     = IndexSet(findindex(A[row, col], "Site"))
                if (side == :left && !right_edge) || (side == :right && !left_edge)
                    AQinds = IndexSet(AQinds, commonindex(findinds(A[row, col], "Link"), findinds(Q[row], "Link"))) # prev_col_ind
                end
                cmb_l[row] = combiner(AQinds, tags="Site,AQ,r$row")[1]
                Ampo[row]  = A[row, col] * cmb_l[row]
                Ampo[row]  = replaceindex!(Ampo[row], next_col_inds[row], dummy_nexts[row])
                Q[row]    *= cmb_l[row]
            end
        end
        @timeit "Q*A -> R" begin
            R           = multMPO(dag(Q), Ampo; kwargs...)
        end
        @timeit "compute overlap" begin
            aqr_overlap = is_gpu ? cuITensor(1.0) : ITensor(1.0)
            a_norm      = is_gpu ? cuITensor(1.0) : ITensor(1.0)
            cis         = [commonindex(Ampo[row], Ampo[row+1]) for row in 1:Ny-1]
            Ampo_       = deepcopy(Ampo)
            for row in 1:Ny
                if row < Ny
                    Ampo_[row] = prime(Ampo_[row], cis[row])
                end
                if row > 1
                    Ampo_[row] = prime(Ampo_[row], cis[row-1])
                end
                aqr_overlap *= Ampo[row] * Q[row] * R[row]
                Q[row]      *= cmb_l[row]
                a_norm      *= dag(Ampo_[row]) * Ampo[row]
            end
            ratio = abs(scalar(collect(aqr_overlap)))/abs(scalar(collect(a_norm)))
        end
        push!(ratio_history, ratio)
        if ratio > best_overlap || iter == 0
            best_Q = deepcopy(Q)
            best_R = deepcopy(R)
            best_overlap = ratio
        end
        ratio > overlap_cutoff && break
        iter += 1
        iter > max_gauge_iter && break
        if (iter > 10 && best_overlap < 0.5) || (iter > 20 && mod(iter, 20) == 0)
            Q_, QR_inds_, next_col_inds_ = initQs(A, col, next_col; kwargs...)
            for row in 1:Ny
                if row < Ny
                    old_u = commonindex(Q[row], Q[row+1])
                    new_u = commonindex(Q_[row], Q_[row+1])
                    replaceindex!(Q_[row], new_u, old_u)
                    replaceindex!(Q_[row+1], new_u, old_u)
                end
                replaceindex!(Q_[row], QR_inds_[row], QR_inds[row])
            end
            for row in 1:Ny
                Q[row]  = Q_[row]
                salt    = is_gpu ? cuITensor(randomITensor(inds(Q[row])))/100.0 : randomITensor(inds(Q[row]))/100.0
                salt_d  = ratio < 0.5 ? norm(salt) : 10.0*norm(salt)
                Q[row] += salt/salt_d
                Q[row] /= sqrt(norm(Q[row]))
            end
        end
    end
    @info "best overlap: ", best_overlap
    return best_Q, best_R, next_col_inds, QR_inds, dummy_nexts
end

function gaugeColumn( A::PEPS, col::Int, side::Symbol; kwargs...)
    Ny, Nx = size(A)

    prev_col_inds = Vector{Index}(undef, Ny)
    next_col_inds = Vector{Index}(undef, Ny)

    next_col   = side == :left ? col - 1 : col + 1
    prev_col   = side == :left ? col + 1 : col - 1
    left_edge  = col == 1
    right_edge = col == Nx
    is_gpu = !(data(store(A[1,1])) isa Array)

    @timeit "gauge QR" begin
        Q, R, next_col_inds, QR_inds, dummy_next_inds = gaugeQR(A, col, side; kwargs...)
    end
    @timeit "update next col" begin
        cmb_r = Vector{ITensor}(undef, Ny)
        cmb_u = Vector{ITensor}(undef, Ny - 1)
        if (side == :left && col > 1 ) || (side == :right && col < Nx)
            next_col_As = MPO(Ny, deepcopy(A[:, next_col]), 0, Ny+1)
            nn_col = side == :left ? next_col - 1 : next_col + 1
            cmb_inds_ = [IndexSet(findindex(A[row, next_col], "Site")) for row in 1:Ny]
            for row in 1:Ny
                cmb_inds_r = cmb_inds_[row]
                if 0 < nn_col < Nx + 1
                    cmb_inds_r = IndexSet(cmb_inds_[row], commonindex(A[row, next_col], A[row, nn_col]))
                end
                cmb_r[row] = combiner(cmb_inds_r, tags="Site,CMB,c$col,r$row")[1]
                next_col_As[row] *= cmb_r[row]
                if (side == :left && !left_edge) || (side == :right && !right_edge)
                    next_col_As[row] = replaceindex!(next_col_As[row], next_col_inds[row], dummy_next_inds[row])
                end
            end
        end
        maxdim::Int  = get(kwargs, :maxdim, 1)
        result       = multMPO(R, next_col_As; kwargs...)
        true_QR_inds = [Index(ITensors.dim(QR_inds[row]), "Link,r,r$row" * (side == :left ? ",c$(col-1)" : ",c$col")) for row in 1:Ny]
        for row in 1:Ny
            A[row, col] = Q[row]
        end
        cUs = [commonindex(A[row, col], A[row+1, col]) for row in 1:Ny-1]
        true_U_inds = [Index(ITensors.dim(cUs[row]), "Link,u,r$row,c$col") for row in 1:Ny-1]
        A[:, col] = [replaceindex!(A[row, col], QR_inds[row], true_QR_inds[row]) for row in 1:Ny]
        A[:, col] = vcat([replaceindex!(A[row, col], cUs[row], true_U_inds[row]) for row in 1:Ny-1], A[Ny, col])
        A[:, col] = vcat(A[1, col], [replaceindex!(A[row, col], cUs[row-1], true_U_inds[row-1]) for row in 2:Ny])

        for row in 1:Ny
            A[row, next_col] = result[row]
        end
        A[:, next_col] = [replaceindex!(A[row, next_col], QR_inds[row], true_QR_inds[row]) for row in 1:Ny]
        cUs            = [commonindex(A[row, next_col], A[row+1, next_col]) for row in 1:Ny-1]
        true_nU_inds   = [Index(ITensors.dim(cUs[row]), "Link,u,r$row,c" * string(next_col)) for row in 1:Ny-1]
        A[:, next_col] = vcat([replaceindex!(A[row, next_col], cUs[row], true_nU_inds[row]) for row in 1:Ny-1], A[Ny, next_col])
        A[:, next_col] = vcat(A[1, next_col], [replaceindex!(A[row, next_col], cUs[row-1], true_nU_inds[row-1]) for row in 2:Ny])
        A[:, next_col] = [A[row, next_col] * cmb_r[row] for row in 1:Ny]
    end
    return A
end

function gaugeColumnForInsert( A::PEPS, col::Int, side::Symbol; kwargs...)
    Ny, Nx = size(A)

    prev_col_inds = Vector{Index}(undef, Ny)
    next_col_inds = Vector{Index}(undef, Ny)

    next_col   = side == :left ? col - 1 : col + 1
    prev_col   = side == :left ? col + 1 : col - 1
    left_edge  = col == 1
    right_edge = col == Nx
    is_gpu = !(data(store(A[1,1])) isa Array)

    @timeit "gauge QR" begin
        Q, R, next_col_inds, QR_inds, dummy_next_inds = gaugeQR(A, col, side; kwargs...)
    end
    return A, Q, R, dummy_next_inds
end


text:
σ(x::Real) = 1.0 / (1.0 + exp(-x))

abstract type LossFunction end

struct Logloss <: LossFunction end

function (logloss::Logloss)(logit, label)
    prob = σ(logit)
    if label == 1
        return -log(prob)
    elseif label == 0
        return -log(1 - prob)
    else
        throw(ArgumentError("label must be 0 or 1."))
    end
end

function grad(loss::Logloss, logit, label)
    return σ(logit) - label
end

struct BPRLoss <: LossFunction end

function (loss::BPRLoss)(pred)
    return -log(σ(pred))
end

function grad(loss::BPRLoss, pred)
    return -1 + σ(pred)
end


text:
using HttpServer, WebSockets, JSON

timestamp(t = time()) = Libc.strftime("%Y-%m-%d %H:%M.%S", t)

function start_http_and_websocket_servers(httpcontent, httpport = 8081, wsport = 8083)
    http = HttpHandler() do req::Request, res::Response
        @show req
        Response( servedcontent )
    end
    httpserver = Server( http )
    println("Starting http server on port $httpport")
    @async run(httpserver, httpport)

    wsh = WebSocketHandler() do req,client
        println("WebSocket opened from $client:\n $req")

        # Send meta-info
        starttime = time()
        ts = timestamp(starttime)
        metainfo = Dict(:info => "meta", :StartTime => ts)
        write(client, JSON.json(metainfo))

        # For sending optimization events we do not set the info key. This saves some time
        # in the json conversion on both sides. We reuse one and the same Dict to save on mem.
        obj = Dict(:ElapsedSeconds => 0, :BestFitness => 1000 + rand()*1000, :NumEvals => 1)

        while true
            sleep(0.1 + rand() * 0.3)
            obj[:ElapsedSeconds] = t = time() - starttime
            obj[:BestFitness] -= (rand() * 100 / (3 + log10(100 + 0.1 * t)))
            obj[:NumEvals] += 1
            try
                write(client, JSON.json(obj))
            catch err
                @show err
                break
            end
        end
    end

    wsserver = Server(wsh)
    println("Starting web socket server on port $wsport")
    @async run(wsserver, wsport)

    return httpserver, wsserver
end

servedcontent = read(open("vega_lite_updating_from_websocket.html"), String)
start_http_and_websocket_servers(servedcontent, 8081, 8083)
sleep(120.0)


text:
# Parse arguments, if any
path = "RUNME.jl"
for arg in ARGS
  if ispath(arg)
    path = arg
  else
    error("Unknown argument: $arg")
  end
end

using Samantha
using Random
using Test

if path == "RUNME.jl"
  @info "Running All Tests"
else
  @info "Running Test $path"
end
include(path)


text:
" Linear solver that uses the built-in linear algebra capabilities "
type JuliaDenseImplicitSolver <: DenseImplicitSolver
    " The System the solve will be done for "
    system::System

    " The 'stiffness matrix' "
    mat::Matrix{Float64}

    " The RHS (force vector) "
    rhs::Vector{Float64}

    " The solution vector "
    solution::Vector{Float64}

    " The ghosted solution vector (for Assembly) "
    ghosted_solution::Vector{Float64}

    " Whether or not this Solver has been initialized "
    initialized::Bool

    JuliaDenseImplicitSolver(system::System) = new(system, Matrix{Float64}(), Vector{Float64}(), Vector{Float64}(), Vector{Float64}(), false)
end

"""
    Solve using the built-in dense matrix/vector types

    `assemble` controls whether or not the system will be assembled automatically
"""
function solve!(solver::JuliaDenseImplicitSolver; assemble=true)
    if !solver.initialized
        initialize!(solver)
    end

    if assemble
        assembleResidualAndJacobian(solver, solver.system)
    end

    # Direct solve
    solver.solution = \(solver.mat, -solver.rhs)
end


text:
using Brown
using Test

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