text:
import LinearAlgebra: ×
import Base: ^, getindex, ndims

export Circle, S¹, Torus
export ×, ndims, info

struct Circle{ℱ}
    a::ℱ
    b::ℱ
end

function Circle(b::FT) where {FT}
    return Circle(FT(0), b)
end

function Circle()
    return Circle(2π)
end

S¹ = Circle

function Base.show(io::IO, Ω::Circle)
    printstyled(io, "[", color = 226)
    a = @sprintf("%.3f", Ω.a)
    b = @sprintf("%.3f", Ω.b)
    printstyled("$a, $b", color = 7)
    printstyled(io, ")", color = 226)
end

struct Torus{DT}
    domains::DT
end

function Base.show(io::IO, Ω::Torus)
    for (i, domain) in enumerate(Ω.domains)
        print(domain)
        if i != length(Ω.domains)
            printstyled(io, "×", color = 118)
        end
    end
end

function ndims(Ω::Circle)
    return 1
end

function ndims(Ω::Torus)
    return length(Ω.domains)
end

getindex(t::Torus, i) = t.domains[i]

# Algebra
×(arg1::Circle, arg2::Circle) = Torus((arg1, arg2))
×(args::Torus, arg2::Circle) = Torus((args.domains..., arg2))
×(arg1::Circle, args::Torus) = Torus((arg1, args.domains...))
×(arg1::Torus, args::Torus) = Torus((arg1.domains..., args.domains...))
×(args::Torus) = Torus(args...)
# Exponentiation
^(circle::Circle, n::Int) = ^(circle, Val(n))
^(circle::Circle, ::Val{1}) = circle
^(circle::Circle, ::Val{2}) = circle × circle
^(circle::Circle, ::Val{n}) where {n} = circle × ^(circle, Val(n - 1))

function info(Ω::Torus)
    println("This is a ", ndims(Ω), "-dimensional Torus.")
    print("The domain is ")
    println(Ω, ".")
    for (i, domain) in enumerate(Ω.domains)
        domain_string = "periodic"
        length = @sprintf("%.2f ", domain.b - domain.a)
        println("The dimension $i domain is ", domain_string, " with length ≈ ", length)

    end
    return nothing
end

ext: jl

max_stars_repo_name: sandreza/FourierCore

max_issues_repo_name: sandreza/FourierCore

max_forks_repo_name: sandreza/FourierCore


text:
function simulate_snv(seq, chr, start, profile_snv)
    seqstr = seq.seq
    for snv in profile_snv
        if chr == snv["chrom"] && start <= snv["pos"] && start+length(seqstr) > snv["pos"]
            seqstr = mutate(seqstr, start, snv)
            # println(seqstr)
        end
    end
    return dna(seqstr)
end

function mutate(seqstr, start, snv)
    # simulate the mutation rate
    if rand() > snv["rate"]
        return seqstr
    end
    offset = snv["pos"] - start
    left = seqstr[1:offset]
    center = ASCIIString(snv["alt"])
    right = seqstr[offset+length(snv["ref"])+1:end]
    seqstr = left * center * right
    return seqstr
end


text:
"""
`struct PIEvaluator` : specifies a PIEvaluator, which is basically defined
through a PIBasis and its coefficients
"""
mutable struct PIEvaluator{T, TPI <: PIBasis}
   pibasis::TPI
   coeffs::Vector{T}
end

==(V1::PIEvaluator, V2::PIEvaluator) =
      (V1.pibasis == V2.pibasis) && (V1.coeffs == V2.coeffs)

# ------------------------------------------------------------
#   Initialisation and Parameter manipulation code
# ------------------------------------------------------------

function set_params!(ev::PIEvaluator, coeffs::AbstractVector)
   ev.coeffs[:] .= coeffs
   return ev
end

function set_params!(ev::PIEvaluator, basis::PIBasis, coeffs::AbstractVector)
   @assert ev.pibasis === basis
   set_params!(ev, coeffs)
end

_get_eff_coeffs!(ceff, basis::SymmetricBasis, c::AbstractVector) =
      genmul!(ceff, transpose(basis.A2Bmap), c, *)

function _get_eff_coeffs(basis::SymmetricBasis, c::AbstractVector{<: Number})
   ceff = zeros(eltype(basis.A2Bmap), size(basis.A2Bmap, 2))
   return _get_eff_coeffs!(ceff, basis, c)
end

set_params!(ev::PIEvaluator, basis::SymmetricBasis, coeffs::AbstractVector) =
      set_params!(ev, basis.pibasis, _get_eff_coeffs(basis, coeffs))

PIEvaluator(basis::SymmetricBasis, c::AbstractVector) =
      PIEvaluator(basis.pibasis, _get_eff_coeffs(basis, c))

# ------------------------------------------------------------
#   FIO code
# ------------------------------------------------------------

write_dict(V::PIEvaluator) = Dict(
      "__id__" => "ACE_PIEvaluator",
     "pibasis" => write_dict(V.pibasis),
      "coeffs" => write_dict(V.coeffs) )

read_dict(::Val{:ACE_PIEvaluator}, D::Dict) =
   PIEvaluator( read_dict(D["pibasis"]),
                read_dict(D["coeffs"]) )

# ------------------------------------------------------------
#   Standard Evaluation code
# ------------------------------------------------------------

# TODO: generalise the R, Z, allocation
alloc_temp(V::PIEvaluator{T}, args...) where {T} =
   (
      tmp_pibasis = alloc_temp(V.pibasis),
   )

evaluate!(tmp, ::LinearACEModel, V::PIEvaluator, cfg::AbstractConfiguration) =
      evaluate!(tmp, V::PIEvaluator, cfg)

# compute one site energy
function evaluate!(tmp, V::PIEvaluator, cfg::AbstractConfiguration)
   A = evaluate!(tmp.tmp_pibasis.A, tmp.tmp_pibasis.tmp1p,
                 V.pibasis.basis1p, cfg)
   spec = V.pibasis.spec
   _real = V.pibasis.real
   # initialize output with a sensible type
   val = zero(eltype(V.coeffs)) * _real(zero(eltype(A)))
   @inbounds for iAA = 1:length(spec)
      aa = A[spec.iAA2iA[iAA, 1]]
      for t = 2:spec.orders[iAA]
         aa *= A[spec.iAA2iA[iAA, t]]
      end
      val += _real(aa) * V.coeffs[iAA]
   end
   return val
end

alloc_temp_d(V::PIEvaluator{T}, N::Integer, args...) where {T} =
      (
       dAco = zeros( complex(eltype(V.coeffs)),
                    length(V.pibasis.basis1p) ),
       tmpd_pibasis = alloc_temp_d(V.pibasis, N),
       dAAdA = (@MVector zeros( fltype(V.pibasis.basis1p),
                               maximum(V.pibasis.spec.orders) )),
      )

grad_config!(g, tmpd, ::LinearACEModel, V::PIEvaluator, cfg::AbstractConfiguration) =
      evaluate_d!(g, tmpd, V::PIEvaluator, cfg)

# compute one site energy
function evaluate_d!(g, tmpd, V::PIEvaluator, cfg::AbstractConfiguration)
   basis1p = V.pibasis.basis1p
   tmpd_1p = tmpd.tmpd_pibasis.tmpd_basis1p
   A = tmpd.tmpd_pibasis.A
   dA = tmpd.tmpd_pibasis.dA
   _real = V.pibasis.real
   dAAdA = tmpd.dAAdA

   # stage 1: precompute all the A values
   evaluate_ed!(A, dA, tmpd_1p, basis1p, cfg)

   # stage 2: compute the coefficients for the ∇A_{klm} = ∇ϕ_{klm}
   dAco = tmpd.dAco
   c = V.coeffs
   spec = V.pibasis.spec
   fill!(dAco, zero(eltype(dAco)))
   @inbounds for iAA = 1:length(spec)
      _AA_local_adjoints!(dAAdA, A, spec.iAA2iA, iAA, spec.orders[iAA], _real)
      @fastmath for t = 1:spec.orders[iAA]
         dAco[spec.iAA2iA[iAA, t]] += dAAdA[t] * complex(c[iAA])
      end
   end

   # stage 3: get the gradients
   fill!(g, zero(eltype(g)))
   for iX = 1:length(cfg)
      @inbounds @fastmath for iA = 1:length(basis1p)
         g[iX] += _real.(dAco[iA] * dA[iA, iX])
      end
   end

   return g
end


text:
include("../../Computation/Computation.jl")
using Main.Computation

activation_func(x) = -x
sum_el(temp, idx, input_shape, activ) = activ(sum(view(temp, idx, 1:input_shape)))

function test_device_dense_logic(device)
    #Length of the array of data that the dense layer takes in
    input_shape = (5)
    #Output Size
    size = 10

    #Creates an array 1->size to be used in parallelized operations. Will be precalculated beforehand for actual layer
    indexes = alloc(device, collect(1:size))

    #Creates a matrix size:input_length that holds all the constants
    constants = alloc(device, fill(2, size, input_shape[1]))

    #Create an intermediate matrix that holds the output of the multiplication step
    temp = alloc(device, Int, size, input_shape[1])

    #The output array
    outs = alloc(device, Int, size)

    #The input array with dummy values
    ins = alloc(device, collect(1:size))

    #Select every row of constants and multiplty every element the corresponding element in inputs (IE: Weight the input).
    #Store the results in temp
    broadcast!(*, temp, constants, ins)

    #Take in the temp array, get the output row from the step above then sum the row and calculate activation value on the sum
    #Store the results in outs array
    broadcast!(sum_el, outs, Ref(temp), indexes, input_shape[1], activation_func)

    return Array(outs)
end

function test()
    println("Starting Device Dense Layer....")
    true_array = [-10, -20, -30, -40, -50, -60, -70, -80, -90, -100]
    for device in get_devices()
        println("\tChecking $device")
        @assert test_device_dense_logic(device) == true_array "Failed $device Dense Layer Test"
        println("\tSuccesfully Completed $device")
    end
    println("Successfully Compeleted Device Dense Layer Test!")
end

test()


text:
include("direct.jl")
include("rts.jl")
include("tv.jl")
include("nltv.jl")


text:
using XLNet
using Test

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


text:
using Pkg
Pkg.add("SpecialFunctions")
Pkg.add("BenchmarkTools")

using SpecialFunctions
using BenchmarkTools
#%%
x = range(0,100, length=10000)

#%%
using Base.Threads
Threads.nthreads()

results = zeros(length(x))

@btime results .= besselj1.(x)

@btime for i in 1:length(x)
    results[i] = besselj1(x[i])
end

@btime @threads for i in 1:length(x)
    results[i] = besselj1(x[i])
end

#%%
function slow_func(x)
    sleep(0.005)
    return x
end

y = 1:0.1:10
res = zeros(length(y))

@btime res .= slow_func.(y)

@btime for i in 1:length(y)
    res[i] = slow_func(y[i])
end

#%%
import Base.Threads.@spawn
using BenchmarkTools

function fib(n::Int)
    if n < 2
        return n
    end
    t = fib(n - 2)
    return fib(n - 1) + t
end

function fib_threads(n::Int)
    if n < 2
        return n
    end
    t = @spawn fib_threads(n - 2)
    return fib_threads(n - 1) + fetch(t)
end

fib(5)
fib_threads(5)
@btime fib(30)
@btime fib_threads(30)

#%%
import Base.Threads.@spawn
using BenchmarkTools

function slow_func(x)
    sleep(0.005)
    return x
end

@btime let
    x = 1:100
    a = @spawn slow_func(2)
    b = @spawn slow_func(4)
    c = @spawn slow_func(42)
    d = @spawn slow_func(12)
    res = fetch(a) .+ fetch(b) .* fetch(c) ./ fetch(d)
end

@btime let
    x = 1:100
    a = slow_func(2)
    b = slow_func(4)
    c = slow_func(42)
    d = slow_func(12)
    res = a .+ b .* c ./ d
end

#%%
@btime let
    x = 1:100
    a = @spawn sin(2)
    b = @spawn sin(4)
    c = @spawn sin(42)
    d = @spawn sin(12)
    res = fetch(a) .+ fetch(b) .* fetch(c) ./ fetch(d)
end

@btime let
    x = 1:100
    a = sin(2)
    b = sin(4)
    c = sin(42)
    d = sin(12)
    res = a .+ b .* c ./ d
end

# %%
using Pkg
Pkg.add("Distributed")

using Distributed
addprocs(4)

nprocs()

fetch(@spawn myid())
fetch(@spawnat 3 myid())

@everywhere using SpecialFunctions

fetch(@spawn gamma(5))

#%%
@everywhere function my_func(x)
    return x^3*cos(x)
end

fetch(@spawnat 2 my_func(0.42))

@everywhere using SharedArrays
res = SharedArray(zeros(10))
@distributed for x in 1:10
    res[x] = my_func(x)
end

res

#%%
using Distributed
using BenchmarkTools
addprocs(4)

@everywhere using SpecialFunctions
@everywhere function my_func(x)
    arg = repeat([x], 1000)
    return besselj1.(arg)
end
# %%
@btime my_func.(1:1000)

@btime my_func.(1:100)

@btime pmap(my_func, 1:1000, batch_size=1)
@btime pmap(my_func, 1:1000, batch_size=100)

max_stars_repo_name: JeffreySarnoff/techytok-examples

max_issues_repo_name: JeffreySarnoff/techytok-examples

max_forks_repo_name: JeffreySarnoff/techytok-examples


text:
# "converts nan values to 0.0"
# _convert_nan(x) = isnan(x) ? 0.0 : x
# _replace_nan(v) = map(x -> isnan(x) ? zero(x) : x, v)

# "Returns true if PowerModelsGraph `graph` has a `property` on an edge or a node `obj`"
# function hasprop(graph::PowerModelsGraph{T}, obj::Union{Int,LightGraphs.AbstractEdge}, property::Symbol) where T <: LightGraphs.AbstractGraph
#     if haskey(graph.metadata, obj)
#         return haskey(graph.metadata[obj], property)
#     else
#         return false
#     end
# end

# "Sets a `property` in the metadata at `key` of `graph` on `obj`"
# function set_property!(graph::PowerModelsGraph{T}, obj::Union{Int,LightGraphs.AbstractEdge}, key::Symbol, property::Any) where T <: LightGraphs.AbstractGraph
#     if !haskey(graph.metadata, obj)
#         graph.metadata[obj] = Dict{Symbol,Any}()
#     end

#     graph.metadata[obj][key] = property
# end

