text:
# Much of this is event-handling to support interactivity

using Gtk.GConstants: GDK_KEY_Left, GDK_KEY_Right, GDK_KEY_Up, GDK_KEY_Down
using Gtk.GConstants.GdkEventMask: KEY_PRESS, SCROLL

abstract type CairoUnit <: Real end

Base.:+(x::U, y::U) where {U<:CairoUnit} = U(x.val + y.val)
Base.:-(x::U, y::U) where {U<:CairoUnit} = U(x.val - y.val)
Base.:<(x::U, y::U) where {U<:CairoUnit} = Bool(x.val < y.val)
Base.:>(x::U, y::U) where {U<:CairoUnit} = Bool(x.val > y.val)
Base.abs(x::U) where {U<:CairoUnit} = U(abs(x.val))
Base.min(x::U, y::U) where {U<:CairoUnit} = U(min(x.val, y.val))
Base.max(x::U, y::U) where {U<:CairoUnit} = U(max(x.val, y.val))
Base.isapprox(x::U, y::U; kwargs...) where U<:CairoUnit =
    isapprox(x.val, y.val; kwargs...)

Base.convert(::Type{T}, x::T) where {T<:CairoUnit} = x
(::Type{T})(x::CairoUnit) where T<:Real = T(x.val)

# Ambiguity resolution
Bool(x::CairoUnit) = Bool(x.val)

# The next three are for ambiguity resolution
Base.promote_rule(::Type{Bool}, ::Type{U}) where {U<:CairoUnit} = Float64
Base.promote_rule(::Type{BigFloat}, ::Type{U}) where {U<:CairoUnit} = BigFloat
Base.promote_rule(::Type{T}, ::Type{U}) where {T<:Irrational,U<:CairoUnit} = promote_type(T, Float64)
Base.promote_rule(::Type{T}, ::Type{U}) where {T<:Real,U<:CairoUnit} = promote_type(T, Float64)

"""
    DeviceUnit(x)

Represent a number `x` as having "device" units (aka, screen
pixels). See the Cairo documentation.
"""
struct DeviceUnit <: CairoUnit
    val::Float64
end

"""
    UserUnit(x)

Represent a number `x` as having "user" units, i.e., whatever units
have been established with calls that affect the transformation
matrix, e.g., [`Graphics.set_coordinates`](@ref) or
[`Cairo.set_matrix`](@ref).
"""
struct UserUnit <: CairoUnit
    val::Float64
end
UserUnit(x::UserUnit) = x

showtype(::Type{UserUnit}) = "UserUnit"
showtype(::Type{DeviceUnit}) = "DeviceUnit"

Base.show(io::IO, x::CairoUnit) = print(io, showtype(typeof(x)), '(', x.val, ')')

Base.promote_rule(::Type{U}, ::Type{D}) where {U<:UserUnit,D<:DeviceUnit} =
    error("UserUnit and DeviceUnit are incompatible, promotion not defined")

function convertunits(::Type{UserUnit}, c, x::DeviceUnit, y::DeviceUnit)
    xu, yu = device_to_user(getgc(c), x.val, y.val)
    UserUnit(xu), UserUnit(yu)
end
function convertunits(::Type{UserUnit}, c, x::UserUnit, y::UserUnit)
    x, y
end
function convertunits(::Type{DeviceUnit}, c, x::DeviceUnit, y::DeviceUnit)
    x, y
end
function convertunits(::Type{DeviceUnit}, c, x::UserUnit, y::UserUnit)
    xd, yd = user_to_device(getgc(c), x.val, y.val)
    DeviceUnit(xd), DeviceUnit(yd)
end

"""
    XY(x, y)

A type to hold `x` (horizontal), `y` (vertical) coordinates, where the
number increases to the right and downward. If used to encode mouse
pointer positions, the units of `x` and `y` are either
[`DeviceUnit`](@ref) or [`UserUnit`](@ref).
"""
struct XY{T}
    x::T
    y::T

    XY{T}(x, y) where {T} = new{T}(x, y)
    XY{U}(x::U, y::U) where {U<:CairoUnit} = new{U}(x, y)
    XY{U}(x::Real, y::Real) where {U<:CairoUnit} = new{U}(U(x), U(y))
end
XY(x::T, y::T) where {T} = XY{T}(x, y)
XY(x, y) = XY(promote(x, y)...)

function XY{U}(w::GtkCanvas, evt::Gtk.GdkEvent) where U<:CairoUnit
    XY{U}(convertunits(U, w, DeviceUnit(evt.x), DeviceUnit(evt.y))...)
end

function Base.show(io::IO, xy::XY{T}) where T<:CairoUnit
    print(io, "XY{$(showtype(T))}(", convert(Float64, xy.x), ", ", convert(Float64, xy.y), ')')
end
Base.show(io::IO, xy::XY) = print(io, "XY(", xy.x, ", ", xy.y, ')')

Base.convert(::Type{XY{T}}, xy::XY{T}) where {T} = xy
Base.convert(::Type{XY{T}}, xy::XY) where {T} = XY(T(xy.x), T(xy.y))

Base.:+(xy1::XY{T}, xy2::XY{T}) where {T} = XY{T}(xy1.x+xy2.x,xy1.y+xy2.y)
Base.:-(xy1::XY{T}, xy2::XY{T}) where {T} = XY{T}(xy1.x-xy2.x,xy1.y-xy2.y)

"""
    MouseButton(position, button, clicktype, modifiers)

A type to hold information about a mouse button event (e.g., a
click). `position` is the canvas position of the pointer (see
[`XY`](@ref)). `button` is an integer identifying the
button, where 1=left button, 2=middle button, 3=right
button. `clicktype` may be `BUTTON_PRESS` or
`DOUBLE_BUTTON_PRESS`. `modifiers` indicates whether any keys were
held down during the click; they may be any combination of `SHIFT`,
`CONTROL`, or `MOD1` stored as a bitfield (test with `btn.modifiers &
SHIFT`).

The fieldnames are the same as the argument names above.

    MouseButton{UserUnit}()
    MouseButton{DeviceUnit}()

Create a "dummy" MouseButton event. Often useful for the fallback to
Reactive's `filterwhen`.
"""
struct MouseButton{U<:CairoUnit}
    position::XY{U}
    button::UInt32
    clicktype::typeof(BUTTON_PRESS)
    modifiers::typeof(SHIFT)
    gtkevent
end
function MouseButton(pos::XY{U}, button::Integer, clicktype::Integer, modifiers::Integer, gtkevent=nothing) where U
    MouseButton{U}(pos, UInt32(button), oftype(BUTTON_PRESS, clicktype), oftype(SHIFT, modifiers), gtkevent)
end
function MouseButton{U}(w::GtkCanvas, evt::Gtk.GdkEvent) where U
    MouseButton{U}(XY{U}(w, evt), evt.button, evt.event_type, evt.state, evt)
end
function MouseButton{U}() where U
    MouseButton(XY(U(-1), U(-1)), 0, 0, 0, nothing)
end

"""
    MouseScroll(position, direction, modifiers)

A type to hold information about a mouse wheel scroll. `position` is the
canvas position of the pointer (see
[`XY`](@ref)). `direction` may be `UP`, `DOWN`, `LEFT`, or
`RIGHT`. `modifiers` indicates whether any keys were held down during
the click; they may be 0 (no modifiers) or any combination of `SHIFT`,
`CONTROL`, or `MOD1` stored as a bitfield.

    MouseScroll{UserUnit}()
    MouseScroll{DeviceUnit}()

Create a "dummy" MouseScroll event. Often useful for the fallback to
Reactive's `filterwhen`.
"""
struct MouseScroll{U<:CairoUnit}
    position::XY{U}
    direction::typeof(UP)
    modifiers::typeof(SHIFT)
end
function MouseScroll(pos::XY{U}, direction::Integer, modifiers::Integer) where U
    MouseScroll{U}(pos, oftype(UP, direction), oftype(SHIFT, modifiers))
end
function MouseScroll{U}(w::GtkCanvas, evt::Gtk.GdkEvent) where U
    MouseScroll{U}(XY{U}(w, evt), evt.direction, evt.state)
end
function MouseScroll{U}() where U
    MouseScroll(XY(U(-1), U(-1)), 0, 0)
end

# immutable KeyEvent
#     keyval
# end

"""
    MouseHandler{U<:CairoUnit}

A type with `Signal` fields for which you can `map` callback actions. The fields are:
  - `buttonpress` for clicks (of type [`MouseButton`](@ref));
  - `buttonrelease` for release events (of type [`MouseButton`](@ref));
  - `motion` for move and drag events (of type [`MouseButton`](@ref));
  - `scroll` for wheelmouse or track-pad actions (of type [`MouseScroll`](@ref));

`U` should be either [`DeviceUnit`](@ref) or [`UserUnit`](@ref) and
determines the coordinate system used for reporting mouse positions.
"""
struct MouseHandler{U<:CairoUnit}
    buttonpress::Signal{MouseButton{U}}
    buttonrelease::Signal{MouseButton{U}}
    motion::Signal{MouseButton{U}}
    scroll::Signal{MouseScroll{U}}
    ids::Vector{Culong}   # for disabling any of these callbacks
    widget::GtkCanvas

    function MouseHandler{U}(canvas::GtkCanvas) where U<:CairoUnit
        pos = XY(U(-1), U(-1))
        btn = MouseButton(pos, 0, BUTTON_PRESS, SHIFT)
        scroll = MouseScroll(pos, UP, SHIFT)
        ids = Vector{Culong}(undef, 0)
        handler = new{U}(Signal(btn), Signal(btn), Signal(btn), Signal(scroll), ids, canvas)
        # Create the callbacks
        push!(ids, Gtk.on_signal_button_press(mousedown_cb, canvas, false, handler))
        push!(ids, Gtk.on_signal_button_release(mouseup_cb, canvas, false, handler))
        push!(ids, Gtk.on_signal_motion(mousemove_cb, canvas, 0, 0, false, handler))
        push!(ids, Gtk.on_signal_scroll(mousescroll_cb, canvas, false, handler))
        handler
    end
end

"""
    GtkReactive.Canvas{U}(w=-1, h=-1, own=true)

Create a canvas for drawing and interaction. The relevant fields are:
  - `canvas`: the "raw" Gtk widget (from Gtk.jl)
  - `mouse`: the [`MouseHandler{U}`](@ref) for this canvas.

See also [`canvas`](@ref).
"""
struct Canvas{U}
    widget::GtkCanvas
    mouse::MouseHandler{U}
    preserved::Vector{Any}

    function Canvas{U}(w::Integer=-1, h::Integer=-1; own::Bool=true) where U
        gtkcanvas = GtkCanvas(w, h)
        # Delete the Gtk handlers
        for id in gtkcanvas.mouse.ids
            signal_handler_disconnect(gtkcanvas, id)
        end
        empty!(gtkcanvas.mouse.ids)
        # Initialize our own handlers
        mouse = MouseHandler{U}(gtkcanvas)
        set_gtk_property!(gtkcanvas, :is_focus, true)
        preserved = []
        canvas = new{U}(gtkcanvas, mouse, preserved)
        gc_preserve(gtkcanvas, canvas)
        canvas
    end
end

"""
    canvas(U=DeviceUnit, w=-1, h=-1) - c::GtkReactive.Canvas

Create a canvas for drawing and interaction. Optionally specify the
width `w` and height `h`. `U` refers to the units for the canvas (for
both drawing and reporting mouse pointer positions), see
[`DeviceUnit`](@ref) and [`UserUnit`](@ref). See also [`GtkReactive.Canvas`](@ref).
"""
canvas(::Type{U}=DeviceUnit, w::Integer=-1, h::Integer=-1) where {U<:CairoUnit} = Canvas{U}(w, h)
canvas(w::Integer, h::Integer) = canvas(DeviceUnit, w, h)

"""
    draw(f, c::GtkReactive.Canvas, signals...)

Supply a draw function `f` for `c`. This will be called whenever the
canvas is resized or whenever any of the input `signals` update. `f`
should be of the form `f(cnvs, sigs...)`, where the number of
arguments is equal to 1 + `length(signals)`.

`f` can be defined as a named function, an anonymous function, or
using `do`-block notation:

    using Graphics, Colors

    draw(c, imgsig, xsig, ysig) do cnvs, img, x, y
        copy!(cnvs, img)
        ctx = getgc(cnvs)
        set_source(ctx, colorant"red")
        set_line_width(ctx, 2)
        circle(ctx, x, y, 5)
        stroke(ctx)
    end

This would paint an image-Signal `imgsig` onto the canvas and then
draw a red circle centered on `xsig`, `ysig`.
"""
function Gtk.draw(drawfun::Function, c::Canvas, signals::Signal...)
    @guarded draw(c.widget) do widget
        # This used to have a `yield` in it to allow the Gtk event queue to run,
        # but that caused
        # https://github.com/JuliaGraphics/Gtk.jl/issues/368
        # and the bizarre failures in
        # https://github.com/JuliaImages/ImageView.jl/pull/153
        drawfun(widget, map(value, signals)...)
    end
    drawsig = map((values...)->draw(c.widget), signals...)
    push!(c.preserved, drawsig)
    drawsig
end

# Painting an image to a canvas
function Base.copy!(ctx::GraphicsContext, img::AbstractArray{C}) where C<:Union{Colorant,Number}
    save(ctx)
    reset_transform(ctx)
    Cairo.image(ctx, image_surface(img), 0, 0, Graphics.width(ctx), Graphics.height(ctx))
    restore(ctx)
