text:
using NiLangCore

const GLOBAL_PHASE = Ref(0.0)
const GLOBAL_AMPLITUDE = Ref(1.0)

export qphase
qphase() = GLOBAL_PHASE[]
qphase(θ) = GLOBAL_PHASE[] = θ
qamplitude() = GLOBAL_AMPLITUDE[]
qamplitude(x) = GLOBAL_AMPLITUDE[] = x

"""
Rx = cos(θ/2) + im*sin(θ/2)*X
2×2 Array{Basic,2}:
    cos((1/2)*θ)  -I*sin((1/2)*θ)
 -I*sin((1/2)*θ)     cos((1/2)*θ)
"""
function Rx!(bit::Bool, θ)
    s = sin(θ/2)
    if rand() < s^2
        qphase(qphase()+(s > 0 ? π/2 : -π/2))
        !bit, θ
    else
        bit, θ
    end
end
iRx!(bit::Bool, θ) = Rx!(bit, -θ)
@dual Rx! iRx!

"""
Ry = cos(θ/2) + im*sin(θ/2)*Y
2×2 Array{Basic,2}:
 cos((1/2)*θ)  -sin((1/2)*θ)
 sin((1/2)*θ)   cos((1/2)*θ)
"""
function Ry!(bit::Bool, θ)
    s = sin(θ/2)
    if rand() < s^2
        qphase(qphase()+(bit ⊻ (s < 0) ? π : 0.0))
        !bit, θ
    else
        bit, θ
    end
end
iRy!(bit::Bool, θ) = Ry!(bit, -θ)
@dual Ry! iRy!

"""
Rz = cos(θ/2) + im*sin(θ/2)*Z
2×2 Array{Basic,2}:
 -I*sin((1/2)*θ) + cos((1/2)*θ)                              0
                              0  I*sin((1/2)*θ) + cos((1/2)*θ)
"""
function Rz!(bit::Bool, θ)
    s = sin(θ/2)
    if bit
        qphase(qphase()+θ/2)
    else
        qphase(qphase()-θ/2)
    end
    bit, θ
end
iRz!(bit::Bool, θ) = Rz!(bit, -θ)
@dual Rz! iRz!

using Test, Random
@testset "Rx, Ry, Rz" begin
    Random.seed!(2)
    for i=1:10
        b = true
        @instr Rx!(b, π)
        @test b == false
        @instr Rx!(b, 2π)
        @test b == false
    end
    qphase(0.0)
    b = true
    @instr Rx!(b, π)
    @test qphase() == π/2

    for i=1:10
        b = true
        @instr Ry!(b, π)
        @test b == false
        @instr Ry!(b, 2π)
        @test b == false
    end
    qphase(0.0)
    b = true
    @instr Ry!(b, π)
    @test qphase() ≈ π

    for i=1:10
        b = true
        @instr Rz!(b, π)
        @test b == true
        @instr Rz!(b, 2π)
        @test b == true
    end
    qphase(0.0)
    b = true
    @instr Rz!(b, π/4)
    @test qphase() ≈ π/8

    qphase(0.0)
    b = false
    @instr Rz!(b, π/4)
    @test qphase() ≈ -π/8
end

@testset "quantum reversibility" begin
    Random.seed!(2)
    cum = [0, 0.0im]
    θ = π/4
    for i=1:100
        # compute <0|Rx(-θ)Rx(θ)|0> and <1|Rx(-θ)Rx(θ)|0>
        b = false
        qphase(0.0)
        @instr Rx!(b, θ)
        @instr (~Rx!)(b, θ)
        cum[Int(b)+1] += exp(im*qphase())
    end
    normalize!(cum)
    @test isapprox(cum[1], 1; atol=1e-1)
    @test isapprox(cum[2], 0; atol=1e-1)
end


text:
"""
# Summary

	abstract type Program <: Any

# Subtypes

	CmdProgram
	JuliaProgram
"""
abstract type Program end

"""
	infer_outputs(p::Program, inputs::Dict{String})

Infer the default outputs from `p::Program` and `inputs::Dict{String}`.
"""
function infer_outputs(p::Program, inputs::Dict{String})
	p.infer_outputs(inputs)
end

function check_keywords(p::Program, inputs::Dict{String}, outputs::Dict{String})
	Set(p.inputs) == Set(keys(inputs)) || error("ProgramInputError: $(p.name): inputs provided and demanded not identical.")
	Set(p.outputs) == Set(keys(outputs)) || error("ProgramOutputError: $(p.name): outputs provided and demanded not identical.")
end
function check_outputs_keywords(p::Program, outputs::Dict{String})
	Set(p.outputs) == Set(keys(outputs)) || error("ProgramOutputError: $(p.name): outputs provided and demanded not identical.")
end

"""
	check_dependency(p::Program)

Check dependencies listed in `p.cmd_dependencies`.
"""
function check_dependency(p::Program)
	foreach(check_dependency, p.cmd_dependencies)
end

function generate_run_uuid(inputs::Dict{String}, outputs::Dict{String})
	all_cmd = ``
	foreach(x -> append!(all_cmd.exec, to_cmd(x).exec), values(inputs))
	foreach(x -> append!(all_cmd.exec, to_cmd(x).exec), values(outputs))
	uuid5(UUID4, string(all_cmd))
end

"""
	inputs_completion(p::Program, inputs::Dict{String})

complete missing keywords in `inputs` and `outputs`.
"""
function inputs_completion(p::Program, inputs::Dict{String})
	value_type = fieldtypes(eltype(inputs))[2]  # value type of inputs
	for (i, keyword) in enumerate(p.inputs)
		if haskey(inputs, keyword)
			# check types
			value = convert_data_type(inputs[keyword], p.input_types[i])
			if !(value isa value_type)
				inputs = convert(Dict{String,Any}, inputs)
			end  # it is ok to replace in for-loop
			inputs[keyword] = value
		else
			# not provided, check default
			default = p.default_inputs[i]
			if isnothing(default)
				# no default
				throw(ErrorException("ArgumentError: Program '$(p.name)' requires '$keyword' in inputs, but it is not provided."))
			else
				if !(default isa value_type)
					inputs = convert(Dict{String,Any}, inputs)
				end  # it is ok to replace in for-loop
				inputs[keyword] = default
			end
		end
	end
	inputs
end
function outputs_completion(p::Program, outputs::Dict{String})
	value_type = fieldtypes(eltype(outputs))[2]  # value type of outputs
	for (i, keyword) in enumerate(p.outputs)
		if haskey(outputs, keyword)
			# check types
			value = convert_data_type(outputs[keyword], p.output_types[i])
			if !(value isa value_type)
				outputs = convert(Dict{String,Any}, outputs)
			end  # it is ok to replace in for-loop
			outputs[keyword] = value
		else
			# not provided, check default
			default = p.default_outputs[i]
			if isnothing(default)
				# no default
				throw(ErrorException("ArgumentError: Program '$(p.name)' requires $keyword in outputs, but it is not provided."))
			else
				if !(default isa value_type)
					outputs = convert(Dict{String,Any}, outputs)
				end  # it is ok to replace in for-loop
				outputs[keyword] = default
			end
		end
	end
	outputs
end

"""
	keyword_interpolation(inputs::Dict{String}, outputs::Dict{String})

Interpolate <keyword> in `String`.
"""
function keyword_interpolation(inputs::Dict{String}, outputs::Dict{String})
	allputs = Dict(union(inputs, outputs))
	keyword_interpolation(allputs::Dict{String})
	for key in keys(inputs)
		inputs[key] = allputs[key]
	end
	for key in keys(outputs)
		outputs[key] = allputs[key]
	end
	inputs, outputs
end
function keyword_interpolation(allputs::Dict{String})
	for key in keys(allputs)
		keyword_interpolation(allputs, key, 1)
	end
	allputs
end
function keyword_interpolation(allputs::Dict{String}, key::String, n_recursion::Int)
	value = allputs[key]
	keywords = find_keywords(value)
	isempty(keywords) && return allputs # no need to interpolate
	for keyword in keywords
		keyword_value = get(allputs, keyword, nothing)
		isnothing(keyword_value) && continue # keyword not match, ignore
		if isempty(find_keywords(keyword_value))
			allputs[key] = replace(allputs[key], "<$keyword>" => str(keyword_value))
		else
			n_recursion > 20 && throw(ErrorException("ProgramArgumentError: too many recursion occurs in keyword interpolation."))
			keyword_interpolation(allputs::Dict{String}, keyword::String, n_recursion + 1)
		end
	end
	allputs
end

function find_keywords(value::T) where T <: AbstractString
	m = eachmatch(r"<([^<>]+)>", value)
	String[i.captures[1] for i in m]
end
find_keywords(not_string) = []

"""
	xxputs_completion_and_check(p::Program, inputs::Dict{String}, outputs::Dict{String})

1. Check and complete `inputs` using types and values stored in `p`.

2. If `outputs` is empty, run `p.infer_outputs`.

3. Check and complete `outputs` using types and values stored in `p`.

4. Check keyword consistency using `p`.

5. Interpolate <keyword> in String in completed `inputs` and `outputs`.

6. Return inputs and outputs.
"""
function xxputs_completion_and_check(p::Program, inputs::Dict{String}, outputs::Dict{String})
	inputs = inputs_completion(p::Program, inputs::Dict{String})

	if isempty(p.outputs)
		# do nothing to outputs
	else
		if isempty(outputs) && p.infer_outputs !== do_nothing
			outputs = to_xxput_dict(p.infer_outputs(inputs))
		end
		outputs = outputs_completion(p::Program, outputs::Dict{String})
	end

	# check keyword consistency (keys in inputs and outputs compatible with the function)
	check_keywords(p, inputs, outputs)

	# parse <keyword> in String
	inputs, outputs = keyword_interpolation(inputs, outputs)
end

function xxputs_completion_and_check(p::Program, inputs, outputs)
	xxputs_completion_and_check(p, to_xxput_dict(inputs), to_xxput_dict(outputs))
end

function Base.run(p::Program, inputs, outputs; kwarg...)
	run(p; inputs=inputs, outputs=outputs, kwarg...)
end

function Base.run(p::Program, inputs; kwarg...)
	run(p; inputs=inputs, kwarg...)
end

function Base.run(p::Program;
	dir::AbstractString = "",
	stdout = nothing, stderr = nothing, stdlog = stderr, append::Bool = false,
	kwarg...
)
	dir_backup = pwd()
	dir == "" || cd(dir) # go to working directory

	res = redirect_to_files(stdout, stderr, stdlog; mode = append ? "a+" : "w+") do
		_run(p; kwarg...)
	end

	dir == "" || cd(dir_backup)
	res
end

ext: jl


text:
# © 2016 Massachusetts Institute of Technology.  See LICENSE file for details.
abstract StoragePool{T}

atype{T}(::StoragePool{T}) = T
Base.eltype{T}(::StoragePool{T}) = eltype(T)

immutable NullPool{T} <: StoragePool{T}
  token::Vector{Symbol}
  arrays::Vector{T}
end

NullPool(T::Type) = NullPool{T}([:nothing], T[])

allocate{T}(n::NullPool{T}, siz::Int...) = let a = T(siz...); push!(n.arrays, a); a end
deallocate{T}(::NullPool{T}, ::T) = nothing

Base.close{T<:Array}(p::NullPool{T}) = empty!(p.arrays)
@cuda function Base.close{T<:CudaArray}(p::NullPool{T})
  for a in p.arrays
    free(a)
  end
  empty!(p.arrays)
end

type PoolCounter
  allocated::Int
  deallocated::Int
  wasted::Int
  unwasted::Int

  maxAllocated::Int
  maxWasted::Int

  PoolCounter() = new(0,0,0,0,0,0)
end