# "Sets multiple `properties` in the metadata of `graph` on `obj` at `key`"
# function set_properties!(graph::PowerModelsGraph{T}, obj::Union{Int,LightGraphs.AbstractEdge}, properties::Dict{Symbol,<:Any}) where T <: LightGraphs.AbstractGraph
#     if !haskey(graph.metadata, obj)
#         graph.metadata[obj] = Dict{Symbol,Any}()
#     end

#     merge!(graph.metadata[obj], properties)
# end

# "Gets the property in the metadata of `graph` on `obj` at `key`. If property doesn't exist, returns `default`"
# function get_property(graph::PowerModelsGraph{T}, obj::Union{Int,LightGraphs.AbstractEdge}, key::Symbol, default::Any) where T <: LightGraphs.AbstractGraph
#     return get(get(graph.metadata, obj, Dict{Symbol,Any}()), key, default)
# end

# "Adds an edge defined by `i` & `j` to `graph`"
# function add_edge!(graph::PowerModelsGraph{T}, i::Int, j::Int) where T <: LightGraphs.AbstractGraph
#     LightGraphs.add_edge!(graph.graph, i, j)
# end

# "Add vertex to the graph"
# function add_vertex!(graph::PowerModelsGraph{T}) where T <: LightGraphs.AbstractGraph
#     LightGraphs.add_vertex!(graph.graph)
# end

# function adjacency_matrix(graph::PowerModelsGraph{T}) where T <: LightGraphs.AbstractGraph
#     return LightGraphs.adjacency_matrix(graph.graph)
# end

# function  dijkstra_shortest_paths(graph::PowerModelsGraph{T}, i) where T <: LightGraphs.AbstractGraph
#     return LightGraphs.dijkstra_shortest_paths(graph.graph,i)
# end

# "Returns an iterator of all of the nodes/vertices in `graph`"
# function vertices(graph::PowerModelsGraph{T}) where T <: LightGraphs.AbstractGraph
#     return LightGraphs.vertices(graph.graph)
# end

# function nv(graph::PowerModelsGraph{T}) where T <: LightGraphs.AbstractGraph
#     return LightGraphs.nv(graph.graph)
# end

# "Returns an iterator of all the edges in `graph`"
# function edges(graph::PowerModelsGraph{T}) where T <: LightGraphs.AbstractGraph
#     return LightGraphs.edges(graph.graph)
# end

# "Returns all of the metadata for `obj` in `graph`"
# function properties(graph::PowerModelsGraph{T}, obj::Union{Int,LightGraphs.AbstractEdge}) where T <: LightGraphs.AbstractGraph
#     return get(graph.metadata, obj)
# end

# "get node/edge data from graph struct"
# function get_data(graph::PowerModelsGraph{T}, obj::Union{Int,LightGraphs.AbstractEdge}) where T <: LightGraphs.AbstractGraph
#     return get_property(graph, obj, :data, Dict{String,Any}())
# end

# "Return label information for `obj` in graph.annotationdata[\"label\"]"
# function get_label(graph::PowerModelsGraph{T}, obj::Union{Int,LightGraphs.AbstractEdge}, default::Any) where T <: LightGraphs.AbstractGraph
#     return get(get(graph.annotationdata, "label", Dict{Union{Int,LightGraphs.AbstractEdge},Any}()), obj, default)
# end

# "Update properties dictionary recursively."
# function update_properties!(default::Dict, new_data::Dict)
#     for (key,v_new) in new_data
#         if haskey(default, key)
#             v_default = default[key]
#             if isa(v_default, Dict) && isa(v_new, Dict)
#                 update_properties!(v_default, v_new)
#             else
#                 default[key] = v_new
#             end
#         else
#             default[key] = v_new
#         end
#     end
# end


text:
"""Julia program to implement Ternary Search algorithm.
Ternary Search is a Divide and Conquer strategy based algorithm similar to binary search.
Here we divide the entire array into three parts and recursively search for the required element
"""

function ternary_search_recursive(arr, low, high, ele)
    while (high >= low)
        # FInd the two required midpoints to divide the array into three.
        midA = low + (high - low) ÷ 3
        midB = high - (high - low) ÷ 3
        # Search, if the element is present in the midpoint indices
        if (arr[midA] == ele)
            return true
        elseif (arr[midB] == ele)
            return true
        end

        # Identify the part in which the required element belongs to and recursively call the function with limited length
        if (ele < arr[midA])
            return ternary_search_recursive(arr, low, midA - 1, ele)
        elseif (ele > arr[midB])
            return ternary_search_recursive(arr, midB + 1, high, ele)
        else
            return ternary_search_recursive(arr, midA + 1, midB - 1, ele)
        end
    end
    return false
end

print("How many numbers are present in the array? ")
n = readline()
n = parse(Int, n)
if (n <= 0)
    println("Array is Empty!!!")
    exit()
end
arr = Int[]
print("Enter the numbers: ")
arr = [parse(Int, num) for num in split(readline())]
print("Which number do you want to search in the array? ")
ele = readline()
ele = parse(Int, ele)
# Sort the array in ascending order
arr = sort(arr)
res = ternary_search_recursive(arr, 1, n, ele)
if (res == 0)
    print("The number $ele is not present in the array")
else
    print("The number $ele is present in the array.")
end

"""
Time Complexity - O(log(n)), where 'n' is the size of the array
Space Complexity - O(n)

SAMPLE INPUT AND OUTPUT

SAMPLE I

How many numbers are present in the array? 5
Enter the numbers: 1 2 3 4 5
Which number do you want to search in the array? 6
The number 6 is not present in the array

SAMPLE II

How many numbers are present in the array? 3
Enter the numbers: 3 1 2
Which number do you want to search in the array? 2
The number 2 is present in the array.

"""

max_stars_repo_name: TechSpiritSS/NeoAlgo


text:
# Warn only once, using the message `msg`.
# `funcsym` is a symbol that determines if the warning has been
# called before (so there is only one warning per `funcsym`).
function warn_once(msg, funcsym; force=true, stacktrace=true)
  if stacktrace
    io = IOBuffer()
    Base.show_backtrace(io, backtrace())
    backtrace_string = String(take!(io))
    backtrace_string *= "\n"
    msg *= backtrace_string
  end
  Base.depwarn(msg, funcsym; force)
  return nothing
end

# Directory helper functions (useful for
# running examples)
src_dir() = dirname(pathof(@__MODULE__))
pkg_dir() = joinpath(src_dir(), "..")
examples_dir() = joinpath(pkg_dir(), "examples")

# Determine version and uuid of the package
function _parse_project_toml(field::String)
  return Pkg.TOML.parsefile(joinpath(pkg_dir(), "Project.toml"))[field]
end
version() = VersionNumber(_parse_project_toml("version"))
uuid() = Base.UUID(_parse_project_toml("uuid"))


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "KAShim"
version = v"0.1.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/vchuravy/KAShim.git", "53a435a7bc2ab1eb1a92c784d79acf2546b4727e")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/KAShim/
make CFLAGS="-I${prefix}/include -fPIC" LDFLAGS="-L${prefix}/lib -luv" install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libkashim", :libkashim)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="LibUV_jll", uuid="183b4373-6708-53ba-ad28-60e28bb38547"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
######################################################################
# max_trajectory.jl
# tries to maximize phi(x) collected
# TODO: this is a wreck. clean it so it looks like clerc_trajectory.
######################################################################

function max_trajectory(em::ErgodicManager, tm::TrajectoryManager, mu::VF, Sigma::MF)
	return max_trajectory(em, tm, [mu], [Sigma])
end

function max_trajectory(em::ErgodicManager, tm::TrajectoryManager, mus::VVF, Sigmas::VMF, weights::VF=ones(length(mus)); max_iters=100)

	# initialize trajectory
	xd, ud = initialize(tm.initializer, em, tm)

	# TODO: really a termination condition
	for i = 1:max_iters

		ad, bd = compute_gradients(tm, xd, ud, mus, Sigmas, weights)
		A, B = linearize(tm.dynamics, xd, ud, tm.h)
		K, C = LQ(A, B, ad, bd, tm.Qn, tm.Rn, tm.N)

		#zd, vd = convex_descent(A, B, ad, bd, N)
		zd, vd = apply_LQ_gains(A, B, K, C)
		step_size = .15 / sqrt(i)
		#step_size = .25 / sqrt(i)

		# printing statistics for testing
		println("i = ",i)
		dd = directional_derivative(ad, bd, zd, vd)
		sdd = scaled_dd(ad, bd, zd, vd)
		#es = ergodic_score(em, xd)
		#ts = total_score(em, xd, ud, T)
		#cs = control_score(ud, 0.01*eye(2), T)
		#println("es = ", es)
		#println("ts = ", ts)
		#println("cs = ", cs)
		println("dd = ", dd)
		println("scaled_dd = ", sdd)
		#println("alpha = ", step_size)
		println("##################################")

		descend!(xd, ud, zd, vd, step_size)
	end

	return xd, ud
end

export max_trajectory

# TODO: should really take in a bunch of gaussians here
# returns a and b, each of which is an array of vectors
function compute_gradients(tm::TrajectoryManager, xd::VVF, ud::VVF, mus::VVF, Sigmas::VMF, weights::VF)
	#a = Array(VF, tm.N+1)
	#b = Array(VF, tm.N)
	a = zeros(tm.dynamics.n, tm.N+1)
	b = zeros(tm.dynamics.m, tm.N)

	for ni = 1:tm.N
		a[1:2, ni] = compute_an(mus, Sigmas, weights, xd[ni], tm.h, tm.q)
		b[:, ni] = tm.h * tm.R * ud[ni]
	end
	a[1:2,tm.N+1] = compute_an(mus, Sigmas, weights, xd[tm.N+1], tm.h, tm.q)

	return a, b
end

function compute_an(mus::VVF, Sigmas::VMF, weights::VF, x::VF, h::Float64, q::Float64)
	num_gauss = length(weights)
	an = zeros(2)
	for i = 1:num_gauss
		mu = mus[i]
		Sigma = Sigmas[i]
		xmu = x[1:2] - mu
		inv_Sigma = inv(Sigma)
		c = q*h*exp(-.5dot(xmu, inv_Sigma*xmu)) / (2*pi*sqrt(det(Sigma)))
		an += weights[i]*c*inv_Sigma*xmu
	end
	return an
end


text:
using JuMP
using Ipopt
using Hiop
using Printf
using DelimitedFiles
using acopf

function solve(opfmodel)
  optimize!(opfmodel)
  status = termination_status(opfmodel)
  if status != MOI.LOCALLY_SOLVED
      println("Could not solve the model to optimality.")
  end
  return opfmodel, status
end

function model(opf_data; max_iter=100, solver="Ipopt")
  #shortcuts for compactness
  Pg0,Qg0,Vm0,Va0 = acopf.initialPt_IPOPT(opf_data)
  lines = opf_data.lines; buses = opf_data.buses; generators = opf_data.generators; baseMVA = opf_data.baseMVA
  busIdx = opf_data.BusIdx; FromLines = opf_data.FromLines; ToLines = opf_data.ToLines; BusGeners = opf_data.BusGenerators;

  nbus  = length(buses); nline = length(lines); ngen  = length(generators)

  #branch admitances
  YffR,YffI,YttR,YttI,YftR,YftI,YtfR,YtfI,YshR,YshI = acopf.computeAdmitances(lines, buses, baseMVA)

  #
  # JuMP model now
  #
  if solver == "Hiop"
    opfmodel = Model(optimizer_with_attributes(Hiop.Optimizer))
  else
    opfmodel = Model(optimizer_with_attributes(Ipopt.Optimizer, "max_iter" => max_iter))
  end

  @variable(opfmodel, generators[i].Pmin <= Pg[i=1:ngen] <= generators[i].Pmax, start = Pg0[i])
  @variable(opfmodel, generators[i].Qmin <= Qg[i=1:ngen] <= generators[i].Qmax, start = Qg0[i])

  @variable(opfmodel, Va[i=1:nbus], start = Va0[i])
  @variable(opfmodel, buses[i].Vmin <= Vm[i=1:nbus] <= buses[i].Vmax, start = Vm0[i])
  # fix the voltage angle at the reference bus
  set_lower_bound(Va[opf_data.bus_ref], buses[opf_data.bus_ref].Va)
  set_upper_bound(Va[opf_data.bus_ref], buses[opf_data.bus_ref].Va)

  # minimize active power
  coeff0 = Vector{Float64}(undef, ngen)
  for (i,v) in enumerate(generators)
      coeff0[i] = v.coeff[v.n]
  end
  coeff1 = Vector{Float64}(undef, ngen)
  for (i,v) in enumerate(generators)
      coeff1[i] = v.coeff[v.n - 1]
  end
  coeff2 = Vector{Float64}(undef, ngen)
  for (i,v) in enumerate(generators)
      coeff2[i] = v.coeff[v.n - 2]
  end

  @NLobjective(opfmodel, Min, sum( coeff2[i]*(baseMVA*Pg[i])^2
                                  + coeff1[i]*(baseMVA*Pg[i])
                                  + coeff0[i] for i=1:ngen))

  #
  # power flow balance
  #

  for b in 1:nbus
      # real part
      @NLconstraint(
      opfmodel,
      ( sum( YffR[l] for l in FromLines[b]) + sum( YttR[l] for l in ToLines[b]) + YshR[b] ) * Vm[b]^2
      + sum( Vm[b]*Vm[busIdx[lines[l].to]]  *( YftR[l]*cos(Va[b]-Va[busIdx[lines[l].to]]  ) + YftI[l]*sin(Va[b]-Va[busIdx[lines[l].to]]  )) for l in FromLines[b] )
      + sum( Vm[b]*Vm[busIdx[lines[l].from]]*( YtfR[l]*cos(Va[b]-Va[busIdx[lines[l].from]]) + YtfI[l]*sin(Va[b]-Va[busIdx[lines[l].from]])) for l in ToLines[b]   )
      - ( sum(baseMVA*Pg[g] for g in BusGeners[b]) - buses[b].Pd ) / baseMVA      # Sbus part
      ==0)
      # imaginary part
  end
  for b in 1:nbus
      @NLconstraint(
      opfmodel,
      ( sum(-YffI[l] for l in FromLines[b]) + sum(-YttI[l] for l in ToLines[b]) - YshI[b] ) * Vm[b]^2
      + sum( Vm[b]*Vm[busIdx[lines[l].to]]  *(-YftI[l]*cos(Va[b]-Va[busIdx[lines[l].to]]  ) + YftR[l]*sin(Va[b]-Va[busIdx[lines[l].to]]  )) for l in FromLines[b] )
      + sum( Vm[b]*Vm[busIdx[lines[l].from]]*(-YtfI[l]*cos(Va[b]-Va[busIdx[lines[l].from]]) + YtfR[l]*sin(Va[b]-Va[busIdx[lines[l].from]])) for l in ToLines[b]   )
      - ( sum(baseMVA*Qg[g] for g in BusGeners[b]) - buses[b].Qd ) / baseMVA      #Sbus part
      ==0)
  end
  #
  # branch/lines flow limits
  #
  nlinelim=0
  for l in 1:nline
      if lines[l].rateA!=0 && lines[l].rateA<1.0e10
      nlinelim += 1
      flowmax=(lines[l].rateA/baseMVA)^2

      #branch apparent power limits (from bus)
      Yff_abs2=YffR[l]^2+YffI[l]^2;        Yft_abs2=YftR[l]^2+YftI[l]^2
      Yre=YffR[l]*YftR[l]+YffI[l]*YftI[l]; Yim=-YffR[l]*YftI[l]+YffI[l]*YftR[l]
      @NLconstraint(
          opfmodel,
          Vm[busIdx[lines[l].from]]^2 *
          ( Yff_abs2*Vm[busIdx[lines[l].from]]^2 + Yft_abs2*Vm[busIdx[lines[l].to]]^2
          + 2*Vm[busIdx[lines[l].from]]*Vm[busIdx[lines[l].to]]
              *(Yre*cos(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]])
              -Yim*sin(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]])
              )
          )
          - flowmax <=0)

      #branch apparent power limits (to bus)
      Ytf_abs2=YtfR[l]^2+YtfI[l]^2; Ytt_abs2=YttR[l]^2+YttI[l]^2
      Yre=YtfR[l]*YttR[l]+YtfI[l]*YttI[l]; Yim=-YtfR[l]*YttI[l]+YtfI[l]*YttR[l]
      @NLconstraint(
          opfmodel,
          Vm[busIdx[lines[l].to]]^2 *
          ( Ytf_abs2*Vm[busIdx[lines[l].from]]^2 + Ytt_abs2*Vm[busIdx[lines[l].to]]^2
          + 2*Vm[busIdx[lines[l].from]]*Vm[busIdx[lines[l].to]]
              *(Yre*cos(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]])
              -Yim*sin(Va[busIdx[lines[l].from]]-Va[busIdx[lines[l].to]])
              )
          )
          - flowmax <=0)
      end
  end

  @printf("Buses: %d  Lines: %d  Generators: %d\n", nbus, nline, ngen)
  println("Lines with limits  ", nlinelim)

  return opfmodel, Pg, Qg, Va, Vm