end
Base.copy!(c::Union{GtkCanvas,Canvas}, img) = copy!(getgc(c), img)
function Base.fill!(c::Union{GtkCanvas,Canvas}, color::Colorant)
    ctx = getgc(c)
    w, h = Graphics.width(c), Graphics.height(c)
    rectangle(ctx, 0, 0, w, h)
    set_source(ctx, color)
    fill(ctx)
end

# TODO: Remove Matrix once the CairoImageSurface constructor has been
# generalized. Cairo.jl issue #252.
image_surface(img::Matrix{Gray24}) =
    Cairo.CairoImageSurface(Matrix(reinterpret(UInt32, img)), Cairo.FORMAT_RGB24)
image_surface(img::Matrix{RGB24})  =
    Cairo.CairoImageSurface(Matrix(reinterpret(UInt32, img)), Cairo.FORMAT_RGB24)
image_surface(img::Matrix{ARGB32}) =
    Cairo.CairoImageSurface(Matrix(reinterpret(UInt32, img)), Cairo.FORMAT_ARGB32)

image_surface(img::AbstractArray{T}) where {T<:Number} =
    image_surface(convert(Matrix{Gray24}, img))
image_surface(img::AbstractArray{T}) where {T<:ColorTypes.AbstractGray} =
    image_surface(convert(Matrix{Gray24}, img))
image_surface(img::AbstractArray{C}) where {C<:Color} =
    image_surface(convert(Matrix{RGB24}, img))
image_surface(img::AbstractArray{C}) where {C<:Colorant} =
    image_surface(convert(Matrix{ARGB32}, img))

# Coordiantes could be AbstractFloat without an implied step, so let's
# use intervals instead of ranges
struct ZoomRegion{T}
    fullview::XY{ClosedInterval{T}}
    currentview::XY{ClosedInterval{T}}
end

"""
    ZoomRegion(fullinds) -> zr
    ZoomRegion(fullinds, currentinds) -> zr
    ZoomRegion(img::AbstractMatrix) -> zr

Create a `ZoomRegion` object `zr` for selecting a rectangular
region-of-interest for zooming and panning. `fullinds` should be a
pair `(yrange, xrange)` of indices, an [`XY`](@ref) object, or pass a
matrix `img` from which the indices will be taken.

`zr.currentview` holds the currently-active region of
interest. `zr.fullview` stores the original `fullinds` from which `zr` was
constructed; these are used to reset to the original limits and to
confine `zr.currentview`.
"""
function ZoomRegion(inds::Tuple{AbstractUnitRange{I},AbstractUnitRange{I}}) where I<:Integer
    ci = ClosedInterval{RInt}.(inds)
    fullview = XY(ci[2], ci[1])
    ZoomRegion(fullview, fullview)
end
function ZoomRegion(fullinds::Tuple{AbstractUnitRange{I},AbstractUnitRange{I}},
                    curinds::Tuple{AbstractUnitRange{I},AbstractUnitRange{I}}) where I<:Integer
    fi = ClosedInterval{RInt}.(fullinds)
    ci = ClosedInterval{RInt}.(curinds)
    ZoomRegion(XY(fi[2], fi[1]), XY(ci[2], ci[1]))
end
ZoomRegion(img) = ZoomRegion(axes(img))
function ZoomRegion(fullview::XY, bb::BoundingBox)
    xview = oftype(fullview.x, bb.xmin..bb.xmax)
    yview = oftype(fullview.y, bb.ymin..bb.ymax)
    ZoomRegion(fullview, XY(xview, yview))
end

reset(zr::ZoomRegion) = ZoomRegion(zr.fullview, zr.fullview)

Base.axes(zr::ZoomRegion) = UnitRange.((zr.currentview.y, zr.currentview.x))

function interior(iv::ClosedInterval, limits::AbstractInterval)
    imin, imax = minimum(iv), maximum(iv)
    lmin, lmax = minimum(limits), maximum(limits)
    if imin < lmin
        imin = lmin
        imax = imin + IntervalSets.width(iv)
    elseif imax > lmax
        imax = lmax
        imin = imax - IntervalSets.width(iv)
    end
    oftype(limits, (imin..imax) ∩ limits)
end

function pan(iv::ClosedInterval, frac::Real, limits)
    s = frac*IntervalSets.width(iv)
    interior(minimum(iv)+s..maximum(iv)+s, limits)
end

"""
    pan_x(zr::ZoomRegion, frac) -> zr_new

Pan the x-axis by a fraction `frac` of the current x-view. `frac>0` means
that the coordinates shift right, which corresponds to a leftward
shift of objects.
"""
pan_x(zr::ZoomRegion, s) =
    ZoomRegion(zr.fullview, XY(pan(zr.currentview.x, s, zr.fullview.x), zr.currentview.y))

"""
    pan_y(zr::ZoomRegion, frac) -> zr_new

Pan the y-axis by a fraction `frac` of the current x-view. `frac>0` means
that the coordinates shift downward, which corresponds to an upward
shift of objects.
"""
pan_y(zr::ZoomRegion, s) =
    ZoomRegion(zr.fullview, XY(zr.currentview.x, pan(zr.currentview.y, s, zr.fullview.y)))

function zoom(iv::ClosedInterval, s::Real, limits)
    dw = 0.5*(s - 1)*IntervalSets.width(iv)
    interior(minimum(iv)-dw..maximum(iv)+dw, limits)
end

"""
    zoom(zr::ZoomRegion, scaleview, pos::XY) -> zr_new

Zooms in (`scaleview` < 1) or out (`scaleview` > 1) by a scaling
factor `scaleview`, in a manner centered on `pos`.
"""
function zoom(zr::ZoomRegion, s, pos::XY)
    xview, yview = zr.currentview.x, zr.currentview.y
    xviewlimits, yviewlimits = zr.fullview.x, zr.fullview.y
    centerx, centery = pos.x.val, pos.y.val
    w, h = IntervalSets.width(xview), IntervalSets.width(yview)
    fx, fy = (centerx-minimum(xview))/w, (centery-minimum(yview))/h
    wbb, hbb = s*w, s*h
    xview = interior(ClosedInterval(centerx-fx*wbb,centerx+(1-fx)*wbb), xviewlimits)
    yview = interior(ClosedInterval(centery-fy*hbb,centery+(1-fy)*hbb), yviewlimits)
    ZoomRegion(zr.fullview, XY(xview, yview))
end

"""
    zoom(zr::ZoomRegion, scaleview)

Zooms in (`scaleview` < 1) or out (`scaleview` > 1) by a scaling
factor `scaleview`, in a manner centered around the current view
region.
"""
function zoom(zr::ZoomRegion, s)
    xview, yview = zr.currentview.x, zr.currentview.y
    xviewlimits, yviewlimits = zr.fullview.x, zr.fullview.y
    xview = zoom(xview, s, xviewlimits)
    yview = zoom(yview, s, yviewlimits)
    ZoomRegion(zr.fullview, XY(xview, yview))
end

"""
    signals = init_pan_scroll(canvas::GtkReactive.Canvas,
                              zr::Signal{ZoomRegion},
                              filter_x::Function = evt->evt.modifiers == SHIFT || event.direction == LEFT || event.direction == RIGHT,
                              filter_y::Function = evt->evt.modifiers == 0 || event.direction == UP || event.direction == DOWN,
                              xpanflip = false,
                              ypanflip  = false)

Initialize panning-by-mouse-scroll for `canvas` and update
`zr`. `signals` is a dictionary holding the Reactive.jl signals needed
for scroll-panning; you can push `true/false` to `signals["enabled"]`
to turn scroll-panning on and off, respectively. Your application is
responsible for making sure that `signals` does not get
garbage-collected (which would turn off scroll-panning).

`filter_x` and `filter_y` are functions that return `true` when the
conditions for x- and y-scrolling are met; the argument is a
[`MouseScroll`](@ref) event. The defaults are that vertical scrolling
is triggered with an unmodified scroll, whereas horizontal scrolling
is triggered by scrolling while holding down the SHIFT key.

You can flip the direction of either pan operation with `xpanflip` and
`ypanflip`, respectively.
"""
function init_pan_scroll(canvas::Canvas{U},
                         zr::Signal{ZoomRegion{T}},
                         filter_x::Function = evt->(evt.modifiers & 0x0f) == SHIFT || evt.direction == LEFT || evt.direction == RIGHT,
                         filter_y::Function = evt->(evt.modifiers & 0x0f) == 0 && (evt.direction == UP || evt.direction == DOWN),
                         xpanflip = false,
                         ypanflip  = false) where {U,T}
    enabled = Signal(true)
    dummyscroll = MouseScroll{U}()
    pan = map(filterwhen(enabled, dummyscroll, canvas.mouse.scroll)) do event
        s = 0.1*scrollpm(event.direction)
        if filter_x(event)
            # println("pan_x: ", event)
            push!(zr, pan_x(value(zr), s))
        elseif filter_y(event)
            # println("pan_y: ", event)
            push!(zr, pan_y(value(zr), s))
        end
        nothing
    end
    Dict("enabled"=>enabled, "pan"=>pan)
end

"""
    signals = init_pan_drag(canvas::GtkReactive.Canvas,
                            zr::Signal{ZoomRegion},
                            initiate = btn->(btn.button == 1 && btn.clicktype == BUTTON_PRESS && btn.modifiers == 0))

Initialize click-drag panning that updates `zr`. `signals` is a
dictionary holding the Reactive.jl signals needed for pan-drag; you
can push `true/false` to `signals["enabled"]` to turn it on and off,
respectively. Your application is responsible for making sure that
`signals` does not get garbage-collected (which would turn off
pan-dragging).

`initiate(btn)` returns `true` when the condition for starting
click-drag panning has been met (by default, clicking mouse button
1). The argument `btn` is a [`MouseButton`](@ref) event.
"""
function init_pan_drag(canvas::Canvas{U},
                       zr::Signal{ZoomRegion{T}},
                       initiate::Function = pandrag_init_default) where {U,T}
    enabled = Signal(true)
    active = Signal(false)
    dummybtn = MouseButton{U}()
    local pos1, zr1, mtrx
    init = map(filterwhen(enabled, dummybtn, canvas.mouse.buttonpress)) do btn
        if initiate(btn)
            push!(active, true)
            # Because the user coordinates will change during panning,
            # convert to absolute position
            pos1 = XY(convertunits(DeviceUnit, canvas, btn.position.x, btn.position.y)...)
            zr1 = value(zr).currentview
            m = Cairo.get_matrix(getgc(canvas))
            mtrx = inv([m.xx m.xy 0; m.yx m.yy 0; m.x0 m.y0 1])
        end
        nothing
    end
    drag = map(filterwhen(active, dummybtn, canvas.mouse.motion)) do btn
        btn.button == 0 && return nothing
        xd, yd = convertunits(DeviceUnit, canvas, btn.position.x, btn.position.y)
        dx, dy, _ = mtrx*[xd-pos1.x, yd-pos1.y, 1]
        fv = value(zr).fullview
        cv = XY(interior(minimum(zr1.x)-dx..maximum(zr1.x)-dx, fv.x),
                interior(minimum(zr1.y)-dy..maximum(zr1.y)-dy, fv.y))
        if cv != value(zr).currentview
            push!(zr, ZoomRegion(fv, cv))
        end
    end
    finish = map(filterwhen(active, dummybtn, canvas.mouse.buttonrelease)) do btn
        btn.button == 0 && return nothing
        push!(active, false)
    end
    Dict("enabled"=>enabled, "active"=>active, "init"=>init, "drag"=>drag, "finish"=>finish)
end
pandrag_button(btn) = btn.button == 1 && (btn.modifiers & 0x0f) == 0
pandrag_init_default(btn) = btn.clicktype == BUTTON_PRESS && pandrag_button(btn)

"""
    signals = init_zoom_scroll(canvas::GtkReactive.Canvas,
                               zr::Signal{ZoomRegion},
                               filter::Function = evt->evt.modifiers == CONTROL,
                               focus::Symbol = :pointer,
                               factor = 2.0,
                               flip = false)

Initialize zooming-by-mouse-scroll for `canvas` and update
`zr`. `signals` is a dictionary holding the Reactive.jl signals needed
for scroll-zooming; you can push `true/false` to `signals["enabled"]`
to turn scroll-zooming on and off, respectively. Your application is
responsible for making sure that `signals` does not get
garbage-collected (which would turn off scroll-zooming).

`filter` is a function that returns `true` when the conditions for
scroll-zooming are met; the argument is a [`MouseScroll`](@ref)
event. The default is to hold down the CONTROL key while scrolling the
mouse.

The `focus` keyword controls how the zooming progresses as you scroll
the mouse wheel. `:pointer` means that whatever feature of the canvas
is under the pointer will stay there as you zoom in or out. The other
choice, `:center`, keeps the canvas centered on its current location.

You can change the amount of zooming via `factor` and the direction of
zooming with `flip`.
"""
function init_zoom_scroll(canvas::Canvas{U},
                          zr::Signal{ZoomRegion{T}},
                          filter::Function = evt->(evt.modifiers & 0x0f) == CONTROL,
                          focus::Symbol = :pointer,
                          factor = 2.0,
                          flip = false) where {U,T}
    focus == :pointer || focus == :center || error("focus must be :pointer or :center")
    enabled = Signal(true)
    dummyscroll = MouseScroll{U}()
    zm = map(filterwhen(enabled, dummyscroll, canvas.mouse.scroll)) do event
        if filter(event)
            # println("zoom scroll: ", event)
            s = factor
            if event.direction == UP
                s = 1/s
            end
            if flip
                s = 1/s
            end
            if focus == :pointer
                # println("zoom focus: ", event)
                push!(zr, zoom(value(zr), s, event.position))
            else
                # println("zoom center: ", event)
                push!(zr, zoom(value(zr), s))
            end
        end
    end
    Dict("enabled"=>enabled, "zoom"=>zm)
