text:
using Revise
using Printf
using PyCall
using PyPlot
using Statistics
using Bem2d

"""
    nearsurfaceerror()

Compare BEM errors from constant and quadratic elements.  The reference
model is an Okada approximation of a 2d thrust fault.  BEM errors are
evaluated not on any boundary but rather at a small distance away from
the free surface boundary.
"""
function nearsurfaceerror()
    mu = 30e9
    nu = 0.25

    # Free surface
    els = Elements(Int(1e5))
    # nfreesurf = 20
    # x1, y1, x2, y2 = discretizedline(-50e3, 0, 50e3, 0, nfreesurf)
    nfreesurf = 200
    x1, y1, x2, y2 = discretizedline(-100e3, 0, 100e3, 0, nfreesurf)
    for i in 1:length(x1)
        els.x1[els.endidx + i] = x1[i]
        els.y1[els.endidx + i] = y1[i]
        els.x2[els.endidx + i] = x2[i]
        els.y2[els.endidx + i] = y2[i]
        els.name[els.endidx + i] = "freesurf"
    end
    standardize_elements!(els)

    @show "element width"
    @show x2[1] - x1[1]

    x1, y1, x2, y2 = discretizedline(-100e3, 0, 100e3, 0, 3 * nfreesurf)
    for i in 1:length(x1)
        els.x1[els.endidx + i] = x1[i]
        els.y1[els.endidx + i] = y1[i]
        els.x2[els.endidx + i] = x2[i]
        els.y2[els.endidx + i] = y2[i]
        els.name[els.endidx + i] = "freesurf2"
    end
    standardize_elements!(els)

    # 45 degree dipping fault
    nfault = 1
    x1, y1, x2, y2 = discretizedline(-10e3, -10e3, 0, 0, nfault)
    y1 = y1 .- 5e3
    y2 = y2 .- 5e3

    for i in 1:length(x1)
        els.x1[els.endidx + i] = x1[i]
        els.y1[els.endidx + i] = y1[i]
        els.x2[els.endidx + i] = x2[i]
        els.y2[els.endidx + i] = y2[i]
        els.name[els.endidx + i] = "fault"
    end
    standardize_elements!(els)

    # Convience data structures
    idx = getidxdict(els)
    partialsconst = initpartials(els)
    partialsconst2 = initpartials(els)
    partialsquad = initpartials(els)

    # Constant slip fault
    partialsconst["disp"]["fault"]["freesurf"], _, partialsconst["trac"]["fault"]["freesurf"] = partialsconstdispstress(slip2dispstress, els, idx["fault"], idx["freesurf"], mu, nu)
    partialsconst["disp"]["freesurf"]["freesurf"], _, partialsconst["trac"]["freesurf"]["freesurf"] = partialsconstdispstress(slip2dispstress, els, idx["freesurf"], idx["freesurf"], mu, nu)
    faultslipconst = sqrt(2) / 2 * [1 ; 1]
    dispfullspaceconst = partialsconst["disp"]["fault"]["freesurf"] * faultslipconst
    dispfreesurfaceconst = inv(partialsconst["trac"]["freesurf"]["freesurf"]) * (partialsconst["trac"]["fault"]["freesurf"] * faultslipconst)
    xplotconst = els.xcenter[idx["freesurf"]]

    # HIGH RES: Constant slip fault
    partialsconst2["disp"]["fault"]["freesurf2"], _, partialsconst2["trac"]["fault"]["freesurf2"] = partialsconstdispstress(slip2dispstress, els, idx["fault"], idx["freesurf2"], mu, nu)
    partialsconst2["disp"]["freesurf2"]["freesurf2"], _, partialsconst2["trac"]["freesurf2"]["freesurf2"] = partialsconstdispstress(slip2dispstress, els, idx["freesurf2"], idx["freesurf2"], mu, nu)
    faultslipconst2 = sqrt(2) / 2 * [1 ; 1]
    dispfullspaceconst2 = partialsconst2["disp"]["fault"]["freesurf2"] * faultslipconst2
    dispfreesurfaceconst2 = inv(partialsconst2["trac"]["freesurf2"]["freesurf2"]) * (partialsconst2["trac"]["fault"]["freesurf2"] * faultslipconst2)
    xplotconst2 = els.xcenter[idx["freesurf2"]]

    # Quadratic slip fault
    partialsquad["disp"]["fault"]["freesurf"], _, partialsquad["trac"]["fault"]["freesurf"] = partialsquaddispstress(slip2dispstress, els, idx["fault"], idx["freesurf"], mu, nu)
    partialsquad["disp"]["freesurf"]["freesurf"], _, partialsquad["trac"]["freesurf"]["freesurf"] = partialsquaddispstress(slip2dispstress, els, idx["freesurf"], idx["freesurf"], mu, nu)
    xplotquad = sort(els.xnodes[idx["freesurf"], :][:])
    faultslipquad = sqrt(2) ./ 2 * ones(6)
    dispfullspacequad = partialsquad["disp"]["fault"]["freesurf"] * faultslipquad
    dispfreesurfacequad = inv(partialsquad["trac"]["freesurf"]["freesurf"]) * (partialsquad["trac"]["fault"]["freesurf"] * faultslipquad)

    # Okada solution
    ow = pyimport("okada_wrapper")# from okada_wrapper import dc3dwrapper
    xokada = collect(LinRange(-2e3, 2e3, 10000))
    offset = abs(els.xnodes[1, 2] - els.xnodes[1, 1])
    offset = 100
    @show offset
    yokada = -offset * ones(size(xokada))
    dispxokada = zeros(length(xokada))
    dispyokada = zeros(length(xokada))
    stressxxokada = zeros(length(xokada))
    stressyyokada = zeros(length(xokada))
    stressxyokada = zeros(length(xokada))

    for i in 1:length(xokada)
        # Fault dipping at 45 degrees
        _, u, s = ow.dc3dwrapper(
            2.0 / 3.0,
            [0, xokada[i] + 5e3, yokada[i]],
            # 5e3,
            10e3,
            45,
            [-1e7, 1e7],
            [-10e3*sqrt(2) / 2, 10e3*sqrt(2) / 2],
            [0.0, 1.0, 0.0],
        )
        dispxokada[i] = u[2]
        dispyokada[i] = u[3]
        dgtxx = s[2, 2]
        dgtyy = s[3, 3]
        dgtxy = s[3, 2]
        dgtyx = s[2, 3]
        exx = dgtxx
        eyy = dgtyy
        exy = 0.5 * (dgtyx + dgtxy)
        sxx = mu * (exx + eyy) + 2 * mu * exx
        syy = mu * (exx + eyy) + 2 * mu * eyy
        sxy = 2 * mu * exy
        stressxxokada[i] = sxx
        stressyyokada[i] = syy
        stressxyokada[i] = sxy
    end

    # Off-fault displacements and stresses
    dispfaultconstvol, stressfaultconstvol = constdispstress(slip2dispstress, xokada, yokada, els, idx["fault"], faultslipconst[1:2:end], faultslipconst[2:2:end], mu, nu)
    dispfreesurfaceconstvol, stressfreesurfaceconstvol = constdispstress(slip2dispstress, xokada, yokada, els, idx["freesurf"], dispfreesurfaceconst[1:2:end], dispfreesurfaceconst[2:2:end], mu, nu)
    dispconst = dispfaultconstvol - dispfreesurfaceconstvol # Note negative sign
    stressconst = stressfaultconstvol - stressfreesurfaceconstvol # Note negative sign

    # HIGH RES: Off-fault displacements and stresses
    dispfaultconstvol2, stressfaultconstvol2 = constdispstress(slip2dispstress, xokada, yokada, els, idx["fault"], faultslipconst[1:2:end], faultslipconst[2:2:end], mu, nu)
    dispfreesurfaceconstvol2, stressfreesurfaceconstvol2 = constdispstress(slip2dispstress, xokada, yokada, els, idx["freesurf2"], dispfreesurfaceconst2[1:2:end], dispfreesurfaceconst2[2:2:end], mu, nu)
    dispconst2 = dispfaultconstvol2 - dispfreesurfaceconstvol2 # Note negative sign
    stressconst2 = stressfaultconstvol2 - stressfreesurfaceconstvol2 # Note negative sign

    qux = transpose(reshape(dispfreesurfacequad[1:2:end], 3, nfreesurf))
    quy = transpose(reshape(dispfreesurfacequad[2:2:end], 3, nfreesurf))
    dispfaultquadvol, stressfaultquadvol = quaddispstress(slip2dispstress, xokada, yokada, els, idx["fault"], transpose(faultslipquad[1:2:end]), transpose(faultslipquad[2:2:end]), mu, nu)
    dispfreesurfacequadvol, stressfreesurfacequadvol = quaddispstress(slip2dispstress, xokada, yokada, els, idx["freesurf"], qux, quy, mu, nu)
    dispquad = dispfaultquadvol - dispfreesurfacequadvol # Note negative sign
    stressquad = stressfaultquadvol - stressfreesurfacequadvol # Note negative sign

    # Plot ux and uy profiles
    xmin = -2e3
    xmax = 2e3
    fontsize = 24
    fontsizelegend = 20
    markersize = 15
    linewidth = 2.0
    close("all")
    figure(figsize = (30, 20))

    ax = subplot(3, 3, 1)
    plot(xokada, log10.(abs.(stressconst[:, 1])), "-c", linewidth=linewidth, label="CS BEM")
    plot(xokada, log10.(abs.(stressconst2[:, 1])), "-m", linewidth=linewidth, label="3 x CS BEM")
    plot(xokada, log10.(abs.(stressquad[:, 1])), "-r", linewidth=linewidth, label="3QN BEM")
    plot(xokada, log10.(abs.(stressxxokada)), "--k", linewidth=3.0, label="analytic")
    gca().set_xlim([xmin, xmax])
    gca().set_ylim([2, 8])
    gca().set_xticks([xmin, 0, xmax])
    # gca().set_yticks([0, 4, 8])
    legend(fontsize=fontsizelegend, frameon=true, facecolor="white", framealpha=1.0, loc=1)
    ax.tick_params("both", labelsize = fontsize)
    ylabel(L"$\log \, |\sigma_{xx}|$ (Pa)", fontsize=fontsize)

    ax = subplot(3, 3, 2)
    plot(xokada, log10.(abs.(stressconst[:, 2])), "-c", linewidth=linewidth, label="CS BEM")
    plot(xokada, log10.(abs.(stressconst2[:, 2])), "-m", linewidth=linewidth, label="3 x CS BEM")
    plot(xokada, log10.(abs.(stressquad[:, 2])), "-r", linewidth=linewidth, label="3QN BEM")
    plot(xokada, log10.(abs.(stressyyokada[:, 1])), "--k", linewidth=3.0, label="analytic")
    gca().set_xlim([xmin, xmax])
    gca().set_ylim([2, 8])
    gca().set_xticks([xmin, 0, xmax])
    # gca().set_yticks([0, 4, 8])
    legend(fontsize=fontsizelegend, frameon=true, facecolor="white", framealpha=1.0, loc=1)
    ax.tick_params("both", labelsize = fontsize)
    ylabel(L"$\log \, |\sigma_{yy}|$ (Pa)", fontsize=fontsize)

    ax = subplot(3, 3, 3)
    plot(xokada, log10.(abs.(stressconst[:, 3])), "-c", linewidth=linewidth, label="CS BEM")
    plot(xokada, log10.(abs.(stressconst2[:, 3])), "-m", linewidth=linewidth, label="3 x CS BEM")
    plot(xokada, log10.(abs.(stressquad[:, 3])), "-r", linewidth=linewidth, label="3QN BEM")
    plot(xokada, log10.(abs.(stressxyokada[:, 1])), "--k", linewidth=3.0, label="analytic")
    gca().set_xlim([xmin, xmax])
    gca().set_ylim([2, 8])
    gca().set_xticks([xmin, 0, xmax])
    # gca().set_yticks([0, 4, 8])
    legend(fontsize=fontsizelegend,frameon=true, facecolor="white", framealpha=1.0, loc=1)
    ax.tick_params("both", labelsize = fontsize)
    ylabel(L"$\log \, |\sigma_{xy}|$ (Pa)", fontsize=fontsize)

    # Absolute error
    ax = subplot(3, 3, 4)
    stressxxconsterror = abs.(stressconst[:, 1] - stressxxokada[:, 1])
    stressxxconsterror2 = abs.(stressconst2[:, 1] - stressxxokada[:, 1])
    stressxxquaderror = abs.(stressquad[:, 1] - stressxxokada[:, 1])
    plot(xokada, log10.(abs.(stressxxconsterror)), "-c", linewidth=linewidth, label=@sprintf "CS BEM mean error = %05.2f" mean(abs.(stressxxconsterror)))
    plot(xokada, log10.(abs.(stressxxconsterror2)), "-m", linewidth=linewidth, label=@sprintf "3 x CS BEM mean error = %05.2f" mean(abs.(stressxxconsterror2)))
    plot(xokada, log10.(abs.(stressxxquaderror)), "-r", linewidth=linewidth, label=@sprintf "3QN BEM mean error = %05.2f" mean(abs.(stressxxquaderror)))
    gca().set_xlim([xmin, xmax])
    gca().set_ylim([2, 8])
    gca().set_xticks([xmin, 0, xmax])
    # gca().set_yticks([-3, 0, 3, 6])
    legend(fontsize=fontsizelegend, frameon=true, facecolor="white", framealpha=1.0, loc=1)
    ax.tick_params("both", labelsize = fontsize)
    ylabel(L"$|\sigma_{xx}$ error|", fontsize=fontsize)

    ax = subplot(3, 3, 5)
    stressyyconsterror = abs.(stressconst[:, 2] - stressyyokada[:, 1])
    stressyyconsterror2 = abs.(stressconst2[:, 2] - stressyyokada[:, 1])
    stressyyquaderror = abs.(stressquad[:, 2] - stressyyokada[:, 1])
    plot(xokada, log10.(abs.(stressyyconsterror)), "-c", linewidth=linewidth, label=@sprintf "CS BEM mean error = %05.2f" mean(abs.(stressyyconsterror)))
    plot(xokada, log10.(abs.(stressyyconsterror2)), "-m", linewidth=linewidth, label=@sprintf "3 x CS BEM mean error = %05.2f" mean(abs.(stressyyconsterror2)))
    plot(xokada, log10.(abs.(stressyyquaderror)), "-r", linewidth=linewidth, label=@sprintf "3QN BEM mean error = %05.2f" mean(abs.(stressyyquaderror)))
    gca().set_xlim([xmin, xmax])
    gca().set_ylim([2, 8])
    gca().set_xticks([xmin, 0, xmax])
    # gca().set_yticks([-3, 0, 3, 6])
    lh = legend(fontsize=fontsizelegend, frameon=true, facecolor="white", framealpha=1.0, loc=1)
    ax.tick_params("both", labelsize = fontsize)
    ylabel(L"$|\sigma_{yy}$ error|", fontsize=fontsize)

    ax = subplot(3, 3, 6)
    stressxyconsterror = abs.(stressconst[:, 3] - stressxyokada[:, 1])
    stressxyconsterror2 = abs.(stressconst2[:, 3] - stressxyokada[:, 1])
    stressxyquaderror = abs.(stressquad[:, 3] - stressxyokada[:, 1])
    plot(xokada, log10.(abs.(stressxyconsterror)), "-c", linewidth=linewidth, label=@sprintf "CS BEM mean error = %05.2f" mean(abs.(stressxyconsterror)))
    plot(xokada, log10.(abs.(stressxyconsterror2)), "-m", linewidth=linewidth, label=@sprintf "3 x CS BEM mean error = %05.2f" mean(abs.(stressxyconsterror2)))
    plot(xokada, log10.(abs.(stressxyquaderror)), "-r", linewidth=linewidth, label=@sprintf "3QN BEM mean error = %05.2f" mean(abs.(stressxyquaderror)))
    gca().set_xlim([xmin, xmax])
    gca().set_ylim([2, 8])
    gca().set_xticks([xmin, 0, xmax])
    # gca().set_yticks([-3, 0, 3, 6])
    legend(fontsize=fontsizelegend, frameon=true, facecolor="white", framealpha=1.0, loc=1)
    ax.tick_params("both", labelsize = fontsize)
    ylabel(L"$|\sigma_{xy}$ error|", fontsize=fontsize)

    ax = subplot(3, 3, 7)
    stressxxconsterror = 100 * (stressconst[:, 1] - stressxxokada[:, 1]) ./ stressxxokada[:, 1]
    stressxxconsterror2 = 100 * (stressconst2[:, 1] - stressxxokada[:, 1]) ./ stressxxokada[:, 1]
    stressxxquaderror = 100 * (stressquad[:, 1] - stressxxokada[:, 1]) ./ stressxxokada[:, 1]
    plot(xokada, log10.(abs.(stressxxconsterror)), "-c", linewidth=linewidth, label=@sprintf "CS BEM median %% error = %05.2f" median(abs.(stressxxconsterror)))
    plot(xokada, log10.(abs.(stressxxconsterror2)), "-m", linewidth=linewidth, label=@sprintf "3 x CS BEM median %% error = %05.2f" median(abs.(stressxxconsterror2)))
    plot(xokada, log10.(abs.(stressxxquaderror)), "-r", linewidth=linewidth, label=@sprintf "3QN BEM median %% error = %05.2f" median(abs.(stressxxquaderror)))
    gca().set_xlim([xmin, xmax])
    gca().set_ylim([-3, 6])
    gca().set_xticks([xmin, 0, xmax])
    gca().set_yticks([-3, 0, 3, 6])
    legend(fontsize=fontsizelegend, frameon=true, facecolor="white", framealpha=1.0, loc=1)
    ax.tick_params("both", labelsize = fontsize)
    xlabel(L"$x$ (m)", fontsize=fontsize)
    ylabel(L"$\log \, |\sigma_{xx} \, \%$ error|", fontsize=fontsize)

    ax = subplot(3, 3, 8)
    stressyyconsterror = 100 * (stressconst[:, 2] - stressyyokada[:, 1]) ./ stressyyokada[:, 1]
    stressyyconsterror2 = 100 * (stressconst2[:, 2] - stressyyokada[:, 1]) ./ stressyyokada[:, 1]
    stressyyquaderror = 100 * (stressquad[:, 2] - stressyyokada[:, 1]) ./ stressyyokada[:, 1]
    plot(xokada, log10.(abs.(stressyyconsterror)), "-c", linewidth=linewidth, label=@sprintf "CS BEM median %% error = %05.2f" median(abs.(stressyyconsterror)))
    plot(xokada, log10.(abs.(stressyyconsterror2)), "-m", linewidth=linewidth, label=@sprintf "3 x CS BEM median %% error = %05.2f" median(abs.(stressyyconsterror2)))
    plot(xokada, log10.(abs.(stressyyquaderror)), "-r", linewidth=linewidth, label=@sprintf "3QN BEM median %% error = %05.2f" median(abs.(stressyyquaderror)))
    gca().set_xlim([xmin, xmax])
    gca().set_ylim([-3, 6])
    gca().set_xticks([xmin, 0, xmax])
    gca().set_yticks([-3, 0, 3, 6])
    lh = legend(fontsize=fontsizelegend, frameon=true, facecolor="white", framealpha=1.0, loc=1)
    ax.tick_params("both", labelsize = fontsize)
    xlabel(L"$x$ (m)", fontsize=fontsize);
    ylabel(L"$\log \, |\sigma_{yy} \, \%$ error|", fontsize=fontsize)

    ax = subplot(3, 3, 9)
    stressxyconsterror = 100 * (stressconst[:, 3] - stressxyokada[:, 1]) ./ stressxyokada[:, 1]
    stressxyconsterror2 = 100 * (stressconst2[:, 3] - stressxyokada[:, 1]) ./ stressxyokada[:, 1]
    stressxyquaderror = 100 * (stressquad[:, 3] - stressxyokada[:, 1]) ./ stressxyokada[:, 1]
    plot(xokada, log10.(abs.(stressxyconsterror)), "-c", linewidth=linewidth, label=@sprintf "CS BEM median %% error = %05.2f" median(abs.(stressxyconsterror)))
    plot(xokada, log10.(abs.(stressxyconsterror2)), "-m", linewidth=linewidth, label=@sprintf "3 x CS BEM median %% error = %05.2f" median(abs.(stressxyconsterror2)))
    plot(xokada, log10.(abs.(stressxyquaderror)), "-r", linewidth=linewidth, label=@sprintf "3QN BEM median %% error = %05.2f" median(abs.(stressxyquaderror)))
    gca().set_xlim([xmin, xmax])
    gca().set_ylim([-3, 6])
    gca().set_xticks([xmin, 0, xmax])
    gca().set_yticks([-3, 0, 3, 6])
    legend(fontsize=fontsizelegend, frameon=true, facecolor="white", framealpha=1.0, loc=1)
    ax.tick_params("both", labelsize = fontsize)
    xlabel(L"$x$ (m)", fontsize=fontsize);
    ylabel(L"$\log \, |\sigma_{xy} \, \%$ error|", fontsize=fontsize)

    tight_layout(pad=5.0)

    show()
