text:
@doc raw"""
    gausslaguerre(n::Integer) -> x, w  # nodes, weights

Return nodes `x` and weights `w` of [Gauss-Laguerre quadrature](https://en.wikipedia.org/wiki/Gauss%E2%80%93Laguerre_quadrature).

```math
\int_{0}^{+\infty} f(x) e^{-x} dx \approx \sum_{i=1}^{n} w_i f(x_i)
```

# Examples
```jldoctest
julia> x, w = gausslaguerre(3);

julia> f(x) = x^4;

julia> I = dot(w, f.(x));

julia> I ≈ 24
true
```
"""
function gausslaguerre(n::Integer)
    return gausslaguerre(n, 0.0)
end

@doc raw"""
    gausslaguerre(n::Integer, α::Real) -> x, w  # nodes, weights

Return nodes `x` and weights `w` of generalized [Gauss-Laguerre quadrature](https://en.wikipedia.org/wiki/Gauss%E2%80%93Laguerre_quadrature).

```math
\int_{0}^{+\infty} f(x) x^\alpha e^{-x} dx \approx \sum_{i=1}^{n} w_i f(x_i)
```

# Examples
```jldoctest
julia> x, w = gausslaguerre(3, 1.0);

julia> f(x) = x^4;

julia> I = dot(w, f.(x));

julia> I ≈ 120
true
```

Optionally, a reduced quadrature rule can be computed. In that case, only those
points and weights are computed for which the weight does not underflow in the
floating point precision type. Supply the optional argument `reduced = true`.

Though the code is generic, heuristical choices on the choice of the algorithm
are based on achieving machine precision accuracy only for `Float64` type. In
case the default choice of algorithm does not achieve the desired accuracy, the
user can manually invoke the following routines:
- `gausslaguerre_GW`: computation based on Golub-Welsch
- `gausslaguerre_rec`: computation based on Newton iterations applied to evaluation
   using the recurrence relation
- `gausslaguerre_asy`: the asymptotic expansions
"""
function gausslaguerre(n::Integer, α::Real; reduced = false)
    if α ≤ -1
        throw(DomainError(α, "The Laguerre parameter α ≤ -1 corresponds to a nonintegrable weight function"))
    end
    if n < 0
        throw(DomainError(n, "gausslaguerre($n,$α) not defined: n must be positive."))
    end

    # Guess the numerical type from the supplied type of α
    # Although the code is generic, the heuristics are derived for Float64 precision
    T = typeof(float(α))
    if n == 0
        T[],T[]
    elseif n == 1
        [1+α], [gamma(1+α)]
    elseif n == 2
        [α + 2-sqrt(α+2),α+2+sqrt(α+2)],
        [((α-sqrt(α+2)+2)*gamma(α+2))/(2*(α+2)*(sqrt(α+2)-1)^2),
         ((α+sqrt(α+2)+2)*gamma(α+2))/(2*(α+2)*(sqrt(α+2)+1)^2)]
    elseif n < 15
        # Use Golub-Welsch for small n
        gausslaguerre_GW(n, α)
    elseif n < 128
        # Use the recurrence relation for moderate n
        gausslaguerre_rec(n, α)
    else
        # Use explicit asymptotic expansions for larger n
        # The restriction to α comes from the restriction on ν in besselroots
        if α < 5
            gausslaguerre_asy(n, α, reduced=reduced, T=-1, recompute=true)
        else
            gausslaguerre_rec(n, α)
        end
    end
end

# Our threshold for deciding on underflow
underflow_threshold(x) = underflow_threshold(typeof(x))
underflow_threshold(::Type{T}) where {T <: AbstractFloat} = 10floatmin(T)

"""
Compute the Gauss-Laguerre rule using explicit asymptotic expansions for the nodes and weights.
Optional parameters are:
- `reduced`: compute a reduced quadrature rule, discarding all points and weights as soon as the weights underflow
- `T`: the order of the expansion. Set `T=-1` to determine the order adaptively depending on the size of the terms in the expansion
- `recompute`: if a crude measure of the error is larger than a tolerance, the point and weight are recomputed using the (slower) recursion+newton approach, yielding more reliable accurate results.
"""
function gausslaguerre_asy(n::Integer, α;
    reduced = false,
    T = max(1, ceil(Int, 50/log(n))),  # Heuristic for number of terms
    recompute = false)

    if α^2/n > 1
        @warn "A large value of α may lead to inaccurate results."
    end

    ELT = typeof(float(α))

    n_alloc = reduced ? 0 : n
    x = zeros(ELT, n_alloc)
    w = zeros(ELT, n_alloc)

    # The expansions are given in powers of 1/(4n+2α+2)
    d = one(ELT)/(4n+2α+2)

    # Heuristical indices for Bessel and Airy regions
    k_bessel = max(ceil(Int, sqrt(n) ), 7)
    k_airy = floor(Int, 0.9*n)

    # The Bessel region
    # First compute the roots of the Bessel function of order α
    jak_vector = approx_besselroots(α, k_bessel)

    bessel_wins = true
    k = 0
    while bessel_wins && k < n
        k += 1
        # We iterate until the estimated error of the bulk expansion is smaller
        # than the one of the Bessel expansion
        jak = (k < k_bessel) ? jak_vector[k] : McMahon(α, k)

        xk, wk, δ_bessel = gausslaguerre_asy_bessel(n, α, jak, d, T)
        xkb, wkb, δ_bulk = gausslaguerre_asy_bulk(n, α, k, d, T)
        if δ_bulk < δ_bessel
            bessel_wins = false
            xk = xkb
            wk = wkb
        end
        if recompute
            δ = min(δ_bessel,δ_bulk)
            if δ > 1e-13
                xk_rec, wk_rec = gl_rec_newton(xk, n, α)
                if abs(xk_rec-xk) < 100δ
                    xk = xk_rec
                    wk = wk_rec
                end
            end
        end
        if reduced
            if abs(wk) < underflow_threshold(ELT)
                return x, w
            else
                push!(x, xk); push!(w, wk)
            end
        else
            x[k] = xk; w[k] = wk
        end
    end

    # The bulk region
    # - First we go from where we left of to our heuristic
    while k < k_airy-1
        k += 1
        xk, wk, δ_bulk = gausslaguerre_asy_bulk(n, α, k, d, T)
        if recompute
            if δ_bulk > 1e-13
                xk_rec, wk_rec = gl_rec_newton(xk, n, α)
                if abs(xk_rec-xk) < 100δ_bulk
                    xk = xk_rec
                    wk = wk_rec
                end
            end
        end
        if reduced
            if abs(wk) < underflow_threshold(ELT)
                return x, w
            else
                push!(x, xk); push!(w, wk)
            end
        else
            x[k] = xk; w[k] = wk
        end
    end

    # - Then we compare to Airy until it wins, and then we switch to just Airy
    bulk_wins = true
    while bulk_wins && k < n
        k += 1
        xk, wk, δ_bulk = gausslaguerre_asy_bulk(n, α, k, d, T)
        xka, wka, δ_airy = gausslaguerre_asy_airy(n, α, k, d, T)
        if δ_airy < δ_bulk
            bulk_wins = false
            xk = xka
            wk = wka
        end
        if recompute
            δ = min(δ_airy,δ_bulk)
            if δ > 1e-13
                xk_rec, wk_rec = gl_rec_newton(xk, n, α)
                if abs(xk_rec-xk) < 100δ
                    xk = xk_rec
                    wk = wk_rec
                end
            end
        end
        if reduced
            if abs(wk) < underflow_threshold(ELT)
                return x, w
            else
                push!(x, xk); push!(w, wk)
            end
        else
            x[k] = xk; w[k] = wk
        end
    end

    # The Airy region
    while k < n
        k += 1
        xk, wk, δ_airy = gausslaguerre_asy_airy(n, α, k, d, T)
        if recompute
            if δ_airy > 1e-13
                xk_rec, wk_rec = gl_rec_newton(xk, n, α)
                if abs(xk_rec-xk) < 100δ_airy
                    xk = xk_rec
                    wk = wk_rec
                end
            end
        end
        if reduced
            if abs(wk) < underflow_threshold(ELT)
                return x, w
            else
                push!(x, xk); push!(w, wk)
            end
        else
            x[k] = xk; w[k] = wk
        end
    end

    # Sanity check
    if ( minimum(x) < 0.0 ) || ( maximum(x) > 4*n + 2*α + 2 ) ||  ( minimum(diff(x)) ≤ 0.0 ) || (minimum(w) < 0.0)
        @warn "Unexpected inconsistency in the computation of nodes and weights"
    end

    return x, w
end

## Expansion coefficients
# These are explicit formulas of the coefficients, up to a simple postprocessing
# that is common to all factors and not included here (see below).
#
# General expressions are given in terms of α, more specific expressions
# follow for the special case α = 0.

## The bulk

# Note: there is always one division by an integer, placed such that it preserves the type of `d`
gl_bulk_x3(t, d, α) = -(12*α^2 + 5*(1-t)^(-2) - 4*(1-t)^(-1) - 4) * d / 12
gl_bulk_x5(t, d, α) = d^3*(1-t)/t/720*(1600*(1-t)^(-6) - 3815*(1-t)^(-5)
    + 480*α^4 +2814*(1-t)^(-4) - 576*(1-t)^(-3) - 960*α^2 - 48*(15*α^4 - 30*α^2 + 7)*(1-t)^(-1)
    -16*(1-t)^(-2) + 224)
gl_bulk_x7(t, d, α) = -d^5/181440*(1-t)^2/t^2*(10797500*(1-t)^(-10) - 43122800*(1-t)^(-9)
    + 66424575*(1-t)^(-8) -48469876*(1-t)^(-7) + 193536*α^6 + 16131880*(1-t)^(-6)
    + 80*(315*α^4 - 630*α^2 -221)*(1-t)^(-4) - 1727136*(1-t)^(-5)
    - 967680*α^4 - 320*(63*α^4 - 126*α^2 +43)*(1-t)^(-3)
    + 384*(945*α^6 - 4620*α^4 + 6405*α^2 - 1346)*(1-t)^(-2)
    + 1354752*α^2   - 23040*(21*α^6 - 105*α^4 + 147*α^2 - 31)*(1-t)^(-1) -285696)
gl_bulk_x9(t, d, α) = d^7/10886400*(1-t)^3/t^3*(43222750000*(1-t)^(-14) - 241928673000*(1-t)^(-13)
    + 566519158800*(1-t)^(-12) -714465642135*(1-t)^(-11) + 518401904799*(1-t)^(-10)
    + 672*(12000*α^4 - 24000*α^2 +64957561)*(1-t)^(-8) - 212307298152*(1-t)^(-9)
    + 24883200*α^8 - 192*(103425*α^4 -206850*α^2 + 15948182)*(1-t)^(-7)
    + 3360*(4521*α^4 - 9042*α^2 - 7823)*(1-t)^(-6) -232243200*α^6
    - 1792*(3375*α^6 - 13905*α^4  + 17685*α^2 - 1598)*(1-t)^(-5)
    + 16128*(450*α^6 - 2155*α^4 + 2960*α^2 - 641)*(1-t)^(-4)
    + 812851200*α^4 -768*(70875*α^8 - 631260*α^6 + 2163630*α^4
    - 2716980*α^2 +555239)*(1-t)^(-3)  + 768*(143325*α^8 - 1324260*α^6
    + 4613070*α^4 -5826660*α^2 + 1193053)*(1-t)^(-2) - 1028505600*α^2
    - 5806080*(15*α^8 -140*α^6 + 490*α^4 - 620*α^2 + 127)*(1-t)^(-1) + 210677760)

gl_bulk_w3(t, d, α) = d^2/6*(2*t + 3)/(t-1)^3
gl_bulk_w5(t, d, α) = (1-t)^2/720/t^2*d^4*(8000*(1-t)^(-8) - 24860*(1-t)^(-7) + 27517*(1-t)^(-6)
    - 12408*(1-t)^(-5) + 1712*(1-t)^(-4) +16*(15*α^4 - 30*α^2 + 7)*(1-t)^(-2) + 32*(1-t)^(-3))
gl_bulk_w7(t, d, α) = -(1-t)^3/90720/t^3*d^6*(43190000*(1-t)^(-12) -204917300*(1-t)^(-11)
    + 393326325*(1-t)^(-10) - 386872990*(1-t)^(-9) + 201908326*(1-t)^(-8)
    + 80*(315*α^4 - 630*α^2 + 53752)*(1-t)^(-6)  - 50986344*(1-t)^(-7)
    - 320*(189*α^4 -378*α^2 - 89)*(1-t)^(-5) + 480*(63*α^4 - 126*α^2
    + 43)*(1-t)^(-4)  -384*(315*α^6 - 1470*α^4 + 1995*α^2 - 416)*(1-t)^(-3)
    + 2304*(21*α^6 -105*α^4 + 147*α^2 - 31)*(1-t)^(-2) )

# And for α = 0
gl_bulk_x3(t, d) = -d/12*(5*(1-t)^(-2) - 4*(1-t)^(-1) - 4)
gl_bulk_x5(t, d) = d^3*(1-t)/t/720*(1600*(1-t)^(-6) - 3815*(1-t)^(-5) + 2814*(1-t)^(-4)
    - 576*(1-t)^(-3) - 48*7*(1-t)^(-1) -16*(1-t)^(-2) + 224)
gl_bulk_x7(t, d) = -d^5/181440*(1-t)^2/t^2*(10797500*(1-t)^(-10)
    - 43122800*(1-t)^(-9) + 66424575*(1-t)^(-8) -48469876*(1-t)^(-7)
    + 16131880*(1-t)^(-6) - 80*221*(1-t)^(-4) - 1727136*(1-t)^(-5) - 320*43*(1-t)^(-3)
    - 384*1346*(1-t)^(-2) + 23040*31*(1-t)^(-1) -285696)
gl_bulk_x9(t, d) = d^7/10886400*(1-t)^3/t^3*(43222750000*(1-t)^(-14)
    - 241928673000*(1-t)^(-13) + 566519158800*(1-t)^(-12) -714465642135*(1-t)^(-11)
    + 518401904799*(1-t)^(-10) + 672*64957561*(1-t)^(-8)   - 212307298152*(1-t)^(-9)
    - 192*15948182*(1-t)^(-7)  - 3360*7823*(1-t)^(-6) + 1792*1598*(1-t)^(-5)
    + 16128*(- 641)*(1-t)^(-4)  -768*555239*(1-t)^(-3)  + 768*1193053*(1-t)^(-2)
    - 5806080*127*(1-t)^(-1) + 210677760)

gl_bulk_w3(t, d) = d^2/6*(2*t + 3)/(t-1)^3
gl_bulk_w5(t, d) = (1-t)^2/720/t^2*d^4*(8000*(1-t)^(-8) - 24860*(1-t)^(-7)
    + 27517*(1-t)^(-6) - 12408*(1-t)^(-5) + 1712*(1-t)^(-4) +16*7*(1-t)^(-2)
    + 32*(1-t)^(-3))