function count_allocated(c::PoolCounter, a::Int, w::Int)
  c.allocated += a
  c.wasted += w
  if c.allocated - c.deallocated > c.maxAllocated
    c.maxAllocated = c.allocated - c.deallocated
    c.maxWasted = c.wasted - c.unwasted
  end
end

function count_deallocated(c::PoolCounter, a::Int, w::Int)
  c.deallocated += a
  c.unwasted += w
end

immutable BuddyPool{T} <: StoragePool{T}
  pool::T
  free::IntSet
  psize::Int
  levels::Int
  map::Dict{Int,Tuple{Int,Int,Int}}
  counter::PoolCounter

  token::Vector{Symbol}

  @compat (::Type{BuddyPool})(pool::Array, N=16) = let NN = min(N, maxlevels(length(pool)))
    init_freelist(new{Array{eltype(pool)}}(pool, IntSet(), div(nextpow2(length(pool)), 1<<(NN-1)), NN, Dict{Int,Int}(), PoolCounter(), [:nothing])) end

  @cuda @compat (::Type{BuddyPool})(pool::CudaArray, N=16) = let NN = min(N, maxlevels(length(pool)))
    init_freelist(new{CudaArray{eltype(pool)}}(pool, IntSet(), div(nextpow2(length(pool)), 1<<(NN-1)), NN, Dict{Int,Int}(), PoolCounter(), [:nothing])) end
end

Base.close{T<:Array}(p::BuddyPool{T}) = nothing
@cuda Base.close{T<:CudaArray}(p::BuddyPool{T}) = free(p.pool)

maxlevels(len) = 64-leading_zeros(UInt64(nextpow2(len)))-4

function init_freelist(b::BuddyPool)
  sizehint!(b.free, 1<<b.levels)
  for i in 1:div(length(b.pool), b.psize)
     deallocate_block(b, (1<<(b.levels-1)) + i - 1)
  end
  b
end

hsize(b::BuddyPool, i::Int) = hsize(sizeof(eltype(b))*i)
hsize(i::Int) = i < 2^10 ? string(i, " B") :
  i < 2^20 ? @sprintf("%4.2f KiB", i/2^10) :
  i < 2^30 ? @sprintf("%4.2f MiB", i/2^20) :
             @sprintf("%4.2f GiB", i/2^30)

counter(b::BuddyPool, i::Int) = "$i ($(hsize(b, i)))"

function print_stats(io::IO, b::BuddyPool)
  len = length(b.pool)
  usable = b.psize*div(len, b.psize)
  println(io, "BuddyPool{$(atype(b))}, $(len) elements ($(hsize(b, len))), $(usable) usable ($(hsize(b, usable)))")
  println(io, "Allocated: $(counter(b, b.counter.allocated)) total, $(counter(b, b.counter.maxAllocated)) max, $(counter(b, b.counter.allocated-b.counter.deallocated)) current")
  println(io, "Wasted: $(counter(b, b.counter.wasted)) total, $(counter(b, b.counter.maxWasted)) max, $(counter(b, b.counter.wasted-b.counter.unwasted)) current")
end

# order 0 is block 1
# order 1 is blocks 2-3
# order 2 is blocks 4-7 (2^order, 2^(order+1)]

function allocate_block(b::BuddyPool, order::Int)
  i,_ = next(b.free, 1<<order)
  if i < (1<<(order+1))
    delete!(b.free, i)
    return i
  end

  if order == 0
    error("memory pool is exhausted")
  end

  a = allocate_block(b, order-1)<<1
  push!(b.free, a+1)
  return a
end

function deallocate_block(b::BuddyPool, i::Int)
  if i > 1 && ((i$1) in b.free)
    delete!(b.free, i$1)
    deallocate_block(b, i>>1)
  else
    push!(b.free, i)
  end
end

function print_free_list(b::BuddyPool, start::Int, depth::Int)
  for i in 1:depth
    for j in 1:(2^(i-1))
      k = (start<<(i-1)) + j - 1
      local ch
      if k in b.free
        ch = "F"
      else
        ch = "x"
        kk = k
        while kk > 0
          if kk in b.free
            ch = " "
            break
          end
          kk >>= 1
        end
      end
      print(ch^(2^(depth-i)))
    end
    println()
  end
end

function allocate{T}(b::BuddyPool{T}, siz::Int...)
  len = prod(siz)
  len = div(len+b.psize-1, b.psize)
  len = nextpow2(len)
  order = leading_zeros(UInt64(len)) + b.levels - 64
  if order < 0
    error("allocation $siz is larger than the pool")
  end
  bl = allocate_block(b, order)
  count_allocated(b.counter, prod(siz), b.psize*len-prod(siz))
  pa = (((1<<order)-1)&bl) * (1<<(b.levels-order-1))
  b.map[pa] = (bl, prod(siz), len)
  block_to_array(b, pa, siz)
end

function deallocate{T}(b::BuddyPool{T}, p::T)
  pa = div(pointer(p) - pointer(b.pool), b.psize*sizeof(eltype(b.pool)))
  bl,siz,len = b.map[pa]
  delete!(b.map, pa)
  deallocate_block(b, bl)
  count_deallocated(b.counter, siz, b.psize*len-siz)
  nothing
end

block_to_array{T}(b::BuddyPool{Array{T}}, p, siz) = pointer_to_array(pointer(b.pool) + b.psize*p*sizeof(eltype(b.pool)), siz)
@cuda block_to_array{T}(b::BuddyPool{CudaArray{T}}, p, siz) = CudaArray(pointer(b.pool) + b.psize*p*sizeof(eltype(b.pool)), siz, device(b.pool))


text:
# Author: AJ Tyler
# Date: 24/01/22
module Classify # This module contains functions which classify configurations based on their symmetries.
using Graphs # This module allows graphical representation of configurations.
export findSymmetries # This function identifies symmetries in a configuration.

"""
This function identifies symmetries in a configuration. First the core atoms of the configuration are tried to be identified. If a core
is found, the caps are attempted to be identified. If no core was found, then symmetric components throughout the configurations are
investigated.

Inputs:
shells: (Vector{Vector{Int64}}) A vector of shell numbers of each atom in each configuration.
atomicProfiles: Array{Dict{String,Int}} An array containing the CNA profile of each configuration.
bondGraphs: Vector{SimpleGraph} A vector of graphical representation of each configuration. Empty if comparing.
L: (Int64) The number of configurations in the file.

Outputs:
classifications: Vector{Dict{Vector{Int64},String}} Vector of dictionaries mapping a list of atoms to their symmetry.
capSymmetries: Vector{Dict{Dict{Int64,Vector{Int64}},String}} Vector of dictionaries mapping list of atoms forming a cap to their symmetry.
"""
function findSymmetries(shells,atomicProfiles,bondGraphs,L)
	classifications = Vector{Dict{Vector{Int64},String}}(undef,L) # Initialise vector of classifications
	capSymmetries = Vector{Dict{Dict{Int64,Vector{Int64}},String}}(undef,L) # Initialise vector of capSymmetries
	for i in 1:L # For each of the configurations in the file
		hasCore, classifications[i] = classifyCore(shells[i],atomicProfiles[1,i,:],bondGraphs[i]) # Try classify cluster core
		if (hasCore) # Check if cluster has identifiable core
			capSymmetries[i] = classifyCaps(shells[i],atomicProfiles[1,i,:],bondGraphs[i]) # Identify caps
		else
			classifications[i] = classifySymmetricComponents(atomicProfiles[1,i,:],bondGraphs[i]) # Classify symmetric components
			capSymmetries[i] = Dict{Dict{Int64,Vector{Int64}},String}() # No caps identifiable
		end
	end

	return classifications,capSymmetries # Return symmetries
end

"""
This function determines the symmetry contained in a CNA profile.

Inputs:
profile: (Vector{Dict{String,Int}}) 1D vector containing atomic CNA profiles.

Ouputs:
classification (String) Detected symmetry
"""
function classifySymmetry(profile)
	n421,n422,n444,n666,n555 = 0,0,0,0,0 # Initialise counts of CNA profiles
	for bond in keys(profile) # For each CNA profile of the current atom
		if bond == "(4,2,1)" # If has the (4,2,1) profile
			n421 += get!(profile,bond,0) # Add to the counter
		elseif bond == "(4,2,2)" # etc
			n422 += get!(profile,bond,0)
		elseif bond == "(4,4,4)"
			n444 += get!(profile,bond,0)
		elseif bond == "(6,6,6)"
			n666 += get!(profile,bond,0)
		elseif bond == "(5,5,5)"
			n555 += get!(profile,bond,0)
		end
	end

	numBonds = sum([n421,n422,n444,n666,n555]) # Compute how many bonds the atom has
	if numBonds == 12 # If has coordination number of 12
		if n421 == 12 # If has 12 (4,2,1) bonds
			return "FCC" # The atom's bonds have FCC symmetry
		elseif (n421 == 6 && n422 == 6) # If has 6 (4,2,1) and 6 (4,2,2) bonds
			return "HCP" # The atom's bonds have HCP symmetry
		elseif n555 == 12 # If has 12 (5,5,5) bonds
			return "ICO" # The atom's bonds have icosahedral symmetry
		else
			return "OTHER"
		end
	elseif (numBonds == 14 && n444 == 6 && n666 ==8) # If has coordination number of 14 and 6 (4,4,4) and 8 (6,6,6) bonds
		return "BCC" # The atom's bonds have BCC symmetry
	else # If doesn't have any of the above symmetries
		return "OTHER" # Then classify as OTHER
	end
end

"""
This function first determines which shell numbers contain the core and then checks if all those atoms have the same CNA profile. If they
do, then a core has been identified. Otherwise, non-cental symmetric components are attempted to be identified instead.

Inputs:
shells: (Vector{Int64}) Shell numbers of each atom in the configuration.
atomicProfiles: Dict{String,Int} A dictionary containing the CNA profile of the configuration.
bondGraph: SimpleGraph A graphical representation of the configuration.

Outputs:
hasCore: (Bool) Whether or not a core has been identified.
classification: (Dict) Dictionary mapping a list of atoms to their symmetry.
"""
function classifyCore(shells,atomicProfiles,bondGraph)
	coreAtoms = findall(x -> x<=floor(0.2*(maximum(shells))),shells) # Determine which atoms may form a core
	coreProfile = atomicProfiles[coreAtoms] # Obtain CNA profiles of the core atoms
	for profile in coreProfile # For each of the core atom CNA profiles
		if (profile != coreProfile[1]) # If the CNA profiles arn't all the same
			return false, Dict{Vector{Int64},String}() # Then core cannot be identified.
		end
	end
	classification = Dict(coreAtoms => classifySymmetry(coreProfile[1])*"(Core)") # Construct core classification

	return true,classification # Found core and return classification
end

"""
This function first finds the atoms in the outer and second most outer shells. If there are more atoms in the second outer shell, then the
outermost shell atoms are cap atoms and all their neighbours are capped atoms.

Inputs:
shells: (Vector{Int64}) Shell numbers of each atom in the configuration.
atomicProfiles: Dict{String,Int} A dictionary containing the CNA profile of the configuration.
bondGraph: SimpleGraph A graphical representation of the configuration.

capSymmetries: (Dict) Dictionary mapping list of atoms forming a cap to their symmetry. The keys are themselves dictionaries, mapping
the cap atom to a vector of the atoms it caps.
"""
function classifyCaps(shells,atomicProfiles,bondGraph)
	capSymmetries = Dict{Dict{Int64,Vector{Int64}},String}() # Initialise dictionary of cap symmetries
	maxShell = maximum(shells) # Find outer shell number
	outerShell = findall(x -> x == maxShell,shells) # Find outer shell atoms
	secondOuterShell = findall(x -> x == maxShell-1,shells) # Find 2nd most outer shell atoms
	if length(outerShell) < length(secondOuterShell) # If there are more atoms in the second outer shell than the outer shell
		for capAtom in outerShell # Each of the outer shell atoms are cap atoms
			cappedAtoms = neighbors(bondGraph,capAtom) # The cap atom caps all of its neighbours
			for cappedAtom in cappedAtoms # For each of the capped atoms
				symmetry = classifySymmetry(atomicProfiles[cappedAtom]) # Identify the capped atom's symmetry
				if symmetry != "OTHER" # If the symmetry is not OTHER
					push!(capSymmetries,Dict(capAtom=>cappedAtoms)=>symmetry) # Then the cap has a known symmetry
					break
				end
			end
			push!(capSymmetries,Dict(capAtom=>cappedAtoms)=>"OTHER") # The cap has an OTHER symmetry
		end
	end

	return capSymmetries # Return cap symmetries
