text:
# GUI event loops and toolkit integration for Python, most importantly
# to support plotting in a window without blocking.  Currently, we
# support wxWidgets, Qt5, Qt4, and GTK+.

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

# global variable to specify default GUI toolkit to use
gui = :default # one of :default, :wx, :qt, :tk, or :gtk

pyexists(mod) = try
    pyimport(mod)
    true
  catch
    false
  end

# Tkinter was renamed to tkinter in Python 3
function tkinter_name()
    return pyversion.major < 3 ? "Tkinter" : "tkinter"
end

pygui_works(gui::Symbol) = gui == :default ||
    ((gui == :wx && pyexists("wx")) ||
     (gui == :gtk && pyexists("gtk")) ||
     (gui == :gtk3 && pyexists("gi")) ||
     (gui == :tk && pyexists(tkinter_name())) ||
     ((gui == :qt5 || gui == :qt_pyqt5) && pyexists("PyQt5")) ||
     (gui == :qt_pyqt4 && pyexists("PyQt4")) ||
     (gui == :qt_pyside && pyexists("PySide")) ||
     (gui == :qt4 && (pyexists("PyQt4") || pyexists("PySide"))) ||
     (gui == :qt && (pyexists("PyQt5") || pyexists("PyQt4") || pyexists("PySide"))))

# get or set the default GUI; doesn't affect running GUI
"""
    pygui()

Return the current GUI toolkit as a symbol.
"""
function pygui()
    global gui
    if gui::Symbol != :default && pygui_works(gui::Symbol)
        return gui::Symbol
    else
        for g in (:tk, :qt, :wx, :gtk, :gtk3)
            if pygui_works(g)
                gui = g
                return gui::Symbol
            end
        end
        error("No supported Python GUI toolkit is installed.")
    end
end
function pygui(g::Symbol)
    global gui
    if g != gui::Symbol
        if !(g in (:wx,:gtk,:gtk3,:tk,:qt,:qt4,:qt5,:qt_pyqt5,:qt_pyqt4,:qt_pyside,:default))
            throw(ArgumentError("invalid gui $g"))
        elseif !pygui_works(g)
            error("Python GUI toolkit for $g is not installed.")
        end
        gui = g
    end
    return g::Symbol
end

############################################################################
# Event loops for various toolkits.

# call doevent(status) every sec seconds
function install_doevent(doevent, sec::Real)
    return Base.Timer(doevent, sec, sec)
end

# For PyPlot issue #181: recent pygobject releases emit a warning
# if we don't specify which version we want:
function gtk_requireversion(gtkmodule::AbstractString, vers::VersionNumber=v"3.0")
    if startswith(gtkmodule, "gi.")
        gi = pyimport("gi")
        if gi[:get_required_version]("Gtk") === nothing
            gi[:require_version]("Gtk", string(vers))
        end
    end
end

# GTK (either GTK2 or GTK3, depending on gtkmodule):
function gtk_eventloop(gtkmodule::AbstractString, sec::Real=50e-3)
    gtk_requireversion(gtkmodule)
    gtk = pyimport(gtkmodule)
    events_pending = gtk["events_pending"]
    main_iteration = gtk["main_iteration"]
    install_doevent(sec) do async
        # handle all pending
        while PyReCall(events_pending, Bool)
            PyReCall(main_iteration, PyObject)
        end
    end
end

# As discussed in PyPlot.jl#278, Qt looks for a file qt.conf in
# the same path as the running executable, which tells it where
# to find plugins etcetera.  For Python's Qt, however, this will
# be in the path of the python executable, not julia.  Furthermore,
# we can't copy it to the location of julia (even if that location
# is writable) because that would assume that all julia programs
# use the same version of Qt (e.g. via the same Python), which
# is not necessarily the case, and even then it wouldn't work
# for other programs linking libjulia.  Unfortunately, there
# seems to be no way to change this.  However, we can at least
# use set QT_PLUGIN_PATH by parsing qt.conf ourselves, and
# this seems to fix some of the path-related problems on Windows.
# ... unfortunately, it seems fixqtpath has to be called before
# the Qt library is loaded.
function fixqtpath(qtconf=joinpath(dirname(pyprogramname),"qt.conf"))
    haskey(ENV, "QT_PLUGIN_PATH") && return false
    if isfile(qtconf)
        for line in eachline(qtconf)
            m = match(r"^\s*prefix\s*=(.*)$"i, line)
            if m !== nothing
                dir = strip(m.captures[1])
                if startswith(dir, '"') && endswith(dir, '"')
                    dir = dir[2:end-1]
                end
                plugin_path = joinpath(dir, "plugins")
                if isdir(plugin_path)
                    # for some reason I don't understand,
                    # if libpython has already been loaded, then
                    # we need to use Python's setenv rather than Julia's
                    PyDict(pyimport("os")["environ"])["QT_PLUGIN_PATH"] = realpath(plugin_path)
                    return true
                end
            end
        end
    end
    return false
end

# Qt: (PyQt5, PyQt4, or PySide module)
function qt_eventloop(QtCore::PyObject, sec::Real=50e-3)
    fixqtpath()
    instance = QtCore["QCoreApplication"]["instance"]
    AllEvents = QtCore["QEventLoop"]["AllEvents"]
    processEvents = QtCore["QCoreApplication"]["processEvents"]
    pop!(ENV, "QT_PLUGIN_PATH", "") # clean up environment
    maxtime = PyObject(50)
    install_doevent(sec) do async
        app = PyReCall(instance, PyObject)
        if app.o != pynothing[]
            app["_in_event_loop"] = true
            PyReCall(processEvents, PyObject, AllEvents, maxtime)
        end
    end
end

qt_eventloop(QtModule::AbstractString, sec::Real=50e-3) =
    qt_eventloop(pyimport("$QtModule.QtCore"), sec)

function qt_eventloop(sec::Real=50e-3)
    for QtModule in ("PyQt5", "PyQt4", "PySide")
        try
            return qt_eventloop(QtModule, sec)
        end
    end
    error("no Qt module found")
end

# wx:  (based on IPython/lib/inputhookwx.py, which is 3-clause BSD-licensed)
function wx_eventloop(sec::Real=50e-3)
    wx = pyimport("wx")
    GetApp = wx["GetApp"]
    EventLoop = wx["EventLoop"]
    EventLoopActivator = wx["EventLoopActivator"]
    install_doevent(sec) do async
        app = PyReCall(GetApp, PyObject)
        if app.o != pynothing[]
            app["_in_event_loop"] = true
            evtloop = PyReCall(EventLoop, PyObject)
            ea = PyReCall(EventLoopActivator, PyObject, evtloop)
            Pending = evtloop["Pending"]
            Dispatch = evtloop["Dispatch"]
            while PyReCall(Pending, Bool)
                PyReCall(Dispatch, PyObject)
            end
            pydecref(ea) # deactivate event loop
            PyReCall(app["ProcessIdle"], PyObject)
        end
    end
end

# Tk: (Tkinter/tkinter module)
function Tk_eventloop(sec::Real=50e-3)
    Tk = pyimport(tkinter_name())
    install_doevent(sec) do async
        root = Tk["_default_root"]
        if root.o != pynothing[]
            PyReCall(root["update"], PyObject)
        end
    end
end
# cache running event loops (so that we don't start any more than once)
const eventloops = Dict{Symbol,Timer}()

"""
    pygui_start(gui::Symbol = pygui())

Start the event loop of a certain toolkit.

The argument `gui` defaults to the current default GUI, but it could be `:wx`, `:gtk`, `:gtk3`, `:tk`, or `:qt`.

"""
function pygui_start(gui::Symbol=pygui(), sec::Real=50e-3)
    pygui(gui)
    if !haskey(eventloops, gui)
        if gui == :wx
            eventloops[gui] = wx_eventloop(sec)
        elseif gui == :gtk
            eventloops[gui] = gtk_eventloop("gtk", sec)
        elseif gui == :gtk3
            eventloops[gui] = gtk_eventloop("gi.repository.Gtk", sec)
        elseif gui == :tk
            eventloops[gui] = Tk_eventloop(sec)
        elseif gui == :qt_pyqt5 || gui == :qt5
            eventloops[gui] = qt_eventloop("PyQt5", sec)
        elseif gui == :qt_pyqt4
            eventloops[gui] = qt_eventloop("PyQt4", sec)
        elseif gui == :qt_pyside
            eventloops[gui] = qt_eventloop("PySide", sec)
        elseif gui == :qt4
            try
                eventloops[gui] = qt_eventloop("PyQt4", sec)
            catch
                eventloops[gui] = qt_eventloop("PySide", sec)
            end
        elseif gui == :qt
            eventloops[gui] = qt_eventloop(sec)
        else
            throw(ArgumentError("unsupported GUI type $gui"))
        end
    end
    gui
end

"""
    pygui_stop(gui::Symbol = pygui())

Stop any running event loop for gui. The `gui` argument defaults to current default GUI.

"""
function pygui_stop(gui::Symbol=pygui())
    if haskey(eventloops, gui)
        Base.close(pop!(eventloops, gui))
        true
    else
        false
    end
end

pygui_stop_all() = for gui in keys(eventloops); pygui_stop(gui); end

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


text:
d = SparseCat([1,2], [0.5, 0.5])
@test sprint(showdistribution, d) == "     SparseCat{Array{Int64,1},Array{Float64,1}} distribution\n     ┌                                        ┐ \n   1 ┤■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.5   \n   2 ┤■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.5   \n     └                                        ┘ "

d = SparseCat(1:50, fill(1/50, 50))
iob = IOBuffer()
io = IOContext(iob, :limit=>true, :displaysize=>(10, 7))
showdistribution(io, d)
@test String(take!(iob)) == "                     SparseCat{UnitRange{Int64},Array{Float64,1}} distribution\n                     ┌                                        ┐ \n                   1 ┤■ 0.02                                    \n                   2 ┤■ 0.02                                    \n                   3 ┤■ 0.02                                    \n   <everything else> ┤■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 0.94   \n                     └                                        ┘ "

# test that it doesn't print <everything else> when there are enough lines
d = SparseCat([:a], 1.0)
iob = IOBuffer()
io = IOContext(iob, :limit=>true, :displaysize=>(10, 7))
showdistribution(io, d)
@test String(take!(iob)) == "      SparseCat{Array{Symbol,1},Float64} distribution\n      ┌                                        ┐ \n   :a ┤■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■ 1.0   \n      └                                        ┘ "

ext: jl


text:
using BinDeps
using Compat

@BinDeps.setup

opencv_version = "3.1.0"

ignore_paths = split(strip(get(ENV, "OPENCVJL_LIBRARY_IGNORE_PATH", "")), ':')

# default validate function
default_validate = function(libpath, handle)
    for path in ignore_paths
        isempty(path) && continue
        ismatch(Regex("^$(path)"), libpath) && return false
    end
    return true
end

function cv2_library_dependency(s, validate=default_validate)
    library_dependency(s, validate=validate)
end

opencv_core = cv2_library_dependency("libopencv_core")

### Source provider ###
github_root = "https://github.com/Itseez/opencv"
provides(Sources,
         URI("$(github_root)/archive/$(opencv_version).tar.gz"),
         [opencv_core],
         unpacked_dir="opencv-$(opencv_version)")

prefix = joinpath(BinDeps.depsdir(opencv_core), "usr")
srcdir = joinpath(BinDeps.depsdir(opencv_core), "src", "opencv-$(opencv_version)")

cmake_options = [
    "-DCMAKE_INSTALL_PREFIX=$prefix",
    "-DBUILD_SHARED_LIBS=ON",
    "-DBUILD_TIFF=ON",
    "-DWITH_CUDA=OFF",
    "-DENABLE_AVX=ON",
    "-DWITH_OPENGL=ON",
    "-DWITH_OPENCL=OFF", # TODO
    "-DWITH_IPP=ON",
    "-DWITH_TBB=ON",
    "-DWITH_EIGEN=ON",
    "-DWITH_V4L=ON",
    "-DBUILD_TESTS=OFF",
    "-DBUILD_PERF_TESTS=OFF",
    "-DBUILD_EXAMPLES=OFF",
    "-DCMAKE_BUILD_TYPE=RELEASE",
    "-DBUILD_opencv_java=OFF",
    "-DBUILD_opencv_python=OFF",
    "-DBUILD_opencv_nonfree=OFF",
    "-DBUILD_opencv_viz=OFF",
    "-DBUILD_DOCS=OFF",
    "-DWITH_VTK=OFF",
    "-DWITH_MATLAB=OFF",
    "-DWITH_FFMPEG=OFF",
    "-DENABLE_PRECOMPILED_HEADERS=OFF",
#    "-DWITH_OPENMP=ON",
#     "-DWITH_QT=ON",
]

### Build opencv from source ###
provides(SimpleBuild,
          (@build_steps begin
              GetSources(opencv_core)
              @build_steps begin
                  ChangeDirectory(srcdir)
                  `mkdir -p build`
                  @build_steps begin
                      ChangeDirectory(joinpath(srcdir, "build"))
                      `rm -f CMakeCache.txt`
                      `cmake $cmake_options ..`
                      `make -j4`
                      `make install`
                  end
                end
          end), opencv_core, os = :Unix)

@BinDeps.install Dict(:libopencv_core => :libopencv_core)