end
# Compute initial point for IPOPT based on the values provided in the case data
function initialPt_IPOPT(opfdata)
Pg=zeros(length(opfdata.generators)); Qg=zeros(length(opfdata.generators)); i=1
for g in opfdata.generators
    # set the power levels in in between the bounds as suggested by matpower
    # (case data also contains initial values in .Pg and .Qg - not used with IPOPT)
    Pg[i]=0.5*(g.Pmax+g.Pmin)
    Qg[i]=0.5*(g.Qmax+g.Qmin)
    i=i+1
end
@assert i-1==length(opfdata.generators)

Vm=zeros(length(opfdata.buses)); i=1;
for b in opfdata.buses
    # set the ini val for voltage magnitude in between the bounds
    # (case data contains initials values in Vm - not used with IPOPT)
    Vm[i]=0.5*(b.Vmax+b.Vmin);
    i=i+1
end
@assert i-1==length(opfdata.buses)

# set all angles to the angle of the reference bus
Va = opfdata.buses[opfdata.bus_ref].Va * ones(length(opfdata.buses))

return Pg,Qg,Vm,Va
end

function outputAll(opfmodel, opf_data, Pg, Qg, Va, Vm)
  #shortcuts for compactness
  lines = opf_data.lines; buses = opf_data.buses; generators = opf_data.generators; baseMVA = opf_data.baseMVA
  busIdx = opf_data.BusIdx; FromLines = opf_data.FromLines; ToLines = opf_data.ToLines; BusGeners = opf_data.BusGenerators;

  nbus  = length(buses); nline = length(lines); ngen  = length(generators)

  # OUTPUTING
  println("Objective value: ", objective_value(opfmodel), "USD/hr")
  VM = Array{Float64}(undef, size(Vm, 1))
  VM=value.(Vm); VA=value.(Va)
  PG=value.(Pg); QG=value.(Qg)

  println("============================= BUSES ==================================")
  println("  BUS    Vm     Va   |   Pg (MW)    Qg(MVAr) ")   # |    P (MW)     Q (MVAr)")  #|         (load)   ")

  println("                     |     (generation)      ")
  println("----------------------------------------------------------------------")
  for i in 1:nbus
    @printf("%4d | %6.2f  %6.2f | %s  | \n",
	    buses[i].bus_i, VM[i], VA[i]*180/pi,
	    length(BusGeners[i])==0 ? "   --          --  " : @sprintf("%7.2f     %7.2f", baseMVA*PG[BusGeners[i][1]], baseMVA*QG[BusGeners[i][1]]))
  end
  println("\n")

  within=20 # percentage close to the limits

  nflowlim=0
  for l in 1:nline
    if lines[l].rateA!=0 && lines[l].rateA<1.0e10
      nflowlim += 1
    end
  end

  if nflowlim>0
    println("Number of lines with flow limits: ", nflowlim)

    optvec=zeros(2*nbus+2*ngen)
    optvec[1:ngen]=PG
    optvec[ngen+1:2*ngen]=QG
    optvec[2*ngen+1:2*ngen+nbus]=VM
    optvec[2*ngen+nbus+1:2*ngen+2*nbus]=VA

    d = JuMP.NLPEvaluator(opfmodel)
    MOI.initialize(d, [:Jac])

    consRhs = zeros(2*nbus+2*nflowlim)
    MOI.eval_constraint(d, consRhs, optvec)

    #println(consRhs)

    @printf("================ Lines within %d %s of flow capacity ===================\n", within, "%")
    println("Line   From Bus    To Bus    At capacity")

    nlim=1
    for l in 1:nline
      if lines[l].rateA!=0 && lines[l].rateA<1.0e10
        flowmax=(lines[l].rateA/baseMVA)^2
        idx = 2*nbus+nlim

        if( (consRhs[idx]+flowmax)  >= (1-within/100)^2*flowmax )
          @printf("%3d      %3d      %3d        %5.3f%s\n", l, lines[l].from, lines[l].to, 100*sqrt((consRhs[idx]+flowmax)/flowmax), "%" )
          # @printf("%7.4f   %7.4f    %7.4f \n", consRhs[idx], consRhs[idx]+flowmax,  flowmax)
        end
        nlim += 1
      end
    end
  end

  #println(getvalue(Vm))
  #println(getvalue(Va)*180/pi)

  #println(getvalue(Pg))
  #println(getvalue(Qg))

  return
end

function main()
  if length(ARGS) < 1
    println("Usage: julia scopf_main.jl case_name")
    println("Cases are in 'data' directory: case9 case30 case57 case118 case300 case1354pegase case2383wp case2736sp case2737sop case2746wop case2869pegase case3012wp  case3120sp case3375wp case9241pegase")
    return
  end
  max_iter = 100
  opfdata = acopf.opf_loaddata(ARGS[1])
  opfmodel, Pg, Qg, Va, Vm = model(opfdata; max_iter = max_iter)
  opfmodel,status = solve(opfmodel)
  if status == MOI.LOCALLY_SOLVED
    outputAll(opfmodel,opfdata, Pg, Qg, Va, Vm)
  end
  objective_value(opfmodel) ≈ 5.2966862e+03
end

main()


text:
export format!, format_conditional!, format_color_gradient!

"""
    format!(client::GoogleSheetsClient, range::CellIndexRange2D, format::CellFormat)::Dict{Any,Any}

Formats a range of cells.

# Arguments
- `client::GoogleSheetsClient`: client.
- `range::CellIndexRange2D`: cell index range.
- `format::CellFormat`: cell format.
"""
function format!(client::GoogleSheetsClient, range::CellIndexRange2D, format::CellFormat)::Dict{Any,Any}
    body = Dict(
        "requests" => [
            Dict(
                "repeatCell" => Dict(
                    "range" => gsheet_json(range),
                    "cell" => Dict(
                        "userEnteredFormat" => gsheet_json(format),
                    ),
                    # "fields" => "userEnteredFormat.numberFormat",
                    "fields" => "userEnteredFormat",
                ),
            ),
        ],
    )

    return batch_update!(client, range.sheet.spreadsheet, body)
end

"""
    format_conditional!(client::GoogleSheetsClient, range::CellIndexRange2D, format::CellFormat,
        condition_type::ConditionType, values...)::Dict{Any,Any}

Formats a range of cells if a condition is met.

# Arguments
- `client::GoogleSheetsClient`: client.
- `range::CellIndexRange2D`: cell index range.
- `format::CellFormat`: cell format.
- `condition_type::ConditionType`: type of condition.
- `values...`: values for the condition.
"""
function format_conditional!(client::GoogleSheetsClient, range::CellIndexRange2D, format::CellFormat,
    condition_type::ConditionType, values...)::Dict{Any,Any}

    body = Dict(
        "requests" => [
            Dict(
                "addConditionalFormatRule" => Dict(
                    "rule" => Dict(
                        "ranges" => [gsheet_json(range)],
                        "booleanRule" => Dict(
                            "condition" => Dict(
                                "type" => gsheet_string(condition_type),
                                "values" => [ Dict("userEnteredValue" => "$value") for value in values ],
                            ),
                            "format" => gsheet_json(format),
                        ),
                    ),
                ),
            ),
        ],
    )

    return batch_update!(client, range.sheet.spreadsheet, body)
end

"""
    format_color_gradient!(client::GoogleSheetsClient, range::CellIndexRange2D;
        min_color::Colorant=colorant"salmon", min_value_type::ValueType=VALUE_TYPE_MIN, min_value::Union{Nothing,Number}=nothing,
        mid_color::Union{Nothing,Colorant}=nothing, mid_value_type::Union{Nothing,ValueType}=nothing, mid_value::Union{Nothing,Number}=nothing,
        max_color::Colorant=colorant"springgreen", max_value_type::ValueType=VALUE_TYPE_MAX, max_value::Union{Nothing,Number}=nothing)::Dict{Any,Any}

Sets color gradient formatting.

# Arguments
- `client::GoogleSheetsClient`: client.
- `range::CellIndexRange2D`: cell index range.

- `min_color::Colorant=colorant"salmon"`: color for the minimum value.
- `min_value_type::ValueType=VALUE_TYPE_MIN`: minimum value type.
- `min_value::Union{Nothing,Number}=nothing`: minimum value.
- `mid_color::Union{Nothing,Colorant}=nothing`: color for the mid value.
- `mid_value_type::Union{Nothing,ValueType}=nothing`: mid value type.
- `mid_value::Union{Nothing,Number}=nothing`: mid value.
- `max_color::Colorant=colorant"springgreen"`: color for the maximum value.
- `max_value_type::ValueType=VALUE_TYPE_MAX`: maximum value type.
- `max_value::Union{Nothing,Number}=nothing`: maximum value.
"""
function format_color_gradient!(client::GoogleSheetsClient, range::CellIndexRange2D;
    min_color::Colorant=colorant"salmon", min_value_type::ValueType=VALUE_TYPE_MIN, min_value::Union{Nothing,Number}=nothing,
    mid_color::Union{Nothing,Colorant}=nothing, mid_value_type::Union{Nothing,ValueType}=nothing, mid_value::Union{Nothing,Number}=nothing,
    max_color::Colorant=colorant"springgreen", max_value_type::ValueType=VALUE_TYPE_MAX, max_value::Union{Nothing,Number}=nothing)::Dict{Any,Any}

    function point(color, type, value)
        rst = Dict{Any,Any}(
                "color" => gsheet_json(color),
        )

        if !isnothing(type)
            rst["type"] = type
        end

        if !isnothing(value)
            rst["value"] = "$value"
        end

        return rst
    end

    value_type(x) =return isnothing(x) ? nothing : gsheet_string(x)

    function gradientRule()
        value_type_min = value_type(min_value_type)
        value_type_max = value_type(max_value_type)
        value_type_mid = value_type(mid_value_type)

        rst = Dict(
                "minpoint" => point(min_color, value_type_min, min_value),
                "maxpoint" => point(max_color, value_type_max, max_value),
        )

        if !isnothing(mid_color)
            rst["midpoint"] = point(mid_color, value_type_mid, mid_value)
        end

        return rst
    end

    body = Dict(
        "requests" => [
            Dict(
                "addConditionalFormatRule" => Dict(
                    "rule" => Dict(
                        "ranges" => [gsheet_json(range)],
                        "gradientRule" => gradientRule(),
                    ),
                ),
            ),
        ],
    )

    return batch_update!(client, range.sheet.spreadsheet, body)
end


text:
"""
    AbstractMutation{T}

Abstract type defining a mutation operation on entities of type `T`.

Implementations are expected to be callable using an entity of type `T` as only input.

May also implement a callable accepting a `AbstractVector{<:T}` if it is useful to work on
all items to mutate at once.
"""
abstract type AbstractMutation{T} end

(m::AbstractMutation{T})(es::AbstractVector{<:T}) where T = m.(es)

"""
    AbstractCrossover{T}

Type alias for `AbstractMutation{Tuple{T,T}}` defining a crossover of two entities of type `T`.

Implementations are expected to be callable using a tuple of two type `T` as only input.
"""
const AbstractCrossover{T} = AbstractMutation{Tuple{T,T}}

"""
    DecoratingMutation{T}

Abstract type indicating that the type itself does not perform any mutation but wraps a type which might do.

Must either implement callable method for `AbstractVector{<:T}` or accept keyword arguments `next=wrapped(m)` and
`noop=identity` along with a single `T`.
"""
abstract type DecoratingMutation{T} <: AbstractMutation{T} end
wrapped(m::DecoratingMutation) = m.m

mutationleaves(m::DecoratingMutation) = (mutationleaves(wrapped(m))...,)
mutationleaves(tm::Tuple) = mapreduce(mutationleaves, (t1,t2) -> (t1...,t2...), tm)
mutationleaves(m) = tuple(m)