end

"""
This function classifies components of the bondgraph which have the same symmetry via a breadth-first-search.

Inputs:
atomicProfiles: Dict{String,Int} A dictionary containing the CNA profile of the configuration.
bondGraph: SimpleGraph A graphical representation of the configuration.

Outputs:
Classification: (Dict) Dictionary mapping a list of atoms to their symmetry.
"""
function classifySymmetricComponents(atomicProfiles,bondGraph)
	atomsToProcess = Set([1]) # Start search from first atom
	atomsProcessed = Set{Int64}() # Initialise set of processed atoms
	classification = Dict{Vector{Int64},String}() # Initialise found symmetries
	while !isempty(atomsToProcess) # While have atoms left to process
		atom = pop!(atomsToProcess) # Pop off atom to process
		push!(atomsProcessed,atom) # Add atom to list of processed atoms
		symmetry = classifySymmetry(atomicProfiles[atom]) # Classify symmetry of atom
		if (symmetry != "OTHER") # If symmetry is not OTHER
			key = Vector{Int64}([atom]) # Initialise list of atoms that have the same symmetry
			# Obtain list of neighbours of the atom that have not already been searched
			neighboursToSearch = setdiff(Set(neighbors(bondGraph,atom)),atomsProcessed)
			while (!isempty(neighboursToSearch)) # While have neighbours left to search
				neighbour = pop!(neighboursToSearch) # Pop off atom to process
				neighbourSymmetry = classifySymmetry(atomicProfiles[neighbour]) # Classify the neighbour's symmetry
				if (neighbourSymmetry == symmetry) # If neighbour has the same symmetry
					push!(key,neighbour) # Add neighbour to key
					# Add the neighbour's neighbours (that haven't already been processed) to search
					neighboursToSearch = union(neighboursToSearch,setdiff(Set(neighbors(bondGraph,neighbour)),atomsProcessed))
					push!(atomsProcessed,neighbour) # Add neighbour to processed atoms
					delete!(atomsToProcess,neighbour) # Remove neighbour from atoms to process
				else # If neighbour has different symmetry
					push!(atomsToProcess,neighbour) # Add neighbour to atomsToProcess
				end
			end

			push!(classification, key => symmetry) # Add component to classifications
		else # If atom's symmetry is OTHER
			# Add the atom's neighbours (that haven't already been processed) to search
			atomsToProcess= union(atomsToProcess,setdiff(Set(neighbors(bondGraph,atom)),atomsProcessed))
		end
	end

	return classification # Return classifications
end

end


text:
export Edge

"""
An `Edge` joins two interfaces (half-edges) `a` and `b`.
"""
mutable struct Edge <: AbstractEdge
    variable::AbstractVariable
    a::Union{Interface, Nothing}
    b::Union{Interface, Nothing}

    function Edge(var::AbstractVariable, a::Interface)
        current_graph = currentGraph()
        (a.partner == nothing && a.edge == nothing) || error("Interface is already connected to an Edge.")
        hasNode(current_graph, a.node) || error("Node $(a.node.id) does not belong to the current graph.")
        hasVariable(current_graph, var) || error("Variable $(var.id) does not belong to the current graph.")

        self = new(var, a, nothing)
        a.edge = self

        push!(var.edges, self)
        push!(current_graph.edges, self)

        return self
    end
end

Edge(var::AbstractVariable, a::Interface, b::Interface) = connect!(Edge(var, a), b)

function degree(edge::Edge)
    deg = 2
    if (edge.a == nothing) || isDeltaConstrained(edge.a)
        deg -= 1
    end
    if (edge.b == nothing) || isDeltaConstrained(edge.b)
        deg -= 1
    end

    return deg
end

"""
Connect loose end of edge to interface b.
"""
function connect!(edge::Edge, b::Interface)
    (edge.b == nothing) || error("The edge does not have a loose end.")
    (b.partner == nothing && b.edge == nothing) || error("Interface is already connected to an Edge.")
    hasNode(currentGraph(), b.node) || error("Node $(b.node.id) does not belong to the current graph.")

    edge.b = b
    b.edge = edge
    edge.a.partner = b
    b.partner = edge.a

    return edge
end

"""
Disconnect edge from interface.
"""
function disconnect!(edge::Edge, interface::Interface)
    if (edge.a !== interface) && (edge.b !== interface)
        error("Cannot disconnect from an interface that is not connected.")
    end

    edge.a.partner = nothing
    edge.b.partner = nothing
    if edge.a === interface
        edge.a.edge = nothing
        edge.a = edge.b
    else
        edge.b.edge = nothing
    end
    edge.b = nothing

    return edge
end

function show(io::IO, edge::Edge)
    if edge.a != nothing
        a = "$(edge.a.node.id)."
        if (a_handle = handle(edge.a)) != ""
            a *= ((typeof(a_handle)==Symbol) ? "i[:$(a_handle)]" : "i[$(a_handle)]")
        else
            a *= "interfaces[$(findfirst(isequal(edge.a), edge.a.node.interface))]"
        end
    else
        a = "NONE"
    end

    if edge.b != nothing
        b = "$(edge.b.node.id)."
        if (b_handle = handle(edge.b)) != ""
            b *= ((typeof(b_handle)==Symbol) ? "i[:$(b_handle)]" : "i[$(b_handle)]")
        else
            b *= "interfaces[$(findfirst(isequal(edge.b), edge.b.node.interfaces))]"
        end
    else
        b = "NONE"
    end

    println(io, "Edge belonging to variable $(edge.variable.id): ( $(a) )----( $(b) ).")
end

Base.isless(e1::Edge, e2::Edge) = isless(name(e1), name(e2))

name(edge::Edge) = name(edge.a)*name(edge.b)

function show(io::IO, edges::Union{Vector{Edge}, Set{Edge}})
    println(io, "Edges:")
    for edge in edges
        show(io, edge)
    end
end


text:
# ## Preliminary steps
#
# Let's generate a `DataFrame` with some dummy regression data, let's also load the good old ridge regressor.

using MLJ, DataFrames, Random
MLJ.color_off() # hide
@load RidgeRegressor pkg=MultivariateStats

Random.seed!(5) # for reproducibility
x1 = rand(300)
x2 = rand(300)
x3 = rand(300)
y = exp.(x1 - x2 -2x3 + 0.1*rand(300))
X = DataFrame(x1=x1, x2=x2, x3=x3)
first(X, 3) |> pretty

# Let's also prepare the train and test split which will be useful later on.

test, train = partition(eachindex(y), 0.8);

# ## Defining a learning network
#
# In MLJ, a *learning network* is a directed acyclic graph (DAG) whose *nodes* apply trained or untrained operations such as a `predict` or `transform` (trained) or `+`, `vcat` etc. (untrained).
# Learning networks can be seen as pipelines on steroids.
#
# Let's consider the following simple DAG:
#
# ![Operation DAG](/assets/diagrams/composite1.svg)
#
# It corresponds to a fairly standard regression workflow: the data is standardized, the target is transformed using a Box-Cox transformation, a ridge regression is applied and the result is converted back by inverting the transform.
#
# **Note**: actually  this DAG is simple enough that it could also have been done with a pipeline.
#
# ### Sources and nodes
#
# In MLJ a learning network starts at **source** nodes and flows through nodes (`X` and `y`) defining operations/transformations (`W`, `z`, `ẑ`, `ŷ`).
# To define the source nodes, use the `source` function, you should specify whether it's a target:

Xs = source(X)
ys = source(y, kind=:target)

# To define an "trained-operation" node, you must simply create a machine wrapping a model and another node (the data) and indicate which operation should be performed (e.g. `transform`):

stand = machine(Standardizer(), Xs)
W = transform(stand, Xs)

# You can `fit!` a trained-operation node at any point, MLJ will fit whatever it needs that is upstream of that node.
# In this case, there is just a source node upstream of `W` so fitting `W` will just fit the standardizer:

fit!(W, rows=train);

# If you want to get the transformed data, you can then call the node speciying on which part of the data the operation should be performed:

W()             # transforms all data
W(rows=test, )  # transforms only test data
W(X[3:4, :])    # transforms specific data

# Let's now define the other nodes:

box_model = UnivariateBoxCoxTransformer()
box = machine(box_model, ys)
z = transform(box, ys)

ridge_model = RidgeRegressor(lambda=0.1)
ridge = machine(ridge_model, W, z)
ẑ = predict(ridge, W)

ŷ = inverse_transform(box, ẑ)

# Note that we have not yet done any training, but if we now call `fit!` on `ŷ`, it will fit all nodes upstream of `ŷ` that need to be re-trained:

fit!(ŷ, rows=train);

# Now that `ŷ` has been fitted, you can apply the full graph on test data (or any compatible data). For instance, let's get the `rms` between the ground truth and the predicted values:

rms(y[test], ŷ(rows=test))

# ### Modifying hyperparameters
#
# Hyperparameters can be accessed using the dot syntax as usual.
# Let's modify the regularisation parameter of the ridge regression:

ridge_model.lambda = 5.0;

# Since the node `ẑ` corresponds to a machine that wraps `ridge_model`, that node has effectively changed and will be retrained:

fit!(ŷ, rows=train)
rms(y[test], ŷ(rows=test))

# ## "Arrow" syntax
# **Important**: for this to work, you need to be using **Julia ≥ 1.3**:
#
# The syntax to define nodes etc. is a bit verbose. MLJ supports a shorter syntax which abstracts away some of the steps. We will refer to it as the "arrow" syntax as it makes use of the `|>` operator which can be interpreted as "data flow".
#
# Let's start with `W` and `z` (the "first layer"):

W = X |> Standardizer()
z = y |> UnivariateBoxCoxTransformer()

# Note that we feed `X` and `y` directly into models. In the background, MLJ will create source nodes and assumes that the operation is a `transform` given the models are unsupervised.
#
# For a node that corresponds to a supervised model, you can feed a tuple where the first element corresponds to the input (here `W`) and the second corresponds to the target (here `z`), MLJ will assume the operation is a `predict`:

ẑ = (W, z) |> RidgeRegressor(lambda=0.1);

# Finally we need to apply the inverse of the transform encapsulated in the node `z`, for this:

ŷ = ẑ |> inverse_transform(z);

# That's it! You can now fit the network as before:

fit!(ŷ, rows=train)
rms(y[test], ŷ(rows=test))

# To *manually* modify hyperparameters on a node, you can access them like so:

ẑ[:lambda] = 5.0;

# Here remember that `ẑ` is a node with a machine that wraps around a ridge regression with a parameter `lambda` so the syntax above is equivalent to

ẑ.machine.model.lambda = 5.0;

# which is relevant if you want to tune the hyperparameter using a `TunedModel`.

fit!(ŷ, rows=train)
rms(y[test], ŷ(rows=test))

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
using JLD2, FileIO
using DataFrames
using Unicode

###### structures #################

