text:
mutable struct Position
    x::Int
    y::Int
end

struct Size
    width::Int
    height::Int
end

s1 = Spaceship(Position(0,0), Size(30,5), Missile);
s2 = Spaceship(Position(10,0), Size(30,5), Laser);
a1 = Asteroid(Position(20,0), Size(20,20));
a2 = Asteroid(Position(0,20), Size(20,20));

struct Rectangle
    top::Int
    left::Int
    bottom::Int
    right::Int
    # return the upper-left and lower-right points
    Rectangle(p::Position, s::Size) =
        new(p.y + s.height, p.x, p.y, p.x+s.width)
end

# check if the 2 rectangles (A & B) overlap
function overlap(A::Rectangle, B::Rectangle)
    return A.left < B.right && A.right > B.left &&
           A.top > B.bottom && A.bottom < B.top
end

function collide(A::Thing, B::Thing)
    println("Checking collision of thing vs. thing")
    rectA = Rectangle(position(A), size(A))
    rectB = Rectangle(position(B), size(B))
    return overlap(rectA, rectB)
end

function collide(A::Spaceship, B::Spaceship)
    println("Checking collision of spaceship vs. spaceship")
    return true # just a test
end

# Randomly pick two things and check
function check_randomly(things)
    for i in 1:5
        two = rand(things, 2)
        collide(two...)
    end
end


text:
module GeometricIntegratorsDiffEq

using Reexport
@reexport using DiffEqBase

using GeometricIntegrators

const warnkeywords =
    (:save_idxs, :d_discontinuities, :unstable_check, :save_everystep,
     :save_end, :initialize_save, :adaptive, :abstol, :reltol, :dtmax,
     :dtmin, :force_dtmin, :internalnorm, :gamma, :beta1, :beta2,
     :qmax, :qmin, :qsteady_min, :qsteady_max, :qoldinit, :failfactor,
     :maxiters, :isoutofdomain, :unstable_check,
     :calck, :progress, :tstops, :saveat, :dense)

function __init__()
    global warnlist = Set(warnkeywords)
end

include("algorithms.jl")
include("solve.jl")

export GeometricIntegratorAlgorithm, GIEuler, GIMidpoint, GIHeun2, GIHeun3,
       GIRalston2, GIRalston3, GIRunge, GIKutta, GIRK4, GIRK416, GIRK438, GISSPRK3,
       GICrankNicolson, GIKraaijevangerSpijker, GIQinZhang, GICrouzeix,
       GIImplicitEuler, GIImplicitMidpoint, GISRK3,
       GIGLRK, GILobattoIIIA, GILobattoIIIB, GILobattoIIIC, GILobattoIIIC̄,
       GILobattoIIID, GILobattoIIIE, GILobattoIIIF, GIRadauIA, GIRadauIIA,
       GISymplecticEulerA, GISymplecticEulerB, GILobattoIIIAIIIB2, GILobattoIIIBIIIA2

end # module

ext: jl


text:
#= basic support for sparse matrices, more Magma style than Julia
  a sparse matrix is a collection (Array) of sparse rows (SRow)
  A sparse matrix cannot have 0 rows!
  A SRow is two Arrays, one (pos) containing the columns, values
  contains the values.
  A[i,j] = A.rows[i].values[A.rows[i].pos[j]]
  to be formal

  The upper_triangular stuff performs very elementary transformations
  until the matrix becomes dense. At this point, the dense bit is extraced and
  converted to an fmpz_mat which is then hnf'ed.

  Missing:
   full HNF, Howell, modular and not
   better elemination strategy

  TODO: sort out the various upper_triangular implementations.
  One (+one with trafo) is probably enough
=#

import Base.push!, Base.max, Nemo.nbits, Base.sparse, Base.Array,
       Base.endof, Base.start, Base.done, Base.next, Base.hcat,
       Base.vcat, Base.max, Base.min

export upper_triangular, vcat!, show, sub, SMat, SRow, random_SMatSLP,
       fmpz_mat, rows, cols, copy, push!, mul, mul!, toNemo, sparse,
       valence_mc, swap_rows!, endof, start, done, next, elementary_divisors,
       randrow, hcat, hcat!, vcat, vcat!, mod!, mod_sym!

function SLP_AddRow(i::Int, j::Int, v::T) where T
  assert(v != 0)
  slp = TrafoAddScaled(i, j, v)
  return slp
end

function SLP_SwapRows(i::Int, j::Int)
  slp = TrafoSwap{fmpz}(i, j)
  return slp
end

# function mul!{T}(a::Array{T, 1}, s::SMatSLP_swap_row)
#   t = a[s.row]
#   a[s.row] = a[s.col]
#   a[s.col] = t
# end
#
# function mul!{T}(a::Array{T, 1}, s::SMatSLP_add_row{T})
#   a[s.row] = a[s.row]*s.val + a[s.col]
# end
#
# function mul_t!{T}(a::Array{T, 1}, s::SMatSLP_swap_row)
#   t = a[s.row]
#   a[s.row] = a[s.col]
#   a[s.col] = t
# end
#
# function mul_t!{T}(a::Array{T, 1}, s::SMatSLP_add_row{T})
#   a[s.col] = a[s.col]*s.val + a[s.row]
# end
#
#
# function apply!{T}(a::Array{T}, b::Array{SMatSLP_add_row{T}, 1})
#   for i=length(b):-1:1
#     mul!(a, b[i])
#   end
#   return a
# end
#
# function apply_t!{T}(a::Array{T}, b::Array{SMatSLP_add_row{T}, 1})
#   for i=1:length(b)
#     mul_t!(a, b[i])
#   end
#   return a
# end
#
function random_SMatSLP(A::SMat{T}, i::Int, v::UnitRange) where T
  a = Array{SMatSLP_add_row{Int}}(i)
  for j=1:i
    c = rand(v)
    while c==0
      c = rand(v)
    end
    i1 = rand(1:rows(A))
    i2 = rand(1:rows(A))
    while i1==i2
      i1 = rand(1:rows(A))
      i2 = rand(1:rows(A))
    end
    a[j] = SLP_AddRow(i1, i2, Int(c))
  end
  return a
end


text:
using MacroTools
using IRTools

using IRTools: branches, block, empty!, evalir, func, branch!, block, IR, @dynamo, xcall

function _mark(label, ex)
    label isa Symbol || error("label has to be a Symbol")
    return Expr(
        :block,
        Expr(:meta, :begin_optional, label),
        esc(ex),
        Expr(:meta, :end_optional, label),
    )
end

macro mark(label, ex)
    _mark(label, ex)
end

foo(x) = bar(baz(x))

function bar(x)
    @mark print iseven(x) && println("The input is even.")
    x
end

function baz(x)
    @mark print x<0 && println("The input is negative.")
    x
end

isbegin(e::Expr) = Meta.isexpr(e,:meta) && e.args[1]===:begin_optional
isend(e::Expr) = Meta.isexpr(e,:meta) && e.args[1]===:end_optional

skip(f::Core.IntrinsicFunction, args...) = f(args...)
skip(f::Core.Builtin, args...) = f(args...)

@dynamo function skip(args...)
    ir = IR(args...)
    delete_line = false
    local orig

    for (x,st) in ir
        is_begin = isbegin(st.expr)
        is_end   = isend(st.expr)

        if is_begin
            delete_line = true
        end

        if is_begin
            orig = block(ir,x)
        elseif is_end
            dest = block(ir,x)
            if orig != dest
                empty!(branches(orig))
                branch!(orig,dest)
            end
        end

        if delete_line
            delete!(ir,x)
        end

        if is_end
            delete_line = false
        end

        if haskey(ir,x) && Meta.isexpr(st.expr,:call)
            ir[x] = IRTools.xcall(skip, st.expr.args...)
        end
    end
    return ir
end

function skip(ex::Expr)
end

macro skip(ex)
    ex.head == :call || error("Input expression has to be a `:call`.")
    return xcall(skip, ex.args...)
end

display(@code_ir foo(-2))
display(@code_ir skip(foo,-2))
display(foo(-2))
@skip foo(-2)


text:
using KernelAbstractions
using Test

include("testsuite.jl")

backend = get(ENV, "KERNELABSTRACTIONS_TEST_BACKEND", "CPU")

if backend != "CPU"
    @info "CPU backend not selected"
    exit()
end

struct CPUDeviceArray{T,N,A} end # Fake and unused
Testsuite.testsuite(CPU, backend, Base, Array, CPUDeviceArray)


text:
module TestDataFrameRow
    using Compat, Compat.Test, DataFrames

    df = DataFrame(a=Union{Int, Missing}[1, 2, 3, 1, 2, 2],
                   b=[2.0, missing, 1.2, 2.0, missing, missing],
                   c=["A", "B", "C", "A", "B", missing],
                   d=CategoricalArray([:A, missing, :C, :A, missing, :C]))
    df2 = DataFrame(a = [1, 2, 3])

    #
    # Equality
    #
    @test_throws ArgumentError DataFrameRow(df, 1) == DataFrameRow(df2, 1)
    @test DataFrameRow(df, 1) != DataFrameRow(df, 2)
    @test DataFrameRow(df, 1) != DataFrameRow(df, 3)
    @test DataFrameRow(df, 1) == DataFrameRow(df, 4)
    @test DataFrameRow(df, 2) == DataFrameRow(df, 5)
    @test DataFrameRow(df, 2) != DataFrameRow(df, 6)

    # isless()
    df4 = DataFrame(a=[1, 1, 2, 2, 2, 2, missing, missing],
                    b=Union{Float64, Missing}[2.0, 3.0, 1.0, 2.0, 2.0, 2.0, 2.0, 3.0],
                    c=[:B, missing, :A, :C, :D, :D, :A, :A])
    @test isless(DataFrameRow(df4, 1), DataFrameRow(df4, 2))
    @test !isless(DataFrameRow(df4, 2), DataFrameRow(df4, 1))
    @test !isless(DataFrameRow(df4, 1), DataFrameRow(df4, 1))
    @test isless(DataFrameRow(df4, 1), DataFrameRow(df4, 3))
    @test !isless(DataFrameRow(df4, 3), DataFrameRow(df4, 1))
    @test isless(DataFrameRow(df4, 3), DataFrameRow(df4, 4))
    @test !isless(DataFrameRow(df4, 4), DataFrameRow(df4, 3))
    @test isless(DataFrameRow(df4, 4), DataFrameRow(df4, 5))
    @test !isless(DataFrameRow(df4, 5), DataFrameRow(df4, 4))
    @test !isless(DataFrameRow(df4, 6), DataFrameRow(df4, 5))
    @test !isless(DataFrameRow(df4, 5), DataFrameRow(df4, 6))
    @test isless(DataFrameRow(df4, 7), DataFrameRow(df4, 8))
    @test !isless(DataFrameRow(df4, 8), DataFrameRow(df4, 7))

    # hashing
    @test hash(DataFrameRow(df, 1)) != hash(DataFrameRow(df, 2))
    @test hash(DataFrameRow(df, 1)) != hash(DataFrameRow(df, 3))
    @test hash(DataFrameRow(df, 1)) == hash(DataFrameRow(df, 4))
    @test hash(DataFrameRow(df, 2)) == hash(DataFrameRow(df, 5))
    @test hash(DataFrameRow(df, 2)) != hash(DataFrameRow(df, 6))

    # check that hashrows() function generates the same hashes as DataFrameRow
    df_rowhashes, _ = DataFrames.hashrows(df, false)
    @test df_rowhashes == [hash(dr) for dr in eachrow(df)]

    # test incompatible frames
    @test_throws UndefVarError DataFrameRow(df, 1) == DataFrameRow(df3, 1)

    # test RowGroupDict
    N = 20
    d1 = rand(map(Int64, 1:2), N)
    df5 = DataFrame(Any[d1], [:d1])
    df6 = DataFrame(d1 = [2,3])

    # test_group("groupby")
    gd = DataFrames.group_rows(df5)
    @test gd.ngroups == 2

    # getting groups for the rows of the other frames
    @test length(gd[DataFrameRow(df6, 1)]) > 0
    @test_throws KeyError gd[DataFrameRow(df6, 2)]
    @test isempty(DataFrames.findrows(gd, df6, (gd.df[1],), (df6[1],), 2))

    # grouping empty frame
    gd = DataFrames.group_rows(DataFrame(x=Int[]))
    @test gd.ngroups == 0

    # grouping single row
    gd = DataFrames.group_rows(df5[1,:])
    @test gd.ngroups == 1
end