text:
# initalize an array
function Init!(X,L)
    size=length(X)
    h=L/size
    for i=0:size-1
        if i>floor(size//8) && i<floor(size//2)+floor(size//8)
            X[i+1]=1.0-2*(i-floor(size/8))*h/L;
        else
            X[i+1]=0.0
        end
    end
end

# proc[1-5]! functions make the same computation, with different programing
# styles.

function proc1!(A::Array{Float64,1},B::Array{Float64,1},
                C::Array{Float64,1},D::Array{Float64,1},niter::Int64)
    for it=1:niter
        @. A=1.7*B-0.8*C-0.9*D
        A,D = D,A
    end
end
function proc2!(A::Array{Float64,1},B::Array{Float64,1},
                C::Array{Float64,1},D::Array{Float64,1},niter::Int64)
    for it=1:niter
        @. A[:]=1.7*B[:]-0.8*C[:]-0.9*D[:]
        A,D = D,A
    end
end
function proc3!(A::Array{Float64,1},B::Array{Float64,1},
                C::Array{Float64,1},D::Array{Float64,1},niter::Int64)
    size=length(A)
    for it=1:niter
        @simd for i in 1:size
            A[i]=1.7*B[i]-0.8*C[i]-0.9*D[i]
        end
        A,D = D,A
    end
end
function proc4!(A::Array{Float64,1},B::Array{Float64,1},
                C::Array{Float64,1},D::Array{Float64,1},niter::Int64)
    for it=1:niter
        A[:]=1.7*B[:]
        A[:]-=0.8*C[:]
        A[:]-=0.9*D[:]
        A,D = D,A
    end
end
function proc5!(A::Array{Float64,1},B::Array{Float64,1},
                C::Array{Float64,1},D::Array{Float64,1},niter::Int64)
    for it=1:niter
        @. A=1.7*B
        @. A-=0.8*C
        @. A-=0.9*D
        A,D = D,A
    end

end

# measure as precisely as possible the execution time for a given
# p= proc[1-5]!
function test(p,A::Array{Float64,1},B::Array{Float64,1},
              C::Array{Float64,1},D::Array{Float64,1},nit::Int)

    niter=nit
    # be sure to run once before actually running the benchmark!
    Init!(A,1.0)
    Init!(B,1.0)
    Init!(C,1.0)
    Init!(D,1.0)
    p(A,B,C,D,niter)
    #
    T=0.0
    while true
        Init!(A,1.0)
        Init!(B,1.0)
        Init!(C,1.0)
        Init!(D,1.0)
        t1 = time_ns()
        p(A,B,C,D,niter)
        t = (time_ns() -t1)/niter
        if abs(t-T)/t<0.05
            break
    else
            T=t
            if niter> 1500
                error("too many iterations",niter)
            else
                niter*=2
            end
    end

end
    T,niter
end
DD=Dict("proc1!"=>"Naïve vectorisation               ",
        "proc2!"=>"Vectorisation with braces         ",
        "proc3!"=>"Unrolled loop                     ",
        "proc4!"=>"Vectorisation with braces by steps",
        "proc5!"=>"Vectorisation with @. by steps    "
        )
# computation starts here:
size=32
sizemax=10^6
const niter=2
fw=open("RunningOn"*gethostname()*"_cl","w")
while size<sizemax
    println("size: ",size,":\n")
    A=Array{Float64}(undef,size)
    B=Array{Float64}(undef,size)
    C=Array{Float64}(undef,size)
    D=Array{Float64}(undef,size)
    tbest=10.0^29
    best=0
    t=0.0
    for p= [proc1!,proc2!,proc3!,proc4!,proc5!]
        t,it=test(p,A,B,C,D,niter)
        if t<tbest
            tbest=t
            best=p
        end
        t*=10.0^(-9)
        println(DD[string(p)]," : t= ",t," seconds")
    end
    nflops=size*5
    flops=nflops/tbest
    println("\nbest: ",DD[string(best)])

    write(fw,string(size)," ",string(tbest*10.0^(-9)),"\n")

    println("nb. flops: ",nflops, ", Gflops/s (best): ",flops)
    println("-------")

    global size*=2
    println()

end
close(fw)

max_stars_repo_name: fluiddyn/BenchmarksPythonJuliaAndCo

max_issues_repo_name: fluiddyn/BenchmarksPythonJuliaAndCo

max_forks_repo_name: fluiddyn/BenchmarksPythonJuliaAndCo


text:
using NLopt
using Optim
using ADCME
using Statistics
using LinearAlgebra
using PyCall
using SparseArrays
using Test
using Random

if has_gpu()
    use_gpu(0)
end
doctor()
sess = Session()

include("layers.jl")
include("sparse.jl")
include("random.jl")
include("io.jl")
include("variable.jl")
include("ops.jl")
include("core.jl")
include("extra.jl")
include("newton.jl")
include("optim.jl")
include("optimizers.jl")
include("ot.jl")
include("ode.jl")
include("flow.jl")


text:
using SFML

music = Music(joinpath(dirname(@__FILE__),"music.ogg"))

set_loop(music, true)
play(music)
playing = true

println("Music is $(as_seconds(get_duration(music))) seconds long")

println("Press 'p' to pause and resume")
println("Press escape to quit")

while true
    sleep(0.1)

    if is_key_pressed(KeyCode.ESCAPE)
        stop(music)
        exit(0)
    end

    if is_key_pressed(KeyCode.P)
        if !playing
            play(music)
        else
            pause(music)
        end
        playing = !playing
    end
end


text:
let
	# Test Optim.nelder_mead for all functions except Large Polynomials in Optim.UnconstrainedProblems.examples
	for (name, prob) in Optim.UnconstrainedProblems.examples
		f_prob = prob.f
		res = Optim.optimize(f_prob, prob.initial_x, NelderMead(), OptimizationOptions(iterations = 10000))
		if name == "Powell"
			res = Optim.optimize(f_prob, prob.initial_x, method=NelderMead(), g_tol = 1e-12)
		elseif name == "Large Polynomial"
			res = Optim.optimize(f_prob, prob.initial_x, method=NelderMead(initial_simplex = Optim.AffineSimplexer(1.,1.)), iterations = 500_000)
		end
		@assert norm(res.minimum - prob.solutions) < 1e-2
	end
end

# Test that deprecated syntax runs
let
	dep_prob = Optim.UnconstrainedProblems.examples["Rosenbrock"]
	optimize(dep_prob.f, dep_prob.initial_x, NelderMead(a = 1.0))
	optimize(dep_prob.f, dep_prob.initial_x, NelderMead(g = 2.0))
	optimize(dep_prob.f, dep_prob.initial_x, NelderMead(b = 0.5))
	optimize(dep_prob.f, dep_prob.initial_x, NelderMead(initial_simplex = Optim.AffineSimplexer()))
	optimize(dep_prob.f, dep_prob.initial_x, NelderMead(parameters = Optim.AdaptiveParameters()))
end


text:
using CausalityTools
using RecipesBase

@testset "Recipes for subpackages" begin

    # Test plot recipes by calling RecipesBase.apply_recipe with empty dict.
    # It should return a vector of RecipesBase.RecipeData
    d = Dict{Symbol,Any}()

    # Create an example embedding
    E = StateSpaceReconstruction.embed([diff(rand(30)) for i = 1:3], [1, 2, 3], [1, 0, -1])
    E_invariant = invariantize(E)
    invm = RectangularInvariantMeasure(E, [4, 5, 6])

    #@testset "StateSpaceReconstruction.AbstractEmbedding" begin
    #    @test typeof(RecipesBase.apply_recipe(d, E)) == Array{RecipesBase.RecipeData,1}
    #end

    @testset "PerronFrobenius.RectangularPartitionTransferOperator" begin
        @test typeof(RecipesBase.apply_recipe(d, invm.transfermatrix)) == Array{RecipesBase.RecipeData,1}
    end

    @testset "Invariant distribution" begin
        @test typeof(RecipesBase.apply_recipe(d, invm.measure)) == Array{RecipesBase.RecipeData,1}
    end

    @testset "RectangularInvariantMeasure" begin
        boxfillfactor = 3
        # Third argument is linesegment, which can be true or false. It determines
        # whether box edges around the partition elements are drawn.
        @test typeof(RecipesBase.apply_recipe(d, (invm, 3, true))) == Array{RecipesBase.RecipeData,1}
        @test typeof(RecipesBase.apply_recipe(d, (invm, 3, false))) == Array{RecipesBase.RecipeData,1}
    end

	@testset "Vizualizing triangulation and simplices" begin
		pts = rand(3, 15)
		n_pts = size(pts, 2)

		# Embed the point
		E = StateSpaceReconstruction.embed(pts)

		# Make sure last point is inside the convex hull of the previous points.
		# If it is not, move it towards the center of the embedding until it is.
		E_linearlyinvariant = invariantize(E)

		# Pick out points and their images (so that they can
		# be indexed with the same simplex index)
		originalpts = E.points[:, 1:end-1]
		forwardpts  = E.points[:, 2:end]

		# Re-create the embedding, but exclude the last point.
		E = StateSpaceReconstruction.embed(originalpts)

		# Triangulate all points but the last point.
		DT = delaunay(E);

		@test typeof(RecipesBase.apply_recipe(d, (pts, E, DT))) == Array{RecipesBase.RecipeData,1}
	end

end


text:
using JuMP

# Creates the equality constraints matrix and vector for quadprog.
# Constrained vertices include - all persistence features (their
# paired minima and maxima), the global minimum and the domain
# edges. The reconstructed function is guaranteed to preserve their
# location and their value - no new minima or maxima are created.
function eq_constraints(p::Persistence)
    n = length(p.data)

    idx = sort([p.minima;p.maxima;p.gminindex])
    sum(idx.==1) != 1 && (idx = [1;idx])
    sum(idx.==n) != 1 && push!(idx, n)

    m = length(idx)

    Aeq = sparse(1:m, idx, ones(m), m, n)
    beq = p.data[idx]

    Aeq,beq
end

# Builds a monotonicity graph based on the indices of the persistent
# extrema features which are used to reconstruct the data.
#
# Assumes that the indices of the local minima and maxima are matched.
#
# For 1D data, a monotonicity matrix is a sparse matrix, with values
# != 0 only on two diagonals, each one 1-removed from the main
# diagonal.  It is used by in_eq_constraints to create inequality
# constraints matrix and vector for quadprog.
function monotonicity_graph(p::Persistence)
    minima = sort([p.minima;p.gminindex])
    maxima = sort(p.maxima)

    first_local_min = if sum(maxima.==1) == 1
        false
    elseif sum(minima.==1) == 1
        true
    elseif minima[1] < maxima[1]
        maxima = [1;maxima]
        false
    else
        minima = [1;minima]
        true
    end

    N = length(p.data)
    if sum(maxima.==N) != 1 && sum(minima.==N) != 1
        if minima[end] > maxima[end]
            maxima = [maxima;N]
        else
            minima = [minima;N]
        end
    end

    pgrad = zeros(N);
    ngrad = zeros(N);

    all_ind = [minima; maxima];
    all_ind = sort(all_ind);

    if first_local_min
        for i=1:2:length(all_ind)-2
            pgrad[all_ind[i]:all_ind[i+1]-1] = 1;
            ngrad[all_ind[i+1]:all_ind[i+2]-1] = 1;
        end
    else
        for i=1:2:length(all_ind)-2
            ngrad[all_ind[i]:all_ind[i+1]-1] = 1;
            pgrad[all_ind[i+1]:all_ind[i+2]-1] = 1;
        end
    end

    # if it divides by 2 and short - need to add 1s to the last section
    if mod(length(all_ind),2) == 0 # even number of indices
        if first_local_min         # mins are at odd, maxs are at even, last index is max
            pgrad[all_ind[length(all_ind)-1]:all_ind[length(all_ind)]] = 1
        else                         # maxs are at odd, mins are at even, last index is min
            ngrad[all_ind[length(all_ind)-1]:all_ind[length(all_ind)]] = 1
        end
    end

    sparse([1:N-1;2:N], [2:N;1:N-1], [pgrad[1:N-1];ngrad[1:N-1]], N, N)
end

# Converts adjacency matrix graphs into
# linear inequality constraints to use with Matlab's quadprog.
# Inequality constraints are used to localize the local extrema
# created in the interpolation process.
#
# The inequality constraints are as follows:
#
# For each vertex which is not a minimum nor a maximum, it will have
# one neighbor smaller than itself and one neighbor greater than
# itself.  Local minima are constrained to be smaller than both of
# their neighbors.  Local maxima are constrained to be greater than
# both of their neighbors.  Within each component between a local
# minimum and a local maximum, monotonicity is guaranteed.
#
# Code is adapted from Alec Jacobson, 2011 (jacobson@inf.ethz.ch)
#
# Jacobson, A., Baran, I., Jovan Popovi\'c, & Sorkine,
# O. (2011). Bounded biharmonic weights for real-time
# deformation. {ACM} Transactions on Graphics, 30(4),
# 1. http://dx.doi.org/10.1145/2010324.1964973
function in_eq_constraints(p::Persistence)
    A = monotonicity_graph(p)
    n = size(A,1)
    I,J = findn(A)
    nleq = length(I)
    m = max(n,nleq)
    Aleq = sparse([1:nleq;1:nleq], [I;J],
                  [ones(nleq);-ones(nleq)],
                  m, m)
    b = zeros(length(p.data))
    Aleq,b
end

# https://en.wikipedia.org/wiki/Quadratic_programming
function quadprog(H,f,A,b,Aeq,beq)
    m = Model()
    M,N = size(H)

    @variable(m, x[1:M])
    @objective(m, Min, (1/2*x'*H*x + f'*x)[1])

    @constraint(m, A*x .<= b)
    @constraint(m, Aeq*x .== beq)

    solve(m)
    getvalue(x)
end

function reconstruct_persistent(data::Vector{T}, threshold::Real,
                                smoothness::Integer = 2, weight = 1e-6) where {T<:Real}
    smoothness > 0 || error("Invalid smoothness requested: $(smoothness)")
    0 <= weight && weight <= 1 || error("Weight $(weight) not in the interval [0,1]")

    p = filter(find_persistence(data), threshold)
    Aeq,beq = eq_constraints(p)
    A,b = in_eq_constraints(p)

    N = length(p.data)
    laplace = spdiagm((-ones(N-1), [1;2ones(N-2);1], -ones(N-1)),
                      (-1, 0, 1))

    op = laplace^smoothness
    op = (op + op')/2

    f = zeros(N)

    if weight != 0
        op += spdiagm(weight*2ones(N))
        f = -2weight*p.data
    end

    quadprog(op, f, A, b, Aeq, beq)
end

export reconstruct_persistent


text:
# Factorize two different matrices (one real, one complex) and solve for multiple rhs and free memory
using Distributed
if nworkers()<3
    addprocs(3)
end
@everywhere using MUMPSjInv
using Test

include("getDivGrad.jl");

A  = getDivGrad(24,20,22);
A2 = getDivGrad(30,31,33);

n  = size(A,1);
n2 = size(A2,1);
A  = A + im*spdiagm(0=>rand(n));

nrhs = 10;
rhs = randn(n,nrhs) + im*randn(n,nrhs);
rhs2 = randn(n2,nrhs);

@test nworkers()>=3

println("Factorize complex matrix on worker 1")
@time begin
F1 = remotecall_fetch(factorMUMPS,workers()[1], A,1)
end
dump(F1)

println("Factorize real matrix on worker 2")
@time begin
	F2 = remotecall_fetch(factorMUMPS,workers()[2],A2,1)
end
dump(F2)

println("Solve complex system on worker 1")
@time begin
	x = remotecall_fetch(applyMUMPS,workers()[3],F1,rhs)
end
err = zeros(nrhs)
for i=1:nrhs
        err[i] =  norm(A*x[:,i]-rhs[:,i]) / norm(rhs[:,i])
end
@test maximum(err) < 1e-14

println("Solve real system on worker 2")
@time begin
	x2 = remotecall_fetch(applyMUMPS,workers()[2],F2,rhs2)
end
err = zeros(nrhs)
for i=1:nrhs
        err[i] =  norm(A2*x2[:,i]-rhs2[:,i]) / norm(rhs2[:,i]);
end
@test maximum(err) < 1e-14

println("Free memory on worker 1")
destroyMUMPS(F1)
# remotecall_fetch(destroyMUMPS,workers()[1],F1);

println("Free memory on worker 2")
destroyMUMPS(F2)
# remotecall_fetch(destroyMUMPS,workers()[2],F2);

println("DONE!")


text:
using DataDeps

register(DataDep(
  "hepatitis",
  "http://archive.ics.uci.edu/ml/datasets/Hepatitis",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/hepatitis/hepatitis.data",
  "9971f9d18abb4fc9dd6749eefb0f864b36fa7bdd59f54188c0b148a5c28214bc",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=1,
        feature_indices=2:20,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
#   http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing,
# software distributed under the License is distributed on an
# "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied.  See the License for the
# specific language governing permissions and limitations
# under the License.

"""
    SymbolicNode

SymbolicNode is the basic building block of the symbolic graph in MXNet.jl.

    (self :: SymbolicNode)(args :: SymbolicNode...)
    (self :: SymbolicNode)(; kwargs...)

Make a new node by composing `self` with `args`. Or the arguments
can be specified using keyword arguments.
"""
mutable struct SymbolicNode
  handle::MX_SymbolHandle
end

const SymbolicNodeOrReal = Union{SymbolicNode, Real}

@unfuse SymbolicNode  # for broadcasting

Base.unsafe_convert(::Type{MX_handle}, obj::SymbolicNode) =
  Base.unsafe_convert(MX_handle, obj.handle)
Base.convert(t::Type{MX_handle}, obj::SymbolicNode) = Base.unsafe_convert(t, obj)
Base.cconvert(t::Type{MX_handle}, obj::SymbolicNode) = Base.unsafe_convert(t, obj)

"""
    deepcopy(self :: SymbolicNode)

Make a deep copy of a SymbolicNode.
"""
function Base.deepcopy(self :: SymbolicNode)
  ref_hdr = Ref{MX_handle}(0)
  @mxcall(:MXSymbolCopy, (MX_handle, Ref{MX_handle}), self, ref_hdr)
  return SymbolicNode(MX_SymbolHandle(ref_hdr[]))
end

"""
    copy(self :: SymbolicNode)

Make a copy of a SymbolicNode. The same as making a deep copy.
"""
function Base.copy(self :: SymbolicNode)
  Base.deepcopy(self)
end

function (self::SymbolicNode)(args :: SymbolicNode...)
  s = deepcopy(self)
  _compose!(s, args...)
end
function (self::SymbolicNode)(;kwargs...)
  s = deepcopy(self)
  _compose!(s; kwargs...)
end

macro _list_symbol_info(self, func_name)
  quote
    ref_sz    = Ref{MX_uint}(0)
    ref_names = Ref{char_pp}(0)
    @mxcall($func_name, (MX_handle, Ref{MX_uint}, Ref{char_pp}),
            $(esc(self)), ref_sz, ref_names)
    narg = ref_sz[]
    names = unsafe_wrap(Array, ref_names[], narg)
    names = [Symbol(unsafe_string(x)) for x in names]
    return names
  end
end

"""
    list_arguments(self :: SymbolicNode)

List all the arguments of this node. The argument for a node contains both
the inputs and parameters. For example, a `FullyConnected` node will
have both data and weights in its arguments. A composed node (e.g. a MLP) will
list all the arguments for intermediate nodes.

Returns a list of symbols indicating the names of the arguments.
"""
function list_arguments(self :: SymbolicNode)
  @_list_symbol_info(self, :MXSymbolListArguments)
end

"""
    list_outputs(self :: SymbolicNode)

List all the outputs of this node.

Returns a list of symbols indicating the names of the outputs.
"""
function list_outputs(self :: SymbolicNode)
  @_list_symbol_info(self, :MXSymbolListOutputs)
end

"""
    list_auxiliary_states(self :: SymbolicNode)

List all auxiliary states in the symbool.

Auxiliary states are special states of symbols that do not corresponds to an argument,
and do not have gradient. But still be useful for the specific operations.
A common example of auxiliary state is the moving_mean and moving_variance in BatchNorm.
Most operators do not have Auxiliary states.

Returns a list of symbols indicating the names of the auxiliary states.
"""
function list_auxiliary_states(self :: SymbolicNode)
  @_list_symbol_info(self, :MXSymbolListAuxiliaryStates)
end

"""
    get_internals(self :: SymbolicNode)

Get a new grouped `SymbolicNode` whose output contains all the internal outputs of
this `SymbolicNode`.
"""
function get_internals(self :: SymbolicNode)
  ref_hdr = Ref{MX_handle}(0)
  @mxcall(:MXSymbolGetInternals, (MX_handle, Ref{MX_handle}), self, ref_hdr)
  return SymbolicNode(MX_SymbolHandle(ref_hdr[]))
end

"""
    get_children(x::SymbolicNode)

Gets a new grouped `SymbolicNode` whose output contains inputs to output
nodes of the original symbol.

```julia
julia> x = mx.Variable(:x)
MXNet.mx.SymbolicNode x

julia> y = mx.Variable(:y)
MXNet.mx.SymbolicNode y

julia> z = x + y
MXNet.mx.SymbolicNode _plus1

julia> a |> mx.get_children |> mx.list_outputs
2-element Array{Symbol,1}:
 :x
 :y
```
"""
function get_children(x::SymbolicNode)
  hdl = Ref{MX_handle}(C_NULL)
  @mxcall(:MXSymbolGetChildren, (MX_handle, Ref{MX_handle}), x, hdl)
  sym = hdl[] |> MX_SymbolHandle |> SymbolicNode
  isempty(list_outputs(sym)) ? nothing : sym
end

"""
    get_attr(self :: SymbolicNode, key :: Symbol)

Get attribute attached to this `SymbolicNode` belonging to key.

Returns the value belonging to key as a `Nullable`.
"""
function get_attr(self :: SymbolicNode, key :: Symbol)
  key_s = string(key)
  ref_out = Ref{Cstring}()
  ref_success = Ref{Cint}(-1)
  @mxcall(:MXSymbolGetAttr, (MX_handle, Cstring, Ref{Cstring}, Ref{Cint}),
          self, key_s, ref_out, ref_success)
  if ref_success[] == 1
    return Nullable{String}(unsafe_string(ref_out[]))
  else
    return Nullable{String}()
  end
end

"""
    list_attr(self :: SymbolicNode)

Get all attributes from a symbol.

Returns a dictionary of attributes.
"""
function list_attr(self :: SymbolicNode)
  ref_sz    = Ref{MX_uint}(0)
  ref_strings = Ref{char_pp}(0)
  @mxcall(:MXSymbolListAttrShallow, (MX_handle, Ref{MX_uint}, Ref{char_pp}),
            self, ref_sz, ref_strings)
  narg = 2*ref_sz[]
  strings = unsafe_wrap(Array, ref_strings[], narg)
  out = Dict{Symbol, String}()
  for i in 1:2:narg
    key = Symbol(unsafe_string(strings[i]))
    value = unsafe_string(strings[i+1]) # Creates a copy of string
    out[key] = value
  end
  return out
end

"""
    list_all_attr(self :: SymbolicNode)

Get all attributes from the symbol graph.

Returns a dictionary of attributes.
"""
function list_all_attr(self :: SymbolicNode)
  ref_sz    = Ref{MX_uint}(0)
  ref_strings = Ref{char_pp}(0)
  @mxcall(:MXSymbolListAttr, (MX_handle, Ref{MX_uint}, Ref{char_pp}),
            self, ref_sz, ref_strings)
  narg = 2*ref_sz[]
  strings = unsafe_wrap(Array, ref_strings[], narg)
  out = Dict{Symbol, String}()
  for i in 1:2:narg
    key = Symbol(unsafe_string(strings[i]))
    value = unsafe_string(strings[i+1])
    out[key] = value
  end
  return out
end

"""
    set_attr(self:: SymbolicNode, key :: Symbol, value :: AbstractString)

Set the attribute key to value for this `SymbolicNode`.

!!! note
    It is encouraged not to call this function directly, unless you know exactly what you are doing. The
    recommended way of setting attributes is when creating the `SymbolicNode`. Changing
    the attributes of a `SymbolicNode` that is already been used somewhere else might
    cause unexpected behavior and inconsistency.
"""
function set_attr(self :: SymbolicNode, key :: Symbol, value :: AbstractString)
  key_s = string(key)
  value_s = String(value)

  @mxcall(:MXSymbolSetAttr, (MX_handle, Cstring, Cstring), self, key_s, value_s)
end

"""
    get_name(self :: SymbolicNode)

Get the name of the symbol.

    julia> x = mx.Variable(:data)
    julia> mx.get_name(x)
    :data

    julia> y = mx.FullyConnected(x, num_hidden = 128)
    julia> mx.get_name(y)
    :fullyconnected0
"""
function get_name(self :: mx.SymbolicNode)
    name = Ref{mx.char_p}(0)
    success = Ref(0)
    @mxcall(:MXSymbolGetName, (MX_handle, Ref{char_p}, Ref{Int}), self.handle.value, name, success)
    @assert success[] != -1

    str = name[]
    if str == C_NULL  # e.g. the symbol returned via get_internals
        string(self.handle.value)
    else
        Symbol(unsafe_string(str))
    end
end

Base.show(io::IO, sym::SymbolicNode) =
  print(io, "$(typeof(sym)) $(get_name(sym))")

import Base: print

function print(io::IO, sym::SymbolicNode)
  out = Ref{mx.char_p}(C_NULL)
  @mx.mxcall(:MXSymbolPrint, (mx.MX_SymbolHandle, Ref{mx.char_p}), sym.handle, out)
  print(io, unsafe_string(out[]))
end

print(sym::SymbolicNode) = print(STDOUT, sym)

"""
    print([io::IO], sym::SymbolicNode)

Print the content of symbol, used for debug.

```julia
julia> layer = @mx.chain mx.Variable(:data)           =>
         mx.FullyConnected(name=:fc1, num_hidden=128) =>
         mx.Activation(name=:relu1, act_type=:relu)
MXNet.mx.SymbolicNode(MXNet.mx.MX_SymbolHandle(Ptr{Void} @0x000055b29b9c3520))

julia> print(layer)
Symbol Outputs:
        output[0]=relu1(0)
Variable:data
Variable:fc1_weight
Variable:fc1_bias
--------------------
Op:FullyConnected, Name=fc1
Inputs:
        arg[0]=data(0) version=0
        arg[1]=fc1_weight(0) version=0
        arg[2]=fc1_bias(0) version=0
Attrs:
        num_hidden=128
--------------------
Op:Activation, Name=relu1
Inputs:
        arg[0]=fc1(0)
Attrs:
        act_type=relu
```
"""
print

"""
    grad(self :: SymbolicNode, wrt :: Vector{SymbolicNode})

Get the autodiff gradient of the current `SymbolicNode`. This function can
only be used if the current symbol is a loss function.

# Arguments:
* `self::SymbolicNode`: current node.
* `wrt::Vector{Symbol}`: the names of the arguments to the gradient.

Returns a gradient symbol of the corresponding gradient.
"""
function grad(self :: SymbolicNode, wrt :: Vector{Symbol})
  hdr_ref = Ref{MX_handle}(0)
  keys = String[string(key) for key in wrt]

  @mxcall(:MXSymbolGrad, (MX_handle, MX_uint, char_pp, Ptr{MX_handle}), self, length(keys), keys, hdr_ref)
  return SymbolicNode(MX_SymbolHandle(hdr_ref[]))
end

"""
    Variable(name :: Union{Symbol, AbstractString})

Create a symbolic variable with the given name. This is typically used as a placeholder.
For example, the data node, acting as the starting point of a network architecture.

# Arguments
* Dict{Symbol, AbstractString} attrs: The attributes associated with this `Variable`.
"""
function Variable(name :: Union{Symbol, AbstractString}; attrs = Dict())
  attrs = convert(Dict{Symbol, AbstractString}, attrs)
  hdr_ref = Ref{MX_handle}(0)
  @mxcall(:MXSymbolCreateVariable, (char_p, Ref{MX_handle}), name, hdr_ref)
  node = SymbolicNode(MX_SymbolHandle(hdr_ref[]))
  for (k, v) in attrs
    set_attr(node, k, v)
  end
  node
end

"""
    @var <symbols>...

A handy macro for creating `mx.Variable`.

```julia
julia> x = @mx.var x
MXNet.mx.SymbolicNode x

julia> x, y, z = @mx.var x y z
(MXNet.mx.SymbolicNode x, MXNet.mx.SymbolicNode y, MXNet.mx.SymbolicNode z)
```
"""
macro var(n::Symbol)
  Expr(:call, :Variable, QuoteNode(n))
end

macro var(names::Symbol...)
  Expr(:tuple, map(n -> Expr(:call, :Variable, QuoteNode(n)), names)...)
end

"""
    Group(nodes :: SymbolicNode...)

Create a `SymbolicNode` by grouping nodes together.
"""
function Group(nodes :: SymbolicNode...)
  handles = MX_handle[nodes...]
  ref_hdr = Ref{MX_handle}(0)
  @mxcall(:MXSymbolCreateGroup, (MX_uint, Ptr{MX_handle}, Ref{MX_handle}),
          length(handles), handles, ref_hdr)
  SymbolicNode(MX_SymbolHandle(ref_hdr[]))
end

function _build_shapes(shape_size::MX_uint, shape_ndim::Ptr{MX_uint}, shape_data::Ptr{Ptr{MX_uint}})
  shape_ndim = unsafe_wrap(Array, shape_ndim, shape_size)
  shape_data = unsafe_wrap(Array, shape_data, shape_size)
  shapes = map(1:shape_size) do i
    my_shape = unsafe_wrap(Array, shape_data[i], shape_ndim[i])
    tuple(flipdim(Int[my_shape...],1)...)
  end
  convert(Vector{Tuple}, shapes)
end

function _infer_shape(self, keys, indptr, sdata)
  ref_arg_shape_size = Ref{MX_uint}(0)
  ref_arg_shape_ndim = Ref{Ptr{MX_uint}}(0)
  ref_arg_shape_data = Ref{Ptr{Ptr{MX_uint}}}(0)
  ref_out_shape_size = Ref{MX_uint}(0)
  ref_out_shape_ndim = Ref{Ptr{MX_uint}}(0)
  ref_out_shape_data = Ref{Ptr{Ptr{MX_uint}}}(0)
  ref_aux_shape_size = Ref{MX_uint}(0)
  ref_aux_shape_ndim = Ref{Ptr{MX_uint}}(0)
  ref_aux_shape_data = Ref{Ptr{Ptr{MX_uint}}}(0)
  ref_complete       = Ref{Cint}(0)
  @mxcall(:MXSymbolInferShape,
          (MX_handle, MX_uint, char_pp, Ptr{MX_uint}, Ptr{MX_uint},
           Ref{MX_uint}, Ref{Ptr{MX_uint}}, Ref{Ptr{Ptr{MX_uint}}},
           Ref{MX_uint}, Ref{Ptr{MX_uint}}, Ref{Ptr{Ptr{MX_uint}}},
           Ref{MX_uint}, Ref{Ptr{MX_uint}}, Ref{Ptr{Ptr{MX_uint}}},
           Ref{Cint}),
          self, length(indptr)-1, keys, indptr, sdata,
          ref_arg_shape_size, ref_arg_shape_ndim, ref_arg_shape_data,
          ref_out_shape_size, ref_out_shape_ndim, ref_out_shape_data,
          ref_aux_shape_size, ref_aux_shape_ndim, ref_aux_shape_data,
          ref_complete)
  if ref_complete[] == 0
    return (nothing, nothing, nothing)
  else
    return (
      _build_shapes(ref_arg_shape_size[], ref_arg_shape_ndim[], ref_arg_shape_data[]),
      _build_shapes(ref_out_shape_size[], ref_out_shape_ndim[], ref_out_shape_data[]),
      _build_shapes(ref_aux_shape_size[], ref_aux_shape_ndim[], ref_aux_shape_data[])
    )
  end
end

"""
    infer_shape(self :: SymbolicNode, args...)
    infer_shape(self :: SymbolicNode; kwargs...)

Do shape inference according to the input shapes. The input shapes could be provided
as a list of shapes, which should specify the shapes of inputs in the same order as
the arguments returned by [`list_arguments`](@ref). Alternatively, the shape information
could be specified via keyword arguments.

Returns a 3-tuple containing shapes of all the arguments, shapes of all the outputs and
shapes of all the auxiliary variables. If shape inference failed due to incomplete
or incompatible inputs, the return value will be `(nothing, nothing, nothing)`.
"""
function infer_shape(self :: SymbolicNode; kwargs...)
  sdata  = MX_uint[]
  indptr = MX_uint[0]
  for (k,v) in kwargs
    append!(sdata, flipdim([v...],1))
    push!(indptr, length(sdata))
  end
  keys = AbstractString[string(x[1]) for x in kwargs]
  _infer_shape(self, keys, indptr, sdata)
end
function infer_shape(self :: SymbolicNode, args :: Union{Tuple, Void}...)
  sdata  = MX_uint[]
  indptr = MX_uint[0]
  for arg in args
    if isa(arg, Void); continue; end
    append!(sdata, flipdim([arg...],1))
    push!(indptr, length(sdata))
  end
  keys = Ptr{char_p}(0)
  _infer_shape(self, keys, indptr, sdata)
end

function _infer_type(self, keys, arg_type_data)
  ref_in_type_size  = Ref{MX_uint}()
  ref_in_type_data  = Ref{Ptr{Cint}}()
  ref_out_type_size = Ref{MX_uint}()
  ref_out_type_data = Ref{Ptr{Cint}}()
  ref_aux_type_size = Ref{MX_uint}()
  ref_aux_type_data = Ref{Ptr{Cint}}()
  ref_complete      = Ref{Cint}()

  @mxcall(:MXSymbolInferType,
          (MX_handle, MX_uint, char_pp, Ptr{Cint},
           Ref{MX_uint}, Ref{Ptr{Cint}},
           Ref{MX_uint}, Ref{Ptr{Cint}},
           Ref{MX_uint}, Ref{Ptr{Cint}},
           Ref{Cint}),
          self, length(arg_type_data)-1, keys, arg_type_data,
          ref_in_type_size, ref_in_type_data,
          ref_out_type_size, ref_out_type_data,
          ref_aux_type_size, ref_aux_type_data,
          ref_complete)

  if ref_complete[] == 0
    return (nothing, nothing, nothing)
  else
    in_type = unsafe_wrap(Array, ref_in_type_data[], ref_in_type_size[])
    out_type = unsafe_wrap(Array, ref_out_type_data[], ref_out_type_size[])
    aux_type = unsafe_wrap(Array, ref_aux_type_data[], ref_aux_type_size[])
    return ([fromTypeFlag(TypeFlag(t)) for t in in_type],
            [fromTypeFlag(TypeFlag(t)) for t in out_type],
            [fromTypeFlag(TypeFlag(t)) for t in aux_type])
  end
end

"""
    infer_type(self :: SymbolicNode; kwargs...)
    infer_type(self :: SymbolicNode, args...)

Do type inference according to the input types. The input types could be provided
as a list of types, which should specify the types of inputs in the same order as
the arguments returned by [`list_arguments`](@ref). Alternatively, the type information
could be specified via keyword arguments.

Returns a 3-tuple containing types of all the arguments, types of all the outputs and
types of all the auxiliary variables. If type inference failed due to incomplete
or incompatible inputs, the return value will be `(nothing, nothing, nothing)`.
"""
function infer_type(self :: SymbolicNode; kwargs...)
  types = Cint[toTypeFlag(x[2]) for x in kwargs]
  keys = AbstractString[string(x[1]) for x in kwargs]
  _infer_type(self, keys, types)
end

function infer_type(self :: SymbolicNode, args :: Union{Tuple, Void}...)
  types = Cint[]
  keys = Ptr{char_p}(0)

  for arg in args
    if isa(arg, Void); continue; end
    push!(types, toTypeFlag(arg))
  end
  _infer_type(self, keys, types)
end

"""
    getindex(self :: SymbolicNode, idx :: Union{Int, Base.Symbol, AbstractString})

Get a node representing the specified output of this node. The index could be
a symbol or string indicating the name of the output, or a 1-based integer
indicating the index, as in the list of [`list_outputs`](@ref).
"""
function Base.getindex(self :: SymbolicNode, idx :: Union{Base.Symbol, AbstractString})
  idx   = Symbol(idx)
  i_idx = find(idx .== list_outputs(self))
  @assert(length(i_idx) > 0, "Cannot find output with name '$idx'")
  @assert(length(i_idx) < 2, "Found duplicated output with name '$idx'")
  Base.getindex(self, i_idx[1])
end
function Base.getindex(self :: SymbolicNode, idx :: Int)
  ref_hdr = Ref{MX_handle}(0)
  # note Julia is 1-based, while MXNet is 0-based
  @mxcall(:MXSymbolGetOutput, (MX_handle, MX_uint, Ref{MX_handle}), self, idx-1, ref_hdr)
  return SymbolicNode(MX_SymbolHandle(ref_hdr[]))
end

import Base: +

"""
    +(args...)
    .+(args...)

Elementwise summation of `SymbolicNode`.
"""
function +(x::SymbolicNode, ys::SymbolicNodeOrReal...)
  ret = x
  for y ∈ ys
    if y isa SymbolicNode
      ret = _plus(ret, y)
    else
      ret = _plus_scalar(ret, scalar=MX_float(y))
    end
  end
  ret
end

+(s::Real, x::SymbolicNode, ys::SymbolicNodeOrReal...) = +(x + s, ys...)

broadcast_(::typeof(+), x::SymbolicNode, ys::SymbolicNodeOrReal...) = +(x, ys...)
broadcast_(::typeof(+), s::Real, x::SymbolicNode, ys::SymbolicNodeOrReal...) = +(x + s, ys...)

import Base: -

"""
    -(x, y)
    .-(x, y)

Elementwise substraction of `SymbolicNode`.
Operating with `Real` is available.
"""
x::SymbolicNode - y::SymbolicNode = _minus(x, y)
x::SymbolicNode - s::Real         = _minus_scalar(x,  scalar=MX_float(s))
s::Real         - x::SymbolicNode = _rminus_scalar(x, scalar=MX_float(s))

-(x::SymbolicNode) = 0 - x

broadcast_(::typeof(-), x::SymbolicNode, y::SymbolicNodeOrReal) = x - y
broadcast_(::typeof(-), s::Real, x::SymbolicNode) = s - x

import Base: *

"""
    .*(x, y)

Elementwise multiplication of `SymbolicNode`.
"""
x::SymbolicNode * s::Real = _mul_scalar(x, scalar=MX_float(s))
s::Real * x::SymbolicNode = _mul_scalar(x, scalar=MX_float(s))

function broadcast_(::typeof(*), x::SymbolicNode, ys::SymbolicNodeOrReal...)
  ret = x
  for y in ys
    if y isa SymbolicNode
      ret = _mul(ret, y)
    else
      ret = _mul_scalar(ret, scalar=MX_float(y))
    end
  end
  ret
end

broadcast_(::typeof(*), s::Real, x::SymbolicNode, ys::SymbolicNodeOrReal...) =
  broadcast_(*, x * s, ys...)

import Base: /

"""
    ./(x, y)

* Elementwise dividing a `SymbolicNode` by a scalar or another `SymbolicNode`
of the same shape.

* Elementwise divide a scalar by an `SymbolicNode`.

* Matrix division (solving linear systems) is not implemented yet.
"""
x::SymbolicNode / s::Real = _DivScalar(x, scalar=MX_float(s))

broadcast_(::typeof(/), x::SymbolicNode, y::SymbolicNode) = _div(x, y)
broadcast_(::typeof(/), x::SymbolicNode, s::Real) = _div_scalar(x,  scalar=MX_float(s))
broadcast_(::typeof(/), s::Real, x::SymbolicNode) = _rdiv_scalar(x, scalar=MX_float(s))

import Base: ^

"""
    .^(x, y)

Elementwise power of `SymbolicNode` and `NDArray`.
Operating with `Real` is available.
"""
^

broadcast_(::typeof(^), x::SymbolicNode, y::SymbolicNode) = _power(x, y)
broadcast_(::typeof(^), x::SymbolicNode, s::Real) = _power_scalar(x,  scalar=MX_float(s))
broadcast_(::typeof(^), s::Real, x::SymbolicNode) = _rpower_scalar(x, scalar=MX_float(s))

broadcast_(::typeof(^), ::Irrational{:e}, x::SymbolicNode) = exp(x)
broadcast_(::typeof(^), x::SymbolicNode, s::Irrational) =
  _power_scalar(x, scalar=MX_float(s))
broadcast_(::typeof(^), s::Irrational, x::SymbolicNode) =
  _rpower_scalar(x, scalar=MX_float(s))

function _compose!(node :: SymbolicNode; kwargs...)
  name     = char_p(0)
  arg_keys = AbstractString[]
  arg_vals = MX_handle[]

  for (k,v) in kwargs
    if k == :name
      name = string(v)
    else
      @assert(isa(v, SymbolicNode), "Compose expect `SymbolicNode` as arguments")
      push!(arg_keys, string(k))
      push!(arg_vals, v)
    end
  end

  @mxcall(:MXSymbolCompose,
          (MX_handle, char_p, MX_uint, Ptr{char_p}, Ptr{MX_handle}),
          node, name, length(arg_keys), arg_keys, arg_vals)
  return node
end
function _compose!(node :: SymbolicNode, args::SymbolicNode...)
  _compose!(node, char_p(0), args...)
end
function _compose!(node :: SymbolicNode, name :: Union{Base.Symbol, char_p}, args::SymbolicNode...)
  if isa(name, Base.Symbol); name = string(name); end
  arg_keys = Ptr{char_p}(0)
  arg_vals = MX_handle[args...]

  @mxcall(:MXSymbolCompose,
          (MX_handle, char_p, MX_uint, Ptr{char_p}, Ptr{MX_handle}),
          node, name, length(arg_vals), arg_keys, arg_vals)
  return node
end

"""
    to_json(self :: SymbolicNode)

Convert a `SymbolicNode` into a JSON string.
"""
function to_json(self :: SymbolicNode)
  ref_json = Ref{char_p}(0)
  @mxcall(:MXSymbolSaveToJSON, (MX_handle, Ref{char_p}), self, ref_json)
  return unsafe_string(ref_json[])
end

"""
    from_json(repr :: AbstractString, ::Type{SymbolicNode})

Load a `SymbolicNode` from a JSON string representation.
"""
function from_json(repr :: AbstractString, ::Type{SymbolicNode})
  ref_hdr = Ref{MX_handle}(0)
  @mxcall(:MXSymbolCreateFromJSON, (char_p, Ref{MX_handle}), repr, ref_hdr)
  return SymbolicNode(MX_SymbolHandle(ref_hdr[]))
end

"""
    load(filename :: AbstractString, ::Type{SymbolicNode})

Load a `SymbolicNode` from a JSON file.
"""
function load(filename :: AbstractString, ::Type{SymbolicNode})
  ref_hdr = Ref{MX_handle}(0)
  @mxcall(:MXSymbolCreateFromFile, (char_p, Ref{MX_handle}), filename, ref_hdr)
  return SymbolicNode(MX_SymbolHandle(ref_hdr[]))
end

"""
    save(filename :: AbstractString, node :: SymbolicNode)

Save a `SymbolicNode` to a JSON file.
"""
function save(filename :: AbstractString, node :: SymbolicNode)
  @mxcall(:MXSymbolSaveToFile, (MX_handle, char_p), node, filename)
end

import Base: reshape

"""
    reshape(sym::SymbolicNode, dim; reverse=false, name)
    reshape(sym::SymbolicNode, dim...; reverse=false, name)

Reshape SymbolicNode operator

Some dimensions of the shape can take special values from the set
{0, -1, -2, -3, -4}.
The significance of each is explained below:

- `0`  copy this dimension from the input to the output shape.

  Example:

  - input shape = (2,3,4), shape = (4,0,2), output shape = (4,3,2)
  - input shape = (2,3,4), shape = (2,0,0), output shape = (2,3,4)

- `-1` infers the dimension of the output shape by using the remainder of the
  input dimensions keeping the size of the new array same as that of the input
  array. At most one dimension of shape can be -1.

  Example:

  - input shape = (2,3,4), shape = (6,1,-1), output shape = (6,1,4)
  - input shape = (2,3,4), shape = (3,-1,8), output shape = (3,1,8)
  - input shape = (2,3,4), shape=(-1,), output shape = (24,)

- `-2` copy all/remainder of the input dimensions to the output shape.

  Example:

  - input shape = (2,3,4), shape = (-2,), output shape = (2,3,4)
  - input shape = (2,3,4), shape = (2,-2), output shape = (2,3,4)
  - input shape = (2,3,4), shape = (-2,1,1), output shape = (2,3,4,1,1)

- `-3` use the product of two consecutive dimensions of the input shape as the
  output dimension.

  Example:

  - input shape = (2,3,4), shape = (-3,4), output shape = (6,4)
  - input shape = (2,3,4,5), shape = (-3,-3), output shape = (6,20)
  - input shape = (2,3,4), shape = (0,-3), output shape = (2,12)
  - input shape = (2,3,4), shape = (-3,-2), output shape = (6,4)

- `-4` split one dimension of the input into two dimensions passed subsequent
  to -4 in shape (can contain -1).

  Example:

  - input shape = (2,3,4), shape = (-4,1,2,-2), output shape = (1,2,3,4)
  - input shape = (2,3,4), shape = (2,-4,-1,3,-2), output shape = (2,1,3,4)

If the argument `reverse` is set to `1`, then the special values are inferred
from right to left.

  Example:

  - with `reverse=false`, for input shape = (10,5,4), shape = (-1,0),
    output shape would be (40,5)
  - with `reverse=true`, output shape will be (50,4).
"""
reshape(sym::SymbolicNode, dim::NTuple{N, Integer}; kwargs...) where {N} =
  _reshape(sym, dim; kwargs...)
reshape(sym::SymbolicNode, dim::Integer...; kwargs...) =
  _reshape(sym, dim; kwargs...)

@inline function _reshape(sym::SymbolicNode, dim::NTuple{N, Integer};
                          reverse::Bool=false, name::String="") where N
  op = _get_cached_libmx_op_handle("reshape")
  node = _create_atomic_symbol(op.value, ["shape", "reverse"],
                               [dump_mx_param(dim), dump_mx_param(!reverse)])
  name = get!(DEFAULT_NAME_MANAGER, name, "reshape")
  _compose!(node, name=name, data=sym)
end

################################################################################
# Atomic SymbolicNode functions dynamically imported from libmxnet
################################################################################
@inline function _create_atomic_symbol(creator::MX_handle, keys::Vector{String},
                                       vals::Vector{String})
  ref_sym_hdr = Ref{MX_handle}(C_NULL)
  @mxcall(:MXSymbolCreateAtomicSymbol,
          (MX_handle, MX_uint, Ptr{char_p}, Ptr{char_p}, Ref{MX_handle}),
          creator, length(keys), keys, vals, ref_sym_hdr)
  SymbolicNode(MX_SymbolHandle(ref_sym_hdr[]))
end

@inline function _create_atomic_symbol(creator::MX_handle, keys::Vector{String},
                                       vals::Vector{String},
                                       attrs::Dict{Symbol, String})
  node = _create_atomic_symbol(creator, keys, vals)
  # set attrs
  for (k, v) in attrs
    set_attr(node, k, v)
  end
  node
end

function _define_atomic_symbol_creator(name :: String)
  handle = _get_libmx_op_handle(name)
  f_desc, key_narg = _get_libmx_op_description(name, handle)

  f_desc *= "* `name::Symbol`: The name of the `SymbolicNode`. (e.g. `:my_symbol`), optional.\n"
  f_desc *= "* `attrs::Dict{Symbol, AbstractString}`: The attributes associated with this `SymbolicNode`.\n\n"

  func_name = Symbol(name)
  func_def = quote
  function $func_name(::Type{SymbolicNode}, args::SymbolicNode...; kwargs...)
    idx = findfirst(x -> x[1] == :name, kwargs)
    if idx > 0
      name = kwargs[idx][2]
    else
      name = ""
    end

    # XXX: hacky way of solving the problem that the arguments of `dot` should be swapped
    # See https://github.com/dmlc/MXNet.jl/issues/55
    if $name == "dot"
      args = reverse(args)
    end

    # XXX: hacky way of solving the semantic difference of the axes parameter in Julia
    # and in libmxnet.
    # See https://github.com/dmlc/MXNet.jl/pull/123
    if $name == "transpose"
      kwargs = Any[key != :axes ? (key, arg) : (key, reverse(map(i->length(arg)-i, arg))) for (key, arg) in kwargs]
    end

    param_keys = String[]
    param_vals = String[]
    symbol_kws = Dict{Symbol, SymbolicNode}()
    attrs = Dict{Symbol, String}()

    $(if key_narg != ""
      quote
        if !in($key_narg, param_keys)
          push!(param_keys, $key_narg)
          push!(param_vals, string(length(args)))
        end
      end
    end)

    for (k,v) in kwargs
      if k == :name; continue; end
      if isa(v, SymbolicNode)
        symbol_kws[k] = v
      elseif k == :attrs
        if isa(v, Dict)
          attrs = convert(Dict{Symbol, String}, v)
        else
          throw(ArgumentError("attrs needs to be a Dictionary"))
        end
      else
        push!(param_keys, string(k))
        push!(param_vals, dump_mx_param(v))
      end
    end

    if length(args) > 1 && length(symbol_kws) != 0
      @assert(false, $name * " only accepts SymbolicNode either as positional or keyword arguments with optional positional `data` argument, not both.")
    end
    $(if key_narg != ""
      quote
        if length(symbol_kws) > 0
          @assert(false, $name * " takes variable number of SymbolicNode arguments, " *
                         "please pass input Symbols via positional arguments, instead of keyword arguments.")
        end
      end
    end)

    local op = _get_cached_libmx_op_handle($name)
    node = _create_atomic_symbol(op.value, param_keys, param_vals, attrs)

    # generate a new name for the new symbol if user not provided in kwargs
    hint = lowercase($name)
    name = get!(DEFAULT_NAME_MANAGER, name, hint)

    if length(symbol_kws) == 0
      _compose!(node, name, args...)
    elseif length(args) == 1
      _compose!(node; name=name, data=args[1], symbol_kws...)
    else
      _compose!(node; name=name, symbol_kws...)
    end

    return node
  end # function
  end # quote

  func_def2 = quote
  @doc $f_desc ->
  function $func_name(args::SymbolicNode...; kwargs...)
    $func_name(SymbolicNode, args...; kwargs...)
  end # function
  end # quote

  return quote
    $func_def
    $func_def2
  end
end

macro _import_atomic_symbol_creators()
  # XXX: those are operators defined for NDArray, we exclude them here
  # because the calling convention for the type signature is not strong
  # enough to disambiguate the method for NDArray and SymbolicNode
  const ignored_ops = ["_set_value", "reshape"]  # in lowercase

  op_names = _get_libmx_op_names()
  func_exprs = map(op_names) do name
    if lowercase(name) ∉ ignored_ops
      expr = _define_atomic_symbol_creator(name)
    end
  end

  esc(quote
    $(func_exprs...)
  end)
end

@_import_atomic_symbol_creators()

################################################################################
# Utility macros to chain up symbols
################################################################################
macro chain(layers)
    exprs = []
    last_layer = nothing

    function _chain_layer(layer, last_layer)
        if isa(last_layer, Void)
            return esc(layer)
        else
            if @capture(layer, f_(x__))
                x′ = esc.(x)
                return :($f($last_layer, $(x′...)))
            else
                throw(AssertionError("$layer is not a valid function call and cannot be chained."))
            end
        end
    end

    while true
        if @capture(layers, l1_=>l2_)
            new_layer = gensym()
            push!(exprs, :($new_layer = $(_chain_layer(l1, last_layer))))
            last_layer = new_layer
            layers = l2
        else
            push!(exprs, _chain_layer(layers, last_layer))
            break
        end
    end
    Expr(:block, exprs...)
end

max_stars_repo_name: Vikas-kum/incubator-mxnet

max_issues_repo_name: Vikas-kum/incubator-mxnet

max_forks_repo_name: Vikas-kum/incubator-mxnet


text:
import Base.LinAlg: axpy!

export svdvals_tr

"""
Matrix of the form
 d1          a_1
    d2       a_2
      ...    ...
         d_l a_l
             d_l+1 e_l+1
                   ...  e_k-1
                        d_k
"""
type BrokenArrowBidiagonal{T} <: AbstractMatrix{T}
    dv::Vector{T}
    av::Vector{T}
    ev::Vector{T}
end

Base.size(B::BrokenArrowBidiagonal) = (n=length(B.dv); (n, n))
function Base.size(B::BrokenArrowBidiagonal, n::Int)
    if n==1 || n==2
        return length(B.dv)
    else
        throw(ArgumentError("invalid dimension $n"))
    end
end

function Base.full{T}(B::BrokenArrowBidiagonal{T})
    n = size(B, 1)
    k = length(B.av)
    M = zeros(T, n, n)
    for i=1:n
        M[i, i] = B.dv[i]
    end
    for i=1:k
        M[i,k+1] = B.av[i]
    end
    for i=k+1:n-1
        M[i, i+1] = B.ev[i-k]
    end
    M
end

Base.svdfact(B::BrokenArrowBidiagonal) = svdfact(full(B)) #XXX This can be much faster

"""
Partial factorization object which is an approximation of a matrix

    A ≈ P * [B 0; 0 β] * Q
"""
type PartialFactorization{T, Tr} <: Factorization{T}
    P :: Matrix{T}
    Q :: Matrix{T}
    B :: AbstractMatrix{Tr}
    β :: Tr
end

"""
The thick-restarted variant of Golub-Kahan-Lanczos bidiagonalization

# Inputs

- `A` : The matrix or matrixlike object whose singular values are desired
- `q` : The starting guess vector in the range of `A`.
        The length of `q` should be the number of columns in `A`.
        Default: A random unit vector.
- `l` : The number of singular values requested.
        Default: 6
- `k` : The maximum number of Lanczos vectors to compute before restarting.
        Default: `2*l`
- `j` : The number of vectors to keep at the end of the restart.
        Default: `l`. We don't recommend j < l.

# Keyword inputs

- `maxiter`: Maximum number of iterations to run
             Default: `minimum(size(A))`
- `tol`    : Maximum absolute error in each desired singular value.
             Default: `√eps()`
- `reltol` : Maximum error in each desired singular value relative to the
             estimated norm of the input matrix.  Default: `√eps()`
- `restart`: Which restarting algorithm to use. Valid choices are:
             - `:ritz`: Thick restart with Ritz values [Wu2000] (default)
             - `:harmonic`: Restart with harmonic Ritz values [Baglama2005]

# Output

- `Σ`: A list of the desired singular values
- `L`: The computed partial factorization of A

# Implementation notes

The implementation of thick restarting follows closely that of SLEPc as
described in [Hernandez2008].

# References

@article{Wu2000,
    author = {Wu, Kesheng and Simon, Horst},
    journal = {SIAM Journal on Matrix Analysis and Applications},
    number = 2,
    pages = {602--616},
    title = {Thick-Restart {L}anczos Method for Large Symmetric Eigenvalue Problems},
    volume = 22,
    year = 2000
}

@article{Baglama2005,
    author = {Baglama, James and Reichel, Lothar},
    doi = {10.1137/04060593X},
    journal = {SIAM Journal on Scientific Computing},
    number = 1,
    pages = {19--42},
    title = {Augmented Implicitly Restarted {L}anczos Bidiagonalization Methods},
    volume = 27,
    year = 2005
}

@article{Hernandez2008,
    author = {Hern\'{a}ndez, Vicente and Rom\'{a}n, Jos\'{e} E and Tom\'{a}s,
    Andr\'{e}s},
    journal = {Electronic Transactions on Numerical Analysis},
    pages = {68--85},
    title = {A Robust and Efficient Parallel {SVD} Solver based on Restarted
        {L}anczos Bidiagonalization},
    url = {http://etna.mcs.kent.edu/volumes/2001-2010/vol31/abstract.php?vol=31\&pages=68-85},
    volume = 31,
    year = 2008
}

"""
function svdvals_tr(A, q::AbstractVector, l::Int=6, k::Int=2l,
    j::Int=l;
    maxiter::Int=minimum(size(A)), tol::Real=√eps(), reltol::Real=√eps(),
    method::Symbol=:ritz)

    @assert k>l
    L = build(A, q, k)

    local F
    for i=1:maxiter
        #@assert size(L.B) == (k, k)
        F = svdfact(L.B)
        i==1 && @assert eltype(F)==eltype(q)
        if method == :ritz
            thickrestart!(A, L, F, j)
        elseif method == :harmonic
            harmonicrestart!(A, L, F, j)
        else
            throw(ArgumentError("Unknown restart method $method"))
        end
        extend!(A, L, k)
        isconverged(L, F, l, tol, reltol) && break
    end
    F[:S][1:l], L
end

"""

# References

The simple error bound dates back at least to Wilkinson's classic book
[Wilkinson1965:Ch.3 §53 p.170]

@book{Wilkinson1965,
    address = {Oxford, UK},
    author = {J H Wilkinson},
    publisher = {Oxford},
    title = {The Algebraic Eigenvalue Problem},
    year = 1965
}

The Rayleigh-Ritz bounds were presented in [Wilkinson1965:Ch.3 §54-55 p.173,
    Yamamoto1980, Ortega1990]

@book{Ortega1990,
    address = {Philadelphia, PA},
    author = {Ortega, James M},
    doi = {10.1137/1.9781611971323},
    edition = {2},
    publisher = {SIAM},
    series = {Classics in Applied Mathematics},
    title = {Numerical Analysis: A Second Course},
    url = {http://epubs.siam.org/doi/book/10.1137/1.9781611971323},
    year = 1990
}

@article{Yamamoto1980,
    author = {Yamamoto, Tetsuro},
    doi = {10.1007/BF01396059},
    journal = {Numerische Mathematik},
    number = 2,
    pages = {189--199},
    title = {Error bounds for computed eigenvalues and eigenvectors},
    volume = {34},
    year = 1980
}

"""
function isconverged(L::PartialFactorization,
        F::Base.LinAlg.SVD, k::Int, tol::Real, reltol::Real, verbose::Bool=false)

    @assert tol ≥ 0

    σ = F[:S][1:k]
    Δσ= L.β*abs(F[:U][end, 1:k])

    #Best available eigenvalue bounds
    δσ = copy(Δσ)

    #Update better error bounds from Rayleigh-Ritz
    #Reference: Ortega, 1972
    #In theory these only apply if we know all the values
    #But so long as the gap between the converged Ritz values and all the
    #others is larger than d then we're fine.
    #Can also get some eigenvector statistics!!
    let
        d = Inf
        for i in eachindex(σ), j=1:i-1
            d = min(d, abs(σ[i] - σ[j]))
        end
        verbose && println("Smallest empirical spectral gap: ", d)
        #Chatelein 1993 - normwise backward error associated with approximate invariant subspace
        #Use the largest singular (Ritz) value to estimate the 2-norm of the matrix
        verbose && println("Normwise backward error associated with subspace: ", L.β/σ[1])
        for i in eachindex(Δσ)
            α = Δσ[i]

            #Simple error bound
            #if 2α ≤ d
            verbose && println("Ritz value ", i, ": ", σ[i])
            verbose && println("Simple error bound on eigenvalue: ", Δσ[i])

            if 2α ≤ d
                #Rayleigh-Ritz bounds
                x = α/(d-α)*√(1+(α/(d-α))^2)
                x=abs(x)
                verbose && println("Rayleigh-Ritz error bound on eigenvector: $x")
                #2α ≤ d && (δσ[i] = min(δσ[i], x))

                y = α^2/d #[Wilkinson:Ch.3 Appendix (4), p.188]
                verbose && println("Rayleigh-Ritz error bound on eigenvalue: $y")
                2α ≤ d && (δσ[i] = min(δσ[i], y))
            end

            #Estimate of the normwise backward error [Deif 1989]
            ##Use the largest singular (Ritz) value to estimate the 2-norm of the matrix
            verbose && println("Normwise backward error estimate: ", α/σ[1])

            #Geurts, 1982 - componentwise backward error also known
            #A. J. Geurts, (1982), A contribution to the theory of condition,
            #Numer.  Math., #39, 85-96.

        end
    end

    #Estimate condition number and see if two-sided reorthog is needed
    if verbose && (F[:S][1]/F[:S][end]) > 1/√eps()
        warn("Two-sided reorthogonalization should be used but is not implemented")
    end

    all(δσ[1:k] .< max(tol, reltol*σ[1]))
end

#Hernandez2008
function build{T}(A, q::AbstractVector{T}, k::Int)
    @assert eltype(A) == T
    m, n = size(A)
    Tr = typeof(real(one(T)))

    P = Array(T, m, k)
    Q = Array(T, n, k+1)
    αs= Array(Tr, k)
    βs= Array(Tr, k-1)
    Q[:, 1] = q
    local β
    p = Array(T, m)
    for j=1:k
        #p = A*q
        A_mul_B!(p, A, q)
        if j>1
            βs[j-1] = β
            axpy!(-β, sub(P, :, j-1), p) #p -= β*P[:,j-1]
        end
        α = convert(Tr, norm(p))
        scale!(p, inv(α))

        αs[j] = α
        P[:, j] = p

        Ac_mul_B!(q, A, p) #q = A'p
        axpy!(-1.0, sub(Q, :, 1:j)*(sub(Q, :, 1:j)'q), q) #orthogonalize

        β = convert(Tr, norm(q))
        scale!(q, inv(β))
        Q[:,j+1] = q
    end
    PartialFactorization{T,Tr}(P, Q, Bidiagonal(αs, βs, true), β)
end

"""
Thick restart (with ordinary Ritz values)

# Reference

[Hernandez2008]
"""
function thickrestart!{T,Tr}(A, L::PartialFactorization{T,Tr},
        F::Base.LinAlg.SVD{Tr,Tr}, l::Int)

    k = size(F[:V], 1)
    m, n = size(A)
    #@assert size(L.P) == (m, k)
    #@assert size(L.Q) == (n, k+1)

    Q = sub(L.Q, :,1:k)*sub(F[:V], :,1:l)
    L.Q = [Q sub(L.Q, :, k+1)]
    #Be pedantic about ensuring normalization
    #L.Q = qr(L.Q)[1]
    #@assert all([norm(L.Q[:,i]) ≈ 1 for i=1:size(L.Q,2)])

    f = A*sub(L.Q, :, l+1)
    ρ = L.β * reshape(F[:U][end, 1:l], l)
    L.P = sub(L.P, :, 1:k)*sub(F[:U], :, 1:l)

    #@assert ρ[i] ≈ f⋅L.P[:, i]
    f -= L.P*ρ
    α = convert(Tr, norm(f))
    scale!(f, inv(α))
    L.P = [L.P f]

    g = A'f - α*L.Q[:, end]
    L.β = β = norm(g)
    L.B = BrokenArrowBidiagonal([F[:S][1:l]; α], ρ, typeof(β)[])
    #@assert size(L.P, 2) == size(L.Q, 2) == size(L.B, 2)
    L
end

"""
Thick restart with harmonic Ritz values

# Reference

[Baglama2005] - note that they have P and Q swapped relative to our notation,
which follows that of [Hernandez2008]
"""
function harmonicrestart!{T,Tr}(A, L::PartialFactorization{T,Tr},
        F::Base.LinAlg.SVD{Tr,Tr}, k::Int)

    m = size(L.B, 1)::Int
    #@assert size(L.P,2)==m==size(L.Q,2)-1

    F0 = F# svdfact(L.B)
    ρ = L.β*F0[:U][end,:] #Residuals of singular values

    F2 = svdfact!([diagm(F0[:S]) ρ'], thin=false)

    #Take k largest triplets
    Σ = (F2[:S]::Vector{Tr})[1:k]
    U = F0[:U]*sub(F2[:U],:,1:k)
    M = eye(T, m+1)
    M[1:m, 1:m] = F0[:V]::Matrix{T}
    M = M * F2[:V]
    Mend = M[end, 1:k]
    #Compute scaled residual from the harmonic Ritz problem
    r0 = zeros(Tr, m)
    r0[end] = 1
    isa(L.B, Bidiagonal) && @assert L.B.isupper
    r = try
        #(L.B\r0)
        A_ldiv_B!(L.B, r0)
    catch exc
        if isa(exc, LinAlg.LAPACKException) ||
                isa(exc, LinAlg.SingularException) #B\r is singular

            pinv(full(L.B))*r0
        else rethrow(exc) end
    end::Vector{Tr}
    scale!(r, L.β)
    M::Matrix{T} = sub(M,1:m, :) + r*sub(M,m+1,:)

    M2 = zeros(T, m+1, k+1)
    M2[1:m, 1:k] = M[:,1:k]
    M2[1:m, k+1] = -r
    M2[m+1, k+1] = 1
    Q, R = qr(M2)

    Q = L.Q*Q
    P = L.P*sub(U,:,1:k)

    R = sub(R,1:k+1,1:k) + sub(R,:,k+1)*Mend

    f = A*sub(Q,:,k+1)
    f -= P*(P'f)
    α = convert(Tr, norm(f))
    scale!(f, inv(α))
    P = [P f]
    B = UpperTriangular{Tr,Matrix{Tr}}([(Diagonal(Σ)*triu(R')); zeros(Tr,1,k) α])
    #@assert size(P, 2) == size(B, 1) == size(Q, 2)
    g = A'f
    q = sub(Q,:,k+1)
    #g-= (g⋅q)*q
    axpy!(-(g⋅q), q, g)
    β = convert(Tr, norm(g))
    scale!(g, inv(β))
    #@assert size(P, 2) == size(Q, 2) == size(B, 2)
    L.P = P
    L.Q = Q
    L.B = B
    L.β = β
    L
end

#Hernandez2008
function extend!{T,Tr}(A, L::PartialFactorization{T, Tr}, k::Int)
    l = size(L.B, 2)::Int-1
    p = L.P[:,l+1]

    m, n = size(A)
    q = zeros(T, n)
    #@assert l+1 == size(L.B, 1) == size(L.B, 2) == size(L.P, 2)

    if isa(L.B, UpperTriangular) #Cannot be appended
        B = zeros(Tr, k, k)
        B[1:size(L.B,1), 1:size(L.B,2)] = L.B
        L.B = UpperTriangular(B)
        #@assert size(L.B) == (k, k)
    end

    β = L.β
    for j=l+1:k
        Ac_mul_B!(q, A, p) #q = A'p
        q -= L.Q*(L.Q'q)   #orthogonalize
        β = norm(q)
        scale!(q, inv(β))

        L.Q = [L.Q q]
        j==k && break

        #p = A*q - β*p
        A_mul_B!(p, A, q)
        Base.LinAlg.axpy!(-β, sub(L.P, :, j), p)

        α = norm(p)
        scale!(p, inv(α))
        if isa(L.B, Bidiagonal) || isa(L.B, BrokenArrowBidiagonal)
            push!(L.B.dv, α)
            push!(L.B.ev, β)
        else
            L.B[j+1, j+1] = α
            L.B[j  , j+1] = β
        end
        L.P = [L.P p]
    end
    L.β = β
    L
end

let
    B = BrokenArrowBidiagonal([1, 2, 3], [1, 2], Int[])
    @assert full(B) == [1 0 1; 0 2 2; 0 0 3]
end


text:
# ---
# title: 930. Binary Subarrays With Sum
# id: problem930
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Hash Table, Two Pointers
# link: <https://leetcode.com/problems/binary-subarrays-with-sum/description/>
# hidden: true
# ---
#
# In an array `A` of `0`s and `1`s, how many **non-empty** subarrays have sum
# `S`?
#
#
#
# **Example 1:**
#
#
#
#     Input: A = [1,0,1,0,1], S = 2
#     Output: 4
#     Explanation:
#     The 4 subarrays are bolded below:
#     [ **1,0,1** ,0,1]
#     [ **1,0,1,0** ,1]
#     [1, **0,1,0,1** ]
#     [1,0, **1,0,1** ]
#
#
#
#
# **Note:**
#
#   1. `A.length <= 30000`
#   2. `0 <= S <= A.length`
#   3. `A[i]` is either `0` or `1`.
#
#
## @lc code=start
using LeetCode

function num_subarrays_with_sum(nums::Vector{Int}, goal::Int)
    cnt = OffsetArray(fill(0, length(nums) + 1), -1)
    s, res = 0, 0
    for num in nums
        cnt[s] += 1
        s += num
        res += get(cnt, s - goal, 0)
    end
    res
end
## @lc code=end


text:
##########################
# FunctionCell Interface #
##########################
"""
	kmeans(k, d=Distances.Euclidean() [;kwargs])

Constructs an untrained cell that when piped data inside, clusters the data using the
k-means algorithm where `k` is the desired number of clusters and `d` is the distance
that will be used to calculate distances from new samples to the clustering generated
from the untrained function cell.

# Keyword arguments (same as in `Clustering`)
  * `init` initialization algorithm; can be `:rand` for random, `:kmpp` for K-medoids++ or `:kmcen` for choosing samples with highest centrality
  or a vector of size 'k' providing the indexes of the initial seeds (default `:kmpp`)
  * `maxiter` maximum number of iterations (default `100`)
  * `tol` tolerable change in objective at convergence (default `1.0e-6`)
  * `weights` can be `nothing` for unit weight or a vector of length `n` for individual sample weights (default `nothing`)
  * `display` Level of information; can be `:none`, `:final`, `:iter` (default `:none`)

Read the `Clustering.jl` documentation for more information.
"""
kmeans(k::Int, d::Distances.PreMetric=Distances.Euclidean(); kwargs...) = FunctionCell(kmeans, (k,d), ModelProperties(), kwtitle("K-means ($k clusters)", kwargs); kwargs...)

############################
# DataCell/Array Interface #
############################
"""
	kmeans(x, k, d=Distances.Euclidean() [;kwargs])

Performs k-means clustering on the dataset `x` using `k` number of clusters and returns the clustering; when piping
data to the resulting clustering (e.g. trained cell), the distances from input data to the center clusters are
returned, using `d` as distance.
"""
# Training
kmeans(x::T where T<:CellData, k::Int, d::Distances.PreMetric=Distances.Euclidean(); kwargs...) = kmeans(getx!(x), k, d; kwargs...)
kmeans(x::T where T<:AbstractVector, k::Int, d::Distances.PreMetric=Distances.Euclidean(); kwargs...) = kmeans(mat(x, LearnBase.ObsDim.Constant{2}()), k, d; kwargs...)
kmeans(x::T where T<:AbstractMatrix, k::Int, d::Distances.PreMetric=Distances.Euclidean(); kwargs...) = begin

	# Perform clustering
	kmeansdata = Clustering.kmeans(x, k; kwargs...)

	@assert size(kmeansdata.centers,2) == k "[kmeans] The number of resulting clusters differes from the desired one k=$k."

	# Build model properties
	modelprops = ModelProperties(size(kmeansdata.centers,1),k) # output dimension is 'k' i.e. the number of cluster centers is 'k'

	FunctionCell(kmeans, Model((d,kmeansdata), modelprops), kwtitle("K-means ($k clusters)", kwargs) );
end

# Execution
kmeans(x::T where T<:CellData, model::Model{<:Tuple{<:Distances.PreMetric, <:Clustering.KmeansResult}}) =
	datacell(kmeans(getx!(x), model), gety(x))
kmeans(x::T where T<:AbstractVector, model::Model{<:Tuple{<:Distances.PreMetric, <:Clustering.KmeansResult}}) =
	kmeans(mat(x, LearnBase.ObsDim.Constant{2}()), model)
kmeans(x::T where T<:AbstractMatrix, model::Model{<:Tuple{<:Distances.PreMetric, <:Clustering.KmeansResult}}) =
	# Return distances from each input sample to clustering centers
	Distances.pairwise(model.data[1], model.data[2].centers, x) # d(centers,x)

##########################
# FunctionCell Interface #
##########################
"""
	kmeans!(centers, d=Distances.Euclidean() [;kwargs])

Constructs an untrained cell that when piped data inside, clusters the data using the
k-means algorithm where `centers` contains the initial cluster centers (will be updated
inplace) and `d` is the distance that will be used to calculate distances from new samples
to the clustering generated from the untrained function cell.

# Keyword arguments (same as in `Clustering`)
  * `maxiter` maximum number of iterations (default `100`)
  * `tol` tolerable change in objective at convergence (default `1.0e-6`)
  * `weights` can be `nothing` for unit weight or a vector of length `n` for individual sample weights (default `nothing`)
  * `display` Level of information; can be `:none`, `:final`, `:iter` (default `:none`)

Read the `Clustering.jl` documentation for more information.
"""
kmeans!(centers::T where T<:AbstractMatrix, d::Distances.PreMetric=Distances.Euclidean(); kwargs...) =
	FunctionCell(kmeans!, (centers,d), ModelProperties(), kwtitle("K-means! ($(size(centers,2)) clusters)", kwargs); kwargs...)

############################
# DataCell/Array Interface #
############################
"""
	kmeans!(x, centers, d=Distances.Euclidean() [;kwargs])

Performs k-means clustering on the dataset `x` using `centers` as initial cluster centers
and returns the clustering; when piping data to the resulting clustering (e.g. trained cell),
the distances from input data to the new cluster centers are returned, using `d` as distance.
"""
# Training
kmeans!(x::T where T<:CellData, centers::S where S<:AbstractMatrix, d::Distances.PreMetric=Distances.Euclidean(); kwargs...) = kmeans!(getx!(x), centers, d; kwargs...)
kmeans!(x::T where T<:AbstractVector, centers::S where S<:AbstractMatrix, d::Distances.PreMetric=Distances.Euclidean(); kwargs...) = kmeans!(mat(x, LearnBase.ObsDim.Constant{2}()), centers, d; kwargs...)
kmeans!(x::T where T<:AbstractMatrix, centers::S where S<:AbstractMatrix, d::Distances.PreMetric=Distances.Euclidean(); kwargs...) = begin

	# Perform clustering
	kmeansdata = Clustering.kmeans!(x, centers; kwargs...)

	# Build model properties
	modelprops = ModelProperties(size(kmeansdata.centers,1), size(kmeansdata.centers,2), nothing, Dict("eval_distance" => d))

	# Return a trained cell that uses `kmeans` and not `kmeans!` as execution function
	# since the only difference is in training (e.g. clustering) and not in
	# assigning new data to clusters.
	FunctionCell(kmeans, Model((d,kmeansdata), modelprops), kwtitle("K-means ($(size(kmeansdata.centers,2)) clusters)", kwargs) );
end


text:
using PascalsTriangle
using BenchmarkTools

function binomialrow(n)
    return binomial.(n, 0:n)
end

function ptrow(n)
    return Row(n)
end

function bin_all_rows(n)
    a = Vector(undef, n+1)
    for i ∈ 0:n
        a[1:i+1] = binomial.(i,0:i)
    end
end

function pt_all_rows(n)
    r = Row(Integer,0,n)
    for i ∈ 1:n
        next!(r)
    end
end

rnums = rand(0:50, 100)

@btime binomialrow.(rnums)
@btime ptrow.(rnums)

@btime bin_all_rows(50)
@btime pt_all_rows(50)


text:
if !isdefined(Main, :Groebner)
    import Groebner
end

import AbstractAlgebra
using BenchmarkTools
using Logging
global_logger(ConsoleLogger(stderr, Logging.Error))

BenchmarkTools.DEFAULT_PARAMETERS.seconds = 10
# BenchmarkTools.DEFAULT_PARAMETERS.samples = 2

function benchmark_system_my(system)
    system = Groebner.change_ordering(system, :degrevlex)
    Groebner.groebner([system[1]])

    # println("length = $(length(system)), nvars = $(AbstractAlgebra.nvars(AbstractAlgebra.parent(system[1])))")
    # println("maxdeg = $(maximum(AbstractAlgebra.total_degree ∘ AbstractAlgebra.leading_monomial, system))")

    @btime Groebner.groebner($system, reduced=true)
end

function run_f4_ff_degrevlex_benchmarks(ground)
    systems = [
    ("cyclic 10", Groebner.rootn(10, ground=ground)),
    ("cyclic 11", Groebner.rootn(11, ground=ground)),
    ("cyclic 12", Groebner.rootn(12, ground=ground)),
    ("cyclic 13", Groebner.rootn(13, ground=ground)),
        ("katsura 6", Groebner.katsura6(ground=ground)),
        ("katsura 9", Groebner.katsura9(ground=ground)),
        ("katsura 10", Groebner.katsura10(ground=ground)),
        ("noon 4"    ,Groebner.noonn(4, ground=ground)),
        ("noon 5"    ,Groebner.noonn(5, ground=ground)),
        ("noon 6"    ,Groebner.noonn(6, ground=ground)),
        ("noon 7"    ,Groebner.noonn(7, ground=ground)),
        ("eco 7"    ,Groebner.eco7(ground=ground)),
        ("eco 10"    ,Groebner.eco10(ground=ground)),
        ("eco 11"    ,Groebner.eco11(ground=ground)),
        ("ku 10"    ,Groebner.ku10(ground=ground)),
        ("s9_1"    ,Groebner.s9_1(ground=ground)),
        ("kinema"    ,Groebner.kinema(ground=ground)),
        ("ojika4_d1R2_d2R5"    ,Groebner.ojika4_d1R2_d2R5(ground=ground))
    ]

    println()

    for (name, system) in systems
        println("$name")
        benchmark_system_my(system)
    end
end

function big_run_f4_ff_degrevlex_benchmarks(ground)
    systems = [
        #("cyclic 14", Groebner.rootn(14, ground=ground)),
        #("cyclic 15", Groebner.rootn(15, ground=ground)),
        #("katsura 11", Groebner.katsura11(ground=ground)),
        #("katsura 12", Groebner.katsura12(ground=ground)),
        #("katsura 13", Groebner.katsura13(ground=ground)),
        ("noon 7"    ,Groebner.noonn(7, ground=ground)),
        ("noon 8"    ,Groebner.noonn(8, ground=ground)),
        ("eco 11"    ,Groebner.eco11(ground=ground)),
        ("eco 12"    ,Groebner.eco12(ground=ground)),
        #("eco 13"    ,Groebner.eco13(ground=ground)),
    ]

    println()

    for (name, system) in systems
        println("$name")
        benchmark_system_my(system, tablesize)
    end
end

println()
ground = AbstractAlgebra.GF(2^31 - 1)
run_f4_ff_degrevlex_benchmarks(ground)

#=
2^6:
cyclic 10
  25.030 ms (61954 allocations: 10.70 MiB)
cyclic 11
  84.327 ms (156817 allocations: 25.69 MiB)
cyclic 12
  352.658 ms (401131 allocations: 70.27 MiB)
cyclic 13
  1.492 s (1032877 allocations: 186.29 MiB)
katsura 6
  3.615 ms (8778 allocations: 2.07 MiB)
katsura 9
  112.370 ms (56824 allocations: 18.90 MiB)
katsura 10
  765.672 ms (152269 allocations: 62.73 MiB)
noon 4
  1.133 ms (4518 allocations: 831.12 KiB)
noon 5
  5.000 ms (19465 allocations: 3.80 MiB)
noon 6
  36.483 ms (89537 allocations: 20.19 MiB)
noon 7
  253.887 ms (426593 allocations: 93.01 MiB)
eco 7
  2.122 ms (8617 allocations: 2.13 MiB)
eco 10
  196.016 ms (139553 allocations: 42.71 MiB)
eco 11
  1.302 s (394404 allocations: 136.15 MiB)
ku 10
  1.146 ms (5884 allocations: 2.53 MiB)
s9_1
  495.600 μs (1901 allocations: 1.82 MiB)
kinema
  2.961 ms (7192 allocations: 3.36 MiB)
ojika4_d1R2_d2R5
  311.500 μs (2154 allocations: 354.58 KiB)

2^4:

2^7:

=#


text:
@testset "Outer face" begin
    tess = Tessellation()

    u, v, w = Vertex(1, 0), Vertex(0, 1), Vertex(-1, 0)
    a1, a2, b1, b2, c1, c2 = [Edge() for _ in 1:6]
    f1, f2 = Face(), Face()

    push!(tess.vertices, u, v, w)
    push!(tess.edges, a1, a2, b1, b2, c1, c2)
    push!(tess.faces, f1, f2)

    u.edge, v.edge, w.edge = a1, b1, c1

    a1.dest, b1.dest, c1.dest = u, v, w
    a2.dest, b2.dest, c2.dest = w, u, v

    a1.twin, a2.twin = a2, a1
    b1.twin, b2.twin = b2, b1
    c1.twin, c2.twin = c2, c1

    a1.next, b1.next, c1.next = b1, c1, a1
    a1.prev, b1.prev, c1.prev = c1, a1, b1

    a2.next, b2.next, c2.next = c2, a2, b2
    a2.prev, b2.prev, c2.prev = b2, c2, a2

    a1.face = b1.face = c1.face = f1
    a2.face = b2.face = c2.face = f2

    f1.edge, f2.edge = a1, a2
    f = find_outer_face(tess)
    @test f == f2
end

@testset "Fix faces" begin
    @testset "Error" begin
        tess = Tessellation()
        f = Face()
        push!(tess.faces, f)
        @test_throws GeometryError fix_faces!(tess)
    end

    @testset "Results" begin
        tess = Tessellation()

        u, v, w = Vertex(1, 0), Vertex(0, 1), Vertex(-1, 0)
        a1, a2, b1, b2, c1, c2 = [Edge() for _ in 1:6]

        push!(tess.vertices, u, v, w)
        push!(tess.edges, a1, a2, b1, b2, c1, c2)

        u.edge, v.edge, w.edge = a1, b1, c1

        a1.dest, b1.dest, c1.dest = u, v, w
        a2.dest, b2.dest, c2.dest = w, u, v

        a1.twin, a2.twin = a2, a1
        b1.twin, b2.twin = b2, b1
        c1.twin, c2.twin = c2, c1

        a1.next, b1.next, c1.next = b1, c1, a1
        a1.prev, b1.prev, c1.prev = c1, a1, b1

        a2.next, b2.next, c2.next = c2, a2, b2
        a2.prev, b2.prev, c2.prev = b2, c2, a2

        fix_faces!(tess)
        @test a1.face == b1.face == c1.face
        @test a2.face == b2.face == c2.face
    end
end


text:
"""
$(SIGNATURES)

Calculates the code phase error in chips.
"""
function dll_disc(
    ::Type{S},
    correlator,
    early_late_sample_shift,
    code_phase_delta
) where S <: AbstractGNSSSystem
    E = abs(get_early(correlator))
    L = abs(get_late(correlator))
    distance_between_early_and_late = 2 * early_late_sample_shift * code_phase_delta
    (E - L) / (E + L) / (2 * (2 - distance_between_early_and_late))
end

"""
$(SIGNATURES)

Calculates the carrier phase error in radians.
"""
function pll_disc(::Type{S}, correlator) where S <: AbstractGNSSSystem
    p = get_prompt(correlator)
    atan(imag(p) / real(p))
end


text:
using ExactD, LinearAlgebra, Test

include("./test_basis_states.jl")
include("./test_expectation_values.jl")
include("./test_many_body_operators.jl")
include("./test_exponential_ops.jl")


text:
using Flux
using Flux: params
using Flux.Data: DataLoader

using Random
rng = MersenneTwister(2345)

@testset "Computer runs without error" for B in 1:2
    R = 1
    N = 5
    W = 10
    X = 2
    Y = 2
    controut = 10
    inputsize = DNC.inputsize(X, R, W)
    x = rand(Float32, X, B)
    y = rand(Float32, Y, B)
    controller = Chain(Dense(inputsize, R*W*2, relu),
                       Dense(R*W*2, controut, relu))
    dnc1 = Dnc(controller, X, Y, controut, N, W, R, B)
    out1 = dnc1(x)
    @test size(out1) == (Y, B)
    @test eltype(out1) == Float32
    dnc2 = Dnc(X, Y, controut, N, W, R, B)
    out2 = dnc2(x)
    @test size(out2) == (Y, B)
    @test eltype(out2) == Float32
    # Test clip value
    clipvalue = 0.1
    dnc3 = Dnc(X, Y, controut, N, W, R, B; clipvalue=clipvalue)
    @test all(dnc3(x) .<= clipvalue)
    @test all(dnc3(x) .>= -clipvalue)
end

R, N, W, X, Y, B = 2, 10, 16, 2, 2, 4
controllerout = 64
dnc = Dnc(X, Y, controllerout, N, W, R, B)
@show(dnc)

function generatedata(in, out, n=100)
    w = rand(rng, Float32, out, in)
    b = rand(rng, Float32, out)
    Xs = rand(rng, Float32, in, n)
    Ys = w*Xs .+ b
    Xs, Ys
end

xs, ys = generatedata(X, Y)
batcheddata = DataLoader(xs, ys, batchsize=B)
loss(x, y) = Flux.mse(dnc(x), y)
opt = ADAM(0.01)
onebatch = Base.iterate(batcheddata, 0)[1]

@testset "Learn" begin
    Flux.train!(loss, params(dnc), batcheddata, opt)
    @test eltype(dnc(onebatch[1])) == Float32
    @test eltype(dnc.cell.memoryaccess.M) == Float32
    @test loss(onebatch...) < 1.0*B
end


text:
using AbstractDomains
using Base.Test
import AbstractDomains: isintersect, subsumes, ⊔
import AbstractDomains: t, f, tf

@test t & f === f
@test tf & f === f
@test t & t === t
@test tf & tf === tf
@test !t === f
@test !tf === tf
@test !f === t
@test tf | t === t
@test tf | f === tf
@test f | t === t
@test t | tf === t

# Lifted equality tests
@test (t == t) === t
@test (t == tf) === tf
@test (tf == t) === tf
@test (f == f) === t
@test (t == f) === f
@test (f == t) === f

# isintersect
@test isintersect(t,f) == false
@test isintersect(f,t) == false
@test isintersect(t,t) == true
@test isintersect(f,f) == true
@test isintersect(tf,t) == true
@test isintersect(t,tf) == true
@test isintersect(tf,f) == true
@test isintersect(f,tf) == true
@test isintersect(tf,tf) == true

# Subsumes
@test subsumes(tf,f) == true
@test subsumes(tf,t) == true
@test subsumes(tf,tf) == true
@test subsumes(f,f) == subsumes(t,t) == true
@test subsumes(f,t) == subsumes(t,f) == false
@test subsumes(f,tf) == subsumes(t,tf) == false

# Join
@test ⊔(t,f) === ⊔(f,t) === tf
@test ⊔(t,t) === t
@test ⊔(f,f) === f
@test ⊔(tf,t) === ⊔(t,tf) === tf
@test ⊔(tf,f) === ⊔(f,tf) === tf
@test ⊔(tf,tf) === tf

#ifelse
@test ifelse(tf,f,t) === tf
@test ifelse(t,f,t) === f
@test ifelse(f,f,t) === t


text:
# A cache type for GenericMERA.
# To be included in MERAKit.jl.

# The storage formats for density matrices, operators, and environments are a little
# different:
#
# For density matrices we always store a `Vector` of the same length as there are
# layers. If a density matrix is not in store, `nothing` is kept as the corresponding
# element.
#
# `operators` is a dictionary with operators as keys, and each element is a
# `Vector` that lists the ascended versions of that operator, starting from the physical
# one, which is nothing but the original operator, possibly wrapped in some container type.
# No place-holder `nothing`s are stored, the `Vector` ends when storage ends.
#
# `environments` is a `Dict` just like `operators`, but its values are `Vectors` that hold
# either `Layer`s or `Nothing`s, where each `Layer` has instead of tensors of the layer,
# environments for the tensors of a layer.
#
# To understand the reason for these differences, note that density matrices are naturally
# generated from the top, and there will never be more than the number of layers of them.
# Similarly it doesn't make sense to have more than one environment for each layer.
# Operators are generated from the bottom, and they may go arbitrarily high in the MERA
# (into the scale invariant part).
#
# In addition we keep `previous_fixedpoint_densitymatrix` and `previous_operatorsum`, to use
# as initial guesses when solving for the fixed point density matrix or scale invariant
# operator sum.

"""
    MERACache{N, LT <: Layer, OT}

A cache for `GenericMERA`. The type parameters are the same as for `GenericMERA`.

See also: [`GenericMERA`](@ref)
"""
mutable struct MERACache{N, LT <: Layer, OT}
    # LT stands for Layer Type, OT for Operator Type.
    # TODO Should we use NTuple{N}s instead for densitymatrices and environments?
    densitymatrices::Vector{Union{Nothing, OT}}
    # probably a better choice to use IdDicts here
    operators::Dict{Any, Vector{OT}}
    environments::Dict{Any, Vector{Union{Nothing, LT}}}
    previous_fixedpoint_densitymatrix::Union{Nothing, OT}
    previous_operatorsum::Union{Nothing, OT}

    function MERACache{N, LT, OT}() where {N, LT, OT}
        @assert OT === operatortype(LT)
        densitymatrices = Vector{Union{Nothing, OT}}(fill(nothing, N))
        operators = Dict{Any, Vector{OT}}()
        environments = Dict{Any, Vector{Union{LT}}}()
        previous_fixedpoint_densitymatrix = nothing
        previous_operatorsum = nothing
        new{N, LT, OT}(
            densitymatrices,
            operators,
            environments,
            previous_fixedpoint_densitymatrix,
            previous_operatorsum
        )
    end
end

function MERACache{N, LT}() where {N, LT}
    OT = operatortype(LT)
    return MERACache{N, LT, OT}()
end

operatortype(::Type{MERACache{N, LT, OT}}) where {N, LT, OT} = OT
layertype(::Type{MERACache{N, LT, OT}}) where {N, LT, OT} = LT
baselayertype(::Type{MERACache{N, LT, OT}}) where {N, LT, OT} = baselayertype(LT)
causal_cone_width(C::Type{<:MERACache}) = causal_cone_width(layertype(C))

function Base.copy!(dst::MERACache, src::MERACache)
    dst.densitymatrices = copy(src.densitymatrices)
    dst.operators = Dict()
    dst.environments = Dict()
    for (k, v) in src.operators
        dst.operators[k] = copy(v)
    end
    for (k, v) in src.environments
        dst.environments[k] = copy(v)
    end
    dst.previous_operatorsum = src.previous_operatorsum
    dst.previous_fixedpoint_densitymatrix = src.previous_fixedpoint_densitymatrix
    return dst
end

Base.copy(c::MERACache) = copy!(typeof(c)(), c)

"""
    replace_layer(c::MERACache, depth)

Create a new `MERACache` that has all the stored pieces removed that are invalidated by
changing the layer at `depth`.
"""
replace_layer(c::MERACache, depth) = replace_layer!(copy(c), depth)

function replace_layer!(c::MERACache{N}, depth) where N
    c = copy(c)
    depth = Int(min(N, depth))
    c.densitymatrices[1:depth] .= nothing
    for v in values(c.operators)
        last_index = min(depth, length(v))
        deleteat!(v, last_index+1:length(v))
    end
    # Changing anything always invalidates all environments, since they depend on things
    # both above and below.
    reset_environment_storage!(c)
    return c
end

function release_transitionlayer(c::MERACache{N, LT}) where {N, LT}
    new_c = MERACache{N+1, LT}()
    copy!(new_c, c)
    density_matrix = new_c.densitymatrices[end]
    push!(new_c.densitymatrices, density_matrix)
    for (k, v) in new_c.environments
        push!(v, nothing)
    end
    return new_c
end

"""
    has_densitymatrix_stored(c::MERACache, depth)

Return whether the density matrix at given depth is already in the cache.

See also: [`set_stored_densitymatrix!`](@ref), [`get_stored_densitymatrix`](@ref)
"""
function has_densitymatrix_stored(c::MERACache{N}, depth) where N
    depth = Int(min(N, depth))
    rho = c.densitymatrices[depth]
    return rho !== nothing
end

"""
    get_stored_densitymatrix(c::MERACache, depth)

Return the density matrix at given depth, assuming it's in the cache.

See also: [`set_stored_densitymatrix!`](@ref), [`has_densitymatrix_stored`](@ref)
"""
function get_stored_densitymatrix(c::MERACache{N}, depth) where N
    depth = Int(min(N, depth))
    rho = c.densitymatrices[depth]
    @assert rho !== nothing
    return rho
end

"""
    set_stored_densitymatrix!(c::MERACache, density_matrix, depth)

Cache the density matrix at given depth.

See also: [`get_stored_densitymatrix`](@ref), [`has_densitymatrix_stored`](@ref)
"""
function set_stored_densitymatrix!(c::MERACache, density_matrix, depth)
    c.densitymatrices[depth] = density_matrix
    return c
end

"""
    operator_storage!(c::MERACache, op)

Return the cached ascended versions of `op`. Initialize the cache for `op` if necessary.
"""
function operator_storage!(c::MERACache{N, LT, OT}, op) where {N, LT, OT}
    if !haskey(c.operators, op)
        op_conv = convert(OT, expand_support(op, causal_cone_width(c)))
        c.operators[op] = [op_conv]
    end
    return c.operators[op]
end

"""
    has_operator_stored(c::MERACache, op, depth)

Return whether the operator `op` ascended to `depth` is already in the cache.

See also: [`set_stored_operator!`](@ref), [`get_stored_operator`](@ref)
"""
function has_operator_stored(c::MERACache, op, depth)
    storage = operator_storage!(c, op)
    return length(storage) >= depth
end

"""
    get_stored_operator(c::MERACache, op, depth)

Return the operator `op` ascended to a given depth, assuming it's in the cache.

See also: [`set_stored_operator!`](@ref), [`has_operator_stored`](@ref)
"""
function get_stored_operator(c::MERACache, op, depth)
    storage = operator_storage!(c, op)
    return storage[depth]
end

"""
    set_stored_operator!(c::MERACache, opasc, op, depth)

Cache `opasc`, the operator `op` ascended to a given depth.

See also: [`get_stored_operator`](@ref), [`has_operator_stored`](@ref)
"""
function set_stored_operator!(c::MERACache, opasc, op, depth)
    storage = operator_storage!(c, op)
    if length(storage) < depth-1
        msg = "Can't store an ascended operator if the lower versions of it aren't in storage already."
        throw(ArgumentError(msg))
    elseif length(storage) == depth-1
        push!(storage, opasc)
    else
        storage[depth] = opasc
    end
    return c
end

"""
    reset_operator_storage!(c::MERACache, op)

Reset cache for a given operator, removing its ascended memory and forcing recomputation.
"""
function reset_operator_storage!(c::MERACache, op)
    delete!(c.operators, op)
    return c
end

"""
    reset_operator_storage!(c::MERACache)

Reset storage for all operators.
"""
function reset_operator_storage!(c::MERACache)
    empty!(c.operators)
    return c
end

"""
    environment_storage!(c::MERACache, op)

Return the environments related to `op`. Initialize the storage `Vector` if necessary.
"""
function environment_storage!(c::MERACache{N}, op) where N
    if !haskey(c.environments, op)
        c.environments[op] = repeat(Union{Nothing, operatortype(c)}[nothing], N)
    end
    return c.environments[op]
end

"""
    has_environment_stored(c::MERACache, op, depth)

Return whether the environment related to `op` at `depth` is already in store.

See also: [`set_stored_environment!`](@ref), [`get_stored_environment`](@ref)
"""
function has_environment_stored(c::MERACache, op, depth)
    storage = environment_storage!(c, op)
    return storage[depth] !== nothing
end

"""
    get_stored_environment(c::MERACache, op, depth)

Return the environment related to `op` at `depth`. A type assertion error is thrown if the
requested environment is not in storage.

See also: [`set_stored_environment!`](@ref), [`has_environment_stored`](@ref)
"""
function get_stored_environment(c::MERACache, op, depth)
    storage = environment_storage!(c, op)
    env = storage[depth]
    @assert env !== nothing
    return env
end

"""
    set_stored_environment!(c::MERACache, env, op, depth)

Store `env`, the environments related to `op` at `depth`.

See also: [`get_stored_environment`](@ref), [`has_environment_stored`](@ref)
"""
function set_stored_environment!(c::MERACache, env, op, depth)
    storage = environment_storage!(c, op)
    storage[depth] = env
    return c
end

"""
    reset_environment_storage!(c::MERACache)

Reset storage for environments.
"""
function reset_environment_storage!(c::MERACache)
    empty!(c.environments)
    return c
end

function expand_bonddim!(c::MERACache{N, LT}, depth, V) where {N, LT}
    # The replacement of layers has already caused the removal of of all necessary stored
    # operators. This only needs to take care of the stuff related to the scale invariant
    # layer.
    depth < N - 1 && return c
    width = causal_cone_width(LT)
    # Pad the stored scale invariant initial guesses.
    old_rho = c.previous_fixedpoint_densitymatrix
    if old_rho !== nothing
        for i in 1:width
            old_rho = pad_with_zeros_to(old_rho, i => V, (i+width) => V')
        end
        c.previous_fixedpoint_densitymatrix = old_rho
    end
    old_opsum = c.previous_operatorsum
    if old_opsum !== nothing
        for i in 1:width
            old_opsum = pad_with_zeros_to(old_opsum, i => V, (i+width) => V')
        end
        c.previous_operatorsum = old_opsum
    end
    return c
end