text:
using DelimitedFiles

@sync begin
  # 2D, flux of Shima et al.
  @async run(`$(Base.julia_cmd()) --check-bounds=no --threads=1 -e 'include("run_benchmarks.jl"); run_benchmarks(3:13, 2, flux_shima_etal_notinlined, flux_shima_etal)'`)
  @async run(`$(Base.julia_cmd()) --check-bounds=no --threads=1 -e 'include("run_benchmarks.jl"); run_benchmarks(14:15, 2, flux_shima_etal_notinlined, flux_shima_etal)'`)

  # 3D, flux of Shima et al.
  @async run(`$(Base.julia_cmd()) --check-bounds=no --threads=1 -e 'include("run_benchmarks.jl"); run_benchmarks(3:13, 3, flux_shima_etal_notinlined, flux_shima_etal)'`)
  @async run(`$(Base.julia_cmd()) --check-bounds=no --threads=1 -e 'include("run_benchmarks.jl"); run_benchmarks(14:15, 3, flux_shima_etal_notinlined, flux_shima_etal)'`)

  # 2D, flux of Ranocha
  @async run(`$(Base.julia_cmd()) --check-bounds=no --threads=1 -e 'include("run_benchmarks.jl"); run_benchmarks(3:13, 2, flux_ranocha_notinlined, flux_ranocha)'`)
  @async run(`$(Base.julia_cmd()) --check-bounds=no --threads=1 -e 'include("run_benchmarks.jl"); run_benchmarks(14:15, 2, flux_ranocha_notinlined, flux_ranocha)'`)

  # 3D, flux of Ranocha
  @async run(`$(Base.julia_cmd()) --check-bounds=no --threads=1 -e 'include("run_benchmarks.jl"); run_benchmarks(3:13, 3, flux_ranocha_notinlined, flux_ranocha)'`)
  @async run(`$(Base.julia_cmd()) --check-bounds=no --threads=1 -e 'include("run_benchmarks.jl"); run_benchmarks(14:15, 3, flux_ranocha_notinlined, flux_ranocha)'`)
end


text:
function axislines!(scene, rect, spinewidth, topspinevisible, rightspinevisible,
    leftspinevisible, bottomspinevisible, topspinecolor, leftspinecolor,
    rightspinecolor, bottomspinecolor)

    bottomline = lift(rect, spinewidth) do r, sw
        y = bottom(r) - 0.5f0 * sw
        p1 = Point2(left(r) - sw, y)
        p2 = Point2(right(r) + sw, y)
        [p1, p2]
    end

    leftline = lift(rect, spinewidth) do r, sw
        x = left(r) - 0.5f0 * sw
        p1 = Point2(x, bottom(r) - sw)
        p2 = Point2(x, top(r) + sw)
        [p1, p2]
    end

    topline = lift(rect, spinewidth) do r, sw
        y = top(r) + 0.5f0 * sw
        p1 = Point2(left(r) - sw, y)
        p2 = Point2(right(r) + sw, y)
        [p1, p2]
    end

    rightline = lift(rect, spinewidth) do r, sw
        x = right(r) + 0.5f0 * sw
        p1 = Point2(x, bottom(r) - sw)
        p2 = Point2(x, top(r) + sw)
        [p1, p2]
    end

    lines!(scene, bottomline, linewidth = spinewidth, show_axis = false,
        visible = bottomspinevisible, color = bottomspinecolor)
    lines!(scene, leftline, linewidth = spinewidth, show_axis = false,
        visible = leftspinevisible, color = leftspinecolor)
    lines!(scene, rightline, linewidth = spinewidth, show_axis = false,
        visible = rightspinevisible, color = rightspinecolor)
    lines!(scene, topline, linewidth = spinewidth, show_axis = false,
        visible = topspinevisible, color = topspinecolor)
end

function interleave_vectors(vec1::Vector{T}, vec2::Vector{T}) where T
    n = length(vec1)
    @assert n == length(vec2)

    vec = Vector{T}(undef, 2 * n)
    @inbounds for i in 1:n
        k = 2(i - 1)
        vec[k + 1] = vec1[i]
        vec[k + 2] = vec2[i]
    end
    vec
end

function shrinkbymargin(rect, margin)
    IRect((rect.origin .+ margin), (rect.widths .- 2 .* margin))
end

function limits(r::Rect{N, T}) where {N, T}
    ows = r.origin .+ r.widths
    ntuple(i -> (r.origin[i], ows[i]), N)
    # tuple(zip(r.origin, ows)...)
end

function limits(r::Rect{N, T}, dim::Int) where {N, T}
    o = r.origin[dim]
    w = r.widths[dim]
    (o, o + w)
end

xlimits(r::Rect{2}) = limits(r, 1)
ylimits(r::Rect{2}) = limits(r, 2)

ext: jl


text:
# Basic way of constructing Neighborhoods.

struct KDECalibrator{K} <: LinearEstimator where K<:ContinuousUnivariateDistribution
    kernel::K
    Q::BinnedCalibrator
    m::Int64
end

function KDECalibrator(::Type{T}, m::Int64, marginal_grid;
    h= default_bandwidth(T, m)) where T<:ContinuousUnivariateDistribution
    kernel = T(h)
    Qc = BinnedCalibrator(marginal_grid, pdf.(kernel, marginal_grid))
    KDECalibrator(kernel, Qc, m)
end

# TODO: Make more general maybe
inference_target(a::KDECalibrator) = MarginalDensityTarget(0.0)

struct SincKernel <: ContinuousUnivariateDistribution
    h::Float64 #Bandwidth
end

pretty_label(kd::KDECalibrator{SincKernel})  = "Sinc Kernel"

default_bandwidth(a::Type{SincKernel}, m) = 1/sqrt(log(m))

cf(a::SincKernel, t) = one(Float64)*(-1/a.h <= t <= 1/a.h)

function pdf(a::SincKernel, t)
    if t==zero(Float64)
        return(one(Float64)/pi/a.h)
    else
        return(sin(t/a.h)/pi/t)
    end
end

struct DeLaValleePoussinKernel <: ContinuousUnivariateDistribution
    h::Float64 #Bandwidth
end

pretty_label(kd::KDECalibrator{DeLaValleePoussinKernel})  = "DLVP Kernel"

default_bandwidth(a::Type{DeLaValleePoussinKernel}, m) = 1/sqrt(log(m))

function cf(a::DeLaValleePoussinKernel, t)
    if abs(t * a.h) <= 1
        return(one(Float64))
    elseif abs(t * a.h) <= 2
        return(2*one(Float64) - abs(t * a.h))
    else
        return(zero(Float64))
    end
end

function pdf(a::DeLaValleePoussinKernel, t)
    if t==zero(Float64)
        return(3*one(Float64)/2/pi/a.h)
    else
        return(a.h*(cos(t/a.h)-cos(2*t/a.h))/pi/t^2)
    end
end

function sinc_kde(Xs, marginal_grid, ::Type{T};
     ws=KernelDensity.UniformWeights(length(Xs))) where T<:ContinuousUnivariateDistribution
    m = length(Xs)
    h = default_bandwidth(T,m)
    ker = T(h)
    (grid_min, grid_max) = extrema(marginal_grid)
    # hack for now to avoid weird floating point tricks wherein the conversion from
    # vector to range fails.. Should switch to ranges everywhere though
    marginal_h = marginal_grid[2] - marginal_grid[1] - 1e-10
    f_marginal = kde(Xs, ws,
         grid_min:marginal_h:grid_max, ker);
    f_marginal
end

function sinc_kde(Xs, marginal_grid; kwargs...)
    sinc_kde(Xs, marginal_grid, DeLaValleePoussinKernel; kwargs... )
end

mutable struct BinnedMarginalDensityNeighborhood
    f::BinnedMarginalDensity
    f_kde::UnivariateKDE
    C_bias::Float64
    C_std::Float64
    η_infl::Float64
end

function BinnedMarginalDensityNeighborhood(f, f_kde, C_bias,C_std)
    BinnedMarginalDensityNeighborhood(f, f_kde, C_bias, C_std, 0.01)
end
# should be a fit function
function fit(::Type{BinnedMarginalDensityNeighborhood}, Xs,
                marginal_grid::Vector{Float64}, ::Type{T};
                nboot=101) where T<:ContinuousUnivariateDistribution

    marginal_h = marginal_grid[2] - marginal_grid[1]
    f_kde = sinc_kde(Xs, marginal_grid, T)

    m = length(Xs)
    C_stds = Vector{Float64}(undef, nboot)

    for k =1:nboot
        # Poisson bootstrap to estimate certainty band
        Z_pois = rand(Poisson(1), m)
        ws =  Weights(Z_pois/sum(Z_pois))
        f_kde_pois =  sinc_kde(Xs, marginal_grid, T; ws=ws)
        C_stds[k] = maximum(abs.(f_kde.density .- f_kde_pois.density))
    end

    C_std = median(C_stds)

    f = BinnedMarginalDensity(f_kde.density * marginal_h, marginal_grid, marginal_h)

    BinnedMarginalDensityNeighborhood(f, f_kde, 0.0, C_std)
end

function fit(::Type{BinnedMarginalDensityNeighborhood}, Xs,
       marginal_grid::Vector{Float64};
       kwargs...)

       fit(BinnedMarginalDensityNeighborhood, Xs,
              marginal_grid, DeLaValleePoussinKernel;
              kwargs...)
end

function fit(::Type{BinnedMarginalDensityNeighborhood}, Xs,
       ds::MixingNormalConvolutionProblem, ::Type{T};
       kwargs...) where T<:ContinuousUnivariateDistribution

    m = length(Xs)
    marginal_grid = ds.marginal_grid

    f = fit(BinnedMarginalDensityNeighborhood, Xs, marginal_grid, T; kwargs...)

    #cb = ComteButucea(MarginalDensityTarget(0.0), m, marginal_grid)
    cb = KDECalibrator(T, m, marginal_grid)

    max_bias = check_bias(cb, ds; maximization=true)
    f.C_bias = max_bias
    f
end

function fit(::Type{BinnedMarginalDensityNeighborhood}, Xs,
       ds::MixingNormalConvolutionProblem;
       kwargs...)

       fit(BinnedMarginalDensityNeighborhood, Xs,
              ds::MixingNormalConvolutionProblem, DeLaValleePoussinKernel;
              kwargs...)
end

function check_bias(cb::Union{ComteButucea, KDECalibrator},
                  ds::MixingNormalConvolutionProblem;
                  maximization=true)

        f = BinnedMarginalDensity([],[],0)

        check_bias(cb.Q, ds, f, cb.m, inference_target(cb);
                         C=Inf, maximization=maximization)
end


text:
program(pgm::Type{Val{:linear}}, t0::Float64, a::Float64=1.0) = k -> t0 / (1+a*k)
program(pgm::Type{Val{:quadratic}}, t0::Float64, a::Float64=1.0) = k -> t0 / (1+a*k^2)
program(pgm::Type{Val{:exp}}, t0::Float64, a::Float64=0.85) = k -> t0*a^k
program(pgm::Type{Val{:log}}, t0::Float64, a::Float64=1.1) = k -> t0 / (1+a*log(1+k))

program(pgm::Symbol, t0::Float64) = program(Val{pgm}, t0)
program(pgm::Symbol, t0::Float64, a::Float64) = program(Val{pgm}, t0, a)

doc"""
Cooling Scheduler and Temperature Management
"""
type CoolingScheduler
    T::Float64
    cycle::Int64
    process::Function
    cooling::Function
    isfrozen::Function

    function CoolingScheduler(pgm::Type{Val{Symbol}}, tmax=500::Float64, tmin=1::Float64, alpha=nothing)
        this = new()
        (this.T, this.cycle) = (tmax, 1)
        this.process = (alpha == nothing)? program(pgm, tmax) : program(pgm, tmax, this.alpha)
        this.isfrozen = () -> this.T <= tmin

        function cooling()
            this.T = this.process(this.cycle)
            cycle += 1
        end
        this.cooling = cooling

        return this
    end
end

CoolingScheduler(pgm::Symbol) = CoolingScheduler(Val{pgm})
CoolingScheduler(pgm::Symbol, tmax::Float64) = CoolingScheduler(Val{pgm}, tmax)
CoolingScheduler(pgm::Symbol, tmax::Float64, tmin::Float64) = CoolingScheduler(Val{pgm}, tmax, tmin)
CoolingScheduler(pgm::Symbol, tmax::Float64, tmin::Float64, alpha::Float64) = CoolingScheduler(Val{pgm}, tmax, tmin, alpha)


text:
abstract Solver
abstract SolverParams

max_forks_repo_name: haroldsoh/RecurrentNN


text:
abstract type AbstractRule end

"""
    Rule{L}

The struct for identifying different rules.

Rule for `ZXDiagram`s:
* `Rule{:f}()`: rule f
* `Rule{:h}()`: rule h
* `Rule{:i1}()`: rule i1
* `Rule{:i2}()`: rule i2
* `Rule{:pi}()`: rule π
* `Rule{:c}()`: rule c

Rule for `ZXGraph`s:
* `Rule{:lc}()`: local complementary rule
* `Rule{:p1}()`: pivoting rule
* `Rule{:pab}()`: rule for removing Pauli spiders adjancent to boundary spiders
* `Rule{:p2}()`: rule p2
* `Rule{:p3}()`: rule p3
* `Rule{:id}()`: rule id
* `Rule{:gf}()`: gadget fushion rule
"""
struct Rule{L} <: AbstractRule end
Rule(r::Symbol) = Rule{r}()