# Apart from being overengineered this helps protecting against forgetting to handle arrays in a DecoratingMutation
# The next/noop happens to work with most existing DecoratingMutations, but it is a bit arbitrary and in some cases
# one must implement both the per-element and the vector of elements versions.
function (m::DecoratingMutation{T})(es::AbstractVector{<:T}) where T
    cnt = Ref(1)
    fornext = Int[]
    next = function(e)
        push!(fornext, cnt[])
        cnt[] += 1
        e
    end
    noop = function(e)
        cnt[] += 1
        e
    end

    allres = m.(es; next, noop)
    mres = wrapped(m)(es[fornext])

    # Mutation might accidentally widen the type compared to allres and then we can't insert mres into allres.
    # Lets fix that if it happens
    RT = typejoin(eltype(allres), eltype(mres))
    res = RT === eltype(allres) ? allres :  convert(Vector{RT}, allres)

    res[fornext] = mres
    return res
end

"""
    MutationProbability{T} <: DecoratingMutation{T}
    MutationProbability(m::AbstractMutation{T}, p::Probability)
    MutationProbability(m::AbstractMutation{T}, p::Number)

Applies `m` with probability `p`.
"""
struct MutationProbability{T, P<:Probability} <: DecoratingMutation{T}
    m::AbstractMutation{T}
    p::P
end
MutationProbability(m::AbstractMutation{T}, p::Number) where T = MutationProbability(m, Probability(p))
(m::MutationProbability{T})(e::T; next=m.m, noop=identity) where T = apply(() -> next(e), m.p, () -> noop(e))

"""
    WeightedMutationProbability{T,F} <: DecoratingMutation{T}
    WeightedMutationProbability(m::AbstractMutation::T, pfun::F)

Applies `m` to an entity `e` with a probability `pfun(e)`.
"""
struct WeightedMutationProbability{T,F} <: DecoratingMutation{T}
    m::AbstractMutation{T}
    pfun::F
end
(m::WeightedMutationProbability{T})(e::T; next=m.m, noop=identity) where T = apply(() -> next(e), m.pfun(e), () -> noop(e))

"""
    HighUtilityMutationProbability(m::AbstractMutation{T}, pbase::Real, rng=rng_default; spread=0.5)

Return a `WeightedMutationProbability` which applies `m` to vertices with an (approximately) average probability of `pbase` and where high `neuronutility` compared to other vertices in same graph means higher probability.

Parameter `spread` can be used to control how much the difference in probability is between high and low utlity. High spread means high difference while low spread means low difference.
"""
HighUtilityMutationProbability(m::AbstractMutation{T}, pbase::Real, rng=rng_default;spread=0.5) where T <: AbstractVertex = WeightedMutationProbability(m, weighted_neuronutility_high(pbase, rng,spread=spread))

"""
    LowUtilityMutationProbability(m::AbstractMutation{T}, pbase::Real, rng=rng_default; spread=2)

Return a `WeightedMutationProbability` which applies `m` to vertices with an (approximately) average probability of `pbase` and where low `neuronutility` compared to other vertices in same graph means higher probability.

Parameter `spread` can be used to control how much the difference in probability is between high and low utlity. High spread means high difference while low spread means low difference.
"""
LowUtilityMutationProbability(m::AbstractMutation{T}, pbase::Real, rng=rng_default;spread=2) where T <: AbstractVertex = WeightedMutationProbability(m, weighted_neuronutility_low(pbase, rng, spread=spread))

weighted_neuronutility_high(pbase, rng=rng_default; spread=0.5) = function(v::AbstractVertex)
    ismissing(NaiveNASflux.neuronutility(v)) && return pbase
    return Probability(fixnan(pbase ^ normexp(v, spread), pbase), rng)
end

weighted_neuronutility_low(pbase, rng=rng_default;spread=2) = function(v::AbstractVertex)
    ismissing(NaiveNASflux.neuronutility(v)) && return pbase
    return Probability(fixnan(pbase ^ (1/normexp(v, 1/spread)), pbase), rng)
end

fixnan(x, rep) = isnan(x) ? rep : clamp(x, 0.0, 1.0)

# This is pretty hacky and arbitrary. Change to something better
function normexp(v::AbstractVertex, s)
    allvertices = filter(allow_mutation, all_in_graph(v))
    allvalues = map(vi -> NaiveNASflux.neuronutility(vi), allvertices)
    meanvalues = map(mean, skipmissing(allvalues))
    meanvalue = mean(meanvalues)
    maxvalue = maximum(meanvalues)
    utlity = mean(NaiveNASflux.neuronutility(v))
    # Basic idea: maxvalue - utlity means the (to be) exponent is <= 0 while the division seems to normalize so that
    # average of pbase ^ normexp across allvertices is near pbase (no proof!). The factor 2 is just to prevent
    # probability of vertex with maxvalue to be 1.
    return (2maxvalue^s - utlity^s) / (2maxvalue^s - meanvalue^s)
end

"""
    MutationChain{T} <: DecoratingMutation{T}
    MutationChain(m::AbstractMutation{T}...)

Chains multiple `AbstractMutation{T}`s after each other.

Input entities will be mutated by the first `AbstractMutation{T}` in the chain and the output will be fed into the next `AbstractMutation{T}` in the chain and so on. The output from the last `AbstractMutation{T}` is returned.
"""
struct MutationChain{T} <: DecoratingMutation{T}
    m::Tuple{Vararg{AbstractMutation{T}}}
end
MutationChain(m::AbstractMutation{T}...) where T = MutationChain(m)
# Identical, but can't use Union due to ambiguity
(m::MutationChain{T})(es::AbstractVector{<:T}) where T = foldl((ei, mi) -> mi(ei), m.m; init=es)
(m::MutationChain{T})(e::T) where T = foldl((ei, mi) -> mi(ei), m.m; init=e)

"""
    RecordMutation{T} <: DecoratingMutation{T}
    RecordMutation(m::AbstractMutation{T})

Records all mutated entities.

Intended use case is to be able to do parameter selection on mutated vertices.
"""
struct RecordMutation{T} <: DecoratingMutation{T}
    m::AbstractMutation{T}
    mutated::Vector{T}
end
RecordMutation(m::AbstractMutation{T}) where T = RecordMutation(m, T[])
function (m::RecordMutation{T})(e::T; next=m.m, noop=identity) where T
    em = next(e)
    push!(m.mutated, em)
    return em
end
function fetchmutated!(m::RecordMutation)
    mutated = copy(m.mutated)
    deleteat!(m.mutated, eachindex(m.mutated))
    return mutated
end

"""
    LogMutation{T} < :DecoratingMutation{T}
    LogMutation(strfun, m::AbstractMutation{T}; level = Logging.Info, nextlogfun=e -> PrefixLogger("   "))
    LogMutation(strfun, level::LogLevel, nextlogfun, m::AbstractMutation{T})

Logs all mutation operations.

Argument `strfun` maps the mutated entity to the logged string.

Calling `nextlogfun(e)` where `e` is the entity to mutate produces an `AbstractLogger` which will be used when applying `m(e)`.

By default, this is used to add a level of indentation to subsequent logging calls which makes logs of hierarchical mutations (e.g. mutate a CompGraph by applying mutations to some of its vertices) easier to read. Set `nextlogfun = e -> current_logger()` to remove this behaviour.
"""
struct LogMutation{T,F,L<:LogLevel,LF} <: DecoratingMutation{T}
    strfun::F
    level::L
    nextlogfun::LF
    m::AbstractMutation{T}
end
LogMutation(strfun, m::AbstractMutation{T}; level = Logging.Info, nextlogfun=e -> PrefixLogger("   ")) where T = LogMutation(strfun, level, nextlogfun, m)
function (m::LogMutation{T})(e::T; next=m.m, noop=identity) where T
    @logmsg m.level m.strfun(e)
    return with_logger(() -> next(e), m.nextlogfun(e))
end

"""
    MutationFilter{T} <: DecoratingMutation{T}
    MutationFilter(predicate, m)

Applies mutation `m` only for entities `e` for which `predicate(e)` returns true.
"""
struct MutationFilter{T,P} <: DecoratingMutation{T}
    predicate::P
    m::AbstractMutation{T}
end
function (m::MutationFilter{T})(e::T; next=m.m, noop=identity) where T
    m.predicate(e) && return next(e)
    return noop(e)
end

"""
    VertexMutation <: DecoratingMutation{CompGraph}
    VertexMutation(m::AbstractMutation{AbstractVertex}, s::AbstractVertexSelection)
    VertexMutation(m::AbstractMutation{AbstractVertex})

Applies a wrapped `AbstractMutation{AbstractVertex}` to each selected vertex in a `CompGraph`.

Vertices to select is determined by the configured `AbstractVertexSelection`.
"""
struct VertexMutation{S<:AbstractVertexSelection} <: DecoratingMutation{CompGraph}
    m::AbstractMutation{AbstractVertex}
    s::S
end
VertexMutation(m::AbstractMutation{AbstractVertex}) = VertexMutation(m, FilterMutationAllowed())
function (m::VertexMutation)(g::CompGraph)
    m.m(select(m.s, g, m))
    return g
end

"""
    NoutMutation <:AbstractMutation{AbstractVertex}
    NoutMutation(l1::Real,l2::Real, rng::AbstractRNG)
    NoutMutation(limit, rng::AbstractRNG=rng_default)
    NoutMutation(l1,l2)

Mutate the out size of a vertex or vector of vertices.

Size is changed by `x * nout(v)` rounded away from from zero where `x` is drawn from `U(minrel, maxrel)` where
`minrel` and `maxrel` are `l1` and `l2` if `l1 < l2` and `l2` and `l1` otherwise.
"""
struct NoutMutation{R<:Real, RNG<:AbstractRNG} <:AbstractMutation{AbstractVertex}
    minrel::R
    maxrel::R
    rng::RNG
    function NoutMutation(l1::R1, l2::R2, rng::RNG) where {R1, R2, RNG}
        R = promote_type(R1, R2)
        return l1 < l2 ? new{R, RNG}(promote(l1, l2)..., rng) : new{R, RNG}(promote(l2, l1)..., rng)
    end
end
NoutMutation(limit, rng::AbstractRNG=rng_default) = NoutMutation(0, limit, rng)
NoutMutation(l1,l2) = NoutMutation(l1,l2, rng_default)
(m::NoutMutation)(v::AbstractVertex) = first(m([v]))
function (m::NoutMutation)(vs::AbstractVector{<:AbstractVertex})

    Δs = Dict{AbstractVertex, Int}()
    shift = m.minrel
    scale = m.maxrel - m.minrel

    for v in vs
        terminputs = findterminating(v, inputs)

        # We are basically just searching for Immutable vertices here, allow_mutation(trait(v)) happens to do just that
        any(tv -> allow_mutation(trait(tv)), terminputs) || continue

        Δfloat = rand(m.rng) * scale + shift

        Δ = ceil(Int, abs(Δfloat) * nout(v)) *  sign(Δfloat)
        minsize = minimum(nout.(terminputs))
        # Or else we might increase the size despite Δ being negative which would be surprising to a user who has specified
        # strictly negative size changes
        minsize + Δ <= 0 && continue

        Δs[v] = Δ
    end

    if !isempty(Δs)
        failmsg = (args...) -> "Could not change nout of $(join(NaiveNASlib.nameorrepr.(keys(Δs)), ", ", " and ")) by $(join(values(Δs), ", ", " and ")). No change!"

        strategy = TimeOutAction(;base=ΔNoutRelaxed(Δs), fallback=LogΔSizeExec(failmsg, Logging.Warn, ΔSizeFailNoOp()))

        Δsize!(strategy)
    end
    return vs
end

"""
    AddVertexMutation <:AbstractMutation{AbstractVertex}
    AddVertexMutation(s::AbstractArchSpace, outselect::Function, WeightInit::AbstractWeightInit, rng::AbstractRNG)
    AddVertexMutation(s, outselect::Function=identity)
    AddVertexMutation(s, rng::AbstractRNG)
    AddVertexMutation(s, wi::AbstractWeightInit)

Insert a vertex from the wrapped `AbstractArchSpace` `s` after a given vertex `v`.

The function `outselect` takes an `AbstractVector{AbstractVertex}` representing the output of `v` and returns an `AbstractVector{AbstractVertex}` which shall be reconnected to the vertex `v'` returned by `s`. Defaults to `identity` meaning all outputs of `v` are reconnected to `v'`.
"""
struct AddVertexMutation{S<:AbstractArchSpace, F, WI<:AbstractWeightInit, RNG<:AbstractRNG} <:AbstractMutation{AbstractVertex}
    s::S
    outselect::F
    weightinit::WI
    rng::RNG
end
AddVertexMutation(s, outselect::Function=identity) = AddVertexMutation(s, outselect, IdentityWeightInit(), rng_default)
AddVertexMutation(s, rng::AbstractRNG) = AddVertexMutation(s, identity, IdentityWeightInit(), rng)
AddVertexMutation(s, wi::AbstractWeightInit) = AddVertexMutation(s, identity, wi, rng_default)

function (m::AddVertexMutation)(v::AbstractVertex)
    insert!(v, vi -> m.s(name(vi), vi, m.rng, outsize=nout(vi), wi=m.weightinit), m.outselect)
    return v
end

"""
    RemoveVertexMutation <:AbstractMutation{AbstractVertex}
    RemoveVertexMutation(s::RemoveStrategy)
    RemoveVertexMutation()

Remove the given vertex `v` using the configured `RemoveStrategy`.

Default size align strategy is `IncreaseSmaller -> DecreaseBigger -> AlignSizeBoth -> FailAlignSizeWarn -> FailAlignSizeRevert`.

Default reconnect strategy is `ConnectAll`.

Note: High likelyhood of large accuracy degradation after applying this mutation.
"""
struct RemoveVertexMutation{S<:RemoveStrategy} <:AbstractMutation{AbstractVertex}
    s::S
end
function RemoveVertexMutation()
    alignstrat = IncreaseSmaller(fallback=DecreaseBigger(fallback=AlignSizeBoth(fallback=FailAlignSizeWarn(msgfun = (vin,vout) -> "Can not remove vertex $(name(vin))! Could not align sizes of neighbours!"))))
    return RemoveVertexMutation(RemoveStrategy(CheckAligned(CheckNoSizeCycle(alignstrat, FailAlignSizeWarn(msgfun = (vin,vout) -> "Can not remove vertex $(name(vin))! Size cycle detected!")))))
end

function (m::RemoveVertexMutation)(v::AbstractVertex)
    remove!(v, m.s)
    return v
end

default_neuronselect(args...) = NaiveNASlib.defaultutility(args...)

