text:
## Code that is Gtk specific

## We have two functions for a given type:
## * gtk_widget(widget::T) to create a Gtk widget for the type T
##
## and to make an output widget
##
## * GtkInteract.show_outwidget(w::GtkInteract.MainWindow, x::T)
##
## For example, this creates a custom label:
##
## function GtkInteract.show_outwidget(w::GtkInteract.MainWindow, x::T)
##     if w.out == nothing
##         w.out = label()
##         push!(w, grow(w.out))
##         display(w)
##     end
##     push!(w.out, x.x)
## end
##
## The key is setting w.out.

## Plotting code is package dependent
## Plots should dispatch to underlying code for Immerse, Winston or PyPlot

## Plots package

## Helper to add canvas to main window if none there
function make_canvas(w::MainWindow, x)
    if w.out == nothing
        w.out = cairographic()
        push!(w, grow(w.out))
        display(w)
    end
    push!(w.out, x)
end

Requires.@require Plots begin
    ## In general use a canvas to display a plot in the main window
    show_outwidget(w::GtkInteract.MainWindow, x::Plots.Plot) = make_canvas(w, x)

    function Base.push!(obj::CairoGraphic, p::Plots.Plot)
        if obj.obj != nothing
            push!(obj.obj, p)
        end
    end
    Base.push!(c::GtkCanvas, p::Plots.Plot) = push!(c, p.o[2])

    ## show_outwidget makes a display widget for a plot

    ## ## For unicode plots we try a label, but doesn't work...
    ## function show_outwidget(w::GtkInteract.MainWindow, p::Plots.Plot{Plots.UnicodePlotsPackage})
    ##      .....
    ##     ## Fails!
    ##     Plots.rebuildUnicodePlot!(p)
    ##     out = sprint(io -> show(io, "text/plain", p.o))
    ##     setproperty!(w.label, :label, out)
    ## end
end

## Immerse
## XXX This has issues, as the canvas doen't get refreshed between draws
Requires.@require Immerse begin
    eval(Expr(:using, :Gadfly))
    eval(Expr(:using, :Compose))

    function gtk_widget(widget::ImmerseFigure)
        widget.obj != nothing && return(widget.obj)

        box, toolbar, cnv = Immerse.createPlotGuiComponents()
        Gtk.G_.size_request(box, 480, 400)
        widget.obj = box
        widget.toolbar = toolbar
        widget.cnv = cnv

        ## Add figure, But should I close the next one??
        i = Immerse.nextfig(Immerse._display)
        f = Immerse.Figure(cnv)
        Immerse.initialize_toolbar_callbacks(f)
        Immerse.addfig(Immerse._display, i, f)

        box
    end

    ## same as Plots.Plot, make canvas window
    function show_outwidget(w::GtkInteract.MainWindow, x::Gadfly.Plot)
        if w.out == nothing
            w.out = immersefigure()
            push!(w, grow(w.out))
            display(w)
        end
        display(Immerse._display, x)
    end

    ## this is used by Plots+Immerse
    function Base.push!(c::GtkCanvas, p::Gadfly.Plot)
        c.back = Gtk.cairo_surface_for(c)
        display(c, Immerse.Figure(c, p))
    end

    function Base.push!(obj::CairoGraphic, p::Gadfly.Plot)
        if obj.obj != nothing
            push!(obj.obj, p)
        end
    end
end

## Winston
## Winston is a problem if `ENV["WINSTON_OUTPUT"] = :gtk` is not set *beforehand*
Requires.@require Winston begin
    ENV["WINSTON_OUTPUT"] = :gtk

    show_outwidget(w::GtkInteract.MainWindow, x::Winston.PlotContainer) = make_canvas(w, x)

    function Base.push!(obj::CairoGraphic, pc::Winston.PlotContainer)
        if obj.obj != nothing
            Winston.display(obj.obj, pc)
        end
    end
    Base.push!(c::GtkCanvas, pc::Winston.PlotContainer) = Winston.display(c, pc)
end

## Gadfly
Requires.@require Gadfly begin
    info("Gadfly support is through the Immerse package.")
end

Requires.@require PyPlot begin
    info("PyPlot support is very buggy!")

    PyPlot.pygui(false)

"""

Overwrite `withfig` from `PyPlot` to work here. This code snippet is
found in PyPlot, save the comment. It is licensed under the MIT
license.

We use this as with `Interact`:

```
f = figure()
@manipulate for n in 1:5
    GtkInteract.withfig(f) do
        xs = linspace(0, n*pi)
        PyPlot.plot(xs, map(sin, xs))
    end
end
```
"""
function withfig(actions::Function, f::PyPlot.Figure; clear=true)
    ax_save = PyPlot.gca()
    PyPlot.figure(f[:number])
    finalizer(f, close)
    try
        if clear && !isempty(f)
            PyPlot.clf()
        end
        actions()
    catch
        rethrow()
    finally
        try
            PyPlot.sca(ax_save) # may fail if axes were overwritten
        end
        ##Main.IJulia.undisplay(f) ## IJulia display queue
    end
    return f
end
export withfig

    " How to show a pyplot figure "
    function show_outwidget(w::GtkInteract.MainWindow, x::PyPlot.Figure)
        error("XXX This is broken XXX")
        if w.out == nothing
            w.out = image()
            push(w, grow(w.out))
            display(w)
        end

        f = tempname() * ".png"
        io = open(f, "w")
        show(io, "image/png", x)
        close(io)
        push!(w.out, f)
        rm(f)
        x[:clear]()
        nothing
    end
end

##################################################
## Controls

## link up widget and its gtk proxy

## button
##
## button("label") is constructor
##

button_cb = function(btnpr::Ptr, user_data)
    w, o = user_data
    push!(w.signal, nothing)
    nothing
end

function gtk_widget(widget::Button)
    obj = @GtkButton()
    Gtk.G_.label(obj, widget.label)

    ## widget -> signal
    id = signal_connect(button_cb, obj, "clicked", Void, (),  false, (widget, obj))

#    id = signal_connect(obj, :clicked) do obj, args...
#        push!(widget, nothing)
#    end
    obj
end

## checkbox
function checkbox_cb(cbptr::Ptr,  user_data)
    w, o = user_data
    val = getproperty(o, :active, Bool)
    push!(w.signal, getproperty(o, :active, Bool))
    nothing
end
function gtk_widget(widget::Checkbox)
    obj = @GtkCheckButton()
    setproperty!(obj, :active, Reactive.value(Interact.signal(widget)))

    ## widget -> signal
    id = signal_connect(checkbox_cb, obj, "toggled", Void, (),  false, (widget, obj))

   Reactive.foreach(widget.signal) do val
       signal_handler_block(obj, id)
       curval = getproperty(obj, :active, Bool)
       val != curval && setproperty!(obj, :active, val)
       signal_handler_unblock(obj, id)
   end

    obj
end

## slider
## following
## https://github.com/timholy/GtkUtilities.jl/blob/b36edfca6b4f0c0b6351f1c1409e4e2d04ca4f8f/src/link.jl#L243-L253
## strangely, this can not be defined within gtk_widget...
function scale_cb(scaleptr::Ptr, user_data)
    w, o = user_data
    val = Gtk.G_.value(o)
    push!(w.signal, val)
    nothing
end

function gtk_widget(widget::Slider)
    obj = @GtkScale(false, first(widget.range), last(widget.range), step(widget.range))
    Gtk.G_.size_request(obj, 200, -1)
    Gtk.G_.value(obj, widget.value)

    ## widget -> signal
    ## https://github.com/JuliaLang/Gtk.jl/blob/master/doc/more_signals.md
    id = signal_connect(scale_cb, obj, "value-changed", Void, (), false, (widget, obj))

    ## ## this might be an issue with #161
    ## id = signal_connect(obj, :value_changed) do obj, args...
    ##     val = Gtk.G_.value(obj)
    ##     push!(widget.signal, val)
    ## end

    ##
    Reactive.foreach(widget.signal) do val
        signal_handler_block(obj, id)
        curval = Gtk.G_.value(obj)
        curval != val && Gtk.G_.value(obj, val)
        signal_handler_unblock(obj, id)
    end

    obj
end

## togglebutton (single one. XXX is label on button or a label? XXX)
function togglebtn_cb(toggleptr::Ptr, user_data)
    w, o = user_data
    value = getproperty(o, :active, Bool)
    push!(w.signal, value)
    nothing
end

function gtk_widget(widget::ToggleButton)
    obj = @GtkToggleButton(string(widget.label))
    setproperty!(obj, :active, widget.value)

    ## widget -> signal
    id = signal_connect(togglebtn_cb, obj, "toggled", Void, (), false, (widget, obj))

    ## signal -> widget
    Reactive.foreach(widget.signal) do val
        signal_handler_block(obj, id)
        setproperty!(obj, :active, val)
        signal_handler_unblock(obj, id)
    end

    obj
end

## textbox
function gtk_widget(widget::Textbox)
    obj = @GtkEntry
    setproperty!(obj, :text, string(widget.signal.value))

    id = signal_connect(obj, :activate) do o
        try
            val = Interact.parse_msg(widget, getproperty(o, :text, AbstractString))
            push!(widget.signal, val)
        catch
            str = getproperty(o, :text, AbstractString)
            messagebox("\"$str\" cannot be parsed as a $(eltype(signal(widget)))", :error)
        end
        false
    end

    ## signal -> widget
    Reactive.foreach(widget.signal) do val
        signal_handler_block(obj, id)
        setproperty!(obj, :text, string(val))
        signal_handler_unblock(obj, id)
    end

    obj
end

## dropdown
function combobox_cb(o::Ptr, user_data)
    w, o = user_data
    index = getproperty(o, :active, Int) + 1
    push!(w.signal, collect(values(w.options))[index])
    nothing
end

function gtk_widget(widget::Options{:Dropdown})
    obj = @GtkComboBoxText(false)
    for key in keys(widget.options)
        push!(obj, key)
    end
    index = findfirst(collect(keys(widget.options)), widget.value_label)
    setproperty!(obj, :active, index - 1)

    ## widget -> signal
    id = signal_connect(combobox_cb, obj, "changed", Void, (), false, (widget, obj))
#    id = signal_connect(obj, :changed) do obj, args...
#        index = getproperty(obj, :active, Int) + 1
#        push!(widget.signal, collect(values(widget.options))[index])
#    end

    ## signal -> widget
    Reactive.foreach(widget.signal) do val
        signal_handler_block(obj, id)
        index = getproperty(obj, :active, Int) + 1
        val = findfirst(collect(values(widget.options)), val)
        if val != index
            setproperty!(obj, :active, val - 1)
        end
        signal_handler_unblock(obj, id)
    end

    obj

end

## radiobuttons
function radiobtn_cb(r::Ptr, user_data)
    w, o = user_data
    if getproperty(o, :active, Bool)
        label = getproperty(o, :label, AbstractString)
        push!(w.signal, w.options[label])
    end
    nothing
end
function gtk_widget(widget::Options{:RadioButtons})
    obj = @GtkBox(false)
    choices = collect(keys(widget.options))
    btns = [@GtkRadioButton(shift!(choices))]
    while length(choices) > 0
        push!(btns, @GtkRadioButton(btns[1], shift!(choices)))
    end
    map(u->push!(obj, u), btns)

    selected = findfirst(collect(values(widget.options)), widget.value)
    setproperty!(btns[selected], :active, true)

    ## widget -> signal
    ids = Dict()
    for btn in btns
        ids[btn] = signal_connect(radiobtn_cb, btn, "toggled", Void, (), false, (widget, btn))
        ## ids[btn] = signal_connect(btn, :toggled) do obj, args...
        ##     if getproperty(obj, :active, Bool)
        ##         label = getproperty(obj, :label, AbstractString)
        ##         push!(widget.signal, widget.options[label])
        ##     end
        ## end
    end
    setproperty!(obj, :visible, true)
    showall(obj)

    ## signal -> widget
    Reactive.foreach(widget.signal) do val
        [signal_handler_block(btn, id) for (btn,id) in ids]
        selected = findfirst(collect(values(widget.options)), val)
        setproperty!(btns[selected], :active, true)
        [signal_handler_unblock(btn, id) for (btn,id) in ids]
    end

    obj

end

## toggle buttons. Exclusive like a radio button
## was seeing crashing, so try this.
function togglebtn_press_event_cb(btnptr::Ptr, evt::Ptr, user_data)
    btn, btns, w, vals, labs = user_data
    val =  getproperty(btn, :active, Bool)
    if !val
        ## set button state
        for b in btns
            setproperty!(b, :active, b==btn)
        end
        ## set widget state
        push!(w.signal, vals[findfirst(labs, getproperty(btn, :label, AbstractString))])
    end
    true
end

function gtk_widget(widget::Options{:ToggleButtons})
    labs = collect(keys(widget.options))
    vals = collect(values(widget.options))

    block = @GtkBox(false)
    function make_button(lab)
        btn =  Gtk.@GtkToggleButton(lab)
        setproperty!(btn, :active, lab == widget.value_label)
        push!(block, btn)
        btn
    end
    btns = map(make_button, labs)

    ## widget -> signal
    ids = Dict()
    for btn in btns
        ids[btn] = signal_connect(togglebtn_press_event_cb, btn,
                                  "button-press-event", Bool, (Ptr{Gtk.GdkEventButton},), false,
                                  (btn, btns, widget, vals, labs))
        ## ids[btn] = signal_connect(btn, :button_press_event) do _,__
        ##     val =  getproperty(btn, :active, Bool)
        ##     if !val
        ##         ## set button state
        ##         for b in btns
        ##             setproperty!(b, :active, b==btn)
        ##         end
        ##         ## set widget state
        ##         push!(widget.signal, vals[findfirst(labs, getproperty(btn, :label, AbstractString))])
        ##     end
        ##     true                # impt: stop evennt propogation
        ## end
    end

    ## ## signal -> widget
    Reactive.foreach(widget.signal) do val
        ## get index from val
        index = findfirst(vals, val)
        lab = labs[index]
        for btn in btns
            signal_handler_block(btn, ids[btn])
            setproperty!(btn, :active, getproperty(btn, :label, AbstractString) == lab)
            signal_handler_unblock(btn, ids[btn])
        end
    end

    block
end

## buttongroup -- non exclusive

