text:
using LazyArrays, FillArrays, ArrayLayouts, StaticArrays, Test
import LazyArrays: CachedArray, CachedMatrix, CachedVector, PaddedLayout, CachedLayout, resizedata!, zero!,
                    CachedAbstractArray, CachedAbstractVector, CachedAbstractMatrix

@testset "Cache" begin
    @testset "basics" begin
        A = 1:10
        C = cache(A)
        @test size(C) == (10,)
        @test axes(C) == (Base.OneTo(10),)
        @test all(Vector(C) .=== Vector(A))
        @test cache(C) isa CachedArray{Int,1,Vector{Int},UnitRange{Int}}
        C2 = cache(C)
        @test C2.data !== C.data
        @test C[:] == C[Base.Slice(Base.OneTo(10))] == C[1:10] == C[collect(1:10)] == C
        @test C[CartesianIndex(2)] == 2

        @test cache(A)[2,1] == 2
        @test_throws BoundsError cache(A)[2,2]

        A = reshape(1:10^2, 10,10)
        C = cache(A)
        @test size(C) == (10,10)
        @test axes(C) == (Base.OneTo(10),Base.OneTo(10))
        @test all(Array(C) .=== Array(A))
        @test cache(C) isa CachedArray{Int,2,Matrix{Int},typeof(A)}
        C2 = cache(C)
        @test C2.data !== C.data
        @test C[1:2,1:2] == [1 11; 2 12]

        @test C[1,:] == C[1,1:10] == C[1,collect(1:10)] == A[1,1:10]
        @test C[:,1] == C[1:10,1] == C[collect(1:10),1] == A[1:10,1]
        @test C[1:10,1:10] == C[collect(1:10),collect(1:10)] == C[:,:] == C[1:10,:] ==
                C[collect(1:10),:] == C[:,1:10] == C[:,collect(1:10)] == A

        A = reshape(1:10^3, 10,10,10)
        C = cache(A)
        @test size(C) == (10,10,10)
        @test axes(C) == (Base.OneTo(10),Base.OneTo(10),Base.OneTo(10))
        @test all(Array(C) .=== Array(A))

        A = reshape(1:10^3, 10,10,10)
        C = cache(A)
        LazyArrays.resizedata!(C,5,5,5)
        LazyArrays.resizedata!(C,8,8,8)
        @test all(C.data .=== Array(A)[1:8,1:8,1:8])

        @test C[1:3,1,:] == A[1:3,1,:]
    end

    @testset "convert" begin
        C = cache(1:5);
        C[1] = 3;
        @test convert(AbstractVector{Int}, C) ≡ convert(AbstractArray{Int}, C) ≡ C
        @test convert(AbstractVector{Float64}, C) isa CachedVector
        @test convert(AbstractArray{Float64}, C) isa CachedVector
        @test convert(AbstractVector{Float64}, C) == convert(AbstractArray{Float64}, C) == C
    end

    @testset "Matrix cache" begin
        A = collect(1:5)
        C = cache(A)
        @test C isa Vector{Int}
        C[1] = 2
        @test A[1] ≠ 2
        @test resizedata!(C,6) === C

        A = cache(Matrix(reshape(1:6,2,3)));
        C = cache(A);
        @test C isa Matrix{Int}
        C[1,1] = 2
        @test A[1,1] ≠ 2
        C[1] = 3
        @test A[1,1] ≠ 3
        @test resizedata!(C,2,3) === C
    end

    @testset "setindex!" begin
        A = (1:5)
        C = cache(A)
        C[1] = 2
        @test C[1] == 2

        A = cache(reshape(1:6,2,3))
        C = cache(A)
        C[1,1] = 2
        @test A[1,1] ≠ 2
        @test C[1,1] == 2
        C[1] = 3
        @test C[1,1] == 3

        @test_throws BoundsError C[3,1]
        @test_throws BoundsError C[7]
    end

    @testset "colsupport past size" begin
        C = cache(Zeros(5,5)); C[5,1];
        @test colsupport(C,1) == Base.OneTo(5)
        @test colsupport(C,3) == 1:0
        @test rowsupport(C,1) == Base.OneTo(1)
    end

    @testset "broadcast" begin
        x = CachedArray([1,2,3],1:8);
        y = 1:8;
        f = (x,y) -> cos(x*y)
        @test f.(x,y) isa CachedArray
        @test @inferred(broadcast(f,x,y)) == f.(Vector(x), Vector(y))

        @test (x + y) isa CachedArray
        @test (x + y).array isa AbstractRange
        @test (x + y) == Vector(x) + Vector(y)

        @testset "Padded" begin
            z = CachedArray([1,4],Zeros{Int}(8));
            @test (x .+ z) isa CachedArray
            @test (x + z) isa CachedArray
            @test Vector( x .+ z) == Vector( x + z) == Vector(x) + Vector(z)
        end

        @testset "Lazy mixed with Static treats as Lazy" begin
            s = SVector(1,2,3,4,5,6,7,8)
            @test f.(x , s) isa CachedArray
            @test f.(x , s) == f.(Vector(x), Vector(s))
        end

        @testset "sub-matrix" begin
            A = cache(Zeros(8,8));
            b = BroadcastVector(exp,randn(8))
            @test view(A, :, 1) .+ x == x .+ view(A, :, 1) == A[:,1] .+ x
            @test view(A, :, 1) .+ b == b .+ view(A, :, 1) == A[:,1] .+ b

            A = cache(Zeros(8,8));
            view(A, :, 2)[3] = 4
            @test A[3,2] == 4
        end
    end

    @testset "padded CachedVector getindex" begin
        v = CachedArray([1,2,3],Zeros{Int}(1000))
        @test v[3:100] == [3; zeros(Int,97)]
        @test v[3:100] isa CachedArray{Int,1,Vector{Int},Zeros{Int,1,Tuple{Base.OneTo{Int}}}}
        @test v[4:end] isa CachedArray{Int,1,Vector{Int},Zeros{Int,1,Tuple{Base.OneTo{Int}}}}
        @test MemoryLayout(v) isa PaddedLayout
        @test all(iszero,v[4:end])
        @test isempty(v[4:0])
        @test norm(v) == norm(Array(v)) == norm(v,2)
        @test norm(v,1) == norm(Array(v),1)
        @test norm(v,Inf) == norm(Array(v),Inf)
        @test norm(v,3) == norm(Array(v),3)
        v = CachedArray([1,2,3],Fill{Int}(1,1000))
        @test v[3:100] == [3; ones(97)]
        @test norm(v) == norm(Array(v)) == norm(v,2)
        @test norm(v,1) == norm(Array(v),1)
        @test norm(v,Inf) == norm(Array(v),Inf)
        @test norm(v,3) ≈ norm(Array(v),3)
    end

    @testset "ambiguity broadcast" begin
        c = cache(1:100)
        v = Vcat([1,2,3],0:96)
        z = Zeros(100)
        @test v .+ c == c .+ v == Array(c) + Array(v)
        @test z .+ c == c .+ z == Array(c)
        @test c .* z ≡ z .* c ≡ z
    end

    @testset "Fill" begin
        c = cache(Fill(1,10))
        @test c[2:3] isa LazyArrays.CachedVector{Int,Vector{Int},Fill{Int,1,Tuple{Base.OneTo{Int}}}}
        @test c[[2,4,6]] isa LazyArrays.CachedVector{Int,Vector{Int},Fill{Int,1,Tuple{Base.OneTo{Int}}}}
    end

    @testset "linalg" begin
        c = cache(Fill(3,3,3))
        @test fill(2,1,3) * c == fill(18,1,3)
        @test ApplyMatrix(exp,fill(3,3,3)) * c == exp(fill(3,3,3)) * fill(3,3,3)
        @test BroadcastMatrix(exp,fill(3,3,3)) * c == exp.(fill(3,3,3)) * fill(3,3,3)
        @test fill(2,3)' * c == fill(18,1,3)
        @test fill(2,3,1)' * c == fill(18,1,3)
    end

    @testset "broadcast" begin
        c = cache(Fill(3.0,10)); c[1] = 2;
        @test exp.(c) isa typeof(c)
        @test exp.(c) == exp.(Vector(c))
        @test c .+ 1 isa typeof(c)
        @test c .+ 1 == Vector(c) .+ 1
        @test 1 .+ c isa typeof(c)
        @test 1 .+ c == 1 .+ Vector(c)
        @test c .+ Ref(1) isa typeof(c)
        @test Ref(1) .+ c isa typeof(c)
    end

    @testset "*" begin
        A = randn(2,2)
        B = ApplyMatrix(exp,A)
        C = BroadcastMatrix(exp,A)
        D = Diagonal(randn(2))
        x = cache(Fill(3,2))
        @test A*x ≈ A*Vector(x)
        @test B*x ≈ Matrix(B)*Vector(x)
        @test C*x ≈ Matrix(C)*Vector(x)
        @test D*x ≈ Matrix(D)*Vector(x)
        @test A'x ≈ Matrix(A)'Vector(x)
        @test B'x ≈ Matrix(B)'Vector(x)
        @test C'x ≈ Matrix(C)'Vector(x)

        @testset "padded" begin
            z = cache(Zeros(2));
            @test MemoryLayout(z) isa PaddedLayout
            @test A*z ≈ A*Vector(z)
            @test B*z ≈ Matrix(B)*Vector(z)
            @test C*z ≈ Matrix(C)*Vector(z)
            @test D*z ≈ Matrix(D)*Vector(z)
            @test A'z ≈ Matrix(A)'Vector(z)
            @test B'z ≈ Matrix(B)'Vector(z)
            @test C'z ≈ Matrix(C)'Vector(z)

            p = Vcat([1,2],Zeros(3)) + cache(Zeros(5))
            @test MemoryLayout(p) isa PaddedLayout
            @test p isa CachedVector{Float64,Vector{Float64},<:Zeros}
            @test p == [1; 2; zeros(3)]
        end
    end

    @testset "copyto!" begin
        a = CachedArray([1,2,3], Zeros{Int}(8));
        b = CachedArray(Int[], Zeros{Int}(8));
        c = CachedArray(Float64[], Zeros{Float64}(8));
        @test copyto!(b, a) == a == b
        @test copyto!(c, a) == a == c

        @test copyto!(a, Zeros{Int}(8)) == zeros(8)
        a = CachedArray([1,2,3], Zeros{Int}(8));
        copyto!(view(a,3:8), Zeros{Int}(6))
        @test a == [1; 2; zeros(6)]

        a = CachedArray([3,missing], Zeros{Union{Int,Missing}}(4));
        b = CachedArray(Union{Int,Missing}[], Zeros{Union{Int,Missing}}(4));
        @test all(copyto!(b, a) .=== a .=== b)

        a = CachedArray([1,2,3], Zeros{Int}(8));
        copyto!(view(a,1:2), [4,5])
        @test a[1:2] == [4,5]
        a = CachedArray(Int[], Zeros{Int}(8));
        copyto!(view(a,1:2), view(Vcat([4,5],Zeros(6)),1:2))
    end

    @testset "fill!/lmul!/rmul!" begin
        a = CachedArray(Zeros{Float64}(100_000_000));
        a[3] = 4;
        fill!(a, 0.0);
        @test a.datasize[1] == 3
        @test a[1:3] == zeros(3)
        a[3] = 4;
        rmul!(a, 2.0);
        @test a.datasize[1] == 3
        @test a[1:3] == [0,0,8]
        lmul!(2.0, a);
        @test a.datasize[1] == 3
        @test a[1:3] == [0,0,16]

        @test_throws ArgumentError fill!(a, 1.0)
        @test_throws ArgumentError rmul!(a, Inf)
        @test_throws ArgumentError lmul!(Inf, a)

        a = CachedArray(Ones{Float64}(100_000_000));
        lmul!(5, view(a, 1:3))
        @test a[1:5] == [fill(5,3); 1; 1]
        rmul!(view(a, 4:5), 6)
        @test a[1:7] == [fill(5,3); fill(6,2); fill(1,2)]
    end

    @testset "Padded broadcast" begin
        a = CachedArray([1,2,3], Zeros{Int}(8));
        r = a .- a;
        @test MemoryLayout(r) isa PaddedLayout
        @test r.datasize[1] == 3
        @test r == Vector(a) - Vector(a)

        a = CachedArray([1,2,3], Zeros{Int}(8));
        b = Fill(2,8);
        r = a .- b;
        @test MemoryLayout(r) isa CachedLayout{DenseColumnMajor,FillLayout}
        @test r.datasize[1] == 3
        @test r == Vector(a) - Vector(b)

        a = CachedArray([1,2,3], Zeros{Int}(8));
        b = Fill(2,8);
        r = b .- a;
        @test MemoryLayout(r) isa CachedLayout{DenseColumnMajor,FillLayout}
        @test r.datasize[1] == 3
        @test r == Vector(b) - Vector(a)

        a = CachedArray([1,2,3], Zeros{Int}(8));
        b = CachedArray([1,2],Fill(2,8));
        r = a .- b;
        @test MemoryLayout(r) isa CachedLayout{DenseColumnMajor,FillLayout}
        @test r.datasize[1] == 3
        @test r == Vector(a) - Vector(b)

        b = CachedArray([1,2,3], Zeros{Int}(8));
        a = CachedArray([1,2],Fill(2,8));
        r = a .- b;
        @test MemoryLayout(r) isa CachedLayout{DenseColumnMajor,FillLayout}
        @test r.datasize[1] == 3
        @test r == Vector(a) - Vector(b)

        a = CachedArray([1,2,3], Zeros{Int}(8));
        @test_throws DimensionMismatch a .+ CachedArray([1,2],Fill(2,6))
    end

    @testset "BigFloat" begin
        a = cache(Zeros{BigFloat}(10));
        resizedata!(a,3);
        resizedata!(a,10);
        @test a == Zeros{BigFloat}(10)
    end

    @testset "zero!" begin
        a = cache(Zeros(10));
        @test zero!(a) ≡ a
        @test a.datasize == (0,)
    end

    @testset "minimum/maximum" begin
        c = cache(Fill(2,4));
        c[1] = 1;
        @test maximum(c) == 2
        @test minimum(c) == 1
        @test sum(c) == 7
        c[1] = 3;
        @test maximum(c) == 3
        @test minimum(c) == 2
        c[1:4] .= 1;
        @test maximum(c) == 1
        c[1:4] .= 3;
        @test minimum(c) == 3
    end

    @testset "permutedims" begin
        a = cache(Zeros(10));
        a[3] = 2
        a.data[4:end] .= NaN # fill with noise
        @test permutedims(a) isa CachedMatrix
        @test permutedims(a).datasize == (1,3)
        @test permutedims(a)[4] == permutedims(a)[1,4] == 0.0
        @test permutedims(a) == [0 0 2 zeros(1,7)]

        A = cache(Zeros(5,6));
        A[2,3] = 2;
        @test permutedims(A) isa CachedMatrix
        @test permutedims(A).datasize == (3,2)
    end

    @testset "CachedAbstract" begin
        A = CachedAbstractArray(1:5);
        A[2] = 3
        @test -A isa CachedAbstractVector
        @test A .+ 1 isa CachedAbstractVector
        @test A .+ Ref(1) isa CachedAbstractVector
        @test 1 .+ A isa CachedAbstractVector
        @test Ref(1) .+ A isa CachedAbstractVector
        @test A .+ A isa CachedAbstractVector

        @test CachedAbstractVector(1:5) isa CachedAbstractVector
        @test CachedAbstractMatrix(reshape(1:6,2,3)) isa CachedAbstractMatrix
    end