"""
    AddEdgeMutation <: AbstractMutation{AbstractVertex}
    AddEdgeMutation(p; rng=rng_default, mergefun=default_mergefun(rng=rng), filtfun=no_shapechange, utilityfun=default_neuronselect)
    AddEdgeMutation(p::Probability; rng=rng_default, mergefun=default_mergefun(rng=rng), filtfun=no_shapechange, utilityfun=default_neuronselect)

Add an edge from a vertex `vi` to another vertex `vo` randomly selected from `vs = filtfun(vi)`.

Higher values of `p` will give more preference to earlier vertices of `vs`.

If `vo` is not capable of having multiple inputs (determined by `singleinput(v) == true`), `vm = mergefun(voi)` where `voi` is a randomly selected input to `vo` will be used instead of `vo` and `vo` will be added as the output of `vm`.

When selecting neurons/outputs after any eventual size change the output of `utilityfun(v)` will be used to determine the utlity of each output in vertex `v`. Note that `length(utilityfun(v)) == nout(v)` must hold.

Note: High likelyhood of large accuracy degradation after applying this mutation.
"""
struct AddEdgeMutation{F1, F2, F3, P<:Probability, RNG} <: AbstractMutation{AbstractVertex}
    mergefun::F1
    filtfun::F2
    utilityfun::F3
    p::P
    rng::RNG
end
AddEdgeMutation(p; rng=rng_default, mergefun=default_mergefun(rng=rng), filtfun=no_shapechange, utilityfun=default_neuronselect) = AddEdgeMutation(Probability(p, rng), rng=rng, mergefun=mergefun, filtfun=filtfun, utilityfun=utilityfun)
AddEdgeMutation(p::Probability; rng=rng_default, mergefun=default_mergefun(rng=rng), filtfun=no_shapechange, utilityfun=default_neuronselect) = AddEdgeMutation(mergefun, filtfun, utilityfun, p, rng)

default_mergefun(pconc = 0.5; rng=rng_default, traitfun = MutationShield ∘ RemoveIfSingleInput ∘ validated() ∘ default_logging(), layerfun = ActivationContribution) = function(vin)
    if rand(rng) > pconc
        return invariantvertex(layerfun(+), vin, traitdecoration=traitfun ∘ named(name(vin) * ".add"))
    end
    return concat(vin, traitfun = traitfun ∘ named(name(vin) * ".cat"), layerfun=layerfun)
end

function no_shapechange(vi)
    # all_in_graph is not sorted, and we want some kind of topoligical order here so that earlier indices are closer to vi
    allsorted = mapreduce(ancestors, vcat, filter(v -> isempty(outputs(v)), all_in_graph(vi))) |> unique

    # Vertices which have the same input as vi and are singleinput
    #   Reason is that this will cause a new vertex to be added between the target output vertex vo
    #   and the input vertex to vi (vii) and this is detected as a size cycle which causes
    #   try_add_edge to fail.
    inouts = filter(singleinput, mapreduce(outputs, vcat, inputs(vi); init=[]))
    # All vertices which are after vi in the topology
    vsafter = setdiff(allsorted, ancestors(vi), outputs(vi), inouts)

    vitrace = shapetrace(vi)
    viorder = allΔshapetypes(vitrace)
    viΔshape = squashshapes(vitrace; order=viorder)

    return filter(vsafter) do vafter
        all(inputs(vafter)) do v
        t = shapetrace(v)
        vΔshape = squashshapes(t; order=union(viorder, allΔshapetypes(t)))
        return viΔshape == vΔshape
        end
    end
end

function (m::AddEdgeMutation)(vi::AbstractVertex)
    # All vertices for which it is allowed to add vi as an input
    allverts = filter(allow_mutation, m.filtfun(vi))
    isempty(allverts) && return vi

    # Higher probability to select a vertex close to v is desired behaviour
    # One line less than a for loop => FP wins!!
    selfun(::Nothing, vc) = apply(m.p) ? vc : nothing
    selfun(vs, vd) = vs
    vo = foldl(selfun, allverts, init=nothing)
    vo = isnothing(vo) ? rand(m.rng, allverts) : vo

    try_add_edge(vi, vo, m.mergefun, m.utilityfun)
    return vi
end

function try_add_edge(vi, vo, mergefun, utilityfun=default_neuronselect)

    # Need to add a vertex which can handle multiple inputs if vo is single input only
    # For cleaning up added vertex if the whole operation fails
    cleanup_failed = () -> nothing
    if singleinput(vo)
        voi = inputs(vo)[1]
        # If the input to vo is capable of multi input we don't need to create a new vertex
        # We must also check that this input does not happen to be an input to vi as this would create a cycle in the graph
        if singleinput(voi) || voi in ancestors(vi)
            vm = mergefun(voi)
            # Insert vm between voi and vo, i.e voi -> vo turns into voi -> vm -> vo
            # vs -> [vo] means only add the new vertex between voi and vo as voi could have other outputs
            insert!(voi, vv -> vm, vs -> [vo])
            cleanup_failed = function()
                length(inputs(vm)) > 1 && return
                remove!(vm, RemoveStrategy(NoSizeChange()))
            end
            vo = vm # vm is the one we shall add an edge to
            @debug "Create new vertex for merging $(name(vo))"
        else
            vo = voi
        end
    end
    # This is mainly because FailAlignSizeRevert does not work when the same vertex is input more than once, but it also seems kinda redundant.
    vi in inputs(vo) && return
    @debug "Create edge between $(name(vi)) and $(name(vo))"
    create_edge!(vi, vo, strategy = create_edge_strat(vo, utilityfun))
    cleanup_failed()
end
# Need to override this one for strange types e.g. layers which support exactly 2 inputs or something.
singleinput(v) = isempty(inputs(v)) || length(inputs(v)) == 1

create_edge_strat(v::AbstractVertex, utilityfun) = create_edge_strat(trait(v), utilityfun)
create_edge_strat(d::DecoratingTrait, utilityfun) = create_edge_strat(base(d), utilityfun)
function create_edge_strat(::SizeInvariant, utilityfun)
    warnfailalign = FailAlignSizeWarn(msgfun = (vin,vout) -> "Could not align sizes of $(name(vin)) and $(name(vout))!")
    alignstrat = AlignSizeBoth(;mapstrat=WithUtilityFun(utilityfun), fallback = warnfailalign)
    # Tricky failure case: It is possible that CheckCreateEdgeNoSizeCycle does not detect any size cycle until after the edge has been created?
    sizecyclewarn = FailAlignSizeWarn(msgfun = (vin,vout) -> "Could not align sizes of $(name(vin)) and $(name(vout))! Size cycle detected!")

    return CheckCreateEdgeNoSizeCycle(ifok=alignstrat, ifnok=sizecyclewarn)
end
function create_edge_strat(::SizeStack, utilityfun)
    warnfailalign = FailAlignSizeWarn(msgfun = (vin,vout) -> "Could not align sizes of $(name(vin)) and $(name(vout))!")
    alignstrat = PostAlign(TruncateInIndsToValid(WithUtilityFun(utilityfun, AlignNinToNout(;fallback=ΔSizeFailNoOp()))), fallback=warnfailalign)

    sizecyclewarn = FailAlignSizeWarn(msgfun = (vin,vout) -> "Could not align sizes of $(name(vin)) and $(name(vout))! Size cycle detected!")
    return CheckCreateEdgeNoSizeCycle(ifok=alignstrat, ifnok=sizecyclewarn)
end

"""
    RemoveEdgeMutation <: AbstractMutation{AbstractVertex}
    RemoveEdgeMutation(;utilityfun=default_neuronselect, rng=rng_default)

Remove an edge from a vertex `vi` to another vertex `vo` randomly selected from `outputs(vi)`.

Vertex `vi` must have more than one output and vertex `vo` must have more than one output for the edge to be removed. Otherwise no change is made.

If there are multiple edges between `vi` and `vo` no change will be made due to NaiveNASlib not being able to revert a failed operation in this case..

When selecting neurons/outputs after any eventual size change the output of `utilityfun(v)` will be used to determine the utlity of each output in vertex `v`. Note that `length(utilityfun(v)) == nout(v)` must hold.

Note: High likelyhood of large accuracy degradation after applying this mutation.
"""
struct RemoveEdgeMutation{F, RNG<:AbstractRNG} <: AbstractMutation{AbstractVertex}
    utilityfun::F
    rng::RNG
end
RemoveEdgeMutation(;utilityfun=default_neuronselect, rng=rng_default) = RemoveEdgeMutation(utilityfun, rng)

function (m::RemoveEdgeMutation)(vi::AbstractVertex)
    length(outputs(vi)) < 2 && return vi

    allverts = filter(vo -> length(inputs(vo)) > 1, outputs(vi))

    isempty(allverts) && return vi

    vo = rand(m.rng, allverts)
    sum(inputs(vo) .== vi) > 1 && return vi# Not implemented in NaiveNASlib

    @debug "Remove edge between $(name(vi)) and $(name(vo))"
    remove_edge!(vi, vo, strategy=remove_edge_strat(vo, m.utilityfun))
    return vi
end

remove_edge_strat(v::AbstractVertex, utilityfun) = remove_edge_strat(trait(v), utilityfun)
remove_edge_strat(d::DecoratingTrait, utilityfun) = remove_edge_strat(base(d), utilityfun)
remove_edge_strat(::SizeInvariant, utilityfun) = NoSizeChange()
remove_edge_strat(t::SizeStack, utilityfun) = create_edge_strat(t, utilityfun)

"""
    KernelSizeMutation{N} <: AbstractMutation{AbstractVertex}
    KernelSizeMutation(Δsizespace::AbstractParSpace{N, Int}; maxsize, pad, rng)
    KernelSizeMutation2D(absΔ::Integer;maxsize, pad, rng)
    KernelSizeMutation(absΔ::Integer...;maxsize, pad, rng)

Mutate the size of filter kernels of convolutional layers.

Note: High likelyhood of large accuracy degradation after applying this mutation.

`KernelSizeMutation2D` is a convenience constructor for `KernelSizeMutation(absΔ, absΔ;...)`.
"""
struct KernelSizeMutation{N,F,P} <: AbstractMutation{AbstractVertex}
    Δsizespace::AbstractParSpace{N, Int}
    maxsize::F
    pad::P
    rng::AbstractRNG
end
KernelSizeMutation(Δsizespace::AbstractParSpace{N, Int}; maxsize = v -> ntuple(i->Inf,N), pad=SamePad(), rng=rng_default) where N = KernelSizeMutation(Δsizespace, maxsize, pad, rng)
KernelSizeMutation2D(absΔ::Integer;maxsize = v -> (Inf,Inf), pad=SamePad(), rng=rng_default) = KernelSizeMutation(absΔ, absΔ, maxsize = maxsize, pad=pad, rng=rng)
KernelSizeMutation(absΔ::Integer...;maxsize = v -> ntuple(i->Inf, length(absΔ)), pad=SamePad(), rng=rng_default) = KernelSizeMutation(ParSpace(UnitRange.(.-absΔ, absΔ));maxsize = maxsize, pad=pad, rng=rng)

function (m::KernelSizeMutation{N})(v::AbstractVertex) where N
    layertype(v) isa FluxConvolutional{N} || return
    l = layer(v)

    currsize = size(NaiveNASflux.weights(l))[1:N]
    Δsize = Int.(clamp.(m.Δsizespace(m.rng), 1 .- currsize, m.maxsize(v) .- currsize)) # ensure new size is > 0 and < maxsize
    # This will eventually boil down to Setfield doing its thing, and that won't be using any convenience constructors
    pad = Flux.calc_padding(typeof(l), m.pad, currsize .+ Δsize, dilation(l), stride(l))
    KernelSizeAligned(Δsize, pad)(v)
    return v
end
dilation(l) = l.dilation
stride(l) = l.stride

"""
    ActivationFunctionMutation{T,R} <: AbstractMutation{AbstractVertex} where {T <: AbstractParSpace{1}, R <: AbstractRNG}
    ActivationFunctionMutation(actspace::AbstractParSpace{1}, rng::AbstractRNG)
    ActivationFunctionMutation(acts...;rng=rng_default)
    ActivationFunctionMutation(acts::AbstractVector;rng=rng_default)

Mutate the activation function of layers which have an activation function.

Note: High likelyhood of large accuracy degradation after applying this mutation.
"""
struct ActivationFunctionMutation{T,RNG} <: AbstractMutation{AbstractVertex} where {T <: AbstractParSpace{1}, R <: AbstractRNG}
    actspace::T
    rng::RNG
end
ActivationFunctionMutation(acts...;rng=rng_default) = ActivationFunctionMutation(collect(acts), rng=rng)
ActivationFunctionMutation(acts::AbstractVector;rng=rng_default) = ActivationFunctionMutation(ParSpace(acts), rng)

function (m::ActivationFunctionMutation)(v::AbstractVertex)
    m(layertype(v), v)
    return v
end
function (m::ActivationFunctionMutation)(t, v) end
(m::ActivationFunctionMutation)(::Union{FluxDense, FluxConvolutional}, v) = NaiveNASflux.setlayer!(v, (σ = m.actspace(m.rng),))
(m::ActivationFunctionMutation)(::FluxParNorm, v) = NaiveNASflux.setlayer!(v, (λ = m.actspace(m.rng),))
function (m::ActivationFunctionMutation)(::FluxRnn, v)
    newcell = setproperties(layer(v).cell, (σ = m.actspace(m.rng),))
    NaiveNASflux.setlayer!(v, (cell = newcell,))
end

"""
    PostMutation{T} <: DecoratingMutation{T}
    PostMutation(actions, m::AbstractMutation{T})
    PostMutation(m::AbstractMutation{T}, actions...)

Performs a set of actions after a wrapped `AbstractMutation` is applied.

Actions will be invoked with arguments (m::PostMutation{T}, e::T) where m is the enclosing `PostMutation` and `e` is the mutated entity of type `T`.
"""
struct PostMutation{T,A} <: DecoratingMutation{T}
    actions::A
    m::AbstractMutation{T}
end
PostMutation(m::AbstractMutation{T}, actions...) where T = PostMutation(actions, m)
PostMutation(action::Function, m::AbstractMutation{T}) where T = PostMutation(m, action)
function (m::PostMutation{T})(e::T; next=m.m, noop=identity) where T
    eout = next(e)
    foreach(a -> a(m, eout), m.actions)
    return eout
end