function buttongroup_cb(o::Ptr, user_data)
    widget, btn, labs, vals = user_data
    val =  getproperty(btn, :active, Bool)
    values = widget.signal.value
    lab = getproperty(btn, :label, AbstractString)
    i = findfirst(labs, lab)
            if val
                !(vals[i] in values) && push!(values, vals[i])
            else
                (vals[i] in values) && (values = filter(x -> vals[i] != x, values))
            end
            push!(widget.signal, values)
    nothing
end
function gtk_widget(widget::VectorOptions{:ButtonGroup})
    labs = collect(keys(widget.options))
    vals = collect(values(widget.options))

    block = @GtkBox(false)

    btns = Dict()
    for (lab, val) in zip(labs, vals)
        btn =  Gtk.@GtkToggleButton(lab)
        btns[lab] = btn
        push!(block, btn)
    end

    ## widget -> signal
    ids = Dict()
    for (lab, btn) in btns
        ids[btn] = signal_connect(buttongroup_cb, btn, "toggled", Void, (), false, (widget, btn, labs, vals))
        ## ids[btn] = signal_connect(btn, :toggled) do btn, xs...
        ##     val =  getproperty(btn, :active, Bool)
        ##     values = widget.signal.value
        ##     lab = getproperty(btn, :label, AbstractString)
        ##     i = findfirst(labs, lab)
        ##     if val
        ##         !(vals[i] in values) && push!(values, vals[i])
        ##     else
        ##         (vals[i] in values) && (values = filter(x -> vals[i] != x, values))
        ##     end
        ##     push!(widget.signal, values)
        ## end
    end

    ## signal -> widget
    Reactive.foreach(widget.signal) do values

        indices = [findfirst(vals, v) for v in values]
        selectedlabs = labs[indices]

        for (lab, btn) in btns
            signal_handler_block(btn, ids[btn])
            setproperty!(btn, :active, getproperty(btn, :label, AbstractString) in selectedlabs)
            signal_handler_unblock(btn, ids[btn])
        end

    end

    ## set initial *after* foreach call
    for (lab, val) in zip(labs, vals)
        btn =  btns[lab]
        setproperty!(btn, :active, val in widget.values)
    end

    block
end

## select -- a grid
function gtk_widget(widget::Options{:Select})

    labs = collect(map(AbstractString,keys(widget.options)))
    vals = collect(map(AbstractString,values(widget.options)))

    m = @GtkListStore(eltype(labs))
    block = @GtkScrolledWindow()
    obj = @GtkTreeView()
    [setproperty!(obj, x, true) for  x in [:hexpand, :vexpand]]
    push!(block, obj)

    Gtk.G_.model(obj, m)
    for lab in labs
        push!(m, (lab,))
    end

    cr = @GtkCellRendererText()
    col = @GtkTreeViewColumn(widget.label, cr, Dict([("text",0)]))
    push!(obj, col)

    ## initial choice
    index = findfirst(labs, widget.value_label)
    selection = Gtk.G_.selection(obj)
    store = getproperty(obj, :model, Gtk.GtkListStoreLeaf)
    iter = Gtk.iter_from_index(store, index)
    Gtk.select!(selection, iter)

    ## set up callback widget -> signal
    id = signal_connect(selection, :changed) do args...
        iter = selected(selection)
        j = Gtk.get_string_from_iter(Gtk.GtkTreeModel(store), iter)
        i = parse(Int, j) + 1
        #map(int, split(get_string_from_iter(treeModel, iter), ":")) + 1
#        i = Gtk.index_from_iter(store, iter)
        push!(widget.signal, vals[i])
    end

    ## push! -> update UI
    Reactive.foreach(widget.signal) do val
        signal_handler_block(selection, id)
        index = findfirst(vals, val)
        iter = Gtk.iter_from_index(store, index)
        Gtk.select!(selection, iter)
        signal_handler_unblock(selection, id)
    end

    block

end

##################################################
## Output widgets

## CairoGraphic
function gtk_widget(widget::CairoGraphic)
    if widget.obj == nothing

        obj = @GtkCanvas(widget.width, widget.height)
        ## how to make winston draw here? Here we store canvas in obj and override push!
        ## is there a more natural way??
        widget.obj = obj
    end

    widget.obj
end

## CairoImageSurface
function gtk_widget(widget::CairoImageSurface)
    if widget.obj == nothing
        obj = @GtkCanvas(widget.width, widget.height)
        widget.obj = obj
        Gtk.draw(obj) do canvas
            ctx = Cairo.getgc(canvas)
            Cairo.save(ctx)
            Cairo.reset_transform(ctx)
            Cairo.image(ctx, widget.surf, 0, 0, Cairo.width(ctx), Cairo.height(ctx))
        end
    end

    widget.obj
end

## Text area.
## unfortunately, setting the font doesn't seem to work.
function gtk_widget(widget::Textarea)
    if widget.obj == nothing
        obj = @GtkTextView()
        block = @GtkScrolledWindow()
        [setproperty!(obj, x, true) for  x in [:hexpand, :vexpand]]
        push!(block, obj)
        setproperty!(obj, :editable, false)

        if widget.buffer == nothing
            widget.buffer = getproperty(obj, :buffer, GtkTextBuffer)
        else
            setproperty!(obj, :buffer, widget.buffer)
        end
        setproperty!(widget.buffer, :text, widget.value)

        widget.obj = block
    end
    widget.obj
end

## change text in view
function Base.push!(obj::Textarea, value::AbstractString)
    if obj.buffer != nothing
        setproperty!(obj.buffer, :text, value)
        obj.value = value
    end
    value
end

## Label
function gtk_widget(widget::Label)
    if widget.obj == nothing
        obj = @GtkLabel(widget.value)
        setproperty!(obj, :selectable, true)
        setproperty!(obj, :use_markup, true)

        widget.obj = obj
    end

    widget.obj
end

function Base.push!(widget::Label, value::AbstractString)
    if widget.obj != nothing
        setproperty!(widget.obj, :label, value)
        setproperty!(widget.obj, :use_markup, true)
        widget.value = value
    end
    value
end

## Catch all `push!` for non AbstractString for label or textarea
typealias TextOrLabel Union{Textarea, Label}
Base.push!{T <: AbstractString}(obj::TextOrLabel, value::Vector{T}) = push!(obj, join(value, "\n"))
Base.push!(obj::TextOrLabel, value::Reactive.Signal) = push!(obj, Reactive.value(value))
Base.push!(obj::TextOrLabel, value) = push!(obj, to_string(value))

## Image
function gtk_widget(widget::Image)
    if widget.obj == nothing
        obj = @GtkImage()
        widget.obj = obj
    end
    if widget.value != nothing
        Gtk.G_.file(obj, widget.value)
    end
    widget.obj
end
function Base.push!(widget::Image, val::AbstractString)
    if widget.obj != nothing
        widget.value = val
        Gtk.G_.file(widget.obj, val)
    end
end
##################################################
## Decorative
## icon (no obj property) as we don't push onto these
function gtk_widget(widget::Icon)
    obj = @GtkImage()
    setproperty!(obj, :icon_name, widget.stock_id)

    if widget.tile != nothing
        child = gtk_widget(widget.tile)
        setproperty!(child, :always_show_image, true)
        Gtk.G_.image(child, obj)
        obj = child
    end

    obj
end

function gtk_widget(widget::Tooltip)
    obj = gtk_widget(widget.tile)
    setproperty!(obj, :tooltip_text, widget.text)
    obj
end

## Progress bar
function gtk_widget(widget::Progress)
    obj = @GtkProgressBar()
    widget.obj = obj

    push!(widget, widget.value)
    obj
end

## push value in range of obj.range
Base.push!(obj::Progress, value::Reactive.Signal) = push!(obj, Reactive.value(value))
function Base.push!(widget::Progress, value)
    if widget.obj != nothing
        frac = clamp((value - first(widget.range)) / (last(widget.range) - first(widget.range)), 0, 1)
        setproperty!(widget.obj, :fraction, frac)
    end
end

##################################################
## Layouts

## Attributes
function gtk_widget(widget::Size)
    # set size of tile
    obj = gtk_widget(widget.tile)

    Gtk.G_.size_request(obj, widget.w_px.value, widget.h_px.value)

    obj
end

function gtk_widget(widget::Grow)
    obj = gtk_widget(widget.tile)
    ## how to use factor? It is in [0,1]?
    widget.factor > 0 && :horizontal in widget.direction && Gtk.G_.hexpand(obj, true)
    widget.factor > 0 && :vertical in widget.direction && Gtk.G_.vexpand(obj, true)

    obj
end

function gtk_widget(widget::Shrink)
    obj = gtk_widget(widget.tile)
    ## how to use factor?
    if widget.factor > 0
        Gtk.G_.hexpand(obj, false)
        Gtk.G_.vexpand(obj, false)
    end

    obj
end

function gtk_widget(widget::Pad)
    # set size of tile
    obj = gtk_widget(widget.tile)

    ## group on keywaor...
    :left in widget.sides   && Gtk.G_.margin_left(obj, widget.len.value)
    :right in widget.sides && Gtk.G_.margin_right(obj, widget.len.value)
    :top in widget.sides    && Gtk.G_.margin_top(obj, widget.len.value)
    :bottom in widget.sides && Gtk.G_.margin_bottom(obj, widget.len.value)

    obj
end

function gtk_widget(widget::Align)
    # set alignment
    obj = gtk_widget(widget.tile)

    Gtk.G_.halign(obj, widget.halign)
    Gtk.G_.valign(obj, widget.valign)

    obj
end

## containers
function gtk_widget(widget::FlowContainer)
    obj = widget.obj =  @GtkBox(widget.direction == "vertical") # GtkFlowBox...
    for child in widget.children
        push!(obj, gtk_widget(child))
    end
    obj
end

function gtk_widget(widget::Separator)
    obj = @GtkLabel("* * * * *")        # poor man's separator...
    Gtk.G_.margin_top(obj, 5)
    Gtk.G_.margin_bottom(obj, 5)
    obj
end

function gtk_widget(widget::Tabs)
    obj = @GtkNotebook()
    for (label, child) in zip(widget.labels, widget.children)
        gchild = gtk_widget(child)
        push!(obj, gchild, label)
        show(gchild)
    end
    Gtk.G_.current_page(obj, widget.initial-1)
    obj
end

formlabel(widget::Widget) = :label in fieldnames(widget) ? widget.label : ""
formlabel(widget::Button) = ""
formlabel(widget) = ""

function gtk_widget(widget::FormLayout)
    obj = @GtkGrid()
    setproperty!(obj, :hexpand, true)
    setproperty!(obj, :row_spacing, 5)
    setproperty!(obj, :column_spacing, 5)

    for (row, child) in enumerate(widget.children)
        al = @GtkAlignment(1.0, 0.0, 0.0, 0.0)
        setproperty!(al, :right_padding, 5)
        setproperty!(al, :left_padding, 5)

        child_widget = gtk_widget(child)
        setproperty!(child_widget, :hexpand, true)

        push!(al, @GtkLabel(formlabel(child)))

        obj[1, row] = al
        obj[2, row] = child_widget
    end

    obj
end

## Toolbar
function gtk_widget(widget::Toolbar)
    obj = @GtkToolbar()
    Gtk.G_.style(obj, Gtk.GConstants.GtkToolbarStyle.GTK_TOOLBAR_BOTH)

    for child in widget.children
        tbchild = gtk_toolbar_widget(child)
        push!(obj, tbchild)
    end

    obj
end

function gtk_toolbar_widget(widget::Button)
    obj = @GtkToolButton("")
    Gtk.G_.label(obj, widget.label)

    ## widget -> signal
    id = signal_connect(button_cb, obj, "clicked", Void, (),  false, (widget, obj))
#    id = signal_connect(obj, :clicked) do obj, args...
#        push!(widget, nothing)
#    end

    obj
end
function gtk_toolbar_widget(widget::ToggleButton)
    obj = @GtkToggleToolButton(widget.label)
    Gtk.G_.label(obj, widget.label)
    setproperty!(obj, :active, widget.value)

    ## widget -> signal
    id = signal_connect(obj, :toggled) do obj, args...
        val = getproperty(obj, :active, Bool)
        push!(widget.signal, val) # call
    end

    obj
end

## XXX This isn't working XXX
function gtk_toolbar_widget(widget::MenuButton)
    obj = @GtkMenuToolButton(widget.label)
    Gtk.G_.label(obj, widget.label)

    m = @GtkMenu()
    Gtk.G_.menu(obj, m)
    push!(m , gtk_menu_widget(menu(widget.children...)))
    showall(obj)
    obj
end

function gtk_toolbar_widget(child::Separator)
    obj = @GtkSeparatorToolItem()
    obj
end

function gtk_toolbar_widget(widget::Icon)
    img = @GtkImage()
    setproperty!(img, :icon_name, widget.stock_id)

    obj = gtk_toolbar_widget(widget.tile) # must be non-empty or an error
    Gtk.G_.icon_widget(obj, img)

    obj
end

function gtk_toolbar_widget(widget::Tooltip)
    obj = gtk_toolbar_widget(widget.tile)
    setproperty!(obj, :tooltip_text, widget.text)
    obj
end

## How to put in a spacer in a toolbar?
## We use vskip() to place a spring in a toolbar
function gtk_toolbar_widget(widget::Size)
    obj = @GtkSeparatorToolItem()
    Gtk.G_.draw(obj, false)
    Gtk.G_.expand(obj, true)
    obj
end

## menu
function gtk_widget(widget::Menu)
    obj = Gtk.@GtkMenuBar()
    Gtk.G_.hexpand(obj, true)
    for child in widget.children
        push!(obj, gtk_menu_widget(child))
    end

    obj
end

function gtk_menu_widget(widget::Button)
    obj = @GtkMenuItem(widget.label)

    id = signal_connect(button_cb, obj, "activate", Void, (),  false, (widget, obj))
    ## widget -> signal
#    id = signal_connect(obj, :activate) do obj, args...
#        push!(widget.signal, Reactive.value(Interact.signal(widget))) # call
#    end

    obj

end

function gtk_menu_widget(widget::ToggleButton)
    error("This menu item is not implemented")
    ## Hack to work around no GtkCheckMenuItem
    ## XXX Doesn't work
    ## obj = @GtkMenuItem()
    ## ## remove
    ## child = Gtk.G_.child(obj)
    ## ccall((:gtk_container_remove, Gtk.libgtk), Void, (Ptr{Gtk.GObject}, Ptr{Gtk.GObject}), obj, child)

    ## ## add back
    ## child = @GtkToggleToolButton(widget.label)
    ## Gtk.G_.label(child, widget.label)
    ## push!(obj, child)

    ## show(child)

    ## ## widget -> signal
    ## id = signal_connect(child, :toggled) do obj, args...
    ##     val = getproperty(child, :active, Bool)
    ##     push!(widget.signal, val) # call
    ## end

    ## obj