end


text:
# This file provides versions of solvers based on communication graphs

using Statistics
using StatsBase
using LinearAlgebra

# Methods for evaluating communication statistics
#
# WARNING: Call after solving a problem
#
# * communication_span: "time" span for communications.
# * communication_messages: Total number of messages sent during execution
#   including each hop
# * communication_volume: Sum of messages and hops multiplied by the number of
#   decisions in each message

###############################
# Adjacenty matrix and plotting
###############################

function make_adjacency_matrix(problem::ExplicitPartitionProblem, range)
  num_agents = get_num_agents(problem)

  adjacency = zeros(Int64, num_agents, num_agents)

  for ii in 1:num_agents, jj in 1:ii-1
    center_ii = get_center(get_agent(problem, ii))
    center_jj = get_center(get_agent(problem, jj))

    if norm(center_ii - center_jj) <= range
      adjacency[ii,jj] = 1
      adjacency[jj,ii] = 1
    end
  end

  adjacency
end

# Ugly lines to propagate the matrix a few hops
function make_hop_adjacency(adjacency::Matrix{Int64}, hops::Int64)
  mat = (adjacency + I)^hops
  map(x -> (x > 0 ? 1 : 0), mat)
  for ii = 1:size(mat, 1)
    mat[ii, ii] = 0
  end

  mat
end

neighbors(adjacency::Matrix{<:Integer}, index::Integer) =
  findall(x -> x > 0, adjacency[index,:])

# Ugly implementation of Dijkstra
function shortest_path(adjacency::Matrix, a::Integer, b::Integer)
  dist = fill(Inf, size(adjacency, 1))
  dist[a] = 0
  prev = fill(0, size(adjacency, 1))

  inds = collect(1:size(adjacency, 1))

  sort_inds() = sort!(inds, by=x->dist[x], rev = true)

  while !isempty(inds)
    sort_inds()
    least = pop!(inds)

    if least == b
      break
    end

    # Update
    for index in neighbors(adjacency, least)
      new_dist = dist[least] + 1
      if dist[index] > new_dist
        dist[index] = new_dist
        prev[index] = least
      end
    end
  end

  path = [b]
  while prev[path[end]] != 0
    push!(path, prev[path[end]])
  end

  (dist[b], path)
end

path_distance(x...) = shortest_path(x...)[1]

# Check connectivity via matrix exponentiation
function is_connected(adjacency::Matrix)
  paths = (adjacency + I) ^ (size(adjacency, 1) - 1)

  all(x -> x > 0, paths)
end

function sample_circle()
  r = sqrt(rand())
  theta = 2*pi*rand()
  r * [cos(theta);sin(theta)]
end

# We need the problem to be connected to avoid bad things
function generate_connected_problem(agent_specification, num_agents;
                                    communication_radius,
                                    objective)
  # Sample the agents
  agents = [make_agent(agent_specification)]
  for ii in 1:num_agents-1
    seed = get_center(sample(agents))

    # sample centers until one is in the unit square
    while true
      center = seed + communication_radius * sample_circle()

      if all(x -> 0 <= x <= 1, center)
        push!(agents, make_agent(agent_specification, center))
        break
      end
    end
  end

  problem = ExplicitPartitionProblem(objective, agents)
end

function plot_adjacency(problem::ExplicitPartitionProblem, range::Real)
  adjacency = make_adjacency_matrix(problem, range)

  plot_adjacency(problem, adjacency)
end
function plot_adjacency(problem::ExplicitPartitionProblem, adjacency::Array)
  num_agents = get_num_agents(problem)

  for ii in 1:num_agents, jj in 1:ii-1
    center_ii = get_center(problem, ii)
    center_jj = get_center(problem, jj)

    if adjacency[ii,jj] > 0
      x = map(first, [center_ii, center_jj])
      y = map(last, [center_ii, center_jj])
      plot(x, y, color="black")
    end
  end
end

# Plots the shortest path.
# * Note that ret[1].remove() will remove the path
function plot_shortest_path(problem::ExplicitPartitionProblem, adjacency, a, b)
  (dist, path) = shortest_path(adjacency, a, b)

  x = Float64[]
  y = Float64[]
  for ii in path
    center = get_center(problem, ii)

    push!(x , first(center))
    push!(y, last(center))
  end

  plot(x, y, color="red")
end

##################
# Multi-hop solver
##################