"""
    RemoveZeroNout()
    RemoveZeroNout(fallback)

Search for vertices with zero output size and remove them and all of their input vertices if possible to do so witout removing an input or output vertex.

Removal is only possible if a vertex is inside a parallel path which will later be concatenated.
"""
struct RemoveZeroNout
    fallback
end
RemoveZeroNout() = RemoveZeroNout(IncreaseZeroNout())
struct IncreaseZeroNout end

(r::RemoveZeroNout)(m, e) = r(e)
function (r::RemoveZeroNout)(g::CompGraph)
    topoligical_order = vertices(g)
    for v in topoligical_order
        nout(v) == 0 || continue

        # Beware! This turned out to be a lot harder than I first thought.
        # I'm not sure the algorithm works (or realizes that it won't work) for all possible cases.

        # This obviously only works if v is inside a parallel path which will later be concatenated

        # To make sure it is, we look ahead in forward and backward direction.
        # If we don't see 1) an input vertex 2) a vertex without outputs (i.e an output vertex) we are good to go
        fseen, fpoints = findforkpoint(v, topoligical_order, inputs, outputs)
        isempty(fpoints) && continue

        bseen, bpoints = findforkpoint(v, topoligical_order, outputs, inputs)
        isempty(bpoints) && continue

        # Ok, fpoints are all vertices where forks with v in them join and bpoints are all vertices where forks with v in them begin

        # If we start removing all seen vertices which are input to an fpoint until we hit a bpoint we should have removed the fork with v in it, right?
        seen = union(fseen, bseen)
        to_rm = intersect(vcat(inputs.(fpoints)...), seen)
        foreach(fpoint -> remove_all_inputs(fpoint, seen, bpoints), to_rm)
    end
end

function findforkpoint(v::AbstractVertex, topoligical_order, f1=inputs, f2=outputs, seen = vcat(v, f1(v)), points = AbstractVertex[])

    if all(x -> x in seen, f1(v))

        # Check if we came across a vertex which previously was thought to be a forkpoint and remove it if so
        if v in points
            deleteat!(points, indexin([v], points))
        end
        # Always visit in reverse topoligical order to mitigate chasing down paths just because we haven't explored them yet
        nextverts = reverse(topoligical_order[unique(indexin(f2(v), topoligical_order))])
        push!(seen, filter(v2 -> !in(v2, seen), nextverts)...)
        foreach(v2 -> findforkpoint(v2, topoligical_order, f1, f2, seen, points), nextverts)
    elseif !(v in points)
        push!(points, v)
    end

    return seen, points
end

function remove_all_inputs(v, seen, stop)
    v in stop && return
    foreach(vrm -> remove_all_inputs(vrm, seen, stop), intersect(inputs(v), seen))
    remove!(v, RemoveStrategy(ConnectNone(), NoSizeChange()))
end

"""
    struct OptimizerMutation{F} <: AbstractMutation{FluxOptimizer}
    OptimizerMutation(optfun)
    OptimizerMutation(os::Union{Tuple, <:AbstractArray})

Mutatates optimizers not wrapped in `ShieldedOpt` through `optfun`.

Invoked recursively for `Flux.Optimiser`s.
"""
struct OptimizerMutation{F} <: AbstractMutation{FluxOptimizer}
    optfun::F
end
OptimizerMutation(os::Union{Tuple, <:AbstractArray}, rng=rng_default) = OptimizerMutation(o -> rand(rng, os)(learningrate(o)))

"""
    LearningRateMutation(rng=rng_default)

Return an `OptimizerMutation` which mutates the learning rate of optimizers.
"""
LearningRateMutation(rng=rng_default) = OptimizerMutation(o -> nudgelr(o, rng))

(m::OptimizerMutation)(opt::Flux.Optimiser) = Flux.Optimiser(m.(opt.os))
(m::OptimizerMutation)(o::ShieldedOpt) = o;
(m::OptimizerMutation)(o::FluxOptimizer) = m.optfun(o)

nudgelr(o, rng=rng_default) = sameopt(o, nudgelr(learningrate(o), rng))
nudgelr(lr::Number, rng=rng_default) = clamp(lr + (rand(rng) - 0.5) * lr * 0.3, 1e-6, 1.0)

learningrate(o::Flux.Optimiser) = prod(learningrate.(o.os))
learningrate(o::ShieldedOpt) = learningrate(o.opt)
learningrate(o) = o.eta

newlr(o, lrf = nudgelr) = sameopt(o, lrf(learningrate(o)))
sameopt(o, lr) = @set o.eta = lr

"""
    AddOptimizerMutation{F} <: AbstractMutation{FluxOptimizer}

Adds optimizer generated by `optgen(os)` to the set of optimizers where `os` is the existing set.

An attempt to merge optimizers of the same type is made using `mergeopts`.
"""
struct AddOptimizerMutation{F} <: AbstractMutation{FluxOptimizer}
    optgen::F
end
(m::AddOptimizerMutation)(o::ShieldedOpt) = o;
(m::AddOptimizerMutation)(o::FluxOptimizer) = m(Flux.Optimiser([o]))
function (m::AddOptimizerMutation)(opt::Flux.Optimiser)
    newopt = m.optgen(opt)
    return Flux.Optimiser(mergeopts(typeof(newopt), newopt, opt.os...))
end


text:
# src/SubTypes.jl
# Custom Subtyping in Julia

__precompile__(true)

"""```
SubTypes
```
Custom Subtypes in Julia.
See https://github.com/bhgomes/SubTypes.jl for more details.
"""
module SubTypes

import Base: convert, eltype

export SubType, check_predicate, predicate, context, Constrained, ConstrainedSymbol, support

"""```
SubType{T, P, Ctx}
```
Subtype structure with underlying type `T`, predicate `P`, and context `Ctx`.
"""
struct SubType{T,P,Ctx}
    value::T
    SubType{T}(x::T) where {T} = new{T,Any,nothing}(x)
    SubType{T,Any}(x::T) where {T} = new{T,Any,nothing}(x)
    function SubType{T,P,Ctx}(x::T) where {T,P,Ctx}
        check_predicate(P, Val(Ctx), x)
        return new{T,P,Ctx}(x)
    end
end

"""```
check_predicate(P, Val(Ctx), x)
```
Check that `x` satisfies the predicate `P` in context `Ctx`.
"""
function check_predicate(P, ::Val{Ctx}, x) where {Ctx}
    throw(ArgumentError("`$x` does not match the predicate `$P` in context `$Ctx`"))
end

"""```
convert(::Type{T}, x::SubType{T}) === x.value
```
Convert a subtype value `x` to its underlying value `x.value`.
"""
convert(::Type{T}, x::SubType{T}) where {T} = x.value

"""```
eltype(::Type{SubType{T}}) === T
```
Return the underlying type of the `SubType` type.
"""
eltype(::Type{SubType{T}}) where {T} = T

"""```
eltype(x::SubType) === eltype(typeof(x))
```
Return the underlying type of the `SubType` value `x`.
"""
eltype(x::SubType) = eltype(typeof(x))

"""```
predicate(::Type{SubType{T, P}}) === P
```
Return the `predicate` of the `SubType` type.
"""
predicate(::Type{SubType{T,P}}) where {T,P} = P

"""```
predicate(x::SubType) === predicate(typeof(x))
```
Return the `predicate` of the `SubType` value `x`.
"""
predicate(x::SubType) = predicate(typeof(x))

"""```
context(::Type{SubType{T, P, Ctx}}) === Ctx
```
Return the `context` of the `SubType` type.
"""
context(::Type{SubType{T,P,Ctx}}) where {T,P,Ctx} = Ctx

"""```
context(x::SubType) === context(typeof(x))
```
Return the `context` of the `SubType` value `x`.
"""
context(x::SubType) = context(typeof(x))

"""```
Constrained{T, S} === SubType{T, S, :Constrained}
```
Constrained types emulate subset inclusion types. That is
`x::Constrained{T,S} <=> x.value::T in S`.
"""
const Constrained{T,S} = SubType{T,S,:Constrained}

"""```
ConstrainedSymbol{S} === Constrained{Symbol, S}
```
Constrained symbols are the constrained type of `Symbol`s.
"""
const ConstrainedSymbol{S} = Constrained{Symbol,S}

"""```
check_predicate(S, ::Val{:Constrained}, x::T)
```
Predicate check implementation for `Constrained` types.
"""
function check_predicate(S, ::Val{:Constrained}, x::T) where {T}
    if !(x in S)
        throw("$x not in support $S")
    end
end

"""```
support(::Type{Constrained{T, S}}) = S
```
Return the set which defines the predicate, i.e. the support of the
indicator function on `S`.
"""
support(::Type{Constrained{T,S}}) where {T,S} = S

"""```
support(constrained::Constrained) === support(typeof(constrained))
```
Return the set which defines the predicate, i.e. the support of the
indicator function on `S`.
"""
support(constrained::Constrained) = support(typeof(constrained))

end  # module SubTypes


text:
module FixedPointNumbers

import Base: ==, <, <=, -, +, *, /, ~, isapprox,
             convert, promote_rule, show, bitstring, abs, decompose,
             isnan, isinf, isfinite, isinteger,
             zero, oneunit, one, typemin, typemax, floatmin, floatmax, eps, reinterpret,
             big, rationalize, float, trunc, round, floor, ceil, bswap, clamp,
             div, fld, cld, rem, mod, mod1, fld1, min, max, minmax,
             signed, unsigned, copysign, flipsign, signbit,
             length

import Statistics   # for _mean_promote
import Random: Random, AbstractRNG, SamplerType, rand!

import Base.Checked: checked_neg, checked_abs, checked_add, checked_sub, checked_mul,
                     checked_div, checked_fld, checked_cld, checked_rem, checked_mod

using Base: @pure

"""
    FixedPoint{T <: Integer, f} <: Real

Supertype of the two fixed-point number types: `Fixed{T, f}` and `Normed{T, f}`.

The parameter `T` is the underlying machine representation and `f` is the number
of fraction bits.
"""
abstract type FixedPoint{T <: Integer, f} <: Real end

export
    FixedPoint,
    Fixed,
    Normed,
    floattype,
# "special" typealiases
    # Q and N typealiases are exported in separate source files
# Functions
    scaledual,
    wrapping_neg, wrapping_abs, wrapping_add, wrapping_sub, wrapping_mul,
    wrapping_div, wrapping_fld, wrapping_cld, wrapping_rem, wrapping_mod,
    saturating_neg, saturating_abs, saturating_add, saturating_sub, saturating_mul,
    saturating_div, saturating_fld, saturating_cld, saturating_rem, saturating_mod,
    wrapping_fdiv, saturating_fdiv, checked_fdiv

include("utilities.jl")

# reinterpretation
reinterpret(x::FixedPoint) = x.i
reinterpret(::Type{T}, x::FixedPoint{T,f}) where {T,f} = x.i
reinterpret(::Type{X}, x::T) where {T <: Integer, X <: FixedPoint{T}} = X(x, 0)

# static parameters
nbitsfrac(::Type{X}) where {T, f, X <: FixedPoint{T,f}} = f
rawtype(::Type{X}) where {T, X <: FixedPoint{T}} = T

# traits based on static parameters
signbits(::Type{X}) where {T, X <: FixedPoint{T}} = T <: Unsigned ? 0 : 1
nbitsint(::Type{X}) where {X <: FixedPoint} = bitwidth(X) - nbitsfrac(X) - signbits(X)

# construction using the (approximate) intended value, i.e., N0f8
*(x::Real, ::Type{X}) where {X <: FixedPoint} = _convert(X, x)
wrapping_mul(x::Real, ::Type{X}) where {X <: FixedPoint} = x % X
saturating_mul(x::Real, ::Type{X}) where {X <: FixedPoint} = clamp(x, X)
checked_mul(x::Real, ::Type{X}) where {X <: FixedPoint} = _convert(X, x)

# type modulus
rem(x::Real, ::Type{X}) where {X <: FixedPoint} = _rem(x, X)
wrapping_rem(x::Real, ::Type{X}) where {X <: FixedPoint} = _rem(x, X)
saturating_rem(x::Real, ::Type{X}) where {X <: FixedPoint} = _rem(x, X)
checked_rem(x::Real, ::Type{X}) where {X <: FixedPoint} = _rem(x, X)

# constructor-style conversions
(::Type{X})(x::X) where {X <: FixedPoint}      = x
(::Type{X})(x::Number) where {X <: FixedPoint} = _convert(X, x)

function (::Type{<:FixedPoint})(x::AbstractChar)
    throw(ArgumentError("FixedPoint (Fixed or Normed) cannot be constructed from a Char"))
end
(::Type{X})(x::Complex) where {X <: FixedPoint} = X(convert(real(typeof(x)), x))
function (::Type{X})(x::Base.TwicePrecision) where {X <: FixedPoint}
    floattype(X) === BigFloat ? X(big(x)) : X(convert(floattype(X), x))
end

# conversions
function Base.Bool(x::FixedPoint)
    x == zero(x) ? false : x == oneunit(x) ? true : throw(InexactError(:Bool, Bool, x))
end
function (::Type{Ti})(x::FixedPoint) where {Ti <: Integer}
    isinteger(x) || throw(InexactError(:Integer, typeof(x), x))
    floor(Ti, x)
end
Base.Rational{Ti}(x::FixedPoint) where {Ti <: Integer} = Rational{Ti}(Rational(x))

big(::Type{<:FixedPoint}) = BigFloat
big(x::FixedPoint) = convert(BigFloat, x)

rationalize(x::FixedPoint; tol::Real=eps(x)) = rationalize(Int, x, tol=tol)
function rationalize(::Type{Ti}, x::FixedPoint; tol::Real=eps(x)) where Ti <: Integer
    tol <= eps(x) ? Rational{Ti}(x) : rationalize(Ti, float(x), tol)
end

"""
    isapprox(x::FixedPoint, y::FixedPoint; rtol=0, atol=max(eps(x), eps(y)))

For FixedPoint numbers, the default criterion is that `x` and `y` differ by no
more than `eps`, the separation between adjacent fixed-point numbers.
"""
@inline function isapprox(x::X, y::X; rtol=0, atol=eps(X)) where {X <: FixedPoint}
    n, m = minmax(x, y) # m >= n
    if rtol == zero(rtol)
        _isapprox_atol(m, n, atol)
    elseif atol == zero(atol)
        _isapprox_rtol(m, n, rtol)
    else
        _isapprox_atol(m, n, atol) | _isapprox_rtol(m, n, rtol)
    end
end
function isapprox(x::FixedPoint, y::FixedPoint; rtol=0, atol=max(eps(x), eps(y)))
    isapprox(promote(x, y)...; rtol=rtol, atol=atol)