gl_bulk_w7(t, d) = -(1-t)^3/90720/t^3*d^6*(43190000*(1-t)^(-12)
    - 204917300*(1-t)^(-11) + 393326325*(1-t)^(-10) - 386872990*(1-t)^(-9)
    + 201908326*(1-t)^(-8) +80*53752*(1-t)^(-6)
    - 50986344*(1-t)^(-7) + 320*89*(1-t)^(-5)
    + 480*43*(1-t)^(-4) + 384*416*(1-t)^(-3) - 2304*31*(1-t)^(-2) )

## The hard edge (Bessel region)

gl_bessel_x3(jak, d, α) = (jak^2 + 2*α^2 - 2)*d^2 / 3
gl_bessel_x5(jak, d, α) = (11*jak^4 +3*jak^2*(11*α^2-19) +46*α^4 -140*α^2 +94)*d^4 / 45
gl_bessel_x7(jak, d, α) = (657*jak^6 +36*jak^4*(73*α^2-181) +2*jak^2*(2459*α^4 -10750*α^2 +14051)
    + 4*(1493*α^6 -9303*α^4 +19887*α^2 - 12077) )*d^6 / 2835
gl_bessel_x9(jak, d, α) = (10644*jak^8 + 60*(887*α^2 - 2879)*jak^6 + (125671*α^4 -729422*α^2 + 1456807)*jak^4
    + 3*(63299*α^6 - 507801*α^4 + 1678761*α^2 - 2201939)*jak^2 + 2*(107959*α^8
    - 1146220*α^6 + 5095482*α^4 -10087180*α^2 + 6029959) )*d^8 / 42525

gl_bessel_w3(jak, d, α) = (α^2 + jak^2 -1)*2*d^2 / 3
gl_bessel_w5(jak, d, α) = (46*α^4 + 33*jak^4 +6*jak^2*(11*α^2 -19) -140*α^2 +94)*d^4 / 45
gl_bessel_w7(jak, d, α) = (1493*α^6 + 657*jak^6 + 27*(73*α^2 - 181)*jak^4 - 9303*α^4
    + (2459*α^4 -10750*α^2 + 14051)*jak^2 + 19887*α^2 - 12077)*4*d^6 / 2835
gl_bessel_w9(jak, d, α) = (215918*α^8 + 53220*jak^8 + 240*(887*α^2 - 2879)*jak^6 -2292440*α^6 +
    3*(125671*α^4 - 729422*α^2 + 1456807)*jak^4 + 10190964*α^4  +
    6*(63299*α^6 - 507801*α^4 + 1678761*α^2 -2201939)*jak^2 -
    20174360*α^2 + 12059918)*d^8 / 42525

# And for α = 0:
gl_bessel_x3(jak, d) = (jak^2 - 2)*d^2 / 3
gl_bessel_x5(jak, d) = (11*jak^4 - 57*jak^2 + 94)d^4 / 45
gl_bessel_x7(jak, d) = (657*jak^6 - 6516*jak^4 + 28102*jak^2 - 48308)*d^6 / 2835
gl_bessel_x9(jak, d) = (10644*jak^8 - 172740*jak^6 + 1456807*jak^4 -  6605817*jak^2 + 12059918)*d^8 / 42525
gl_bessel_x11(jak, d) = (410649*jak^10 -  9908262*jak^8 + 138902061*jak^6 - 1248722004*jak^4 + 6028914206*jak^2 - 11427291076)*d^10 / 1403325

gl_bessel_w3(jak, d) = (jak^2 - 1)*2*d^2 / 3
gl_bessel_w5(jak, d) = (33*jak^4 -114*jak^2 + 94)*d^4 / 45
gl_bessel_w7(jak, d) = (657*jak^6 - 4887*jak^4 + 14051*jak^2 - 12077)*4*d^6 / 2835
gl_bessel_w9(jak, d) = (53220*jak^8 - 690960*jak^6 + 4370421*jak^4 - 13211634*jak^2 + 12059918)*d^8 / 42525
gl_bessel_w11(jak, d) = (1231947*jak^10 - 24770655*jak^8 + 277804122*jak^6 - 1873083006*jak^4 + 6028914206*jak^2 - 5713645538)*2*d^10 / 1403325

## The soft edge (Airy region)

gl_airy_x1(ak, d, α) = 1/d + ak*(d/4)^(-1/3)
gl_airy_x3(ak, d, α) = ak^2*(d*16)^(1/3)/5 + (11/35-α^2-12/175*ak^3)*d + (16/1575*ak+92/7875*ak^4)*2^(2/3)*d^(5/3)
gl_airy_x5(ak, d, α) = -(15152/3031875*ak^5+1088/121275*ak^2)*2^(1/3)*d^(7/3)

gl_airy_x1(ak, d) = 1/d + ak*(d/4)^(-1/3)
gl_airy_x3(ak, d) = ak^2*(d*16)^(1/3)/5 + (11/35-12/175*ak^3)*d + (16/1575*ak+92/7875*ak^4)*2^(2/3)*d^(5/3)
gl_airy_x5(ak, d) = -(15152/3031875*ak^5+1088/121275*ak^2)*2^(1/3)*d^(7/3)

# Sum the first Q elements of vals, and return the absolute value of the next
# value in the list (or of the last value in the list)
function sum_explicit(vals, Q)
    T = eltype(vals[1])
    z = zero(T)
    for q = min(Q,length(vals)):-1:1
        z += vals[q]
    end
    if Q < length(vals)
        delta = abs(vals[Q+1])
    else
        delta = abs(vals[end])
    end
    z, delta
end

function sum_adaptive(vals)
    z = vals[1]
    i = 1
    while (i < length(vals)) && (abs(vals[i+1]) < abs(vals[i]))
        i += 1
        z += vals[i]
    end
    delta = abs(vals[min(i+1,length(vals))])
    z, delta
end

function gl_bulk_solve_t(n, k, d)
    T = typeof(d)
    pt = (4n-4k+3)*d
    t = T(π)^2/16*(pt-1)^2
    diff = 100
    iter = 0
    maxiter = 20
    while (abs(diff) > 100eps(T)) && (iter < maxiter)
        iter += 1
        diff = (pt*π +2*sqrt(t-t^2) -acos(2*t-1) )*sqrt(t/(1-t))/2
        t -= diff
    end
    if iter == maxiter
        @warn "Maximal number of iterations reached in the computation of t for the bulk"
    end
    t
end

function gausslaguerre_asy_bulk(n, α, k, d, T)
    if α == 0
        return gausslaguerre_asy0_bulk(n, k, d, T)
    end

    t = gl_bulk_solve_t(n, k, d)
    x3 = gl_bulk_x3(t, d, α)
    x5 = gl_bulk_x5(t, d, α)
    x7 = gl_bulk_x7(t, d, α)
    x9 = gl_bulk_x9(t, d, α)
    w3 = gl_bulk_w3(t, d, α)
    w5 = gl_bulk_w5(t, d, α)
    w7 = gl_bulk_w7(t, d, α)

    xs = (x3, x5, x7, x9)
    ws = (w3, w5, w7)

    xk, xdelta = (T > 0) ? sum_explicit(xs, (T-1)>>1) : sum_adaptive(xs)
    wk, wdelta = (T > 0) ? sum_explicit(ws, (T-1)>>1) : sum_adaptive(ws)

    xk += t/d

    wfactor = xk^α * exp(-xk) * 2π * sqrt(t/(1-t))
    wk = wfactor * (1+wk)
    wdelta *= wfactor

    xk, wk, max(xdelta,wdelta)
end

function gausslaguerre_asy0_bulk(n, k, d, T)
    t = gl_bulk_solve_t(n, k, d)
    x3 = gl_bulk_x3(t, d)
    x5 = gl_bulk_x5(t, d)
    x7 = gl_bulk_x7(t, d)
    x9 = gl_bulk_x9(t, d)
    w3 = gl_bulk_w3(t, d)
    w5 = gl_bulk_w5(t, d)
    w7 = gl_bulk_w7(t, d)

    xs = (x3, x5, x7, x9)
    ws = (w3, w5, w7)

    xk, xdelta = (T > 0) ? sum_explicit(xs, (T-1)>>1) : sum_adaptive(xs)
    wk, wdelta = (T > 0) ? sum_explicit(ws, (T-1)>>1) : sum_adaptive(ws)

    xk += t/d

    wfactor = exp(-xk) * 2π * sqrt(t/(1-t))
    wk = wfactor * (1+wk)
    wdelta *= wfactor

    xk, wk, max(xdelta,wdelta)
end

function gausslaguerre_asy_bessel(n, α, jak, d, T)
    if α == 0
        return gausslaguerre_asy0_bessel(n, jak, d, T)
    end
    x3 = gl_bessel_x3(jak, d, α)
    x5 = gl_bessel_x5(jak, d, α)
    x7 = gl_bessel_x7(jak, d, α)
    x9 = gl_bessel_x9(jak, d, α)
    w3 = gl_bessel_w3(jak, d, α)
    w5 = gl_bessel_w5(jak, d, α)
    w7 = gl_bessel_w7(jak, d, α)
    w9 = gl_bessel_w9(jak, d, α)

    xs = (x3, x5, x7, x9)
    ws = (w3, w5, w7, w9)

    xk, xdelta = (T > 0) ? sum_explicit(xs, (T-1)>>1) : sum_adaptive(xs)
    wk, wdelta = (T > 0) ? sum_explicit(ws, (T-1)>>1) : sum_adaptive(ws)

    xfactor = jak^2 * d
    xk = xfactor * (1 + xk)
    xdelta *= xfactor

    # Invoking the besselj function below is the cause of memory
    # allocation of this routine
    wfactor = 4d * xk^α * exp(-xk) / besselj(α-1, jak)^2
    wk = wfactor * (1 + wk)
    wdelta *= wfactor

    return xk, wk, max(xdelta,wdelta)
end

function gausslaguerre_asy0_bessel(n, jak, d, T)
    x3 = gl_bessel_x3(jak, d)
    x5 = gl_bessel_x5(jak, d)
    x7 = gl_bessel_x7(jak, d)
    x9 = gl_bessel_x9(jak, d)
    x11 = gl_bessel_x11(jak, d)
    w3 = gl_bessel_w3(jak, d)
    w5 = gl_bessel_w5(jak, d)
    w7 = gl_bessel_w7(jak, d)
    w9 = gl_bessel_w9(jak, d)
    w11 = gl_bessel_w11(jak, d)

    xs = (x3, x5, x7, x9, x11)
    ws = (w3, w5, w7, w9, w11)

    xk, xdelta = (T > 0) ? sum_explicit(xs, (T-1)>>1) : sum_adaptive(xs)
    wk, wdelta = (T > 0) ? sum_explicit(ws, (T-1)>>1) : sum_adaptive(ws)

    xfactor = jak^2 * d
    xk = xfactor * (1 + xk)
    xdelta *= xfactor

    wfactor = 4d * exp(-xk) / besselj(-1, jak)^2
    wk = wfactor * (1 + wk)
    wdelta *= wfactor

    return xk, wk, max(xdelta,wdelta)
end

function compute_airy_root(n, k)
    index = n-k+1
    if index ≤ 11
        ak = AIRY_ROOTS[index]
    else
        t = 3 * π/2 * (index-0.25)
        ak = -t^(2/3)*(1 + 5/48/t^2 - 5/36/t^4 + 77125/82944/t^6 -10856875/6967296/t^8)
    end
    ak
end

function gausslaguerre_asy_airy(n, α, k, d, T)
    if α == 0
        return gausslaguerre_asy0_airy(n, k, d, T)
    end

    ak = compute_airy_root(n, k)
    x1 = gl_airy_x1(ak, d, α)
    x3 = gl_airy_x3(ak, d, α)
    x5 = gl_airy_x5(ak, d, α)

    xs = (x1, x3, x5)

    xk, xdelta = (T > 0) ? sum_explicit(xs, (T+1)>>1) : sum_adaptive(xs)

    wk = 4^(1/3)*xk^(α+1/3)*exp(-xk)/(airyaiprime(ak))^2
    wdelta = abs(wk)

    return xk, wk, max(xdelta,wdelta)
end

function gausslaguerre_asy0_airy(n, k, d, T)
    ak = compute_airy_root(n, k)

    x1 = gl_airy_x1(ak, d)
    x3 = gl_airy_x3(ak, d)
    x5 = gl_airy_x5(ak, d)

    xs = (x1, x3, x5)

    xk, xdelta = (T > 0) ? sum_explicit(xs, (T+1)>>1) : sum_adaptive(xs)

    wk = 4^(1/3) * xk^(1/3) * exp(-xk) / (airyaiprime(ak))^2
    wdelta = abs(wk)

    return xk, wk, max(xdelta,wdelta)
end

"""
Calculate Gauss-Laguerre nodes and weights from the eigenvalue decomposition of
the Jacobi matrix.
"""
function gausslaguerre_GW(n, α)
    _α = 2*(1:n) .+ (α-1)  # 3-term recurrence coeffs a and b
    β = sqrt.( (1:n-1).*(α .+ (1:n-1)) )
    T = SymTridiagonal(promote(collect(_α), β)...)
    x, V = eigen(T)  # eigenvalue decomposition
    w = gamma(α+1)*V[1,:].^2  # Quadrature weights
    return x, vec(w)
end

########################## Routines for the forward recurrence ##########################

function gl_rec_newton(x0, n, α; maxiter = 20, computeweight = true)
    T = eltype(x0)
    step = x0
    iter = 0
    xk = x0

    xk_prev = xk
    pn_prev = floatmax(T)
    pn_deriv = zero(T)
    while (abs(step) > 40eps(T)*xk) && (iter < maxiter)
        iter += 1
        pn, pn_deriv = evalLaguerreRec(n, α, xk)
        if abs(pn) ≥ abs(pn_prev)*(1-50eps(T))
            # The function values do not decrease enough any more due to roundoff errors.
            xk = xk_prev # Set to the previous value and quit.
            break
        end
        step = pn / pn_deriv
        xk_prev = xk
        xk -= step
        pn_prev = pn
    end
    if ( xk < 0 ) || ( xk > 4n + 2α + 2 ) || ( iter == maxiter )
        @warn "Newton method may not have converged in gausslaguerre_rec($n,$α)."
    end
    wk = oftype(xk, 0)
    if computeweight
        pn_min1, ~ = evalLaguerreRec(n-1, α, xk)
        wk = (n^2 +α*n)^(-1/2)/pn_min1/pn_deriv
    end

    return xk, wk
end

"Compute Gauss-Laguerre rule based on the recurrence relation, using Newton iterations on an initial guess."
function gausslaguerre_rec(n, α; reduced = false)
    T = typeof(float(α))

    n_alloc = reduced ? 0 : n
    w = zeros(T, n_alloc)
    x = zeros(T, n_alloc)

    # We compute up to 7 starting values for the Newton iterations
    n_pre = min(n, 7)

    ν = 4n + 2α + 2
    x_pre = T.(approx_besselroots(α, n_pre)).^2 / ν # this is a lower bound by [DLMF 18.16.10]

    noUnderflow = true  # this flag turns false once the weights start to underflow
    for k in 1:n
        # Use sextic extrapolation for a new initial guess
        xk = (k ≤ n_pre) ? x_pre[k] : 7*x[k-1] -21*x[k-2] +35*x[k-3] -35*x[k-4] +21*x[k-5] -7*x[k-6] +x[k-7]

        xk, wk = gl_rec_newton(xk, n, α, maxiter = 20, computeweight = noUnderflow)
        if noUnderflow && abs(wk) < underflow_threshold(T)
            noUnderflow = false
        end

        if reduced
            if !noUnderflow
                return x, w
            else
                push!(x, xk); push!(w, wk)
            end
        else
            x[k] = xk
            w[k] = wk
        end
    end

    return x, w