mutable struct DataConfig
  paths::DefaultDict{AbstractString,String}
  sources::DefaultDict{AbstractString,String}
  loaders::DefaultDict{AbstractString,String}
  data_dirs::DefaultDict{AbstractString,String}
  args_filename_builder::DefaultDict{AbstractString,Any}
  args_data_loader::DefaultDict{AbstractString,Any}
  args_population_loader::DefaultDict{AbstractString,Any}
  data_maps::DefaultDict{AbstractString,Matrix{<:Real}}
  outbreak_conditions::DefaultDict{AbstractString,Function}
    function DataConfig()
      new(
        DefaultDict{AbstractString,String}(k -> k; passkey=true), # paths
        DefaultDict{AbstractString,String}("native"), # sources
        DefaultDict{AbstractString,String}("native"), # loaders
        DefaultDict{AbstractString,String}(""), # data_dirs
        DefaultDict{AbstractString,Any}(k -> k; passkey=true), # args_filename_builder
        DefaultDict{AbstractString,Any}(k -> k; passkey=true), # args_data_loader
        DefaultDict{AbstractString,Any}(k -> k; passkey=true), # args_population_loader
        DefaultDict{AbstractString,Matrix{<:Real}}(Diagonal(ones(Int,4))), # data_maps
        DefaultDict{AbstractString,Function}(c -> findfirst(values(c[:,1].>0))), # outbreak_conditions
      )
    end
end

struct DataStruct{T<:Integer,D<:TimeType}
    cases::TimeArray{T,2,D}
    outbreakdate::Union{D,Missing}
    lastdate::D
    population::T
    map::Matrix{<:Real}
    source::String
    function DataStruct(cases::TimeArray{T,2,D},
                        outbreakdate::Union{D,Missing},
                        lastdate::D,
                        population::T,
                        datamap::Matrix{<:Real},
                        source::AbstractString) where {T<:Integer,D<:TimeType}
        new{T,D}(cases,outbreakdate,lastdate,population,datamap,source)
    end
end

###### filename builder #################
function build_filename(config::DataConfig,region::AbstractString)
    _build_filename(config.args_filename_builder[region])
end

function _build_filename(region::NTuple{N,AbstractString}) where {N}
    fn = _build_filename(region[1])
    for i=2:N
        fn = string(_build_filename(region[i]),"-",fn)
    end
    fn
end

_build_filename(region::NTuple{2,AbstractString}) =
        _build_filename(region[2])*"-"*_build_filename(region[1])

function _build_filename(region::AbstractString)
    char_map = Dict('ä' => "ae", 'ö' => "oe", 'ü' => "ue", 'ß' => "ss", ' ' => "_")
    fn = region
    fn = Unicode.normalize(fn,compat=true,stripcc=true,rejectna=true)
    fn = join(haskey(char_map,c) ? char_map[c] : c for c in fn)
    fn = titlecase(fn,strict=false)
    fn = Unicode.normalize(fn,stripmark=true)
    fn = replace(fn, r"[\W]" => "")
    return fn
end

###### data loaders #################

function load_data(config::DataConfig,region::AbstractString)
    source = config.sources[region]
    path = config.paths["raw"]
    loader = config.loaders[source]
    cases = _load_data(Val(Symbol(loader)),config.args_data_loader[region],
                        joinpath(path,config.data_dirs[source]))
    outbreakdate = get_outbreak_date(config,cases,loader)
    lastdate = timestamp(cases)[end]
    population = _load_population(config.args_population_loader[region],path)
    datamap = config.data_maps[loader]
    DataStruct(cases,outbreakdate,lastdate,population,datamap,source)
end

include("data_loaders/native.jl")
include("data_loaders/population.jl")
include("data_loaders/JHU.jl")

###### outbreak date #################
function get_outbreak_date(config::DataConfig,
                           cases::TimeArray,
                           loader::AbstractString)
    outbreakindex = config.outbreak_conditions[loader](cases)
    if isnothing(outbreakindex)
        outbreakdate = missing
    else
        outbreakdate = timestamp(cases)[outbreakindex]
    end
    return outbreakdate
end

###### save/load results #################

function save(config::DataConfig, region::AbstractString, da::DA;
              interactive=false::Bool, kwargs...)
    path = joinpath(config.paths["results"],build_filename(config,region))
    mkpath(path)
    filename = joinpath(path,"da.jld2")
    if interactive && isfile(filename)
        response = ""
        while isempty(response)
            response = input("overwrite `$filename'? ")
        end
        response == "y" || return
    end
    _save(filename, da; kwargs...)
end

function _save(filename::AbstractString, da::DA;
               data::Union{Missing,DataStruct}=missing)
    dat = ["da", da]
    if !ismissing(data)
        push!(dat, "data", data)
    end
    FileIO.save( filename, dat...)
end

function load_model_params(config::DataConfig,region::AbstractString)
    filename = joinpath(config.paths["results"],build_filename(config,region),"da.jld2")
    da = FileIO.load(filename,"da")
    da.p
end


text:
function normal_clusters(nx, ny; dim=2, σ=0.1)
    x0 = poisson(nx; dim)
    λ = ny / nx
    y0 = clusters(x0, () -> normal(rand(Poisson(λ)); dim) * σ)
    x0, y0
end

function circle_clusters(nx, ny; radius=0.05)
    λ = ny / nx
    x0 = poisson(nx)
    y0 = clusters(x0, () -> radius * circle_poisson(rand(Poisson(λ))))
    x0, y0
end

function anarchy(n)
    z = poisson(n)
    ix = map(w -> w[1] > w[2], eachcol(z))
    x0 = z[:, ix]
    y0 = z[:, .!ix]
    x0, y0
end

function poisson(nx, ny; dim=2)
    x0 = poisson(nx; dim)
    y0 = poisson(ny; dim)
    x0, y0
end

function colored_clusters(nsup, npp; σ=0.1)
    dim = 2
    z = poisson(4nsup; dim)
    d = pairwise(PeriodicEuclidean(ones(dim)), z; dims=2)
    mix = map(i -> minimum(d[:, (1:end) .!= i]), 1:4nsup)
    ix = sortperm(mix)
    z = z[:, 2nsup + 1:end]
    xs = z[:, 1:2:end]
    ys = z[:, 2:2:end]
    λ = npp / nsup
    x0 = clusters(xs, () -> normal(rand(Poisson(λ)); dim) * σ)
    y0 = clusters(ys, () -> normal(rand(Poisson(λ)); dim) * σ)
    x0, y0
end

function target_circles(n=1500, radius=0.125)
    dim = 2
    z = poisson(n; dim)
    center = reshape([0.5, 0.5], 2, 1)
    d = pairwise(PeriodicEuclidean(ones(2)), z, center; dims=dim)[:]
    ix = mod.(div.(d, radius), 2) .== 1
    x0 = z[:, ix]
    y0 = z[:, .!ix]
    x0, y0
end