struct MultiHopSolver <: DAGSolver
  adjacency::Array
  hop_adjacency::Array
  hops::Integer
  nominal_solver::DAGSolver

  function MultiHopSolver(problem::PartitionProblem;
                          solver::DAGSolver,
                          communication_range,
                          num_hops)

    adjacency = make_adjacency_matrix(problem, communication_range)

    hop_adjacency = make_hop_adjacency(adjacency, num_hops)

    new(adjacency, hop_adjacency, num_hops, solver)
  end
end

sequence(x::MultiHopSolver) = sequence(x.nominal_solver)
partitions(x::MultiHopSolver) = partitions(x.nominal_solver)

# Keep only neighbors that are within the communication range
# Note: in-neighbors for the DAG solver
function in_neighbors(x::MultiHopSolver, agent_index)
  nominal_neighbors = in_neighbors(x.nominal_solver, agent_index)

  neighbors = Int64[]
  for neighbor in nominal_neighbors
    if x.hop_adjacency[agent_index, neighbor] > 0
      push!(neighbors, neighbor)
    end
  end

  neighbors
end

# Multi-hop communications

solver_rank(x::MultiHopSolver) = solver_rank(x.nominal_solver)
communication_span(x::MultiHopSolver) =
  x.hops * (communication_span(x.nominal_solver) - 1)
# Number of messages multiplied by the number of hops the message travels
function communication_messages(x::MultiHopSolver)
  num_agents = solver_rank(x)

  sum(1:num_agents) do a
    mapreduce(+, in_neighbors(x, a), init=0) do b
      path_distance(x.adjacency, a, b)
    end
  end
end
# Robots send a single decision at a time
communication_volume(x::MultiHopSolver) = communication_messages(x)

#################################################
# Sequential solver with communication statistics
#################################################

struct SequentialCommunicationSolver
  adjacency::Array
end

function SequentialCommunicationSolver(problem::PartitionProblem,
                                       communication_range)
  adjacency = make_adjacency_matrix(problem, communication_range)

  SequentialCommunicationSolver(adjacency)
end

solver_rank(x::SequentialCommunicationSolver) = size(x.adjacency, 1)

# Simply solve using the standard sequential solver
function solve_problem(::SequentialCommunicationSolver,
                       p::PartitionProblem; kwargs...)
  solve_sequential(p)
end
# Each robot except the first receives a single message from the previous
# containing all prior decisions. Include path lengths
function communication_messages(x::SequentialCommunicationSolver)
  num_agents = solver_rank(x)

  sum(1:num_agents-1) do index
    path_distance(x.adjacency, index, index + 1)
  end
end
# Sum of message hops times message size
# The message size is equal to the index of the sender
function communication_volume(x::SequentialCommunicationSolver)
  num_agents = solver_rank(x)

  sum(1:num_agents-1) do index
    index * path_distance(x.adjacency, index, index + 1)
  end
end
# Span is the distance from each agent to the next
# This will be the same as the communication_messages
function communication_span(x::SequentialCommunicationSolver)
  num_agents = solver_rank(x)

  sum(1:num_agents - 1) do index
    path_distance(x.adjacency, index, index + 1)
  end
end

################
# Auction solver
################

# Based on the distributed planner in
# by Choi et al. Transactions on Robots, 2009
# as well as the work of
# Luo et al. IROS, 2016
#
# Modifications:
# * Each agent has its own set of actions as for the definition of the
#   PartitionProblem type
# * Messages are maps from ids to values rather than vectors containing all
#   assignments

# Abstract auctions type for brevity when defining the local auctions type
abstract type AbstractAuctionSolver end
abstract type AbstractAuctionAgent end

# We will evaluate message counts and such
mutable struct AuctionSolver <: AbstractAuctionSolver
  adjacency::Matrix

  # We will allow for convergence before reaching the maximum number of steps
  # (equivalent span plus 1)
  nominal_steps::Integer

  span::Integer
  messages::Integer
  volume::Integer

  solved::Bool

  AuctionSolver(adjacency, span) = new(adjacency, span, 0, 0, 0, false)
end

# Second term is the upper bound on convergence time
# (longest possible path times number of assignments)
AuctionSolver(adjacency) = AuctionSolver(adjacency, size(adjacency, 1)^2)

function AuctionSolver(problem::PartitionProblem, communication_range, x...)
  adjacency = make_adjacency_matrix(problem, communication_range)

  AuctionSolver(adjacency, x...)
end

# Contains details for a given agent (its assignments and index)
struct AuctionAgent <: AbstractAuctionAgent
  index::Integer
  assignments::Vector{ExplicitSolutionElement}
end

AuctionAgent(index::Integer) = AuctionAgent(index, ExplicitSolutionElement[])
# Constructor for incrementing the assignments
AuctionAgent(a::AuctionAgent, s::Vector) = AuctionAgent(a.index, s)

get_index(x::AbstractAuctionAgent) = x.index
# Returns the agent's current list of assignments to all agents
get_assignments(x::AbstractAuctionAgent) = x.assignments
# Returns the assignment to this agent
get_assignment(agent::AbstractAuctionAgent) =
  first(filter(x->first(x) == get_index(agent), get_assignments(agent)))

num_assignments(agent::AbstractAuctionAgent) = length(get_assignments(agent))

neighbors(adjacency::Matrix, a::AbstractAuctionAgent) =
  neighbors(adjacency, a.index)

# Get the set of neighbors by adjacency and filtering
function neighbors(adjacency::Matrix, agents::Vector{<:AbstractAuctionAgent},
                   a::AbstractAuctionAgent)
  ns = neighbors(adjacency, a)

  filter(x -> in(x.index,ns), agents)
end

is_solved(x::AbstractAuctionSolver) = x.solved
assert_solved(x::AbstractAuctionSolver) =
  is_solved(x) ? true : error("Solve problem before calling")

communication_span(x::AbstractAuctionSolver) = assert_solved(x) && x.span
communication_messages(x::AbstractAuctionSolver) =
  assert_solved(x) && x.messages
communication_volume(x::AbstractAuctionSolver) = assert_solved(x) && x.volume

function converged(agents::Vector{<:AbstractAuctionAgent})
  # As a proxy for requiring maximal assignments, use the fact that the rank is
  # equal to the number of agents
  ( all(x -> length(x.assignments) == length(agents), agents)
  # Assert that all agents have the same assignments
  && all(x -> Set(x.assignments) == Set(first(agents).assignments), agents))
end

# Priority goes to agents with lower indices
priority(a::Int64, b::Int64) = a < b
priority(a::Tuple{Int64,Int64}, b::Tuple{Int64,Int64}) =
  priority(first(a), first(b)) ||
  (first(a) == first(b) && priority(last(a), last(b)))

# Run general greedy assignment
#
# Being careful to
# * Break ties deterministically, even within proposed assignments
#   (to ensure convergence)
# * Track independence
function greedy_assignment(problem::ExplicitPartitionProblem,
                           agent::AuctionAgent,
                           assignments::Vector;
                           prior_selections = ExplicitSolutionElement[]
                          )
  block = get_element_indices(problem, get_index(agent))

  # Construct a local ground with the current block and incoming assignments
  G = convert(Vector{ExplicitSolutionElement}, union(block,
                                                     get_assignments(agent),
                                                     assignments...))

  # Assignments
  X = prior_selections

  filter!(x->can_augment(X, x), G)

  # Greedy selection process
  while !isempty(G)
    values = map(G) do x
      marginal_gain(y->objective(problem, y), x, X)
    end

    max_value = maximum(values)
    maximizers = G[findall(x->x==max_value, values)]

    # Priority goes to least index
    max_priority = first(sort(maximizers; lt=priority))

    # Conditionally add the new solution element
    if independent(union(X, [max_priority]))
      X = union(X, [max_priority])
    end

    # Element either assigned or cannot be added to this or future solutions
    setdiff!(G, [max_priority])

    # Remove elements that cannot be added to the solution
    filter!(x->can_augment(X, x), G)
  end

  X
end

function solve_problem(solver::AuctionSolver, problem::PartitionProblem;
                       kwargs...)
  num_agents = get_num_agents(problem)
  adjacency = solver.adjacency
  agents = map(x->AuctionAgent(x), 1:num_agents)

  for ii in 1:solver.nominal_steps
    # Update span
    if ii != 1
      # Skip the first round (which we can complete without communication)
      solver.span += 1
    end

    updated_agents = map(agents) do agent
      ns = neighbors(adjacency, agents, agent)

      assignments = map(get_assignments, ns)

      # Update messages and volume
      if ii != 1
        solver.messages += length(ns)
        solver.volume += mapreduce(num_assignments, +, ns, init=0)
      end

      # Identify decisions that have reached consensus
      least_difference = minimum(ns) do neighbor
        first_difference(agent, neighbor, include_unassigned=true)
      end
      prior_selections = get_assignments(agent)[1:least_difference-1]

      # Construct the new set of assigments
      new_assignments = greedy_assignment(problem, agent, assignments,
                                          prior_selections=prior_selections)

      AuctionAgent(agent, new_assignments)
    end

    agents[:] = updated_agents

    if converged(agents)
      break
    end
  end

  solver.solved = true

  assignments = map(get_assignment, agents)

  # Return the solution
  evaluate_solution(problem, assignments)
end

######################
# Local auction solver
######################

# This is a similar auction solver but does not assume that agents can
# re-evaluate others' objective values

# We will evaluate message counts and such
mutable struct LocalAuctionSolver <: AbstractAuctionSolver
  adjacency::Matrix

  # We will allow for convergence before reaching the maximum number of steps
  # (equivalent span plus 1)
  nominal_steps::Integer

  span::Integer
  messages::Integer
  volume::Integer

  solved::Bool

  LocalAuctionSolver(adjacency, span) = new(adjacency, span, 0, 0, 0, false)
end

# Second term is the upper bound on convergence time
# (longest possible path times number of assignments)
LocalAuctionSolver(adjacency) =
  LocalAuctionSolver(adjacency, size(adjacency, 1)^2)