"""
    Match{T<:Integer}

A struct for saving matched vertices.
"""
struct Match{T<:Integer}
    vertices::Vector{T}
end

"""
    match(r, zxd)

Returns all matched vertices, which will be store in sturct `Match`, for rule `r`
in a ZX-diagram `zxd`.
"""
Base.match(::AbstractRule, zxd::AbstractZXDiagram{T, P}) where {T, P} = Match{T}[]

function Base.match(::Rule{:f}, zxd::ZXDiagram{T, P}) where {T, P}
    matches = Match{T}[]
    for v1 in spiders(zxd)
        if spider_type(zxd, v1) == SpiderType.Z || spider_type(zxd, v1) == SpiderType.X
            for v2 in neighbors(zxd, v1)
                if spider_type(zxd, v1) == spider_type(zxd, v2) && v2 >= v1
                    push!(matches, Match{T}([v1, v2]))
                end
            end
        end
    end
    return matches
end

function Base.match(::Rule{:h}, zxd::ZXDiagram{T, P}) where {T, P}
    matches = Match{T}[]
    for v1 in spiders(zxd)
        if spider_type(zxd, v1) == SpiderType.X
            push!(matches, Match{T}([v1]))
        end
    end
    return matches
end

function Base.match(::Rule{:i1}, zxd::ZXDiagram{T, P}) where {T, P}
    matches = Match{T}[]
    for v1 in spiders(zxd)
        if spider_type(zxd, v1) == SpiderType.Z || spider_type(zxd, v1) == SpiderType.X
            if phase(zxd, v1) == 0 && (degree(zxd, v1)) == 2
                push!(matches, Match{T}([v1]))
            end
        end
    end
    return matches
end

function Base.match(::Rule{:i2}, zxd::ZXDiagram{T, P}) where {T, P}
    matches = Match{T}[]
    for v1 in spiders(zxd)
        if spider_type(zxd, v1) == SpiderType.H && (degree(zxd, v1)) == 2
            for v2 in neighbors(zxd, v1)
                if spider_type(zxd, v2) == SpiderType.H && (degree(zxd, v2)) == 2
                    v2 >= v1 && push!(matches, Match{T}([v1, v2]))
                end
            end
        end
    end
    return matches
end

function Base.match(::Rule{:pi}, zxd::ZXDiagram{T, P}) where {T, P}
    matches = Match{T}[]
    for v1 in spiders(zxd)
        if spider_type(zxd, v1) == SpiderType.X && phase(zxd, v1) == one(P) && (degree(zxd, v1)) == 2
            for v2 in neighbors(zxd, v1)
                if spider_type(zxd, v2) == SpiderType.Z
                    push!(matches, Match{T}([v1, v2]))
                end
            end
        end
    end
    return matches
end

function Base.match(::Rule{:c}, zxd::ZXDiagram{T, P}) where {T, P}
    matches = Match{T}[]
    for v1 in spiders(zxd)
        if spider_type(zxd, v1) == SpiderType.X && phase(zxd, v1) == zero(P) && (degree(zxd, v1)) == 1
            for v2 in neighbors(zxd, v1)
                if spider_type(zxd, v2) == SpiderType.Z
                    push!(matches, Match{T}([v1, v2]))
                end
            end
        end
    end
    return matches
end

function Base.match(::Rule{:b}, zxd::ZXDiagram{T, P}) where {T, P}
    matches = Match{T}[]
    for v1 in spiders(zxd)
        if spider_type(zxd, v1) == SpiderType.X && phase(zxd, v1) == zero(P) && (degree(zxd, v1)) == 3
            for v2 in neighbors(zxd, v1)
                if spider_type(zxd, v2) == SpiderType.Z && phase(zxd, v2) == zero(P) && (degree(zxd, v2)) == 3 && mul(zxd.mg, v1, v2) == 1
                    push!(matches, Match{T}([v1, v2]))
                end
            end
        end
    end
    return matches
end

