text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module GaussFactorials

using Nemo, IterTools
using RecordSearch, NumberTheory, Products

export ModuleGaussFactorials, GaussFactorial
export V216919, V001783, V055634, V232980, V232981, V232982, V124441, V124442, V066570
export I193338, F193338, L193338, V193338
export I193339, F193339, L193339, V193339

"""
The Gauß factorial is ``∏_{1 ≤ j ≤ N, j ⊥ n} j``, the product of the positive integers which are ``≤ N`` and are prime to ``n``.

* GaussFactorial, I193338, F193338, L193338, V193338, I193339, F193339, L193339, V193339, V216919, V001783, V055634, V232980, V232981, V232982, V124441, V124442, V066570
"""
const ModuleGaussFactorials = ""

"""
Return ``∏_{1 ≤ j ≤ N, j ⊥ n} j``, the product of the positive integers which are ``≤ N`` and are prime to ``n``.
"""
GaussFactorial(N, n) = ∏([j for j in 1:Int(N) if ⊥(j, n)])

"""
Return ``∏_{1 ≤ j ≤ n, j ⊥ n} j``, the product of the positive integers which are ``≤ n`` and are prime to ``n``.
"""
GaussFactorial(n) = GaussFactorial(n, n)

# V000142(n) = n! = GaussFactorial(n, 1).
# V001147(n) = GaussFactorial(2*n, 2), double factorial of odd numbers.

"""
Return GaussFactorial(N, n).
"""
V216919(N, n) = GaussFactorial(N, n)

"""
Return ``∏_{1 ≤ j ≤ n, j ⊥ n} j``, the product of the positive integers which are ``≤ N`` and are prime to ``n``, a.k.a. the phi-torial of n.
"""
V001783(n) = GaussFactorial(n, n)

"""
Return the 2-adic factorial of n.
"""
V055634(n) = GaussFactorial(n, 2)*(-1)^n

"""
Return GaussFactorial(n, 3).
"""
V232980(n) = GaussFactorial(n, 3)

"""
Return GaussFactorial(n, 5).
"""
V232981(n) = GaussFactorial(n, 5)

"""
Return GaussFactorial(n, 6).
"""
V232982(n) = GaussFactorial(n, 6)

"""
Return GaussFactorial(div(n, 2), n).
"""
V124441(n) = GaussFactorial(div(n, 2), n)

"""
Return GaussFactorial(n, n)/GaussFactorial(div(n, 2), n).
"""
V124442(n) = div(GaussFactorial(n, n), GaussFactorial(div(n, 2), n))

"""
Return GaussFactorial(n, 1)/GaussFactorial(n, n).
"""
V066570(n) = div(GaussFactorial(n, 1), GaussFactorial(n, n))

# -----------------------------------------------------------------
# Indices!

"""
Iterate over the indices of the first ``n`` record values of the Gauß factorial.
"""
I193339(n) = Records(GaussFactorial, n, true, true)

"""
Iterate over indices of the record values of the Gauß factorial which do not exceed ``n`` (``1 ≤ i ≤ n``).
"""
F193339(n) = Records(GaussFactorial, n, false, true)

"""
Return the indices of the first ``n`` record values of the Gauß factorial as an array.
"""
L193339(n) = collect(I193339(n))

"""
Return the index of the ``n``-th record value of the Gauß factorial.
"""
V193339(n) = nth(I193339(n), n)

# -----------------------------------------------------------------
# Values!

"""
Iterate over the first ``n`` record values of the Gauß factorial (``1 ≤ r``).
"""
I193338(n) = Records(GaussFactorial, n, true, false)

"""
Iterate over the record values of the Gauß factorial which do not exceed ``n`` (``1 ≤ i ≤ n``).
"""
F193338(n) = Records(GaussFactorial, n, false, false)

"""
Return the first ``n`` record values of the Gauß factorial as an array.
"""
L193338(n) = collect(I193338(n))

"""
Return the (value of the) ``n``-th record of the Gauß factorial.
"""
V193338(n) = nth(I193338(n), n)

#START-TEST-########################################################

using Test, SeqTests

function test()

    @testset "GaussFactorial" begin

        #@test
        if is_oeis_installed()

            V = [V001783, V124441, V124442, V066570]
            for v in V SeqTest(v, 'V', 1) end

            W = [V055634, V232980, V232981, V232982]
            for w in W SeqTest(w, 'V', 0) end
        end
    end
end

function demo()
    println("\n Iterate over the indices of the first 10 record values of the Gauß factorial.")
    for r in I193339(10) print(r, ", ") end; println("...")

    println("\n Iterate over indices of the record values of the Gauß factorial which do not exceed 10.")
    for r in F193339(10) print(r, ", ") end; println("...")

    println("\n Return the indices of the first 10 record values of the Gauß factorial as an array.")
    println(L193339(10))

    println("\n Return the index of the 9-th record value of the Gauß factorial.")
    println("9 -> ", V193339(9))

    println("\n Iterate over the first 10 record values of the Gauß factorial.")
    for r in I193338(10) print(r, ", ") end; println("...")

    println("\n Iterate over the record values of Gauß factorial which do not exceed 10.")
    for r in F193338(10) print(r, ", ") end; println("...")

    println("\n Return the first 10 record values of the Gauß factorial as an array.")
    println(L193338(10))

    println("\n Return the (value of the) 9-th record of the Gauß factorial.")
    println("9 -> ", V193338(9))
end

function perf()
end

function main()
    test()
    demo()
    perf()
end

main()

end # module

#=
Iterate over the indices of the first 10 record values ofthe Gauß factorial.
1, 3, 4, 5, 7, 9, 11, 13, 17, ...

Iterate over indices of the record values of the Gauß factorial which do not exceed 10.
1, 3, 4, 5, 7, 9, ...

Return the indices of the first 10 record values of the Gauß factorial as an array.
Nemo.fmpz[1, 3, 4, 5, 7, 9, 11, 13, 17]

Return the index of the 9-th record value of the Gauß factorial.
9 -> 13

Iterate over the first 10 record values of the Gauß factorial.
1, 2, 3, 24, 720, 2240, 3628800, 479001600, 20922789888000, ...

Iterate over the record values of Gauß factorial which do
not exceed 10.
1, 2, 3, 24, 720, 2240, ...

Return the first 10 record values of the Gauß factorial as an array.
Nemo.fmpz[1, 2, 3, 24, 720, 2240, 3628800, 479001600, 20922789888000]

Return the (value of the) 9-th record of the Gauß factorial.
9 -> 479001600
=#


text:
export WaveletOp

mutable struct WaveletOp{T,F1<:FuncOrNothing,F2<:FuncOrNothing,F3<:FuncOrNothing} <: AbstractLinearOperator{T,F1,F2,F3}
  nrow   :: Int
  ncol   :: Int
  symmetric :: Bool
  hermitian :: Bool
  prod   :: F1 # apply the operator to a vector
  tprod  :: F2 # apply the transpose operator to a vector
  ctprod :: F3 # apply the transpose conjugate operator to a vector
end

"""
  WaveletOp(shape, wt=wavelet(WT.db2))

returns a `ẀaveletOp <: AbstractLinearOperator`, which performs a Wavelet transform on
a given input array.

# Arguments

* `shape`                 - size of the Array to transform
* (`wt=wavelet(WT.db2)`)  - Wavelet to apply
"""
function WaveletOp(shape, wt=wavelet(WT.db2))
  return WaveletOp{Float64,Function,Nothing,Function}(maximum(shape)^2, prod(shape), false, false
            , x->waveletProd(x,shape,wt)
            , nothing
            , y->waveletCTProd(y,shape,wt) )