function LocalAuctionSolver(problem::PartitionProblem, communication_range,
                            x...)
  adjacency = make_adjacency_matrix(problem, communication_range)

  LocalAuctionSolver(adjacency, x...)
end

# Contains details for a given agent (its assignments and index)
struct LocalAuctionAgent <: AbstractAuctionAgent
  index::Integer
  assignments::Vector{ExplicitSolutionElement}
  values::Vector{Float64}
end

LocalAuctionAgent(index::Integer) =
  LocalAuctionAgent(index, ExplicitSolutionElement[], Float64[])
# Constructor for incrementing the assignments
LocalAuctionAgent(a::LocalAuctionAgent, s::Vector{ExplicitSolutionElement},
                  v::Vector{Float64}) = LocalAuctionAgent(a.index, s, v)

get_values(x::LocalAuctionAgent) = x.values

get_value_assignment(x::LocalAuctionAgent, index::Integer) =
  (get_values(x)[index], get_assignments(x)[index])

# Returns the first difference between a and b
#
# If including unassigned values, assignments past the length of the shorter are
# considered. If both arrays are the same up to that point, the first difference
# is at length(shorter) + 1
#
# Else return nothing if there is no difference up to the length of the shorter
function first_difference(a::AbstractAuctionAgent, b::AbstractAuctionAgent;
                          include_unassigned = false
                         )
  len = min(num_assignments(a), num_assignments(b))

  index = findfirst(x-> x[1] != x[2],
                    collect(zip(get_assignments(a)[1:len],
                                get_assignments(b)[1:len])))

  if isnothing(index) && include_unassigned
    min(num_assignments(a), num_assignments(b)) + 1
  else
    index
  end
end

function dominates(a::Tuple{Float64,ExplicitSolutionElement},
                   b::Tuple{Float64,ExplicitSolutionElement})
  first(a) > first(b) ||
    (first(a) == first(b) && priority(last(a), last(b)))
end
# Return true if any solution element in A beats a solution in B
# before B beats A
function dominates(a::LocalAuctionAgent, b::LocalAuctionAgent)
  index = first_difference(a, b)

  if isnothing(index)
    # Default: solutions are identical or one has additional elements
    num_assignments(a) > num_assignments(b)
  else
    dominates(get_value_assignment(a, index), get_value_assignment(b, index))
  end
end

# Local maximization when updating another agents' solution.
# Chop off the incoming solution after winning because we cannot re-evaluate
# later values
#
# Proceed from the first difference. Maximize until winning or no other
# solution. Chop the incoming assignments and append new solution element
function local_agent_maximize(problem::ExplicitPartitionProblem,
                              agent::LocalAuctionAgent,
                              other_agent::LocalAuctionAgent)
  difference_index = first_difference(agent, other_agent)

  if isnothing(difference_index)
    # We already know that other_agent is not dominated
    difference_index = num_assignments(agent) + 1
  end

  for ii in difference_index:num_assignments(other_agent)+1
    prev_assignments = get_assignments(other_agent)[1:ii-1]
    prev_values = get_values(other_agent)[1:ii-1]

    x = solve_block(problem, get_index(agent), prev_assignments)

    value = marginal_gain(y->objective(problem, y), x, prev_assignments)

    if num_assignments(other_agent) < ii ||
      dominates((value, x), get_value_assignment(other_agent, ii))

      return LocalAuctionAgent(agent, vcat(prev_assignments, x),
                               vcat(prev_values, value))
    end
  end

  error("Failed to update solution")
end

# Update the agent using only local information
#
# * Compare solutions, solution with the first winning selection dominates
#   * Return the identity if the local agent wins
# * If agent already has a solution element in the new solution, return this
#   solution (both will have the same solution element for the local agent by
#   monotonicity of solutions)
# * Proceed from the first difference. Maximize until winning or no other
#   solution. Chop the incoming assignments and append new solution element
#
function local_agent_update(problem::ExplicitPartitionProblem,
                            agent::LocalAuctionAgent,
                            other_agent::LocalAuctionAgent)
  if dominates(agent, other_agent)
    agent
  elseif !can_augment(get_assignments(other_agent), get_index(agent))
    LocalAuctionAgent(agent, get_assignments(other_agent),
                      get_values(other_agent))
  else
    local_agent_maximize(problem, agent, other_agent)
  end
end

function solve_problem(solver::LocalAuctionSolver, problem::PartitionProblem;
                       kwargs...)
  num_agents = get_num_agents(problem)
  adjacency = solver.adjacency
  agents = map(x->LocalAuctionAgent(x), 1:num_agents)

  for ii in 1:solver.nominal_steps
    # Update span
    if ii != 1
      # Skip the first round (which we can complete without communication)
      solver.span += 1
    end

    updated_agents = map(agents) do agent
      ns = neighbors(adjacency, agents, agent)

      # Update messages and volume
      if ii != 1
        solver.messages += length(ns)
        solver.volume += mapreduce(num_assignments, +, ns, init=0)
      end

      # Try maximizing on all neighbors in sequence.
      #
      # Alternatively, I could identify which solution dominates
      dominating_agent = foldl(ns, init=agent) do best, next
        dominates(next, best) ? next : best
      end

      local_agent_update(problem, agent, dominating_agent)
    end

    agents[:] = updated_agents

    if converged(agents)
      break
    end
  end

  solver.solved = true

  assignments = map(get_assignment, agents)

  # Return the solution
  evaluate_solution(problem, assignments)
end

ext: jl

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
# Adopted from https://discourse.julialang.org/t/get-cpu-usage/24468/2
function get_cpu_stats()
  st = read("/proc/stat", String)
  lines = split(st, "\n")
  m = map(line -> match(r"^cpu(\d) (.*)", line), lines)
  matches = filter(m -> m !== nothing, m)
  Dict(parse(Int, match.captures[1]) => parse.(Int, split(match.captures[2], ' ')) for match in matches)
end

cpu_percentages = Dict(i => (0, 0) for i in 0:3)
function compute_percentage(fields, cpu_idx)
  lastidle, lasttotal = cpu_percentages[cpu_idx]
  idle, total = fields[4], sum(fields)
  Δidle, Δtotal = idle - lastidle, total - lasttotal
  utilization = 100 * (1 - Δidle / Δtotal)
  cpu_percentages[cpu_idx] = idle, total
  utilization
end


text:
module ExtendedVariantModeToggleDashDirectionTrigger

using ..Ahorn, Maple

@mapdef Trigger "ExtendedVariantMode/ToggleDashDirectionTrigger" ToggleDashDirectionTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16,
	dashDirection::Integer=0b1000000000, enable::Bool=true, revertOnLeave::Bool=false, revertOnDeath::Bool=true)

const placements = Ahorn.PlacementDict(
	"Toggle Dash Direction (Extended Variant Mode)" => Ahorn.EntityPlacement(
		ToggleDashDirectionTrigger,
		"rectangle"
	)
)

Ahorn.editingOptions(trigger::ToggleDashDirectionTrigger) = Dict{String, Any}(
	"dashDirection" => Dict{String, Any}(
		"Top" =>          0b1000000000,
		"Top-Right" =>    0b0100000000,
		"Right" =>        0b0010000000,
		"Bottom-Right" => 0b0001000000,
		"Bottom" =>       0b0000100000,
		"Bottom-Left" =>  0b0000010000,
		"Left" =>         0b0000001000,
		"Top-Left" =>     0b0000000100
	)
)

end


text:
import TimeZones: Transition
import TimeZones.Olson: ZoneDict, RuleDict, zoneparse, ruleparse, resolve, parsedate, order_rules
import Base.Dates: Hour, Minute, Second

### parsedate ###

# Variations of until dates
@test parsedate("1945") == (DateTime(1945), 'w')
@test parsedate("1945 Aug") == (DateTime(1945,8), 'w')
@test parsedate("1945 Aug 2") == (DateTime(1945,8,2), 'w')
@test parsedate("1945 Aug 2 12") == (DateTime(1945,8,2,12), 'w')  # Doesn't actually occur
@test parsedate("1945 Aug 2 12:34") == (DateTime(1945,8,2,12,34), 'w')
@test parsedate("1945 Aug 2 12:34:56") == (DateTime(1945,8,2,12,34,56), 'w')

# Make sure parsing can handle additional spaces.
@test parsedate("1945  Aug") == (DateTime(1945,8), 'w')
@test parsedate("1945  Aug  2") == (DateTime(1945,8,2), 'w')
@test parsedate("1945  Aug  2  12") == (DateTime(1945,8,2,12), 'w')
@test parsedate("1945  Aug  2  12:34") == (DateTime(1945,8,2,12,34), 'w')
@test parsedate("1945  Aug  2  12:34:56") == (DateTime(1945,8,2,12,34,56), 'w')

# Explicit zone "local wall time"
@test_throws ArgumentError parsedate("1945w")
@test_throws KeyError parsedate("1945 Augw")
@test_throws ArgumentError parsedate("1945 Aug 2w")
@test parsedate("1945 Aug 2 12w") == (DateTime(1945,8,2,12), 'w')
@test parsedate("1945 Aug 2 12:34w") == (DateTime(1945,8,2,12,34), 'w')
@test parsedate("1945 Aug 2 12:34:56w") == (DateTime(1945,8,2,12,34,56), 'w')

# Explicit zone "UTC time"
@test_throws ArgumentError parsedate("1945u")
@test_throws KeyError parsedate("1945 Augu")
@test_throws ArgumentError parsedate("1945 Aug 2u")
@test parsedate("1945 Aug 2 12u") == (DateTime(1945,8,2,12), 'u')
@test parsedate("1945 Aug 2 12:34u") == (DateTime(1945,8,2,12,34), 'u')
@test parsedate("1945 Aug 2 12:34:56u") == (DateTime(1945,8,2,12,34,56), 'u')