end
nearsurfaceerror()


text:
struct AhoCorasickAutomaton{T <: Unsigned}
    base::Vector{T}
    from::Vector{T}
    ikey::Vector{T}
    deep::Vector{T}
    back::Vector{T}
    arcs::Vector{Vector{UInt8}}
end

function AhoCorasickAutomaton{T}() where T
    base = T[1]
    from = T[1]
    ikey = T[0]
    deep = T[0]
    back = T[0]
    arcs = [UInt8[]]
    return AhoCorasickAutomaton{T}(base, from, ikey, deep, back, arcs)
end

function AhoCorasickAutomaton{T}(keys::Vector, sort::Bool) where T
        obj = AhoCorasickAutomaton{T}()
        if sort && !issorted(keys) Base.sort!(keys) end
        for (key, i) in keys
                addkey!(obj, codeunits(key), T(i))
        end
        shrink!(obj)
        @inbounds fillback!(obj)
        validate(obj)
        resize!(obj.arcs, 0)
        return obj
end

function AhoCorasickAutomaton{T}(keys::AbstractDict{String, Ti}; sort::Bool = true) where {T, Ti}
        return AhoCorasickAutomaton{T}(collect(keys), sort)
end

function AhoCorasickAutomaton{T}(keys::Vector{String}; sort::Bool = true) where T
    return AhoCorasickAutomaton{T}(collect(zip(keys, 1:length(keys))), sort)
end

function shrink!(obj::AhoCorasickAutomaton{T})::T where T
    actlen = findlast(!iszero, obj.from)
    if (actlen < length(obj.from))
        resize!(obj.base, actlen)
        resize!(obj.from, actlen)
        resize!(obj.ikey, actlen)
        resize!(obj.deep, actlen)
        resize!(obj.back, actlen)
        resize!(obj.arcs, actlen)
    end
    return actlen
end

function enlarge!(obj::AhoCorasickAutomaton{T}, newlen::T)::T where T
    base = obj.base; from = obj.from; deep = obj.deep; back = obj.back; ikey = obj.ikey; arcs = obj.arcs;
    oldlen::T = length(obj.base)
    newlen2 = oldlen
    while newlen2 < newlen newlen2 *= 2 end
    if (oldlen < newlen2)
        resize!(base, newlen2)
        resize!(from, newlen2)
        resize!(ikey, newlen2)
        resize!(deep, newlen2)
        resize!(back, newlen2)
        resize!(arcs, newlen2)
        # for i = oldlen + 1:newlen2 base[i] = i end
        base[oldlen + 1:newlen2] .= 0
        from[oldlen + 1:newlen2] .= 0
        ikey[oldlen + 1:newlen2] .= 0
        deep[oldlen + 1:newlen2] .= 0
        back[oldlen + 1:newlen2] .= 0
        for i in oldlen + 1:newlen2
            arcs[i] = UInt8[]
        end
    end
    return newlen2
end

function addkey!(obj::AhoCorasickAutomaton{T}, code::Base.CodeUnits{UInt8,String}, icode::T)::Nothing where T
    base = obj.base; from = obj.from; deep = obj.deep; back = obj.back; ikey = obj.ikey; arcs = obj.arcs;
    cur::T = 1
    nxt::T = 0
    for c in code
        nxt = base[cur] + c
        enlarge!(obj, nxt)
        if (from[nxt] == 0)
            from[nxt] = cur
            push!(arcs[cur], c)
            deep[nxt] = deep[cur] + 1
            cur = nxt
        elseif (from[nxt] == cur)
            cur = nxt
        else # from[nxt] != cur
            push!(arcs[cur], c)
            if length(arcs[cur]) <= length(arcs[from[nxt]]) || from[nxt] == from[cur]
                rebase!(obj, cur)
                nxt = base[cur] + c
            else
                rebase!(obj, from[nxt])
            end
            from[nxt] = cur
            deep[nxt] = deep[cur] + 1
            cur = nxt
        end
    end
    ikey[cur] = icode
    return nothing
end

function rebase!(obj::AhoCorasickAutomaton{T}, cur::T)::Nothing where T
    base = obj.base; from = obj.from; deep = obj.deep; back = obj.back; ikey = obj.ikey; arcs = obj.arcs;
    oldbase = base[cur]
    @assert length(arcs[cur]) > 0 string(cur)
    newbase = findbase(obj, cur)
    enlarge!(obj, newbase + maximum(arcs[cur]))
    for i = eachindex(arcs[cur])
        # arc = arcs[cur][i]
        newson = newbase + arcs[cur][i]
        from[newson] = cur
        oldson = oldbase + arcs[cur][i]
        if (from[oldson] != cur) continue end
        base[newson] = base[oldson]
        ikey[newson] = ikey[oldson]
        deep[newson] = deep[oldson]
        z = arcs[newson]; arcs[newson] = arcs[oldson]; arcs[oldson] = z;
        # grandsons
        for arc in arcs[newson]
            from[base[newson] + arc] = newson
        end
        # oldson
        base[oldson] = from[oldson] = ikey[oldson] = deep[oldson] = 0
    end
    base[cur] = newbase
    return nothing
end

function findbase(obj::AhoCorasickAutomaton{T}, cur::T)::T where T
    base = obj.base; from = obj.from; deep = obj.deep; back = obj.back; ikey = obj.ikey; arcs = obj.arcs;
    n::T = length(from)
    for b = max(cur + 1, base[cur]):n
        ok = true
        for arc in arcs[cur]
            @inbounds ok &= arc + b > n || from[arc + b] == 0
        end
        if (ok)
            return b
        end
    end
    return T(n + 1)
end