end

function waveletProd(x::Vector{T},shape, wt) where T
  if shape[1] != shape[2]
    xSquare = zeros(T, maximum(shape), maximum(shape))
    xSquare[1:shape[1],1:shape[2]] = reshape(x,shape)
  else
    xSquare = reshape(x,shape)
  end
  return vec( dwt(xSquare, wt) )
end

function waveletCTProd(y::Vector{T},shape, wt) where T
  squareSize = (maximum(shape), maximum(shape))
  x = idwt( reshape(y, squareSize), wt)
  return vec( x[1:shape[1],1:shape[2]] )
end

\(A::WaveletOp, x::Vector) = adjoint(A) * x

ext: jl


text:
# functions for building various standard types of graphs

export Complete, Path, Cycle, RandomGraph, RandomRegular, RandomSBM
export RandomTree, code_to_tree
export Grid, Wheel, Cube, BuckyBall, Johnson, Doyle
export Petersen, Kneser, Paley, Knight, Frucht, Hoffman, HoffmanSingleton, Spindle, Golomb, is_unit_distance

"""
`Complete(n)` returns a complete graph with `n` vertices `1:n`.

`Complete(n,m)` returns a complete bipartite graph with `n` vertices
in one part and `m` vertices in the other.

`Complete([n1,n2,...,nt])` returns a complete multipartite graph with
parts of size `n1`, `n2`, ..., `nt`.
"""
function Complete(n::Int)
    G = IntGraph(n)

    for k = 1:n-1
        for j = k+1:n
            add!(G, j, k)
        end
    end
    name(G, "Complete")
    return G
end

# Create a complete bipartite graph
function Complete(n::Int, m::Int)
    G = IntGraph(n + m)
    for u = 1:n
        for v = n+1:n+m
            add!(G, u, v)
        end
    end
    name(G, "Complete($n,$m)")
    return G
end

# Create the complete multipartite graph with given part sizes
function Complete(parts::Array{Int,1})
    # check all part sizes are positive
    for p in parts
        if p < 1
            error("All part sizes must be positive")
        end
    end

    n = sum(parts)
    G = IntGraph(n)

    np = length(parts)
    if np < 2
        return G
    end

    # create table of part ranges
    ranges = Array{Int}(undef, np, 2)   # old ranges = Array(Int,np,2)
    ranges[1, 1] = 1
    ranges[1, 2] = parts[1]
    for k = 2:np
        ranges[k, 1] = ranges[k-1, 2] + 1
        ranges[k, 2] = ranges[k, 1] + parts[k] - 1
    end

    # Add all edges between all parts
    for i = 1:np-1
        for j = i+1:np
            for u = ranges[i, 1]:ranges[i, 2]
                for v = ranges[j, 1]:ranges[j, 2]
                    add!(G, u, v)
                end
            end
        end
    end
    name(G, "Complete($parts)")
    return G
end

# Create a path graph on n vertices
"""
`Path(n)` creates a path graph with `n` vertices named `1:n`.

`Path(array)` creates a path graph with vertices `array[1]`,
`array[2]`, etc.
"""
function Path(n::Int)
    G = IntGraph(n)
    for v = 1:n-1
        add!(G, v, v + 1)
    end
    name(G, "Path")
    set_rot(G)
    d = Dict{Int,Vector{Float64}}()
    for v = 1:n
        d[v] = [v - n / 2, 0]
    end
    embed(G, d)
    return G
end

# Create a path graph from a list of vertices
function Path(verts::Array{T}) where {T}
    G = SimpleGraph{T}()
    n = length(verts)

    if n == 1
        add!(G, verts[1])
    end
    for k = 1:n-1
        add!(G, verts[k], verts[k+1])
    end
    name(G, "Path")
    set_rot(G)
    return G
end

# Create a cycle graph on n vertices
"""
`Cycle(n)` creates a cycle with vertex set `1:n`.
"""
function Cycle(n::Int)
    if n < 3
        error("Cycle requires 3 or more vertices")
    end
    G = Path(n)
    add!(G, 1, n)
    name(G, "Cycle")
    set_rot(G)
    embed(G)
    return G
end

# Create the wheel graph on n vertices: a cycle on n-1 vertices plus
# an additional vertex adjacent to all the vertices on the wheel.
"""
`Wheel(n)` creates a wheel graph with `n` vertices. That is, a cycle
with `n-1` vertices `1:(n-1)` all adjacent to a common single vertex,
`n`.
"""
function Wheel(n::Int)
    if n < 4
        error("Wheel graphs must have at least 4 vertices")
    end
    G = Cycle(n - 1)
    for k = 1:n-1
        add!(G, k, n)
    end
    name(G, "Wheel")
    embed(G, :tutte, outside = collect(1:n-1))
    embed_rot(G)
    return G
end

# Create a grid graph
"""
`Grid(n,m)` creates an `n`-by-`m` grid graph. For other grids, we
suggest `Path(n1)*Path(n2)*Path(n3)` optionally wrapped in
`relabel`. See also: `Cube`.
"""
function Grid(n::Int, m::Int)
    G = SimpleGraph{Tuple{Int,Int}}()

    # add the vertices
    for u = 1:n
        for v = 1:m
            add!(G, (u, v))
        end
    end

    #horizontal edges
    for u = 1:n
        for v = 1:m-1
            add!(G, (u, v), (u, v + 1))
        end
    end

    # vertical edges
    for v = 1:m
        for u = 1:n-1
            add!(G, (u, v), (u + 1, v))
        end
    end

    # embedding
    d = Dict{Any,Array{Float64,1}}()
    for v in G.V
        x, y = v
        d[v] = [Float64(x), Float64(y)]
    end
    embed(G, d)
    recenter(G)
    embed_rot(G)

    name(G, "Grid($n,$m)")
    return G
end

# Create an Erdos-Renyi random graph
"""
`RandomGraph(n,p=0.5)` creates an Erdos-Renyi random graph with `n`
vertices and edge probability `p`.
"""
function RandomGraph(n::Int, p::Real = 0.5)
    G = IntGraph(n)

    # guess the size of the edge set to preallocate storage
    m = round(Int, n * n * p) + 1

    # generate the edges
    for v = 1:n-1
        for w = v+1:n
            if (rand() < p)
                add!(G, v, w)
            end
        end
    end
    name(G, "Erdos-Renyi, p=$p")
    return G
end

include("prufer.jl")

# Generate a random tree on vertex set 1:n. All n^(n-2) trees are
# equally likely.

"""
`RandomTree(n)` creates a random tree on `n` vertices each with
probability `1/n^(n-2)`.
"""
function RandomTree(n::Int)
    if n < 0   # but we allow n==0 to give empty graph
        error("Number of vertices cannot be negative")
    end

    if n < 2
        return IntGraph(n)
    end

    code = [mod(rand(Int), n) + 1 for _ = 1:n-2]
    G = prufer_restore(code)
    name(G, "Tree")
    return G
end

# # This is a helper function for RandomTree that converts a Prufer code
# # to a tree. No checks are done on the array fed into this function.
# function code_to_tree(code::Array{Int,1})
#     n = length(code) + 2
#     G = IntGraph(n)
#     degree = ones(Int, n)  # initially all 1s