# Explicit zone "standard time"
@test_throws ArgumentError parsedate("1945s")
@test_throws KeyError parsedate("1945 Augs")
@test_throws ArgumentError parsedate("1945 Aug 2s")
@test parsedate("1945 Aug 2 12s") == (DateTime(1945,8,2,12), 's')
@test parsedate("1945 Aug 2 12:34s") == (DateTime(1945,8,2,12,34), 's')
@test parsedate("1945 Aug 2 12:34:56s") == (DateTime(1945,8,2,12,34,56), 's')

# Invalid zone
@test_throws ArgumentError parsedate("1945 Aug 2 12:34i")

# Actual until date found in Zone "Pacific/Apia"
@test parsedate("2011 Dec 29 24:00") == (DateTime(2011,12,30), 'w')

### order_rules ###

# Rule    Poland  1918    1919    -   Sep 16  2:00s   0       -
# Rule    Poland  1919    only    -   Apr 15  2:00s   1:00    S
# Rule    Poland  1944    only    -   Apr  3  2:00s   1:00    S
rule_a = ruleparse("1918", "1919", "-", "Sep", "16", "2:00s", "0", "-")
rule_b = ruleparse("1919", "only", "-", "Apr", "15", "2:00s", "1:00", "S")
rule_c = ruleparse("1944", "only", "-", "Apr", "3", "2:00s", "1:00", "S")

# Note: Alternatively we could be using `permutations` here.
for rules in ([rule_a, rule_b, rule_c], [rule_c, rule_b, rule_a], [rule_a, rule_c, rule_b])
    dates, ordered = order_rules(rules)

    @test dates == [Date(1918, 9, 16), Date(1919, 4, 15), Date(1919, 9, 16), Date(1944, 4, 3)]
    @test ordered == [rule_a, rule_b, rule_a, rule_c]
end

# ignore rules starting after the cutoff
dates, ordered = order_rules([rule_a, rule_b, rule_c], max_year=1940)
@test dates == [Date(1918, 9, 16), Date(1919, 4, 15), Date(1919, 9, 16)]
@test ordered == [rule_a, rule_b, rule_a]

# truncate rules ending after the cutoff
rule_pre = ruleparse("1999", "only", "-", "Jun", "7", "2:00s", "0", "P" )
rule_overlap = ruleparse("1999", "2001", "-", "Jan", "1", "0:00s", "0", "-")
rule_endless = ruleparse("1993", "max", "-", "Feb", "2", "6:00s", "0", "G")
rule_post = ruleparse("2002", "only", "-", "Jan", "1", "0:00s", "0", "IP")

dates, ordered = order_rules([rule_post, rule_endless, rule_overlap, rule_pre], max_year=2000)
@test dates == [
    Date(1993, 2, 2),
    Date(1994, 2, 2),
    Date(1995, 2, 2),
    Date(1996, 2, 2),
    Date(1997, 2, 2),
    Date(1998, 2, 2),
    Date(1999, 1, 1),
    Date(1999, 2, 2),
    Date(1999, 6, 7),
    Date(2000, 1, 1),
    Date(2000, 2, 2),
]
# Equality check based on reference
@test ordered == [
    rule_endless,
    rule_endless,
    rule_endless,
    rule_endless,
    rule_endless,
    rule_endless,
    rule_overlap,
    rule_endless,
    rule_pre,
    rule_overlap,
    rule_endless,
]

### resolve ###

warsaw = resolve("Europe/Warsaw", tzdata["europe"]...)

# Europe/Warsaw time zone has a combination of factors that requires computing
# the abbreviation to be done in a specific way.
@test warsaw.transitions[1].zone.name == :LMT
@test warsaw.transitions[2].zone.name == :WMT
@test warsaw.transitions[3].zone.name == :CET   # Standard time
@test warsaw.transitions[4].zone.name == :CEST  # Daylight saving time
@test issorted(warsaw.transitions)

zone = Dict{AbstractString,FixedTimeZone}()
zone["LMT"] = FixedTimeZone("LMT", 5040, 0)
zone["WMT"] = FixedTimeZone("WMT", 5040, 0)
zone["CET"] = FixedTimeZone("CET", 3600, 0)
zone["CEST"] = FixedTimeZone("CEST", 3600, 3600)
zone["EET"] = FixedTimeZone("EET", 7200, 0)
zone["EEST"] = FixedTimeZone("EEST", 7200, 3600)

@test warsaw.transitions[1] == Transition(typemin(DateTime), zone["LMT"])  # Ideally -Inf
@test warsaw.transitions[2] == Transition(DateTime(1879,12,31,22,36), zone["WMT"])
@test warsaw.transitions[3] == Transition(DateTime(1915,8,4,22,36), zone["CET"])
@test warsaw.transitions[4] == Transition(DateTime(1916,4,30,22,0), zone["CEST"])
@test warsaw.transitions[5] == Transition(DateTime(1916,9,30,23,0), zone["CET"])
@test warsaw.transitions[6] == Transition(DateTime(1917,4,16,1,0), zone["CEST"])
@test warsaw.transitions[7] == Transition(DateTime(1917,9,17,1,0), zone["CET"])
@test warsaw.transitions[8] == Transition(DateTime(1918,4,15,1,0), zone["CEST"])
@test warsaw.transitions[9] == Transition(DateTime(1918,9,16,1,0), zone["EET"]) #
@test warsaw.transitions[10] == Transition(DateTime(1919,4,15,0,0), zone["EEST"])
@test warsaw.transitions[11] == Transition(DateTime(1919,9,16,0,0), zone["EET"])
@test warsaw.transitions[12] == Transition(DateTime(1922,5,31,22,0), zone["CET"]) #
@test warsaw.transitions[13] == Transition(DateTime(1940,6,23,1,0), zone["CEST"])

@test warsaw.transitions[14] == Transition(DateTime(1942,11,2,1,0), zone["CET"])
@test warsaw.transitions[15] == Transition(DateTime(1943,3,29,1,0), zone["CEST"])
@test warsaw.transitions[16] == Transition(DateTime(1943,10,4,1,0), zone["CET"])
@test warsaw.transitions[17] == Transition(DateTime(1944,4,3,1,0), zone["CEST"])
@test warsaw.transitions[18] == Transition(DateTime(1944,10,4,0,0), zone["CET"])

# Zone Pacific/Honolulu contains the following properties which make it good for testing:
# - Zone's contain save in rules field
# - Zone abbreviation redefined: HST
# - Is not cutoff

honolulu = resolve("Pacific/Honolulu", tzdata["northamerica"]...)

zone = Dict{AbstractString,FixedTimeZone}()
zone["LMT"] = FixedTimeZone("LMT", -37886, 0)
zone["HST"] = FixedTimeZone("HST", -37800, 0)
zone["HDT"] = FixedTimeZone("HDT", -37800, 3600)
zone["HST_NEW"] = FixedTimeZone("HST", -36000, 0)

@test honolulu.transitions[1] == Transition(typemin(DateTime), zone["LMT"])
@test honolulu.transitions[2] == Transition(DateTime(1896,1,13,22,31,26), zone["HST"])
@test honolulu.transitions[3] == Transition(DateTime(1933,4,30,12,30), zone["HDT"])
@test honolulu.transitions[4] == Transition(DateTime(1933,5,21,21,30), zone["HST"])
@test honolulu.transitions[5] == Transition(DateTime(1942,2,9,12,30), zone["HDT"])
@test honolulu.transitions[6] == Transition(DateTime(1945,9,30,11,30), zone["HST"])
@test honolulu.transitions[7] == Transition(DateTime(1947,6,8,12,30), zone["HST_NEW"])

@test length(honolulu.transitions) == 7
@test isnull(honolulu.cutoff)

# Zone Pacific/Apia contains the following properties which make it good for testing:
# - Offset switch from -11:00 to 13:00
# - Rules interaction with a large negative offset
# - Rules interaction with a large positive offset
# - Includes a DateTime Julia could consider invalid: "2011 Dec 29 24:00"
# - Changed zone format while in a non-standard transition
# - Zone abbreviation redefined: LMT, WSST

apia = resolve("Pacific/Apia", tzdata["australasia"]...)

zone = Dict{AbstractString,FixedTimeZone}()
zone["LMT_OLD"] = FixedTimeZone("LMT", 45184, 0)
zone["LMT"] = FixedTimeZone("LMT", -41216, 0)
zone["WSST_OLD"] = FixedTimeZone("WSST", -41400, 0)
zone["SST"] = FixedTimeZone("SST", -39600, 0)
zone["SDT"] = FixedTimeZone("SDT", -39600, 3600)
zone["WSST"] = FixedTimeZone("WSST", 46800, 0)
zone["WSDT"] = FixedTimeZone("WSDT", 46800, 3600)

@test apia.transitions[1] == Transition(typemin(DateTime), zone["LMT_OLD"])
@test apia.transitions[2] == Transition(DateTime(1879,7,4,11,26,56), zone["LMT"])
@test apia.transitions[3] == Transition(DateTime(1911,1,1,11,26,56), zone["WSST_OLD"])
@test apia.transitions[4] == Transition(DateTime(1950,1,1,11,30), zone["SST"])
@test apia.transitions[5] == Transition(DateTime(2010,9,26,11), zone["SDT"])
@test apia.transitions[6] == Transition(DateTime(2011,4,2,14), zone["SST"])
@test apia.transitions[7] == Transition(DateTime(2011,9,24,14), zone["SDT"])
@test apia.transitions[8] == Transition(DateTime(2011,12,30,10), zone["WSDT"])
@test apia.transitions[9] == Transition(DateTime(2012,3,31,14), zone["WSST"])
@test apia.transitions[10] == Transition(DateTime(2012,9,29,14), zone["WSDT"])

# Zone Europe/Madrid contains the following properties which make it good for testing:
# - Observed midsummer time
# - End of midsummer time also switches both the UTC offset and the saving time
# - In 1979-01-01 switches from "Spain" to "EU" rules which could create a redundant entry
madrid = resolve("Europe/Madrid", tzdata["europe"]...)