function fillback!(obj::AhoCorasickAutomaton{T})::Nothing where T
    base = obj.base; from = obj.from; deep = obj.deep; back = obj.back; ikey = obj.ikey; arcs = obj.arcs;
    #println(arcs)
    n::T = length(arcs)
    root::T = 1
    que = similar(base); head::T = 1; tail::T = 2;
    que[1] = root; back[root] = root;
    while head < tail
        cur = que[head]; head += 1;
        for arc in arcs[cur]
            chd = base[cur] + arc
            chdback = root
            if (cur != root)
                chdback = back[cur]
                while chdback != root && (base[chdback] + arc > n || from[base[chdback] + arc] != chdback)
                    chdback = back[chdback]
                end
                if base[chdback] + arc <= n && from[base[chdback] + arc] == chdback
                    chdback = base[chdback] + arc
                end
            end
            back[chd] = chdback
            que[tail] = chd; tail += 1;
        end
    end
    return nothing
end

function validate(obj::AhoCorasickAutomaton{T})::Nothing where T
    base = obj.base; from = obj.from; deep = obj.deep; back = obj.back; ikey = obj.ikey; arcs = obj.arcs;
    root = 1
    que = similar(base); head = 1; tail = 2;
    que[1] = root;
    while head < tail
        cur = que[head]; head += 1;
        for arc in arcs[cur]
            chd = base[cur] + arc
            # @assert from[chd] == cur && back[chd] != chd && back[chd] != 0 string(chd, " fa=", cur, " from=", from[chd], " back=", back[chd])
            @assert from[chd] == cur string("cur=", cur, " chd=", chd, " from=", from[chd])
            que[tail] = chd; tail += 1;
        end
    end
    return nothing
end

"""
ACMatch has 3 fields:

    1. s : start of match
    2. t : stop of match, [s, t), using str[s:prevind(str, t)] to get matched patterns
    3. i : index of the key in *obj*, which is the original insertion order of keys to *obj*

The field *i* may be use as index of external property arrays, i.e., the AhoCorasickAutomaton
can act as a `Map{String, Any}`.
"""
struct ACMatch
    s::Int
    t::Int
    i::Int
end

import Base.length
length(x::ACMatch) = x.t - x.s

function isless(x::ACMatch, y::ACMatch)::Bool
    return x.s < y.s || x.s == y.s && x.t < y.t || x.s == y.s && x.t == y.t && x.i < y.i
end

function eachmatch(obj::AhoCorasickAutomaton{T}, text::AbstractString)::Vector{ACMatch} where T
    return eachmatch(obj, codeunits(text))
end

function eachmatch(obj::AhoCorasickAutomaton{T}, codes::DenseVector{UInt8})::Vector{ACMatch} where T
    base = obj.base; from = obj.from; deep = obj.deep; back = obj.back; ikey = obj.ikey; arcs = obj.arcs;
    n = length(base)
    root = cur = T(1)
    res = ACMatch[]
    for i = 1:length(codes)
        c = codes[i]
        while cur != root && (base[cur] + c > n || from[base[cur] + c ] != cur)
            cur = back[cur]
        end
        if (base[cur] + c <= n && from[base[cur] + c] == cur)
            cur = base[cur] + c
        end
        # if (ikey[cur] > 0)
            node = cur
            while node != root
                if (ikey[node] > 0) push!(res, ACMatch(i + 1 - deep[node], i + 1, ikey[node])) end
                node = back[node]
            end
        # end
    end
    return res
end

#### utils
import Base.getindex
getindex(xs::String, match::ACMatch) = xs[match.s:prevind(xs, match.t)]

function get(obj::AhoCorasickAutomaton{T}, key::DenseVector{UInt8}, default::T)::T where T
    cur::T = 1
    n::T = length(obj.from)
    for c in key
        nxt = obj.base[cur] + c
        if (nxt <= n && obj.from[nxt] == cur)
            cur = nxt
        else
            return default
        end
    end
    return obj.ikey[cur]
end

function get(obj::AhoCorasickAutomaton{T}, key::AbstractString, default::T)::T where T
    return get(obj, codeunits(key), default)
end

getindex(aca::AhoCorasickAutomaton{T}, key::AbstractString) where {T} = get(aca, key, T(0))

function in(key::AbstractString, obj::AhoCorasickAutomaton{T})::Bool where T
    return get(obj, key, T(0)) > 0
end

function in(key::DenseVector{UInt8}, obj::AhoCorasickAutomaton{T})::Bool where T
    return get(obj, key, T(0)) > 0
end

function length(obj::AhoCorasickAutomaton{T}) where T
    return count(!iszero, obj.ikey)
end

function collect(obj::AhoCorasickAutomaton{T}) where T
    base = obj.base
    from = obj.from
    ikey = obj.ikey
    res = Pair{String, Int}[]
    for i = 1:length(ikey)
        if ikey[i] == 0 continue end
        codes = UInt8[]
        j = i
        while j > 1
            c = j - base[from[j]]
            push!(codes, c)
            j = from[j]
        end
        push!(res, String(reverse!(codes)) => ikey[i])
    end
    return res
end

function keys(obj::AhoCorasickAutomaton{T}) where T
    return map(first, collect(obj))
end

function values(obj::AhoCorasickAutomaton{T}) where T
    return filter(!iszero, obj.ikey)
end

function ==(x::AhoCorasickAutomaton, y::AhoCorasickAutomaton)
    return x.base == y.base && x.from == y.from && x.ikey == y.ikey && x.deep == y.deep && x.back == y.back
end

function show(io::IO, o::AhoCorasickAutomaton{T}) where T
    print(io, "AhoCorasickAutomaton{$T}", "(keys=", length(o), ")")
end

ext: jl


text:
#! /usr/bin/julia

# Rosetta Code, Sorting algorithms/Bubble sort

function bubblesort{T}(a::AbstractArray{T,1})
    b = copy(a)
    isordered = false
    span = length(b)
    while !isordered && span > 1
        isordered = true
        for i in 2:span
            if b[i] < b[i-1]
                t = b[i]
                b[i] = b[i-1]
                b[i-1] = t
                isordered = false
            end
        end
        span -= 1
    end
    return b
end

a = [rand(-100:100) for i in 1:20]
println("Before bubblesort:")
println(a)
a = bubblesort(a)
println("\nAfter bubblesort:")
println(a)

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using InvariantMeasures
using Test

@testset "InvariantMeasures.jl" begin

    include("TestContractors.jl")
    include("TestDynamic.jl")
    include("TestHat.jl")
    include("TestUlam.jl")
    include("TestAssemble.jl")
    include("TestAssembleHat.jl")
    include("TestEstimate.jl")
    include("TestNormOfPowers.jl")
    include("TestPreimages.jl")
    include("TestChebyshev.jl")
end


text:
struct BandedBlockBandedSizes <: AbstractBlockSizes{2}
    block_sizes::DefaultBlockSizes{2}
    data_block_sizes::DefaultBlockSizes{2}
    l::Int
    u::Int
    λ::Int
    μ::Int
end

function BandedBlockBandedSizes(bs::BlockSizes{2}, l, u, λ, μ)
    # data matrix has row blocks all the same size but same column blocks
    # we access the cumul vec directly to reuse bs.cumul_sizes[2]
    d_bs = BlockSizes((BlockArrays._cumul_vec(fill(max(λ+μ+1,0), max(l+u+1,0))),bs.cumul_sizes[2]))
    BandedBlockBandedSizes(bs, d_bs, l, u, λ, μ)
end

BandedBlockBandedSizes(bs::AbstractBlockSizes{2}, l, u, λ, μ) =
    BandedBlockBandedSizes(convert(BlockSizes{2},bs), l, u, λ, μ)

BandedBlockBandedSizes(rows::AbstractVector{Int}, cols::AbstractVector{Int}, l, u, λ, μ) =
    BandedBlockBandedSizes(BlockSizes(rows,cols), l, u, λ, μ)

cumulsizes(B::BandedBlockBandedSizes) = cumulsizes(B.block_sizes)

convert(::Type{BlockSkylineSizes}, B::BandedBlockBandedSizes) =
    BlockBandedSizes(B.block_sizes, B.l, B.u)
convert(::Type{BlockBandedSizes}, B::BandedBlockBandedSizes) =
    BlockBandedSizes(B.block_sizes, B.l, B.u)

convert(::Type{BlockSizes}, B::BandedBlockBandedSizes) = B.block_sizes
convert(::Type{BlockSizes{2}}, B::BandedBlockBandedSizes) = B.block_sizes
convert(::Type{DefaultBlockSizes{2}}, B::BandedBlockBandedSizes) = B.block_sizes

const BandedOrBlockBandedSizes = Union{BandedBlockBandedSizes,BlockBandedSizes}

==(A::BandedOrBlockBandedSizes, B::BandedOrBlockBandedSizes) = A.block_sizes == B.block_sizes
==(A::BandedOrBlockBandedSizes, B::BlockSkylineSizes) = A.block_sizes == B.block_sizes
==(A::BlockSkylineSizes, B::BandedOrBlockBandedSizes) = A.block_sizes == B.block_sizes
==(A::BandedOrBlockBandedSizes, B::AbstractBlockSizes) = A.block_sizes == B
==(A::AbstractBlockSizes, B::BandedOrBlockBandedSizes) = A == B.block_sizes

BlockBandedSizes(B::BandedBlockBandedSizes) = convert(BlockBandedSizes, B)
BlockSkylineSizes(B::BandedBlockBandedSizes) = convert(BlockSkylineSizes, B)

function check_data_sizes(data::AbstractBlockMatrix, B::BandedBlockBandedSizes)
    bs = data.block_sizes
    c_rows, c_cols = bs.cumul_sizes
    if length(c_rows) ≠ B.l + B.u + 2 && !(length(c_rows) == 1 && -B.l > B.u)
        throw(ArgumentError("Data matrix must have number of row blocks equal to number of block bands"))
    end
    for k = 1:length(c_rows)-1
        if c_rows[k+1]-c_rows[k] ≠ B.λ + B.μ + 1 && !(c_rows[k+1] == c_rows[k]  && -B.λ > B.μ)
            throw(ArgumentError("Data matrix must have row block sizes equal to number of subblock bands"))
        end
    end
    if c_cols ≠ B.block_sizes.cumul_sizes[2]
        throw(ArgumentError("Data matrix must have same column blocks as matrix"))
    end
end

function _BandedBlockBandedMatrix end

# Represents a block banded matrix with banded blocks
#   similar to BandedMatrix{BandedMatrix{T,Matrix{T}}}
# Here the data is stored by blocks, in a way that is consistent with
# BandedMatrix
#

struct BandedBlockBandedMatrix{T, BLOCKS} <: AbstractBlockBandedMatrix{T}
    data::BLOCKS
    block_sizes::BandedBlockBandedSizes

    l::Int  # block lower bandwidth
    u::Int  # block upper bandwidth
    λ::Int  # sub lower bandwidth
    μ::Int  # sub upper bandwidth

    global function _BandedBlockBandedMatrix(data::AbstractBlockMatrix,
                                             block_sizes::BandedBlockBandedSizes)
      @boundscheck check_data_sizes(data, block_sizes)
      new{eltype(data), typeof(data)}(data, block_sizes, block_sizes.l,
                                      block_sizes.u, block_sizes.λ,
                                      block_sizes.μ)
    end
end

const DefaultBandedBlockBandedMatrix{T} = BandedBlockBandedMatrix{T, PseudoBlockMatrix{T, Matrix{T}}}

@inline _BandedBlockBandedMatrix(data::AbstractMatrix, block_sizes::BandedBlockBandedSizes) =
    _BandedBlockBandedMatrix(PseudoBlockArray(data, block_sizes.data_block_sizes), block_sizes)

BandedBlockBandedMatrix{T, B}(::UndefInitializer, block_sizes::BandedBlockBandedSizes) where {T, B} =
    _BandedBlockBandedMatrix(B(undef, block_sizes.data_block_sizes), block_sizes)

"""
    BandedBlockBandedMatrix{T}(undef, (rows, cols), (l, u), (λ, μ))

returns an undef `sum(rows)`×`sum(cols)` banded-block-banded matrix `A`
of type `T` with block-bandwidths `(l,u)` and where `A[Block(K,J)]`
is a `BandedMatrix{T}` of size `rows[K]`×`cols[J]` with bandwidths `(λ,μ)`.
"""
BandedBlockBandedMatrix

BandedBlockBandedMatrix{T, B}(::UndefInitializer, dims::NTuple{2, AbstractVector{Int}},
                              lu::NTuple{2, Int}, λμ::NTuple{2, Int}) where {T, B} =
    BandedBlockBandedMatrix{T, B}(undef, BandedBlockBandedSizes(dims..., lu..., λμ...))

BandedBlockBandedMatrix{T, B}(::UndefInitializer, bs::AbstractBlockSizes,
                              lu::NTuple{2, Int}, λμ::NTuple{2, Int}) where {T, B} =
    BandedBlockBandedMatrix{T, B}(undef, BandedBlockBandedSizes(bs, lu..., λμ...))

BandedBlockBandedMatrix{T}(::UndefInitializer, args...) where T =
    DefaultBandedBlockBandedMatrix{T}(undef, args...)

# Auxiliary outer constructors
@inline _BandedBlockBandedMatrix(data::AbstractMatrix, dims::NTuple{2, AbstractVector{Int}},
                                         lu::NTuple{2, Int}, λμ::NTuple{2, Int}) =
    _BandedBlockBandedMatrix(data, BandedBlockBandedSizes(dims..., lu..., λμ...))