end

function gtk_menu_widget(widget::Separator)
    obj = @GtkSeparatorMenuItem()
    obj
end

function gtk_menu_widget(widget::Icon)
    ## images don't work?
    obj = gtk_menu_widget(widget.tile)  # will haver error if no tile

    img = @GtkImage()
    setproperty!(img, :icon_name, widget.stock_id)

    Gtk.G_.icon_widget(obj, img)
end

function gtk_menu_widget(widget::Menu)
    obj = @GtkMenuItem(widget.label)
    submenu = @GtkMenu(obj)
    for child in widget.children
         push!(submenu, gtk_menu_widget(child))
    end

    obj
end

## XXX This isn't working XXX
function gtk_widget(widget::MenuButton)
    error("GtkMenuButton is not (yet) in Gtk.jl")
    #obj = @GtkMenuButton()              # error, not (yet) part of Gtk.jl
    for child in widget.chldren
        push!(obj, gtk_menu_widget(child))
    end
    obj
end

##
function gtk_widget(widget::Window)
    obj = @GtkWindow(title=widget.title)
    widget.obj = obj
    widget.width > 0 && widget.height > 0 && resize!(obj, widget.width, widget.height)

    ## interior packing box...
    box = @GtkBox(true)
    Gtk.G_.hexpand(box, true)
    Gtk.G_.vexpand(box, true)
    push!(obj, box)

    for child in widget.children
        push!(box, gtk_widget(child))
    end
    showall(obj)
    obj
end
function Base.display(widget::Window)
    showall(gtk_widget(widget))
end
Base.display(widget::Widget) = nothing

## Main window
function gtk_widget(widget::MainWindow)
    obj =  @GtkWindow(title=widget.title)
    resize!(obj, widget.width, widget.height)
    widget.window = obj
    signal_connect(obj, :destroy) do win
        closerefs!(widget.refs)
    end

    push!(obj, gtk_widget(formlayout(widget.children...)))

    obj
end

## for displaying an @manipulate object, we need this
function Base.display(x::ManipulateWidget)
    Reactive.foreach(a -> show_outwidget(x.w, a), x.a)
    ## hack to get output widgets passed in
    tmp = filter(u->!isa(u, GtkInteract.OutputWidget), x.w.children)[1]
    push!(tmp, value(tmp))
    nothing
end

## Catch all for showing outwidget
function show_outwidget(w::GtkInteract.MainWindow, x)

    if w.out == nothing
        w.out = label("")
        push!(w, grow(w.out))
        display(w)
    end
    x != nothing &&  push!(w.out, to_string(x))
end

## convert object to string for display through label
to_string(x::AbstractString) = x
to_string(x) = sprint(io -> show(io, "text/plain", x))

## Dialogs
function gtk_widget(widget::MessageBox)
    if widget.style == :warnt
        fn = Gtk.warn_dialog
    elseif widget.style == :error
        fn = Gtk.error_dialog
    else
        fn = Gtk.info_dialog
    end

    fn(widget.msg)

end

function gtk_widget(widget::ConfirmBox)
    Gtk.ask_dialog(widget.msg)
end

function gtk_widget(widget::InputBox)
    ret, val = Gtk.input_dialog(widget.msg, widget.default)
    if ret == 0
        val = utf8("")
    end
    val
end

function gtk_widget(widget::OpenFile)
    open_dialog(widget.title)
end

function gtk_widget(widget::SaveFile)
    save_dialog(widget.title)
end

function gtk_widget(widget::SelectDir)
    dlg = @GtkFileChooserDialog(widget.title, Gtk.GtkNullContainer(),
                                Gtk.GConstants.GtkFileChooserAction.SELECT_FOLDER,
                                (("_Cancel", Gtk.GConstants.GtkResponseType.CANCEL),
                                 ("_Save",   Gtk.GConstants.GtkResponseType.ACCEPT))
                                )
    dlgp = GtkFileChooser(dlg)
    response = run(dlg)
    if response == Gtk.GConstants.GtkResponseType.ACCEPT
        selection = bytestring(GAccessor.filename(dlgp))
    else
        selection = utf8("")
    end
    destroy(dlg)
    selection
end

## Typeography
gtk_widget(obj::Bold) = gtk_widget(label("<b>$(obj.label)</b>"))
gtk_widget(obj::Emph) = gtk_widget(label("<i>$(obj.label)</i>"))
gtk_widget(obj::Code) = gtk_widget(label("<tt>$(obj.label)</tt>"))


text:
# Based on https://gist.github.com/simonbyrne/5443843

function IBM2Float64(x::Array{Uint8, 1})
    ibm = reinterpret(Uint64,x)[1]
    if ibm == 0
        println("ibm 0")
        return zero(Float64)
    end

    mantissa = ibm & 0x00ff_ffff_ffff_ffff
    exponent = ibm & 0x7f00_0000_0000_0000
    if mantissa == 0
        # if exponent == 0x5f00_0000_0000_0000
        #     #._
        # elseif  exponent == 0x2e00_0000_0000_0000
        #     #.
        # elseif  exponent == 0x4100_0000_0000_0000
        #     #.A
        # elseif exponent == 0x4200_0000_0000_0000
        #     #.B
        #     #.....
        # elseif exponent == 0x5a00_0000_0000_0000
        #     #.Z
        # end
        return NA
    end

    # normalise mantissa to base 2.
    offset = 1
    while mantissa & 0x0080_0000_0000_0000 == 0
        offset += 1
        mantissa <<= 1
    end
    # drop leading 1 (8th bit), and shift remainder to 12th bit.
    mantissa = (mantissa $ 0x0080_0000_0000_0000) >> 3

    sign = ibm & 0x8000_0000_0000_0000

    # convert from base 16 to base 2, adjust for normalisation.
    exponent = ((((exponent >> 56) - 64) << 2) - offset + 1023) << 52

    return reinterpret(Float64, sign | exponent | mantissa)
end

ext: jl


text:
@fastmath @views function ltc!(du,u,p,t, n_in, n_neurons, w_sens, w_syns)
  v = u
  stim, C, G_L, E_L, G_sens, s_sens, h_sens, E_sens, G_syns, s_syns, h_syns, E_syns = ltc_unpack_p(n_in,n_neurons,p)
  #@unpack stim, C, G_L, E_L, G_sens, s_sens, h_sens, E_sens, G_syns, s_syns, h_syns, E_syns = p

  #I_senss = reshape([G_sens[s,n] * NNlib.σ((stim[s] - h_sens[s,n]) * s_sens[s,n]) * (v[n] - E_sens[s,n]) for s in 1:size(G_sens,1), n in 1:size(G_sens,2)], size(G_sens,1), :)
  #I_synss = reshape([G_syns[s,n] * NNlib.σ((v[s] - h_syns[s,n]) * s_syns[s,n]) * (v[n] - E_syns[s,n]) for s in 1:size(G_syns,1), n in 1:size(G_syns,2)], size(G_syns,1), :)

  #I_senss = zeros(eltype(stim), n_in, n_neurons)
  #I_synss = zeros(eltype(stim), n_neurons, n_neurons)
  I_senss = Zygote.Buffer(stim, n_in, n_neurons)
  I_synss = Zygote.Buffer(stim, n_neurons, n_neurons)

  #I_sens = zeros(eltype(stim), n_neurons)
  #I_syns = zeros(eltype(stim), n_neurons)
  # TODO: tullio?
  for n in 1:n_neurons
    for s in 1:n_in
      I_senss[s,n] = w_sens[s,n] * G_sens[s,n] * mysigm((stim[s] - h_sens[s,n]) * s_sens[s,n]) * (v[n] - E_sens[s,n])
      #I_sens[n] += w_sens[s,n] * G_sens[s,n] * mysigm((stim[s] - h_sens[s,n]) * s_sens[s,n]) * (v[n] - E_sens[s,n])
    end
  end
  for n in 1:n_neurons
    for s in 1:n_neurons
      I_synss[s,n] = w_syns[s,n] * G_syns[s,n] * mysigm((v[s] - h_syns[s,n]) * s_syns[s,n]) * (v[n] - E_syns[s,n])
      #I_syns[n] += w_syns[s,n] * G_syns[s,n] * mysigm((v[s] - h_syns[s,n]) * s_syns[s,n]) * (v[n] - E_syns[s,n])
    end
  end
  I_senss = copy(I_senss)
  I_synss = copy(I_synss)

  #I_senss = reshape([G_sens[s,n] * mysigm((stim[s] - h_sens[s,n]) * s_sens[s,n]) * (v[n] - E_sens[s,n]) for s in 1:size(G_sens,1), n in 1:size(G_sens,2)], size(G_sens,1), :)
  #I_synss = reshape([G_syns[s,n] * mysigm((v[s] - h_syns[s,n]) * s_syns[s,n]) * (v[n] - E_syns[s,n]) for s in 1:size(G_syns,1), n in 1:size(G_syns,2)], size(G_syns,1), :)

  I_sens = reshape(sum(I_senss, dims=1), :)
  I_syns = reshape(sum(I_synss, dims=1), :)

  # I_sens = zeros(eltype(stim), n_neurons)
  # I_syns = zeros(eltype(stim), n_neurons)
  # @inbounds for n in 1:n_neurons
  #   for s in 1:n_in
  #     #w_sens[s,n] == 0 && continue
  #     I_sens[n] += w_sens[s,n] * G_sens[s,n] * mysigm((stim[s] - h_sens[s,n]) * s_sens[s,n]) * (v[n] - E_sens[s,n])
  #   end
  #   # I_sens[n] += sum(G_sens[:,n] .* mysigm.((stim .- h_sens[:,n]) .* s_sens[:,n]) .* (v[n] .- E_sens[:,n]))
  #   for s in 1:n_neurons
  #     #w_sens[s,n] == 0 && continue
  #     I_syns[n] += w_syns[s,n] * G_syns[s,n] * mysigm((v[s] - h_syns[s,n]) * s_syns[s,n]) * (v[n] - E_syns[s,n])
  #   end
  #   # I_syns[n] += sum(G_syns[:,n] .* mysigm.((v .- h_syns[:,n]) .* s_syns[:,n]) .* (v[n] .- E_syns[:,n]))
  # end

  #du .= (-((G_L .* (v .- E_L)) .+ I_sens .+ I_syns) ./ C) #.* dZ
  @.. du = (-((G_L * (v - E_L)) + I_sens + I_syns) / C) #.* dZ

  nothing
end

function ltc_oop!(u,p,t, n_in, n_neurons, w_sens, w_syns)
  v = u
  #stim, C, G_L, E_L, G_sens, s_sens, h_sens, E_sens, G_syns, s_syns, h_syns, E_syns = ltc_unpack_p(n_in,n_neurons,p)
  @unpack stim, C, G_L, E_L, G_sens, s_sens, h_sens, E_sens, G_syns, s_syns, h_syns, E_syns = p

  I_sens = zeros(eltype(stim), n_neurons)
  I_syns = zeros(eltype(stim), n_neurons)

  # TODO: tullio?
  for n in 1:n_neurons
    for s in 1:n_in
      #I_senss[s,n] = w_sens[s,n] * G_sens[s,n] * mysigm((stim[s] - h_sens[s,n]) * s_sens[s,n]) * (v[n] - E_sens[s,n])
      I_sens[n] += w_sens[s,n] * G_sens[s,n] * mysigm((stim[s] - h_sens[s,n]) * s_sens[s,n]) * (v[n] - E_sens[s,n])
    end
  end

  for n in 1:n_neurons
    for s in 1:n_neurons
      #I_synss[s,n] = w_syns[s,n] * G_syns[s,n] * mysigm((v[s] - h_syns[s,n]) * s_syns[s,n]) * (v[n] - E_syns[s,n])
      I_syns[n] += w_syns[s,n] * G_syns[s,n] * mysigm((v[s] - h_syns[s,n]) * s_syns[s,n]) * (v[n] - E_syns[s,n])
    end
  end
  #I_sens = reshape(sum(I_senss, dims=1), :)
  #I_syns = reshape(sum(I_synss, dims=1), :)

  du = (-((G_L * (v - E_L)) + I_sens + I_syns) / C) #.* dZ

  return du
end

function ltc_init_u0p(n_in, n_neurons; T=Float32, wiring=nothing)
  v = rand_uniform(T, 0.01, 0.1, n_neurons) #.± 0.1f0

  stim   = fill(T(0), n_in)
  C      = rand_uniform(T,  1.0,     1.002,   n_neurons)
  G_L    = rand_uniform(T,  0.0001,  0.1,     n_neurons)
  E_L    = rand_uniform(T, -0.3,     0.3,     n_neurons) #.± 0.1f0
  G_sens = rand_uniform(T,  0.0001,  0.1,     n_in, n_neurons)
  s_sens = rand_uniform(T,  3,       8,       n_in, n_neurons)
  h_sens = rand_uniform(T,  0.3,     0.8,     n_in, n_neurons)
  E_sens = rand_uniform(T, -0.9,     0.9,     n_in, n_neurons)
  G_syns = rand_uniform(T,  0.0001,  0.1,     n_neurons, n_neurons)
  s_syns = rand_uniform(T,  3,       8,       n_neurons, n_neurons)
  h_syns = rand_uniform(T,  0.3,     0.8,     n_neurons, n_neurons)
  E_syns = rand_uniform(T, -0.9,     0.9,     n_neurons, n_neurons)

  if wiring !== nothing
    E_sens = wiring.matrices[:E_sens]
    E_syns = wiring.matrices[:E_syns]
  end

  u0 = v
  p = vcat( stim,
    C, G_L, E_L,
    vec(G_sens), vec(s_sens), vec(h_sens), vec(E_sens),
    vec(G_syns), vec(s_syns), vec(h_syns), vec(E_syns),
  )
  return u0, p
end