text:
@testset "Size mutations" begin

    inpt(size, id="in") = inputvertex(id, size)
    nt(name) = t -> NamedTrait(t, name)
    tf(name) = t -> nt(name)(SizeChangeValidation(t))
    cc(in...; name="cc") = conc(in..., dims=1, traitdecoration = tf(name))
    ea(ins...; name="ea") = +(traitconf(tf(name)) >> ins[1], ins[2:end]...)

    @testset "Basic tests" begin

        av(outsize, in, name="av") = absorbvertex(identity, outsize, in, traitdecoration=tf(name))

        @testset "AbsorbVertex" begin
            iv = av(2, inpt(2), "iv")
            v1 = av(3, iv, "v1")

            @test nin(v1) == [2]
            @test nout(v1) == 3

            @test outputs.(inputs(v1)) == [[v1]]

            Δnin(v1, 2)
            @test nin(v1) == [4]

            Δnin(v1, -3)
            @test nin(v1) == [1]

            Δnout(v1, -2)
            @test nout(v1) == 1

            # Add one vertex and see that change propagates
            v2 = av(4, v1, "v2")
            @test nout(v1) == nin(v2)[1]
            @test outputs(v1) == [v2]
            @test inputs(v2) == [v1]

            Δnin(v2, 4)
            @test nout(v1) == nin(v2)[1] == 5

            Δnout(v1, -2)
            @test nout(v1) == nin(v2)[1] == 3

            # Fork of v1 into a new vertex
            v3 = av(2, v1, "v3")
            @test outputs(v1) == [v2, v3]
            @test inputs(v3) == inputs(v2) == [v1]

            Δnout(v1, -2)
            @test nout(v1) == nin(v2)[1] == nin(v3)[1] == 1

            Δnin(v3, 3)
            @test nout(v1) == nin(v2)[1] == nin(v3)[1] == 4

            Δnin(v2, -2)
            @test nout(v1) == nin(v2)[1] == nin(v3)[1] == 2

            ivc = clone(iv, inputs(iv)[1])
            v1c = clone(v1, ivc)
            v2c = clone(v2, v1c)
            v3c = clone(v3, v1c)

            @test issame(iv, ivc)
            @test issame(v1, v1c)
            @test issame(v2, v2c)
            @test issame(v3, v3c)
        end

        @testset "StackingVertex" begin

            @testset "StackingVertex 1 to 1" begin

                iv = av(2, inpt(2), "iv")
                tv = cc(iv)
                io = av(2, tv, "io")

                @test outputs.(inputs(io)) == [[io]]
                @test outputs.(inputs(tv)) == [[tv]]
                @test outputs(iv) == [tv]
                @test neighbours(Both(), tv) == [iv, io]

                Δnout(iv, 3)
                @test [nout(iv)] == nin(tv) == [nout(tv)] == nin(io) == [5]

                Δnin(io, -2)
                @test [nout(iv)] == nin(tv) == [nout(tv)] == nin(io) == [3]

                Δnin(tv, +1)
                @test [nout(iv)] == nin(tv) == [nout(tv)] == nin(io) == [4]

                Δnout(tv, -1)
                @test [nout(iv)] == nin(tv) == [nout(tv)] == nin(io) == [3]

                ivc = clone(iv, inputs(iv)[])
                tvc = clone(tv, ivc)
                ioc = clone(io, tvc)

                @test issame(iv, ivc)
                @test issame(tv, tvc)
                @test issame(io, ioc)
            end

            @testset "StackingVertex 2 inputs" begin
                # Try with two inputs to StackingVertex
                iv1 = av(2, inpt(2,"in1"), "iv1")
                iv2 = av(3, inpt(3,"in2"), "iv2")
                tv = cc(iv1, iv2)
                io1 = av(5, tv, "io1")

                @test inputs(tv) == [iv1, iv2]
                @test outputs.(inputs(tv)) == [[tv], [tv]]
                @test outputs(iv1) == [tv]
                @test outputs(iv2) == [tv]
                @test neighbours(Both(), tv) == [iv1, iv2, io1]

                @test [nout(iv1) + nout(iv2)] == [sum(nin(tv))] == [nout(tv)] == nin(io1) == [5]

                Δnout(iv2, -2)
                @test nin(io1) == [nout(tv)] == [sum(nin(tv))] == [nout(iv1) + nout(iv2)]  ==  [3]

                Δnin(io1, 3)
                @test [nout(iv1) + nout(iv2)] == [sum(nin(tv))] == [nout(tv)] == nin(io1) == [6]

                Δnout(iv1, 4)
                Δnin(io1, -8)
                @test [nout(iv1) + nout(iv2)] == [sum(nin(tv))] == [nout(tv)] == nin(io1) == [2]
                @test nout(iv1) == nout(iv2) == 1

                #Add another output
                io2 = av(2, tv, "io2")

                @test [nout(iv1) + nout(iv2)] == [sum(nin(tv))] == [nout(tv)] == nin(io1) == nin(io2) == [2]

                Δnout(iv1, 3)
                @test [nout(iv1) + nout(iv2)] == [sum(nin(tv))] == [nout(tv)] == nin(io1) == nin(io2) == [5]

                Δnin(io2, -2)
                @test [nout(iv1) + nout(iv2)] == [sum(nin(tv))] == [nout(tv)] == nin(io1) == nin(io2) == [3]

                Δnin(io1, 3)
                @test [nout(iv1) + nout(iv2)] == [sum(nin(tv))] == [nout(tv)] == nin(io1) == nin(io2) == [6]

                Δnin(tv, -1, missing)
                @test [nout(iv1) + nout(iv2)] == [sum(nin(tv))] == [nout(tv)] == nin(io1) == nin(io2) == [5]

                Δnin(tv, missing, 2)
                @test [nout(iv1) + nout(iv2)] == [sum(nin(tv))] == [nout(tv)] == nin(io1) == nin(io2) == [7]

                Δnout(tv, -1)
                @test [nout(iv1) + nout(iv2)] == [sum(nin(tv))] == [nout(tv)] == nin(io1) == nin(io2) == [6]

                iv1c = clone(iv1, inputs(iv1)[])
                iv2c = clone(iv2, inputs(iv2)[])
                tvc = clone(tv, iv1c, iv2c)
                io1c = clone(io1, tvc)
                io2c = clone(io2, tvc)

                @test issame(iv1, iv1c)
                @test issame(iv2, iv2c)
                @test issame(tv, tvc)
                @test issame(io1, io1c)
                @test issame(io2, io2c)
            end

        end
        @testset "InvariantVertex" begin

            @testset "InvariantVertex 1 to 1" begin
                #Behaviour is identical to StackingVertex in this case
                iv = av(2, inpt(2), "iv")
                inv = ea(iv)
                io = av(2, inv, "io")

                @test outputs.(inputs(io)) == [[io]]
                @test outputs.(inputs(inv)) == [[inv]]
                @test outputs(iv) == [inv]

                Δnout(iv, 3)
                @test [nout(iv)] == nin(inv) == [nout(inv)] == nin(io) == [5]

                Δnin(io, -2)
                @test [nout(iv)] == nin(inv) == [nout(inv)] == nin(io) == [3]

                ivc = clone(iv, inputs(iv)[])
                invc = clone(inv, ivc)
                ioc = clone(io, invc)

                @test issame(iv, ivc)
                @test issame(inv, invc)
                @test issame(io, ioc)
            end

            @testset "invariantvertex 2 inputs" begin
                # Try with two inputs a eavariantVertex
                iv1 = av(3, inpt(2, "in1"), "iv1")
                iv2 = av(3, inpt(2, "in2"), "iv2")
                eav = ea(iv1, iv2)
                io1 = av(5, eav, "io1")

                @test inputs(eav) == [iv1, iv2]
                @test outputs.(inputs(eav)) == [[eav], [eav]]
                @test outputs(iv1) == [eav]
                @test outputs(iv2) == [eav]

                @test nout(iv1) == nout(iv2) == nin(eav)[1] == nin(eav)[2] == nout(eav) == nin(io1)[] == 3

                Δnout(iv2, -2)
                @test nin(io1)[] == nout(eav) == nin(eav)[1] == nin(eav)[2]  == nout(iv1) == nout(iv2)  ==  1

                Δnin(io1, 3)
                @test nout(iv1) == nout(iv2) == nin(eav)[1] == nin(eav)[2] == nout(eav) == nin(io1)[] == 4

                #Add another output
                io2 = av(4, eav, "io2")

                @test nout(iv1) == nout(iv2) == nin(eav)[1] == nin(eav)[2] == nout(eav) == nin(io1)[] == nin(io2)[] == 4

                Δnout(iv1, -3)
                @test nout(iv1) == nout(iv2) == nin(eav)[1] == nin(eav)[2] == nout(eav) == nin(io1)[] == nin(io2)[] == 1

                Δnin(io2, 2)
                @test nout(iv1) == nout(iv2) == nin(eav)[1] == nin(eav)[2] == nout(eav) == nin(io1)[] == nin(io2)[] == 3

                Δnout(iv2, 3)
                @test nout(iv1) == nout(iv2) == nin(eav)[1] == nin(eav)[2] == nout(eav) == nin(io1)[] == nin(io2)[] == 6

                iv1c = clone(iv1, inputs(iv1)[])
                iv2c = clone(iv2, inputs(iv2)[])
                eavc = clone(eav, iv1c, iv2c)
                io1c = clone(io1, eavc)
                io2c = clone(io2, eavc)

                @test issame(iv1, iv1c)
                @test issame(iv2, iv2c)
                @test issame(eav, eavc)
                @test issame(io1, io1c)
                @test issame(io2, io2c)
            end
        end

        @testset "SizeChangeValidation" begin

            # Mock for creating invalid size changes
            import NaiveNASlib:OnlyFor
            struct Ignore <: MutationSizeTrait end
            function NaiveNASlib.Δnin(::OnlyFor, ::Ignore, v, Δ::Union{Missing, T}...) where T end
            function NaiveNASlib.Δnout(::OnlyFor, ::Ignore, v, Δ) end
            NaiveNASlib.all_in_Δsize_graph(::Ignore, d::Direction, v, visited) = all_in_Δsize_graph(SizeAbsorb(), d, v, visited)
            function NaiveNASlib.ninconstraint!(s, ::Ignore, v, data) end

            v1 = inpt(3, "v1")
            v2 = vertex(identity, 5, Ignore() |> tf("v2"), v1)
            v3 = av(7, v2, "v3")
            v4 = vertex(identity, 5, Ignore() |> tf("v4"), v3)

            @test_throws ArgumentError Δnout(v2, 2)
            @test_throws ArgumentError Δnin(v3, -3)
            @test_throws ArgumentError Δnin(v4, 5)
            @test_throws ArgumentError Δnout(v3, -7)

            # Too many Δs!
            @test_throws AssertionError Δnin(v3, 1, 1)
        end

        @testset "Fail mutate immutable" begin
            v0 = inpt(3)
            v1 = av(5, v0, "v1")

            @test_throws ErrorException Δsize(ΔNout{Exact}(v0, 2, ΔSizeFailError("Success?!")), [v0, v1])
            @test [nout(v0)] == nin(v1) == [3]

            @test_throws ErrorException Δsize(ΔNin{Exact}(v1, [-2], ΔSizeFailError("Success?!")), [v1, v0])
            @test [nout(v0)] == nin(v1) == [3]
        end
    end

    @testset "Mutate tricky structures" begin

        ## Helper functions
        rb(start, residual,name="add") = traitconf(tf(name)) >> start + residual
        iv(in; name="iv") = invariantvertex(identity, in, traitdecoration=tf(name))
        concd1(paths...;name="conc") = conc(paths..., dims=1, traitdecoration=tf(name))
        mm(nin, nout) = x -> x * reshape(collect(1:nin*nout), nin, nout)
        av(in, outsize, name="comp") = absorbvertex(mm(nout(in), outsize), outsize, in, traitdecoration = tf(name))
        function stack(start, nouts...; bname = "stack")
            # Can be done on one line with mapfoldl, but it is not pretty...
            next = start
            for i in 1:length(nouts)
                next = av(next, nouts[i], "$(bname)_$i")
            end
            return next
        end

        @testset "Residual fork block" begin
            start = av(inputvertex("in", 3), 9, "start")
            p1 = stack(start, 3,4, bname = "p1")
            p2 = stack(start, 4,5, bname = "p2")
            resout = rb(start, concd1(p1, p2))
            out = av(resout, 9, "out")

            @test nout(resout) == 9

            # Propagates to out, start outputs of start and also to p1 as objective is minimized by increasing p1 by 1
            Δnout(p2, -2)
            @test nout(p2) == 3
            @test nin(out) == [nout(start)] == [8]
            @test nout(p2) + nout(p1) == 8
            #outputs(start) = first vertex in p1, p2 and resout (which has two inputs)
            @test foldl(vcat, nin.(outputs(start))) == [8, 8, 8, 8]

            Δnin(out, +2)
            @test nin(out) == [nout(start)] == [10]
            @test foldl(vcat, nin.(outputs(start))) == [10, 10, 10, 10]
        end

        @testset "Half transparent residual fork block" begin
            start = av(inputvertex("in", 3), 8, "start")
            split = av(start, 4, "split")
            p1 = iv(split, name="p1") #Just an identity vertex
            p2 = stack(split, 3,2,4, bname="p2")
            resout = rb(start, concd1(p1, p2))
            out = av(resout, 3, "out")

            @test nout(resout) == 8

            # Propagates to input of first vertex of p2, input of out and start
            # via p1 and resout as well as to input of split
            Δnout(split, -1)
            @test nin(out) == [nout(start)] == nin(split) == [8]
            @test foldl(vcat, nin.(outputs(split))) == [3, 3]

            # Should basically undo the previous mutation
            Δnin(out, +1)
            @test nin(out) == [nout(start)] == nin(split) == [9]
        end

        @testset "Transparent fork block" begin
            start = av(inputvertex("in", 3), 4, "start")
            p1 = iv(start, name="p1")
            p2 = iv(start, name="p2")
            joined = concd1(p1, p2, name="join")
            out = av(joined, 3, "out")

            @test nout(joined) == 8

            # Evil action: This will propagate to both p1 and p2 which are in
            # turn both input to the conc before resout. Simple dfs will
            # fail as one will hit the conc through p1 before having
            # resolved the path through p2.
            Δnout(start, -1)
            @test nin(out) == [2nout(start)] == [6]

            # Should basically undo the previous mutation
            @test minΔninfactor(out) == 2
            Δnin(out, +2)
            @test nin(out) == [2nout(start)] == [8]
        end

        @testset "Transparent residual fork block" begin
            start = av(inputvertex("in", 3), 8, "start")
            split = av(start, 4, "split")
            p1 = iv(split, name="p1")
            p2 = iv(split, name="p2")
            resout = rb(start, concd1(p1, p2, name="join"))
            out = av(resout, 3, "out")

            @test nout(resout) == 8

            # Evil action: This will propagate to both p1 and p2 which are in
            # turn both input to the conc before resout. Simple dfs will
            # fail as one will hit the conc through p1 before having
            # resolved the path through p2.
            Δnout(split, -1)
            @test nin(out) == [nout(start)] == nin(split) == [6]

            # Should basically undo the previous mutation
            Δnin(out, +2)
            @test nin(out) == [nout(start)] == nin(split) == [8]

            @test minΔninfactor(out) == 2
            Δnout(start, -2)
            @test nin(out) == [nout(start)] == [6]
            @test nout(split) == 3
        end

        @testset "Transparent residual fork block with single absorbing path" begin
            start = av(inputvertex("in", 3), 8, "start")
            split = av(start, 3, "split")
            p1 = iv(split, name="p1")
            p2 = iv(split, name="p2")
            p3 = av(split, 2, "p3")
            resout = rb(start, concd1(p1, p2, p3, name="join"), "add")
            out = av(resout, 3, "out")

            @test nout(resout) == 8

            # Evil action: This will propagate to both p1 and p2 which are in
            # turn both input to the conc before resout. Simple dfs will
            # fail as one will hit the conc through p1 before having
            # resolved the path through p2.
            Δnout(split, -1)
            @test nin(out) == [nout(start)] == nin(split) == [6]

            # Should basically undo the previous mutation
            Δnin(out, +2)
            @test nin(out) == [nout(start)] == nin(split) == [8]

            @test minΔninfactor(out) == 2
            Δnout(start, -2)
            @test nin(out) == [nout(start)] == [2*nout(split) + nout(p3)] == [6]
            @test nout(split) == 2
            @test nout(p3) == 2
        end
    end

    @testset "Size Mutation possibilities" begin

        # Helpers
        struct SizeConstraint constraint; end
        NaiveNASlib.minΔnoutfactor(c::SizeConstraint) = c.constraint
        NaiveNASlib.minΔninfactor(c::SizeConstraint) = c.constraint
        function NaiveNASlib.compconstraint!(s, c::SizeConstraint, data)
            fv_out = JuMP.@variable(data.model, integer=true)
            JuMP.@constraint(data.model, c.constraint * fv_out ==  nout(data.vertex) - data.noutdict[data.vertex])

            ins = filter(vin -> vin in keys(data.noutdict), inputs(data.vertex))
            fv_in = JuMP.@variable(data.model, [1:length(ins)], integer=true)
            JuMP.@constraint(data.model, [i=1:length(ins)], c.constraint * fv_in[i] ==  nout(ins[i]) - data.noutdict[ins[i]])
        end
        av(size, csize, in... ;name = "av") = absorbvertex(SizeConstraint(csize), size, in..., traitdecoration=tf(name))

        @testset "SizeStack multi inputs" begin
            v1 = av(6,3, inpt(3), name="v1")
            v2 = av(7,2, inpt(3), name="v2")
            v3 = av(8,2, inpt(3), name="v3")

            cc1 = cc(v1, v2, name="cc1")
            cc2 = cc(v3, v2, v1, v2, name="cc2")
            @test minΔnoutfactor(cc1) == 6
            @test minΔnoutfactor(cc2) == 12

            # Expect only v1 to change as size change is not compatible with v2
            Δnout(cc1, -3)
            @test nin(cc1) == [nout(v1), nout(v2)] == [3, 7]
            @test nin(cc2) == [nout(v3), nout(v2), nout(v1), nout(v2)] == [8, 7, 3, 7]
            @test nout(cc1) == sum(nin(cc1))
            @test nout(cc2) == sum(nin(cc2))

            # v1 can't change as it is too small already
            Δnout(cc2, -6)
            @test nin(cc1) == [nout(v1), nout(v2)] == [3, 5]
            @test nin(cc2) == [nout(v3), nout(v2), nout(v1), nout(v2)] == [6, 5, 3, 5]
            @test nout(cc1) == sum(nin(cc1))
            @test nout(cc2) == sum(nin(cc2))

            # Evil action! Must have understanding that change will propagate to cc1 from
            # both v1 and v2
            Δnout(cc2, +9)
            @test nin(cc1) == [nout(v1), nout(v2)] == [6, 7]
            @test nin(cc2) == [nout(v3), nout(v2), nout(v1), nout(v2)] == [8, 7, 6, 7]
            @test nout(cc1) == sum(nin(cc1))
            @test nout(cc2) == sum(nin(cc2))

        end

        @testset "Stacked SizeStacks" begin
            v1 = av(100,1, inpt(3), name="v1")
            v2 = av(100,2, inpt(3), name="v2")
            v3 = av(100,3, inpt(3), name="v3")
            v4 = av(100,5, inpt(3), name="v4")

            cc1 = cc(v2, v3, name="cc1")
            cc2 = cc(cc1, v1, name="cc2")
            cc3 = cc(cc2, v4, v2, name="cc3")
            @test minΔnoutfactor(cc1) == minΔninfactor(cc1) == 6
            @test minΔnoutfactor(cc2) == minΔninfactor(cc2)== 6
            @test minΔnoutfactor(cc3) == minΔninfactor(cc3) == 60 # v2 is input twice through sc2->cc1

            Δnout(cc3, -60)
            @test nin(cc1) == nout.(inputs(cc1)) == [86, 91]
            @test nin(cc2) == nout.(inputs(cc2)) == [177, 87]
            @test nin(cc3) == nout.(inputs(cc3)) == [264, 90, 86]
            @test nout(cc1) == nout(v2) + nout(v3) == 177
            @test nout(cc2) == nout(cc1) + nout(v1) == 264
            @test nout(cc3) == nout(cc2) + nout(v4) + nout(v2) == 440

            v5 = av(10, 3, cc3, name="v5")
            @test minΔnoutfactor(cc1) == minΔninfactor(cc1) == 6
            @test minΔnoutfactor(cc2) == minΔninfactor(cc2)== 6
            @test minΔnoutfactor(cc3) == minΔninfactor(cc3) == 60

            Δnout(v1, 3)
            @test nin(cc1) == nout.(inputs(cc1)) == [86, 91]
            @test nin(cc2) == nout.(inputs(cc2)) == [177, 90]
            @test nin(cc3) == nout.(inputs(cc3)) == [267, 90, 86]
            @test nout(cc1) == nout(v2) + nout(v3) == 177
            @test nout(cc2) == nout(cc1) + nout(v1) == 267
            @test [nout(cc3)] == [(nout(cc2) + nout(v4) + nout(v2))] == nin(v5) == [443]

            # Evil action! Must have understanding that the change in v2 will propagate
            # to cc3 input 1 through cc1 and cc2
            Δnout(v2, -12)
            @test nin(cc1) == nout.(inputs(cc1)) == [74, 103]
            @test nin(cc2) == nout.(inputs(cc2)) == [177, 90]
            @test nin(cc3) == nout.(inputs(cc3)) == [267, 90, 74]
            @test nout(cc1) == nout(v2) + nout(v3) == 177
            @test nout(cc2) == nout(cc1) + nout(v1) == 267
            @test [nout(cc3)] == [(nout(cc2) + nout(v4) + nout(v2))] == nin(v5) == [431]
        end

        @testset "SizeInvariant multi input" begin
            v1 = av(10,1, inpt(3), name="v1")

            ea1 = ea(v1, name="ea1")
            ea2 = ea(ea1, v1, name="ea2")

            Δnout(ea2, -2)

            @test nout(ea2) == nout(ea1) == nout(v1) == 8
            @test nin(ea2) == [nout(ea1), nout(v1)] == [8, 8]
            @test nin(ea1) == [nout(ea1)] == [8]
        end

        @testset "SizeInvariant multi SizeStack input" begin
            v1 = av(7,1, inpt(3), name="v1")
            v2 = av(13,1, inpt(3), name="v2")

            cc1 = cc(v1, v2, name="cc1")
            cc2 = cc(v1, v2, name="cc2")

            ea1 = ea(cc1, cc2, name="ea1")
            ea2 = ea(cc1, cc2, name="ea2")
            ea3 = ea(ea1, ea2, name="ea3")

            Δnout(ea3, -4)

            @test nout(ea3) == nout(ea2) == nout(ea1) == nout(cc2) == nout(cc1) == nout(v2) + nout(v1) == 16
            @test nin(ea3) == nin(ea2) == nin(ea3) == [nout(cc1), nout(cc2)] == [16, 16]
            @test nin(cc2) == nin(cc1) == [nout(v1), nout(v2)] == [6, 10]
        end

        @testset "Stacked SizeInvariant" begin
            v1 = av(100,1, inpt(3), name="v1")
            v2 = av(100,2, inpt(3), name="v2")
            v3 = av(100,3, inpt(3), name="v3")
            v4 = av(100,5, inpt(3), name="v4")

            ea1 = ea(v2, v3, name="ea1")
            ea2 = ea(ea1, v1, name="ea2")
            ea3 = ea(ea2, v4, v2, name="ea3")

            # Everything touches everything in this setup
            @test minΔnoutfactor(ea1) == minΔninfactor(ea1) == 1*2*3*5
            @test minΔnoutfactor(ea2) == minΔninfactor(ea2) == 1*2*3*5
            @test minΔnoutfactor(ea3) == minΔninfactor(ea3) == 1*2*3*5

            Δnout(ea3, -30)
            @test nout(v1) == nout(v2) == nout(v3) == nout(v4) == 70
            @test nout(ea1) == nout(ea2) == nout(ea3) == 70

            v5 = av(10, 3, ea3, name="v5")

            Δnout(v1, 60)
            @test nout(v1) == nout(v2) == nout(v3) == nout(v4) == 130
            @test [nout(ea1)] == [nout(ea2)] == [nout(ea3)] == nin(v5) == [130]

            # Evil action! Must have understanding that the change in v2 will propagate
            # to ea3 input 1 through ea1 and ea2 and hold off updating it through input 3
            Δnout(v2, -30)
            @test nout(v1) == nout(v2) == nout(v3) == nout(v4) == 100
            @test [nout(ea1)] == [nout(ea2)] == [nout(ea3)] == nin(v5) == [100]
        end

        @testset "Stacked input to SizeInvariant" begin
            v0 = inpt(3, "in")
            v1 = av(20, 2, v0, name="v1")
            v2 = av(30, 3, v0, name="v2")
            v3 = cc(v1,v2,v1, name="v3")
            v4 = av(70, 5, v0, name="v4")
            v5 = ea(v4,v3, name="v5")
            v6 = av(10, 7, v5, name="v6")

            # Evilness: Invariant vertex must change all its inputs and therefore it must take their minΔnoutfactors into account when computing minΔninfactor.
            @test minΔnoutfactor(v4) == 2*2*3*5*7

            Δnout(v5, 2*2*3*5*7)

            @test nout(v5) == nout(v4) == nout(v3) == 490
            @test nin(v5) == [nout(v4), nout(v3)] == [490, 490]
            @test nin(v3) == [nout(v1), nout(v2), nout(v1)] == [140, 210, 140]
        end

        @testset "Invariant input to SizeInvariant" begin
            v0 = inpt(3, "in")
            v1 = av(20, 2, v0, name="v1")
            v2 = av(20, 3, v0, name="v2")
            v3 = ea(v1,v2, name="v3")
            v4 = av(20, 5, v0, name="v4")
            v5 = ea(v4,v3, name="v5")
            v6 = av(10, 7, v5, name="v6")

            # Evilness: Infinite recursion without memoization
            @test minΔnoutfactor(v4) == 2*3*5*7

            Δnout(v5, 2*3*5*7)

            @test nout(v5) == nout(v4) == nout(v3) == 230
            @test nin(v5) == [nout(v4), nout(v3)] == nin(v3) == [nout(v1), nout(v2)] == [230, 230]
        end

        @testset "SizeStack duplicate through SizeInvariant" begin
            v0 = inpt(5, "in")
            v1 = av(2, 2, v0, name="v1")
            v2 = ea(v1, name="v2")
            v3 = ea(v1, name="v3")
            v4 = cc(v2,v3, name="v4")
            v5 = ea(v4, name="v5")
            v6 = cc(v5, v2, name="v6")

            @test minΔnoutfactor(v6) == 3*2
            Δnout(v6, 12)

            @test nout(v6) == sum(nin(v6)) == 18
            @test nout(v5) == nin(v5)[] == 12
            @test nout(v4) == sum(nin(v4)) == 12
            @test nout(v1) == 6
        end

        @testset "SizeInvariant zig-zag" begin
                v0 = inpt(5, "in")
                v1 = av(2, 2, v0, name="v1")
                v2 = av(3, 3, v0, name="v2")
                v3 = cc(v1,v2, name="v3")
                function zigzag(vin1, sc, vin2=v0;name="zig")
                    vnew = av(nout(vin1), sc, vin2, name=name*"_new")
                    vout = ea(vnew, name=name*"_ivB")
                    vcon = ea(vout, vin1, name=name*"_ivA")
                    return vout
                end
                v4 = zigzag(v3, 5, name="z1")
                v5 = zigzag(v4, 7, name="z2")
                v6 = av(11, 11, v5, name="v6")

                expectedΔf = 2*3*5*7*11
                @test minΔnoutfactor(v1) == expectedΔf / 3 # 3 is size constraint for v2
                @test minΔnoutfactor(v2) == expectedΔf / 2 # 2 is size constraint for v1
                @test minΔnoutfactor(v3) == expectedΔf

                Δnout(v2, expectedΔf)
                @test nout(v3) == nout(v4) == nout(v5) == expectedΔf + 2+3
                @test nin(v4) == nout.(inputs(v4)) == [expectedΔf + 2+3]
                @test nin(v5) == nout.(inputs(v5)) == [expectedΔf + 2+3]
        end

        @testset "SizeStack duplicate SizeInvariant mini-zig-zag" begin
            v0 = inpt(5, "in")
            v1 = av(2, 2, v0, name="v1")
            v2 = av(2, 3, v0, name="v2")
            v3 = ea(v1, name="v3")
            v4 = ea(v3, v2, name="v4")
            v5 = av(3, 5, v3, name="v5")
            v6 = cc(v4,v4, name="v6")
            v7 = ea(v6, v6, name="v7")
            v8 = av(4, 7, v7, name="v8")

            expectedΔf = 2*2*3*5*7
            @test minΔnoutfactor(v6) == expectedΔf

            Δnout(v6, expectedΔf)
            @test nin(v8) == [nout(v7)] == [nout(v6)] == [4+expectedΔf]

            @test nin(v6) == [2 + expectedΔf ÷ 2, 2 + expectedΔf ÷ 2]
            @test nout(v1) == nout(v2) == 2 + expectedΔf ÷ 2
            @test nin(v5) == [2 + expectedΔf ÷ 2]
        end

        @testset "Deep SizeStack" begin
            v0 = inpt(3, "in")
            v1 = av(8, 1, v0, name="v1")
            v2 = av(4, 1, v0, name="v2")
            v3 = cc(v1,v2, name= "v3")
            pa1 = ea(v3, name="pa1")
            pb1 = ea(v3, name="pb1")
            pc1 = ea(v3, name="pc1")
            pd1 = av(5, 1, v3, name="pd1")
            pa1pa1 = ea(pa1, name="pa1pa1")
            pa1pb1 = ea(pa1, name="pa1pb1")
            pa2 = cc(pa1pa1, pa1pb1, name = "pa2")
            v4 = cc(pa2, pb1, pc1, pd1, name = "v4")

            @test minΔnoutfactor(v4) == 4

            Δnout(v4, 8)

            @test nout(v4) == 61
            @test nin(v4) == nout.(inputs(v4)) == [28, 14, 14, 5]
            @test nin(pa2) == nout.(inputs(pa2)) == [14, 14]
            @test nin(pa1pa1) == nin(pa1pb1) == [14]
            @test nin(pa1) == nin(pb1) == nin(pc1) == [nout(v3)] == [14]
            @test nin(v3) == [10, 4]
        end

        @testset "Fail invalid size change" begin
            v0 = inpt(3)
            v1 = av(3,3, v0, name="v1")
            v2 = av(5,2, v1, name="v2")
            import NaiveNASlib:Exact, Relaxed

            @test_throws ErrorException newsizes(ΔNout{Exact}(v1, 2, ΔSizeFailError("")), all_in_graph(v1))
            @test_throws ErrorException newsizes(ΔNout{Exact}(v1, 3, ΔSizeFailError("")), all_in_graph(v1))

            @test_throws ErrorException newsizes(ΔNin{Exact}(v2, [2], ΔSizeFailError("")), all_in_graph(v1))
            @test_throws ErrorException newsizes(ΔNin{Exact}(v2, [3], ΔSizeFailError("")), all_in_graph(v1))

            @test newsizes(ΔNout{Exact}(v1, 2, ΔSizeFailNoOp()), all_in_graph(v1))[1] == false
            @test newsizes(ΔNin{Exact}(v2, [2], ΔSizeFailNoOp()), all_in_graph(v2))[1] == false

            using Logging
            @test (@test_logs (:info, "Giving up") newsizes(ΔNout{Exact}(v1, 2, LogΔSizeExec(Logging.Info, "Giving up")), all_in_graph(v1)))[1] == false

            @test_logs (:warn, r"Could not change nout of .* by 2! Relaxing constraints...") Δnout(v1, 2)
            @test [nout(v1)] == nin(v2) == [9]

            @test_logs (:warn, r"Could not change nin of .* by -2! Relaxing constraints...") Δnin(v2, -2)
            @test [nout(v1)] == nin(v2) == [3]

            v3 = cc(v1,v2, name="v3")

            @test_logs (:warn, r"Could not change nin of .* by 2, 6! Relaxing constraints...") Δnin(v3, 2, 6)
            @test [nout(v1)] == nin(v2) == [9]
            @test nin(v3) == [nout(v1), nout(v2)] == [9, 11]

            @test_logs (:warn, r"Could not change nout of .* by -3! Relaxing constraints...") Δnout(v3, -3)
            @test [nout(v1)] == nin(v2) == [9]
            @test nin(v3) == [nout(v1), nout(v2)] == [9, 9]

        end

    end

    @testset "SizeChangeLogger" begin

        traitfun(name) = t -> SizeChangeLogger(NameInfoStr(), NamedTrait(t, name))
        av(in, size ;name = "av") = absorbvertex(identity, size, in, traitdecoration=traitfun(name))

        @testset "Log size change" begin
                v1 = inpt(3, "v1")
                v2 = av(v1, 10, name="v2")
                v3 = av(v2, 4, name="v3")

                @test_logs (:info, "Change nin of v3 by 3") (:info, "Change nout of v2 by 3") Δnin(v3, 3)

                @test_logs (:info, "Change nout of v2 by -3") (:info, "Change nin of v3 by -3") Δnout(v2, -3)

                v4 = av(v1, nout(v3), name="v4")
                v5 = traitconf(traitfun("v5")) >> v3 + v4

                @test_logs  (:info, "Change nin of v5 by 30, 30") (:info, "Change nout of v5 by 30") (:info, "Change nout of v3 by 30") (:info, "Change nout of v4 by 30") Δnout(v5, 30)

                @test_logs  (:info, "Change nin of v5 by [1, 2, 3, 4, -1×30], [1, 2, 3, 4, -1×30]") (:info, "Change nout of v5 by [1, 2, 3, 4, -1×30]") (:info, "Change nout of v3 by [1, 2, 3, 4, -1×30]") (:info, "Change nout of v4 by [1, 2, 3, 4, -1×30]") Δoutputs(v5, v -> 1:nout_org(v))
        end

        @testset "Clone" begin
            slg = SizeChangeLogger(NameInfoStr(), NamedTrait(SizeAbsorb(), "test"))

            @test slg == clone(slg)

            rename(x;cf=rename) = clone(x, cf=cf)
            rename(x::String; cf) = "success"

            @test name(rename(slg), nothing) == "success"
        end
    end