zone = Dict{AbstractString,FixedTimeZone}()
zone["WET"] = FixedTimeZone("WET", 0, 0)
zone["WEST"] = FixedTimeZone("WEST", 0, 3600)
zone["WEMT"] = FixedTimeZone("WEMT", 0, 7200)
zone["CET"] = FixedTimeZone("CET", 3600, 0)
zone["CEST"] = FixedTimeZone("CEST", 3600, 3600)

@test madrid.transitions[23] == Transition(DateTime(1939,4,15,23), zone["WEST"])
@test madrid.transitions[24] == Transition(DateTime(1939,10,7,23), zone["WET"])
@test madrid.transitions[25] == Transition(DateTime(1940,3,16,23), zone["WEST"])
@test madrid.transitions[26] == Transition(DateTime(1942,5,2,22), zone["WEMT"])

@test madrid.transitions[33] == Transition(DateTime(1945,9,29,23), zone["WEST"])
@test madrid.transitions[34] == Transition(DateTime(1946,4,13,22), zone["WEMT"])
@test madrid.transitions[35] == Transition(DateTime(1946,9,29,22), zone["CET"])
@test madrid.transitions[36] == Transition(DateTime(1949,4,30,22), zone["CEST"])

# Redundant transition would be around 1979-01-01T00:00:00 as CET
@test madrid.transitions[47] == Transition(DateTime(1978,9,30,23), zone["CET"])
@test madrid.transitions[48] == Transition(DateTime(1979,4,1,1), zone["CEST"])

# Zone America/Anchorage contains the following properties which make it good for testing:
# - Uses a format containing a slash which indicates the abbreviations for STD/DST
#   Alternatives include: Europe/London, Europe/Dublin, and Europe/Moscow.
# - Observed war/peace time
anchorage = resolve("America/Anchorage", tzdata["northamerica"]...)

zone = Dict{AbstractString,FixedTimeZone}()
zone["CAT"] = FixedTimeZone("CAT", -36000, 0)
zone["CAWT"] = FixedTimeZone("CAWT", -36000, 3600)
zone["CAPT"] = FixedTimeZone("CAPT", -36000, 3600)
zone["AHST"] = FixedTimeZone("AHST", -36000, 0)
zone["AHDT"] = FixedTimeZone("AHDT", -36000, 3600)

@test anchorage.transitions[3] == Transition(DateTime(1900,8,20,21,59,36), zone["CAT"])
@test anchorage.transitions[4] == Transition(DateTime(1942,2,9,12), zone["CAWT"])
@test anchorage.transitions[5] == Transition(DateTime(1945,8,14,23), zone["CAPT"])
@test anchorage.transitions[6] == Transition(DateTime(1945,9,30,11), zone["CAT"])
@test anchorage.transitions[7] == Transition(DateTime(1967,4,1,10), zone["AHST"])
@test anchorage.transitions[8] == Transition(DateTime(1969,4,27,12), zone["AHDT"])

# Zone Europe/Ulyanovsk contains the following properties which make it good for testing:
# - With the exception of LMT all Zone and Rule abbreviations are UTC offsets which should
#   be treated as NULL.
ulyanovsk = resolve("Europe/Ulyanovsk", tzdata["europe"]...)
@test all(t -> string(t.zone.name) == "", ulyanovsk.transitions[2:end])

# Fake Zone Pacific/Cutoff contains the following properties which make it good for testing:
# - Having a single transition on the first year allows us to test the special case where we
#   need to include a cutoff while only having a single transition
# - Having no rules ensures that cutoff is calculated correctly with only zones

# Manually generate zones and rules as if we had read them from a file.
zones = ZoneDict()
rules = RuleDict()

zones["Pacific/Cutoff"] = [
    zoneparse("-10:00", "-", "CUT-1", "1933 Apr 1 2:00s"),
    zoneparse("-11:00", "-", "CUT-2", ""),
]

cutoff_tz = resolve("Pacific/Cutoff", zones, rules)

zone = Dict{AbstractString,FixedTimeZone}()
zone["CUT-1"] = FixedTimeZone("CUT-1", -36000)
zone["CUT-2"] = FixedTimeZone("CUT-2", -39600)

@test cutoff_tz.transitions[1] == Transition(typemin(DateTime), zone["CUT-1"])
@test cutoff_tz.transitions[2] == Transition(DateTime(1933,4,1,12), zone["CUT-2"])
@test isnull(cutoff_tz.cutoff)

cutoff_tz = resolve("Pacific/Cutoff", zones, rules, max_year=1900)

# Normally resolve TimeZones with a single transition are returned as a FixedTimeZone except
# when a cutoff is included.
@test isa(cutoff_tz, VariableTimeZone)
@test length(cutoff_tz.transitions) == 1
@test get(cutoff_tz.cutoff) == DateTime(1933,4,1,12)

# Behaviour of mixing "RULES" as a String and as a Time. In reality this behaviour has never
# been observed.

# Manually generate zones and rules as if we had read them from a file.
zones = ZoneDict()
rules = RuleDict()

zones["Pacific/Test"] = [
    zoneparse("-10:00", "-", "TST-1", "1933 Apr 1 2:00s"),
    zoneparse("-10:00", "1:00", "TDT-2", "1933 Sep 1 2:00s"),
    zoneparse("-10:00", "Testing", "T%sT-3", "1934 Sep 1 3:00s"),
    zoneparse("-10:00", "1:00", "TDT-4", "1935 Sep 1 3:00s"),
    zoneparse("-10:00", "Testing", "T%sT-5", ""),
]
rules["Testing"] = [
    ruleparse("1934", "1935", "-", "Apr", "1", "3:00s", "1", "D"),
    ruleparse("1934", "1935", "-", "Sep", "1", "3:00s", "0", "S"),
]

test = resolve("Pacific/Test", zones, rules)

zone = Dict{AbstractString,FixedTimeZone}()
zone["TST-1"] = FixedTimeZone("TST-1", -36000, 0)
zone["TDT-2"] = FixedTimeZone("TDT-2", -36000, 3600)
zone["TST-3"] = FixedTimeZone("TST-3", -36000, 0)
zone["TDT-3"] = FixedTimeZone("TDT-3", -36000, 3600)
zone["TDT-4"] = FixedTimeZone("TDT-4", -36000, 3600)
zone["TST-5"] = FixedTimeZone("TST-5", -36000, 0)
zone["TDT-5"] = FixedTimeZone("TDT-5", -36000, 3600)

@test test.transitions[1] == Transition(typemin(DateTime), zone["TST-1"])
@test test.transitions[2] == Transition(DateTime(1933,4,1,12), zone["TDT-2"]) # -09:00
@test test.transitions[3] == Transition(DateTime(1933,9,1,12), zone["TST-3"])
@test test.transitions[4] == Transition(DateTime(1934,4,1,13), zone["TDT-3"])
@test test.transitions[5] == Transition(DateTime(1934,9,1,13), zone["TDT-4"])
@test test.transitions[6] == Transition(DateTime(1935,9,1,13), zone["TST-5"])

# Note: Due to how the the zones/rules were written a redundant transition could be created
# such that `test.transitions[6] == test.transitions[7]`. The TimeZone code can safely
# handle redundant transitions but ideally they should be eliminated.
@test length(test.transitions) == 6

# Make sure that we can deal with Links. Take note that the current implementation converts
# links into zones which makes it hard to explicitly test for a link. We expect that the
# following link exists:
#
# Link  Europe/Oslo  Arctic/Longyearbyen

# Make sure that that the link time zone was parsed.
zone_names = keys(tzdata["europe"][1])
@test "Arctic/Longyearbyen" in zone_names

oslo = resolve("Europe/Oslo", tzdata["europe"]...)
longyearbyen = resolve("Arctic/Longyearbyen", tzdata["europe"]...)

@test oslo.transitions == longyearbyen.transitions

# Zones that don't include multiple lines and no rules should be treated as a FixedTimeZone.
mst = resolve("MST", tzdata["northamerica"]...)
@test isa(mst, FixedTimeZone)


text:
using EmpiricalRisks
using Base.Test
import EmpiricalRisks: gets

## auxiliary testing facilities

function _val_and_addgrad(rm::SupervisedRiskModel,
                          β::Float64, g0::StridedVecOrMat, α::Float64, θ::StridedVecOrMat, X, y)
    value_and_addgrad!(rm, β, copy(g0), α, θ, X, y)
end

function verify_risk(pm::PredictionModel, loss::Loss,
                     θ::VecOrMat{Float64}, X::Matrix{Float64}, y::VecOrMat)

    rm = riskmodel(pm, loss)

    # produce ground-truth
    n = size(X, 2)
    @assert ninputs(pm, X) == n
    @assert size(y, ndims(y)) == n

    R0 = zeros(n)
    if ndims(θ) == 1
        d = length(θ)
        G0 = zeros(d, n)
        for i = 1:n
            r, g = _risk_and_grad(pm, loss, θ, gets(X,i), gets(y,i))
            R0[i] = r
            G0[:,i] = g
        end
    else
        d1, d2 = size(θ)
        G0 = zeros(d1, d2, n)
        for i = 1:n
            r, g = _risk_and_grad(pm, loss, θ, gets(X,i), gets(y,i))
            R0[i] = r
            G0[:,:,i] = g
        end
    end

    # perform verification

    # over individual samples
    g0 = rand(size(θ))

    for i = 1:n
        x_i = gets(X, i)
        y_i = gets(y, i)
        g_i = gets(G0, i)
        v_i = R0[i]

        @test_approx_eq v_i value(rm, θ, x_i, y_i)

        v, g = value_and_grad(rm, θ, x_i, y_i)
        @test_approx_eq v_i v
        @test_approx_eq g_i g

        for β in [0.0, 0.5, 1.0], α in [0.5, 1.0, 2.0]
            v, g = _val_and_addgrad(rm, β, g0, α, θ, x_i, y_i)
            @test_approx_eq α * v_i v
            @test_approx_eq β * g0 + α * g_i g
        end
    end

    # over sample batch
    rv = sum(R0)
    rg = gets(sum(G0, ndims(G0)), 1)

    @test_approx_eq rv value(rm, θ, X,y)

    buffer = ndims(θ) == 2 ? zeros(size(θ,1), size(X,2)) : zeros(size(X,2))
    @test_approx_eq rv value!(buffer, rm, θ, X, y)

    v, g = value_and_grad(rm, θ, X, y)
    @test_approx_eq rv v
    @test_approx_eq rg g

    for β in [0.0, 0.5, 1.0], α in [0.5, 1.0, 2.0]
        v, g = _val_and_addgrad(rm, β, g0, α, θ, X, y)
        @test_approx_eq α * rv v
        @test_approx_eq β * g0 + α * rg g
    end