end

"""
Evaluate the orthonormal associated Laguerre polynomial with positive leading coefficient,
as well as its derivative, in the point x using the recurrence relation.
"""
function evalLaguerreRec(n, α, x)
    T = typeof(α)
    pnprev = zero(T)
    pn = 1/sqrt(gamma(α+1))
    pndprev = zero(T)
    pnd = zero(T)
    for k in 1:n
        pnold = pn
        pn = (x -2*k -α+1)/sqrt(k*(α+k))*pnold-sqrt((k-1+α)*(k-1)/k/(k+α))*pnprev
        pnprev = pnold
        pndold = pnd
        pnd = (pnold+(x-2*k-α+1)*pndold)/sqrt(k*(α+k)) -sqrt((k-1+α)*(k-1)/k/(α+k))*pndprev
        pndprev = pndold
    end

    return pn, pnd
end


text:
module GeometricIntegrators

    using Reexport

    @reexport using GeometricBase
    @reexport using GeometricBase.Config
    @reexport using GeometricBase.Utils
    @reexport using GeometricEquations

    include("Solutions.jl")
    @reexport using .Solutions
    include("Discontinuities.jl")
    @reexport using .Discontinuities
    include("Integrators.jl")
    @reexport using .Integrators
    include("Simulations.jl")
    @reexport using .Simulations
    include("Tableaus.jl")
    @reexport using .Tableaus

    function __init__()
        add_config(:verbosity, 1)
    end

end

ext: jl


text:
function HyperDiffusionVec!(FUCG,v1CG,v2CG,wCG,ThCG,RhoCG,CG,Param)
OP=CG.OrdPoly+1;
NF=Param.Grid.NumFaces;
nz=Param.Grid.nz;
FRTemp1 = Param.CacheC5
FRTemp2 = Param.CacheC6
FRotCurl2VecDSS!(FRTemp1,FRTemp2,v1CG,v2CG,CG,Param);
FRotCurl2Vec!(view(FUCG,:,:,:,:,Param.uPos:Param.vPos),FRTemp1,FRTemp2,CG,Param)

FGradDiv2VecDSS!(FRTemp1,FRTemp2,v1CG,v2CG,CG,Param);
FGradDiv2Vec!(view(FUCG,:,:,:,:,Param.uPos:Param.vPos),FRTemp1,FRTemp2,CG,Param)

# FTempW=FDivGrad2Vec(0.5*(W(:,1:nz)+W(:,2:nz+1)),CG,Param);
# FTempW=-Param.HyperDDiv*FDivGrad2Vec(FTempW,CG,Param);
#FW(:,2:nz)=0.5*(FTempW(:,1:nz-1)+FTempW(:,2:nz));
#
# FTempT=FDivGrad2VecDSS(ThCG./RhoCG,CG,Param);
# FUCG(:,:,:,:,4)=-Param.HyperDDiv*FDivRhoGrad2Vec(FTempT,RhoCG,CG,Param);
if Param.HyperDDiv>0
  FDivGrad2VecDSS!(FRTemp1,ThCG./RhoCG,CG,Param)
  FDivRhoGrad2Vec!(view(FUCG,:,:,:,:,Param.ThPos),FRTemp1,RhoCG,CG,Param)
end
end

function HyperDiffusionVec(v1CG,v2CG,wCG,ThCG,RhoCG,CG,Param)
OP=CG.OrdPoly+1;
NF=Param.Grid.NumFaces;
nz=Param.Grid.nz;
HyperDCurl = Param.HyperDCurl
HyperDGrad = Param.HyperDGrad
HyperDDiv = Param.HyperDDiv
FUCG=zeros(OP,OP,NF,nz,4);
(FRTemp1,FRTemp2)=FRotCurl2VecDSS(v1CG,v2CG,CG,Param);
(FGTemp1,FGTemp2)=FGradDiv2VecDSS(v1CG,v2CG,CG,Param);
#FRTemp1=FRTemp1-FGTemp1;
#FRTemp2=FRTemp2-FGTemp2;
FUCG[:,:,:,:,1:2]=FUCG[:,:,:,:,1:2]-HyperDCurl*FRotCurl2Vec(FRTemp1,FRTemp2,CG,Param);
FUCG[:,:,:,:,1:2]=FUCG[:,:,:,:,1:2]+
  HyperDGrad*FGradDiv2Vec(-FGTemp1,-FGTemp2,CG,Param);

# FTempW=FDivGrad2Vec(0.5*(W(:,1:nz)+W(:,2:nz+1)),CG,Param);
# FTempW=-Param.HyperDDiv*FDivGrad2Vec(FTempW,CG,Param);
#FW(:,2:nz)=0.5*(FTempW(:,1:nz-1)+FTempW(:,2:nz));
#
# FTempT=FDivGrad2VecDSS(ThCG./RhoCG,CG,Param);
# FUCG(:,:,:,:,4)=-Param.HyperDDiv*FDivRhoGrad2Vec(FTempT,RhoCG,CG,Param);
if HyperDDiv>0
  FTempT=FDivGrad2VecDSS(ThCG./RhoCG,CG,Param);
  FUCG[:,:,:,:,4]=-HyperDDiv*FDivRhoGrad2Vec(FTempT,RhoCG,CG,Param);
end
return FUCG
end


text:
"""
# script test_illusions

- Julia version: 1.4
- Author: niallcullinane
- Date: 2020-08-20

Script to run set of images with LAMINART and plot layers and activation.
"""

using DrWatson
@quickactivate "MEngProject"
using MEngProject,
    CUDA,
    DifferentialEquations,
    PyPlot,
    NNlib,
    ImageFiltering,
    Images,
    ImageIO,
    MEngProject.LaminartKernels,
    MEngProject.LaminartInitFunc,
    MEngProject.Utils,
    BenchmarkTools,
    Test

using OrdinaryDiffEq,
    ParameterizedFunctions, LSODA, Sundials, DiffEqDevTools, Noise

batch = 1002

# files = readdir(datadir("img"))

files = [
    "Iine_gap_1_100_gs.png",
    "Iine_gap_2_100_gs.png",
    "Iine_gap_3_100_gs.png",
    "Iine_gap_4_100_gs.png",
    "Iines_gaps_100_gs.png",
    "diag_dots_100_gs.png",
    "diag_gap_100_gs.png",
    "kan_sq_cont.png",
    "kan_sq_cont_l.png",
    #      "mo05709.png",
    "stairs_100gs.png",
    "stairs_200gs.png",
    "viper00187.png",
    "viper00661.jpg",
    "viper00715.png",
    "viper00717.png",
    "viper00720.png",
    "viper00721.png",
    "viper00842.jpg",
    "viper00891.jpg",
    "viper00904.jpg",
    "viper00921.jpg",
    "viper01006.jpg",
    "viper01333.jpg",
]

batch_ = string(batch, "_", rand(1000:9999))
mkdir(plotsdir(string("illusions", batch_)))
mkdir(plotsdir(string("illusions_nocb", batch_)))

tspan = (0.0f0, 800f0)

for file in files
    try
        p = LaminartInitFunc.parameterInit_conv_gpu(
            datadir("img", file),
            Parameters.parameters_f32,
        )

        u0 = cu(reshape(
            zeros(Float32, p.dim_i, p.dim_j * (5 * p.K + 2)),
            p.dim_i,
            p.dim_j,
            5 * p.K + 2,
            1,
        ))

        arr1 = similar(@view u0[:, :, 1:2, :])
        arr2 = similar(@view u0[:, :, 1:1, :])

        f = LaminartFunc.LamFunction(
            arr1, #x
            similar(arr1), #m
            similar(arr1), #s
            arr2, #x_lgn,
            similar(arr1), #C,
            similar(arr1), #H_z,
            similar(arr1), # dy_temp,
            similar(arr1), # dm_temp,
            similar(arr1), # dz_temp,
            similar(arr1), # ds_temp,
            similar(arr2), # dv_temp,
            similar(arr1), # H_z_temp,
            similar(arr2), #  V_temp_1,
            similar(arr2), #  V_temp_2,
            similar(arr1), #  Q_temp,
            similar(arr1), #   P_temp
        )

        prob = ODEProblem(f, u0, tspan, p)
        # 	@benchmark sol = solve(prob)
        sol = solve(prob, BS3())

        # plots
        for t ∈ [25, 50, 100, 200, 400, 800]
            Utils.plot_k2_vsep1(sol, t, "illusions", batch_, file)
            Utils.plot_k2_vsep1(
                sol,
                t,
                "illusions_nocb",
                batch_,
                file,
                cb = false,
            )
        end

        Utils.plot_t_act(sol, "illusions", batch_, file)
    catch err
        print(err)
    finally
        u0 = nothing
        p = nothing
        arr1 = nothing
        arr2 = nothing
        f = nothing
        prob = nothing
        sol = nothing
        CUDA.reclaim()
    end
end


text:
using ReTest, Random, AdvancedHMC, ForwardDiff
include("common.jl")

using OrdinaryDiffEq
using LinearAlgebra: dot
using Statistics: mean

@testset "Integrator" begin
    ϵ = 0.01
    lf = Leapfrog(ϵ)

    θ_init = randn(D)
    h = Hamiltonian(UnitEuclideanMetric(D), ℓπ, ∂ℓπ∂θ)
    r_init = AdvancedHMC.rand(h.metric)

    n_steps = 10

    @testset "step" begin
        z = AdvancedHMC.phasepoint(h, copy(θ_init), copy(r_init))
        z_step_loop = z

        t_step_loop = @elapsed for i = 1:n_steps
            z_step_loop = AdvancedHMC.step(lf, h, z_step_loop)
        end

        t_step = @elapsed z_step = AdvancedHMC.step(lf, h, z, n_steps)

        @info "Performance of loop of step() v.s. step()" n_steps t_step_loop t_step t_step_loop / t_step

        @test z_step_loop.θ ≈ z_step.θ atol=DETATOL
        @test z_step_loop.r ≈ z_step.r atol=DETATOL
    end

    @testset "jitter" begin
        @testset "Leapfrog" begin
            ϵ0 = 0.1
            lf = Leapfrog(ϵ0)
            @test lf.ϵ == ϵ0
            @test AdvancedHMC.nom_step_size(lf) == ϵ0
            @test AdvancedHMC.step_size(lf) == ϵ0

            lf2 = AdvancedHMC.jitter(Random.GLOBAL_RNG, lf)
            @test lf2 === lf
            @test AdvancedHMC.nom_step_size(lf2) == ϵ0
            @test AdvancedHMC.step_size(lf2) == ϵ0
        end

        @testset "JitteredLeapfrog" begin
            ϵ0 = 0.1
            lf = JitteredLeapfrog(ϵ0, 0.5)
            @test lf.ϵ0 == ϵ0
            @test AdvancedHMC.nom_step_size(lf) == ϵ0
            @test lf.ϵ == ϵ0
            @test AdvancedHMC.step_size(lf) == lf.ϵ

            lf2 = AdvancedHMC.jitter(Random.GLOBAL_RNG, lf)
            @test lf2.ϵ0 == ϵ0
            @test AdvancedHMC.nom_step_size(lf2) == ϵ0
            @test lf2.ϵ != ϵ0
            @test AdvancedHMC.step_size(lf2) == lf2.ϵ
        end
    end

    @testset "update_nom_step_size" begin
        @testset "Leapfrog" begin
            ϵ0 = 0.1
            lf = Leapfrog(ϵ0)
            @test AdvancedHMC.nom_step_size(lf) == ϵ0

            lf2 = AdvancedHMC.update_nom_step_size(lf, 0.5)
            @test lf2 !== lf
            @test AdvancedHMC.nom_step_size(lf2) == 0.5
            @test AdvancedHMC.step_size(lf2) == 0.5
        end

        @testset "JitteredLeapfrog" begin
            ϵ0 = 0.1
            lf = JitteredLeapfrog(ϵ0, 0.5)
            @test AdvancedHMC.nom_step_size(lf) == ϵ0

            lf2 = AdvancedHMC.update_nom_step_size(lf, 0.2)
            @test lf2 !== lf
            @test AdvancedHMC.nom_step_size(lf2) == 0.2
            # check that we've only updated nominal step size
            @test AdvancedHMC.step_size(lf2) == ϵ0
        end
    end

    @testset "temper" begin
        αsqrt = 2.0
        lf = TemperedLeapfrog(ϵ, αsqrt ^ 2)
        r = ones(5)
        r1 = AdvancedHMC.temper(lf, r, (i=1, is_half=true), 3)
        r2 = AdvancedHMC.temper(lf, r, (i=1, is_half=false), 3)
        r3 = AdvancedHMC.temper(lf, r, (i=2, is_half=true), 3)
        r4 = AdvancedHMC.temper(lf, r, (i=2, is_half=false), 3)
        r5 = AdvancedHMC.temper(lf, r, (i=3, is_half=true), 3)
        r6 = AdvancedHMC.temper(lf, r, (i=3, is_half=false), 3)
        @test r1 == αsqrt * ones(5)
        @test r2 == αsqrt * ones(5)
        @test r3 == αsqrt * ones(5)
        @test r4 == inv(αsqrt) * ones(5)
        @test r5 == inv(αsqrt) * ones(5)
        @test r6 == inv(αsqrt) * ones(5)
        @test_throws BoundsError AdvancedHMC.temper(lf, r, (i=4, is_half=false), 3)

    end

    @testset "Analyitcal solution to Eq (2.11) of Neal (2011)" begin
        negU = q -> -dot(q, q) / 2

        ϵ = 0.01
        for lf in [
            Leapfrog(ϵ),
            DiffEqIntegrator(ϵ, VerletLeapfrog())
        ]
            q_init = randn(1)
            h = Hamiltonian(UnitEuclideanMetric(1), negU, ForwardDiff)
            p_init = AdvancedHMC.rand(h.metric)

            q, p = copy(q_init), copy(p_init)
            z = AdvancedHMC.phasepoint(h, q, p)

            n_steps = 10_000
            qs = zeros(n_steps)
            ps = zeros(n_steps)
            Hs = zeros(n_steps)
            for i = 1:n_steps
                z = AdvancedHMC.step(lf, h, z)
                qs[i] = z.θ[1]
                ps[i] = z.r[1]
                Hs[i] = -AdvancedHMC.neg_energy(z)
            end

            # Throw first 1_000 steps
            qs = qs[1_000:end]
            ps = ps[1_000:end]
            Hs = Hs[1_000:end]

            # Check if all points located at a cirle centered at the origin
            rs = sqrt.(qs.^2 + ps.^2)
            @test all(x-> abs(x - mean(rs)) < 2e-3, rs)

            # Check if the Hamiltonian energy is stable
            @test all(x-> abs(x - mean(Hs)) < 2e-3, Hs)
        end
    end

end


text:
using Random
using LightGraphs
using BlossomV

using ChpSim
include("BiMaps.jl"); using .BiMaps