end
function _isapprox_atol(m::X, n::X, atol::X) where {X <: FixedPoint}
    unsigned(m.i - n.i) <= unsigned(max(atol, zero(X)).i)
end
function _isapprox_atol(m::X, n::X, atol) where {X <: FixedPoint}
    unsigned(m.i - n.i) <= div(atol, eps(X))
end
function _isapprox_rtol(m::X, n::X, rtol) where {X <: FixedPoint}
    unsigned(m.i - n.i) <= rtol * max(unsigned(abs(m.i)), unsigned(abs(n.i)))
end

# predicates
isinteger(x::FixedPoint) = x == trunc(x) # TODO: use floor(x) when dropping support for Fixed{Int8,8}
isfinite(x::FixedPoint) = true
isnan(x::FixedPoint) = false
isinf(x::FixedPoint) = false

# identities
zero(::Type{X}) where {X <: FixedPoint} = X(zero(rawtype(X)), 0)
oneunit(::Type{X}) where {X <: FixedPoint} = X(rawone(X), 0)
one(::Type{X}) where {X <: FixedPoint} = oneunit(X)

# for Julia v1.0, which does not fold `div_float` before inlining
inv_rawone(x) = (@generated) ? (y = 1.0 / rawone(x); :($y)) : 1.0 / rawone(x)

# traits
eps(::Type{X}) where {X <: FixedPoint} = X(oneunit(rawtype(X)), 0)
typemax(::Type{T}) where {T <: FixedPoint} = T(typemax(rawtype(T)), 0)
typemin(::Type{T}) where {T <: FixedPoint} = T(typemin(rawtype(T)), 0)
floatmin(::Type{T}) where {T <: FixedPoint} = eps(T)
floatmax(::Type{T}) where {T <: FixedPoint} = typemax(T)

"""
    floattype(::Type{T})::Type{<:AbstractFloat}

Return a minimal type suitable for performing computations with instances of type `T` without integer overflow.

The fallback definition of `floattype(T)` applies only to `T<:AbstractFloat`.
However, it is permissible to extend `floattype` to return types that are not subtypes of
`AbstractFloat`; the key characteristic is that the return type should support computation without integer overflow.

In general the returned type should have the minimum bitwidth needed to encode the full precision of the input type.
however, a priority should be placed on computational efficiency; consequently, types like `Float16` should be avoided
except in scenarios where they are guaranteed to have hardware support.

# Examples

A classic usage is to avoid overflow behavior by promoting `FixedPoint` to `AbstractFloat`

```jldoctest; setup = :(using FixedPointNumbers)
julia> x = N0f8(1.0)
1.0N0f8

julia> x + x # overflow
0.996N0f8

julia> T = floattype(x)
Float32

julia> T(x) + T(x)
2.0f0
```

The following represents a valid extension of `floattype` to non-AbstractFloats:

```julia
julia> using FixedPointNumbers, ColorTypes

julia> floattype(RGB{N0f8})
RGB{Float32}
```

`RGB` itself is not a subtype of `AbstractFloat`, but unlike `RGB{N0f8}` operations with `RGB{Float32}` are not subject to integer overflow.
"""
floattype(::Type{T}) where {T <: AbstractFloat} = T # fallback (we want a MethodError if no method producing AbstractFloat is defined)
floattype(::Type{T}) where {T <: Union{ShortInts, Bool}} = Float32
floattype(::Type{T}) where {T <: Integer} = Float64
floattype(::Type{T}) where {T <: LongInts} = BigFloat
floattype(::Type{T}) where {I <: Integer, T <: Rational{I}} = typeof(zero(I)/oneunit(I))
floattype(::Type{<:AbstractIrrational}) = Float64
floattype(::Type{X}) where {T <: ShortInts, X <: FixedPoint{T}} = Float32
floattype(::Type{X}) where {T <: Integer, X <: FixedPoint{T}} = Float64
floattype(::Type{X}) where {T <: LongInts, X <: FixedPoint{T}} = BigFloat

# Non-Real types
floattype(::Type{Complex{T}}) where T = Complex{floattype(T)}
floattype(::Type{Base.TwicePrecision{Float64}}) = Float64    # wider would be nice, but hardware support is paramount
floattype(::Type{Base.TwicePrecision{T}}) where T<:Union{Float16,Float32} = widen(T)

float(x::FixedPoint) = convert(floattype(x), x)

# wrapping arithmetic
wrapping_neg(x::X) where {X <: FixedPoint} = X(-x.i, 0)
wrapping_abs(x::X) where {X <: FixedPoint} = X(abs(x.i), 0)
wrapping_add(x::X, y::X) where {X <: FixedPoint} = X(x.i + y.i, 0)
wrapping_sub(x::X, y::X) where {X <: FixedPoint} = X(x.i - y.i, 0)
wrapping_mul(x::X, y::X) where {X <: FixedPoint} = (float(x) * float(y)) % X
function wrapping_fdiv(x::X, y::X) where {X <: FixedPoint}
    z = floattype(X)(x.i) / floattype(X)(y.i)
    isfinite(z) ? z % X : zero(X)
end
function wrapping_div(x::X, y::X, r::RoundingMode = RoundToZero) where {T, X <: FixedPoint{T}}
    z = round(floattype(X)(x.i) / floattype(X)(y.i), r)
    isfinite(z) || return zero(T)
    if T <: Unsigned
        _unsafe_trunc(T, z)
    else
        z > typemax(T) ? typemin(T) : _unsafe_trunc(T, z)
    end
end
wrapping_fld(x::X, y::X) where {X <: FixedPoint} = wrapping_div(x, y, RoundDown)
wrapping_cld(x::X, y::X) where {X <: FixedPoint} = wrapping_div(x, y, RoundUp)
wrapping_rem(x::X, y::X, r::RoundingMode = RoundToZero) where {T, X <: FixedPoint{T}} =
    X(x.i - wrapping_div(x, y, r) * y.i, 0)
wrapping_mod(x::X, y::X) where {X <: FixedPoint} = wrapping_rem(x, y, RoundDown)

# saturating arithmetic
saturating_neg(x::X) where {X <: FixedPoint} = X(~min(x.i - true, x.i), 0)
saturating_neg(x::X) where {X <: FixedPoint{<:Unsigned}} = zero(X)

saturating_abs(x::X) where {X <: FixedPoint} =
    X(ifelse(signbit(abs(x.i)), typemax(x.i), abs(x.i)), 0)

saturating_add(x::X, y::X) where {X <: FixedPoint} =
    X(x.i + ifelse(x.i < 0, max(y.i, typemin(x.i) - x.i), min(y.i, typemax(x.i) - x.i)), 0)
saturating_add(x::X, y::X) where {X <: FixedPoint{<:Unsigned}} = X(x.i + min(~x.i, y.i), 0)

saturating_sub(x::X, y::X) where {X <: FixedPoint} =
    X(x.i - ifelse(x.i < 0, min(y.i, x.i - typemin(x.i)), max(y.i, x.i - typemax(x.i))), 0)
saturating_sub(x::X, y::X) where {X <: FixedPoint{<:Unsigned}} = X(x.i - min(x.i, y.i), 0)

saturating_mul(x::X, y::X) where {X <: FixedPoint} = clamp(float(x) * float(y), X)

saturating_fdiv(x::X, y::X) where {X <: FixedPoint} =
    clamp(floattype(X)(x.i) / floattype(X)(y.i), X)

function saturating_div(x::X, y::X, r::RoundingMode = RoundToZero) where {T, X <: FixedPoint{T}}
    z = round(floattype(X)(x.i) / floattype(X)(y.i), r)
    isnan(z) && return zero(T)
    if T <: Unsigned
        isfinite(z) ? _unsafe_trunc(T, z) : typemax(T)
    else
        _unsafe_trunc(T, clamp(z, typemin(T), typemax(T)))
    end
end
saturating_fld(x::X, y::X) where {X <: FixedPoint} = saturating_div(x, y, RoundDown)
saturating_cld(x::X, y::X) where {X <: FixedPoint} = saturating_div(x, y, RoundUp)
function saturating_rem(x::X, y::X, r::RoundingMode = RoundToZero) where {T, X <: FixedPoint{T}}
    T <: Unsigned && r isa RoundingMode{:Up} && return zero(X)
    X(x.i - saturating_div(x, y, r) * y.i, 0)
end
saturating_mod(x::X, y::X) where {X <: FixedPoint} = saturating_rem(x, y, RoundDown)

# checked arithmetic
checked_neg(x::X) where {X <: FixedPoint} = checked_sub(zero(X), x)
function checked_abs(x::X) where {X <: FixedPoint}
    abs(x.i) >= 0 || throw_overflowerror_abs(x)
    X(abs(x.i), 0)
end
function checked_add(x::X, y::X) where {X <: FixedPoint}
    r, f = Base.Checked.add_with_overflow(x.i, y.i)
    z = X(r, 0) # store first
    f && throw_overflowerror(:+, x, y)
    z
end
function checked_sub(x::X, y::X) where {X <: FixedPoint}
    r, f = Base.Checked.sub_with_overflow(x.i, y.i)
    z = X(r, 0) # store first
    f && throw_overflowerror(:-, x, y)
    z
end
function checked_mul(x::X, y::X) where {X <: FixedPoint}
    z = float(x) * float(y)
    typemin(X) - eps(X)/2 <= z < typemax(X) + eps(X)/2 || throw_overflowerror(:*, x, y)
    z % X
end
function checked_fdiv(x::X, y::X) where {T, X <: FixedPoint{T}}
    y === zero(X) && throw(DivideError())
    z = floattype(X)(x.i) / floattype(X)(y.i)
    if T <: Unsigned
        z < typemax(X) + eps(X)/2 || throw_overflowerror(:/, x, y)
    else
        typemin(X) - eps(X)/2 <= z < typemax(X) + eps(X)/2 || throw_overflowerror(:/, x, y)
    end
    z % X
end
function checked_div(x::X, y::X, r::RoundingMode = RoundToZero) where {T, X <: FixedPoint{T}}
    y === zero(X) && throw(DivideError())
    z = round(floattype(X)(x.i) / floattype(X)(y.i), r)
    if T <: Signed
        z <= typemax(T) || throw_overflowerror_div(r, x, y)
    end
    _unsafe_trunc(T, z)
end
checked_fld(x::X, y::X) where {X <: FixedPoint} = checked_div(x, y, RoundDown)
checked_cld(x::X, y::X) where {X <: FixedPoint} = checked_div(x, y, RoundUp)
function checked_rem(x::X, y::X, r::RoundingMode = RoundToZero) where {T, X <: FixedPoint{T}}
    y === zero(X) && throw(DivideError())
    fx, fy = floattype(X)(x.i), floattype(X)(y.i)
    z = fx - round(fx / fy, r) * fy
    if T <: Unsigned && r isa RoundingMode{:Up}
        z >= zero(z) || throw_overflowerror_rem(r, x, y)
    end
    X(_unsafe_trunc(T, z), 0)
end
checked_mod(x::X, y::X) where {X <: FixedPoint} = checked_rem(x, y, RoundDown)

# default arithmetic
const DEFAULT_ARITHMETIC = :wrapping

for (op, name) in ((:-, :neg), (:abs, :abs))
    f = Symbol(DEFAULT_ARITHMETIC, :_, name)
    @eval begin
        $op(x::X) where {X <: FixedPoint} = $f(x)
    end
end
for (op, name) in ((:+, :add), (:-, :sub), (:*, :mul))
    f = Symbol(DEFAULT_ARITHMETIC, :_, name)
    @eval begin
        $op(x::X, y::X) where {X <: FixedPoint} = $f(x, y)
    end
end
# force checked arithmetic
/(x::X, y::X) where {X <: FixedPoint} = checked_fdiv(x, y)
div(x::X, y::X, r::RoundingMode = RoundToZero) where {X <: FixedPoint} = checked_div(x, y, r)
fld(x::X, y::X) where {X <: FixedPoint} = checked_div(x, y, RoundDown)
cld(x::X, y::X) where {X <: FixedPoint} = checked_div(x, y, RoundUp)
rem(x::X, y::X) where {X <: FixedPoint} = checked_rem(x, y, RoundToZero)
rem(x::X, y::X, ::RoundingMode{:Down}) where {X <: FixedPoint} = checked_rem(x, y, RoundDown)
rem(x::X, y::X, ::RoundingMode{:Up}) where {X <: FixedPoint} = checked_rem(x, y, RoundUp)
mod(x::X, y::X) where {X <: FixedPoint} = checked_rem(x, y, RoundDown)

function minmax(x::X, y::X) where {X <: FixedPoint}
    a, b = minmax(reinterpret(x), reinterpret(y))
    X(a,0), X(b,0)
end

bitstring(x::FixedPoint) = bitstring(x.i)

bswap(x::X) where {X <: FixedPoint} = sizeof(X) == 1 ? x : X(bswap(x.i), 0)

# At least on Julia v1.5.0 or earlier, the following specialization helps the
# SIMD vectorization. (cf. PR #194)
clamp(x::X, lo::X, hi::X) where {X <: FixedPoint} = X(clamp(x.i, lo.i, hi.i), 0)

clamp(x, ::Type{X}) where {X <: FixedPoint} = clamp(x, typemin(X), typemax(X)) % X

# Since `FixedPoint` is not an integer type, it is not clear in what type
# `signed` and `unsigned` for `FixedPoint` should return values. They should
# currently throw errors in case we support "unsigned Fixed" or "signed Normed"
# in the future. The following "incomplete" code is necessary for Julia v1.0
# etc. to prevent accidental conversion to an integer type.
signed(x::X) where {X <: FixedPoint} = signed(X)(signed(x.i), 0)
unsigned(x::X) where {X <: FixedPoint} = unsigned(X)(unsigned(x.i), 0)

function copysign(x::X, y::Real) where {T, X <: FixedPoint{T}}
    T <: Signed ? X(copysign(x.i, y), 0) : throw_not_a_signed_number_error(x)
end
function flipsign(x::X, y::Real) where {T, X <: FixedPoint{T}}
    T <: Signed ? X(flipsign(x.i, y), 0) : throw_not_a_signed_number_error(x)
end
if copysign(-1, 0x1) !== 1 # for Julia v1.0 and v1.1 (julia #30748)
    copysign(x::X, y::Unsigned) where {T, X <: FixedPoint{T}} = copysign(x, signed(y))
    flipsign(x::X, y::Unsigned) where {T, X <: FixedPoint{T}} = flipsign(x, signed(y))