end


text:
module SlackSDK
using Swagger

export WebClient

include("web/Web.jl")
using .Web

global const apiNames = Ref{Array{String}}()

# Struct encapsulating a connection to the Slack API.
struct WebClient
    token::String
    swagger_client::Swagger.Client
    signature::String

    function WebClient(token=get(ENV, "SLACK_API_TOKEN",
	""); base_url="https://www.slack.com/api/", signature="")
        if isempty(token)
            error("Token is empty!\nA Slack API token is required to perform API operations. Set an environment variable `SLACK_API_TOKEN` or provide directly to the `WebClient` constructor" )
        else
            hdrs = Dict{String,String}("Authorization"=>"Bearer $token")
            return new(token, Swagger.Client(base_url, headers=hdrs), signature)
        end
    end
end

Base.show(io::IO, wc::WebClient) = print(io, "Slack Web API client, with base url '$(wc.swagger_client.root)' and token '$(wc.token[1:8])....$(wc.token[end-1:end])'")

function Base.getproperty(wc::WebClient, method::Symbol)
    if hasfield(WebClient, method)
        return getfield(wc, method)
    else
        nm = apiNames[]
        name_list = split(String(method), "_"; keepempty=false)
        api_object_name = join(map(x->titlecase(x, strict=false), name_list[1:end-1])) * "Api"
        if findfirst(x -> x==api_object_name, nm) === nothing
            throw(UndefVarError(method))
        end
        api_object_type = getfield(SlackSDK.Web, Symbol(api_object_name))
        api_object = api_object_type(wc.swagger_client)
        method_name = name_list[1] * join(map(x->titlecase(x, strict=false), name_list[2:end]))
        try
            method_fn = getfield(SlackSDK.Web, Symbol(method_name))
            return (x...; y...) -> begin
                    method_fn(api_object, x...; y...)
                end
        catch e
            if typeof(e) == UndefVarError
                rethrow(UndefVarError(method))
            else
                rethrow()
            end
        end
    end