function ltc_u0p_ca(wiring; T=Float32)
  n_in = wiring.s_in
  n_neurons = wiring.n_total
  v = rand_uniform(T, 0.01, 0.1, n_neurons) #.± 0.1f0

  stim   = fill(T(0), n_in)

  C      = rand_uniform(T,  1.0,     1.002,   n_neurons)
  G_L    = rand_uniform(T,  0.0001,  0.1,     n_neurons)
  E_L    = rand_uniform(T, -0.3,     0.3,     n_neurons) #.± 0.1f0
  G_sens = rand_uniform(T,  0.0001,  0.1,     n_in, n_neurons)
  s_sens = rand_uniform(T,  3,       8,       n_in, n_neurons)
  h_sens = rand_uniform(T,  0.3,     0.8,     n_in, n_neurons)
  E_sens = rand_uniform(T, -0.9,     0.9,     n_in, n_neurons)
  G_syns = rand_uniform(T,  0.0001,  0.1,     n_neurons, n_neurons)
  s_syns = rand_uniform(T,  3,       8,       n_neurons, n_neurons)
  h_syns = rand_uniform(T,  0.3,     0.8,     n_neurons, n_neurons)
  E_syns = rand_uniform(T, -0.9,     0.9,     n_neurons, n_neurons)

  #E_sens = wiring.matrices[:E_sens]
  #E_syns = wiring.matrices[:E_syns]

  #u0 = ComponentArray{T}(v=v)
  u0 = v
  p = ComponentArray{T}( stim=stim,
    C=C, G_L=G_L, E_L=E_L,
    G_sens=G_sens, s_sens=s_sens, h_sens=h_sens, E_sens=E_sens,
    G_syns=G_syns, s_syns=s_syns, h_syns=h_syns, E_syns=E_syns)

  u0, p
end

function ltc_bounds(n_in, n_neurons; T=Float32)
  lb = vcat(
    [1.0 for _ in 1:n_neurons],               # C
    [0 for _ in 1:n_neurons],                 # G_L
    [-1 for _ in 1:n_neurons],                # E_L
    [0 for _ in 1:n_in*n_neurons],            # G_sens
    [2 for _ in 1:n_in*n_neurons],            # s_sens
    [-0.3 for _ in 1:n_in*n_neurons],          # h_sens
    [-1 for _ in 1:n_in*n_neurons],           # E_sens
    [0 for _ in 1:n_neurons*n_neurons],       # G_syns
    [2 for _ in 1:n_neurons*n_neurons],       # s_syns
    [-0.3 for _ in 1:n_neurons*n_neurons],     # h_syns
    [-1 for _ in 1:n_neurons*n_neurons],      # E_syns

    [-0.4 for _ in 1:n_neurons])             # v

  ub = vcat(
    [2.0 for _ in 1:n_neurons],
    [0.4 for _ in 1:n_neurons],
    [1 for _ in 1:n_neurons],
    [1 for _ in 1:n_in*n_neurons],
    [9 for _ in 1:n_in*n_neurons],
    [0.9 for _ in 1:n_in*n_neurons],
    [1 for _ in 1:n_in*n_neurons],
    [1 for _ in 1:n_neurons*n_neurons],
    [9 for _ in 1:n_neurons*n_neurons],
    [0.9 for _ in 1:n_neurons*n_neurons],
    [1 for _ in 1:n_neurons*n_neurons],

    [0.8 for _ in 1:n_neurons]) # v)

  T.(lb), T.(ub)
end

function ltc_unpack_p(n_in,n_neurons,p)
  stiml = n_in
  Cl = n_neurons
  G_Ll = n_neurons
  E_Ll = n_neurons
  sensl = n_in*n_neurons
  synsl = n_neurons*n_neurons
  G_sensl = sensl
  s_sensl = sensl
  h_sensl = sensl
  E_sensl = sensl
  G_synsl = synsl
  s_synsl = synsl
  h_synsl = synsl
  E_synsl = synsl

  s = 1
  stim = @view p[s:s+stiml-1]
  s += stiml
  C = @view p[s:s+Cl-1]
  s += Cl
  G_L = @view p[s:s+G_Ll-1]
  s += G_Ll
  E_L = @view p[s:s+E_Ll-1]
  s += E_Ll
  G_sens = @view p[reshape(s:s+G_sensl-1, n_in, n_neurons)]
  s += G_sensl
  s_sens = @view p[reshape(s:s+s_sensl-1, n_in, n_neurons)]
  s += s_sensl
  h_sens = @view p[reshape(s:s+h_sensl-1, n_in, n_neurons)]
  s += h_sensl
  E_sens = @view p[reshape(s:s+E_sensl-1, n_in, n_neurons)]
  s += E_sensl

  G_syns = @view p[reshape(s:s+G_synsl-1, n_neurons, n_neurons)]
  s += G_synsl
  s_syns = @view p[reshape(s:s+s_synsl-1, n_neurons, n_neurons)]
  s += s_synsl
  h_syns = @view p[reshape(s:s+h_synsl-1, n_neurons, n_neurons)]
  s += h_synsl
  E_syns = @view p[reshape(s:s+E_synsl-1, n_neurons, n_neurons)]
  s += E_synsl
  return stim, C, G_L, E_L, G_sens, s_sens, h_sens, E_sens, G_syns, s_syns, h_syns, E_syns
end

function ltc_unpack_p_ca(n_in,n_neurons,p)
end

# function LTC(n_in, n_neurons, solver, sensealg;
#   T=Float32, tspan=T.((0, 1)), n_sens=n_neurons, n_out=n_neurons,
#   w_sens=ones(T, n_in, n_neurons),
#   w_syns=ones(T, n_neurons, n_neurons),
#   wiring=nothing,
#   mtkize=true, gen_jac=false, kwargs...)

#   if wiring !== nothing
#     w_sens = wiring.matrices[:w_sens]
#     w_syns = wiring.matrices[:w_syns]
#   end

#   lb, ub = ltc_bounds(n_in, n_neurons; T)
#   u0, p = ltc_init_u0p(n_in, n_neurons; T, wiring)
#   #u0, p = ltc_sys_u0p_ca(n_in, n_neurons; T)
#   dudt!(du,u,p,t) = ltc!(du,u,p,t, n_in, n_neurons, w_sens, w_syns)

#   rnncell = BCTRNNCell(n_in, n_sens, n_neurons, n_out, solver, sensealg, dudt!, u0, tspan, p, lb, ub; mtkize, gen_jac, kwargs...)
#   MyRecur(rnncell)
# end

function LTC(wiring, solver, sensealg;
  T=Float32, tspan=T.((0, 1)),
  #wiring=nothing,
  mtkize=true, gen_jac=false, kwargs...)

  w_sens = wiring.matrices[:w_sens]
  w_syns = wiring.matrices[:w_syns]

  s_in = wiring.s_in
  n_neurons = wiring.n_total

  lb, ub = ltc_bounds(s_in, n_neurons; T)
  u0, p = ltc_init_u0p(s_in, n_neurons; T, wiring)
  #u0, p = ltc_u0p_ca(wiring; T)
  dudt!(du,u,p,t) = ltc!(du,u,p,t, s_in, n_neurons, w_sens, w_syns)

  rnncell = BCTRNNCell(wiring, solver, sensealg, dudt!, u0, tspan, p, lb, ub; mtkize, gen_jac, kwargs...)
  MyRecur(rnncell)
end

# function bcde_sys!(du,u,p,t, c_mapperp_re)
#   v, c = u
#   stim, c_mapperp_p, C, G_L, E_L, G_sens, s_sens, h_sens, E_sens, G_syns, s_syns, h_syns, E_syns = p

#   I_leak = G_L .* (v .- E_L)
#   I_senss = zeros(eltype(v), n_stim, size(v,1))
#   I_synss = zeros(eltype(v), size(v,1), size(v,1))

#   for n in 1:size(G_sens,2)
#     for s in 1:size(G_sens,1)
#       I_senss[s,n] = G_sens[s,n] * NNlib.σ((stim[s] - h_sens[s,n]) * s_sens[s,n]) * (v[n] - E_sens[s,n])
#     end
#   end
#   for n in 1:size(G_syns,2)
#     for s in 1:size(G_syns,1)
#       I_synss[s,n] = G_syns[s,n] * NNlib.σ((v[s] - h_syns[s,n]) * s_syns[s,n]) * (v[n] - E_syns[s,n])
#     end
#   end
#   I_sens = vec(sum(I_senss, dims=1))
#   I_syns = vec(sum(I_synss, dims=1))

#   I_tot = I_leak .+ I_sens .+ I_syns

#   du[1:size(v,1)] .= (-(I_tot .+ c) ./ C) #.* dZ

#   cmin = [vcat(stim, c[i], v[i], I_tot[i], t) for i in 1:size(v, 1)]
#   du[size(v,1)+1:2*size(v,1)] .= c_mapperp_re(c_mapperp_p)(cmin...) #.* dZ

#   nothing
# end


text:
using Plots
sphere(x) = sum(x.^2)

dpi = 600
x = y = range(-3, stop=3, length=50)

l = Vector{Float64}()
for i in x
   for j in y
       v = [i j]
       append!(l, 1/sphere(v))
   end
end

surface(x,y,l, dpi=dpi)
png("./inverse_sphere_surface.png")
contour(x,y,l, dpi=dpi)
png("./inverse_sphere_contour.png")


text:
## Notation
"""
    ket([T = Float64], i::Integer, d::Integer) -> SparseVector{Complex{T}}

Create a vector representing the `i`th computational basis vector in dimension `d`.

## Example

```jldoctest
julia> ket(1,2)
2-element SparseVector{Complex{Float64},Int64} with 1 stored entry:
  [1]  =  1.0+0.0im

julia> collect(ans)
2-element Array{Complex{Float64},1}:
 1.0 + 0.0im
 0.0 + 0.0im
```
"""
ket(::Type{T}, i::Integer, d::Integer) where {T} = sparsevec([i], [one(complex(T))], d)

"""
    bra([T = Float64], i::Integer, d::Integer) -> SparseVector{Complex{T}}'

Create a dual vector representing the bra associated to `i`th computational basis vector in dimension `d`.

## Example

```jldoctest
julia> bra(1,2)
1×2 LinearAlgebra.Adjoint{Complex{Float64},SparseVector{Complex{Float64},Int64}}:
 1.0-0.0im  0.0-0.0im

julia> collect(ans)
1×2 Array{Complex{Float64},2}:
 1.0-0.0im  0.0-0.0im
```
"""
bra(::Type{T}, i::Integer, d::Integer) where {T} = ket(T, i, d)'

ket(i::Integer, d::Integer) = ket(Float64, i, d)
bra(i::Integer, d::Integer) = bra(Float64, i, d)

"""
    dm(ψ::AbstractVector)

Creates the density matrix version of a pure state `ψ` via the outer product.
"""
dm(ψ::AbstractVector) = ψ' ⊗ ψ

if VERSION < v"1.2-pre"
    (I::UniformScaling)(n::Integer) = Diagonal(fill(I.λ, n))
end

const ⊗ = kron

"""
    BB84_states([T::Type = Float64])

Generates the BB84 states `|0⟩`, `|1⟩`, `|-⟩`, and `|+⟩`, for use in [`guesswork`](@ref) or other functions. The
numeric type can be optionally specified by the argument.
"""
function BB84_states end

BB84_states() = BB84_states(Float64)

function BB84_states(::Type{T}) where {T<:Number}
    ketplus = (ket(T, 1, 2) + ket(T, 2, 2)) / sqrt(T(2))
    ketminus = (ket(T, 1, 2) - ket(T, 2, 2)) / sqrt(T(2))
    ketzero = ket(T, 1, 2)
    ketone = ket(T, 2, 2)
    return dm.([ketzero, ketone, ketminus, ketplus])
end

"""
    iid_copies(ρBs::AbstractVector{<:AbstractMatrix}, n::Integer) -> Vector{Matrix}

Create a vector of all states of the form ``ρ_1 \\otimes \\dotsm \\otimes ρ_n``
where the ``ρ_i`` range over the set `ρBs`.
"""
function iid_copies(ρBs::AbstractVector{<:AbstractMatrix}, n::Integer)
    inds = cartesian_power(length(ρBs), n)
    [foldl(⊗, ρBs[collect(I)]) for I in inds]
end

## Random states

"""
    simplexpt(unif)

Takes a vector of length `d-1` of numbers between `0.0` and `1.0` and converts it a point on the standard `d` dimensional simplex.
"""
function simplexpt(unif)
    d = length(unif) + 1
    T = eltype(unif)
    w = zeros(T, d + 1)
    w[2:d] .= sort(unif)
    w[d+1] = one(T)
    diff(w)
end

"""
    randprobvec([rng, T=Float64], d)

Generates points of type `T`, uniformly at random on the standard `d-1` dimensional simplex using an algorithm by [Smith and Tromble](http://www.cs.cmu.edu/~nasmith/papers/smith+tromble.tr04.pdf).

## Example

```julia
julia> randprobvec(3)
3-element Array{Float64,1}:
 0.24815974900033688
 0.17199716455672287
 0.5798430864429402

```
"""
randprobvec(rng::AbstractRNG, ::Type{T}, d::Integer) where T = simplexpt(rand(rng, T, d - 1))
randprobvec(::Type{T}, d::Integer) where T = randprobvec(Random.GLOBAL_RNG, T, d)
randprobvec(d::Integer) = randprobvec(Float64, d)

"""
    randunitary([rng, T=Float64], d)

Generates a unitary matrix of dimension `d` at random according to the Haar measure, using an algorithm described by Maris Ozols in ["How to generate a random unitary matrix"](http://home.lu.lv/~sd20008/papers/essays/Random%20unitary%20%5Bpaper%5D.pdf).
"""
function randunitary(rng::AbstractRNG, ::Type{T}, d::Integer) where T
    # `randn` for `BigFloat`'s isn't supported yet
    # (https://github.com/JuliaLang/julia/issues/17629)
    # so we convert afterwards instead (losing some randomness).
    rg1 = T.(randn(rng, d, d))
    rg2 = T.(randn(rng, d, d))
    RG = rg1 + im * rg2
    Q, R = qr(RG)
    r = diag(R)
    L = diagm(0 => r ./ abs.(r))
    return Q * L
end

randunitary(::Type{T}, d::Integer) where {T} = randunitary(Random.GLOBAL_RNG, T, d)
randunitary(d::Integer) = randunitary(Float64, d)