end

scrollpm(direction::Integer) =
    direction == UP ? -1 :
    direction == DOWN ? 1 :
    direction == RIGHT ? 1 :
    direction == LEFT ? -1 : error("Direction ", direction, " not recognized")

##### Callbacks #####
function mousedown_cb(ptr::Ptr, eventp::Ptr, handler::MouseHandler{U}) where U
    evt = unsafe_load(eventp)
    push!(handler.buttonpress, MouseButton{U}(handler.widget, evt))
    Int32(false)
end
function mouseup_cb(ptr::Ptr, eventp::Ptr, handler::MouseHandler{U}) where U
    evt = unsafe_load(eventp)
    push!(handler.buttonrelease, MouseButton{U}(handler.widget, evt))
    Int32(false)
end
function mousemove_cb(ptr::Ptr, eventp::Ptr, handler::MouseHandler{U}) where U
    evt = unsafe_load(eventp)
    pos = XY{U}(handler.widget, evt)
    # This doesn't support multi-button moves well, but those are rare in most GUIs and
    # users can examine `modifiers` directly.
    button = 0
    if evt.state & Gtk.GdkModifierType.BUTTON1 != 0
        button = 1
    elseif evt.state & Gtk.GdkModifierType.BUTTON2 != 0
        button = 2
    elseif evt.state & Gtk.GdkModifierType.BUTTON3 != 0
        button = 3
    end
    push!(handler.motion, MouseButton(pos, button, evt.event_type, evt.state))
    Int32(false)
end
function mousescroll_cb(ptr::Ptr, eventp::Ptr, handler::MouseHandler{U}) where U
    evt = unsafe_load(eventp)
    push!(handler.scroll, MouseScroll{U}(handler.widget, evt))
    Int32(false)
end


text:
export norm_relation

function _add_sunits_from_norm_relation!(c, UZK, N)
  cp = sort!(collect(Set(minimum(x) for x = c.FB.ideals)))
  K = N.K
  for i = 1:length(N)
    k, mk = subfield(N, i)
    @vprint :NormRelation 1 "Computing maximal order ..."
    zk = maximal_order(k)
    @vprint :NormRelation 1 "Computing lll basis ... "
    zk = lll(zk)
    @vprint :NormRelation 1 "Computing class group of $k... "
    class_group(zk)
    @vprint :NormRelation 1 "done"
    lpk = NfOrdIdl[ P[1] for p in cp for P = prime_decomposition(zk, p)]
    @vprint :NormRelation 1 "Now computing the S-unit group for lp of length $(length(lpk))"
    @assert length(lpk) > 0
    Szk, mS = Hecke.sunit_mod_units_group_fac_elem(lpk)

    D = Dict{nf_elem, nf_elem}()
    function N_mk(x, D, i)
      if haskey(D, x)
        return D[x]
      else
        y = N(x, i)
        D[x] = y
        return y
      end
    end

    cc = Vector{Tuple{Int, fmpz}}[]
    for j in 1:length(N.coefficients_gen[i])
      @vprint :NormRelation 1 "Inducing the action ... "
      z = induce_action(N, i, j, lpk, c.FB, cc)

      @vprint :NormRelation 1 "Feeding in the S-units of the small field ... "

      for l=1:ngens(Szk)
        u = mS(Szk[l])  #do compact rep here???
        valofnewelement = mul(mS.valuations[l], z)
        el_to_add = FacElem(Dict{nf_elem, AnticNumberField}((N(x, i, j), v) for (x,v) = u.fac))
        Hecke.class_group_add_relation(c, el_to_add, valofnewelement)#, always = false)
      end
    end

    # Skipping the units
    U, mU = unit_group_fac_elem(zk)
    for j=2:ngens(U) # I cannot add a torsion unit. He will hang forever.
      u = mU(U[j])
      Hecke.add_unit!(UZK, FacElem(Dict((N_mk(x, D, i), v) for (x,v) = u.fac)))
    end
    UZK.units = Hecke.reduce(UZK.units, UZK.tors_prec)
  end
end

function _compute_sunit_and_unit_group!(c, U, N, saturate = true)
  cp = sort!(collect(Set(minimum(x) for x = c.FB.ideals)))
  K = N.K
  skipped_units = FacElem{nf_elem, AnticNumberField}[]
  autos = automorphisms(field(N), copy = false)

  for i = 1:length(N)
    k, mk = subfield(N, i)
    zk = lll(maximal_order(k))
    print("Computing class group of $k... ")
    @time class_group(zk, redo = true, use_aut = true)
    println("done")
    lpk = NfOrdIdl[ P[1] for p in cp for P = prime_decomposition(zk, p)]
    println("Now computing the S-unit group for lp of length $(length(lpk))")
    @assert length(lpk) > 0
    @time Szk, mS = Hecke.sunit_mod_units_group_fac_elem(lpk)
    @show length(N.coefficients_gen[i])

    #D = Dict{nf_elem, nf_elem}()
    cc = Vector{Tuple{Int, fmpz}}[]
    induced = induce_action_from_subfield(N, i, lpk, c.FB, cc)
    for l in 1:ngens(Szk)
      u = mS(Szk[l])
      for j in 1:length(induced)
        aut = autos[j]
        valofnewelement = mul(mS.valuations[l], induced[j])
        el_to_add = FacElem(Dict{nf_elem, AnticNumberField}((aut(embedding(N, i)(x)), v) for (x,v) = u.fac))
        Hecke.class_group_add_relation(c, el_to_add, valofnewelement)#, always = false)
      end
    end
    println("done")

    u, mu = unit_group_fac_elem(zk)
    for n=2:ngens(u) # I cannot add a torsion unit. He will hang forever.
      uelem = mu(u[n])
      for j in 1:length(induced)
        aut = autos[j]
        lifted_unit = FacElem(Dict((aut(embedding(N, i)(x)), v) for (x,v) = uelem.fac))
        bb = Hecke.add_unit!(U, lifted_unit)
        @show bb
        if !bb
          push!(skipped_units, lifted_unit)
        end
      end
      U.units = Hecke.reduce(U.units, U.tors_prec)
    end
  end

  @show length(skipped_units)

  for lifted_unit in skipped_units
    @show Hecke.regulator(U.units, 64)
    Hecke._add_dependent_unit(U, lifted_unit)
  end

  U.tentative_regulator = Hecke.regulator(U.units, 64)

  println("Now saturating ... ")

  if saturate
    for (p, e) in factor(index(N))
      b = Hecke.saturate!(c, U, Int(p))
      while b
        b = Hecke.saturate!(c, U, Int(p))
      end
    end
  end
  return c, U
end

################################################################################
#
#  Redo the S-unit computation for Brauer norm relation
#
################################################################################

function _embed(N::NormRelation, i::Int, a::nf_elem)
  b = get(N.embed_cache_triv[i], a, nothing)
  if b === nothing
    _, mk = subfield(N, i)
    c = mk(a)
    N.embed_cache_triv[i][a] = c
    return c
  else
    return b
  end
end

function _get_sunits(N, i, lp)
  k = subfield(N, i)[1]
  if degree(k) > 10
    Gk, mGk = automorphism_group(k)
    if has_useful_brauer_relation(Gk)
      Szk, mS = _sunit_group_fac_elem_via_brauer(k, lp, true, index(N) == 1 ? 0 : index(N))
    else
      Szk, mS = Hecke.sunit_group_fac_elem(lp)
    end
  else
    Szk, mS = Hecke.sunit_group_fac_elem(lp)
  end
  return Szk, mS
end

# pure

function _add_sunits_from_brauer_relation!(c, UZK, N; invariant::Bool = false, compact::Int = 0, saturate_units::Bool = false)
  # I am assuming that c.FB.ideals is invariant under the action of the automorphism group used by N
  cp = sort!(collect(Set(minimum(x) for x = c.FB.ideals)))
  K = N.K
  add_unit_later = FacElem{nf_elem, AnticNumberField}[]
  @vprint :NormRelation 1 "Adding trivial relations\n"
  for I in c.FB.ideals
    a = I.gen_one
    Hecke.class_group_add_relation(c, K(a), fmpq(abs(a)^degree(K)), fmpz(1), orbit = false)
    b = I.gen_two.elem_in_nf
    bn = Hecke.norm_div(b, fmpz(1), 600)
    if nbits(numerator(bn)) < 550
      Hecke.class_group_add_relation(c, b, abs(bn), fmpz(1), orbit = false)
    end
  end
  for i = 1:length(N)
    if isdefined(N, :nonredundant) && !(i in N.nonredundant)
      continue
    end
    k, mk = subfield(N, i)
    @vprint :NormRelation 1 "Looking at the subfield $i / $(length(N)) with defining equation $(k.pol) \n"
    @vprint :NormRelation 1 "Computing lll basis of maximal order ...\n"
    zk = maximal_order(k)
    zk = lll(zk)
    lpk = NfOrdIdl[]
    for p in cp
      append!(lpk, prime_ideals_over(zk, p))
    end
    @assert length(lpk) > 0
    @vprint :NormRelation 1 "Computing sunit group for set of size $(length(lpk)) ... \n"
    Szk, mS = _get_sunits(N, i, lpk)
    @vprint :NormRelation 1 "Coercing the sunits into the big field ...\n"
    z = induce_action_just_from_subfield(N, i, lpk, c.FB, invariant)

    #found_torsion = false
    #CAREFUL! I AM ASSUMING THAT THE FIRST ELEMENT IS A TORSION UNIT!
    for l=2:ngens(Szk)
      @vprint :NormRelation 3 "Sunits $l / $(ngens(Szk))\n"
      u = mS(Szk[l])  #do compact rep here???
      if iszero(mS.valuations[l])
        if UZK.finished
          continue
        end
        if compact != 0
          @vprint :NormRelation 3 "  Compact presentation ...\n"
          @vtime :NormRelation 4 u = Hecke.compact_presentation(u, compact, decom = Dict{NfOrdIdl, Int}())
        elseif saturate_units
          @vprint :NormRelation 3 "  Compact presentation ...\n"
          @vtime :NormRelation 4 u = Hecke.compact_presentation(u, ispower(index(N))[2], decom = Dict{NfOrdIdl, Int}())
        end
        @vtime :NormRelation 4 img_u = FacElem(Dict{nf_elem, fmpz}((_embed(N, i, x), v) for (x,v) = u.fac))
        #=
        if !found_torsion
          fl = Hecke.istorsion_unit(img_u, false, 16)[1]
          if fl
            found_torsion = true
            continue
          end
        end
        =#
        has_full_rk = Hecke.has_full_rank(UZK)
        @vtime :NormRelation 4 ff = Hecke.add_unit!(UZK, img_u)
        if !has_full_rk && !ff
          push!(add_unit_later, img_u)
        end
      else
        valofnewelement = mul(mS.valuations[l], z)
        #=
        if isdefined(c.M, :basis)
          push!(deb_rr, (deepcopy(c.M.basis), deepcopy(valofnewelement)))
          rr = Hecke.reduce_right!(c.M.basis, deepcopy(valofnewelement))
          MM = matrix(c.M.basis)
          vv = zero_matrix(FlintZZ, 1, ncols(MM))
          for (jj, vvv) in valofnewelement
            vv[1, jj] = vvv
          end
          Hecke.reduce_mod_hnf_ur!(vv, MM)
          @assert iszero(rr) == iszero(vv)
          if iszero(rr)
            continue
          end
        end
        =#
        if compact != 0
          sup = Dict{NfOrdIdl, Int}((lpk[w], t) for (w, t) in mS.valuations[l])
          @vprint :NormRelation 3 "  Compact presentation ...\n"
          @vtime :NormRelation 4 u = Hecke.compact_presentation(u, compact, decom = sup)
        end
        @vtime :NormRelation 4 img_u = FacElem(Dict{nf_elem, fmpz}((_embed(N, i, x), v) for (x, v) = u.fac if !iszero(v)))
        @hassert :NormRelation 1 sparse_row(FlintZZ, [ (j, valuation(img_u, p)) for (j, p) in enumerate(c.FB.ideals) if valuation(img_u, p) != 0]) == valofnewelement
        @vtime :NormRelation 4 Hecke.class_group_add_relation(c, img_u, valofnewelement)
        #=
        if rank(c.M) == length(c.FB.ideals)
          h, piv = Hecke.class_group_get_pivot_info(c)
        end
        =#
      end
    end
    @vprint :NormRelation 4 "Reducing the units\n"
    @vtime :NormRelation 4 UZK.units = Hecke.reduce(UZK.units, UZK.tors_prec)
  end

  if length(add_unit_later) > 0
    for uu in add_unit_later
      Hecke.add_unit!(UZK, uu)
    end
    UZK.units = Hecke.reduce(UZK.units, UZK.tors_prec)
  end
  return nothing
end