end

function Base.propertynames(wc::WebClient)
       return Symbol.(methodNames)
end

function __init__()
    nm = String.(names(SlackSDK.Web))
    apiNames[] = copy(filter(x->endswith(x, "Api"), nm))
end

#=
Generated with the following code. Not doing this at runtime, since reading non-source files are problematic if packages are compiled to sysimages.
Regenerate if the spec file changes.
     using JSON3
     file = read("specs/slack_web_openapi_v2_without_examples.json", String);
     j = JSON3.read(file);
     a = map(x->replace(x[2:end], "."=>"_"), string.(keys(j.paths)));
     print(a)
=#
const methodNames = [
    "admin_apps_approve",
	"admin_apps_approved_list",
	"admin_apps_requests_list",
	"admin_apps_restrict",
	"admin_apps_restricted_list",
	"admin_conversations_archive",
	"admin_conversations_convertToPrivate",
	"admin_conversations_create",
	"admin_conversations_delete",
	"admin_conversations_disconnectShared",
	"admin_conversations_ekm_listOriginalConnectedChannelInfo",
	"admin_conversations_getConversationPrefs",
	"admin_conversations_getTeams",
	"admin_conversations_invite",
	"admin_conversations_rename",
	"admin_conversations_restrictAccess_addGroup",
	"admin_conversations_restrictAccess_listGroups",
	"admin_conversations_restrictAccess_removeGroup",
	"admin_conversations_search",
	"admin_conversations_setConversationPrefs",
	"admin_conversations_setTeams",
	"admin_conversations_unarchive",
	"admin_emoji_add",
	"admin_emoji_addAlias",
	"admin_emoji_list",
	"admin_emoji_remove",
	"admin_emoji_rename",
	"admin_inviteRequests_approve",
	"admin_inviteRequests_approved_list",
	"admin_inviteRequests_denied_list",
	"admin_inviteRequests_deny",
	"admin_inviteRequests_list",
	"admin_teams_admins_list",
	"admin_teams_create",
	"admin_teams_list",
	"admin_teams_owners_list",
	"admin_teams_settings_info",
	"admin_teams_settings_setDefaultChannels",
	"admin_teams_settings_setDescription",
	"admin_teams_settings_setDiscoverability",
	"admin_teams_settings_setIcon",
	"admin_teams_settings_setName",
	"admin_usergroups_addChannels",
	"admin_usergroups_addTeams",
	"admin_usergroups_listChannels",
	"admin_usergroups_removeChannels",
	"admin_users_assign",
	"admin_users_invite",
	"admin_users_list",
	"admin_users_remove",
	"admin_users_session_invalidate",
	"admin_users_session_reset",
	"admin_users_setAdmin",
	"admin_users_setExpiration",
	"admin_users_setOwner",
	"admin_users_setRegular",
	"api_test",
	"apps_event_authorizations_list",
	"apps_permissions_info",
	"apps_permissions_request",
	"apps_permissions_resources_list",
	"apps_permissions_scopes_list",
	"apps_permissions_users_list",
	"apps_permissions_users_request",
	"apps_uninstall",
	"auth_revoke",
	"auth_test",
	"bots_info",
	"calls_add",
	"calls_end",
	"calls_info",
	"calls_participants_add",
	"calls_participants_remove",
	"calls_update",
	"chat_delete",
	"chat_deleteScheduledMessage",
	"chat_getPermalink",
	"chat_meMessage",
	"chat_postEphemeral",
	"chat_postMessage",
	"chat_scheduleMessage",
	"chat_scheduledMessages_list",
	"chat_unfurl",
	"chat_update",
	"conversations_archive",
	"conversations_close",
	"conversations_create",
	"conversations_history",
	"conversations_info",
	"conversations_invite",
	"conversations_join",
	"conversations_kick",
	"conversations_leave",
	"conversations_list",
	"conversations_mark",
	"conversations_members",
	"conversations_open",
	"conversations_rename",
	"conversations_replies",
	"conversations_setPurpose",
	"conversations_setTopic",
	"conversations_unarchive",
	"dialog_open",
	"dnd_endDnd",
	"dnd_endSnooze",
	"dnd_info",
	"dnd_setSnooze",
	"dnd_teamInfo",
	"emoji_list",
	"files_comments_delete",
	"files_delete",
	"files_info",
	"files_list",
	"files_remote_add",
	"files_remote_info",
	"files_remote_list",
	"files_remote_remove",
	"files_remote_share",
	"files_remote_update",
	"files_revokePublicURL",
	"files_sharedPublicURL",
	"files_upload",
	"migration_exchange",
	"oauth_access",
	"oauth_token",
	"oauth_v2_access",
	"pins_add",
	"pins_list",
	"pins_remove",
	"reactions_add",
	"reactions_get",
	"reactions_list",
	"reactions_remove",
	"reminders_add",
	"reminders_complete",
	"reminders_delete",
	"reminders_info",
	"reminders_list",
	"rtm_connect",
	"search_messages",
	"stars_add",
	"stars_list",
	"stars_remove",
	"team_accessLogs",
	"team_billableInfo",
	"team_info",
	"team_integrationLogs",
	"team_profile_get",
	"usergroups_create",
	"usergroups_disable",
	"usergroups_enable",
	"usergroups_list",
	"usergroups_update",
	"usergroups_users_list",
	"usergroups_users_update",
	"users_conversations",
	"users_deletePhoto",
	"users_getPresence",
	"users_identity",
	"users_info",
	"users_list",
	"users_lookupByEmail",
	"users_profile_get",
	"users_profile_set",
	"users_setActive",
	"users_setPhoto",
	"users_setPresence",
	"views_open",
	"views_publish",
	"views_push",
	"views_update",
	"workflows_stepCompleted",
	"workflows_stepFailed",
	"workflows_updateStep"
]