end

## (safe) reference implementation to compare with

function _risk_and_grad(::LinearPred, ::SqrLoss, w::StridedVector, x::StridedVector, y::Real)
    r = dot(w, x) - y
    (r^2 / 2, x * r)
end

function _risk_and_grad(pm::AffinePred, ::SqrLoss, wa::StridedVector, x::StridedVector, y::Real)
    x_ = [x; pm.bias]
    r = dot(wa, x_) - y
    (r^2 / 2, x_ * r)
end

function _risk_and_grad(::MvLinearPred, ::SumSqrLoss, W::StridedMatrix, x::StridedVector, y::StridedVector)
    r = W * x - y
    (sumabs2(r) / 2, r * x')
end

function _risk_and_grad(pm::MvAffinePred, ::SumSqrLoss, Wa::StridedMatrix, x::StridedVector, y::StridedVector)
    d = length(x)
    x_ = [x; pm.bias]
    p = Wa * x_
    r = p - y
    (sumabs2(r) / 2, r * x_')
end

function _risk_and_grad(pm::MvLinearPred, ::MultiLogisticLoss, W::StridedMatrix, x::StridedVector, y::Int)
    p = W * x
    ep = exp(p)
    q = ep ./ sum(ep)
    q[y] -= 1
    (log(sum(ep)) - p[y], q * x')
end

### Univariate prediction + Univariate Loss

d = 5
n = 8
a = randn()
bias = 2.5
w = randn(d)
wa = [w; a]
X = randn(d, n)
y = randn(n)

verify_risk(LinearPred(d), SqrLoss(), w, X, y)
verify_risk(AffinePred(d), SqrLoss(), wa, X, y)
verify_risk(AffinePred(d, bias), SqrLoss(), wa, X, y)

### Multivariate prediction + Univariate Loss

k = 3
W = randn(k, d)
a = randn(k)
Wa = [W a]
X = randn(d, n)
y = randn(k, n)

verify_risk(MvLinearPred(d, k), SumSqrLoss(), W, X, y)
verify_risk(MvAffinePred(d, k), SumSqrLoss(), Wa, X, y)
verify_risk(MvAffinePred(d, k, bias), SumSqrLoss(), Wa, X, y)

y = [1, 2, 3, 3, 2, 1, 2, 3]
@assert length(y) == n
verify_risk(MvLinearPred(d, k), MultiLogisticLoss(), W, X, y)


text:
using DataStructures

export fitness, lsfits, fitlevs, levcounts, neighbors, neutralnet, neutralnets, netcounts

@doc """lsfits(ls::Landscape)

Return a vector containing the fitnesses of all possible genotypes for the
given landscape.
"""
function lsfits(ls::Landscape)
  fits = zeros(Float64, ls.a^ls.n)
  for i = 0:(ls.a^ls.n - 1)
    fits[i+1] = fitness(Genotype(i, ls))
  end
  return fits
end

@doc """fitlevs(ls::Landscape, intervals::Int64, fits::Vector{Float64}, epsilon::Float64)

Return a vector of fitness levels of each genotype possible on the given
landscape. The fitness level is the index of the bin into which it would fall
assuming the interval [0,1] was divided into `intervals` bins.
"""
function fitlevs(ls::Landscape, intervals::Int64, fits::Vector{Float64}, epsilon::Float64)
  levels = zeros(Int64, ls.a^ls.n)
  for i = 0:(ls.a^ls.n - 1)
    levels[i + 1] = floor(Int64, (fits[i + 1] + epsilon) * intervals)
  end
  return levels
end

@doc """fitlevs(ls::NKqLandscape, intervals::Int64, fits::Vector{Float64})
"""
function fitlevs(ls::NKqLandscape, intervals::Int64, fits::Vector{Float64})
  fitlevs(ls, intervals, fits, eps())
end

@doc """fitlevs(ls::Landscape, intervals::Int64, fits::Vector{Float64})
"""
function fitlevs(ls::NKLandscape, intervals::Int64, fits::Vector{Float64})
  fitlevs(ls, intervals, fits, 0.0)
end

@doc """fitlevs(ls::Landscape, intervals::Int64)
"""
function fitlevs(ls::Landscape, intervals::Int64)
  fitlevs(ls, intervals, lsfits(ls))
end

@doc """levcounts(ls::Landscape, intervals::Int64, levels::Vector{Int64})

Return a vector such that the i'th element of the vector contains the number of
genotypes that fall into the i'th fitness level bin.
"""
function levcounts(ls::Landscape, intervals::Int64, levels::Vector{Int64})
  counts = zeros(Int64, intervals)
  for i = 0:(ls.a^ls.n - 1)
    counts[levels[i + 1] + 1] += 1
  end
  return counts
end

@doc """levcounts(ls::Landscape, intervals::Int64, fits::Vector{Float64})
"""
function levcounts(ls::Landscape, intervals::Int64, fits::Vector{Float64})
  levcounts(ls, intervals, fitlevs(ls, intervals, fits))
end

@doc """levcounts(ls::Landscape, intervals::Int64)
"""
function levcounts(ls::Landscape, intervals::Int64)
  levcounts(ls, intervals, lsfits(ls))
end

@doc """ function neutralnet(g::Genotype, levels::Vector{Int64})

Returns the alleles of the genotypes of the connected neutral network of which the given genotype
is a part as a vector of `AlleleString`s. Two genotypes are considered to be
neighbors if they have the same fitness level. Two genotypes are considered to
be connected if they are neighbors (they differ at one locus).

`levels` is a vector of fitness levels, or bins, as returned by the `fitlevs`
function.
"""
function neutralnet(g::Genotype, levels::Vector{Int64})
  f0 = levels[g.alleles + 1]
  closed = Set{AlleleString}()
  stack = Stack(Genotype)
  push!(stack, g)
  while !isempty(stack)
    ng = pop!(stack)
    if !in(ng.alleles, closed)
      push!(closed, ng.alleles)
      for nbr = filter(nbr -> levels[nbr.alleles + 1] == f0, all_neighbors(ng))
        push!(stack, nbr)
      end
    end
  end
  return collect(closed)
end

@doc """neutralnets(ls::Landscape, levels::Vector{Int64})

Returns an instance of `IntDisjointSet` that contains, as disjoint sets, the
connected neutral networks present in the given landscape.
"""
function neutralnets(ls::Landscape, levels::Vector{Int64})
  nets = IntDisjointSets(ls.a^ls.n)
  for i = 0:(ls.a^ls.n - 1)
    g = Genotype(i,ls)
    if (find_root(nets, i + 1) != i + 1) || (nets.ranks[i + 1] > 0)
      continue   # g has already been accounted for
    end
    net = neutralnet(g, levels)
    if length(net) > 0
      push!(net, i)  # add i to net
    else
      net = AlleleString[i]   # Array containing only i
    end
    for j = 2:length(net)
      union!(nets, net[j - 1] + 1, net[j] + 1)
    end
  end
  return nets
end

@doc """netcounts(nets::IntDisjointSets, levels::Vector{Int64})

Returns a vector of triples, each of which corresponds to a connected neutral
network within the given `IntDisjointSets` instance, which will most likely
come from a call to `neutralnets(...)`. Each triple consists of the following
data:

  1. Integer genotype representative of the network
  2. Number of genotypes in the network
  3. Integer fitness level of the network

The resulting vector is sorted by the number of genotypes in each network.

`levels` is the fitness levels vector that was provided to `neutralnets`.
"""
function netcounts(nets::IntDisjointSets, levels::Vector{Int64})
  c = counter(nets.parents)
  nn_list = map(x -> (x, c[x], levels[x]), keys(c))
  sort!(nn_list, by=x -> x[2]) # Sort by network size
  return nn_list
end


text:
using JuMP
using Cbc

LP=Model(Cbc.Optimizer)
@variable(LP, x>=5)
@variable(LP, y>=6)
@constraint(LP, x+y==18)
@constraint(LP, 20x+30y<=470)
@objective(LP, Max, 2x+3y)
status= JuMP.optimize!(LP)
println("X: $(JuMP.value(x))")
println("Y : $(JuMP.value(y))")
println("zzz")

max_stars_repo_name: valenmgama/optimization-tec

max_issues_repo_name: valenmgama/optimization-tec

max_forks_repo_name: valenmgama/optimization-tec


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Mortar2D.jl/blob/master/LICENSE

using PyPlot
using Mortar2D: calculate_normals, project_from_master_to_slave, project_from_slave_to_master, calculate_segments, calculate_mortar_matrices, calculate_mortar_assembly

coords = Dict(1 => [8.0, 10.0],
         2 => [7.0, 7.0],
         3 => [4.0, 3.0],
         4 => [0.0, 0.0],
         5 => [-3.0, 0.0],
         6 => [12.0, 10.0],
         7 => [10.0, 4.0],
         8 => [7.0, 2.0],
         9 => [4.0, -2.0],
         10 => [0.0, -3.0],
         11 => [-4.0, -3.0])

elements = Dict(
    1 => [1, 2],
    2 => [2, 3],
    3 => [3, 4],
    4 => [4, 5],
    5 => [6, 7],
    6 => [7, 8],
    7 => [8, 9],
    8 => [9, 10],
    9 => [10, 11])

slave_element_ids = [1, 2, 3, 4]
slave_elements = Dict(i => elements[i] for i in slave_element_ids)
master_element_ids = [5, 6, 7, 8, 9]
element_types = Dict(i => :Seg2 for i=1:length(elements))
normals = calculate_normals(slave_elements, element_types, coords)
segmentation = calculate_segments(slave_element_ids, master_element_ids,
                                  elements, element_types, coords, normals)

function plot1(; plot_element_normals=false, plot_nodal_normals=false, plot_segmentation=false)

    figure(figsize=(4, 4))

    for i in slave_element_ids
        x1,y1 = coords[elements[i][1]]
        x2,y2 = coords[elements[i][2]]
        xmid = 1/2*(x1+x2)
        ymid = 1/2*(y1+y2)
        n = [y1-y2, x2-x1]
        n /= norm(n)
        plot([x1,x2], [y1,y2], "-bo")
        if plot_element_normals
            arrow(x1, y1, n[1], n[2], head_width=0.1, head_length=0.2, fc="b", ec="b")
            arrow(x2, y2, n[1], n[2], head_width=0.1, head_length=0.2, fc="b", ec="b")
        end
    end

    for i in master_element_ids
        x1,y1 = coords[elements[i][1]]
        x2,y2 = coords[elements[i][2]]
        plot([x1,x2], [y1,y2], "-ro")
    end

    if plot_nodal_normals
        for i in keys(normals)
            x = coords[i]
            n = normals[i]
            arrow(x[1], x[2], n[1], n[2], head_width=0.1, head_length=0.2, fc="k", ec="k")
        end
    end

    if plot_segmentation
        for (sid, seg) in segmentation
            scon = elements[sid]
            xs1 = coords[scon[1]]
            xs2 = coords[scon[2]]
            ns1 = normals[scon[1]]
            ns2 = normals[scon[2]]
            for (mid, xi) in seg
                mcon = elements[mid]
                xm1 = coords[mcon[1]]
                xm2 = coords[mcon[2]]
                for xi_s in xi
                    #xi_s = (1-xi)/2*p + (1+xi)/2*p
                    N1 = [(1-xi_s)/2 (1+xi_s)/2]
                    n_s = N1[1]*ns1 + N1[2]*ns2
                    x_g = N1[1]*xs1 + N1[2]*xs2
                    xi_m = project_from_slave_to_master(Val{:Seg2}, x_g, n_s, xm1, xm2)
                    N2 = [(1-xi_m)/2 (1+xi_m)/2]
                    x_m = N2[1]*xm1 + N2[2]*xm2
                    plot([x_g[1], x_m[1]], [x_g[2], x_m[2]], "--k")
                end
            end
        end
    end

    axis("equal")
    annotate("\$\\Gamma_{\\mathrm{c}}^{\\left(1\\right)}\$", xy=(2.5, 5.0))
    annotate("\$\\Gamma_{\\mathrm{c}}^{\\left(2\\right)}\$", xy=(8.0, 0.0))
    axis("off")
end


text:
insts = OrderedDict()
insts["minimal"] = [
    ((2, 3, 4, false, 0, 0, 0, 0, 0),),
    ((5, 3, 4, false, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((5, 3, 4, false, 0.1, 0.1, 0.1, 0.2, 0.2), :SOCExpPSD),
    ((5, 4, 4, true, 0.1, 0, 0, 0, 0),),
    ((5, 3),),
    ((5, 3), :SOCExpPSD),
    ]
insts["fast"] = [
    ((3, 4, 5, false, 0, 0, 0, 0, 0),),
    ((3, 4, 5, false, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((3, 4, 5, false, 0, 0.1, 0.1, 0, 0),),
    ((20, 5),),
    ((6, 5), :SOCExpPSD),
    ((10, 20, 20, false, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((20, 10, 20, true, 0.1, 0, 0, 0, 0),),
    ((50, 8, 12, false, 0, 0, 0, 0, 0),),
    ((50, 8, 12, false, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((50, 8, 12, false, 0.1, 0.1, 0.1, 0.2, 0.2), :SOCExpPSD),
    ((15, 8, 12, true, 0.1, 0, 0, 0, 0),),
    ((15, 8, 8, true, 0, 0, 0.1, 0, 0), :SOCExpPSD),
    ]
insts["various"] = [
    ((3, 2, 80, true, 0.1, 0.1, 0.1, 0.2, 0.2), :SOCExpPSD),
    ((3, 2, 80, true, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((3, 2, 80, false, 0.1, 0.1, 0.1, 0.2, 0.2), :SOCExpPSD),
    ((3, 2, 80, false, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((6, 4, 80, true, 0.1, 0.1, 0.1, 0.2, 0.2), :SOCExpPSD),
    ((6, 4, 80, true, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((6, 4, 80, false, 0.1, 0.1, 0.1, 0.2, 0.2), :SOCExpPSD),
    ((6, 4, 80, false, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((12, 8, 80, true, 0.1, 0.1, 0.1, 0.2, 0.2), :SOCExpPSD),
    ((12, 8, 80, true, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ((12, 8, 80, false, 0.1, 0.1, 0.1, 0.2, 0.2),),
    ]
return (MatrixRegressionJuMP, insts)


text:
# https://adventofcode.com/2021/day/11
using AdventOfCode

input = split.(readlines("2021/data/day_11.txt"),"") |> Base.Fix1(reduce,hcat) .|> parse

function part_1(input)
    oct=copy(input)
    CIs=CartesianIndices(oct)
    fI, lI = first(CIs), last(CIs)
    I1 = oneunit(fI)
    flashes=0
    for t in 1:100
        oct.+=1
        while(any(oct.>9))
            excitedI=CIs[oct .> 9]
            for eI in excitedI
                neighbors=@view (oct[max(fI,eI-I1):min(lI,eI+I1)])
                neighbors.+=sign.(neighbors)
            end
            oct[excitedI].=0
        end
        flashes+=sum(==(0),oct)
    end
    flashes
end

function part_2(input)
    oct=copy(input)
    CIs=CartesianIndices(oct)
    fI, lI = first(CIs), last(CIs)
    I1 = oneunit(fI)
    t=0
    while true
        t+=1
        oct.+=1
        while(any(oct.>9))
            excitedI=CIs[oct .> 9]
            for eI in excitedI
                neighbors=@view (oct[max(fI,eI-I1):min(lI,eI+I1)])
                neighbors.+=sign.(neighbors)
            end
            oct[excitedI].=0
        end
    if all(==(0),oct) return(t) end
    end
end


text:
@testset "arrays can be diffed" begin
    a1 = [1, 2, 3, 4]

    # one number changed
    a2 = [1, 7, 3, 4]
    d = deepdiff(a1, a2)
    @test added(d) == [2]
    @test removed(d) == [2]
    @test changed(d) == []
    @test before(d) == a1
    @test after(d) == a2

    # removed from middle
    d = deepdiff(a1, [1, 3, 4])
    @test removed(d) == [2]
    @test added(d) == []

    # removed from beginning
    d = deepdiff(a1, [2, 3, 4])
    @test removed(d) == [1]
    @test added(d) == []

    # removed from end
    d = deepdiff(a1, [1, 2, 3])
    @test removed(d) == [4]
    @test added(d) == []

    # added to end
    d = deepdiff(a1, [1, 2, 3, 4, 5])
    @test removed(d) == []
    @test added(d) == [5]

    # added to beginning
    d = deepdiff(a1, [0, 1, 2, 3, 4])
    @test removed(d) == []
    @test added(d) == [1]

    # two additions
    d = deepdiff(a1, [1, 4, 2, 5, 3, 4])
    @test removed(d) == []
    @test added(d) == [2, 4]

    # two removals
    d = deepdiff(a1, [2, 4])
    @test removed(d) == [1, 3]
    @test added(d) == []
end


text:
@testset "$TEST $G" begin

g5 = DG(4)
add_edge!(g5,1,2); add_edge!(g5,2,3); add_edge!(g5,1,3); add_edge!(g5,3,4)
@test isapprox(pagerank(g5)[3], 0.318,  atol=0.001)
@test_throws ErrorException pagerank(g5, 2)
@test_throws ErrorException pagerank(g5, 0.85, 2)

end # testset


text:
# Loop Perforation for Approximate Computation
#
# ############################

# Loops over an expression block or function and finds the types of loops
# which we are looking to modify and applies the modification function - `replacementfn`
function LoopPerforation( ex, looptype, replacementfn )
    for i in 1:length(ex.args)
        if(typeof(ex.args[i]) == Expr && length(ex.args[i].args) > 1)
            LoopPerforation(ex.args[i], looptype, replacementfn)
        end

        if(typeof(ex.args[i]) == Expr && ex.args[i].head == :for)
            res = EvaluateForLoop(ex.args[i].args[1], looptype, replacementfn)
            ex.args[i].args[1].args[2] = res
        end
    end
end

# When a for loop is found this function tests to see if it is the correct type
# and passes the evaluated Expr to the replacement function so that the replacement
# will be returned. If nothing can be done, the input is returned.
function EvaluateForLoop(ex,  looptype, replacementfn)
    c = 0
    try
        c  = eval(ex)
    catch
        println( "Unable to evaluate" )
    end

    if( isa(c, looptype) )
        return replacementfn(c)
    end

     return ex
end

# Simple provided perforation functions:

function ClipFrontAndBack(c)
    b = c.start+1
    e = c.stop-1
    return :($(b):$(e))
end

function OnlyOddIterations(c)
    b = c.start
    e = c.stop
    return :($(b):2:$(e))
end

function OnlyEvenIterations(c)
    b = c.start+1
    e = c.stop
    return :($(b):2:$(e))
end

function OnlyFirstHalf(c)
    b = c.start
    e = c.stop
    return :($(b):$(e)/2)
end

function OnlySecondHalf(c)
    b = c.start
    e = c.stop
    return :($(e)/2:$(e))
end


text:
macro grab_output(ex)
    quote
        mktemp() do fname, fout
            ret = nothing
            open(fname, "w") do fout
                redirect_stdout(fout) do
                    ret = $(esc(ex))
                end
            end
            ret, read(fname, String)
        end
    end
end