function convert(::Type{<:BandedBlockBandedMatrix}, B::BandedMatrix)
    if isdiag(B)
        _BandedBlockBandedMatrix(copy(B.data), (fill(1,size(B,1)),fill(1,size(B,2))), (0,0), (0,0))
    else
        _BandedBlockBandedMatrix(copy(B.data), [size(B,1)], [size(B,2)], (0,0), (B.l,B.u))
    end
end

function BandedBlockBandedMatrix{T, B}(Z::Zeros, dims::NTuple{2,AbstractVector{Int}},
                                       lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T, B}
   if size(Z) ≠ sum.(dims)
       throw(DimensionMismatch())
   end

   bs = BandedBlockBandedSizes(dims..., lu..., λμ...)
   d_bs = bs.data_block_sizes
   blocks = fill!(B(undef, d_bs), zero(T))
   _BandedBlockBandedMatrix(blocks, bs)
end

function BandedBlockBandedMatrix{T, B}(E::Eye, dims::NTuple{2,AbstractVector{Int}},
                                       lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T, B}
    if size(E) ≠ sum.(dims)
        throw(DimensionMismatch())
    end
    ret = BandedBlockBandedMatrix{T, B}(Zeros{T}(E), dims, lu, λμ)
    ret[diagind(ret)] .= one(T)
    ret
end

function BandedBlockBandedMatrix{T, B}(A::UniformScaling, dims::NTuple{2, AbstractVector{Int}},
                                       lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T, B}
    ret = BandedBlockBandedMatrix{T, B}(Zeros{T}(sum.(dims)), dims, lu, λμ)
    ret[diagind(ret)] .= convert(T, A.λ)
    ret
end