end


text:
module TestUnique

include("preamble.jl")

@testset begin
    @test unique!!([3, 2, 1, 2, 3]) == [3, 2, 1]
    @test unique!!(x for x in [3, 2, 1, 2, 3]) == [3, 2, 1]
end

@testset "mutation" begin
    @testset "Vector" begin
        x = [3, 2, 1, 2, 3]
        y = unique(x)
        @test unique!!(x) === x == y
    end
    @testset "Set" begin
        x = Set([0])
        @test unique!!(x) === x
    end
end

end  # module


text:
export islin, coord, lift

mydot(a, b) = sum(map(*, a, b))
mydot(a::AbstractVector, b::AbstractVector) = dot(f, v)
mydot(f::FixedVector, v::FixedVector) = dot(f, v)
#mydot{T<:Union{Ray,Line}}(a::T, b::T) = mydot(a.r, b.r)
#mydot{T<:Union{Ray,Line}}(a, b::T) = mydot(a, b.r)
#mydot{T<:Union{Ray,Line}}(a::T, b) = mydot(a.r, b)

vecconv{T}(::Type{T}, a::AbstractVector) = AbstractVector{T}(a)
vecconv{T}(::Type{T}, a::FixedVector) = FixedVector{T}(a)

abstract HRepElement{N,T}

# <a, x> <= β
type HalfSpace{N,T} <: HRepElement{N,T}
  a
  β
end

type HyperPlane{N,T} <: HRepElement{N,T}
  a
  β
end
HalfSpace(a, β) = HalfSpace{fulldim(a), eltype(a)}(a, eltype(a)(β))
HyperPlane(a, β) = HyperPlane{fulldim(a), eltype(a)}(a, eltype(a)(β))

Base.convert{N,Tin,Tout}(::Type{HalfSpace{N,Tout}}, h::HalfSpace{N,Tin}) = HalfSpace{N,Tout}(vecconv(Tout, h.a), Tout(h.β))
Base.convert{N,Tin,Tout}(::Type{HyperPlane{N,Tout}}, h::HyperPlane{N,Tin}) = HyperPlane{N,Tout}(vecconv(Tout, h.a), Tout(h.β))

islin(v::HalfSpace) = false
islin(v::HyperPlane) = true

function (*){ElemT<:HRepElement}(h::ElemT, P::Matrix)
  Tout = mypromote_type(eltype(ElemT), eltype(P))
  ElemTout = changeboth(ElemT, size(P, 2), Tout)
  ElemTout(h.a' * P, h.β)
end
function zeropad{ElemT<:HRepElement}(h::ElemT, n::Integer)
  if n == 0
    h
  else
    ElemTout = changefulldim(ElemT, fulldim(h) + abs(n))
    T = eltype(ElemT)
    if n < 0
      aout = [spzeros(T, n); h.a]
    else
      aout = [h.a; spzeros(T, n)]
    end
    ElemTout(aout, h.β)
  end
end
function Base.round{ElemT<:HRepElement}(h::ElemT, P::Matrix)
  ElemT(round(h.a), round(h.β))
end

# Point:
# Point{N, T} or AbstractVector{T}
# Linear Point:
# SymPoint{N,T}
# Ray:
# Vec{N, T} or Ray{N, T}
# Linear Ray:
# Line{N, T}

type SymPoint{N, T}
  a
end

type Ray{N, T}
  a
end

type Line{N,T}
  a
end

typealias FixedVRepElement{N,T} Union{Point{N,T},Vec{N,T},Ray{N,T},Line{N,T},SymPoint{N,T}}
typealias VRepElement{N,T} Union{FixedVRepElement{N,T},AbstractVector{T}}
typealias RepElement{N,T} Union{HRepElement{N,T}, VRepElement{N,T}}
typealias FixedRepElement{N,T} Union{HRepElement{N,T}, FixedVRepElement{N,T}}

fulldim{N}(e::FixedRepElement{N}) = N
eltype{N,T}(e::FixedRepElement{N, T}) = T
fulldim{T<:FixedRepElement}(::Type{T}) = T.parameters[1]
eltype{T<:FixedRepElement}(::Type{T}) = T.parameters[2]

fulldim(v::AbstractVector) = length(v)
fulldim{T<:AbstractVector}(::Type{T}) = T.parameters[1]

islin{T<:Union{Point,AbstractVector,Vec,Ray}}(v::T) = false
islin{T<:Union{SymPoint,Line}}(v::T) = true
isray{T<:Union{Point,AbstractVector,SymPoint}}(v::T) = false
isray{T<:Union{Vec,Ray,Line}}(v::T) = true

coord{ElemT<:Union{Point,AbstractVector,Vec}}(v::ElemT) = v
coord{ElemT<:Union{HRepElement,SymPoint,Ray,Line}}(v::ElemT) = v.a

typealias VRepElementContainer Union{SymPoint, Ray, Line}

@generated function (*){ElemT<:VRepElement}(P::Matrix, v::ElemT)
  if ElemT <: VRepElementContainer
    Tout = mypromote_type(eltype(ElemT), eltype(P))
    ElemTout = changeboth(ElemT, size(P, 1), Tout)
    return :(:ElemTout(P * v.a))
  else
    return :(P * v.a)
  end
end
function zeropad{ElemT<:VRepElement}(v::ElemT, n::Integer)
  if n == 0
    v
  else
    ElemTout = changefulldim(ElemT, fulldim(h) + abs(n))
    T = eltype(ElemT)
    if n < 0
      aout = [spzeros(T, n); coef(v)]
    else
      aout = [coef(v); spzeros(T, n)]
    end
    ElemTout(aout)
  end
end

for ElemT in [:HalfSpace, :HyperPlane, :Point, :SymPoint, :Vec, :Ray, :Line]
  @eval begin
    changeeltype{N,T,Tout}(::Type{$ElemT{N,T}}, ::Type{Tout}) = $ElemT{N,Tout}
    changefulldim{N,T}(::Type{$ElemT{N,T}}, Nout) = $ElemT{Nout,T}
    changeboth{N,T,Tout}(::Type{$ElemT{N,T}}, Nout, ::Type{Tout}) = $ElemT{Nout,Tout}
  end
end

changeeltype{VecT<:AbstractVector,Tout}(::Type{VecT}, ::Type{Tout}) = AbstractVector{Tout}
changefulldim{VecT<:AbstractVector}(::Type{VecT}, Nout) = VecT
changeboth{VecT<:AbstractVector,Tout}(::Type{VecT}, Nout, ::Type{Tout}) = AbstractVector{Tout}

Base.in{N}(h::HalfSpace{N}, r::Vec{N}) = mynonpos(mydot(h.a, r))
Base.in{N}(h::HalfSpace{N}, r::Ray{N}) = mynonpos(mydot(h.a, r.r))
Base.in{N}(h::HalfSpace{N}, l::Line{N}) = mynonpos(mydot(h.a, l.r))
Base.in{N}(h::HalfSpace{N}, p::Point{N}) = myleq(mydot(h.a, p), β)
Base.in{N}(h::HalfSpace{N}, p::AbstractVector) = myleq(mydot(h.a, p), β)
Base.in{N}(h::HalfSpace{N}, p::SymPoint{N}) = myleq(mydot(h.a, p.p), β)

Base.in{N}(h::HyperPlane{N}, r::Vec{N}) = myeqzero(mydot(h.a, r))
Base.in{N}(h::HyperPlane{N}, r::Ray{N}) = myeqzero(mydot(h.a, r.r))
Base.in{N}(h::HyperPlane{N}, l::Line{N}) = myeqzero(mydot(h.a, l.r))
Base.in{N}(h::HyperPlane{N}, p::Point{N}) = myeq(mydot(h.a, p), β)
Base.in{N}(h::HyperPlane{N}, p::AbstractVector) = myeq(mydot(h.a, p), β)
Base.in{N}(h::HyperPlane{N}, p::SymPoint{N}) = myeq(mydot(h.a, p.p), β)

import Base.vec
function Base.vec{N,T}(x::FixedVector{N,T})
  y = Vector{T}(N)
  for i in 1:N
    y[i] = x[i]
  end
  y
end
#Base.vec{ElemT::VRepElementContainer}(x::ElemT) = ElemT(vec(x.a))

pushbefore(a::AbstractVector, β) = [β; a]
function pushbefore{ElemT<:FixedVector}(a::ElemT, β, ElemTout = changefulldim(ElemT, fulldim(ElemT)+1))
  ElemTout([β; vec(a)])
end

function lift{ElemT <: HRepElement}(h::ElemT)
  ElemT(pushbefore(h.a, h.β), zero(eltype(ElemT)))
end
lift{N,T}(h::Vec{N,T}) = pushbefore(h, zero(T))
lift{N,T}(h::Ray{N,T}) = Ray{N+1,T}(pushbefore(h.a, zero(T)))
lift{N,T}(h::Line{N,T}) = Line{N+1,T}(pushbefore(h.a, zero(T)))
lift{N,T}(h::Point{N,T}) = pushbefore(h, one(T)), Vec{N+1,T}
lift{T}(h::AbstractVector{T}) = Ray{length(h)+1,T}(pushbefore(h.a, one(T)))
lift{N,T}(h::SymPoint{N,T}) = Line{N+1,T}(pushbefore(h.a, one(T)))