function Base.match(::Rule{:lc}, zxg::ZXGraph{T, P}) where {T, P}
    matches = Match{T}[]
    vs = spiders(zxg)
    vB = [get_inputs(zxg); get_outputs(zxg)]
    for i = 1:length(vB)
        push!(vB, neighbors(zxg, vB[i])[1])
    end
    sort!(vB)
    for v in vs
        if spider_type(zxg, v) == SpiderType.Z &&
            (phase(zxg, v) in (1//2, 3//2))
            if length(searchsorted(vB, v)) == 0
                if degree(zxg, v) == 1
                    # rewrite phase gadgets first
                    pushfirst!(matches, Match{T}([neighbors(zxg, v)[]]))
                    pushfirst!(matches, Match{T}([v]))
                else
                    push!(matches, Match{T}([v]))
                end
            end
        end
    end
    return matches
end

function Base.match(::Rule{:p1}, zxg::ZXGraph{T, P}) where {T, P}
    matches = Match{T}[]
    vs = spiders(zxg)
    vB = [get_inputs(zxg); get_outputs(zxg)]
    for i = 1:length(vB)
        push!(vB, neighbors(zxg, vB[i])[1])
    end
    sort!(vB)
    for v1 in vs
        if spider_type(zxg, v1) == SpiderType.Z && length(searchsorted(vB, v1)) == 0 &&
            (phase(zxg, v1) in (0, 1))
            for v2 in neighbors(zxg, v1)
                if spider_type(zxg, v2) == SpiderType.Z && length(searchsorted(vB, v2)) == 0 &&
                    (phase(zxg, v2) in (0, 1)) && v2 > v1
                    push!(matches, Match{T}([v1, v2]))
                end
            end
        end
    end
    return matches
end

function Base.match(::Rule{:pab}, zxg::ZXGraph{T, P}) where {T, P}
    matches = Match{T}[]
    vs = spiders(zxg)
    vB = [get_inputs(zxg); get_outputs(zxg)]
    for i = 1:length(vB)
        push!(vB, neighbors(zxg, vB[i])[1])
    end
    sort!(vB)
    for v2 in vB
        if spider_type(zxg, v2) == SpiderType.Z && length(neighbors(zxg, v2)) > 2
            for v1 in neighbors(zxg, v2)
                if spider_type(zxg, v1) == SpiderType.Z && length(searchsorted(vB, v1)) == 0 &&
                    (phase(zxg, v1) in (0, 1))
                    push!(matches, Match{T}([v1, v2]))
                end
            end
        end
    end
    return matches
end

function Base.match(::Rule{:p2}, zxg::ZXGraph{T, P}) where {T, P}
    matches = Match{T}[]
    vs = spiders(zxg)
    vB = [get_inputs(zxg); get_outputs(zxg)]
    for i = 1:length(vB)
        push!(vB, neighbors(zxg, vB[i])[1])
    end
    sort!(vB)
    gadgets = T[]
    for v in vs
        if spider_type(zxg, v) == SpiderType.Z && length(neighbors(zxg, v)) == 1
            push!(gadgets, v, neighbors(zxg, v)[1])
        end
    end
    sort!(gadgets)

    v_matched = T[]

    for v1 in vs
        if spider_type(zxg, v1) == SpiderType.Z && length(searchsorted(vB, v1)) == 0 &&
            (degree(zxg, v1)) > 1 && (rem(phase(zxg, v1), 1//2) != 0) &&
            length(neighbors(zxg, v1)) > 1 && v1 ∉ v_matched
            for v2 in neighbors(zxg, v1)
                if spider_type(zxg, v2) == SpiderType.Z &&
                    length(searchsorted(vB, v2)) == 0 &&
                    (phase(zxg, v2) in (0, 1))
                    if length(searchsorted(gadgets, v2)) == 0  && v2 ∉ v_matched
                        push!(matches, Match{T}([v1, v2]))
                        push!(v_matched, v1, v2)
                    end
                end
            end
        end
    end
    return matches
end

function Base.match(::Rule{:p3}, zxg::ZXGraph{T, P}) where {T, P}
    matches = Match{T}[]
    vs = spiders(zxg)
    vB = [get_inputs(zxg); get_outputs(zxg)]
    for i = 1:length(vB)
        push!(vB, neighbors(zxg, vB[i])[1])
    end
    sort!(vB)
    gadgets = T[]
    for v in vs
        if spider_type(zxg, v) == SpiderType.Z && length(neighbors(zxg, v)) == 1
            push!(gadgets, v, neighbors(zxg, v)[1])
        end
    end
    sort!(gadgets)

    v_matched = T[]

    for v1 in vB
        if spider_type(zxg, v1) == SpiderType.Z && length(searchsorted(vB, v1)) > 0 &&
            (rem(phase(zxg, v1), 1//2) != 0) && length(neighbors(zxg, v1)) > 1 &&
            v1 ∉ v_matched
            for v2 in neighbors(zxg, v1)
                if spider_type(zxg, v2) == SpiderType.Z && length(searchsorted(vB, v2)) == 0 &&
                    (phase(zxg, v2) in (0, 1)) && length(searchsorted(gadgets, v2)) == 0 && v2 ∉ v_matched
                    push!(matches, Match{T}([v1, v2]))
                    push!(v_matched, v1, v2)
                end
            end
        end
    end
    return matches
end

function Base.match(::Rule{:id}, zxg::ZXGraph{T, P}) where {T, P}
    matches = Match{T}[]
    for v2 in spiders(zxg)
        nb2 = neighbors(zxg, v2)
        if spider_type(zxg, v2) == SpiderType.Z && length(nb2) == 2
            v1, v3 = nb2
            if phase(zxg, v2) == 0
                if spider_type(zxg, v1) == SpiderType.Z && spider_type(zxg, v3) == SpiderType.Z
                    push!(matches, Match{T}([v1, v2, v3]))
                end
            else phase(zxg, v2) == 1
                if spider_type(zxg, v1) == SpiderType.Z && spider_type(zxg, v3) == SpiderType.Z
                    if degree(zxg, v1) == 1
                        push!(matches, Match{T}([v1, v2, v3]))
                    elseif degree(zxg, v3) == 1
                        push!(matches, Match{T}([v1, v2, v3]))
                    end
                end
            end
        end
    end
    return matches
end

function Base.match(::Rule{:gf}, zxg::ZXGraph{T, P}) where {T, P}
    matches = Match{T}[]
    vs = spiders(zxg)
    gad_ids = vs[[spider_type(zxg, v) == SpiderType.Z && (degree(zxg, v)) == 1 for v in vs]]
    gads = [(v, neighbors(zxg, v)[1], setdiff(neighbors(zxg, neighbors(zxg, v)[1]), [v])) for v in gad_ids]

    for i = 1:length(gads)
        v1, v2, gad_v = gads[i]
        for j in (i+1):length(gads)
            u1, u2, gad_u = gads[j]
            if gad_u == gad_v && phase(zxg, v2) in (zero(P), one(P)) && phase(zxg, u2) in (zero(P), one(P))
                push!(matches, Match{T}([v1, v2, u1, u2]))
            end
        end
    end
    return matches
end

function Base.match(::Rule{:scalar}, zxg::Union{ZXGraph{T, P}, ZXDiagram{T, P}}) where {T, P}
    matches = Match{T}[]
    vs = spiders(zxg)
    for v in vs
        if degree(zxg, v) == 0
            if spider_type(zxg, v) in (SpiderType.Z, SpiderType.X)
                push!(matches, Match{T}([v]))
            end
        end
    end
    return matches
end

"""
    rewrite!(r, zxd, matches)

Rewrite a ZX-diagram `zxd` with rule `r` for all vertices in `matches`. `matches`
can be a vector of `Match` or just an instance of `Match`.
"""
function rewrite!(r::AbstractRule, zxd::AbstractZXDiagram{T, P}, matches::Vector{Match{T}}) where {T, P}
    for each in matches
        rewrite!(r, zxd, each)
    end
    return zxd
end

function rewrite!(r::AbstractRule, zxd::AbstractZXDiagram{T, P}, matched::Match{T}) where {T, P}
    vs = matched.vertices
    if check_rule(r, zxd, vs)
        rewrite!(r, zxd, vs)
    end
    return zxd
end

function check_rule(r::Rule{:f}, zxd::ZXDiagram{T, P}, vs::Vector{T}) where {T, P}
    v1, v2 = vs
    (has_vertex(zxd.mg, v1) && has_vertex(zxd.mg, v2)) || return false
    if spider_type(zxd, v1) == SpiderType.Z || spider_type(zxd, v1) == SpiderType.X
        if v2 in neighbors(zxd, v1)
            if spider_type(zxd, v1) == spider_type(zxd, v2) && v2 >= v1
                return true
            end
        end
    end
    return false
end

function rewrite!(r::Rule{:f}, zxd::ZXDiagram{T, P}, vs::Vector{T}) where {T, P}
    v1, v2 = vs
    for v3 in neighbors(zxd, v2)
        if v3 != v1
            add_edge!(zxd, v1, v3)
        end
    end
    set_phase!(zxd, v1, phase(zxd, v1)+phase(zxd, v2))
    rem_spider!(zxd, v2)
    return zxd
end

function check_rule(r::Rule{:h}, zxd::ZXDiagram{T, P}, vs::Vector{T}) where {T, P}
    @inbounds v1 = vs[1]
    has_vertex(zxd.mg, v1) || return false
    if spider_type(zxd, v1) == SpiderType.X
        return true
    end
    return false
end

function rewrite!(r::Rule{:h}, zxd::ZXDiagram{T, P}, vs::Vector{T}) where {T, P}
    @inbounds v1 = vs[1]
    for v2 in neighbors(zxd, v1)
        if v2 != v1
            insert_spider!(zxd, v1, v2, SpiderType.H)
        end
    end
    zxd.st[v1] = SpiderType.Z
    return zxd
end

function check_rule(r::Rule{:i1}, zxd::ZXDiagram{T, P}, vs::Vector{T}) where {T, P}
    @inbounds v1 = vs[1]
    has_vertex(zxd.mg, v1) || return false
    if spider_type(zxd, v1) == SpiderType.Z || spider_type(zxd, v1) == SpiderType.X
        if phase(zxd, v1) == 0 && (degree(zxd, v1)) == 2
            return true
        end
    end
    return false
end

function rewrite!(r::Rule{:i1}, zxd::ZXDiagram{T, P}, vs::Vector{T}) where {T, P}
    @inbounds v1 = vs[1]
    v2, v3 = neighbors(zxd, v1, count_mul = true)
    add_edge!(zxd, v2, v3)
    rem_spider!(zxd, v1)
    return zxd
end

function check_rule(r::Rule{:i2}, zxd::ZXDiagram{T, P}, vs::Vector{T}) where {T, P}
    @inbounds v1, v2 = vs
    (has_vertex(zxd.mg, v1) && has_vertex(zxd.mg, v2)) || return false
    if spider_type(zxd, v1) == SpiderType.H && spider_type(zxd, v2) == SpiderType.H && has_edge(zxd.mg, v1, v2)
        if (degree(zxd, v1)) == 2 && (degree(zxd, v2)) == 2
            return true
        end
    end
    return false
end

function rewrite!(r::Rule{:i2}, zxd::ZXDiagram{T, P}, vs::Vector{T}) where {T, P}
    @inbounds v1, v2 = vs
    nb1 = neighbors(zxd, v1, count_mul = true)
    nb2 = neighbors(zxd, v2, count_mul = true)
    @inbounds v3 = (nb1[1] == v2 ? nb1[2] : nb1[1])
    @inbounds v4 = (nb2[1] == v1 ? nb2[2] : nb2[1])
    add_edge!(zxd, v3, v4)
    rem_spiders!(zxd, [v1, v2])
    return zxd
end

function check_rule(r::Rule{:pi}, zxd::ZXDiagram{T, P}, vs::Vector{T}) where {T, P}
    v1, v2 = vs
    (has_vertex(zxd.mg, v1) && has_vertex(zxd.mg, v2)) || return false
    if spider_type(zxd, v1) == SpiderType.X && phase(zxd, v1) == one(phase(zxd, v1)) &&
            (degree(zxd, v1)) == 2
        if v2 in neighbors(zxd, v1)
            if spider_type(zxd, v2) == SpiderType.Z
                return true
            end
        end
    end
    return false
end

function rewrite!(r::Rule{:pi}, zxd::ZXDiagram{T, P}, vs::Vector{T}) where {T, P}
    v1, v2 = vs
    add_global_phase!(zxd, phase(zxd, v2))
    set_phase!(zxd, v2, -phase(zxd, v2))
    nb = neighbors(zxd, v2, count_mul = true)
    for v3 in nb
        # TODO
        v3 != v1 && insert_spider!(zxd, v2, v3, SpiderType.X, phase(zxd, v1))
    end
    if neighbors(zxd, v1) != [v2]
        add_edge!(zxd, neighbors(zxd, v1))
        rem_spider!(zxd, v1)
    end
    return zxd
end

function check_rule(r::Rule{:c}, zxd::ZXDiagram{T, P}, vs::Vector{T}) where {T, P}
    v1, v2 = vs
    (has_vertex(zxd.mg, v1) && has_vertex(zxd.mg, v1)) || return false
    if spider_type(zxd, v1) == SpiderType.X && phase(zxd, v1) == zero(phase(zxd, v1)) && (degree(zxd, v1)) == 1
        if v2 in neighbors(zxd, v1)
            if spider_type(zxd, v2) == SpiderType.Z
                return true
            end
        end
    end
    return false
end

function rewrite!(r::Rule{:c}, zxd::ZXDiagram{T, P}, vs::Vector{T}) where {T, P}
    v1, v2 = vs
    ph = phase(zxd, v1)
    rem_spider!(zxd, v1)
    add_power!(zxd, 1)
    nb = neighbors(zxd, v2, count_mul = true)
    for v3 in nb
        add_spider!(zxd, SpiderType.X, ph, [v3])
        add_power!(zxd, -1)
    end
    rem_spider!(zxd, v2)
    return zxd
end

function check_rule(r::Rule{:b}, zxd::ZXDiagram{T, P}, vs::Vector{T}) where {T, P}
    v1, v2 = vs
    (has_vertex(zxd.mg, v1) && has_vertex(zxd.mg, v2)) || return false
    if spider_type(zxd, v1) == SpiderType.X && phase(zxd, v1) == 0 && (degree(zxd, v1)) == 3
        if v2 in neighbors(zxd, v1)
            if spider_type(zxd, v2) == SpiderType.Z && phase(zxd, v2) == 0 && (degree(zxd, v2)) == 3 && mul(zxd.mg, v1, v2) == 1
                return true
            end
        end
    end
    return false
end

function rewrite!(r::Rule{:b}, zxd::ZXDiagram{T, P}, vs::Vector{T}) where {T, P}
    v1, v2 = vs
    nb1 = neighbors(zxd, v1)
    nb2 = neighbors(zxd, v2)
    v3, v4 = nb1[nb1 .!= v2]
    v5, v6 = nb2[nb2 .!= v1]

    # TODO
    a1 = insert_spider!(zxd, v1, v3, SpiderType.Z)[1]
    a2 = insert_spider!(zxd, v1, v4, SpiderType.Z)[1]
    a3 = insert_spider!(zxd, v2, v5, SpiderType.X)[1]
    a4 = insert_spider!(zxd, v2, v6, SpiderType.X)[1]
    rem_spiders!(zxd, [v1, v2])

    add_edge!(zxd, a1, a3)
    add_edge!(zxd, a1, a4)
    add_edge!(zxd, a2, a3)
    add_edge!(zxd, a2, a4)
    add_power!(zxd, 1)
    return zxd
end

function check_rule(::Rule{:lc}, zxg::ZXGraph{T, P}, vs::Vector{T}) where {T, P}
    @inbounds v = vs[1]
    has_vertex(zxg.mg, v) || return false
    if has_vertex(zxg.mg, v)
        if spider_type(zxg, v) == SpiderType.Z &&
            (phase(zxg, v) in (1//2, 3//2))
            if is_interior(zxg, v)
                return true
            end
        end
    end
    return false
end

function rewrite!(r::Rule{:lc}, zxg::ZXGraph{T, P}, vs::Vector{T}) where {T, P}
    @inbounds v = vs[1]
    phase_v = phase(zxg, v)
    if phase_v == 1//2
        add_global_phase!(zxg, P(1//4))
    else
        add_global_phase!(zxg, P(-1//4))
    end
    nb = neighbors(zxg, v)
    n = length(nb)
    add_power!(zxg, (n-1)*(n-2)//2)
    rem_spider!(zxg, v)
    for u1 in nb, u2 in nb
        if u2 > u1
            add_edge!(zxg, u1, u2, EdgeType.HAD)
        end
    end
    for u in nb
        set_phase!(zxg, u, phase(zxg, u)-phase_v)
    end
    return zxg
end

function check_rule(::Rule{:p1}, zxg::ZXGraph{T, P}, vs::Vector{T}) where {T, P}
    v1, v2 = vs
    (has_vertex(zxg.mg, v1) && has_vertex(zxg.mg, v2)) || return false
    if has_vertex(zxg.mg, v1)
        if spider_type(zxg, v1) == SpiderType.Z && is_interior(zxg, v1) &&
            (phase(zxg, v1) in (0, 1))
            if v2 in neighbors(zxg, v1)
                if spider_type(zxg, v2) == SpiderType.Z && is_interior(zxg, v2) &&
                    (phase(zxg, v2) in (0, 1)) && v2 > v1
                    return true
                end
            end
        end
    end
    return false
end

function rewrite!(::Rule{:p1}, zxg::ZXGraph{T, P}, vs::Vector{T}) where {T, P}
    u, v = vs
    phase_u = phase(zxg, u)
    phase_v = phase(zxg, v)
    add_global_phase!(zxg, phase_u*phase_v)
    nb_u = setdiff(neighbors(zxg, u), [v])
    nb_v = setdiff(neighbors(zxg, v), [u])

    U = setdiff(nb_u, nb_v)
    V = setdiff(nb_v, nb_u)
    W = intersect(nb_u, nb_v)
    add_power!(zxg, (length(U)+length(V)-2)*length(W) + (length(U)-1)*(length(V)-1))

    rem_spiders!(zxg, vs)
    for u0 in U, v0 in V
        add_edge!(zxg, u0, v0)
    end
    for u0 in U, w0 in W
        add_edge!(zxg, u0, w0)
    end
    for v0 in V, w0 in W
        add_edge!(zxg, v0, w0)
    end
    for u0 in U
        set_phase!(zxg, u0, phase(zxg, u0)+phase_v)
    end
    for v0 in V
        set_phase!(zxg, v0, phase(zxg, v0)+phase_u)
    end
    for w0 in W
        set_phase!(zxg, w0, phase(zxg, w0)+phase_u+phase_v+1)
    end
    return zxg
end

function check_rule(::Rule{:pab}, zxg::ZXGraph{T, P}, vs::Vector{T}) where {T, P}
    v1, v2 = vs
    (has_vertex(zxg.mg, v1) && has_vertex(zxg.mg, v2)) || return false
    if has_vertex(zxg.mg, v1)
        if spider_type(zxg, v1) == SpiderType.Z && is_interior(zxg, v1) &&
            (phase(zxg, v1) in (0, 1))
            if v2 in neighbors(zxg, v1)
                if spider_type(zxg, v2) == SpiderType.Z && !is_interior(zxg, v2) &&
                    length(neighbors(zxg, v2)) > 2
                    return true
                end
            end
        end
    end
    return false
end

function rewrite!(::Rule{:pab}, zxg::ZXGraph{T, P}, vs::Vector{T}) where {T, P}
    u, v = vs
    phase_v = phase(zxg, v)
    nb_v = neighbors(zxg, v)
    v_bound = zero(T)
    for v0 in nb_v
        if spider_type(zxg, v0) != SpiderType.Z
            v_bound = v0
            break
        end
    end
    @inbounds if is_hadamard(zxg, v, v_bound)
        # TODO
        w = insert_spider!(zxg, v, v_bound)[1]

        zxg.et[(v_bound, w)] = EdgeType.SIM
        v_bound_master = v_bound
        v_master = neighbors(zxg.master, v_bound_master)[1]
        w_master = insert_spider!(zxg.master, v_bound_master, v_master,  SpiderType.Z)[1]
        # @show w, w_master
        zxg.phase_ids[w] = (w_master, 1)

        # set_phase!(zxg, w, phase(zxg, v))
        # zxg.phase_ids[w] = zxg.phase_ids[v]
        # set_phase!(zxg, v, zero(P))
        # zxg.phase_ids[v] = (v, 1)
    else
        # TODO
        w = insert_spider!(zxg, v, v_bound)[1]
        # insert_spider!(zxg, w, v_bound, phase_v)
        # w = neighbors(zxg, v_bound)[1]
        # set_phase!(zxg, w, phase(zxg, v))
        # zxg.phase_ids[w] = zxg.phase_ids[v]

        v_bound_master = v_bound
        v_master = neighbors(zxg.master, v_bound_master)[1]
        w_master = insert_spider!(zxg.master, v_bound_master, v_master,  SpiderType.X)[1]
        # @show w, w_master
        zxg.phase_ids[w] = (w_master, 1)

        # set_phase!(zxg, v, zero(P))
        # zxg.phase_ids[v] = (v, 1)
        # rem_edge!(zxg, w, v_bound)
        # add_edge!(zxg, w, v_bound, EdgeType.SIM)
    end
    return rewrite!(Rule{:p1}(), zxg, Match{T}([u, v]))
end

function check_rule(::Rule{:p2}, zxg::ZXGraph{T, P}, vs::Vector{T}) where {T, P}
    v1, v2 = vs
    if has_vertex(zxg.mg, v1)
        if spider_type(zxg, v1) == SpiderType.Z && is_interior(zxg, v1) &&
            (degree(zxg, v1)) > 1 && (rem(phase(zxg, v1), 1//2) != 0) &&
            length(neighbors(zxg, v1)) > 1
            if v2 in neighbors(zxg, v1)
                if spider_type(zxg, v2) == SpiderType.Z && is_interior(zxg, v2) &&
                    (phase(zxg, v2) in (0, 1))
                    if all(length(neighbors(zxg, u)) > 1 for u in neighbors(zxg, v2))
                        return true
                    end
                end
            end
        end
    end
    return false
end

function rewrite!(::Rule{:p2}, zxg::ZXGraph{T, P}, vs::Vector{T}) where {T, P}
    # u has non-Clifford phase
    u, v = vs
    phase_u = phase(zxg, u)
    phase_v = phase(zxg, v)
    nb_u = setdiff(neighbors(zxg, u), [v])
    nb_v = setdiff(neighbors(zxg, v), [u])

    U = setdiff(nb_u, nb_v)
    V = setdiff(nb_v, nb_u)
    W = intersect(nb_u, nb_v)
    add_power!(zxg, (length(U)+length(V)-1)*length(W) + length(U)*(length(V)-1))

    phase_id_u = zxg.phase_ids[u]
    sgn_phase_v = iseven(Phase(phase_v)) ? 1 : -1

    if sgn_phase_v < 0
        zxg.phase_ids[u] = (phase_id_u[1], -phase_id_u[2])
        phase_id_u = zxg.phase_ids[u]
    end
    rem_spider!(zxg, u)
    for u0 in U, v0 in V
        add_edge!(zxg, u0, v0)
    end
    for u0 in U, w0 in W
        add_edge!(zxg, u0, w0)
    end
    for v0 in V, w0 in W
        add_edge!(zxg, v0, w0)
    end
    for u0 in U
        set_phase!(zxg, u0, phase(zxg, u0)+phase_v)
    end
    for w0 in W
        set_phase!(zxg, w0, phase(zxg, w0)+phase_v+1)
    end
    gad = add_spider!(zxg, SpiderType.Z, P(sgn_phase_v*phase_u))
    add_edge!(zxg, v, gad)
    set_phase!(zxg, v, zero(P))
    zxg.phase_ids[gad] = phase_id_u
    zxg.phase_ids[v] = (v, 1)

    rem_vertex!(zxg.layout, v)
    return zxg
end

function check_rule(::Rule{:p3}, zxg::ZXGraph{T, P}, vs::Vector{T}) where {T, P}
    v1, v2 = vs
    if has_vertex(zxg.mg, v1)
        if spider_type(zxg, v1) == SpiderType.Z && !is_interior(zxg, v1) &&
            (rem(phase(zxg, v1), 1//2) != 0) && length(neighbors(zxg, v1)) > 1
            if v2 in neighbors(zxg, v1)
                if spider_type(zxg, v2) == SpiderType.Z && is_interior(zxg, v2) &&
                    (phase(zxg, v2) in (0, 1))
                    if all(length(neighbors(zxg, u)) > 1 for u in neighbors(zxg, v2))
                        return true
                    end
                end
            end
        end
    end
    return false
end

function rewrite!(::Rule{:p3}, zxg::ZXGraph{T, P}, vs::Vector{T}) where {T, P}
    # u is the boundary spider
    u, v = vs
    phase_u = phase(zxg, u)
    phase_v = phase(zxg, v)
    nb_u = setdiff(neighbors(zxg, u), [v])
    nb_v = setdiff(neighbors(zxg, v), [u])
    bd_u = nb_u[findfirst([u0 in get_inputs(zxg) || u0 in get_outputs(zxg) for u0 in nb_u])]
    setdiff!(nb_u, [bd_u])

    U = setdiff(nb_u, nb_v)
    V = setdiff(nb_v, nb_u)
    W = intersect(nb_u, nb_v)
    add_power!(zxg, (length(U)+length(V))*length(W) + (length(U)+1)*(length(V)-1))

    phase_id_u = zxg.phase_ids[u]
    sgn_phase_v = iseven(Phase(phase_v)) ? 1 : -1

    if sgn_phase_v < 0
        zxg.phase_ids[u] = (phase_id_u[1], -phase_id_u[2])
        phase_id_u = zxg.phase_ids[u]
    end
    phase_id_v = zxg.phase_ids[v]
    rem_edge!(zxg, u, v)
    for u0 in U, v0 in V
        add_edge!(zxg, u0, v0)
    end
    for u0 in U, w0 in W
        add_edge!(zxg, u0, w0)
    end
    for v0 in V, w0 in W
        add_edge!(zxg, v0, w0)
    end
    for u0 in U
        rem_edge!(zxg, u, u0)
        set_phase!(zxg, u0, phase(zxg, u0)+phase_v)
    end
    for v0 in V
        add_edge!(zxg, u, v0)
    end
    for w0 in W
        set_phase!(zxg, w0, phase(zxg, w0)+phase_v+1)
    end
    set_phase!(zxg, v, zero(P))
    set_phase!(zxg, u, phase_v)
    gad = add_spider!(zxg, SpiderType.Z, P(sgn_phase_v*phase_u))
    add_edge!(zxg, v, gad)
    zxg.phase_ids[gad] = phase_id_u
    zxg.phase_ids[u] = phase_id_v
    zxg.phase_ids[v] = (v, 1)

    rem_vertex!(zxg.layout, v)

    if is_hadamard(zxg, u, bd_u)
        rem_edge!(zxg, u, bd_u)
        add_edge!(zxg, u, bd_u, EdgeType.SIM)
    else
        rem_edge!(zxg, u, bd_u)
        add_edge!(zxg, u, bd_u)
    end
    return zxg
end

function rewrite!(::Rule{:pivot}, zxg::ZXGraph{T, P}, vs::Vector{T}) where {T, P}
    # This rule should be only used in circuit extraction.
    # This rule will do pivoting on u, v but preserve u, v.
    # And the scalars are not considered in this rule.
    # gadget_u is the non-Clifford spider
    u, gadget_u, v = vs
    phase_u = phase(zxg, u)
    phase_v = phase(zxg, v)

    for v0 in neighbors(zxg, v)
        if spider_type(zxg, v0) != SpiderType.Z
            if is_hadamard(zxg, v0, v)
                v1 = insert_spider!(zxg, v0, v)
                insert_spider!(zxg, v0, v1)
            else
                insert_spider!(zxg, v0, v)
            end
            break
        end
    end

    nb_u = setdiff(neighbors(zxg, u), [v, gadget_u])
    nb_v = setdiff(neighbors(zxg, v), [u])

    U = setdiff(nb_u, nb_v)
    V = setdiff(nb_v, nb_u)
    W = intersect(nb_u, nb_v)
    add_power!(zxg, length(U)*length(V) + length(V)*length(W) + length(W)*length(U))

    phase_id_gadget_u = zxg.phase_ids[gadget_u]
    phase_gadget_u = phase(zxg, gadget_u)
    if !iseven(Phase(phase_u))
        zxg.phase_ids[gadget_u] = (phase_id_gadget_u[1], -phase_id_gadget_u[2])
        phase_id_gadget_u = zxg.phase_ids[gadget_u]
        phase_gadget_u = -phase(zxg, gadget_u)
    end

    for u0 in U, v0 in V
        add_edge!(zxg, u0, v0)
    end
    for u0 in U, w0 in W
        add_edge!(zxg, u0, w0)
    end
    for v0 in V, w0 in W
        add_edge!(zxg, v0, w0)
    end

    for w0 in W
        set_phase!(zxg, w0, phase(zxg, w0)+1)
    end

    set_phase!(zxg, v, phase_gadget_u)
    zxg.phase_ids[v] = phase_id_gadget_u
    zxg.phase_ids[u] = (u, 1)

    rem_spider!(zxg, gadget_u)
    return zxg
end

function check_rule(::Rule{:id}, zxg::ZXGraph{T, P}, vs::Vector{T}) where {T, P}
    v1, v2, v3 = vs
    if has_vertex(zxg.mg, v2)
        nb2 = neighbors(zxg, v2)
        if spider_type(zxg, v2) == SpiderType.Z && length(nb2) == 2
            (v1 in nb2 && v3 in nb2) || return false
            if phase(zxg, v2) == 0
                if spider_type(zxg, v1) == SpiderType.Z && spider_type(zxg, v3) == SpiderType.Z
                    return true
                end
            else phase(zxg, v2) == 1
                if spider_type(zxg, v1) == SpiderType.Z && spider_type(zxg, v3) == SpiderType.Z
                    return degree(zxg, v1) == 1 || degree(zxg, v3) == 1
                end
            end
        end
    end
    return false
end

function rewrite!(::Rule{:id}, zxg::ZXGraph{T, P}, vs::Vector{T}) where {T, P}
    v1, v2, v3 = vs
    if phase(zxg, v2) == 1
        set_phase!(zxg, v2, zero(P))
        set_phase!(zxg, v1, -phase(zxg, v1))
        zxg.phase_ids[v1] = (zxg.phase_ids[v1][1], -zxg.phase_ids[v1][2])
    end
    set_phase!(zxg, v3, phase(zxg, v3)+phase(zxg, v1))
    id1, mul1 = zxg.phase_ids[v1]
    id3, mul3 = zxg.phase_ids[v3]
    set_phase!(zxg.master, id3, (mul3 * phase(zxg.master, id3) + mul1 * phase(zxg.master, id1)) * mul3)
    set_phase!(zxg.master, id1, zero(P))
    for v in neighbors(zxg, v1)
        v == v2 && continue
        add_edge!(zxg, v, v3, is_hadamard(zxg, v, v1) ? EdgeType.HAD : EdgeType.SIM)
    end
    rem_spiders!(zxg, [v1, v2])
    return zxg
end

function check_rule(::Rule{:gf}, zxg::ZXGraph{T, P}, vs::Vector{T}) where {T, P}
    @inbounds if all(has_vertex(zxg.mg, v) for v in vs)
        v1, v2, u1, u2 = vs
        if spider_type(zxg, v1) == SpiderType.Z && (degree(zxg, v1)) == 1 &&
            spider_type(zxg, u1) == SpiderType.Z && (degree(zxg, u1)) == 1
            if v2 == neighbors(zxg, v1)[1] && u2 == neighbors(zxg, u1)[1]
                gad_v = setdiff(neighbors(zxg, v2), [v1])
                gad_u = setdiff(neighbors(zxg, u2), [u1])
                if gad_u == gad_v && phase(zxg, v2) in (zero(P), one(P)) && phase(zxg, u2) in (zero(P), one(P))
                    return true
                end
            end
        end
    end
    return false
end

function rewrite!(::Rule{:gf}, zxg::ZXGraph{T, P}, vs::Vector{T}) where {T, P}
    v1, v2, u1, u2 = vs
    if phase(zxg, v2) == 1
        add_global_phase!(zxg, phase(zxg, v1))
        set_phase!(zxg, v2, zero(P))
        set_phase!(zxg, v1, -phase(zxg, v1))
        zxg.phase_ids[v1] = (zxg.phase_ids[v1][1], -zxg.phase_ids[v1][2])
    end
    if phase(zxg, u2) == 1
        add_global_phase!(zxg, phase(zxg, u1))
        set_phase!(zxg, u2, zero(P))
        set_phase!(zxg, u1, -phase(zxg, u1))
        zxg.phase_ids[u1] = (zxg.phase_ids[u1][1], -zxg.phase_ids[u1][2])
    end

    set_phase!(zxg, v1, phase(zxg, v1)+phase(zxg, u1))

    idv, mulv = zxg.phase_ids[v1]
    idu, mulu = zxg.phase_ids[u1]
    set_phase!(zxg.master, idv, (mulv * phase(zxg.master,idv) + mulu * phase(zxg.master,idu)) * mulv)
    set_phase!(zxg.master, idu, zero(P))

    add_power!(zxg, degree(zxg, v2)-2)
    rem_spiders!(zxg, [u1, u2])
    return zxg
end

function check_rule(::Rule{:scalar}, zxg::Union{ZXGraph{T, P}, ZXDiagram{T, P}}, vs::Vector{T}) where {T, P}
    @inbounds v = vs[1]
    if has_vertex(zxg.mg, v)
        if degree(zxg, v) == 0
            if spider_type(zxg, v) in (SpiderType.Z, SpiderType.X)
                return true
            end
        end
    end
    return false
end

function rewrite!(::Rule{:scalar}, zxg::Union{ZXGraph{T, P}, ZXDiagram{T, P}}, vs::Vector{T}) where {T, P}
    @inbounds v = vs[1]
    rem_spider!(zxg, v)
    return zxg
end


text:
"""
    all_backedges(mi::MethodInstance)

Return a list of all backedges (direct and indirect) of `mi`.
"""
function all_backedges(mi::MethodInstance)
    backedges = Set{MethodInstance}()
    visit_backedges(x->(push!(backedges, x); true), mi)
    delete!(backedges, mi)
    return collect(backedges)
end

"""
    with_all_backedges(itr)

Return all MethodInstances detected when iterating through items in `itr` and any
their backedges. The result includes both MethodTable and MethodInstance backedges.
"""
function with_all_backedges(iter)
    backedges = Set{MethodInstance}()
    visited = IdSet{MethodInstance}()
    for item in iter
        visit_backedges(item, visited) do edge
            if isa(edge, MethodInstance)
                push!(backedges, edge)
            end
            true
        end
    end
    return collect(backedges)
end

"""
    terminal_backedges(mi::MethodInstance)

Obtain the "ultimate callers" of `mi`, i.e., the reason(s) `mi` was compiled.
"""
function terminal_backedges(mi::MethodInstance)
    backedges = Set{MethodInstance}()
    visit_backedges(mi) do x
        if !isdefined(x, :backedges) || isempty(x.backedges)
            push!(backedges, x)
        end
        true
    end
    delete!(backedges, mi)
    return collect(backedges)
end

"""
    direct_backedges(f::Function; skip=true)

Collect all backedges for a function `f` as pairs `instance=>caller` or `sig=>caller` pairs.
The latter occur for MethodTable backedges.
If `skip` is `true`, any `caller` listed in a MethodTable backedge is omitted from the instance backedges.
"""
function direct_backedges(f::Union{Method,Callable}; skip::Bool=true)
    bes = []
    _skip = Set{MethodInstance}()
    visit(f) do item
        if isa(item, MethodTable)
            mt = item::MethodTable
            if isdefined(mt, :backedges)
                sigmis = mt.backedges::Vector{Any}
                for i = 1:2:length(sigmis)
                    sig, mi = sigmis[i], sigmis[i+1]
                    push!(bes, Pair{Any,MethodInstance}(sig, mi))
                    push!(_skip, mi)
                end
            end
            return false
        elseif isa(item, MethodInstance)
            callee = item::MethodInstance
            if isdefined(callee, :backedges)
                for caller in callee.backedges
                    skip && caller ∈ _skip && continue
                    push!(bes, callee=>caller)
                end
            end
            return false
        end
        true
    end
    return bes
end

"""
    direct_backedges(mi::MethodInstance)

A vector of all direct backedges of `mi`. This is equivalent to `mi.backedges` except that it's "safe,"
meaning it returns an empty list even when `mi.backedges` is not defined.
"""
function direct_backedges(mi::MethodInstance)
    out = MethodInstance[]
    if isdefined(mi, :backedges)
        append!(out, mi.backedges)
    end
    return out
end


text:
using DSGE, FileIO, JLD2, ModelConstructors, Test, Random, Dates
path = dirname(@__FILE__)

# Set up arguments
m = AnSchorfheide(testing = true)
m <= Setting(:date_forecast_start, quartertodate("2015-Q4"))
m <= Setting(:forecast_horizons, 1)

system, histshocks = JLD2.jldopen("$path/../reference/forecast_args.jld2","r") do file
    read(file, "system"), read(file, "histshocks")
end

# Read expected output
exp_states, exp_obs, exp_pseudo =
    JLD2.jldopen("$path/../reference/shock_decompositions_out.jld2", "r") do file
        read(file, "exp_states"), read(file, "exp_obs"), read(file, "exp_pseudo")
    end

# With shockdec_startdate not null
states, obs, pseudo = shock_decompositions(m, system, histshocks)

@testset "Test shockdec with non-null startdate" begin
    @test @test_matrix_approx_eq exp_states[:startdate] states
    @test @test_matrix_approx_eq exp_obs[:startdate]    obs
    @test @test_matrix_approx_eq exp_pseudo[:startdate] pseudo
end

# With shockdec_startdate null
#m <= Setting(:shockdec_startdate, Nullable{Date}())
m <= Setting(:shockdec_startdate, nothing)
states, obs, pseudo = shock_decompositions(m, system, histshocks)

@testset "Test shockdec with null startdate" begin
    @test @test_matrix_approx_eq exp_states[:no_startdate] states
    @test @test_matrix_approx_eq exp_obs[:no_startdate]    obs
    @test @test_matrix_approx_eq exp_pseudo[:no_startdate] pseudo
end

@testset "Deterministic trends" begin
    dettrend = DSGE.deterministic_trends(m, system, zeros(8))
    @test @test_matrix_approx_eq dettrend[1] zeros(size(dettrend[1]))
    @test @test_matrix_approx_eq dettrend[2] zeros(size(dettrend[2]))
    @test @test_matrix_approx_eq dettrend[3] zeros(size(dettrend[3]))
end

@testset "Trends" begin
    out = DSGE.trends(system)
    @test @test_matrix_approx_eq out[1] system[:CCC]
    @test @test_matrix_approx_eq out[2] system[:ZZ] * system[:CCC] + system[:DD]
    @test @test_matrix_approx_eq out[3] system[:ZZ_pseudo] * system[:CCC] + system[:DD_pseudo]
end

# Check trends for "fake" regime-switching system
reg_sys = RegimeSwitchingSystem(System{Float64}[system, system])
@testset "Regime-switching trends" begin
    out = DSGE.trends(reg_sys)
    for i in 1:2
        @test @test_matrix_approx_eq out[1][:, i] reg_sys[i, :CCC]
        @test @test_matrix_approx_eq out[2][:, i] reg_sys[i, :ZZ] * reg_sys[i, :CCC] + reg_sys[i, :DD]
        @test @test_matrix_approx_eq out[3][:, i] reg_sys[i, :ZZ_pseudo] * reg_sys[i, :CCC] + reg_sys[i, :DD_pseudo]
    end
end

m.settings[:regime_dates] = Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2010, 3, 31)))
m.test_settings[:regime_dates] = Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m), 2 => Date(2010, 3, 31)))
m.settings[:regime_switching] = Setting(:regime_switching, true)
m.test_settings[:regime_switching] = Setting(:regime_switching, true)
setup_regime_switching_inds!(m)
@testset "Regime-switching deterministic trends" begin
    dettrend = DSGE.deterministic_trends(m, reg_sys, zeros(8))
    @test @test_matrix_approx_eq dettrend[1] zeros(size(dettrend[1]))
    @test @test_matrix_approx_eq dettrend[2] zeros(size(dettrend[2]))
    @test @test_matrix_approx_eq dettrend[3] zeros(size(dettrend[3]))
end

# Check shock decompositions for regime-switching system
out_shockdec1 = shock_decompositions(m, reg_sys, histshocks,
                                    date_presample_start(m), date_forecast_start(m)) # check we only return the shockdecs for requested period
out_shockdec2 = shock_decompositions(m, reg_sys, histshocks[:, 1:end - 1]) # check default

@testset "Shock decompositions with regime-switching" begin
    @test @test_matrix_approx_eq out_shockdec1[1] states
    @test @test_matrix_approx_eq out_shockdec1[2] obs
    @test @test_matrix_approx_eq out_shockdec1[3] pseudo
    @test @test_matrix_approx_eq out_shockdec2[1] states
    @test @test_matrix_approx_eq out_shockdec2[2] obs
    @test @test_matrix_approx_eq out_shockdec2[3] pseudo
end

## Shock decompositions with time-varying CCC
# Set up
m = Model1002("ss10"; custom_settings = Dict{Symbol, Setting}(:add_altpolicy_pgap => Setting(:add_altpolicy_pgap, true),
                                                              :add_altpolicy_ygap => Setting(:add_altpolicy_ygap, true)))
m <= Setting(:regime_switching, true)
m <= Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m),
                                            2 => Date(2020, 6, 30),
                                            3 => Date(2020, 9, 30),
                                            4 => Date(2020, 12, 31),
                                            5 => Date(2021, 3, 31),
                                            6 => Date(2021, 6, 30)))
m <= Setting(:date_forecast_start, Date(2020, 6, 30))
m <= Setting(:date_conditional_end, Date(2020, 6, 30))
m <= Setting(:tvis_information_set, [1:1, 2:2, 3:6, 4:6, 5:6, 6:6])
m <= Setting(:replace_eqcond, true)
m <= Setting(:gensys2, true)
zlb_rule_eqcond = DSGE.EqcondEntry(DSGE.zlb_rule(), [1., 0.])
m <= Setting(:regime_eqcond_info, Dict(3 => deepcopy(zlb_rule_eqcond),
                                             4 => deepcopy(zlb_rule_eqcond),
                                             5 => deepcopy(zlb_rule_eqcond),
                                             6 => DSGE.EqcondEntry(DSGE.flexible_ait(), [1., 0.])))
m <= Setting(:temporary_altpolicy_names, [:zlb_rule])
setup_regime_switching_inds!(m; cond_type = :full)
df = load(joinpath(path, "..", "reference", "regime_switch_data.jld2"), "regime_switch_df_full")
sys = compute_system(m; tvis = true)
_, histshocks, _, init_states = smooth(m, df, sys; cond_type = :full)
output = DSGE.forecast_one_draw(m, :mode, :full, [:forecastobs, :histpseudo, :forecastpseudo,
                                                  :histstates, :forecaststates],
                                [x.value for x in m.parameters], df, regime_switching = true,
                                n_regimes = get_setting(m, :n_regimes))
shockstates, shockobs, shockpseudo = shock_decompositions(m, sys, histshocks,
                                                          date_mainsample_start(m), date_conditional_end(m), :full)
dettrendstates, dettrendobs, dettrendpseudo = deterministic_trends(m, sys, init_states,
                                                                   date_mainsample_start(m),
                                                                   date_conditional_end(m), :full)
trendstates, trendobs, trendpseudo = trends(m, sys, date_mainsample_start(m), date_conditional_end(m), :full)

states = hcat(output[:histstates], output[:forecaststates])[:, index_shockdec_start(m):end]
obs = hcat(df_to_matrix(m, df; cond_type = :full, include_presample = false)[:, 1:end - 1],
           output[:forecastobs])[:, index_shockdec_start(m):end]
pseudo = hcat(output[:histpseudo], output[:forecastpseudo])[:, index_shockdec_start(m):end]

impl_states = dropdims(sum(shockstates, dims = 3), dims = 3) + dettrendstates + trendstates
impl_obs = dropdims(sum(shockobs, dims = 3), dims = 3) + dettrendobs + trendobs
impl_pseudo = dropdims(sum(shockpseudo, dims = 3), dims = 3) + dettrendpseudo + trendpseudo

@test states ≈ impl_states
@test obs[vcat(1:9, 13), :] ≈ impl_obs[vcat(1:9, 13), :]
@test pseudo ≈ impl_pseudo

df_states = DSGE.prepare_means_table_trend_nostates(m, :full, :state,
                                                    date_mainsample_start(m), date_conditional_end(m), apply_altpolicy = true,
                                                    annualize = false)
df_obs    = DSGE.prepare_means_table_trend_nostates(m, :full, :obs,
                                                    date_mainsample_start(m), date_conditional_end(m), apply_altpolicy = true,
                                                    annualize = false)
df_pseudo = DSGE.prepare_means_table_trend_nostates(m, :full, :pseudo,
                                                    date_mainsample_start(m), date_conditional_end(m), apply_altpolicy = true,
                                                    annualize = false)

@test all(Matrix(df_states[:, 2:end]) .≈ 0.)
@test all((df_to_matrix(m, df_obs) .- sys[n_regimes(sys), :DD]) .≈ 0.)
pseudo_obs = vcat(1:13, 19:n_pseudo_observables(m)) .+ 1 # Remove forward-looking pseudo-obs
@test all((Matrix(df_pseudo[:, pseudo_obs])' .- sys[n_regimes(sys), :DD_pseudo][pseudo_obs .- 1]) .≈ 0.)

nothing


text:
using Pkg
Pkg.activate(pwd())

# # Linear regression
# ## Loading and preparing data

using Plots
using StatsPlots
using RDatasets

iris = dataset("datasets", "iris")
iris[1:5,:]

# ### Exercise:
# We will simplify the goal and estimate the dependence of petal width on petal length. # Create the data $X$ (do not forget to add the bias) and the labels $y$.
#
# Make a graph of the dependence of petal width on petal length.
#
# ---
# ### Solution:

y = iris.PetalWidth
X = hcat(iris.PetalLength, ones(length(y)))

@df iris scatter(
    :PetalLength,
    :PetalWidth;
    label="",
    xlabel = "Petal length",
    ylabel = "Petal width"
)

# ---
#
# ## Training the classifier
# ### Exercise:
# Use the closed-form formula to get the coefficients $w$ for the linear regression. Then
# use the `optim` method derived in the previous lecture to solve the optimization problem
# via gradient descent. The results should be identical.
#
#---
#### Solution:

abstract type Step end

struct GD <: Step
    α::Real
end

optim_step(s::GD, f, g, x) = -s.α*g(x)

function optim(f, g, x, s::Step; max_iter=100)
    for i in 1:max_iter
        x += optim_step(s, f, g, x)
    end
    return x
end

#+

using LinearAlgebra

w = (X'*X) \ (X'*y)
g(w) = X'*(X*w-y)
w2 = optim([], g, zeros(size(X,2)), GD(1e-4); max_iter=10000)

norm(w - w2)

# ---
#
# ### Exercise:
# Write the dependence on the petal width on the petal length. Plot it in the previous
# graph.
#
# ---
# ### Solution:

f_pred(x::Real, w) = w[1]*x + w[2]
x_lims = extrema(iris.PetalLength) .+ [-0.1, 0.1]

@df iris scatter(
    :PetalLength,
    :PetalWidth;
    xlabel = "Petal length",
    ylabel = "Petal width",
    label = "",
    legend = :topleft,
)

plot!(x_lims, x -> f_pred(x,w); label = "Prediction", line = (:black,3))


text:
function jt02()
#
#  test polynomial suggested by Jenkins and Traub
#
p = poly([(1.0:17)...])
    z = [(1.0:17)...]
    z = [z ones(17)]
    p, PolyZeros(z)
end


text:
type VIPData
    F::Array{JuMP.NonlinearExpression,1}
    var::Array{JuMP.Variable,1}
    relation::Dict{JuMP.Variable, JuMP.NonlinearExpression}
end

function VIPModel(solver=Ipopt.IpoptSolver(print_level=0))
    m = Model(solver=solver)
    m.ext[:VIP] = VIPData(Array(JuMP.NonlinearExpression,0), Array(JuMP.Variable,0), Dict{JuMP.Variable, JuMP.NonlinearExpression}() )
    return m
end

function getVIPData(m::Model)
    if haskey(m.ext, :VIP)
        return m.ext[:VIP]::VIPData
    else
        error("The 'getVIPData' function is only for VIP models as in VariationalInequality.jl")
    end
end

################################################################################
# correspond interface
# The most basic/important one. All other interfaces call this method.
function correspond(m::Model, expression::JuMP.NonlinearExpression, variable::JuMP.Variable)
    data = getVIPData(m)
    data.relation[variable] = expression
end
# Alternative
function correspond(m::Model, variable::JuMP.Variable, expression::JuMP.NonlinearExpression)
    correspond(m, expression, variable)
end

# Do we need the below?
# function correspond(m::Model, variables::Array{JuMP.Variable,1}, expressions::Array{JuMP.NonlinearExpression,1})
#     correspond(m, expressions, variables)
# end
# function correspond(m::Model, expressions::Array{JuMP.NonlinearExpression,1}, variables::Array{JuMP.Variable,1})
#     @assert length(expressions) == length(variables)
#     for i in 1:length(variables)
#         correspond(m, expressions[i], variables[i])
#     end
# end

function correspond(m::Model, expressions::Array{JuMP.NonlinearExpression}, variables::Array{JuMP.Variable})
    expressions = collect(expressions)
    variables = collect(variables)
    @assert length(expressions) == length(variables)
    for i in 1:length(variables)
        correspond(m, expressions[i], variables[i])
    end
end
function correspond(m::Model, variables::Array{JuMP.Variable}, expressions::Array{JuMP.NonlinearExpression})
    correspond(m, expressions, variables)
end

# Do we need the below?
# function correspond(m::Model, expressions::JuMP.JuMPArray, variables::JuMP.JuMPArray)
#     variables = collect(variables.innerArray)
#     expressions = collect(expressions.innerArray)
#     @assert length(expressions) == length(variables)
#     for i in 1:length(variables)
#         correspond(m, expressions[i], variables[i])
#     end
# end
################################################################################

# function setVIP(m::Model, expressions::Array{JuMP.NonlinearExpression,1}, variables::Array{JuMP.Variable,1})
#     @assert length(expressions) == length(variables)
#
#     data = getVIPData(m)
#     data.F = expressions
#     data.var = variables
#
#     return m
# end

# function solveVIP(m::Model    ; step_size=0.01,
#                                 algorithm=:fixed_point,
#                                 tolerance=1e-10,
#                                 max_iter=1000            )
#
#     F = getVIPData(m).F
#     var = getVIPData(m).var
#     @assert length(F) == length(var)
#     dim = length(F)
#
#     x0 = zeros(dim)
#
#     return solveVIP(m, x0, step_size=step_size, algorithm=algorithm, tolerance=tolerance, max_iter=max_iter)
#
# end

# Copied from JuMP.jl
# internal method that doesn't print a warning if the value is NaN
_getvalue(v::JuMP.Variable) = v.m.colVal[v.col]

function clearValues(m)
    relation = getVIPData(m).relation
    for variable in keys(relation)
        setvalue(variable, NaN)
    end
end

function initial_projection(m::JuMP.Model)
    relation = getVIPData(m).relation

    initial_values = Dict{JuMP.Variable, Float64}()
    for variable in keys(relation)
        val = _getvalue(variable)
        if isnan(val)
            val = 0.0
        end
        initial_values[variable] = val
    end

    @NLobjective(m, Min, sum{ ( variable - initial_values[variable] )^2, variable in keys(relation)} )
    status = solve(m)
    @assert status==:Optimal
end

function gap_function(m)
    relation = getVIPData(m).relation
    var = getVariables(relation)
    y = getCurrentX(relation)
    Fy = getCurrentF(relation)

    @NLobjective(m, Max,
        sum{ Fy[j] * ( y[j] - var[j] ),    j=1:length(var)}
    )
    solve(m)
    return getobjectivevalue(m)
end

function saveSolution(m)
    relation = getVIPData(m).relation
    var = getVariables(relation)
    x = getCurrentX(relation)
    F = getCurrentF(relation)

    solution = Dict(zip(var, x))
    F_value = Dict(zip(var, F))
    gap = gap_function(m)

    return solution, F_value, gap
end

function solveVIP(m::Model;    step_size=0.01,
                                algorithm=:fixed_point,
                                tolerance=1e-6,
                                max_iter=1000            )

    initial_projection(m)

    if algorithm==:fixed_point
        _fixed_point(m, step_size, tolerance, max_iter)
    elseif algorithm==:extra_gradient
        _extra_gradient(m, step_size, tolerance, max_iter)
    # elseif algorithm==:hyperplane
    #     _hyperplane(m, step_size, tolerance, max_iter)
    end

    # sol, Fval, gap = saveSolution(m)
    return saveSolution(m)
end


text:
# spectra normalization functions, copied from an old version of RvSpectML.EchelleInstruments to here due to the most up-to-date version giving errors
# TODO: push changes to EchelleInstruments to fix these errors, and update RvLineList to use the official version of EchelleInstruments

""" Normalize spectrum based on continuum model from FITS file. """
function continuum_normalize_spectrum!(spectrum::ST) where { ST<:AbstractSpectra }
    @assert haskey(spectrum.metadata,:continuum)
    if spectrum.metadata[:normalization] == :raw
        spectrum.flux ./= spectrum.metadata[:continuum] .* spectrum.metadata[:blaze]
        spectrum.var ./= (spectrum.metadata[:continuum].* spectrum.metadata[:blaze] ) .^2
        spectrum.metadata[:normalization] = :continuum
    elseif spectrum.metadata[:normalization] == :blaze
        spectrum.flux ./= spectrum.metadata[:continuum]
        spectrum.var ./= (spectrum.metadata[:continuum] ) .^2
        spectrum.metadata[:normalization] = :continuum
    elseif spectrum.metadata[:normalization] == :continuum
        # do nothing
    else
        @error "Normalizing from " * string(spectrum.metadata.normalization) * " to continuum isn't implemented yet"
    end
    return spectrum
end

""" Normalize each spectrum based on continuum model from FITS files. """
function continuum_normalize_spectra!(spectra::AS) where { ST<:AbstractSpectra, AS<:AbstractArray{ST} }
    for spectrum in spectra
        continuum_normalize_spectrum!(spectrum)
    end
    return spectra
end

max_stars_repo_name: alexander-wise/RvLineList

max_issues_repo_name: alexander-wise/RvLineList

max_forks_repo_name: alexander-wise/RvLineList


text:
"""
    LogisticLikelihood()

Bernoulli likelihood with a logistic link for the Bernoulli likelihood
```math
    p(y|f) = \\sigma(yf) = \\frac{1}{1+\\exp(-yf)},
```
(for more info see : [wiki page](https://en.wikipedia.org/wiki/Logistic_function))

---

For the analytic version the likelihood, it is augmented via:
```math
    p(y|f,ω) = exp(0.5(yf - (yf)^2 ω))
```
where ``ω ~ PG(ω | 1, 0)``, and `PG` is the Polya-Gamma distribution
See paper : [Efficient Gaussian Process Classification Using Polya-Gamma Data Augmentation](https://arxiv.org/abs/1802.06383)
"""
struct LogisticLikelihood{T<:Real,A<:AbstractVector{T}} <: ClassificationLikelihood{T}
    c::A
    θ::A
    function LogisticLikelihood{T}() where {T<:Real}
        return new{T,Vector{T}}()
    end
    function LogisticLikelihood{T}(c::A, θ::A) where {T<:Real,A<:AbstractVector{T}}
        return new{T,A}(c, θ)
    end
end

function LogisticLikelihood()
    return LogisticLikelihood{Float64}()
end

function implemented(
    ::LogisticLikelihood, ::Union{<:AnalyticVI,<:QuadratureVI,<:GibbsSampling}
)
    return true
end

function init_likelihood(
    ::LogisticLikelihood{T},
    inference::AbstractInference{T},
    nLatent::Int,
    nSamplesUsed::Int,
) where {T}
    if inference isa AnalyticVI || inference isa GibbsSampling
        LogisticLikelihood{T}(abs.(rand(T, nSamplesUsed)), zeros(T, nSamplesUsed))
    else
        LogisticLikelihood{T}()
    end
end

function (::LogisticLikelihood)(y::Real, f::Real)
    return logistic(y * f)
end

function Distributions.loglikelihood(::LogisticLikelihood{T}, y::Real, f::Real) where {T}
    return -log(one(T) + exp(-y * f))
end

function Base.show(io::IO, ::LogisticLikelihood)
    return print(io, "Bernoulli Likelihood with Logistic Link")
end

function compute_proba(l::LogisticLikelihood{T}, f::Real) where {T<:Real}
    return l(1, f)
end

function compute_proba(
    ::LogisticLikelihood{T}, μ::AbstractVector{<:Real}, σ²::AbstractVector{<:Real}
) where {T<:Real}
    N = length(μ)
    pred = zeros(T, N)
    σ²_pred = zeros(T, N)
    for i in 1:N
        x = pred_nodes .* sqrt(max(σ²[i], zero(T))) .+ μ[i]
        pred[i] = dot(pred_weights, logistic.(x))
        σ²_pred[i] = max(dot(pred_weights, logistic.(x) .^ 2) - pred[i]^2, zero(T))
    end
    return pred, σ²_pred
end

### Local Updates Section ###

function local_updates!(
    l::LogisticLikelihood{T}, ::AbstractVector, μ::AbstractVector, diagΣ::AbstractVector
) where {T}
    @. l.c = sqrt(diagΣ + abs2(μ))
    @. l.θ = 0.5 * tanh(0.5 * l.c) / l.c
end

function sample_local!(l::LogisticLikelihood, ::AbstractVector, f::AbstractVector)
    return set_ω!(l, rand.(PolyaGamma.(1, abs.(f))))
end

### Natural Gradient Section ###

∇E_μ(::LogisticLikelihood, ::AOptimizer, y::AbstractVector) = (0.5 * y,)
∇E_Σ(l::LogisticLikelihood, ::AOptimizer, ::AbstractVector) = (0.5 * l.θ,)

### ELBO Section ###

function expec_loglikelihood(
    l::LogisticLikelihood{T},
    ::AnalyticVI,
    y::AbstractVector,
    μ::AbstractVector,
    diag_cov::AbstractVector,
) where {T}
    tot = -(0.5 * length(y) * logtwo)
    tot += 0.5 .* (dot(μ, y) - dot(l.θ, diag_cov) - dot(l.θ, μ))
    return tot
end

AugmentedKL(l::LogisticLikelihood, ::AbstractVector) = PolyaGammaKL(l)

function PolyaGammaKL(l::LogisticLikelihood{T}) where {T}
    return sum(broadcast(PolyaGammaKL, ones(T, length(l.c)), l.c, l.θ))
end

### Gradient Section ###

∇loglikehood(::LogisticLikelihood{T}, y::Real, f::Real) where {T} = y * logistic(-y * f)

function hessloglikehood(::LogisticLikelihood{T}, y::Real, f::Real) where {T<:Real}
    return -exp(y * f) / logistic(-y * f)^2
end


text:
function isdf_error( errmsg, path )  #=

prints an error message given as the first argument "errmsg" and the path to get the place where the error occurs. The path is given by the second argument "path". After the error message and the path have been printed into the stdout, the execution is immediately terminated.

Author: Taku Noda
Started on Nov. 16, 2019
Update on June 29, 2021

Copyright 2021 Central Research Institute of Electric Power Industry

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

=#

println( "" )
println( "iSDF error:" )
println( errmsg )
println( "" )
println( "The error ocurred at:" )
println( path )

# exit( 1 )  # currently commented out to avoid the termination of Julia environment when an error occurs. This should be decommented when the program is released for practical use.

end

max_stars_repo_name: takunoda/omegasim

max_issues_repo_name: takunoda/omegasim

max_forks_repo_name: takunoda/omegasim


text:
"""
`print_map_statistics` prints basic map file properties and characteristics
of road network (OpenStreetMapX.MapData type) built on top of that

**Fields**
* `datapath` : path to map file
* `mapfile` : map file name
* `roadclass` : highest road classe available on the map

"""
function print_map_statistics(datapath::String, mapfile::String, roadclass::Int)
    #Load map
    map_data = get_map_data(datapath, mapfile, use_cache=false; road_levels = Set(1:roadclass));
    #Map file info
    println("Map name: $mapfile")
    println("Size: $(filesize(datapath*mapfile)) B")
    filehash=hash(read(datapath*mapfile))
    println("Map file hash generated with Base.hash function: $(string(filehash,base=16))")
    #Map load statement
    println("\nMap statistics generated for map loaded with statement:")
    println("OpenStretMapX.get_map_data('$(datapath)', '$(mapfile)', use_cache=false; road_levels = Set(1:$(roadclass)))")

    #Bounds
    println("\n---Map bounds (LLA coordinates)---")
    println("MinY = $(map_data.bounds.min_y)");
    println("MaxY = $(map_data.bounds.max_y)");
    println("MinX = $(map_data.bounds.min_x)");
    println("MaxY = $(map_data.bounds.max_y)");
    #Graph
    println("\n---Graph characteristics---")
    print("No. of vertices: "); println(length(map_data.v))
    print("No. of edges: "); println(length(map_data.e))

    #Transportation network
    println("\n---Edges characteristics---") #2176
    get_velocities = OpenStreetMapX.get_velocities(map_data)
    velo_dict = sort(countmap(get_velocities.nzval));
    println("\nNumber of edges with given speed limits:")
    for (k,v) in velo_dict
        println("$(round(k*3.6, digits=2)) km/h: $v edges")
    end

    println("\nNumber of edges with given lanes:") #2176
    roads_lanes = Dict{Int,Int}()
    for r in map_data.roadways
        OpenStreetMapX.haslanes(r) ? lanes = OpenStreetMapX.getlanes(r) : lanes = 1
        roads_lanes[r.id] = lanes
    end
    segments = OpenStreetMapX.find_segments(map_data.nodes, map_data.roadways, map_data.intersections)
    filtered_segments = [[s.node0, s.node1, roads_lanes[s.parent]] for s in segments]
    unique_segments = unique(filtered_segments)
    lanes = [s[3] for s in unique_segments]
    lanes_dict = sort(countmap(lanes))
    for (k,v) in lanes_dict
        if k==1
            println("$k lane: $v edges")
        else
            println("$k lanes: $v edges")
        end
    end

    println("\n---Edges length distribution---") #2176
    edge_len = map_data.w.nzval
    length_quan = quantile(edge_len,0:0.1:1)
    println()
    print(summarystats(edge_len))
    return true
end


text:
# Checks BS runs

println("Running BS tests...")

n = 1000;        # Number of samples
λ = 100;         # Frequencey of changes

########################
# Normal mean segments #
########################
μ, σ = Normal(0.0, 10.0), 1.0
sample, cps = @changepoint_sampler n λ Normal(μ, σ)
seg_costs = NormalMeanSegment(sample);
BS(seg_costs, n);
@BS sample Normal(?,σ)

#######################
# Normal var segments #
#######################
μ, σ = 1.0, Uniform(2.0, 15.0)
sample, cps = @changepoint_sampler n λ Normal(μ, σ)
seg_costs = NormalVarSegment(sample, μ)
BS(seg_costs, n)
@BS sample Normal(μ, ?)

############################
# Exponential changepoints #
############################
μ = Uniform(0.0, 10.0)
sample, cps = @changepoint_sampler n λ Exponential(μ)
seg_costs = ExponentialSegment(sample)
BS(seg_costs, n)
@BS sample Exponential(?)

########################
# Poisson changepoints #
########################
μ = Uniform(0.0, 10.0)
sample, cps = @changepoint_sampler n λ Poisson(μ)
seg_costs = PoissonSegment(sample)
BS(seg_costs, n)
@BS sample Poisson(?)


text:
function KiteGrid=GridToKiteGrid(Grid,OrientFace)
KiteGrid.nBar=[ 0  1   0   -1
                1  0  -1    0];
KiteGrid.Dim=3;
KiteGrid.Type='Quad';

NumNodes=Grid.NumNodes+Grid.NumEdges+Grid.NumFaces;
Nodes(1:NumNodes)=Node([0 0 0],0);
NodeNumber=1;
for iN=1:Grid.NumNodes
  Nodes(NodeNumber)=Node(Grid.Nodes(iN).P,NodeNumber);
  NodeNumber=NodeNumber+1;
end
for iE=1:Grid.NumEdges
  P=0.5*(Grid.Nodes(Grid.Edges(iE).N(1)).P...
    +Grid.Nodes(Grid.Edges(iE).N(2)).P);
  Nodes(NodeNumber)=Node(P,NodeNumber);
  NodeNumber=NodeNumber+1;
end
for iF=1:Grid.NumFaces
  P=[0 0 0];
  for iN=1:size(Grid.Faces(iF).N,2)
    P=P+Grid.Nodes(Grid.Faces(iF).N(iN)).P;
  end
  P=(1/size(Grid.Faces(iF).N,2))*P;
  Nodes(NodeNumber)=Node(P,NodeNumber);
  NodeNumber=NodeNumber+1;
end
KiteGrid.Nodes=Nodes;
KiteGrid.NumNodes=NumNodes;

NumEdges=2*Grid.NumEdges;
for iF=1:Grid.NumFaces
  NumEdges=NumEdges+size(Grid.Faces(iF).N,2);
end

Edges(1:NumEdges)=Edge([1 2],KiteGrid,0,0);
EdgeNumber=1;
for iE=1:Grid.NumEdges
  N1=Grid.Edges(iE).N(1);
  N2=iE+Grid.NumNodes;
  N3=Grid.Edges(iE).N(2);
  Edges(EdgeNumber)=Edge([N1 N2],KiteGrid,EdgeNumber,EdgeNumber);
  EdgeNumber=EdgeNumber+1;
  Edges(EdgeNumber)=Edge([N2 N3],KiteGrid,EdgeNumber,EdgeNumber);
  EdgeNumber=EdgeNumber+1;
end
NumInnerNode=Grid.NumNodes+Grid.NumEdges+1;
for iF=1:Grid.NumFaces
  N2=NumInnerNode;
  for iE=1:size(Grid.Faces(iF).E,2)
    N1=Grid.NumNodes+Grid.Edges(Grid.Faces(iF).E(iE)).E;
    Edges(EdgeNumber)=Edge([N1 N2],KiteGrid,EdgeNumber,EdgeNumber);
    EdgeNumber=EdgeNumber+1;
  end
  NumInnerNode=NumInnerNode+1;
end
KiteGrid.Edges=Edges;
KiteGrid.NumEdges=NumEdges;
KiteGrid.NumEdgesI=NumEdges;
KiteGrid.NumEdgesB=0;
NumFaces=0;
for iF=1:Grid.NumFaces
  NumFaces=NumFaces+size(Grid.Faces(iF).N,2);
end
Faces(1:NumFaces)=Face([0 0 0 0],KiteGrid,0,'x');
FaceNumber=1;
NumInnerEdge=2*Grid.NumEdges+1;
Type='o';
for iF=1:Grid.NumFaces
  for iE=1:size(Grid.Faces(iF).E,2)-1
    if Grid.Edges(Grid.Faces(iF).E(iE)).N(1)...
        ==Grid.Faces(iF).N(iE+1)
      E1=2*Grid.Edges(Grid.Faces(iF).E(iE)).E-1;
    else
      E1=2*Grid.Edges(Grid.Faces(iF).E(iE)).E;
    end
    if Grid.Edges(Grid.Faces(iF).E(iE+1)).N(1)...
        ==Grid.Faces(iF).N(iE+1)
      E2=2*Grid.Edges(Grid.Faces(iF).E(iE+1)).E-1;
    else
      E2=2*Grid.Edges(Grid.Faces(iF).E(iE+1)).E;
    end
    E3=NumInnerEdge+1;
    E4=NumInnerEdge;
    [Faces(FaceNumber),KiteGrid]=Face([E1 E2 E3 E4],KiteGrid,...
      FaceNumber,Type,OrientFace);
    FaceNumber=FaceNumber+1;
    NumInnerEdge=NumInnerEdge+1;
  end
  iE=size(Grid.Faces(iF).E,2);
  if Grid.Edges(Grid.Faces(iF).E(iE)).N(1)...
      ==Grid.Faces(iF).N(1)
    E1=2*Grid.Edges(Grid.Faces(iF).E(iE)).E-1;
  else
    E1=2*Grid.Edges(Grid.Faces(iF).E(iE)).E;
  end
  if Grid.Edges(Grid.Faces(iF).E(1)).N(1)...
      ==Grid.Faces(iF).N(1)
    E2=2*Grid.Edges(Grid.Faces(iF).E(1)).E-1;
  else
    E2=2*Grid.Edges(Grid.Faces(iF).E(1)).E;
  end
  E3=NumInnerEdge+1-size(Grid.Faces(iF).E,2);
  E4=NumInnerEdge;
    [Faces(FaceNumber),KiteGrid]=Face([E1 E2 E3 E4],KiteGrid,...
      FaceNumber,Type,OrientFace);
    FaceNumber=FaceNumber+1;
    NumInnerEdge=NumInnerEdge+1;
end
KiteGrid.Faces=Faces;
KiteGrid.NumFaces=NumFaces;

KiteGrid=Orientation(KiteGrid);
KiteGrid=Renumbering(KiteGrid);
end


text:
using DocStringExtensions
using Test
import Markdown
import LibGit2

include("tests.jl")


text:
"""
    InitialState

A type used for some AbstractTrees.jl iterators to indicate that iteration is in its initial state.
Typically this is used for wrapper types to indicate that the `iterate` function has not yet  been called
on the wrapped object.
"""
struct InitialState end

"""
    TreeCursor{P,N}

Abstract type for tree cursors which when constructed from a node can be used to
navigate the entire tree descended from that node.

Tree cursors satisfy the abstract tree interface with a few additional guarantees:
- Tree cursors all have the [`StoredParents`](@ref) and [`StoredSiblings`](@ref) traits.
- All functions acting on a cursor which returns a tree node is guaranteed to return another `TreeCursor`.
    For example, `children`, `parent` and `nextsiblin` all return a `TreeCursor` of the same type as
    the argument.

Tree nodes which define `children` and have the traits [`StoredParents`](@ref) and
[`StoredSiblings`](@ref) satisfy the `TreeCursor` interface, but calling `TreeCursor(node)` on such
a node wraps them in a [`TrivialCursor`](@ref) to maintain a consistent interface.

## Constructors
All `TreeCursor`s possess (at least) the following constructors
- `T(node)`
- `T(parent, node)`

In the former case the `TreeCursor` is constructed for the tree of which `node` is the root.
"""
abstract type TreeCursor{P,N} end

"""
    nodevaluetype(csr::TreeCursor)

Get the type of the wrapped node.  This should match the return type of [`nodevalue`](@ref).
"""
nodevaluetype(::Type{<:TreeCursor{P,N}}) where {P,N} = N
nodevaluetype(csr::TreeCursor) = nodevaluetype(typeof(csr))

"""
    parenttype(csr::TreeCursor)

The return type of `parent(csr)`.  For properly constructed `TreeCursor`s this is guaranteed to be another
`TreeCursor`.
"""
parenttype(::Type{<:TreeCursor{P,N}}) where {P,N} = P
parenttype(csr::TreeCursor) = parenttype(typeof(csr))

# this is a fallback and may not always be the case
Base.IteratorSize(::Type{<:TreeCursor{P,N}}) where {P,N} = IteratorSize(childtype(N))

Base.length(tc::TreeCursor) = (length ∘ children ∘ nodevalue)(csr)

Base.IteratorEltype(::Type{<:TreeCursor}) = EltypeUnknown()

ParentLinks(::Type{<:TreeCursor}) = StoredParents()

SiblingLinks(::Type{<:TreeCursor}) = StoredSiblings()

# all TreeCursor give children on iteration
children(tc::TreeCursor) = tc

ChildIndexing(tc::TreeCursor) = ChildIndexing(nodevalue(tc))

nodevalue(tc::TreeCursor) = tc.node

parent(tc::TreeCursor) = tc.parent

"""
    TrivialCursor{P,N} <: TreeCursor{P,N}

A [`TreeCursor`](@ref) which matches the functionality of the underlying node.  Tree nodes wrapped by this
cursor themselves have most of the functionality required of a `TreeCursor`, this type exists entirely
for the sake of maintaining a fully consistent interface with other `TreeCursor` objects.
"""
struct TrivialCursor{P,N} <: TreeCursor{P,N}
    parent::P
    node::N
end

parent(csr::TrivialCursor) = parent(csr.node)

TrivialCursor(node) = TrivialCursor(parent(node), node)

function Base.iterate(csr::TrivialCursor, s=InitialState())
    cs = (children ∘ nodevalue)(csr)
    r = s isa InitialState ? iterate(cs) : iterate(cs, s)
    isnothing(r) && return nothing
    (n′, s′) = r
    o = TrivialCursor(n′)
    (o, (o, s′))
end

function nextsibling(csr::TrivialCursor)
    n = nextsibling(csr.node)
    isnothing(n) ? nothing : TrivialCursor(csr.parent, n)
end

function prevsibling(csr::TrivialCursor)
    p = prevsibling(csr.node)
    isnothing(p) ? nothing : TrivialCursor(csr.parent, p)
end

"""
    ImplicitCursor{P,N,S} <: TreeCursor{P,N}

A [`TreeCursor`](@ref) which wraps nodes which cannot efficiently access either their parents or siblings directly.
This should be thought of as a "worst case scenario" tree cursor.  In particular, `ImplicitCursor`s store the
child iteration state of type `S` and for any of `ImplicitCursor`s method to be type-stable it must be possible
to infer the child iteration state type, see [`childstatetype`](@ref).
"""
struct ImplicitCursor{P,N,S} <: TreeCursor{P,N}
    parent::P
    node::N
    sibling_state::S

    ImplicitCursor(p::Union{Nothing,ImplicitCursor}, n, s=InitialState()) = new{typeof(p),typeof(n),typeof(s)}(p, n, s)
end

ImplicitCursor(node) = ImplicitCursor(nothing, node)

Base.IteratorEltype(::Type{<:ImplicitCursor}) = HasEltype()

function Base.eltype(::Type{ImplicitCursor{P,N,S}}) where {P,N,S}
    cst = (childstatetype ∘ nodevalueeltype)(P)
    P′ = ImplicitCursor{P,N,S}
    ImplicitCursor{P′,childtype(N),cst}
end

function Base.eltype(csr::ImplicitCursor)
    cst = (childstatetype ∘ parent ∘ nodevalue)(csr)
    ImplicitCursor{typeof(csr),childtype(nodevalue(csr)),cst}
end

function Base.iterate(csr::ImplicitCursor, s=InitialState())
    cs = (children ∘ nodevalue)(csr)
    # do NOT just write an iterate(x, ::InitialState) method, it's an ambiguity nightmare
    r = s isa InitialState ? iterate(cs) : iterate(cs, s)
    isnothing(r) && return nothing
    (n′, s′) = r
    o = ImplicitCursor(csr, n′, s′)
    (o, s′)
end

function nextsibling(csr::ImplicitCursor)
    isroot(csr) && return nothing
    cs = (children ∘ nodevalue ∘ parent)(csr)
    # do NOT just write an iterate(x, ::InitialState) method, it's an ambiguity nightmare
    r = csr.sibling_state isa InitialState ? iterate(cs) : iterate(cs, csr.sibling_state)
    isnothing(r) && return nothing
    (n′, s′) = r
    ImplicitCursor(parent(csr), n′, s′)
end

"""
    IndexedCursor{P,N} <: TreeCursor{P,N}

A [`TreeCursor`](@ref) for tree nodes with the [`IndexedChildren`](@ref) trait but for which parents and siblings
are not directly accessible.

This type is very similar to [`ImplicitCursor`](@ref) except that it is free to assume that the child iteration
state is an integer starting at `1` which drastially simplifies type inference and slightly simplifies the
iteration methods.
"""
struct IndexedCursor{P,N} <: TreeCursor{P,N}
    parent::P
    node::N
    index::Int

    IndexedCursor(p::Union{Nothing,IndexedCursor}, n, idx::Integer=1) = new{typeof(p),typeof(n)}(p, n, idx)
end

IndexedCursor(node) = IndexedCursor(nothing, node)

Base.IteratorSize(::Type{<:IndexedCursor}) = HasLength()

function Base.eltype(::Type{IndexedCursor{P,N}}) where {P,N}
    P′ = IndexedCursor{P,N}
    IndexedCursor{P′,childtype(N)}
end
Base.eltype(csr::IndexedCursor) = IndexedCursor{typeof(csr),childtype(nodevalue(csr))}
Base.length(csr::IndexedCursor) = (length ∘ children ∘ nodevalue)(csr)

function Base.getindex(csr::IndexedCursor, idx)
    cs = (children ∘ nodevalue)(csr)
    IndexedCursor(csr, cs[idx], idx)
end

function Base.iterate(csr::IndexedCursor, idx=1)
    idx > length(csr) && return nothing
    (csr[idx], idx+1)
end

function nextsibling(csr::IndexedCursor)
    p = parent(csr)
    isnothing(p) && return nothing
    idx = csr.index + 1
    idx > length(p) && return nothing
    p[idx]
end

function prevsibling(csr::IndexedCursor)
    idx = csr.index - 1
    idx < 1 && return nothing
    parent(csr)[idx]
end

"""
    SiblingCursor{P,N} <: TreeCursor{P,N}

A [`TreeCursor`](@ref) for trees with the [`StoredSiblings`](@ref) trait.
"""
struct SiblingCursor{P,N} <: TreeCursor{P,N}
    parent::P
    node::N

    SiblingCursor(p::Union{Nothing,SiblingCursor}, n) = new{typeof(p),typeof(n)}(p, n)
end

SiblingCursor(node) = SiblingCursor(nothing, node)

Base.IteratorSize(::Type{SiblingCursor{P,N}}) where {P,N} = IteratorSize(childtype(N))

Base.IteratorEltype(::Type{<:SiblingCursor}) = HasEltype()

function Base.eltype(::Type{SiblingCursor{P,N}}) where {P,N}
    cst = (childstatetype ∘ nodevaluetype)(P)
    P′ = SiblingCursor{P,N}
    SiblingCursor{P′,childtype(N)}
end

Base.eltype(csr::SiblingCursor) = SiblingCursor{typeof(csr),childtype(nodevalue(csr))}

function Base.iterate(csr::SiblingCursor, (c, s)=(nothing, InitialState()))
    cs = (children ∘ nodevalue)(csr)
    r = s isa InitialState ? iterate(cs) : iterate(cs, s)
    isnothing(r) && return nothing
    (n′, s′) = r
    o = SiblingCursor(csr, n′)
    (o, (o, s′))
end

function nextsibling(csr::SiblingCursor)
    n = nextsibling(nodevalue(csr))
    isnothing(n) ? nothing : SiblingCursor(parent(csr), n)
end

function prevsibling(csr::SiblingCursor)
    p = prevsibling(nodevalue(csr))
    isnothing(n) ? nothing : SiblingCursor(parent(csr), p)
end

TreeCursor(::ChildIndexing, ::StoredParents, ::StoredSiblings, node) = TrivialCursor(node)

TreeCursor(::ChildIndexing, ::ImplicitParents, ::StoredSiblings, node) = SiblingCursor(node)

TreeCursor(::NonIndexedChildren, ::ParentLinks, ::ImplicitSiblings, node) = ImplicitCursor(node)
TreeCursor(::IndexedChildren, ::ParentLinks, ::ImplicitSiblings, node) = IndexedCursor(node)

TreeCursor(node) = TreeCursor(ChildIndexing(node), ParentLinks(node), SiblingLinks(node), node)


text:
abstract AbstractLotkaVolterra

#abstract
"""The per capita growth rate of each speicies in the ecosystem"""
specificrate(Z::AbstractLotkaVolterra,x::Vector{Float64}) = intrinsicrate(Z)+communitymatrix(Z)*x
"""The total growth rate of each species in the ecosystem"""
rate(Z::AbstractLotkaVolterra,x::Vector{Float64}) = x.*specificrate(Z,x)

"""The Jacobian matrix of the rate function at the specified density x"""
jacobian(Z::AbstractLotkaVolterra,x::Vector{Float64}) = diagm(specificrate(Z,x))+ diagm(x)*communitymatrix(Z)
"""The Jacobian matrix of the rate function at the interior fixed point"""
jacobian(Z::AbstractLotkaVolterra) = Diagonal(fixedpoint(Z))*communitymatrix(Z)
"""The jacobian matrix of the rate function at the fixed point of the specified sector"""
function jacobian(Z::AbstractLotkaVolterra,sector::Vector{Int64})
    x = fixedpoint(Z,sector)
    diagterm = specificrate(Z,copy(x))
    diagterm[sector] = 0.
    squareterm = Diagonal(x)*communitymatrix(Z)
    return Diagonal(diagterm) + squareterm
end

#sector check and generation functions
"""The set of all sector vectors in the ecosystem with a specified number of species (dim)"""
sectors(Z::AbstractLotkaVolterra,dim::Int64) = subsets(collect(1:nspecies(Z)),dim)
"""The set of all sector vectors in the ecosystem"""
sectors(Z::AbstractLotkaVolterra) = chain(Any[Int64[]],[sectors(Z,dim) for dim = 1:nspecies(Z)]...)
"""The set of all sector vectors of a given kind.
The supported kinds include: :feasible,:infeasible,:stable,:unstable,:viable,:unviable  """
function sectors(Z::AbstractLotkaVolterra,kind::Symbol)
  allsectors = filter(sector->isfixedpoint(Z,sector),sectors(Z))#this needs to be pulled out, just here for now
  @match kind begin
    :feasible => filter(sector->isfeasible(Z,sector),allsectors)
    :infeasible => filter(sector->!isfeasible(Z,sector),allsectors)
    :stable => filter(sector->isstable(Z,sector),allsectors)
    :unstable => filter(sector->!isstable(Z,sector),allsectors)
    :viable => filter(sector->isviable(Z,sector),allsectors)
    :unviable => filter(sector->!isviable(Z,sector),allsectors)
    _ => error("not a valid sector filter kind")
  end
end
"""Tests whether a given vector is a valid sector for the ecosystem"""
function issector(Z::AbstractLotkaVolterra,sector::Vector{Int64})
  N = length(sector)
  usector = unique(sector)
  return length(usector) == N && all(0.< usector .<= nspecies(Z))
end

#fixed point checks and generation
"""Checks for the existence of a fixed point in the interior sector of an ecosystem.
The test may fail due to numerical illposedness of the community matrix"""
function isfixedpoint(Z::AbstractLotkaVolterra)
  try
    fixedpoint(Z)
    return true
  catch
    return false
  end
end
"""Checks for the existence of a fixed point in the specified sector of an ecosystem.
The test may fail due to numerical illposedness of the community matrix"""
function isfixedpoint(Z::AbstractLotkaVolterra,sector::Vector{Int64})
  try
    fixedpoint(Z,sector)
    return true
  catch
    return false
  end
end
"""Returns the interior fixed point of the ecosystem or throws an error if it does nto exist"""
fixedpoint(Z::AbstractLotkaVolterra) = -communitymatrix(Z)\intrinsicrate(Z)
"""Returns the fixed point of the ecosystem in the specified sector or throws an error if it does to exist"""
function fixedpoint(Z::AbstractLotkaVolterra,sector::Vector{Int64})
    A = communitymatrix(Z)[sector,sector]
    r = intrinsicrate(Z)[sector]
    x = zeros(nspecies(Z))
    x[sector] = -A\r
    return x
end
"""Returns all the fixed points of an ecosystem"""
fixedpoints(Z::AbstractLotkaVolterra) = imap(sector->fixedpoint(Z,sector),filter(sector->isfixedpoint(Z,sector),sectors(Z)))
"""Returns all the fixed points of an ecosystem of a specified kind"""
function fixedpoints(Z::AbstractLotkaVolterra,kind::Symbol)
  return imap(sector->fixedpoint(Z,sector),filter(sector->isfixedpoint(Z,sector),sectors(Z,kind)))
end

isfeasible(Z::AbstractLotkaVolterra,sector::Vector{Int64}) = length(sector)>0 && all(fixedpoint(Z,sector)[sector] .> 0.)
isfeasible(Z::AbstractLotkaVolterra) = all(fixedpoint(Z) .> 0.)

isstable(Z::AbstractLotkaVolterra) = real(eigmaxreal(jacobian(Z))) < 0.
isstable(Z::AbstractLotkaVolterra,sector::Vector{Int64}) = real(eigmaxreal(jacobian(Z,sector))) < 0.

isviable(Z::AbstractLotkaVolterra) = isfeasible(Z) & isstable(Z)
isviable(Z::AbstractLotkaVolterra,sector::Vector{Int64}) = isfeasible(Z,sector) & isstable(Z,sector)

feasibility(Z::AbstractLotkaVolterra) = isfeasible(Z)? 1 : 0
stability(Z::AbstractLotkaVolterra) = isstable(Z) < 0 ? 1 : 0
viability(Z::AbstractLotkaVolterra) = isviable(Z) ? 1 : 0

"""Solve the initial value problem for the ecosystem desnity of the ecosystem ODE. API subject to change."""
function odeint(Z::AbstractLotkaVolterra,x0::Vector{Float64},times)
    f(t::Float64,x::Vector{Float64}) = rate(Z,x)
    return ode45(f,x0,times)
end


text:
using CorrNoise
using Test

# RNG tests

ref = [0.9854981268, 0.2536861135,
       0.8791850018, 0.8541532028,
       0.4161281283, 0.6491611481,
       0.8336713058, 0.2541507778,
       0.0673349826, 0.7994472017,
       0.2600378725, 0.5134736516,
       0.4544143085, 0.9879117254,
       0.1875971474, 0.1698732623,
       0.9844831470, 0.2436698098,
       0.2365035815, 0.9442295933]

flatrng = initflatrng128()
for (idx, sample) in enumerate([rand(flatrng) for i = 1:length(ref)])
    @assert ref[idx] ≈ sample
end

ref = [1.8051588123, -1.0150233503,
       -0.5138818361, 0.6974503387,
       0.0959640650, -1.7090943033,
       0.2825155379, -0.0263954840,
       -0.4486580098, -0.4245704031,
       -1.6567426564, -1.2158240845,
       0.3230971021, -1.0730251963,
       -0.6426576020, 1.0481948256,
       -0.0527193353, -0.2020275077,
       0.4885079587, 0.2111337896]

normrng = GaussRNG(initflatrng128())
for (idx, sample) in enumerate([randn(normrng) for i = 1:length(ref)])
    @assert ref[idx] ≈ sample
end

ref = [2.0886370365, -0.6074844495,
       -0.3464683356, 0.8936789422,
       0.4167745672, -1.6416294545,
       0.1259486772, -0.1427304943,
       -0.6395859847, -0.7526144784,
       -2.3116079175, -2.3217436788,
       -0.9229344940, -2.4367338428,
       -2.2756948201, -0.5210397131,
       -1.4655131788, -1.6547241628,
       -0.9190954646, -1.0864977062]

normrng = GaussRNG(initflatrng128())
oof2rng = Oof2RNG(normrng, 1.15e-5, 0.05, 1.0)
for (idx, sample) in enumerate([randoof2(oof2rng) for i = 1:length(ref)])
    @assert ref[idx] ≈ sample
end

ref = [2.0400958027, -0.6813881047,
       -0.3848671621, 0.84863765610,
       0.34838982440, -1.6703162833,
       0.13588215630, -0.1362109745,
       -0.6176813736, -0.7037434318,
       -2.2015645097, -2.1249009577,
       -0.6866971017, -2.1640479726,
       -1.9400718236, -0.1768193407,
       -1.1308790625, -1.2995275170,
       -0.5582322798, -0.7326164277]

normrng = GaussRNG(initflatrng128())
oofrng = OofRNG(normrng, -1.7, 1.15e-5, 0.05, 1.0)
for (idx, sample) in enumerate([randoof(oofrng) for i = 1:length(ref)])
    @assert ref[idx] ≈ sample
end