BandedBlockBandedMatrix{T}(m::Union{Zeros, Eye, UniformScaling},
                           dims::NTuple{2,AbstractVector{Int}},
                           lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where T =
  DefaultBandedBlockBandedMatrix{T}(m, dims, lu, λμ)

function BandedBlockBandedMatrix{T, B}(Z::Zeros, block_sizes::BandedBlockBandedSizes) where {T, B}
   if size(Z) ≠ size(block_sizes)
       throw(DimensionMismatch())
   end

   d_bs = block_sizes.data_block_sizes
   data = fill!(B(undef, d_bs), zero(T))
   _BandedBlockBandedMatrix(data, block_sizes)
end

BandedBlockBandedMatrix{T}(Z::Zeros, block_sizes::BandedBlockBandedSizes,
                           lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where T =
  DefaultBandedBlockBandedMatrix{T}(Z, block_sizes)

function BandedBlockBandedMatrix{T, Blocks}(A::AbstractMatrix, block_sizes::BandedBlockBandedSizes) where {T, Blocks}
    ret = BandedBlockBandedMatrix{T, Blocks}(Zeros{T}(size(A)), block_sizes)
    L,M = block_sizes.λ, block_sizes.μ
    for J = Block.(1:nblocks(ret, 2)), K = blockcolrange(ret, Int(J))
        kr, jr = globalrange(block_sizes, (Int(K), Int(J)))

        # We have the correct block - now we need to only add the entries from
        # the correct bands
        B = view(A, kr, jr)
        R = view(ret, K, J)
        rows, cols = size(B)
        for λ = 1:L, j = 1:min(rows-λ, cols)
            view(R, j+λ, j) .= view(B, j+λ, j)
        end
        for i = 1:min(rows, cols)
            view(R, i, i) .= view(B, i, i)
        end
        for μ = 1:M, k = 1:min(rows, cols-μ)
            view(R, k, k+μ) .= view(B, k, k+μ)
        end
    end
    ret
end

BandedBlockBandedMatrix{T}(A::AbstractMatrix, block_sizes::BandedBlockBandedSizes) where T =
    DefaultBandedBlockBandedMatrix{T}(A, block_sizes)

BandedBlockBandedMatrix{T}(A::AbstractMatrix, block_sizes::BandedBlockBandedSizes,
                           lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where T =
  DefaultBandedBlockBandedMatrix{T}(A, block_sizes, lu, λμ)

BandedBlockBandedMatrix(A::Union{AbstractMatrix,UniformScaling},
                        block_sizes::BandedBlockBandedSizes) =
    BandedBlockBandedMatrix{eltype(A)}(A, block_sizes)

BandedBlockBandedMatrix(A::Union{AbstractMatrix,UniformScaling},
                        block_sizes::AbstractBlockSizes, lu::NTuple{2,Int},
                        λμ::NTuple{2,Int}) =
                        BandedBlockBandedMatrix(A, BandedBlockBandedSizes(block_sizes, lu..., λμ...))

BandedBlockBandedMatrix{T}(A::Union{AbstractMatrix,UniformScaling},
                           dims::NTuple{2,AbstractVector{Int}}, lu::NTuple{2,Int},
                           λμ::NTuple{2,Int}) where T =
    BandedBlockBandedMatrix{T}(A, BandedBlockBandedSizes(dims..., lu..., λμ...))

BandedBlockBandedMatrix{T,B}(A::Union{AbstractMatrix,UniformScaling},
                           dims::NTuple{2,AbstractVector{Int}}, lu::NTuple{2,Int},
                           λμ::NTuple{2,Int}) where {T,B} =
    BandedBlockBandedMatrix{T,B}(A, BandedBlockBandedSizes(dims..., lu..., λμ...))

BandedBlockBandedMatrix(A::Union{AbstractMatrix,UniformScaling},
                        dims::NTuple{2, AbstractVector{Int}}, lu::NTuple{2,Int},
                        λμ::NTuple{2,Int}) =
    BandedBlockBandedMatrix{eltype(A)}(A, dims, lu, λμ)

BandedBlockBandedMatrix(A::AbstractMatrix) =
    BandedBlockBandedMatrix(A, blocksizes(A), blockbandwidths(A), subblockbandwidths(A))

similar(A::BandedBlockBandedMatrix, T::Type=eltype(A), bs::BandedBlockBandedSizes=blocksizes(A)) =
      BandedBlockBandedMatrix{T}(undef, bs)

function sparse(A::BandedBlockBandedMatrix)
    i = Vector{Int}()
    j = Vector{Int}()
    z = Vector{eltype(A)}()
    for J = Block.(1:nblocks(A,2)), K = blockcolrange(A, J)
        B = view(A, K, J)
        ĩ = _banded_rowval(B)
        j̃ = _banded_colval(B)
        z̃ = _banded_nzval(B)
        ĩ .+= cumulsizes(A, 1, Int(K))-1
        j̃ .+= cumulsizes(A, 2, Int(J))-1
        append!(i, ĩ)
        append!(j, j̃)
        append!(z, z̃)
    end
    sparse(i, j, z, size(A)...)
end

################################
# BandedBlockBandedMatrix Interface #
################################

MemoryLayout(::Type{<:BandedBlockBandedMatrix}) = BandedBlockBandedColumnMajor()
BroadcastStyle(::Type{<:BandedBlockBandedMatrix}) = BandedBlockBandedStyle()

isbandedblockbanded(_) = false
isbandedblockbanded(::BandedBlockBandedMatrix) = true

blockbandwidths(A::BandedBlockBandedMatrix) = (A.l, A.u)

"""
    subblockbandwidths(A)

returns the sub-block bandwidths of `A`, where `A` is a banded-block-banded
matrix. In other words, `A[Block(K,J)]` will return a `BandedMatrix` with
bandwidths given by `subblockbandwidths(A)`.
"""
subblockbandwidths(A::BandedBlockBandedMatrix) = (A.λ, A.μ)

# default is to use bandwidth
subblockbandwidths(A::AbstractMatrix) = bandwidths(A)

"""
    subblockbandwidth(A, i)

returns the sub-block lower (`i == 1`) or upper (`i == 2`) bandwidth of `A`,
where `A` is a banded-block-banded matrix. In other words, `A[Block(K,J)]` will
return a `BandedMatrix` with the returned lower/upper bandwidth.
"""
subblockbandwidth(A::AbstractArray, k::Integer) = subblockbandwidths(A)[k]

isdiag(A::BandedBlockBandedMatrix) = A.λ == A.μ == A.l == A.u

################################
# AbstractBlockArray Interface #
################################

@inline blocksizes(block_array::BandedBlockBandedMatrix) = block_array.block_sizes

zeroblock(A::BandedBlockBandedMatrix, K::Int, J::Int) =
    BandedMatrix(Zeros{eltype(A)}(blocksize(A, (K, J))), (A.λ, A.μ))

@inline function getblock(A::BandedBlockBandedMatrix, K::Int, J::Int)
    @boundscheck blockcheckbounds(A, K, J)
    if -A.l ≤ J - K ≤ A.u
        convert(BandedMatrix, view(A, Block(K, J)))
    else
        zeroblock(A, K, J)
    end
end

# @inline function Base.getindex(block_arr::BlockArray{T,N}, blockindex::BlockIndex{N}) where {T,N}
#     @boundscheck checkbounds(block_arr.blocks, blockindex.I...)
#     @inbounds block = block_arr.blocks[blockindex.I...]
#     @boundscheck checkbounds(block, blockindex.α...)
#     @inbounds v = block[blockindex.α...]
#     return v
# end

###########################
# AbstractArray Interface #
###########################

# @inline function Base.similar(block_array::BandedBlockBandedMatrix{T}, ::Type{T2}) where {T,N,T2}
#     BlockArray(similar(block_array.blocks, Array{T2, N}), copy(block_array.block_sizes))
# end

Base.size(arr::BandedBlockBandedMatrix) =
    @inbounds return size(arr.block_sizes)

@inline function getindex(A::BandedBlockBandedMatrix{T}, i::Int, j::Int) where T
    @boundscheck checkbounds(A, i, j)
    bi = global2blockindex(A.block_sizes, (i, j))
    @inbounds v = view(A, Block(bi.I))[bi.α...]
    return v
end

@inline function setindex!(A::BandedBlockBandedMatrix{T}, v, i::Int, j::Int) where T
    @boundscheck checkbounds(A, i, j)
    bi = global2blockindex(A.block_sizes, (i, j))
    if -A.l ≤ bi.I[2] - bi.I[1] ≤ A.u
        V = view(A, Block(bi.I))
        @inbounds V[bi.α...] = convert(T, v)::T
    elseif !iszero(v)
        throw(BandError(A))
    end
    return v
end

## structured matrix methods ##
function Base.replace_in_print_matrix(A::BandedBlockBandedMatrix, i::Integer, j::Integer, s::AbstractString)
    bi = global2blockindex(A.block_sizes, (i, j))
    I,J = bi.I
    i,j = bi.α
    -A.l ≤ J-I ≤ A.u && -A.λ ≤ j-i ≤ A.μ ? s : Base.replace_with_centered_mark(s)
end

############
# Indexing #
############

# function _check_setblock!(block_arr::BlockArray{T, N}, v, block::NTuple{N, Int}) where {T,N}
#     for i in 1:N
#         if size(v, i) != blocksize(block_arr.block_sizes, (i, block[i]))
#             throw(DimensionMismatch(string("tried to assign $(size(v)) array to ", blocksize(block_arr, block), " block")))
#         end
#     end
# end
#
#
@inline function setblock!(A::BandedBlockBandedMatrix, v, K::Int, J::Int)
    @boundscheck blockcheckbounds(A, K, J)

    @boundscheck (bandwidth(v, 1) > A.λ || bandwidth(v, 2) > A.μ) && throw(BandError())
    V = view(A, Block(K), Block(J))
    V .= v
    return A
end
#
# @propagate_inbounds function Base.setindex!(block_array::BlockArray{T, N}, v, block_index::BlockIndex{N}) where {T,N}
#     getblock(block_array, block_index.I...)[block_index.α...] = v
# end

########
# Misc #
########

# @generated function Base.Array(block_array::BlockArray{T, N, R}) where {T,N,R}
#     # TODO: This will fail for empty block array
#     return quote
#         block_sizes = block_array.block_sizes
#         arr = similar(block_array.blocks[1], size(block_array)...)
#         @nloops $N i i->(1:nblocks(block_sizes, i)) begin
#             block_index = @ntuple $N i
#             indices = globalrange(block_sizes, block_index)
#             arr[indices...] = getblock(block_array, block_index...)
#         end
#
#         return arr
#     end
# end
#
# @generated function Base.copyto!(block_array::BlockArray{T, N, R}, arr::R) where {T,N,R <: AbstractArray}
#     return quote
#         block_sizes = block_array.block_sizes
#
#         @nloops $N i i->(1:nblocks(block_sizes, i)) begin
#             block_index = @ntuple $N i
#             indices = globalrange(block_sizes, block_index)
#             copyto!(getblock(block_array, block_index...), arr[indices...])
#         end
#
#         return block_array
#     end
# end
#
# function Base.fill!(block_array::BlockArray, v)
#     for block in block_array.blocks
#         fill!(block, v)
#     end
# end

##################
# BandedBlockBandedBlock
#
#   views of the blocks satisfy the BandedMatrix interface, and are memory-compatible
#   with BLASBandedMatrix.
##################

const BandedBlockBandedBlock{T, BLOCKS} = SubArray{T,2,BandedBlockBandedMatrix{T, BLOCKS},Tuple{BlockSlice1,BlockSlice1},false}

isbanded(::BandedBlockBandedBlock) = true
MemoryLayout(::BandedBlockBandedBlock) = BandedColumnMajor()
BroadcastStyle(::Type{<: BandedBlockBandedBlock}) = BandedStyle()

function inblockbands(V::BandedBlockBandedBlock)
    A = parent(V)
    K_sl, J_sl = parentindices(V)
    K, J = K_sl.block, J_sl.block
    -A.l ≤ Int(J-K) ≤ A.u
end

######################################
# BandedMatrix interface  for Blocks #
######################################
@inline function bandwidths(V::BandedBlockBandedBlock)
    inblockbands(V) && return subblockbandwidths(parent(V))
    (-720,-720)
end

# gives the columns of parent(V).data that encode the block
blocks(V::BandedBlockBandedBlock)::Tuple{Int,Int} = Int(first(parentindices(V)).block),
                                                    Int(last(parentindices(V)).block)

function bandeddata(V::BandedBlockBandedBlock{T}) where T
    inblockbands(V) || return Array{T}(undef, 0, size(V,2))
    A = parent(V)
    u = A.u
    K_sl, J_sl = parentindices(V)
    K, J = K_sl.block, J_sl.block
    view(A.data, u + K - J + 1, J)
end

@inline function inbands_getindex(V::BandedBlockBandedBlock, k::Int, j::Int)
    A = parent(V)
    banded_getindex(bandeddata(V), A.λ, A.μ, k, j)
end

@inline function inbands_setindex!(V::BandedBlockBandedBlock, v, k::Int, j::Int)
    A = parent(V)
    banded_setindex!(bandeddata(V), A.λ, A.μ, v, k, j)
end

# we need getindex as we use block-wise get index above
@propagate_inbounds function getindex(V::BandedBlockBandedBlock, k::Int, j::Int)
    @boundscheck checkbounds(V, k, j)
    A = parent(V)
    K,J = blocks(V)
    if -A.l ≤ J-K ≤ A.u && -A.λ ≤ j-k ≤ A.μ
        inbands_getindex(V, k, j)
    else
        zero(eltype(V))
    end
end

@propagate_inbounds function setindex!(V::BandedBlockBandedBlock, v, k::Int, j::Int)
    @boundscheck checkbounds(V, k, j)
    A = parent(V)
    K,J = blocks(V)
    if -A.l ≤ J-K ≤ A.u && -A.λ ≤ j-k ≤ A.μ
        inbands_setindex!(V, v, k, j)
    elseif iszero(v) # allow setindex for 0 datya
        v
    else
        throw(BandError(V, J-K))
    end
end

BLAS.axpy!(a::T, A::BandedBlockBandedMatrix{T}, B::BandedBlockBandedMatrix{T}) where T =
    B .= a .* A .+ B

lmul!(x::Number, A::BandedBlockBandedMatrix) = (lmul!(x, A.data); A)
rmul!(A::BandedBlockBandedMatrix, x::Number) = (rmul!(A.data, x); A)
*(x::Number, A::BandedBlockBandedMatrix) = _BandedBlockBandedMatrix(x*A.data, A.block_sizes)
*(A::BandedBlockBandedMatrix, x::Number) = _BandedBlockBandedMatrix(A.data*x, A.block_sizes)


text:
module XLA
    using ProtoBuf
    using TensorFlow

    import Base: run

    # Load generated xla protobuf definitions
    include(joinpath(dirname(@__FILE__), "gen","tensorflow.jl"))
    include(joinpath(dirname(@__FILE__), "gen","xla.jl"))
    include(joinpath(dirname(@__FILE__), "gen","xrt.jl"))
    using .xrt
    using .xla
    using .tensorflow.tpu
    const Shape = ShapeProto
    const ProgramShape = ProgramShapeProto

    include("literal.jl")
    include("xrt.jl")
    include("topology.jl")
    include("multicore.jl")
    include("hlo.jl")
    include("shape_infer.jl")
    include("execute.jl")
    include("emulate.jl")
    include("linalg.jl")
    include("explain.jl")
    include("compiler.jl")
    include("compiler_passes.jl")
    include("compiler_interface.jl")
    include("utils.jl")
    include("grad.jl")
    include("infeed.jl")
    include("flux.jl")

    function regen_proto()
        tf_path = "/home/keno/tensorflow/tensorflow/"
        xrt_proto = tf_path * "compiler/xrt/xrt.proto"
        topo_proto = tf_path * "contrib/tpu/proto/topology.proto"
        cd(@__DIR__) do
            run(ProtoBuf.protoc(`-I=/home/keno/tensorflow/ --julia_out=gen $xrt_proto $topo_proto`))
        end
    end

    export XRTArray, GenericHloOp, HloConstant, HloDot, HloSlice, ImmutableChain, XRTAllocation, XRTRemoteStruct

    # TPU Sessions/devices
    export TPUSession, all_tpu_devices

    export @tpu_compile, @tpu_dump, @tpu

    function __init__()
        TensorFlow.import_op("XRTCompile")
        TensorFlow.import_op("XRTExecute")
        TensorFlow.import_op("XRTAllocate")
        TensorFlow.import_op("XRTReadLiteral")
        TensorFlow.import_op("XRTReadLiteralAndRelease")
        TensorFlow.import_op("XRTReleaseAllocationHandle")
        TensorFlow.import_op("XRTReleaseCompilationHandle")
        TensorFlow.import_op("ConfigureDistributedTPU")
        TensorFlow.import_op("ShutdownDistributedTPU")
    end

end


text:
type Set_Data{T1,T2,T3}
    X :: Vector{T1}
    Y :: Vector{T2}
    Z :: Vector{T3}
    var_idx :: Dict{Tuple{Int32,Int32,Int32},Int32}
    Set_Data() = new([],[],[],Dict{Tuple{T1,T2,T3},Int32}())
end

to_yz{I<:Integer              }(d::Set_Data, y::I, z::I)         :: Int32    = Int32( y*d.length(Z) + z )
varidx{I₁<:Integer,I₂<:Integer}(d::Set_Data, x::I₁,yz::I₂)       :: Int32    = get( d.var_idx, (x,yz), 0)

function Set_Data{T1,T2,T3,FLOAT}(pdf::Dict{Tuple{T1,T2,T3},FLOAT})
    local var_list::Vector{ Pair{Tuple{Int32,Int32},Int32} }
    local X_set::Set{T1}
    local Y_set::Set{T2}
    local Z_set::Set{T3}

    for (xyz,val) ∈ pdf
        if val > 0
            x,y,z = xyz
            push!(X_set,x)
            push!(Y_set,y)
            push!(Z_set,z)
        end
    end

    d = Set_Data{T1,T2,T3}
    d.X = Vector{T1}( collect(X_set) )
    d.Y = Vector{T1}( collect(Y_set) )
    d.Z = Vector{T1}( collect(Z_set) )

    CONTINUE HERE
    m_XY::Array{FLOAT} = [  sum(view(pdf[x,y,:]))   for x ∈ d.X, y ∈ d.Y  ]
    m_XZ::Array{FLOAT} = [  sum(view(pdf[x,:,z]))   for x ∈ d.X, z ∈ d.Z  ]

    ...................................................................................................

    local counter::Int = 0
            counter += 1

    sizehint!(d.var_idx, counter)

    sort!(d.X)
    sort!(d.Y)
    sort!(d.Z)

    # We want an X x (YxZ) in column major:
    for i_y,y in enumerate(d.Y):
        for i_z,z in enumerate(d.Z):
            for i_x,x in enumerate(d.X):
                if get(pdf,(x,y,z),-1) > 0
                    push!(var_list, (Int32(i_x),to_yz(i_y,i_z)) => Int32(1+length(var_list)) )
                end
            end
        end
    end

    d.var_idx = Dict{Tuple{Int32,Int32},Int32}( var_list )

    return d
    ;
end #^ Set_Data() [constructor]

max_stars_repo_name: dot-at/newpid

max_issues_repo_name: dot-at/newpid

max_forks_repo_name: dot-at/newpid


text:
# Pulse Parameters
vmax = 5  # volts
fc = 5e9  # GHz
bw = 50  # percent
v0 = 0.00001  # Maximum voltage allowed at t = 0
ptype = "normal"

# Time vector
t = linspace(0, 5e-9, 2001);

# Compute sigma (control's pulse width)
sigma = 1/(2*pi*(fc*bw/100)/2.354820045030949);

# Compute mu (time offset)
# Note that this is only approximate for the Differentiated Pulse.
mu = sigma*sqrt(-2*log(v0/vmax));

if ptype == "normal"
    pulse = vmax*exp.(-(t - mu).^2/(2*sigma^2));
elseif ptype == "diff"
    pulse = -vmax/sigma*(t - mu).*exp.(0.5-(t - mu).^2/(2*sigma^2));
elseif ptype == "sinmod"
    pulse = vmax*exp.(-(t - mu).^2/(2*sigma^2)).*cos.(2*pi*fc*(t - mu));
else
    error("Unsupported pulse type")
end

using PyPlot
plot(t*1e9, pulse, linewidth=2.0)
xlabel("Time (ns)")
ylabel("Amplitude (Volts)")
grid(true)

max_stars_repo_name: engineeringterminal/gaussian-pulses

max_issues_repo_name: engineeringterminal/gaussian-pulses

max_forks_repo_name: engineeringterminal/gaussian-pulses


text:
using myexample
using Test

@testset "myexample.jl" begin
    @test my_f(2,3)==13
    @test my_f(1,2)==8

end

@testset "Derivative Tests" begin
    @test my_derivative_f(2,1)==2
end


text:
ia> [RpcApi.SimpleOrganism(0x000000000004a6e9,:(function (o::Creature.Organism,)
     function func_32(var_30::Int8=-88,var_31::Bool=true)
         local var_34::String = "HmoFFpog"
         return var_30
     end
     function func_19(var_17::String="A06aRuJT",var_18::Int64=-4494005246402814635)
         local var_35::Int8 = -121
         return var_17
     end
     local var_1::Int8 = 121
     for i = 1:var_1
         Creature.stepRight(o)
         Creature.stepRight(o)
     end
     Creature.stepDown(o)
     Creature.eatUp(o,Int(var_1))
     for i = 1:var_1
         local var_15::Bool = true
         local var_9::Bool = true
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatLeft(o,Int(var_6))
         Creature.stepRight(o)
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.stepRight(o)
         var_9 = var_9 | var_9
         Creature.eatRight(o,Int(var_6))
         Creature.stepLeft(o)
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatUp(o,Int(var_6))
         Creature.stepRight(o)
         Creature.stepRight(o)
         Creature.eatRight(o,Int(var_6))
         Creature.eatDown(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatUp(o,Int(var_6))
     end
     Creature.stepUp(o)
     Creature.eatDown(o,Int(var_1))
     for i = 1:var_1
         Creature.stepLeft(o)
     end
     for i = 1:var_1
         Creature.stepLeft(o)
     end
     for i = 1:var_1
         Creature.stepLeft(o)
         Creature.stepLeft(o)
         Creature.stepLeft(o)
     end
     Creature.eatDown(o,Int(var_1))
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     for i = 1:var_1
         local var_25::Bool = true
     end
     for i = 1:var_1
         local var_29::Int64 = -5975265535146971507
         Creature.stepLeft(o)
     end
     return true
 end),[100,300,95,0,1,1,1],1,6745,1,5281227,0x011a2250,Dict(23832=>23832),[414,333],51),RpcApi.SimpleOrganism(0x000000000004a714,:(function (o::Creature.Organism,)
     function func_19(var_17::String="A06aRuJT",var_18::Int64=-4494005246402814635)
         Creature.stepLeft(o)
         return var_17
     end
     local var_1::Int8 = 121
     for i = 1:var_1
         Creature.stepRight(o)
         Creature.stepRight(o)
     end
     Creature.eatUp(o,Int(var_1))
     for i = 1:var_1
         local var_15::Bool = true
         local var_9::Bool = true
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         Creature.eatRight(o,Int(var_6))
         Creature.eatLeft(o,Int(var_6))
         Creature.stepRight(o)
         Creature.eatRight(o,Int(var_6))
         var_9 = var_9 & var_9
         Creature.eatRight(o,Int(var_6))
         Creature.stepRight(o)
         var_9 = var_9 | var_9
         Creature.eatRight(o,Int(var_6))
         Creature.stepLeft(o)
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatLeft(o,Int(var_6))
         Creature.stepRight(o)
         Creature.eatRight(o,Int(var_6))
         Creature.eatDown(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatUp(o,Int(var_6))
     end
     Creature.stepUp(o)
     for i = 1:var_1
         Creature.stepLeft(o)
         Creature.stepLeft(o)
     end
     for i = 1:var_1
         Creature.stepLeft(o)
         Creature.stepLeft(o)
     end
     Creature.eatDown(o,Int(var_1))
     Creature.stepLeft(o)
     for i = 1:var_1
         local var_25::Bool = true
     end
     for i = 1:var_1
         local var_29::Int64 = -5975265535146971507
         Creature.stepLeft(o)
     end
     return true
 end),[100,300,95,0,1,1,1],1,6745,1,5033218,0x01186d70,Dict(23832=>23832),[2,1],43),RpcApi.SimpleOrganism(0x00000000000476de,:(function (o::Creature.Organism,)
     function func_19(var_17::String="A06aRuJT",var_18::Int64=-4494005246402814635)
         Creature.stepLeft(o)
         return var_17
     end
     local var_1::Int8 = 121
     for i = 1:var_1
         Creature.stepRight(o)
         Creature.stepRight(o)
     end
     Creature.eatUp(o,Int(var_1))
     for i = 1:var_1
         local var_15::Bool = true
         local var_9::Bool = true
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         Creature.eatRight(o,Int(var_6))
         Creature.eatLeft(o,Int(var_6))
         Creature.stepRight(o)
         Creature.eatRight(o,Int(var_6))
         var_9 = var_9 & var_9
         Creature.eatRight(o,Int(var_6))
         Creature.stepRight(o)
         var_9 = var_9 | var_9
         Creature.eatRight(o,Int(var_6))
         Creature.stepLeft(o)
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatLeft(o,Int(var_6))
         Creature.stepRight(o)
         Creature.eatRight(o,Int(var_6))
         Creature.eatDown(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatUp(o,Int(var_6))
     end
     Creature.stepUp(o)
     for i = 1:var_1
         Creature.stepLeft(o)
         Creature.stepLeft(o)
     end
     for i = 1:var_1
         Creature.stepLeft(o)
         Creature.stepLeft(o)
     end
     Creature.eatDown(o,Int(var_1))
     Creature.stepLeft(o)
     for i = 1:var_1
         local var_25::Bool = true
     end
     for i = 1:var_1
         local var_29::Int64 = -5975265535146971507
         Creature.stepLeft(o)
     end
     return true
 end),[100,300,95,0,1,1,1],1,6745,1,4852283,0x01186d70,Dict(23832=>23832),[1,1],43),RpcApi.SimpleOrganism(0x000000000004a6df,:(function (o::Creature.Organism,)
     function func_32(var_30::Int8=-88,var_31::Bool=true)
         local var_34::String = "HmoFFpog"
         return var_30
     end
     function func_19(var_17::String="A06aRuJT",var_18::Int64=-4494005246402814635)
         local var_36::String = "jdoTH62s"
         return var_17
     end
     local var_1::Int8 = 121
     for i = 1:var_1
         Creature.stepRight(o)
         Creature.stepRight(o)
     end
     Creature.stepDown(o)
     Creature.eatUp(o,Int(var_1))
     for i = 1:var_1
         local var_15::Bool = true
         local var_9::Bool = true
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatLeft(o,Int(var_6))
         Creature.stepRight(o)
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.stepRight(o)
         var_9 = var_9 | var_9
         Creature.eatRight(o,Int(var_6))
         Creature.stepLeft(o)
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatUp(o,Int(var_6))
         Creature.stepRight(o)
         Creature.stepRight(o)
         Creature.eatRight(o,Int(var_6))
         Creature.eatDown(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatUp(o,Int(var_6))
     end
     Creature.stepUp(o)
     for i = 1:var_1
         Creature.stepLeft(o)
     end
     for i = 1:var_1
         local var_35::Bool = true
         Creature.stepLeft(o)
     end
     for i = 1:var_1
         Creature.stepLeft(o)
         Creature.stepLeft(o)
         Creature.stepLeft(o)
     end
     Creature.eatDown(o,Int(var_1))
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     for i = 1:var_1
         local var_25::Bool = true
     end
     for i = 1:var_1
         local var_29::Int64 = -5975265535146971507
         Creature.stepLeft(o)
     end
     return true
 end),[100,300,95,0,1,1,1],1,6745,1,4130906,0x011a2250,Dict(23832=>23832),[412,333],51),RpcApi.SimpleOrganism(0x000000000004a6a8,:(function (o::Creature.Organism,)
     function func_32(var_30::Int8=-88,var_31::Bool=true)
         for i = 1:var_30
         end
         return var_30
     end
     function func_19(var_17::String="A06aRuJT",var_18::Int64=-4494005246402814635)
         return var_17
     end
     local var_1::Int8 = 121
     Creature.stepDown(o)
     Creature.eatUp(o,Int(var_1))
     for i = 1:var_1
         local var_15::Bool = true
         local var_9::Bool = true
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatLeft(o,Int(var_6))
         Creature.stepRight(o)
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.stepRight(o)
         Creature.stepDown(o)
         Creature.eatRight(o,Int(var_6))
         Creature.eatDown(o,Int(var_6))
         Creature.stepUp(o)
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatUp(o,Int(var_6))
         Creature.stepRight(o)
         Creature.stepRight(o)
         Creature.eatRight(o,Int(var_6))
         o.mem[var_8] = var_8
         Creature.eatRight(o,Int(var_6))
         Creature.eatUp(o,Int(var_6))
     end
     Creature.stepUp(o)
     for i = 1:var_1
     end
     for i = 1:var_1
         Creature.stepLeft(o)
     end
     for i = 1:var_1
         Creature.stepLeft(o)
         Creature.stepLeft(o)
         Creature.stepLeft(o)
     end
     Creature.eatUp(o,Int(var_1))
     Creature.eatRight(o,Int(var_1))
     Creature.eatDown(o,Int(var_1))
     for i = 1:var_1
         local var_25::Bool = true
         Creature.stepRight(o)
     end
     return true
 end),[100,300,95,0,1,1,1],2,1581,1,4076260,0x01201b60,Dict(23832=>23832),[275,316],43),RpcApi.SimpleOrganism(0x000000000004a726,:(function (o::Creature.Organism,)
     function func_19(var_17::String="A06aRuJT",var_18::Int64=-4494005246402814635)
         Creature.stepLeft(o)
         return var_17
     end
     local var_1::Int8 = 121
     for i = 1:var_1
         Creature.stepRight(o)
         Creature.stepRight(o)
     end
     Creature.eatUp(o,Int(var_1))
     for i = 1:var_1
         local var_15::Bool = true
         local var_9::Bool = true
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         Creature.eatRight(o,Int(var_6))
         Creature.eatLeft(o,Int(var_6))
         Creature.stepRight(o)
         Creature.eatRight(o,Int(var_6))
         var_9 = var_9 & var_9
         Creature.eatRight(o,Int(var_6))
         Creature.stepRight(o)
         var_9 = var_9 | var_9
         Creature.eatRight(o,Int(var_6))
         var_8 = var_8 - var_8
         Creature.stepLeft(o)
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatLeft(o,Int(var_6))
         Creature.stepRight(o)
         Creature.eatRight(o,Int(var_6))
         Creature.eatDown(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatUp(o,Int(var_6))
     end
     Creature.stepUp(o)
     for i = 1:var_1
         Creature.stepLeft(o)
         Creature.stepLeft(o)
     end
     for i = 1:var_1
         Creature.stepLeft(o)
         Creature.stepLeft(o)
     end
     Creature.eatDown(o,Int(var_1))
     Creature.stepLeft(o)
     for i = 1:var_1
         local var_25::Bool = true
     end
     for i = 1:var_1
         local var_29::Int64 = -5975265535146971507
         Creature.stepLeft(o)
     end
     return true
 end),[100,300,95,0,1,1,1],1,6745,1,4039362,0x01188348,Dict(23832=>23832),[4,1],44),RpcApi.SimpleOrganism(0x000000000004a5ff,:(function (o::Creature.Organism,)
     function func_32(var_30::Int8=-88,var_31::Bool=true)
         local var_34::String = "HmoFFpog"
         return var_30
     end
     function func_19(var_17::String="A06aRuJT",var_18::Int64=-4494005246402814635)
         local var_35::Int8 = -121
         return var_17
     end
     local var_1::Int8 = 121
     for i = 1:var_1
         Creature.stepRight(o)
         Creature.stepRight(o)
     end
     Creature.stepDown(o)
     Creature.eatUp(o,Int(var_1))
     for i = 1:var_1
         local var_15::Bool = true
         local var_9::Bool = true
         local var_8::Int16 = 23832
         local var_6::Int8 = 66
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatLeft(o,Int(var_6))
         Creature.stepRight(o)
         Creature.eatRight(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.stepRight(o)
         var_9 = var_9 | var_9
         Creature.eatRight(o,Int(var_6))
         Creature.stepLeft(o)
         var_8 = if haskey(o.mem,var_8)
                 o.mem[var_8]
             else
                 var_8
             end
         Creature.eatUp(o,Int(var_6))
         Creature.stepRight(o)
         Creature.stepRight(o)
         Creature.eatRight(o,Int(var_6))
         Creature.eatDown(o,Int(var_6))
         Creature.eatRight(o,Int(var_6))
         Creature.eatUp(o,Int(var_6))
     end
     Creature.stepUp(o)
     Creature.eatDown(o,Int(var_1))
     for i = 1:var_1
         Creature.stepLeft(o)
     end
     for i = 1:var_1
         Creature.stepLeft(o)
     end
     for i = 1:var_1
         Creature.stepLeft(o)
         Creature.stepLeft(o)
         Creature.stepLeft(o)
     end
     Creature.eatDown(o,Int(var_1))
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     Creature.stepLeft(o)
     for i = 1:var_1
         local var_25::Bool = true
     end
     for i = 1:var_1
         local var_29::Int64 = -5975265535146971507
         Creature.stepLeft(o)
     end
     return true
 end),[100,300,95,0,1,1,1],1,6745,1,3829111,0x011a2250,Dict(23832=>23832),[422,333],51)]

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
@testset "GPS L1" begin

    @test @inferred(get_center_frequency(GPSL1)) == 1.57542e9Hz
    @test @inferred(get_code_length(GPSL1)) == 1023
    @test @inferred(get_secondary_code_length(GPSL1)) == 1
    @test @inferred(get_code(GPSL1, 0, 1)) == 1
    @test @inferred(get_code(GPSL1, 0.0, 1)) == 1
    @test @inferred(get_code_unsafe(GPSL1, 0.0, 1)) == 1
    @test @inferred(get_data_frequency(GPSL1)) == 50Hz
    @test @inferred(get_code_frequency(GPSL1)) == 1023e3Hz
    @test get_code.(GPSL1, 0:1022, 1) == L1_SAT1_CODE

end


text:
using Oceananigans
using CairoMakie
using LazyGrids
# using UnicodePlots

grid = RegularRectilinearGrid(size=(128, 128), x=(-5, 5), y=(-5, 5),
topology=(Periodic, Periodic, Flat))
model = NonhydrostaticModel(grid=grid, tracers=:c, buoyancy=nothing)

initial_c(x, y, z) = exp.(-x.^2 - y.^2)
set!(model, u=1, c=initial_c)

xs =  LinRange(-5, 5, 128)
ys =  LinRange(-5, 5, 128)

global i = 400

while i in 700

    simulation = Simulation(model, Δt=1e-2, stop_iteration = i)
    run!(simulation)

    fig = Figure(resolution=(700, 450), fontsize=18, font="sans")
    ax = fig[1, 1] = Axis(fig, xlabel="x", ylabel="y")
    cMatrixTrimmed = interior(model.tracers.c)[:, :, 1]'
    CairoMakie.heatmap!(xs, ys, cMatrixTrimmed, colormap = :deep)

    fig

    i += 100

end


text:
using NBodySimulator, Dates

const T = 370 # °K
const T0 = 275 # °K
const kb = 8.3144598e-3 # kJ/(K*mol)
const ϵOO = 0.1554253*4.184 # kJ
const σOO = 0.3165492 # nm
const ρ = 997/1.6747# Da/nm^3
const mO = 15.999 # Da
const mH = 1.00794 # Da
const mH2O = mO+2*mH
const N = 216#floor(Int, ρ * L^3 / m)
const L = (mH2O*N/ρ)^(1/3)
const R = 0.9 # ~3*σOO
const Rel = 0.49*L
const v_dev = sqrt(kb * T /mH2O)
const τ = 0.5e-3 # ps
const t1 = 0τ
const t2 = 5τ # ps
const k_bond = 1059.162*4.184*1e2 # kJ/(mol*nm^2)
const k_angle = 75.90*4.184 # kJ/(mol*rad^2)
const rOH = 0.1012 # nm
const ∠HOH = 113.24*pi/180 # rad
const qH = 0.41
const qO = -0.82
const k = 138.935458 #
bodies = generate_bodies_in_cell_nodes(N, mH2O, v_dev, L)
jl_parameters = LennardJonesParameters(ϵOO, σOO, R)
e_parameters = ElectrostaticParameters(k, Rel)
spc_paramters = SPCFwParameters(rOH, ∠HOH, k_bond, k_angle)
pbc = CubicPeriodicBoundaryConditions(L)
water = WaterSPCFw(bodies, mH, mO, qH, qO,  jl_parameters, e_parameters, spc_paramters);
simulation = NBodySimulation(water, (t1, t2), pbc, kb);
result = @time run_simulation(simulation, VelocityVerlet(), dt=τ)

#=
time_now = Dates.format(now(), "yyyy_mm_dd_HH_MM_SS")
Nactual = length(bodies)
timesteps = round(length(result.solution.t))

ts = t1:τ:result.solution.t[end-1]

for t in ts
    coordinates = get_position(result, t)
end
=#


text:
# This file is imported by most Augmentor.Operations
# documentation pages. The purpose of this file is to generate
# one or more example images for what the operation does.
# Additionally the result is returned as a markdown table
# showing the input image and the corresponding output image.
#
# Example use for single output image as png:
#
# ```@eval
# include("optable.jl")
# @optable Rotate(15)
# ```
#
# Example use for 8 output images compiled into a gif:
#
# ```@eval
# include("optable.jl")
# @optable 8 => Rotate(-15:15)
# ```

using Augmentor, Images, Colors
using Reel, PaddedViews, OffsetArrays
Reel.set_output_type("gif")

srand(1337)

if !isfile("../assets/testpattern.png")
    pattern = imresize(testpattern(), (240, 320))
    save("../assets/testpattern.png", pattern)
end

pattern = load("../assets/testpattern.png")
pattern_noalpha = ((1 .- alpha.(pattern)) .* colorant"#F3F6F6") .+ (alpha.(pattern) .* color.(pattern))

function drawborder!(img, col)
    img[1:end,   1] .= fill(col, size(img,1))
    img[1:end, end] .= fill(col, size(img,1))
    img[1,   1:end] .= fill(col, size(img,2))
    img[end, 1:end] .= fill(col, size(img,2))
    img
end

centered(img) = OffsetArray(img, convert(Tuple, 1 .- round.(Int, ImageTransformations.center(img))))

macro optable(expr)
    if expr.args[1] == :(=>) && expr.args[2] isa Int
        n = expr.args[2]
        nexpr = expr.args[3]
        name = string(nexpr.args[1])
        descr = string(nexpr)
        :(optable($(esc(nexpr)), $name, $descr, $n))
    elseif expr.args[1] == :(=>) && expr.args[2] isa String
        name = expr.args[2]
        nexpr = expr.args[3]
        descr = string(nexpr)
        :(optable($(esc(nexpr)), $name, $descr))
    else
        name = string(expr.args[1])
        descr = string(expr)
        :(optable($(esc(expr)), $name, $descr))
    end
end

function optable(op, name, descr)
    fname = joinpath("..", "assets", string(name, ".png"))
    i = 2
    while isfile(fname)
        fname = joinpath("..", "assets", string(name, i, ".png"))
        i = i + 1
    end
    out = augment(pattern, op)
    save(fname, out)
    header = length(descr) < 20 ? "Output for `$descr`" : "`$descr`"
    tbl = string(
        "Input | $header\n",
        "------|--------\n",
        "![input](../assets/testpattern.png) | ![output]($fname)\n"
    )
    Markdown.parse(tbl)
end

function optable(op, name, descr, n)
    fname = joinpath("..", "assets", string(name, ".gif"))
    i = 2
    while isfile(fname)
        fname = joinpath("..", "assets", string(name, i, ".gif"))
        i = i + 1
    end
    raw_imgs = [centered(drawborder!(augment(pattern_noalpha, op), colorant"pink")) for i in 1:n]
    imgs = map(parent, map(copy, [paddedviews(colorant"#F3F6F6", raw_imgs...)...]))
    insert!(imgs, 1, first(imgs)) # otherwise loop isn't smooth
    film = roll(imgs, fps = 2)
    write(fname, film)
    header = length(descr) < 20 ? "Samples for `$descr`" : "`$descr`"
    tbl = string(
        "Input | $header\n",
        "------|--------\n",
        "![input](../assets/testpattern.png) | ![output]($fname)\n"
    )
    Markdown.parse(tbl)
end


text:
#!/usr/bin/env julia
###### AC-OPF using Optim ######
#
# implementation reference: https://julianlsolvers.github.io/Optim.jl/stable/#examples/generated/ipnewton_basics/
# currently does not converge to a feasible point, root cause in unclear
# `debug/optim-debug.jl` can be used to confirm it will converge if given a suitable starting point
#

import PowerModels
import Optim

function solve_opf(file_name)
    time_data_start = time()

    data = PowerModels.parse_file(file_name)
    PowerModels.standardize_cost_terms!(data, order=2)
    PowerModels.calc_thermal_limits!(data)
    ref = PowerModels.build_ref(data)[:it][:pm][:nw][0]

    data_load_time = time() - time_data_start

    time_model_start = time()

    bus_pd = Dict(i => 0.0 for (i,bus) in ref[:bus])
    bus_qd = Dict(i => 0.0 for (i,bus) in ref[:bus])

    bus_gs = Dict(i => 0.0 for (i,bus) in ref[:bus])
    bus_bs = Dict(i => 0.0 for (i,bus) in ref[:bus])

    for (i,bus) in ref[:bus]
        if length(ref[:bus_loads][i]) > 0
            bus_pd[i] = sum(ref[:load][l]["pd"] for l in ref[:bus_loads][i])
            bus_qd[i] = sum(ref[:load][l]["qd"] for l in ref[:bus_loads][i])
        end

        if length(ref[:bus_shunts][i]) > 0
            bus_gs[i] = sum(ref[:shunt][s]["gs"] for s in ref[:bus_shunts][i])
            bus_bs[i] = sum(ref[:shunt][s]["bs"] for s in ref[:bus_shunts][i])
        end
    end

    br_g = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_b = Dict(i => 0.0 for (i,branch) in ref[:branch])

    br_tr = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_ti = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_tm = Dict(i => 0.0 for (i,branch) in ref[:branch])

    br_g_fr = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_b_fr = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_g_to = Dict(i => 0.0 for (i,branch) in ref[:branch])
    br_b_to = Dict(i => 0.0 for (i,branch) in ref[:branch])

    for (i,branch) in ref[:branch]
        g, b = PowerModels.calc_branch_y(branch)
        tr, ti = PowerModels.calc_branch_t(branch)

        br_g[i] = g
        br_b[i] = b

        br_tr[i] = tr
        br_ti[i] = ti
        br_tm[i] = tr^2 + ti^2

        br_g_fr[i] = branch["g_fr"]
        br_b_fr[i] = branch["b_fr"]
        br_g_to[i] = branch["g_to"]
        br_b_to[i] = branch["b_to"]
    end

    var_lookup = Dict{String,Int}()

    var_init = Float64[]
    var_lb = Float64[]
    var_ub = Float64[]

    var_idx = 1
    for (i,bus) in ref[:bus]
        push!(var_init, 0.0) #va
        push!(var_lb, -Inf)
        push!(var_ub, Inf)
        var_lookup["va_$(i)"] = var_idx
        var_idx += 1

        push!(var_init, 1.0) #vm
        push!(var_lb, bus["vmin"])
        push!(var_ub, bus["vmax"])
        var_lookup["vm_$(i)"] = var_idx
        var_idx += 1
    end

    for (i,gen) in ref[:gen]
        #push!(var_init, 0.0) #pg
        push!(var_init, (gen["pmax"]+gen["pmin"])/2) # non-standard start
        push!(var_lb, gen["pmin"])
        push!(var_ub, gen["pmax"])
        var_lookup["pg_$(i)"] = var_idx
        var_idx += 1

        #push!(var_init, 0.0) #qg
        push!(var_init, (gen["qmax"]+gen["qmin"])/2) # non-standard start
        push!(var_lb, gen["qmin"])
        push!(var_ub, gen["qmax"])
        var_lookup["qg_$(i)"] = var_idx
        var_idx += 1
    end

    for (l,i,j) in ref[:arcs]
        branch = ref[:branch][l]

        push!(var_init, 0.0) #p
        push!(var_lb, -branch["rate_a"])
        push!(var_ub,  branch["rate_a"])
        var_lookup["p_$(l)_$(i)_$(j)"] = var_idx
        var_idx += 1

        push!(var_init, 0.0) #q
        push!(var_lb, -branch["rate_a"])
        push!(var_ub,  branch["rate_a"])
        var_lookup["q_$(l)_$(i)_$(j)"] = var_idx
        var_idx += 1
    end

    @assert var_idx == length(var_init)+1
    #total_callback_time = 0.0
    function opf_objective(x)
        #start = time()
        cost = 0.0
        for (i,gen) in ref[:gen]
            pg = x[var_lookup["pg_$(i)"]]
            cost += gen["cost"][1]*pg^2 + gen["cost"][2]*pg + gen["cost"][3]
        end
        #total_callback_time += time() - start
        return cost
    end

    function opf_constraints(c,x)
        #start = time()
        va = Dict(i => x[var_lookup["va_$(i)"]] for (i,bus) in ref[:bus])
        vm = Dict(i => x[var_lookup["vm_$(i)"]] for (i,bus) in ref[:bus])

        pg = Dict(i => x[var_lookup["pg_$(i)"]] for (i,gen) in ref[:gen])
        qg = Dict(i => x[var_lookup["qg_$(i)"]] for (i,gen) in ref[:gen])

        p = Dict((l,i,j) => x[var_lookup["p_$(l)_$(i)_$(j)"]] for (l,i,j) in ref[:arcs])
        q = Dict((l,i,j) => x[var_lookup["q_$(l)_$(i)_$(j)"]] for (l,i,j) in ref[:arcs])

        vm_fr = Dict(l => vm[branch["f_bus"]] for (l,branch) in ref[:branch])
        vm_to = Dict(l => vm[branch["t_bus"]] for (l,branch) in ref[:branch])
        va_fr = Dict(l => va[branch["f_bus"]] for (l,branch) in ref[:branch])
        va_to = Dict(l => va[branch["t_bus"]] for (l,branch) in ref[:branch])

        va_con = [va[i] for (i,bus) in ref[:ref_buses]]

        #     @constraint(model,
        #         sum(p[a] for a in ref[:bus_arcs][i]) ==
        #         sum(pg[g] for g in ref[:bus_gens][i]) -
        #         sum(load["pd"] for load in bus_loads) -
        #         sum(shunt["gs"] for shunt in bus_shunts)*vm[i]^2
        #     )
        power_balance_p_con = [
           sum(pg[j] for j in ref[:bus_gens][i]; init=0.0) -
           bus_pd[i] -
           bus_gs[i]*vm[i]^2 -
           sum(p[a] for a in ref[:bus_arcs][i])
           for (i,bus) in ref[:bus]
        ]

        #     @constraint(model,
        #         sum(q[a] for a in ref[:bus_arcs][i]) ==
        #         sum(qg[g] for g in ref[:bus_gens][i]) -
        #         sum(load["qd"] for load in bus_loads) +
        #         sum(shunt["bs"] for shunt in bus_shunts)*vm[i]^2
        #     )
        power_balance_q_con = [
           sum(qg[j] for j in ref[:bus_gens][i]; init=0.0) -
           bus_qd[i] +
           bus_bs[i]*vm[i]^2 -
           sum(q[a] for a in ref[:bus_arcs][i])
           for (i,bus) in ref[:bus]
        ]

        # @NLconstraint(model, p_fr ==  (g+g_fr)/tm*vm_fr^2 + (-g*tr+b*ti)/tm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/tm*(vm_fr*vm_to*sin(va_fr-va_to)) )
        power_flow_p_from_con = [
           (br_g[l]+br_g_fr[l])/br_tm[l]*vm_fr[l]^2 +
           (-br_g[l]*br_tr[l]+br_b[l]*br_ti[l])/br_tm[l]*(vm_fr[l]*vm_to[l]*cos(va_fr[l]-va_to[l])) +
           (-br_b[l]*br_tr[l]-br_g[l]*br_ti[l])/br_tm[l]*(vm_fr[l]*vm_to[l]*sin(va_fr[l]-va_to[l])) -
           p[(l,i,j)]
           for (l,i,j) in ref[:arcs_from]
        ]

        # @NLconstraint(model, p_to ==  (g+g_to)*vm_to^2 + (-g*tr-b*ti)/tm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-b*tr+g*ti)/tm*(vm_to*vm_fr*sin(va_to-va_fr)) )
        power_flow_p_to_con = [
           (br_g[l]+br_g_to[l])*vm_to[l]^2 +
           (-br_g[l]*br_tr[l]-br_b[l]*br_ti[l])/br_tm[l]*(vm_to[l]*vm_fr[l]*cos(va_to[l]-va_fr[l])) +
           (-br_b[l]*br_tr[l]+br_g[l]*br_ti[l])/br_tm[l]*(vm_to[l]*vm_fr[l]*sin(va_to[l]-va_fr[l])) -
           p[(l,i,j)]
           for (l,i,j) in ref[:arcs_to]
        ]

        # @NLconstraint(model, q_fr == -(b+b_fr)/tm*vm_fr^2 - (-b*tr-g*ti)/tm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/tm*(vm_fr*vm_to*sin(va_fr-va_to)) )
        power_flow_q_from_con = [
           -(br_b[l]+br_b_fr[l])/br_tm[l]*vm_fr[l]^2 -
           (-br_b[l]*br_tr[l]-br_g[l]*br_ti[l])/br_tm[l]*(vm_fr[l]*vm_to[l]*cos(va_fr[l]-va_to[l])) +
           (-br_g[l]*br_tr[l]+br_b[l]*br_ti[l])/br_tm[l]*(vm_fr[l]*vm_to[l]*sin(va_fr[l]-va_to[l])) -
           q[(l,i,j)]
           for (l,i,j) in ref[:arcs_from]
        ]

        # @NLconstraint(model, q_to == -(b+b_to)*vm_to^2 - (-b*tr+g*ti)/tm*(vm_to*vm_fr*cos(va_fr-va_to)) + (-g*tr-b*ti)/tm*(vm_to*vm_fr*sin(va_to-va_fr)) )
        power_flow_q_to_con = [
           -(br_b[l]+br_b_to[l])*vm_to[l]^2 -
           (-br_b[l]*br_tr[l]+br_g[l]*br_ti[l])/br_tm[l]*(vm_to[l]*vm_fr[l]*cos(va_fr[l]-va_to[l])) +
           (-br_g[l]*br_tr[l]-br_b[l]*br_ti[l])/br_tm[l]*(vm_to[l]*vm_fr[l]*sin(va_to[l]-va_fr[l])) -
           q[(l,i,j)]
           for (l,i,j) in ref[:arcs_to]
        ]

        # @constraint(model, va_fr - va_to <= branch["angmax"])
        # @constraint(model, va_fr - va_to >= branch["angmin"])
        power_flow_vad_con = [
           va_fr[l] - va_to[l]
           for (l,i,j) in ref[:arcs_from]
        ]

        # @constraint(model, p_fr^2 + q_fr^2 <= branch["rate_a"]^2)
        power_flow_mva_from_con = [
           p[(l,i,j)]^2 + q[(l,i,j)]^2
           for (l,i,j) in ref[:arcs_from]
        ]

        # @constraint(model, p_to^2 + q_to^2 <= branch["rate_a"]^2)
        power_flow_mva_to_con = [
           p[(l,i,j)]^2 + q[(l,i,j)]^2
           for (l,i,j) in ref[:arcs_to]
        ]

        c .= [
            va_con...,
            power_balance_p_con...,
            power_balance_q_con...,
            power_flow_p_from_con...,
            power_flow_p_to_con...,
            power_flow_q_from_con...,
            power_flow_q_to_con...,
            power_flow_vad_con...,
            power_flow_mva_from_con...,
            power_flow_mva_to_con...,
        ]
        #total_callback_time += time() - start
        return c
    end

    con_lbs = Float64[]
    con_ubs = Float64[]

    #@constraint(model, va[i] == 0)
    for (i,bus) in ref[:ref_buses]
        push!(con_lbs, 0.0)
        push!(con_ubs, 0.0)
    end

    #power_balance_p_con
    for (i,bus) in ref[:bus]
        push!(con_lbs, 0.0)
        push!(con_ubs, 0.0)
        #push!(con_lbs, -Inf)
        #push!(con_ubs, Inf)
    end

    #power_balance_q_con
    for (i,bus) in ref[:bus]
        push!(con_lbs, 0.0)
        push!(con_ubs, 0.0)
        #push!(con_lbs, -Inf)
        #push!(con_ubs, Inf)
    end

    #power_flow_p_from_con
    for (l,i,j) in ref[:arcs_from]
        push!(con_lbs, 0.0)
        push!(con_ubs, 0.0)
    end

    #power_flow_p_to_con
    for (l,i,j) in ref[:arcs_to]
        push!(con_lbs, 0.0)
        push!(con_ubs, 0.0)
    end

    #power_flow_q_from_con
    for (l,i,j) in ref[:arcs_from]
        push!(con_lbs, 0.0)
        push!(con_ubs, 0.0)
    end

    #power_flow_q_to_con
    for (l,i,j) in ref[:arcs_to]
        push!(con_lbs, 0.0)
        push!(con_ubs, 0.0)
    end

    #power_flow_vad_con
    for (l,i,j) in ref[:arcs_from]
        branch = ref[:branch][l]
        push!(con_lbs, branch["angmin"])
        push!(con_ubs, branch["angmax"])
    end

    #power_flow_mva_from_con
    for (l,i,j) in ref[:arcs_from]
        branch = ref[:branch][l]
        push!(con_lbs, -Inf)
        push!(con_ubs, branch["rate_a"]^2)
    end

    #power_flow_mva_to_con
    for (l,i,j) in ref[:arcs_to]
        branch = ref[:branch][l]
        push!(con_lbs, -Inf)
        push!(con_ubs, branch["rate_a"]^2)
    end

    model_variables = length(var_init)
    model_constraints = length(opf_constraints(zeros(length(con_lbs)), var_init))
    println("variables: $(model_variables), $(length(var_lb)), $(length(var_ub))")
    println("constraints: $(model_constraints), $(length(con_lbs)), $(length(con_ubs))")

    df = Optim.TwiceDifferentiable(opf_objective, var_init)
    dfc = Optim.TwiceDifferentiableConstraints(opf_constraints, var_lb, var_ub, con_lbs, con_ubs)

    model_build_time = time() - time_model_start

    time_solve_start = time()

    options = Optim.Options(show_trace=true)

    # NOTE: had to change initial guess to be an interior point, otherwise getting NaN values
    res = Optim.optimize(df, dfc, var_init, Optim.IPNewton(), options)
    #res = Optim.optimize(df, dfc, var_init, Optim.LBFGS(), options) #  StackOverflowError:
    #res = Optim.optimize(df, dfc, var_init, Optim.NelderMead(), options) #  StackOverflowError:
    display(res)

    sol = res.minimizer
    cost = res.minimum

    solve_time = time() - time_solve_start
    total_time = time() - time_data_start

    # NOTE: confirmed these constraint violations can be eliminated
    # if a better starting point is used
    sol_eval = opf_constraints(zeros(length(con_lbs)), sol)
    vio_lb = [max(v,0) for v in (con_lbs .- sol_eval)]
    vio_ub = [max(v,0) for v in (sol_eval .- con_ubs)]
    const_vio = vio_lb .+ vio_ub
    #println(const_vio)
    println("total constraint violation: $(sum(const_vio))")
    constraint_tol = 1e-6
    feasible = (sum(const_vio) <= constraint_tol)

    if !feasible
        @warn "Optim optimize failed to satify the problem constraints"
    end

    println("")
    println("\033[1mSummary\033[0m")
    println("   case........: $(file_name)")
    println("   variables...: $(model_variables)")
    println("   constraints.: $(model_constraints)")
    println("   feasible....: $(feasible)")
    println("   cost........: $(round(Int, cost))")
    println("   total time..: $(total_time)")
    println("     data time.: $(data_load_time)")
    println("     build time: $(model_build_time)")
    println("     solve time: $(solve_time)")
    # println("      callbacks: $(total_callback_time)")
    println("")

    return Dict(
        "case" => file_name,
        "variables" => model_variables,
        "constraints" => model_constraints,
        "feasible" => feasible,
        "cost" => cost,
        "time_total" => total_time,
        "time_data" => data_load_time,
        "time_build" => model_build_time,
        "time_solve" => solve_time,
        #"time_callbacks" => TBD,
    )
end

if isinteractive() == false
    solve_opf("$(@__DIR__)/data/pglib_opf_case5_pjm.m")
end

max_stars_repo_name: lanl-ansi/rosetta-opf

max_issues_repo_name: lanl-ansi/rosetta-opf

max_forks_repo_name: lanl-ansi/rosetta-opf


text:
#    Problem : GROUP A
#    *********
#    This is a simplified version of problem NONDIA.
#
#    Original SIF Source:
#    G. Li,
#    "The secant/finite difference algorithm for solving sparse
#    nonlinear systems of equations",
#    SIAM Journal on Optimization, (to appear), 1990.
#
#    Implementation translated from Source:
#    http://eprints.tsu.ge/234/14/Tests%20collection-K-F.pdf
#
#    LIARWHD.SIF classification  SUR2-AN-V-0
#
#    Number of variables (at least 2)
#
# Daniel Henderson, 08/2021

f = (x) -> begin
    fx = 0.0
    for i in firstindex(x):lastindex(x)
        α = 2.0(x[i]^2 - x[1])
        γ = x[i]-1
        fx += α^2 + γ^2
    end
    return fx
end

g! = (g, x) -> begin
    for i in firstindex(x):lastindex(x)
        α = 2.0(x[i]^2 - x[1])
        γ = x[i]-1
        g[i] = 8.0x[i]α + 2.0γ
        g[1] -= 4.0α
    end
    return g
end

fg! = (g, x) -> begin
    fx = 0.0
    for i in firstindex(x):lastindex(x)
        α = 2.0(x[i]^2 - x[1])
        γ = x[i]-1
        fx += α^2 + γ^2
        g[i] = 8.0x[i]α + 2.0γ
        g[1] -= 4.0α
    end
    return fx, g
end

init = (n::Int=5000) -> begin
	x0 = 4.0*ones(n)
    return n, x0
end

TestSet["LIARWHD"] = UncProgram("LIARWHD", f, g!, fg!, init)


text:
using Compat
import Compat.read

export PcapFileHeader, PcapRec, PcapOffline,
       pcap_get_record

mutable struct PcapFileHeader
    magic_number::UInt32
    version_major::UInt16
    version_minor::UInt16
    thiszone::Int32
    sigfigs::UInt32
    snaplen::UInt32
    network::UInt32
    PcapFileHeader() = new(0,0,0,0,0,0,0)
end # PcapFileHeader structure

mutable struct PcapRec
    ts_sec::UInt32
    ts_usec::UInt32
    incl_len::UInt32
    orig_len::UInt32
    payload::Vector{UInt8}
    PcapRec() = new(0,0,0,0, Vector{UInt8}(undef, 0))
end # PcapRec structure

mutable struct PcapOffline
    filename::AbstractString
    file::IO
    filehdr::PcapFileHeader
    record::PcapRec
    is_big::Bool
    function PcapOffline(fn::AbstractString)
        filename = fn
        file = open(fn, "r+")
        filehdr, is_big = decode_hdr(file)
        record = PcapRec()
        new(filename, file, filehdr, record, is_big)
    end # constructor
end # PcapOffline structure

#----------
# decode PCap file format header
#----------
function decode_hdr(file::Any)
    filehdr = PcapFileHeader()
    filehdr.magic_number  = read(file, UInt32)
    big_endian = false
    if filehdr.magic_number in Set([0xa1b23c4d, 0xa1b2c3d4]) # added support for nanosecond-resolution files
        big_endian = true
    end
    filehdr.version_major = big_endian ? read(file, UInt16) : ntoh(read(file, UInt16))
    filehdr.version_minor = big_endian ? read(file, UInt16) : ntoh(read(file, UInt16))
    filehdr.thiszone      = read(file, Int32)
    filehdr.sigfigs       = big_endian ? read(file, UInt32) : ntoh(read(file, UInt32))
    filehdr.snaplen       = big_endian ? read(file, UInt32) : ntoh(read(file, UInt32))
    filehdr.network       = big_endian ? read(file, UInt32) : ntoh(read(file, UInt32))
    return [filehdr, big_endian]
end # function decode_hdr

#----------
# decode next record in PCap file
#----------
function pcap_get_record(s::PcapOffline)
    rec = PcapRec()
    if (!eof(s.file))
        rec.ts_sec   = s.is_big ? read(s.file, UInt32) : ntoh(read(s.file, UInt32))
        rec.ts_usec  = s.is_big ? read(s.file, UInt32) : ntoh(read(s.file, UInt32))
        rec.incl_len = s.is_big ? read(s.file, UInt32) : ntoh(read(s.file, UInt32))
        rec.orig_len = s.is_big ? read(s.file, UInt32) : ntoh(read(s.file, UInt32))
        rec.payload  = read(s.file, rec.incl_len)
        return rec
    end
    nothing
end # function pcap_get_record


text:
function comment_using_revise(filename)
    lines = readlines(filename)
    idx = findfirst(lines .== "using Revise")
    if isnothing(idx)
        @warn "Could not find \"using Revise\" in file $filename, doing nothing."
    else
        lines[idx] = "# using Revise"
        open(filename, "w") do f
            for l in lines
                println(f, l)
            end
        end
    end
end

function process_folder(folderpath)
    testfiles = filter(x -> startswith(x, "test_"), readdir(folderpath))
    testfiles = (folderpath*"/") .* testfiles
    for f in testfiles
        comment_using_revise(f)
    end
end

foldernames = filter(x->isdir("test/"*x), readdir("test"))
folderpaths = "test/" .* foldernames

for path in folderpaths
    process_folder(path)
end


text:
function mpccmodel_test_testvectors_evalpt_basic(model_cfg::MPCCModelConfig, model::MPCCModel)

    ptevalreq = MPCCPointEvalReq()

    nt = length(model_cfg.testvectors)
    test_results = Vector{Bool}(undef, nt)

    for lp_t=1:nt

        # Get the reference point and parameters for this test
        x0 = model_cfg.testvectors[lp_t].x_val
        pr0 = model_cfg.testvectors[lp_t].pr_val
        ps0 = model_cfg.testvectors[lp_t].ps_val

        # Calculate point using model
        evalpt = mpccmodel_pointeval_basic(model, ptevalreq, x0, pr0, ps0)

        # Do the test
        test_results[lp_t] = model_cfg.testvectors[lp_t].eval_test_pt ≈ evalpt

    end

    return test_results
end


text:
import FrankWolfe
using LinearAlgebra
using Random
using Test
using SparseArrays

n = Int(1e4)
k = 3000

s = 41
Random.seed!(s)

xpi = rand(n);
total = sum(xpi);
xp = xpi # ./ total;

f(x) = norm(x - xp)^2
function grad!(storage, x)
    @. storage = 2 * (x - xp)
end

lmo = FrankWolfe.KSparseLMO(100, 1.0)

x0 = FrankWolfe.compute_extreme_point(lmo, spzeros(size(xp)...))

x, v, primal, dual_gap, _ = FrankWolfe.blended_conditional_gradient(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    L=2,
    verbose=true,
    trajectory=false,
    K=1.00,
    weight_purge_threshold=1e-10,
    epsilon=1e-9,
)

@test dual_gap ≤ 5e-4
@test f(x0) - f(x) ≥ 180

x0 = FrankWolfe.compute_extreme_point(lmo, spzeros(size(xp)...))

x, v, primal_cut, dual_gap, _ = FrankWolfe.blended_conditional_gradient(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    L=2,
    verbose=true,
    trajectory=false,
    K=1.00,
    weight_purge_threshold=1e-10,
    epsilon=1e-9,
    timeout=3.0,
)

@test primal ≤ primal_cut