export x_gate!, z_gate!, reset!, measure_reset!, rz_gate!,
    NoiseModel, SyndromeCircuit, BasicSyndrome,
    CodeDistanceSim, CodeDistanceRun, LogicalOpSim, LogicalOpRun,
    MatchingGraphWeights, apply_sim_error!, do_single_run, do_n_runs,
    do_single_logical_op_run, do_n_logical_op_runs

const rng = Random.GLOBAL_RNG

function Base.append!(dict::Dict, iter)
    for (k, v) in iter
        dict[k] = v
    end
end

# Extra ChpSim methods
function x_gate!(state::ChpState, qubit::Int)
    hadamard!(state, qubit)
    phase!(state, qubit)
    phase!(state, qubit)
    hadamard!(state, qubit)
end
function z_gate!(state::ChpState, qubit::Int)
    phase!(state, qubit)
    phase!(state, qubit)
end
function reset!(state::ChpState, qubit::Int)
    if measure!(state, qubit, bias=0).value
        x_gate!(state, qubit)
    end
    nothing
end
function measure_reset!(state::ChpState, qubit::Int)
    meas = measure!(state, qubit)
    if meas.value
        x_gate!(state, qubit)
    end
    meas
end

# New non-cliford gate, returns a list of sampled bits of length dist.
# rz is translated as a stochastic Z flips along the logical operators
function rz_gate!(state::ChpState, theta::Float64, dist::Int, z_qubits::Vector{Int})
    r = rand(rng, Float64, dist)
    pI = cos(theta/2)^2
    pZ = sin(theta/2)^2
    projected=BitArray(undef, dist)
    for i in 1:dist
        if r[i] < pI
            projected[i]=0
        end
        if r[i] >= pI
            projected[i]=1
            z_gate!(state, z_qubits[i])
        end
    end
    pI, pZ, projected
end

"""
    NoiseModel(...)

Holds probabilities of various quantum error sources
"""
struct NoiseModel
    errors::Dict{Symbol, Float64}
end
function Base.getproperty(noise::NoiseModel, s::Symbol)
    s == :errors && return getfield(noise, :errors)
    getfield(noise, :errors)[s]
end

"""
    SyndromeCircuit

Parent type representing a style of syndrome measurement circuit.
"""
abstract type SyndromeCircuit end
"""
    BasicSyndrome

Basic syndrome measurement with simple errors
"""
struct BasicSyndrome <: SyndromeCircuit end

NodeT = Tuple{Int, Symbol, Int, Int}
PlaqInfoT = NamedTuple{(:ancilla, :data), Tuple{Int, Vector{Int}}}

"""
    CodeDistanceSim(z_dist[, x_dist, m_dist], syndrome_circuit, noise_model)

Stores the configuration and precomputed data used for error correction.
"""
struct CodeDistanceSim
    z_dist::Int
    x_dist::Int
    m_dist::Int
    syndrome_circuit::SyndromeCircuit
    noise_model::NoiseModel
    num_qubits::Int
    anc_qubits::Vector{Int}
    data_qubits::Vector{Int}
    z_plaqs::Vector{Tuple{Int, Int}}
    z_plaq_info::Vector{PlaqInfoT}
    z_space_boundary::Vector{Tuple{Int, Int}}
    z_doubled_boundary::Vector{Tuple{Int, Int}}
    z_graph_nodes::BiMap{NodeT, Int}
    z_graph::Graph{Int}
    z_costs::Matrix{Float64}
    z_bpaths::Set{Tuple{Int, Int}}
    z_path_lengths::Dict{Tuple{Int, Int}, Int}
    x_plaqs::Vector{Tuple{Int, Int}}
    x_plaq_info::Vector{PlaqInfoT}
    x_space_boundary::Vector{Tuple{Int, Int}}
    x_doubled_boundary::Vector{Tuple{Int, Int}}
    x_graph_nodes::BiMap{NodeT, Int}
    x_graph::Graph{Int}
    x_costs::Matrix{Float64}
    x_bpaths::Set{Tuple{Int, Int}}
    x_path_lengths::Dict{Tuple{Int, Int}, Int}
end

function CodeDistanceSim(dist::Int, syndrome_circuit::SyndromeCircuit,
                         noise_model::NoiseModel)
    CodeDistanceSim(dist, dist, dist, syndrome_circuit, noise_model)
end
function CodeDistanceSim(z_dist::Int, x_dist::Int, m_dist::Int,
                         syndrome_circuit::SyndromeCircuit, noise_model::NoiseModel)
    # Graphs
    (z_plaqs, z_space_boundary, z_doubled_boundary), (
        x_plaqs, x_space_boundary, x_doubled_boundary) = (
        make_plaqs(z_dist, x_dist))
    num_qubits, anc_qubits, data_qubits, z_plaq_info, x_plaq_info, _, _, _ = (
        make_qubit_assignments(z_dist, z_plaqs, x_dist, x_plaqs))
    z_graph_nodes, z_graph = (
        make_graph(m_dist+1, z_plaqs, z_space_boundary, false, false))
    x_graph_nodes, x_graph = (
        make_graph(m_dist+1, x_plaqs, x_space_boundary, false, false))

    z_costs, z_bpaths, z_path_lengths = constuct_graph_costs(
        z_graph_nodes, z_graph, z_doubled_boundary,
        noise_model, syndrome_circuit, false, z_dist, m_dist)
    x_costs, x_bpaths, x_path_lengths = constuct_graph_costs(
        x_graph_nodes, x_graph, x_doubled_boundary,
        noise_model, syndrome_circuit, true, x_dist, m_dist)

    CodeDistanceSim(
        z_dist, x_dist, m_dist, syndrome_circuit, noise_model,
        num_qubits, anc_qubits, data_qubits,
        z_plaqs, z_plaq_info, z_space_boundary, z_doubled_boundary,
            z_graph_nodes, z_graph, z_costs, z_bpaths, z_path_lengths,
        x_plaqs, x_plaq_info, x_space_boundary, x_doubled_boundary,
            x_graph_nodes, x_graph, x_costs, x_bpaths, x_path_lengths,
    )
end

#####################

struct LogicalOpSim
    z_dist::Int
    x_dist::Int
    m_dist::Int
    theta::Float64
    syndrome_circuit::SyndromeCircuit
    noise_model::NoiseModel
    num_qubits::Int
    anc_qubits::Vector{Int}
    z_anc_qubits::Vector{Int}
    x_anc_qubits::Vector{Int}
    data_qubits::Vector{Int}
    logical_z_qubits::Vector{Int}
    diagonal_x_ancilla::Vector{Int}
    #logical_op_projection::AbstractArray{Bool}
    z_plaqs::Vector{Tuple{Int, Int}}
    z_plaq_info::Vector{PlaqInfoT}
    z_space_boundary::Vector{Tuple{Int, Int}}
    z_doubled_boundary::Vector{Tuple{Int, Int}}
    z_graph_nodes::BiMap{NodeT, Int}
    z_graph::Graph{Int}
    z_costs::Matrix{Float64}
    z_bpaths::Set{Tuple{Int, Int}}
    z_path_lengths::Dict{Tuple{Int, Int}, Int}
    x_plaqs::Vector{Tuple{Int, Int}}
    x_plaq_info::Vector{PlaqInfoT}
    x_space_boundary::Vector{Tuple{Int, Int}}
    x_doubled_boundary::Vector{Tuple{Int, Int}}
    x_graph_nodes::BiMap{NodeT, Int}
    x_graph::Graph{Int}
    x_costs::Matrix{Float64}
    x_bpaths::Set{Tuple{Int, Int}}
    x_path_lengths::Dict{Tuple{Int, Int}, Int}
end

function LogicalOpSim(dist::Int, theta::Float64, syndrome_circuit::SyndromeCircuit,
                         noise_model::NoiseModel)
    LogicalOpSim(dist, dist, dist, theta, syndrome_circuit, noise_model)
end
function LogicalOpSim(z_dist::Int, x_dist::Int, m_dist::Int, theta::Float64,
                         syndrome_circuit::SyndromeCircuit, noise_model::NoiseModel)
    # Graphs
    (z_plaqs, z_space_boundary, z_doubled_boundary), (
        x_plaqs, x_space_boundary, x_doubled_boundary) = (
        make_plaqs(z_dist, x_dist))
    # Qubits
    num_qubits, anc_qubits, data_qubits, z_plaq_info, x_plaq_info, lz_qubits, z_anc_qubits, x_anc_qubits = (
        make_qubit_assignments(z_dist, z_plaqs, x_dist, x_plaqs))
    z_graph_nodes, z_graph = (
        make_graph(m_dist+1, z_plaqs, z_space_boundary, false, false))
    x_graph_nodes, x_graph = (
        make_graph(m_dist+1, x_plaqs, x_space_boundary, false, false))

    z_costs, z_bpaths, z_path_lengths = constuct_graph_costs(
        z_graph_nodes, z_graph, z_doubled_boundary,
        noise_model, syndrome_circuit, false, z_dist, m_dist)
    x_costs, x_bpaths, x_path_lengths = constuct_graph_costs(
        x_graph_nodes, x_graph, x_doubled_boundary,
        noise_model, syndrome_circuit, true, x_dist, m_dist)
    #lz_proj = BitArray(undef, z_dist) # dummy value for now

    diagonal_anc = diagonal_x_plaqs(x_plaq_info, lz_qubits)

    LogicalOpSim(
        z_dist, x_dist, m_dist, theta, syndrome_circuit, noise_model,
        num_qubits, anc_qubits, z_anc_qubits, x_anc_qubits,
        data_qubits, lz_qubits, diagonal_anc,
        z_plaqs, z_plaq_info, z_space_boundary, z_doubled_boundary,
            z_graph_nodes, z_graph, z_costs, z_bpaths, z_path_lengths,
        x_plaqs, x_plaq_info, x_space_boundary, x_doubled_boundary,
            x_graph_nodes, x_graph, x_costs, x_bpaths, x_path_lengths,
    )
end

#####################

function make_plaqs(z_dist, x_dist)
    z_plaqs = Tuple{Int, Int}[
        (x, y)
        for x in 1:z_dist-1
        for y in 0:x_dist
        if mod(x+y, 2) == 0
    ]
    z_space_boundary = Tuple{Int, Int}[
        (x, y)
        for (x, y) in z_plaqs
        if x == 1 || x == z_dist-1
    ]
    z_doubled_boundary = Tuple{Int, Int}[  # Plaqs with two boundary qubits
        (x, y)
        for (x, y) in z_space_boundary
        if 1 <= y <= x_dist-1
    ]
    x_plaqs = Tuple{Int, Int}[
        (x, y)
        for x in 0:z_dist
        for y in 1:x_dist-1
        if mod(x+y, 2) == 1
    ]
    x_space_boundary = Tuple{Int, Int}[
        (x, y)
        for (x, y) in x_plaqs
        if y == 1 || y == x_dist-1
    ]
    x_doubled_boundary = Tuple{Int, Int}[  # Plaqs with two boundary qubits
        (x, y)
        for (x, y) in x_space_boundary
        if 1 <= x <= z_dist-1
    ]
    (z_plaqs, z_space_boundary, z_doubled_boundary), (
        x_plaqs, x_space_boundary, x_doubled_boundary)
end
function make_qubit_assignments(z_dist, z_plaqs, x_dist, x_plaqs)
    counter = Iterators.Stateful(Iterators.countfrom(1))
    data_qubits = Dict{Tuple{Int, Int}, Int}(
        (x, y) => popfirst!(counter)
        for x in 1:z_dist
        for y in 1:x_dist
    )
    lz_qubits = Dict{Tuple{Int, Int}, Int}(
        (z_dist-lz+1, lz) => data_qubits[(z_dist-lz+1, lz)]
        for lz in 1:z_dist
    )
    anc_qubits = Dict{Tuple{Int, Int}, Int}()
    append!(anc_qubits, xy => popfirst!(counter) for xy in z_plaqs)
    z_anc_qubits = [anc_qubits[xy] for xy in z_plaqs]
    append!(anc_qubits, xy => popfirst!(counter) for xy in x_plaqs)
    x_anc_qubits = [anc_qubits[xy] for xy in x_plaqs]
    num_qubits = length(data_qubits) + length(anc_qubits)
    make_plaq_info(plaqs)::Vector{PlaqInfoT} = [
        begin
            qubits = Int[
                data_qubits[(xx, yy)]
                for (xx, yy) in [(x, y), (x, y+1), (x+1, y), (x+1, y+1)]
                if (xx, yy) in keys(data_qubits)
            ]
            @assert(length(qubits) in (2, 4), "Invalid plaquettes")
            (ancilla=anc_qubits[(x, y)], data=qubits)
        end
        for (x, y) in plaqs
    ]
    z_plaq_info = make_plaq_info(z_plaqs)
    x_plaq_info = make_plaq_info(x_plaqs)
    (num_qubits, collect(values(anc_qubits)), collect(values(data_qubits)),
        z_plaq_info, x_plaq_info, collect(values(lz_qubits)), z_anc_qubits, x_anc_qubits)
end
function make_graph(m_dist, plaqs, space_boundary,
                    start_boundary::Bool, end_boundary::Bool)
    graph_nodes = BiMap{NodeT, Int}()
    counter = Iterators.Stateful(Iterators.countfrom(1))
    # Assign node ids
    if start_boundary
        for (x, y) in plaqs
            graph_nodes[(0, :tboundary, x, y)] = popfirst!(counter)
        end
    end
    for t in 1:m_dist
        for (x, y) in plaqs
            graph_nodes[(t, :plaq, x, y)] = popfirst!(counter)
        end
        for (x, y) in space_boundary
            graph_nodes[(t, :sboundary, x, y)] = popfirst!(counter)
        end
    end
    if end_boundary
        for (x, y) in plaqs
            graph_nodes[(m_dist+1, :tboundary, x, y)] = popfirst!(counter)
        end
    end
    # Make graph
    plaq_set = Set(plaqs)
    graph = Graph{Int}(length(graph_nodes))
    # Space edges
    for t in 1:m_dist
        for (x, y) in plaqs
            for (x2, y2) in [(x-1, y-1), (x-1, y+1), (x+1, y-1), (x+1, y+1)]
                (x2, y2) in plaq_set || continue
                add_edge!(graph,
                          graph_nodes[(t, :plaq, x, y)],
                          graph_nodes[(t, :plaq, x2, y2)])
            end
        end
        for (x, y) in space_boundary
            add_edge!(graph,
                      graph_nodes[(t, :plaq, x, y)],
                      graph_nodes[(t, :sboundary, x, y)])
            # Connect boundaries
            if (x, y) != space_boundary[1]
                add_edge!(graph,
                          graph_nodes[(t, :sboundary, x, y)],
                          graph_nodes[(t, :sboundary, space_boundary[1]...)])
            end
        end
    end
    # Time edges
    for t in 1-start_boundary:m_dist-1+end_boundary
        type1 = t == 0 ? :tboundary : :plaq
        type2 = t == m_dist ? :tboundary : :plaq
        for (x, y) in plaqs
            add_edge!(graph,
                      graph_nodes[(t, type1, x, y)],
                      graph_nodes[(t+1, type2, x, y)])
        end
        # Connect boundaries
        if 0 < t < m_dist
            add_edge!(graph,
                      graph_nodes[(t, :sboundary, space_boundary[1]...)],
                      graph_nodes[(t+1, :sboundary, space_boundary[1]...)])
        end
    end
    if start_boundary
        # Connect boundaries
        for (x, y) in plaqs
            add_edge!(graph,
                      graph_nodes[(0, :tboundary, x, y)],
                      graph_nodes[(1, :sboundary, space_boundary[1]...)])
        end
    end
    if end_boundary
        # Connect boundaries
        for (x, y) in plaqs
            add_edge!(graph,
                      graph_nodes[(m_dist+1, :tboundary, x, y)],
                      graph_nodes[(m_dist, :sboundary, space_boundary[1]...)])
        end
    end
    graph_nodes, graph