#     #every time a vertex appears in code[], up its degree by 1
#     for c in code
#         degree[c] += 1
#     end

#     for u in code
#         for v = 1:n
#             if degree[v] == 1
#                 add!(G, u, v)
#                 degree[u] -= 1
#                 degree[v] -= 1
#                 break
#             end
#         end
#     end

#     last = findall(x -> x != 0, degree)
#     add!(G, last[1], last[2])

#     return G
# end

# Create the Cube graph with 2^n vertices
"""
`Cube(n)` creates the `n`-dimensional cube graph. This graph has `2^n`
vertices named by all possible length-`n` strings of 0s and 1s. Two
vertices are adjacent iff they differ in exactly one position.
"""
function Cube(n::Integer = 3)
    G = StringGraph()
    for u = 0:2^n-1
        for shift = 0:n-1
            v = xor((1 << shift), u)
            add!(G, string(u, base = 2, pad = n), string(v, base = 2, pad = n))
        end
    end
    name(G, "Cube($n)")

    if n < 3
        set_rot(G)
    end

    if n == 3
        F = ["000", "001", "011", "010"]
        embed(G, :tutte, outside = F)
        embed_rot(G)
    end

    return G
end

# Create the BuckyBall graph

"""
`BuckyBall()` returns the Bucky ball graph.
"""
function BuckyBall()
    G = IntGraph()
    edges = [
        (1, 3),
        (1, 49),
        (1, 60),
        (2, 4),
        (2, 10),
        (2, 59),
        (3, 4),
        (3, 37),
        (4, 18),
        (5, 7),
        (5, 9),
        (5, 13),
        (6, 8),
        (6, 10),
        (6, 17),
        (7, 8),
        (7, 21),
        (8, 22),
        (9, 10),
        (9, 57),
        (11, 12),
        (11, 13),
        (11, 21),
        (12, 28),
        (12, 48),
        (13, 14),
        (14, 47),
        (14, 55),
        (15, 16),
        (15, 17),
        (15, 22),
        (16, 26),
        (16, 42),
        (17, 18),
        (18, 41),
        (19, 20),
        (19, 21),
        (19, 27),
        (20, 22),
        (20, 25),
        (23, 24),
        (23, 32),
        (23, 35),
        (24, 26),
        (24, 39),
        (25, 26),
        (25, 31),
        (27, 28),
        (27, 31),
        (28, 30),
        (29, 30),
        (29, 32),
        (29, 36),
        (30, 45),
        (31, 32),
        (33, 35),
        (33, 40),
        (33, 51),
        (34, 36),
        (34, 46),
        (34, 52),
        (35, 36),
        (37, 38),
        (37, 41),
        (38, 40),
        (38, 53),
        (39, 40),
        (39, 42),
        (41, 42),
        (43, 44),
        (43, 47),
        (43, 56),
        (44, 46),
        (44, 54),
        (45, 46),
        (45, 48),
        (47, 48),
        (49, 50),
        (49, 53),
        (50, 54),
        (50, 58),
        (51, 52),
        (51, 53),
        (52, 54),
        (55, 56),
        (55, 57),
        (56, 58),
        (57, 59),
        (58, 60),
        (59, 60),
    ]
    for e in edges
        add!(G, e[1], e[2])
    end
    name(G, "Buckyball")
    embed(G, :tutte, outside = [15, 16, 42, 41, 18, 17])
    embed_rot(G)
    return G
end

# The Kneser graph Kneser(n,k) has C(n,k) vertices that are the
# k-element subsets of 1:n in which two vertices are adjacent if (as
# sets) they are disjoint. The Petersen graph is Kneser(5,2).

# import IterTools.subsets

"""
`Kneser(n,m)` creates the Kneser graph whose vertices are all the
`m`-element subsets of `1:n` in which two vertices are adjacent iff
they are disjoint.
"""
function Kneser(n::Int, k::Int)
    A = collect(1:n)
    vtcs = [Set(v) for v in IterTools.subsets(A, k)]
    G = SimpleGraph{Set{Int}}()

    for v in vtcs
        add!(G, v)
    end

    nn = length(vtcs)
    for i = 1:nn-1
        u = vtcs[i]
        for j = i+1:nn
            v = vtcs[j]
            if length(intersect(u, v)) == 0
                add!(G, u, v)
            end
        end
    end
    name(G, "Kneser($n,$k)")
    return G
end

"""
`Johnson(n,k)` creates the Johnson graph whose vertices
and the `k`-element subsets of `{1,2,...,n}`. Vertices `v`
and `w` are adjacent if their intersection has size `k-1`.
"""
function Johnson(n::Int, k::Int)
    if k < 0 || n < 0
        error("n,k must be nonnegative")
    end

    if k > n
        error("k must not be greater than n")
    end

    A = collect(1:n)
    vtcs = [Set{Int}(v) for v in IterTools.subsets(A, k)]
    G = SimpleGraph{Set{Int}}()

    for v in vtcs
        add!(G, v)

        for j = 1:n
            if !in(j, v)
                for x in v
                    w = deepcopy(v)
                    push!(w, j)
                    pop!(w, x)
                    add!(G, v, w)
                end
            end
        end
    end

    name(G, "Johnson($n,$k)")
    return G
end

# Create the Petersen graph.
"""
`Petersen()` returns the Petersen graph. The vertices are labeled as
the 2-element subsets of `1:5`. Wrap in `relabel` to have vertices
named `1:10`. See also: `Kneser`.
"""
function Petersen()
    G = Kneser(5, 2)
    name(G, "Petersen")
    embed(
        G,
        :tutte,
        outside = [Set([1, 4]), Set([3, 5]), Set([2, 4]), Set([1, 3]), Set([2, 5])],
    )
    embed_rot(G)
    return G
end

# Create Paley graphs

"""
`Paley(p)` creates the Paley graph with `p` vertices named
`0:(p-1)`. Here `p` must be a prime with `p%4==1`. Vertices `u` and
`v` are adjacent iff `u-v` is a quadratic residue (perfect square)
modulo `p`.
"""
function Paley(p::Int)
    if mod(p, 4) != 1 || ~isprime(p)
        error("p must be a prime congruent to 1 mod 4")
    end

    # Quadratic residues mod p
    qrlist = unique([mod(k * k, p) for k = 1:p])

    G = IntGraph()
    for u = 0:p-1
        for k in qrlist
            v = mod(u + k, p)
            add!(G, u, v)
        end
    end
    name(G, "Paley")
    return G
end

"""
`Frucht()` returns the Frucht graph: A 12-vertex, 3-regular
graph with no non-nontrivial automorphisms.
"""
function Frucht()
    G = Cycle(12)
    more_edges = [(1, 6), (2, 4), (3, 11), (5, 7), (8, 10), (9, 12)]
    add_edges!(G, more_edges)
    F = [3, 11, 10, 8, 7, 5, 4]
    embed(G, :tutte, outside = F)
    embed_rot(G)
    name(G, "Frucht")
    return G
end

# Called by RandomRegular ... one step
function RandomRegularBuilder(n::Int, d::Int)
    # caller has already checked the values of n,d are legit
    vlist = randperm(n * d)
    G = IntGraph(n * d)
    for v = 1:2:n*d
        add!(G, vlist[v], vlist[v+1])
    end

    for v = n:-1:1
        mushlist = collect(d*(v-1)+1:v*d)
        for k = d-1:-1:1
            contract!(G, mushlist[k], mushlist[k+1])
        end
    end
    return relabel(G)