# TODO: Encode that the i-th field is normal over Q
function induce_action_just_from_subfield(N::NormRelation, i, s, FB, invariant = false)
  S = FB.ideals
  ZK = order(S[1])

  z = zero_matrix(SMat, FlintZZ, 0, length(S))

  mk = embedding(N, i)
  zk = order(s[1])

  @assert mod(degree(ZK), degree(zk)) == 0
  reldeg = divexact(degree(ZK), degree(zk))

  for l in 1:length(s)
    v = Tuple{Int, fmpz}[]
    P = s[l]
    genofsl = elem_in_nf(_embed(N, i, P.gen_two.elem_in_nf))
    pmin = minimum(P, copy = false)
    # Compute the number of ideals
    # We use this to speed things up if L/K and L/Q are normal.
    # We are checking this below.
    local numb_ideal::Int
    if isnormal(number_field(ZK))
      rele = divexact(ramification_index((FB.fb[pmin].lp)[1][2]), ramification_index(P))
      relf = divexact(degree((FB.fb[pmin].lp)[1][2]), degree(P))
      @assert mod(reldeg, rele * relf) == 0
      numb_ideal = divexact(reldeg, rele * relf)
    else
      numb_ideal = -1
    end
    found = 0
    for k in 1:length(S)
      Q = S[k]
      if minimum(Q, copy = false) == pmin
        if genofsl in Q
          found += 1
          @assert mod(ramification_index(Q), ramification_index(s[l])) == 0
          push!(v, (k, divexact(ramification_index(Q), ramification_index(s[l]))))
        end
      end
      if invariant && N.isnormal[i] && isnormal(number_field(ZK))
				if found == numb_ideal
					break
				end
			end
    end
    sort!(v, by = x -> x[1])
    push!(z, sparse_row(FlintZZ, v))
  end
  return z
end

function norm_relation(K::AnticNumberField, coprime::Int = 0; small_degree = true, cached = true)
  local N
  if cached
    if has_attribute(K, :norm_relation)
      N = get_attribute(K, :norm_relation)::Vector{NormRelation{Int}}
      if coprime == 0
        return true, N[1]::NormRelation{Int}
      else
        for i in 1:length(N)
          if coprime == N[i].isoptimal
            return true, N[i]::NormRelation{Int}
          end
        end
        fl, M = has_coprime_norm_relation(K, fmpz(coprime))
        if fl
          push!(N, M)
          return true, M
        end
        return false, NormRelation{Int}()
      end
    end
  end
  if coprime == 0
    if !has_useful_brauer_relation(automorphism_group(K)[1])
      return false, NormRelation{Int}()
    end
    M = _norm_relation_setup_generic(K, pure = true, small_degree = small_degree)
    set_attribute!(K, :norm_relation, NormRelation{Int}[M])
    return true, M::NormRelation{Int}
  else
    fl, M = has_coprime_norm_relation(K, fmpz(coprime))
    if fl
      set_attribute!(K, :norm_relation, NormRelation{Int}[M])
      return true, M::NormRelation{Int}
    end
    return false, NormRelation{Int}()
  end
end

function _sunit_group_fac_elem_quo_via_brauer(K::AnticNumberField, S, n::Int, invariant::Bool = false; saturate_units::Bool = false)
  @vprint :NormRelation 1 "Setting up the norm relation context ... \n"
  fl, N = norm_relation(K, n)
  if !fl
    fl, N = norm_relation(K, 0, small_degree = false)
  end
  @assert fl
  @vprint :NormRelation 1 "Using norm relation $N\n"
  compact = 0
  g = gcd(index(N), n)
  compact = 0
  if !isone(g)
    compact = ispower(g)[2]
  end
  return __sunit_group_fac_elem_quo_via_brauer(N, S, n, invariant, compact, saturate_units = saturate_units)::Tuple{GrpAbFinGen, Hecke.MapSUnitGrpFacElem}
end

function _sunit_group_fac_elem_via_brauer(K::AnticNumberField, S::Vector{NfOrdIdl}, invariant::Bool = false, compact::Int = 0)
  @vprint :NormRelation 1 "Setting up the norm relation context ... \n"
  fl, N = norm_relation(K, 0, small_degree = false)
  @assert fl
  @vprint :NormRelation 1 "Using norm relation $N\n"
  return __sunit_group_fac_elem_quo_via_brauer(N, S, 0, invariant, compact)::Tuple{GrpAbFinGen, Hecke.MapSUnitGrpFacElem}
end

function sunit_group_fac_elem_via_brauer(N::NormRelation, S, invariant::Bool = false)
  return _sunit_group_fac_elem_quo_via_brauer(N, S, 0, invariant)::Tuple{GrpAbFinGen, Hecke.MapSUnitGrpFacElem}
end

function _setup_for_norm_relation_fun(K, S = prime_ideals_up_to(maximal_order(K), Hecke.factor_base_bound_grh(maximal_order(K))))
  ZK = order(S[1])
  FB = NfFactorBase(ZK, S)
  c = Hecke.class_group_init(FB)
  U1 = Hecke._get_UnitGrpCtx_of_order(ZK, false)
  if U1 == nothing
    UZK = Hecke.UnitGrpCtx{FacElem{nf_elem, AnticNumberField}}(ZK)
    Hecke._set_UnitGrpCtx_of_order(ZK, UZK)
  else
    UZK = U1
  end
  return c, UZK
end

function _find_perm(v::Vector{NfOrdIdl}, w::Vector{NfOrdIdl})
  p = Dict{Int, Int}()
  for i = 1:length(v)
    pi = findfirst(isequal(v[i]), w)
    @assert pi != nothing
    p[pi] = i
  end
  return p
end

function __sunit_group_fac_elem_quo_via_brauer(N::NormRelation, S::Vector{NfOrdIdl}, n::Int, invariant::Bool = false, compact::Int = 0; saturate_units::Bool = false)
  O = order(S[1])
  K = N.K
  if invariant
    c, UZK = _setup_for_norm_relation_fun(K, S)
    _add_sunits_from_brauer_relation!(c, UZK, N, invariant = true, compact = compact, saturate_units = saturate_units)
  else
    cp = sort!(collect(Set(minimum(x) for x = S)))
    Sclosed = NfOrdIdl[]
    # If K is non-normal and N does not use the full automorphism group, we
    # might be enlarging S too much.
    for p in cp
      append!(Sclosed, prime_ideals_over(O, p))
    end
    if length(Sclosed) == length(S)
      invariant = true
      Sclosed = S
    end
    if !invariant
      @vprint :NormRelation 1 "I am not Galois invariant. Working with S of size $(length(Sclosed))\n"
    end
    c, UZK = _setup_for_norm_relation_fun(K, Sclosed)
    _add_sunits_from_brauer_relation!(c, UZK, N, invariant = true, compact = compact, saturate_units = saturate_units)
  end

  if gcd(index(N), n) != 1
    # I need to saturate
    for (p, e) in factor(index(N))
      @vprint :NormRelation 1 "Saturating at $p \n"
      b = Hecke.saturate!(c, UZK, Int(p), 3.5, easy_root = true, use_LLL = true)
      while b
        b = Hecke.saturate!(c, UZK, Int(p), 3.5, easy_root = true, use_LLL = true)
      end
    end
    UZK.finished = true
  end

  if index(N) == 1
    UZK.finished = true
  end

  if saturate_units && !UZK.finished
    for (p, e) in factor(index(N))
      @vprint :NormRelation 1 "Saturating at $p \n"
      b = Hecke.saturate!(UZK, Int(p), 3.5, easy_root = true, use_LLL = true)
      while b
        b = Hecke.saturate!(UZK, Int(p), 3.5, easy_root = true, use_LLL = true)
      end
    end
    UZK.finished = true
  end

  # This makes c.R.gen a basis of the S-units (modulo torsion)
  c = Hecke.RelSaturate.simplify(c, UZK, use_LLL = true)
  perm_ideals = _find_perm(S, c.FB.ideals)
  if invariant
    sunitsmodunits = FacElem{nf_elem, AnticNumberField}[x for x in c.R_gen] # These are generators for the S-units (mod units, mod n)
    valuations_sunitsmodunits = Vector{SRow{fmpz}}(undef, length(S))
    for i = 1:length(sunitsmodunits)
      r = Tuple{Int, fmpz}[(perm_ideals[j], v) for (j, v) in c.M.bas_gens[i]]
      sort!(r, lt = (a,b) -> a[1] < b[1])
      valuations_sunitsmodunits[i] = sparse_row(FlintZZ, r)
    end
  else
    # I need to extract the S-units from the Sclosed-units
    # Now I need to find the correct indices in the c.FB.ideals
    sunitsmodunits = FacElem{nf_elem, AnticNumberField}[]
    valuations_sunitsmodunits = SRow{fmpz}[]
    ind = Int[]
    for P in S
      for i in 1:length(c.FB.ideals)
        if P == c.FB.ideals[i]
          push!(ind, i)
          break
        end
      end
    end
    sort!(ind)
    # ind = indices of S inside c.FB.ideals
    @assert length(Sclosed) == length(c.FB.ideals)
    @assert length(ind) == length(S)
    z = zero_matrix(FlintZZ, length(c.R_gen), length(Sclosed) - length(S))
    for i in 1:length(c.R_gen)
      k = 1
      for j in 1:length(Sclosed)
        if !(j in ind)
          z[i, k] = c.M.bas_gens[i, j]
          if k == ncols(z)
            break
          end
          k = k + 1
        end
      end
    end
    k, K = left_kernel(z)
    for i in 1:nrows(K)
      if iszero_row(K, i)
        continue
      end
      push!(sunitsmodunits, FacElem(c.R_gen, fmpz[K[i, j] for j in 1:ncols(K)]))
      v_c = sum(SRow{fmpz}[K[i, j]*c.M.bas_gens[j] for j = 1:ncols(K)])
      r = Tuple{Int, fmpz}[(perm_ideals[j], v) for (j, v) in v_c]
      sort!(r, lt = (a,b) -> a[1] < b[1])
      push!(valuations_sunitsmodunits, sparse_row(FlintZZ, r))
    end
  end

  unitsmodtorsion = UZK.units # These are generators for the units (mod n)
  T, mT = torsion_unit_group(O)
  Q, mQ = quo(T, n, false)
  # Q need not be in SNF, but we need that it is generated by one element.
  @assert ngens(Q) == 1
  m = order(Q)
  if !isone(m)
    units = Vector{FacElem{nf_elem, AnticNumberField}}(undef, length(unitsmodtorsion)+1)
    units[1] = FacElem(elem_in_nf(mT(mQ\Q[1])))
    for i = 2:length(units)
      units[i] = unitsmodtorsion[i-1]
    end
    res_group = abelian_group(append!(fmpz[m], [fmpz(n) for i in 1:(length(sunitsmodunits) + length(unitsmodtorsion))]))
  else
    units = unitsmodtorsion
    res_group = abelian_group([fmpz(n) for i in 1:(length(sunitsmodunits) + length(unitsmodtorsion))])
  end
  local exp
  let res_group = res_group, units = units
    function exp(a::GrpAbFinGenElem)
      @assert parent(a) == res_group
      z = prod(units[i]^a[i] for i = 1:length(units))
      if !isempty(sunitsmodunits)
        zz = prod(sunitsmodunits[i]^a[length(units) + i] for i in 1:length(sunitsmodunits))
        mul!(z, z, zz)
      end

      for (k, v) in z.fac
        if iszero(v)
          delete!(z.fac, k)
        end
      end
      return z
    end
  end

  disclog = function(a)
    throw(NotImplemented())
  end

  r = Hecke.MapSUnitGrpFacElem()
  r.valuations = Vector{SRow{fmpz}}(undef, ngens(res_group))
  for i = 1:length(units)
    r.valuations[i] = sparse_row(FlintZZ)
  end
  for i = 1:length(sunitsmodunits)
    r.valuations[i+length(units)] = valuations_sunitsmodunits[i]
  end
  r.idl = S
  r.isquotientmap = n

  r.header = Hecke.MapHeader(res_group, FacElemMon(nf(O)), exp, disclog)
  @hassert :NormRelation 9001 begin
    _S, _mS = Hecke.sunit_group_fac_elem(S)
    @show _S
    @show res_group
    _Q, _mQ = quo(_S, n)
    V = quo(_Q, [_mQ(_mS\(r(res_group[i]))) for i in 1:ngens(res_group)])
    order(_Q) == order(res_group) && order(V[1]) == 1
  end
  @hassert :NormRelation 9000 begin
    fl = true
    for i = 1:ngens(res_group)
      el = r(res_group[i])
      if sparse_row(FlintZZ, [ (j, valuation(el, S[j])) for j = 1:length(S) if valuation(el, S[j]) != 0]) != r.valuations[i]
        fl = false
        break
      end
    end
    fl
  end
  return (res_group, r)::Tuple{GrpAbFinGen, Hecke.MapSUnitGrpFacElem}
end

ext: jl


text:
# define exceptions for stopping conditions in algorithm
type OptimalPointFound <: Exception end
type UnboundedFunctionException <: Exception end

function get_min(x, d, assts)
    n = length(assts)
    e = zeros(n)
    e[assts] = d
    ratios = []
    for i in 1:n
        if (assts[i]) & (e[i] > 0)
            push!(ratios, (x[i]/e[i], i))
        end
    end

    try
        return minimum(ratios)
    catch
        println(ratios)
    end
end

function get_q(s_N, assts)
    n = length(assts)
    e = zeros(n)
    e[!assts] = s_N

    try
        return indmax(e.<0)
    catch
            print e
    end
end