end
function diagonal_x_plaqs(x_plaq_info, logical_z_qubits)
    # check if neighbor has lz_qubits
    plaqs = x_plaq_info
    lz_qubits = logical_z_qubits
    res = Vector{Int}()
    for p in plaqs
        #println(p[:data])
        if size(p[:data]) == (4,)
            for d in p[:data]
                if d in lz_qubits
                    append!(res, p[:ancilla])
                    break
                end
            end
        end
    end
    res
end

"""
Used by construct_graph_costs().
"""
struct MatchingGraphWeights{SyndromeCircuitT} <: AbstractMatrix{Float64}
    graph_nodes::BiMap{NodeT, Int}
    doubled_boundary::Set{Tuple{Int, Int}}
    noise_model::NoiseModel
    syndrome_circuit::SyndromeCircuitT
    space_edge1::Float64
    time_edge1::Float64
    space_edge::Float64
    time_edge::Float64
end
Base.size(w::MatchingGraphWeights) = (l=length(w.graph_nodes); (l, l))
function Base.getindex(w::MatchingGraphWeights, i, j)::Float64
    r = rev(w.graph_nodes)
    n, m = r[i], r[j]
    # Inter-boundary
    (n[2] != :plaq && m[2] != :plaq) && return 0.0
    # Time edge (including boundary)
    n[1] != m[1] && return (min(n[1], m[1]) == 0 ? w.time_edge1 : w.time_edge)
    # Space edge (including boundary)
    return (n[1] == 1 ? w.space_edge1 : w.space_edge)
    #n[2] != :plaq && ((n, m) = (m, n))  # Make it plaq->boundary or plaq->plaq
    #return (m[2] != :plaq && (m[3], m[4]) in w.doubled_boundary
    #    # Plaq with two boundary qubits
    #    ? w.space_double_edge
    #    # Edge error probability dependent on only one qubit
    #    : w.space_edge
    #)
end
function matching_space_edge(::SyndromeCircuit, noise::NoiseModel,
                             is_x::Bool, m_dist::Int, fisrt_layer::Bool)
    -log(noise.errors[:uniform_data])
end
function matching_time_edge(::SyndromeCircuit, noise::NoiseModel,
                            is_x::Bool, m_dist::Int, first_layer::Bool)
    -log(noise.errors[:uniform_anc])
end

function constuct_graph_costs(graph_nodes::BiMap{NodeT, Int}, graph::Graph{Int},
                              doubled_boundary::Vector{Tuple{Int, Int}},
                              noise_model::NoiseModel,
                              syndrome_circuit::SyndromeCircuit,
                              is_x::Bool, dist::Int, m_dist::Int)
    rev_graph_nodes = rev(graph_nodes)
    weights = MatchingGraphWeights(
        graph_nodes, Set(doubled_boundary), noise_model, syndrome_circuit,
        matching_space_edge(syndrome_circuit, noise_model, is_x, m_dist, true),
        matching_time_edge(syndrome_circuit, noise_model, is_x, m_dist, true),
        matching_space_edge(syndrome_circuit, noise_model, is_x, m_dist, false),
        matching_time_edge(syndrome_circuit, noise_model, is_x, m_dist, false),
    )
    paths = floyd_warshall_shortest_paths(graph, weights)
    boundary_ids = Set{Int}(
        id
        for ((_, kind, _, _), id) in graph_nodes
        if kind != :plaq
    )
    boundary_paths = Set{Tuple{Int, Int}}()
    path_parities = Dict{Tuple{Int, Int}, Int}()
    for i in 1:length(graph_nodes)-1
        for j in i+1:length(graph_nodes)
            #i >= j && continue
            if i in boundary_ids && j in boundary_ids
                push!(boundary_paths, (i, j))
                push!(boundary_paths, (j, i))
                continue
            end
            jj = j
            non_boundary_count = 0
            hits_boundary = false
            while jj != 0
                in_boundary = jj in boundary_ids
                jj2 = paths.parents[i, jj]
                is_time_edge = (
                    !in_boundary && jj2 != 0 && !(jj2 in boundary_ids)
                    && rev_graph_nodes[jj][1] != rev_graph_nodes[jj2][1]
                )
                non_boundary_count += !in_boundary && !is_time_edge
                hits_boundary |= in_boundary
                jj = jj2
            end
            if hits_boundary
                push!(boundary_paths, (i, j))
                push!(boundary_paths, (j, i))
            end
            path_parity = non_boundary_count - !hits_boundary
            if !(i in boundary_ids && j in boundary_ids)
                path_parities[(i, j)] = path_parities[(j, i)] = path_parity
            end
        end
    end
    paths.dists, boundary_paths, path_parities
end

"""
    CodeDistanceRun(ctx)

Contains the per error correction run state.
"""
struct CodeDistanceRun
    ctx::CodeDistanceSim
    state::ChpState
    zx_error_counts::Vector{Int}
    zx_meas_error_counts::Vector{Int}
    z_prev::Vector{Bool}
    x_prev::Vector{Bool}
    z_syndromes::Matrix{Bool}
    x_syndromes::Matrix{Bool}
    sim_noise_params::Dict{Symbol, Float64}
end
function CodeDistanceRun(ctx::CodeDistanceSim)
    state = ChpState(ctx.num_qubits, bitpack=false)
    zx_error_counts = Int[0, 0]
    zx_meas_error_counts = Int[0, 0]
    z_prev = zeros(Bool, length(ctx.z_plaqs))
    x_prev = zeros(Bool, length(ctx.x_plaqs))
    z_syndromes = Matrix{Bool}(undef, length(z_prev), ctx.m_dist+1)
    x_syndromes = Matrix{Bool}(undef, length(x_prev), ctx.m_dist+1)
    sim_noise_params = simulation_noise_parameters(ctx.syndrome_circuit,
                                                   ctx.noise_model, ctx)
    CodeDistanceRun(
        ctx, state, zx_error_counts, zx_meas_error_counts,
        z_prev, x_prev, z_syndromes, x_syndromes,
        sim_noise_params,
    )
end

struct LogicalOpRun
    ctx::LogicalOpSim
    state::ChpState
    zx_error_counts::Vector{Int}
    zx_meas_error_counts::Vector{Int}
    z_prev::Vector{Bool}
    x_prev::Vector{Bool}
    z_syndromes::Matrix{Bool}
    x_syndromes::Matrix{Bool}
    sim_noise_params::Dict{Symbol, Float64}
end
function LogicalOpRun(ctx::LogicalOpSim)
    state = ChpState(ctx.num_qubits, bitpack=false)
    zx_error_counts = Int[0, 0]
    zx_meas_error_counts = Int[0, 0]
    z_prev = zeros(Bool, length(ctx.z_plaqs))
    x_prev = zeros(Bool, length(ctx.x_plaqs))
    z_syndromes = Matrix{Bool}(undef, length(z_prev), ctx.m_dist+1)
    x_syndromes = Matrix{Bool}(undef, length(x_prev), ctx.m_dist+1)
    sim_noise_params = simulation_noise_parameters(ctx.syndrome_circuit,
                                                   ctx.noise_model, ctx)
    #sim_noise_params = nothing
    LogicalOpRun(
        ctx, state, zx_error_counts, zx_meas_error_counts,
        z_prev, x_prev, z_syndromes, x_syndromes,
        sim_noise_params,
    )
end

function simulation_noise_parameters(::SyndromeCircuit, model::NoiseModel,
                                     ctx::CodeDistanceSim)
    Dict{Symbol, Float64}(
        p_data => model.uniform_data,
        p_anc_z => model.uniform_anc,
        p_anc_x => model.uniform_anc,
        p_cnot1 => 0,
        p_cnot => 0,
    )
end
function simulation_noise_parameters(::SyndromeCircuit, model::NoiseModel,
                                     ctx::LogicalOpSim)
    Dict{Symbol, Float64}(
        p_data => model.uniform_data,
        p_anc_z => model.uniform_anc,
        p_anc_x => model.uniform_anc,
        p_cnot1 => 0,
        p_cnot => 0,
    )
end

function apply_sim_error!(noise_model::Nothing,
                          state::ChpState, zx_counts_out::Vector{Int},
                          kind::Symbol, qubits::NTuple{N, Int} where N)
    # Apply no error when no noise model
end
function apply_sim_error!(noise_model::NoiseModel,
                          state::ChpState, zx_counts_out::Vector{Int},
                          kind::Symbol, qubits::NTuple{N, Int} where N)
    p = noise_model.errors[kind]
    # TODO: rng
    z_count = x_count = 0
    for q in qubits
        r = rand(rng, Float64)
        r >= p && return
        # Apply X, Y, or Z with (1/3)p probability
        if r < (2/3)*p
            z_gate!(state, q)
            z_count += 1
        end
        if r >= (1/3)*p
            x_gate!(state, q)
            x_count += 1
        end
    end
    zx_counts_out[2] += z_count  # Z errors cause X syndromes
    zx_counts_out[1] += x_count
    nothing
end

function exec_syndrome_layer(noise_model::Union{NoiseModel, Nothing},
                             syndrome_circuit::BasicSyndrome,
                             run::CodeDistanceRun,
                             layer_i::Int)
    ctx = run.ctx
    state = run.state
    # Apply errors
    for q in ctx.data_qubits
        apply_sim_error!(noise_model, state, run.zx_error_counts,
                         :uniform_data, (q,))
    end
    # Run circuit
    for info in ctx.z_plaq_info
        anc = info.ancilla
        for dat in info.data
            cnot!(state, dat, anc)
        end
    end
    for info in ctx.x_plaq_info
        anc = info.ancilla
        hadamard!(state, anc)
        for dat in info.data
            cnot!(state, anc, dat)
        end
        hadamard!(state, anc)
    end
    # Apply errors
    for q in ctx.anc_qubits
        apply_sim_error!(noise_model, state, run.zx_meas_error_counts,
                         :uniform_anc, (q,))
    end
    # Measure
    for (i, info) in enumerate(ctx.z_plaq_info)
        meas = measure_reset!(state, info.ancilla).value  # TODO: rng
        run.z_syndromes[i, layer_i] = run.z_prev[i] ⊻ meas
        run.z_prev[i] = meas
    end
    for (i, info) in enumerate(ctx.x_plaq_info)
        meas = measure_reset!(state, info.ancilla).value  # TODO: rng
        run.x_syndromes[i, layer_i] = run.x_prev[i] ⊻ meas
        run.x_prev[i] = meas
    end
    nothing
end

function reset_all_qubits!(run::CodeDistanceRun)
    run.state.z[1:run.ctx.num_qubits] .= false  # Clean start
    run.state.x[1:run.ctx.num_qubits] .= false
    run.zx_error_counts .= 0
    run.zx_meas_error_counts .= 0
end

function reset_all_qubits!(run::LogicalOpRun)
    run.state.z[1:run.ctx.num_qubits] .= false  # Clean start
    run.state.x[1:run.ctx.num_qubits] .= false
    run.zx_error_counts .= 0
    run.zx_meas_error_counts .= 0
end

"""
    simulate_syndrome_run(run)

Simulate m_dist rounds of syndrome measurement (with start and end boundaries).
"""
function simulate_syndrome_run(run::CodeDistanceRun)
    simulate_syndrome_run(run.ctx.syndrome_circuit, run)
end
function simulate_syndrome_run(syndrome_circuit::SyndromeCircuit, run::CodeDistanceRun)
    ctx = run.ctx
    reset_all_qubits!(run)  # Clean start
    run.zx_error_counts .= 0
    run.zx_meas_error_counts .= 0
    exec_syndrome_layer(nothing, syndrome_circuit, run, 1)
    for i in 1:ctx.m_dist+1
        noise = i == ctx.m_dist+1 ? nothing : ctx.noise_model  # Noise-free end layer
        exec_syndrome_layer(noise, syndrome_circuit, run, i)
    end
    run.z_syndromes, run.x_syndromes
end

function simulate_logical_op_syndrome_run(run::LogicalOpRun)
    simulate_logical_op_syndrome_run(run.ctx.syndrome_circuit, run)
end
function simulate_logical_op_syndrome_run(syndrome_circuit::SyndromeCircuit, run::LogicalOpRun)
    ctx = run.ctx
    reset_all_qubits!(run)  # Clean start
    state = run.state
    run.zx_error_counts .= 0
    run.zx_meas_error_counts .= 0
    # Logical operation
    for q in ctx.data_qubits
        hadamard!(state, q) # transversal hadamard
    end
    #ctx.logical_op_projection = projected

    exec_syndrome_layer(nothing, syndrome_circuit, run, 1)

    _, _, projected = rz_gate!(state, ctx.theta, ctx.z_dist, ctx.logical_z_qubits)
    for i in 1:ctx.m_dist+1
        noise = i == ctx.m_dist+1 ? nothing : ctx.noise_model  # Noise-free end layer
        exec_syndrome_layer(noise, syndrome_circuit, run, i)
    end
    run.z_syndromes, run.x_syndromes, projected
end

function syndromes_to_error_ids(syndromes, plaqs, graph_nodes)
    @assert size(syndromes)[1] == length(plaqs)
    Int[
        graph_nodes[(t, :plaq, plaqs[i]...)]
        for t in axes(syndromes, 2)
        for i in axes(syndromes, 1)
        if syndromes[i, t]
    ]
end

function syndromes_to_all_errors(syndromes, anc_ids)
    @assert size(syndromes)[1] == length(anc_ids) # same order as plaqs
    Int[
        anc_ids[i]
        for t in axes(syndromes, 2)
        for i in axes(syndromes, 1)
        if syndromes[i, t]
    ]
end

function construct_matching_graph(graph_nodes, costs, extra_boundary::NodeT,
                                  error_ids::Vector{Int})
    node_count = div(length(error_ids) + 1, 2) * 2  # Round up to even
    edge_count = div(node_count * (node_count-1), 2)
    matching = Matching(Float64, node_count, edge_count)
    if node_count > length(error_ids)  # Round up to even
        push!(error_ids, graph_nodes[extra_boundary])
    end
    for i in 1:node_count-1#length(error_locations)-1
        for j in i+1:node_count#length(error_ids)
            c = costs[error_ids[i], error_ids[j]]
            add_edge(matching, i-1, j-1, c)
        end
    end
    matching
end

function count_corrected_errors(matching, error_ids, path_lengths)
    count = 0
    for i in 1:length(error_ids)-1
        j = get_match(matching, i-1) + 1  # get_match is zero-indexed
        i < j || continue  # Count each pair once
        count += path_lengths[error_ids[i], error_ids[j]]
    end
    count
end