end

"""
`RandomRegular(n,d)` creates a random `d`-regular graph on `n`
vertices. This can take a while especially if the arguments are
large. Call with an optional third argument to activate verbose
progress reports: `RandomRegular(n,p,true)`.
"""
function RandomRegular(n::Int, d::Int, verbose::Bool = false)
    # sanity checks
    if n < 1 || d < 1 || (n * d) % 2 == 1
        error("n,d must be positive integers and n*d even")
    end
    if verbose
        println("Trying to build ", d, "-regular graph on ", n, " vertices")
        count::Int = 0
    end

    while true
        if verbose
            count += 1
            println("Attempt ", count)
            tic()
        end
        g = RandomRegularBuilder(n, d)
        if verbose
            toc()
        end
        dlist = deg(g)
        if dlist[1] == dlist[n]
            if verbose
                println("Success")
            end
            name(g, "Regular, d=$d")
            return g
        end
        if verbose
            println("Failed; trying again")
        end
    end
end

"""
`RandomSBM(bmap,pmat)` creates a random stochastic block model random graph.
The vector `bmap` is a list of `n` positive integers giving the block number
of vertices `1:n`. The `i,j`-entry of the matrix `pmat` gives the probability
of an edge from a vertex in block `i` to a vertex in block `j`.

`RandomSBM(n,pvec,pmat)` creates such a graph with `n` vertices. The vector
`pvec` gives the probabilities that vertices fall into a given block.
"""
function RandomSBM(bmap::Vector{Int}, pmat::Array{S,2}) where {S<:Real}
    n = length(bmap)    # no of vertices
    b = maximum(bmap)   # no of blocks

    @assert minimum(bmap) > 0 "Block numbers must be positive"
    @assert pmat == pmat' "Probability matrix must be square and symmetrical"
    r = size(pmat, 1)
    @assert b <= r "Insufficient rows/cols in probability matrix"
    @assert minimum(pmat) >= 0 && maximum(pmat) <= 1 "Entries in probability matrix out of range"

    G = IntGraph(n)
    for u = 1:n-1
        bu = bmap[u]
        for v = u+1:n
            bv = bmap[v]
            p = pmat[bu, bv]
            if rand() <= p
                add!(G, u, v)
            end
        end
    end
    name(G, "Stochastic Block Model")
    return G
end

function RandomSBM(n::Int, pvec::Vector{S}, pmat::Array{T,2}) where {S<:Real,T<:Real}
    @assert minimum(pvec) >= 0 "Entries in pvec must be nonnegative"
    @assert sum(pvec) == 1 "Entries in pvec must sum to 1"
    bmap = [random_choice(pvec) for _ = 1:n]
    return RandomSBM(bmap, pmat)
end

"""
`Knight(r::Int=8,c::Int=8)` creates a Knight's Moves graph on a
`r`-by-`c` grid. That is, the vertices of this graph are the squares
of an `r`-by-`c` chess board. Two vertices are adjacent if a Knight
can go from one of these squares to the other in a single move.
"""
function Knight(r::Int = 8, c::Int = 8)
    vtcs = collect(Base.Iterators.product(1:r, 1:c))
    G = SimpleGraph{Tuple{Int64,Int64}}()
    d = Dict{Tuple{Int,Int},Vector{Float64}}()

    for v in vtcs
        add!(G, v)
        d[v] = collect(v)
    end

    for v in vtcs
        for w in vtcs
            xv = collect(v)
            xw = collect(w)
            z = sort(map(abs, xv - xw))
            if z == [1, 2]
                add!(G, v, w)
            end
        end
    end
    name(G, "Knight($r,$c)")
    embed(G, d)

    return G
end

"""
`HoffmanSingleton()` creates the Hoffman-Singleton graph. This is
a 7-regular graph whose diameter is 2 and whose girth is 5.
"""
function HoffmanSingleton()
    G = StringGraph()
    # P-pentagons
    for i = 0:4
        for j = 0:4
            jj = mod(j + 1, 5)
            v = "P$i$j"
            w = "P$i$jj"
            add!(G, v, w)
        end
    end

    # Q-pentagrams
    for i = 0:4
        for j = 0:4
            jj = mod(j + 2, 5)
            v = "Q$i$j"
            w = "Q$i$jj"
            add!(G, v, w)
        end
    end

    # Connections
    for i = 0:4
        for j = 0:4
            for k = 0:4
                v = "P$i$j"
                x = mod(i * k + j, 5)
                w = "Q$k$x"
                add!(G, v, w)
            end
        end
    end
    name(G, "Hoffman-Singleton")
    return G
end