# a single iteration of the simplex algorithm
function simplex_iteration(assts, x, A, b, c)
    """
    We have a linear program of the form: min c^T x s.t. Ax = b

    arguments:
        assts: true indicates the element belongs to the basic set and false to the nonbasic set
        x: the vector x as given above in the LP
        A: constraint matrix, of the form Ax = b
        b: solutions of constraints
        c: coefficients of the LP function
    """

    B = A[:,assts]
    N = A[:,!assts]

    x_B = x[assts]
    x_N = x[!assts]

    # solve B^T lambda = c_B for lambda
    lambda = B' \ c[assts]

    # compute s_N = c_N - N^\top \lambda
    s_N = c[!assts] - N' * lambda

    # check if all elements of s_N are >= 0
    if all(s_N .>= 0)
        throw(OptimalPointFound())
    end

    # select q \in N with s_q < 0 as entering index
    q = get_q(s_N, assts) # here we always choose the first negative element

    # solve Bd = A_q for d
    d = B \ A[:,q]

    # check d
    if all(s_N .>= 0)
        throw(UnboundedFunctionException())
    end

    # calculate x_q+
    x_q, p = get_min(x, d, assts)

    # update x
    x[assts] = x[assts] - d * x_q
    x[q] = x_q

    # update the assts
    assts[p] = !assts[p]
    assts[q] = !assts[q]

    # returns new assts, x
    return assts, x

end

function run_simplex(assts, x, A, b, c, num_iters, verbose=false)

    for i in 1:num_iters
        #println("Starting iteration: ", i)

        try
            assts, x = simplex_iteration(assts, x, A, b, c)

            if verbose
                println("Resulting x, assts: \n\t", (x,assts))
            end

        catch err
            if isa(err, UnboundedFunctionException)

                println("Unbounded function exception!")
                return assts, x

            elseif isa(err, OptimalPointFound)

                println("Optimal point found!")

                if verbose
                    println("Optimal x, assts: \n\t", (x,assts))
                end

                return assts, x
            else
                println(err)

                throw(err)
            end
        end

    end
    return assts, x
end