function match_and_evaluate_syndromes(plaqs, graph_nodes, space_boundary,
                                      costs, path_lengths,
                                      syndromes, error_count)
    error_ids = syndromes_to_error_ids(syndromes, plaqs, graph_nodes)
    extra_boundary = (1, :sboundary, space_boundary[1]...)  # Doesn't matter which
    matching = construct_matching_graph(graph_nodes, costs, extra_boundary, error_ids)
    solve(matching)
    corrected = count_corrected_errors(matching, error_ids, path_lengths)
    finalize(matching)
    failed = mod(error_count + corrected, 2) != 0
    failed
end

"""
    do_single_run(run, z_only=false)

Simulate and error-correct once.  Returns if there was a Z error, or X error,
and inserted error counts.
"""
function do_single_run(run::CodeDistanceRun, z_only::Bool=false)
    z_syndromes, x_syndromes = simulate_syndrome_run(run)
    z_errors, x_errors = run.zx_error_counts

    ctx = run.ctx
    z_fail = match_and_evaluate_syndromes(
        ctx.z_plaqs, ctx.z_graph_nodes, ctx.z_space_boundary, ctx.z_costs,
        ctx.z_path_lengths, z_syndromes, z_errors)
    if z_only
        x_fail = false
    else
        x_fail = match_and_evaluate_syndromes(
            ctx.x_plaqs, ctx.x_graph_nodes, ctx.x_space_boundary, ctx.x_costs,
            ctx.x_path_lengths, x_syndromes, x_errors)
    end
    z_fail, x_fail, (z_errors, x_errors, run.zx_meas_error_counts...,)
end

"""
    do_n_runs(run, n, z_only=false)

Simulate and error-correct n times.  Returns the error probability.
"""
function do_n_runs(run::CodeDistanceRun, n::Int, z_only::Bool=false)
    fail_count = 0
    if z_only
        for _ in 1:n
            z_fail, = do_single_run(run, true)
            fail_count += z_fail
        end
    else
        for _ in 1:n
            z_fail, x_fail = do_single_run(run)
            fail_count += z_fail || x_fail
        end
    end
    #println(fail_count)
    return fail_count / n
end

#############################################

function do_single_logical_op_run(run::LogicalOpRun, z_only::Bool=false)
    z_syndromes, x_syndromes, projected = simulate_logical_op_syndrome_run(run)
    # if syndromes along the diagonal region, discard: todo
    discarded = false
    logical_err = false
    #z_errors, x_errors = run.zx_error_counts

    ctx = run.ctx
    x_ancs = ctx.x_anc_qubits
    z_ancs = ctx.z_anc_qubits
    z_error_ids = syndromes_to_all_errors(x_syndromes, x_ancs)
    x_error_ids = syndromes_to_all_errors(z_syndromes, z_ancs)
    #for err in error_ids
    #    if err in ctx.diagonal_x_ancilla
    #        discarded = true
    #    end
    #end
    if length(z_error_ids) > 0 || length(x_error_ids) > 0 # discard if any error detected
        discarded = true
    end
    if discarded == false
        theta = ctx.theta
        amp_I = cos(theta/2)
        amp_Z = sin(theta/2)
        dist = ctx.z_dist
        weight = sum(projected)
        if 0 < weight < dist
            logical_err = true
            #println(projected)
        end
        # Actual logical space Z expected value sin(theta_L/2)
        actual_op = (amp_Z^(dist-weight)*amp_I^(weight)) / (sqrt(amp_I^(2*(dist-weight))*amp_Z^(2*weight) + amp_Z^(2*(dist-weight))*amp_I^(2*weight)))
        # Actual angle is theta_L
        actual_angle = 2 * asin(actual_op) # allowed domain for sin is -pi/2 to pi/2
    else
        actual_angle = 0.0
    end

    discarded, logical_err, actual_angle
end

function do_n_logical_op_runs(run::LogicalOpRun, n::Int, z_only::Bool=false)
    fail_count = 0
    logical_err_count = 0
    ctx = run.ctx
    diagonal_anc = ctx.diagonal_x_ancilla
    theta = ctx.theta
    amp_I = cos(theta/2)
    amp_Z = sin(theta/2)
    dist = ctx.z_dist
    # Predicted logical space Z expected value for sin(theta_L/2)
    ideal_op = (amp_Z^dist)/(sqrt(amp_I^(2*dist) + amp_Z^(2*dist)))
    # Predict angle is theta_L
    ideal_angle = 2 * asin(ideal_op) # allowed domain for sin is -pi/2 to pi/2
    actual_sum = 0.0
    weighted_err_sum = 0.0
    for _ in 1:n
        discarded, logical_err, actual_angle = do_single_logical_op_run(run)
        if discarded
            fail_count += 1
        elseif logical_err
            logical_err_count += 1
            weighted_err_sum += 1 * (sin((ideal_angle-actual_angle) / 2))^2
            actual_sum += actual_angle
        else
            actual_sum += actual_angle
        end
    end
    #println(logical_err_count, " ", n, " ", fail_count)
    actual_avg = n == fail_count ? 0.0 : actual_sum / (n-fail_count)
    logical_err_rate = n == fail_count ? 0.0 : logical_err_count / (n-fail_count)
    logical_infidelity = n == fail_count ? 0.0 : weighted_err_sum / (n-fail_count)
    fail_rate = fail_count / n
    #println(actual_avg, " ", actual_sum, " ", fail_count)
    return fail_rate, ideal_angle, actual_avg, logical_infidelity
end

max_stars_repo_name: yongshanding/VLQ

max_issues_repo_name: yongshanding/VLQ

max_forks_repo_name: yongshanding/VLQ


text:
function MakeModel(A,b,C,XLB,XUB,VariableDef,ObjInd;Cuts = [])
    m,n = size(A);
    ObjNum = size(C)[1];
    MyModel = Model();
    @variable(MyModel,XLB[i] <= X[i=1:n] <= XUB[i],category = VariableType(VariableDef,i))#
    @objective(MyModel, Min, sum(C[ObjInd,i]*X[i] for i=1:n))
    for i = 1:m
        @constraint(MyModel,sum(A[i,j]*X[j] for j = 1:n) >= b[i] )
    end
    return MyModel;
end

function MakeModel(A,b,C,XLB,XUB,ObjInd;Cuts = [])
    m,n = size(A);
    ObjNum = size(C)[1];
    MyModel = Model();
    @variable(MyModel, XLB[i] <= X[i=1:n] <= XUB[i])
    @objective(MyModel, Min, sum(C[ObjInd,i]*X[i] for i=1:n))
    for i = 1:m
        @constraint(MyModel,sum(A[i,j]*X[j] for j = 1:n) >= b[i] )
    end
    return MyModel;
end


text:
# Write your tests here.
graph = [
    [2, 3],
    [1], # list of nodes that can be reached from node 2
    [1], # ", node 3
    [5], # ", node 4
    [4]] # ", node 5
dir_graph = [
[2],
[3,5], # list of nodes that can be reached from node 2
[4], # ", node 3
[1,2], # ", node 4
[4]] # ", node 5
# @test TDD.get_neighbors(graph, 1) == [1,2,3] # note the order.
# @test TDD.get_neighbors(graph, 2) == [1,2]
# @test TDD.get_neighbors(graph, 5) == [4,5]
# @test TDD.get_neighbors(graph, 5) == [4,5] # test if graph will change after one evaluation of function.
# @test TDD.get_neighbors(dir_graph, 2) == [2,3,5]
# @test TDD.get_neighbors(dir_graph, 5) == [4,5]
# Test for the 1st problem, neighbors
@test all(i in TDD.get_neighbors(graph, 1) for i in graph[1])
@test all(i in TDD.get_neighbors(graph, 5) for i in graph[5])
@test all(i in TDD.get_neighbors(dir_graph, 3) for i in dir_graph[3])
@test all(i in TDD.get_neighbors(dir_graph, 4) for i in dir_graph[4])
for i in 1:length(dir_graph)
@test i in TDD.get_neighbors(dir_graph, i) # include the starting index as neighbor
end
# Test for the 2nd problem
@test all(i in TDD.get_connected(graph, 1) for i in [1,2,3])
@test all(i in TDD.get_connected(graph, 4) for i in [4,5])
@test all(i in TDD.get_connected(dir_graph, 1) for i in 1:5)
# Test for the 3rd problem
@test all(comp in TDD.connected_components(graph) for comp in [Set([1,2,3]),Set([4,5])])
@test Set([1:5...]) in TDD.connected_components(dir_graph)

A = Bool[
1 1 1 0 0;
1 1 0 0 0;
1 0 1 0 0;
0 0 0 1 1;
0 0 0 1 1;]

Adir = Bool[
0 1 0 0 0;
0 0 1 0 1;
0 0 0 1 0;
1 1 0 1 0;
0 0 0 1 0;]
@test all(i in TDD.get_neighbors(A, 1) for i in graph[1])
@test all(i in TDD.get_neighbors(A, 5) for i in graph[5])
@test all(i in TDD.get_neighbors(Adir, 3) for i in dir_graph[3])
@test all(i in TDD.get_neighbors(Adir, 4) for i in dir_graph[4])
for i in 1:TDD.n_nodes(Adir)
@test i in TDD.get_neighbors(Adir, i) # include the starting index as neighbor
end


text:
# -*- coding: utf-8 -*-
# lib_fnc_defs

VERSION >= v"0.4.0-dev+6521" && __precompile__()
module DLL_Loader

import Relocator

export load, unload

# must load :freetype before :d3dxfreetype2 (or place to current directory)
# must load :dx9adl before :d3dxglyph (or place to current directory)
const _dlls = [:d3d9, :d3dx9, :d3dxconsole, :freetype, :d3dxfreetype2,
  :d3dxtexturebmp, :dx9adl, :d3dxglyph]

function load(bp::AbstractString="")
  Relocator._init(_dlls, bp)
end

function unload()
  Relocator._close()
end

end # DLL_Loader

macro juliaobj(ptr)
  quote
    unsafe_pointer_to_objref($ptr)
  end
end

macro ptr(obj)
  quote
    pointer_from_objref($obj)
  end
end

macro ptr_as(typ, obj)
  quote
    convert(Ptr{$typ}, pointer_from_objref($obj))
  end
end

# typ must be bitstype (TypStructBits <- TypStruct)
macro as_bits(typ, obj)
  quote
    pointer_to_array((@ptr_as $typ $obj), 1)[]
  end
end

# typ must be immutable (ImmuStruct <- TypStructBits)
macro bits_to(typ, obj)
  quote
    pointer_to_array((@ptr_as $typ $obj), 1)[]
  end
end

# typ must be immutable (ImmuStruct <- Array{TypStructBits,1}[idx])
macro bits_elem(typ, ary, idx)
  quote
    pointer_to_array(convert(Ptr{$typ}, pointer($ary)), idx)[idx]
  end
end

macro as_array(m)
  quote
    pointer_to_array((@ptr_as Float32 $m), (4, 4))
  end
end

macro q2m(mq, m0, m1) # transposed matrix reversed multiply
  quote
    array_to($mq, (@as_array $m1) * (@as_array $m0))
  end
end

macro q3m(mq, m0, m1, m2) # transposed matrix reversed multiply
  quote
    array_to($mq, (@as_array $m2) * (@as_array $m1) * (@as_array $m0))
  end
end

import Relocator: _mf, @mf, @cf, @wf

@wf kernel32 UInt32 GetModuleHandleA (Ptr{Void},)

@cf Ptr{Void} memcpy (Ptr{Void}, Ptr{Void}, UInt32,)

@wf d3d9 Ptr{Void} Direct3DCreate9 (UInt32,)

@wf d3dx9 Ptr{Void} D3DXVec3Project (Ptr{Void}, Ptr{Void}, Ptr{Void},
                      Ptr{Void}, Ptr{Void}, Ptr{Void},)
@wf d3dx9 UInt32 D3DXCreateFontIndirectA (Ptr{Void},
                   Ptr{Void}, Ptr{Ptr{Void}},)
@wf d3dx9 UInt32 D3DXCreateFontA (Ptr{Void},
                   UInt32, UInt32, UInt32, UInt32, UInt32,
                   UInt32, UInt32, UInt32, UInt32,
                   Ptr{Void}, Ptr{Ptr{Void}},)
@wf d3dx9 UInt32 D3DXCreateSprite (Ptr{Void}, Ptr{Ptr{Void}},)
@wf d3dx9 UInt32 D3DXCreateTextureFromFileA (Ptr{Void},
                   Ptr{Void}, Ptr{Ptr{Void}},)
@wf d3dx9 Ptr{Void} D3DXMatrixRotationX (Ptr{Void}, Float32,)
@wf d3dx9 Ptr{Void} D3DXMatrixRotationY (Ptr{Void}, Float32,)
@wf d3dx9 Ptr{Void} D3DXMatrixRotationZ (Ptr{Void}, Float32,)
@wf d3dx9 Ptr{Void} D3DXMatrixLookAtLH (Ptr{Void},
                      Ptr{Void}, Ptr{Void}, Ptr{Void},)
@wf d3dx9 Ptr{Void} D3DXMatrixPerspectiveFovLH (Ptr{Void},
                      Float32, Float32, Float32, Float32,)

@mf d3dxconsole Void debugalloc ()
@mf d3dxconsole Void debugout (Ptr{UInt8},)
@mf d3dxconsole Void debugout (Ptr{UInt8}, UInt32,)
@mf d3dxconsole Void debugout (Ptr{UInt8}, UInt32, UInt32,)
@mf d3dxconsole Void debugfree ()

@mf d3dxfreetype2 UInt32 D3DXFT2_GlyphAlloc (Ptr{Void}, # GLYPH_TBL
                           UInt32, UInt32,)
@mf d3dxfreetype2 UInt32 D3DXFT2_GlyphFree (Ptr{Void},) # GLYPH_TBL
@mf d3dxfreetype2 UInt32 D3DXFT2_GlyphBmp (Ptr{Void}, # GLYPH_TBL
                           Ptr{Void}, Int32, Int32,) # FT_Bitmap
@mf d3dxfreetype2 UInt32 D3DXFT2_GlyphShowBmp (Ptr{Void},) # GLYPH_TBL
@mf d3dxfreetype2 UInt32 D3DXFT2_GlyphOutline (Ptr{Void},) # GLYPH_TBL
@mf d3dxfreetype2 UInt32 D3DXFT2_Status (Ptr{Void},) # GLYPH_TBL

@mf d3dxtexturebmp UInt32 D3DXTXB_CreateTexture (Ptr{Void}, # LPDIRECT3DDEVICE9
                            UInt32, UInt32,
                            Ptr{Ptr{Void}},) # LPDIRECT3DTEXTURE9 *
@mf d3dxtexturebmp UInt32 D3DXTXB_LockRect ( # LPDIRECT3DTEXTURE9 *
                            Ptr{Ptr{Void}}, Ptr{Void},) # D3DLOCKED_RECT *
@mf d3dxtexturebmp UInt32 D3DXTXB_UnlockRect ( # LPDIRECT3DTEXTURE9 *
                            Ptr{Ptr{Void}},)
@mf d3dxtexturebmp UInt32 D3DXTXB_RenderFontGlyph (Ptr{Void}, # GLYPH_TBL
                            Ptr{Void}, Int32, Int32,) # FT_Bitmap