"""
    randdm([rng, T = Float64], d)

Generates a density matrix with numeric type `Complex{T}`, of dimension `d` at
random.

## Example

```julia
julia> randdm(2)
2×2 Array{Complex{Float64},2}:
 0.477118+0.0im        0.119848-0.0371569im
 0.119848+0.0371569im  0.522882+0.0im
```
"""
function randdm(rng::AbstractRNG, ::Type{T}, d::Integer) where T
    eigs = diagm(0 => randprobvec(rng, T, d))
    U = randunitary(rng, T, d)
    ρ = U * eigs * (U')
    return Matrix((ρ + ρ') / 2)
end

randdm(::Type{T}, d::Integer) where T = randdm(Random.GLOBAL_RNG, T, d)
randdm(d) = randdm(Float64, d)


text:
function buildNet(numInputs::Integer, numOutputs::Integer, hiddenStructure::AVec{Int};
                  hiddenMapping::Mapping = SigmoidMapping(),
                  finalMapping::Mapping = SigmoidMapping(),
                  params = NetParams(),
                  solverParams = SolverParams(),
                  inputTransformer::Transformer = IdentityTransformer(),
                  wgt::Weight = EqualWeight(),
                  LAYER::DataType = Layer
                  )
  layers = LAYER[]
  nin = numInputs

  # the first layer will get the "input" dropout probability
  pDropout = getDropoutProb(params, true)

  for nout in hiddenStructure

    # push the hidden layer
    push!(layers, LAYER(nin,
                        nout,
                        hiddenMapping,
                        params.updater,
                        pDropout;
                        wgt = wgt,
                        weightInit = params.weightInit
                       ))
    nin = nout

    # next layers will get the "hidden" dropout probability
    pDropout = getDropoutProb(params, false)
  end

  # push the output layer
  push!(layers, LAYER(nin,
                      numOutputs,
                      finalMapping,
                      params.updater,
                      pDropout;
                      wgt = wgt,
                      weightInit = params.weightInit
                     ))

  NeuralNet{LAYER}(layers, params, solverParams, inputTransformer)
end

buildClassificationNet(args...; kwargs...) = buildNet(args...; kwargs..., finalMapping = SigmoidMapping())
buildTanhClassificationNet(args...; kwargs...) = buildNet(args...; kwargs..., hiddenMapping = TanhMapping(), finalMapping = TanhMapping())
buildRegressionNet(args...; kwargs...) = buildNet(args...; kwargs..., finalMapping = IdentityMapping())

buildNet{LAYER <: NeuralNetLayer}(::Type{LAYER}, args...; kw...) = buildNet(args...; kw..., LAYER = LAYER)


text:
import Pkg
Pkg.activate("./")

using PDMesh

obj = Cuboid([-10 10; 0 3; -2 2])

f = out -> begin
x=out[1];
mask = (x[1, :] .< 0.1) .& (x[1, :] .> -0.1) .& (x[2, :] .> 2)
mask
end

obj = delete(obj,  f)

x, v, y, vol, type = create(obj, resolution=0.1, rand_=0.0, type=1)

write_data("./examples/ovito_ref.data", x, type, vol)

write_data("./examples/ovito.data", y, v, type, vol)

max_stars_repo_name: ravinderbhattoo/PDMesh

max_issues_repo_name: ravinderbhattoo/PDMesh

max_forks_repo_name: ravinderbhattoo/PDMesh


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.eyeM)),NamedTuple{(:like,),Tuple{Array{Int8,2}}},typeof(eyeM),Int64,Int64})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.freqspaceM)),NamedTuple{(:dim,),Tuple{Int64}},typeof(freqspaceM),Int64})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.ndgridM)),NamedTuple{(:dim,),Tuple{Int64}},typeof(ndgridM),StepRange{Int64,Int64}})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.onesM)),NamedTuple{(:like,),Tuple{Array{Float64,2}}},typeof(onesM),Array{Int64,1}})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.onesM)),NamedTuple{(:like,),Tuple{Array{Int32,2}}},typeof(onesM),Symbol,Int64})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.onesM)),NamedTuple{(:like,),Tuple{Array{Int64,2}}},typeof(onesM),Int64})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.randM)),NamedTuple{(:like,),Tuple{Array{Float64,2}}},typeof(randM),Array{Int64,1}})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.randM)),NamedTuple{(:like,),Tuple{Array{Int32,2}}},typeof(randM),Symbol,Int64})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.randM)),NamedTuple{(:like,),Tuple{Array{Int64,2}}},typeof(randM),Int64})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.sortM)),NamedTuple{(:I,),Tuple{Bool}},typeof(sortM),Array{Int64,1}})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.sortM)),NamedTuple{(:I,),Tuple{Bool}},typeof(sortM),Array{Int64,2},Int64,Symbol})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.sortperm)),NamedTuple{(:dims,),Tuple{Int64}},typeof(sortperm),Array{Int64,2}})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.zerosM)),NamedTuple{(:like,),Tuple{Array{Float64,2}}},typeof(zerosM),Array{Int64,1}})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.zerosM)),NamedTuple{(:like,),Tuple{Array{Int32,2}}},typeof(zerosM),Symbol,Int64})
    Base.precompile(Tuple{Core.kwftype(typeof(MatLang.zerosM)),NamedTuple{(:like,),Tuple{Array{Int64,2}}},typeof(zerosM),Int64})
    Base.precompile(Tuple{typeof(blkdiagM),Array{Float64,2},Array{Float64,2},Vararg{Array{Float64,2},N} where N})
    Base.precompile(Tuple{typeof(blkdiagM),Symbol,Array{Float64,2},Vararg{Array{Float64,2},N} where N})
    Base.precompile(Tuple{typeof(catM),Int64,Array{Float64,2},Vararg{Array{Float64,2},N} where N})
    Base.precompile(Tuple{typeof(catM),Int64,Array{Float64,3},Vararg{Array{Float64,3},N} where N})
    Base.precompile(Tuple{typeof(clcM)})
    Base.precompile(Tuple{typeof(diagM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(diagM),Array{Int64,2}})
    Base.precompile(Tuple{typeof(doubleM),Float64})
    Base.precompile(Tuple{typeof(eyeM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(eyeM),Int64,Int64})
    Base.precompile(Tuple{typeof(eyeM),Int64})
    Base.precompile(Tuple{typeof(eyeM),Tuple{Int64,Int64}})
    Base.precompile(Tuple{typeof(eyeM),Type,Array{Int64,1}})
    Base.precompile(Tuple{typeof(eyeM),Type,Int64,Int64})
    Base.precompile(Tuple{typeof(eyeM),Type,Tuple{Int64,Int64}})
    Base.precompile(Tuple{typeof(falseM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(falseM),Int64})
    Base.precompile(Tuple{typeof(falseM),Symbol,Int64})
    Base.precompile(Tuple{typeof(falseM),Tuple{Int64,Int64}})
    Base.precompile(Tuple{typeof(flipM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(flipM),Array{Int64,2}})
    Base.precompile(Tuple{typeof(flipM),String})
    Base.precompile(Tuple{typeof(int16M),Int64})
    Base.precompile(Tuple{typeof(int32M),Int64})
    Base.precompile(Tuple{typeof(int64M),Int64})
    Base.precompile(Tuple{typeof(int8M),Int64})
    Base.precompile(Tuple{typeof(iscolumnM),Array{Float64,1}})
    Base.precompile(Tuple{typeof(iscolumnM),Array{Float64,2}})
    Base.precompile(Tuple{typeof(iscolumnM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(iscolumnM),Array{Int64,2}})
    Base.precompile(Tuple{typeof(iscolumnM),Array{String,1}})
    Base.precompile(Tuple{typeof(iscolumnM),Array{String,2}})
    Base.precompile(Tuple{typeof(iscolumnM),Int64})
    Base.precompile(Tuple{typeof(iscolumnM),String})
    Base.precompile(Tuple{typeof(iscolumnM),Symbol,Array{Float64,1}})
    Base.precompile(Tuple{typeof(iscolumnM),Symbol,Array{Float64,2}})
    Base.precompile(Tuple{typeof(iscolumnM),Symbol,Array{Int64,1}})
    Base.precompile(Tuple{typeof(iscolumnM),Symbol,Array{Int64,2}})
    Base.precompile(Tuple{typeof(iscolumnM),Symbol,Array{String,1}})
    Base.precompile(Tuple{typeof(iscolumnM),Symbol,Int64})
    Base.precompile(Tuple{typeof(iscolumnM),Symbol,String})
    Base.precompile(Tuple{typeof(ismatrixM),Array{Float64,1}})
    Base.precompile(Tuple{typeof(ismatrixM),Array{Float64,2}})
    Base.precompile(Tuple{typeof(ismatrixM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(ismatrixM),Array{Int64,2}})
    Base.precompile(Tuple{typeof(ismatrixM),Array{String,1}})
    Base.precompile(Tuple{typeof(ismatrixM),Array{String,2}})
    Base.precompile(Tuple{typeof(ismatrixM),String})
    Base.precompile(Tuple{typeof(ismatrixM),Symbol,Array{Float64,1}})
    Base.precompile(Tuple{typeof(ismatrixM),Symbol,Array{Float64,2}})
    Base.precompile(Tuple{typeof(ismatrixM),Symbol,Array{Int64,1}})
    Base.precompile(Tuple{typeof(ismatrixM),Symbol,Array{Int64,2}})
    Base.precompile(Tuple{typeof(ismatrixM),Symbol,Array{String,1}})
    Base.precompile(Tuple{typeof(ismatrixM),Symbol,String})
    Base.precompile(Tuple{typeof(isrowM),Array{Float64,1}})
    Base.precompile(Tuple{typeof(isrowM),Array{Float64,2}})
    Base.precompile(Tuple{typeof(isrowM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(isrowM),Array{Int64,2}})
    Base.precompile(Tuple{typeof(isrowM),Array{String,1}})
    Base.precompile(Tuple{typeof(isrowM),Array{String,2}})
    Base.precompile(Tuple{typeof(isrowM),Int64})
    Base.precompile(Tuple{typeof(isrowM),String})
    Base.precompile(Tuple{typeof(isrowM),Symbol,Array{Float64,1}})
    Base.precompile(Tuple{typeof(isrowM),Symbol,Array{Float64,2}})
    Base.precompile(Tuple{typeof(isrowM),Symbol,Array{Int64,1}})
    Base.precompile(Tuple{typeof(isrowM),Symbol,Array{Int64,2}})
    Base.precompile(Tuple{typeof(isrowM),Symbol,Array{String,1}})
    Base.precompile(Tuple{typeof(isrowM),Symbol,Int64})
    Base.precompile(Tuple{typeof(isrowM),Symbol,String})
    Base.precompile(Tuple{typeof(isscalarM),Array{Float64,3}})
    Base.precompile(Tuple{typeof(isscalarM),Array{String,1}})
    Base.precompile(Tuple{typeof(isscalarM),Array{String,2}})
    Base.precompile(Tuple{typeof(isscalarM),Bool})
    Base.precompile(Tuple{typeof(isscalarM),Char})
    Base.precompile(Tuple{typeof(isscalarM),Int64})
    Base.precompile(Tuple{typeof(isscalarM),Symbol,Array{Float64,3}})
    Base.precompile(Tuple{typeof(isscalarM),Symbol,Array{Int64,1}})
    Base.precompile(Tuple{typeof(isscalarM),Symbol,Array{String,1}})
    Base.precompile(Tuple{typeof(isscalarM),Symbol,Array{String,2}})
    Base.precompile(Tuple{typeof(isscalarM),Symbol,Bool})
    Base.precompile(Tuple{typeof(isscalarM),Symbol,Char})
    Base.precompile(Tuple{typeof(isscalarM),Symbol,Int64})
    Base.precompile(Tuple{typeof(isscalarM),Symbol,String})
    Base.precompile(Tuple{typeof(isvectorM),Array{Float64,1}})
    Base.precompile(Tuple{typeof(isvectorM),Array{Float64,2}})
    Base.precompile(Tuple{typeof(isvectorM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(isvectorM),Array{Int64,2}})
    Base.precompile(Tuple{typeof(isvectorM),Array{String,1}})
    Base.precompile(Tuple{typeof(isvectorM),Array{String,2}})
    Base.precompile(Tuple{typeof(isvectorM),Int64})
    Base.precompile(Tuple{typeof(isvectorM),String})
    Base.precompile(Tuple{typeof(isvectorM),Symbol,Array{Float64,1}})
    Base.precompile(Tuple{typeof(isvectorM),Symbol,Array{Float64,2}})
    Base.precompile(Tuple{typeof(isvectorM),Symbol,Array{Int64,1}})
    Base.precompile(Tuple{typeof(isvectorM),Symbol,Array{Int64,2}})
    Base.precompile(Tuple{typeof(isvectorM),Symbol,Array{String,1}})
    Base.precompile(Tuple{typeof(isvectorM),Symbol,Int64})
    Base.precompile(Tuple{typeof(isvectorM),Symbol,String})
    Base.precompile(Tuple{typeof(lengthM),Array{Any,2}})
    Base.precompile(Tuple{typeof(lengthM),Array{Int64,2}})
    Base.precompile(Tuple{typeof(linspaceM),Int64,Int64,Int64})
    Base.precompile(Tuple{typeof(linspaceM),Int64,Int64})
    Base.precompile(Tuple{typeof(linspaceM),Symbol,Int64,Int64,Int64})
    Base.precompile(Tuple{typeof(linspaceM),Symbol,Int64,Int64})
    Base.precompile(Tuple{typeof(logspaceM),Int64,Int64,Int64})
    Base.precompile(Tuple{typeof(logspaceM),Int64,Int64,Symbol})
    Base.precompile(Tuple{typeof(logspaceM),Int64,Int64})
    Base.precompile(Tuple{typeof(meshgridM),StepRange{Int64,Int64}})
    Base.precompile(Tuple{typeof(meshgridM),UnitRange{Int64},UnitRange{Int64},UnitRange{Int64}})
    Base.precompile(Tuple{typeof(ndgridM),UnitRange{Int64},UnitRange{Int64},UnitRange{Int64}})
    Base.precompile(Tuple{typeof(onesM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(onesM),Int64})
    Base.precompile(Tuple{typeof(onesM),Symbol,Int64})
    Base.precompile(Tuple{typeof(onesM),Tuple{Int64,Int64}})
    Base.precompile(Tuple{typeof(onesM),Type,Int64,Int64})
    Base.precompile(Tuple{typeof(onesM),Type,Tuple{Int64,Int64}})
    Base.precompile(Tuple{typeof(onesM),Type{Int32},Array{Int64,1}})
    Base.precompile(Tuple{typeof(randM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(randM),Int64})
    Base.precompile(Tuple{typeof(randM),Symbol,Int64})
    Base.precompile(Tuple{typeof(randM),Tuple{Int64,Int64}})
    Base.precompile(Tuple{typeof(randM),Type,Int64,Int64})
    Base.precompile(Tuple{typeof(randM),Type,Tuple{Int64,Int64}})
    Base.precompile(Tuple{typeof(randM),Type{Int32},Array{Int64,1}})
    Base.precompile(Tuple{typeof(repelemM),Array{Int64,2},Int64})
    Base.precompile(Tuple{typeof(repmatM),Array{Int64,1},Array{Int64,1}})
    Base.precompile(Tuple{typeof(repmatM),Array{Int64,1},Int64,Int64})
    Base.precompile(Tuple{typeof(repmatM),Array{Int64,1},Int64})
    Base.precompile(Tuple{typeof(repmatM),Array{Int64,1},Tuple{Int64,Int64}})
    Base.precompile(Tuple{typeof(repmatM),Int64,Int64,Int64})
    Base.precompile(Tuple{typeof(singleM),Float64})
    Base.precompile(Tuple{typeof(sizeM),Array{Float64,3},Int64,Int64})
    Base.precompile(Tuple{typeof(sizeM),Array{Int64,2},Int64})
    Base.precompile(Tuple{typeof(sizeM),Symbol,Array{Float64,3},Int64,Int64})
    Base.precompile(Tuple{typeof(sizeM),Symbol,Array{Int64,2},Int64})
    Base.precompile(Tuple{typeof(sizeM),Symbol,Array{Int64,2}})
    Base.precompile(Tuple{typeof(sortM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(sortM),Array{Int64,2}})
    Base.precompile(Tuple{typeof(sortM),Array{Integer,3},Int64})
    Base.precompile(Tuple{typeof(squeezeM),Array{Float64,3}})
    Base.precompile(Tuple{typeof(transposeM),Array{Int64,2}})
    Base.precompile(Tuple{typeof(transposeM),Symbol,Array{Int64,2}})
    Base.precompile(Tuple{typeof(trueM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(trueM),Int64})
    Base.precompile(Tuple{typeof(trueM),Symbol,Int64})
    Base.precompile(Tuple{typeof(trueM),Tuple{Int64,Int64}})
    Base.precompile(Tuple{typeof(uint16M),Int64})
    Base.precompile(Tuple{typeof(uint32M),Int64})
    Base.precompile(Tuple{typeof(uint64M),Int64})
    Base.precompile(Tuple{typeof(uint8M),Int64})
    Base.precompile(Tuple{typeof(zerosM),Array{Int64,1}})
    Base.precompile(Tuple{typeof(zerosM),Int64})
    Base.precompile(Tuple{typeof(zerosM),Symbol,Int64})
    Base.precompile(Tuple{typeof(zerosM),Tuple{Int64,Int64}})
    Base.precompile(Tuple{typeof(zerosM),Type,Int64,Int64,Int64})
    Base.precompile(Tuple{typeof(zerosM),Type,Int64,Int64})
    Base.precompile(Tuple{typeof(zerosM),Type,Tuple{Int64,Int64}})
    Base.precompile(Tuple{typeof(zerosM),Type{Int32},Array{Int64,1}})
    isdefined(MatLang, Symbol("##32#33")) && Base.precompile(Tuple{getfield(MatLang, Symbol("##32#33")),Array{Int64,1}})
    isdefined(MatLang, Symbol("##sortperm#31")) && Base.precompile(Tuple{getfield(MatLang, Symbol("##sortperm#31")),Int64,Bool,typeof(sortperm),Array{Int64,2}})
end

max_stars_repo_name: juliamatlab/MatLang

max_issues_repo_name: juliamatlab/MatLang

max_forks_repo_name: juliamatlab/MatLang


text:
#################### Posterior Statistics ####################
"""
    autocor(c::AbstractChains; lags::Vector=[1, 5, 10, 50],
             relative::Bool=true)

Compute lag-k autocorrelations for MCMC sampler output.

Returns a `ChainSummary` type object with model parameters indexed by the first dimension of `value`, lag-autocorrelations by the second, and chains by the third.

* `c` : sampler output on which to perform calculations.

* `lags` : lags at which to compute autocorrelations.

* `relative` : whether the lags are relative to the thinning interval of the output (`true`) or relative to the absolute iteration numbers (`false`).
"""
function autocor(c::AbstractChains; lags::Vector = [1, 5, 10, 50], relative::Bool = true)
    if relative
        lags *= step(c)
    elseif any(lags .% step(c) .!= 0)
        throw(ArgumentError("lags do not correspond to thinning interval"))
    end
    labels = map(x -> "Lag " * string(x), lags)
    vals = mapslices(x -> autocor(x, lags)', c.value, dims = [1, 2])
    ChainSummary(vals, c.names, labels, header(c))
end
"""
    cor(c::AbstractChains)

Compute cross-correlations for MCMC sampler output.

Returns a `ChainSummary` type object with the first and second dimensions of the `value` field indexing the model parameters between which correlations. Results are for all chains combined.

* `c` : sampler output on which to perform calculations.
"""
function cor(c::AbstractChains)
    ChainSummary(cor(combine(c)), c.names, c.names, header(c))
end
"""
    changerate(c::AbstractChains)

Estimate the probability, or rate per iteration, ``\\Pr(\\theta^i \\ne \\theta^{i-1})`` of a state space change for iterations ``i = 2, \\ldots, N`` in MCMC sampler output. Estimation is performed for each parameter univariately as well as for the full parameter vector multivariately. For continuous output generated from samplers, like Metropolis-Hastings, whose algorithms conditionally accept candidate draws, the probability can be viewed as the acceptance rate.

Returns a `ChainSummary` type object with parameters in the rows of the `value` field, and the estimated rates in the column. Results are for all chains combined.

* `c` : sampler output on which to perform calculations.
"""
function changerate(c::AbstractChains)
    n, p, m = size(c.value)
    r = zeros(Float64, p, 1, 1)
    r_mv = 0.0
    delta = Array{Bool}(undef, p)
    for k = 1:m
        prev = c.value[1, :, k]
        for i = 2:n
            for j = 1:p
                x = c.value[i, j, k]
                dx = x != prev[j]
                r[j] += dx
                delta[j] = dx
                prev[j] = x
            end
            r_mv += any(delta)
        end
    end
    vals = round.([r; r_mv] / (m * (n - 1)), digits = 3)
    ChainSummary(vals, [c.names; "Multivariate"], ["Change Rate"], header(c))
end

describe(c::AbstractChains; args...) = describe(stdout, c; args...)
"""
    describe(io::IO, c::AbstractChains;
                  q::Vector=[0.025, 0.25, 0.5, 0.75, 0.975], etype=:bm, args...)

Compute summary statistics for MCMC sampler output.

Returns results from calls to `summarystats(c, etype, args...)` and `quantile(c, q)` are printed for all chains combined, and a value of `nothing` is returned.

* `c` : sampler output on which to perform calculations.

* `q` : probabilities at which to calculate quantiles.

* `etype` : method for computing Monte Carlo standard errors. See `mcse()` for options.

* `args...` : additional arguments to be passed to the `etype` method.
"""
function describe(
    io::IO,
    c::AbstractChains;
    q::Vector = [0.025, 0.25, 0.5, 0.75, 0.975],
    etype = :bm,
    args...,
)
    ps_stats = summarystats(c; etype = etype, args...)
    ps_quantiles = quantile(c, q = q)
    println(io, ps_stats.header)
    print(io, "Empirical Posterior Estimates:\n")
    show(io, ps_stats)
    print(io, "Quantiles:\n")
    show(io, ps_quantiles)
end
"""
    hpd(x::Vector{T}; alpha::Real=0.05) where {T<:Real}

Compute highest posterior density (HPD) intervals of Chen and Shao [16] for MCMC sampler output. HPD intervals have the desirable property of being the smallest intervals that contain a given probability. However, their calculation assumes unimodal marginal posterior distributions, and they are not invariant to transformations of parameters like central (quantile-based) posterior intervals.

Returns a `ChainSummary` type object with parameters contained in the rows of the value field, and lower and upper intervals in the first and second columns. Results are for all chains combined.

* `c` : sampler output on which to perform calculations.

* `alpha` : the `100 * (1 - alpha)`% interval to compute.

* `x` : vector on which to perform calculations.
"""
function hpd(x::Vector{T}; alpha::Real = 0.05) where {T<:Real}
    n = length(x)
    m = max(1, ceil(Int, alpha * n))

    y = sort(x)
    a = y[1:m]
    b = y[(n-m+1):n]
    _, i = findmin(b - a)

    [a[i], b[i]]
end
"""
    hpd(c::AbstractChains; alpha::Real=0.05)
"""
function hpd(c::AbstractChains; alpha::Real = 0.05)
    pct = first(showoff([100.0 * (1.0 - alpha)]))
    labels = ["$(pct)% Lower", "$(pct)% Upper"]
    vals = permutedims(
        mapslices(x -> hpd(vec(x), alpha = alpha), c.value, dims = [1, 3]),
        [2, 1, 3],
    )
    ChainSummary(vals, c.names, labels, header(c))
end
"""
    quantile(c::AbstractChains; q::Vector=[0.025, 0.25, 0.5, 0.75, 0.975])

Compute posterior quantiles for MCMC sampler output.

Returns a `ChainSummary` type object with parameters contained in the rows of the `value` field, and quantiles in the columns. Results are for all chains combined.

* `c` : sampler output on which to perform calculations.

* `q` : probabilities at which to compute quantiles.
"""
function quantile(c::AbstractChains; q::Vector = [0.025, 0.25, 0.5, 0.75, 0.975])
    labels = map(x -> string(100 * x) * "%", q)
    vals =
        permutedims(mapslices(x -> quantile(vec(x), q), c.value, dims = [1, 3]), [2, 1, 3])
    ChainSummary(vals, c.names, labels, header(c))
end
"""
    summarystats(c::AbstractChains; etype=:bm, args...)

Compute posterior summary statistics for MCMC sampler output.

Returns a `ChainSummary` type object with parameters in the rows of the `value` field; and the sample mean, standard deviation, standard error, Monte Carlo standard error, and effective sample size in the columns. Results are for all chains combined.

* `c` : sampler output on which to perform calculations.

* `etype` : method for computing Monte Carlo standard errors. See `mcse()` for options.

* `args...` : additional arguments to be passed to the `etype` method.
"""
function summarystats(c::AbstractChains; etype = :bm, args...)
    f = x -> [mean(x), std(x), sem(x), mcse(vec(x), etype; args...)]
    labels = ["Mean", "SD", "Naive SE", "MCSE", "ESS"]
    vals = permutedims(mapslices(x -> f(x), c.value, dims = [1, 3]), [2, 1, 3])
    stats = [vals min.((vals[:, 2] ./ vals[:, 4]) .^ 2, size(c.value, 1))]
    ChainSummary(stats, c.names, labels, header(c))
end

max_stars_repo_name: erathorn/MCPhylo

max_issues_repo_name: erathorn/MCPhylo

max_forks_repo_name: erathorn/MCPhylo


text:
using Quantica: nvertices, nedges, nsimplices

@test begin
    mesh = marchingmesh(0:0.1:1, 1:0.1:2)
    nvertices(mesh) == 121 && nedges(mesh) == 320 && nsimplices(mesh) == 200
end


text:
# irflp.m
# IRF Estimates and stdev of IRF by Local Projections
# Luca Brugnolini 2014
# Based on Kilian

function irflp(z, paic, H, A0inv, COVsig)

# Input : data z, lag order paic (Hx1 vector), horizon for IRF H,
#         A0inv- relationship btw sf & rf shocks estimated from VAR
#         COVsig- variance due to A0inv, estimated from VAR
# Output: IRF, Std error of IRF, reduced form IRF

(T,K) = size(z);

# IRF estimate and STdev of IRF for Horizon 0
A1 = vec(A0inv');                                # IRF for horizon 0
A2 = vec(eye(K));                                # Reduced form IRF
COVA1 = zeros(K^2,1);                            # Std(IRF0)=0

# IRF estimate and Stdev of IRF for Horizon 1~H
 for h = 1:H
     ph = paic[h];                     # lag-order for horizon h
     ph = itrunc(ph);
     ys = z[ph+h:T,:];

    # RHS variables: lags
     yt = z[ph:T-h,:];                 # RHS variable of interest: y(t)
     x = ones(T-ph-h+1,1);             # constant term

        for i = 2:ph
            x = [x z[ph+1-i:T-i-h+1,:]];  # other RHS variables: y(t-1)~y(t-p+1)

        end

     Mx = eye(size(x,1))-x/(x'*x)*x';         # annhilation matrix
     B1 = (yt'*Mx*yt)\(yt'*Mx*ys);            # IRF by Local Projection
     A1 = [A1 vec(A0inv'*B1)];                # structural IRF
     A2 = [A2 vec(B1)];                       # reduced form IRF

     # Newey-West HAC Estimator of Stdev (truncation point is h)
     u = Mx*ys - Mx*yt*B1;                                           # residual from LPs
     meanu = [mean(u[:,1]) mean(u[:,2]) mean(u[:,3]) mean(u[:,4])];
     u0 = u-kron(ones(size(u,1),1),meanu);
     R0 = (u0'*u0)./(size(u,1));         # here was length(u)                            # R(0): smpl cov(u)
     SIGu = R0;
     M = h;                                                          # truncation point

        for j = 1:M
            Rj = (u0[1:end-j,:]'*u0[j+1:end,:])./(size(u,1));   # R(j)
            SIGu = SIGu + (1-j/(M+1))*(Rj+Rj');                 # varcov(u)

        end

     invytMxyt  = inv(yt'*Mx*yt);
     SIGB1 = kron(SIGu, invytMxyt);                               # var(vec(B1))
     SIGA1 = kron(eye(K), A0inv')*SIGB1*kron(eye(K), A0inv')';     # var(vec(A0inv*B1))
     COVA1 = [COVA1 reshape(diag(SIGA1), K^2, 1)];

end

    stdA1 = sqrt(COVA1+COVsig);                                       # STD2 is variance due to

     return A1, stdA1, A2

end


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

module TLEAliases

using ..TLE
using ..Empty
using ..SK
using ..EA
using ..SAT

export Graph0TLE, GraphSKTLE, GraphEATLE, GraphSATTLE

const Graph0TLE{M,γT,λT} = GraphTopologicalLocalEntropy{M,γT,λT,GraphEmpty}

# """
#     Graph0TLE(...)
#
# TODO
# """
Graph0TLE(Nk::Integer, M::Integer, γ::Float64, λ::Float64, β::Float64) = GraphTopologicalLocalEntropy(Nk, M, γ, λ, β, GraphEmpty, Nk)

const GraphSKTLE{M,γT,λT} = GraphTopologicalLocalEntropy{M,γT,λT,GraphSK}

# """
#     GraphSKTLE(N::Integer, M::Integer, γ::Float64, λ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
GraphSKTLE(Nk::Integer, M::Integer, γ::Float64, λ::Float64, β::Float64) = GraphTopologicalLocalEntropy(Nk, M, γ, λ, β, GraphSK, SK.gen_J(Nk))

const GraphEATLE{M,γT,λT,twoD} = GraphTopologicalLocalEntropy{M,γT,λT,GraphEANormal{twoD}}

# """
#     GraphEATLE(L::Integer, D::Integer, M::Integer, γ::Float64, λ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
function GraphEATLE(L::Integer, D::Integer, M::Integer, γ::Float64, λ::Float64, β::Float64)
    D ≥ 1 || throw(ArgumentError("D must be ≥ 0, given: $D"))
    A = EA.gen_EA(L, D)
    N = length(A)
    J = EA.gen_J(Float64, N, A) do
        4*rand() - 2
    end

    GraphTopologicalLocalEntropy(N, M, γ, λ, β, GraphEANormal{2D}, L, A, J)
end

function GraphEATLE(fname::AbstractString, M::Integer, γ::Float64, λ::Float64, β::Float64)
    L, D, A, J = EA.gen_AJ(fname)
    N = length(A)
    @assert N == L^D
    GraphTopologicalLocalEntropy(N, M, γ, λ, β, GraphEANormal{2D}, L, A, J)
end

function GraphEATLE(X::GraphEANormal{twoD}, M::Integer, γ::Float64, λ::Float64, β::Float64) where {twoD}
    @assert iseven(twoD)
    D = twoD ÷ 2
    N = X.N
    L = round(Int, N^(1/D))
    @assert L^D == N
    GraphTopologicalLocalEntropy(N, M, γ, λ, β, GraphEANormal{twoD}, L, X.A, X.J)
end

const GraphSATTLE{M,γT,λT} = GraphTopologicalLocalEntropy{M,γT,λT,GraphSAT}

# """
#     GraphSATTLE(N::Integer, K::Integer, α::Real, M::Integer, γ::Float64, λ::Float64, β::Float64) <: DoubleGraph
#
# TODO
# """
function GraphSATTLE(N::Integer, K::Integer, α::Real, M::Integer, γ::Float64, λ::Float64, β::Float64)
    A, J = SAT.gen_randomKSAT(N, K, α)
    GraphTopologicalLocalEntropy(N, M, γ, λ, β, GraphSAT, N, A, J)
end

function GraphSATTLE(X::GraphSAT, M::Integer, γ::Float64, λ::Float64, β::Float64)
    N = X.N
    GraphTopologicalLocalEntropy(N, M, γ, λ, β, GraphSAT, N, X.A, X.J)
end

end # module


text:
function assert_allow_step(state)
    if state.broke_on_error
        printstyled(stderr, "Cannot step after breaking on error\n"; color=:red)
        return false
    end
    if state.level != 1
        printstyled(stderr, "Cannot step in a non leaf frame\n"; color=:red)
        return false
    end
    return true
end

function show_breakpoint(io::IO, bp::BreakpointRef)
    outbuf = IOContext(IOBuffer(), io)
    if bp.err === nothing
        print(outbuf, "Hit breakpoint: ")
    else
        print(outbuf, "Breaking on error: ")
    end
    if checkbounds(Bool, bp.framecode.breakpoints, bp.stmtidx)
        lineno = linenumber(bp.framecode, bp.stmtidx)
        print(outbuf, bp.framecode.scope, ", line ", lineno)
    else
        print(outbuf, bp.framecode.scope, ", %", bp.stmtidx)
    end
    if bp.err !== nothing
        print(outbuf, ", ", bp.err)
    end
    print(io, String(take!(outbuf.io)))
    println(io)
end

function execute_command(state::DebuggerState, ::Union{Val{:c},Val{:nc},Val{:n},Val{:se},Val{:s},Val{:si},Val{:sg},Val{:so}}, cmd::AbstractString)
    assert_allow_step(state) || return false
    cmd == "so" && (cmd = "finish")
    ret = debug_command(state.frame, cmd)
    if ret === nothing
        state.overall_result = get_return(root(state.frame))
        state.frame = nothing
        return false
    else
        state.frame, pc = ret
        if pc isa BreakpointRef
            if pc.stmtidx != 0 # This is the dummy breakpoint to stop just after entering a call
                if state.terminal !== nothing # fix this, it happens when a test hits this and hasnt set a terminal
                    show_breakpoint(Base.pipe_writer(state.terminal), pc)
                end
            end
            if pc.err !== nothing
                state.broke_on_error = true
            end
        end
        JuliaInterpreter.maybe_next_call!(state.frame)
        return true
    end
end

function execute_command(state::DebuggerState, ::Val{:bt}, cmd)
    num = 0
    frame = state.frame
    while frame !== nothing
        num += 1
        print_frame(Base.pipe_writer(state.terminal), num, frame)
        frame = caller(frame)
    end
    println()
    return false
end

function stacklength(frame)
    s = 0
    JuliaInterpreter.traverse(fr -> (s += 1; JuliaInterpreter.caller(fr)), JuliaInterpreter.leaf(frame))
    return s
end

function execute_command(state::DebuggerState, ::Union{Val{:f}, Val{:fr}}, cmd)
    subcmds = split(cmd,' ')[2:end]
    if isempty(subcmds) || subcmds[1] == "v"
        print_frame(Base.pipe_writer(state.terminal), state.level, active_frame(state))
        return false
    else
        new_level = parse(Int, subcmds[1])
        if new_level > stacklength(state.frame) || new_level < 1
            printstyled(stderr, "Not a valid frame index\n"; color=:red)
            return false
        end
        state.level = new_level
        return true
    end
end

function execute_command(state::DebuggerState, _, cmd)
    println("Unknown command `$cmd`. Executing `?` to obtain help.")
    execute_command(state, Val{Symbol("?")}(), "?")
end

function execute_command(state::DebuggerState, ::Val{:?}, cmd::AbstractString)
    display(
            @md_str """
    Basic Commands:\\
    - `n`: step to the next line\\
    - `s`: step into the next call\\
    - `so`: step out of the current call\\
    - `c`: continue execution until a breakpoint is hit\\
    - `bt`: show a simple backtrace\\
    - ``` `stuff ```: run `stuff` in the current function's context\\
    - `fr [v::Int]`: show all variables in the current frame, `v` defaults to `1`\\
    - `f [n::Int]`: go to the `n`-th frame\\
    - `q`: quit the debugger, returning `nothing`\\
    Advanced commands:\\
    - `nc`: step to the next call\\
    - `se`: step one expression step\\
    - `si`: same as `se` but step into a call if a call is the next expression\\
    - `sg`: step into a generated function\\
    """)
    return false
end


text:
function get_gm(obj::DSSTNewtonianAttractionContext)
    return jcall(obj, "getGM", jdouble, ())
end


text:
function cuMPS(psi::MPS)
    phi = copy(psi)
    for site in 1:length(psi)
        phi.data[site] = cuITensor(psi.data[site])
    end
    return phi
end

cu(ψ::MPS) = cuMPS(ψ)

cuMPS() = MPS()
function cuMPS(sites) # random MPS
  N = length(sites)
  v = Vector{ITensor}(undef, N)
  l = [Index(1, "Link,l=$ii") for ii=1:N-1]
  for ii in eachindex(sites)
    s = sites[ii]
    if ii == 1
      v[ii] = cuITensor(l[ii], s)
    elseif ii == N
      v[ii] = cuITensor(l[ii-1], s)
    else
      v[ii] = cuITensor(l[ii-1],s,l[ii])
    end
  end
  return MPS(v,0,N+1)
end

function randomCuMPS(sites)
  M = cuMPS(sites)
  for i in eachindex(sites)
    randn!(M[i])
    normalize!(M[i])
  end
  M.llim = 1
  M.rlim = length(M)
  return M
end

function productCuMPS(::Type{T}, ivals::Vector{<:IndexVal}) where {T<:Number}
  N     = length(ivals)
  As    = Vector{ITensor}(undef,N)
  links = Vector{Index}(undef,N)
  for n=1:N
    s = ITensors.ind(ivals[n])
    links[n] = Index(1,"Link,l=$n")
    if n == 1
      A = ITensor(T, s, links[n])
      A[ivals[n],links[n](1)] = 1.0
    elseif n == N
      A = ITensor(T, links[n-1], s)
      A[links[n-1](1),ivals[n]] = 1.0
    else
      A = ITensor(T, links[n-1], s, links[n])
      A[links[n-1](1),ivals[n],links[n](1)] = 1.0
    end
    As[n] = cuITensor(A)
  end
  return MPS(As,0,2)
end
productCuMPS(ivals::Vector{<:IndexVal}) = productCuMPS(Float64, ivals)

function productCuMPS(::Type{T}, sites,
                      states) where {T<:Number}
  if length(sites) != length(states)
    throw(DimensionMismatch("Number of sites and and initial states don't match"))
  end
  ivals = [state(sites[n],states[n]) for n=1:length(sites)]
  return productCuMPS(T, ivals)
end

productCuMPS(sites, states) = productCuMPS(Float64, sites, states)


text:
using JuMP
using Clp
using Cbc
using GLPK

include("../OptMaps.jl")
using .OptMaps

function solve_model(solver, model_log, solver_options=Dict())
    model = Model(solver)

    for (k, v) in solver_options
        set_optimizer_attribute(model, k, v)
    end

    #sets
    xs = [1, 2, 3, 4]
    ys = [:skilled, :unskilled]
    resources = [:mach1, :mach2]

    #parameters
    unit_prices = [300 260 220 180]
    labor_cost_hour = [8 6]

    rsc_per_prod = [11 7 6 5;
                    4 6 5 4]

    hours_per_prod = [8 5 5 6;
                      7 8 7 4]

    resource_limits = [700 500]
    hour_limits = [600 650]

    n_products = length(xs)
    n_labor_types = length(ys)
    n_resources = length(resources)

    n_variables = n_products + n_labor_types
    n_constraints = n_resources

    ###############################################################
    #decision variables
    @variable(model, x[1:n_products] >= 0)
    @variable(model, y[1:n_labor_types] >= 0)

    starting_point = zeros(Float64, n_variables)
    solution_point = zeros(Float64, n_variables)

    ###############################################################
    #objective

    #regular function
    total_profit_f(v) = sum([v[i] * unit_prices[i] for i in 1:n_products]) - sum(v[l+4] * labor_cost_hour[l] for l in 1:n_labor_types)

    #jump expression
    total_profit = @expression(model, total_profit, sum(x[p] * unit_prices[p] for p in 1:n_products) - sum(y[l] * labor_cost_hour[l] for l in 1:n_labor_types))

    @objective(model, Max, total_profit)

    ###############################################################
    #constraints:

    #regular functions
    constraints_f = []

    for r in 1:n_resources
        push!(constraints_f, (v) -> sum([v[p] * rsc_per_prod[r,p] for p in 1:n_products]) <= resource_limits[r])
    end

    for l in 1:n_labor_types
        push!(constraints_f, (v) -> sum([v[p] * hours_per_prod[l,p] for p in 1:n_products]) <= v[l+4])
    end

    for l in 1:n_labor_types
        push!(constraints_f, (v) -> v[l+4] <= hour_limits[l])
    end

    #jump expressions
    @constraint(model, rsc[r = 1:n_resources], sum((x[p]*rsc_per_prod[r,p]  for p in 1:n_products)) <= resource_limits[r])
    @constraint(model, hours[l = 1:n_labor_types], sum((x[p]*hours_per_prod[l,p]  for p in 1:n_products)) <= y[l])
    @constraint(model, hours_lim[l = 1:n_labor_types], y[l] <= hour_limits[l])

    ###############################################################
    #save model textual description
    open(model_log, "w") do f
        print(f, model)

        ###############################################################
        #solve model
        elapsed_time = @elapsed optimize!(model)

        println(f, "-----------------------------------------------")
        println(f, "Elapsed time: $elapsed_time")
        println(f, "-----------------------------------------------")
        println(f, "Status:    $(termination_status(model))")
        println(f, "Objective: $(objective_value(model))")
        println(f, "-----------------------------------------------")
        println(f, "Variables:")

        for p in 1:n_products
            println(f, "$(xs[p]): $(value(x[p]))")
            solution_point[p] = value(x[p])
        end

        for l in 1:n_labor_types
            println(f, "$(ys[l]): $(value(y[l]))")
            solution_point[l+4] = value(y[l])
        end

        # println(f, "-----------------------------------------------")
        # println(f, "Shadow prices:")
        # for r in 1:n_resources
        #     println(f, "$r - $(resources[r]):  $(shadow_price(dmd[r]))")
        # end
    end

    variable_names = vcat(["x[$(string(p))]" for p in xs], ["y[$(string(l))]" for l in ys])
    constraint_names = ["mach1", "mach2", "skilled", "unskilled", "skilled_tot", "unskilled_tot"]

    obj_val = objective_value(model)

    return total_profit_f, constraints_f, variable_names, constraint_names, starting_point, solution_point, obj_val
end

#The datatype used in the tuples is used to determine if the variable is "real" or integer
variable_domains = [(-100.0,100.0), (-100.0,100.0), (-100.0,100.0), (-100.0,100.0), (0.0,1000.0), (0.0,1000.0)]

optimizers = [Clp.Optimizer, GLPK.Optimizer, Cbc.Optimizer]
opt_names = ["clp", "glpk", "cbc"]

problems = Dict()

for (opt, opt_name) in zip(optimizers, opt_names)
    objective_f, constraints_f, variable_names, constraint_names, starting_point, solution_point, obj_val = solve_model(opt, "prod_sched_$(opt_name).log")
    problems[opt_name] = (objective_f, constraints_f, variable_names, constraint_names, starting_point, solution_point, obj_val)
end

#All functions are equal, get first to initialize OptMap
objective_f, constraints_f, variable_names, constraint_names, _, _, obj_val = problems[opt_names[1]]
O = OptMap(variable_names, variable_domains, objective_f; constraint_names=constraint_names, constraints_f=constraints_f, verbose=true)
create_map!(O)

for opt_name in opt_names
    _, _, _, _, starting_point, solution_point, _ = problems[opt_name]

    path_to_solution = zeros(Float64, length(variable_domains), 2)
    path_to_solution[:,1] = starting_point
    path_to_solution[:,2] = solution_point

    draw_solution!(O, path_to_solution)
    save_all(O, "prod_sched_$(opt_name)", ".")
end

max_stars_repo_name: nmaterise/optmap

max_issues_repo_name: nmaterise/optmap

max_forks_repo_name: nmaterise/optmap


text:
struct OneBlockMatrix{T,Block<:AbstractMatrix{T}} <: AbstractMatrix{T}
    block::Block
    m::Int
    n::Int
    function OneBlockMatrix(block::Block, m, n) where {T,Block<:AbstractMatrix{T}}
        require_one_based_indexing(block)
        mb,nb = size(block)
        mb ≤ m && nb ≤ n ||
            throw(DimensionMismatch("Size of subblock ($mb,$nb) must be smaller than or equal to size of matrix ($m,$n)"))
        new{T,Block}(block, m, n)
    end
end

Base.size(A::OneBlockMatrix) = (A.m, A.n)
Base.size(A::OneBlockMatrix, i) = size(A)[i]

Base.similar(A::OneBlockMatrix, ::Type{T}, m::Integer, n::Integer) where T =
    OneBlockMatrix(similar(A.block, T), m, n)

Base.similar(A::OneBlockMatrix, ::Type{T}, (m,n)::Tuple{<:Integer,<:Tuple}) where T =
    similar(A, T, m, n)

Base.similar(A::OneBlockMatrix, ::Type{T}) where T =
    OneBlockMatrix(similar(A.block, T), A.m, A.n)

# Apparently ambiguous
Base.similar(A::OneBlockMatrix{T}, dims::Union{Integer, AbstractUnitRange}...) where T=
    similar(A, T, dims...)

function Base.getindex(A::OneBlockMatrix{T}, i::Integer, j::Integer) where T
    m,n = size(A.block)
    i ≤ m && j ≤ n && return A.block[i,j]
    zero(T)
end

function Base.getindex(A::OneBlockMatrix{T}, i::AbstractVector, j::AbstractVector) where T
    b = zeros(T, length(i), length(j))
    for (k,i) in enumerate(i)
        for (l,j) in enumerate(j)
            b[k,l] = A[i,j]
        end
    end
    b
end

function Base.replace_in_print_matrix(A::OneBlockMatrix,i::Integer,j::Integer,s::AbstractString)
    m,n = size(A.block)
    i ≤ m && j ≤ n && return s
    Base.replace_with_centered_mark(s)
end

for T in [Integer,Real]
    @eval begin
        LinearAlgebra.:(^)(A::OneBlockMatrix, p::$T) =
            OneBlockMatrix(A.block^p, A.m, A.n)
    end
end

Base.:(-)(A::OneBlockMatrix) =
    OneBlockMatrix(-A.block, A.m, A.n)

Base.:(*)(α::Number, A::OneBlockMatrix) =
    OneBlockMatrix(α*A.block, A.m, A.n)

Base.:(*)(A::OneBlockMatrix, α::Number) =
    OneBlockMatrix(α*A.block, A.m, A.n)

Base.:(/)(A::OneBlockMatrix, α::Number) = A * inv(α)
Base.:(//)(A::OneBlockMatrix, α::Number) = A * (1//α)

for f in [:exp, :sin, :cos, :tan, :sqrt]
    @eval begin
        function Base.$f(A::OneBlockMatrix)
            M,N = size(A)
            M == N || throw(DimensionMismatch("matrix is not square: dimensions are ($(M), $(N))"))

            b = A.block
            m,n = size(b)

            block = if m == n
                $f(b)
            else
                ee = eigen(A)
                Q = first(ee.vectors.blocks)
                Q*Diagonal($f.(ee.values[1:m]))*inv(Q)
            end

            fz = $f(zero(eltype(block)))
            D = Diagonal(fill(fz, M))

            BlockSparseDiagonal(D, [block-fz*I], [1])
        end
    end
end

# * Matrix arithmetic

for op in [:(+), :(-)]
    @eval begin
        function $op(A::OneBlockMatrix, B::OneBlockMatrix)
            size(A) == size(B) || throw(DimensionMismatch("Cannot add matrices of different sizes"))

            mA,nA = size(A.block)
            mB,nB = size(B.block)

            T = Base.promote_op($op, eltype(A), eltype(B))
            b = zeros(T, max(mA,mB), max(nA,nB))

            view(b, 1:mA, 1:nA) .+= A.block
            view(b, 1:mB, 1:nB) .+= $(op).(B.block)

            OneBlockMatrix(b, size(A)...)
        end

        function $op(I::UniformScaling, A::OneBlockMatrix)
            M,N = size(A)
            M == N || throw(DimensionMismatch("matrix is not square: dimensions are ($(M), $(N))"))

            T = Base.promote_op($op, eltype(A), eltype(I))

            D = Diagonal{T}(I, M)
            BlockSparseDiagonal(D, [$op(A.block)], [1])
        end

        function $op(A::OneBlockMatrix, I::UniformScaling)
            M,N = size(A)
            M == N || throw(DimensionMismatch("matrix is not square: dimensions are ($(M), $(N))"))

            T = Base.promote_op($op, eltype(A), eltype(I))

            D = $op(Diagonal{T}(I, M))
            BlockSparseDiagonal(D, [A.block], [1])
        end
    end
end

# * Matrix products

# ** Arbitrary matrices

function _mul!(y::AbstractArray{T,N}, A::OneBlockMatrix, x::AbstractArray{T,N},
               α::Number=true, β::Number=false) where {T,N}
    mY = size(y,1)
    nY = size(y,2)
    mA = size(A,1)
    nX = size(x,2)
    mY == mA && nY == nX ||
        throw(DimensionMismatch("Dimensions of y ($(mY),$(nY)) must agree with first dimension of A ($(mA)), and second dimension of x ($(nX)), respectively"))

    if iszero(β)
        y .= false
    elseif !isone(β)
        lmul!(β, y)
    end

    iszero(α) && return y

    m,n = size(A.block)
    mul!(view(y, 1:m, :), A.block, view(x, 1:n, :),
         α, true)

    y
end

LinearAlgebra.mul!(y::AbstractVector, A::OneBlockMatrix, x::AbstractVector,
                   α::Number=true, β::Number=false) = _mul!(y, A, x, α, β)

LinearAlgebra.mul!(y::AbstractMatrix, A::OneBlockMatrix, x::AbstractMatrix,
                   α::Number=true, β::Number=false) = _mul!(y, A, x, α, β)

# ** OneBlockMatrices

function LinearAlgebra.mul!(y::OneBlockMatrix, A::OneBlockMatrix, x::OneBlockMatrix,
                            α::Number=true, β::Number=false)
    mY,nY = size(y)
    mA = size(A,1)
    nX = size(x,2)
    mY == mA && nY == nX ||
        throw(DimensionMismatch("Dimensions of y ($(mY),$(nY)) must agree with first dimension of A ($(mA)), and second dimension of x ($(nX)), respectively"))

    mY,nY = size(y.block)
    mA,nA = size(A.block)
    mX,nX = size(x.block)

    mY == mA && nY == nX ||
        throw(DimensionMismatch("Dimensions of y subblock ($(mY),$(nY)) must agree with first dimension of A subblock ($(mA)), and second dimension of x subblock ($(nX)), respectively"))

    if iszero(β)
        y.block .= false
    elseif !isone(β)
        lmul!(β, y.block)
    end

    iszero(α) && return y

    interval = 1:min(nA,mX)
    mul!(y.block, view(A.block, :, interval), view(x.block, interval, :), α, true)

    y
end

function *(A::OneBlockMatrix, B::OneBlockMatrix)
    T = Base.promote_op(*, eltype(A), eltype(B))
    m = size(A.block,1)
    n = size(B.block,2)
    block = Matrix{T}(undef, m, n)
    C = OneBlockMatrix(block, size(A, 1), size(B, 2))
    mul!(C, A, B)
end

# ** Diagonal matrices

function mul!(y::OneBlockMatrix, A::OneBlockMatrix, D::Diagonal,
              α::Number=true, β::Number=false)
    mY,nY = size(y)
    mA = size(A,1)
    nD = size(D,2)
    mY == mA && nY == nD ||
        throw(DimensionMismatch("Dimensions of y ($(mY),$(nY)) must agree with first dimension of A ($(mA)), and second dimension of D ($(nD)), respectively"))

    m = size(A.block,1)
    mul!(y.block, A.block, Diagonal(view(D.diag, 1:m)), α, β)
    y
end

*(A::OneBlockMatrix, D::Diagonal) =
    mul!(similar(A), A, D)

function mul!(y::OneBlockMatrix, D::Diagonal, A::OneBlockMatrix,
              α::Number=true, β::Number=false)
    mY,nY = size(y)
    mD = size(D,1)
    nA = size(A,2)
    mY == mD && nY == nA ||
        throw(DimensionMismatch("Dimensions of y ($(mY),$(nY)) must agree with first dimension of D ($(mD)), and second dimension of A ($(nA)), respectively"))

    m = size(A.block,1)
    mul!(y.block, Diagonal(view(D.diag, 1:m)), A.block, α, β)
    y
end

*(D::Diagonal, A::OneBlockMatrix) =
    mul!(similar(A), D, A)

# ** Various banded matrices

const VariousBandedMatrices = Union{Tridiagonal,SymTridiagonal,Bidiagonal,BandedMatrix,SkewTridiagonal}

function _mul!(y::OneBlockMatrix, T::VariousBandedMatrices, A::OneBlockMatrix,
               α::Number=true, β::Number=false)
    mY,nY = size(y)
    mT = size(T,1)
    nA = size(A,2)
    mY == mT && nY == nA ||
        throw(DimensionMismatch("Dimensions of y ($(mY),$(nY)) must agree with first dimension of T ($(mT)), and second dimension of A ($(nA)), respectively"))

    b = y.block
    mY,nY = size(b)
    mA,nA = size(A.block)

    nrows = bandwidths(T)[1]
    mn = min(mA+nrows,size(T,1))
    mY == mn && nY == nA ||
        throw(DimensionMismatch("Dimensions of y subblock ($(mY),$(nY)) must agree with the dimensions of A subblock ($(mA)+nrows,$(nA))"))

    mul!(b, view(BandedMatrix(T), 1:mn, 1:mA), A.block)
    y
end

mul!(y::OneBlockMatrix, T::Union{Bidiagonal, SymTridiagonal, Tridiagonal}, A::OneBlockMatrix, args...) =
    _mul!(y, T, A, args...)

mul!(y::OneBlockMatrix, T::Union{BandedMatrix, SkewTridiagonal}, A::OneBlockMatrix, args...) =
    _mul!(y, T, A, args...)

function *(T::VariousBandedMatrices, A::OneBlockMatrix)
    m,n = size(A.block)
    U = Base.promote_op(*, eltype(T), eltype(A))
    M = size(T,1)
    nrows = bandwidths(T)[1]
    b = zeros(U, min(m+nrows,M), n)
    mul!(OneBlockMatrix(b, M, A.n), T, A)
end

function _mul!(y::OneBlockMatrix, A::OneBlockMatrix, T::VariousBandedMatrices,
               α::Number=true, β::Number=false)
    mY,nY = size(y)
    mA = size(A,1)
    nT = size(T,2)
    mY == mA && nY == nT ||
        throw(DimensionMismatch("Dimensions of y ($(mY),$(nY)) must agree with first dimension of A ($(mA)), and second dimension of T ($(nT)), respectively"))

    b = y.block
    mY,nY = size(b)
    mA,nA = size(A.block)

    ncols = bandwidths(T)[2]
    mn = min(nA+ncols,size(A,2))
    mY == mA && nY == mn ||
        throw(DimensionMismatch("Dimensions of y subblock ($(mY),$(nY)) must agree with the dimensions of A subblock ($(mA),$(nA)+1)"))

    mul!(b, A.block, view(BandedMatrix(T), 1:nA, 1:mn))
    y
end

mul!(y::OneBlockMatrix, A::OneBlockMatrix, T::Union{Bidiagonal, SymTridiagonal, Tridiagonal}, args...) =
    _mul!(y, A, T, args...)

mul!(y::OneBlockMatrix, A::OneBlockMatrix, T::Union{BandedMatrix, SkewTridiagonal}, args...) =
    _mul!(y, A, T, args...)

function *(A::OneBlockMatrix, T::VariousBandedMatrices)
    m,n = size(A.block)
    U = Base.promote_op(*, eltype(A), eltype(T))
    N = size(T,2)
    ncols = bandwidths(T)[2]
    b = zeros(U, m, min(n+ncols,N))
    mul!(OneBlockMatrix(b, A.m, N), A, T)
end

export OneBlockMatrix

# * Eigenfactorization

function LinearAlgebra.eigen(A::OneBlockMatrix{T}) where T
    M,N = size(A)
    M == N || throw(DimensionMismatch("matrix is not square: dimensions are ($(M), $(N))"))
    m,n = size(A.block)

    mn = max(m,n)
    Mmn = M-mn

    b = zeros(T, mn, mn)
    copyto!(view(b, 1:m, 1:n), A.block)

    ee = eigen(b)
    U = eltype(ee)
    λ = vcat(ee.values, zeros(U, Mmn))
    Φ = BlockSparseDiagonal(Diagonal(ones(U, M)), [ee.vectors - I], [1])

    Eigen(λ, Φ)
end


text:
using Test
#using MATLAB
#using SSST
using LinearAlgebra

# work around until added to General Registry
include("../src/SparseSymmetricTensors.jl")
using Main.SparseSymmetricTensors

function set_up(n,ord,nnz)
    x = rand(n)
    indices = rand(1:n,nnz,ord)
    for i = 1:nnz
        indices[i,:] = indices[i,sortperm(indices[i, :])]
    end
    vals = rand(nnz)
    x, indices, vals
end

# define common variables usable in all tests here

n = 5
valid_edges = [([1,2,3],1.0),([1,2,2],-1.0),([3,3,3],1)]
unordered_edge = [([2,3,1],1.0)]
negative_index = [([-1,2,3],1.0)]

# ------------      used in CoreOperations    ------------    #
 ord = 5
 tol = 1e-9
 nnz = 10
tempath = "tmp_tensor.ssten"

#include("SparseSymmetricTensors_Tests.jl")
#include("CoreOperations_Tests.jl")
include("Helpers_Tests.jl")
#include("HyperGraphAlgos_Tests.jl")
#include("NumericalRoutines_Tests.jl")