text:
function SimpRulT(A,B,Q,R,U,bet,n1,n2,F,SigmaXX,x10,cutoff=1.0)
#SimpRulT   Solving RE model for given simple rule u(t)  = -F * x(t),
#           where   x(t)  = [x1(t)' x2(t)']'.
#
#
#
#  Usage:    (M,C,J0) = SimpRulT(A,B,Q,R,U,bet,n1,n2,F,SigmaXX,x10,cutoff)
#
#  Input:     A,B,Q,R,U,bet,n1,n2: see ComItAlg
#             F          kxn matrix, decision rule
#             SigmaXX    n1xn1 covariance matrix of shocks
#             x10        n1 vector, initial values for predetermined
#                        variables. Not n1x1 array.
#             cutoff     scalar, max modulus of backward looking variables
#
#  Output:    M          n1xn1 matrix, x1(t+1) = M*x1(t) + e(t+1)
#             C          n2xn1 matrix, x2(t) = C*x1(t)
#             J0         scalar, value of loss function
#
#
#
#  Paul Soderlind, Paul.Soderlind@unisg.ch, Aug 2000, to Julia Jan 2016
#-----------------------------------------------------------------------

  Q = (Q + Q')/2                  #to make symmetric
  R = (R + R')/2

  A = A - B*F                     #subst for u(t) in evolution
  n = n1 + n2

  GS = schur(Matrix(1.0I,n,n),A)
  logconA = abs.(GS.beta) .<= (abs.(GS.alpha)*cutoff) #selecting stable eigenvalues
  ordschur!(GS,logconA)                                 #reordering, stable first

  if sum(logconA) < n1
    @warn("Too few stable roots: no stable solution")
    return NaN,NaN,NaN
  elseif sum(logconA) > n1
    @warn("Too many stable roots: inifite number of stable solutions")
    return NaN,NaN,NaN
  end

  Stt = GS.S[1:n1,1:n1]
  Zkt = GS.Z[1:n1,1:n1]
  Zlt = GS.Z[n1+1:n,1:n1]
  Ttt = GS.T[1:n1,1:n1]

  if cond(Zkt) > 1e+14
    @warn("Zkt is singular: rank condition for solution not satisfied")
    return NaN,NaN,NaN
  end

  Zkt_1 = inv(Zkt)         #inverting
  Stt_1 = inv(Stt)

  M = real(Zkt*Stt_1*Ttt*Zkt_1)       #x1(t+1) = M*x1(t) + e(t+1)
  C = real(Zlt*Zkt_1)                 #x2(t) = C*x1(t)

  P = [ I;
        C;
        (-F*[I;C]) ]
  QUUR = [ Q   U;
           U'  R ]
  PQUURP = P'QUUR*P

  vecV = (I - kron(M',bet*M'))\vec(PQUURP)
  V    = reshape(vecV,n1,n1)                 #solves V = W + bet*M'V*M

  J0 =  x10'V*x10 + (bet/(1-bet)) * tr(V*SigmaXX)

  return M,C,J0

end
#----------------------------------------------------------------------

max_stars_repo_name: PaulSoderlind/PaulSoderlindCode

max_issues_repo_name: PaulSoderlind/PaulSoderlindCode

max_forks_repo_name: PaulSoderlind/PaulSoderlindCode


text:
@testset "Default reward" begin
    @testset "set_reward!(DecisionPhase)" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        lh = SeaPearl.LearnedHeuristic(agent)
        SeaPearl.update_with_cpmodel!(lh, model)

        lh.reward.value = 0
        SeaPearl.set_reward!(SeaPearl.DecisionPhase, lh, model)
        @test lh.reward.value == -0.025f0
    end

    @testset "set_reward!(EndingPhase)" begin
        trailer = SeaPearl.Trailer()
        model = SeaPearl.CPModel(trailer)

        lh = SeaPearl.LearnedHeuristic(agent)
        SeaPearl.update_with_cpmodel!(lh, model)

        lh.reward.value = 5
        model.statistics.numberOfNodes = 30
        SeaPearl.set_reward!(SeaPearl.EndingPhase, lh, model, nothing)
        @test lh.reward.value == 6
    end
end


text:
using Colors, GLVisualize
using GLVisualize.ComposeBackend, Gadfly, DataFrames, RDatasets

if !isdefined(:runtests)
    window = glscreen()
    composebackend = ComposeBackend.GLVisualizeBackend(window)
end
const static_example = true

p = plot(dataset("car", "SLID"), x="Wages", color="Language", Geom.histogram)

draw(composebackend, p)

if !isdefined(:runtests)
	renderloop(window)
end


text:
# testing some special cases that can get missed in more complicated tests
# these are more checks for crashes than accuracy
# tests zeroth order mass action reactions and passing just one mass action reaction
# tests FW and tuple based mix of mass action jump / constant jump

using DiffEqBase, DiffEqJump
using Test

doprint = false
#using Plots; plotlyjs()
doplot = false

methods = (RDirect(), RSSACR(), Direct(), DirectFW(), FRM(), FRMFW(), SortingDirect(), NRM(), RSSA(), DirectCR())

# one reaction case, mass action jump, vector of data
rate = [2.0]
rs = [[0 => 1]]
ns = [[1 => 1]]
jump = MassActionJump(rate, rs, ns)
prob = DiscreteProblem([100],(0.,100.))
jump_prob = JumpProblem(prob, Direct(), jump)
sol = solve(jump_prob, SSAStepper())
if doprint
    println("mass act jump using vectors of data: last val = ", sol[end, end])
end
if doplot
    plothand = plot(sol, lab="vec, 0=>1")
end
@test sol[end,end] > 200

# one reaction case, mass action jump, data as scalars
rate = 2.0
rs = [0 => 3]                   # stoich power should be ignored
ns = [1 => 1]
jump = MassActionJump(rate, rs, ns)
jump_prob = JumpProblem(prob, Direct(), jump)
sol = solve(jump_prob, SSAStepper())
if doprint
    println("mass act jump using scalar data: last val = ", sol[end, end])
end
if doplot
    plot!(plothand, sol, lab="scalar, 0=>1")
end
@test sol[end,end] > 200

# one reaction case, mass action jump, vector of empty data
rate = [2.0]
rs = [Vector{Pair{Int,Int}}()]
ns = [[1 => 1]]
jump = MassActionJump(rate, rs, ns)
prob = DiscreteProblem([100],(0.,100.))
jump_prob = JumpProblem(prob, Direct(), jump)
sol = solve(jump_prob, SSAStepper())
if doprint
    println("mass act jump using vector of Pair{Int,Int}: last val = ", sol[end, end])
end
if doplot
    plot!(plothand, sol, lab="vec, ()")
end
@test sol[end,end] > 200

# one reaction case, mass action jump, data as scalars
rate = 2.0
rs = Vector{Pair{Int,Int}}()
ns = [1 => 1]
jump = MassActionJump(rate, rs, ns)
jump_prob = JumpProblem(prob, Direct(), jump)
sol = solve(jump_prob, SSAStepper())
if doprint
    println("mass act jump using scalar Pair{Int,Int}: last val = ", sol[end, end])
end
if doplot
    plot!(plothand, sol, lab="scalar, ()")
end
@test sol[end,end] > 200

# mix two rx types
# 0 -> A and A -> 0
rate = 600.0
rs = [0 => 3]                   # stoich power should be ignored
ns = [1 => 1]
jump = MassActionJump(rate, rs, ns)
ratefun = (u,p,t) -> 2.0*u[1]
affect! = function(integrator)
    integrator.u[1] -= 1
end
jump2 = ConstantRateJump(ratefun, affect!)
if doplot
    plothand2 = plot(reuse=false)
end

dep_graph = [
    [1, 2],
    [1, 2]
]
spec_to_dep_jumps = [[2]]
jump_to_dep_specs = [[1],[1]]
namedpars = (dep_graph=dep_graph, vartojumps_map=spec_to_dep_jumps, jumptovars_map=jump_to_dep_specs)

for method in methods
    jump_prob = JumpProblem(prob, method, jump, jump2; namedpars...)
    sol = solve(jump_prob, SSAStepper())

    if doplot
        plot!(plothand2, sol, label=("A <-> 0, " * string(method)))
    end

    if doprint
        println("Mix of constant and mass action jumps, method = ", typeof(method), ", sol[end] = ", sol[end,end])
    end
    @test sol[end,end] > 200
end

if doplot
    display(plothand)
    display(plothand2)
end


text:
get_char(::Nothing) = '⋅'
get_color(::Nothing) = :white

get_char(::Agent) = '☻'
get_char(::Agent, ::Up) = '↑'
get_char(::Agent, ::Down) = '↓'
get_char(::Agent, ::Left) = '←'
get_char(::Agent, ::Right) = '→'
get_color(::Agent) = :red

get_char(::Wall) = '█'
get_color(::Wall) = :white

get_char(::Goal) = '♥'
get_color(::Goal) = :red

get_char(::Door) = '⩎'
get_color(::Door) = :yellow

get_char(::Key) = '⚷'
get_color(::Key) = :yellow

get_char(::Gem) = '♦'
get_color(::Gem) = :magenta

get_char(::Obstacle) = '⊗'
get_color(::Obstacle) = :yellow

get_char(::Box) = '▒'
get_color(::Box) = :yellow

get_char(::Target) = '✖'
get_color(::Target) = :red

get_char(::Target1) = '✖'
get_color(::Target1) = :yellow

get_char(::Target2) = '✖'
get_color(::Target2) = :green

get_char(::Body) = '∘'
get_color(::Body) = :green

get_char(::Food) = '♦'
get_color(::Food) = :yellow

get_char(::Basket) = '⊔'
get_color(::Basket) = :red

get_char(::Ball) = '∘'
get_color(::Ball) = :yellow

#####
# GridWorldBase
#####

function get_first_object(grid::AbstractArray{Bool, 3}, objects, pos::CartesianIndex{2})
    idx = findfirst(grid[:, pos])
    if isnothing(idx)
        return nothing
    else
        return objects[idx]
    end
end

function get_render_data(world::GridWorldBase, pos::CartesianIndex{2})
    grid = get_grid(world)
    objects = get_objects(world)

    object = get_first_object(grid, objects, pos)
    background = :black
    foreground = get_color(object)
    char = get_char(object)

    return background, foreground, char
end

function get_render_data(world::GridWorldBase, pos::CartesianIndex{2}, layer)
    if world[layer, pos]
        object = get_objects(world)[layer]
        return :black, get_color(object), get_char(object)
    else
        return :black, get_color(nothing), get_char(nothing)
    end
end

function render(io::IO, ::MIME"text/plain", world)
    for i in 1:get_height(world)
        for j in 1:get_width(world)
            pos = CartesianIndex(i, j)
            background, foreground, char = get_render_data(world, pos)
            print(io, Crayons.Crayon(background = background, foreground = foreground, bold = true, reset = true), char)
        end
        println(io, Crayons.Crayon(reset = true))
    end
end

function render_layers(io::IO, ::MIME"text/plain", world)
    objects = get_objects(world)
    for (layer, object) in enumerate(objects)
        println("layer = $layer, object = $object")
        for i in 1:get_height(world)
            for j in 1:get_width(world)
                pos = CartesianIndex(i, j)
                background, foreground, char = get_render_data(world, pos, layer)
                print(io, Crayons.Crayon(background = background, foreground = foreground, bold = true, reset = true), char)
            end
            println(io, Crayons.Crayon(reset = true))
        end
    end
end

#####
# AbstractGridWorld
#####

function get_render_data(env::E, pos::CartesianIndex{2}) where {E<:AbstractGridWorld}
    grid = get_grid(env)
    objects = get_objects(env)

    object = get_first_object(grid, objects, pos)
    background = :black
    foreground = get_color(object)
    char = get_char(object)

    if hasfield(E, :agent_pos)
        agent_pos = get_agent_pos(env)
        if hasfield(E, :agent_dir)
            agent_dir = get_agent_dir(env)
            if isa(RLBase.StateStyle(env), RLBase.Observation)
                if pos in get_grid_inds(agent_pos.I, get_half_size(env), agent_dir)
                    background = :dark_gray
                end
            end
            if object === AGENT
                char = get_char(object, agent_dir)
            end
        else
            if isa(RLBase.StateStyle(env), RLBase.Observation)
                if pos in get_grid_inds(agent_pos.I, get_half_size(env))
                    background = :dark_gray
                end
            end
        end
    end

    return background, foreground, char
end

function get_render_data(env::E, pos::CartesianIndex{2}, layer) where {E<:AbstractGridWorld}
    world = get_world(env)
    grid = get_grid(world)

    object = nothing
    background = :black
    foreground = get_color(object)
    char = get_char(object)
    if world[layer, pos]
        object = get_objects(world)[layer]
        foreground = get_color(object)
        char = get_char(object)
    end

    if hasfield(E, :agent_pos)
        agent_pos = get_agent_pos(env)
        if hasfield(E, :agent_dir)
            agent_dir = get_agent_dir(env)
            if isa(RLBase.StateStyle(env), RLBase.Observation)
                if pos in get_grid_inds(agent_pos.I, get_half_size(env), agent_dir)
                    background = :dark_gray
                end
            end
            if object === AGENT
                char = get_char(object, agent_dir)
            end
        else
            if isa(RLBase.StateStyle(env), RLBase.Observation)
                if pos in get_grid_inds(agent_pos.I, get_half_size(env))
                    background = :dark_gray
                end
            end
        end
    end

    return background, foreground, char
end

get_render_data(env::Snake, pos::CartesianIndex{2}) = get_render_data(get_world(env), pos)
get_render_data(env::Snake, pos::CartesianIndex{2}, layer) = get_render_data(get_world(env), pos, layer)


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    imfilter_fft(img, kern)

Perform filtering on `img` with kernel `kern` using the FFT algorithm.
"""
imfilter_fft(img::AbstractArray{T,N}, kern::AbstractArray{K,N}) where {T<:Real,K<:Real,N} =
  imfilter(img, centered(kern), Inner(), Algorithm.FFT())

"""
    convdist(img, patch, weights)

Compute distance between all patches of `img` and a single `patch`
using `weights` for each pixel in the `patch`.
"""
function convdist(img::AbstractArray, patch::AbstractArray, weights::AbstractArray)
  wpatch = weights.*patch

  A² = imfilter_fft(img.^2, weights)
  AB = imfilter_fft(img, wpatch)
  B² = sum(wpatch .* patch)

  abs.(A² .- 2AB .+ B²)
end

"""
    selectpatch(img, patchsize, center)

Helper function to extract a patch of size
`patchsize` from `img` centered at Cartesian
index `center`.
"""
function selectpatch(img, patchsize, center)
  start  = CartesianIndex(@. center.I - (patchsize-1)÷2)
  finish = CartesianIndex(@. center.I + (patchsize  )÷2)
  view(img, start:finish)
end


text:
module Todo

using Redux

# actions
struct AddTodo{T<:AbstractString} <: AbstractSyncAction
    text::T
end

struct DeleteTodo <: AbstractSyncAction
    id::Int
end

struct EditTodo{T<:AbstractString} <: AbstractSyncAction
    id::Int
    text::T
end

struct CompleteTodo <: AbstractSyncAction
    id::Int
end

struct CompleteAllTodos <: AbstractSyncAction end

struct ClearCompleted <: AbstractSyncAction end

# state
struct State <: AbstractImmutableState
    id::Int
    completed::Bool
    text::String
end

const INITIAL_STATE = State[State(0, false, "Use Redux")]

# reducers
todo(state::AbstractState, action::AbstractAction) = state
todo(state::Vector{<:AbstractState}, action::AbstractAction) = state
todo(s::Vector{State}, a::AddTodo) = State[s..., State(isempty(s) ? 1 : s[end].id + 1, false, a.text)]
todo(s::Vector{State}, a::DeleteTodo) = filter(s -> s.id !== a.id, s)
todo(s::Vector{State}, a::EditTodo) = map(s) do s
    s.id === a.id ? State(s.id, s.completed, a.text) : s
end
todo(s::Vector{State}, a::CompleteTodo) = map(s) do s
    s.id === a.id ? State(s.id, !s.completed, s.text) : s
end
function todo(s::Vector{State}, a::CompleteAllTodos)
    are_all_marked = all(x->x.completed, s)
    map(x->State(x.id, !are_all_marked, x.text), s)
end
todo(s::Vector{State}, a::ClearCompleted) = filter(s -> s.completed == false, s)

end # module


text:
using Pkg, Random
Pkg.activate(".")

module EnzymeMFABC
include("../example.jl")
end

using FileIO, JLD2

Random.seed!(888)
BSL_stats = EnzymeMFABC.do_bsl.([2_500, 5_000, 10_000]; N_batch=500)

c = map(BSL_stats) do x
    x[3]
end
v = map(BSL_stats) do x
    x[2]
end

save("BSL_stats.jld2", Dict("c" => c, "v" => v))

max_stars_repo_name: tpprescott/mf-lf

max_issues_repo_name: tpprescott/mf-lf

max_forks_repo_name: tpprescott/mf-lf


text:
using ParticleMethods
using Test

using StaticArrays: MVector

@testset "ParticleMethods.jl" begin
    @testset "Particles" begin
        include("particle_tests.jl")
    end

    @testset "Particle Lists" begin
        include("particle_list_tests.jl")
    end
end


text:
using LinearAlgebra
using CompScienceMeshes
using SauterSchwab3D
using StaticArrays

const pI   = point(1,0,0)
const pII  = point(0,1,0)
const pIII = point(0,0,0)
const pIV  = point(0,0,1)

const qIV  = point(0,0,-1)

const P = simplex(pI,pII,pIV,pIII)
const Q = simplex(pI,pII,qIV,pIII)

const sing = SauterSchwab3D.singularity_detection(P,Q)

Accuracy2 = 18
cf_ref = CommonFace6D(sing,SauterSchwab3D._legendre(Accuracy2,0.0,1.0))

function integrand(x,y)
      return ((x-pI)'*(y-pIV))*exp(-im*1*norm(x-y))/(4pi*norm(x-y))
end

function INTEGRAND(u,v)
   j1 = volume(P) * factorial(dimension(P))
   j2 = volume(Q) * factorial(dimension(Q))
   x = barytocart(P,u)
   y = barytocart(Q,v)
   output = integrand(x,y)*j1*j2

   return output
end

print("Ref: ")
ref = sauterschwab_parameterized(INTEGRAND, cf_ref)
println(ref)
println()

res_tp =[]
res_sp =[]
res_gm =[]
n1 = []
n2 = []
n3 = []
for i in 2:1:15
   Accuracy = i
   cf = CommonFace6D(sing,SauterSchwab3D._legendre(Accuracy,0.0,1.0))

   int_tp = sauterschwab_parameterized(INTEGRAND, cf)

   num_pts = 15*length(cf.qps)^6
   push!(n1,num_pts)

   push!(res_tp,int_tp)

end

for i in 2:1:7
   Accuracy = i
   cf_s = CommonFace6D_S(sing,(SauterSchwab3D._legendre(Accuracy,0.0,1.0),
                          SauterSchwab3D._shunnham2D(Accuracy),
                          SauterSchwab3D._shunnham3D(Accuracy)))

   int_sp = sauterschwab_parameterized(INTEGRAND, cf_s)

   num_pts = 3*length(cf_s.qps[1])*length(cf_s.qps[2])*length(cf_s.qps[3])+
            12*length(cf_s.qps[3])^2
   push!(n2,num_pts)
   push!(res_sp,int_sp)
end

for i in 2:1:12
   Accuracy = i
   cf_gm = CommonFace6D_S(sing,(SauterSchwab3D._legendre(Accuracy,0.0,1.0),
                          SauterSchwab3D._grundmannMoeller2D(2*Accuracy-1),
                          SauterSchwab3D._grundmannMoeller3D(2*Accuracy-1)))

   int_gm = sauterschwab_parameterized(INTEGRAND, cf_gm)

   num_pts = 3*length(cf_gm.qps[1])*length(cf_gm.qps[2])*length(cf_gm.qps[3])+
            12*length(cf_gm.qps[3])^2
   push!(n3,num_pts)
   push!(res_gm,int_gm)
end

err_tp = norm.(res_tp.-ref)/norm(ref)
err_sp = norm.(res_sp.-ref)/norm(ref)
err_gm = norm.(res_gm.-ref)/norm(ref)

println(err_tp)
println(err_sp)
println(err_gm)

using Plots
plot( yaxis=:log, xaxis=:log, fontfamily="Times")
plot!(n1,err_tp, label="Gauss Tensor-Product",markershape=:circle)
plot!(n2,err_sp, label="Simplex Tensor-Product",markershape=:rect)
#plot!(n3,err_gm, label="Simplex-Product GM",markershape=:x)
plot!(xlims=(1e2,1e7),ylims=(1e-8,1))
plot!(xlabel="#Quad. pts/Func. evals", ylabel="Rel. Error.", title="Common Face 6D",legend=:bottomleft)

using BenchmarkTools

ref = 0.0011204365593010864 - 0.000870131729342847im

cf = CommonFace6D(sing,SauterSchwab3D._legendre(4,0.0,1.0))

int_tp = sauterschwab_parameterized(INTEGRAND, cf)

num_pts = 15*length(cf.qps)^6
println("#Pts: ",num_pts)
err_tp = norm.(int_tp.-ref)/norm(ref)
println("Rel. Err: ",err_tp)
@time for i in 1:100 sauterschwab_parameterized(INTEGRAND, cf) end

cf_s = CommonFace6D_S(sing,(SauterSchwab3D._legendre(5,0.0,1.0),
                         SauterSchwab3D._shunnham2D(5),
                         SauterSchwab3D._shunnham3D(5)))

int_sp = sauterschwab_parameterized(INTEGRAND, cf_s)

num_pts = 3*length(cf_s.qps[1])*length(cf_s.qps[2])*length(cf_s.qps[3])+
            12*length(cf_s.qps[3])^2
println("#Pts: ",num_pts)
err_sp = norm.(int_sp.-ref)/norm(ref)
println("Rel Err: ",err_sp)
@time for i in 1:100 sauterschwab_parameterized(INTEGRAND, cf_s) end
#=
cf_gm = CommonFace6D_S(sing,(SauterSchwab3D._legendre(5,0.0,1.0),
                         SauterSchwab3D._grundmannMoeller2D(2*5-1),
                         SauterSchwab3D._grundmannMoeller3D(2*5-1)))

int_gm = sauterschwab_parameterized(INTEGRAND, cf_gm)

num_pts = 3*length(cf_gm.qps[1])*length(cf_gm.qps[2])*length(cf_gm.qps[3])+
            12*length(cf_gm.qps[3])^2
println("#Pts: ",num_pts)
err_gm = norm.(int_gm.-ref)/norm(ref)
println("Rel Err: ",err_gm)

=#


text:
###
### LongSequence
###
###
### A general purpose biological sequence representation.
###
### It is cheap to create a subsequence from a sequence because sequences can
### share underlying data: creating a subsequence is just allocating a new
### sequence object defining a part of the original sequence without copy.
### Destructive operations create a copy of the underlying data if and only if it
### is shared between (sub)sequences. This is often called as copy-on-write
### strategy in computer science and should be transparent to the user.
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

# About Internals
# ---------------
#
# The `data` field of a `LongSequence{A}` object contains binary representation
# of a biological character sequence. Each character is encoded with an encoder
# corresponding to the alphabet `A` and compactly packed into `data`. To extract
# a character from a sequence, you should decode this binary sequence with a
# decoder that is a pair of the encoder. The length of encoded binary bits is
# fixed, and hence a character at arbitrary position can be extracted in a
# constant time. To know the exact location of a character at a position, you
# can use the `bitindex(seq, i)` function, which returns a pair of element's index
# containing binary bits and bits' offset. As a whole, character extraction
# `seq[i]` can be written as:
#
#     j = bitindex(seq, i)
#     decode(A, (seq.data[index(j)] >> offset(j)) & mask(A))
#
#  index :        index(j) - 1       index(j)       index(j) + 1
#   data : ....|xxxxx...........|xxXxxxxxxxxxxxxx|............xxxx|....
# offset :                          |<-offset(j)-|
#  width :      |<---- 64 ---->| |<---- 64 ---->| |<---- 64 ---->|
#
#  * '.' : unused (4 bits/char)
#  * 'x' : used
#  * 'X' : used and pointed by index `i`

"""
Biological sequence data structure indexed by an alphabet type `A`.
"""
mutable struct LongSequence{A <: Alphabet} <: BioSequence{A}
    data::Vector{UInt64}  # encoded character sequence data
    part::UnitRange{Int}  # interval within `data` defining the (sub)sequence
    shared::Bool          # true if and only if `data` is shared between sequences

    function LongSequence{A}(data::Vector{UInt64},
                            part::UnitRange{Int},
                            shared::Bool) where A
        return new(data, part, shared)
    end
end

const LongNucleotideSequence = LongSequence{<:NucleicAcidAlphabet}
const LongDNASeq       = LongSequence{DNAAlphabet{4}}
const LongRNASeq       = LongSequence{RNAAlphabet{4}}
const LongAminoAcidSeq = LongSequence{AminoAcidAlphabet}
const LongCharSeq      = LongSequence{CharAlphabet}

###
### Required type traits and methods
###

"Gets the alphabet encoding of a given BioSequence."
BioSymbols.alphabet(::Type{LongSequence{A}}) where {A} = alphabet(A)
Alphabet(::Type{LongSequence{A}}) where {A <: Alphabet} = A()
Base.length(seq::LongSequence) = last(seq.part) - first(seq.part) + 1
bindata(seq::LongSequence) = seq.data
Base.eltype(::Type{LongSequence{A}}) where {A} = eltype(A)

@inline function seq_data_len(::Type{A}, len::Integer) where A <: Alphabet
    return cld(len, div(64, bits_per_symbol(A())))
end

@inline function encoded_data(seq::LongSequence)
    return seq.data
end

# Replace a LongSequence's data with a copy, copying only what's needed.
# The user should never need to call this, as it has no outward effect on the
# sequence.
function orphan!(seq::LongSequence,
		 size::Integer = length(seq),
		 force::Bool = false)
	if !seq.shared & !force
	    return seq
	end
	return _orphan!(seq, size)
end

function _orphan!(seq::LongSequence{A},
		 size::Integer = length(seq)) where {A}

    j, r = bitindex(seq, 1)
    data = Vector{UInt64}(undef, seq_data_len(A, size))

    @inbounds if !isempty(seq) & !isempty(data)
        x = seq.data[j] >> r
        m = index(bitindex(seq, lastindex(seq))) - j + 1
        l = min(lastindex(data), m)
        @simd for i in 1:l-1
            y = seq.data[j + i]
            data[i] = x | y << (64 - r)
            x = y >> (r & 63)
        end
        if m <= l
            data[l] = x
        else
            y = seq.data[j + l]
            data[l] = x | y << (64 - r)
        end
    end

    seq.data = data
    seq.part = 1:size
    seq.shared = false
    return seq
end

include("indexing.jl")
include("constructors.jl")
include("printing.jl")
include("copying.jl")
include("conversion.jl")
include("stringliterals.jl")
include("transformations.jl")
include("operators.jl")
include("counting.jl")


text:
using DifferentialEquations
using StaticArrays
using FastClosures
#the Pauli matrices
function pauli_matrices(; type=Float64)
	σ1 = SMatrix{2, 2, Complex{type}}([0 1; 1 0])
	σ2 = SMatrix{2, 2, Complex{type}}([0 -1im; 1im 0])
	σ3 = SMatrix{2, 2, Complex{type}}([1 0; 0 -1])
	σ4 = SMatrix{2, 2, Complex{type}}([1 0; 0 1])
	σ = (σ1, σ2, σ3, σ4)
	return σ
end

function fidelity(U, U_ideal, dims)
	U_ideal_dag = adjoint(U_ideal)
	M = U_ideal_dag[1: dims, 1: dims]*U[1: dims, 1: dims]
	M_dag = adjoint(M)
	fidelity = (tr(M*M_dag) + abs(tr(M))^2)/(dims*(dims + 1))
	return real(fidelity)
end
#Decompose a generic 2x2 Hermitian matrix M in the Pauli basis.
function pauli_decomp(M)
	c = zeros(4)
	c[1] = real(M[1,2])
	c[2] = imag(M[2,1])
	c[3] = real(M[1,1] - M[2,2])/2
	c[4] = real(M[1,1] + M[2,2])/2
	return c
end
#Definite integral of the envelope function f.
function c(t, tf)
	theta = pi/2
	x = t/tf
	c_t = (x - sin(2*pi*x)/(2*pi))/2
	return theta*c_t
end
#Envelope function.
function f(t, tf; derivative::Bool = false)
	theta = pi/2
	x = t/tf
	if derivative == false
		return theta*(1 - cos(2*pi*x))/tf
	else
		return theta*sin(2*pi*x)*2*pi/(tf*tf)
	end
end
#Hamiltonian of an isolated 2 level system in the rotating frame.
function H0(t, tf, cr_terms)
	H0t = zeros(4)
	ft = f(t, tf)/2
	if cr_terms == false
		H0t[1] = ft
	else
		H0t[1] = ft*(1 + cos(2*t))
		H0t[2] = -ft*sin(2*t)
	end
	return H0t
end
#Differential equation function of the evolution operator.
function dU!(dU, U, args, t)
	Omega = args[1]
	σ = args[2]
	H = args[3]
	Omega_t = args[4]

	Omega_t = Omega(t)
	#Hamiltonian
	H .= zero(t)
	for i = 1: 4
		H .= H .+ Omega_t[i]*σ[i]
	end
	H .= -1im*H
	mul!(dU,H,U)
end
#Commutator of operators A and B.
function comm!(commAB, A, B)
	commAB[1] = (A[2]*B[3] - A[3]*B[2])*2
	commAB[2] = (A[3]*B[1] - A[1]*B[3])*2
	commAB[3] = (A[1]*B[2] - A[2]*B[1])*2
	return nothing
end
#Hamiltonian of a 2 level system in the interaction picture.
function V_I(t, tf)
	ft = f(t, tf)/2		#note the division by 2!!!
	c_t = c(t, tf)

	V_It = zeros(4)
	V_It[1] = ft*cos(2*t)
	V_It[2] = -ft*sin(2*t)*cos(2*c_t)
	V_It[3] = ft*sin(2*t)*sin(2*c_t)
	return V_It
end
#The lindbladian superoperator.
#ρ is the density operator. c is the collapse opperator.
function lindbladian!(lρ, c, ρ)
	cdag = adjoint(c)
	m = cdag*c
	lρ .= ρ*cdag
	lρ .= c*lρ .- (m*ρ .+ ρ*m)./2
	return nothing
end
#The master equations for the two level system.
#H is the Hamiltonian. c_list is a tuple with the collapse operators.
function dρ!(dρ, ρ, args, t)
	H = args[1]
	clist = args[2]
	ncollapse = arg[3]
	Ht = args[3]
	Hmatrix = args[4]
	lρ = args[5]

	Ht .= H(t)
	Hmatrix .= Ht[1].*σ[1]
	for i = 2: 4
		Hmatrix .= Hmatrix .+ Ht[i].*σ[i]
	end
	#unitary evolution
	dρ .= 1im.*(ρ*Hmatrix .- Hmatrix*ρ)
	#non-unitary evolution
	for i = 1: ncollapse
		lindbladian!(lρ, c1, ρ)
		dρ .= dρ .+ lρ
	end
end
#
function solveU(H, tf, σ)
	U0 = Matrix{Complex{Float64}}(I, 2, 2)
	tspan = (zero(tf), tf)
	args = (H, σ, zeros(Complex, 2, 2), zeros(4))
	prob = ODEProblem(dU!, U0, tspan, args)
	sol = solve(prob, Vern9(), saveat=[tf], reltol=1e-10, abstol=1e-14)
	return sol
end
#
function U(tf, σ; ideal_dynamics=false)
	cr_terms = !ideal_dynamics
	H = @closure t-> H0(t, tf, cr_terms)
	sol = solveU(H, tf, σ)
	return sol(tf)
end
#
function solveρ(tf, H, clist, ρ0)
	ρ0 = convert(Array{Complex{Float64}}, ρ0)
	#clist = (sqrt(γ1)*[[0 0];[1.0 0]], sqrt(γ2)*σ3)
	tspan = (zero(tf), tf)
	args = (H, clist, length(clist), zeros(4), zeros(Complex{Float64}, 2, 2),
			zeros(Complex{Float64}, 2, 2), zeros(Complex{Float64}, 2, 2))
	prob = ODEProblem(dρ!, ρ0, tspan, args)
	sol = solve(prob, Vern9(), saveat=[tf], reltol=1e-10, abstol=1e-14)
	return sol(tf)
end
#This function has the Bloch vector a as input and returns the density operator.
function ρbloch(a)
	u, v, w = a[1], a[2], a[3]
	ρ11, ρ12 = 1 + w, u - 1im*v
	ρ21, ρ22 = u + 1im*v, 1 - w
	#ρ = [[1+w u-1im*v];[u+1im*v 1-w]]/2
	ρ = [ρ11 ρ12; ρ21 ρ22]/2
	return ρ
end
#This function calculates the average fidelity in a two level system.
#Hideal is the ideal Hamiltonian; H is the nonideal Hamiltonian; γ1 and γ2 are
#the dissipation and dephasing constants.
function average_fidelity(tf, H, Hideal, clist)
	U = solveU(Hideal, tf, σ)
	dict = Dict("+x"=> (1,0,0), "-x"=> (-1,0,0),
				"+y"=> (0,1,0), "-y"=> (0,-1,0),
				"+z"=> (0,0,1), "-z"=> (0,0,-1))

	F = zero(tf)
	for tag in ("+x", "-x", "+y", "-y", "+z", "-z")
		a = dict[tag]
		ρ0 = ρbloch(a)
		ρf = solveρ(tf, H, clist, ρ0)
		ρfideal = U*ρ0*adjoint(U)
		F += tr(ρfideal*ρf)/6
	end
	if imag(F) > 10.0^(-10.0)
		error("Complex fidelity!")
		return -1
	else
		return real(F)
	end
end

max_stars_repo_name: thalesfr/MagnusCorrection

max_issues_repo_name: thalesfr/MagnusCorrection

max_forks_repo_name: thalesfr/MagnusCorrection


text:
# ---
# title: 646. Maximum Length of Pair Chain
# id: problem646
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/maximum-length-of-pair-chain/description/>
# hidden: true
# ---
#
# You are given `n` pairs of numbers. In every pair, the first number is always
# smaller than the second number.
#
# Now, we define a pair `(c, d)` can follow another pair `(a, b)` if and only if
# `b < c`. Chain of pairs can be formed in this fashion.
#
# Given a set of pairs, find the length longest chain which can be formed. You
# needn't use up all the given pairs. You can select pairs in any order.
#
# **Example 1:**
#
#
#
#     Input: [[1,2], [2,3], [3,4]]
#     Output: 2
#     Explanation: The longest chain is [1,2] -> [3,4]
#
#
# **Note:**
#
#   1. The number of given pairs will be in the range [1, 1000].
#
#
## @lc code=start
using LeetCode

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


text:
import Base: mean, std, repr, length
import IterTools: take

function set_rng!(solver::Solver, rng::AbstractRNG)
    solver.rng = rng
end

function test_run(eval_problem::NoCrashMDP, initial_state::MLState, solver_problem::NoCrashMDP, solver::Solver, rng_seed::Integer, max_steps=10000)
    set_rng!(solver, MersenneTwister(rng_seed))
    sim = POMDPSimulators.HistoryRecorder(rng=MersenneTwister(rng_seed), max_steps=max_steps, capture_exception=true)
    terminal_problem = deepcopy(eval_problem)
    terminal_problem.dmodel.lane_terminate = true
    r = simulate(sim, terminal_problem, solve(solver,solver_problem), initial_state)
    return sim
end

function test_run_return_policy(eval_problem::NoCrashMDP, initial_state::MLState, solver_problem::NoCrashMDP, solver::Solver, rng_seed::Integer, max_steps=10000)
    set_rng!(solver, MersenneTwister(rng_seed))
    sim = POMDPSimulators.HistoryRecorder(rng=MersenneTwister(rng_seed), max_steps=max_steps)
    policy = solve(solver, solver_problem)
    terminal_problem = deepcopy(eval_problem)
    terminal_problem.dmodel.lane_terminate = true
    r = simulate(sim, terminal_problem, policy, initial_state)
    return sim, policy
end

function run_simulations(stats::DataFrame,
                         objects::Dict;
                         parallel=true,
                         max_steps=10000,
                         progress=true,
                         desc="Progress: ")

    nb_sims = nrow(stats)

    all_solvers = Vector{Any}(undef, nb_sims)
    all_problems = Vector{Any}(undef, nb_sims)
    all_solver_problems = Vector{Any}(undef, nb_sims)
    all_initial = Vector{Any}(undef, nb_sims)

    solvers = objects["solvers"]
    problems = objects["problems"]
    initial_states = objects["initial_states"]

    for i in 1:nrow(stats)
        all_solvers[i] = deepcopy(solvers[stats[i,:solver_key]])
        all_problems[i] = problems[stats[i,:problem_key]]
        all_solver_problems[i] = problems[stats[i,:solver_problem_key]]
        all_initial[i] = initial_states[stats[i,:initial_key]]
    end

    return run_simulations(all_problems,
                           all_initial,
                           all_solver_problems,
                           all_solvers,
                           stats[:rng_seed],
                           parallel=parallel,
                           max_steps=max_steps,
                           progress=progress,
                           desc=desc)
end

function run_simulations(eval_problems::AbstractVector,
                         initial_states::AbstractVector,
                         solver_problems::AbstractVector,
                         solvers::AbstractVector,
                         rng_seeds::AbstractVector=collect(1:length(eval_problems));
                         parallel=true,
                         max_steps=10000,
                         progress=true,
                         desc="Progress: ")

    N = length(eval_problems)
    prog = ProgressMeter.Progress( N, dt=0.1, barlen=30, output=STDERR, desc=desc)
    if parallel
        if progress
            sims = pmap(test_run,
                        prog,
                        eval_problems,
                        initial_states,
                        solver_problems,
                        solvers,
                        rng_seeds,
                        max_steps*ones(Int,N)
                        )
        else
            sims = pmap(test_run,
                        eval_problems,
                        initial_states,
                        solver_problems,
                        solvers,
                        rng_seeds,
                        max_steps*ones(Int,N)
                        )
        end
    else
        sims = Vector{HistoryRecorder}(undef, length(eval_problems))
        if progress
            @showprogress for j in 1:length(eval_problems)
                sims[j] = test_run(eval_problems[j],
                                   initial_states[j],
                                   solver_problems[j],
                                   solvers[j],
                                   rng_seeds[j],
                                   max_steps)
            end
        else
            for j in 1:length(eval_problems)
                sims[j] = test_run(eval_problems[j],
                                   initial_states[j],
                                   solver_problems[j],
                                   solvers[j],
                                   rng_seeds[j],
                                   max_steps)
            end
        end
    end

    for (i,r) in enumerate(sims)
        if isa(r, RemoteException) || nothing !== (r.exception)
            println("Exception in simulation $(i)!")
            println("===============================")
            if isa(r, RemoteException)
                Base.showerror(STDOUT, r.captured.ex)
                println()
                # @show r.captured.processed_bt
                Base.show_backtrace(STDOUT, r.captured.processed_bt)
                println()
            else
                Base.showerror(STDOUT, r.exception)
                println()
                Base.show_backtrace(STDOUT, r.backtrace)
                println()
            end
            println("===============================")
            println("")
        end
    end

    # return a vector of history recorders
    return sims
end

get_brake_terminate_thresh(p::NoCrashProblem) = p.dmodel.brake_terminate_thresh
get_brake_terminate_thresh(p::Union{SuccessMDP, SuccessPOMDP}) = NA
get_lambda(p::NoCrashProblem) = p.rmodel.cost_dangerous_brake/p.rmodel.reward_in_target_lane
get_lambda(p::Union{SuccessMDP, SuccessPOMDP}) = NA

function fill_stats!(stats::DataFrame, objects::Dict, sims::Vector;
                     metrics::AbstractVector=get(objects, "metrics", []))
    @assert nrow(stats) == length(sims)
    problems = objects["problems"]
    nb_sims = length(sims)
    eval_problems = [problems[key] for key in stats[:problem_key]]

    # add new columns
    stats[:reward] = DataArray(Float64, nb_sims)
    stats[:brake_penalty_thresh] = DataArray(Float64, length(eval_problems))
    stats[:lambda] = DataArray(Float64, length(eval_problems))
    for i in 1:length(eval_problems)
        p = eval_problems[i]
        stats[:brake_penalty_thresh][i] = get_brake_terminate_thresh(p)
        stats[:lambda][i] = get_lambda(p)
    end
    stats[:brake_terminate_thresh] = [p.dmodel.brake_terminate_thresh for p in eval_problems]
    stats[:nb_brakes] = DataArray(Int, nb_sims)
    stats[:steps_to_lane] = DataArray(Int, nb_sims)
    stats[:time_to_lane] = DataArray(Float64, nb_sims)
    stats[:steps_in_lane] = DataArray(Int, nb_sims)
    stats[:steps] = Int[length(s.action_hist) for s in sims]
    stats[:crash] = DataArray(Bool, nb_sims)
    stats[:term] = DataArray(String, nb_sims)

    for m in metrics
        stats[key(m)] = DataArray(datatype(m), nb_sims)
    end

    for i in 1:nb_sims
        r = 0.0
        dt = eval_problems[i].dmodel.phys_param.dt
        steps_to_lane = nothing
        steps_in_lane = 0
        nb_brakes = 0
        crashed = false

        for (k,a) in enumerate(sims[i].action_hist)
            s = sims[i].state_hist[k]
            sp = sims[i].state_hist[k+1]
            t = (k-1)*dt
            r += reward(eval_problems[i], s, a, sp)

            nb_brakes += detect_braking(eval_problems[i], s, sp)

            if sp.cars[1].y == eval_problems[i].rmodel.target_lane
                steps_in_lane += 1
            end
            if s.cars[1].y == eval_problems[i].rmodel.target_lane
                if nothing === (steps_to_lane)
                    steps_to_lane = k-1
                end
            end

            if is_crash(eval_problems[i], s, sp)
                crashed = true
            end
        end

        if sims[i].state_hist[end].cars[1].y == eval_problems[i].rmodel.target_lane
            if nothing === (steps_to_lane)
                steps_to_lane = length(sims[i].state_hist)-1
            end
        end

        if nothing === (sims[i].state_hist[end].terminal)
            stats[:term][i] = ""
        else
            stats[:term][i] = string(sims[i].state_hist[end].terminal)
        end

        stats[:reward][i] = r
        stats[:nb_brakes][i] = nb_brakes
        stats[:steps_to_lane][i] = nothing === (steps_to_lane) ? NA : steps_to_lane
        stats[:time_to_lane][i] = stats[:steps_to_lane][i]*dt
        stats[:steps_in_lane][i] = steps_in_lane
        stats[:crash][i] = crashed

        for m in metrics
            stats[key(m)][i] = calculate(m, eval_problems[i], sims[i])
        end
    end
    return stats
end

function assign_keys(problems::Vector, initial_states::Vector; rng=MersenneTwister(rand(UInt32)))
    p_keys = String[randstring(rng) for p in problems]
    is_keys = String[randstring(rng) for is in initial_states]
    return Dict{String, Any}(
        "eval_problems"=>Dict{String,Any}(p_keys[i]=>problems[i] for i in 1:length(problems)),
        "initial_states"=>Dict{String,Any}(is_keys[i]=>initial_states[i] for i in 1:length(initial_states)),
    )
end

function assign_keys(problems::AbstractVector; rng=MersenneTwister(rand(UInt32)))
    p_keys = String[randstring(rng) for p in problems]
    return Dict{String,Any}([(p_keys[i], problems[i]) for i in 1:length(problems)])
end

"""
Checks that all common keys have equal values in both dictionaries and adds keys not in common
"""
function careful_merge!(d1::Dict, d2::Dict)
    for k in keys(d2)
        if haskey(d1, k)
            @assert d1[k] == d2[k]
        else
            d1[k] = d2[k]
        end
    end
    return d1
end

function merge_results!(r1::Dict{S1, Any}, r2::Dict{S2, Any}; careful=true) where {S1<:AbstractString, S2<:AbstractString}
    merge!(r1["solvers"], r2["solvers"])
    merge!(r1["problems"], r2["problems"])
    if careful
        careful_merge!(r1["initial_states"], r2["initial_states"])
    else
        merge!(r1["initial_states"], r2["initial_states"])
    end
    merge!(r1["behaviors"], r2["behaviors"])
    r1["param_table"] = join(r1["param_table"], r2["param_table"], on=names(r1["param_table"]), kind=:outer)
    if haskey(r1, "stats")
        # len = nrow(r1["stats"])
        append!(r1["stats"], r2["stats"])
        # r1["stats"][:id][end-len+1:end] += len
        r1["stats"][:id]=1:nrow(r1["stats"])
    else
        r1["stats"] = r2["stats"]
    end
    if haskey(r2, "histories")
        if haskey(r1, "histories")
            if r1["histories"] !== nothing && r2["histories"] !== nothing
                append!(r1["histories"], r2["histories"])
            end
        else
        r1["histories"] = r2["histories"]
        end
    end

    r1ips = r1["initial_physical_states"]
    if length(r2["initial_physical_states"]) > length(r1ips)
        r1["initial_physical_states"] = r2["initial_physical_states"]
    end
    if careful
        for i in 1:min(length(r1ips), length(r2["initial_physical_states"]))
            @assert r1ips[i] == r2["initial_physical_states"][i]
        end
    end

    merge_tests!(r1["tests"], r2["tests"])
    return r1
end

function rerun(results::Dict{S, Any}, id; enforce_match=true) where S<:AbstractString
    stats = results["stats"]
    @assert stats[:id][id] == id
    problem = results["problems"][stats[:problem_key][id]]
    is = results["initial_states"][stats[:initial_key][id]]
    solver = results["solvers"][stats[:solver_key][id]]
    solver_problem = results["problems"][stats[:solver_problem_key][id]]
    rng_seed = stats[:rng_seed][id]
    if enforce_match
        steps = stats[:steps][id]
    else
        steps = 10000
    end
    sim, policy = test_run_return_policy(problem, is, solver_problem, solver, rng_seed, steps)
    r = sum([reward(problem, sim.state_hist[i], sim.action_hist[i], sim.state_hist[i+1]) for i in 1:length(sim.action_hist)])
    if enforce_match
        @assert r == stats[:reward][id]
    end
    return problem, sim, policy
end


text:
@rep_op "a | n > = sqrt(n)| n-1 >" 0:4
@assert d" act_on(a, | 1,2,3 > + 2| 3,3,1 >, 2) == sqrt(2)| 1,1,3 > + 2*sqrt(3)| 3,2,1 > "
@assert d" act_on(a', | 1,2,3 > + 2| 3,2,1 >, 2) == sqrt(3)| 1,3,3 > + 2*sqrt(3)| 3,3,1 > "

k3 = sum(i->d" (i + (i * im)) * | i > ", 1:2)^3

@rep_op "a_2 | x,y,z > = sqrt(y) * | x, y-1, z >" 0:3 0:3 0:3
@rep_op "< x,y,z | a_2_c = sqrt(y+1) * < x, y+1, z |" 0:2 0:2 0:2

@def_op "lowf_2 | x,y,z > = sqrt(y) * | x, y-1, z >"
@def_op "< x,y,z | lowf_2 = sqrt(y + 1) * < x, y+1, z |"

filt3s(k,c) = ! (in(3, klabel(k)) || in(3, blabel(k)))

@assert filter(filt3s, a_2) == filter(filt3s, a_2_c)
@assert a_2 * k3 == act_on(a, k3, 2) == lowf_2 * k3  == lower(k3, 2)
@assert a_2' * k3 == act_on(a', k3, 2) == lowf_2' * k3 == raise(k3, 2)
@assert k3' * a_2  ==  k3' * lowf_2  == raise(k3', 2)
@assert k3' * a_2'  ==  k3' * lowf_2'  == lower(k3', 2)


text:
module LOLTools

include("helpers.jl")

# This API is being deprecated on Monday, September 13th, 2021.
# Please use match-v5 as a replacement.
export MatchV4
include("MatchV4.jl")

export MatchV5
include("MatchV5.jl")

export SummonerV4
include("SummonerV4.jl")

export SpectatorV4
include("SpectatorV4.jl")

export LeagueV4
include("LeagueV4.jl")

export LeagueExpV4
include("LeagueExpV4.jl")

export ChampionMasteryV4
include("ChampionMasteryV4.jl")

export ChampionV3
include("ChampionV3.jl")

export LoLStatusV3
include("LoLStatusV3.jl")

export GameClient
include("GameClient.jl")

end # module LOLTools


text:
abstract type AbstractMsgFormat end

"""
Intermediate format based on JSON.
A JSON object with `cmd` (string), `args` (array), `vargs` (dict).
"""
struct JSONMsgFormat <: AbstractMsgFormat
end

wireformat(fmt::JSONMsgFormat, cmd::String, args...; data...) = JSON.json(_dict_fmt(cmd, args...; data...))
wireformat(fmt::JSONMsgFormat, code::Int, headers::Dict{String,String}, resp, id=nothing) = JSON.json(_dict_fmt(code, headers, resp, id))
juliaformat(fmt::JSONMsgFormat, msgstr) = JSON.parse(msgstr)

"""
Intermediate format based on Julia serialization.
A dict with `cmd` (string), `args` (array), `vargs` (dict).
"""
struct SerializedMsgFormat <: AbstractMsgFormat
end

wireformat(fmt::SerializedMsgFormat, cmd::String, args...; data...) = _dict_ser(_dict_fmt(cmd, args...; data...))
wireformat(fmt::SerializedMsgFormat, code::Int, headers::Dict{String,String}, resp, id=nothing) = _dict_ser(_dict_fmt(code, headers, resp, id))
juliaformat(fmt::SerializedMsgFormat, msgstr) = _dict_dser(msgstr)

"""
Intermediate format that is just a Dict. No serialization.
A Dict with `cmd` (string), `args` (array), `vargs` (dict).
"""
struct DictMsgFormat <: AbstractMsgFormat
end

wireformat(fmt::DictMsgFormat, cmd::String, args...; data...) = _dict_fmt(cmd, args...; data...)
wireformat(fmt::DictMsgFormat, code::Int, headers::Dict{String,String}, resp, id=nothing) = _dict_fmt(code, headers, resp, id)
juliaformat(fmt::DictMsgFormat, msg) = msg

##############################################
# common interface methods for message formats
##############################################
cmd(fmt::AbstractMsgFormat, msg) = get(msg, "cmd", "")
args(fmt::AbstractMsgFormat, msg) = get(msg, "args", [])
function data(fmt::AbstractMsgFormat, msg)
    vargs = get(msg, "vargs", Dict{Symbol,Any}())
    isa(vargs, Dict{Symbol,Any}) && (return vargs)
    dict = Dict{Symbol,Any}()
    for (n,v) in vargs
        dict[Symbol(n)] = v
    end
    dict
end

"""
extract and return the response data as a direct function call would have returned
but throw error if the call was not successful.
"""
fnresponse(fmt::AbstractMsgFormat, resp) = _dict_fnresponse(resp)

"""construct an HTTP Response object from the API response"""
httpresponse(fmt::AbstractMsgFormat, resp) = _dict_httpresponse(resp)

##############################################
# common utility methods for message formats
##############################################
function data_dict(data)
    d = Dict{Symbol,Any}()
    for (n,v) in data
        d[n] = v
    end
    d
end

function _dict_ser(d)
    io = IOBuffer()
    serialize(io, d)
    take!(io)
end
_dict_dser(b) = deserialize(IOBuffer(b))

function _dict_fmt(cmd::String, args...; data...)
    req = Dict{String,Any}()

    req["cmd"] = cmd
    isempty(args) || (req["args"] = args)
    isempty(data) || (req["vargs"] = data_dict(data))
    req
end

function _dict_fmt(code::Int, headers::Dict{String,String}, resp, id=nothing)
    msg = Dict{String,Any}()

    (id == nothing) || (msg["nid"] = id)

    if !isempty(headers)
        msg["hdrs"] = headers
        @debug("sending headers", headers)
    end

    msg["code"] = code
    msg["data"] = resp
    msg
end

function _dict_httpresponse(resp)
    hdrs = Dict{String,String}()
    if "hdrs" in keys(resp)
        for (k,v) in resp["hdrs"]
            hdrs[k] = v
        end
    end
    data = get(resp, "data", "")
    respdata = isa(data, Array) ? convert(Array{UInt8}, data) :
               isa(data, Dict) ? JSON.json(data) :
               string(data)

    HTTP.Response(resp["code"], hdrs; body=respdata)
end

function _dict_fnresponse(resp)
    data = get(resp, "data", "")
    (resp["code"] == ERR_CODES[:success][1]) || error("API error: $data")
    data
end


text:
module Mining

ccall(:jl_exit_on_sigint, Cvoid, (Cint,), 0)

include("model.jl")

using LightGraphs
using StatsBase

"""
Decompose graph into components to serve as seeds for structure candidates
"""
function decompose_graph(G, minimum_component_size::Int64=10, neighborhood_size::Int64=1, return_slashed::Bool=false)
    components = []
    slashed = []

    ccs = sort!(connected_components(G), by=length, rev=true)
    lcc = ccs[1]
    lccG, vmap = induced_subgraph(G, lcc)
    largest_degree = Δ(lccG)

    while length(lcc) >= minimum_component_size && largest_degree >= minimum_component_size
        new_hub = sort!(filter(x->x[end] == largest_degree, collect(zip(vertices(lccG),degree(lccG)))), by=x->x[1])[1][1]
        new_hub_original = vmap[new_hub]
        push!(slashed, new_hub_original)
        hub_neighbors = neighborhood(G, new_hub_original, neighborhood_size)
        push!(components, hub_neighbors)
        for edge in filter(x -> src(x) == new_hub_original || dst(x) == new_hub_original, collect(edges(G)))
            rem_edge!(G, edge)
        end
        ccs = sort!(connected_components(G), by=length, rev=true)
        lcc = ccs[1]
        lccG, vmap = induced_subgraph(G, lcc)
        largest_degree = Δ(lccG)
    end

    filter!(x->length(x) >= minimum_component_size, ccs)
    append!(components, ccs)
    append!(slashed, repeat([nothing],length(ccs)))
    if !return_slashed
        return components
    else
        return components, slashed
    end
end

"""
Create a star candidate from a component
"""
function find_star_from_component(component, G)
    cG, vmap = induced_subgraph(G, component)
    largest_degree = Δ(cG)
    top_nodes = sort!(filter(x->x[end] == largest_degree, collect(zip(vertices(cG),degree(cG)))), by=x->x[1])
    if length(top_nodes) > 1
        return nothing, [], nothing
    else
        new_hub = top_nodes[1][1]
        new_hub_original = vmap[new_hub]
        sG, vmap = induced_subgraph(G, neighborhood(G, new_hub_original, 1))
        nodes_original = [vmap[x] for x in vertices(sG)]
        n_spokes = nv(sG) - 1
        fraction = 0.1
        candidates = sort!(filter(x -> n_spokes - 1 > x[end] - 1 > 0.05 * n_spokes, collect(zip(vertices(sG),degree(sG)))), by=x->x[end], rev=true)
        cut_point = ceil(fraction * length(candidates))
        nodes_to_cut = candidates[1:convert(Int, cut_point)]
        while !isempty(nodes_to_cut) && nv(sG) >= 10
            nodes_to_cut_original = [vmap[x[1]] for x in nodes_to_cut]
            nodes_to_keep = filter(x -> x ∉ nodes_to_cut_original, nodes_original)
            sG, vmap = induced_subgraph(G, nodes_to_keep)
            nodes_original = [vmap[x] for x in vertices(sG)]
            n_spokes = nv(sG) - 1
            fraction = min(fraction + 0.01, 1.0)
            candidates = sort!(filter(x -> n_spokes - 1 > x[end] - 1 > 0.05 * n_spokes, collect(zip(vertices(sG),degree(sG)))), by=x->x[end], rev=true)
            cut_point = ceil(fraction * length(candidates))
            nodes_to_cut = candidates[1:convert(Int, cut_point)]
        end
        if nv(sG) < 10
            return nothing, [], nothing
        else
            return sG, [vmap[x] for x in vertices(sG) if vmap[x] != new_hub_original], new_hub_original
        end
    end
end

"""
Create a clique candidate from a component
"""
function find_clique_from_component(component, G, min_connectivity_fraction=0.5)
    cG, vmap = induced_subgraph(G, component)

    maxcliques = sort!(maximal_cliques(cG), by=x->length(x), rev=true)
    nodes_original = [vmap[x] for x in maxcliques[1]]
    cG, vmap = induced_subgraph(G, nodes_original)

    while true
        cutoff = min_connectivity_fraction * length(nodes_original)
        nodes_to_add = sort!([x[1] for x in filter(x -> x[end] >= cutoff && x[1] ∉ nodes_original,
                                collect(countmap(vcat([neighbors(G, x) for x in nodes_original]...))))], by=x->degree(G,x), rev=true)
        if !isempty(nodes_to_add)
            node_to_add = nodes_to_add[1]
        else
            break
        end
        union!(nodes_original, [node_to_add])
        cG, vmap = induced_subgraph(G, nodes_original)
    end

    if nv(cG) < 10 || !is_connected(cG)
        return nothing, []
    else
        return cG, nodes_original
    end
end

"""
Create a biclique candidate from a component
"""
function find_biclique_from_component(component, G)
    max_intra_connectivity = 0.05
    min_inter_connectivity = 0.5
    cG, vmap = induced_subgraph(G, component)

    right = independent_set(cG, MaximalIndependentSet(), seed=1234)
    original_right = [vmap[x] for x in right]
    sort!(original_right, by=x->degree(G, x), rev=true)
    original_right = original_right[1:min(length(original_right),5)]
    cG, vmap = induced_subgraph(G, original_right)

    left_candidates = [x[1] for x in filter(x -> x[end] >= min_inter_connectivity * length(original_right) && x[1] ∉ original_right, collect(countmap(vcat([neighbors(G, x) for x in original_right]...))))]
    lG, vmap2 = induced_subgraph(G, left_candidates)
    left = independent_set(lG, MaximalIndependentSet(), seed=1234)
    original_left = [vmap2[x] for x in left]
    sort!(original_left, by=x->degree(G, x), rev=true)
    original_left = original_left[1:min(length(original_left),5)]
    left_to_add = [x[1] for x in filter(x -> x[end] <= max_intra_connectivity * length(original_left) && x in left_candidates, collect(countmap(vcat([neighbors(G, x) for x in original_left]...))))]
    if !isempty(left_to_add)
        union!(original_left,left_to_add)
    end

    if length(original_left) < 3 || length(original_right) < 5
        return nothing, [], []
    end

    while true
        added_left = filter(
                    x ->
                    x[1] ∉ original_left
                    && x[1] ∉ original_right
                    && x[end] >= min_inter_connectivity * length(original_right)
                    && length(intersect(neighbors(G, x[1]), original_left)) <= max_intra_connectivity * length(original_left),
                        collect(countmap(vcat([neighbors(G, x) for x in original_right]...)))
                    )
        if !isempty(added_left)
            added_left = map(x->x[1],sort!(added_left, by=x->x[end], rev=true))
            append!(original_left, [added_left[1]])
        end

        added_right = filter(
                    x ->
                    x[1] ∉ original_left
                    && x[1] ∉ original_right
                    && x[end] >= min_inter_connectivity * length(original_left)
                    && length(intersect(neighbors(G, x[1]), original_right)) <= max_intra_connectivity * length(original_right),
                        collect(countmap(vcat([neighbors(G, x) for x in original_left]...)))
                    )
        if !isempty(added_right)
            added_right = map(x->x[1],sort!(added_right, by=x->x[end], rev=true))
            append!(original_right, [added_right[1]])
        end

        if isempty(added_left) && isempty(added_right)
            break
        end
    end

    cG, vmap = induced_subgraph(G, union(original_left, original_right))
    isolates = [vmap[n] for n in vertices(cG) if degree(cG, n) == 0]
    setdiff!(original_left, isolates)
    setdiff!(original_right, isolates)
    cG, vmap = induced_subgraph(G, union(original_left, original_right))

    if nv(cG) < 10 || !is_connected(cG) #|| length(original_left) < 3 || length(original_right) < 3
        return nothing, [], []
    else
        return cG, original_left, original_right
    end
end

"""
Create a starclique candidate from a component
"""
function find_starclique_from_component(component, G)
    max_intra_connectivity = 0.05
    min_intra_connectivity = 0.5
    min_inter_connectivity = 0.5
    cG, vmap = induced_subgraph(G, component)

    maxcliques = sort!(maximal_cliques(cG), by=x->length(x), rev=true)
    original_left = [vmap[x] for x in maxcliques[1]]
    cG, vmap = induced_subgraph(G, original_left)

    right_candidates = [x[1] for x in filter(x -> x[end] >= min_inter_connectivity * length(original_left) && x[1] ∉ original_left, collect(countmap(vcat([neighbors(G, x) for x in original_left]...))))]
    rG, rvmap = induced_subgraph(G, right_candidates)
    right = independent_set(rG, MaximalIndependentSet(), seed=1234)
    original_right = [rvmap[x] for x in right]

    while true
        added_left = filter(
                    x ->
                    x[1] ∉ original_left
                    && x[1] ∉ original_right
                    && x[end] >= min_inter_connectivity * length(original_right)
                    && length(intersect(neighbors(G, x[1]), original_left)) >= min_intra_connectivity * length(original_left),
                        collect(countmap(vcat([neighbors(G, x) for x in original_right]...)))
                    )
        if !isempty(added_left)
            added_left = map(x->x[1],sort!(added_left, by=x->x[end], rev=true))
            append!(original_left, [added_left[1]])
        end

        added_right = filter(
                    x ->
                    x[1] ∉ original_left
                    && x[1] ∉ original_right
                    && x[end] >= min_inter_connectivity * length(original_left)
                    && length(intersect(neighbors(G, x[1]), original_right)) <= max_intra_connectivity * length(original_right),
                        collect(countmap(vcat([neighbors(G, x) for x in original_left]...)))
                    )
        if !isempty(added_right)
            added_right = map(x->x[1],sort!(added_right, by=x->x[end], rev=true))
            append!(original_right, [added_right[1]])
        end

        if isempty(added_left) && isempty(added_right)
            break
        end
    end

    cG, vmap = induced_subgraph(G, union(original_left, original_right))
    isolates = [vmap[n] for n in vertices(cG) if degree(cG, n) == 0]
    setdiff!(original_left, isolates)
    setdiff!(original_right, isolates)
    cG, vmap = induced_subgraph(G, union(original_left, original_right))
    if nv(cG) < 10 || !is_connected(cG) || length(original_left) < 3 || length(original_right) < 3
        return nothing, [], []
    else
        return cG, original_left, original_right
    end
end

"""
Helper for sorting candidates
"""
function get_structure_type_priority(structure_type)
    if structure_type == "biclique"
        return 4
    elseif structure_type == "starclique"
        return 3
    elseif structure_type == "star"
        return 2
    else # "clique"
        return 1
    end
end

"""
Model building function which does the heavy lifting
"""
function build_decomposition_based_model!(M::Model.GraphModel, size_threshold::Float64=10., stop_threshold::Float64=10., stop_after_n_structures::Float64=Inf, sort_heuristic::String="plain")
    minimum_component_size = convert(Int64, ceil(size_threshold))

    G = deepcopy(M.G)
    Model.log_progress!(M, "Created graph with $(nv(G)) nodes and $(ne(G)) edges...")
    components = decompose_graph(G, minimum_component_size, 1, false)
    Model.log_progress!(M, "Decomposed graph")

    stars, star_spokes_original, hubs = generate_star_component_candidates(M, minimum_component_size, components)
    Model.log_progress!(M, "Generated star candidates")
    cliques, clique_nodes_original = generate_clique_component_candidates(M, minimum_component_size, components)
    Model.log_progress!(M, "Generated clique candidates")
    bicliques, left_nodes_original, right_nodes_original = generate_biclique_component_candidates(M, minimum_component_size, components)
    Model.log_progress!(M, "Generated biclique candidates")
    starcliques, s_left_nodes_original, s_right_nodes_original = generate_starclique_component_candidates(M, minimum_component_size, components)
    Model.log_progress!(M, "Generated starclique candidates")

    star_structures = generate_all_star_structures(hubs, star_spokes_original, stars)
    clique_structures = generate_all_clique_structures(clique_nodes_original, cliques)
    biclique_structures = generate_all_biclique_structures(left_nodes_original, right_nodes_original, bicliques, M.G, "biclique")
    starclique_structures = generate_all_biclique_structures(s_left_nodes_original, s_right_nodes_original, starcliques, M.G, "starclique")
    all_structures = vcat(clique_structures, star_structures, biclique_structures, starclique_structures)
    sort!(all_structures, by=x -> (x.n_nodes_total, x.n_edges_total, get_structure_type_priority(x.structure_type)), rev=true)
    Model.log_progress!(M, "List of candidate structures, ordered by size:")
    for structure in all_structures
        Model.log_progress!(M, "$(structure)")
    end
    Model.log_progress!(M, "")

    successive_failures = 0
    total_failures = 0
    while !isempty(all_structures)
        Model.log_progress!(M, "$(length(all_structures)) structures remaining")
        structure = popfirst!(all_structures)
        if Model.test_adding_structure!(M, structure, size_threshold)
            Model.log_progress!(M, "added: $(structure)")
            Model.log_progress!(M, "number of equivalence classes: $(length(M.equivalence_class_edge_probabilities))")
            Model.log_progress!(M, "number of lambdas: $(length(M.lambdas))")
            successive_failures = 0
            union!(M.covered_nodes,Mining.Model.Structures.get_node_set(structure))
            if sort_heuristic == "uncovered_nodes"
                sort!(all_structures, by=x ->  (get_uncovered_score(x, M), x.n_nodes_total, x.n_edges_total), rev=true)
            end
        else
            successive_failures += 1
            total_failures += 1
            if total_failures >= stop_threshold
                Model.log_progress!(M, "Stopping model computation - $(total_failures) successive failures...")
                break
            end
        end
        if length(M.macro_structures) - 2 >= stop_after_n_structures
            Model.log_progress!(M, "Stopping model computation - $(stop_after_n_structures) structures found...")
            break
        end
    end
end

"""
Helper for sort heuristic 'uncovered_nodes'
"""
function get_uncovered_score(structure, M)
    nodes = Mining.Model.Structures.get_node_set(structure)
    n_uncovered_nodes = length(setdiff(nodes, M.covered_nodes))
    n_covered_nodes = length(nodes) - n_uncovered_nodes
    return n_uncovered_nodes * (n_uncovered_nodes - 1) / 2 + n_uncovered_nodes * n_covered_nodes # this can be optimistic or pessimistic depending on the structure
end

"""
Create *all* star candidates (raw)
"""
function generate_star_component_candidates(M::Model.GraphModel, minimum_component_size::Int64=10, components=[])
    stars = []
    star_spokes_original = []
    hubs = []
    for c in components
        star, original_nodes, hub = find_star_from_component(c, M.G)
        if star !== nothing
            push!(stars,star)
            push!(star_spokes_original, original_nodes)
            push!(hubs, hub)
        end
    end
    return stars, star_spokes_original, hubs
end

"""
Create *all* clique candidates (raw)
"""
function generate_clique_component_candidates(M::Model.GraphModel, minimum_component_size::Int64=10, components=[])
    cliques = []
    clique_nodes_original = []
    for c in components
        clique, original_nodes = find_clique_from_component(c, M.G)
        clique_length = length(original_nodes)
        found = false
        for (idx,other_nodes) in enumerate(clique_nodes_original)
            if length(intersect(original_nodes, other_nodes))/min(clique_length,length(other_nodes)) >= 0.9 # was: 0.95
                clique_nodes_original[idx] = union(original_nodes, other_nodes)
                cliques[idx], _ = induced_subgraph(M.G, clique_nodes_original[idx])
                found = true
                break
            end
        end
        if !found && clique !== nothing
            push!(cliques, clique)
            push!(clique_nodes_original, original_nodes)
        end
    end
    return cliques, clique_nodes_original
end

"""
Create *all* biclique candidates (raw)
"""
function generate_biclique_component_candidates(M::Model.GraphModel, minimum_component_size::Int64=10, components=[], bipartite_type="biclique")
    bicliques = []
    biclique_left_nodes = []
    biclique_right_nodes = []
    for c in components
        biclique, original_nodes_left, original_nodes_right = find_biclique_from_component(c, M.G)
        left_length = length(original_nodes_left)
        right_length = length(original_nodes_right)
        found = false
        for (idx,(left,right)) in enumerate(zip(biclique_left_nodes, biclique_right_nodes))
            if length(intersect(original_nodes_left, left))/min(left_length,length(left)) >= 0.90 && length(intersect(original_nodes_right, right))/min(right_length,length(right)) >= 0.90 # was: 0.9
                biclique_left_nodes[idx] = union(original_nodes_left, left)
                biclique_right_nodes[idx] = setdiff(union(original_nodes_right, right), biclique_left_nodes[idx])
                bicliques[idx], _ = induced_subgraph(M.G, union(biclique_left_nodes[idx], biclique_right_nodes[idx]))
                found = true
                break
            end
        end
        if !found && biclique !== nothing
            push!(bicliques, biclique)
            push!(biclique_left_nodes, original_nodes_left)
            push!(biclique_right_nodes, original_nodes_right)
        end
    end
    return bicliques, biclique_left_nodes, biclique_right_nodes
end

"""
Create *all* starclique candidates (raw)
"""
function generate_starclique_component_candidates(M::Model.GraphModel, minimum_component_size::Int64=10, components=[])
    starcliques = []
    starclique_left_nodes = []
    starclique_right_nodes = []
    for c in components
        starclique, original_nodes_left, original_nodes_right = find_starclique_from_component(c, M.G)
        left_length = length(original_nodes_left)
        right_length = length(original_nodes_right)
        found = false
        for (idx,(left,right)) in enumerate(zip(starclique_left_nodes, starclique_right_nodes))
            if length(intersect(original_nodes_left, left))/min(left_length,length(left)) >= 0.9 && length(intersect(original_nodes_right, right))/min(right_length,length(right)) >= 0.9
                starclique_left_nodes[idx] = union(original_nodes_left, left)
                starclique_right_nodes[idx] = setdiff(union(original_nodes_right, right), starclique_left_nodes[idx])
                starcliques[idx], _ = induced_subgraph(M.G, union(starclique_left_nodes[idx], starclique_right_nodes[idx]))
                found = true
                break
            end
        end
        if !found && starclique !== nothing
            push!(starcliques, starclique)
            push!(starclique_left_nodes, original_nodes_left)
            push!(starclique_right_nodes, original_nodes_right)
        end
    end
    return starcliques, starclique_left_nodes, starclique_right_nodes
end

"""
Create *all* star candidates (as objects)
"""
function generate_all_star_structures(hubs, star_spokes_original, stars)
    star_structures = []
    for (hub, spokes, sG) in zip(hubs, star_spokes_original, stars)
        if !isempty(spokes)
            s = Mining.Model.Structures.Star(ne(sG), hub, spokes)
            push!(star_structures, s)
        end
    end
    return star_structures
end

"""
Create *all* clique candidates (as objects)
"""
function generate_all_clique_structures(clique_nodes_original, cliques)
    clique_structures = []
    clique_nodes = []
    for (nodes, cG) in zip(clique_nodes_original, cliques)
        for (idx,elem) in enumerate(map(elem -> elem ⊆ nodes, clique_nodes))
            if elem
                clique_nodes[idx] = nodes
                clique_structures[idx] = Mining.Model.Structures.Clique(ne(cG), nodes)
            end
        end
        if !any(map(elem -> nodes ⊆ elem, clique_nodes))
            if !isempty(nodes)
                c = Mining.Model.Structures.Clique(ne(cG), nodes)
                push!(clique_structures, c)
                push!(clique_nodes, Set(nodes))
            end
        end
    end
    return clique_structures
end

"""
Create *all* biclique or starclique candidates (as objects)
"""
function generate_all_biclique_structures(biclique_left_nodes, biclique_right_nodes, bicliques, G, bipartite_type="biclique")
    biclique_structures = []
    for (l, r, bG) in zip(biclique_left_nodes, biclique_right_nodes, bicliques)
        if !isempty(l)
            n_edges_in_left = ne(induced_subgraph(G, l)[1])
            n_edges_in_right = ne(induced_subgraph(G, r)[1])
            if bipartite_type == "biclique"
                b = Mining.Model.Structures.Biclique(ne(bG), l, r, n_edges_in_left, n_edges_in_right)
            else # assume "starclique"
                b = Mining.Model.Structures.Starclique(ne(bG), l, r, n_edges_in_left, n_edges_in_right)
            end
            push!(biclique_structures, b)
        end
    end
    return biclique_structures
end

"""
Model building entry function
"""
function compute_anytime_model(dataset::String, surprise_threshold::Float64, size_threshold::Float64, stop_threshold::Float64, stop_after_n_structures::Float64, sort_heuristic::String; overwrite::Bool=false, debug::Bool=false)
    dataset_core = splitdir(splitdir(dataset)[1])[end]*"/"*splitext(splitdir(dataset)[end])[1]
    logdir, logfile = splitdir("../results/$(dataset_core)_size-$(size_threshold)_max-$(stop_after_n_structures)_sort-$(sort_heuristic).log")
    # we don't want to accidentally overwrite logs when there is a corresponding model already
    if !overwrite && isfile(logdir * "/" * logfile[1:end-3] * "json")
        println("Skipping $(logfile) because there exists a corresponding model and overwrite is false...")
        return
    end
    mkpath(logdir)
    open("$(logdir)/$(logfile)", "w") do io
        M = Model.GraphModel(dataset, io)
        M.parameters["surprise_threshold"] = surprise_threshold
        M.parameters["size_threshold"] = size_threshold
        M.parameters["stop_after_n_structures"] = stop_after_n_structures
        if debug
            Model.log_progress!(M, "starting description length: $(M.start_description_length_maxent)\n")
            build_decomposition_based_model!(M, size_threshold, stop_threshold, stop_after_n_structures, sort_heuristic)
            Model.save_model(M, "../results"; overwrite=overwrite)
        else
            try
                Model.log_progress!(M, "starting description length: $(M.start_description_length_maxent)\n")
                build_decomposition_based_model!(M, size_threshold, stop_threshold, stop_after_n_structures, sort_heuristic)
            catch InterruptException
                println("\nCatching interrupt exception...")
            finally
                println("Saving state...")
                Model.save_model(M, "../results"; overwrite=overwrite)
                println("Done.")
            end
        end
    end
end

end

max_stars_repo_name: dataspider/momo

max_issues_repo_name: dataspider/momo

max_forks_repo_name: dataspider/momo


text:
module Spglib
using AffineTransforms

# Load library after build
if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
  include(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
else
  error("Spglib not properly installed. Please run Pkg.build(\"Spglib\")")
end

" version number of the underlying C-API "
const version = VersionNumber(
  ccall((:spg_get_major_version, spglib), Cint, ()),
  ccall((:spg_get_minor_version, spglib), Cint, ()),
  ccall((:spg_get_micro_version, spglib), Cint, ()),
)

function symmetry_operations(lattice::AbstractMatrix,
                             positions::AbstractMatrix,
                             types::AbstractVector;
                             symprec::Real=1e-8)
  if size(positions, 2) != length(types)
    error("Number of positions and types do not match")
  end
  if size(positions, 1) != 3
    error("Operating in 3D here")
  end
  const maxsize::Integer = 52
  rotations = Array{Cint}((3, 3, maxsize))
  translations = Array{Cdouble}((3, maxsize))

  unique_types = unique(types)
  type_indices = Cint[findfirst(unique_types, u) for u in types]
  clattice = convert(Matrix{Cdouble}, lattice)
  cpositions = convert(Matrix{Cdouble}, positions)
  numops = ccall(
      (:spg_get_symmetry, spglib),
      Cint,
      (Ptr{Cint}, Ptr{Cdouble}, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint},
       Cint, Cdouble),
    rotations, translations, maxsize, clattice, cpositions, type_indices, length(type_indices), symprec
  )
  if numops == 0
    error("Could not determine symmetries")
  end
  [AffineTransform(transpose(rotations[:, :, i]), translations[:, i]) for i in 1:numops]
end

function international_symbol(lattice::AbstractMatrix,
                              positions::AbstractMatrix,
                              types::AbstractVector;
                              symprec::Real=1e-8)
  result = zeros(Cchar, 11)

  unique_types = unique(types)
  type_indices = Cint[findfirst(unique_types, u) for u in types]
  clattice = convert(Matrix{Cdouble}, lattice)
  cpositions = convert(Matrix{Cdouble}, positions)
  numops = ccall(
      (:spg_get_international, spglib),
      Cint,
      (Ptr{Cchar}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint},
       Cint, Cdouble),
    result, clattice, cpositions, type_indices, length(type_indices), symprec
  )

  if numops == 0
    error("Could not determine internation symbol")
  end
  join(convert(Array{Char}, result[1:findfirst(result, 0) - 1]))
end

function schoenflies_symbol(lattice::AbstractMatrix,
                           positions::AbstractMatrix,
                           types::AbstractVector;
                           symprec::Real=1e-8)
  result = zeros(Cchar, 11)

  unique_types = unique(types)
  type_indices = Cint[findfirst(unique_types, u) for u in types]
  clattice = convert(Matrix{Cdouble}, lattice)
  cpositions = convert(Matrix{Cdouble}, positions)
  numops = ccall(
      (:spg_get_schoenflies, spglib),
      Cint,
      (Ptr{Cchar}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cint},
       Cint, Cdouble),
    result, clattice, cpositions, type_indices, length(type_indices), symprec
  )

  if numops == 0
    error("Could not determine Schoenflies symbol")
  end
  join(convert(Array{Char}, result[1:findfirst(result, 0) - 1]))
end

export symmetry_operations, international_symbol, schoenflies_symbol
end # module


text:
"""
Init basic C libraries
"""
function InitBase()
  """
  Base.__init__()
  Sys.__init__() #fix https://github.com/JuliaLang/julia/issues/30479
  """
end

"""
# Initialize REPL module for Docs
"""
function InitREPL()
  """
  using REPL
  Base.REPL_MODULE_REF[] = REPL
  """
end

function Include(path)
  """
  M = Module()
  # Include into anonymous module to not polute namespace
  @eval(M, (Base.include(\$M, $(repr(path)))))
  """
end

"""
Exit hooks can get serialized and therefore end up in weird behaviour
When incrementally compiling
"""
function ExitHooksStart()
    """
    atexit_hook_copy = copy(Base.atexit_hooks) # make backup
    # clean state so that any package we use can carelessly call atexit
    empty!(Base.atexit_hooks)
    """
end

function ExitHooksEnd()
    """
    Base._atexit() # run all exit hooks we registered during precompile
    empty!(Base.atexit_hooks) # don't serialize the exit hooks we run + added
    # atexit_hook_copy should be empty, but who knows what base will do in the future
    append!(Base.atexit_hooks, atexit_hook_copy)
    """
end

"""
The command to pass to julia --output-o, that runs the julia code in `path` during compilation.
"""
function PrecompileCommand(path)
    ExitHooksStart() * InitBase() * InitREPL() * Include(path) * ExitHooksEnd()
end

"""
    compile_incremental(
        toml_path::String, snoopfile::String;
        force = false, precompile_file = nothing, verbose = true,
        debug = false, cc_flags = nothing
    )

    Extract all calls from `snoopfile` and ahead of time compiles them
    incrementally into the current system image.
    `force = true` will replace the old system image with the new one.
    The argument `toml_path` should contain a project file of the packages that `snoopfile` explicitly uses.
    Implicitly used packages & modules don't need to be contained!

    To compile just a single package, see the simpler version  `compile_incremental(package::Symbol)`:
"""
function compile_incremental(
        toml_path::Union{String, Nothing}, precompiles::String;
        force = false, verbose = true,
        debug = false, cc_flags = nothing
    )
    systemp = sysimg_folder("sys.a")
    sysout = sysimg_folder("sys.$(Libdl.dlext)")
    code = PrecompileCommand(precompiles)
    run_julia(
        code, O = 3, output_o = systemp, g = 1,
        track_allocation = "none", startup_file = "no", code_coverage = "none"
    )
    build_shared(sysout, systemp, false, sysimg_folder(), verbose, "3", debug, system_compiler, cc_flags)
    curr_syso = current_systemimage()
    force && cp(sysout, curr_syso, force = true)
    return sysout, curr_syso
end

"""
    compile_incremental(
        packages::Symbol...;
        force = false, reuse = false, verbose = true,
        debug = false, cc_flags = nothing,
        blacklist::Vector = []
    )

    Incrementally compile `package` into the current system image.
    `force = true` will replace the old system image with the new one.
    This process requires a script that julia will run in order to determine
    which functions to compile. A package may define a script called `snoopfile.jl`
    for this purpose. If this file cannot be found the package's test script
    `Package/test/runtests.jl` will be used. `compile_incremental` will search
    for `snoopfile.jl` in the package's root directory and in the folders
    `Package/src` and `Package/snoop`. For a more explicit version of compile_incremental,
    see: `compile_incremental(toml_path::String, snoopfile::String)`

    Not all packages can currently be compiled into the system image. By default,
    `compile_incremental(:Package) will also compile all of Package's dependencies.
    It can still be desirable to compile packages with dependencies that cannot be
    compiled. For this reason `compile_incremental` offers
    the ability for the user to pass a list of blacklisted packages
    that will be ignored during the compilation process. These are passed as a
    vector of package names (defined as either strings or symbols) using the
    `blacklist keyword argument`
"""
function compile_incremental(pkg::Symbol, packages::Symbol...;
                             blacklist::Vector=[], kw...)
    toml, precompile = snoop_packages(pkg, packages...; blacklist=blacklist)
    compile_incremental(toml, precompile; kw...)
end


text:
using MagneticReadHead
using Test

@testset "Should be able to add and remove breakpoints" begin
    @test length(set_breakpoint!(Test)) == 1
    @test length(rm_breakpoint!(Test)) == 0

    @test length(set_uninstrumented!(Test)) == 1
    @test length(rm_uninstrumented!(Test)) == 0
end


text:
@doc """
# Description
A function for optimising the GP hyperparameters based on type II maximum likelihood estimation. This function performs gradient based optimisation using the Optim pacakge to which the user is referred to for further details.

# Arguments:
* `gp::GP`: Predefined Gaussian process type
* `noise::Bool`: Noise hyperparameters should be optmized
* `mean::Bool`: Mean function hyperparameters should be optmized
* `kern::Bool`: Kernel function hyperparameters should be optmized
* `kwargs`: Keyword arguments for the optimize function from the Optim package

# Return:
* `::Optim.MultivariateOptimizationResults{Float64,1}`: optimization results object
""" ->
function optimize!(gp::GP; noise::Bool=true, mean::Bool=true, kern::Bool=true, method::Optim.Optimizer=ConjugateGradient(), kwargs...)
    function mll(hyp::Vector{Float64})
        set_params!(gp, hyp; noise=noise, mean=mean, kern=kern)
        update_mll!(gp)
        return -gp.mLL
    end
    function dmll!(hyp::Vector{Float64}, grad::Vector{Float64})
        set_params!(gp, hyp; noise=noise, mean=mean, kern=kern)
        update_mll_and_dmll!(gp; noise=noise, mean=mean, kern=kern)
        grad[:] = -gp.dmLL
    end
    function mll_and_dmll!(hyp::Vector{Float64}, grad::Vector{Float64})
        set_params!(gp, hyp; noise=noise, mean=mean, kern=kern)
        update_mll_and_dmll!(gp; noise=noise, mean=mean, kern=kern)
        grad[:] = -gp.dmLL
        return -gp.mLL
    end

    func = DifferentiableFunction(mll, dmll!, mll_and_dmll!)
    init = get_params(gp;  noise=noise, mean=mean, kern=kern)  # Initial hyperparameter values
    optimize(func, init, method, kwargs...)                    # Run optimizer
end


text:
###
### Basic MA operations
###

MA.mutability(::Type{<:DGAE}) = MA.IsMutable()

function MA.mutable_copy(expr::DGAE{C,V,P}) where {C,V,P}
    return DGAE{C,V,P}(MA.mutable_copy(expr.v), MA.mutable_copy(expr.p))
end

###
### MA.promote_operation
###

# Parameter -- Variable

function MA.promote_operation(
    ::Union{typeof(+),typeof(-)},
    ::Type{VariableRef},
    ::Type{ParameterRef},
)
    return PAE{Float64}
end

function MA.promote_operation(
    ::Union{typeof(+),typeof(-)},
    ::Type{ParameterRef},
    ::Type{VariableRef},
)
    return PAE{Float64}
end

# GAEp -- Variable

function MA.promote_operation(
    ::Union{typeof(+),typeof(-)},
    ::Type{VariableRef},
    ::Type{GAEp{C}},
) where {C}
    return PAE{C}
end

function MA.promote_operation(
    ::Union{typeof(+),typeof(-)},
    ::Type{GAEp{C}},
    ::Type{VariableRef},
) where {C}
    return PAE{C}
end

# GAEp -- GAE

function MA.promote_operation(
    op::Union{typeof(+),typeof(-)},
    ::Type{GAE{C1,ParameterRef}},
    ::Type{GAE{C2,V}},
) where {C1,C2,V}
    C = MA.promote_operation(op, C1, C2)
    return DGAE{C,V,ParameterRef}
end

function MA.promote_operation(
    op::Union{typeof(+),typeof(-)},
    ::Type{GAE{C1,V}},
    ::Type{GAE{C2,ParameterRef}},
) where {C1,C2,V}
    C = MA.promote_operation(op, C1, C2)
    return DGAE{C,V,ParameterRef}
end

# DGAE -- Number

function MA.promote_operation(
    op::Union{typeof(*),typeof(+),typeof(-)},
    ::Type{C1},
    ::Type{DoubleGenericAffExpr{C2,V,P}},
) where {C1<:Number,C2,V,P}
    C = MA.promote_operation(op, C1, C2)
    return DoubleGenericAffExpr{C,V,P}
end

function MA.promote_operation(
    op::Union{typeof(*),typeof(+),typeof(-)},
    ::Type{DoubleGenericAffExpr{C2,V,P}},
    ::Type{C1},
) where {C1<:Number,C2,V,P}
    C = MA.promote_operation(op, C1, C2)
    return DoubleGenericAffExpr{C,V,P}
end

# DGAE -- Variable

function MA.promote_operation(
    ::Union{typeof(+),typeof(-)},
    ::Type{VariableRef},
    ::Type{DGAE{C,VariableRef,P}},
) where {C,P}
    return DGAE{C,VariableRef,P}
end

function MA.promote_operation(
    ::Union{typeof(+),typeof(-)},
    ::Type{DGAE{C,VariableRef,P}},
    ::Type{VariableRef},
) where {C,P}
    return DGAE{C,VariableRef,P}
end

# DGAE -- Parameter

function MA.promote_operation(
    ::Union{typeof(+),typeof(-)},
    ::Type{ParameterRef},
    ::Type{DGAE{C,V,ParameterRef}},
) where {C,V}
    return DGAE{C,V,ParameterRef}
end

function MA.promote_operation(
    ::Union{typeof(+),typeof(-)},
    ::Type{DGAE{C,V,ParameterRef}},
    ::Type{ParameterRef},
) where {C,V}
    return DGAE{C,V,ParameterRef}
end

# DGAE -- GAE

function MA.promote_operation(
    op::Union{typeof(+),typeof(-)},
    ::Type{GAE{C1,V1}},
    ::Type{DGAE{C2,V1,V2}},
) where {C1,C2,V1,V2}
    C = MA.promote_operation(op, C1, C2)
    return DGAE{C,V1,V2}
end

function MA.promote_operation(
    op::Union{typeof(+),typeof(-)},
    ::Type{DGAE{C1,V1,V2}},
    ::Type{GAE{C2,V1}},
) where {C1,C2,V1,V2}
    C = MA.promote_operation(op, C1, C2)
    return DGAE{C,V1,V2}
end

function MA.promote_operation(
    op::Union{typeof(+),typeof(-)},
    ::Type{GAE{C1,V2}},
    ::Type{DGAE{C2,V1,V2}},
) where {C1,C2,V1,V2}
    C = MA.promote_operation(op, C1, C2)
    return DGAE{C,V1,V2}
end

function MA.promote_operation(
    op::Union{typeof(+),typeof(-)},
    ::Type{DGAE{C1,V1,V2}},
    ::Type{GAE{C2,V2}},
) where {C1,C2,V1,V2}
    C = MA.promote_operation(op, C1, C2)
    return DGAE{C,V1,V2}
end

# DGAE -- DGAE

function MA.promote_operation(
    op::Union{typeof(+),typeof(-)},
    ::Type{DGAE{C1,V1,V2}},
    ::Type{DGAE{C2,V1,V2}},
) where {C1,C2,V1,V2}
    C = MA.promote_operation(op, C1, C2)
    return DGAE{C,V1,V2}
end

###
### MA.operate!(zero/one, ...)
###

function MA.operate!(::typeof(zero), aff::DGAE)
    MA.operate!(zero, aff.v)
    MA.operate!(zero, aff.p)
    return aff
end

function MA.operate!(::typeof(one), aff::DGAE)
    MA.operate!(one, aff.v)
    MA.operate!(zero, aff.p)
    return aff
end

###
### MA.operate!(op::MA.AddSubMul, ...)
###

# 4-argument functions

# DGAE - V

function MA.operate!(
    op::MA.AddSubMul,
    aff::DGAE{C,V,P},
    c::Number,
    x::Union{V,GAE{C,V}},
) where {C,V,P}
    if !iszero(c)
        MA.operate!(op, aff.v, c, x)
    end
    return aff
end

function MA.operate!(
    op::MA.AddSubMul,
    aff::DGAE{C,V,P},
    x::Union{V,GAE{C,V}},
    c::Number = 1,
) where {C,V,P}
    return MA.operate!(op, aff, c, x)
end

# DGAE - P

function MA.operate!(
    op::MA.AddSubMul,
    aff::DGAE{C,V,P},
    c::Number,
    x::Union{P,GAE{C,P}},
) where {C,V,P}
    if !iszero(c)
        MA.operate!(op, aff.p, c, x)
    end
    return aff
end

function MA.operate!(
    op::MA.AddSubMul,
    aff::DGAE{C,V,P},
    x::Union{P,GAE{C,P}},
    c::Number = 1,
) where {C,V,P}
    return MA.operate!(op, aff, c, x)
end

# DGAE - DGAE

function MA.operate!(
    op::MA.AddSubMul,
    aff::DGAE{C,V,P},
    c::Number,
    x::DGAE{C,V,P},
) where {C,V,P}
    if !iszero(c)
        MA.operate!(op, aff.p, c, x.p)
        MA.operate!(op, aff.v, c, x.v)
    end
    return aff
end

function MA.operate!(
    op::MA.AddSubMul,
    aff::DGAE{C,V,P},
    x::DGAE{C,V,P},
    c::Number = 1,
) where {C,V,P}
    return MA.operate!(op, aff, c, x)
end

# DGAE - Number

function MA.operate!(
    op::MA.AddSubMul,
    aff::DGAE,
    c::Number,
    x::Number = 1,
)
    if !iszero(c) && !iszero(x)
        MA.operate!(op, aff.v, c, x)
    end
    return aff
end

# n-argument functions

@generated function _add_sub_mul_reorder!(
    op::MA.AddSubMul,
    expr::DGAE,
    args::Vararg{Any,N},
) where {N}
    n = length(args)
    @assert n ≥ 3
    varidx = findall(t -> !(t <: Number), collect(args))
    allscalar = all(t -> (t <: Number), args[setdiff(1:n, varidx)])
    # We need to get down to only two factors
    # If there are only constants and one JuMP expressions, then we multiply
    # the constants together. Otherwise we multiply all factors except the
    # last one, there may be a better thing to do here.
    idx = (allscalar && length(varidx) == 1) ? varidx[1] : n
    coef = Expr(:call, :*, [:(args[$i]) for i in setdiff(1:n, idx)]...)
    return :(MA.operate!(op, expr, $coef, args[$idx]))
end

function MA.operate!(
    op::MA.AddSubMul,
    expr::DGAE,
    x,
    y,
    z,
    other_args::Vararg{Any,N},
) where {N}
    return _add_sub_mul_reorder!(op, expr, x, y, z, other_args...)
end

function MA.operate!(::typeof(*), expr::DGAE, x::Number)
    MA.operate!(*, expr.v, x)
    MA.operate!(*, expr.p, x)
    return expr
end

function MA.operate!(::typeof(+), expr::DGAE, x)
    return JuMP.add_to_expression!(expr, x)
end

function MA.operate!(::typeof(-), expr::DGAE, x)
    return JuMP.add_to_expression!(expr, -1, x)
end

###
### JuMP.add_to_expression
###

# 2-argument functions

function JuMP.add_to_expression!(
    lhs::DGAE{C1,V,P},
    rhs::Union{Number,V,GAE{C2,V}},
) where {C1,C2,V,P}
    JuMP.add_to_expression!(lhs.v, rhs)
    return lhs
end

function JuMP.add_to_expression!(
    lhs::DGAE{C1,V,P},
    rhs::Union{P,GAE{C2,P}},
) where {C1,C2,V,P}
    JuMP.add_to_expression!(lhs.p, rhs)
    return lhs
end

function JuMP.add_to_expression!(lhs::DGAE, rhs::DGAE)
    JuMP.add_to_expression!(lhs.p, rhs.p)
    JuMP.add_to_expression!(lhs.v, rhs.v)
    return lhs
end

# 3-argument functions

# V -- Number

function JuMP.add_to_expression!(
    lhs::DGAE{C1,V,P},
    x::Union{Number,V,GAE{C2,V}},
    y::Number,
) where {C1,C2,V,P}
    if !iszero(y)
        add_to_expression!(lhs.v, x, y)
    end
    return lhs
end

function JuMP.add_to_expression!(
    lhs::DGAE{C1,V,P},
    y::Number,
    x::Union{V,GAE{C2,V}},
) where {C1,C2,V,P}
    return add_to_expression!(lhs, x, y)
end

# P -- Number

function JuMP.add_to_expression!(
    lhs::DGAE{C1,V,P},
    x::Union{P,GAE{C2,P}},
    y::Number,
) where {C1,C2,V,P}
    if !iszero(y)
        JuMP.add_to_expression!(lhs.p, x, y)
    end
    return lhs
end

function JuMP.add_to_expression!(
    lhs::DGAE{C1,V,P},
    y::Number,
    x::Union{P,GAE{C2,P}},
) where {C1,C2,V,P}
    return add_to_expression!(lhs, x, y)
end

# DGAE -- Real

function JuMP.add_to_expression!(lhs::DGAE, x::DGAE, y::Real)
    if iszero(y)
        return lhs
    end
    if !iszero(x.v)
        JuMP.add_to_expression!(lhs.v, x.v, y)
    end
    if !iszero(x.p)
        JuMP.add_to_expression!(lhs.p, x.p, y)
    end
    return lhs
end

function JuMP.add_to_expression!(lhs::DGAE, y::Real, x::DGAE)
    return add_to_expression!(lhs, x, y)
end


text:
using Gadfly
using RDatasets

plot(data("datasets", "iris"),
     x="Sepal.Length", y="Sepal.Width", color="Species",
     Geom.point)