@mf d3dxtexturebmp UInt32 D3DXTXB_CopyTexture ( # D3DLOCKED_RECT *
                            Ptr{Void}, Ptr{Void}, Int32,)
@mf d3dxtexturebmp UInt32 D3DXTXB_RewriteTexture ( # LPDIRECT3DTEXTURE9 *
                            Ptr{Ptr{Void}}, Ptr{Ptr{Void}},)

@mf dx9adl UInt32 ReleaseNil (Ptr{Ptr{Void}},)
@mf dx9adl Ptr{Ptr{Void}} PtrPtrS (Ptr{Void}, UInt32,) # RenderD3DItemsState
@mf dx9adl Ptr{Void} PtrSO (Ptr{Void}, UInt32,) # RenderD3DItemsState
@mf dx9adl Ptr{Ptr{Void}} PtrPtrU (Ptr{Void}, UInt32,) # RenderD3DItemsState
@mf dx9adl Ptr{Void} PtrUO (Ptr{Void}, UInt32,) # RenderD3DItemsState
@mf dx9adl UInt32 SetupMatrices (Ptr{Void},) # RenderD3DItemsState
@mf dx9adl UInt32 DrawAxis (Ptr{Void},) # RenderD3DItemsState
@mf dx9adl UInt32 DX9_CreateVertexBuffer (Ptr{Void}, # LPDIRECT3DDEVICE9
                    UInt32, Ptr{Ptr{Void}},) # LPDIRECT3DVERTEXBUFFER9 *
@mf dx9adl UInt32 DX9_Lock (Ptr{Ptr{Void}}, # LPDIRECT3DVERTEXBUFFER9 *
                    UInt32, Ptr{Ptr{Void}},) # VOID **
@mf dx9adl UInt32 DX9_Unlock (Ptr{Ptr{Void}},) # LPDIRECT3DVERTEXBUFFER9 *
@mf dx9adl UInt32 DX9_StoreVertices (
                    Ptr{Ptr{Void}}, # LPDIRECT3DVERTEXBUFFER9 *
                    UInt32, Ptr{Void},) # CUSTOMVERTEX *
@mf dx9adl UInt32 DX9_LoadVertices (
                    Ptr{Void}, # CUSTOMVERTEX *
                    UInt32, Ptr{Ptr{Void}},) # LPDIRECT3DVERTEXBUFFER9 *
@mf dx9adl UInt32 DX9_SetTexture (Ptr{Void}, # LPDIRECT3DDEVICE9
                    Ptr{Ptr{Void}}, # LPDIRECT3DTEXTURE9 *
                    UInt32,)
@mf dx9adl UInt32 DX9_RenderVertices (Ptr{Void}, # LPDIRECT3DDEVICE9
                    Ptr{Ptr{Void}}, # LPDIRECT3DVERTEXBUFFER9 *
                    UInt32, UInt32, UInt32,)
@mf dx9adl UInt32 DrawString (Ptr{Void}, # RenderD3DItemsState
                    UInt32, Ptr{Cchar}, UInt32,
                    Float32, Float32, Float32, Float32, Float32, Float32,)
@mf dx9adl UInt32 BltString (Ptr{Void}, # RenderD3DItemsState
                    UInt32, Ptr{Cchar}, UInt32, UInt32, UInt32, Float32,)
@mf dx9adl UInt32 BltTexture (Ptr{Void}, # RenderD3DItemsState
                    UInt32, Ptr{Ptr{Void}}, UInt32, UInt32, UInt32, UInt32,
                    Float64, Float64, Float64, Float64, Float64, Float64,)
@mf dx9adl UInt32 InitD3DApp (UInt32, UInt32, Ptr{UInt8}, Ptr{UInt8},
                    Ptr{Void}, # RenderD3DItemsState
                    Ptr{Void}, Ptr{Void}, Ptr{Void},) # cfunc, cfunc, cfunc
@mf dx9adl UInt32 MsgLoop (Ptr{Void},) # RenderD3DItemsState
@mf dx9adl UInt32 DX9_Mode (Ptr{Void},) # RenderD3DItemsState
@mf dx9adl UInt32 DX9_Status (Ptr{Void},) # RenderD3DItemsState

@mf d3dxglyph UInt32 D3DXGLP_InitFont ()
@mf d3dxglyph UInt32 D3DXGLP_CleanupFont ()
@mf d3dxglyph UInt32 D3DXGLP_GlyphContours (Ptr{Void},) # GLYPH_TBL
@mf d3dxglyph UInt32 D3DXGLP_DrawGlyph (Ptr{Void},) # GLYPH_TBL


text:
@testset "Testing markov_approx.jl" begin
    rough_kwargs = Dict(:atol => 1e-8, :rtol => 1e-8)

    # set up

    ρ, σ_u = rand(2)
    μ = 0.0

    for n=3:25
        mc = rouwenhorst(n, ρ, σ_u, μ)

        @test isapprox(sum(mc.state_values) , 0.0; rough_kwargs...)

        for m=1:4
            mc = tauchen(n, ρ, σ_u, μ, m)
            @test isapprox(sum(mc.state_values), 0.0; rough_kwargs...)
        end
    end
end  # @testset


text:
import LinearAlgebra: lowrankupdate!, lowrankdowndate!

# Matlab style identity matrix convenience function
eye(FT, d) = diagm(0 => fill(one(FT),d))
eye(d) = eye(Float64, d)

# Helper to calculate z = sc_y*y  + sc*A*x (in place), where A assumed
# lower triangular. Argument y can be same as z.
# NB: No checks about dimensions!
@inline function lowerTriInplaceMultiplyAdd!(z::zT,
    sc::FT, A::AT, x::xT, sc_y::FT, y::yT) ::Nothing where {
    FT <: AbstractFloat,
    zT <: AbstractVector{FT},
    AT <: AbstractMatrix{FT},
    xT <: AbstractVector{FT},
    yT <: AbstractVector{FT}
    }
    n::Int64 = length(x)
    if (sc_y == zero(FT))
        @inbounds @simd for i = 1:n
            u::FT = zero(FT)
            for j = 1:i
                u += A[i,j]*x[j]
            end
            z[i] = sc*u
        end
    else
        @inbounds @simd for i = 1:n
            u::FT = zero(FT)
            for j = 1:i
                u += A[i,j]*x[j]
            end
            z[i] = sc_y*y[i] + sc*u
        end
    end
    nothing
end

# These are from cholesky.jl: just to make these work with MVector...
function lowrankupdate!(C::Cholesky{FT,<:MMatrix{n,n}},
    v::MVector{n,FT}) ::Nothing where {n, FT<:AbstractFloat}

    A = C.factors
    if C.uplo == 'U'
        conj!(v)
    end

    @inbounds for i = 1:n

        # Compute Givens rotation
        c, s, r = LinearAlgebra.givensAlgorithm(A[i,i], v[i])

        # Store new diagonal element
        A[i,i] = r

        # Update remaining elements in row/column
        if C.uplo == 'U'
            @simd for j = i + 1:n
                Aij = A[i,j]
                vj  = v[j]
                A[i,j]  =   c*Aij + s*vj
                v[j]    = -s'*Aij + c*vj
            end
        else
            @simd for j = i + 1:n
                Aji = A[j,i]
                vj  = v[j]
                A[j,i]  =   c*Aji + s*vj
                v[j]    = -s'*Aji + c*vj
            end
        end
    end
    nothing
end

function lowrankdowndate!(C::Cholesky{FT,<:MMatrix{n,n}},
    v::MVector{n,FT}) ::Nothing where {n, FT<:AbstractFloat}
    A = C.factors
    if C.uplo == 'U'
        conj!(v)
    end

    @inbounds for i = 1:n

        Aii = A[i,i]

        # Compute Givens rotation
        s = conj(v[i]/Aii)
        s2 = abs2(s)
        if s2 > one(FT)
            throw(LinearAlgebra.PosDefException(i))
        end
        c = sqrt(one(FT) - abs2(s))

        # Store new diagonal element
        A[i,i] = c*Aii

        # Update remaining elements in row/column
        if C.uplo == 'U'
            @simd for j = i + 1:n
                vj = v[j]
                Aij = (A[i,j] - s*vj)/c
                A[i,j] = Aij
                v[j] = -s'*Aij + c*vj
            end
        else
            @simd for j = i + 1:n
                vj = v[j]
                Aji = (A[j,i] - s*vj)/c
                A[j,i] = Aji
                v[j] = -s'*Aji + c*vj
            end
        end
    end
    nothing
end


text:
# test 1: compare SDDP and extensive formulation
using Distributed;
addprocs(30);
@everywhere include("loadMod.jl");
@everywhere const GUROBI_ENV = Gurobi.Env();
#pmap(i -> importIpopt(),1:30);

τ = 2;
T = 12;
Δt = 0.25;
N = 5;
###########################################################################
data = [];
ωSet = [5,(2,9),(8,12)];

for T in [8,12,16]
    for j in procs()
        remotecall_fetch(readInData_old,j,T,ωSet,10000);
    end

    # global mExt = Model(solver = IpoptSolver(linear_solver = "ma27"));
    global mExt = Model(solver = GurobiSolver());
    startTE = time();
    global mExt = extForm(1, 0, [[],bData.bInv,[]], 1, τ, Δt, T, fData, bData, dData, pDistr, true);
    solve(mExt);
    elapsedTE = time() - startTE;
    mExtObj = getobjectivevalue(mExt);

    startT = time();
    #cutDict,LBHist,UBHist,UBuHist,UBlHist = solveMain(τ, T, Δt, fData, pDistr, bData, dData, N, false, 20, 20);
    cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(τ, T, Δt, N, false, false, 100, 100);
    elapsedT = time() - startT;
    push!(data,(T,elapsedTE,mExtObj,elapsedT,LBHist));
end
save("extensiveSDDP.jld","data",data);

τ = 2;
T = 8;
Δt = 0.25;
N = 5;
disAdd = "testProbRead_96.csv";
pDistr = readDisruption(disAdd,"csv");
ωSet0 = [i for i in keys(pDistr.ωDistrn)];
ωSetTot = [[5,(2,9),(8,12)],
        [2,3,5,(2,3),(2,9),(8,12),(9,12),(10,13)],
        ωSet0];

for i in 1:3
    for j in procs()
        remotecall_fetch(readInData_old,j,T,ωSetTot[i]);
    end

    global mExt = Model(solver = GurobiSolver(NumericFocus = 3));
    startTE = time();
    global mExt = extForm(1, 0, [[],bData.bInv,[]], 1, τ, Δt, T, fData, bData, dData, pDistr, true);
    solve(mExt);
    elapsedTE = time() - startTE;
    mExtObj = getobjectivevalue(mExt);

    startT = time();
    cutDict,LBHist,UBHist,UBuHist,UBlHist,timeHist = solveMain(τ, T, Δt, N, false, false, 100, 100);
    elapsedT = time() - startT;
    push!(data,(T,elapsedT,cutDict,LBHist,UBHist,UBuHist,UBlHist));
end


text:
"""
cohorts as 3 matrices: year_b - the initial year,j0 - youngest age, j1 - oldest age, jmax - maximal age,  x - abundance, y - catch, w - weight (obsolete)
"""
type cohort1
    year_b
    j0
    j1
    jmax
     x::Array{Float64,2}
     y::Array{Float64,2}
     w::Array{Float64,2}
 end

"""
cohorts as 3 matrices: year_b - the initial year,   x - abundance, y - catch, w - weight (standard)
"""

type cohort
     year_b

     x::Array{Float64,2}
     y::Array{Float64,2}
     w::Array{Float64,2}
 end

"""
list of cohorts (obsolete)
"""
type  cohorts
     year_b
     tmax
     jmax
     cs::Array{cohort,1}
end

"""
cohorts as 3 matrices: year_b - the initial year, tmax - number of seasons, jmax - maximal age,  x - abundance, y - catch, w - weight (obsolete)
"""

type  cohorts1
     year_b
     tmax
     jmax
     cs::Array{cohort1,1}
end

"""
cohorts as 3 matrices: year_b - the initial year, x - abundance, y - catch, w - weight (standard)
"""
type  cohorts_m
     year_b
    x
     y
    w

end

"""
cohorts as 3 matrices: year_b - initial year, year_e - last year,  x - abundance, y - catch, w - weight (obsolete)
"""
type  cohorts_t
     year_b
    year_e
    x
     y
    w

end

"""
cohort as model parameters: year_b - initial year, j_max - maximum age,  x1j - first row of x,  w - weight, k_x, k_w - coefficients for x and y, p_m - natural mortality parameters, h - H matrix,  p_s - - fishing mortality parameters,  ft
"""
type cohort_mp
     year_b
     j_max
     x1j::Array{Float64,2}# row
     w::Array{Float64,2}
     k_x::Float64
     k_w::Float64
     p_m
     h::Array{Float64,2}
     p_s
    ft::Array{Float64,1}# column
end

"""
transforms cohort from cohort_mp to  cohort_m
"""
function cmp2cm(cmp)

    jmax=cmp.j_max
     x1j=cmp.x1j
          kx=cmp.k_x
     kw=cmp.k_w
     a1m1,a2m1,b1m1,b2m1,cm1,d1m1,d2m1=cmp.p_m
     Hb=cmp.h
     as1,bs1=cmp.p_s
    ftb=cmp.ft

    tmax=length(ftb)

    x,y=gm_burago(tmax,jmax,Mb,Fb,W,kx,kw,x1j)#

    return cohorts_m(cmp.year_b,x,y,cmp.w)

end

# from cohort_m to cohort_mp !

max_stars_repo_name: picestcode/Fishmetica

max_issues_repo_name: picestcode/Fishmetica

max_forks_repo_name: picestcode/Fishmetica


text:
#
# Download index and extra data specific from suite sparse at TAMU
#

using MAT
using DataFrames
using Base.Filesystem

#const SS_SITE = "https://sparse.tamu.edu"
#const SS_FILES = "/files/ssstats.csv"
#const SS_INDEX = "/files/ss_index.mat"
#const SS_GRAPHICS = "/files/" # to be followed by /$Group/$(Name)$ending
# where endings indicate graphics about the main matrix of different kinds
# ".png"|"_thumb.png"               - sparsity structure
# "_graph.gif"|"_graph_thumb.gif"   - iconnectivity graph
# "_svd.png"                        - singular values plot

const SS_SVDDIR = "/svd/" # to be followed by "$Group/$(Name)_SVD{.mat|_piroband.mat}

# the files "$SS_SITE/$Group/$Name" contain html-formatted further information
# the directories "/mat" "/MM" "/RB"
# contain data files in $Group/$Name.mat or $Group/Name.tar.gz" in different formats
#
# the SVD-based statistics data are only available on the "further information html file"
# They may be derived from the complete list of SVD values from "/SVD" if available
#
# file "/statistics" contains html with legend of data contained in "ss_index"
# see also a copy at the end of this file
#

const MFIELDNAMES = [:Group, :Name, :nrows, :ncols, :nnz, :nzero, :nnzdiag,
                    :pattern_symmetry, :numerical_symmetry,
                    :isBinary, :isReal, :isND, :isGraph, :posdef, :cholcand,
                    :amd_lnz, :amd_vnz, :amd_rnz,:amd_flops,
                    :nblocks, :ncc, :sprank, :RBtype,
                    :lowerbandwidth, :upperbandwidth,
                    :rcm_lowerbandwidth, :rcm_upperbandwidth,
                    :xmin, :xmax]