function get_starting_values(A, b)
    """
    given a LP min c'x s.t. Ax = b, we want to create the new LP e'z s.t. Ax + Ez = b, where x=0, e=ones(.)

    A: the original constraint matrix
    b: the original b

    returns new A, x, c, assts to pass into simplex algorithm
    """
    nrows, ncols = size(A)
    # init x, c
    x = zeros(ncols + nrows)
    c = zeros(ncols + nrows)

    # add z to x
    x[ncols+1:ncols+nrows,:] = abs(b')
    c[ncols+1:ncols+nrows,:] = 1.

    # assignments are non-zero elements of x
    assts = (x .!= 0)

    E = zeros(nrows, nrows)
    vals = b .> 0

    for i in 1:nrows
        if vals[i]
            E[i,i] = 1
        else
            E[i,i] = -1
        end
    end

    return [A E], x, c, assts

end

function run_two_stage_simplex(A, b, c, num_iters)
    time1 = time()
    A0, x0, c0, assts0 = get_starting_values(A, b)

    nrows, ncols = size(A)
    l = length(x0)

    println("\nRunning Phase I of simplex algorithm\n")
    # run phase I
    assts_start, x_start = run_simplex(assts0,x0,A0,b,c0,num_iters)

    # run phase II
    try
        @assert x_start[l] == 0
    catch err
        println("No feasible solution for m=", size(A)[2])
    end

    # discard z
    xnew = x_start[1:l-nrows]
    assts_new = assts_start[1:l-nrows]

    println("\nRunning Phase II of simplex algorithm \n")

    x, assts = run_simplex(assts_new,xnew,A,b,c,num_iters)

    time2 = time()

    @printf("\nTime: %.5f (s)", time2-time1)

    return x, assts, (time2-time1)

end

#run_two_stage_simplex(A, b, c, 10)

function solve_subproblem(assts, x_k, G, c, A, b)
    g_k = G * x_k + c

    ## if all assignments are false
    if all(assts .== false)
        return G \ -g_k

    else
        dimx = length(x_k)
        A_w = A[assts,:]
        k,k = size(G)
        m,n = size(A_w)

        # construct matrix
        mat = [G A_w'; A_w zeros(k+m - n,k+m - n)]

        # pad with 0's
        d, = size(g_k)
        answ = [-g_k; zeros(k+m - d)]

        # compute solution and get x
        sol = mat \ answ
        x = sol[1:dimx,:]

        return x
    end
end

function compute_multipliers(assts, x_k, G, c, A)
    g = G * x_k + c
    lambda = A[assts,:]' \ g

    n = length(assts)
    return lambda, (1:n)[assts]
end

function compute_step_length(A, b, x_k, p_k, assts)
    """
    Returns step length alpha_k and blocking constraint of choice -- if the latter is 0, then alpha_k = 1
    """
    n = length(b)

    # by default, return 1 with no blocking constraints
    min = 1
    blocking = 0

    for i in 1:n
        # if i is in the working set, continue
        if assts[i]
            continue
        end

        a_i = A[i,:] ## this is a row vector so no need to transpose
        prod = a_i * p_k
        @assert length(prod) == 1

        # if a_i^ p_k >= 0, continue
        if (prod)[1] >= 0
            continue
        end

        val = (b[i] - a_i * x_k) / (prod)
        @assert length(val) == 1
        if val[1] < min
            min = val[1]
            blocking = i
        end
    end

    return min, blocking
end

function active_set_iteration(assts, x_k, G, c, A, b)
    """
    assts: assignments of constraint indices to the working set -- 1's are in the working set, 0's are not
    x_k: the current value of x_k obtained from the previous iteration
    G: the matrix in the quadratic function, i.e., x^T G x + x^T c
    c: the c from the quadratic function above
    A: matrix, where each row is the a_i from the constraint
    b: vector where each index is the b_i corresponding to constraint i
    """
    # solve subproblem to find p_k
    p_k = solve_subproblem(assts, x_k, G, c, A, b)

    #if all(p_k == 0) sufficient close to 0
    if all(abs(p_k-0) .<= 1e-6)

        if all(assts .== false)
            throw(OptimalPointFound())
        else

            # compute lagrange multipliers lambda_i that satisfy
            lambda, inds = compute_multipliers(assts, x_k, G, c, A)
        end

        # if lambda_i >= 0 for all i in W_k \cap I
        if all(lambda .>= 0)
            # stop with solution x* = x_k
            throw(OptimalPointFound())

        # (not all lambda_i >= 0)
        else
            # index of most negative multiplier lambda_j
            j = inds[indmin(lambda)]

            # remove j from working set
            assts[j] = false
            x_new = x_k
        end

    # (p_k != 0)
    else
        # compute step length alpha_k
        alpha_k, blocking = compute_step_length(A, b, x_k, p_k, assts)

        # update x_{k} += alpha_k * p_k
        x_new = x_k + alpha_k * p_k

        # if there are blocking constraints, update W_k (assts) by adding one of them to W_k
        if blocking != 0
            assts[blocking] = true
        end
    end

    return assts, x_new

end

function run_active_set(num_iters, assts, x_0, G, c, A, b, verbose=true)
    x_k = x_0
    for i in 1:num_iters
        if verbose
            println("Iteration:\t", i)
        end
        try
            assts, x_k = active_set_iteration(assts, x_k, G, c, A, b)
            if verbose
                println("Resulting x, assts: \n\t", (x_k,assts))
            end
        catch err
            if isa(err, OptimalPointFound)

                println("Optimal point found!")

                if verbose
                    println("Optimal x, assts: \n\t", (x_k,assts))
                end

                return assts, x_k
            else
                println(err)
                throw(err)
            end
        end

    end
    return assts, x_k
end

function run_active_set_full(num_iters, assts, x_0, G, c, A, b, verbose=true)
    x_k = x_0
    vals = []
    for i in 1:num_iters
        if verbose
            println("Iteration:\t", i)
        end
        try
            assts, x_k = active_set_iteration(assts, x_k, G, c, A, b)

            val = 1/2 * x_k' * G * x_k + c' * x_k
            push!(vals, val)

            if verbose
                println("Resulting x, assts: \n\t", (x_k,assts))
            end
        catch err
            if isa(err, OptimalPointFound)

                println("Optimal point found!")

                if verbose
                    println("Optimal x, assts: \n\t", (x_k,assts))
                end

                return assts, x_k, vals
            else
                println(err)
                throw(err)
            end
        end

    end
    return assts, x_k, vals
end


text:
using JuLIP
using Pkg
# @show pathof(JuLIP)

@info("nthreads = $(Threads.nthreads())")

function runtests(at, usethreads)
   JuLIP.usethreads!(usethreads)
   energy(at); forces(at); virial(at)
   tE = minimum( (@elapsed energy(at)) for _=1:5 )
   tF = minimum( (@elapsed forces(at)) for _=1:5 )
   tV = minimum( (@elapsed virial(at)) for _=1:5 )
   tf = usethreads ? "t" : "f"
   println(tf * " -> E: $tE; F: $tF; V: $tV")
   nothing
end

for mul in [5, 10, 20]
   at = bulk(:Cu, cubic=true) * mul
   set_calculator!(at, JuLIP.Potentials.EMT())
   @info("Run timings with $(length(at)) atoms")
   runtests(at, true)
   runtests(at, false)
end


text:
mutable struct IntrusiveListNode{T}
    value::T
    next::IntrusiveListNode{T}
    prev::IntrusiveListNode{T}

    function IntrusiveListNode{T}(value::T) where {T}
        x = new{T}(value)
        x.next = x
        x.prev = x
        x
    end
end

IntrusiveListNode(value::T) where {T} = IntrusiveListNode{T}(value)

Base.show(io::IO, x::IntrusiveListNode) = show(io, x.value)


text:
using Pkg

Pkg.add(
    url="https://github.com/terasakisatoshi/HelloCExecutable_jll.jl.git",
    rev="HelloCExecutable-v0.1.0+0",
)


text:
using ISA #References: Sandoval, Steven, and Phillip L. De Leon. "The Instantaneous Spectrum: A General Framework for Time-Frequency Analysis." IEEE Transactions on Signal Processing 66.21 (2018): 5679-5693.
using Plots

fs=500
t = collect(0.0:1/fs:2.0)

a₀(t) = exp(-t^2)
ω₀(t) = 50
φ₀ = 0
𝐶₀ = Tuple([a₀,ω₀,φ₀])

a₁(t) = exp(-(t-1.0)^2)
ω₁(t) = 2*π*100 + 150*sin(2t)
φ₁ = π
𝐶₁ = Tuple([a₁,ω₁,φ₁])

a₂(t) = 0.8*cos(2t)
ω₂(t) = 2*π*(-100 + 50t)
φ₂ = π
𝐶₂ = Tuple([a₂,ω₂,φ₂])

𝑆 = [𝐶₀,𝐶₁,𝐶₂]
z = AMFMmodel(𝑆)

p1 = isaPlot3d(z, t, FreqUnits="Hz")

windowFunction = ones(Float64,25)
frameAdvance = 1
zₖ = frameSignal(z(t), windowFunction, frameAdvance)
Z,freqs = STFT(z(t), windowFunction, frameAdvance,fs=fs)

p2 = contour(t[1:frameAdvance:end], freqs, abs.(reduce(hcat,Z)), fill=:true, seriescolor=cgrad(ISA.cmap),levels=5 )

plot(p1,p2,layout=(2,1))

max_stars_repo_name: kricheso/ISA

max_issues_repo_name: kricheso/ISA

max_forks_repo_name: kricheso/ISA


text:
using VCFTools

truth  = "target.full.vcf.gz"
mendel = "mendel.imputed.vcf.gz"
beagle = "beagle.imputed.maf0.05.vcf.gz"
#mmac4  = "minimac4.result.dose.vcf.gz"
imp5   = "impute5.result.vcf.gz"

Xtrue = convert_gt(UInt8, truth, msg="importing...")
Xmendel = convert_gt(UInt8, mendel, msg="importing...")
Xbeagle = convert_gt(UInt8, beagle, msg="importing...")
#Xmmac4 = convert_gt(UInt8, mmac4, msg="importing...")
Ximp5  = convert_gt(UInt8, imp5, msg="importing...")
n, p = size(Xtrue)

println("error mendel = $(sum(Xmendel .!= Xtrue) / n / p)")
println("error beagle = $(sum(Xbeagle .!= Xtrue) / n / p)")
#println("error minimac4 = $(sum(Xmmac4 .!= Xtrue) / n / p)")
println("error impute5 = $(sum(Ximp5 .!= Xtrue) / n / p)")

max_stars_repo_name: Hua-Zhou/MendelImpute


text:
module finiteMPS

# import
using TensorOperations, KrylovKit, LinearAlgebra
using MPOmodule

# export
export randMPS
export rightCanMPS, leftCanMPS, leftCanSite, rightCanSite, mixedCanMPS
export applyTM_OP, applyTM_MPO
export evo_sweep, evo_sweep_2Site, vmps_sweep, iter_applyMPO
export getOverlap, getExpValue, singleSiteExpValue
export applyHeff, applyHCeff

#include("sturctMPS.jl")
include("initial.jl")
include("canForm.jl")
include("TransferM.jl")
include("Heff.jl")
include("sweepSchemes.jl")
include("essentials.jl")

end # module

max_stars_repo_name: Gawatz/finiteMPS

max_issues_repo_name: Gawatz/finiteMPS

max_forks_repo_name: Gawatz/finiteMPS


text:
"""
    NegativeBinomialExpert(n, p)

PMF:

```math
P(X = k) = \\frac{\\Gamma(k+r)}{k! \\Gamma(r)} p^r (1 - p)^k, \\quad \\text{for } k = 0,1,2,\\ldots.
```

See also: [Negative Binomial Distribution](https://reference.wolfram.com/language/ref/NegativeBinomialDistribution.html) (Wolfram)

"""
struct NegativeBinomialExpert{T<:Real} <: NonZIDiscreteExpert
    n::T
    p::T
    NegativeBinomialExpert{T}(n, p) where {T<:Real} = new{T}(n, p)
end

function NegativeBinomialExpert(n::T, p::T; check_args=true) where {T <: Real}
    check_args && @check_args(NegativeBinomialExpert, 0 <= p <= 1 && n > zero(n))
    return NegativeBinomialExpert{T}(n, p)
end

## Outer constructors
NegativeBinomialExpert(n::Real, p::Real) = NegativeBinomialExpert(promote(n, p)...)
NegativeBinomialExpert(n::Integer, p::Integer) = NegativeBinomialExpert(float(n), float(p))
NegativeBinomialExpert() = NegativeBinomialExpert(1, 0.5)

## Conversion
function convert(::Type{NegativeBinomialExpert{T}}, n::S, p::S) where {T <: Real, S <: Real}
    NegativeBinomialExpert(T(n), T(p))
end
function convert(::Type{NegativeBinomialExpert{T}}, d::NegativeBinomialExpert{S}) where {T <: Real, S <: Real}
    NegativeBinomialExpert(T(d.n), T(d.p), check_args=false)
end
copy(d::NegativeBinomialExpert) = NegativeBinomialExpert(d.n, d.p, check_args=false)

## Loglikelihood of Expoert
logpdf(d::NegativeBinomialExpert, x...) = isinf(x...) ? -Inf : Distributions.logpdf.(Distributions.NegativeBinomial(d.n, d.p), x...)
pdf(d::NegativeBinomialExpert, x...) = isinf(x...) ? 0.0 : Distributions.pdf.(Distributions.NegativeBinomial(d.n, d.p), x...)
logcdf(d::NegativeBinomialExpert, x...) = isinf(x...) ? 0.0 : Distributions.logcdf.(Distributions.NegativeBinomial(d.n, d.p), x...)
cdf(d::NegativeBinomialExpert, x...) = isinf(x...) ? 1.0 : Distributions.cdf.(Distributions.NegativeBinomial(d.n, d.p), x...)

## expert_ll, etc
expert_ll_exact(d::NegativeBinomialExpert, x::Real) = LRMoE.logpdf(LRMoE.NegativeBinomialExpert(d.n, d.p), x)
function expert_ll(d::NegativeBinomialExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    # d_exp = NegativeBinomialExpert(d.n*exposure, d.p)
    expert_ll = (yl == yu) ? logpdf.(d, yl) : logcdf.(d, yu) + log1mexp.(logcdf.(d, ceil.(yl) .- 1) - logcdf.(d, yu))
    return expert_ll
end
function expert_tn(d::NegativeBinomialExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    # d_exp = NegativeBinomialExpert(d.n*exposure, d.p)
    expert_tn = (tl == tu) ? logpdf.(d, tl) : logcdf.(d, tu) + log1mexp.(logcdf.(d, ceil.(tl) .- 1) - logcdf.(d, tu))
    return expert_tn
end
function expert_tn_bar(d::NegativeBinomialExpert, tl::Real, yl::Real, yu::Real, tu::Real)
    # d_exp = NegativeBinomialExpert(d.n*exposure, d.p)
    expert_tn_bar = (tl == tu) ? log1mexp.(logpdf.(d, tl)) : log1mexp.(logcdf.(d, tu) + log1mexp.(logcdf.(d, ceil.(tl) .- 1) - logcdf.(d, tu)))
    return expert_tn_bar
end

exposurize_expert(d::NegativeBinomialExpert; exposure = 1) = NegativeBinomialExpert(d.n*exposure, d.p)

## Parameters
params(d::NegativeBinomialExpert) = (d.n, d.p)
function params_init(y, d::NegativeBinomialExpert)
    μ, σ2 = mean(y), var(y)
    p_init = μ / σ2
    n_init = μ*p_init/(1-p_init)
    try
        NegativeBinomialExpert(n_init, p_init)
    catch;
        NegativeBinomialExpert()
    end
end

## Simululation
sim_expert(d::NegativeBinomialExpert) = Distributions.rand(Distributions.NegativeBinomial(d.n, d.p), 1)[1]

## penalty
penalty_init(d::NegativeBinomialExpert) = [2.0 10.0]
no_penalty_init(d::NegativeBinomialExpert) = [1.0 Inf]
penalize(d::NegativeBinomialExpert, p) = (p[1]-1)*log(d.n) - d.n/p[2]

## statistics
mean(d::NegativeBinomialExpert) = mean(Distributions.NegativeBinomial(d.n, d.p))
var(d::NegativeBinomialExpert) = var(Distributions.NegativeBinomial(d.n, d.p))
quantile(d::NegativeBinomialExpert, p) = quantile(Distributions.NegativeBinomial(d.n, d.p), p)

## Misc functions for E-Step

function _sum_densy_series(d::NegativeBinomialExpert, yl, yu)
    if isinf(yu)
        series = 0:(max(yl-1, 0))
        return d.n*(1-d.p)/d.p - sum(pdf.(d, series) .* series)[1]
    else
        series = yl:yu
        return sum(pdf.(d, series) .* series)[1]
    end
end

function _int_obs_Y_raw(d::NegativeBinomialExpert, yl, yu)
    return _sum_densy_series(d, yl, yu)
end

function _int_lat_Y_raw(d::NegativeBinomialExpert, tl, tu)
    return d.n*(1-d.p)/d.p - _sum_densy_series(d, tl, tu)
end

function _sum_denslogy_series(n_new, d::NegativeBinomialExpert, yl, yu, exposure)
    upper_finite = isinf(yu) ? quantile(d, 1-1e-8) : yu
    series = yl:(max(yl, min(yu, upper_finite+1)))
    # series = yl:yu
    return sum(pdf.(d, series) .* loggamma.(series .+ n_new*exposure))[1]
end

function _int_obs_logY_raw(n_new, d::NegativeBinomialExpert, yl, yu, exposure)
    if yl == yu
        return pdf.(d, yl) .* loggamma.(yl .+ n_new*exposure)
    else
        return _sum_denslogy_series.(n_new, d, yl, yu, exposure)
    end
end

function _int_lat_logY_raw(n_new, d::NegativeBinomialExpert, tl, tu, exposure)
    # return (tl==0 ? 0.0 : _sum_denslogy_series.(n_new, d, 0, tl)) + (isinf(tu) ? 0.0 : _sum_denslogy_series.(n_new, d, tu, Inf))
    return (tl==0 ? 0.0 : _sum_denslogy_series.(n_new, d, 0, ceil(tl)-1, exposure)) +
        (isinf(tu) ? 0.0 : _sum_denslogy_series.(n_new, d, floor(tu)+1, Inf, exposure))
end

function _negativebinomial_n_to_p(n, sum_term_zkz, sum_term_zkzy;
                        penalty = true, pen_pararms_jk = [])
    return 1 - sum_term_zkzy/(sum_term_zkz + sum_term_zkzy)
end

function _negativebinomial_optim_n(logn,
                        d_old,
                        tl, yl, yu, tu,
                        exposure,
                        z_e_obs, z_e_lat, k_e,
                        Y_e_obs, Y_e_lat;
                        penalty = true, pen_pararms_jk = [])
    # Optimization in log scale for unconstrained computation
    n_tmp = exp(logn)

    # Further E-Step
    # yl_yu_unique = unique_bounds(yl, yu)
    # int_obs_logY_tmp = _int_obs_logY_raw.(n_tmp, d_old, yl_yu_unique[:,1], yl_yu_unique[:,2])
    # logY_e_obs = exp.(-expert_ll_pos) .* int_obs_logY_tmp[match_unique_bounds(hcat(vec(yl), vec(yu)), yl_yu_unique)]
    # nan2num(logY_e_obs, 0.0) # get rid of NaN

    # tl_tu_unique = unique_bounds(tl, tu)
    # int_lat_logY_tmp = _int_lat_logY_raw.(n_tmp, d_old, tl_tu_unique[:,1], tl_tu_unique[:,2])
    # logY_e_lat = exp.(-expert_tn_bar_pos) .* int_lat_logY_tmp[match_unique_bounds(hcat(vec(tl), vec(tu)), tl_tu_unique)]
    # nan2num(logY_e_lat, 0.0) # get rid of NaN

    logY_e_obs = fill(0.0, length(yl))
    logY_e_lat = fill(0.0, length(yl))

    for i in 1:length(yl)
        d_expo = exposurize_expert(d_old, exposure = exposure[i])
        expert_ll_pos = expert_ll.(d_expo, tl[i], yl[i], yu[i], tu[i])
        expert_tn_bar_pos = expert_tn_bar.(d_expo, tl[i], yl[i], yu[i], tu[i])
        logY_e_obs[i] = exp(-expert_ll_pos) * _int_obs_logY_raw(n_tmp, d_expo, yl[i], yu[i], exposure[i])
        logY_e_lat[i] = exp(-expert_tn_bar_pos) * _int_lat_logY_raw(n_tmp, d_expo, tl[i], tu[i], exposure[i])
    end

    nan2num(logY_e_obs, 0.0) # get rid of NaN
    nan2num(logY_e_lat, 0.0) # get rid of NaN
    inf2num(logY_e_obs, 0.0) # get rid of inf
    inf2num(logY_e_lat, 0.0) # get rid of inf

    term_zkz = (z_e_obs .* n_tmp .* exposure) .+ (z_e_lat .* k_e .* n_tmp .* exposure)
    term_zkz_Y = (z_e_obs .* Y_e_obs) .+ (z_e_lat .* k_e .* Y_e_lat)
    term_zkz_logY = (z_e_obs .* logY_e_obs) .+ (z_e_lat .* k_e .* logY_e_lat)

    sum_term_zkz = sum(term_zkz)[1]
    sum_term_zkzy = sum(term_zkz_Y)[1]
    sum_term_zkzlogy = sum(term_zkz_logY)[1]

    p_tmp = _negativebinomial_n_to_p(n_tmp, sum_term_zkz, sum_term_zkzy, penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    obj = sum_term_zkzlogy - sum(z_e_obs .* loggamma.(n_tmp .* exposure))[1] + sum_term_zkz*log(p_tmp) + sum_term_zkzy*log(1-p_tmp)
    p = penalty ? (pen_pararms_jk[1]-1)*log(n_tmp) - n_tmp/pen_pararms_jk[2] : 0.0
    return (obj + p) * (-1.0)
end

## EM: M-Step
function EM_M_expert(d::NegativeBinomialExpert,
                    tl, yl, yu, tu,
                    exposure,
                    z_e_obs, z_e_lat, k_e;
                    penalty = true, pen_pararms_jk = [2.0 1.0])

    # Further E-Step
    # yl_yu_unique = unique_bounds(yl, yu)
    # int_obs_Y_tmp = _int_obs_Y_raw.(d, yl_yu_unique[:,1], yl_yu_unique[:,2])
    # Y_e_obs = exp.(-expert_ll_pos) .* int_obs_Y_tmp[match_unique_bounds(hcat(vec(yl), vec(yu)), yl_yu_unique)]
    # nan2num(Y_e_obs, 0.0) # get rid of NaN

    # tl_tu_unique = unique_bounds(tl, tu)
    # int_lat_Y_tmp = _int_lat_Y_raw.(d, tl_tu_unique[:,1], tl_tu_unique[:,2])
    # Y_e_lat = exp.(-expert_tn_bar_pos) .* int_lat_Y_tmp[match_unique_bounds(hcat(vec(tl), vec(tu)), tl_tu_unique)]
    # nan2num(Y_e_lat, 0.0) # get rid of NaN
    Y_e_obs = fill(0.0, length(yl))
    Y_e_lat = fill(0.0, length(yl))

    for i in 1:length(yl)
        d_expo = exposurize_expert(d, exposure = exposure[i])
        expert_ll_pos = expert_ll.(d_expo, tl[i], yl[i], yu[i], tu[i])
        expert_tn_bar_pos = expert_tn_bar.(d_expo, tl[i], yl[i], yu[i], tu[i])
        Y_e_obs[i] = exp(-expert_ll_pos) * _int_obs_Y_raw(d_expo, yl[i], yu[i])
        Y_e_lat[i] = exp(-expert_tn_bar_pos) * _int_lat_Y_raw(d_expo, tl[i], tu[i])
    end

    nan2num(Y_e_obs, 0.0) # get rid of NaN
    nan2num(Y_e_lat, 0.0) # get rid of NaN
    inf2num(Y_e_obs, 0.0) # get rid of inf
    inf2num(Y_e_lat, 0.0) # get rid of inf

    # Update parameters
    logn_new = Optim.minimizer( Optim.optimize(x -> _negativebinomial_optim_n(x, d,
                                                tl, yl, yu, tu,
                                                exposure,
                                                z_e_obs, z_e_lat, k_e,
                                                Y_e_obs, Y_e_lat,
                                                penalty = penalty, pen_pararms_jk = pen_pararms_jk),
                                                log(d.n)-0.5, log(d.n)+0.5 )) # ,
                                                # max(log(d.n)-2.0, 0.0), log(d.n)+2.0 )) # ,
                                                # log(d.k)-2.0, log(d.k)+2.0 )) # ,
                                                # GoldenSection() ))
                                                # , rel_tol = 1e-8) )
    n_new = exp(logn_new)

    term_zkz = (z_e_obs .* n_new .* exposure ) .+ (z_e_lat .* k_e  .* n_new .* exposure)
    term_zkz_Y = (z_e_obs .* Y_e_obs) .+ (z_e_lat .* k_e .* Y_e_lat)

    sum_term_zkz = sum(term_zkz)[1]
    sum_term_zkzy = sum(term_zkz_Y)[1]

    p_new = _negativebinomial_n_to_p(n_new, sum_term_zkz, sum_term_zkzy, penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    # Deal with zero mass
    if (p_new^n_new > 0.999999) || (isnan(p_new^n_new))
        p_new, n_new = d.p, d.n
    end

    return NegativeBinomialExpert(n_new, p_new)

end

## EM: M-Step, exact observations
function _negativebinomial_optim_n_exact(logn, d_old,
                        ye, exposure,
                        z_e_obs;
                        penalty = true, pen_pararms_jk = [])
    # Optimization in log scale for unconstrained computation
    n_tmp = exp(logn)

    Y_e_obs = ye

    # Further E-Step
    logY_e_obs = loggamma.(ye .+ n_tmp .* exposure)
    nan2num(logY_e_obs, 0.0) # get rid of NaN

    term_zkz = z_e_obs .* n_tmp .* exposure
    term_zkz_Y = z_e_obs .* Y_e_obs
    term_zkz_logY = z_e_obs .* logY_e_obs

    sum_term_zkz = sum(term_zkz)[1]
    sum_term_zkzy = sum(term_zkz_Y)[1]
    sum_term_zkzlogy = sum(term_zkz_logY)[1]

    p_tmp = _negativebinomial_n_to_p(n_tmp, sum_term_zkz, sum_term_zkzy, penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    obj = sum_term_zkzlogy - sum(z_e_obs .* loggamma.(n_tmp .* exposure))[1] + sum_term_zkz*log(p_tmp) + sum_term_zkzy*log(1-p_tmp)
    p = penalty ? (pen_pararms_jk[1]-1)*log(n_tmp) - n_tmp/pen_pararms_jk[2] : 0.0
    return (obj + p) * (-1.0)
end

function EM_M_expert_exact(d::NegativeBinomialExpert,
                    ye, exposure,
                    # expert_ll_pos,
                    z_e_obs;
                    penalty = true, pen_pararms_jk = [1.0 Inf])

    # Update parameters
    logn_new = Optim.minimizer( Optim.optimize(x -> _negativebinomial_optim_n_exact(x, d,
                                                ye, exposure,
                                                z_e_obs,
                                                penalty = penalty, pen_pararms_jk = pen_pararms_jk),
                                                log(d.n)-0.5, log(d.n)+0.5 ))
    n_new = exp(logn_new)

    term_zkz = z_e_obs .* n_new .* exposure
    term_zkz_Y = z_e_obs .* ye

    sum_term_zkz = sum(term_zkz)[1]
    sum_term_zkzy = sum(term_zkz_Y)[1]

    p_new = _negativebinomial_n_to_p(n_new, sum_term_zkz, sum_term_zkzy, penalty = penalty, pen_pararms_jk = pen_pararms_jk)

    # Deal with zero mass
    if (p_new^n_new > 0.999999) || (isnan(p_new^n_new))
        p_new, n_new = d.p, d.n
    end

    return NegativeBinomialExpert(n_new, p_new)

end

max_stars_repo_name: sparktseung/JuliaLRMoE

max_issues_repo_name: sparktseung/JuliaLRMoE

max_forks_repo_name: sparktseung/JuliaLRMoE


text:
module JsonLinkedData

function ActiveContext(json::Dict{String, Any}, base_iri::String)
  ActiveContext(json, base_iri)
end

# Need some types for algorithm implementation
struct TermDefinition
  iri_mapping::Dict{String, Any}
  reverse_property::Bool
  type_mapping::Nullable{} # What sort of type?
  language_mapping::Nullable{} # What sort of type?
  context::Nullable{} # What is this for? What type?
  prefix_flag::Bool # What is this for?
  container_mapping::Dict{String, Any} # What is this for?
end

# Lots of type enforcement stuff can be handled implicitly
struct ActiveContext
  term_definitions::Dict{String, TermDefinition}
  base_iri::String
  vocabulary_mapping::Dict{String, Any}
  default_language::String
end

jsonld_keywords = Set{String}([
  ":",
  "@base",
  "@container",
  "@context",
  "@graph",
  "@id",
  "@index",
  "@language",
  "@list",
  "@nest",
  "@none",
  "@prefix",
  "@reverse",
  "@set",
  "@type",
  "@value",
  "@version",
  "@vocab"
])

# package code goes here
include("iri_expansion.jl")
include("create_term_definition.jl")
include("context_processing.jl")

end # module


text:
@testset "MOI_wrapper.jl" begin
    for opt in dual_linear_optimizer
        linear_config = MOIT.TestConfig(atol = 1e-6, rtol = 1e-6)
        linear_cache =
            MOIU.UniversalFallback(Dualization.DualizableModel{Float64}())
        MOI.empty!(opt)
        linear_cached = MOIU.CachingOptimizer(linear_cache, opt)
        linear_bridged = MOIB.full_bridge_optimizer(linear_cached, Float64)

        @testset "linear test" begin
            MOIT.contlineartest(
                linear_bridged,
                linear_config,
                [
                    "linear8b", # Asks for infeasibility ray
                    "linear8c", # Asks for infeasibility ray
                    "linear12", # Asks for infeasibility ray
                    "linear13", # Feasibility problem
                    "linear15",  # Feasibility when written in the canonical form
                ],
            )
        end
    end

    for opt in dual_conic_optimizer
        conic_config = MOIT.TestConfig(atol = 1e-4, rtol = 1e-4)
        conic_cache =
            MOIU.UniversalFallback(Dualization.DualizableModel{Float64}())
        conic_cached = MOIU.CachingOptimizer(conic_cache, opt)
        conic_bridged = MOIB.full_bridge_optimizer(conic_cached, Float64)

        @testset "conic linear, soc, rsoc and sdp test" begin
            MOIT.contconictest(
                conic_bridged,
                conic_config,
                [
                    "lin3", # Feasibility problem
                    "lin4", # Feasibility problem
                    "geomean3f",
                    "geomean3v", # CSDP does not converge after https://github.com/jump-dev/Dualization.jl/pull/86
                    "normone2", # Feasibility problem
                    "norminf2", # Feasibility problem
                    "soc3", # Feasibility problem
                    "rotatedsoc2", # Feasibility problem
                    "rotatedsoc4", # Termination status is `ALMOST_OPTIMAL`
                    "exp", # Tested in exp and power cone test
                    "dualexp", # Tested in exp and power cone test
                    "pow", # Tested in exp and power cone test
                    "dualpow", # Tested in exp and power cone test
                    "rootdet", # Dual not defined in MOI
                    "logdet", # Dual not defined in MOI
                    "relentr", # Dual not defined in MOI
                    "psdt3", # atol should be 1e-2
                    "psds3", # atol should be 1e-2
                ],
            )
        end

        @testset "quadratically constrained" begin
            MOIT.contquadratictest(
                conic_bridged,
                conic_config,
                ["qp", "ncqcp", "socp"],
            )
        end
    end

    for opt in dual_power_cone_optimizer
        power_cone_config = MOIT.TestConfig(atol = 1e-3, rtol = 1e-3)
        power_cone_cache =
            MOIU.UniversalFallback(Dualization.DualizableModel{Float64}())
        power_cone_cached = MOIU.CachingOptimizer(power_cone_cache, opt)
        power_cone_bridged =
            MOIB.full_bridge_optimizer(power_cone_cached, Float64)

        @testset "power cone test" begin
            MOIT.contconictest(
                power_cone_bridged,
                power_cone_config,
                [
                    "lin", # Tested in coninc linear, soc, rsoc and sdp test
                    "normone", # Tested in coninc linear, soc, rsoc and sdp test
                    "norminf", # Tested in coninc linear, soc, rsoc and sdp test
                    "soc", # Tested in coninc linear, soc, rsoc and sdp test
                    "rsoc", # Tested in coninc linear, soc, rsoc and sdp test
                    "geomean", # Tested in coninc linear, soc, rsoc and sdp test
                    "sdp", # Tested in coninc linear, soc, rsoc and sdp test
                    "rootdet", # Dual not defined in MOI
                    "logdet", # Dual not defined in MOI
                    "relentr", # Dual not defined in MOI
                ],
            )
        end
    end

    @testset "attributes" begin
        for optimizer in [dual_conic_optimizer; dual_linear_optimizer]
            before = MOI.get(optimizer, MOI.Silent())
            MOI.set(optimizer, MOI.Silent(), !before)
            @test MOI.get(optimizer, MOI.Silent()) == !before
            MOI.set(optimizer, MOI.Silent(), before)
            @test MOI.get(optimizer, MOI.Silent()) == before
        end
        for i in eachindex(dual_conic_optimizer)
            @test MOI.get(dual_conic_optimizer[i], MOI.SolverName()) ==
                  "Dual model with $(MOI.get(primal_conic_optimizer[i], MOI.SolverName())) attached"
        end
        for i in eachindex(dual_linear_optimizer)
            @test MOI.get(dual_linear_optimizer[i], MOI.SolverName()) ==
                  "Dual model with $(MOI.get(primal_linear_optimizer[i], MOI.SolverName())) attached"
        end
    end

    @testset "support" begin
        for opt in dual_linear_optimizer
            @test !MOI.supports_constraint(opt, SVF, MOI.Integer)
            @test MOI.supports(opt, MOI.ObjectiveSense())
        end
        for opt in dual_conic_optimizer
            @test MOI.supports_constraint(
                opt,
                VVF,
                MOI.PositiveSemidefiniteConeTriangle,
            )
        end
    end

    @testset "dual_status" begin
        @test Dualization.dual_status(MOI.INFEASIBLE) == MOI.DUAL_INFEASIBLE
        @test Dualization.dual_status(MOI.DUAL_INFEASIBLE) == MOI.INFEASIBLE
        @test Dualization.dual_status(MOI.ALMOST_INFEASIBLE) ==
              MOI.ALMOST_DUAL_INFEASIBLE
        @test Dualization.dual_status(MOI.ALMOST_DUAL_INFEASIBLE) ==
              MOI.ALMOST_INFEASIBLE
    end

    @testset "DualOptimizer" begin
        for opt in primal_linear_optimizer
            err = ErrorException("DualOptimizer must have a solver attached")
            @test_throws err DualOptimizer()
            dual_opt_f32 = Dualization.DualOptimizer{Float32}(opt)
            Caching_OptimizerType = typeof(dual_opt_f32.dual_problem.dual_model)
            @test typeof(dual_opt_f32) ==
                  DualOptimizer{Float32,Caching_OptimizerType}
        end
    end
end


text:
# Check overloaded access functions
m = DSGEVAR(AnSchorfheide(), [:z_sh])

@testset "Access functions" begin
    # Number of anticipated policy shocks
    @test n_anticipated_shocks(m) == get_setting(m, :n_anticipated_shocks)
    @test n_anticipated_shocks_padding(m) == get_setting(m, :n_anticipated_shocks_padding)

    # Dates, indices, number of periods for each regime
    @test date_presample_start(m) == get_setting(m, :date_presample_start)
    @test date_mainsample_start(m) == get_setting(m, :date_mainsample_start)
    @test date_zlb_start(m) == get_setting(m, :date_zlb_start)
    # date_zlb_end(m) == get_setting(m, :date_zlb_end)

    @test date_presample_end(m) == Dates.lastdayofquarter(get_setting(m, :date_mainsample_start) - Dates.Month(3))
    @test date_prezlb_end(m) == Dates.lastdayofquarter(get_setting(m, :date_zlb_start) - Dates.Month(3))
    @test date_mainsample_end(m) == Dates.lastdayofquarter(get_setting(m, :date_forecast_start) - Dates.Month(3))
    @test date_conditional_end(m) == get_setting(m, :date_conditional_end)

    @test index_presample_start(m) == 1
    @test index_mainsample_start(m) == subtract_quarters(date_mainsample_start(m), date_presample_start(m)) + 1
    @test index_zlb_start(m) == subtract_quarters(date_zlb_start(m), date_presample_start(m)) + 1
    @test index_forecast_start(m) == subtract_quarters(date_forecast_start(m), date_presample_start(m)) + 1

    # Interface for fields and size
    @test DSGE.n_lags(m)        == m.lags
    @test DSGE.get_lags(m)      == m.lags
    @test DSGE.get_shocks(m)    == m.shocks
    @test DSGE.get_λ(m)         == m.λ
    @test DSGE.get_dsge(m)      == m.dsge
    @test DSGE.n_observables(m) == length(m.observables)
    @test DSGE.n_shocks(m)      == length(m.shocks)

    # Interface for accessing parameters
    @test DSGE.get_parameters(m) == m.dsge.parameters
    @test n_parameters(m)   == n_parameters(m.dsge)
    @test n_parameters_steady_state(m) == n_parameters_steady_state(m.dsge)
    @test n_parameters_free(m) == n_parameters_free(m.dsge)

    # Interface for accessing rng
    @test DSGE.get_rng(m) == m.dsge.rng

    # Interface for accessing settings dictionary
    @test DSGE.get_settings(m) == m.dsge.settings

    # Interface for accessing observables dictionary
    @test DSGE.get_observables(m) == m.observables

    # Interface for data
    @test cond_vintage(m)    == get_setting(m, :cond_vintage)
    @test cond_id(m)         == get_setting(m, :cond_id)
    @test cond_full_names(m) == get_setting(m, :cond_full_names)
    @test cond_semi_names(m) == get_setting(m, :cond_semi_names)
    @test use_population_forecast(m) == get_setting(m, :use_population_forecast)
    @test DSGE.hpfilter_population(m)     == get_setting(m, :hpfilter_population)

    # Interface for general computation settings
    @test use_parallel_workers(m)    == get_setting(m, :use_parallel_workers)

    # Interface for estimation settings
    @test reoptimize(m)          == get_setting(m, :reoptimize)
    @test calculate_hessian(m) == get_setting(m, :calculate_hessian)
    @test_throws KeyError DSGE.hessian_path(m)      == get_setting(m, :hessian_path)
    m <= Setting(:hessian_path, "")
    @test DSGE.hessian_path(m)      == get_setting(m, :hessian_path)

    @test n_hessian_test_params(m) == get_setting(m, :n_hessian_test_params)

    # Interface for Metropolis-Hastings settings
    @test n_mh_blocks(m)       == get_setting(m, :n_mh_blocks)
    @test DSGE.n_mh_param_blocks(m) == get_setting(m, :n_mh_param_blocks)
    @test n_mh_simulations(m)  == get_setting(m, :n_mh_simulations)
    @test n_mh_burn(m)         == get_setting(m, :n_mh_burn)
    @test mh_thin(m)           == get_setting(m, :mh_thin)

    # Interface for forecast settings
    @test date_forecast_start(m)   == get_setting(m, :date_forecast_start)
    @test forecast_block_size(m)   == get_setting(m, :forecast_block_size)
    @test forecast_start_block(m)  == get_setting(m, :forecast_start_block)
    @test forecast_input_file_overrides(m) == get_setting(m, :forecast_input_file_overrides)
    @test isnull(forecast_uncertainty_override(m)) # == get_setting(m, :forecast_uncertainty_override)
    @test forecast_smoother(m)     == get_setting(m, :forecast_smoother)
    @test forecast_tdist_df_val(m) == get_setting(m, :forecast_tdist_df_val)
    @test forecast_tdist_shocks(m) == get_setting(m, :forecast_tdist_shocks)
    @test forecast_zlb_value(m)    == get_setting(m, :forecast_zlb_value)
    @test impulse_response_horizons(m) == get_setting(m, :impulse_response_horizons)
    DSGE.get_forecast_input_file(m, :mode) == DSGE.get_forecast_input_file(m.dsge, :mode)

    # Interface for alternative policy settings
    @test alternative_policy(m) == get_setting(m, :alternative_policy)

    @test date_forecast_end(m) == Dates.lastdayofquarter(date_forecast_start(m) + Dates.Month(3 * (forecast_horizons(m)-1)))

    @test DSGE.forecast_horizons(m) == get_setting(m, :forecast_horizons)
    @test DSGE.forecast_horizons(m; cond_type = :full) == get_setting(m, :forecast_horizons) - DSGE.n_conditional_periods(m.dsge)
end

@testset "Updating parameters and settings, transforming parameters" begin
    dsge = Model1002("ss10", testing = true)
    m = DSGEVAR(dsge, [:ztil_sh])
    oldvals = map(x -> x.value, DSGE.get_parameters(m))
    newvals = oldvals + randn(n_parameters(m)) .* 1e-3
    DSGE.update!(dsge, oldvals)
    DSGE.transform_to_model_space!(m, ModelConstructors.transform_to_real_line(DSGE.get_parameters(m),
                                                                               newvals))
    DSGE.transform_to_model_space!(dsge, ModelConstructors.transform_to_real_line(DSGE.get_parameters(dsge),
                                                                                newvals))
    @test map(x -> x.value, DSGE.get_parameters(m)) ≈ map(x -> x.value, DSGE.get_parameters(dsge))
    @test map(x -> x.value, DSGE.get_dsge(m).steady_state) ≈ map(x -> x.value, dsge.steady_state)

    DSGE.update!(dsge, oldvals)
    DSGE.update!(m, oldvals)
    @test map(x -> x.value, DSGE.get_parameters(m)) ≈ map(x -> x.value, DSGE.get_parameters(dsge))
    @test map(x -> x.value, DSGE.get_dsge(m).steady_state) ≈ map(x -> x.value, dsge.steady_state)

    newvals = oldvals + randn(n_parameters(m)) .* 1e-3
    oldpvec = dsge.parameters
    DSGE.update!(dsge, newvals)
    DSGE.update!(m, newvals)
    DSGE.update!(dsge, oldpvec)
    DSGE.update!(m, oldpvec)
    @test map(x -> x.value, DSGE.get_parameters(m)) ≈ map(x -> x.value, DSGE.get_parameters(dsge))
    @test map(x -> x.value, DSGE.get_dsge(m).steady_state) ≈ map(x -> x.value, dsge.steady_state)

    @test DSGE.get_parameters(m)[DSGE.get_dsge(m).keys[:α]].value != 0.5
    m <= parameter(:α, 0.5)
    @test DSGE.get_parameters(m)[DSGE.get_dsge(m).keys[:α]].value == 0.5

    @test DSGE.get_dsge(m).steady_state[DSGE.get_dsge(m).keys[:z_star] - n_parameters(m)].value != 0.5
    m <= SteadyStateParameter(:z_star, 0.5)
    @test DSGE.get_dsge(m).steady_state[DSGE.get_dsge(m).keys[:z_star] - n_parameters(m)].value == 0.5

    @test !haskey(DSGE.get_settings(m), :test_setting) && !haskey(DSGE.get_dsge(m).test_settings, :test_setting)
    m <= Setting(:test_setting, true)
    @test get_setting(m, :test_setting) && haskey(DSGE.get_dsge(m).test_settings, :test_setting)
    dsge = Model1002("ss10")
    m = DSGEVAR(dsge, [:ztil_sh])
    @test !haskey(DSGE.get_settings(m), :test_setting) && !haskey(DSGE.get_dsge(m).test_settings, :test_setting)
    m <= Setting(:test_setting, true)
    @test get_setting(m, :test_setting) && haskey(DSGE.get_settings(m), :test_setting)

    @test prior(m) == prior(DSGE.get_dsge(m))
end

@testset "File paths" begin
    m = DSGEVAR(AnSchorfheide("ss0", testing = true), [:z_sh])
    for folder in ["raw", "work", "figures", "tables", "log"]
        @test eval(Symbol(folder, "path"))(m, "estimate") ==
            joinpath(saveroot(m), "output_data", "an_schorfheide", "ss0", "dsgevar_an_schorfheide", "ss0", "estimate", folder)
    end
end


text:
using NbodyGradient, LinearAlgebra

export State

#include("../src/ttv.jl")
#include("/Users/ericagol/Software/TRAPPIST1_Spitzer/src/NbodyGradient/src/ttv.jl")
# Specify the initial conditions for the outer solar
# system
#n=6
n=5
xout = zeros(3,n)
# Positions at time September 5, 1994 at 0h00 in days (from Hairer, Lubich & Wanner
# 2006, Geometric Numerical Integration, 2nd Edition, Springer, pp. 13-14):

xout .= transpose([-2.079997415328555E-04  7.127853194812450E-03 -1.352450694676177E-05;
        -3.502576700516146E+00 -4.111754741095586E+00  9.546978009906396E-02;
         9.075323061767737E+00 -3.443060862268533E+00 -3.008002403885198E-01;
         8.309900066449559E+00 -1.782348877489204E+01 -1.738826162402036E-01;
         1.147049510166812E+01 -2.790203169301273E+01  3.102324955757055E-01]) #;
        #-1.553841709421204E+01 -2.440295115792555E+01  7.105854443660053E+00])
vout = zeros(3,n)
vout .= transpose([-6.227982601533108E-06  2.641634501527718E-06  1.564697381040213E-07;
         5.647185656190083E-03 -4.540768041260330E-03 -1.077099720398784E-04;
         1.677252499111402E-03  5.205044577942047E-03 -1.577215030049337E-04;
         3.535508197097127E-03  1.479452678720917E-03 -4.019422185567764E-05;
         2.882592399188369E-03  1.211095412047072E-03 -9.118527716949448E-05]) #;
         #2.754640676017983E-03 -2.105690992946069E-03 -5.607958889969929E-04]);
# Units of velocity are AU/day

# Specify masses, including terrestrial planets in the Sun:
m = [1.00000597682,0.000954786104043,0.000285583733151,
         0.0000437273164546,0.0000517759138449] #,6.58086572e-9];

# Compute the center-of-mass:
vcm = zeros(3);
xcm = zeros(3);
for j=1:n
    vcm .+= m[j]*vout[:,j];
    xcm .+= m[j]*xout[:,j];
end
vcm ./= sum(m);
xcm ./= sum(m);
# Adjust so CoM is stationary
for j=1:n
    vout[:,j] .-= vcm[:];
    xout[:,j] .-= xcm[:];
end

struct CartesianElements{T} <: NbodyGradient.InitialConditions{T}
    x::Matrix{T}
    v::Matrix{T}
    m::Vector{T}
    nbody::Int64
end

ic = CartesianElements(xout,vout,m,5);

function NbodyGradient.State(ic::NbodyGradient.InitialConditions{T}) where T<:AbstractFloat
    n = ic.nbody
    x = copy(ic.x)
    v = copy(ic.v)
    jac_init = zeros(7*n,7*n)
    xerror = zeros(T,size(x))
    verror = zeros(T,size(v))
    jac_step = Matrix{T}(I,7*n,7*n)
    dqdt = zeros(T,7*n)
    dqdt_error = zeros(T,size(dqdt))
    jac_error = zeros(T,size(jac_step))

    rij = zeros(T,3)
    a = zeros(T,3,n)
    aij = zeros(T,3)
    x0 = zeros(T,3)
    v0 = zeros(T,3)
    input = zeros(T,8)
    delxv = zeros(T,6)
    rtmp = zeros(T,3)
    return State(x,v,[0.0],copy(ic.m),jac_step,dqdt,jac_init,xerror,verror,dqdt_error,jac_error,n,
    rij,a,aij,x0,v0,input,delxv,rtmp)
end

function compute_energy(m::Array{T,1},x::Array{T,2},v::Array{T,2},n::Int64) where {T <: Real}
  KE = 0.0
  for j=1:n
    KE += 0.5*m[j]*(v[1,j]^2+v[2,j]^2+v[3,j]^2)
  end
  PE = 0.0
  for j=1:n-1
    for k=j+1:n
       PE += -NbodyGradient.GNEWT*m[j]*m[k]/norm(x[:,j] .- x[:,k])
    end
  end
  ang_mom = zeros(3)
  for j=1:n
    ang_mom .+= m[j]*cross(x[:,j],v[:,j])
  end
  return KE,PE,ang_mom
end

# Now, integrate this forward in time:
hgrid = [1.5625,3.125,6.25,12.5,25.0,50.0,100.0,200.0]
nstepgrid = [32000000,16000000,8000000,4000000,2000000,1000000,500000,250000]
nskip = [128,64,32,16,8,4,2,1]

#h = 200.0 # 200-day time-step chosen to be <1/20 of the orbital period of Jupiter
#h = 100.0 # 100-day time-step chosen to check conservation of energy/angular momentum with time step
#h = 50.0 # 50-day time-step chosen to check conservation of energy/angular momentum with time step
h = 25.0 # 25-day time-step chosen to check conservation of energy/angular momentum with time step
#h = 12.5 # 12.5-day time-step chosen to check conservation of energy/angular momentum with time step
#h = 6.25 # 6.25-day time-step chosen to check conservation of energy/angular momentum with time step
#h = 3.125 # 3.125-day time-step chosen to check conservation of energy/angular momentum with time step
#h = 1.5625 # 1.5625-day time-step chosen to check conservation of energy/angular momentum with time step

# 50 days x 1e6 time steps ~ 137,000 yr (takes about 15 seconds to run)
grad = false
nstep = maximum(nstepgrid)
ngrid = 8
xsave = zeros(3,n,nstep,ngrid)
vsave = zeros(3,n,nstep,ngrid)
PE = zeros(nstep,ngrid); KE=zeros(nstep,ngrid); ang_mom = zeros(3,nstep,ngrid)
telapse = zeros(ngrid)
for j=1:8
  h = hgrid[j]
  nstep = nstepgrid[j];
  s = State(ic)
  pair = zeros(Bool,s.n,s.n)
  if grad; d = Derivatives(T,s.n); end

  # Set up array to save the state as a function of time:
  # Save the potential & kinetic energy, as well as angular momentum:
  # Time the integration:
  tstart = time()
  # Carry out the integration:
  for i=1:nstep
    if grad
      ahl21!(s,d,h,pair)
    else
      ahl21!(s,h,pair)
    end
    xsave[:,:,i,j] .= s.x
    vsave[:,:,i,j] .= s.v
    KE_step,PE_step,ang_mom_step=compute_energy(s.m,s.x,s.v,n)
    KE[i,j] = KE_step
    PE[i,j] = PE_step
    ang_mom[:,i,j] = ang_mom_step
  end
  s.t[1] = h*nstep
  telapse[j] = time()- tstart
  println("h: ",h," nstep: ",nstep," time: ",telapse[j])
end

using PyPlot,Statistics

for j=8:-1:1
  plot(collect(1:1:nstepgrid[j])*hgrid[j],KE[1:nstepgrid[j],j] .+ PE[1:nstepgrid[j],j]); println(j," ",hgrid[j]," ",std(KE[1:nstepgrid[j],j] .+ PE[1:nstepgrid[j],j])); Emean[j] = mean(KE[1:nstepgrid[j],j] .+ PE[1:nst
end

read(stdin,Char)

clf()
loglog(hgrid,(Emean .- minimum(Emean)) ./abs(minimum(Emean)) ,"o")
loglog(hgrid,(Emean[8] .- minimum(Emean) ) ./abs(minimum(Emean)) .* (hgrid ./ hgrid[8]).^4 )
xlabel("Time step [d]")
ylabel("Fractional change in mean energy")