"""
`Hoffman()` creates the Hoffman graph which is cospectral with,
but not isomorphic to, `Cube(4)`.
"""
function Hoffman()
    D = [
        1 1 1 1 0 0 0 0
        1 1 1 0 1 0 0 0
        1 0 0 1 0 1 1 0
        0 1 0 1 0 1 0 1
        0 0 1 1 0 0 1 1
        1 0 0 0 1 1 1 0
        0 1 0 0 1 1 0 1
        0 0 1 0 1 0 1 1
    ]

    A = [0*D D; D' 0*D]

    G = SimpleGraph(A)
    name(G, "Hoffman")
    return G
end

"""
`Doyle()` creates the Doyle/Holt graph.
See article on [Mathworld](http://mathworld.wolfram.com/DoyleGraph.html)
"""
function Doyle()
    T = Tuple{Int,Int}
    G = SimpleGraph{T}()
    for a = 0:8
        for b = 0:2
            v = (a, b)
            aa = mod(4a + 1, 9)
            bb = mod(b - 1, 3)
            w = (aa, bb)
            add!(G, v, w)
            aa = mod(4a - 1, 9)
            w = (aa, bb)
            add!(G, v, w)
        end
    end
    name(G, "Doyle")
    return G
end

"""
`Spindle()` returns the Moser spindle graph. This is a seven-vertex
unit distance graph with chromatic number equal to 4.
"""
function Spindle()
    G = IntGraph(7)
    edges = [
        (1, 2),
        (1, 3),
        (2, 3),
        (2, 4),
        (3, 4),
        (1, 5),
        (1, 6),
        (5, 6),
        (5, 7),
        (6, 7),
        (4, 7),
    ]
    add_edges!(G, edges)

    # compute a rotation system that's planar
    F = [1, 3, 4, 7, 5]
    embed(G, :tutte, outside = F)
    embed_rot(G)

    # but then give a unit-distance embedding
    d = Dict{Int,Vector}()
    a = sqrt(3) / 2

    pts = [0 1/2 -1/2 0; 0 a a 2a]

    theta = acos(5 / 6) / 2
    R = [cos(theta) -sin(theta); sin(theta) cos(theta)]

    p1 = R * pts
    for k = 1:4
        d[k] = p1[:, k]
    end

    p2 = R' * pts
    for k = 5:7
        d[k] = p2[:, k-3]
    end
    embed(G, d)
    SimpleGraphs.name(G, "Moser Spindle")
    return G
end

"""
    Golomb()
Create the Golomb graph. This is a unit-distance graph with chromatic number 4.
It has 10 vertices and 18 edges.
"""
function Golomb()::SimpleGraph{Int}
    G = Cycle(6)
    for v = 1:6
        add!(G, 0, v)
    end

    add!(G, 1, 7)
    add!(G, 3, 8)
    add!(G, 5, 9)
    add!(G, 7, 8)
    add!(G, 7, 9)
    add!(G, 8, 9)

    xy = Dict{Int,Vector{Float64}}()

    xy[0] = [0, 0]
    for k = 1:6
        x, y = reim(exp((k - 1) * im * 2 * π / 6))
        xy[k] = [x, y]
    end

    r = sqrt(3) / 3
    θ = acos(r / 2)

    for k = 0:2
        x, y = reim(r * exp((θ + 2 * k * π / 3) * im))
        xy[k+7] = [x, y]
    end

    embed(G, xy)
    name(G, "Golomb")

    return G
end

"""
    is_unit_distance(G, tol=1e-10)
Check to see if the embedded graph `G` is a unit-distance graph. That is, two vertices of `G`
should be distance 1 apart if and only if they are adjacent. The optional `tol` gives some tolerance
to this assessment to deal with roundoff.
"""
function is_unit_distance(G::SimpleGraph{T}, tol = 1e-10)::Bool where {T}
    xy = getxy(G)
    VV = vlist(G)
    n = NV(G)
    for i = 1:n-1
        u = VV[i]
        xyu = xy[u]
        for j = u+1:n
            v = VV[j]
            xyv = xy[v]
            unit_check = abs(norm(xyu - xyv) - 1) <= tol

            if has(G, u, v)
                if !unit_check
                    @info "$u and $v are adjacent, but not distance 1 (within tolerance $tol)"
                    return false
                end
            else
                if unit_check
                    @info "$u and $v are not adjacent, but are distance 1 (within tolerance $tol)"
                    return false
                end
            end

        end
    end
    return true

end


text:
struct BrownianMotion{F<:Number} <: Motion
    noise::F
    persistence_time::F
    persistence_length::F
    effective_diffusion::F
    friction::F
end

function BrownianMotion( noise::F, friction::F ) where F<:Number
    persistence_time = 1/friction
    persistence_length = sqrt( noise / friction )
    effective_diffusion = noise
    BrownianMotion{F}( noise, persistence_time, persistence_length, effective_diffusion, friction )
end # function

function BrownianMotion(; noise = 0.0, friction = 0.0 )
    BrownianMotion( noise, friction )
end # function

function step( motion::BrownianMotion, position::P, velocity::V, dt ) where {P,V}
    γ = motion.friction
    ξ = motion.noise
    new_velocity = velocity - V(
        γ * velocity[1] * dt - sqrt( 2ξ * γ * dt ) * randn(),
        γ * velocity[2] * dt - sqrt( 2ξ * γ * dt ) * randn(),
        )
    new_position = position + P(
        new_velocity[1] * dt,
        new_velocity[2] * dt,
        )
    return new_position, new_velocity
end # function


text:
"""
Returns true if point, pt, is contained by the polygon defined by points in poly

Pulled from:
https://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html

Based on Jordan curve theorem.

Draw a semi-infinite ray from the point (x infinite, y fixed) & count how many
edges that ray crosses (or actually, how many edges it is not parallel to). If
the ray crosses an odd number of edges when its slope is larger than the slope
of the edge, its contained. Only consider edges for which the point falls within
its y-range.

The semi-infinite ray method is inconsistent on boundaries, so include a test
for the border at the very end.
"""
function pt_in_poly(pt, poly)
  is_contained = false
  N = size(poly,1)
  # only count the distinct vertices
  if poly[1,:] == poly[N,:]
    N = N-1
  end
  j = N
  for i = 1:N
    i_vert = poly[i,:][:]
    j_vert = poly[j,:][:]
    # is point is on same side of edges endpoints in y?
    if (pt[2] > i_vert[2]) != (pt[2] > j_vert[2])
      # is the line from pt->i_vert parallel to j_vert->i_vert?
      if pt[1] < (j_vert[1]-i_vert[1])*(pt[2]-i_vert[2])/(j_vert[2]-i_vert[2]) + i_vert[1]
        is_contained = !is_contained
      end
    end
    j = i
  end

  if !is_contained
    j = N
    for i = 1:N
      i_vert = poly[i,:][:]
      j_vert = poly[j,:][:]
      on_line = pt_on_line_segment(pt, (i_vert, j_vert))
      if on_line
        return true
      end
    end
  end

  return is_contained
end

"""
If point exists along the line segment, the triangle it defines should be flat.
"""
function pt_on_line_segment(pt, line)
  epsilon = 1e-4
  a, b = line
  return -epsilon < norm(a-b) - (norm(a-pt) + norm(b-pt)) < epsilon
end

"""
Returns true if polyA contains polyB
"""
function poly_contains_poly(polyA, polyB)
  for i in 1:size(polyB, 1)
    pt = polyB[i,:][:]
    if pt_in_poly(pt, polyA)
      return true
    end
  end
  return false
end

"""
Returns true if polyA & polyB intersect
"""
function poly_intersects(polyA, polyB)
  return poly_contains_poly(polyA, polyB) || poly_contains_poly(polyB, polyA)
end

"""
Sutherland-Hodgman

* Assumes clip polygon is convex

https://en.wikipedia.org/wiki/Sutherland%E2%80%93Hodgman_algorithm
"""
function clip_polygon(subject, clip)
    output = [[subject[i,:]...] for i in 1:size(subject,1)]
    clip_edges = hcat(clip[1:end-1,:], clip[2:end,:])
    edge1 = line_to_vector(clip_edges[1,:])
    edge2 = line_to_vector(clip_edges[2,:])
    is_counterclockwise = cross2(edge1, edge2) > 0
    for i in 1:size(clip_edges,1)
        edge = clip_edges[i,:]
        input = output
        output = []
        start_pt = input[end]
        for j in 1:size(input,1)
            end_pt = input[j]
            if pt_on_left(end_pt, edge) == is_counterclockwise
                if pt_on_left(start_pt, edge) != is_counterclockwise
                    mid_pt = find_line_intersection([start_pt..., end_pt...], edge)
                    push!(output, mid_pt)
                end
                push!(output, end_pt)
            else
                if pt_on_left(start_pt, edge) == is_counterclockwise
                    mid_pt = find_line_intersection([start_pt..., end_pt...], edge)
                    push!(output, mid_pt)
                end
            end
            start_pt = end_pt
        end
    end
    if length(output) > 1
      push!(output, output[1])
    end
    return hcat(output...)'
end

"""
Determine if point is to the left of an edge

"""
function pt_on_left(point, line)
    pt_vector = line_to_vector([line[1:2], point...])
    line_vector = line_to_vector(line)
    return cross2(line_vector, pt_vector) > 0
end

"""
Find intersecting point of two line segments

Inputs:
    * a: first line segment (array of two points)
    * b: second line segment (array of two points)

Outputs:
    * if the line segments intersect: point
    * if the line segments don't intersect: nothing

http://stackoverflow.com/questions/563198/how-do-you-detect-where-two-line-segments-intersect
"""
function find_line_segment_intersection(a, b)
    p = a[1:2]
    r = line_to_vector(a)
    q = b[1:2]
    s = line_to_vector(b)
    if cross2(r, s) == 0
        if cross2(q-p, r) != 0
            # parallel, non-intersecting
            return nothing
        else
            # collinear
            t0 = (q-p)' * r/(r'*r)
            t1 = t0 + s' * r/(r'*r)
            if (s'*r)[1] < 0
                return p + t1[1]*r
            else
                return p + t0[1]*r
            end
        end
    else
        t = cross2(q-p, s / cross2(r, s))
        u = cross2(q-p, r / cross2(r, s))
        if (0 <= t <=1) & (0 <= u <= 1)
            # intersect @ one point
            return p + t*r
        else
            # do not intersect
            return nothing
        end
    end
end

"""
Find intersection of two lines - no protection
"""
function find_line_intersection(a, b)
    p = a[1:2]
    r = line_to_vector(a)
    q = b[1:2]
    s = line_to_vector(b)
    if cross2(r, s) == 0
      # collinear
      t0 = (q-p)' * r/(r'*r)
      t1 = t0 + s' * r/(r'*r)
      if (s'*r)[1] < 0
          return p + t1[1]*r
      else
          return p + t0[1]*r
      end
    else
      t = cross2(q-p, s / cross2(r, s))
      # u = cross2(q-p, r / cross2(r, s))
      return p + t*r
    end
end

"""
Convert 4-element line segment to 2-element vector
"""
function line_to_vector(a)
    return [a[3]-a[1], a[4]-a[2]]
end

"""
Given two 2-element vectors, return their determinant
"""
function cross2(a, b)
    return a[1]*b[2] - a[2]*b[1]
end

"""
Given an Nx2 list of convex polygon vertices, calculate polygon area
"""
function poly_area(a)
  N = size(a,1)
  if a[1,:] == a[N,:]
    N -= 1
  end
  area = 0
  for i = 1:N
    area += abs(cross2(a[i,:], a[N,:]))
  end
  return area / 2
end


text:
using RobustDMD, Random, Printf

# Generate DMD Synthetic Data
# ------------------------------------------------------------------------------
# dimensions
m = 500;    # temporal dimension
n = 1000;    # spatial dimension
k = 3;      # number of modes
T = Float64;
sigma = 1e-3;
# generate data
X, t, at, Bt = simDMD(m, n, k, T; seed=123456);
Xdat = X + sigma * Random.randn!(copy(X))
# specify the loss
function lossFunc(r)
	return sum(abs2, r)
end
#
function lossGrad(r)
	conj!(r)
end
#
# Create object
# ------------------------------------------------------------------------------
params = DMDParams(k, Xdat, t, lossFunc, lossGrad, inner_directl2=true);
#
# Initial Guess
# ------------------------------------------------------------------------------
# create initial guess using exact dmd
a0, B0 = dmdexactestimate(m, n, k, Xdat, t);
#
# Create Object
# ------------------------------------------------------------------------------
params = DMDParams(k, Xdat, t, lossFunc, lossGrad);
copyto!(params.a, a0);
copyto!(params.B, B0);
# Random.randn!(params.ar);
#
# Apply Solver
# ------------------------------------------------------------------------------
tau = 10;
eta = 1e1;
dms = 500;
itm = 1000;
tol = 1e-6;
ptf = 10;

function prox_stab(ar)
    # project for numerical stability
    k = length(ar) >> 1;
    @inbounds for i = 1:k
        ar[i << 1 - 1] = min(2.0, ar[i << 1 - 1]);
    end
end

opts = DMD_SPG_Options(tau, eta, dms=dms, itm=itm, tol=tol, ptf=ptf,
                       prox=prox_stab);

@time obj_his, err_his = solveDMD_withSPG(params, opts);

@printf "%7.4e relative residual after SPG\n" sqrt(dmdobjective(params) / lossFunc(X))

params2 = DMDParams(k, Xdat, t, lossFunc, lossGrad);
copyto!(params2.a,at)
@printf "%7.4e relative residual with true evals\n" sqrt(dmdobjective(params2) / lossFunc(X))
copyto!(params2.a,a0)
@printf "%7.4e relative residual with initial guess\n" sqrt(dmdobjective(params2) / lossFunc(X))

@printf "%7.4e abs err evals after SPG\n" besterrperm(params.a, at)
@printf "%7.4e abs err evals initial guess\n" besterrperm(a0, at)


text:
# usage example for package developer:
#
#     import HistogramAdjustmentAPI: AbstractHistogramAdjustmentAlgorithm,
#                             adjust_histogram, adjust_histogram!

"""
    AbstractHistogramAdjustmentAlgorithm <: AbstractImageFilter

The root type for `ImageContrastAdjustment` package.

Any concrete histogram adjustment algorithm shall subtype it to support
[`adjust_histogram`](@ref) and [`adjust_histogram!`](@ref) APIs.

# Examples

All histogram adjustment algorithms in ImageContrastAdjustment are called in the
following pattern:

```julia
# first generate an algorithm instance
f = LinearStretching()

# then pass the algorithm to `adjust_histogram`
img_adjusted = adjust_histogram(img, f)

# or use in-place version `adjust_histogram!`
img_adjusted = similar(img)
adjust_histogram!(img_adjusted, img, f)
```

Some algorithms also receive additional information as an argument,
e.g., `nbins` of `Equalization`.

```julia
# you can explicit specify the parameters
f = Equalization(nbins = 32)
```

For more examples, please check [`adjust_histogram`](@ref),
[`adjust_histogram!`](@ref) and concrete algorithms.
"""
abstract type AbstractHistogramAdjustmentAlgorithm <: AbstractImageFilter end

adjust_histogram!(out::Union{GenericGrayImage, AbstractArray{<:Color3}},
          img,
          f::AbstractHistogramAdjustmentAlgorithm,
          args...; kwargs...) =
    f(out, img, args...; kwargs...)

# TODO: Relax this to all color types
function adjust_histogram!(img::Union{GenericGrayImage, AbstractArray{<:Color3}},
                   f::AbstractHistogramAdjustmentAlgorithm,
                   args...; kwargs...)
    tmp = copy(img)
    f(img, tmp, args...; kwargs...)
    return img
end

function adjust_histogram(::Type{T},
                  img,
                  f::AbstractHistogramAdjustmentAlgorithm,
                  args...; kwargs...) where T
    out = similar(Array{T}, axes(img))
    adjust_histogram!(out, img, f, args...; kwargs...)
    return out
end

adjust_histogram(img::AbstractArray{T},
                 f::AbstractHistogramAdjustmentAlgorithm,
                 args...; kwargs...) where T <: Colorant =
         adjust_histogram(T, img, f, args...; kwargs...)

# Do not promote Number to Gray{<:Number}
adjust_histogram(img::AbstractArray{T},
                 f::AbstractHistogramAdjustmentAlgorithm,
                 args...; kwargs...) where T <: Number =
        adjust_histogram(T, img, f, args...; kwargs...)

# Handle instance where the input is a sequence of images.
adjust_histogram!(out_sequence::Vector{T},
          img_sequence,
          f::AbstractHistogramAdjustmentAlgorithm,
          args...; kwargs...) where T <: Union{GenericGrayImage, AbstractArray{<:Color3}} =
    f(out_sequence, img_sequence, args...; kwargs...)

# TODO: Relax this to all color types
function adjust_histogram!(img_sequence::Vector{T},
                   f::AbstractHistogramAdjustmentAlgorithm,
                   args...; kwargs...) where T <: Union{GenericGrayImage, AbstractArray{<:Color3}}
    tmp = copy(img_sequence)
    f(img_sequence, tmp, args...; kwargs...)
    return img_sequence
end

function adjust_histogram(::Type{T},
                  img_sequence::Vector{<:AbstractArray},
                  f::AbstractHistogramAdjustmentAlgorithm,
                  args...; kwargs...) where T
    N  = length(img_sequence)
    out_sequence = [similar(Array{T}, axes(img_sequence[n])) for n = 1:N]
    adjust_histogram!(out_sequence, img_sequence, f, args...; kwargs...)
    return out_sequence
end

adjust_histogram(img_sequence::Vector{<:AbstractArray{T}},
                 f::AbstractHistogramAdjustmentAlgorithm,
                 args...; kwargs...) where T <: Colorant =
         adjust_histogram(T, img_sequence, f, args...; kwargs...)

# Do not promote Number to Gray{<:Number}
adjust_histogram(img_sequence::Vector{<:AbstractArray{T}},
                 f::AbstractHistogramAdjustmentAlgorithm,
                 args...; kwargs...) where T <: Number =
        adjust_histogram(T, img_sequence, f, args...; kwargs...)

### Docstrings

"""
    adjust_histogram!([out,] img, f::AbstractHistogramAdjustmentAlgorithm, args...; kwargs...)

Adjust histogram of `img` using algorithm `f`.

# Output

If `out` is specified, it will be changed in place. Otherwise `img` will be changed in place.

# Examples

Just simply pass an algorithm to `adjust_histogram!`:

```julia
img_adjusted = similar(img)
adjust_histogram!(img_adjusted, img, f)
```

For cases you just want to change `img` in place, you don't necessarily need to manually
allocate `img_adjusted`; just use the convenient method:

```julia
adjust_histogram!(img, f)
```

See also: [`adjust_histogram`](@ref)
"""
adjust_histogram!

"""
    adjust_histogram([T::Type,] img, f::AbstractHistogramAdjustmentAlgorithm, args...; kwargs...)

Adjust histogram of `img` using algorithm `f`.

# Output

The return image `img_adjusted` is an `Array{T}`.

If `T` is not specified, then it's inferred.
# Examples

Just simply pass the input image and algorithm to `adjust_histogram`

```julia
img_adjusted = adjust_histogram(img, f)
```

This reads as "`adjust_histogram` of image `img` using algorithm `f`".

You can also explicitly specify the return type:

```julia
img_adjusted_float32 = adjust_histogram(Gray{Float32}, img, f)
```

See also [`adjust_histogram!`](@ref) for in-place histogram adjustment.
"""
adjust_histogram


text:
module SimpleWeightedGraphs

using LightGraphs
using LinearAlgebra
using Markdown
using SparseArrays

import Base:
    convert, eltype, show, ==, Pair, Tuple, copy, length, issubset, zero

import LightGraphs:
    _NI, AbstractGraph, AbstractEdge, AbstractEdgeIter,
    src, dst, edgetype, nv, ne, vertices, edges, is_directed,
    add_vertex!, add_edge!, rem_vertex!, rem_edge!,
    has_vertex, has_edge, inneighbors, outneighbors,
    indegree, outdegree, degree, has_self_loops, num_self_loops,

    add_vertices!, adjacency_matrix, weights, connected_components, cartesian_product,

    AbstractGraphFormat, loadgraph, loadgraphs, savegraph,
    pagerank, induced_subgraph

export
    AbstractSimpleWeightedGraph,
    AbstractSimpleWeightedEdge,
    SimpleWeightedEdge,
    SimpleWeightedGraph,
    SimpleWeightedGraphEdge,
    SimpleWeightedDiGraph,
    SimpleWeightedDiGraphEdge,
    weight,
    weighttype,
    get_weight,
    WGraph,
    WDiGraph,
    SWGFormat

include("simpleweightededge.jl")

"""
    AbstractSimpleWeightedGraph

An abstract type representing a simple graph structure.
AbstractSimpleWeightedGraphs must have the following elements:
- weightmx::AbstractSparseMatrix{Real}
"""
abstract type AbstractSimpleWeightedGraph{T<:Integer,U<:Real} <: AbstractGraph{T} end

function show(io::IO, g::AbstractSimpleWeightedGraph{T, U}) where T where U
    dir = is_directed(g) ? "directed" : "undirected"
    print(io, "{$(nv(g)), $(ne(g))} $dir simple $T graph with $U weights")
end

# conversion to SparseMatrixCSC
convert(::Type{SparseMatrixCSC{T, U}}, g::AbstractSimpleWeightedGraph) where T<:Real where U<:Integer = SparseMatrixCSC{T, U}(g.weights)

### INTERFACE

nv(g::AbstractSimpleWeightedGraph{T, U}) where T where U = T(size(g.weights, 1))
vertices(g::AbstractSimpleWeightedGraph{T, U}) where T where U = one(T):nv(g)
eltype(x::AbstractSimpleWeightedGraph{T, U}) where T where U = T
weighttype(x::AbstractSimpleWeightedGraph{T, U}) where T where U = U

has_edge(g::AbstractSimpleWeightedGraph{T, U}, e::AbstractSimpleWeightedEdge) where T where U =
    g.weights[dst(e), src(e)] != zero(U)

# handles single-argument edge constructors such as pairs and tuples
has_edge(g::AbstractSimpleWeightedGraph{T, U}, x) where T where U = has_edge(g, edgetype(g)(x))
add_edge!(g::AbstractSimpleWeightedGraph{T, U}, x) where T where U = add_edge!(g, edgetype(g)(x))

# handles two-argument edge constructors like src,dst
has_edge(g::AbstractSimpleWeightedGraph, x, y) = has_edge(g, edgetype(g)(x, y, 0))
add_edge!(g::AbstractSimpleWeightedGraph, x, y) = add_edge!(g, edgetype(g)(x, y, 1))
add_edge!(g::AbstractSimpleWeightedGraph, x, y, z) = add_edge!(g, edgetype(g)(x, y, z))

function issubset(g::T, h::T) where T<:AbstractSimpleWeightedGraph
    (gmin, gmax) = extrema(vertices(g))
    (hmin, hmax) = extrema(vertices(h))
    return (hmin <= gmin <= gmax <= hmax) && issubset(edges(g), edges(h))
end

has_vertex(g::AbstractSimpleWeightedGraph, v::Integer) = v in vertices(g)

function rem_edge!(g::AbstractSimpleWeightedGraph{T, U}, u::Integer, v::Integer) where T where U
    rem_edge!(g, edgetype(g)(T(u), T(v), one(U)))
end

@doc_str """
    rem_vertex!(g::AbstractSimpleWeightedGraph, v)

Remove the vertex `v` from graph `g`. Return false if removal fails
(e.g., if vertex is not in the graph); true otherwise.

### Implementation Notes
This operation has to be performed carefully if one keeps external
data structures indexed by edges or vertices in the graph, since
internally the removal results in all vertices with indices greater than `v`
being shifted down one.
"""
function rem_vertex!(g::AbstractSimpleWeightedGraph, v::Integer)
    v in vertices(g) || return false
    n = nv(g)

    newweights = g.weights[1:nv(g) .!= v, :]
    newweights = newweights[:, 1:nv(g) .!= v]

    g.weights = newweights
    return true
end

function outneighbors(g::AbstractSimpleWeightedGraph, v::Integer)
    mat = g.weights
    return mat.rowval[mat.colptr[v]:mat.colptr[v+1]-1]
end

get_weight(g::AbstractSimpleWeightedGraph, u::Integer, v::Integer) = g.weights[v, u]

zero(g::T) where T<:AbstractSimpleWeightedGraph = T()

# TODO: manipulte SparseMatrixCSC directly
add_vertex!(g::AbstractSimpleWeightedGraph) = add_vertices!(g, 1)

copy(g::T) where T <: AbstractSimpleWeightedGraph =  T(copy(g.weights))

const SimpleWeightedGraphEdge = SimpleWeightedEdge
const SimpleWeightedDiGraphEdge = SimpleWeightedEdge
include("simpleweighteddigraph.jl")
include("simpleweightedgraph.jl")
include("overrides.jl")
include("persistence.jl")

const WGraph = SimpleWeightedGraph
const WDiGraph = SimpleWeightedDiGraph

SimpleWeightedDiGraph(g::SimpleWeightedGraph) = SimpleWeightedDiGraph(g.weights)
SimpleWeightedDiGraph{T,U}(g::SimpleWeightedGraph) where T<:Integer where U<:Real =
    SimpleWeightedDiGraph(SparseMatrixCSC{U, T}(g.weights))

SimpleWeightedGraph(g::SimpleWeightedDiGraph) = SimpleWeightedGraph(g.weights .+ g.weights')
SimpleWeightedGraph{T,U}(g::SimpleWeightedDiGraph) where T<:Integer where U<:Real =
    SimpleWeightedGraph(SparseMatrixCSC{U, T}(g.weights .+ g.weights'))

end # module


text:
# Copyright 2020 Google LLC

# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at

#     https://www.apache.org/licenses/LICENSE-2.0

# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

struct SetInstruction
    addr::Int
    value::Int
end

struct MaskInstrunction
    newmask::AbstractString
end

Instruction = Union{SetInstruction, MaskInstrunction}

mutable struct ComputerState
    mem::Dict{Int, Int}
    mask::AbstractString
end

function applymask(mask, n::Number)::Number
    result = 0
    rev_mask = reverse(mask)
    for i ∈ reverse(1:36)
        bit = nothing
        if rev_mask[i] == 'X'
            bit = (n ÷ 2^(i-1)) % 2
        elseif rev_mask[i] == '0'
            bit = 0
        elseif rev_mask[i] == '1'
            bit = 1
        else
            throw(ErrorException("Unexpected mask character"))
        end
        result += bit * 2^(i-1)
    end
    return result
end

function parseinstruction(line::AbstractString)::Instruction
    if startswith(line, "mask = ")
        return MaskInstrunction(split(line)[3])
    else
        m = match(r"^mem\[([0-9]+)\] = ([0-9]+)$", line)
        if m === nothing
            throw(ErrorException("Unexpected instruction: " * line))
        end
        addr = parse(Int, m.captures[1])
        value = parse(Int, m.captures[2])
        return SetInstruction(addr, value)
    end
end

function execute!(state::ComputerState, instr::Instruction)
    if typeof(instr) == MaskInstrunction
        state.mask = instr.newmask
    elseif typeof(instr) == SetInstruction
        state.mem[instr.addr] = applymask(state.mask, instr.value)
    end
end

function part1(filename = "day14input.txt")
    instructions = parseinstruction.(readlines(filename))
    mask = 'X'^36
    mem = Dict{Int, Int}()
    state = ComputerState(mem, mask)
    for instr ∈ instructions
        execute!(state, instr)
    end
    return sum(values(mem))
end

function applyfloatingtozero(mask, n::Number)::Number
    result = 0
    rev_mask = reverse(mask)
    for i ∈ reverse(1:36)
        bit = nothing
        if rev_mask[i] == '0'
            bit = (n ÷ 2^(i-1)) % 2
        elseif rev_mask[i] == 'x'
            bit = 0
        elseif rev_mask[i] == '1'
            bit = 1
        else
            throw(ErrorException("Unexpected mask character"))
        end
        result += bit * 2^(i-1)
    end
    return result
end

function getmaskedaddresses(mask, addr::Number)
    if contains(mask, 'X')
        ind = findfirst('X', mask)
        place = 2^(length(mask) - ind)
        newmask = mask[1:ind-1] * "x" * mask[ind+1:end]
        addrs = getmaskedaddresses(newmask, addr)
        to_add = Int[]
        for a ∈ addrs
            push!(to_add, place + a)
        end
        union!(addrs, to_add)
    else
        res = applyfloatingtozero(mask, addr)
        addrs = Set{Int}()
        push!(addrs, res)
        return addrs
    end
end

function executev2!(state::ComputerState, instr::Instruction)
    if typeof(instr) == MaskInstrunction
        state.mask = instr.newmask
    elseif typeof(instr) == SetInstruction
        addrs = getmaskedaddresses(state.mask, instr.addr)
        for addr ∈ addrs
            state.mem[addr] = instr.value
        end
    end
end

function part2(filename = "day14input.txt")
    instructions = parseinstruction.(readlines(filename))
    mask = 'X'^36
    mem = Dict{Int, Int}()
    state = ComputerState(mem, mask)
    for instr ∈ instructions
        executev2!(state, instr)
    end
    return sum(values(mem))
end


text:
reference_solver = GLPKSolverLP()

regularizers = [DontRegularize(),
                RegularizedDecomposition(penaltyterm = Linearized()),
                TrustRegion(),
                LevelSet(penaltyterm = InfNorm(), projectionsolver = reference_solver)]

aggregators = [DontAggregate(),
               PartialAggregate(2),
               Aggregate(),
               DynamicAggregate(2, SelectUniform(2))]

consolidators = [Consolidate(), DontConsolidate()]

@testset "Structured Solvers" begin
    @info "Running L-shaped tests..."
    @testset "L-shaped: simple problems" begin
        @testset "$(solverstr(ls)): $name" for ls in [LShapedSolver(reference_solver,
                                                                    crash = Crash.EVP(),
                                                                    regularize = regularizer,
                                                                    aggregate = aggregator,
                                                                    consolidate = consolidator,
                                                                    log = false)
                                                      for regularizer in regularizers, aggregator in aggregators, consolidator in consolidators], (sp,res,name) in problems
            tol = 1e-5
            optimize!(sp, solver=reference_solver)
            x̄ = optimal_decision(sp)
            Q̄ = optimal_value(sp)
            if name == "Infeasible"
                with_logger(NullLogger()) do
                    @test optimize!(infeasible, solver=ls) == :Infeasible
                    add_params!(ls, feasibility_cuts = true)
                end
            end
            @test optimize!(sp, solver=ls) == :Optimal
            @test isapprox(optimal_value(sp), Q̄, rtol = tol)
            @test isapprox(optimal_decision(sp), x̄, rtol = sqrt(tol))
        end
    end
    @info "Running progressive-hedging tests..."
    @testset "Progressive-hedging: simple problems" begin
        @testset "Progressive-hedging: $name" for (sp,res,name) in problems
            tol = 1e-2
            ph = ProgressiveHedgingSolver(reference_solver, penaltyterm = Linearized(nbreakpoints=30), log = false)
            optimize!(sp, solver=reference_solver)
            x̄ = optimal_decision(sp)
            Q̄ = optimal_value(sp)
            @test optimize!(sp, solver=ph) == :Optimal
            @test isapprox(optimal_value(sp), Q̄, rtol = tol)
            @test isapprox(optimal_decision(sp), x̄, rtol = sqrt(tol))
        end
    end
end