end
@noinline function throw_not_a_signed_number_error(x)
    throw(ArgumentError("$x is not a signed number."))
end

signbit(x::X) where {X <: FixedPoint} = signbit(x.i)

~(x::X) where {X <: FixedPoint} = X(~x.i, 0)

function round(x::FixedPoint, r::RoundingMode=RoundNearest;
               digits::Union{Nothing, Integer}=nothing,
               sigdigits::Union{Nothing, Integer}=nothing,
               base::Union{Nothing, Integer}=nothing)

    d = digits === nothing ? 0 : Int(digits)

    if base !== nothing && base != 10
        throw(ArgumentError("`base` numbers other than 10 are not supported."))
    elseif sigdigits !== nothing
        throw(ArgumentError("`sigdigits` is not supported."))
    elseif d < 0
        throw(ArgumentError("negative `digits` is not supported."))
    end
    d === 0 ? _round_digits0(x, r) : _round_digits(x, r, d)
end
function _round_digits(x::X, r::RoundingMode, d::Int) where {T, f, X <: FixedPoint{T,f}}
    log10_2 = 0.3010299956639812
    d > floor(Int, ((f + 1) * log10_2)) && return x
    r = round(float(x), r, digits=d)
    typemin(X) - eps(X)/2 <= r < typemax(X) + eps(X)/2 || throw_converterror(X, r)
    clamp(r, X)
end

trunc(x::X) where {X <: FixedPoint{<:Unsigned}} = floor(x)
trunc(::Type{Ti}, x::X) where {X <: FixedPoint{<:Unsigned}, Ti <: Integer} = floor(Ti, x)

for f in (:zero, :oneunit, :one, :eps, :rawone, :rawtype, :floattype)
    @eval begin
        $f(x::FixedPoint) = $f(typeof(x))
    end
end
for f in (:(==), :<, :<=, :fld1)
    @eval begin
        $f(x::X, y::X) where {X <: FixedPoint} = $f(x.i, y.i)
    end
end
for f in (:mod1, :min, :max)
    @eval begin
        $f(x::X, y::X) where {X <: FixedPoint} = X($f(x.i, y.i), 0)
    end
end
for (m, f) in ((:(:Nearest), :round),
               (:(:ToZero), :trunc),
               (:(:Up), :ceil),
               (:(:Down), :floor))
    @eval begin
        _round_digits0(x::FixedPoint, ::RoundingMode{$m}) = $f(x)
        round(::Type{Ti}, x::FixedPoint, ::RoundingMode{$m}) where {Ti <: Integer} = $f(Ti, x)
    end
end

function length(r::StepRange{X,X}) where {X <: FixedPoint{<:ShorterThanInt}}
    start, step, stop = Int(reinterpret(r.start)), Int(reinterpret(r.step)), Int(reinterpret(r.stop))
    return div((stop - start) + step, step)
end
function length(r::StepRange{X,X}) where {X <: FixedPoint}
    start, step, stop = reinterpret(r.start), reinterpret(r.step), reinterpret(r.stop)
    return checked_div(checked_add(checked_sub(stop, start), step), step)
end
function length(r::StepRange{<:FixedPoint})
    start, step, stop = float(r.start), r.step, float(r.stop)
    return div((stop - start) + step, step)
end

hasalias(::Type) = false
hasalias(::Type{X}) where {T<:NotBiggerThanInt64, f, X<:FixedPoint{T,f}} = f isa Int

# `alias_symbol` is used to define type aliases, so we need this before we
# include "src/fixed.jl" / "src/normed.jl".
function alias_symbol(@nospecialize(X))
    Symbol(type_prefix(X), nbitsint(X), 'f', nbitsfrac(X))
end

function _alias_symbol(::Type{X}) where {X <: FixedPoint}
    if @generated
        sym = string(alias_symbol(X))
        return :(Symbol($sym))
    else
        return alias_symbol(X)
    end
end

@inline function showtype(io::IO, ::Type{X}) where {X <: FixedPoint}
    if hasalias(X)
        write(io, _alias_symbol(X))
    else
        print(io, X)
    end
    io
end

function show(io::IO, x::FixedPoint{T,f}) where {T,f}
    compact = get(io, :compact, false)::Bool
    log10_2 = 0.3010299956639812
    digits = min(ceil(Int, f * log10_2), compact ? 6 : typemax(Int))
    val = round(convert(Float64, x), digits=digits)
    if compact || get(io, :typeinfo, Any) === typeof(x)
        show(io, val)
    elseif hasalias(typeof(x))
        show(io, val)
        showtype(io, typeof(x))
    else
        print(io, typeof(x), '(', val, ')')
    end
end

if VERSION < v"1.6.0-DEV.356" # JuliaLang/julia#36107
    function Base.showarg(io::IO, a::Array{X}, toplevel) where {X<:FixedPoint}
        toplevel || print(io, "::")
        print(io, "Array{")
        showtype(io, X)
        print(io, ",$(ndims(a))}")
        toplevel && hasalias(X) && print(io, " with eltype ", X)
    end
end

include("fixed.jl")
include("normed.jl")
include("deprecations.jl")
const UF = (N0f8, N6f10, N4f12, N2f14, N0f16)

# Promotions
promote_rule(::Type{X}, ::Type{Tf}) where {X <: FixedPoint, Tf <: AbstractFloat} =
    promote_type(floattype(X), Tf)

# Note that `Tr` does not always have enough domains.
promote_rule(::Type{X}, ::Type{Tr}) where {X <: FixedPoint, Tr <: Rational} = Tr

promote_rule(::Type{X}, ::Type{Ti}) where {X <: FixedPoint, Ti <: Integer} = floattype(X)

function promote_rule(::Type{X1}, ::Type{X2}) where {T1, f1, X1 <: FixedPoint{T1,f1},
                                                     T2, f2, X2 <: FixedPoint{T2,f2}}
    X = wrapper(X1)
    X !== wrapper(X2) && return promote_type(floattype(X1), floattype(X2))

    f = max(f1, f2)  # ensure we have enough precision
    Tp = promote_type(T1, T2)
    T = (T1 <: Signed || T2 <: Signed) ? signedtype(Tp) : Tp
    # make sure we have enough integer bits
    m = max(nbitsint(X1), nbitsint(X2))
    _widen_rawtype(X{T,f}, m)
end

function _widen_rawtype(::Type{X}, m) where {T, f, X<:FixedPoint{T,f}}
    nbitsint(X) >= m && return X
    Tw = widen1(T)
    T === Tw && return X
    _widen_rawtype(wrapper(X){Tw,f}, m)
end

# Promotions for reductions
const Treduce = Float64
Base.add_sum(x::FixedPoint, y::FixedPoint) = Treduce(x) + Treduce(y)
Base.reduce_empty(::typeof(Base.add_sum), ::Type{F}) where {F<:FixedPoint}  = zero(Treduce)
Base.reduce_first(::typeof(Base.add_sum), x::FixedPoint)   = Treduce(x)
Base.mul_prod(x::FixedPoint, y::FixedPoint) = Treduce(x) * Treduce(y)
Base.reduce_empty(::typeof(Base.mul_prod), ::Type{F}) where {F<:FixedPoint} = one(Treduce)
Base.reduce_first(::typeof(Base.mul_prod), x::FixedPoint)  = Treduce(x)

if isdefined(Statistics, :_mean_promote)
    Statistics._mean_promote(x::Real, y::FixedPoint) = Treduce(y)
end

"""
    sd, ad = scaledual(s::Number, a)

Return `sd` and `ad` such that `sd * ad == s * a`.
When `a` is an array of FixedPoint numbers, `sd*ad` might be faster to compute than `s*a`.
"""
scaledual(b::Number, x::Union{Number,AbstractArray{<:Number}}) = b, x
scaledual(b::Number, x::FixedPoint) = b/rawone(x), reinterpret(x)
scaledual(b::Number, x::AbstractArray{T}) where T <: FixedPoint =
    b/rawone(T), reinterpret(rawtype(T), x)

scaledual(::Type{Tdual}, x::Union{Number,AbstractArray{<:Number}}) where Tdual = oneunit(Tdual), x
scaledual(::Type{Tdual}, x::FixedPoint) where Tdual = convert(Tdual, 1/rawone(x)), reinterpret(x)
scaledual(::Type{Tdual}, x::AbstractArray{T}) where {Tdual, T <: FixedPoint} =
    convert(Tdual, 1/rawone(T)), reinterpret(rawtype(T), x)

@noinline function throw_converterror(::Type{X}, @nospecialize(x)) where X <: FixedPoint
    nbits = bitwidth(rawtype(X))
    io = IOBuffer()
    showtype(io, X)
    print(io, " is ")
    print(io, nbits == 8 ? "an " : "a ", nbits, "-bit type representing ")
    print(io, nbits <= 16 ? string(2^nbits) : "2^$nbits", " values from ")
    print(IOContext(io, :compact=>true), typemin(X), " to ")
    print(IOContext(io, :compact=>true), typemax(X), "; ")
    print(io, "cannot represent ", x)
    throw(ArgumentError(String(take!(io))))
end

@noinline function throw_overflowerror(op::Symbol, @nospecialize(x), @nospecialize(y))
    io = IOBuffer()
    print(io, x, ' ', op, ' ', y, " overflowed for type ")
    showtype(io, typeof(x))
    throw(OverflowError(String(take!(io))))
end
@noinline function throw_overflowerror_abs(@nospecialize(x))
    io = IOBuffer()
    print(io, "abs(", x, ") overflowed for type ")
    showtype(io, typeof(x))
    throw(OverflowError(String(take!(io))))
end
@noinline function throw_overflowerror_div(r::RoundingMode, @nospecialize(x), @nospecialize(y))
    io = IOBuffer()
    op = r === RoundUp ? "cld(" : r === RoundDown ? "fld(" : "div("
    print(io, op, x, ", ", y, ") overflowed for type ", rawtype(x))
    throw(OverflowError(String(take!(io))))
end
@noinline function throw_overflowerror_rem(r::RoundingMode, @nospecialize(x), @nospecialize(y))
    io = IOBuffer()
    print(io, "rem(", x, ", ", y, ", ", r, ") overflowed for type ", typeof(x))
    throw(OverflowError(String(take!(io))))
end

function Random.rand(r::AbstractRNG, ::SamplerType{X}) where X <: FixedPoint
    X(rand(r, rawtype(X)), 0)
end

function rand!(r::AbstractRNG, A::Array{X}, ::SamplerType{X}) where {T, X <: FixedPoint{T}}
    At = unsafe_wrap(Array, reinterpret(Ptr{T}, pointer(A)), size(A))
    Random.rand!(r, At, SamplerType{T}())
    A
end

if VERSION >= v"1.1" # work around https://github.com/JuliaLang/julia/issues/34121
    include("precompile.jl")
    _precompile_()
end

end # module


text:
module RegisterUtilities

export Counter

#### Counter ####
#
# Stolen from Grid.jl. Useful when you want to do more math on the iterator.

struct Counter
    max::Vector{Int}
end
Counter(sz::Tuple) = Counter(Int[sz...])
function Base.iterate(c::Counter)
    N = length(c.max)
    (N == 0 || any(c.max.<=0)) && return nothing
    state = ones(Int,N)
    copy(state), state
end
function Base.iterate(c::Counter, state)
    state[1] += 1
    i = 1
    while state[i] > c.max[i] && i < length(state)
        state[i] = 1
        i += 1
        state[i] += 1
    end
    state[end] > c.max[end] && return nothing
    copy(state), state
end

# Below functions are from RegisterTestUtilities

using RegisterCore, LinearAlgebra

export quadratic, block_center, tighten

function quadratic(m, n, shift, Q)
    A = zeros(m, n)
    c = block_center(m, n)
    cntr = [shift[1]+c[1], shift[2]+c[2]]
    u = zeros(2)
    for j = 1:n, i = 1:m
        u[1], u[2] = i-cntr[1], j-cntr[2]
        A[i,j] = dot(u,Q*u)
    end
    A
end

quadratic(shift, Q, denom::Matrix) = MismatchArray(quadratic(size(denom)..., shift, Q), denom)

function block_center(sz...)
    ntuple(i->sz[i]>>1+1, length(sz))
end

function tighten(A::AbstractArray)
    T = typeof(first(A))
    for a in A
        T = promote_type(T, typeof(a))
    end
    At = similar(A, T)
    copyto!(At, A)
end

end


text:
str = "Hello!"
char = 'a'
println("This is a string: ", str)
println("This is a character: ", char)

# These variable names are not valid:
# more@, 123name, keywords such as struct, const, if, begin

# String and Char operations
# "2" - "1", "eggs" / "easy", "third" + "a charm" are invalid
println("\'a\' * \'b\' = ", 'a' * 'b') # The result is a String
println("\"abc\" * \"def\" = ", "abc" * "def")
println("\"a\"^3 = ", "a"^3)
println("\"abc\"^3 = ", "abc"^3)

# Conversions
println(Int('x'))
println(Char(60))
println(isvalid(Char, 0x110000))
# U+0000 through U+D7FF and U+E000 through U+10FFFF are valid

# Access and indexing
# These return characters
println(str[begin]) # Equivalent to print(str[1])
println(str[3])
println(str[end])
# Strings are immutable, str[2] = 'j' is invalid

# These return new strings
println(str[3:3])
println(str[begin:end-1]) # Inclusve start:end
println(str[1:3])

# This returns a SubString, i.e. a view into a String
substr = SubString(str, 1, 4)
println(substr)
println(typeof(substr))

# String lengths
println(length(str))
println(length(substr))
println(length(char))

# Interpolation
a, b = 2, 5
println("a + b = $(a + b)")
println("$str, $char.")

println("""

This
is
a test

""")
# In the repl, you will get "\nThis\nis\na test\n\n"

# String comparison
println("abracadabra" < "xylophone")
println( "abracadabra" == "xylophone")
println( "Hello, world." != "Goodbye, world.")

# Searching
println(findfirst('o', "xylophone"))
println(findlast('o', "xylophone"))
println(findfirst('z', "xylophone")) # Note: typeof(nothing) is nothing
println(findnext('o', "xylophone", 1))
println(findnext('o', "xylophone", 5))
println(occursin("world", "Hello, world."))
println(occursin("o", "Xylophon"))
println(occursin('o', "Xylophon"))

# Reading input
text = readline()
println("This was the input: ", text)
number = parse(Int64, readline())
println("This was the input: $number of type $(typeof(number))")


text: result(h::Handle) = LibClingo.clingo_solve_handle_get(h)