const MFIELDTYPES = [String, String, Int, Int, Int, Int, Int,
                     Float64, Float64,
                     Bool, Bool, Bool, Bool, Bool, Bool,
                     Int, Int, Int, Int,
                     Int, Int, Int, String,
                     Int, Int, Int, Int,
                     ComplexF64, ComplexF64]
const MFIELDNAMES2 =[:nnzdiag,
                    :pattern_symmetry, :numerical_symmetry,
                    :isND, :isGraph, :posdef, :cholcand,
                    :amd_lnz, :amd_vnz, :amd_rnz,:amd_flops,
                    :nblocks, :ncc, :sprank,
                    :lowerbandwidth, :upperbandwidth,
                    :rcm_lowerbandwidth, :rcm_upperbandwidth,
                    :xmin, :xmax]
const MFIELDTYPES2 = [Int,
                     Float64, Float64,
                     Bool, Bool, Bool, Bool,
                     Int, Int, Int, Int,
                     Int, Int, Int,
                     Int, Int, Int, Int,
                     ComplexF64, ComplexF64]

function readindex(remote::SSRemoteType, db::MatrixDatabase)
    df = read_ss_index()
    for id = 1:size(df, 1)
        mt = copy(df[id,:])
        name = mt.name
        data = get!(db.data, name) do
            info = MetaInfo(mt.nrows, mt.ncols, mt.nnz, 0, "", 0, "", "", "", "", "")
            RemoteMatrixData{typeof(remote)}(name, id, info)
        end
        for (name, T) in zip(MFIELDNAMES2, MFIELDTYPES2)
            setproperty!(data.header, name, (getfield(mt, name)))
        end
    end
end

function MetaInfo(a::Int64, b::Int64, c::Int64, d::Int64, e::String,
                  f::Int64, g::String, h::String, i::String, j::String, k::String)

    MetaInfo(a, b, c, d, e, f, g, h ,i, j, k,
            0, 0.0, 0.0, false, false, false, false, 0, 0, 0, 0, 0, 0, 0,
            0, 0, 0, 0, Complex(0.0), Complex(0.0),
            false, 0.0, 0.0, 0.0, 0, 0, 0.0, "", "", Float64[])
end

function load_ss_index()
    file = localindex(preferred(SSRemoteType))
    if !isfile(file) || Base.Filesystem.stat(file).size == 0
        url = redirect(indexurl(SS_REMOTE))
        println("downloading: ", url)
        downloadfile(url, file)
    end
    file
end

"""
    read_ss_index(d::DataFrame)

Read file "ss_index.mat" from local data directory and add all fields to `d`.
Read file "ssstats.csv" and add field "kinds" to `d`.
As a result, all available metadata for all problems in `"SuiteSparseMatrixCollection"`
as in the data frame `d`.
"""
function read_ss_index(d::DataFrame=DataFrame())
    file = load_ss_index()
    m = matopen(file) do io
        read(io, "ss_index")
    end
    for (fn, T) in zip(MFIELDNAMES, MFIELDTYPES)
        a = m[string(fn)]
        b = vec([from_matdat(T, s) for s in a])
        d[!,fn] = b
    end
    d[!,:nnz] += d[!,:nzero]
    jname(a::String, b::String) = a * '/' * b
    d[!,:name] = jname.(d[!,:Group],d[!,:Name])
    select!(d, DataFrames.Not(:Group))
    d
end

from_matdat(T::Type, s::Any) = T(s)
from_matdat(::Type{String}, s::AbstractString) = s

"""
    loadsvd(data::RemoteMatrixData)

Download the extra information related to SVD from the backing repository.
That is currently the TAMU site for the UFl collection.
The files are in MAT format, and are uncompressed and un-tar-ed if necessary.
"""
function loadsvd(data::RemoteMatrixData{SSRemoteType})
    file = svdfile(data)
    if isfile(file)
        return 0
    end
    dir = dirname(localdir(data))
    url = redirect(datasvdurl(data))

    isdir(dir) || mkpath(dir)

    try
        println("downloading: ", url)
        downloadfile(url, file)
        addsvd!(data)
        1
    catch
        rm(file, force=true)
        write(file, "")
        0
    end
end
loadsvd(data::MatrixData) = 0

function addsvd!(data::RemoteMatrixData{SSRemoteType})
    file = svdfile(data)
    if !isempty(file) && Filesystem.stat(file).size > 0
        matopen(file) do io
            d = read(io, "S")
            pushsvd!(data, d["status"], d["how"], d["s"])
        end
    else
        0
    end
end
addsvd!(::MatrixData) = 0

function datasvdurl(data::RemoteMatrixData)
    string(siteurl(data), '/', "svd", '/', data.name, "_SVD.mat")
end
function svdfile(data::RemoteMatrixData)
    abspath(localdir(data), string(basename(data.name), "_SVD.mat"))
end

function notesfile(data::RemoteMatrixData)
    abspath(localdir(data), string(basename(data.name), "_notes.txt"))
end

function pushsvd!(data::RemoteMatrixData, status::AbstractString, how::AbstractString, sv)
    info = data.header
    status = lowercase(status)
    sv = _vector(sv)
    res = info.svdstatus != status || info.svdhow != how || !isequal(info.sv, sv)
    info.svdok = status == "ok"
    info.svdstatus = status
    info.svdhow = how
    info.sv = sv
    st = svdstatistics(info.sv)
    info.norm = st.maxsv
    info.minsv = st.minsv
    info.cond = st.cond
    info.rank = st.rank
    info.nullspace = st.nsd
    info.svgap = st.gap
    res
end
_vector(sv::AbstractArray) = vec(sv)
_vector(sv::Number) = [sv]

datadetailsurl(data::RemoteMatrixData) = string(siteurl(data), '/', data.name)

"""
    svdstatistics(svd::Vector{<:Real})

Return a names tuple containing several values derivable from the svd vector.
The vector should be sorted in decreasing order.

Return type: `NamedTuple{(:norm, :svdmin, :cond, :rank, :nsd, :gap)}
* maxsv: maximal singular value
* minsv: minimal singular value
* cond: condition number = `maxsv / minsv`
* rank: number of singular values > `tol = (length(sv)) * eps(norm)`
* nsd: nullspace dimension ? `length(sv) - rank`
* gap: relation between smallest `sv > tol` and greatest `sv <= tol`

"""
function svdstatistics(sv::Vector{T}) where T<:Real
    mn = length(sv)
    minsv, maxsv = mn > 0 ? extrema(sv) : (zero(T), zero(T))
    norm = maxsv
    cond = (iszero(minsv) && mn > 0 ? one(T) : maxsv) / minsv
    tol = mn * eps(norm)
    rank = count(x -> x > tol, sv)
    nsd = mn - rank
    mint = maxsv
    maxt = zero(T)
    for s in sv
        if s > tol && s < mint
            mint = s
        elseif s < tol && s > maxt
            maxt = s
        end
    end
    gap = (iszero(maxt) && mn > 0 ? one(T) : mint) / maxt
    (maxsv = norm, minsv = minsv, cond = cond, rank = rank, nsd = nsd, gap = gap)
end

#=

Statistics computed for the SuiteSparse Matrix Collection

LastRevisionDate
    This is a single string kept in the UF_index MATLAB struct that states when the collection or the index was last modified.
DownloadTimeStamp
    The date and time the index that you last downloaded the index.
Group
    A cell array. Group{id} is the group name for the matrix whose serial number is 'id'. Each matrix has a unique id number in the range of 1 to the number of matrices in the collection. Once an id is assigned to a matrix, it never changes.
Name
    Name{id} is the name of the matrix (excluding the Group). Name{id} is not unique. The full name of a matrix should always be given as Group/Name.
nrows
    The number of rows in the matrix.
ncols
    The number of columns in the matrix.
nnz
    The number of numerically nonzero entries in the matrix, or nnz(Problem.A) in MATLAB, where Problem=UFget(id) is a struct containing the MATLAB format of the problem. This statistic can differ from the number of 'entries' explicitly stored in the matrix, however, since some of these entries may be numerically zero. In the MATLAB format, these explicit zero entries are stored in the binary Problem.Zeros matrix, since MATLAB drops all explicit zeros from its sparse matrix storage. The Problem.A matrix in MATLAB has nnz entries in it, with no explicit zeros. In the Matrix Market and Rutherford-Boeing format, a single file holds all entries, both nonzero and the explicit zero entries.
nzero
    The number of explicit entries present in the matrix that are provided by the matrix author but which are numerically zero. nzero(id) is nnz(Problem.Zeros).
pattern_symmetry
    Let S=spones(Problem.A) be the binary pattern of the explicit nonzeros in the matrix. Let pmatched be the number of matched offdiagonal entries, where both S(i,j) and S(j,i) are one, with i not equal to j. Let nzoffdiag be the number of offdiagonal entries in S. Then pattern_symmetry is the ratio of pmatched/nzoffdiag. Note that if S(i,j) and S(j,i) are both one, then this pair of entries is counted twice in both pmatched and nzoffdiag. If the matrix is rectangular, this statistic is zero. If there are no offdiagonal entries, the statistic is 1.
numerical_symmetry
    Let xmatched be the number of matched offdiagonal entries, where A(i,j) is equal to the complex conjugate of A(j,i) and where i and j are not equal. Then numerical_symmetry is the ration xmatched / nzoffdiag (or 1 if nzoffdiag is zero). This statistic is zero for rectangular matrices. Note that this statistic measures how close a matrix is to being Hermitian (A=A' in MATLAB). For complex symmetric matrices (A=A.' in MATLAB), this ratio will be less than one (unless all offdiagonal entries are real).
isBinary
    1 if the matrix is binary (all entries are 0 or 1), 0 otherwise.
isReal
    1 if the matrix is real, 0 if complex.
nnzdiag
    The number of numerically nonzero entries on the diagonal (nnz (diag (Problem.A)) in MATLAB notation). This statistic ignores explicit zero entries (Problem.Zeros in the MATLAB struct).
posdef
    1 if the matrix is known to be symmetric positive definite (or Hermitian positive definite for the complex case), 0 if it is known not to be, -1 if it is symmetric (or Hermitian) but with unknown positive-definiteness. If the statistic is unknown (-1) this may be revised in subsequent versions of the index.
amd_lnz
    Let C=S+S' where S=spones(A) is the binary pattern of A. Then amd_lnz is number of nonzeros in the Cholesky factorization of the matrix C(p,p) (assuming C is positive definite) where p=amd(C) is the AMD fill-reducing ordering. This statistic is -2 for rectangular matrices or for graph problems. It is -1 if it is not computed. This figure gives an estimate of the memory requirements for x=A\b in MATLAB for square matrices.
amd_flops
    The floating-point operation count for computing the Cholesky factorization of C(p,p) (see above).
amd_vnz
    The number of entries in a Householder-vector representation of the Q factor of R (but not the QR in MATLAB), after a COLAMD fill-reducing ordering. This is an upper bound on L for the LU factorization of A.
amd_rnz
    The number of entries in R for the QR factorization of A, after a COLAMD fill-reducing ordering. This is an upper bound on U for the LU factorization of A.
nblocks
    The number of blocks from dmperm (see dmperm in MATLAB).
sprank
    The structural rank of the matrix, which is the number of maximual number of nonzero entries that can be permuted to the diagonal (see dmperm, or sprank in MATLAB). This statistic is not computed for problems that represent graphs, since in those cases the diagonal of the matrix is often not relevant (self-edges are often ignored).
RBtype
    The Rutherford Boeing type of the matrix (ignoring explicit zeros in Problem.Zeros). RBtype is a a 3-letter lower-case string. The first letter is:

    r
        if A is real but not binary
    p
        if A is binary
    c
        if A is complex
    i
        if A is integer

    The second letter:

    r
        if A is rectangular
    u
        if A is square and unsymmetric
    s
        if A is symmetric (nnz(A-A.') is zero in MATLAB)
    h
        if A is Hermitian (nnz(A-A') is zero in MATLAB)
    z
        if A is skew-symmetric (nnz(A+A.') is zero in MATLAB)

    The third letter is always 'a' (for 'assembled'). The RB format allows for unassembled finite-element matrices, but they are converted to assembled format for this collection.
cholcand
    1 if the matrix is symmetric (Hermitian if complex) and if all the diagonal entries are postive and real. Zero otherwise. If 1, the matrix is a candidate for a Cholesky factorization, which MATLAB will first try when computing x=A\b.
ncc
    The number of of strongly-connected components in the graph of A (if A is square) or in the bipartite graph if A is rectangular. The diagonal is ignored.
isND
    1 if the problem comes from a 2D/3D discretization, zero otherwise. This determination is not a property of the matrix, but a qualitative assesment of the kind of problem the matrix represents.
isGraph
    1 if the problem is best considered as a graph rather than a system of equations, zero otherwise. This determination is not a property of the matrix, but a qualitative assesment of the kind of problem the matrix represents.

UFstats.csv

A CSV file is also available with some of this index information (UFstats.csv). The first line of the CSV file gives the number of matrices in the collection, and the second line gives the LastRevisionDate. Line k+2 in the file lists the following statistics for the matrix whose id number is k: Group, Name, nrows, ncols, nnz, isReal, isBinary, isND, posdef, pattern_symmetry, numerical_symmetry, and kind.
SVD-based statistics

The following statistics are not (yet) in the UFindex. They are currently available only on the web page for each matrix. You can also download the singular values at www.cise.ufl.edu/research/sparse/svd. These were typically calculated with s=svd(full(A)) in MATLAB, are are thus only available for modest-sized matrices.

norm(A)
    The 2-norm of A (the largest singular value)
min(svd(A))
    The smallest singular value
cond(A)
    The 2-norm condition number, which is the ratio of the largest over the smallest singular value.
rank(A)
    The rank of the matrix, which is the number of singular values larger than the tolerance of max(m,n)*eps(norm(A)). This tolerance is plotted in green in the figure.
sprank(A)-rank(A)
    sprank(A) (see above) is an upper bound on the rank of A.
null space dimension
    The dimension of the null space (zero if it has full numerical rank). This is simply min(m,n)-rank(A).
full numerical rank?
    'yes' or 'no'.
singular value gap
    If k=rank(A), and the matrix is rank deficient, then this is the ratio s(k)/s(k+1). A red line between the kth and (k+1)st singular value is drawn to illustrate this gap.
singular values
    These can be downloaded as a MATLAB MAT-file. Each file contains a struct with the fields: s (a vector containing the singular values), how (a string stating how the SVD was computed), and status (a string that is either 'ok' or a warning). If the status shows that the SVD did not converge, the singular values are probably not computed accurately.

=#


text:
using ErrorAnalysis
using Base.Test

# write your own tests here
@testset "Generic" begin
    x = 3.1234;
    y = 3.021;
    A = [0.531844 0.717453; 0.552965 0.421109];
    B = [0.273785 0.212329; 0.175248 0.598923];
    @test !iswithinerrorbars(x,y,0.1,true)
    @test iswithinerrorbars(x,y,0.11,true) == true
    @test !iswithinerrorbars(A,B,fill(0.01,size(A)...), true)
    @test iswithinerrorbars(A,B,fill(0.6,size(A)...), true)

end